Sudden change of the thermal contact between two quantum systems

In this paper, we address the issue of the stability of the thermal equilibrium of large quantum systems with respect to variations of the thermal contact between them. We study the Schr\"odinger time evolution of a free bosonic field in two coupled one-dimensional cavities after a sudden change of the contact between the cavities. Though the coupling we consider is thermodynamically small, modifying it has a considerable impact on the two-point correlation functions of the system. We find that they do not return to equilibrium but essentially oscillate with a period proportional to the length of the cavities. We compare this coupled cavities system with the perfect gas which is described by similar expressions but behaves very differently.


Introduction
The physically relevant degrees of freedom of a large isolated system initially out of equilibrium, are expected to relax to their thermal equilibrium values. Obviously, no relaxation behavior can be observed for the complete quantum state of the system but interesting degrees of freedom can evolve irreversibly in the limit of a large system. This issue of thermalisation in isolated quantum systems has been principally investigated by considering composite systems consisting of two distinguishable subsystems of very different sizes. For such systems, the smaller subsystem relaxes to thermal equilibrium if the larger one is assumed to be in an equilibrium mixed state [1]. Recently, it has been shown that a priori thermal averaging is not essential. The small subsystem can thermalise whereas the whole isolated system is in a pure quantum state [2].
An isolated system of identical particles cannot be divided into distinguishable subsystems and the relevant degrees of freedom are, in this case, the n-particle reduced distribution functions such as the particle number density. Relaxation behaviors have been obtained for both integrable [3,4,5,6] and nonintegrable systems [7,8,9]. In these studies, the system is supposed to be initially in the ground state of a given Hamiltonian but its subsequent time evolution is governed by a different Hamiltonian. For example, a gas is considered to be initially confined in a subvolume of a larger accessible space. The single-particle distributions were always found to relax to asymptotic profiles but which are, intriguingly, not all well described by an equilibrium ensemble.
To better apprehend the thermalisation process in isolated many-body systems, it is necessary not to be limited to ground states and to be able to treat also finite initial temperatures. As we are concerned with isolated systems, such an initial condition must be represented by a microcanonical mixed state, according to the postulate of equal a priori probabilities, or, equivalently, by a pure state of macroscopically welldefined energy [2,10,11]. In Ref. [11], the time-dependent particle density of a perfect quantum gas, initially in such a state, has been determined. It was found that this density relaxes from a thermal profile to a non-thermal one. In view of these results, one may wonder how two quantum systems, initially independent of each other but at the same temperature, evolve if they are brought into thermal contact. Do their relevant degrees of freedom remain at equilibrium once they are coupled ?
In this paper, we address this issue by considering a free bosonic field in two onedimensional cavities coupled to each other by one of their ends. Free fields are commonly used to describe the large number of environmental degrees of freedom which give rise to dissipation in a smaller system [1,12]. Here, we are interested in two systems of comparable size, neither of them acts as a heat bath for the other one. In section 2, we present our model Hamiltonian of two quantum systems in thermal contact. This model allows to describe from perfectly coupled to completely uncoupled systems. We show that the considered coupling between the two cavities is thermodynamically small by evaluating the equilibrium total energy and correlation functions. In section 3, we determine how the energy contents of the cavities and the field-field correlations evolve in response to a sudden change of the coupling between the two cavities. Our model is simple enough to obtain exact results for arbitrary values of the coupling strength.
In section 4, we discuss why, though the expressions are formally similar, relaxation behaviors are obtained for a perfect gas but not for the system studied here. Finally, in the last section, we summarize our results and mention a few open questions raised by our work.

