The fate of quantum shock waves at late times

TL;DR

This paper investigates how quantum shock waves in fermionic many-body systems evolve and decay over time, especially when weak interactions are present, extending understanding from free to interacting systems.

Contribution

It generalizes the study of quantum shock waves from free fermions to include weak interactions, analyzing their decay and shape evolution at late times.

Findings

Shock waves decay and change shape due to interactions

Weak interactions cause shock wave broadening

Results extend understanding of hydrodynamics in quantum systems

Abstract

Shock waves are an ubiquitous feature of hydrodynamic theories. Given that fermionic quantum many-body systems admit hydrodynamical descriptions on length scales large compared to the Fermi wavelength, it is natural to ask what the status of shock waves is in such systems. Free fermions provide a solvable yet non-trivial example, and here we generalise to include generic (non-integrable) weak interactions to understand how a shock wave decays and changes its shape well after forming.

Equations

$$\rho (x)=\frac{k^0}{\pi }+\frac{k^m-k^0}{\pi }f\left(\frac{x}{w}\right),$$

$$H=\sum _p\frac{p^2}{2m}{\psi }_p^{†}{\psi }_p^{},$$

$$\begin{gathered} W(x,k,t)\equiv \int \mathrm{d}y{e}^{-iky} \\ ⟨\Psi |e^{iHt}{\psi }^{†}\left(x-\frac{y}{2}\right)\psi \left(x+\frac{y}{2}\right)e^{-iHt}|\Psi ⟩, \end{gathered}$$

$$\rho (x,t)=\int \frac{\mathrm{d}k}{2\pi }W(x,k,t),$$

$$\left({\partial }_t+\frac{k}{m}{\partial }_x\right)W(x,k,t)=0.$$

$$W_0(x,k,0)\equiv \theta (k_F(x)-k)\theta (k_F(x)+k).$$

$$W_0(x,k,0)\approx \theta (k_F(x)-k)\theta (k^0+k).$$

$$k_F(x,t)=k_F(x-k_F(x,t)t/m,0),$$

$$|\Psi ⟩=e^{i\int \mathrm{d}x{\rho }^R(x)\Phi (x)}|0⟩.$$

$$W(x,k,0)=\int \mathrm{d}y{e}^{i(k^0-k)y}\frac{e^{i(\Phi (x+y/2)-\Phi (x-y/2))}}{2\pi i(y+i{0}^{+})}.$$

$$\frac{\delta \rho (x,t)}{k^m-k^0}$$

$$=\frac{t_s}{t\Delta N^{1/3}}[{\mathrm{Ai}}^{\prime }{\left(\frac{x_{+}(t)-x}{\ell (t)}\right)}^2$$

$$-\frac{x_{+}(t)-x}{\ell (t)}{\mathrm{Ai}}^2\left(\frac{x_{+}(t)-x}{\ell (t)}\right)],$$

$$\ell (t)$$

$$\sim \frac{1}{2^{1/3}{\left(\Delta N\right)}^{2/3}}\frac{t}{t_S}w,$$

$$\lambda (x,t)\equiv \frac{x_{+}(t)-x}{x_{+}(t)-x_{-}(t)}.$$

$$\frac{\delta \rho [\lambda ,t]}{\left(k^m-k^0\right)}\approx \frac{\sqrt{\lambda }}{\pi }\frac{t_S}{t}\left(1+\frac{\sin \left(\frac{2}{3}(\lambda /{\lambda }_{\mathrm{cr}}{)}^{3/2}\right)}{2\lambda /{\lambda }_{\mathrm{cr}}}\right),$$

$$\left({\partial }_t+\frac{k}{m}{\partial }_x\right)f(x,k,t)=I\left\{f\right\},$$

$$\left({\partial }_t+\frac{k}{m}{\partial }_x\right)f(x,k,t)=J(x,k,t)-\Gamma (k)f(x,k,t).$$

$$J(x,k,t)$$

$$\Gamma (k)$$

