A coordinate Bethe ansatz approach to the calculation of equilibrium and nonequilibrium correlations of the one-dimensional Bose gas

The Lieb–Liniger model [1] describes a system of N indistinguishable bosons subject to a delta-function interaction potential in a periodic 1D geometry of length L. We work in units such that ${\hslash }=1$ and the particle mass $m=1/2$, and so the Hamiltonian of this system reads

$$\begin{eqnarray}&&\hat{H}=-\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{\partial }^{2}}{\partial {x}_{i}^{2}}+2c\displaystyle \sum _{i\lt j}^{N}\delta ({x}_{i}-{x}_{j}),\end{eqnarray} \tag{ 1 }$$

where c is the interaction strength. The coordinate Bethe ansatz yields eigenstates $| \{{\lambda }_{j}\}\rangle$ of Hamiltonian (1) with spatial representation [17]

$$\begin{eqnarray}{\zeta }_{\{{\lambda }_{j}\}}(\{{x}_{i}\}) & \equiv & \langle \{{x}_{i}\}| \{{\lambda }_{j}\}\rangle \\ & = & \;{A}_{\{{\lambda }_{j}\}}\displaystyle \sum _{\sigma }\ \mathrm{exp}\left[{\rm{i}}\displaystyle \sum _{m=1}^{N}{x}_{m}{\lambda }_{\sigma (m)}\right]\displaystyle \prod _{k\gt l}\left(1-\displaystyle \frac{{\rm{i}}c\;\mathrm{sgn}({x}_{k}-{x}_{l})}{{\lambda }_{\sigma (k)}-{\lambda }_{\sigma (l)}}\right),\end{eqnarray} \tag{ 2 }$$

where the rapidities ${\lambda }_{j}$ (or quasimomenta) are solutions of the Bethe equations

$$\begin{eqnarray}&&{\lambda }_{j}=\displaystyle \frac{2\pi }{L}{m}_{j}-\displaystyle \frac{2}{L}\displaystyle \sum _{k=1}^{N}\mathrm{arctan}\left(\displaystyle \frac{{\lambda }_{j}-{\lambda }_{k}}{c}\right).\end{eqnarray} \tag{ 3 }$$

The quantum numbers m_j are any N distinct integers (half-integers) in the case that N is odd (even) [2], and ${\sum }_{\sigma }$ denotes a sum over all $N!$ permutations $\sigma =\{\sigma (j)\}$ of $\{1,2,...,N\}$. The normalization constant reads [17]

$$\begin{eqnarray}&&{A}_{\{{\lambda }_{j}\}}=\displaystyle \frac{{\displaystyle \prod }_{k\gt l}({\lambda }_{k}-{\lambda }_{l})}{{[N!\mathrm{det}\{{M}_{\{{\lambda }_{j}\}}\}{\displaystyle \prod }_{k\gt l}[{({\lambda }_{k}-{\lambda }_{l})}^{2}+{c}^{2}]]}^{1/2}},\end{eqnarray} \tag{ 4 }$$

where ${M}_{\{{\lambda }_{j}\}}$ is the N × N matrix with elements

$$\begin{eqnarray}{[{M}_{\{{\lambda }_{j}\}}]}_{{kl}} & = & {\delta }_{{kl}}\left(L+\displaystyle \sum _{m=1}^{N}\displaystyle \frac{2c}{{c}^{2}+{({\lambda }_{k}-{\lambda }_{m})}^{2}}\right)-\displaystyle \frac{2c}{{c}^{2}+{({\lambda }_{k}-{\lambda }_{l})}^{2}}.\end{eqnarray} \tag{ 5 }$$

The rapidities determine the total momentum $P={\sum }_{j=1}^{N}{\lambda }_{j}$ and energy $E={\sum }_{j=1}^{N}{\lambda }_{j}^{2}$ of the system in each eigenstate. The ground state of the system corresponds to the set of N rapidities that minimize E and constitute the (pseudo-)Fermi sea of the 1D Bose gas [17]. The Fermi momentum

$$\begin{eqnarray}&&{k}_{{\rm{F}}}=\displaystyle \frac{2\pi }{L}\displaystyle \frac{N-1}{2}\end{eqnarray} \tag{ 6 }$$

is the magnitude of the largest rapidity occurring in the ground state in the Tonks–Girardeau limit of strong interactions [3]. The only parameter of the Lieb–Liniger model in the thermodynamic limit is the dimensionless coupling $\gamma \equiv c/n$, where $n\equiv N/L$ is the 1D density. In finite systems, physical quantities also depend on the particle number N (see, e.g., section 3.3), whereas the length L of our system, and therefore also the density n, are arbitrary. Consequently, in this article we will specify both N and γ. Unless specified otherwise, we measure time in units of ${k}_{{\rm{F}}}^{-2}$, energy in units of ${k}_{{\rm{F}}}^{2}$, and length in units of ${k}_{{\rm{F}}}^{-1}$.

As the eigenstates $| \{{\lambda }_{j}\}\rangle$ form a complete basis [91] for the state space of the Lieb–Liniger model, the expectation value $\langle \hat{O}{\rangle }_{t}=\mathrm{Tr}\{\hat{\rho }(t)\hat{O}\}$ of an arbitrary operator $\hat{O}$ in a Schrödinger-picture density matrix $\hat{\rho }(t)$ can be expressed as a sum of matrix elements of $\hat{O}$ between the states $| \{{\lambda }_{j}\}\rangle$. In particular, in a pure state $| \psi (t)\rangle ={\sum }_{\{{\lambda }_{j}\}}{C}_{\{{\lambda }_{j}\}}(t)| \{{\lambda }_{j}\}\rangle$ we have

$$\begin{eqnarray}\langle \hat{O}{\rangle }_{t} & \equiv & \langle \psi (t)| \hat{O}| \psi (t)\rangle = & \displaystyle \sum _{\{{\lambda }_{j}\}}\displaystyle \sum _{\{{\lambda }_{j}^{\prime }\}}{C}_{\{{\lambda }_{j}^{\prime }\}}^{*}(t){C}_{\{{\lambda }_{j}\}}^{}(t)\langle \{{\lambda }_{j}^{\prime }\}| \hat{O}| \{{\lambda }_{j}\}\rangle ,\end{eqnarray} \tag{ 7 }$$

whereas in a statistical ensemble with density matrix ${\hat{\rho }}_{\mathrm{SE}}={\sum }_{\{{\lambda }_{j}\}}{\rho }_{\{{\lambda }_{j}\}}^{\mathrm{SE}}| \{{\lambda }_{j}\}\rangle \langle \{{\lambda }_{j}\}|$, we find

$$\begin{eqnarray}\langle \hat{O}\rangle & = & \displaystyle \sum _{\{{\lambda }_{j}\}}{\rho }_{\{{\lambda }_{j}\}}^{\mathrm{SE}}\langle \{{\lambda }_{j}\}| \hat{O}| \{{\lambda }_{j}\}\rangle .\end{eqnarray} \tag{ 8 }$$

