Nonequilibrium dynamics of noninteracting fermions in a trap

TL;DR

This paper analyzes the real-time dynamics of large numbers of noninteracting fermions in one-dimensional traps, providing exact solutions and universal scaling behaviors for their phase space distributions, correlations, and relaxation to equilibrium.

Contribution

It offers exact analytical results for the Wigner function and correlation kernels in various trapping potentials, extending equilibrium scaling laws to non-equilibrium dynamics.

Findings

Exact solutions for Wigner function, density, and correlations in harmonic and inverse square traps.

Universal edge scaling functions (Airy and Bessel types) persist during dynamics.

Fermi surf dynamics can be described by equilibrium scaling functions with time-dependent modifications.

Abstract

We consider the real time dynamics of $N$ noninteracting fermions in $d=1$. They evolve in a trapping potential $V(x)$, starting from the equilibrium state in a potential $V_0(x)$. We study the time evolution of the Wigner function $W(x,p,t)$ in the phase space $(x,p)$, and the associated kernel which encodes all correlation functions. At $t=0$ the Wigner function for large $N$ is uniform in phase space inside the Fermi volume, and vanishes at the Fermi surf over a scale $e_N$ being described by a universal scaling function related to the Airy function. We obtain exact solutions for the Wigner function, the density, and the correlations in the case of harmonic and inverse square potentials, for several $V_0(x)$. In the large $N$ limit, near the edges where the density vanishes, we obtain limiting kernels (of the Airy or Bessel types) that retain the form found in equilibrium, up to a time dependent rescaling. For non-harmonic traps the evolution of the Fermi volume is more complex. Nevertheless we show that, for intermediate times, the Fermi surf is still described by the same equilibrium scaling function, with a non-trivial time and space dependent width which we compute analytically. We discuss the multi-time correlations and obtain their explicit scaling forms valid near the edge for the harmonic oscillator. Finally, we address the large time limit where relaxation to the Generalized Gibbs Ensemble (GGE) was found to occur in the "classical" regime $\hbar \sim 1/N$. Using the diagonal ensemble we compute the Wigner function in the quantum case (large $N$, fixed $\hbar$) and show that it agrees with the GGE. We also obtain the higher order (non-local) correlations in the diagonal ensemble.

Figures

Equations

$${\rho }_{sc}(x)=\frac{2}{\pi x_e^2}\sqrt{x_e^2-x^2}$$

$$H_0(x,p)=\frac{p^2}{2}+V_0(x),H(x,p)=\frac{p^2}{2}+V(x)$$

$$W(x,p,t=0)\simeq \frac{1}{2\pi \hslash }\theta (\mu -H_0(x,p))$$

$$W(x,p,t=0)\simeq \frac{1}{2\pi \hslash }\mathcal{W}(a),\mathcal{W}(a)={\int }_{2^{2/3}a}^{+\infty }\mathrm{Ai}(u)du$$

$$a=\frac{H_0(x,p)-\mu }{e(x,p)},e(x,p)=\frac{{\hslash }^{\frac{2}{3}}}{2^{\frac{1}{3}}}[p^2{V}_0^{\prime \prime }(x)+V_0^{\prime }(x{)}^2{]}^{\frac{1}{3}}$$

$$\Psi (x_1,\cdots x_N;t)=\frac{1}{\sqrt{N!}}\underset{1\le k,j\le N}{\det }{\psi }_k(x_j,t)$$

$$K_{\mu }(x,x^{\prime },t)=\sum _{k=1}^N{\psi }_k^{\ast }(x,t){\psi }_k(x^{\prime },t)$$

$$W(x,p,t)=\frac{N}{2\pi \hslash }{\int }_{-\infty }^{+\infty }d{x}_2..d{x}_{N}dy{e}^{\frac{ipy}{\hslash }}$$

$$\times {\Psi }^{\ast }(x+\frac{y}{2},x_2,\cdots x_N;t)\Psi (x-\frac{y}{2},x_2,\cdots x_N;t)$$

$$W(x,p,t)=\frac{1}{2\pi \hslash }{\int }_{-\infty }^{+\infty }dy{e}^{\frac{ipy}{\hslash }}{K}_{\mu }(x+\frac{y}{2},x-\frac{y}{2},t)$$

$$\hslash {\partial }_t{K}_{\mu }=-\frac{i{\hslash }^2}{2}({\partial }_x^2-{\partial }_{x^{\prime }}^2)K_{\mu }+i(V(x)-V(x^{\prime }))K_{\mu }$$

$${\partial }_{t}W=\left(-p{\partial }_x+\frac{i}{\hslash }V(x-\frac{i\hslash }{2}{\partial }_p)-\frac{i}{\hslash }V(x+\frac{i\hslash }{2}{\partial }_p)\right)W$$

$${\partial }_{t}W=\left(-p{\partial }_x+V^{\prime }(x){\partial }_p\right)W$$

$$W(x,p,t)$$

$$=W_0(x\cos (\omega t)-\frac{p}{\omega }\sin (\omega t),p\cos (\omega t)+\omega x\sin (\omega t))$$

$$x_e(t)=\frac{\sqrt{\mu }}{{\omega }_0\omega }\sqrt{{\omega }^2+{\omega }_0^2+({\omega }^2-{\omega }_0^2)\cos (2\omega t)}$$

