Non-Equilibrium Thermodynamics of Harmonically Trapped Bosons

We apply the framework of non-equilibrium quantum thermodynamics to the physics of quenched small-sized bosonic quantum gases in a one-dimensional harmonic trap. We show that dynamical orthogonality can occur in these few-body systems with strong interactions after a quench and we find its occurrence analytically for an infinitely repulsive pair of atoms. We further show this phenomena is related to the fundamental excitations that dictate the dynamics from the spectral function. We establish a clear qualitative link between the amount of (irreversible) work performed on the system and the establishment of entanglement. We extend our analysis to multipartite systems by examining the case of three trapped atoms. We show the initial (pre-quench) interactions play a vital role in determining the dynamical features, while the qualitative features of the two particle case appear to remain valid. Finally, we propose the use of the atomic density profile as a readily accessible indicator of the non-equilibrium properties of the systems in question.


Introduction
Quantum gases offer a valuable platform for the study of quantum phenomena in interacting many-body systems. The availability of high-quality and reliable experimental control techniques and the low-level of external influences makes such systems excellent candidates for the simulation of quantum processes [1] and the exploration of the interplay between these and quantum critical behaviours.
Recent experimental progress has shown the possibility to observe non-equilibrium physics [2], which is quickly leading to the establishment of an experimental ultracoldatom framework for the exploration of complex phenomena, such as many-body localisation [3]. In particular, this framework can be used to test the recently developed ideas for the finite-time thermodynamics of closed quantum systems [4], which include tools for the quantification of thermodynamically relevant quantities, such as work and entropy, after a finite-time, non-equilibrium quantum quench. To date, this powerful formalism of non-equilibrium quantum thermodynamics has found only limited experimental validation and has been mostly applied to nuclear magnetic resonance settings [5]. However, notwithstanding the exquisite control available over such systems, they are hard to scale and offer only few possibilities for the inclusion of many-body effects. On the contrary, ultracold atomic gases offer solutions for both of these issues and we will therefore in the following study the connections between the phenomenology of non-equilibrium quantum gases and finite-time thermodynamics.
For this, we consider small-size gases of interacting bosonic atoms that are perturbed out of their equilibrium configuration. In particular, we focus on the ground state of interacting, harmonically trapped bosonic atoms in one dimension, and subject them to a sudden quench of the Hamiltonian parameters. Such a perturbation, which has been recently explored to characterise the occurrence of Anderson's orthogonality catastrophe [6,7,8], embodies the paradigm of a non-equilibrium process and has been shown to capture perfectly the complexity arising from quantum many-body effects in quantum spin systems [9]. In Ref. [10], the case of a fermionic system was addressed and developed, the study of the statistics of work in bosonic Josephson junctions was presented in Ref. [11], and the study of quenched attractive cold gases has recently been examined in Ref. [12]. Furthermore, Ref. [13] proposed realizing many-particle thermal machines using harmonically trapped bosons where a quantum advantage can be achieved.
Here we explore the non-equilibrium properties of these trapped ultracold atoms. First, we find the analytical expressions for relevant thermodynamical quantities, such as the Loschmidt Echo (LE) and the average work, for a single atom subjected to a quench in the trapping frequency and elucidate the relation between the entropy production and the tendency toward dynamical orthogonality, although we find that the evolved state is never completely orthogonal to the initial one for finite quenches in the trapping frequency. Enlarging the system, we study the role particle interactions play. For the infinitely repulsive Tonks-Girardeu molecule we find analytical expressions for the same figures of merit as in the single atom case. While the qualitative behaviour is consistent with the single atom case, we now see due to the strong interactions, finite sized quenches lead to dynamical orthogonality. For finite interactions, we find numerically that the average in time of the entanglement and the (irreversible) average work are linked. This is a crucial result in this paper, as we find that even in the smallest interacting system, there exists a link between the establishment of correlations (entanglement) and the work done. We relate these features to the spectral function, which gives the fundamental excitations that dictate the dynamics. We complete our work studying the smallest possible mixture of ultracold bosons with both inter-and intra-species interactions: two atoms of type X and one atom of type Y and employ the same figures of merit as in the previous cases. We show a relationship between the LE and the density profile, which is an experimentally measurable quantity. Furthermore this gives information about the correlations established among the atoms and indicates a trend for larger gases subjected to a quench, a situation of great theoretical interest and experimental relevance [2]. Finally, we remark here that, using symmetry arguments, we show that the mixture of three atoms behaves identically to that of the three indistinguishable-boson system if both the inter-and intra-species constants are equal.
Our presentation is organised as follows. In Section 2 we provide a brief introduction to non-equilibrium quantum thermodynamics, focussing in particular on the consequences arising from a sudden quench. This formal framework is then applied in Sections 3-5 to the single atom, trapped molecule and three atom mixture when subjected to the quench of their Hamiltonian parameters. Finally, Section 6 draws our conclusions, while a set of technical considerations and details are presented in the Appendices.