In this article, we focus in particular on the normalized mth-order equal-time correlation functions

$$\begin{eqnarray}&&{g}^{(m)}({x}_{1},...,{x}_{m},{x}_{1}^{\prime },...,{x}_{m}^{\prime };t)\equiv \displaystyle \frac{\langle {\hat{{\rm{\Psi }}}}^{\dagger }({x}_{1})\cdots {\hat{{\rm{\Psi }}}}^{\dagger }({x}_{m})\hat{{\rm{\Psi }}}({x}_{1}^{\prime })\cdots \hat{{\rm{\Psi }}}({x}_{m}^{\prime })\rangle }{{[\langle \hat{n}({x}_{1})\rangle \cdots \langle \hat{n}({x}_{m})\rangle \langle \hat{n}({x}_{1}^{\prime })\rangle \cdots \langle \hat{n}({x}_{m}^{\prime })\rangle ]}^{1/2}},\end{eqnarray} \tag{ 9 }$$

where ${\hat{{\rm{\Psi }}}}^{(\dagger )}(x)$ is the annihilation (creation) operator for the Bose field and $\hat{n}(x)\equiv {\hat{{\rm{\Psi }}}}^{\dagger }(x)\hat{{\rm{\Psi }}}(x)$. Here and in the following we drop the time index t of the state vectors.

Since the Hamiltonian we consider in this article is translationally invariant along the periodic volume of length L, the mean density $\langle \hat{n}(x)\rangle \equiv n$ is constant in both time and space, and ${g}^{(m)}({x}_{1},...,{x}_{m},{x}_{1}^{\prime },...,{x}_{m}^{\prime };t)=\langle {\hat{{\rm{\Psi }}}}^{\dagger }({x}_{1})\cdots {\hat{{\rm{\Psi }}}}^{\dagger }({x}_{m})\hat{{\rm{\Psi }}}({x}_{1}^{\prime })\cdots \hat{{\rm{\Psi }}}({x}_{m}^{\prime })\rangle {/n}^{m}$. The correlation functions ${g}^{(m)}({x}_{1},...,{x}_{m},{x}_{1}^{\prime },...,{x}_{m}^{\prime };t)$ can therefore be expressed as the expectation values of the operators ${\hat{g}}^{(m)}({x}_{1},...,{x}_{m},{x}_{1}^{\prime },...,{x}_{m}^{\prime })\equiv {\hat{{\rm{\Psi }}}}^{\dagger }({x}_{1})\cdots {\hat{{\rm{\Psi }}}}^{\dagger }({x}_{m})\hat{{\rm{\Psi }}}({x}_{1}^{\prime })\cdots \hat{{\rm{\Psi }}}({x}_{m}^{\prime }){/n}^{m}$. We note that for the same reasons as above the matrix elements $\langle \{{\lambda }_{j}^{\prime }\}| {\hat{g}}^{(m)}({x}_{1},...,{x}_{m},{x}_{1}^{\prime },...,{x}_{m}^{\prime })| \{{\lambda }_{j}\}\rangle$ are invariant under global coordinate shifts $x\to x+d$ and thus, without loss of generality, we can set one of the spatial variables to zero. For the first-order correlation function, the matrix elements are

$$\begin{eqnarray}\langle \{{\lambda }_{j}^{\prime }\}| {\hat{g}}^{(1)}(0,x)| \{{\lambda }_{j}\}\rangle & \equiv & \langle \{{\lambda }_{j}^{\prime }\}| {\hat{{\rm{\Psi }}}}^{\dagger }(0)\hat{{\rm{\Psi }}}(x)| \{{\lambda }_{j}\}\rangle \\ & = & \displaystyle \frac{N}{n}\displaystyle \int {\rm{d}}{x}_{1}\cdots {\rm{d}}{x}_{N-1}\;{\zeta }_{\{{\lambda }_{j}^{\prime }\}}^{*}(0,{x}_{1},...,{x}_{N-1}){\zeta }_{\{{\lambda }_{j}\}}(x,{x}_{1},...,{x}_{N-1}).\end{eqnarray} \tag{ 10 }$$

The evaluation of the integral in equation (10) is complicated by the sign function in equation (2) and the associated nonanalyticities in ${\zeta }_{\{{\lambda }_{j}\}}(\{{x}_{i}\})$ where any two particle coordinates x_k and x_l coincide. However, we can use the Bose symmetry of the wave function ${\zeta }_{\{{\lambda }_{j}\}}(\{{x}_{i}\})$ to reexpress this matrix element as a sum of integrals

$$\begin{eqnarray}\langle \{{\lambda }_{j}^{\prime }\}| {\hat{g}}^{(1)}(0,x)| \{{\lambda }_{j}\}\rangle & = & \displaystyle \frac{N!}{n}\displaystyle \sum _{{\ell }=0}^{N-1}{\displaystyle \int }_{{{ \mathcal R }}_{N-1,{\ell }}(x)}{\rm{d}}{x}_{1}\cdots {\rm{d}}{x}_{N-1}\;\\ & & \times {\zeta }_{\{{\lambda }_{j}^{\prime }\}}^{*}(0,{x}_{1},...,{x}_{N-1})\;{\zeta }_{\{{\lambda }_{j}\}}({x}_{1},...,{x}_{{\ell }},x,{x}_{{\ell }+1},...,{x}_{N-1}),\end{eqnarray} \tag{ 11 }$$

over the ordered domains [10]

$$\begin{eqnarray}&&{{ \mathcal R }}_{M,j}(x):\quad 0\leqslant {x}_{1}\lt \cdots \lt {x}_{j}\lt x\lt {x}_{j+1}\lt \cdots \lt {x}_{M}\leqslant L.\end{eqnarray} \tag{ 12 }$$

Substituting the coordinate-space form (equation (2)) of the Lieb–Liniger eigenfunctions, we obtain