$${\psi }_k(x,t)=\frac{e^{i\frac{L^{\prime }(t)}{2L(t)\hslash }{x}^2-i\frac{{ϵ}_k^0\tau (t)}{\hslash }}}{\sqrt{L(t)}}{\varphi }_k^0(\frac{x}{L(t)}),$$

$${\partial }_t^{2}L(t)+\omega (t{)}^{2}L(t)=\frac{{\omega }_0^2}{L(t{)}^3}$$

$$K_{\mu }(x,x^{\prime },t)=\frac{e^{i\frac{L^{\prime }(t)}{2L(t)\hslash }({x^{\prime }}^2-x^2)}}{L(t)}{K}_{\mu }(\frac{x}{L(t)},\frac{x^{\prime }}{L(t)},0)$$

$$W(x,p,t)=W_0(\frac{x}{L(t)},pL(t)-x{L}^{\prime }(t))$$

$$\rho (x,t)=\frac{1}{w_N(t)}{F}_1(\frac{x-x_e(t)}{w_N(t)})$$

$${\rho }_N(x,t)=\frac{2}{\pi x_e}\{\begin{array}{ll}& 0,\sin \theta <-\sin \omega t\\ & \frac{\sin (\omega t+\theta )}{\sin \omega t},|\sin \theta |<\sin \omega t\\ & 2\cos \theta ,\sin \theta >\sin \omega t\end{array}$$

$$H(x,p,t)=\frac{p^2}{2}+\frac{1}{2}\omega (t{)}^2{x}^2+\frac{\alpha (\alpha -1){\hslash }^2}{2x^2},x>0$$

$${\rho }_N(x,t)\simeq \frac{1}{L(t)}\tilde{\rho }(\frac{x}{L(t)},0),\tilde{\rho }(x)=\frac{4\theta (x)}{\pi x_e^2}\sqrt{x_e^2-x^2}$$

$$K_{\mu }(x,y,t)\simeq \frac{2k_F^2\sqrt{x{x}^{\prime }}}{L(t{)}^2}{K}_{\nu }^{\mathrm{Be}}(\frac{k_F^2{x}^2}{L(t{)}^2},\frac{k_F^2{y}^2}{L(t{)}^2})$$

$$K_{\nu }^{\mathrm{Be}}(u,v)=\frac{\sqrt{v}{J}_{\nu }^{\prime }(\sqrt{v})J_{\nu }(\sqrt{u})-\sqrt{u}{J}_{\nu }^{\prime }(\sqrt{u})J_{\nu }(\sqrt{v})}{2(u-v)}$$

$$\big{(}\partial_{t}+p~{}\partial_{x}-V^{\prime}(x)\partial_{p}+\frac{\hbar^{2}}{24}V^{\prime\prime\prime}(x)\partial_{p}^{3}\big{)}W=0$$

$$W(x,p,t)=\frac{1}{2\pi\hbar}\theta\big{(}\mu-\frac{p_{0}(x,p,t)^{2}}{2}-V_{0}(x_{0}(x,p,t))\big{)}$$

$$W(x,p,t)=\frac{1}{2\pi \hslash }\theta (p-p_{-}(x,t))\theta (p_{+}(x,t)-p)$$

$${\partial }_{t}p(x,t)+p(x,t){\partial }_{x}p(x,t)+V^{\prime }(x)=0$$

$$p_{\pm }(x,t)=\pm p_{\mu }^0(x),p_{\mu }^0(x)=\sqrt{2(\mu -V_0(x){)}_{+}}$$

Full text

There was recent progress in the theoretical study of noninteracting fermions in a confining trap at equilibrium [1, 2, 3, 4, 5, 6, 7, 8, 9, 10], motivated in part by progress in manipulating and imaging cold atoms [11, 12, 13, 14, 15]. Although they are noninteracting, because of Pauli exclusion they exhibit non trivial quantum correlations. In particular at zero temperature and in one dimension $d=1$, it was shown that the correlations are universal and related to random matrix theory (RMT) [16, 17]. For instance, for $N$ fermions in a harmonic trap $V(x)=\frac{1}{2}m\omega^{2}x^{2}$, the positions $x_{i}$ are in one-to-one correspondence with the eigenvalues of an $N\times N$ random matrix drawn from the Gaussian unitary ensemble. This implies that the average density (with unit normalization) is given for large $N$ by the Wigner semi-circle $\rho_{N}\simeq\rho_{sc}$,

[TABLE]

which has edges at $x=\pm x_{e}$, $x_{e}=\sqrt{\frac{\hbar}{m\omega}}\sqrt{2N}$. In the bulk, i.e. away from the edges, the density correlations are described by the sine-kernel which can be obtained through the local density approximation (LDA) [2] or in a more controlled way using the connection to RMT [8]. Near the edge, within a scale $w_{N}\sim N^{-1/6}$, the fluctuations are enhanced and are described by the so-called Airy kernel [4, 5, 6]. This result extends to a large class of smooth potentials [8]. In presence of hard edges, similar results were shown in terms of the Bessel kernel [18, 19, 20]. Both kernels are well known to describe soft and hard edges in RMT. Extensions to finite temperature $T>0$, $d>1$ and interactions were studied in [6, 21, 7, 22, 23, 24].

It is natural to ask how these properties behave under real time quantum dynamics. For equilibrium dynamics this was studied in [25], and related to the extended Sine and Airy processes in imaginary time. Here we consider non-equilibrium quantum quenches, first from zero temperature $T=0$, and later from finite temperature) and study the time evolution of the density and correlations. In the $T=0$ quench, the many body system is prepared in the ground state of a given Hamiltonian ${\cal H}_{0}$ and evolves under the Hamiltonian ${\cal H}\neq{\cal H}_{0}$. Such quantum quenches in traps have been studied in a number of works, however most of them focus on interacting bosons, in particular in $d=1$ in the Tonks-Girardeau limit of infinite repulsion [27, 26, 29, 28] and were studied for finite $N$. Although in this limit there is a formal relation to noninteracting fermions [18, 30], the correlations, beyond the density, are different. For instance, the possible edge universality of the time dependent fermion problem at large $N$ has not been addressed. Very recently the evolution in the bulk was investigated in the limit $\hbar\sim 1/N$, $N\to\infty$ (which we call here ”classical”), and at large time thermalization was found to occur for a one-particle observable [31].