$$f(x,k,t)=\{\begin{array}{ll}{f}_{\mathrm{low}}(x,k,t),& k^0<k\le k_{-}\\ f_{\mathrm{high}}(x,k,t),& k_{-}<k\le k_{+}\end{array}.$$

$$\left({\partial }_t+\frac{k}{m}{\partial }_x\right)f_{\mathrm{high}}(x,k,t)=-\Gamma (k)f_{\mathrm{high}}(x,k,t).$$

$$\begin{gathered} f_{\mathrm{high}}(x,k,t)= \\ \theta \left(k-k_{-}(x,t)\right)\theta \left(k_{+}(x,t)-k\right)e^{-t\Gamma (k)}. \end{gathered}$$

$$J(x,k,t)$$

$$={\int }_{k_{-}(x,t)}^{k_{+}(x,t)}\mathrm{d}p{e}^{-t\Gamma (p)}{\mathcal{W}}_{p\to k}.$$

$$\begin{gathered} f_{\mathrm{low}}(x,k,t)= \\ {\int }_0^t\mathrm{d}{t}^{\prime }{e}^{-(t-t^{\prime })\Gamma (k)}J\left(x-k(t-t^{\prime }),k,t^{\prime }\right). \end{gathered}$$

$$\begin{gathered} H=\sum _p\frac{p^2}{2m}{\psi }_p^{†}{\psi }_p \\ +\frac{1}{2L}\sum _{\begin{array}{c}q\ne 0\\ \alpha =R,L\end{array}}\left(V_q{\rho }_q^{\alpha }{\rho }_{-q}^{\alpha }+2U_q{\rho }_q^R{\rho }_q^L\right), \end{gathered}$$

$${\mathcal{W}}_{p\to k}$$

$$\frac{\Gamma (p)}{\Gamma (k^m)}$$

$$J(x,k,t)\approx \left(k_{+}(x,t)-k_{-}(x,t)\right)e^{-t\Gamma (\bar{k})}{\mathcal{W}}_{\bar{k}\to k},$$

$$x-k(t-t^{\prime })\le x_{+}(t^{\prime }),$$

$$t^{\prime }\ge t_1\equiv \frac{x-kt}{k^m-k},$$

Full text

I Introduction

The notion of shock waves is well established in classical hydrodynamicswitham . If the fluid velocity is an increasing function of density, then any smooth density profile with a local maximum will eventually form a shock wave: a physical quantity becomes non-analytic as a function of spatial coordinates. In fermionic systems, there is a tension between this singular behaviour and the dispersive broadening one may expect at the level of single-particle quantum mechanics.

The nature of shock waves in the context of free fermions has been the subject of previous theoretical investigation Protopopov ; BettelheimGlazman ; AbanovWiegmann1 ; AbanovWiegmann2 ; Bettelheim . In a particular classical limit, the formation of a shock wave is exhibited as a non-analyticity in the density $\rho(x)$. Semi-classical corrections modify this by smoothing out the behaviour at the shock front through the introduction of quantum ripples. It should be emphasised that prior work has focussed on times close to the formation of the shock, where the ripples may be significant across the entirety of the structure associated with the shock. In Section II, we recapitulate these results and observe that for times well after shock formation a parametrically large spatial region has only negligible quantum corrections to the density.

Recently, the topic of generalised hydrodynamicsBertiniPRL ; DoyonDubail ; DoyonYoshimura ; CastroAlvaredo has led to significant progress in understanding the dynamics of quantities such as the density for integrable systems. The consequences of generic (i.e. non-integrable) interactions are not clear from this picture, however. In this paper we linearise a Boltzmann equation and use single particle decay rates to describe the effect of interactions. This technique allows us to investigate the shock wave at all spatial scales, excluding only a small spatial region affected by quantum corrections, and is valid for times well beyond that of shock formation.

The kinetic theory developed in Section III allows us to find the deformation of the spatial distribution of the density caused by relaxation. Despite the exponential decay of the number of fermions forming the shock wave, a substantial section of the shock wave retains its profile.

We present the final conclusions in Section IV, where we associate the dissolution of the shock wave with the interplay between the quantum-mechanical dispersion and the quasiparticle kinetics.

II Shock waves for free fermions

In order to present a self-contained discussion, we begin by recapitulating and extending some results from Ref. BettelheimGlazman, . The question we wish to address is the following: for a system of spinless fermions, given an initial density profile

[TABLE]

how does the density evolve as a function of time? In Eq. (1) $k^{0}/\pi$ is a background density corresponding to a uniform Fermi sea, and $\left(k^{m}-k^{0}\right)/\pi$ is the height of an isolated, smooth perturbation with profile $f(s)$. This has a single maximum at $f(0)\sim 1$ and $\lim_{|s|\to\infty}f(s)=0$, varying on the scale $s\sim 1$ (the perturbation is of extent $\sim w$). We restrict to the scenario where the height of the perturbation is small $k^{m}-k^{0}\ll k^{0}$, and the number of particles contained in the perturbation is large $\Delta N\sim\left(k^{m}-k^{0}\right)w\gg 1$. $\Delta N$ will be our large parameter for a semi-classical treatment.

The small height of the perturbation implies that excitations are confined to be particles and holes in the vicinity of the Fermi points. This gives us a well-defined notion of right- and left-movers. We initially consider the case of free fermions with a parabolic dispersion relation and mass $m$ as given by the Hamiltonian

[TABLE]

where $\psi^{\dagger}_{p}$ and $\psi^{\phantom{\dagger}}_{p}$ are fermionic creation/annihilation operators at momentum $p$ and obey the standard anti-commutation relations, and we set $\hbar=1$ throughout.

We are interested in a semi-classical description of the problem, and so introduce the Wigner function, defined by

[TABLE]

where $|\Psi\rangle$ is the initial state at $t=0$. The Wigner function is useful for a number of reasons: it allows us to perform a controlled semi-classical approximation with large parameter $\Delta N$, it provides simple access to the density, given by

[TABLE]

and finally, for $H$ given by Eq. (2), $W(x,k,t)$ obeys the simple linear differential equation

[TABLE]

This results in the time-evolved Wigner function having the form $W(x,k,t)=W\left(x-kt/m,k,0\right)$.

II.1 Classical picture

It is a natural ansatz that, due to the smooth variation of the density on the scale of the Fermi wavelength, the Wigner function of Eq. (3) may be described by a “local Fermi surface” i.e.

[TABLE]

Right- and left-movers separate on a timescale $t_{LR}\sim\frac{mw}{k^{0}}$. We therefore choose to ignore left-movers with no loss of generality, and simplify the above to

[TABLE]

This leads to the implicit equation

[TABLE]

which gives rise to multi-valued solutions on a time-scale $t_{S}\sim\frac{mw}{k^{m}-k^{0}}$. This is consistent with ignoring left-movers, as $t_{LR}\ll t_{S}$. The region where $k_{F}(x,t)$ is multi-valued exists between the front of the shock, which we denote $x_{+}(t)$; and the back of the shock $x_{-}(t)$. Formally $x_{\pm}(t)$ are the two solutions of $\partial_{k}x_{F}(k,t)=0$ where $x_{F}(k,t)$ satisfies $k_{F}\left(x_{F}(k,t),t\right)=k$. At $t=t_{S}$ these two solutions coincide. For $t-t_{S}\ll t_{S}$, $x_{+}(t)-x_{-}(t)\sim w\left(\frac{t-t_{S}}{t_{S}}\right)^{2}$. For $t\gg t_{S}$ the difference between them (i.e. the extent of the shock) grows linearly in time as $x_{+}(t)-x_{-}(t)\sim(k^{m}-k^{0})t$.

Between the points $x_{-}(t)$ and $x_{+}(t)$, $k_{F}(x,t)$ has three branches which we will denote $k_{F}^{(0)}(x,t)\leq k_{-}(x,t)\leq k_{+}(x,t)$. It is evident from Fig. 1 that near $x_{+}(t)$ the density acquires square-root behaviour in $x_{+}(t)-x$, and so within the ansatz of Eq. (7) a non-analyticity in the density arises.

II.2 Semi-classical corrections

The main result of Ref. BettelheimGlazman, is to quantify how, for a specific form of initial state $|\Psi\rangle$, including the leading semi-classical correction rounds off the non-analytic behaviour. We begin from the same point, specifying the initial state as

[TABLE]

Here $\rho^{R}(x)$ is the density associated with right-movers, $|0\rangle$ is the (translationally invariant) ground state with Fermi momentum $k^{0}$, and $\Phi(x)$ is a smooth function corresponding to a density $\rho(x)=k^{0}/\pi+\Phi^{\prime}(x)/(2\pi)$ i.e. $\Phi^{\prime}(x)\leftrightarrow k_{F}(x)-k^{0}$. This state has convenient analytic structure, and is experimentally relevant in terms of being preparable by a sudden large perturbationCobdenMuzykantskii ; AbanovWiegmann2 .

Considering only right-movers, standard bosonisation techniques give an explicit integral representation for the Wigner function at $t=0$ of

[TABLE]

Performing a gradient expansion of $\Phi$ in the exponent, it is clear that retaining only the linear-in-$y$ term leads to the step-function approximation of Eq. (7). This approximation is justified in the $\Delta N\to\infty$ limit, where the ansatz of Eq. (7) as describing Eq. (3) is exact. Keeping the $y^{3}$ term in Eq. (10) amounts to including semi-classical corrections.

The excess density may be expressed in terms of the distance from the front of the shock as

[TABLE]

where the length-scale $\ell(t)$ is given by

[TABLE]

and ${\rm Ai}$ is an Airy functionBettelheimGlazman . Note that the previous work focussed on times shortly after the formation of the classical shock ($t-t_{S}\ll t_{S}$, $x_{+}-x_{-}\ll w$), and identified the small parameter $1/\Delta N^{2/3}$ required for the classical description to be valid. The main result of Ref. BettelheimGlazman, is that the scaling function of Eq. (11) gives a good description of the ripples for the entire interval $x_{+}-x_{-}\lesssim w$ until times $t-t_{S}\sim t_{S}$. In fact, as shown in Appendix A, Eq. (11) continues to give a good description at all times $t\gtrsim t_{S}$ for $x_{+}(t)-x\ll x_{+}(t)-x_{-}(t)$.

Because the extent of the shock grows linearly with time for $t\gg t_{S}$, it will be useful to introduce a dimensionless parameter measuring the distance from the front of the shock:

[TABLE]

In terms of this dimensionless variable, we may express the asymptote of Eq. (11) for $\lambda>0$ as

[TABLE]

where the crossover scale $\lambda_{\rm cr}=\frac{1}{(\Delta N)^{2/3}}$. For the regime $1\gg\lambda\gg\lambda_{\rm cr}$, semi-classical corrections are negligible and the simple step-function of Eq. (7) captures the essential physics. In other words, although the spatial window where quantum corrections are appreciable grows with time, it is parametrically small on the length-scale of the shock.

How this window changes over time can be understood simply: the length-scale for quantum corrections $\ell(t)$ grows linearly in time. At times $t-t_{S}\ll t_{S}$, the extent of the shock is small as $x_{+}(t)-x_{-}(t)\sim w\left(t-t_{S}\right)^{2}/{t_{S}}^{2}$ and the quantum corrections are significant. However, at late times $t-t_{S}\gg t_{S}$, while the length-scale $\ell(t)$ grows linearly in time, so too does the extent of the shock and we have that $\frac{x_{+}(t)-x_{-}(t)}{\ell(t)}\sim\Delta N^{2/3}$ i.e. the ripples are squeezed into a fraction $\Delta N^{-2/3}$ of the shock. Therefore quantum corrections are most significant around the time $t=t_{S}$, when the shock nucleates. For $t-t_{S}\gg t_{S}$, the fraction of the shock smeared by quantum fluctuations $\lambda_{\rm cr}$ remains finite and independent of time.

We wish to add small, generically integrability-breaking interactions to this picture. Having established the regime within which semi-classical corrections are small, restricting to this will allow us to make further analytic progress.

III Adding generic interactions

We wish to understand how adding interactions changes the behaviour at late times, restricting to small interactions such that the shock structure of free fermions can become established before decay processes start taking effect. Well after the formation of the overhanging profile of Fig. 1(b), particles above the Fermi surface will begin to relax towards lower energies. This modifies the Wigner function from the step-like behaviour of Fig. 1, and it will generically be non-zero for $k^{0}<k<k_{-}(x,t)$ and $x_{-}(t)<x<x_{+}(t)$. We will examine how to modify the free fermionic description to account for integrability-breaking interactions, and how this changes the evolution of the density as a function of time.

In one dimension, two-particle collisions do not redistribute energy and momentum. Generic interactions permit 3-particle collisions, which leads to the relaxation of excited statesMatveev2 . We denote the decay rate for single-particle excitations over the Fermi sea with momentum $k>k^{0}$ by $\Gamma(k)$. We wish to incorporate this decay rate, with characteristic magnitude $\Gamma(k^{m})$, into our description of the time-evolution of the shock. By working explicitly in the regime where $\Gamma(k^{m})t_{S}\ll 1$ the shock profile is established before decay processes become important. We also require $t\gg t_{S}$, $1\gg\lambda\gg\lambda_{\rm cr}$ such that we may dispense with ripples. Accordingly, one may view the Wigner function as the distribution function in the classical limit, $f(x,k,t)$. In the absence of integrability, three-particle collisions lead to a redistribution of the occupied states, and this is captured by the kinetic (Boltzmann) equationlandafshitzsp ; Matveev1

[TABLE]

where $I$ is a three-particle collision integral.

We wish to evaluate how the shock structure fades on times $t\gtrsim 1/\Gamma(k^{m})$. At these long times, the variation of the spatial structure is smooth on the scale of the Fermi wavelength. Intuitively, interactions will lead to a decay of $f(x,k,t)$ at “high energies” (i.e. for $k$ between $k_{-}(x,t)$ and $k_{+}(x,t)$), which will act as a source for “low energies” ($k$ between $k^{0}$ and $k_{-}(x,t)$). Linearising the collision integral $I$ gives an equation of the form

[TABLE]

Here $J(x,k,t)$ represents “high energy” particles decaying and acting as a source for $f(x,k,t)$, and $\Gamma(k)$ is the aggregate of decay processes from momentum $k$ to lower energies. In terms of the decay rate from $p$ to the interval $[k,k+\delta k]$, which we denote $\mathcal{W}_{p\to k}\delta k$, $J(x,k,t)$ and $\Gamma(k)$ are given by

[TABLE]

Formally, Eq. (16) is a linear integro-differential equation. We decompose $f(x,k,t)$ into “low-energy” (below $k_{-}(x,t)$) and “high-energy” (between $k_{-}(x,t)$ and $k_{+}(x,t)$) pieces:

[TABLE]

If we focus on the region in $k$ between $k_{-}(x,t)$ and $k_{+}(x,t)$, there is no source for particles: $J(x,k,t)=0$ for $k_{-}(x,t)<k<k_{+}(x,t)$. In the accepted approximation $f_{\rm high}(x,k,t)$ satisfies the equation

[TABLE]

The initial conditions are defined by the free evolution within the time frame $t_{S}\ll t\ll\frac{1}{\Gamma(k^{m})}$. The corresponding solution of Eq. (19) is

[TABLE]

This is simply the result of Eq. (7) augmented with the finite lifetime of fermions above the Fermi surface. Immediately below this region, $k\lesssim k_{-}(x,t)$, the only contribution to the source term $J(x,k,t)$ in Eq. (16) comes from $f_{\rm high}(x,k,t)$. It is therefore appropriate that Eq. (17) may be approximated by

[TABLE]

We will comment upon the consistency of this approximation at the end of this section. This approach means that $J(x,k,t)$ and $\Gamma(k)$ are independent of $f(x,k,t)$, and the solution of Eq. (16) for $k^{0}<k<k_{-}(x,t)$ is easily verified as

[TABLE]

This corresponds to integrating over all contributions from modes which are sourced by the term $J$, and also allows for decay.

Concretely, we now consider a Hamiltonian of the form

[TABLE]

where $\rho^{R/L}_{q}=\sum_{k}{\psi^{R/L}_{k-q}}^{\dagger}{\psi^{R/L}_{k}}$ is the Fourier component of the density operator for right- and left-movers respectively. It is therefore sensible to consider $\mathcal{W}_{p\to k}$ (and accordingly $\Gamma(k)$) as that given by the single-particle decay rate of Ref. PKKG, . In terms of the decay rate $\Gamma(k^{m})$ at $k^{m}$ we can express

[TABLE]

with corrections suppressed by factors of $(p-k^{0})/k^{0}\ll 1$. Here $c=168$ is a normalisation constant such that $\int_{k^{0}}^{k^{m}}\mathrm{d}k\mathcal{W}_{k^{m}\to k}=\Gamma(k^{m})$. In evaluating Eq. (21), we assume that $t\gg t_{S}\sqrt{\lambda}$, and may replace the exponential and the rate $\mathcal{W}_{p\to k}$ by their averages in the interval:

[TABLE]

where $\bar{k}=\frac{1}{2}\left(k_{+}(x,t)+k_{-}(x,t)\right)$. and vanishes otherwise. To ease notation, we set $m=1$ for the remainder of the paper.

Having dispensed with ripples, simple geometric considerations dictate that the integrand in Eq. (22) is only nonzero for

[TABLE]

which determines an inequality for $t^{\prime}$ in Eq. (22)

[TABLE]

and so we can rewrite it as

[TABLE]

Hereinafter, $x$ should be understood as a function of $\lambda$ and $t$. Using these new variables, we may use the explicit expression for $\mathcal{W}$ given by Eq. (24) and discard subleading corrections in $\lambda$ to find

[TABLE]

To determine the behaviour of this integral, it is crucial to know the behaviour of the exponential inside the integrand. To make this clearer, we introduce the dimensionless variable

[TABLE]

In terms of this, we may write

[TABLE]

It will be helpful to also introduce a dimensionless variable $\gamma$ which interpolates in the $k$-direction between the overhanging tip ($\gamma=0$) and the background Fermi sea ($\gamma=1$), defined by

[TABLE]

Rewriting in terms of the variables $\lambda$, $\gamma$, $t$, the leading behaviour of Eq. (28) may be approximated by

[TABLE]

The contribution to the density coming from the “high energy” region i.e. $k_{-}(x,t)<k<k_{+}(x,t)$ behaves as

[TABLE]

The contribution to the density coming from the “low energy” region is given by integrating Eq. (33) over $k^{0}<k<k_{-}(x,t)$, which translates to

[TABLE]

Using Eq. (33) and noticing that the dominant contribution to $\delta\rho_{\rm low}$ comes from $\gamma\sim\mathcal{O}(1)$, Eq. (35) may be evaluated to leading order in $\lambda$, see Appendix B, yielding

[TABLE]

Here $\delta\rho_{\rm high}$ is given by Eq. (34) and we have used the function

[TABLE]

which is plotted in Fig. 2. We note that the form of $F(z)$ appears to be largely insensitive to the particular form of $\mathcal{W}$. This correction is maximal at $\lambda t\Gamma(k^{m})\sim 1$.

We have argued that the ballistic result of Eq. (34) is modified with a contribution from “lower energies”, giving

[TABLE]

which remains a monotonic function of $\lambda$, as shown in Fig. 3.

In dimensionful variables the shock wave preserves its form away from the front over the time-independent scale $x_{+}(t)-x\sim\frac{w}{\Gamma(k^{m})t_{S}}$, while the entire shock wave structure (see Fig. 4) expands linearly in time as $x_{+}(t)-x_{-}(t)\sim w\frac{t}{t_{S}}$.

With our expression for $f_{\rm low}$ in Eq. (33), it is possible to comment on the consistency of the assumption of the region $[k_{-}(x,t),k_{+}(x,t)]$ being the dominant contribution to $J(x,k,t)$, as asserted in Eq. (21). By evaluating the contribution to $J$ for momenta below $k_{-}(x,t)$, as outlined in Appendix C, one finds that the relative correction is not asymptotically small for all times, with behaviour similar to that of Eq. (36). However, the correction is in fact numerically small, on the order of $2\%$ at the largest, and so the approximation of Eq. (21) appears to be consistent.

IV Conclusions

For free fermions, introducing a density disturbance of width $w$ containing $\Delta N\gg 1$ particles leads to the formation of a shock structure on a timescale which may be estimated from classical mechanics as $t_{S}\sim\frac{w^{2}}{\Delta N}$. Semi-classical corrections in the form of quantum ripples exist in a region of the front of the shock for all times. As quantified in Section II, the size of this region grows linearly with time, but with a parametrically small coefficient. In dimensionful variables, the quantum ripples occur on the scale $x_{+}(t)-x\sim\frac{wt}{t_{S}\left(\Delta N\right)^{2/3}}$. The fraction of the shock for which ripples are significant is given by $\lambda_{\rm cr}\sim\left(\Delta N\right)^{-2/3}$, and so most of the shock is well-described classically.

Utilising this picture and disregarding ripples, one may now consider turning on interactions. Generic interactions lead to the decay of excitations high above the Fermi sea at the expense of creating a large number of low-energy excitations. These high-energy excitations have characteristic decay rate $\Gamma(k^{m})$. To establish continuity between the free and interacting pictures, we make the restriction that $t_{S}\Gamma(k^{m})\ll 1$: the profile of the shock is manifestly unchanged for short times, as interactions have barely “turned on”. Decay processes are only pertinent at $t\Gamma(k^{m})\gtrsim 1$, well after the shock has been established.

Our analysis of the linearised Boltzmann equation of Eq. (16) yields the behaviour of the fermionic density at times $t\gg t_{S}$. For any time $t\gtrsim\Gamma(k^{m})^{-1}$ the deviations from a simple exponentially decaying ballistic theory are given by the function $F(\lambda t\Gamma(k^{m}))$ in Eq. (38). This term is significant only for $\lambda\sim 1/\left[t\Gamma(k^{m})\right]$. In dimensionful variables, the kinetic corrections to the density profile become important only at distances $x_{+}(t)-x\gtrsim w/\left[t_{S}\Gamma(k^{m})\right]$ away from the tip of the shock wave, illustrated in Fig. (4). We find that the corrections to the shape are small at $x_{+}(t)-x\ll w/\Gamma(k^{m})t_{S}$, albeit the density of fermions residing in the shock wave is suppressed by the factor $e^{-t\Gamma(k^{m})}$.

As the region of quantum ripples grows linearly with time, it eventually overtakes the scale $w/\Gamma(k^{m})t_{S}$. It is meaningless to make claims on the modified shape of the classical profile once quantum ripples have encroached upon this region. The length-scale for ripples and the length-scale for kinematic corrections are comparable for $t\Gamma(k^{m})\sim\left(\Delta N\right)^{2/3}$. This time-scale therefore determines when Eq. (38) is no longer legitimate, as quantum corrections are important. Until this time, however, the main result of Eq. (38) holds: the shape of the shock is modified and subject to exponential decay of the magnitude with rate $\Gamma(k^{m})$. We note that by the time $t\sim(\Delta N)^{2/3}/\Gamma(k^{m})$ at which the “quantum” distortion meets with the “kinetic” one, destroying the shape as described by Eq. (38), its amplitude becomes exponentially small in the large parameter $(\Delta N)^{2/3}$ of the semi-classical theory.

A priori it is not clear how interactions should modify a shock wave. We have provided a picture motivated by a Boltzmann equation, where we explicitly determine the shape of the propagating shock at times well after the formation of the shock, including generic interactions. Although $\Gamma(k^{m})$ is given by a perturbative evaluation of Fermi’s Golden Rule for generic density-density interactions between spinless fermionsPKKG , the strong momentum-dependence of $\mathcal{W}$ in Eq. (24) is a consequence of the limited phase space for scattering, which remains true even in the case of spinful fermions. We conjecture that our observations should be quite generic, as the form of $F(z)$ in Eq. (37) is not sensitive to the precise details of the rate $\mathcal{W}$. The main result of Eq. (36) exhibits the corrections to the naïve picture of a kinematic shock with exponential decay.

It remains of interest to investigate if there is a direct connection with the non-linear Luttinger liquidSIG picture, in order to investigate this question in a more general, non-perturbative contextidrisovschmidt .

V Acknowledgements

We acknowledge support from the Yale Postdoctoral Prize Fellowship (TV) and NSF DMR Grant No. 1603243 (LG).

Appendix A Evaluating the bulk profile

The approach we take to understand the time-evolved behaviour of Eq. (10) for times $t\gg t_{S}$ is to consider a quadratic profile

[TABLE]

When focussing on the front of the shock only the curvature at the maximum of the initial density perturbation should be important. The time evolution of the Wigner function is given by $x\to x-\frac{\hbar kt}{m}$, and at fixed $x,t$ the roots of $k_{F}\left(x-\frac{kt}{m}\right)-k$ determine the upper and lower boundaries of support of the overhanging section of classical Wigner function. The case of Eq. (39) causes these to be the roots of a quadratic polynomial

[TABLE]

where

[TABLE]

Here we have assumed that $\lambda\gg\frac{1}{\tau}$. To evaluate the density we must plug the time-evolved Eq. (10) into Eq. (4). We wish to evaluate the excess density, defined as

[TABLE]

For the quadratic profile this is equivalent to restricting the $k$-integration of Eq. (4) to the region for which $|x-\frac{kt}{m}|\leq w$. In the new variables of Eq. (41) this corresponds to

[TABLE]

This window is sufficiently large to capture smearing of the Wigner function in the $k$-direction even for non-quadratic profiles. Integrating $\kappa$ between the bounds of Eq. (43) yields

[TABLE]

Rewriting the quadratic polynomial explicitly as in Eq. (40), and introducing $y^{\prime}\equiv y\left(\frac{k^{m}-k^{0}}{4w^{2}}\right)^{1/3}$ and $\kappa^{\prime}\equiv\kappa-\lambda$, gives

[TABLE]

where we have identified $\Delta N=2w(k^{m}-k^{0})$. We may now examine the $\tau\to\infty$ behaviour of Eq. (45). Up to a prefactor, this is in fact a function of one variable: $\lambda/\lambda_{\rm cr}$, where

[TABLE]

This may be seen by defining $p\equiv\kappa^{\prime}\tau(\Delta N)^{1/3}$, using the expressions for $\kappa_{\pm}$ from Eq. (41), and extending the range of the $p$-integration to infinity to find

[TABLE]

Taking spatial derivatives of Eq. (47), and making the substitutions (assuming $\lambda>0$)

[TABLE]

gives the form

[TABLE]

which is valid at any $\lambda/\lambda_{\rm cr}$. Using the definitions of $\lambda$, $\lambda_{\rm cr}$, $\tau$, $\ell(t)$ (of Eq. (12)), and $\Delta N$, this is equivalent to Eq. (11). The leading asymptote for $\lambda/\lambda_{\rm cr}\gg 1$ is

[TABLE]

Appendix B Determining corrections in “low energy” region

We begin from the expression for $f_{\rm low}$ of Eq. (33):

[TABLE]

The contribution to the density, given by integrating over $\gamma$ as in Eq. (35), relative to $\delta\rho_{\rm high}$ is given by

[TABLE]

We will show that this integral is dominated by $\gamma\sim\mathcal{O}(1)$. In this case we may approximate Eq. (52) by setting $\lambda=0$ in the lower integration limit and in $\left[\gamma+\lambda\right]^{11/2}$, and dropping the $(1+\tau)^{13/2}$ denominator, which will be $1+\mathcal{O}(\lambda)$. This then yields the simpler expression

[TABLE]

Eq. (53) depends only on the parameter $\lambda t\Gamma(k^{m})$, and may be simply rewritten as

[TABLE]

We will now justify this procedure. First, we split the $\gamma$ integral into two regions, focussing first on the region $\gamma<\lambda$:

[TABLE]

By changing integration variables it is clear that to leading order in $\lambda$ this is given by

[TABLE]

where $G(k)\sim k$ for $k\ll 1$, and is bounded by a constant for all $k>0$. Turning now to the region $\gamma>\lambda$:

[TABLE]

It is helpful to introduce the change of variables $z\equiv\tau t\Gamma(k^{m})\gamma$, and explicitly include $c$ from Eq. (24), such that

[TABLE]

The leading behaviour in $\lambda$ is given by setting $\lambda=0$ in both the fractional power term and the lower limit of the $\gamma$ integral, as well as discarding the denominator in the $\tau$ integral, giving

[TABLE]

with $F(z)$ given by Eq. (54). The scale for $I_{1}$ to be comparable to $I_{2}$ is $\Gamma t\lambda\sim\lambda^{-10}$. For decay to be significant we require $\Gamma(k^{m})t\gtrsim 1$, and as $\lambda\ll 1$ we are justified in considering only the contribution from $I_{2}$. Putting all this together allows us to write the leading contribution as

[TABLE]

with $\delta\rho_{\rm high}[\lambda,t]$ given by Eq. (34). We note that the form of $\mathcal{W}$ determines the particular numerical coefficients and powers appearing in the above expressions, but the form of the integral of $F$ is insensitive to this. Indeed, the factors of $1/(1+\tau)$ we neglect in the integrand of Eq. (57) come from the $(k-p)$ dependence of $\mathcal{W}$ in Eq. (24). This dependence arises from a combination of matrix elements and density-of-states of low-energy excitations. The density-of-states is small even in the presence of spin, and so it is possible that the form of $F$ survives even in the case of weakly-interacting spin-$\frac{1}{2}$ fermions.

Appendix C Consistency of approximation

We wish to understand if the assumption of Eq. (21) is consistent. In order for this to be the case the contribution to the source term from the “low-energy” region should be small compared to that of the “high-energy” region. This entails examining

[TABLE]

We can rewrite this by using the expression for $\mathcal{W}$ from Eq. (24), the dimensionless variable $\gamma$ introduced in Eq. (32), and the same approximation for the denominator as in Eq. (25) to give

[TABLE]

By applying the same approximation technique as in Appendix B, one may obtain that the leading (in $\lambda$) relative correction to the source term is given by

[TABLE]

where $F(z)$ is the same as in Eq. (54). We observe that although the corrections are $\mathcal{O}(1)$ for $\lambda t\Gamma(k^{m})\sim 1$, they are nonetheless numerically small, with the largest corrections being below $2\%$.