Non-equilibrium thermodynamics of quantum quenches
In the following we will first briefly summarise the key notions of finite-time thermodynamics in closed quantum systems. For this we consider a system whose Hamiltonian, H, depends on an externally controlled, time-dependent work parameter λ t . The system is assumed to be in contact with a bath at inverse temperature β for a time long enough to have reached equilibrium. At t = 0, the system is detached from the reservoir and its energy is changed by modifying the value of the work parameter from λ 0 to λ τ . The evolution is accounted for by the unitary propagator U τ . As the system is detached from the surrounding world, such a change of energy can only be interpreted as work done on/by the system, which can be characterised by introducing the work probability distribution [4] Here E n (E M ) is the n th (M th ) eigenvalue of the associated eigenstate |n (|M ), of the initial (final) Hamiltonian. Moreover, p(n, M ) = Tr[|M M | U τ |n n| ρ s |n n| U † τ ] is the joint probability of finding the system in |n at time t = 0 and in state |M at time τ , after the evolution by the time-propagator U τ . Obviously, such a joint probability can be decomposed as p(n, M ) = p 0 n p τ M |n , where p 0 n is the probability that the system is found in state |n at time t = 0 and p τ M |n is the conditional probability to find the system in |M at time τ if it was initially in |n . Therefore, P (W ) contains information about the statistics of the initial state and the fluctuations arising from quantum dynamics and measurement statistics. The characteristic function of the work probability distribution of P (W ) is defined as [4] with ρ 0 eq being the initial equilibrium state of the system and H(τ ) the Hamiltonian of the system when the work parameter takes the value λ τ .
For a quasistatic process, the change in free energy ∆F of the system is equal to the average work done on/by it. The former can be written as ∆F = ∆E − ∆S/β with ∆E being the change in energy and ∆S the corresponding entropy variation. On the other hand, if the process is fast (i.e. not quasistatic), then the relation W ≥ ∆F holds, accounting for the fact that part of the work performed on/by the system is dissipated due to the abrupt nature of the transformation. By introducing the standard definition of non-equilibrium entropy production where Q is the average heat exchanged with the environment, we find for a closed, unitarily evolving system This allows to quantify the irreversible nature of a given process in terms of the discrepancy between ∆F and W . The definition of Σ allows for the consideration of the so-called irreversible work which will be extensively used in this work. A general approach to irreversible entropy in open quantum systems (including non-equilibrium ones) can be found in Ref. [14], while a different quantifier which is based on the use of adiabatic transformations (rather than the implicit isothermal ones considered here) has been proposed in Ref. [15]. A very useful lower bound to the non-equilibrium entropy production, Σ , can be based on the unitarily invariant Bures angle (see Ref. [16]). For arbitrary density matrices ρ 1,2 , the Bures angle is defined as B = arccos F (ρ 1 , ρ 2 ) with F (ρ 1 , ρ 2 ) the fidelity between the two states. Using this we find where B eq (τ ) = arccos F (ρ τ eq , ρ τ ) is the angle between the non-equilibrium state ρ τ of a closed quantum system, and its equilibrium version ρ τ eq . Eq. (6) defines a thermodynamic distance that is valid arbitrarily far from equilibrium, and can thus be used to characterise the departure from equilibrium following an arbitrary driving process.
Since we are interested in examining the dynamics of a cold atomic system after a sudden quench, we will make use of the fact that we start from the ground state of a given system and that its state remains pure throughout the whole dynamics. We will use the spectral decomposition of the initial and final Hamiltonians of the system H α = j E α j ψ α j ψ α j with α = I (α = F ) denoting the initial (final) Hamiltonian operator. Here E α j is the j th eigenvalue of H α with associated eigenstate ψ α j . A key figure of merit for our system is the Loschmidt echo (LE), which is defined as where we have assumed that the initial state of the system |Ψ 0 coincides with the ground state ψ I 0 of H I . The LE is closely related to the characteristic function of the probability distribution of the work done on/by the system upon subjecting it to the quench considered here. In fact, for a sudden quench we have that U τ = 1 1, with 1 1 the identity operator, and thus χ(u, τ ) ≡ χ(τ ) = Tr[e iuH(τ ) e −iuH(0) ρ 0 eq ]. Here, we are taking ρ 0 eq = ψ I 0 ψ I 0 and, by using the identifications H(0) = H I and H(τ ) = H F , we find L(t) = |χ(t)| 2 with From this the average work is given by while from the definition of irreversible entropy production Eq. (4), we can introduce a quantifier of the dissipated work due to the non-quasistatic nature of the quench.