$$\begin{eqnarray}\langle \{{\lambda }_{j}^{\prime }\}| {\hat{g}}^{(1)}(0,x)| \{{\lambda }_{j}\}\rangle & = & \displaystyle \frac{N!}{n}{A}_{\{{\lambda }_{j}^{}\}}^{}{A}_{\{{\lambda }_{j}^{\prime }\}}^{*}\displaystyle \sum _{\sigma }\displaystyle \sum _{{\sigma }^{\prime }}\displaystyle \prod _{j\gt k}\left(1-\displaystyle \frac{{\rm{i}}c}{{\lambda }_{\sigma (j)}-{\lambda }_{\sigma (k)}}\right)\displaystyle \prod _{{j}^{\prime }\gt {k}^{\prime }}\left(1+\displaystyle \frac{{\rm{i}}c}{{\lambda }_{{\sigma }^{\prime }({j}^{\prime })}^{\prime }-{\lambda }_{{\sigma }^{\prime }({k}^{\prime })}^{\prime }}\right)\\ & & \times \displaystyle \sum _{{\ell }=0}^{N-1}\mathrm{exp}({\rm{i}}{\lambda }_{\sigma ({\ell }+1)}x){\displaystyle \int }_{{{ \mathcal R }}_{N-1,{\ell }}(x)}{\rm{d}}{x}_{1}\cdots {\rm{d}}{x}_{N-1}\mathrm{exp}\left({\rm{i}}\displaystyle \sum _{m=1}^{N-1}({\lambda }_{{\sigma }^{({\ell }+1)}(m)}-{\lambda }_{{\sigma }^{\prime }(m+1)}^{\prime }){x}_{m}\right),\end{eqnarray} \tag{ 13 }$$

where ${\sigma }^{({\ell }+1)}=(\sigma (1),...,\sigma ({\ell }),\sigma ({\ell }+2),...,\sigma (N))$. The matrix elements of the second-order correlation operator ${\hat{g}}^{(2)}(0,x)\equiv {\hat{{\rm{\Psi }}}}^{\dagger }(0){\hat{{\rm{\Psi }}}}^{\dagger }(x)\hat{{\rm{\Psi }}}(x)\hat{{\rm{\Psi }}}(0)/{n}^{2}$ are similarly given by

$$\begin{eqnarray}\langle \{{\lambda }_{j}^{\prime }\}| {\hat{g}}^{(2)}(0,x)| \{{\lambda }_{j}\}\rangle & = & \displaystyle \frac{N!}{{n}^{2}}{A}_{\{{\lambda }_{j}^{}\}}^{}{A}_{\{{\lambda }_{j}^{\prime }\}}^{*}\displaystyle \sum _{\sigma }\displaystyle \sum _{{\sigma }^{\prime }}\displaystyle \prod _{j\gt k}\left(1-\displaystyle \frac{{\rm{i}}c}{{\lambda }_{\sigma (j)}-{\lambda }_{\sigma (k)}}\right)\displaystyle \prod _{{j}^{\prime }\gt {k}^{\prime }}\left(1+\displaystyle \frac{{\rm{i}}c}{{\lambda }_{{\sigma }^{\prime }({j}^{\prime })}^{\prime }-{\lambda }_{{\sigma }^{\prime }({k}^{\prime })}^{\prime }}\right)\\ & & \times \displaystyle \sum _{{\ell }=0}^{N-2}\mathrm{exp}({\rm{i}}({\lambda }_{\sigma ({\ell }+2)}-{\lambda }_{{\sigma }^{\prime }({\ell }+2)}^{\prime })x){\displaystyle \int }_{{{ \mathcal R }}_{N-2,{\ell }}(x)}{\rm{d}}{x}_{1}\cdots {\rm{d}}{x}_{N-2}\\ & & \times \mathrm{exp}\left({\rm{i}}\displaystyle \sum _{m=1}^{N-2}({\lambda }_{{\sigma }^{(1,{\ell }+2)}(m)}-{\lambda }_{{{\sigma }^{\prime }}^{(1,{\ell }+2)}(m)}^{\prime }){x}_{m}\right),\end{eqnarray} \tag{ 14 }$$

where ${\sigma }^{(1,{\ell }+2)}=(\sigma (2),...,\sigma ({\ell }+1),\sigma ({\ell }+3),...,\sigma (N))$ and ${{\sigma }^{\prime }}^{(1,{\ell }+2)}$ is defined analogously in terms of the elements of ${\sigma }^{\prime }$. In the limit $x\to 0$ this expression simplifies somewhat, and in general the matrix elements of the local mth-order correlation operator ${\hat{g}}^{(m)}(0)\equiv {[{\hat{{\rm{\Psi }}}}^{\dagger }(0)]}^{m}{[\hat{{\rm{\Psi }}}(0)]}^{m}/{n}^{m}$ are given by the expression

$$\begin{eqnarray}\langle \{{\lambda }_{j}^{\prime }\}| {\hat{g}}^{(m)}(0)| \{{\lambda }_{j}\}\rangle & = & \displaystyle \frac{N!}{{n}^{m}}{A}_{\{{\lambda }_{j}^{}\}}^{}{A}_{\{{\lambda }_{j}^{\prime }\}}^{*}\displaystyle \sum _{\sigma }\displaystyle \sum _{{\sigma }^{\prime }}\displaystyle \prod _{j\gt k}\left(1-\displaystyle \frac{{\rm{i}}c}{{\lambda }_{\sigma (j)}-{\lambda }_{\sigma (k)}}\right)\displaystyle \prod _{{j}^{\prime }\gt {k}^{\prime }}\left(1+\displaystyle \frac{{\rm{i}}c}{{\lambda }_{{\sigma }^{\prime }({j}^{\prime })}^{\prime }-{\lambda }_{{\sigma }^{\prime }({k}^{\prime })}^{\prime }}\right)\\ & & \times \displaystyle \sum _{{\ell }=0}^{N-m}{\displaystyle \int }_{{{ \mathcal R }}_{N-m}}{\rm{d}}{x}_{1}\cdots {\rm{d}}{x}_{N-m}\mathrm{exp}\left({\rm{i}}\displaystyle \sum _{n=1}^{N-m}({\lambda }_{\sigma (m+n)}-{\lambda }_{{\sigma }^{\prime }(m+n)}^{\prime }){x}_{n}\right),\end{eqnarray} \tag{ 15 }$$

where the domain ${{ \mathcal R }}_{M}:0\leqslant {x}_{1}\lt {x}_{2}\lt \cdots \lt {x}_{M}\leqslant L$. We note, moreover, that equations (13)–(15) include as degenerate cases the diagonal matrix elements (see [10]) appropriate to the calculation of correlations in the ground state (section 3) and in statistical ensembles (section 4).

The calculation of correlation functions from equations (13)–(15) involves the evaluation of integrals of the general form