We consider noninteracting fermions in the quantum regime for fixed $\hbar$. The initial and evolution many body Hamiltonians thus take the forms ${\cal H}_{0}=\sum_{i=1}^{N}H_{0}(x_{i},p_{i})$ and ${\cal H}=\sum_{i=1}^{N}H(x_{i},p_{i})$, in terms of the single particle Hamiltonians (below we work in units where $m=1$)

[TABLE]

A convenient way to study noninteracting fermions [32, 33] is via the Wigner function $W(x,p,t)$ (see definition below) which plays the role of a (quasi-)probability density in the phase space $(x,p)$ [34]. In particular $\int dp~{}W(x,p,t)=\tilde{\rho}(x,t)=N\rho_{N}(x,t)$, the time dependent number density. We choose the initial condition $W(x,p,t=0)$ as the Wigner function associated to the ground state of ${\cal H}_{0}$, studied in [35, 36, 37] for various potentials $V_{0}$. It was shown that in the large $N$ limit it takes the form in the bulk

[TABLE]

with $\theta(x)$ the Heaviside step function. It depends on $N$ only through $\mu$, the Fermi energy, i.e. the energy of the highest occupied single particle state, which is an increasing function of $N$. The curve $(x_{e},p_{e})$ where $W$ in (3) vanishes, i.e. $\mu=H_{0}(x_{e},p_{e})$, defines an edge in phase space (the Fermi surf $S_{0}$, enclosing the Fermi volume $\Omega_{0}$). At large $N$, i.e. for large $\mu$ the step has a finite width $e_{N}=e(x_{e},p_{e})\ll\mu$, i.e. near the edge and for smooth potentials $V_{0}(x)$, the Wigner function takes the universal scaling form

[TABLE]

in terms of the scaled distance $a$ to the Fermi surf

[TABLE]

Our goal in this paper is to study the time dependent Wigner function starting from this initial condition, e.g. (3), (4), and to understand how its bulk and edge properties evolve with time, in particular whether their universality is preserved under nonequilibrium dynamics. We obtain exact solutions for time dependent harmonic and inverse square potentials, and calculate the Wigner function and the density for various $V_{0}(x)$. For large $N$, near the edges where the density vanishes, we obtain limiting kernels (of the Airy or Bessel types) that retain the form found in the statics, up to a time dependent rescaling. It implies the occurrence of the Tracy-Widom distribution [38] for the position of the rightmost fermion [6]. For generic non-harmonic traps, at large $N$ and intermediate times, the Wigner function is still uniform on the Fermi volume, whose evolution can be more complex. We show that at generic points the Fermi surf is still described by the same universal function as in equilibrium (4), (33), with a space and time dependent width that we calculate. We express multi-point and multi-time correlations in terms of the Wigner function and its associated kernel. In the infinite time limit, we obtain exact formula for the Wigner function under a thermalization hypothesis, which we show coincide with the prediction of the Generalized Gibbs Ensemble (GGE) [39, 40].

The time evolved $N$ body wave function keeps the form of a Slater determinant

[TABLE]

where the $\psi_{k}(x,t)$ are the solutions of $i\hbar\partial_{t}\psi=H\psi$ with initial condition $\psi_{k}(x,t=0)=\phi^{0}_{k}(x)$, where $\phi^{0}_{k}(x)$, $k\geq 1$, are the single particles eigenstates of $H_{0}$ in increasing order of energies $\epsilon^{0}_{k}$. We are interested in the quantum probability $|\Psi(x_{1},\cdots x_{N};t)|^{2}=\frac{1}{N!}\det_{1\leq i,j\leq N}K_{\mu}(x_{i},x_{j},t)$, described by the time dependent kernel

[TABLE]

where the $\psi_{k}(x,t)$’s form an orthonormal basis for all $t$. All $m$-point correlations $R_{m}(x_{1},\dots x_{m};t)\mathrel{\mathop{\mathchar 58\relax}}=\frac{N!}{(N-m)!}\prod_{j=m+1}^{N}\int dx_{j}|\Psi(x_{1},\cdots x_{N};t)|^{2}$ are determinants, $R_{m}(x_{1},\dots x_{m};t)=\det_{1\leq i,j\leq m}K_{\mu}(x_{i},x_{j},t)$, i.e. the $x_{i}$’s form a determinantal point process [41, 42]. One defines the $N$-body Wigner function