Derivation of the Hamiltonian
To describe two coupled one-dimensional bosonic systems characterised by length L and velocity c, we start from the Hamiltonian where the fields Π and φ are canonically conjugate to each other. The position-dependent refractive index n(x) is equal to n > 1 for −d < x < d and to 1 elsewhere. Diverse physical systems can be represented by this Hamiltonian. The fields Π and ∂ x φ can be interpreted as the electric and magnetic components of the electromagnetic field, or as the charge and current distributions of an LC transmission line [1,13]. We write these fields as where the operators P and a q satisfy the commutation relations [P, a q ] = [a q , a q ′ ] = 0 and [a q , a † q ′ ] = δ qq ′ , and Λ is a cut-off length. Throughout this paper, we use units in which = k B = 1. The eigenmodes ϕ q of H ′ are solutions of the differential equation The even eigenmodes are proportional to cos[q(|x| − L)] for |x| > d, with q given by tan[q(L − d)] = −n tan(2nqd).
To obtain the Hamiltonian of two cavities coupled to each other by one of their ends, we consider the limits n ≫ 1 and d ≪ L with the length kept fixed. In this regime, equation (4) simplifies to tan(qL) = 2/Dq which can be rewritten as This second form will be useful in the following. For the even eigenmodes, the wavenumbers become q = pπ/L where p is a positive integer. It can then be shown, using (2) and (3), that the Hamiltonian (1) turns into where H right is given by this expression with −L and 0 replaced by 0 and L, respectively, and In the following, we study the Hamiltonian (7) which allows to describe from perfectly coupled to completely uncoupled cavities. Its eigenmodes ϕ q obey ∂ /D which leads to even ϕ q (x) = L −1/2 cos(qx) independent of the coupling characteristic length D. The odd eigenmodes read where A q = [L + (D/2)(1 + (Dq/2) 2 ) −1 ] −1/2 and q is a positive solution of (6). In the large D limit, the two cavities decouple from each other as can be seen from the fact that the boundary conditions at the ends 0 − and 0 + become similar to that at the ends −L and L. In this limit, the solutions q of the equation (6) are the multiples of π/L and the corresponding eigenmodes ϕ q simplify to L −1/2 sgn(x) cos(qx) for q > 0 and to (2L) −1/2 sgn(x) for the solution q → 0. In the opposite limit, D = 0, the odd eigenmodes are given by ϕ q (x) = L −1/2 sin(qx) with the wavenumbers q = (p+1/2)π/L where p ∈ N. The transmission at x = 0 is perfect in this case.

Thermal contact
The coupling (9) is thermodynamically small, i.e., its contributions to the thermodynamic functions of the complete system described by (7) (2) and (3) gives As there is only one solution to equation (6) in a given interval [pπ/L, (p + 1)π/L] where p ∈ N, the above sum over q can be approximated by an integral, which leads to The first two extensive terms stem from the integral approximation, see Appendix A.
The energy E is the difference between this approximation and the sum (11) and is thus expected to be negligible in the limit L ≫ c/T . For large and vanishing D, the wavenumbers q are equally spaced and it is then possible to evaluate E, see Appendix A. We find Note that the energy H T for large D is not exactly twice the energy of an isolated cavity of length L, the coupling between the two cavities contributes an energy T /2. In the limit L ≫ c/T , the energy (13) reaches T /2 and (14) vanishes. For an arbitrary characteristic coupling length D, E/T converges also to a finite value in this limit, as shown by Fig. 1, and hence the energy of the complete system H T ≃ Lc/πΛ 2 + πLT 2 /3c is essentially independent of D. We remark that, since the solutions of (6) are q ≃ pπ/L + 2/πpD where p is a positive integer, for large D, E first increases with decreasing D. The expressions (13) and (14) are valid for any LT /c. We see that, for L ≪ c/T , the contribution of the coupling (9) cannot be neglected as E (∞, T L/c) and E (0, T L/c) diverge differently in this limit. We obtain this diverging behavior of the energy E for other values of D as well, see Fig. 1.
1L and 0.2L, and D = 0.1L, L and 10L. The dotted line corresponds to the thermodynamic limit expression (18). For the above parameters, the agreement with this approximation is excellent for x ′ > x.

Equilibrium correlation functions
It is also instructive to consider the thermal two-point correlation functions of the system, for example where p runs over all integers and the second sum over the negative and positive solutions of (6). For vanishing and large D, this expression can be simplified by using the relation (A.1) and the function (A.3). We find, for D = 0 and D → ∞, respectively, where p runs over Z and ǫ over {−1, 1}. For L ≫ c/T , the only terms which contribute significantly to the sum (16) are (ǫ, p) = (−1, 0), (1, −1) and (1, 0). The last two are important only for (x, x ′ ) in the vicinity of (−L, −L) and (L, L). The vanishing of (17) for xx ′ < 0 shows that there is no correlation between the two cavities in the large D limit. The comparison of (16) and (17) reveals that, in the large D case, the thermal correlation function (15) is identical to that of two isolated cavities of length L. Consequently, in these two limiting cases, for x and x ′ not too close to the cavities ends. More precisely, Π(x)Π(x ′ ) T deviates notably from this approximate expression if x or x ′ is at a distance smaller than c/T from −L, 0 or L. For other values of D, the field-field correlations (15) are also well described, in the large L limit, by (18). This can be seen as follows. For x ≃ x ′ not too close to the ends −L, 0 and L, the terms exp[iq(x + x ′ )] and exp[iq(|x| + |x ′ |)] in (15) vary a lot from one value of q to the next and hence give vanishing contributions. Moreover, the normalisation factor A q is close to L −1/2 except possibly for the first few roots of (6), and this equation has only one solution in a given interval [pπ/L, (p+1)π/L] where p ∈ N. As a result, the two remaining sums in (15) are accurately approximated by the same integral which leads to (18). The thermal correlation function (15) differs significantly from this approximation only for x and x ′ close to the cavities ends as shown in Fig. 2.
For the field c∂ x φ, we obtain the expression where p runs over Z and the second sum over the solutions of (6). For D = 0 and D → ∞, this correlation function is given by (16) and (17), respectively, with an extra factor −ǫ in the summand. For other values of D, the above arguments can be used to simplify (19). Thus, the approximate expression (18) applies to the field c∂ x φ as well.
For the correlations between the fields Π and ∂ x φ, the completeness of the basis {ϕ q } implies, for any D, In summary, in the thermodynamic limit, the thermal correlation functions are essentially independent of the characteristic coupling length D.