$$\begin{eqnarray}{\displaystyle \int }_{{{ \mathcal R }}_{M,{\ell }}(x)}{\rm{d}}{x}_{1}\cdots {\rm{d}}{x}_{M}\mathrm{exp}\left({\rm{i}}\displaystyle \sum _{m=1}^{M}{\kappa }_{m}{x}_{m}\right) & = & {\displaystyle \int }_{x}^{L}{\rm{d}}{x}_{M}{{\rm{e}}}^{{\rm{i}}{\kappa }_{M}{x}_{M}}{\displaystyle \int }_{x}^{{x}_{M}}{\rm{d}}{x}_{M-1}{{\rm{e}}}^{{\rm{i}}{\kappa }_{M-1}{x}_{M-1}}\cdots {\displaystyle \int }_{x}^{{x}_{{\ell }+2}}{\rm{d}}{x}_{{\ell }+1}{{\rm{e}}}^{{\rm{i}}{\kappa }_{{\ell }+1}{x}_{{\ell }+1}}\\ & & \times {\displaystyle \int }_{0}^{x}{\rm{d}}{x}_{{\ell }}{{\rm{e}}}^{{\rm{i}}{\kappa }_{{\ell }}{x}_{{\ell }}}{\displaystyle \int }_{0}^{{x}_{{\ell }}}{\rm{d}}{x}_{{\ell }-1}{{\rm{e}}}^{{\rm{i}}{\kappa }_{{\ell }-1}{x}_{{\ell }-1}}\cdots {\displaystyle \int }_{0}^{{x}_{2}}{\rm{d}}{x}_{1}{{\rm{e}}}^{{\rm{i}}{\kappa }_{1}{x}_{1}},\end{eqnarray} \tag{ 16 }$$

where (for the repulsive interactions $c\gt 0$ considered in this article) the ${\kappa }_{m}$ are real numbers. A single closed form for this integral does not exist, as in general one or more ${\kappa }_{m}$ may vanish, and this must be handled separately from the case of ${\kappa }_{m}\ne 0$. However, given knowledge of the particular sets of rapidities $\{{\lambda }_{j}\}$ and $\{{\lambda }_{j}^{\prime }\}$ (and permutations σ and ${\sigma }^{\prime }$), and thus of the locations of zero exponents ${\kappa }_{m}=0$ in equation (16), each individual integral of this form can be reduced to an algebraic expression in terms of $\{{\kappa }_{m}\}$. More specifically, each successive integration $\int {\rm{d}}{x}_{m}$ yields a term (involving, in general, ${x}_{m+1}$) arising from the primitive integral [92]

$$\begin{eqnarray}\displaystyle \int {\rm{d}}x\ {x}^{p}{{\rm{e}}}^{{\rm{i}}{kx}} & = & -{({\rm{i}}/k)}^{p+1}{\rm{\Gamma }}(p+1,-{\rm{i}}{kx})\\ & = & -p!{({\rm{i}}/k)}^{p+1}{{\rm{e}}}^{{\rm{i}}{kx}}\displaystyle \sum _{s=0}^{p}\displaystyle \frac{{(-{\rm{i}}{kx})}^{s}}{s!},\end{eqnarray} \tag{ 17 }$$

in the case that ${\kappa }_{m}$ is nonzero, or from $\int {\rm{d}}x\;{x}^{p}$ otherwise. In our calculations, the construction of algebraic expressions for the integrals occurring in equations (13)–(15) in terms of the rapidities ${\lambda }_{j}$ is efficiently performed by a simple computer algorithm that accounts for and combines the symbolic terms that arise from these successive reductions. We note that, e.g., each matrix element $\langle \{{\lambda }_{j}^{\prime }\}| {\hat{g}}^{(1)}(0,x)| \{{\lambda }_{j}\}\rangle$ is a sum of N integrals over $(N-1)$-dimensional domains and that the integrand in each case comprises ${(N!)}^{2}$ terms [10], illustrating the dramatically increasing computational cost of evaluating correlation functions with increasing N. Nevertheless, the explicit closed-form expression for the integral produced by our algorithm can be evaluated to obtain a numerically exact result by substituting in the values of the rapidities. The latter are obtained by solving equation (3) numerically using Newton's method, starting in the Tonks–Girardeau regime of strong interactions $\gamma \gg 1$ and iteratively progressing to smaller values of γ using initial guesses given by linear extrapolation of the solutions at stronger interaction strengths.

We note that this algorithmic approach also provides for the efficient and accurate calculation of the overlaps $\langle \{{\lambda }_{j}\}| \{{\mu }_{j}\}\rangle$ between eigenstates of Hamiltonian (1) corresponding to different values of γ, which we make use of in our analysis of nonequilibrium dynamics in section 4. In particular, the overlap between an arbitrary eigenstate $| \{{\lambda }_{j}\}\rangle$ of $\hat{H}$ at a finite interaction strength $\gamma \gt 0$ and the noninteracting ground state $| 0\rangle$, with constant spatial representation $\langle \{{x}_{i}\}| 0\rangle ={L}^{-N/2}$, is simply given by

$$\begin{eqnarray}&&\langle \{{\lambda }_{j}\}| 0\rangle =\displaystyle \frac{N!}{{L}^{N/2}}{A}_{\{{\lambda }_{j}\}}\displaystyle \sum _{\sigma }\displaystyle \prod _{j\gt k}\left(1+\displaystyle \frac{{\rm{i}}c}{{\lambda }_{\sigma (j)}-{\lambda }_{\sigma (k)}}\right){\displaystyle \int }_{{{ \mathcal R }}_{N}}{\rm{d}}{x}_{1}\cdots {\rm{d}}{x}_{N}\mathrm{exp}\left(-{\rm{i}}\displaystyle \sum _{n=1}^{N}{\lambda }_{\sigma (n)}{x}_{n}\right),\end{eqnarray} \tag{ 18 }$$

which can easily be evaluated semi-analytically using our algorithm. In practice we find that the results we obtain for the overlaps from our evaluation of equation (18) agree with the recently derived closed-form expressions for these quantities [84, 93–95], which imply in particular that $\langle \{{\lambda }_{j}\}| 0\rangle \propto 1/{\lambda }_{j}^{2}$ as any ${\lambda }_{j}\to \infty$.

We begin by considering the first-order correlation function ${g}^{(1)}(x)\equiv {g}^{(1)}(0,x)$ in the ground state of the Lieb–Liniger model. In figure 1(a) we plot ${g}^{(1)}(x)$ for N = 7 particles for a range of interaction strengths γ, which exhibits the expected decrease in spatial phase coherence with increasing γ [16]. As is well known, true long-range order, ${\mathrm{lim}}_{x\to \infty }\;{g}^{(1)}(x)={n}_{0}\gt 0$ [98, 99], is prohibited in an interacting homogeneous 1D Bose gas in the thermodynamic limit, even at zero temperature (see [97] and references therein). Indeed the Lieb–Liniger system is quantum critical at zero temperature, and the asymptotic long-range behavior of ${g}^{(1)}(x)$ is a power-law decay (so-called quasi-long-range order) [17].

Zoom In Zoom Out Reset image size

This power-law scaling of ${g}^{(1)}(x)$ is only expected to be realized at separations x large compared to the healing length $\xi =1/\sqrt{\gamma }$ and, in a finite periodic geometry such as we consider here, is curtailed by the finite extent L of the system (see, e.g., [16]). Indeed, for $\gamma =0.1$, the power-law decay is not visible in our finite-sized calculation, although as the interaction strength γ increases ${g}^{(1)}(x)$ exhibits behavior consistent with power-law decay over an increasingly large range of x, see figure 1(a). In particular, for $\gamma \gtrsim 10$, our results for ${g}^{(1)}(x)$ seem to converge toward the asymptotic scaling of the Tonks–Girardeau limit (black dotted–dashed line) with increasing γ.