[TABLE]

which is related to the kernel via [37]

[TABLE]

Using the Schrödinger equation for the $\psi_{k}(x^{\prime},t)$ in (7) we see that the kernel $K_{\mu}(x,x^{\prime},t)$ satisfies

[TABLE]

which, via Fourier transformation, leads to the evolution equation for $W(x,p,t)$, i.e. the Wigner equation (WE)

[TABLE]

Although derived here for any $N$, this equation does not explicitly depend on $N$, hence it is identical to the equation for $N=1$ [34]. It was also obtained in second quantized form in [43, 31]. The crucial difference thus lies in the initial condition, discussed above in the large $N$ limit, see (3)-(4)-(5) at $T=0$. In the classical limit, i.e. setting $\hbar=0$ in the WE, one obtains the Liouville equation (LE)

[TABLE]

It verifies $dW(x(t),p(t),t)/dt=0$ along the classical trajectories $\dot{x}=p$ and $\dot{p}=-V^{\prime}(x)$. Denoting $x_{0}(x,p,t)$, $p_{0}(x,p,t)$ the initial conditions (at $t=0$) as functions of the final ones (at $t$), the general solution of the LE can thus be written as $W(x,p,t)=W_{0}(x_{0}(x,p,t),p_{0}(x,p,t))$ where $W_{0}(x,p)=W(x,p,0)$ is the initial Wigner function.

Let us start with the simplest solvable case of the harmonic oscillator, $V(x)=\frac{1}{2}\omega^{2}x^{2}$. Since $V^{\prime\prime\prime}(x)=0$, the WE (11), reduces to the LE (12). Hence the solution is

[TABLE]

which is a simple rotation in phase space of the initial Wigner function and is valid for any $N$. Let us specify the initial Hamiltonian $H_{0}$ to be a harmonic oscillator, i.e. $V_{0}(x)=\frac{1}{2}\omega_{0}^{2}x^{2}$ in (2). Hence the initial density $\rho_{N}(x,t=0)$ is the semi-circle (1) for large $N$ with $x_{e}=\sqrt{2\mu}/\omega_{0}=\sqrt{2N\hbar/\omega_{0}}$. Substituting the initial condition (3) into (13) one immediately finds that for large $N$ the density remains a semicircle with $x_{e}\to x_{e}(t)$

[TABLE]

and $\mu=N\hbar\omega_{0}$. In fact a solution of the general time dependent harmonic oscillator, $H(t)=\frac{p^{2}}{2}+\frac{1}{2}\omega(t)^{2}x^{2}$, can be obtained for the same initial condition using the ‘rescaling’ method [44, 45]. The wave functions are given by

[TABLE]

where $\tau(t)=\int_{0}^{t}dt^{\prime}/L^{2}(t^{\prime})$ and $L(t)$ satisfies Ermakov’s equation

[TABLE]

with $L(0)=1$, $L^{\prime}(0)=0$. The kernel and the density are thus given in terms of their initial values, for any $N$ by

[TABLE]

and $\rho_{N}(x,t)=K_{\mu}(x,x,t)/N=\frac{1}{L(t)}\rho_{N}(\frac{x}{L(t)},0)$. In addition we obtain [46] the exact formula for

[TABLE]

which reduces to (13) for $\omega(t)=\omega$. Hence we see that for large $N$ the density remains a semi-circle with an edge at $x_{e}(t)=L(t)x_{e}$, formula (14) being recovered in the case $\omega(t)=\omega$. For the kernel, the scaling forms obtained at equilibrium, namely the sine-kernel (in the bulk) and the Airy kernel (at the edge) [8], are preserved, up to a change in the scale. In particular the density at the edge reads

[TABLE]

with $F_{1}(y)={\rm Ai}^{\prime}(y)^{2}-y{\rm Ai}(y)^{2}$ and a width $w_{N}(t)=w_{N}L(t)$ and $w_{N}=\sqrt{\frac{\hbar}{2\omega_{0}}}N^{-1/6}\sim\mu^{-1/6}$.

The Wigner function method allows to study a variety of initial conditions. One example is $H_{0}$ given by the $\omega_{0}$ harmonic oscillator restricted to $x>0$ with an impenetrable wall at $x=0$. The initial kernel, $K^{+}_{\mu}$, is obtained via the image method $K^{+}_{\mu}(x,x^{\prime})=(K_{\mu}(x,x^{\prime})-K_{\mu}(x,-x^{\prime}))\theta(x)\theta(x^{\prime})$. The associated Wigner function in the bulk is unity on the half-ellipse $H_{0}(x,p)<\mu$ with $x>0$. Under evolution with $\omega(t)=\omega$, this half ellipse rotates clockwise (see Fig. 1). Consider for simplicity $\omega=\omega_{0}$. Defining $x=x_{e}\sin\theta$, $\theta\in[-\frac{\pi}{2},\frac{\pi}{2}]$, where $x_{e}=\sqrt{2\mu}/\omega_{0}$, the density reads for $0<\omega t<\pi/2$

[TABLE]

and in addition satisfies $\rho_{N}(x,t)=\rho_{N}(-x,\frac{\pi}{\omega}-t)$ and is time periodic of period $2\pi/\omega$. Hence it is not a semi-circle, and there is now a moving front (e.g. at $x=-x_{e}\sin\omega t$ for $\omega t<\frac{\pi}{2}$) where the density vanishes linearly (see [46] for details). A formula for $W$ can also be obtained for any $\omega(t)$ and arbitrary initial condition (see [46]).