Sudden change of the coupling strength
We consider here that the system is initially at equilibrium with temperature T and characteristic coupling length D 0 and then evolves under the Hamiltonian (7) with D = D 0 . We study the time-dependent two-point correlation functions and the time evolution of the energy contained in one cavity. Since, for these expectation values, a microcanonical mixed state or a pure state of macroscopically well-defined energy are, in the thermodynamic limit, equivalent to a canonical ensemble [2], we evaluate canonical averages in the following.

Time-evolved field operators
To obtain the field-field correlations at any time t, we write the time-evolved field operators in terms of the creation operators a † k corresponding to D 0 as The coefficients Θ k and Ω k are given by are the eigenmodes corresponding to D and D 0 , respectively. If ϕ (0) k is an even function of x then only the even ϕ q contribute to the sums (22) and (23). Moreover, in this case, since the even eigenmodes do not depend on the characteristic coupling length, these expressions simplify to Θ k = −(k/L) 1/2 sin(kx) exp(ickt) and Ω k = (k/L) 1/2 cos(kx) exp(ickt). For odd ϕ (0) k , it is instructive to first consider the two limiting cases of vanishing and large D.
3.1.1. Perfectly coupled cavities For D = 0, the odd eigenmodes are ϕ q (x) = L −1/2 sin(qx) with the wavenumbers q = (p + 1/2)π/L where p ∈ N. By evaluating (k|q) and using the relation (A.1), we obtain where p runs over all integers, τ = ct/2L,x = x/2L, ⌊τ ⌋ denotes the largest integer smaller than τ , and ǫ(τ ) = τ − ⌊τ ⌋. Inspection of this expression shows that the field Π is the superposition of two waves travelling with velocity c and reflected with no phase shift at x = ±L. For a given x, |Ω k | takes two different values as time goes on, see Fig  3. We find where . . do not depend on the wavenumber k. As we will see below, this independence plays an essential role in the behavior of physical properties such as the field-field correlations or the cavities energies. For the field ∂ x φ, the coefficient (22) is given by a similar expression.

Uncoupled cavities
For large D, we obtain Here, the field Π is the superposition of two waves reflected at x = 0 and L (−L) in the right (left) cavity. The expression (25) becomes where p = 0 in A/B, 1 in A ′ /B ′ , . . . . Contrary to the above case, the field Π is always discontinuous at x = 0. The field ∂ x φ satisfies ∂ x φ(0, t) = 0, i.e., the reflection is perfect at x = 0 + and 0 − .