Single trapped atom
Let us start by considering the simplest possible scenario of a harmonically trapped single atom in one dimension. The Hamiltonian of the system reads with m the mass of the atom and ω 1 the frequency of the trapping potential. In the following we will consider a quench in the trapping potential frequency ω 1 → ω 2 and, in order to simplify our notation, we rescale the position of the atom with respect to the ground state length a ho = /mω 2 , and its energy with respect to ω 2 . The dimensionless initial HamiltonianH = H/( ω 2 ) of the system then reads wherex = x/a ho and = ω 2 /ω 1 . The results presented in this section are closely related to those presented in Refs. [10,16]. However an explicit re-examination of these calculations will be useful for understanding the upcoming sections. To calculate any of the above quantities requires determining the overlap between the initial (ground) state and the eigenstates Eq. (9)] which in this case is done using the known wavefunctions with the associated energies E I 0 = 1/(2 ) and E F n = (n + 1/2), and where H n (y) is the Hermite polynomial of order n and argument y. Exploiting the fact that which is valid only for even values of n and directly leads to Note that the average work is dimensionless in our chosen units. As our system is pure, the free energy difference is simply the difference between the initial and final ground state energies, ∆F = 1 2 ( − 1), and thus The lower bound of the entropy produced dynamically is given by and in Fig. 1 we show the behaviour of these quantities for different representative values of the quench. Examining panels (a) and (b) we see an oscillating pattern stemming from the harmonic oscillator dynamics and find that the behaviour of the lower bound on the irreversible entropy is strongly correlated with the behaviour of the LE. Its value at a given time grows with the strength of the quench as a consequence of the fact that, as grows, the ground state ofH F becomes increasingly different from ψ I 0 . When examined against time, we find that the maximum entropy production is achieved in correspondence with the minimum value of L(t). At this time, the state of the system is as different as possible from the the initial one, which coincides with the maximum irreversible entropy generated. By inspection of Eq. (14), we note that full dynamical orthogonality never occurs for finite quenches of the trapping frequency. In Fig. 1 (c) we show the behavior of the average work, irreversible work and the free energy against the strength of the quench. While, naturally, all quantities grow with increasingly large quenching strengths, W irr grows much more significantly than ∆F , precisely inline with the increasingly large maxima attained in the entropy produced for larger quenches.

Trapped molecule
We now move to a more complex multi-particle system and examine the effect of atomic interactions on the quantities discussed above. For this we consider two atoms of equal mass m, jointly trapped in a harmonic potential of frequency ω 1 and then quench the trap frequency to ω 2 , as before. The Hamiltonian model of the system scaled in the same way as in Eq. (11) reads and corresponds to the initial Hamiltonian for k = 1 the the final one for k = 2. The parameter g is the coupling constant, which characterises the boson-boson interaction in the limit of low temperatures. We will assume the interactions to be repulsive throughout this work and therefore g to be positive. Note that the rescaling leads to a coupling constant in units of a ho ω 2 and this again allows to define the parameter = ω 2 /ω 1 . It is important to note that while in principle we have the freedom to quench the interaction strength or the trapping frequency, quenching ω will implicitly lead to a quench in the interaction strength, if it was initially finite. We will show below that this can be dealt with cleanly for two-particle systems, however for larger system sizes more care must be taken (see Sec. 5).