It is possible to extend the exact solution to the case of an additional $1/x^{2}$ wall at $x=0$, following [47]. Consider

[TABLE]

and denote $\omega_{0}=\omega(0)$. The initial kernel has an explicit expression [46] in terms of Laguerre polynomials for any $N$. At large $N$, the density in the bulk is a time-dependent half semi-circle

[TABLE]

where $L(t)$ is the same solution of (16), e.g. $L(t)=x_{e}(t)/x_{e}$ for $\omega(t>0)=\omega$. Near the wall the kernel takes the form

[TABLE]

up to a phase factor, where $\mu=\frac{\hbar^{2}}{2}k_{F}^{2}$, $\nu=\alpha+\frac{1}{2}$, in terms of the Bessel kernel

[TABLE]

characteristic of the ”hard-edge” universality class of RMT [22], thus preserved by the dynamics in this case.

Consider now generic smooth confining potentials $V(x)$, $V_{0}(x)$, with $V^{\prime\prime\prime}(x)\neq 0$. Our strategy is as follows. We first argue that, for large $N$ (or equivalently large $\mu$), $W$ obeys Liouville equation (12) (and at $T=0$ is a theta function, i.e. takes value in $\{0,1\}$). We then study the solution to the LE, leading to Burgers equation and fermion hydrodynamics. Finally, we obtain the leading (quantum) corrections to that picture, of foremost importance near the Fermi surf.

Let us expand (11) to the first ”quantum” order as

[TABLE]

To check that the quantum terms are subdominant, we can evaluate each term at $t=0$ using the initial condition (4)-(5). The derivatives are non-zero only near the Fermi-surf: consider a generic point there, $(x_{e},p_{e})$, such that $p_{e}^{2}\sim V_{0}(x_{e})\sim\mu$, leading to $e_{N}\sim(\hbar\mu/x_{e})^{2/3}$. As $t$ increases, we assume that $V(x)$ and its derivatives still scale the same way as $V_{0}(x)$, i.e. $V^{(n)}(x_{e})\sim\mu/x_{e}^{n}$. Each order in $\hbar$ brings $\partial_{x}$ acting on $V$ and $\partial_{p}$ acting on $W$, hence, using the scaling form (4)-(5), a factor either $\sim\hbar/(x_{e}p_{e})$ or $\sim\hbar p_{e}/(x_{e}e_{N})$. The latter dominates, but is $\sim\hbar/(x_{e}^{1/3}\mu^{1/6})\ll 1$ at large $\mu$.

Under the Liouville equation the Fermi volume $\Omega_{t}$ and surf $S_{t}$ are simply transported by the (phase-space area preserving) classical equation of motion, i.e.

[TABLE]

There are several parametrizations to describe the transported (and possibly wildly deformed) Fermi volume $\Omega_{t}$ [46]. In the simplest case one can write

[TABLE]

for $x\in{\cal I}=[x^{-}_{e}(t),x^{+}_{e}(t)]$ the (single) support of the number density $\tilde{\rho}(x,t)=\frac{1}{2\pi\hbar}(p_{+}(x,t)-p_{-}(x,t))$. For $x\in{\cal I}$, each $p_{\pm}$ must then satisfy the Burgers equation [48]

[TABLE]

It recovers the standard free fermion hydrodynamics with a velocity field $v=\frac{1}{2}(p_{+}+p_{-})$ with $\tilde{\rho}v=\int dp~{}pW$, and the local kinetic energy density $\int dp~{}\frac{p^{2}}{2}W=\frac{1}{2}\tilde{\rho}v^{2}+\frac{\hbar^{2}}{6}\pi^{2}\tilde{\rho}^{3}$ [46]. The equilibrium is recovered for $V(x)=V_{0}(x)$ with

[TABLE]

and $W$ is time independent. The form (27) using (28) automatically satisfies Liouville equation.

We now study the behavior of the Wigner function near the transported Fermi-surf $S_{t}$. We know that at $t=0$ it is given by (4)-(5) and has a thickness $e_{N}=e(x_{e},p_{e})$. We thus look for a solution of (25) for $p\approx p_{+}(x,t)$ of the form

[TABLE]

with a time-dependent thickness $D(x,t)$, $F(+\infty)=0$ and $F(-\infty)=1$. Remarkably, plugging this form into Eq.(25) allows to determine uniquely [46]

[TABLE]