General case
For a finite value of D, the fields Π and ∂ x φ are also superpositions of propagating waves but which are, contrary to the above studied special cases, both transmitted and reflected at x = 0. We write the coefficients Ω k and Θ k as where the sum runs over the positive and negative solutions of (6). To evaluate the functions Γ ± k , we first note that, since q is a solution of (6), the summand in (30) can be multiplied by a factor [(Dq + 2i) exp(−2iqL)/(Dq − 2i)] m where m is any integer. We then rewrite Γ ± k as where X(z) = (z + 2i) exp(−2izL/D)/(z − 2i) and C is a contour enclosing the real axis.
Expressing the above integral in terms of the residues of the integrand poles on the real axis leads to (30). For times t between t ± m and t ± m+1 where t ± m = ∓|x|/c + 2mL/c, the upper (lower) part of the contour C can be closed in the upper (lower) half plane. Moreover, for positive m, the integrand in (31) has poles only on the real axis and at z = 2i. Thus, for t > 0, the first term of the expression (31) is determined by the residue of this last pole if m is chosen as the largest integer smaller than (ct ± |x|)/2L. It can be shown that the functional dependence of this term on τ = c(t−t ± m )/D is exp(−2τ )P ± m (τ ) where P + m (P − m ) is a polynomial of order m − 1 (m), see Appendix B. The functions Ω k and Θ k change discontinuously at t = t ± m . For D ≪ L, they then reach the values given by the second term of (31) in a time of order D/c. For D = 0, one retrieves the Dirac delta functions of (24) and the second term of (31) gives the expression (25). In the opposite limit, D ≫ L, the first term of (31) is negligible with respect to the second one which leads to the expression (27). Contrary to the limiting cases of vanishing and large D studied above, since the successive reflections and transmissions at x = 0 deform the propagating waves, |Ω k | and |Θ k | are not strictly periodic here.

Correlation functions
With the notations introduced in (20) and (21), the two-point correlation functions of the cavities read Since the functions Θ k and Ω k present discontinuities at positions independent of k, as discussed above, the field-field correlations at a given time change abruptly at x = this term and using the properties X(Dk) * = X(Dk) −1 = X(−Dk), we obtain where the second sum runs over the positive and negative solutions of (6) with the characteristic coupling length D 0 . The coefficients α k and β k are given by where X = X(Dk) and the regions A, B, A ′ , . . . are presented in Fig. 3. Therefore, the approximate correlation function (35), for given x and x ′ , is periodic with period 2L/c. When x and x ′ are both in a region of type A, the expression (35) is identical to the initial condition (15). As discussed after equation (18), the second and fourth terms of (35) are negligible except for x and x ′ near a cavity end, and hence (35) is practically equal to (15) for x and x ′ both in a region of type B, see Fig. 4. When x and x ′ are in regions of different types, (35) can differ considerably from (15). For example, for cavities perfectly coupled (D = 0) but initially uncoupled (D 0 → ∞), α k = β k = 0 in this case and thus the spatial correlations of the field Π are reduced by a factor of two with respect to the initial condition (18). It is interesting to note that even for uncoupled cavities (D → ∞), the sudden change of the characteristic coupling length  For the other correlation functions, we find where α k and β k are given by (36), and, for x = x ′ , where the sum runs over the solutions of (6), the upper sign is for |x| > |x ′ |, and Contrary to the correlations (35) and (37) of a field with itself, the correlations between the fields Π and ∂ x φ given by (38) vanish when x and x ′ are in the same region, see Fig. 4. For x and x ′ not too close to a cavity end, 0) )/c where the sign depends on the relative positions of 0, x and x ′ . We have seen in the previous section that there is no correlation between the two fields and that the correlations of a field with itself are essentially the same everywhere at equilibrium. Switching the characteristic coupling length to another value induces correlations between the two fields and affects the whole system even far away from the interface between the two cavities.

Energy of the cavities
We study here the time evolution of the energy content of the cavities. For that purpose, we first note that the energy density is even with respect to x. This symmetry property leads to The time dependence of the cavities energies is thus simply related to that of the coupling energy. As the eigenmodes of the total Hamiltonian H satisfy /D, the time-evolved coupling Hamiltonian (9) is proportional to ∂ x φ(0, t) and hence can be expressed in terms of the functions (22). We obtain The even eigenmodes do not contribute to this sum since their derivatives vanish at x = 0. For the odd eigenmodes, we have seen above that, for D ≪ L, the first term of the expression (31) contributes only for t close to 2mL/c where m is an integer. The second term of (31) gives, using |X(Dk)| = 1, where the sum runs over the wavenumbers k corresponding to odd ϕ (0) k . For a small but finite D, the coupling energy deviates from this value only for short periods of time when t is much smaller than a characteristic time t 0 which increases with decreasing D. In this short time regime, the energy density (40) is essentially constant except in two regions wich propagate with velocity c and are reflected and transmitted at x = 0, L and −L. The multiple peak structure shown in left part of Fig. 6 results from the preceding reflections and transmissions at x = 0. The temporal evolution of H coupling (t) becomes more complicated for t > t 0 , see Fig. 6. For D ≫ L, the contributions of the first and second terms of (31) to the functions Θ k (0, t) are of the same order and hence the approximate expression (43) does not apply. For D → ∞, ∂ x φ(0, t) vanishes, as mentioned after (27). However, because of the factor D in (42), some care must be taken in the evaluation of H coupling (t) in this limit. For large D, the wavenumbers q are practically the multiples of π/L except the lowest one q ≃ (2/DL) 1/2 . The terms of the sum (22) are thus essentially equal to their infinite D values except the first one and varies sinusoidally with a finite amplitude.