Tonks-Girardeau pair of atoms
The coupling constant can range from zero to infinity, which is the so-called Tonks-Girardeau (TG) limit. Since in the non-interacting case all the results of the previous section still apply, we will in the following start by carefully studying the TG limit, where the atoms behave as hardcore bosons and are readily amenable to analytic treatment. The wavefunction of the system can be split into its centre-of-mass (COM) and relative (REL) coordinates, |ψ α n = |η α n |ϕ α n , where |η α n (|ϕ α n ) refers to the COM (REL) degree of freedom. The LE and characteristic function depend on the overlap between the initial and final wavefunctions, as before. In fact, the wavefunctions of the COM terms are precisely the same as in the single atom problem and therefore the overlap is given by Eq. (13). However, due to the infinite interaction each even REL state becomes degenerate with the next higher lying odd state, such that it is sufficient to work only with the odd states. The required initial and final eigenstates are wherex =x 1 −x 2 , with the associated energies A I 1 = 3/2 and A F 2n+1 = (2n + 3/2). We can then express the average work and the LE in terms of these functions as with the overlap between the REL wavefunctions given by Using this expression we can calculate the LE using Eq. (20) and Σ B using Eq. (6). The infinite sums can be evaluated explicitly, and the final expression for the average work turns out to be exactly four times the single-atom work showing that in this regime work is an extensive quantity as we have performed a quench of the trap frequency for two atoms [compare to Eq. (15)] ‡. This can be further understood from the fact that in the TG limit the bosons behave as non-interacting fermions [10]. In Fig. 2 we show the behaviour of L(t) and Σ B . While the qualitative behaviour is consistent with the single atom case, we see that the effect of the interactions is to magnify these features as we now must account for many-body effects. In particular we also see the system periodically evolves into orthogonal states as the interacting twobody system can be moved further out of equilibrium, as evidenced by the vanishing values of LE which are already less than 10 −2 for = 6. ‡ In the TG limit we find the average work scales as N 2 times the single atom average work.