and $F^{\prime\prime\prime}(a)=4aF^{\prime}(a)$, which is obeyed by $F(a)={\cal W}(a)={\rm Ai}_{1}(2^{2/3}a)$. It also provides the correct matching with the initial condition, since $-(\frac{\hbar^{2}}{2}\partial^{2}_{x}p_{\mu}^{0}(x))^{1/3}=\frac{e(x,p)}{p}|_{p=p_{\mu}^{0}(x)}$ on the Fermi surf. The solution corresponding to $\theta(p-p_{-})$ is obtained similarly with $\tilde{a}\to-\tilde{a}$. Furthermore we show [46] that for smooth potentials, adding the missing higher order terms in $\hbar$ of the WE in Eq.(25) does not change the result. Thus derivatives up to and including $V^{\prime\prime\prime}$ determine the dynamical width, while only $V^{\prime\prime}$ and $V^{\prime}$ enter the static width $e_{N}$. Hence we have found the universal scaling form near a generic point of the transported Fermi-surf $S_{t}$. The width (31) can be calculated for any $V(x)$ solving the Burgers equation (28) with initial condition $p_{+}(x,t=0)=p_{\mu}^{0}(x)$. Note that related results were obtained [49, 50, 51, 52] for a single particle in the semi-classical limit, and for fermions at equilibrium [35, 36]. Other expressions can be obtained [46] at non-generic points, e.g. when $\partial^{2}_{x}p(x,t)=0$ [51, 53, 37]. We assumed here that the Fermi surf maintained its integrity, which in some cases requires finite time. More generally, the section $\Omega_{t}(x)$ of $\Omega_{t}$ (such that $p\in\Omega_{t}(x)$ iff $(x,p)\in\Omega_{t}$) can contain multiple intervals $\Omega_{t}(x)=\cup_{j=1}^{n(x)}[p^{j}_{-}(x,t),p^{j}_{+}(x,t)]$ which merge or disappear as $x$ varies. Locally, however, as in the statics as long as one can find a smooth local parametrization of the Fermi surf, similar scaling form as (30) should hold. This dynamical scenario remains valid over a time scale when the Fermi surf does not change drastically from the initial condition. Beyond this time scale, the dynamics becomes quite different, as discussed later. However, it is hard to estimate precisely this time scale.

The above results can be extended to a quench from an initial system prepared at finite temperature $T=1/\beta$, in the grand canonical (GC) ensemble with chemical potential $\tilde{\mu}$. The equal time correlations are again determinantal with the GC kernel (overbars denote GC averages)

[TABLE]

with $\bar{n}^{0}_{k}=1/(1+e^{\beta(\epsilon^{0}_{k}-\tilde{\mu})})$ the mean occupation number of the energy level $\epsilon^{0}_{k}$ of $H_{0}$. The GC Wigner function $W_{\tilde{\mu}}$ (defined in [37]) is again the Fourier transform of the GC kernel, as in (9). It is clear from (32) that it obeys (10) hence $W$ still the solution of (11), albeit with a different initial condition. For large $N$, in the bulk it reads

[TABLE]

Introducing the important dimensionless parameter $b=e_{N}/T$, the (same) Fermi surf is still well defined at low $T=e_{N}/b\ll\mu$ (with $\tilde{\mu}\simeq\mu$) and Eq. (4) still applies with

[TABLE]

a universal function depending only on the dimensionless parameter $b$ (with ${\cal W}_{+\infty}(a)={\cal W}(a)$ in the $T=0$ limit). Formula (13), (17), (18), readily apply with, for large $N,\tilde{\mu}$, these new initial conditions. Eqs. (20), (24), (19) admit $T>0$ extensions [46].

From the linearity of the W-equation (11) one can express its solution at $T>0$ in terms of the $T=0$ solution as

[TABLE]

since it holds at $t=0$ (equilibrium) as shown in [37]. At low $T$ near the Fermi surf, it allows to generalize Eq. (30)

[TABLE]

where the dependence in $\mu$ of $p_{+}(x,t,\mu)$ is indicated explicitly, and the parameter $b$ acquires a dependence in $x,t$, which however disappears for the harmonic oscillator for which $b(x,t)=b=\beta(\hbar\omega_{0})^{2/3}\mu^{1/3}$, as for equilibrium.

We now study the multi-time quantum correlations. As for equilibrium dynamics [25] one shows from the Eynard-Mehta theorem [42] that they are determinants based on the non-equilibrium space-time extended kernel [54]

[TABLE]

with $\theta_{+}(x)=1$ for $x>0$ and $\theta_{+}(x)=0$ for $x\leq 0$. For instance the density-density correlation for $t_{1}<t_{2}$ reads $\langle\hat{\rho}(x_{1},t_{1})\hat{\rho}(x_{2},t_{2})\rangle=\det_{1\leq i,j\leq 2}K_{\tilde{\mu}}(x_{i},t_{i};x_{j}t_{j})$. Consider now $H$ to be the $\omega(t)$-harmonic oscillator and use the rescaling method. If $H_{0}$ is the $\omega_{0}$-harmonic oscillator, one easily obtains, up to a phase factor [46]

[TABLE]

where $L(t)$ obeys (16), and $\tau(t)=\int_{0}^{t}dt^{\prime}/L(t^{\prime})^{2}$. For $\omega(t)=\omega$, then $L(t)=x_{e}(t)/x_{e}$ as given in (14) and $\tan\omega_{0}\tau(t)=\frac{\omega_{0}}{\omega}\tan\omega t$. Here $K_{\tilde{\mu}}^{\rm eq}(x,t;x^{\prime},t^{\prime})$ is the equilibrium two time kernel of the $\omega_{0}$-harmonic oscillator studied in [25], given by (36) with $\psi_{k}^{*}(x,t)\psi_{k}(x^{\prime},t^{\prime})\to{\phi^{0}_{k}}^{*}(x)\phi^{0}_{k}(x^{\prime})e^{-i\omega_{0}k(t^{\prime}-t)}$. This relation between $K_{\tilde{\mu}}$ and $K_{\tilde{\mu}}^{\rm eq}$ is special to the harmonic oscillator. $K_{\tilde{\mu}}^{\rm eq}$ is periodic in both times of period $2\pi/\omega_{0}$, and $K_{\tilde{\mu}}$ of period $2\pi/\omega$. Within a period $K_{\tilde{\mu}}^{\rm eq}$ decays very fast in $t-t^{\prime}$ with a time scale $\sim 1/N$ in the bulk and $\sim 1/N^{1/3}$ at the edge [25]. One can thus expand near $t\simeq t^{\prime}$ and one obtains the following scaling forms for $K_{\tilde{\mu}}$. In the bulk

