Generalized hydrodynamics in the one-dimensional Bose gas: theory and experiments

It is instructive to start with the case of N = 2 particles on an infinite line. In first quantization, using center-of-mass and relative coordinates X = (x₁ + x₂)/2 and Y = x₁ − x₂, the Hamiltonian (6) splits into a sum of two independent one-body problems,

The eigenstates of the center-of-mass Hamiltonian are plane waves, and the Hamiltonian for the relative coordinate Y is that of a particle of mass 1/2 in the presence of a delta potential at Y = 0. Because of that delta potential, the first derivative of the wavefunction φ(Y) must have a discontinuity at Y = 0: φ′(0⁺) − φ′(0⁻) − c φ(0) = 0. Coming back to the original coordinates, one sees that the two-body wavefunction satisfies

The same condition holds for x₁ exchanged with x₂, since the wavefunction is symmetric. Thus the eigenstates of (8) are

corresponding to the eigenvalues . For θ₁ > θ₂, the two terms and correspond to the incoming and outgoing pairs of particles in a two-body scattering process. The ratio of their amplitudes is the two-body scattering phase,

An equivalent expression for that phase, often used in the literature and which we also use below, is ϕ(θ) = 2 arctan(θ/c) ∈ [−π, π].

It was pointed out by Eisenbud (1948) and by Wigner (1955) that the scattering phase may be viewed semiclassically as a ‘time delay’. Let us briefly sketch the argument of Wigner (1955). First, we note that, for a single particle, a simple substitute for a wavepacket is a superposition of two plane waves with momenta θ and θ + δθ,

Such a superposition evolves in time as e^{ixθ−itɛ(θ)} + e^{ix(θ+δθ)−itɛ(θ+δθ)}, where ɛ(θ) = θ²/2 is the energy. The center of this ‘wavepacket’ is at the position where the phases of the two terms coincide, namely the point where xδθ − t[ɛ(θ + δθ) − ɛ(θ)] = 0, which gives x ≃ vt with the group velocity v = dɛ/dθ = θ. So this is indeed a ‘wavepacket’ moving at speed θ. Next, consider two incoming particles in a state such that the center of mass (x₁ + x₂)/2 has momentum θ₁ + θ₂, while the relative coordinate x₁ − x₂ is in a ‘wavepacket’ moving at velocity (θ₁ − θ₂)/2,

According to equations (10) and (11), the corresponding outgoing state would be

Then, repeating the previous argument of phase stationarity, one finds that the relative coordinate is at position at time t. Since the center of mass is not affected by the collision and moves at the group velocity (θ₁ + θ₂)/2, we see that the position of the two semiclassical particles after the collision will be

where the scattering shift Δ(θ) is given by the derivative of the scattering phase,

The two particles are delayed: their position after the collision is the same as if they were late by a time δt₁ = Δ(θ₂ − θ₁)/v₁ and δt₂ = Δ(θ₂ − θ₁)/v₂ respectively.

For more particles, the eigenstates of the Hamiltonian (6) on the infinite line are Bethe states labeled by a set of N numbers , called the rapidities. In the domain x₁ < x₂ < ⋯ < x_N , the wavefunction is (Gaudin 2014, Korepin et al 1997, Lieb and Liniger 1963)

and it is extended to other domains by symmetry x_i ↔ x_j . Here the sum runs over all permutations σ of N elements (so there are N! terms) and (−1)^|σ| is the signature of the permutation. The momentum and energy of the eigenstate (17) are

The rapidities θ_a are conveniently thought of as the asymptotic momenta in an N-body scattering process. For θ₁ > θ₂ > ⋯ > θ_N , the combination of two terms

that appears in (17) can be viewed as the sum of incoming () and outgoing states () in an N-body scattering process, see figure 2. Their respective amplitude is a many-body phase, which depends on all incoming rapidities. Crucially, this many-body phase factorizes into a product of two-body scattering phases (11); this is a central property of all quantum integrable systems.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The fact that the rapidities are the asymptotic momenta in a scattering process implies that they can be measured by letting the bosons expand freely along the infinite line (Bolech et al 2012, Buljan et al 2008, Campbell et al 2015, Caux et al 2019, Del Campo 2008, Jukić et al 2008, Malvania et al 2021, Mei et al 2016, Minguzzi and Gangardt 2005, Rigol and Muramatsu 2005, Wilson et al 2020). Here we follow the argument of Campbell et al (2015).

Consider a state of N bosons confined to some interval around the origin at time t. This could be, for instance, the ground state in a trapping potential V(x), or some out-of-equilibrium state produced by some quench protocol, also in a trapping potential to ensure that the bosons are initially confined. This many-body state can be expanded on the basis of Bethe states,

where the Bethe states on the infinite line are normalized such that (assuming that both sets of rapidities are ordered, θ₁ > ⋯ > θ_N and ). The integral is restricted to the domain θ₁ > θ₂ > ⋯ > θ_N in the first line to avoid double-counting. Notice that, with the definition (17), the Bethe states are anti-symmetric under exchange of two rapidities θ_a ↔ θ_b . Then, plugging (17) into (20), and using this anti-symmetry, one obtains (Campbell et al 2015)

for x₁ < x₂ < ⋯ < x_N . This expression is particularly convenient for analyzing the expansion. When the trapping potential V(x) is switched off at time t, and the bosons are left to evolve freely along the infinite line, the probability of finding them at positions x₁, x₂, …, x_N after an expansion time t_exp is

From the second to the third line, we have used the stationary phase approximation, taking the limit t_exp → ∞ while keeping the ratios x₁/t_exp, …, x_N /t_exp fixed. The proportionality factor is fixed by imposing that .

In conclusion, we see from (22) that the joint probability distribution of the positions of the atoms after a large 1D expansion time directly reflects the distribution of rapidities in the state just before the expansion. This is very important because it means that the rapidities can be measured experimentally, by performing such 1D expansions. This was done experimentally for the first time by Wilson et al (2020). This experiment is discussed in section 3 below.

In the two previous subsections, we have focused on a finite number of bosons on an infinite line, corresponding to a vanishing density of particles. But, to understand the thermodynamic properties of the model, one needs to work with a finite density N/L. This can be done by imposing periodic boundary conditions, identifying the points x = 0 and x = L in the system. Imposing periodic boundary conditions on the Bethe wavefunction (17), i.e. , leads to the Bethe equations

where the two-body scattering phase ϕ(θ_a − θ_b ) is defined in equation (11). Taking the logarithm on both sides, one gets the following system of N coupled non-linear equations:

It is convenient to think of the numbers p_a as the momenta of N non-interacting fermions (with periodic or anti-periodic boundary conditions, depending on the parity of N). These fermion momenta have the following interpretation: for fixed N and L, one can adiabatically follow each eigenstate as one varies the repulsion strength c. In the infinite repulsion limit c → +∞, the bosonic wavefunction (17) is, up to multiplication by a sign ∏_a<b  sign(x_b − x_a ), equal to the Slater determinant of N non-interacting fermions, i.e. (Girardeau 1960). In that limit, the rapidities are nothing but the momenta of these non-interacting fermions: θ_a = p_a . Importantly, the fermions in the c → +∞ limit must obey the Pauli exclusion principle, so all momenta should be different: p_a ≠ p_b if a ≠ b. In the following we order both the rapidities and the fermion momenta as

It is natural to wonder what happens in the opposite limit of non-interacting bosons, c → 0. This is easily answered by introducing the ‘boson momenta’ , related to the fermion momenta p_a as

Notice that ; in particular, two or more boson momenta can coincide. Using the fact that , equation (24) is equivalent to

For c > 0, the ‘momenta’ are just another way of parameterizing the solutions of the Bethe equations; they should not be confused with the momenta of the atoms, which would be obtained by computing the momentum distribution , where is the Fourier mode of the creation operator Ψ^†(x). However, in the limit of vanishing repulsion c → 0, the rapidities are nothing but the boson momenta, . Moreover, in that limit, the Bethe wavefunction (17) is nothing but the permanent , i.e. the wavefunction of N non-interacting bosons. So, in that limit, the rapidities coincide with the atom momenta.

Away from these two limits, the rapidities θ_a correspond to an adiabatic interpolation between the non-interacting fermion (c → +∞) and boson (c → 0) momenta, obtained by solving the Bethe equations (24), see figure 3.

In general, the Bethe equations (24) cannot be solved analytically, but they can easily be solved numerically. One efficient way of doing this is to use the Newton–Raphson method.

The eigenstates of the Lieb–Liniger Hamiltonian (6) are Bethe states labeled by their sets of rapidities. This allows us to define a family of charge operators Q[f], diagonal in the eigenbasis and parameterized by functions , such that

Both the momentum operator and the Hamiltonian are of that form, with f(θ) = θ and f(θ) = θ²/2 respectively, see equation (18). It is the integrability of the model, reflected in the structure of the eigenstates (17), which allows us to consider the more general conserved charges (28). By construction, all these operators commute: [Q[f₁], Q[f₂]] = 0. In general, an explicit expression for Q[f] in second-quantized form (like Q[θ²/2] given by equation (6)) is not known, and typically regularization issues appear when one tries to write it (Davies 1990). Nevertheless, even in the absence of such direct expressions for the charges, the conserved charges Q[f] defined formally by equation (28) prove to be very useful. There are other ways to perform calculations with these charges, that do not require one to know their explicit second-quantized form, in particular the algebraic Bethe ansatz, see e.g. (Korepin et al 1997) for an introduction.

From their definition (28), one expects the Q[f] charges to be extensive with N, and to be the integral of a charge density

For f sufficiently regular, the charge density q[f](x) is sufficiently local, meaning that it vanishes quickly away from the point x. (In the hard-core limit g → +∞, this is a consequence of the Paley–Wiener theorem. At finite repulsion strength g, and more generally in interacting integrable models, the locality properties of charge densities are an advanced topic that is beyond the scope of this review, see e.g. (Doyon 2017, Ilievski et al 2016, Palmai and Konik 2018) for introductions.) By definition, the expectation value of the charge density in a Bethe state (normalized as ) is

It is independent of x because the Bethe state is translation invariant.

To the charge density q[f], one associates a current operator j[f] through the continuity equation,

As sketched in the introduction, continuity equations are the basic ingredient of hydrodynamics. To write useful hydrodynamic equations, however, one must be able to evaluate the currents in given stationary states. Until very recently, it was not known how to evaluate the expectation values of the current j[f]. However, thanks to developments in integrability, in particular in form factor techniques and the algebraic Bethe ansatz, a remarkable exact formula has just been discovered by Borsi et al (2020) for the expectation value of j[f] in a Bethe state (see also (Pozsgay 2020a, 2020b) for further developments in the context of spin chains, as well as the review article by Borsi et al (2021) in this volume),

Here ɛ′(θ) = θ is the derivative of ɛ(θ) = θ²/2, and is the Jacobian matrix of the transformation from the p_a to the θ_b defined by equation (24), known as the Gaudin matrix,

(where the p_a in equation (24) are no longer restricted to be in ). The Gaudin matrix is symmetric, , as a consequence of the fact that the scattering phase ϕ depends on the rapidities θ_a and θ_b only through the difference θ_a − θ_b , see equation (16).

Let us mention that the remarkable formula (33) is a particular case of a more general result, also obtained in (Borsi et al 2020, Pozsgay 2020a, 2020b). One can define generalized currents j[h, f](x) through the generalization of the continuity equation (31):

Then the general formula for the expectation value reads

and the above physical current is the special case h(θ) = ɛ(θ) = θ²/2. We note that such generalized currents had also been considered in (Castro-Alvaredo et al 2016) in the thermodynamic limit.

