Thermalization, prethermalization and impact of the temperature in the quench dynamics of two $unequal$ Luttinger liquids

We study the effect of a quantum quench between two tunnel coupled Tomonaga-Luttinger liquids (TLLs) with different speed of sound and interaction parameter. The quench dynamics is induced by switching off the tunnelling and letting the two systems evolve independently. We fully diagonalize the problem within a quadratic approximation for the initial tunnelling. Both the case of zero and finite temperature in the initial state are considered. We focus on correlation functions associated with the antisymmetric and symmetric combinations of the two TLLs (relevant for interference measurements), which turn out to be coupled due to the asymmetry in the two systems' Hamiltonians. The presence of different speeds of sound leads to multiple lightcones separating different decaying regimes. In particular, in the large time limit, we are able to identify a prethermal regime where the two-point correlation functions of symmetric and antisymmetric sector can be characterized by two emerging effective temperatures, eventually drifting towards a thermal regime, where these two correlators become time independent and are characterized by a unique effective temperature. If the initial state is at equilibrium at non-zero temperature $T_0$, all the effective temperatures acquire a linear correction in $T_0$, leading to faster decay of the correlation functions. Such effects can play a crucial role for the correct description of currently running cold atoms experiments.

In the theoretical description, it is often assumed that the two TLLs are identical, meaning they are characterized by equal sound velocities and Luttinger parameters. In this case, the theory consists of a quantum sine-Gordon model and a free boson [50,51], describing respectively the antisymmetric and symmetric combinations of the phase fields (see section II for proper definitions). Importantly, as a consequence of the symmetry between the two TLLs, these two sectors are not coupled g > 0 FIG. 1. Schematic picture of the system studied in this paper. It consists of two unequal Luttinger liquids, with sound velocity ui and Luttinger parameters Ki (i = 1, 2). We quench the system by switching off the tunnelling g, starting from a non zero value (both ground state and finite temperature thermal states are considered as initial states). This corresponds to suddenly raising the barrier of the double well potential separating the two sides. and thus can be treated as isolated systems. In particular, time-dependent correlation functions of the antisymmetric sector (directly related to interference measurements [52]) after a sudden change in the tunnelling strength (a so-called quantum quench [13]) have been widely studied [53]. They have been obtained by relying, for example, on a simple harmonic approximation [54][55][56][57] and, more recently, on a refined selfconsistent version of it [58,59]. Exact results have been further obtained at the Luther-Emery point [55] and by means of techniques such as integrability [51,60,61] and semi-classical methods [62]. A truncated conformal approach was considered in [63], while a combination of analytic (based on Keldish formalism [64]) and numerical methods was used in [65]. Finally, an effective model for the relative degrees of freedom was recently derived in [66]. In these studies the existence of a prethermal regime was demonstrated.
Much less attention has been devoted so far to the effect of introducing an "imbalance" between the two systems. On the theory side such a case is interesting since, due to the presence of two velocities, one can expect multiple lightcones to emerge, separating different decaying regimes (as opposed to the single lightcone effect [13,29] usually observed in systems of identical TLLs [31,57]). Because of the coupling between the modes one can also expect that the prethermal regime evidenced in the antisymmetric sector could now decay into another thermal regime. Whether such a regime could be characterized by a single temperature despite the integrable nature of the underlying model [17,18] is an interesting question. However due to the complexity of such a situation the asymmetric has been much less studied. A noteworthy exception is provided by Ref. [67], where, relying on a phenomenological approach for the quench, the authors consider the quench dynamics of two bosonic tubes characterized by different densities. This studies showed the existence of two light cones and the potential decay of the prethermal regime into a long time thermal one.
Given the importance of the physical effects in the asymmetric situation, it would thus be highly desirable to have: i) a full theoretical derivation of the quench of two different TLLs; ii) to allow for the general case of two different sound velocities u i and Luttinger parameters K i (i = 1, 2), to be able to disentangle the effects of the variation of the velocities from the ones of the Luttinger parameters and to allow for all possible sources for the imbalance. Such a study is the goal of the present paper.
The paper is organized as follows. In section II we introduce the model and the quench dynamics we focus on. Section III and section IV discuss the Bogoliubov transformation which diagonalize the hamiltonian at initial time and introduce the correlations functions of interest, respectively. In section V a detailed analysis of the dynamics when starting from the ground state (i.e., at zero temperature) of the initial hamiltonian is carried out. The same analysis is extended to quenches starting from a thermal states in section VI. A discussion of the results, also in connection with previous literature, is left to section VII. Conclusions and future perspectives are finally collected in section VIII. Details regarding the calculations are reported in the appendices.