[TABLE]

where ${\cal K}^{\rm bulk,r}$ is the real time continuation of the equilibrium extended sine kernel given at $T=0$ e.g. in Eq. (70-71) of [25] and [46]. The width $\xi_{x,t}=2\hbar L(t)^{2}/(\omega_{0}\sqrt{x_{e}(t)^{2}-x^{2}})\sim N^{-1/2}$ and the velocity $v_{xt}=xL^{\prime}(t)/L(t)\sim N^{1/2}$ are new features of the non-equilibrium dynamics ($v_{xt}=0$ at equilibrium). Near the edge it takes the scaling form

[TABLE]

where $v_{e}(t)=x_{e}(t)L^{\prime}(t)$ is the velocity of the edge and ${\cal K}_{b}^{\rm edge,r}$ is the equilibrium extended Airy kernel (continued to real time) given e.g. in Eq. (156) of [25] and [46].

Finally we discuss the large time behavior. In the ”classical” limit $\hbar\sim 1/N$, $N\to\infty$, it was shown recently [31] that for generic confining potentials, e.g. different from the HO, the Wigner function has a large time stationary limit. It was found to coincide with the prediction of the GGE. Here we investigate this question for fixed $\hbar$. Starting from the time-dependent kernel (32), and expanding on the basis of the eigenfunctions $\phi_{\ell}(x)=\langle x|\phi_{\ell}\rangle$ of $H$ (of energies $\epsilon_{\ell}$) we obtain

[TABLE]

For a confining potential the spectrum of $H$ is non-degenerate, hence the only non-oscillating terms are $\ell=\ell^{\prime}$. Keeping only these terms one obtains

[TABLE]

where the $0\leq\nu_{\ell}\leq 1$ are ”effective mean occupation numbers”, with $\sum_{\ell}\nu_{\ell}=\sum_{k}\bar{n}_{k}^{0}=\bar{N}$ (i.e. $N$ at $T=0$), as expected from particle number conservation. This is the so-called diagonal approximation (DA) [39] which, in the present case, can also be obtained as a time average, i.e. $K^{\rm di}_{\tilde{\mu}}(x,y)=\lim_{\tau\to\infty}\frac{1}{\tau}\int_{0}^{\tau}dtK_{\tilde{\mu}}(x,y,t)$. Under certain conditions, including absence of time periodicity, the large time limit exists $K^{\infty}_{\tilde{\mu}}(x,y)=\lim_{t\to\infty}K_{\tilde{\mu}}(x,y,t)$ and then $K_{\tilde{\mu}}^{\infty}=K_{\tilde{\mu}}^{\rm di}$. Let us define a fonction $\nu_{\tilde{\mu}}(\epsilon)$ such that $\nu_{\tilde{\mu}}(\epsilon_{\ell})=\nu_{\ell}$ for all $\ell$. Then one can obtain the diagonal Wigner function

[TABLE]

in terms of the Wigner function, $W^{H}_{\mu}(x,p)$, associated to the ground state of ${\cal H}$ i.e to the kernel $K^{H}_{\mu}(x,y)=\sum_{\ell}\theta(\mu-\epsilon_{\ell})\phi_{\ell}^{*}(x)\phi_{\ell}(y)$ (one first shows that $K_{\tilde{\mu}}^{\rm di}$ and $K_{\mu}^{H}$ are related via (41), which implies (41) for their Wigner functions). Note that if one defines $\hat{\nu}_{\tilde{\mu}}(\epsilon)=\sum_{\ell}\nu_{\ell}\delta(\epsilon-\epsilon_{\ell})$, then $\hat{\nu}_{\tilde{\mu}}(\epsilon)=\nu_{\tilde{\mu}}(\epsilon)\rho(\epsilon)$, where $\rho(\epsilon)=\sum_{\ell}\delta(\epsilon-\epsilon_{\ell})$ is the density of states. The above holds for arbitrary $N,\tilde{\mu}$.

For large $N$, we insert in (41) the known form of the Wigner function from (4), (5), (33), i.e. $W^{H}(x,p)\simeq\frac{1}{2\pi\hbar}{\cal W}_{b}(a)$ and obtain one of our main result

[TABLE]

where at $T=0$, ${\cal W}^{\prime}_{+\infty}(a)=-2^{1/3}{\rm Ai}(2^{2/3}a)$. If one neglects the quantum width of the edge, $\delta\mu^{\prime}\sim e_{N}\sim\hbar^{2/3}$ for fixed $\mu^{\prime}$, it is equivalent to insert instead $W^{H}(x,p)\simeq\frac{1}{2\pi\hbar}\theta(\mu-H(x,p))$ in (41). This leads to

[TABLE]