Finitely interacting pair of atoms
When the frequency of the trap is quenched while the interaction strength is finite, i.e. between the non-interacting and the TG limit, the resulting complexity requires we use a numerical approach to study the system, and we refer the reader to Refs. [17,18] for details of the recipes employed here. We remark that one could equally employ the known analytic solutions in Ref. [19], however as the infinite sums and overlaps must still be performed numerically such an approach is equivalent to the one employed here. The dynamics well characterised by the von Neumann entropy (vNE) of the atomic state, which is given by where ρ is the reduced density matrix (RDM) for one of the atoms. Since we consider the system to be always in a pure state, the vNE is a good measure of the entanglement in the system. For finite values of g, quenching the trap frequency implies also a quench of the interaction strength between the atoms and the evolution of the vNE will therefore be determined by a competition between these two mechanisms. § The resulting behaviour is shown in Fig. 3 for various values of the quench amplitude and the coupling strength g. For small-amplitude quenches and weakly interacting atoms (lower blue curve in Fig. 3(a)) the vNE oscillates, as expected for a quench in a harmonic oscillator, with an amplitude modulation due to the effective interaction quench. At larger strengths of the quench (lower blue curve of panel (b)) this behaviour is strongly modified. The absolute values of the entanglement increase, as the system becomes more strongly correlated, but there is no longer evidence of regular oscillations as the spectrum has become anharmonic due to the interactions. Looking at the spectral function of the out-of-equilibrium state, A(ω) = 2Re χ(t)e iωt dt [20], we can identify § In the TG limit the vNE is constant for all quenches, due to the infinite value of the interaction. We remark S = 0.68275 (two particles) and S = 1.0574 (three particles).    these different excitation frequencies inherent in the evolution, see Fig. 4. The majority of the motion is governed by the quasi particle peak at the ground state energy of the quenched state E F 0 + A F 0 , and smaller contributions stem from combinations of COM and REL even states at higher energies (there are no contributions from the odd states). At larger interaction strengths the high energy peaks in the spectral function approach each other, as the system becomes doubly degenerate in the TG limit, see Fig. 4 (b). This causes larger interference effects that are apparent in the entropy evolution (c.f. the upper red curves in Fig. 3 (a) and (b)). The periodic nature of the revivals is destroyed for the large quench of = 5 due to the broadening of the high energy peaks (see inset of Figs. 4 (c) and (d)) which highlights the chaotic dynamics of the strongly quenched state. A clear signature of interference of the COM and REL states is manifested in the appearance of Fano-resonances in the spectral function, a feature which is typical in systems with two scattering amplitudes which overlap (see all insets in Fig. 4) [21].
To evaluate the finite-time thermodynamics of the system following the quench, we show in Fig. 5 the LE for the same parameters as used in Fig. 3. The periodic nature of the echo is visible for the small quench ( = 2) exhibiting breathing dynamics which are a consequence of the non-trivial energy shifts the system caused by the interaction [22,23,24]. In this case the frequency of the fast oscillations are given by the energy of the dominant quasi-particle peak in the corresponding spectral functions, which are governed by the ground state energy of the quenched state. The slower frequency envelope is a consequence of the finite interactions which cause a splitting of the first excited state into a pair of COM and REL states of comparable energy (see inset of Fig. 4 (a) and (b)). The interference of these two states causes the breathing observed in the LE and the difference in their energy controls the breathing frequency. For the larger quench ( = 5) it is clear that this beating is destroyed resulting in orthogonality and further destructive interference effects from the contributions of higher energy states.
In the remainder of this Section we concentrate on the behaviour of the irreversible work W irr and its connection with the amount of entanglement created between the two atoms and the amount of work dissipated. To this aim, we calculate the mean vNE which is time-averaged over an interval τ that is long enough to include many periods of oscillation. In Fig. 6 we show a clear qualitative link between S and both W and W irr . We find all quantities follow the same qualitative behaviour (predominately linear, although exhibiting a slight systematic curvature) as a function of the interaction strength, thus suggesting a causal link between these figures of merit. This is intriguing as it evokes a possible cost function-role played by the thermodynamic irreversibility in the establishment of quantum correlations within an interacting system. It is worth noting that a number of recent works have managed to establish a rigorous link between the appearance of correlations (both quantum and classical) and the associated thermodynamic cost [25] (albeit applying a separate formalism to the one considered here), that is complementary to our analysis. Indeed, our results suggest a significant role for thermodynamic work in the establishment of quantum correlations between the atoms. Such a connection has previously been highlighted in the context of spring-like coupled bosons, but it was limited to Gaussian states and quadratic evolutions [26]. Our results go significantly beyond these restrictions as they address the case of contact-like bosonic interactions, which are highly experimentally relevant. In the following we will use this study on two coupled bosons as a benchmark for contactlike couplings among multiple atoms.

Three trapped atoms
We now extend the Hamiltonian model addressed so far to the more complex case of a one-dimensional mixture of two identical bosons of the same species X, whose coordinates will be indicated as x 1 and x 2 , and one impurity atom of a different species Y, with coordinate y. We assume that all atoms have the same mass m and are trapped with the same oscillator frequency ω. The interactions are of contact form and characterised by the intra-and inter-species coupling constants g X and g XY . In this situation, the Hamiltonian reads where the dimensionless coordinates are defined by rescaling energies and coupling rates as above. Note that the eigenfunctions of Eq. (26) have to be symmetric with respect to the exchange of the X bosons, but no symmetry restriction for the interchange of the X atoms with the Y atom exists. A detailed study of the forms and properties of the eigenstates Ψ(x 1 , x 2 , y) of Eq. (26), focusing on the degeneracies of the spectrum, is given in the supplementary material. The RDM for the single atoms can be calculated by taking the partial trace of the state of the three-atom system over two of the atoms and is given for one of the X atoms Analogously, for the impurity atom of species Y, we have Here the functions f X,Y k are the natural orbitals that diagonalise the RDMs with natural orbital occupations λ X,Y k . The vNE of the impurity or of one of the atoms of species X can then be calculated as As mentioned previously, for any finite interaction strength, a quench in the trapping frequency results in an effective change of the interaction strength. While for the trapped molecule a quench of either parameter led to the same behaviour, the presence of two coupling constants in the current system complicates the matter. In order to clearly identify the role atomic interactions play in the dynamics, we therefore restrict ourselves in the following to quenching the interaction strengths directly. It is worth noting that, when g X = g XY = 0 or g X = g XY ≈ ∞ , the behavior of the threeatom system is qualitatively the same as the single and two atom cases, respectively.
For our two-component system there are in principle five different strategies for quenching the coupling constants (g X , g XY ), but in the following we will focus on just two strategies as illustrated in Fig. 7, as these encompass the most relevant features of the thermodynamic properties of the system (see Appendix A for a discussions of the other strategies) A: (0, 0) → (g, g), B: (0, 0) → (0, g).
The energy spectra corresponding to the adiabatic version of these quenches are shown in Fig. 8 [27,28,29,30] and the visible differences are directly related to the distinct In our simulations g = 20 is large enough to effectively reach the TG limit. outcomes of the particular quenching protocols. The spectrum for situations adhering to strategy A shows the emergence of three-fold degenerate states in the limit of infinite interaction strengths, while in the same limit strategy B leads to a two-fold degeneracy. The differences among these energy spectra make it clear that each quenching protocol provides distinct thermodynamic and density evolutions, which we will discuss below. We will show that both the entanglement and the thermodynamical quantities of interest are strongly affected by the chosen quenching strategy and we highlight a link between these quantities and the behaviour of the atomic density profiles. We remark that a consistent feature of all results is the periodic nature of these functions, which is due to the gas refocussing in its external harmonic oscillator potential at approximate multiples of the trap frequency.