II. Setting of the quench
We consider two different Luttinger liquids which are initially tunnel-coupled and then evolve independently: this is the simplest quench one can look at, since the evolution is the one of two free (compactified) bosons, while the coupling between the two is only in the initial state. This protocol has also the advantage to be easily implementable in a controlled way in cold atom experiments.
Microscopically, the system corresponds to two interacting 1d Bose gases, represented by bosonic fields Ψ i (i = 1, 2) of mass M i and effective 1d interaction U i . We are going to work with their phase θ i (x) and the fluctuation of the densities n i (x), related to the original field via the bosonization formula [3,48,49]: and ρ i is the average density of the i-th tube. In terms of these variables, the system is supposed to be prepared in the ground state (or in a thermal state) of the (generalized) Sine-Gordon Hamiltonian: where H i are the Luttinger liquid Hamiltonians [49]: and the cosine term originates from the tunnelling (Ψ † 1 Ψ 2 + h.c.), with strength tuned by g. In (3) K i is the Luttinger liquid parameter which encodes the interaction of the system and u i is the speed of sound. They are related to the microscopic parameters via Therefore one can get unequal TLLs in many different settings, depending on the values of M i , U i and ρ i . Hereafter we will set = 1. At time t = 0 the interaction between the two systems is switched off and the final Hamiltonian simply reads: As the study of the initial hamiltonian (2) is particularly involved, we resort to a semiclassical (harmonic) approximation: Note that in our quench the approximation is only on the initial state, while the dynamics can be obtained exactly. Such approximation is expected to hold as long as the cosine term in (2) is highly relevant in a renormalization group (RG) sense (in the case of identical TLLs, this corresponds to K large enough [49], while the same RG analysis is missing for the more generic case considered here; note, however, that in the experiments involving bosons with contact interactions we can safely assume that we are in the relevant regime). Remarkably, for identical TLLs, it has been shown by means of exact calculations that the dynamics starting from (6) is qualitatively the same as from the Luther-Emery point where the full cosine term can be taken into account [55]. The fields θ i and n i admit a decomposition in normal modes: (8) In terms of these bosons the final Hamiltonian is diagonal, namely: where the zero modes (i.e., θ i,0 , n i,0 ) have been neglected. The Hamiltonian H SC initial , instead, is quadratic but needs to be diagonalized via a Bogoliubov transformation (see Section III below).
To highlight the difference with the case of two identical systems it is useful to introduce the symmetric (+) and antisymmetric (−) modes: which satisfy canonical commutation relations. In terms of these variables the final Hamiltonian reads: with u K = 1 2 Therefore we see that in the case of two identical systems the final hamiltonian is decoupled in the symmetric and in the antisymmetric sector and the quench occurs only in the antisymmetric variable. The situation that we consider in this work is more involved as this decoupling is not possible and to study correlation functions of θ − , which are those interesting for experiments, one has to consider the dynamics of θ 1 and θ 2 which are correlated via the initial condition.

III. Bogolioubov transformation for two species of bosons
In order to characterize the state that is evolving we aim to diagonalize the initial Hamiltonian H SC initial and write it as: up to an unimportant overall constant, that we neglect. This transformation amounts to a Bogoliubov rotation of a four component vector, mixing the modes (p, −p) of the two initial species of bosons. Specifically, we introduce the vectors of bosons of the initial and the fi- . These two are related by a matrix multiplication b p = B(φ p )η p with B(φ p ) depending on the set of parametersφ p = {ϕ 1,p , ϕ 2,p , ∆ p , φ p } and parametrized as follows [68]: Details on the derivation are reported in Appendix A. The parameters of the matrix B(φ p ) in (15) have the following interpretation: ϕ 1,p and ϕ 2,p define Bogoliubov rotations associated to the two bosons, separately. φ p is the mixing angle between them. Finally, ∆ p exists only when the Bogoliubov rotation and the mixing of different bosons appear at the same time [68]. Explicitly, they are given by: in terms of the parameters (for i = 1, 2) and the eigenvalues of the hamiltonian (14): Note that, at the leading order in p → 0, the eigenvalues (19) read: with: in terms of the parameters defined in (12). Therefore, they describe a massive and a massless mode. Note also that, in the limit of equal TLLs, they would coincide with the antisymmetric and symmetric modes, respectively.