Comparison to perfect gas
The time-dependent average energies H left (t) , H right (t) and H coupling (t) , and the correlation functions (32)-(34) can be written as whereÔ stands for H left , . . . , the sums run over the positive and negative q corresponding to the characteristic coupling length D, O qq ′ are appropriate coefficients, ǫ q = cq and O(ω) = q,q ′ >0 O qq ′ δ(ω − ǫ q + ǫ q ′ ). We remark that the even eigenmodes ϕ q contribute only to the steady component of Ô (t) , see the discussion after equation (23). The expectation values of the single-particle observables of an isolated perfect gas, its particle number density at a given point for instance, are also given by expressions of the form (45) with the single-particle energies ǫ q [10]. Consider as an example a one-dimensional perfect gas confined in a box of length L. In this case, the single-particle energies are ǫ q = p 2 π 2 /2mL 2 where m is the particle mass and p an integer. The frequencies ǫ q − ǫ q ′ are thus regularly spaced and Ô (t) is periodic with period 4mL 2 /π. However, in the Joule expansion studied in Ref. [10], for times t ≪ mL 2 , the function O in the integral expression (45) of the gas density profile, can be approximated as a continuous function plus a term O ∞ δ(ω) and Ô (t) essentially relaxes from its initial value to O ∞ . The situation is radically different for the free field system considered in this paper. There obviously also exists a δ(ω) contribution to O(ω) but Ô (t) does not relax. First of all, we observe that, in the limiting cases of vanishing and large D, the function O is a sum of equally spaced Dirac delta functions and Ô (t) is strictly periodic with period 2L/c as shown by (25) and (27). For a finite D, the frequencies ω qq ′ = ǫ q − ǫ q ′ = c(q − q ′ ) are not exactly equal to multiples of πc/L but most of them are close to these values, see Fig. 7. As it is clear from equation (6), the spacing between consecutive frequencies ω qq ′ varies and goes to zero for ω qq ′ → pπc/L where p ∈ Z. We see on Fig. 7 that O qq ′ → 0 in this limit. Time regimes may then exist in which O(ω) can be approximated as a sum of smooth functions with finite support. However, these functions clearly vanish at the right end of their support interval but not at the left one. Consequently, their Fourier transforms do not vanish fast enough at long times to lead to a relaxation behavior of Ô (t) . We find similar functions O for the energy H coupling (t) with the noticeable difference that the coefficients O qq ′ corresponding to the lowest frequencies ω qq ′ ≃ ±2c(2/DL) 1/2 are dominant in the large D limit, giving rise to the sinusoidal behavior discussed at the end of the previous section.

Conclusion
In this paper, we have studied how changing the thermal contact between two quantum systems whose elementary excitations are noninteracting bosons, affects these systems. More precisely, our model consists of two one-dimensional cavities of equal length coupled to each other by one of their ends. The coupling we have considered is thermodynamically small, as shown by the fact that the equilibrium total energy and correlation functions are practically independent of its strength, and allows to describe from perfectly coupled to completely uncoupled systems. We have seen that a sudden change of this thermal contact has a considerable impact on the correlation functions of the two cavities. They do not return to equilibrium but essentially oscillate with a period equal to twice the time required for a signal to propagate from one end of a cavity to the other. Their exact time evolution is mainly controlled by the final coupling strength which determines the reflection and transmission coefficients at the interface between the two systems. We found a similar behavior for the energy content of the cavities. For weak coupling between them, it varies sinusoidally with a period set by the coupling strength. The absence of relaxation behavior is not simply a consequence of the fact that the considered Hamiltonian is quadratic. The Joule expansion of a perfect quantum gas has, for example, been obtained in Ref. [10]. As discussed in the previous section, this difference in behavior comes from qualitative differences in the elementary excitation energy spectrum. Clearly, the simple dispersion relation of our one-dimensional continuous field model, plays a crucial role in the temporal quasiperiodicity obtained. It would be thus interesting to study higher-dimensional and lattice free field systems. Are interactions between the bosonic excitations necessary to ensure the stability of the thermal equilibrium ?
The above integral is determined by the residue of the function F at z = 2i. We write With these notations, (B.1) becomes for m > 1 (P + 0 = 0). We find, for example,