Thermodynamic quantities
For strategy A the behaviour of the vNE is shown in Fig. 9 (a). It is worth noting that the behaviour for this quantity is the same whether two X atoms or one X and one Y atom are traced out. This can be understood by realising that when quenching the two coupling constants to the same value, the system behaves as one composed of three indistinguishable atoms. It indicates that only the energy levels that both systems have in common are affected by the quench and there is no difference between the respective RDMs. A more formal argument using group theory is presented in Appendix B. We can also see that the entanglement increases with the amplitude of the quench and that dips appear with a periodicity close to multiples of the trapping frequency. These are more harmonic and narrower when the system is quenched close to the TG regime (g = 20), as it can then be mapped onto non-interacting fermions. In Figs. 9 (b) and (c) we show the evolution of the vNEs for strategy B, which now differ for the atoms of different species. Again, the vNE grows for stronger quenches, but the periodic dips are less pronounced compared to the ones observed with strategy A. This is not surprising, as the spectrum for this situation is denser and therefore more states contribute to the evolution, which makes perfect refocussing less likely.
The LEs following the quenches are shown in Figs. 10 (a) and (c) and, similar to the behaviour of the vNE above, they display regular revivals, whose period moves closer to the harmonic oscillator time scale of multiples of ωt/π for large quench amplitudes (when the system is quenched to the TG regime). The different behaviours for each quench are most clearly visible by looking at the behaviour of the LE around the revivals. Let us first discuss quenches into the TG regime, for which the LE always displays a periodicity. For strategy A the revivals have a period ωt/π, which is due to the evolution being governed by the energy differences between the quenched states and the initial state, ∆E = E n − E 0 = q, where q is an integer as the trapping frequencies do not change. For case B, the LE is more varied and shows periodicities with 4ωt/π and 2ωt/π. This is a consequence of the symmetry considerations discussed in the supplementary material, which results in energy differences (a combination of integer and half integer ∆E) for strategy B. For weaker quenches, g < 20, the LEs show complex temporal patterns which is a sign that the energy spectrum is not as degenerate as in the TG limit.
The correspondence between the LE and the characteristic function of the work distribution allows us to evaluate the average work done on the system as a result of the quench (see Figs. 10 (b) and (d)). Quite interestingly, we find that W depends linearly on the quenching strength in both settings, while the irreversible work W irr produced in such non-quasistatic processes behaves linearly only at larger values of the quench amplitude. This arises from the impossibility of the system to generate increasing amounts of 'useful' work as the strength of the couplings grows. Indeed, the free energy difference ∆F = ∆E = E 0 −E 0 (see insets of Figs. 10 (b) and (d)), levels off at large values of the interaction strengths, showing that the system soon saturates its capabilities to produce work that can be usefully extracted. Large quenching amplitudes are thus associated with an increasing degree of irreversibility. These considerations allow us to identify an optimal configuration of quenching, associated with moderate quenching amplitudes, that without requiring sophisticated control techniques is able to attain the maximum allowed useful work without generating an unbounded amount of thermodynamic irreversibility.
The lower bound of the average irreversible entropy produced as a result of the quenches is shown in Fig. 11 (a)-(b) and it is apparent that Σ B closely follows the same temporal trend as L. This is not surprising, given that for a time-independent Hamiltonian and a system prepared in the ground state of the initial Hamiltonian, the LE coincides with the state fidelity, and the latter directly enters the definition of the Bures angle. In Fig. 11 (c) we show a qualitative link between the behaviour of the LE and the single atom density of the system. The visible regular dependence suggests the possibility that the behaviour of thermodynamically relevant quantities, such as the irreversible entropy production, can be inferred from experimentally accessible figures of merit such as the density profile, which we will explore in the next Section. (a) (b) Figure 12. Evolution of the density profile ρ X (x) for a quench of strategy A for (a) g = 2 and (b) g = 20.