The discovery and proof of formula (34) or its generalization (35) required advanced techniques (Borsi et al 2020, Pozsgay 2020a, 2020b), however the result is simple and its physical interpretation is quite clear. The bth boson, with rapidity θ_b , carries an amount of charge density . In the absence of other particles, it would travel at the single-particle group velocity v_b = ∂ɛ(θ_b )/∂p(θ_b ) = ∂ɛ(θ_b )/∂θ_b (or its generalization , resulting in the current .

In the presence of other particles, the group velocity of the bth boson is modified. To compute it, one can consider a small variation of the fermion momentum p_b → p_b + δp_b in equation (24). This results in a small change of the total momentum δP = δp_b , and of the total energy (more generally, of the total charge ). Thus, the modified group velocity is , resulting in formula (32), or more generally , resulting in (35).

For further discussions of the physical interpretation of equations (32) and (35), see (Borsi et al 2020) and also (Bertini et al 2016, Bonnes et al 2014, Castro-Alvaredo et al 2016, Doyon 2019b, Doyon et al 2018), where similar discussions had been presented previously for the thermodynamic version of these formulas (see equation (46) below). See also the two reviews by Borsi et al (2021) and by Cubero et al (2021) in this volume.

So far we have focused on a finite number of bosons N, first on an infinite line, and then in a periodic box of length L. For hydrodynamics, one needs first to understand the thermodynamic properties of the system. In this subsection, we briefly review the techniques for taking the thermodynamic limit N, L → ∞, keeping the density of bosons n = N/L fixed.

The key idea is to focus on an infinite sequence of eigenstates of the Lieb–Liniger Hamiltonian (6), with L = N/n, such that the limit of the distribution of rapidities

is well defined and is a (piecewise) smooth function of θ. The thermodynamic properties of the system (such as its energy density, pressure, etc) then become particular functionals of that rapidity density ρ(θ), and the goal is to find these functionals and to evaluate them. In what follows, we write ‘lim_therm.’ for this limiting procedure.

For example, consider the expectation values of the charge densities (30): in the thermodynamic limit, these become

In particular, the density of particles, the momentum density and the energy density are, respectively, , and .

1.5.1. Thermodynamic form of the Bethe equations

Crucially, since all the states in the infinite sequence are Bethe states, each set of rapidities must satisfy the Bethe equations (24). To implement that constraint, it is customary to consider the set of fermion momenta associated to the set of rapidities , both of them ordered as in (25), and to define the density of states ρ_s(θ) as (see figure 4)

where the sequence of indices a in the rhs is chosen so that lim_therm. θ_a = θ. Because the fermion momenta p_a must satisfy the Pauli exclusion principle (they must all be different), it is clear that . Also, notice that, by definition, . Consequently, the Fermi occupation ratio

must always satisfy

Moreover, the rapidity density ρ(θ) and the density of states ρ_s(θ) are related by the thermodynamic version of the Bethe equation (24). Plugging equation (24) into the definition (38) leads to the constitutive equation

where Δ(θ − θ′) is the differential two-body scattering shift (16).

In practice, to construct interesting thermodynamic states, one can specify the Fermi occupation ratio ν(θ), and then use the constitutive equation (41) to reconstruct the rapidity density ρ(θ) and the density of states ρ_s(θ). One important example of this is the ground state of the Lieb–Liniger Hamiltonian, which corresponds to an occupation ratio which is a rectangular function: ν(θ) = 1 for θ ∈ [−θ_F, θ_F], and ν(θ) = 0 otherwise. Here θ_F is the Fermi rapidity, which is a function of the density of particles n. In that case, the constitutive equation becomes the Lieb equation (Lieb and Liniger 1963) (also known as the Love equation (Love 1949); for studies of this particular equation see e.g. (Lang et al 2017, Mariño and Reis 2019, Popov 1977, Prolhac 2017, Takahashi 1975). Another important example is that of a thermal equilibrium distribution ρ(θ) obtained by solving the Yang–Yang equation (equation (57) below).

In general, the constitutive equation cannot be solved analytically, however, since it is linear, it is easily solved numerically by discretizing the integral.

1.5.2. The dressing

In thermodynamic manipulations, it turns out that the following operation is ubiquitous: to a function f(θ), one has to associate its ‘dressed’ counterpart f^dr(θ), defined by the integral equation

Although it is not explicit in the notation, f^dr(θ) is always a functional of the rapidity distribution, through its dependence on the Fermi occupation ratio. For instance, with this definition, the constitutive equation (41) is recast as

where 1(θ) = 1 is the constant function.

Another example where the dressing (42) pops out is in manipulations that involve the Gaudin matrix. This is important for this review article, because to establish hydrodynamic equations one needs the thermodynamic limit of the expectation value of the current, see equation (32). The following identity holds:

where, again, the relation between θ in the lhs and the rapidity θ_a in the rhs is θ = lim_therm. θ_a . This identity is easily derived as follows. Using the definition (33),

where the ‘undressing’ is the inverse of the dressing, i.e. . Inverting this formula gives equation (44).

We will see a few more examples of physical quantities whose computation involves the dressing operation below. References where this operation is used extensively include, for example, the original derivation of the GHD equations in integrable quantum field theories (Castro-Alvaredo et al 2016), the calculation of Drude weights and other two-point correlations of charge and currents in the Lieb–Liniger model (Doyon and Spohn 2017a) and the inclusion of force fields (Doyon and Yoshimura 2017), adiabatically varying interactions (Bastianello et al 2019) or diffusive corrections (De Nardis et al 2018, 2019, Gopalakrishnan et al 2018) into the GHD equations.

1.5.3. Expectation values of the currents in the thermodynamic limit

We now present the central ingredient of GHD. The thermodynamic expectation value of the current j[f] (see subsection 1.4) is

where the effective velocity is a functional of the rapidity distribution defined by

with ɛ′(θ) = id(θ) = θ and 1(θ) = 1. The remarkable result (46) was first obtained by Bertini et al (2016), Castro-Alvaredo et al (2016), and it was the key observation that triggered all the later developments of GHD in quantum integrable systems. Bertini et al (2016) relied partially on (Bonnes et al 2014), where the formula for the effective velocity (47) had first appeared in the context of a quantum integrable system. In retrospect, the thermodynamic result (46) can be viewed as a consequence of the finite-size formula (32) of Borsi et al (2020), Pozsgay (2020a, 2020b), using the fact that the dressing is the thermodynamic limit of the Gaudin matrix, see equation (44). Historically though, the thermodynamic result was discovered before its finite-size counterpart. Since 2016, several works have aimed to establish the validity of the thermodynamic formula (46) in various models, by relying on various approaches. Let us mention the form factor approaches of Cubero (2020), Cubero and Panfil (2020), Vu and Yoshimura (2019) for quantum field theories, and arguments based on the symmetry of the charge–current correlations (Yoshimura and Spohn 2020) or exact results in the classical integrable model of the Toda chain (Bulchandani et al 2019, Cao et al 2019, Doyon 2019a, Spohn 2020). We also refer to the two reviews on this topic by Cubero et al (2021) and by Borsi et al (2021) in this volume.

The effective velocity (47) solves the equation

This is analogous to equation (4) in the introduction, which defines the effective velocity in the hard rod gas. The main difference is that the scattering shift Δ(θ − θ′) is now rapidity-dependent, while in the hard rod gas Δ is a constant equal to minus the diameter of the balls. The physical interpretation of equation (48) is analogous to that shown in figure 1: a ‘tracer’ quasiparticle with rapidity (asymptotic momentum) θ, which would normally travel at constant speed θ in a vacuum, finds its velocity modified by the presence of a finite density ρ(θ′) of other quasiparticles. From time t to t + δt, the tracer typically scatters against a number δt × |v^eff[ρ](θ) − v^eff[ρ](θ′)|ρ(θ′) of quasiparticles with rapidity θ′. At each collision, the tracer is shifted backwards by an amount Δ(θ − θ′): this is the physical effect that is encoded by formula (48).

To check that the effective velocity (47) solves equation (48) as claimed, one can use the definition of the dressing and the constitutive relation:

In the previous subsection we illustrated how physical observables, such as the expectation values of charges and currents, become functionals of the rapidity distribution ρ(θ) in the thermodynamic limit. We did not explain how to construct physically meaningful rapidity distributions though (except for the ground state of the Lieb–Liniger Hamiltonian, for which ν(θ) is a rectangular function, see subsection 1.5).

For instance, what is the rapidity distribution corresponding to a thermal equilibrium state at non-zero temperature? This question was answered in the pioneering work of Yang and Yang (1969), which we now briefly review.

First, we observe that there are many different choices of sequences of eigenstates that lead to the same thermodynamic rapidity distribution (36). The description of the system in terms of a rapidity distribution ρ(θ) is only a coarse-grained description: one should think of the rapidity distribution ρ(θ) as characterizing a macrostate of the system, corresponding to a very large number of possible microstates . For thermodynamics, one needs to estimate the number of such microstates.

To estimate that number, one focuses on a small rapidity cell [θ, θ + δθ], which contains Lρ(θ)δθ rapidities. The Bethe equations (24) relate these rapidities to fermion momenta p_a in a momentum cell [p, p + δp], where δp/δθ ≃ 2πρ_s(θ), see equation (38). Importantly, the fermion momenta p_a satisfy the Pauli exclusion principle. Then the number of microstates is evaluated by counting how many configurations of mutually distinct Lρ(θ)δθ fermion momenta can fit into the box [p, p + δp]. Since the minimal spacing between two momenta is , the answer is

The total number of microstates is the product of all such configurations over all the rapidity cells [θ, θ + δθ]. Taking the logarithm, and replacing the sum by an integral over dθ, we obtain the Yang–Yang entropy

The notation indicates that the Yang–Yang entropy is a functional of ρ only, and not of ρ_s; this is because ρ_s must always be obtained from ρ by the constitutive equation (41).

Now let us consider the thermal equilibrium density matrix at temperature T,

where the sum runs over all eigenstates. In fact, a straightforward generalization consists of considering the GGE (Rigol et al 2008, 2007) density matrix

for some function f. We would like to compute expectation values w.r.t. this density matrix, for example

for some observable . When the observable is sufficiently local, it is believed that the expectation value does not depend on the specific microstate of the system, so that it becomes a functional of ρ in the thermodynamic limit,

This assumption is related to a ‘generalized eigenstate thermalization hypothesis’, see e.g. (Cassidy et al 2011, Dymarsky and Pavlenko 2019, He et al 2013, Pozsgay 2011, 2014, Vidmar and Rigol 2016). Under that assumption, one can replace the above sum over all eigenstates by a functional integral over the coarse-grained rapidity distribution ρ,

The functional integral is then dominated by the root distribution which minimizes a (generalized) free energy functional:

Using the definition of the Yang–Yang entropy and the constitutive equation (24), one obtains the following relation between the function f(θ) defining the diagonal density matrix (52) and the rapidity distribution ρ(θ) dominating the functional integral (55):

This equation is known as the (generalized) Yang–Yang equation, or (generalized) thermodynamic Bethe ansatz equation. Again, the term ‘generalized’ refers to the replacement of the thermal equilibrium density matrix by a GGE, (51) → (52), see e.g. (Caux and Konik 2012, Wouters et al 2014). Like most of the equations encountered so far, in general the Yang–Yang equation cannot be solved analytically, but it can be efficiently solved numerically, by iteration. In particular, this allows us to compute the rapidity distribution ρ(θ) at thermal equilibrium.

This is particularly useful in applications discussed later in this review, because in experiments one often assumes that the system is (at least initially) at thermal equilibrium.

For instance, using the Yang–Yang equation, it is possible to tabulate the equilibrium pressure as a function of the particle density n and the energy per particle e. To do this, one needs to first solve numerically equation (57) with f(θ) = (ɛ(θ) − μ)/T, and equation (41), to get ρ(θ) and ρ_s(θ). Then the equilibrium pressure is given by (Korepin et al 1997, Yang and Yang 1969)

where the free energy is F/L = ∫(ɛ(θ) − μ)ρ(θ)dθ − TS_YY[ρ]. This gives the thermodynamic equilibrium pressure at density n = ∫ρ(θ)dθ and energy per particle . Alternatively, the pressure can be identified with the momentum current j_P = j[id] (with id(θ) = θ) as we did in the introduction, see equation (3). Thus, according to equation (46), we must also have

The equivalence between the two formulas (58) and (59) follows from manipulations of the dressing operation (42), which we leave as an exercise to the interested reader. (Hint: with the definition of the effective velocity, formula (59) is equivalent to , while differentiating (57) w.r.t. θ and using the definition of the dressing operation leads to .)

It is not a priori clear whether GGEs are relevant in the context of an isolated Lieb–Liniger gas. This question is linked to the notion of relaxation in isolated many-body quantum systems, which has been at the heart of many studies in recent decades (Polkovnikov et al 2011).

It is now well established that GGEs are relevant to describing locally an isolated Lieb–Liniger gas after relaxation, see for instance the review articles (Caux 2016, Essler and Fagotti 2016, Vidmar and Rigol 2016).

Since this point is essential in GHD theory, we briefly recall the underlying physics. Let us consider a Lieb–Liniger gas, confined in a box-like potential, and let us assume it is initially in an out-of-equilibrium state that is a pure quantum state. This quantum state expands onto many Bethe ansatz eigenstates:

Typically, the rapidity distributions of the Bethe ansatz states involved in this expansion gather around a given averaged rapidity distribution. During the time evolution, the different Bethe ansatz states, each of which evolves with its own energy, will dephase. Because of this dephasing, the contribution of cross terms will vanish at long times when computing the mean value of an observable. Thus, the mean values of the observable will undergo a relaxation and take the asymptotic value

We then invoke the generalized eigenstate thermalization hypothesis which states that expectation values of a local observable do not depend on the specific Bethe ansatz state, but are a smooth functional of the rapidity distribution. Thus, expectation values of local observables are identical for all diagonal ensembles, provided they are peaked around a given rapidity distribution. One can choose the GGE corresponding to the correct rapidity distribution, as was done in equation (54). One then finds that

To compute local observables after relaxation, one can represent the whole isolated system of length L by any diagonal ensemble peaked around the correct rapidity distribution, and the GGE plays a special role in describing a subsystem of length l. If l is both much smaller than L and much larger than microscopic correlation lengths, then the subsystem is correctly described by a GGE, the GGE accounting properly for the fluctuations in the subsystem. This property was used when analyzing local fluctuation measurements (Armijo et al 2010, Jacqmin et al 2011). The physical picture supporting the relevance of the GGE to describe the small system is that the large system acts as a reservoir of rapidities for the subsystem.

In order to show that the GGE indeed describes a system in contact with reservoirs of rapidities, we propose the following picture. Let us first discretize the rapidity space: we split it into intervals [θ_i , θ_i + δθ], where i ∈ Z, θ_i = iδθ and δθ is much smaller than the scale of variation of ρ. We assume now that the system is, for each integer i, in contact with a reservoir of rapidities lying in the interval [θ_i , θ_i + δθ], and we label the reservoir with the integer i. Then, statistical mechanics tells us that the density matrix of the system is

where N_i is the number of rapidities of the state |{θ_a }⟩ lying in the interval [θ_i , θ_i + δθ], and f_i is the temperature parameter associated to the reservoir number i. Equation (63) is nothing other than the GGE ensemble given in equation (52), with a discretized function f: in each segment [θ_i , θ_i + δθ], f takes the constant value f_i .

So far we have seen that the thermodynamics of the Lieb–Liniger model is determined by the rapidity distribution ρ(θ), which parameterizes an infinite family of stationary states. This is in contrast with generic chaotic Galilean invariant gases, where the thermodynamic properties would depend only on the atomic density, the momentum density and the energy density. Despite the infinite-dimensional parameter space of stationary states, the Lieb–Liniger gas possesses a small number of asymptotic regimes where its description simplifies.

In this subsection, we review the three main asymptotic regimes of the Lieb–Liniger gas: the ideal Bose gas, the quasicondensate and the hard-core regimes. These three regimes arise when one compares the typical energy per atom e to two energy scales: the scattering energy mg²/ℏ² and the mean-field interaction energy gn (where n is the atom density):

• e ≫ mg²/ℏ², gn: ideal Bose gas regime
• e ≃ gn ≫ mg²/ℏ²: quasicondensate regime
• mg²/ℏ² ≫ e: hard-core regime.

In the following, we discuss the phenomenology of these three regimes.

1.8.1. Ideal Bose gas regime

When e ≫ mg²/ℏ², gn, the interactions are negligible, and the gas behaves like a gas of non-interacting bosons. The rapidity distribution coincides with the momentum distribution of the bosons, as discussed after equation (27). Alternatively, this can also be understood as a consequence of the fact that the momentum distribution of the ideal gas is preserved during a 1D expansion, and the fact that the rapidities are the asymptotic momenta after such an expansion, see subsection 1.2.

The GGE that describes the local properties of the gas takes the form of a Gaussian density matrix, , for some function h(p). Here Ψ_p = ∫e^−ipx Ψ(x)dx is the Fourier mode of the boson annihilation operator Ψ(x). One of the consequences of that general Gaussian form is that, because of Wick’s theorem, the two-body zero-distance correlation

Thus, the gas exhibits the bosonic bunching phenomenon: whenever a boson is found inside a small interval [x, x + dx], the probability of finding another boson in that same interval is enhanced (i.e. it is larger than n dx).

We stress that the ideal Bose gas regime is not restricted to the classical, or non-degenerate, limit. The population of some bosonic modes can be highly occupied, and thus the ideal Bose gas can be highly degenerate.

This is exemplified by the case of thermal equilibrium at temperature T and chemical potential μ, with μ < 0 for the ideal Bose gas. The distribution of bosons is given by the Bose–Einstein distribution , and there is a crossover between the classical gas (which corresponds to |μ| ≫ k_B T) and the degenerate ideal gas (|μ| ≪ k_B T). In the degenerate regime, for momenta , the momentum distribution is close to a Lorentzian of half-width at half-maximum . The atom density is then , so we can estimate that the gas is degenerate as long as , or, equivalently, as long as

On the other hand, the typical energy per particle is e ≃ |μ|. The above condition e ≫ gn, which ensures that the gas, although it is degenerate, is in the ideal Bose gas regime as opposed to the quasicondensate regime, then reads

As long as both conditions (65) and (66) are fulfilled, the thermal equilibrium gas is not in the quasicondensate regime, even though it is highly degenerate. We note that the condition (66) was first established by Kheruntsyan et al (2003), who estimated the effects of interactions on g⁽²⁾(0) perturbatively, and required that they remain small.

1.8.2. The quasicondensate regime

This regime is reached when and the typical energy per atom e stays close to its value in the ground state, e ≃ gn/2. It is characterized by very small density fluctuations, with

Correlations are weak in this regime: the probability of finding an atom in a small interval [x, x + dx] is barely affected by the presence of another atom in this interval.

A good description of the gas in that regime is provided by Bogoliubov theory, or more precisely the extension of Bogoliubov theory to quasicondensates (Mora and Castin 2003). This approach assumes a phase-density representation of the bosonic field: one writes the atomic field Ψ(x) as where ϕ and δn are the phase and density fluctuation fields, which fulfill [δn(x), ϕ(x′)] = iδ(x − x′). This approach is a coarse-grained approximation, valid for length scales much larger than the interparticle distance. The Bogoliubov approximation assumes small density fluctuations, δn(x) ≪ n, and small phase gradient, ∂ϕ(x)/∂x ≪ n. Inserting this phase–density representation into the Hamiltonian (6), one finds to second order:

This quadratic Hamiltonian allows us to determine the quantum fluctuations around the classical profile, which solves the Gross–Pitaevskii equation, i.e. n = N/L = μ/g where μ > 0 is the chemical potential. It is easily diagonalized by a Bogoliubov transformation. Defining the bosonic mode such that [B(x), B^†(x′)] = δ(x − x′), and its Fourier transform with , one finds that the quadratic Hamiltonian becomes, up to a constant term,

where we have used μ = gn. Then the Bogoliubov transformation

with and , where , gives

with a dispersion relation .

The Bogoliubov model (70) is obviously an integrable model, since it amounts to a collection of independent harmonic modes, and its integrals of motion are the population in each mode. Making the link between the Bogoliubov modes and the rapidities in the Bethe ansatz solution of the Hamiltonian (6) is, however, a difficult task. In his seminal work, Lieb (1963) described this link for states close to the ground state. He identified the so-called ‘Lieb-I excitation’ (or ‘particle excitation’) branch to the Bogoliubov modes. In a more recent investigation, Ristivojevic (2014) found that this holds in fact only for large enough momenta. Thus, to our knowledge, making the connection between rapidities and Bogoliubov modes precise remains an open problem. The difficulty of this problem is related to the difficulty of developing schemes for expansions of the thermodynamic form of the Bethe equations (41) at small γ, see e.g. (Lang et al 2017, Mariño and Reis 2019, Popov 1977, Prolhac 2017, Takahashi 1975).

1.8.3. Hard-core regime

The hard-core regime is reached when the scattering energy mg²/ℏ² is much larger than all other intensive energy scales in the system. This is equivalent to taking g → +∞. Then, in all two-body scattering processes, the scattering phase factor (11) is one, and the scattering shift (16) vanishes. In that regime, two atoms can never be at the same position, which results in

In that regime, when a boson is found in a small interval [x, x + dx], then the probability of finding another one in the same interval is zero. This property reflects the Pauli principle satisfied by non-interacting fermions, which are related to the hard-core bosons by the non-local transformation

This transformation, closely related to the Jordan–Wigner transformation between lattice hard-core bosons and spin chains of spin-1/2, is defined such that the fermion creation/annihilation operators satisfy the canonical anticommutation relation . In terms of these fermions, the Lieb–Liniger Hamiltonian (6) with g → +∞ becomes

This is the Hamiltonian of a non-interacting Fermi gas. The identification of hard-core bosons with non-interacting fermions remains valid in the presence of an external potential V(x).

The Bethe wavefunction is, up to a sign ∏_a<b  sign(x_b − x_a ), the Slater determinant of N non-interacting fermions (Girardeau 1960), and the rapidities are simply the momenta of the underlying non-interacting fermions, as discussed below equation (24). The GGEs, which describe relaxed states and which are characterized by the rapidity distribution, correspond, for the fermionic gas, to GGE states that are obtained as a product of the Gaussian density matrix for each momentum state, , where is the Fourier mode of the fermion creation operator defined equation (72). (To be more precise, for finite L the fermions obey either periodic or anti-periodic boundary conditions depending on the parity of the particle number N, so the GGE is rather of the form , where P_{N even} and P_{N odd} are projectors onto the even and odd sectors. This complication can usually be omitted when one is interested in expectation values of local observables. For an example where it cannot be omitted, see e.g. (Bouchoule et al 2020) where the change of boundary conditions plays a key role in determining the effect of particle losses on the resulting GGE.)

The number of works that have exploited the mapping from hard-core bosons to non-interacting fermions is too large to allow us to review them all here. We simply mention a few such works that are representative in that they illustrate the typical calculations that can be done in that regime. For instance, early studies of the momentum distribution (Lenard 1964, Vaidya and Tracy 1979) revealed the presence of algebraically decaying ground state correlations, as well as the presence of tails in the momentum distribution of the atoms decaying as 1/p⁴ (Minguzzi et al 2002, Rigol and Muramatsu 2004). As is often the case with the 1D Bose gas, these observations remain valid beyond the hard-core regime, see e.g. the review articles (Cazalilla 2004, Cazalilla et al 2011) for long-range correlations, or (Olshanii and Dunjko 2003) about the tails of the momentum distribution. Many advances have been made in out-of-equilibrium quantum dynamics thanks to the study of the hard-core limit, for instance the early works on GGEs (Rigol et al 2007) or, more recently, investigations of trap releases (Collura et al 2013a, 2013b) and equilibration toward a GGE, or Floquet dynamics in harmonic traps (Scopa and Karevski 2017, Scopa et al 2018). For more references on hard-core bosons, we refer the reader to the bibliographies of these papers.

1.8.4. The thermal equilibrium phase diagram

In general, because of its integrability, an isolated Lieb–Liniger gas has no reason to be described by a thermal equilibrium state. Even after relaxation, the system is expected to be described by a GGE parameterized by a whole function (the rapidity distribution, see subsections 1.6 and 1.7), rather than by a Gibbs ensemble parameterized by only two parameters: the atom density n and the energy density e. However, in experiments, weak perturbations violate integrability, for instance the presence of transversely excited states (Li et al 2020, Mazets et al 2008) or the longitudinal potential (Bastianello et al 2020a). The Lieb–Liniger GGE will then be observable on an intermediate time scale, long compared to the relaxation time of the Lieb–Liniger model, but short enough that the effect of integrability breaking perturbations is still negligible. The system on this intermediate time scale is called prethermalized (Berges 2004). At very long times, the integrability breaking mechanisms will induce relaxation toward a thermal equilibrium state. Here we discuss the thermal equilibrium behavior of the 1D Bose gas, following Gangardt and Shlyapnikov (2003a), Petrov et al (2000) and especially Kheruntsyan et al (2003).

Assuming that the homogeneous 1D Bose gas is at thermal equilibrium, its state is characterized by only two dimensionless parameters: the dimensionless interaction strength γ, and the dimensionless temperature t (not to be confused with time in this subsection)

To identify the different asymptotic regimes discussed above in the phase diagram (γ, t), one can rely on the fact that g⁽²⁾(0) allows one to distinguish between them (Kheruntsyan et al 2003), see equations (64), (67) and (71). One can compute g⁽²⁾(0) using the Hellmann–Feynman theorem,

where F is the free energy per unit length, which can be computed using Yang–Yang thermodynamics (Yang and Yang 1969), see subsection 1.6. As imposed by the Hohenberg–Mermin–Wagner theorem, there is no phase transition, however the three aforementioned asymptotic regimes appear in the phase diagram, separated by smooth crossovers, see figure 5.

Zoom In Zoom Out Reset image size

In figure 5, we represent the crossover between the ideal Bose gas regime and the quasicondensate regime, t ≃ γ^−3/2, see the condition (66). The dashed line is the quantum degeneracy condition t ≃ γ⁻², see (65). Above this line, the occupation numbers of single-particle quantum states are small: quantum effects are small and the gas behaves mainly as a classical gas. Below this line, the behavior depends on the regime. In the ideal Bose gas regime, low-energy single-particle states get highly populated. In the hard-core regime, the rapidity distribution, which corresponds to the momentum distribution of the equivalent Fermi gas, becomes close to that of a zero-temperature Fermi sea, namely a rectangular function.

1.8.5. The classical field approximation

We conclude this survey of the regimes of the Lieb–Liniger model by discussing the classical field approach, which, owing to its simplicity and its relevance to the description of the crossover between the quasicondensate and the degenerate ideal Bose gas, is a popular technique. In this approach, the quantization of the atomic field, i.e. the discrete nature of atoms, is ignored, and the physics boils down to that of a classical complex field ψ(x). The energy functional of this field is

and the Lagrangian is L = (iℏ/2)∫dx(ψ*∂ψ/∂t − ∂ψ*/∂tψ) − E[ψ], such that the time evolution of ψ is given by the Gross–Pitaevskii equation

The effect of an external potential is easily taken into account within this approach, by adding a term ∫dx  V(x)|ψ|² (or V(x)ψ x) to the energy functional (or to the Gross–Pitaevskii equation). This approach is expected to be meaningful in the degenerate ideal Bose gas regime, as well as in the high temperature quasicondensate regime, where the population of the modes is high. In particular, it captures the crossover between the ideal Bose gas regime and the quasicondensate regime. This approach has been used, at thermal equilibrium, to compute correlation functions (Bouchoule et al 2012, Castin 2004, Castin et al 2000, Jacqmin et al 2012) and full counting statistics (Arzamasovs and Gangardt 2019), and to investigate non-equilibrium dynamics (Bouchoule et al 2016, Thomas et al 2021). In the absence of an external potential, the Gross–Pitaevskii equation is the non-linear Schrödinger equation, a classical integrable model. The link between the classical integrals of motion and the rapidity distribution of the quantum model has been discussed in (Bettelheim 2020, Del Vecchio et al 2020).

The classical field approach is plagued by an overestimation of the role of high wavevector components of ψ: their thermal mean value scales as mk_B T/(2ℏ) in the classical field approach, instead of the expected Gaussian behavior . This can strongly affect some observables. In higher dimensions, this induces an ultraviolet divergence of the density, which leads to the well-known black body problem. In 1D, the density does not diverge within the classical field approximation, but other observables do, like the energy density. To cure this problem, refined classical field approaches have been developed, that include a cut-off (Blakie et al 2008, Cockburn et al 2011a, Cockburn et al 2011b, Davis et al 2001).

The idea of a hydrodynamic description of the 1D Bose gas does not date back to 2016; it is of course much older. One popular hydrodynamic approach in the atomic gas literature is to start from the Gross–Pitaevskii description of the weakly interacting gas at zero temperature. Writing the wavefunction of the quasicondensate as , and defining the velocity ), one gets the Madelung form of the Gross–Pitaevskii equation (Cazalilla et al 2011, Stringari 1996, 1998),

where is the pressure of the gas in the quasicondensate regime. Clearly, these two equations look like the first two Euler hydrodynamic equations (3) in the introduction, up to the so-called quantum pressure term . This term is beyond the Euler scale though: it involves higher-order derivatives, so in the Euler limit of a slowly varying density n(x), this term vanishes. The reason why there are only two equations in (79), instead of the three Euler equations in the introduction, is that the temperature is zero. Indeed, the third Euler equation in (3) can be recast as a conservation law for the entropy of the fluid, which is automatically satisfied at zero temperature because the entropy identically vanishes.

There have been several attempts at extending this description of the gas beyond the weakly interacting regime, see e.g. (Damski 2006, Kolomeisky et al 2000, Menotti and Stringari 2002). One idea that is often used (Damski 2006, Menotti and Stringari 2002, Öhberg and Santos 2002, Pedri et al 2003, Peotta and Di Ventra 2014, Sarishvili et al 2016) is to replace the pressure of the quasicondensate regime by the true pressure of the Lieb–Liniger model at zero temperature, calculated from equation (59). This gives a closed system of hydrodynamic equations for the gas that can be solved numerically. This approach can also be used at finite temperature (Bouchoule et al 2016, Doyon et al 2017, Schemmer et al 2019): in that case one numerically solves the three Euler hydrodynamic equations from the introduction with the equilibrium pressure at finite temperature.

This ‘conventional’ Euler hydrodynamic approach, which assumes local relaxation of the gas to a thermal equilibrium state, has been successfully applied in states not far from thermal equilibrium, see e.g. (Menotti and Stringari 2002, Öhberg and Santos 2002, Pedri et al 2003). However, it breaks down away from equilibrium. A good illustration of the problems that are typically encountered can be found in (Peotta and Di Ventra 2014), see figure 7. In that reference, the conventional hydrodynamic equations (79) are solved numerically and compared to numerically exact time-dependent density matrix renormalization group (t-DMRG) calculations for small numbers of bosons. The authors find that the hydrodynamic equations can describe the breathing of the atom cloud very well after a quench of the harmonic trapping frequency. However, those hydrodynamic equations are unable to describe a quench from a double-well to a harmonic potential away from the weakly interacting regime.

Zoom In Zoom Out Reset image size

The reason for the failure of this conventional hydrodynamic approach in the latter setup is that it develops a shock after some fraction of the oscillation period. When the atom cloud initially has two well-separated density peaks, some atoms from the left peak move to the right with large velocity, while other atoms from the right peak move to the left with large velocity. Then in the center of the cloud, the fluid is similar to a two-component fluid, with one component moving quickly to the right, the other to the left. Locally, this is a state that is very far from a thermal equilibrium state. The above approach, which enforces local thermal equilibrium, fails to capture that situation. Instead, the Euler-scale hydrodynamic equations (79) develop a shock, which is regulated by higher-order derivative terms, like the quantum pressure term. The solution after the shock depends very strongly on the details of the regularization, so that in general it loses its validity after the first shock. The results of Peotta and Di Ventra (2014) show that the simple insertion of the quantum pressure term into the zero-temperature hydrodynamic equations (79) does not provide the correct regularization for finite repulsion strength. Only in the limit of weak repulsion (g → 0), i.e. when equation (79) is mathematically equivalent to the standard Gross–Pitaevskii equation, does the quantum pressure term provide the right regularization. In that case, the Euler-scale shock is visible in the Gross–Pitaevskii solution through the appearance of strong oscillations of short wavelength in the density profile n(x, t), see e.g. (Simmons et al 2020). These oscillations are beyond the Euler scale, and one can in principle average them over fluid cells to recover the correct Euler-scale description after the shock. Bettelheim (2020) has shown that the Euler-scale dynamics obtained from the Gross–Pitaevskii equation in this way exactly coincides with GHD (at zero temperature and in the limit of weak repulsion).

Away from the Gross–Pitaevskii limit, the analysis carried out in (Doyon et al 2017) shows that the above approach is, in fact, well justified only at zero temperature and before the appearance of the first shock. Moreover, in that case, the conventional hydrodynamics (79) at the Euler scale (i.e. neglecting the quantum pressure term) turns out to be exactly equivalent to GHD. However, in any other situation, it is in principle not applicable, and it leads to quantitatively wrong results. This is illustrated in figure 7.

Contrary to previous hydrodynamic approaches to the 1D Bose gas, GHD is not based on the assumption of local thermal equilibrium, but only on local relaxation to a stationary state which, in general, is a GGE. This allows the description of situations that are very far from thermal equilibrium. Arguably, the most paradigmatic such out-of-equilibrium situation is the quantum Newton cradle setup of Kinoshita et al (2006). There, the atoms, which are initially at equilibrium in a harmonic potential, are suddenly given a large momentum ±q_Bragg by a Bragg pulse. Half of the atoms move to the right, the other half move to the left. Because of the harmonic trapping potential, the two packets of atoms oscillate in the trap, colliding twice during each oscillation cycle. In this subsection we review the theoretical work on this paradigmatic setup. The pioneering experiment of Kinoshita et al (2006), which motivated all these theoretical works, is discussed below in section 4.

We stress that, before the advent of GHD, a direct simulation of the quantum Newton cradle setup with experimentally realistic parameters (in particular, number of atoms N ∼ 10²–10³) was completely out of reach. This is of course because of the exponential growth of the Hilbert space in many-body quantum systems, which makes all direct approaches, such as an exact diagonalization of (a discretized version of) the Lieb–Liniger Hamiltonian (6), numerically intractable for more than a dozen atoms. Thanks to the discovery of GHD in 2016, this situation has now completely changed. Nowadays, it is very easy to model the 1D Bose gas in a quantitatively reliable way.

Bulchandani et al (2017) used GHD to model two packets of atoms colliding against each other on an infinite line. Then, a complete study of the Newton cradle setup, including the trapping potential that gives rise to oscillations of the packets and therefore multiple collisions, was performed by Caux et al (2019), see figure 8. The numerical solution of the GHD equation in that reference was obtained by a classical molecular dynamics simulation of the so-called flea gas model (Doyon et al 2018), which is an extension of the classical hard rod gas that incorporates the scattering shift discussed in subsection 1.1. There are other ways of numerically solving the GHD equations, which have been discussed, to some extent, in (Bastianello et al 2019, Bastianello et al 2020a, Bulchandani et al 2017, 2018, Doyon et al 2017, Møller and Schmiedmayer 2020, Møller et al 2021) (see also the appendix of the review by Bastianello et al (2021) in this volume). Among these works, we draw the reader’s attention in particular to the GHD code ‘iFluid’ of Møller and Schmiedmayer (2020), which is publicly available.

Zoom In Zoom Out Reset image size

In figure 8, one can observe the evolution of the rapidity distribution predicted by GHD in a harmonic potential (with oscillation period τ = 2π/ω), and also in a potential with a small anharmonicity . The initial state is constructed so as to mimic the effect of the Bragg pulse sequence that imparts the initial momentum to the atoms. Before the sequence, the gas is described in a hydrostatic (or local density) approximation by its local distribution of rapidities ρ(x, θ, t < 0) = ρ_thermal(x, θ), obtained by solving the Yang–Yang equation (57) for a thermal equilibrium distribution. Then momentum ±q_Bragg is imparted in a random fashion to all quasiparticles in the system. This results in a distribution of rapidities at t = 0

This simple ansatz for the initial state can be justified using the results of van den Berg et al (2016), which showed that the momentum distribution function of the bosons is affected in this way by a Bragg pulse, and to a good approximation the same holds for the rapidities. The rapidity distribution is evaluated at later times by solving the GHD equation (78). One observes that the two blobs, initially well separated in momentum space for sufficiently large q_Bragg, evolve by performing a deformed rotation-like movement around the origin of phase space. In the harmonic case, over the first two or three oscillation cycles, their evolution is not drastically affected by the collisions. However, at later times the two blobs ultimately merge due to inter-cloud interactions. With a small anharmonicity, the two blobs get deformed much more quickly and the distribution ρ(x, θ, t) gets more and more stirred up after a few periods. This dephasing effect would also be present for the single particle in an anharmonic trap, see e.g. (Bastianello et al 2017). Many-body dephasing is also present: without interactions, the original blobs would disintegrate into long spiraling filaments; instead, here the filaments merge and high-energy (longer-period) tails scatter to lower energies, leading to the reformation of new blobs.

Importantly, the ability to perform GHD modeling of the Newton cradle setup has opened the possibility of studying theoretically one fundamental question raised by the experiment of Kinoshita et al (2006). If one waits long enough, does the gas in the trap ultimately reach thermal equilibrium?

This question is non-trivial because, although the Lieb–Liniger gas is integrable, its integrability is broken by the trapping potential V(x). Yet, at the Euler scale, the potential V(x) varies very slowly compared to microscopic scales, so the breaking of integrability by the external potential is weak. It is not obvious how much of the original conservation laws should be reflected in the stationary state.

Cao et al (2018) studied this question for the classical hard rod gas in a trapping potential, and found that the answer is negative: the gas exhibits ‘incomplete thermalization’. It reaches a stationary state of the GHD equation (78), namely a rapidity distribution ρ_stat.(x, θ) which satisfies

but this distribution does not need to be a thermal equilibrium distribution in the trap. The possibility of GHD-stationary distributions of the form (81) that are not thermal has also been discussed in (Doyon and Yoshimura 2017). Caux et al (2019) arrived at the same conclusion for the 1D Bose gas in the Newton cradle setup, within the framework of the aforementioned flea gas model. In (Caux et al 2019), the existence of non-thermal stationary states was argued to be a consequence of the conservation of certain quantities S[f] under evolution generated by the Euler-scale GHD equation (78), even in the presence of an external potential V(x). The conservation of these quantities is incompatible with convergence toward thermal equilibrium. These quantities can be constructed out of the Fermi occupation ratio ν(x, θ, t), and read

for arbitrary functions f. We stress that these quantities are different from the standard conserved charges of the form (28). Instead, the quantities S[f] look more like generalizations of the Yang–Yang entropy: the Yang–Yang entropy, integrated over space, corresponds to the specific choice f(ν) = −ν  log ν − (1 − ν) log(1 − ν), see equation (49). The fact that

can be checked directly using the GHD equation (78), see (Caux et al 2019). We stress that the Yang–Yang entropy, and more generally the quantities S[f], are conserved only at the Euler scale. When higher-order terms are included into the hydrodynamic equations (78), as discussed in subsection 2.4.3 below, the entropy increases with time, and the other quantities S[f] are no longer constant. In particular, it has been argued recently by Bastianello et al (2020a) that the inclusion of a Navier–Stokes-like higher-order term in (78) does lead to ‘complete thermalization’, see subsection 2.4.3 below.

Let us briefly review some other physically relevant setups that have been investigated with the new toolbox provided by GHD.

De Nardis and Panfil (2018) considered the case of two atom clouds that are prepared at different temperatures T₁ ≠ T₂, put in the same trapping potential. At the junction between the two clouds, the local state displays an edge singularity in its response function and quasilong-range order.

Dubessy et al (2021) studied the Lieb–Liniger gas in an infinite flat box potential (figure 9). Initially the gas is in its ground state. At time t = 0, it is instantaneously boosted by a momentum k₀, a protocol that can be realized experimentally by phase imprinting. Then the particles in the gas start reflecting against the two infinite walls at x = 0 and x = L. This is modeled in GHD by the following boundary condition:

The same boundary condition holds at x = L. Dubessy et al (2021) used a trick to implement easily these boundary conditions: the system can be glued together with its mirror image, to give a periodic system of length 2L. The rapidity distribution in that periodic system of size 2L is related to that in the infinite box potential by

With this trick, Dubessy et al (2021) studied the formation of shock waves that oscillate in the box for several periods, with a period fixed by the sound velocity in the gas, see figure 9.

Zoom In Zoom Out Reset image size

Doyon et al (2017) pointed out that a drastic simplification of the GHD equations occurs when the gas is initially in its ground state. In an external potential V(x), the ground state is modeled by hydrostatics, or equivalently by the local density approximation (LDA), which gives the distribution of rapidities ρ(x, θ, t = 0). In the ground state, the Yang–Yang entropy of the gas vanishes. Then, since entropy is always conserved by Euler-scale hydrodynamic equations, it must vanish at all times. This puts very strong restrictions on the class of local stationary states that are explored by the system under GHD evolution: these local states must be either the ground state itself (up to a Galilean boost), or a ‘split Fermi sea’ (Eliëns 2017, Eliëns and Caux 2016, Fokkema et al 2014). Within the framework of GHD, this is easily understood by using the so-called convective form of the GHD equation (78), which gives the evolution of the Fermi occupation ratio ν(x, θ, t) = ρ(x, θ, t)/ρ_s(x, θ, t),

(It is easy to see that this form of the GHD equation is equivalent to (78). Plugging the constitutive relation (41) into (78), one finds that the first equation (78) satisfied by ρ(x, θ, t) is also satisfied by ρ_s(x, θ, t). This directly leads to (86) for the ratio ρ/ρ_s.)

In the ground state, the Fermi occupation ratio is either zero or one: ν(x, θ) = 1 if θ ∈ [−θ_F(x), θ_F(x)], and ν(x, θ) = 0 otherwise. Here θ_F(x) is a position-dependent Fermi momentum, which depends on the atom density n(x), see subsection 1.5. This specific form of ν(x, θ, t) is preserved under (86): at any time, ν(x, θ, t) is either zero or one. Consequently, the state of the system at time t is parameterized by a contour Γ_t in phase space (figure 9, left, and figure 13, top row), which separates the region where ν = 1 from the one where ν = 0, namely

Writing the contour as Γ_t = {(x_t (s), θ_t (s)); s ∈ [0, 2π)}, and plugging this into equation (86) one finds that its evolution equation reads

This ‘zero-entropy GHD’ (Doyon et al 2017) is very useful for numerical purposes, because it provides a very efficient way of solving the GHD equations: it is easier to compute the evolution of the contour Γ_t , rather than to simulate the evolution of the full distribution ρ(x, θ, t), even though both formulations are equivalent in the end. In figure 13 (top row), one can see the evolution of the contour Γ_t when the gas is suddenly released from its ground state in a double-well potential to a harmonic potential. After a fraction of the period of the harmonic trap frequency, one observes the appearance of local split Fermi seas. Hence, such a setup could not be described by the conventional hydrodynamic approaches of subsection 2.1, because the appearance of such multiple Fermi seas would translate into shock formations in those approaches. GHD, on the other hand, does not have shocks (Bulchandani 2017, Doyon et al 2017, El and Kamchatnov 2005) and remains valid after the appearance of multiple Fermi seas.

So far we have reviewed results obtained on the 1D Bose gas with the Euler-scale GHD equation (78). Since 2016, the original framework has been extended in several directions that could be relevant to the description of existing or future experiments. We now turn to these developments.

2.4.1. Adiabatically varying interactions

The original formulation of GHD of Bertini et al (2016), Castro-Alvaredo et al (2016) was for a time-independent translation invariant Hamiltonian acting on a spatially inhomogeneous state. In particular, no external potential term −(∂_x V)(∂_θ ρ) was included originally. Doyon and Yoshimura (2017) considered the addition of generalized potentials to the Hamiltonian H → H + ∫V_f (x)q[f](x)dx, where q[f](x) is the charge density of the charge operator (28). This includes, in particular, the case of a standard external potential H → H + ∫V(x)Ψ^†(x)Ψ(x)dx, corresponding to the simplest choice q[f] = 1, which results in the form (78) of the GHD equation with the acceleration term −∂_x V∂_θ ρ. This is the form of the GHD equation that is most relevant for the description of existing experiments. However, we stress that the results of Doyon and Yoshimura (2017) are in principle more general.

A further extension of the original GHD equation was obtained by Bastianello et al (2019), who considered the case of a non-uniform repulsion strength g(x, t) through the gas, see figure 10. Under the assumption of slow variation of g(x, t) (or of ) in position and time, they found the following additional term:

where we have dropped the dependence on x, θ and t for better readability, and where the functions F and Λ are

For a derivation of this result, see the original paper (Bastianello et al 2019).

Zoom In Zoom Out Reset image size

2.4.2. Boltzmann kinetic equation for a Bose gas at the 1D–3D crossover

Taking inspiration from the developments of GHD, Møller et al (2021) studied a three-dimensional (3D) Bose gas at the crossover to the 1D regime, and introduced a phenomenological description of the gas dynamics at that crossover. The idea is the following. The atoms lie in the 3D potential V(x) + V_⊥(r_⊥), where is the distance from the axis at y = z = 0. The transverse potential is harmonic, , so that the wavefunction of each atom ψ(x, r_⊥) can be expanded on the eigenstates ψ_n (r_⊥) of the transverse harmonic oscillator, see equation (103) below. Strictly in the 1D regime (i.e. when ℏω_⊥ is much larger than all other energy scales), the transverse ground state (n = 0) is the only state that is occupied. But if the transverse confinement is not strong enough, transversely excited states will also be occupied. Møller et al (2021) consider the first three transverse states (n = 0, 1, 2) (their degeneracy is neglected), and regard the resulting system as a three-component 1D Bose gas, with a Hamiltonian of the form

Here the excitation terms are of the form and correspond, for instance, to a two-body collision where two atoms in the transverse ground state get excited to the first excited state. The coupling constants g_a and g_ab (a, b = 0, 1, 2) are effective 1D coupling constants resulting from the 3D interaction. For a 3D scattering length much smaller than , they can be calculated from the shape of the transverse orbitals, see the discussion in subsection 3.1.2. Møller et al (2021) simply assume that they are equal: g_a = g_ab = g. Under that assumption, the multi-component Bose gas with V(x) = 0 is integrable (Caux et al 2009, Gaudin 2014, Klauser and Caux 2011, Klümper and Pâţu 2011, Pâţu and Klümper 2015, Robinson et al 2016, Robinson and Konik 2017, Yang 1967), however, the description of its thermodynamic limit is considerably more complicated than that given by the Lieb–Liniger model reviewed in section 1. In particular, while the macrostates of the Lieb–Liniger model in the thermodynamic limit are characterized by their distribution of rapidities ρ(θ), those of the multi-component Bose gas are characterized not only by ρ(θ), but also by infinitely many rapidity distributions for pseudo-spin bound states. Møller et al (2021) assume that the population of such pseudo-spin bound states is small and can be neglected, and keep only the distribution of rapidities ρ(θ), whose evolution under the Hamiltonian (90) then coincides with that of the Lieb–Liniger gas, and is therefore nothing but the GHD equation (78).

Finally, the effects of the excitation terms in the Hamiltonian (90) are introduced in the form of a Boltzmann collision integral:

The collision integral is constructed phenomenologically, by considering a simple model for two-body collisions. This phenomenological approach does not incorporate interactions with the other particles, contrary to what usually happens in exact Bethe ansatz calculations (where the interaction effects appear through the dressing of the various quantities that enter all the formulas). Within the framework of that very simple model, the probability of two colliding atoms with momenta p₁ = ℏk₁, p₂ = ℏk₂ initially in the same transverse state getting excited to a different transverse state is estimated as P_↕(k, q) ≃ 4c² kq/[k² q² + c²(k + q)²] where k = |k₁ − k₂| and . Then the authors assume that the atom momenta may be replaced by the rapidities, and arrive at an expression of the form , where ν_n is the probability that an atom is in the nth transversely excited state (assumed to be ≪1 for n = 1, 2), β_n is the number of atoms changing state in a collision (β₂ = 1 and β₁ = 2 in the model of Møller et al (2021)), and

with similar expressions for , and , where the indices ‘p’ and ‘h’ refer to ‘particle’ and ‘hole’, see (Møller et al 2021). Here are the momenta of the two atoms after getting excited to the transverse state in that model of two-body collisions.

The various approximations involved in the construction of the effective Boltzmann equations (91) and (92) are particularly meaningful in the ideal Bose gas regime. There, this phenomenological approach is expected to work well. Notice that in the ideal Bose gas regime, the GHD dynamics reduces to that of free Bosons, with the effective velocity appearing in equation (91) reducing to the bare velocity v^eff(θ) = θ/m. In that regime the aforementioned assumption that all intra- and inter-component coupling strengths are equal is actually no longer necessary for the gas to be integrable. Beyond the ideal Bose gas regime, the accuracy of the description of the 1D–3D crossover by equations (91) and (92) remains to be investigated.

Motivated by the recent experiment in the Newton cradle setup of Li et al (2020), which observed thermalization attributed to excitations of the transverse modes, Møller et al (2021) simulated the Newton cradle with equations (91) and (92) in the ideal Bose gas regime, and found that excitations to the transverse modes do indeed induce thermalization, see figure 11.

Zoom In Zoom Out Reset image size

2.4.3. Beyond the Euler scale: Navier–Stokes diffusive corrections

As mentioned in the introduction, Euler-scale hydrodynamic equations are but the zeroth order in a gradient expansion. More precisely, under the assumption of local relaxation, the expectation values of the currents are functions of the charges and of their derivatives. Schematically, , which is then expanded as . The zeroth order gives the Euler-scale hydrodynamic equations, while the next order gives the hydrodynamic equation at the diffusive scale. These hydrodynamic equations include a diffusive, entropy-producing (or Navier–Stokes) term.

The ‘Navier–Stokes’ GHD equation, with the diffusive term, was obtained in (De Nardis et al 2018). The result reads

where is a kernel (i.e. ) describing a Markov process of random momentum exchanges via two-body collisions (De Nardis et al 2019, Gopalakrishnan et al 2018). It is defined by the relation

with

For a detailed derivation of these equations, see (De Nardis et al 2019).

As mentioned above, interestingly, taking into account the diffusive correction changes the conclusion about thermalization in the quantum Newton cradle setup, see (Bastianello et al 2020a) and figure 12. With the inclusion of the Navier–Stokes term, the quantities (82) are no longer conserved, and it is found that the gas ultimately reaches thermal equilibrium, contrary to what had been observed previously in (Cao et al 2018, Caux et al 2019).

Zoom In Zoom Out Reset image size

For more studies of diffusion in the Lieb–Liniger gas, see e.g. (Medenjak et al 2020, Panfil and Pawełczyk 2019). Diffusion has also been studied very extensively in spin chains; for this topic we refer the reader to the review articles by Bulchandani et al (2021) and by De Nardis et al (2021) in this volume.

2.4.4. Quantum fluctuating hydrodynamics

As emphasized above, Euler-scale hydrodynamic equations are based on separation of scales (figure 6), which allows for a description of the gas as a continuum of independent fluid cells, each of which has relaxed to a stationary state.

In that hydrodynamic picture, a small perturbation of the system at spacetime position (x, t) generates sound waves that propagate through the gas, so that the perturbation may be observable at a different position (x′, t′) (Kadanoff and Martin 1963). Such dynamical correlations have been studied in the framework of GHD in (Doyon 2018, Møller et al 2020); see also the review by De Nardis et al (2021) in this volume.

However, at equal time, the fluid cells at different positions x and x′ are independent. Thus, all equal-time connected correlations vanish at the Euler scale: for any local observables ,

This is somewhat puzzling because, in a quantum system, equal-time correlations are typically non-zero: for instance, in the ground state of the Lieb–Liniger gas they are non-zero, see e.g. the reviews (Cazalilla 2004, Cazalilla et al 2011). However, this is not in contradiction with (95): the reason is simply that such non-zero equal-time correlations are an effect occurring beyond the Euler scale.

A good illustration of this phenomenon is provided by the standard fluid described by the Euler equations (3). At zero temperature, the third Euler equation (3) is automatically satisfied as a consequence of the first two (see the discussion in subsection 2.1), so we are left with

For simplicity, we now consider the ground state of the spatially homogeneous system (V(x) = 0), with n(x, t) = n and u(x, t) = 0. Linearizing the system (96) around that solution leads to the equation of propagation of sound waves (δn(x, t), δu(x, t))

or equivalently

Here v is the sound velocity in the gas, and K is a dimensionless parameter (called the Luttinger parameter, see below) normalized such that K → 1 in the hard-core limit. We see from equation (98) that there are right- and left-moving sound waves, corresponding to specific linear combinations of δn and δu parameterized by K, traveling at velocity ±v.

The sound waves can be used as the basic ingredient in a quantized theory of the fluid described by the Euler equations (96). The basic idea is to view δn(x) and δu(y) as operators in a quantum theory, and to impose the canonical commutation relations of a one-component fluid (Landau 1941),

and . Then, one should find a Hamiltonian H such that the Heisenberg equations and coincide with the equations of motion (97). The simplest choice is:

which is nothing but the Hamiltonian of a Luttinger liquid, see e.g. (Cazalilla 2004, Giamarchi 2003, Tsvelik 2007) for introductions. Thus, we see that quantizing the sound waves of a standard Euler fluid at zero temperature (96) directly leads to a Luttinger liquid. With that observation, one can get the leading behavior of equal-time correlations beyond the Euler scale. For instance, the equal-time density–density correlation in the ground state of the Hamiltonian (100) is (Cazalilla 2004, Giamarchi 2003, Tsvelik 2007):

Such correlation functions can then be propagated in time using the sound wave equation (98). For instance, propagating equation (101) in time leads to a combination of two terms, one coming form the right-moving sound wave, the other from the left-moving one:

Other time-dependent correlations can be obtained in a similar way within the framework of Luttinger liquid theory, see e.g. (Cazalilla 2004, Giamarchi 2003, Tsvelik 2007). These dynamical correlation functions obtained from the propagation of linear sound waves are valid on time scales that are not too long, before non-linear effects kick in. For a review on such non-linear effects, see (Imambekov et al 2012).

It is natural to ask whether such a program can be implemented, replacing the standard Euler equations at zero temperature (96) by those of GHD at zero temperature (88). This was achieved, to some extent, by Ruggiero et al (2020). In this work, the equations for linear sound waves around zero-entropy GHD are derived, and quantized by imposing canonical commutation relations. This procedure results in a time-dependent, spatially inhomogeneous, multi-component Luttinger liquid.

In the ground state of the gas, the Hamiltonian obtained by Ruggiero et al (2020) is that of an inhomogeneous Luttinger liquid, see e.g. (Bastianello et al 2020c, Brun and Dubail 2018, Cazalilla 2004, Dubail et al 2017). At later times it follows the evolution of the system under the zero-entropy GHD equation (88), encoding the propagation of the quantum fluctuations from time t = 0 to times t > 0, see figure 13. For more details we refer the reader to (Collura et al 2020, Ruggiero et al 2019, 2020), and to the review by Alba et al (2021) in this volume.

Zoom In Zoom Out Reset image size

We also mention that there is, in principle, another approach to quantum fluctuations in GHD (Fagotti 2017, 2020), which consists of viewing the GHD equation (78) as the zeroth order in the evolution equation for the Wigner quasiprobability distribution (Bettelheim et al 2006, Bettelheim and Glazman 2012, Bettelheim and Wiegmann 2011, Dean et al 2018, 2019, Doyon et al 2017, Fagotti 2020, Moyal 1949, Protopopov et al 2013, Ruggiero et al 2019). The evolution equation of the Wigner function may be expanded in powers of ℏ (Moyal 1949), with higher-order terms that give the corrections to the Euler-scale GHD equation. In this approach, it is not only the quantum fluctuations at time zero that are propagated in time as above. Corrections also appear dynamically because of the modified evolution equation. So far, this approach has been limited to the hard-core limit g → +∞, or more generally to spin chains that map to non-interacting fermions. To our knowledge, the extension of this approach to the interacting case is an open problem. We refer the reader to the aforementioned review by Alba et al (2021) in this volume for a thorough discussion of this topic.

3.1.1. One-dimensional geometries

Confining potentials for cold atoms can be realized by different means (see for instance the book by Pethick and Smith (2008)). One can use laser fields: the interaction between the induced atomic dipole and the laser field results in a force acting on the center-of-mass motion of the atom, which is conservative if the laser frequency is sufficiently far from any atomic resonance and which is proportional to the laser intensity gradient. One can also rely on a non-vanishing magnetic moment of the atoms to realize potentials using spatially varying magnetic fields. For large enough magnetic fields, adiabatic following of the orientation of the magnetic dipole of the atoms ensures the realization of conservative potentials.

One-dimensional gases are realized in cold atom setups when the atoms are confined in guides with a transverse confinement large enough that the energy gap between the transverse ground state and the first excited state is much larger than the typical energy per atom. Then the transverse degrees of freedom get frozen and the resulting dynamics are effectively 1D. The strong transverse confinement required to reach the 1D regime can be realized using laser fields with large intensity gradients resulting from interference: a two-dimensional (2D) optical lattice formed by lasers far-detuned from the atomic transitions realizes a 2D array of 1D tubes for the atoms where the 1D regime can be accessed (Fabbri et al 2011, Haller et al 2009, Kinoshita et al 2005), see figure 14. Such setups present several advantages. First, they offer the possibility of realizing very strong transverse confinements, since the light intensity is modulated at distances as small as the laser wavelength. Second, cold atoms are loaded into such a lattice adiabatically from a 3D Bose–Einstein condensate (BEC), which permits one to obtain very cold 1D gases. Finally, in this setup the magnetic field is a free parameter that can be tuned to approach Feshbach resonances and thus adjust the interaction strength between atoms (Haller et al 2009). On the other hand, optical lattice setups also have some limitations. The first is that all measurements are ensemble measurements: the value of the measured observable is averaged over several 1D tubes, with parameters such as the number of atoms or the interaction strength varying from tube to tube, which complicates the data analysis. This averaging also prevents the investigation of fluctuations and correlations. Second, the 2D optical lattice has a finite length in the longitudinal direction—the direction along the tubes—since it is produced using lasers with relatively small waists. This imposes restrictions on the spatial extension of the 1D gases: the clouds cannot be too long. However, this limitation does not prevent cloud expansions that are long enough to allow access to the rapidity distribution, see section 3.3.

Zoom In Zoom Out Reset image size

Another type of experiment uses magnetic trapping of atoms to reach the 1D regime. Strong magnetic gradients are obtained in the vicinity of microwires running an electrical current in the so-called atom chip setup first developed in Munich (Reichel et al 1999) and in Vienna (Folman et al 2000) (see (Reichel and Vuletic 2011) for a review). Close enough to the microstructures, the 1D regime can be reached (Armijo et al 2011, Jacqmin et al 2011, van Amerongen et al 2008). The idea of magnetic guiding is simple. Consider a current-carrying wire immersed in a homogeneous magnetic field which is perpendicular to it. The total magnetic field vanishes on a line parallel to the wire, and atoms, when polarized in a low-field seeker magnetic state, will be guided along this line. The transverse homogeneous magnetic field can be replaced by the magnetic field produced by two wires placed on each side of the central wire, so that one has the configuration depicted in figure 14. Strong transverse confinement is obtained close to the wires and, in some experiments, atoms are guided at a distance of as small as 15 μm from the wires (Jacqmin et al 2011). The advantage of atom chip setups is that a single 1D gas is realized and studied. This allows direct comparison with theoretical predictions. This also permits the study of fluctuations, and of their correlations. The guiding of atoms can be realized on very long distances, as the transverse confinement is invariant along the whole microwire. However, atom chips also have some drawbacks. In particular, the 1D gases are prepared and cooled down directly in the 1D geometry, and it turns out that temperatures obtained in these setups are higher than those obtained using 2D optical lattices. The physical phenomena limiting the temperatures that can be reached are not completely elucidated, although the role of atom losses has recently been singled out (Bouchoule and Schemmer 2020, Rauer et al 2016, Schemmer and Bouchoule 2018).

3.1.2. Effective 1D interaction parameter

In experiments, although the transverse confinement is large enough to freeze the transverse degrees of freedom, it is typically small enough that the width of the transverse ground state wavefunction is much larger than the range of the 3D interaction potential. In that case, one can model the 3D interaction by a contact potential characterized by the scattering length a_3D. The relation between the 3D scattering properties and the 1D scattering properties was worked out by Olshanii (1998), and we briefly review the main steps of his argument. A harmonic transverse confinement of frequency ω_⊥ is assumed, and the scattering state of two atoms is explicitly constructed. Because of the harmonic nature of the transverse confinement, the center-of-mass motion decouples from that of the relative coordinate. The wavefunction of the relative coordinate obeys the Schrödinger equation for a particle of mass m/2 in a harmonic transverse confinement of frequency ω_⊥ and with the 3D contact potential at the origin. It can be expanded as a sum of transverse eigenstates, with x-dependent amplitudes. One considers a state whose energy is small enough that, apart from the transverse ground state, the amplitudes are exponentially decreasing functions of x. Then the 3D wavefunction reads,

where the sum runs over the transversely excited states of vanishing angular momentum whose wavefunctions are ψ_n (r_⊥), where r_⊥ is the transverse distance to the origin. The parameters {a_n } and δ are determined in the following way. First, one imposes that ∂ψ/∂x is regular everywhere in the plane x = 0 except at the 3D origin. Then, one imposes that the wavefunction fulfills the 3D boundary condition: for a zero-range 3D interacting potential of scattering length a_3D, the wavefunction diverges as 1/r − 1/a_3D at short distance. This fixes the value of tan(δ), which is found to scale as 1/k at small k (Olshanii 1998). For a pure 1D problem, with a contact potential g_1D δ(x), one would have tan(δ) = −mg_1D/(2ℏ² k). Comparing with the small k limit of the 3D problem, one finds that the effective 1D coupling constant is (Olshanii 1998)

where is the width of the transverse ground state and .

In most experimental realizations, a_3D ≪ a_⊥, so that Olshanii’s relation (104) reduces to g_1D = 2ℏω_⊥ a_3D. This expression can be recovered by simple reasoning, valid for a weakly interacting gas in the quasicondensate regime. Since we assume a_3D ≪ a_⊥, interactions have a 3D nature and their effect is obtained simply by averaging over the transverse density profile. More precisely, the interaction energy per unit length in the longitudinal direction is , where g_3D = 4πℏ² a_3D/m is the 3D interaction strength. On the other hand, considering the 1D problem one obtains . Equating both expressions, and using the Gaussian shape of ψ₀, one recovers g_1D = 2ℏω_⊥ a_3D.

Equation (104) shows a resonance behavior when reaches one, which can be interpreted as a Feshbach resonance involving a bound state related to a transversely excited state (Bergeman et al 2003). Experimentally, approaching this resonance would require very large transverse confinements for standard 3D scattering lengths. However, close to a 3D scattering resonance, one can reach the regime where a_3D becomes close to or larger than a_⊥. This was used by Haller et al (2009) to reach the 1D hard-core regime, and even realize a metastable state with attractive 1D interactions (i.e. g_1D < 0).

3.1.3. Validity of the assumption of separation of scales

The realization of the 1D Bose gas with contact repulsive interactions in experimental cold atom setups has been established by comparisons of experimental data with exact predictions from the Lieb–Liniger model. The theoretical predictions have used either the machinery of integrable systems reviewed in section 1, or numerically exact Monte Carlo calculations, valid for gases at thermal equilibrium. Here we present a few results which, in our view, constitute landmarks in this field. All the results reviewed here are about gases of many atoms (N ∼ 10–10⁴) and are in good agreement with the LDA analysis.

3.2.1. Measurement of the zero-distance correlation function

The zero-distance two-body correlation function can be obtained exactly using the Hellmann–Feynman theorem, see equation (75). For the ground state it gives

where e₀ is the ground state energy density and the derivative is taken at constant linear density n. This quantity can be computed numerically (Lieb 1963). As discussed in subsection 1.8, the zero-distance correlation function characterizes the crossover between the quasicondensate regime (where g⁽²⁾(0) ≃ 1) and the hard-core regime (g⁽²⁾(0) ≃ 0).

Experimentally, the zero-distance correlation function can be measured by probing two-body losses induced by photoassociation (Kinoshita et al 2005), see figure 15. The idea is to shine a laser beam on the cloud that is sufficiently far from the atomic resonance to leave isolated atoms at rest, but whose frequency is chosen to induce, for pairs of atoms that are very close to each other, a transition toward an excited diatomic molecule. The distance required to perform a photoassociation is much smaller than the typical length scale of variation of , such that the rate of production of excited molecules is simply proportional to n² g⁽²⁾(0). When a photoassociated molecule decays, it typically produces atoms whose kinetic energy is much larger than the trap depth; these atoms thus leave the atomic cloud. In the end, one observes a decrease of the total atom number N, and the loss rate dN/dt is proportional to n² g⁽²⁾(0).

Zoom In Zoom Out Reset image size

The experiment of Kinoshita et al (2005) is realized in a 2D optical lattice, so that a 2D array of 1D inhomogeneous tubes is populated. More precisely, the atom density varies within each tube because of the harmonic longitudinal confinement V(x), and the atom number N per tube varies among the tubes. The data analysis uses the LDA. The local loss rate is proportional to the local value of n² g⁽²⁾(0), where g⁽²⁾(0) depends on n via its dependence on the interaction parameter γ = mg/(ℏ² n): . One assumes that is linear in log(γ), which is a good approximation for the data range explored (see figure 15). Then the loss rate is expected to be proportional to , where is the mean linear density seen by atoms and γ_eff is such that log(γ_eff) is the spatial average of log(γ), weighted by n². To evaluate γ_eff and , the distribution of atoms among the tubes is assumed to be that corresponding to the initial shape of the 3D BEC, and, within a 1D tube, the Lieb–Liniger equation of state is used. Finally, the value of deduced from experimental data is shown in figure 15. It is in remarkable agreement with the prediction from the Lieb–Liniger model. The data show the crossover between the quasicondensate regime, where g⁽²⁾(0) ≃ 1, as in a true BEC, and the hard-core regime, where g⁽²⁾(0) ≪ 1, as for a Fermi gas.

Shortly before that measurement of g⁽²⁾(0), the observation of 1D Bose gases in the hard-core regime had been achieved by Kinoshita et al (2004), who measured the total energy of an array of 1D gases by 1D expansion, as well as the length of the trapped atom clouds. The approach to the hard-core regime was also signaled by a measured reduction of the three-body loss rate (Tolra et al 2004) (see section 5 for a discussion of three-body losses).

3.2.2. Analysis of density profiles: Yang–Yang thermodynamics

Remarkably, the results for g⁽²⁾(0) briefly reviewed in the previous subsection are in good agreement with the theoretical prediction for the ground state of the 1D Bose gas. However, in many experimental situations, the gas is not in its ground state. As reviewed in subsection 1.6, exact theory results from Yang and Yang (1969) are available for thermodynamic quantities at thermal equilibrium. In several experimental works, the data were found to be in very good agreement with predictions from Yang–Yang thermodynamics.

The first comparison between experiments and the finite-temperature thermodynamics of Yang and Yang was done by van Amerongen et al (2008), who analyzed the density profiles of trapped 1D gases. This work used an atom chip setup: atoms are confined in magnetic traps realized by microwires deposited on a chip. In contrast with optical trapping, a single 1D cloud is observed, which allows a more direct comparison with theoretical predictions, since no averaging over tubes is required. This permits investigation of the density profile of the 1D gas, confined in a longitudinal potential V(x). As discussed in subsection 3.1.3, the longitudinal trapping is weak enough that the LDA is valid. For a gas at thermal equilibrium, the LDA means that the gas at position x is described by a gas at temperature T and chemical potential μ(x) = μ₀ − V(x). Thus, the density at position x is

where n_YY(T, μ) is the linear density of the homogeneous Lieb–Liniger gas at temperature T and chemical potential μ, referred to as the Yang–Yang equation of state. Thus, if V(x) is known, the experimental density profile can be compared to the theoretical curve, with two adjustable parameters: the temperature T and the chemical potential μ₀. In the experiment of van Amerongen et al (2008), however, the population of transversely excited states was not negligible. This was taken into account in the theoretical model by treating each transversely excited state as an ideal 1D Bose gas, at thermal equilibrium with the 1D gas in the transverse ground state. This treatment results in a modified Yang–Yang equation of state n(T, μ). Figure 16 reproduces the experimental density profiles of van Amerongen et al (2008) compared with the theoretical curves of the modified Yang–Yang equation of state, with μ₀ and T as fitting parameters. The theoretical curves reproduce the measured profiles very well.

Zoom In Zoom Out Reset image size

Since these pioneering results, the analysis of density profiles of 1D gases using the Yang–Yang equation of state has been realized successfully in other experiments (Armijo et al 2011, Vogler et al 2013). In (Vogler et al 2013) the experimental setup uses an optical lattice, and density profiles are acquired using an electron beam propagating perpendicularly to the longitudinal direction x. The resolution attained using electron-beam imaging is sufficiently good to resolve individual 1D tubes. However, the data analysis is complicated by the fact that column-integrated density profiles are acquired, the integration being done over the direction of propagation of the electron beam. Thus, the raw data mix information about different 1D tubes that have different density profiles. The density profile of individual 1D tubes can, however, be extracted deconvoluting for the effect of the column integration, if one assumes invariance by rotation of the tubes’ distribution. The deconvolution technique, that does not require an a priori knowledge of the tube distribution, uses the Abel transformation. The resulting density profiles of individual 1D gases, shown in figure 16, fit remarkably well with the Yang–Yang theory.

3.2.3. Density fluctuations

A more stringent test of Yang–Yang thermodynamics can be made if one investigates not the density profile but its fluctuations. For this purpose, one can investigate the atom-number fluctuations in a pixel whose size Δ is much smaller than the length of the cloud, and much larger than the correlation length of the gas. We note the atom-number fluctuations in a pixel δN = N − ⟨N⟩, where N is the atom number in the pixel.

Then, if one assumes thermal equilibrium, the gas contained in the pixel can be described by a Gibbs ensemble, the rest of the cloud acting as a reservoir of energy and particles. Then atom-number fluctuations fulfill, for any integer K,

where n(μ, T) is the equation of state of the gas. Thus, the measured density fluctuations can be compared to expectation values derived from the Yang–Yang equation of state. Notice that, in contrast with the experiments presented above, no precise knowledge of the longitudinal potential V(x) is required. The longitudinal potential is irrelevant in the data analysis: it is simply useful to sample different densities at a given temperature.

In order to measure density fluctuations, ensemble measurements, such as those realized in setups using 2D optical lattices, are prohibited. Thus, such measurements have been realized only in an atom chip setup. To extract density fluctuations, one records an ensemble of density profiles, all taken with identical experimental parameters, from which a statistical analysis is performed. Measurements of density fluctuations in 1D gases were first performed in (Esteve et al 2006). Comparison with Yang–Yang predictions were performed in (Armijo et al 2011, 2010, Jacqmin et al 2011). In (Armijo et al 2011), the dimensional crossover from 3D to 1D was investigated, and it was shown that it is well accounted for by the aforementioned modified Yang–Yang equation of state, which includes the effect of transversely excited states treated as ideal Bose gases. Figure 17 displays a selection of results from the two latter references.

Zoom In Zoom Out Reset image size

3.2.4. Techniques for measuring the momentum distribution

While all the above results concern observables that involve only the real-space density, information is also contained in other observables, in particular the momentum distribution. The momentum distribution has been measured in several experiments using different techniques, which we briefly review now.

The first technique that was used was the so-called Bragg technique. The idea is to shine two laser fields with wavevector difference q and frequency difference ω onto the atomic cloud. The action on the cloud is described by the addition of a new term to the Hamiltonian, which reads

where each term corresponds to the absorption of a photon in one laser beam and stimulated emission in the other, a process called a Bragg process. For counterpropagating lasers, q is usually very large compared to the typical momentum width of the cloud, and ℏ² q²/m is typically much larger than the energy per atom. In this case, called the Doppler regime, the atom promoted to the momentum state k + q by the Bragg process can be considered as effectively removed from the system, and the final energy is about ℏ²(k + q)²/(2m) ≃ E_rec + ℏ² kq/m, where E_rec = ℏ² q²/m is the recoil energy, and the second term is the Doppler shift. Then the Fermi golden rule leads to a rate of Bragg transfer Γ which obeys

where conservation of energy imposes k = m(ℏω − E_rec)/(ℏ² q). Thus, measuring the loss rate of the atomic cloud as a function of the Bragg detuning ω gives access to the momentum distribution of the atoms. The energy deposited in the system per unit time is nothing but Γℏω, such that one can equivalently deduce the momentum distribution from the measurement of the energy increase rate versus ω. Note that, by choosing ω close to N² E_rec, where N is an integer, one can induce Nth-order Bragg processes: the momentum transfer is then Nq, which allows deeper access inside the Doppler regime. This technique was first used in (Stenger et al 1999) in a 3D BEC, where it confirmed the presence of long-range order in the cloud. Applied to a very elongated BEC, it was used to demonstrate the presence of longitudinal thermally excited phase fluctuations (Richard et al 2003). More recently, it was implemented in 1D gases realized in a 2D optical lattice (Fabbri et al 2011). Note that, going beyond the Doppler limit of large momentum transfer, the Bragg techniques allow the measurement of the dynamic structure factor. Such measurements have been compared to results based on the Lieb–Liniger model in (Fabbri et al 2015) and (Meinert et al 2015).

Bragg momentum spectroscopy suffers from low signal amplitude, since it is based on a perturbative analysis, and, most importantly, it does not allow recording of the whole momentum distribution in a single measurement. This impedes the measurement of correlations in momentum space. We now discuss another technique that enables the recording of the whole momentum distribution of 1D gases in a single measurement. The key point is to be able to switch off the interactions almost instantaneously with respect to the longitudinal motion. This can be achieved, close to Feshbach resonances, by manipulating the interaction strength with a magnetic field (Stewart et al 2010). However, in the special case of 1D gases, this switch-off can be achieved very easily by removing the transverse confinement: after the switch-off of the transverse potential, the transverse wavefunction expands in a typical time of the order of 1/ω_⊥, much shorter than typical times associated with the longitudinal motion, and this amounts to an almost instantaneous vanishing of the effective 1D interaction strength with respect to the longitudinal motion.

Once interactions have been effectively switched off, one is left with the task of measuring the momentum distribution of an ideal gas. This could be done by a simple ballistic expansion, following the sudden switch-off of the longitudinal potential: after an expansion time long enough that the final cloud size is much larger than its initial size, the density distribution becomes homothetic to the momentum distribution. This technique, usually referred to as the time-of-flight technique, is used for instance in (Fabbri et al 2011) and (Malvania et al 2021, Wilson et al 2020). However, in experiments using very long initial clouds and/or gases lying deep in the quasicondensate regime, such as atom chip experiments, this technique typically requires unrealistic expansion times. To overcome this difficulty, one can use the so-called focusing method. This method, first implemented in (Shvarchuck et al 2002), consists of applying a short pulse of a strong longitudinal harmonic potential. Atoms do not have time to move during this pulse, but they acquire a momentum kick δp = −Ax proportional to their distance x from the center. The cloud then undergoes free evolution. After a time equal to the focusing time t_f = m/A, the density distribution is homothetic to the initial momentum distribution. This technique effectively erases the information on the initial cloud longitudinal spread.

3.2.5. Results in momentum space

All measurements reviewed in subsections 3.2.1–3.2.3 probe thermodynamic quantities. On the theory side, those are easily accessible numerically using Yang–Yang thermodynamics for thermal equilibrium states. In contrast, the momentum distribution is not a thermodynamic quantity.

The momentum distribution of 1D gases lying in the quasicondensate regime was measured using Bragg spectroscopy by Fabbri et al (2011). The authors show that the measured momentum distribution is close to a Lorentzian. A Lorentzian shape of full-width at half-maximum mk_B T/(ℏ² n) is expected for homogeneous gases in this regime, for wavevectors lying in the phononic regime (Jacqmin et al 2012). Within the LDA, the total momentum distribution n(p) is obtained by summing the contributions of each small fluid cell. For a harmonically confined quasicondensate, one finds a momentum distribution that stays very close to a Lorentzian (Jacqmin et al 2012).

The focusing method was used for 1D gases trapped on an atom chip in (Davis et al 2012). The measured momentum distribution n(p) was compared to improved classical field calculations, expected to be valid in the quasicondensate and ideal Bose gas regimes. The transversely excited states cannot be neglected in those experiments and they were taken into account assuming they behave as ideal Bose gases. In fact, in this paper, the authors propose a thermometry method. More precisely, they evaluate the total kinetic energy E_K = ∫n(p)p²/(2m)dp from the measured momentum distribution n(p). E_K is a thermodynamic quantity that can be calculated with Yang–Yang thermodynamics, and that can be fitted to extract the temperature. The fact that E_K is a thermodynamic quantity follows from the following argument. E_K can be computed using the LDA as the integral of the kinetic energy density e_K(x) over the cloud. The latter is e_K(x) = e(n(x), T) − e_int(n(x), T) where e(n, T) is the energy density of the Lieb–Liniger gas at density n and temperature T, and e_int(n, T) is its interaction energy. Both are thermodynamic quantities: the interaction energy can be obtained from the Hellmann–Feynman theorem, see equation (75). For data analysis, one adds the contribution of the transversely excited states, considered as ideal Bose gases.

The first comparison between the measured momentum distribution of a 1D Bose gas and exact calculations was done in (Jacqmin et al 2012). The first-order correlation function at thermal equilibrium was computed with a Monte Carlo method, which gave very good agreement with the measured momentum distribution. The measured in situ density fluctuations (subsection 3.2.3) also fitted very well with the temperature fitted from the Monte Carlo simulations. The momentum distribution n(p) was measured across the smooth transition between the ideal Bose gas regime and the quasicondensate regime. No striking modification of n(p) was observed across the transition; in particular, it conserves a Lorentzian-like shape. (We recall that the crossover occurs for a degenerate gas whose momentum distribution, in the ideal gas regime, has a Lorentzian central shape.)

Using the focusing technique, the full momentum distribution of the sample can be recorded in a single measurement. Thus, it is possible not only to extract the mean momentum distribution n(p), but also its fluctuations δn_p . In (Fang et al 2016), the correlations ⟨δn_p δn_p′⟩ are deduced from a statistical analysis of hundreds of images taken in the same experimental conditions. The results are in very good agreement with numerically exact results obtained, for a gas at thermal equilibrium, using a quantum Monte Carlo algorithm, as seen in figure 18. The calculation uses a discretized model and a worm algorithm. The crossover between the ideal Bose gas and the quasicondensate regime has a clear signature in momentum space correlations. The ideal Bose gas regime is characterized by a bunching phenomenon for equal momenta. In the quasicondensate regime, on top of the bunching seen on the diagonal, anti-correlations appear for different momenta. Those anti-correlations ensure small density fluctuations, the latter being inhibited in the quasicondensate regime because of the large interaction energy they would require. The presence of those anti-correlations in the quasicondensate regime was predicted using Bogoliubov theory in (Bouchoule et al 2012).

Zoom In Zoom Out Reset image size

3.2.6. Other higher-order correlation functions

The two-body correlations in momentum space presented in figure 18 are related to the four-point correlation function of the atomic creation/annihilation operator. The four-point correlation function can also be probed by investigating density fluctuations resulting from a short free evolution time. The idea is to switch off the interactions between atoms abruptly, for instance by removing the transverse confinement (subsection 3.2.4), and then to remove the longitudinal confinement and let the gas evolve freely for a short time t_free. In contrast with the time-of-flight method used to measure momentum distribution (subsection 3.2.4), here the free evolution time is assumed to be short, such that the mean density profile of the cloud barely changes. Thus, no information is gained from the mean profile. Instead, one investigates the density fluctuations, also called density ripples. They are characterized by their spectral density ⟨|ρ(q)|²⟩. For wavelengths much shorter than the size of the mean density profile, the latter is related to the four-point correlation function of the atomic field,

where expectation values are taken in the 1D gas before the expansion. For very large q, since there is no long-range order in the 1D gas, the atomic fields at positions separated by qt_free are uncorrelated and the above expression vanishes. On the other hand, for very small q one recovers the spectral density of density fluctuations that were present in the initial gas. This measurement turns out to be particularly relevant in the quasicondensate regime. There, the initial density fluctuations are negligible and ⟨|ρ(q)|²⟩ results from the transformation of phase fluctuations into density fluctuations during the free evolution. Density ripples were first investigated in (Dettmer et al 2001) in a quasi-1D setup. It was then used for thermometry in (Manz et al 2010) in a 1D gas realized in an atom chip setup. For q small enough that ℏqt_free/m is very small compared to the correlation length of the phase fluctuations, one finds from equation (110) that the power spectrum of density fluctuations after free evolution during t_free is simply proportional to the power spectrum of the initial phase fluctuations at the same wavevector. Thus, the analysis of density ripples allows one to access each Bogoliubov mode individually. This feature was used in (Schemmer et al 2018) to monitor the dynamics produced by an interaction quench. Finally, let us stress that the measurement of density ripples is clearly very different from the measurement of the momentum distribution, as the momentum distribution n(p) mixes all Bogoliubov modes.

3.2.7. Validity of the thermodynamic equilibrium

As explained in section 1.2, rapidities are the asymptotic momenta of the atoms after a 1D expansion. This characterization can be viewed as the definition of the rapidities. It also shows that the rapidity distribution is an observable. Owing to the important role of the rapidity distribution, the experimental ability to measure it is a key development for the study of 1D gases.

The first experimental measurement of the rapidity distribution was done by Wilson et al (2020), for gases lying quite deep inside the hard-core regime. In this experiment the trapping potential is the sum of a 2D array of 1D tubes realized by a 2D optical lattice, and a slowly varying trapping potential, which provides a longitudinal confinement along the tubes. Removing the slowly varying potential, a 1D expansion is performed within each tube. After a sufficiently long expansion time, the momentum distribution converges toward the rapidity distribution. The momentum distribution is then measured using the time-of-flight technique, performed by suddenly turning off all confining potentials (including both the 1D longitudinal confinement and the 2D array of tubes). Interactions are effectively almost instantaneously turned off by the rapid transverse expansion of each tube. Thus, the cloud performs a ballistic expansion such that, at long time, the density profile reflects the momentum distribution.

The momentum distribution measured for different 1D expansion times is found to be in very good agreement with theoretical predictions, see figure 19 from (Wilson et al 2020). The calculation of the evolution during the 1D expansion assumes hard-core bosons and is done for a gas initially in the ground state. For long enough expansion times, the momentum distribution converges toward the rapidity distribution.

Zoom In Zoom Out Reset image size

In their famous experiment on non-equilibrium dynamics in a 1D Bose gas, Kinoshita et al (2006) use an array of 1D tubes, in a blue-detuned 2D optical lattice setup, with a longitudinal confinement, approximately harmonic, realized by an additional smooth dipole trap. The interaction parameter γ (see equation (7)) at the center of the trap, averaged over the collection of 1D gases, ranges from 0.6 to 4, depending on the data set.

An out-of-equilibrium initial situation is realized by applying two Bragg pulses on the cloud, such that the zero-momentum state is mostly transferred to a superposition of the two momentum states p = ±2ℏk. (Because of the interactions between atoms, that momentum transfer is not perfect; for a theoretical study of the state generated by the Bragg pulse, see (van den Berg et al 2016).) The cloud then evolves freely inside the trap up to time t. At time t, the longitudinal potential is turned off and atoms undergo a 1D expansion during an expansion time t_exp., large enough so that the final size of the atomic cloud is much larger than its initial size. After the 1D expansion, an image is recorded. The image performs a column integration in a direction perpendicular to the longitudinal direction. Figure 20 shows such images, for different evolution times t that span an oscillation period of the longitudinal potential: τ = (2π)/ω_∥, where ω_∥ is the frequency of the longitudinal potential. At t = 0, one observes two well-separated clouds, corresponding to the two components of different momenta created by the Bragg pulses. Then, one sees a dynamics resembling that which would be found if the atoms were non-interacting. In particular, at time τ/2, one recovers a situation close to the initial one: roughly speaking, the ‘cloud’ of initial momentum 2ℏk has performed half an oscillation in the harmonic trap and its momentum is now −2ℏk.

Zoom In Zoom Out Reset image size

After an evolution of about 10τ the images show small variations on an oscillation period. This is consistent with the dephasing effect due to the spreading of ω_∥ among the 1D tubes. A slow time evolution of the longitudinal profile is observed, see figure 20. This slow evolution is attributed to atom losses and to a small heating rate. Importantly, the measured longitudinal distribution does not evolve toward a thermal equilibrium distribution, at least not on the time scale probed in the experiment.

Theoretical modeling of this experiment is easy only in two asymptotic regimes: the ideal Bose gas regime, and the hard-core regime. In both cases, the dynamics is that of a non-interacting gas, see subsection 1.8. In those limits, one expects to observe undamped oscillations going on for ever, for a single 1D tube and a purely harmonic trap. The data of Kinoshita et al (2006) are compatible with this interpretation: at short times, the evolution is close to that of an ideal gas, and at longer times the observed damping can be attributed to dephasing between 1D gases. At very long times the evolution can be attributed to atom losses. The observed behavior is thus qualitatively similar to that expected for hard-core bosons, and this is because γ is sufficiently large so the 1D gases are quite well inside the hard-core regime. However, a quantitatively accurate modeling of those results, properly taking into account the finite value of γ, was completely out of reach at the time when (Kinoshita et al 2006) was published. In the decade that followed the experiment, no theory was capable of simulating it, taking into account the finite value of the interaction strength.

GHD is the first and, until now, only theory capable of obtaining quantitative predictions for such an experiment, valid for any initial situation (see section 2.2). This illustrates the power of that theory. The physical picture, detailed in section 2.2, is summarized below. The time evolution of the rapidity density ρ(x, θ) develops sharp structures due to the trap anharmonicity (Caux et al 2019). At some point those structures will be so narrow that the large-scale approximation of Euler-scale GHD will fail. One then expects the fine structures to disappear and the rapidity distribution to tend to a rapidity distribution that is a stationary solution of the GHD equations and that is non-thermal (Cao et al 2018, Caux et al 2019). On even longer time scales, under the combined effects of the integrability breaking produced by the longitudinal potential and diffusive terms, which are beyond Euler-scale GHD, the system will eventually drift toward a thermal equilibrium state (Bastianello et al 2020a). Relaxation toward a thermal equilibrium state has also been observed in classical field numerical simulations of the Newton cradle setup (Thomas et al 2021). The classical field, which is expected to describe weakly interacting gases with large mode population (see section 1.8.5), is not restricted to the description of long-wavelength behavior: it is beyond Euler-scale GHD, and this explains why it can lead to thermalization in the presence of an external potential.

In the quantum Newton cradle experiment, effects that are beyond pure 1D physics may also play a role. The population of transversely excited states is negligible since ℏω_⊥ (i.e. the energy gap between the transverse ground state and the first excited state) greatly exceeds the typical longitudinal energy per atom. In particular, thanks to the use of blue-detuned lasers for the realization of the 2D lattice, the potential depth, limited by longitudinal trapping, is smaller than ℏω_⊥, which ensures that the longitudinal energy stays smaller than ℏω_⊥. However, even when they are not populated, transversely excited states can contribute as virtual states in three-body processes (Mazets et al 2008). This phenomenon introduces an effective three-body interaction that breaks the integrability of the Lieb–Liniger model, and it might contribute significantly to the relaxation toward a thermal equilibrium.

Similar quantum Newton cradle experiments have been reproduced in (Li et al 2020, Schemmer et al 2019, Tang et al 2018). Li et al (2020) studied the effect of integrability breaking due to the presence of atoms in transversely excited states, while Tang et al (2018) investigated the effects of dipolar interactions.

A first experimental demonstration of the validity and relevance of GHD was carried out by Schemmer et al (2019). In this experiment, a 1D gas is realized on an atom chip. The initial cloud is at equilibrium in a regime close to the quasicondensate regime. Dynamics is initiated by a quench of the longitudinal potential, and the time evolution of the density profiles is recorded. For all situations considered, the results are in very good agreement with the predictions of the GHD theory.

The data are also compared to predictions of standard hydrodynamics, given by the Euler equations (3) of the introduction (see also subsection 2.1). In contrast with GHD, this standard hydrodynamic approach assumes that the gas is locally at thermal equilibrium. The numerical solution of the Euler equations involves the numerically tabulated pressure which is obtained from the Yang–Yang equation at thermal equilibrium, see subsection 1.6. A clear failure of standard hydrodynamics is found when the 1D gas is prepared at equilibrium in a double-well potential, and the double well is suddenly switched off and the gas is allowed to expand freely in 1D. In that case, the predictions of standard hydrodynamics differ strongly from those of GHD. Figure 21 shows the experimental data, together with GHD calculations and standard hydrodynamics calculations. The experiment clearly discriminates between both theories, and the data are found to be in agreement with GHD but not with standard hydrodynamics.

Zoom In Zoom Out Reset image size

The origin of the failure of standard hydrodynamics in the above scenario is revealed by the following simple picture (figure 21, bottom). During the time evolution, the two clouds that were initially in each of the potential wells spread, the negative rapidities moving to the left and the positive ones to the right. After some expansion time, at the central position, the positive rapidities from the left cloud meet negative rapidities coming from the right cloud. The resulting rapidity distribution is then double-peaked. The standard hydrodynamic theory cannot capture this feature, because it assumes local thermal equilibrium and the rapidity distribution of thermal states are single-peaked. This striking difference between GHD and conventional Euler hydrodynamics is clearly visible in figure 21, where the calculated Fermi occupation ratio ν(x, θ) is displayed for both theories.

The Newton cradle scenario has also been reproduced in (Schemmer et al 2019), by quenching the longitudinal potential from a double-well to a harmonic potential, see figure 22. One then initiates a dynamics similar to that studied in (Kinoshita et al 2006), with two clouds that oscillate and collide in a harmonic potential. The experimental results compare well with the prediction from GHD. For this scenario, conventional Euler hydrodynamics completely fails: it predicts the formation of a very sharp structure in the density distribution, corresponding to a large gradient of the density, that eventually leads to a shock, at times as small as about 30 ms. The agreement between experimental data and GHD is less good than for the scenario of figure 21. One of the reasons for this might be the effect of three-body losses. In the Newton cradle scenario, large peak densities are attained when both clouds superpose, and the three-body recombination process occurs, leading to atom losses. We find experimentally that the total atom number decreases by about 15% on the time evolution shown in figure 22.

Zoom In Zoom Out Reset image size

In the experiment of Schemmer et al (2019), GHD is tested in a very large atom cloud that contains thousands of atoms, whose longitudinal size, of the order of 100 μm, is very large compared to microscopic scales. In these conditions, GHD is clearly expected to be valid. In a more recent experiment by Malvania et al (2021), GHD is tested in a setup that uses a 2D lattice of 1D gases. In this experiment, the typical atom number N is as low as ten to 20 per 1D gas. Moreover, the (quasi-)harmonic longitudinal potential V(x) can be quenched very dramatically: at time t = 0, the amplitude of the trap can be increased by a factor of as large as 100, so the atom cloud is very strongly compressed. The validity of GHD is severely challenged in this situation, yet the experimental data, which are averaged over all the 1D tubes, are still correctly described by GHD over the first few oscillation cycles. We now discuss the results of Malvania et al (2021) in more detail.

The 1D clouds are prepared from a 3D BEC by adiabatically increasing the depth of the 2D lattice. When the 2D lattice depth is large enough, the gas decouples into independent 1D tubes. The temperature is extremely low, so that each 1D gas is close to the ground state of the Lieb–Liniger Hamiltonian in each tube. The dynamics is generated by suddenly increasing the longitudinal potential V(x) felt by the atoms in each tube. The longitudinal potential V(x), realized using an optical beam, has a Gaussian shape and its depth is increased by a large factor, of ten or 100 depending on the experimental data set. In sharp contrast with (Schemmer et al 2019), the initial cloud lies in the hard-core regime, with mean γ, averaged over the distribution of linear densities, as large as nine. Malvania et al (2021) measure the time evolution of the rapidity distribution (figure 23), using the technique of Wilson et al (2020) reviewed in subsection 3.3. This measurement is a global measurement: the rapidity distribution is integrated over positions x within each tube, and averaged over all the tubes.

Zoom In Zoom Out Reset image size

The results are shown in figure 23, for a quench of the 1D trap amplitude by a factor of 100. They are in excellent agreement with GHD predictions, which use the Lieb–Liniger ground state within the LDA as the initial state and the zero-entropy GHD equation (88) discussed in subsection 2.3 for the time evolution (figure 24). In the theoretical calculation it is found that, for an increase of the trap depth of the quasi-harmonic (Gaussian) potential by a factor of 100, the dimensionless repulsion strength γ, taken to be initially of order one in the center of the cloud, drops to a value of order 0.1 at the maximum of the compression: the variation of γ as a function of position and time over one compression cycle is very large. In this setup, the contour Γ_t (see discussion in subsection 2.3 and the caption of figure 24) gets deformed very quickly. Already in the first compression cycle, there are positions in the cloud where the rapidity distribution is no longer that of a thermal equilibrium state. Instead, the local state of the gas consists of a split Fermi sea. The appearance of such exotic rapidity distributions originates both from the non-trivial effect of interactions between atoms—which are stronger at the maximum of the compression, where the atom density is high and the cloud is far from the hard-core regime—and from the anharmonicity of the longitudinal potential. The presence of multiple Fermi seas rules out the possibility of describing the dynamics by standard hydrodynamic approaches (subsection 2.1).

Zoom In Zoom Out Reset image size

The measurement of the rapidity distribution, integrated over the cloud and averaged over all tubes, allows one to monitor the evolution of the rapidity energy. This is defined as E = ∫ρ(x, θ)θ²/(2m)dx dθ for a single 1D cloud, and it is averaged over all clouds in the measurement. The rapidity energy is the total energy, which is conserved, minus the potential energy associated with the confining potential. Since the potential energy oscillates as the cloud breathes, so does the rapidity energy E, as seen in figure 23(b). Moreover, besides the rapidity distribution, Malvania et al (2021) also measure the momentum distribution using a time-of-flight technique (see subsection 3.2.4), from which the kinetic energy E_K can be extracted. The time evolution of the kinetic energy is shown in figure 23(c) (notice that E_K is simply noted as K in the figure). The difference between the rapidity energy E and the kinetic energy E_K is the interaction energy. If the gas remained in the hard-core regime at all times, then one would always have E_K = E: the atoms would never be at the same position, so the contact interaction energy would be zero and the rapidity energy would be entirely in the form of kinetic energy. The narrow dip of E_K at the time when the cloud is most compressed clearly demonstrates that, at this time, the cloud leaves the hard-core regime. (Notice that it is important that it is the second moment of w(p), i.e. the kinetic energy, that is used as a diagnostic for the departure from the hard-core regime. If, instead, one studied the half-width of w(p), then its evolution would show pronounced dips at the time when the cloud is the most compressed, even if the gas stayed in the hard-core regime throughout the evolution (Atas et al 2017).)

To conclude this section, we stress that the good agreement found by Malvania et al (2021) between experimental data with small atom numbers per tube and the predictions of GHD is a topic that deserves further investigation. As discussed in the supplemental material of (Malvania et al 2021), the averaging over 1D tubes results in a smoothing of small atom number effects, which may explain the good agreement with the hydrodynamic theory to some extent. It would be interesting to investigate the question of the accuracy of GHD for single 1D clouds of a few atoms more systematically. On the theory side, this could be done by performing numerical simulations of a single 1D gas for very small atom numbers (similarly to the numerical results shown in figure 13, but for longer times).

Cold atom gases always suffer from losses. Different mechanisms for losses can be present, which are distinguished by the number of atoms K (K = 1, 2, 3, …) involved in each loss event.

• One-body losses (K = 1) can occur due to collisions with hot atoms from the residual gas in the vacuum chamber. This typically imparts a kinetic energy to the (cold) atoms that is sufficiently large that they leave the trap. One-body losses can also result from de-excitation for atoms lying in a metastable state, spin-flips to an untrapped magnetic state in the case of magnetically trapped atoms (Burrows et al 2017), collisions with energetic electrons (Labouvie et al 2016), or coupling to an untrapped state (Bouchoule and Schemmer 2020, Rauer et al 2016).
• Two-body losses (K = 2) occur, for instance, when atoms are not in their internal ground state and exothermic two-body collisions that change the internal state of the atoms are present (Traverso et al 2009, Yamaguchi et al 2008). The collision residues leave the trap because their internal state is not trapped, or because their kinetic energy exceeds the trap depth. In the presence of a laser, one could also have, starting from two nearby atoms, photoassociation toward excited molecules, which, after de-excitation, produce two very energetic atoms that leave the trap (Kinoshita et al 2005).
• Importantly, cold atom experiments always suffer from three-body losses (K = 3). This is due to three-body recombination, where a deeply bound molecule is formed. The binding energy, typically very large, is released in the form of kinetic energy, and the collision residues leave the trap (Söding et al 1999, Tolra et al 2004).
• In principle, loss processes involving more than three atoms also exist. In particular, losses involving K = 4 atoms have been reported in (Ferlaino et al 2009, Gurian et al 2012).

In the following we consider a general K-body loss process, for a fixed positive integer K. We now explain why the natural theoretical framework to model such a K-body loss process is the Lindblad equation (111) below.

Because the lost atoms can be viewed as escaping toward a reservoir of particles whose state is not being monitored, the atoms remaining in the cloud no longer follow a unitary dynamics. Instead, the density matrix of the remaining atoms in the gas (not to be confused with the rapidity distribution ρ(θ)) evolves according to a Lindblad equation (see equation (111) below). More precisely, the evolution of the gas under losses is described by a Lindblad equation if one assumes that the dynamics remains Markovian. This assumption holds if the energy width of the reservoir, E_res, which is the energy width spanned by the states of the continuum that are coupled to the trapped atoms by the loss process, is much larger than the energy width involved in the dynamics of the gas. In temporal terms, it corresponds to the fact that the intrinsic duration of the loss event, equal to ℏ/E_res, is much smaller than all other evolution time scales of the gas.

In experiments involving atoms in their internal ground state, the range of the interaction between atoms is typically much smaller than the typical distance between them. Losses can then be modeled by purely local processes: a loss event can occur only when K atoms are found at the same position.

Under these assumptions, the evolution of the density matrix of a uniform Lieb–Liniger gas of length L under losses is

where H is the Lieb–Liniger Hamiltonian (6). Here Ψ(x)^K is the operator that destroys K bosons at position x, and G is a constant characterizing the loss process, with dimensions of length^K−1  time⁻¹.

Calculating the evolution of the density matrix directly from the Lindblad equation (111) for more than a few atoms is, in general, an intractable task. Numerically, the size of the matrix quickly becomes prohibitively large. Therefore, the role of equation (111) is merely to give a formal definition of the theoretical problem one would like to solve. To make progress, further assumptions are needed in order to simplify the description.

In the context of this review of hydrodynamics, where one focuses on effective hydrodynamic descriptions valid at large scales assuming local relaxation, a natural assumption is to consider the limit of adiabatic losses. When the parameter G is small enough that the dynamics induced by losses is much slower than the relaxation time τ_relax of the system, the gas always remains in a relaxed state on long time scales. Its local properties are entirely described by the rapidity distribution ρ(θ), and the problem then boils down to computing the time evolution of ρ(θ). To lowest order in the small parameter τ_relax Gn^K−1 (where n = N/L = ∫ρ(θ)dθ is the atom density), the evolution of the rapidity distribution must be of the form

where F[ρ](θ) is some functional of ρ at time t, which needs to be determined from equation (111). The functional F[ρ] has been studied in (Bouchoule et al 2020), which we briefly review now.

The idea is to consider the adiabatic evolution of the conserved local charges Q[f], parameterized by some functions f (see equation (28)). To lighten the notation we simply write ‘Q’ for such a generic charge. Under Lindblad evolution (111), the expectation value of the charge changes as .

Because Q is the integral of a local charge density, Q = ∫q(x)dx, and because Ψ^K (x) and Ψ^†K (x) are local operators, the two operators between the brackets are local. Then we know that, after the relaxation time τ_relax, their expectation values relax to their values in a GGE, see subsection 1.6. This GGE is a diagonal density matrix which can be characterized by its distribution of rapidities ρ(θ), see subsection 1.6. Using the fact that , and writing for the expectation value of an observable w.r.t. the GGE density matrix parameterized by the rapidity density ρ, we get the evolution equation for the expectation values of the charges

We have used translation invariance, which implies that the integrand is independent of x, to go from the first to the second line.

Equation (113) determines the time evolution of the expectation values of all charges Q under adiabatic losses. We see that the problem of modeling losses boils down to computing the rhs of (113), namely the expectation value in a GGE of an operator of the form Ψ^†K (0)[Q, Ψ^K (0)], where Q is a generic conserved charge, and Ψ^K (0) is the operator that removes K bosons at the same position.

To connect this simple general result to the evolution of the rapidity distribution (112), we specialize Q to an operator which measures the rapidity distribution. To elaborate, recall that the charges Q[f] (28) are diagonal in the eigenbasis and that their expectation value in a Bethe state is . Formally, we can consider the charge Q[f] corresponding to the choice f(α) = δ(θ − α), which directly measures the distribution of rapidities ρ(θ). However, to ensure that the charge Q[f] has good locality properties, it is safer to work with a regularization of the Dirac delta function, δ_σ , of typical width σ and of total weight ∫δ_σ (θ)dθ = 1, which remains a smooth function of θ for any σ > 0 (for example a Gaussian of width σ). Then the choice f(α) = δ_σ (α − θ) defines a charge Q[f] = Q[δ_σ (. − θ)] which remains sufficiently local so that (113) applies. Thus, we see that the functional F[ρ] entering the evolution equation (112) must be given by

It is essentially equivalent to calculating the rhs of (113) for arbitrary charges Q, or to computing the functional F[ρ](θ) defined by (114). Either way, the difficulty lies in computing the expectation value of a specific local operator in a GGE.

Up to now, the following results have been obtained in connection with this problem.

• The functional F[ρ] can been evaluated numerically for a given ρ by performing a Monte Carlo summation over Bethe states (Bouchoule et al 2020), using exact formulas for the matrix elements of ψ^K (0) between Bethe states, see (Piroli and Calabrese 2015, Pozsgay 2011). The differential equation (112) can then be integrated numerically, see the example shown in figure 25. However, this procedure is numerically heavy: it typically takes a few hours to compute F[ρ] for a given rapidity distribution ρ on a single core (but the procedure can be trivially parallelized).
• The functional F[ρ] is known analytically in the ideal Bose gas regime,
  
  and also in the hard-core regime,
  
  where is the Hilbert transform of the rapidity distribution.The result for the ideal Bose gas is very simple, reflecting the fact that the rapidities are the momenta of the non-interacting bosons in that regime. The combinatorial factor K K! comes from the local K-body correlation (a consequence of Wick’s theorem), and the additional factor K simply comes from the fact that there are K atoms lost in each event.In contrast, the result for the hard-core regime for K = 1 is much more complex, even though it is also related to an underlying model of non-interacting particles, see subsection 1.8. In particular, one sees that F[ρ] is both non-linear in ρ(θ), and non-local in rapidity space (because the Hilbert transform depends on ρ at all values of the rapidity, not just at θ). These properties are thus expected to hold generically for finite repulsion strength.In the hard-core regime, F[ρ](θ) vanishes for K ⩾ 2 because two atoms (or more) can never be found at the same position; thus local K-body processes with K ⩾ 2 are suppressed.
• Instead of aiming directly at the time variation of the full rapidity distribution ρ(θ), another possibility consists of studying the variation of particular conserved charges. For instance, specifying f(θ) = 1 in (113) leads to the evolution of the density of particles,
  
  Similarly, specifying f(θ) = θ²/2 gives the evolution of the energy density,
  
  and so on. General expressions are available for the local K-body correlation g^(K)(0) as a functional of the distribution of rapidities ρ(θ), so that it is possible to compute dn/dt efficiently (Bastianello and Piroli 2018, Bastianello et al 2018, Pozsgay 2011). (The topic of the evaluation of g^(K)(0) in the Lieb–Liniger gas has a long history, see e.g. (Cheianov et al 2006a, 2006b, Gangardt and Shlyapnikov 2003a, Kheruntsyan et al 2003, Kormos et al 2011, 2009).) On the experimental side, the measurement of dn/dt was used by Tolra et al (2004) as a first demonstration that the inferred zero-distance correlation g^(K)(0) can take a value of substantially below one in strongly interacting 1D gases. More recently, a strong dependence of the effective loss constant K_eff = Kg^(K)(0) on the energy of colliding clouds has been reported for strongly interacting gases in (Zundel et al 2019). Such a behavior is compatible with the expected strong dependence of g^(K)(0) with the spread in rapidity space (Gangardt and Shlyapnikov 2003b). This behavior is also recovered by an analysis of the three-body problem (Mehta et al 2007).

Zoom In Zoom Out Reset image size

Until very recently not much was known about the evolution of other charge densities. Hutsalyuk and Pozsgay (2020) focused on the energy density and managed to compute the rhs of (118) as an explicit functional of the rapidity distribution ρ(θ). It is likely that their method can be generalized to some other charge densities, and this could perhaps ultimately lead to the variation of the full distribution of rapidities. At the moment, nothing is known beyond the atom density and the energy density though, and finding expressions similar to those of Bastianello et al (2018), Hutsalyuk and Pozsgay (2020), Pozsgay (2011) for the variation of other charges remains an open problem.

Again, we also refer the reader to the review of Bastianello et al (2021) in this volume for a thorough discussion of related recent results.

One striking effect of losses is found when investigating the evolution of the high-rapidity tails of the rapidity distribution ρ(θ). It was shown in (Bouchoule and Dubail 2020) that, under losses, ρ(θ) develops algebraically decaying tails, ρ(θ) ∼ 1/θ⁴ when |θ| → ∞. This is in contrast with more standard cases of rapidity distributions, e.g. those in thermal equilibrium states, which typically decay exponentially or even as Gaussians.

The physical origin of the development of those 1/θ⁴ tails in the rapidity distribution lies in the cusp singularity of the wavefunction when two atoms are at the same position, see equation (9):

Consider the case of one-body losses (K = 1), for simplicity. At a time t_l , the lth atom is suddenly removed from the system. Immediately after the loss event, the wavefunction of the remaining atoms is , where x_l is the position of the lost atom. This wavefunction, viewed as a function of x_i , still has a cusp singularity at x_i = x_l , even though there is no longer a particle at x_l . On the other hand, the eigenstates of the Lieb–Liniger Hamiltonian (6) for N − 1 atoms are smooth functions of z_i around z_i = z_l (for fixed values of x_j ≠ x_l for j ≠ i). Expanding the wavefunction over the eigenstates for N − 1 particles, one finds that the coefficients of the Bethe states decay as . The rapidity distribution is obtained by averaging w.r.t. to the squared amplitudes of these coefficients, so it must decay as 1/θ⁴.

Rapidity distributions decaying as 1/|θ|⁴ are not very common, however, there is at least one other known physical scenario where they appear: a sudden quench of the repulsion strength g. In (De Nardis et al 2014), the rapidity distribution ρ(θ) after a quench from g = 0 at zero temperature (BEC) to g > 0 is computed exactly, and it is found that ρ(θ) ∼ 1/θ⁴ for large rapidities. Notice that this is consistent with the discussion above based on cusps of the wavefunction: when the repulsion strength g is suddenly changed, the wavefunction immediately after the quench violates the cusp condition (119), resulting in an expansion over Bethe states with coefficients decaying slowly at large rapidities, as in the case of losses.

One important physical consequence of the presence of the 1/θ⁴ rapidity tails, worked out in (Bouchoule and Dubail 2020), is the breakdown of a famous exact relation between the tails of the momentum distribution w(p) in the gas and the ‘contact’ (Minguzzi et al 2002, Olshanii and Dunjko 2003),

Here the momentum distribution w(p) is normalized such that ∫w(p)dp = n. The relation (120) has been extended to higher dimensions and to fermionic gases or general Bose–Fermi mixtures with contact interaction (Tan 2008a, 2008b, 2008c). It is known as ‘Tan’s adiabatic theorem’ or simply ‘Tan’s relation’, and it has been studied extensively over the past 15 years, both theoretically—see e.g. (Barth and Zwerger 2011, Braaten and Platter 2008, Minguzzi et al 2002, Olshanii and Dunjko 2003, Tan 2008a, 2008b, 2008c, Werner and Castin 2012a, 2012b, Yao et al 2018)—and experimentally (Kuhnle et al 2010, Stewart et al 2010, Wild et al 2012).

However, in (Bouchoule and Dubail 2020), it is argued that the relation (120) breaks down when the rapidity distribution decays as 1/θ⁴ at large θ. In that case, the rapidity tail adds to the one caused by the two-body contact term. More precisely, introducing C_r = lim_{θ → ∞} ρ(θ) θ⁴ the amplitude of the 1/θ⁴ rapidity tails, the relation (120) is superseded by the relation

In most known stationary states of the gas, in particular in thermal equilibrium states, C_r = 0, so that Tan’s relation (120) holds. However, in a gas subject to losses, or after a quench of the interaction strength g, C_r > 0 and Tan’s relation is violated. Moreover, in (Bouchoule and Dubail 2020), the amplitude of the term C_r is found to be potentially much larger than . For one-body losses (K = 1) the ratio increases exponentially in time under lossy evolution, while it grows logarithmically for K = 2 and remains bounded (but not necessarily close to one) for K ⩾ 3.

In the asymptotic regime of a quasicondensate, the effect of losses can be investigated using the Bogoliubov description of the gas. The Bogoliubov approach is an approximate description of the gas, that constitutes a trivial integrable model since it amounts to independent bosonic modes: the integrals of motion are nothing other than the population in each mode. The effect of losses within the Bogoliubov approach was first investigated in (Grišins et al 2016, Rauer et al 2016), where the emphasis was put on the long-wavelength modes, the so-called phononic modes. These studies were then extended to all Bogoliubov modes, and also to K-body processes (Bouchoule et al 2018, Johnson et al 2017). The population of the phononic modes is expected to reach a value corresponding to a temperature T_p that fulfills

where μ ≃ gn is the chemical potential of the gas and α is a numerical factor, of order one, which depends on K. This prediction is in agreement with experimental studies made for three-body (Schemmer and Bouchoule 2018) and one-body (Bouchoule and Schemmer 2020) losses. Bogoliubov modes of shorter wavelength are, on the other hand, expected to reach a higher temperature (Grišins et al 2016, Johnson et al 2017). This peculiar state, in which the temperature of the short-wavelength modes is larger than that of the phonons, might be the reason for the experimental observation of long-lived non-thermal states (Johnson et al 2017).

The integrals of motion of the Bogoliubov model are however not the real integrals of motion, which are given by the rapidity distribution. The precise link between the Bogoliubov modes and the rapidities is still an open question. For short-wavelength modes, however, one can identify the population of the Bogoliubov modes with rapidities, as done already by Lieb (1963) in a seminal contribution. When doing this identification, the results of the Bogoliubov theory for the effect of losses on short-wavelength modes coincide with the expected behavior for the tails of the rapidity distribution (Bouchoule and Dubail 2020). For the phonons, on the other hand, there is no one-to-one correspondence with rapidities. Neither the validity of equation (122) at long terms, nor its compatibility with the time evolution of the rapidity distribution, has been established yet.

In (Rossini et al 2020), a lattice Bose gas with onsite repulsion (the Bose–Hubbard model) and onsite two-body losses is studied. In the limit of very fast losses, all configurations with more than one particle per site decay extremely quickly, so the slow dynamics occurs within the restricted subspace of configurations with at most one boson per site. In the effective dynamics of these emerging hard-core lattice bosons, two-body losses are still present, however, they occur on nearest-neighbor sites and they are very slow (García-Ripoll et al 2009). Thus, Rossini et al (2020) effectively work with a lattice hard-core boson model subject to adiabatic two-body losses. The corresponding (lattice) rapidity distribution is evaluated in a way that parallels the above discussion, see equation (112), and results analogous to those of Bouchoule et al (2020) are obtained.

As mentioned above, losses in the Lieb–Liniger gas are but one particular example of a mechanism that breaks the integrability of the underlying model. Other mechanisms have been studied, for instance the coupling between tubes in an array of 1D gases (Caux et al 2019), or dephasing (Bastianello et al 2020b). Let us also mention closely related works on integrable spin chains evolving under Lindblad evolution (Lange et al 2017, 2018, Lenarčič et al 2018), where an adiabatic limit analogous to the one discussed above is implemented using truncated GGEs, or the crossover from ballistic to diffusive transport induced by weak integrability breaking terms (Friedman et al 2020, Žnidarič 2020). These results, and more, are discussed in the review of Bastianello, De Luca and Vasseur in this volume.