Equilibration and nonclassicality of a double-well potential

A double well loaded with bosonic atoms represents an ideal candidate to simulate some of the most interesting aspects in the phenomenology of thermalisation and equilibration. Here we report an exhaustive analysis of the dynamics and steady state properties of such a system locally in contact with different temperature reservoirs. We show that thermalisation only occurs 'accidentally'. We further examine the nonclassical features and energy fluxes implied by the dynamics of the double-well system, thus exploring its finite-time thermodynamics in relation to the settlement of nonclassical correlations between the wells.

Here ˆf  describes the free evolution of the atomic systems in the two wells, each occurring at the rate set by the single-atom energy ω j , and with ,ˆ † a a j j the associated annihilation and creation operators for each well. The Hamiltonian term ˆs i  accounts for the self-interaction (at rate U) between atoms occupying the same well, while ˆt  stands for the tunnelling term, which occurs at rate  . We will focus mostly on the tunnelling-dominated regime associated with U = 0. However, the interaction-dominated regime corresponding to  = 0, and the intermediate regime will also be addressed. The focus of our investigation will be the phenomenology of thermalisation of the system, both at the single and two-well level. We remark that model, Eq. (1), can be realised in a variety of settings including superconducting Josephson junctions 16 , trapped ions 17 , bimodal optical cavities 18 and optomechanical setups 19 .
While important insight will be gathered by addressing the unitary evolution induced by considering  , the overarching goal of this work is the study of the open-system evolution created by the contact of the two wells with their respective reservoirs. We are interested in addressing the dynamics induced by the master equation 20 where we have introduced the overall-system density matrix at a generic time t,  t , and the Lindblad super-operators n a a aa n a a a a 1 1 2 1 2 3 j t j j j t j j j t j j j t j j j t which describe the incoherent particle-exchange process (occurring at rate γ j ) between a well and the respective reservoir (assumed to have a thermal occupation number n j ). Eq. (2) is the key equation in our analysis to follow. We remark that in certain working conditions this description of the dynamics is not always valid. In particular, when the scattering length of the BEC is large, non-Markovian dynamics can play an important role 21 . We therefore assume that the scattering length is sufficiently small to ensure the validity of the Markovian approximation 21 .