Density evolutions
The evolution of the density profile after quench A for two different strengths is shown in Fig. 12 and we note that, as expected, the profiles obtained from the RDM for X and Y are identical. We can see that the system is localised around x = 0 and oscillates at a frequency which depends on the strength of the quench. In fact, these oscillations mirror the appearance of the dips in the vNE and the peaks in the LE (see Figs. 9 (a) and 10 (a)). A larger amplitude of the quench leads to narrower revival peaks, which again corresponds to tighter dips (peaks) in the vNE (LE). For the situation of strategy B, where g X = 0 and g XY is quenched, it was shown above that the LE and the vNE show complex dynamics for small quench values. As the two species are distinguishable in this case, we show in Fig. 13 the density evolutions for an atom of species X and Y separately. For small g XY , the density profile for species X remains localised in the center of the trap, while the impurity can spread to the edges of the distribution for the X atoms, forming a double peaked structure at certain times [27]. The intrincate structure of the LE in this case can be related to the irregular temporal evolution of the Y atom compared to the more periodic oscillations of the X atoms. However for larger g XY , the periodicity of the evolution becomes more regular for both species as the energy structure becomes more degenerate, resulting in complementary trends visible in the corresponding vNE and LE (see Figs. 9 (c) and (d) and 10 (c)).
A connection between the dynamics of the density profile and the evolution of thermodynamically relevant quantities would allow insight into the finite-time thermodynamics of non-equilibrium processes without requiring measurements of hardto-access variables [4]. In particular, Fig. 11 (c) highlights the possibility of postprocessing data acquired on the density profile at the centre of the trap within one period of the evolution to infer the corresponding value taken by the LE.