an approximation valid when $e_{N}$ is much smaller than the width $\delta\mu^{\prime}$ over which the effective occupation number $\nu_{\tilde{\mu}}(\mu^{\prime})$ drops from $1$ (for $\mu^{\prime}\ll\tilde{\mu}$) to [math] (for $\mu^{\prime}\gg\tilde{\mu}$). This is the case in the classical limit $\hbar\sim 1/N$ which we now consider [55]. We note that $\frac{1}{N}\hat{\nu}_{\tilde{\mu}}(\epsilon)$ is the probability density of the energy $H$ of a single particle in the initial state. In the classical limit it is evaluated through a phase space average over the initial Wigner function $W(x,p,0)$, which together with (43) leads to

[TABLE]

where $\rho^{sc}(\epsilon)$ is the semi-classical density of states. This is our result in the classical limit $\hbar\sim 1/N$. In the simplest case considered in [31] (i.e. $T=0$ and a single orbit at each energy $\epsilon$) one has $\rho^{sc}(\epsilon)=\frac{T(\epsilon)}{2\pi\hbar}$ where $T(\epsilon)$ is the period of the classical orbit at energy $\epsilon$. Inserting $W(x_{0},p_{0},0)=\frac{1}{2\pi\hbar}\theta(\mu-H_{0}(x_{0},p_{0}))$, integrating over $p_{0}$ in (44), and putting all together [46], one recovers the large time limit result of [31] (obtained by a quite different method) valid in the classical case. Our formula (41) is thus a good candidate for the large time limit of the Wigner function in the quantum regime.

The GGE thermalization hypothesis [39, 40, 56] states that for local observables $O$, $\lim_{t\to+\infty}\langle O(t)\rangle={\rm Tr}{\cal D}_{\rm GGE}O$, where the density matrix ${\cal D}_{\rm GGE}$ involves an extensive number of conserved quantities. For non-interacting fermions, as discussed in [39, 31], these are the occupation numbers $c^{\dagger}_{\ell}c_{\ell}$ of single particle energy levels $\epsilon_{\ell}$ of $H$. The natural candidate is thus ${\cal D}_{\rm GGE}=\frac{1}{Z_{\rm GGE}}e^{\sum_{\ell}f_{\ell}c^{\dagger}_{\ell}c_{\ell}}$ where the $f_{\ell}$ are determined from the initial condition. One has $\langle c^{\dagger}_{\ell}c_{\ell}\rangle_{t=0}=\nu_{\ell}$, where the $\nu_{\ell}$ are exactly the ones obtained in (40) [46]. Hence one finds $\langle c^{\dagger}_{\ell}c_{\ell}\rangle_{\rm GGE}=\frac{1}{1+e^{-f_{\ell}}}=\nu_{\ell}$. Thus we find that at the level of the single particle observables (the kernel and Wigner function) the result of the GGE is identical to the diagonal approximation.

What about the multi-point time dependent correlations ? Since for any initial state equal to an eigenstate of ${\cal H}_{0}$, labeled by a set ${\bf n}^{0}=\{n^{0}_{k}\}$ of occupation numbers, $n^{0}_{k}=0,1$, the time evolved state is determinantal, correlations in this state have the form

[TABLE]

Expanding on the eigenstates of $H$ now involves a sum over two sets of $m$-fermion eigenstates of ${\cal H}_{m}$. The DA restricts that sum to identical eigenstates in both sets, and one obtains [46] (after a grand canonical average)

[TABLE]

where $\nu_{\ell,\ell^{\prime}}=\langle\phi_{\ell}|(1+e^{\beta(H_{0}-\mu)})^{-1}|\phi_{\ell^{\prime}}\rangle=\langle c_{\ell}^{\dagger}c_{\ell^{\prime}}\rangle_{t=0}$. This DA for $p\geq 2$ differs from the prediction of the simplest GGE ansatz discussed above. Indeed, the latter predicts that the $m$-point correlation $R_{m}$ is equal to the $m\times m$ determinant built from the kernel $K^{\rm GGE}=K^{\rm di}_{\tilde{\mu}}$ given in (40). Agreement $R_{m}^{\rm di}=R_{m}^{GGE}$ would hold only if one could neglect the off-diagonal terms in the determinant in (45), i.e. replace $\det_{p\times p}\nu_{\ell_{i},\ell_{j}}\to\nu_{\ell_{1}}\dots\nu_{\ell_{p}}$, as implied by the Cauchy-Binet identity. This occurs in some models, see a bosonic example in [40] which has translational invariance and the GGE was shown to hold. Here, connecting $R_{m}^{\rm di}$ to the GGE (possibly in an inhomogeneous version involving local Lagrange multipliers $f_{\ell}(x)$) remains open.

In conclusion we studied the evolution of the Wigner function for noninteracting fermions in $d=1$ after a quantum quench. At times such that the integrity of the Fermi surf is retained, we found that it is described by the same universal function as at equilibrium, although with a space time dependent width that we calculated. We obtained an exact candidate formula for the large time limit of the Wigner function for generic (e.g. non-harmonic) confining potentials, and showed that it agrees with the GGE. We calculated the diagonal approximation to the many point correlations. It would be interesting to explore $d>1$ and detect these predictions in cold atom experiments.

Note added: at completion of this work we became aware of the recent work in [57] which studies a related problem.