Due to the translational invariance of our system, the first-order correlations of the Lieb–Liniger ground state are encoded in the momentum distribution

$$\begin{eqnarray}&&\tilde{n}({k}_{j})=n{\displaystyle \int }_{0}^{L}{\rm{d}}x\ {{\rm{e}}}^{-{\rm{i}}{k}_{j}x}{g}^{(1)}(x),\end{eqnarray} \tag{ 19 }$$

which, in our finite periodic geometry, is only defined for discrete momenta ${k}_{j}=2\pi j/L$, with j an integer. In figure 1(b) we plot the momentum distributions $\tilde{n}({k}_{j})$ corresponding to the first-order correlation functions ${g}^{(1)}(x)$ shown in figure 1(a). The first feature that we note in figure 1(b) is that for all interaction strengths, $\tilde{n}(k)$ exhibits a power-law decay $\tilde{n}(k)\propto {k}^{-4}$ (dotted–dashed black line) at high momenta. This is a universal result for delta-function interactions in 1D [21, 89, 100] (and indeed also in higher dimensions [101]). The effects of the finite extent L of the system on the first-order correlations are again evident in this momentum-space representation. For $\gamma =0.1$, no deviation from the ∝k⁻⁴ scaling is observed for the smallest (nonzero) momenta k_j that can be resolved in the periodic geometry. For larger values of the interaction strength, $\tilde{n}(k)$ departs from the ∝k⁻⁴ scaling at increasingly large values of k with increasing γ, and develops a hump at momenta near ${k}_{{\rm{F}}}$ for $\gamma \gtrsim 10$ [89]. We note that although the small-k behavior of $\tilde{n}(k)$ tends towards the $\propto {k}^{-1/2}$ scaling exhibited by the Tonks–Girardeau gas in the thermodynamic limit, the rounding off of the power-law decay of ${g}^{(1)}(x)$ as $x\to L/2$ precludes $\tilde{n}(k)$ from reaching the known asymptotic $k\to 0$ behavior in our finite geometry.

In figure 1(c), we present the nonlocal second-order coherence ${g}^{(2)}(x)\equiv {g}^{(2)}(0,x,x,0)$, which provides a measure of density-density correlations, for N = 7 particles at a range of interaction strengths γ. In the limiting case of an ideal gas ($\gamma =0$), the ground state of the system is a Fock state of N particles in the zero-momentum single-particle mode, and the second-order coherence ${g}_{\gamma =0}^{(2)}(x)=1-{N}^{-1}$ (horizontal dashed line) is therefore independent of x. As the interaction strength γ is increased, the second-order coherence is increasingly suppressed at zero spatial separation and correspondingly enhanced at separations $x\gtrsim 2{k}_{{\rm{F}}}^{-1}$. Oscillations in ${g}^{(2)}(x)$ develop at finite x as the system enters the strongly interacting regime $\gamma \gg 1$ [9, 17] and, in particular, for $\gamma =100$ (dashed cyan line), our numerical results are practically indistinguishable from the exact Tonks–Girardeau limit result (solid black line) [3].

An alternative representation of the second-order correlations of the ground state is given by the static structure factor S(k), which is related to ${g}^{(2)}(x)$ by [11]

$$\begin{eqnarray}&&S({k}_{j})=1+n{\displaystyle \int }_{0}^{L}{\rm{d}}x\ {{\rm{e}}}^{-{\rm{i}}{k}_{j}x}[{g}^{(2)}(x)-1].\end{eqnarray} \tag{ 20 }$$

In figure 1(d) we present the structure factors S(k) corresponding to the correlation functions ${g}^{(2)}(x)$ shown in figure 1(c). For all values of γ, $S(0)=0$ due to particle-number conservation and translational invariance. In the ideal-gas limit (red circles) $S({k}_{j})=1$ for all nonzero k_j. In the opposite limit of a Tonks–Girardeau gas