Conclusions
We have studied the finite-time thermodynamics of a small-sized gas of interacting bosonic atoms subjected to a sudden quench of the Hamiltonian parameters. By first reviewing the out-of-equilibrium dynamics of a quenched single atom we have confirmed that larger quenches lead to larger amounts of entropy produced which implies an increase in the amount of work and irreversible work injected into the system. However, the system only appears to evolve into an orthogonal state when the quench is infinitely large. For the two-atom case we have investigated the interesting role that particle interactions play. In particular, starting from the analytically tractable Tonks-Girardeu regime, we noted that for such infinitely repulsive bosons, the strong interaction enhances the entropy production and the system can now exhibit full orthogonality. We also established the extensive nature of work in this system. For finite interactions between the atoms, we have shown a clear qualitative link between the amount of (irreversible) work performed on the system and the increase in the degree of interatomic entanglement. Moving into the multipartite case, and despite the significant increase in the complexity of the problem (as evidenced by the range of inherently different dynamics and sprectra observable simply by altering the initial interaction strengths), we have highlighted that the qualitative features of the two particle case appear to persist, i.e. the initial interactions strongly dictate the dynamical features. Finally, we have shown that the behaviour of the atomic density profile of a single atom can be a useful tool in exploring the non-equilibrium properties of a system, even in the case of complex multipartite systems.  Figure 14. (a) Sketch of the three other quenching strategies considered for the threeatom system. The quenches D and E could in principle lead to the same final values for g X and g XY . (b) Corresponding energy eigenspectra for these systems with different scattering symmetry requirements. The sprectra relating to quench E fixes the intraspecies coupling strength, g X , and changes the inter-species coupling strength, g XY , while for the spectra C and D the opposite is the case. In situation C only scattering with a symmetry requirement is present and it should therefore be compared with situation B where all symmetry requirements are absent. Finally, in situation D the symmetric scattering is fixed and the scattering without the symmetry condition is varied, which should be compared to the situation E where the opposite is the case. Figure 15. (a) and (b) Entropies S X and S Y following the quenching strategy D, where g XY = 20 and g X is quenched from 0 to 2 or 20 (black and red curves respectively). common with strategy B, with the important distinction that the atoms are initially interacting at t < 0, and therefore many natural orbitals have non-zero occupation even before the quench [27,31,32].
The dynamical behaviour of the vNE for strategy D is shown in Figs. 15 (a) and (b). As the interaction between the two species is already large at the beginning, the initial values of the correlations are now finite and the quench increases them to a similar level as in the cases considered in the main text. Similarly, periodic dips appear again around the refocussing time of the harmonic trap. Figs. 16 (a) shows the LE that, similar to strategy B, exhibits periodicities around 4ωt/π and 2ωt/π due to the energy differences of ∆E = q + 1/2. We see from panel Fig. 16 (b) the behaviour of the average (irreversible) work and free energy are qualitatively the same as for strategy B. Fig. 17 shows the evolution of the density profile for a quench of g X taken from 0 to a value that is either much smaller than g XY , or comparable to it. The phenomenology is strikingly different for both cases. For g X g XY the density profile for the X atoms is peaked at center of the trap, while the impurity Y has a double-peak structure which is localised at the sides of the density of X. This distribution shows only a weak temporal change and the separation is maintained, which corresponds to the flat vNE of the individual species (see Fig. 15 (a) and (b)). The case in which the final value of g X is comparable to the inter-species coupling rate shows more pronounced temporal Here g XY = 20 is constant and we consider g X quenched to g = 2 and g = 20 [top and bottom row of plots, respectively]. Same color-scale as in Fig. 12.
oscillations, which are also seen in the behaviour of the vNE . The atomic species are strongly correlated regardless of the strength of the quench. However, large quenching amplitudes result in dips of the vNE at the refocusing time of the density profile that are much less pronounced than those occurring at small final values of g A .

Appendix B: Comparison with indistinguishable atoms
It is also interesting to compare the evolution of the density, the vNE, and occupation of the natural orbitals shown in Figs. 9 and 12 with that of a system of three indistinguishable atoms. The energy spectrum for such a system is shown in Fig. 8 ( * ), and clear differences from that of two atoms plus a third distinguishable one are visible. In the latter case triple degenerate states occur in the limit g X , g XY → ∞ (see Fig. 8 quench A) and a discrete group theory analysis presented in Refs. [27,33] showed that these three degenerate states belong to different irreducible representations of the group. The discrete group to which these solutions belong is the discrete rotational group of order 2, C 2 , restricted by the bosonic symmetry under interchange of the two indistinguishable atoms. Indeed, for all finite values of the coupling constants, all wave functions can be classified according to the possible irreducible representations of this group. In Ref. [27] it was shown that the ground state for the three indistinguishable atoms was the same as the 2+1 case for all values of g = g X = g XY . This is not too surprising, as these solutions obey all symmetries under interchange of two atoms required by the three-indistinguishable atoms which coincide with the ones required by the corresponding irreducible representation of the group in the 2+1 system. By comparing the dynamical evolution of the density, the vNE, and the occupations of the natural orbitals for a system of three-indistinguishable atoms with the 2+1 setting, we find that they all coincide. The reason for this is that the initial state is a noninteracting Gaussian state with certain symmetries, which corresponds to the absence of a change in the sign of the wave function when any pair of atoms interchanges their coordinates. Therefore it belongs to a definite irreducible representation of the C 2 , restricted by the bosonic symmetry under the interchange of the two X atoms [27]. If the system has this symmetry initially, the dynamical evolution has to conserve it, so only part of the energy spectra in the 2+1 case plays a role in the evolution. This part of the energy spectra is exactly the same as in the case of three indistinguishable atoms.