Results
Exact solutions of the tunneling-dominated regime. In order to gather insight into the basic coherent processes of the system in the case of tunnelling-dominated regimes, we set U = 0 in  and address the unitary evolution first. We define the canonical quadrature operators , , ,ˆˆx p x p { }  )ˆˆˆP x p x p 1 1 2 2 and the expectation value of such vector (calculated over the state of the system), which bear full information on the state of the system. Both are readily gathered as ] 0 the covariance matrix [vector of phase-space displacements] of the initial state of the system and σ ( ) t θ θ θ θ We now concentrate on the situation where the particles in each well are initially at thermal equilibrium with their local reservoirs. The initial covariance matrix will thus be that of a two-mode thermal state a a j j j the mean number of particles in the j th well. For such an initial state, the phase-space displacements are all null and full information on the evolved state is provided by the covariance matrix , i.e. the system does not evolve in time and the two wells remain at their thermal equilibrium, notwithstanding the tunneling. This is a clear interference effect. Moreover, for identical single-atom energy in each well, i.e. Δ = 0, c 1 is null, showing that the position [momentum] In general, such correlations do not imply necessarily the setting of entanglement between the wells 24 . Indeed, the tunneling term of the Hamiltonian,  t in Eq. (1), can generate entanglement only when the state of at least one of the wells is sufficiently non-classical. In the context of our investigation here, this basically implies the preparation of squeezed states of the wells 25 . This can be understood by noticing the formal analogy between ˆt  and the generator of a two-mode mixing transformation and considering, for the sake of argument, the resonant case Δ = 0. Under such conditions, moving to the interaction picture with respect to ˆf  , the time evolution operator would correspond to ( ) U Jt T , which gives rise to no entanglement between the two wells when they are prepared in thermal states (even at different effective temperatures), as demonstrated in ref. 25. However, this does not imply that the dynamics of the two-well system is trivial. In fact, in general, quantum correlations (of a form weaker than entanglement) are set by ( ) U Jt when acting on thermal states with ≠ n n 1 2 . We will address the emergence of discord-like quantum correlations [26][27][28] and its relation to the inter-well exchange process in a later section.
We now move to solving the full dissipative dynamics governed by Eq. (2) for U = 0. The problem can be efficiently solved by using a suitable Gaussian ansatz: We first translate Eq. (2) into the phase space by deriving a differential equation for the symmetrically ordered characteristic function after a lengthy but otherwise straightforward calculation, we find that Eq. (2) takes the form of the Fokker-Planck equation ] and expressing the characteristic function in terms of the entries of the vector of quadrature variables, we can write χ σ y z y z 1 1 2 2 and matrix σ  whose elements we aim at finding, which we do by solving Eq. (14). In Methods we provide the set of differential equations for the elements of X and σ  obtained when evaluating both sides of Eq. (14) and equating them term by term. The explicit solution of the problem at hand leads to a time-evolved covariance matrix of the general block form where c is a 2 × 2 matrix of correlations among the quadrature operators of the system. The diagonal structure of the blocks pertaining to the the individual wells shows that, locally, the system thermalises at temperatures determined by the explicit form of m 1,2 . However, as c is, in general, not null, global thermalisation is not achieved: the overall system never thermalises, notwithstanding an explicitly dissipative evolution. This is clearly seen by looking at the general form of the steady state. Although the analytic form of the non-zero elements is readily achievable for any value of the parameters involved in the problem, they are, in general, too cumbersome to be reported here. However, assuming γ 1 = γ 2 = γ, the steady-state of the system is determined by the covariance matrix Clearly, only for = n n 1 2 the structure of the global covariance matrix takes a thermal-like form. However, this does not preclude the possibility to achieve accidental thermalisation, i.e. situations such that the state of the system either becomes globally/locally thermal, or closely approximates an equilibrium configuration. This will be the focus of the following analysis.
Assessment of dynamical thermalisation. We shall start with the study of the unitary case. Dynamical thermalisation in closed-system dynamics is a topic of vast interest, which has recently attracted considerable attention at both the theoretical and experimental level 29 . Our approach is based on the assessment of the distance between the time-dependent state of the system and a generic (either global or local) thermal state [30][31][32][33] . Quantitatively, as a measure of the distance between two states ρ A,B , we use the Ulhmann fidelity 34 i y y is the two-mode symplectic matrix (σ y being the y-Pauli matrix) and , which is the one associated with the tensor product of locally thermal states (each with mean number of excitations μ j ). For clarity, we have indicated with n the identity matrix of dimension n.
We present the case of global thermalisation first: after calculating the time behaviour of σ σ ( ( ), ) F t u G 2 for various choices of Δ , we have numerically evaluated the value of μ that achieves the maximum of σ σ ( ( ), ) F t u G 2 . In Fig. 1 we show both such value and the corresponding estimate for μ. The state fidelity remains evidently quite large, being only partially depleted by an increasing value of Δ (the dependence on such parameter is quite non-trivial, given that for Δ = 2.5, for instance, values very close to those associated with Δ = 0 can be achieved, at suitable times in the evolution). However, while at small values of Δ the target state changes very little with time, this is not the case for increasing bias: the value of μ corresponding to a non-zero Δ oscillates with a non-negligible amplitude as this parameter grows. In any case, perfect thermalisation is never achieved, a result that is strengthened by the analysis that we will report in the next Section.
The situation is somehow different when locally thermal target states are considered [cf. Fig. 2]: besides the expected times at which a full period of the evolution is achieved, it is possible to identify instants of time at which the state of the double-well system is indeed very close to a locally thermal state , which would suggest the occurrence of accidental dynamical thermalisation.
In the open-system dynamics case, a similar calculation allows us to evaluate the fidelity with both a globally thermal and locally thermal state as shown in Fig. 3, which studies the effects of both the energy bias [panel (a)] and a difference in the damping rates of the two wells [panel (b)]. Quite evidently, both effects spoil the state fidelity, which however achieves values that are either precisely 1 or very close to it. Indeed, focusing on the unbiased case with Δ = 0 and γ 1 = γ 2 , we know that the solution is given in the form of Eq. (15). Moreover, the off-diagonal block matrix c turns out to be anti-diagonal with entries equal in modulus but opposite in sign. Therefore, in order to determine if the system has accidentally thermalised, we need only to determine if, at some value of t, these entries are identically zero. After some manipulation, this condition reduces to the transcendental equation Interestingly, this 'accidental' thermalisation is independent of the temperature of either well and only concerned with the tunnelling strength and the damping rate. For the same parameters taken to obtain the red curve in Fig. 3(b), we find Eq. (20) has two solutions: ω ∼ . , clearly corresponding to the two instances of local thermalisation in Fig. 3(b). Furthermore, we find the thermal occupation numbers of the wells at the first instance of thermalisation are = . , thus suggesting that the two instants of accidental thermalisation correspond to an almost swap of the two local thermal states. Increasing J leads to more instances of accidental thermalisation occurring before the system equilibrates to its steady state.
Assessment of the non-classical nature of the state of the system. Values of fidelity so close to unity should not lead to misinterpretation of the actual nature of the state of the two-well system. In fact, any . assessment of fidelity should be accompanied by the study of problem-specific figures of merit able to provide a more fine-grained characterisation of the state at hand. For the sake of a study on thermalisation, a significant class of such quantifiers is embodied by measures of quantum correlations. In this respect, it is important here to assess the role, if any, various forms of quantum correlations play in the dynamics highlighted above. It is quickly confirmed that, as anticipated before, the system never becomes entangled. While this is expected in light of the nature of the interaction and initial state being considered, nothing prevents the settlement of weaker forms of quantum correlations, such as quantum discord (QD) 28 . QD is the difference between two classically equivalent definitions of mutual information when applied to a quantum system 26,27 . A non-zero degree of QD implies that, in a bipartite system composed of parties A and B, information can be gathered on system A by interrogating party B. For Gaussian states, QD is captured by the Gaussian quantum discord [36][37][38] , which entails that the interrogation of B only involves Gaussian measurements. For a generic covariance matrix  I  I  I  I  I  I I  I  I  I   I  I I  I I  I I   I  I I  I  I  I  I  I I  I I I  where  . In Fig. 4(a) we study QD against the energy bias for the case of the unitary solution Eq. (11). Intuitively we would expect that for Δ = 0, owing to the full symmetry enforced in the system, QD will be maximised. This is indeed the case, as it can be seen in Fig. 4(a). However, an interesting feature appears as we increase the bias. At ∆/ ∼ . J 2 5, = n 1 1 , and = n 5 2 , QD exhibits a plateau, which implies the existence of an 'optimal' value of the bias, dependent on the temperature difference, that helps amplify the non-classicality of the system. Further increase of Δ pushes the systems too far off resonance, and the coherence decays. In Fig. 4(b) we examine this phenomenon closer, for a fixed temperature difference and small biasing, Δ /J = 1 (red), we see the oscillatory behaviour changes and the first zero-point is lifted. At the optimal value of Δ (solid black) the plateau is clearly evident. When we increase the bias further, we see the decay in the non-classicality, as well as a change in the periodicity of the system.
Turning our attention to the dissipative case, in Fig. 5 we compare the unitary dynamics with the dissipative case for unbiased wells [panel (a)], and biased ones [panel (b)] at various differences in temperature. For unbiased wells, we see dissipation quickly suppresses the the oscillations and we reach a steady state with non-vanishing QD. As we increase the temperature between the wells we see an increase in the QD for the unitary case and the steady state QD is larger for increasing temperature difference. When we bias the wells, taking Δ /J = 2.5 for all temperature differences, we see the dissipative dynamics clearly show the enhanced non-classicality. While the time to reach equilibrium appears unaffected, the steady state is significantly more non-classical than in the unbiased situation. This may imply that in this situation the non-classicality plays no role in reaching equilibrium. In Fig. 5(c) we examine the effect that self-interaction has on the dynamics of nonclassicality. In order to do so, we compute the Gaussian discord of the hypothetical Gaussian state having, as covariance matrix, the one achieved by calculating the entries σ ij over the non-Gaussian state resulting from a chosen non-zero values of U. Evidently the larger the self interaction, the more self ordered each well becomes, diminishing the effect of the tunnelling and reducing the amount of nonclassicality present. Also we see the system tends to equilibrate faster.
The nonclassicality of the steady state is delicately dependent on the temperature difference, as well as the tunneling rate and the bias. In Fig. 6 we examine this behaviour closer, fixing n 1 and γ ω γ ω / = / = 1 1 1 2 1 with J = 2 and Δ = 5. The only conditions for which the system does not exhibit nonclassical correlations is the trivial one of = n n 2 1 . As we increase the temperature imbalance we see that QD increases to a maximum value before slowly  decaying [cf. Fig. 6(a,e)]. If we fix the temperature difference such that = n 1 1 and = n 5 2 , we see in panels (b) and (c) that there are optimal values of the remaining parameters that give the largest value of discord. While the reservoirs have been kept at moderately low energies, in panel (f) we significantly increase both n 1 and n 2 and see that large values of QD can still be achieved. An unbiased configuration leads to values of discord of the order of 10 -5 . Increasing the bias, such values are raised by up to one order of magnitude.
Dynamics of the energy flux between the wells. It is important to gather insight into the details of the exchange of energy between the wells of the system, which is at the basis of the process of quasi-thermalisation highlighted so far and takes place in two forms: an exchange of particles between the wells and a similar process occurring at the interface between the double-well system and the reservoirs. The aim of this section is to identify the contribution coming from both such fluxes. We are thus interested in quantifying the flux into/from well j = 1, 2, which we label as  j  , and the total flux   tot . These are given by the quantities is the free evolution of a single well and  j is the density matrix of well j. Conveniently, these quantities can be directly evaluated from the covariance matrix (and we will assume both wells to have the same damping rate, i.e. γ ω γ ω γ / = / = 2 and clearly taking γ = 0 or Δ = 0 recovers the unitary and unbiased limits respectively. However, this behaviour is only for the special case of the wells initially being thermalised with their baths, while taking a different initial state this behaviour no longer holds. Indeed, what is special about our initial state is that it conserves the total energy of the system. This is readily seen given that =  0 tot  and it is easy to confirm that  . (e) Steady-state discord studied against both n 1 and n 2 when J = 2, Δ = 5, and γ ω γ ω / = / = 1  The tunnelling term in Eq. (1) commutes with  i , and when U = 0 the only contribution to the total flux is from the free evolution of each well. Therefore we are interested in calculating In a tedious but otherwise straightforward calculation, we can explicitly evaluate this expression when assuming H . We find that both the terms entering Eq.
(28) are identically zero, thus showing that, in the U = 0 case, the net heat flux is null due to two special circumstances: on one hand our chosen initial state, on the other hand the tunnelling term commutes with the super operators. In Fig. 7(a) we show the dynamics of the various energy fluxes given in Eq. (24). We notice how the flux into the cooler well is proportional to the flux out of the hotter well, which results in a null net flux. Needless to say, the single-well fluxes only account for the net intake/outtake of particles for one of the wells and do not provide information on the actual balance between the contribution due to the coupling to the reservoir and that due to the coherent inter-well interaction.
We can study the intermediate dynamical regime where the self-interaction is non-zero and comparable with the tunnelling by numerically solving Eq. (2) and examining the behaviour of the heat fluxes, of which we illustrate some examples in Fig. 7(b,c) [we refer to the caption for an account of the parameters used in the simulations]. The total flux is now non-zero, and the energy is not conserved. However, the average occupation number is conserved, i.e. ω ω / = − /ˆˆ † † a a a a 1 1 1 22 2 , which follows directly from the previous arguments.

Discussion
The analysis above shows that neither global nor local thermalisation with the reservoirs is achieved. The fidelity between the density matrices of the time-evolved state and the target thermal one (whether globally or locally) connects the closeness of the populations of the energy levels of the former to the statistics of the latter. However, the interaction between the wells establishes strong quantum coherence between the particles of the systems, which in turn results in the generation of a substantive degree of quantum correlations, albeit of a nature weaker than entanglement, which prevents the thermal character of the resulting state. The analysis reported here also has the merit of providing rather deep insight into the phenomenology of quantum correlations between the wells. We have qualitatively and quantitatively examined the dynamics and steady state of a BEC loaded into a double well potential. While the wells remain separable at all times, thus sharing no entanglement, by exploring the behaviour of the quantum discord we find the system to be always non-classical, except under trivial, uninteresting conditions. Furthermore, the degree of nonclassicality of the system is dependent on the energy bias between the two wells. For identical wells, a significant amount of QD is possible, provided that a large temperature imbalance is established. Such nonclassicality can be greatly enhanced by taking a suitable value of tunnelling, which must be a function of the given bias. The transfer of heat in the system is equally complex.

Differential Equations.
Here we provide the complete set of differential equations that describe the dissipative dynamics considered throughout.