IV. Correlation functions after the quench
We will be mostly interested in the correlation functions: where the expectation value · T0 is on the initial state, that we choose to be either the ground state (T 0 = 0) or a finite temperature (T 0 = 0) equilibrium state of the initial hamiltonian H SC initial . Note that in our approach, due to the absence of decoupling between symmetric and antisymmetric variables, θ ± are not anymore the preferable variables to work with (as it was the case in the symmetric quench [56,57]). Instead, we will stick to the initial fields, θ 1 and θ 2 . In terms of those variables, the one (two) point function of the symmetric or antisymmetric fields is recast into a two (four) point function.
We start by defining the parametersˆ p = {u 1 |p|, u 2 |p|} entering in the definition of the time evolution operator and the matrices: where ⊗ denotes the Kronecker product. For a generic quench starting from a thermal state of (6) at temperature T 0 , Eq. (22) takes the compact form (see Appendix B for details): (25) For convenience we have introduced an ultraviolet cutoff α −1 . We further denoted by W ± µν the elements of the matrices Note that only two elements of the whole matrices are needed to fully characterize the correlation functions (25). Moreover, thanks to the quadratic approximation in the initial hamiltonian, those coefficients can be written explicitly (see Eq. (B8) in Appendix B). In order to define effective temperatures for C ± , we are going to compare these post-quench correlations with the equilibrium ones at finite temperature T : which present an exponential decay in space with (inverse) correlation length: V. Quench from the ground state We consider here the quench from the ground state (T 0 = 0) of the Hamiltonian (6) and we defer the solution of the dynamics from a thermal state at temperature T 0 to section VI. In this section, expectation values over the ground state will be simply denoted as · .

