$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\bR}{\mathbb{R}}$ $\newcommand{\bC}{\mathbb{C}}$ $\newcommand{\bZ}{\mathbb{Z}}$ $\newcommand{\const}{\operatorname{const}}$ $\newcommand{\sgn}{\operatorname{sgn}}$ $\newcommand{\rank}{\operatorname{rank}}$
We consider Cauchy problem for Burgers equation \begin{align} &u_t+uu_x=\varepsilon u_{xx},\qquad t>0, -\infty < x<\infty \label{eq-9.1.1}\\ &u|_{t=0}=f(x). \label{eq-9.1.2} \end{align} Smooth global solution $u=u(x,t;\varepsilon)$ exists for any $\varepsilon >0$ (we assume that $f(x)$ is a bounded function). We are interested in its asymptotics as $\varepsilon \to +0$.
The theory of this problem with $\varepsilon =0$ \begin{align} &v_t+vv_x=0,\qquad t>0, -\infty < x<\infty \label{eq-9.1.3}\\ &v|_{t=0}=f(x). \label{eq-9.1.4} \end{align} is studied in Section 12.1 of online PDE textbook. Read it!
If $f(x)$ is a monotone-nondecreasing function then solution to (\ref{eq-9.1.3})-(\ref{eq-9.1.2}) is at least as smooth as $f(x)$ itself (furthermore, discontinuities may dissipate). Solution is defined as \begin{equation} x= z+f(z)t,\qquad v=f(z) \label{eq-9.1.5} \end{equation} or equivalently as an implicit function \begin{equation} v=f(x-vt). \label{eq-9.1.6} \end{equation}
However if $f(x)$ is not a monotone-nondecreasing function solution $v(x,t)$ develops jumps. For continuous $v$ solution satisfies \begin{equation} v_t+(\frac{1}{2}v^2)_x=0,\qquad t>0, -\infty < x<\infty \label{eq-9.1.7} \end{equation} which we consider as a main equation as $v$ is not continuous (then it is understood in the sense of distributions; see Chapter 11 of online PDE textbook. It implies that if $x=\xi(t)$ is a jump then \begin{equation} \frac{d\xi}{dt}= \frac{1}{2} \bigl(v(\xi(t)+0,t)+v(\xi(t)-0,t)\bigr). \label{eq-9.1.8} \end{equation}
For continuous solutions also \begin{equation} (\frac{1}{2}v^2)_t+(\frac{1}{3}v^3)_x=0,\qquad t>0, -\infty < x<\infty \label{eq-9.1.9} \end{equation} but for discontinuous \begin{equation} (\frac{1}{2}v^2)_t+(\frac{1}{3}v^3)_x\le 0,\qquad t>0, -\infty < x<\infty \label{eq-9.1.10} \end{equation} where for distributions this inequality is understood in the following sense: Distribution $U \ge 0$ if for all non-negative test functions $\varphi$ $U(\varphi) \ge 0$.
One can prove that solution of problem (\ref{eq-9.1.1})-(\ref{eq-9.1.2}) $u(x,t;\varepsilon)$ converges (in the sense of distributions) to solution of problem (\ref{eq-9.1.1})-(\ref{eq-9.1.2}) $v(x,t)$ as $\varepsilon\to +0$.
We claim that this solution satisfies (\ref{eq-9.1.10}). Indeed, from (\ref{eq-9.1.1}) we conclude that \begin{equation*} (\frac{1}{2}u^2)_t+(\frac{1}{3}u^3)_x= \varepsilon uu_{xx}= \varepsilon (\frac{1}{2}u^2)_{xx} - \varepsilon u_x^2; \end{equation*} applying it to test function $\varphi \ge 0$ we pass to the limit in the left (requires justification) and in the first term on the right (also requires justification, results in $0$) but the last term on the right is $-\varepsilon u_x^2 (\varphi)$ which results in $\le 0$ (but not necessarily equal).
It was proven about 50 y.a. that under additional restriction (\ref{eq-9.1.10}) problem (\ref{eq-9.1.3})--(\ref{eq-9.1.4}) in the sense of distributions has a unique solution (which is not the case without this restriction).
Luckily with Cole-Hopf transform \begin{equation} u= -2\varepsilon (\log \varphi)_x= -2\varepsilon\frac{\varphi_x}{\varphi} \label{eq-9.1.11} \end{equation} with the inverse transform \begin{equation} \varphi (x,t)=\exp \Bigl( -\int_{-\infty}^x u(y,t)\,dy) \label{eq-9.1.12} \end{equation} we reduce (\ref{eq-9.1.1}) to a heat equation for $\varphi$: \begin{equation} \varphi _t =\varepsilon\varphi_{xx}. \label{eq-9.1.13} \end{equation}
Indeed, plugging (\ref{eq-9.1.11}) into the left-hand expression of (\ref{eq-9.1.1}) we get \begin{equation*} \bigl(-2\varepsilon (\log \varphi)_t + 2 \varepsilon^2 (\log \varphi)_x^2\bigr)_x= \bigl(-2\varepsilon \frac{\varphi_t}{\varphi} + 2 \varepsilon^2 \frac{\varphi_x^2}{\varphi^2}\bigr)_x \end{equation*} plugging (\ref{eq-9.1.11}) into the right-hand expression of (\ref{eq-9.1.1}) we get \begin{equation*} \bigl(-2\varepsilon^2 (\log \varphi) \bigr)_{xxx}= \bigl(\bigl(-2\varepsilon^2 (\log \varphi) \bigr)_{xx}\bigr)_x= \bigl(-2\varepsilon^2 \frac{\varphi_{xx}}{\varphi} + 2 \varepsilon^2 \frac{\varphi_x^2}{\varphi^2} \bigr)_x \end{equation*} which follows from (\ref{eq-9.1.13}).
Since we know the formula for a solution of Cauchy problem for heat equation (\ref{eq-9.1.13}) (see (3.1.14) of PDE textbook) \begin{equation} \varphi (x,t;\varepsilon) = \frac{1}{\sqrt{4\pi \varepsilon t}}\int_{-\infty}^{\infty} e^{-\frac{1}{4\varepsilon t}(x-z)^2}\varphi(y,0;\varepsilon)\,dz; \label{eq-9.1.14} \end{equation} then using (\ref{eq-9.1.12}) and (\ref{eq-9.1.2}) we have \begin{equation} \varphi(x,0;\varepsilon)=e^{-\frac{1}{2\varepsilon}F(x)} \label{eq-9.1.15} \end{equation} with \begin{equation} F(x)=\int_0 ^x f(y)\,dy \label{eq-9.1.16} \end{equation} Then plugging into (\ref{eq-9.1.14}) and the result in (\ref{eq-9.1.11}) and cancelling factors $\frac{1}{\sqrt{4\pi \varepsilon t}}$ in both numerator and denominator we get \begin{equation} u(x,t;\varepsilon)= \frac{\int_{-\infty}^{\infty} e^{-\frac{1}{2\varepsilon} \Phi(x,z,t)} \frac{1}{t} (x-z)\,dz} {\int_{-\infty}^{\infty} e^{-\frac{1}{2\varepsilon} \Phi(x,z,t)}\,dz} \label{eq-9.1.17} \end{equation} with \begin{equation} \Phi(x,z,t)=F(z)+\frac{1}{2 t} (x-z)^2. \label{eq-9.1.18} \end{equation}