$$\begin{eqnarray}{S}_{\gamma =\infty }({k}_{j})=\left\{\begin{array}{ll}\displaystyle \frac{| {k}_{j}| (1-{N}^{-1})}{2{k}_{{\rm{F}}}} & \quad \quad | {k}_{j}| \leqslant 2{k}_{{\rm{F}}}\\ 1 & \quad \quad | {k}_{j}| \gt 2{k}_{{\rm{F}}},\end{array}\right.\end{eqnarray} \tag{ 21 }$$

which tends, in the thermodynamic limit, to the well-known result (see, e.g., [9]) $S(k)=| k| /2{k}_{{\rm{F}}}$ for $| k| \leqslant 2{k}_{{\rm{F}}}$, and $S(k)=1$ for $| k| \gt 2{k}_{{\rm{F}}}$. Just as for ${g}^{(2)}(x)$, we observe that for $\gamma =100$ (cyan plus symbols), our numerical results for S(k) are almost identical to the known exact expression (equation (21)) for the Tonks–Girardeau limit (black crosses). For smaller values of γ our mesoscopic results for S(k) appear consistent with those of [22, 25], obtained using quantum Monte Carlo and algebraic-Bethe ansatz techniques, respectively.

We now focus in more detail on local correlation functions. We note that the local second-order coherence has recently been proposed as a measure of quantum criticality in the 1D boson system [102], while the local third-order correlations have received increasing attention both in theory [103] and experiment [47, 104–106]. The local fourth-order correlations for the Lieb–Liniger model have also been investigated [107]. In figure 2, we plot the local second-order coherence ${g}^{(2)}(0)$ (solid red line), together with the local third-order coherence ${g}^{(3)}(0)=\langle {[{\hat{{\rm{\Psi }}}}^{\dagger }(0)]}^{3}{[\hat{{\rm{\Psi }}}(0)]}^{3}\rangle /{n}^{3}$ (dotted green line), and the local fourth-order coherence ${g}^{(4)}(0)=\langle {[{\hat{{\rm{\Psi }}}}^{\dagger }(0)]}^{4}{[\hat{{\rm{\Psi }}}(0)]}^{4}\rangle /{n}^{4}$ (dashed blue line) for N = 7 particles and a broad range of interaction strengths γ. For comparison, we also plot the asymptotic results obtained in the Bogoliubov limit of weak interactions ($\gamma \to 0$) in the thermodynamic limit [12, 18] (left-hand dotted–dashed lines). The numerical results for small γ are broadly comparable to these thermodynamic-limit results. However, for the small particle numbers considered here, the suppression of ${g}^{(2)}(0)$, ${g}^{(3)}(0)$, and ${g}^{(4)}(0)$ due to interactions in the limit of small γ is overshadowed by the suppression due to the finite population of the system [20]. At larger γ, the effects of interactions dominate, and the numerical results converge closely to the appropriate strong-coupling expressions [18] (right-hand dotted–dashed lines). We note, therefore, that the local correlations of the Lieb–Liniger ground state, and particularly their scaling with γ, appear to be quite insensitive to the infrared cutoff imposed by the finite extent of our system in the strongly interacting regime $\gamma \gg 1$.

Zoom In Zoom Out Reset image size

The results we have obtained so far indicate that, as expected, the small size of our system leads to corrections to correlation functions as compared to their known asymptotic forms in the thermodynamic limit. However, our results also suggest that the effects of finite system size are comparatively less important for local correlations, particularly in the limit of large interaction strengths $\gamma \gg 1$. To further elucidate the potential significance of finite-size effects in our calculations of nonequilibrium dynamics [90], here we give a brief characterization of the dependence of correlation functions of the Lieb–Liniger ground state on the particle number N at a fixed value of the interaction strength γ.

Specifically we consider the case for $\gamma =10$, as this value places the system in the strongly interacting regime $\gamma \gg 1$ (which appears less sensitive to finite-size effects than the weakly interacting regime $\gamma \lesssim 1$), while still exhibiting significant deviations from the Tonks–Girardeau limit (see, e.g., [9]). Whereas elsewhere in this paper we quote momenta (lengths) in units of ${k}_{{\rm{F}}}$ (${k}_{{\rm{F}}}^{-1}$), in comparing results between systems with different particle numbers N we quote momenta (lengths) in units of $\pi n$ ${\rm{[}}{(\pi n)}^{-1}$], so as to avoid a potentially misleading dependence of the unit of length on N (see equation (6)).

In figure 3(a) we plot ${g}^{(1)}(x)$ for particle numbers $N=3,4,5,6,$ and 7. For small x, the curves fall nearly perfectly on one line. The same behavior can be observed for the large-k tail of the corresponding momentum distribution $\tilde{n}(k)$, which we plot in figure 3(b). Indeed, at larger momenta $k\gtrsim 2\pi n$, $\tilde{n}(k)$ appears to exhibit a rapid collapse to a single curve with increasing N [21, 109]. However, the differences in $\tilde{n}(k)$ are so small that they can not be seen in figure 3(b). For small momenta, our choice of units implies an increasing resolution with increasing particle number, specifically ${k}_{1}=2\pi /L\times {(\pi n)}^{-1}=2/N$. However, this lowest resolvable momentum seems to fall on one line for increasing particle number, indicating that the infrared behavior of large systems can be at least partly accessed by our mesoscopic system sizes.

Zoom In Zoom Out Reset image size

Luttinger-liquid theory predicts a long-range power-law decay ${g}^{(1)}(x)\propto | x{| }^{-1/2K}$, where the Luttinger parameter K can be calculated from the thermodynamic limit of the Bethe ansatz solution (see, e.g., [16, 17] and references therein). For our parameters we have K = 1.40, implying an asymptotic scaling ${g}^{(1)}(x)\propto | x{| }^{-0.357}$ (black dotted–dashed line in figure 3(a)). This corresponds to a power-law behavior $\tilde{n}(k)\propto | k{| }^{-1+1/2K}=| k{| }^{-0.643}$ [16] (dotted–dashed line in figure 3(b)) for small momenta. We note that this infrared scaling is a true many-body effect and as such does not show up for N = 2 particles. Indeed, one can show analytically that, for N = 2, the momentum distribution $\tilde{n}(k)\propto {({\lambda }_{1}^{2}-{k}^{2})}^{-2}$ and thus ${k}^{-4}$ is the highest power in the series expansion of $\tilde{n}(k)$.

In figure 3(c) we plot the nonlocal second-order coherence ${g}^{(2)}(x)$ for $\gamma =10$ and $N=3,4,5,6,$ and 7. The corresponding static structure factor S(k) is shown in figure 3(d). In figure 3(d) we also plot (black dotted–dashed line) the form of S(k) resulting from the phenomenological expression proposed in [108] (see also [110]). This expression involves the limiting dispersions and edge exponents of the Lieb–Liniger model, which we obtain by numerically solving the appropriate integral equations [1, 111]. We also plot the corresponding prediction for ${g}^{(2)}(x)$ (black dotted–dashed line) in figure 3(c). We note that the numerical results for our mesoscopic systems are, in general, rather close to the phenomenological thermodynamic-limit expressions even for the relatively small particle numbers considered here.

In figure 4 we plot the time evolution of the local second-order coherence ${g}^{(2)}(0,t)$ for N = 5 particles following quenches of the interaction strength from initial values ${\gamma }_{0}=0$ (red dotted line) and ${\gamma }_{0}=100$ (blue dashed line) to the common final value ${\gamma }^{*}$. For the quench from the noninteracting initial state (${\gamma }_{0}=0$), as time evolves the local second-order coherence decays from its initial value ${g}^{(2)}(0,t=0)=1-{N}^{-1}$ before settling down to fluctuate about the diagonal-ensemble expectation value ${g}_{\mathrm{DE}}^{(2)}(0)$ (horizontal dotted–dashed line). This behavior is consistent with results obtained for similar quenches of the interaction strength from zero to a positive value in [90]. For the quench from ${\gamma }_{0}=100$, the value of ${g}^{(2)}(0)$ in the initial 'fermionized' state is ${g}^{(2)}(0)\approx {10}^{-3}$. In this case ${g}^{(2)}(0,t)$ rises as time progresses, and then exhibits somewhat irregular oscillations about ${g}_{\mathrm{DE}}^{(2)}(0)$ (horizontal solid line). We observe that the decay (growth) of ${g}^{(2)}(0,t)$ to its diagonal-ensemble value and the onset of irregular oscillations about this value occur on comparable time scales in the two quenches.

Zoom In Zoom Out Reset image size

We note that the predictions of the diagonal ensemble for the local second-order coherence ${g}_{\mathrm{DE}}^{(2)}(0)$ are very similar for the two quenches, despite the significant difference between the values of ${g}^{(2)}(0)$ in the two initial states. However, they are clearly distinct—${g}_{\mathrm{DE}}^{(2)}(0)$ for the quench from the noninteracting state is in fact larger than that for the quench from the correlated state by an amount ≈0.0125, demonstrating that the system retains some memory of its initial state in the long time limit as is expected for an integrable system. We analyze this difference in more detail in section 4.3.

We now turn our attention to the time evolution of the full nonlocal second-order correlation function ${g}^{(2)}(x,t)$. In figure 5(a) we show the dependence of ${g}^{(2)}(x,t)$ on separation x for the quench from the noninteracting initial state at four representative times. (Note that the upper limit $x=2\pi {k}_{{\rm{F}}}^{-1}$ of the x axis in figure 5(a) corresponds to $x=L/2$ in the present case of N = 5 particles.) At t = 0 (horizontal solid line), the second-order coherence has the constant form of the noninteracting ground state. At short times (e.g., $t=0.01\;{k}_{{\rm{F}}}^{-2}$, red dashed line) a minimum in ${g}^{(2)}(x)$ develops at zero separation, together with the corresponding maximum required by the conservation of ${\displaystyle \int }_{0}^{L}{\rm{d}}x\ {g}^{(2)}(x,t)$ [66]. As time progresses a wave pattern of maxima and minima develops and propagates away from the origin (e.g., $t=0.1{k}_{{\rm{F}}}^{-2}$, green dotted line). By time $t=1\;{k}_{{\rm{F}}}^{-2}$ (blue dotted–dashed line), the distinct maxima and minima of ${g}^{(2)}(x,t)$ have broadened in such a way that they are no longer clearly distinguishable and the correlation function agrees reasonably well with its diagonal-ensemble form (black dotted–dashed line) for small separations $x\lesssim 0.25\times 2\pi {k}_{{\rm{F}}}^{-1}$. In figure 5(b) we show the full space and time dependence of ${g}^{(2)}(x,t)$ following a quench from ${\gamma }_{0}=0$, which gives a more complete picture of the development of a correlation wave at short length scales and its propagation to larger values of x as time progresses. The correlation wave we observe here is consistent with the results of previous investigations of the dynamics following the sudden introduction of repulsive interactions in an initially noninteracting gas [63, 64, 66, 78, 118].

Zoom In Zoom Out Reset image size

In figure 5(d) we plot the spatial form of ${g}^{(2)}(x,t)$ for the quench from ${\gamma }_{0}=100$ at the same four representative times considered in figure 5(a). Despite the obvious distinction that the initial (t = 0, solid gray line) correlation function is in the fermionized regime with ${g}^{(2)}(0)\ll 1$, the behavior of ${g}^{(2)}(x,t)$ for this quench is qualitatively similar to that observed for the quench from ${\gamma }_{0}=0$, in that at early times (e.g., $t=0.01{k}_{{\rm{F}}}^{-2}$, red dashed line), deviations from ${g}^{(2)}(x,t=0)$ occur only at small separations $x\ll 2\pi {k}_{{\rm{F}}}^{-1}$. Moreover, as time evolves and ${g}^{(2)}(0,t)$ increases towards ${g}_{\mathrm{DE}}^{(2)}(0)$, larger modulations of ${g}^{(2)}(x,t)$ about its initial functional form develop (e.g., $t=0.1{k}_{{\rm{F}}}^{-2}$, green dotted line). At later times (e.g., $t=1{k}_{{\rm{F}}}^{-2}$, blue dotted–dashed line), ${g}^{(2)}(x,t)$ is close to ${g}_{\mathrm{DE}}^{(2)}(x)$ at small separations $x\lesssim 0.25\times 2\pi {k}_{{\rm{F}}}^{-1}$, but exhibits large excursions away from it at larger x. In figure 5(e) we plot the full space and time dependence of ${g}^{(2)}(x,t)$ following the quench from ${\gamma }_{0}=100$. Although the behavior of ${g}^{(2)}(x,t)$ here obviously differs from that following a quench from the noninteracting initial state (figure 5(b)), with the 'fermionic' depression around x = 0 lessening rather than growing in magnitude, a similar pattern of propagating correlation waves in ${g}^{(2)}(x,t)$ can again be seen.

The correlation-wave pattern common to both quenches is more clearly exhibited by the change${g}^{(2)}(x,t)-{g}^{(2)}(x,0)$ in the correlation function following the quench, which we plot in figures 5(c) and (f). This representation of the postquench second-order coherence of the system reveals a remarkably similar pattern of propagating waves in both cases, although the maxima and minima of the two wave patterns are inverted relative to one another. Fitting a power law to the position x(t) of the first propagating extremum of each of the two correlation waves, we find $x\propto {t}^{0.516\pm 0.012}$ for the quench from ${\gamma }_{0}=0$ and $x\propto {t}^{0.496\pm 0.005}$ for the quench from ${\gamma }_{0}=100$, which we indicate by the solid black lines in figures 5(c) and (f). These power-law trajectories are consistent with the 'telescoping' $x\propto {t}^{1/2}$ behavior obtained for a quench $\gamma =0\to \infty$ in [63], and for quenches from finite repulsive interactions to the noninteracting limit in [119] (see also [120]). The small scale features on top of the main propagating extrema differ for the two quenches, with fast oscillations appearing more pronounced for the quench $\gamma =0\to {\gamma }^{*}$ in figure 5(c). Even though hardly visible in figure 5(f), they are still present for the quench from $\gamma =100\to {\gamma }^{*}$, but due to the different distribution of overlaps in the final basis compared to the quench from ${\gamma }_{0}=0$ (see section 4.3), they contain more high-frequency components and therefore the fine structure differs.

We now compare the time-averaged second-order correlation functions following the two quenches with the form of this function that would be obtained if, following the quench, the system relaxed to thermal equilibrium. As in [90] we make use of the canonical ensemble, for which the density matrix is given by

$$\begin{eqnarray}&&{\hat{\rho }}_{\mathrm{CE}}={Z}_{\mathrm{CE}}^{-1}\displaystyle \sum _{\{{\lambda }_{j}\}}{{\rm{e}}}^{-\beta {E}_{\{{\lambda }_{j}\}}}\ | \{{\lambda }_{j}\}\rangle \langle \{{\lambda }_{j}\}| ,\end{eqnarray} \tag{ 27 }$$

where the partition function ${Z}_{\mathrm{CE}}={\sum }_{\{{\lambda }_{j}\}}\mathrm{exp}(-\beta {E}_{\{{\lambda }_{j}\}})$. The inverse temperature β is determined implicitly by fixing the mean energy in the state ${\hat{\rho }}_{\mathrm{CE}}$ to the common postquench energy, i.e., $\mathrm{Tr}\{{\hat{\rho }}_{\mathrm{CE}}\hat{H}({\gamma }^{*})\}={E}_{0\to {\gamma }^{*}}$. The sum in equation (27), like that in equation (26), formally ranges over an infinite number of eigenstates. We therefore truncate this sum by applying a cutoff in energy, as described in appendix A.

In figure 6(a) we plot the second-order correlation function ${g}_{\mathrm{CE}}^{(2)}(x)=\mathrm{Tr}\{{\hat{\rho }}_{\mathrm{CE}}\;{\hat{g}}^{(2)}(0,x)\}$ in the canonical ensemble (black dotted–dashed line), along with the diagonal-ensemble predictions ${g}_{\mathrm{DE}}^{(2)}(x)$ for the quenches from ${\gamma }_{0}=0$ (red solid line) and from ${\gamma }_{0}=100$ (blue dotted line). For comparison we also plot the correlation functions in the initial states with ${\gamma }_{0}=0$ (horizontal line), ${\gamma }_{0}=100$ (gray dashed line), as well as the ground state for $\gamma ={\gamma }^{*}$ (solid black line). For the quench from ${\gamma }_{0}=0$, the time-averaged value ${g}_{\mathrm{DE}}^{(2)}(0)$ is smaller than the corresponding thermal value ${g}_{\mathrm{CE}}^{(2)}(0)$, consistent with the results of [84, 90, 112]. In fact ${g}_{\mathrm{DE}}^{(2)}(x)$ is suppressed below ${g}_{\mathrm{CE}}^{(2)}(x)$ over a range of separations $x\lesssim 0.4\times 2\pi {k}_{{\rm{F}}}^{-1}$. Correspondingly, ${g}_{\mathrm{DE}}^{(2)}(x)\gt {g}_{\mathrm{CE}}^{(2)}(x)$ at larger separations x due to particle number and momentum conservation. For the quench $\gamma =100\to {\gamma }^{*}$, the diagonal-ensemble coherence function ${g}_{\mathrm{DE}}^{(2)}(x)$ is similar in shape to that of the quench from ${\gamma }_{0}=0$. However, it is somewhat smaller at x = 0, and correspondingly larger at large x. This indicates some memory of the initial state preserved by the dynamics of the integrable Lieb–Liniger system [58, 85]. Despite these differences, on the whole both functions ${g}_{\mathrm{DE}}^{(2)}(x)$ are comparable to ${g}_{\mathrm{CE}}^{(2)}(x)$ (see also [66]). We note, however, that they are also both reasonably close to the ground state result for ${g}^{(2)}(x)$ at interaction strength ${\gamma }^{*}$ (solid black line), although the local value ${g}_{\mathrm{DE}}^{(2)}(0)$ for both quenches is much closer to the thermal value than the ground state value.

Zoom In Zoom Out Reset image size

Since the system is in its ground state before the quench for both ${\gamma }_{0}=0$ and ${\gamma }_{0}=100$, and the total momentum operator $\hat{P}$ commutes with the Hamiltonian, the postquench states at ${\gamma }^{*}$ only have support on eigenstates with total momentum P = 0. Furthermore, the spatially structureless initial state at ${\gamma }_{0}=0$ implies additional parity-invariance ($\{{\lambda }_{j}\}=\{-{\lambda }_{j}\}$) in Bethe rapidity space for the postquench eigenstates [93–95]. Thus an interesting question to ask is if we constructed a canonical density matrix (27) restricted to P = 0 states, or one further restricted to parity-invariant states (which are a subset of the P = 0 states), would these yield better agreement with the diagonal ensemble predictions for the quenches? We have performed these constructions with the temperature in both cases fixed via the postquench energy in the same way as for the canonical ensemble, see equation (27) and the following text.

In figure 6(b), we plot the resulting second-order correlation function ${g}_{\mathrm{CE}}^{(2)}(x)=\mathrm{Tr}\{{\hat{\rho }}_{\mathrm{CE}}\;{\hat{g}}^{(2)}(0,x)\}$ for the standard canonical ensemble (black dotted–dashed line), as well as in the restricted P = 0 ensemble (solid black line), and the parity-invariant ensemble (solid gray line). We also include the diagonal-ensemble predictions ${g}_{\mathrm{DE}}^{(2)}(x)$ for the quenches from ${\gamma }_{0}=0$ (red solid line) and from ${\gamma }_{0}=100$ (blue dotted line). It can be seen that the restricted ensembles give results for the correlation function that are quite close to the standard canonical ensemble, and are no closer to the diagonal ensemble results.

The relaxation of the nonlocal correlations ${g}^{(2)}(x,t)$ takes place on a similar time scale to that of the local coherence ${g}^{(2)}(0,t)$ for both of the quenches considered here. This should be contrasted with, e.g., the behavior following a quench from the noninteracting limit to $\gamma =100$ reported in [90], in which ${g}^{(2)}(0,t)$ decays rapidly and the development and propagation of correlation waves occurs over a significantly longer time scale. We identify the absence of a significant separation of the time scales of local and nonlocal evolution here as a consequence of the fact that only a small number of eigenstates contribute significantly to the postquench dynamics (see [90] and references therein). Indeed, we find that the purity ${{\rm{\Gamma }}}_{\mathrm{DE}}\equiv \mathrm{Tr}\{{({\hat{\rho }}_{\mathrm{DE}})}^{2}\}$ of the diagonal-ensemble density matrix takes values ≈0.52 for the quench $\gamma =0\to {\gamma }^{*}$ and ≈0.63 for the quench $\gamma =100\to {\gamma }^{*}$, indicating rather weak participation of the eigenstates $| \{{\lambda }_{j}\}\rangle$ in the dynamics. The difference in the purities can largely be attributed to the somewhat greater occupation of the ground state of $\hat{H}({\gamma }^{*})$ following the quench from the $\gamma =100$ initial state.

To further illustrate the difference in the final states, in figure 7 we plot the occupations of eigenstates with energy ${E}_{\{{\lambda }_{j}\}}$ for the quenches from ${\gamma }_{0}=0$ (red crosses) and ${\gamma }_{0}=100$ (blue squares). For the quench from ${\gamma }_{0}=100$, significantly more eigenstates have occupations above a given threshold than in the case of ${\gamma }_{0}=0$, resulting in a much larger basis size in this case. However, the occupation of the ground state of $\hat{H}({\gamma }^{*})$ is somewhat larger for the quench from ${\gamma }_{0}=100$ than for ${\gamma }_{0}=0$, and the low-lying excited states are comparatively weakly occupied for ${\gamma }_{0}=100$, see figure 7(b). This result is reasonably intuitive, as the ground state for $\gamma ={\gamma }^{*}$ is moderately correlated, and will be more similar to the $\gamma =100$ than the $\gamma =0$ ground state. The distribution of normalization over eigenstates $| \{{\lambda }_{j}\}\rangle$ is thus more sharply 'localized' on the ground state in this case, resulting in the somewhat larger value of the purity ${{\rm{\Gamma }}}_{\mathrm{DE}}$ following this quench.

Zoom In Zoom Out Reset image size

For comparison, we also plot the occupations of the three ensembles introduced in section 4.2 in figure 7. The restrictions lead to a reduction in available eigenstates for any given energy-window, and correspondingly the temperature of the canonical ensemble is smaller than that of the P = 0 ensemble, which is in turn smaller than that of the parity-invariant ensemble. The occupations of eigenstates for the quench from ${\gamma }_{0}=0$ (red crosses) and from ${\gamma }_{0}=100$ (blue squares) are suggestive of power-law decay at high energies. For small energies on the other hand, figure 7(b) shows that the functional form is not incompatible with exponential decay.