A. Eigenmodes dynamics
An important observation is that in the limit T 0 → 0, we have cotanh(λ i,p /(2T 0 )) → 1 in Eq. (25), and it turns out that the leading order as p → 0 of C ± is captured uniquely by the first term, namely by the massive mode. The main contribution is better characterized by introducing the dynamics of the modes of phase and density. In particular, by using the following decomposition one finds that: Ki|p| . The expectation value of the two point function over the ground state simplifies to: for i, j = 1, 2, namely the initial correlations between θ i and n j do not enter in the eigenmodes' dynamics.
By plugging the asymptotic expressions (20), one can further check, that at the leading order in p → 0 it holds: where we defined u1 . Note the initial anticorrelations between the densities of the two systems, which will have a role on the evolution of the phase.
The leading divergence as p → 0 in C ± (x, t) ≡ C ± (x, t, T 0 = 0) (cfr. Eq. (25)) comes from the initial density fluctuation while the part coming from the phase is negligible (this is due to the term α i,p α j,−p ∝ 1/p 2 in the eigenmodes' dynamics (31)). Notice that since the sound velocity a appears only in the phase fluctuations, at this order the massless mode will not play any role in the correlation functions, consistently with what anticipated from Eq. (25).
If we define the building block of the correlations (25) as: (34) such that then from (31) and (33) we have: If we neglect the cutoff, the integrals (36) can be analytically evaluated. They are of the form This shows explicitly the emergence of (sharp) lightcones, associated to each velocity u within the correlation functions. Note that the light-cones are smoothen out (as physically expected) by reintroducing the cutoff. Correlations like those in Eq. (36) appear in the exponent of C ± . Therefore, we expect the approximation (36) (whose integrand behaves as 1/p 2 ) to capture only their exponential decay. A careful analysis should take into account possible power law corrections which come from the next-to-leading order correction (corresponding to an integrand ∝ 1/p). These can be computed explicitly as follows and therefore grow unbounded at large distances. As we are going to discuss, these terms are actually important for the dynamics of C + : in fact, there are regimes where the exponential behavior vanishes and power laws become leading.
B. Two-point function: transient, prethermal and stationary state By looking at the Eq. (36), we can read the leading terms in the two point function (22), which presents a very reach behavior. Without loss of generality, we may assume u 1 > u 2 . Then, we find: We stress that the expression above only captures the exponential decay of C ± , while power-law corrections are not included. In particular, within this approximation Eq. (39a) shows no spatial dependence: this does not mean that it does not decay at all, but that the next to leading term should be taken into account. From a more refined analysis it turns out that Eqs. (39) are actually correct also at the next-to-leading order for C − , whereas the symmetric sector has logarithmic (power-law) corrections at all times, which are particularly important in the first large distance regime that actually reads: This result is easy to understand from a physical point of view since before the quench θ − is a massive mode that does not decay to zero while θ + is massless with leading power-law correlations. In the regime of large distances and short times therefore we find the memory of such initial condition.
Moreover, from (39a) and its refined version (40), making use of the cluster property, lim x→∞ O(x, t)O(0, t) = O(t) 2 , we can read the behavior of the one point function: For the antisymmetric sector we obtain the exponential decay: Note that for the symmetric quench (K 1 = K 2 = K and u 1 = u 2 = u) we obtain the results of [13,69,70] with the scaling dimension of θ − equal to h = 1/(4K). Contrarily, for the symmetric sector, due to the power-law correction in (40), we find a vanishing one-point function at all times, i.e., A + (t) = 0. From Eq. (39e) we see that C ± reaches a stationary state at short length scales (large times). Note also that the next to leading order terms in 1/p, such as integrals of the form (38) do not lead to important time corrections in the limit of large times and that formally one expect the limit of t → ∞ to be really time independent as all the oscillating factors die out. From (39e), one can therefore read off an associated correlation length: equal for both the symmetric and the antisymmetric mode, signaling that the correlations between the system one and two are lost. Comparing with (28) this defines an effective temperature: Then, if u 1 ≈ u 2 , we have u 1 + u 2 ≈ 2u 1 ≈ 2u 2 |u 1 − u 2 |, and from Eq. (39d), we see that one can define a quasi-stationary prethermal state with correlation length and thus effective temperature different for the symmetric and the antisymmetric mode: This is the regime to which the dynamics stops in the limit of u 1 = u 2 and thus in particular for the symmetric quench. One can indeed check that for the symmetric quench (u 1 = u 2 = u and K 1 = K 2 = K) we recover the results T eff − = m 0 /4 and T eff + = 0 as expected from [4,56] and from the decoupling of the modes.
In Fig. 2 we show the logarithm of the correlation functions C ± (x, t) after a quench from the ground state of the Hamiltonian (6) (T 0 = 0) as a function of distance x and time t. The exact expressions (left panels), numerically computed from Eq. (25), are compared with the small momenta approximation derived in Eq. (39) (right panels). The position of the lightcones are also shown. While C − is well approximated by its exponential decay only (according to (39)), for a correct description of C + power-law correction must be included (this becomes more and more visible when looking at smaller scales). Note that for the parameters chosen, we are in the regime u 1 ≈ u 2 . This is the reason why the regimes in (39b), (39c) and (39d), shown as dashed lines in the Figure, are not well separated. Finally, the dot-dashed line corresponds to the last lightcone at x ≈ |u 1 − u 2 |t, separating prethermal and final thermal regime.
We then focus on the spatial decay of C + (x, t) and C − (x, t), for different (fixed) times (this would corresponds to horizontal cuts in Fig. 2). In Fig. 3 we compare this decay with the equilibrium correlation functions at temperatures T eff and T eff + for C + (x, t) and T eff and T eff − for C − (x, t), which captures the first two exponentially decaying regimes. In fact, for both correlations the longest time (short distance) shows the crossover between the thermal and the prethermal regime at distances around |u 1 − u 2 |t. After this decay C + (x, t) is characterized by a non monotonic behavior in the intermediate regimes. At large distances, differently as compared with C − (x, t), it does not saturates but it slowly decreases due to the power law corrections. The shortest times (long distance) of C − (x, t) instead show a light-cone like behavior toward a constant value for large distances. For the choice of parameters in the figure, the first three light cones in (39) are not separately visible in C − (x, t) because they are very close (in C + (x, t) they correspond to the non monotonic behavior). The presence of the cut-off in (25) also tends to smear out the sharp transitions in (39), as anticipated.
In Fig. 4 we compare the spatial decay of C − (x, t) and C + (x, t) in the thermal and in prethermal regime, comparing also with the thermal correlations at temperature T eff , T eff − and T eff + . The plot shows that the correlation length of the two quantities coincides in the first (thermal) regime and is also compatible with the equilibrium decay at temperature T eff . From this analysis therefore, the last regime can be thought (at the leading order) as a thermal regime, at least for the observables considered here. However, as we discuss in section V C the more robust thermodynamic interpretation of the stationary state is in terms of two temperatures, one for the first system and one for the second, which combined give rise coherently to Eq. (44). At larger distances the two correlations C ± depart from the thermal regime and from each other, and agree with an equilibrium-like behavior at temperature T eff − and T eff + , respectively (here we explicitly see a dependence on the observable chosen).

C. Interpretation as a two temperature system
The final Hamiltonian (5) has clearly two extensive and different conserved quantities: the energy of each subsystem. Therefore we expect to be able to define an effective temperature associated to each of them from the expectation values of the energy densities of the modes i,p ≡ u i |p| b † i,p b i,p with i = 1, 2 separately. In the limit of small momenta the expectation of each mode is dominated by a constant term equal for all the modes. Thanks to a classical equipartition approximation this allows us to interpret such constant as the effective temperature of the two systems [57]. In particular we have: Note that: as for the symmetric quench [4,56]. However, contrary to the symmetric limit where all the energy is stored in the antisymmetric sector (being isolated from the symmetric one), in the general case part of it is shared with the symmetric mode as well. Note that if one supposes the two systems equilibrated at different temperatures, the correlation length associated to the decay of C eq ± turns out to be: thus generalizing the expression (28). This is perfectly consistent with the effective temperatures (46) and the post quench correlation length (43). Specifically, the unique effective temperature that we can read off from C ± at large times is related to the ones in (46) through The two-temperature interpretation is sustained also by an FDT (fluctuation-dissipation theorem [71][72][73][74]) argument, analyzing the correlation and the response functions associated to the Green's functions of the two systems -in the limit of small ω (see Appendix C). We close this section with an interesting remark. Knowing that in our quench (when u 1 = u 2 ) the correlations between the two systems are lost in the latest thermal regime, one can expect the thermal state to which the system evolves to coincide with the final state reached by the same system of two bosons but after two independent quenches with initial energies (or initial masses, equivalenty) fixed by (46). As main difference, in this simpler quench, correlations are absent also in the initial state. If fact, since each H i (i = 1, 2) describes a conformally invariant system, one can directly apply the results of [13,69] for the correlation functions to see that at largest times: The expectation value · m1,m2 in the first line is taken on a factorized state characterized by mass m i for the ith system, which, therefore, simply splits in expectation values over the two systems (second line). The results of [13,69] have been applied to each · mi . In the last step, we used the explicit form h i = 1/(4K i ) for conformal dimensions of the (primary) operators e ±iθi(x,t) [75]. We stress, however, that, even if the result (50) is consistent with the last regime with associated effective temperature (49), the transient and prethermal regime are not captured by this simple picture.
VI. Quench from a thermal state: corrections due to the initial temperature If the initial state is prepared at finite temperature T 0 , the full expression for the correlation function is still the one in Eq. (25). Now, however, one sees that, differently from the quench from the ground state, the leading contribution as p → 0 includes a term coming from the massless mode. One can in principle carry a similar analysis as the one of the section V (notice in particular that Eq. (31) remains true also when starting from a thermal state), leading to different regimes during the evolution. In particular in Appendix D we sketch the derivation of the leading order term contributing to C ± (x, t, T 0 ) showing that the same light cones as for T 0 = 0 appear, with different correlation lengths and coherent times. Here however, we focus on the last two regimes (at large times), being the most relevant for the relaxation dynamics. As before, indeed, they allow for a definition of a prethermal and a thermal correlation length, for both the symmetric and the antisymmetric mode. The associated prethermal effective temperatures now read: For the symmetric quench we recover T eff as in [76] and T eff + = T 0 , as expected from the decoupling of the modes. A crucial observation here is that in this symmetric limit, the antisymmetric sector is almost unaffected by the true temperature of the system: in fact T eff − m 0 /4, namely it is independent on T 0 (as long as it is low), while the thermal fluctuations are present only in the symmetric mode, as reflected by its effective temperature. The reason is that, while θ 1 and θ 2 are subject to thermal fluctuations, those cancel out in their difference (namely in θ − ), while remain present in their sum (i.e., in θ + ) [33]. Importantly, this picture completely changes as soon as an asymmetry is induced in the parameters u i , K i associated with the two tubes. In fact, Eq. (51) clearly shows a correction linear in T 0 for the effective temperature. To be more precise, for such linear correction to be present in the antisymmetric mode as well, different sound velocities, i.e., u 1 = u 2 , are needed (while a difference in the Luttinger parameters K i does not seem to play a main role here). In this case, the initial temperature plays a crucial role in the decay of all correlation functions. Specifically, since the term proportional to T 0 in (51) is always positive, it leads to a faster decay of C ± .
The final regime is instead described by: which also shows a term depending linearly on the initial temperature, leading to faster decaying correlations.
We mention that the limit of shallow quench m 0 → 0 (which amounts to do nothing to the system) does not reproduce the equilibrium result T eff = T 0 . This signals that the limit of small momenta p → 0 used in deriving T eff in (52) does not commute with the limit m 0 → 0.
In Fig. 5 we show how the effective temperatures derived above capture the main decay of the correlation functions, also in this thermal quench. In particular, it shows the spatial decay of C − (x, t) starting from a thermal state at temperature T 0 = 0.5 and different times. The correlation lengths in the thermal and in the prethermal regime are compared with the one at equilibrium at temperature T eff − from (51) and T eff from (52). The plateau attained at large distances (short times) is instead a property of the (massive) initial condition.
In Fig. 6 we plot the same correlator as a function of time. The top panel shows the time dependence of C − (x, t) again starting from a thermal state at temperature T 0 = 0.5 and different points in space. In the bottom panel, instead, the inverse correlation length at fixed distance is showed, as obtained from the spatial derivative of the exponent in Eq. (25). We see that in the regime u 1 , u 2 |u 1 − u 2 | there is an intermediate prethermal regime where the correlation length is compatible with the equilibrium one at temperature T eff − . At later times this quantity crosses over towards the asymptotic regime, compatible with the equilibrium one at temperature T eff .

VII. Discussions
Let us make some comments about the results that we have obtained in the previous sections, also in comparison with previous works.
To start with, as a generic dependence on the sound velocities u i and Luttinger parameters K i is considered in our analysis, it is worth stressing the different role that they play in the dynamics. In fact, while we saw that, in the Hamiltonian (5), which governs the evolution after the quench, a coupling between symmetric and antisymmetric mode is present as soon as the system one and two differ in either of those parameters (cfr. Eq. (11)), the consequences of having different K i or different u i separately are not the same. If u 1 = u 2 , then the correlation functions (39) are much simplified and only one lightcone appears, with the dynamics never reaching a thermal regime (39e). This means that the symmetric and the antisymmetric modes always show different effective temperatures. Moreover the linear correction of the effective temperature of the antisymmetric mode (51) due to the initial temperature T 0 vanishes. This suggests that a sort of decoupling between different sectors still exists. In fact, in the final Hamiltonian (5), one could rescale the field θ i by √ K i and n i by 1/ √ K i in such a way to respect canonical commutation relations and end up in a system of effectively identical TLLs, allowing for additional conservation laws than those associated to H 1 and H 2 (as for the symmetric quench). On the other hand having different u i but same K i does not modify the generic picture outlined in (39), which is characterized by the presence of multiple lightcones and regimes. In fact, in this case the difference in the two tubes can not be reabsorbed in a rescaling of the variables. One peculiarity of this limit, however, is the fact that the effective temperature of the symmetric mode in the prethermal regime (45) is zero.
Let us now turn to the final stationary regime reached by the dynamics. As we discuss in section V C and Appendix C, such regime, under the approximation of small momenta, is compatible with an equilibrium-like result associated to the two systems thermalized at temperatures T eff 1 and T eff 2 , in accordance with the classical equipartition theorem and the FDT in its classical (low frequency) approximation. While this appears as a generalization to a two temperature equilibrium state of previous results [13,30,56], it might sound surprising given that the underlying dynamics conserves the energy of each mode. In fact, a generalized Gibbs ensemble [17,18] would rather appear from Eq. (25), if we would take into account the full dependence on the momenta p in the integrals. However, as we have seen, the full dynamics of the asymmetric mode (see the bottom panels of Figure 2) and the stationary part for the symmetric one (see the regime within the first lightcone at short distances in the top panels of Figure 2) are well captured by the leading order term in p → 0 of the integrands which gives the expressions (39) and in particular (39e). The accuracy of the approximation made in (39) actually, is quite good also to describe the prethermal phase of the symmetric mode (as can be seen in the corresponding regime of Figure 2 and in Figure 4).
Moreover, in our discussion we referred to the regime (39d) (at least in the limit of u 1 u 2 ) as a prethermal one, in analogy with the work [30] (there, the dynamics stops at this prethermal regime due to the decoupling between symmetric and antisymmetric modes, which holds in the symmetric setting). More generally, prethermalization has been discussed in many works and it is often associated to a slow evolving intermediate state attained by the system before a complete thermalization takes place, as it happens in integrable systems in presence of a small integrability breaking perturbation [33,[77][78][79][80][81] or in other more exotic scenarios as in [82]. In this sense, the fact of considering different u i can be considered as a symmetry breaking mechanism that removes the degeneracy of the hamiltonian driving the dynamics. In fact, from Figure 6 (particularly if focusing on the inverse correlation length, bottom panel) one clearly sees the presence of a first rapid transient regime, followed by a quasi-stationary one for a relatively large time (divergent in the limit u 1 → u 2 ) and later evolving towards its asymptotic value. Note however that in order for the final state to be reached, the pretharmal plateau cannot be really time-independent and this is in fact clearly visible when looking at the correlation function C − itself (the top panel of the same Figure), which shows a slow ramp towards the final stationary regime. Note that this ramp can be increasing or decreasing according to the sign of T eff − − T eff , which can be tuned upon varying T 0 . About the main experimental implications of our results, one of the most surprising effects of considering two TLLs with different parameters is the (positive) linear correction in T 0 to the effective temperature T eff − of the antisymmetric mode, in contrast to the insensitivity of the same in the symmetric scenario [30]. This implies correlations decaying faster and it might be a non negligible effect in the dynamics, due to the relative high temperature at which experiments are carried out. For example, we would expect a similar correction to take place in the experiment discussed in [83], where a similar analysis can be carried out, while for now a theoretical understanding of the observed "effective" dissipation mechanism is still missing [58,59,84].
Remarkably, the phenomenological description of the unbalanced splitting protocol outlined in the aforementioned Ref. [67] for two bosonic tubes at different densities agrees in many aspects with the overall picture emerging from our general analysis of the quench dynamics in unbalanced TLLs coupled by tunnelling.
The transition from a prethermal to a thermal regime, both characterized by an exponential decay of correlation functions, with a multi light-cone dynamics that signals the sharp transition between different correlation lengths exists in both studies, as well as a linear correction in the initial temperature to the final effective temperature, shared by both the symmetric and the antisymmetric mode. This is due to the fact that their phenomenological input consists in an expression for the initial correlations, whose leading term behaves as p −2 in agreement with our derived expressions (see Eqs. (D1) and (D2)). This behavior leads to the exponential decay and a correction in the initial temperature showing up in the final correlation length. The relation T eff = (T eff − + T eff + )/2 in [67], connecting the prethermal and the thermal effective temperatures, is found to hold in our more general setting (cfr. Eqs. (51) and (52)).
There are however some interesting differences. In particular the case of density imbalance, ρ 1 = ρ 2 studied in Ref. [67] leads to the vanishing of the n ± mixing term in (11) (i.e., Λ n = 0 in (13)), while more general imbalances (for example, from different 1d interactions U i ) would allow for the presence of such a term. Due to the difference in the two protocols (the starting point of [67] is an imbalanced splitting of a single tube, while we start directly from two different tubes with non-zero tunnelling), our temperatures show a different dependence on the density as one can easily check by substituting the parameters (4) in our expressions. In particular, note that, in our protocol, for an imbalance in the densities only and in the limit T 0 → 0, we get that the temperatures of the two systems given in (46) are the same, i.e., T eff 1 = T eff 2 , and therefore they also coincide with the final temperature of the symmetric and antisymmetric modes. This, however, is not the case anymore at finite temperature T 0 = 0, and a linear correction in T 0 appears also to the prethermal temperature T eff − of the asymmetric mode, due to the difference between the two velocities.
It would be very interesting to test the previously highlighted features, displaying strong differences as compared to the equal TLLs scenario. This could be done e.g. in experiments similar to the ones of the Vienna's group [30][31][32]45]. Given the importance of the velocities in the dynamics, the presence of the harmonic confinement potential leading to a spatially dependent velocity is clearly a highly unwanted complication. Fortunately the recent realization of boxlike potentials in such experiments [85] shows great promise that the features analyzed in the present paper could be tested in a near future. Note that although we focussed here mainly on the phase correlation functions our analysis gives a fully diagonalization of the problem, so in principle other correlation functions are also easily accessible depending on the ones that can be accessed in future experiments.

VIII. Conclusions
In this work we have considered a quench in the tunneling strength of two TLLs with different parameters, under a quadratic approximation for the initial tunnelling term.
Our results show that the fact of considering two different systems leads to a much reacher physics than the one observed in presence of decoupling of the modes. This is manifested, for instance, in the emergence of multiple light cones. Moreover under this dynamics the prethermal regime discussed in [30] is followed by a final stationary state, never reached after the symmetric quench, where symmetric and antisymmetric mode display the same effective temperature (spatial decay). Due to the coupling between the modes one observes also an important effect of the initial temperature on the correlation length (effective temperature) measured via the decay of the antisymmetric mode, which otherwise would be only slightly modified in the limit of large initial masses.
Our prediction could be tested in experiments similar to the ones performed [30][31][32]45] for the symmetric quenches.
Beyond the current work the generalized Bogolioubov transformations developed in this paper allow us to address also different settings and a natural sequel of this work would be to consider the opposite quench, namely from a massless (uncoupled) initial condition to a massive (coupled) dynamics [86]. Another interesting direction to pursue is to understand the solution of the dynamics outlined in this work from the perspectives of a conformal field theory (CFT) approach [87], generalizing the ideas of [4,13,69] to the quench of two independent CFTs coupled by a (conformal) initial condition.
We want to diagonalize the hamiltonian (6). To this aim, we go to Fourier space, where it can be decomposed as The problem is thus reduced to the diagonalization of the 4 × 4 matrix H p . This can be achieved via a Bogoliubov transformation [88,89], which is a linear transformation B on the bosons b p . Restricting to real transformations, it has 16 free parameters. However, it has to satisfy some constraints [68]. First of all, the 4 bosonic modes defining b p are not independent, but are instead related (in pairs) by p → −p. This reduces the free parameters to 8, and constrains the corresponding Bogoliubov matrix to be of the form and the same for β, γ, δ. Moreover, we want B to preserve canonical commutation relations, i.e, where µ, ν = {1, 2, 3, 4}. This requirement leads to the condition namely B must be a symplectic matrix. Eq. (A5) is equivalent to If we take (α 2 1 − α 2 2 ) ≥ 0 and the same for β, γ, δ, the solutions can be parametrized by a Bogoliubov matrix of the form given in Eq. (15), with B ≡ B(φ p ) depending on a set of 4 parametersφ p = {ϕ 1,p , ϕ 2,p , ∆ p , φ p }. Finally, their value is uniquely fixed by the requirement for B to diagonalize H p . Note that this is not a standard diagonalization problem, because of the symplectic nature of B. The standard procedure [90,91] amounts to finding the spectrum of H p , by introducing the matrix H p J. This one can now be diagonalized in a standard way, meaning via a unitary transformation T as with Λ p diagonal (the corresponding spectrum in our case is given by Eq. (19) in the main text). Eventually, one imposes B t H p B = Λ p . This fixes the parametersφ p to be of the form given in Eq. (17). We start by considering the logarithm of C ± (x, t, T 0 ) defined in (22), i.e., [θ ± (x, t) − θ ± (0, t)] 2 T0 .