Quenches in initially coupled Tomonaga-Luttinger Liquids: a conformal field theory approach

We study the quantum quench in two coupled Tomonaga-Luttinger Liquids (TLLs), from the off-critical to the critical regime, relying on the conformal field theory approach and the known solutions for single TLLs. We consider a squeezed form of the initial state, whose low energy limit is fixed in a way to describe a massive and a massless mode, and we encode the non-equilibrium dynamics in a proper rescaling of the time. In this way, we compute several correlation functions, which at leading order factorize into multipoint functions evaluated at different times for the two modes. Depending on the observable, the contribution from the massive or from the massless mode can be the dominant one, giving rise to exponential or power-law decay in time, respectively. Our results find a direct application in all the quench problems where, in the scaling limit, there are two independent massless fields: these include the Hubbard model, the Gaudin-Yang gas, and tunnel-coupled tubes in cold atoms experiments.


Introduction
In recent times the theoretical understanding of out-of-equilibrium homogeneous systems in 1D has become central in statistical and condensed matter physics, as counterpart to the enormous experimental advances brought by cold atoms [1,2]. Several aspects have indeed been tackled by a variety of techniques, ranging from numerical methods (with particular reference to TEBD -time evolved block decimation [3]-and to DRMG -density matrix renormalization group- [4,5], and its time dependent extension [6,7]) to field theoretical techniques [8][9][10][11], integrability [12][13][14][15][16][17][18][19] and much more (see, e.g., [20][21][22][23][24][25] as more comprehensive reviews). Due to the complexity of generic out-of-equilibrium protocols, many studies focused on the simplified setup where the system is prepared in the ground state of some hamiltonian H 0 and then let evolve with a different one H: the famous quantum quench. While bringing important simplifications from a theory viewpoint, following the first remarkable example of Ref. [26], quenches have been realized in a variety of cases in cold atomic systems (see e.g. [27][28][29][30][31][32]). The possibility of experimental realizations triggered a corresponding theoretical effort to set a framework in which to study quenches [33][34][35][36][37].
While for free models, both on the lattice and in the continuum, quantum quench problems are most often analytically treatable as, e.g., reviewed in [15], it has been understood that the powerful tools of integrability out-of-equilibrium [13,14] can lead to exact analytic results only for a limited (but very interesting and experimentally relevant) class of initial states compatible with integrability [38]. Consequently, many interesting scenarios can be studied only with the help some approximate methods.
In this respect, when H is at or close to a quantum critical point, a very powerful approach is brought by conformal methods. Specifically, when the initial hamiltonian is massive, the problem can be tackled relying on an imaginary time path integral approach. In particular, in (1+1)D, the problem is mapped to a boundary conformal field theory (BCFT) one. This is the key result of the works by Calabrese and Cardy [33,35]. This description gives rise to exponential decay in time of correlations (with decay rate fixed by the initial mass). In contrast, when also the initial hamiltonian is massless, the correlations are expected to decay algebraically. Such behavior is recovered for generic systems relying on the Tomonaga Luttinger Liquid (TLL) paradigm [8,39,40], where initial and final hamiltonians are fully characterised by two parameters, known as sound velocity u and Luttinger parameter K. In this case, the power-law decays can be related in a simple way to initial and final Luttinger parameters only, as shown by Cazalilla in [34] (see also [9,[41][42][43][44][45][46][47][48][49][50][51][52][53][54][55][56] for some generalizations).
These results are somehow complementary, giving access to the dynamics after a quench starting from different classes of initial states. Still, it is important to keep in mind their range of applicability, especially when aiming at describing quantitatively the non-equilibrium dynamics of realistic microscopical models such as spin chains and quantum gases. In fact, at equilibrium it is very well established that the TLL approximation is quantitatively correct in the low-energy/large-distance regime; conversely because of the instantaneous nature of the quench after a quantum quench the system has an energy in the middle of the many-body spectrum and by no means a low-energy description is justified. Consequently, the extent to which Tomonaga Luttinger liquid theory can be used for the quantitative analysis of quench dynamics is a non trivial question. Of course, if the quench is near instantaneous but slow enough compared to high energy scales of the microscopic model (see e.g. [57]), the field theory description remains valid.
We will assume that we can fully describe the system by an effective field theory. In addition to the practical possibility of finding ramps with the proper speed for the field theory to remain valid, the theoretical study of the non-equilibrium dynamics of these conformal systems has its per se interest and provides very fundamental qualitative features that are difficult to get by other means in such generality (e.g., the previously mentioned exponential and power-law decays of correlations). Furthermore, many works attempted a detailed comparison between the CFT predictions and the actual non-equilibrium dynamics of lattice models, as e.g., [58][59][60][61][62][63][64][65][66][67] and, maybe surprisingly, it turned out that many features of the quench dynamics are not only qualitatively but also quantitatively captured by the TLL approximation. The effect of perturbations away from TLL model has also been analysed in [68][69][70][71] using renormalization group (RG) methods showing that the above mentioned studies provided a useful starting point.
Another crucial point is that the above-mentioned standard TLL approaches are only valid when the system under consideration consists of a single field or, trivially, of several independent (pre-and post-quenches) modes. On the other hand, in realistic (even 1D) systems this is not always the case: fermions usually carry spin, thus doubling the degrees of freedom, like in the celebrated (Fermi-)Hubbard model [8,72] as well as Gaudin-Yang gases [73,74]; more generally speaking, it is very common to end up in situations where the system at low-energy consists few coupled species of quasiparticles like in the experimentally relevant example of two (or more) tubes of interacting cold bosonic gases which are tunnelcoupled [75][76][77][78][79][80][81][82][83]. However, in the presence of more interacting species, the quench problem becomes very challenging both numerically and analytically. As examples, we mention the numerical analysis performed in Refs. [84][85][86][87][88][89][90][91][92][93][94] with a variety of methods. The situation is even more complicated on the analytical side, where only few exact results are available and mainly focused on the characterization of the final steady state [95][96][97] or work in some limits/approximations [98][99][100]. It is thus clear that any semi-quantitative, or even qualitative general picture for such problems would be not only useful, but highly desirable and this is the strategy employed in some related works appeared in the literature [101][102][103][104][105][106][107]. A great simplification occurs when the Hamiltonian has some symmetry and the problem with two degrees of freedom can be studied by introducing a suitable change of variable leading to effectively decoupled modes. In the case of tunnel coupled condensates this has allowed to study the quench of two identical tubes with a mass coupling between the two that is suddenly removed [104,[108][109][110][111][112][113].
The study of the problem is, instead, much more complicated when no such obvious change of variable allows to reduce the problem to two decoupled modes, see for instance [105]. In this direction, few recent studies investigated the quench dynamics in the tunnelling coupling of two TLLs with different sound velocity and/or Luttinger parameter [114][115][116], aiming at understanding the effect of such "imbalance" between them. Using a semiclassical approximation, the problem was solved via a Bogoliubov approach. For those specific situations, this approximation was shown to give access to a very rich phenomenology [116], with (i) the emergence of multiple lightcones, separating different decaying regimes; (ii) a prethermal regime eventually decaying into a quasi-thermal one; (iii) non-trivial effects of a non-zero temperature in the initial state. However, how much the above results depend on the specific quench considered there is not clear.

Goal and main results
The goal of the paper is to understand the quench dynamics of two initially coupled and generically different Luttinger liquids, which are let evolve independently, by relying on conformal methods where (differently from other approaches) general universal aspects are expected to clearly manifest.
In order to do this, we take advantage of the explicit solution of a particular (quadratic) quench problem, i.e., the one studied in Ref. [116], and in view of possible generalization, we focus on the form of the initial state in terms of the post-quench modes b † i,p (i = 1, 2) (rather than looking at it as the vacuum of the pre-quench modes, as one would do in standard Bogoliubov approach). This turns out to have the following squeezed form: where N is a normalization, 0⟩ the vacuum of the post-quench hamiltonian (two uncoupled TLLs), and W p is a two by two matrix containing all information the initial state. In particular, we realize that a crucial information is encoded in its low-energy expansion, i.e., W p=0 , and, specifically, in its eigenvalues: an eigenvalue equal to one will be associated to a massive mode, while an eigenvalue smaller than one to a massless mode. This is consistent with the coefficients of the squeezed initial states emerging in massive and massless quenches in a single TLLs (see Section 3). In order to get meaningful results for the dynamics, however, one must take into account the next-to-leading order term in W p . This is also something that is corroborated by our understanding of the massive quench for single Luttinger liquid (cf. Eq. (13)). For a generic theory of two different Luttinger liquids the next-to-leading order approximation of W p in (1), results in a matrix at linear order in p which presents an off-diagonal term coupling the two diagonal modes.
In this work, focusing on the case of one massive and one massless mode, we show that if one is interested in the leading behavior -the meaning of "leading" will be clear in the following -for large time and space separation, one can safely disregard such coupling term and solve the quench of two decoupled modes. This comes at the price of dealing with multitimes correlation functions, which can be ultimately traced back to the different speed of sounds characterizing the two systems.
The generality of our approach lies in that the state in (1), with W p associated with one massive and one massless mode, can be seen as an effective description of the ground state of two LLs coupled by a generic RG relevant term, in the sense that it reproduces the leading order of equilibrium correlation functions, and is here assumed to be also a good starting point from the out-of-equilibrium problem. The main effect of going from a quadratic to a more complicated coupling (like a cosine term, associated with a truly interacting Hamiltonian) is expected to be a change in the value of some parameters in the matrix at the next to leading order: this is something well known in the context of a single theory and has been widely used for non-equilibrium settings (see Ref. [35] or Section 3 for a brief discussion). While the exact value of such parameters cannot be computed exactly in genuinely interacting theories, our interest here is in the functional dependence of correlation functions.
Altogether our method allows us to understand the underlying structure of correlations functions in terms of the two quasi-independent modes, and compute straightforwardly several of them. Specifically, we first compute the one-and two-point correlation functions of vertex operators of symmetric and antisymmetric sectors (see below for proper definitions), given by Eqs. (36)(37)(38)(39)(40) and Eqs. (42-44) (respectively). Moreover, out of this approach, we can easily get other correlators, e.g., the correlations of density, Eqs. (51)(52), and currents, Eqs. (53)(54), which, in the particular case of a quadratic Hamiltonian as the one studied in [116], were not clearly accessible within a Bogoliouv approach.

Organization of the paper
The paper is organized as follows. In Section 2 we set the problem, including some reminders of the results obtained in Ref. [116] within the semiclassical approximation. In Section 3 we give a brief overview of the two approaches to quantum quenches in conformal field theory for a single field. The problem of two initially coupled TLLs is then studied relying on such results in Section 4 and Section 5. Specifically, Section 4 discusses in details the main ingredients needed for the calculation of correlation functions, which is then carried out in Section 5. We finally conclude in Section 6. In order to keep the paper fluid to read, we chose to collect most of the calculations in four appendices.

Setting of the problem
As mentioned in the introduction we are interested in studying the time evolution after a quench in which the post-quench low-energy physics is captured by two different Tomonaga-Luttinger liquids (in the CFT language, two free compact bosons).
Without loss of generality we can write the post-quench hamiltonian as with (i = 1, 2) The fields θ i , n i in (3) are related to the bosonic creation/annihilation operators b ( †) i,p via: where L is the system size, and we introduced an ultraviolet cutoff Λ in order to ensure convergence of correlation functions. In the rest of the paper we focus on the thermodynamic limit (TDL), namely infinite system size.
The initial state, of the form (1), is assumed to couple the two Hilbert spaces associated to i = 1, 2. Otherwise, the problem factorizes in the initial variables and we go back (at leading order) to the dynamics of two decoupled single fields, which, as mentioned in the introduction, has been already largely addressed. Moreover we will assume that it describes a massive and a massless mode. This is the case of the ground state of considered in [116], or of a generic hamiltonian with a quadratic coupling in the fields θ i (which can be recast in the hamiltonian (7), at the price of having renormalized Luttinger parameters). More generally adding to the hamiltonian of two LLs a coupling term relevant in RG sense (e.g., a cosine) will open a gap in the spectrum, thus giving rise -at equilibrium -to exponentially decaying correlations (associated with an effective mass), plus possible powerlaw corrections (multiplicative and/or additive ones). Such behaviour can be effectively reproduced by a state of the form (1) with a specific low energy expansion (mentioned in the introduction and explicitly given in Eq. (27) below). Crucially, the correlations functions of such a state at large distance are characterised only by a mass which can always be though as an effective parameter in the ground state of two TTLs. What is non-trivial, then, is that this ansatz for the state still correctly describes the dynamics after a quench. In general this is not guaranteed already for a single theory, as subleading corrections might become dominant at late times. In the case of an interacting post-quench hamiltonian, indeed, whether a squeezed state approximation is justified in some regimes is subject of current research [121][122][123][124][125]. If the final hamiltonian is critical, instead, arguments of RG theory of boundary critical behaviour [120] have been used to justify a squeezed state form [35] (see Section 3.1 below). Since we focus on a critical final hamiltonian, we assume that this is still the case for two theories. Whether at some later time, our assumption on the initial state breaks down is an open problem, that, nevertheless, is beyond the aim of this paper.
In the following sections we discuss the dynamics of the symmetric and antisymmetric modes with speed of sound u and TLL parameter K equal and given by θ ± are often relevant in systems with two types of degrees of freedom, as the Hubbard model or the tunnel-coupled condensates. While these variables are obviously the most appropriate ones in the case of two identical systems with coupling in (θ 1 − θ 2 ), as the initial and final Hamiltonian are decoupled in this basis and the quench occurs only in the antisymmetric sector [117,118], this is not true in general when the two systems are different.

Reminders of Bogoliubov approach and some notations
To make contact with Ref. [116], we recall that in the Bogoliubov approach to the quench dynamics [15], the initial state is assumed to be the ground state (or even another eigenstate [119], but we do not consider this case here) of a quadratic hamiltonian, as, e.g., (7). Hence, the standard way to solve the quench dynamics of interest is to exploit its quadratic nature and perform a Bogoliubov transformation to bring it in the following diagonal form so that its ground state is factorized as the product of the vacua of the pre-quench modes: Here we are interested in theories that have a massive (m) and a massless (0) mode, i.e. in which the small momentum behavior of the two dispersions reads For the hamiltonian (7), considered in Ref. [116], we have m 0 = gu K which is the mass of the massive mode, and v = u 1 u 2 K 1 K 2 K the speed of sound of the orthogonal (massless) mode. As mentioned, also in the general case, we can think of the initial state as the ground state of an effective quadratic hamiltonian. In this case m 0 is the effective mass (whose value depends on the interaction details), while v is the sound velocity of the massless mode.

Preliminary results: quenches in a single CFT
Before studying the setting of two coupled CFTs introduced above, we shorty review the results available for the simpler case of a quench in a single CFT. In this case, the postquench Hamiltonian is H u,K [θ, n], while the initial state can have or not have a gap (i.e., an effective mass).

Massive quench
This quench has been solved exploiting the conformal invariance of the problem, considering an imaginary time path integral approach, that we recall in this section. The results of this method have been developed in [33,35], and later clarified and generalised in [126][127][128][129][130][131]. The framework is quite general and applies to quenches starting from a translationally invariant massive state ψ 0 ⟩, namely any state with short-range correlations.
The objects of interest are expectation values of local operators φ j (x j ) after the quench, namely In imaginary time, Eq. (12) can be represented as a path integral over a strip with operator insertions and ψ 0 ⟩ playing the role of boundary condition imposed at initial and final times. A crucial point is that, exploiting the powerful tools of Renormalization Group (RG) theory of boundary critical phenomena [132], a short-ranged initial state ψ 0 ⟩ can be always replaced by the appropriate RG-invariant boundary state B⟩ to which it flows. The distance of the actual boundary state is taken into account (to leading order) introducing an extrapolation length τ 0 , and approximating the state as ψ 0 ⟩ ≃ e −τ 0 H B⟩, with H the post-quench conformal hamiltonian. Note that in terms of the creation operators b † p , it takes the form of a squeezed state, where 0 b ⟩ is the vacuum of bosons of the final Hamiltonian. The extrapolation length τ 0 is expected to be of the order of the inverse mass, e.g., τ 0 ∝ 1 √ g in the case of pre-quench hamiltonian (7), and, more generally speaking, of the inverse gap for a gapped interacting theory [35,132]. Accordingly, Eq. (12) in can be rewritten as where the problem has been mapped to a boundary conformal field theory (BCFT). Namely, Eq. (14) is given by a path integral over a strip of width 2τ 0 with conformally invariant boundary conditions, and operators inserted at τ j (with τ j ∈ [0, 2τ 0 ]). Eq. (14) is often denoted as ⟨φ 1 (x 1 , τ 1 )⋯φ n (x n , τ n )⟩ slab(2τ 0 ) : we will use this convention in Appendix A. One can then rely on standard CFT calculations, based on conformal maps and on the transformation of operators under those, to compute (14) exactly. The imaginary times τ j have to be analytically continued to τ 0 + it j as the final step, to recover the real time evolution. Within this framework, very general results can be obtained for n-point correlation functions, which show exponential decay in time before relaxation, with the appearance of the famous lightcone effect [33,35,[135][136][137][138][139][140]. The steady state shows a finite correlation length typical of a thermal system and the deviations from a thermal state generated by the integrability of the model are small scale details not captured by the too simplistic approximation of the initial state (the modifications necessary to observe the relaxation to a generalized Gibbs ensemble [141] within this approach have been worked out in [126]).

Massless quench
The path integral approach of the previous subsection does not apply to initial critical states, which are long ranged and therefore would be associated to a diverging extrapolation length (i.e., a vanishing gap). This case has instead been considered in Ref. [34] (see also [68,69], that we closely follow in terms of notation, and [9] as review on the subject), where the quench dynamics of a TLL after a sudden change of the TLL parameter, say from K 0 to K f , is studied via a Bogoliubov approach.
In this case, the initial hamiltonian is diagonal in some operator basis η p , and the final one in some other basis b p . They are related by a Bogoliubov transformation Note that this diagonalization also holds when the quench occurs in the sound velocities as well (i.e., for the more general case The ground state of the initial hamiltonian can be written, again, as a squeezed state in the final basis, i.e., where N = ∏ p>0 (1 − W 2 ) −1 2 is the normalization factor. A crucial difference as compared to the boundary state (13) is that here W < 1: this is ultimately responsible for the power-law decay of correlation functions in the state (16) versus the exponential one in (13). The relaxation is towards a genuine non-equilibrium steady state, namely a generalized Gibbs ensemble [141], determined by the underlying integrability of the model. In fact, the late-time spatial decay is power-law and governed by an exponent that is different from the one that governs asymptotic ground state correlations (i.e., K f ). In particular this Luttinger parameter gets renormalized by a function of the ratio K 0 K f [68,69], as one might expect from the transformation in (15).

Initial state and operators' dynamics
In this section we initiate the conformal field theory study of two initially coupled TLLs in the setting of Section 2. We discuss first the low energy properties of the initial state, and then the operators dynamics in the Heisenberg picture.

Leading order for small momentum
A state of the form (1), with generic W p , cannot be directly handled with conformal methods because of the non-trivial dependence of the momentum in W p . However, we anticipated that, invoked RG ideas, we can focus on its low energy limit, i.e., the limit p → 0 of W p . In the basis of the post-quench hamiltonian, a zero-momentum matrix with a massive and a massless mode can be parametrized as follows where I is the identity matrix, and ϕ and ν such that (recall that we can always think of the state as the ground state of a quadratic hamiltonian of the form (10) and the value of the effective mass does not enter in the parameters below) The matrix W 0 can be diagonalized via a rotation, parametrized by an angle ν The two eigenvalues of W 0 are {1, cos 2ϕ}, associated respectively to two orthogonal sectors that we dubbed {A, B}, and the value of these eigenvalues determines the spectrum and the decay of correlations of the two modes, as we saw in the previous sections and can be understood from (104). In the case of two identical systems u 1 = u 2 and K 1 = K 2 , the two modes are associated to the antisymmetric/symmetric fields θ ± introduced above. However such correspondence does not hold in general.
The state (1) in the low-energy approximation is thus factorized in the basis {A, B}. It consists of an infinite mass state (associated to the eigenvalue 1, see Eq. (13)) in the A-sector, and a massless state (associated to cos 2ϕ < 1, see Eq. (16)) in the B-sector. Importantly, these two are the low energy states that characterize the dynamics studied, respectively, by Calabrese-Cardy [33,35], and by Cazalilla [34], as discussed in Sections 3.1 and 3.2. The dynamics in the B-sector can be interpreted as a quench in the TLL parameter. Indeed, for the ground state of the hamiltonian (7), it corresponds to a quench from Γ in the initial state to K + in the post-quench hamiltonian, as follows from identifying in Eq. (16) tanh δ with cos 2ϕ = Γ−K+ Γ+K+ . Crucially, the factorization of the state at this order is, by construction, independent of the momentum p, in such a way that we can define the fields The TLL parameters K A and K B are auxiliary variables which are free variables and, as we will see, they will not enter in the formulas for the dynamics of physical observables. Always to have in mind a specific example, we can plug in the above equations the value of ν in Eq. (18), corresponding to hamiltonian (7), to have This equations shows that the field θ A remains aligned to the massive hamiltonian term (θ 1 − θ 2 ), cfr. Eq. (7).
By inverting Eq. (21) for {θ i , n i }, and plugging them into the post-quench TLL hamiltonian (2), we get with the coupling of the A − B sectors given by so, in general, the two sectors are coupled. Moreover, we have that fix the sound velocities u A B of θ A B . In the case of the ground state of hamiltonian (7), these two velocities are u A = K 1 K 2 K+ u K and u B = 1 K+ uK. We conclude this subsection with two comments. The modes A B allow us to write the initial state as a factorized squeezed state at low-energy with operators acting on the physical vacuum of the post-quench hamiltonian. This is different from writing the state as product of the two pre-quench vacua of (10). As a little detour, we note that the rotation that diagonalizes W 0 in (17) is the same one introduced in Ref. [142] for permeable interfaces in CFT. Indeed, there the scattering matrix is just given by either S + (ν) or S − (ν). Given that W 0 and S − commute, they are diagonalized by the same transformation. This observation is the starting point for a possible connection between permeable interfaces and quench problems that will be investigated in a forthcoming work [143].

Beyond the leading order
The next to leading order in p of the initial state ψ 0 ⟩ in general breaks the factorization in the A B sectors. It is then convenient to write these next-to-leading order terms in p directly in the basis A B, in which the initial state has the general form where the normalization N is reported in Appendix D, see Eq. (109). Clearly τ A , τ B and γ are functions of the initial parameters {K i , u i } (e.g., for the hamiltonian (7) they can be explicitly worked out), but the precise functional dependence is actually not needed. The two velocities u A and u B in W p are defined in (25) and (26), respectively.
is nothing but the extrapolation length of the massive quench introduced ad hoc in the previous section (cfr. Eq. (13)). This length is interpreted as the "distance" from the infinite-mass state. As already mentioned, it is expected to be of the order of the inverse gap m −1 0 . This term is the one generating exponential decay of correlation functions.
The O(p) correction to the BB-component (parametrised by τ B ), only produce subleading corrections, as shown in Appendix D. Hence it is neglected in what follows.
The factorization of the initial state is spoiled by the presence of the off-diagonal matrix element γ ≠ 0. To proceed, however, in the following we are going to assume a diagonal form of the state also at this order (i.e., γ = 0). The consequences of a non-zero value of γ will be discussed for the correlation functions under consideration. As we are going to argue, the role of γ is to renormalize subleading power-law exponents in some cases, or just modify non-universal prefactors in other. Crucially, it will never affect the leading behaviour of the correlation functions of interest.

Decoupling of operator dynamics
We are going to work out the non-equilibrium dynamics in the Heisenberg picture, in which the time dependence is entirely encoded into operators, while the state does not evolve. A suitable rescaling of the times will allow us to always write our observables in a decoupled form with respect to the A and B degrees of freedom.
Let us focus on the field θ 1 , for θ 2 the derivation is identical. Its dynamics is given by In the above equation we can replace with u 0 a common (arbitrary) velocity for the two TLL hamiltonians of θ 1 2 . The presence of H u 0 ,K 2 [θ 2 , n 2 ] does not affect the dynamics of θ 1 because the two commute. Importantly, due to the rescaling of time, now θ 1 2 have the same (auxiliary and fictitious) sound velocity u 0 . Then (cfr. Eq. (23)) where we used that, in these rescaled time, u 1 = u 2 = u 0 implies (see Eqs. (25)(26)) u A = u B = u 0 as well, and λ AB = 0 in (23). This is crucial, because now the hamiltonian in the rhs of (30) acts separately on θ A B . Finally, defining the rescaled times the rhs of (29) can be recast in the form which plugged in (28) gives where the second equation for θ 2 follows from a very similar calculation. In fact, (33) is nothing but the time-dependent version of Eq. (22). The rescaling of time introduced above is particularly important when considering observables which are functions of both θ 1 and θ 2 , such as, for example, of the symmetric and antisymmetric fields θ ± , which decouple in terms of θ A B at any time In summary, the general idea, exploited in the following section, is to use a time rescaling to reabsorb the different velocities of the two initial LLs into the times at which the observables are evaluated. Hence, using rescaled modes with the same sound velocity is enough to ensure an exact decoupling into the time-dependent θ A and θ B at any time. The price to pay is that equal-time observables and correlations become multi-times ones. This is evident in (34), where a single time in the lhs results in two different times in the rhs.

Correlation functions
In the previous section, we achieved the two necessary conditions to compute correlations functions, namely • the factorization of the state in the basis which diagonalizes the fields θ A B (assuming to neglect the coupling γ in Eq. (27)); • the decoupling of the operator dynamics with respect to the same basis.
Using these properties, all the correlation functions can be computed independently in the massive and in the massless sector as multi-point functions at different times. In particular, for the massive sector, one can apply the method developed in [33,35], while the computations in the massless sector are equivalent to those in [34]. In the following, we will see how our conformal approach provides the correct result of the leading decay of the correlation functions.

One-point functions of θ ±
As a first non-trivial example, we consider the exponential one-point functions of θ ± (i.e., vertex operators in CFT language) where ⟪⋅⟫ γ denotes the expectation value over the state (27), and below we consider γ = 0. These correlations of θ ± are the experimentally relevant ones in the context of tunneled-coupled tubes in cold atoms experiments [83]. They are also the most natural also in the Hubbard and Gaudin-Yang models, where they are associated with spin and charge sectors [72]. Making use of the decomposition derived above when γ = 0 (cfr. Eq. (34)), Eq. (35) can be cast in the following form where from now on expectation values of observables which are functions of θ A are understood to be taken on the AA-component of the state (27). Similarly, functions of θ B are evaluated on the BB-component of (27) with τ B = 0 (as already mentioned, it does not contribute up to subleading corrections). Moreover, to lighten the notation, we simply used t i instead of t a i , because the correlation univocally specifies whether a = A or a = B.
Thus, we conclude that one point functions of exponential of θ ± become a product of two-point functions of θ A B . Now massive and massless part can be computed separately.
We start from the massive part. Using the approach of Section 3.1, the object that we need to evaluate is a path integral over a slab of width W = 2τ A with the two vertex operators (which are primary operators in the CFT) inserted at different points. The final result reads (up to a unimportant prefactor) where we used Eq. (60) in Appendix A with the specific values h 1 = 1 K 2 for the conformal weights of the corresponding vertex operators. Eq. (37) reproduces the leading behavior for large t of C − 1,γ (t) obtained in Ref. [116] via a Bogoliubov calculation if we identify which is consistent with the standard interpretation of τ A as inverse initial mass gap [33]. Let us now move to the B mode. In this sector the initial state is massless, so that we can use the results in [34] for the two-point function, but generalized to the case of unequal times. This is, once again, a standard Bogoliubov calculation (that we report in Appendix B). The final result reads (cos(2pu 1 t) + cos(2pu 2 t)) ± cos(p(u 1 + u 2 t)) and K ± = Γ K+ ± K+ Γ (cfr. Eq. (75), specialized to our quench with K 0 K f = Γ K + ). According to the ± sign, this integral has or has not an infrared (small p) divergence. Specifically, such divergence gives C + 1,0 (t) = 0 when plugging Eq. (40) in (39) (with the + sign). This in fact is the correct result for C + 1,γ (t) known also from Ref. [116] (see Eq. (82) in App. C), which remains valid for γ ≠ 0. Conversely, Eq. (39) (together with (40)) leads to a power-law decay for C − 1,0 (t). The exponent is however different from the one found for γ ≠ 0 (cfr. Eq. (82) in App. C). We will discuss this discrepancy in Section 5.3. For now, we just point out that at leading order log C ± 1,γ (t) = log C ± 1,0 (t), while the same is not guaranteed for subleading corrections (the logarithm is important for the correctness of this statement, since the corrections from the massless sector to C ± 1,γ are multiplicative).

Two-point functions of θ ±
Similar results can be found for the (exponential) two-point function of θ ± As before, we start by rewriting it in a factorized form for γ = 0, i.e.
⟩ . (42) Therefore, vertex two-point functions of θ ± are mapped into the product of two four-point functions of θ A B , that we can compute separately. For the massive part, we now have a four-point function to be evaluated in the same strip geometry considered before. The result is given by (see Appendix A, Eq. (65)) where, without loss of generality, we assumed u 1 > u 2 . Again, if we fix τ A as in (38), this reproduces the correct exponential decay of both C ± 2,γ found in [116]. Given that the B sector provides algebraic correlation, the exponential contribution in Eq. (43) represents always the leading decay both in x and t, as already pointed out in Ref. [116]. A possible special case is the short time regime 2u 1 t < x (first case in (43)) where there is no x-dependence. Hence, the possible space dependence is entirely in the subleading powerlaw contributions which we now study. The result for the B-part of this two-point function Note that from (44) one can easily read off all the different regimes (the same as in Eq. (43)), sharply separated by lightcones. Those are nonetheless smoothen out when reintroducing the ultraviolet cutoff [116]. Note also that both K A and K B cancel in the above expressions.
In the aforementioned regime of short time (x ≫ u 1 t), Eq. (44) gives that C − 2,0 is constant (the various exponents sum up to zero) and C + 2,0 ∼ x (K 0 K f ) (K+) = x 1 Γ . In this regime, these correlation functions match exactly the results from the Bogoliubov calculation in [116] (reported, for completeness, in Appendix C).
In the other regimes, the power-law scaling in Eq. (44) have exponents that are, in general, different compared to the ones for γ ≠ 0. As mentioned already for the one-point functions, this disagreement represent the limits of the conformal method, which does not gives access to all power-law contributions. Anyway, the conclusion also for the two-point function, is that the leading term is well captured is all the regimes. We will come back to this issue in Section 5.3.

Derivative operators
We focus here on fluctuations of the initial fields (i = 1, 2) where j i (x, t) the current density associated to θ i (x, t). Density and current correlations can be related to correlators of the derivative operators, which is also a primary operator of the CFT [148]. In fact it holds For definiteness, below we look at density-density correlations, while results for the currents can be similarly derived. Following the same logic used for the vertex operators, we exploit the factorization of the state (27) and the decoupling of observable to get We see that the first term in the above equation is associated to the massive mode, and therefore, according to [33,35], decays exponentially in time. Hence, the leading term is now given by the part associated to the massless mode, giving rise to a power law decay, according to [34]. Such power-law decay comes out very naturally within this approach.
To get this leading term, we compute where we defined s j ≡ s a j = (u j u a )s. This is a two-point function at equal times when i = j and at different times when i ≠ j. In both cases we can evaluate it using (48) and the results in Appendix B. For i ≠ j we get (49) and in the case i = j The density-density fluctuations finally read (i ≠ j) where above we defined c = x t < ∞, and (for i = j) Similarly, for the current-current correlations we get for i ≠ j and, for i = j, As we are going to justify in the following subsection, the leading algebraic decay is not influenced by the inclusion of a γ ≠ 0 term.

Corrections from off-diagonal terms and comparison with Bogoliubov approach
In the previous two subsections we calculated several correlation functions following standard RG ideas which are completely under control at equilibrium. In particular, we focused on the first terms in a small momentum expansion of the initial state. However, it in unclear how these RG reasonings capture the real out-of-equilibrium time evolution of the two generally coupled TLLs. The crucial point is represented by the generic form of the initial state in Eq. (27) which shows that at order O(p) (i.e., for γ ≠ 0) a term breaking the factorization of the initial state appears, modifying qualitatively a few aspects of our approach. It is straightforward to realize that this term generates further algebraic decay in space and time separations, as can be understood from the results derived in Appendix D. As a consequence, for all the correlation functions with a leading exponential behavior, the presence of such off-diagonal coupling only provide subleading corrections to the result we obtained assuming factorization. Conversely, power-law terms are in general affected by the presence of the off-diagonal term, and so are not correctly captured within our approach working within a factorized initial state. Interestingly (and maybe surprisingly), for all the correlations presented above, every time that the leading term is algebraic, such off-diagonal term always leave it untouched. Let us consider as a first example C ± 2,γ (x, t), defined in Eq. (41) of Section 5.1. In this correlation, the leading behaviour with the assumption of factorized initial state is given by the exponential terms in Eq. (43) coming from the massive sector. The massless sector provides only a correction to the power-law multiplicative correction reported in Eq. (44). The exact result (with the correct powers), as obtained from the Bogoliubov approach, is reported in the appendix, cfr. Eq. (97). It is evident that the two are in general different: in fact, while (44) only has two free parameters (i.e., K ± ), in (97) we have four of them (i.e., Θ and {a ij } in (97), defined in Appendix C). Remarkably, they exhibit the same lightcone structure.
Nonetheless, as anticipated in the previous section, it turns out that in short time regime also power laws are correctly captured. This agreement does not come as a surprise, because, in the short-time regime (namely, before the first lightcone), the correlations reduce to the ones in the initial state. Anyhow, this obvious result comes from a non-trivial limit and was worth to test. Moreover, since this is the only regime where power-laws become leading (cf. (43)), the conclusion is that the leading term in both x and t is always correctly captured by our approach. Note that in the intermediate regimes x and t have to scale in the same way (by definition, x t must be finite and within the limits defining the corresponding regime). Therefore in this case, we have effectively just one independent variable, and the leading term is exponential. Coincidently, in the specific case K 1 = K 2 , the exponential decay in x also vanishes in the prethermal regime (i.e., 2u 2 t > x > u 1 − u 2 t), and one can verify that the power law correction is correctly described as well.
The other correlation function of interest is the density-density one that we considered in Section 5.2. The effect of having γ ≠ 0 in this correlation is to add a mixing term ⟪∂θ A ∂θ B ⟫ γ to Eq. (47). However, it is clear that this further term decay with the same power-law as the leading B/B term in Eq. (47): therefore the presence of γ only changes the prefactor in front of the power-law decay.
Note that there is a main difference between the two examples discussed above. For the vertex operators, the correlations of θ A B appear in the exponent. As a consequence, the non-diagonal contribution from γ ≠ 0 are multiplicative ones and therefore they renormalize the power-laws. For the derivative operators, instead, the corrections in γ are additive, and thus they do not change the power-law exponent, but just modify the non-universal prefactor.

Particular limits
In this short subsection, we analyze how our general quench simplifies when the velocities or the TLL parameters of the two species are the same.

Same velocities u 1 = u 2
When u 1 = u 2 (even for K 1 ≠ K 2 ) in the post-quench Hamiltonian (2) H u 1 ,K 1 and H u 1 ,K 2 have the same spectrum. As a result, for any choice of the basis parametrized by ν in Eq. (20), the hamiltonian corresponds to two decoupled TLLs, as can be seen from (23). The initial state then selects a particular angle ν (cf. Eq. (18)) that ensures the decoupling between the massive and the massless degrees of freedom in the pre-quench Hamiltonian (7). In this case therefore, the initial state is exactly factorized in the two sectors and consequently the solution of nonequilibrium dynamics does not require any time rescaling (as a trivial consequence of the equality of two velocities). Hence, one can compute simply equal time correlation functions, and the light cone structure is largely simplified: this is evident, e.g., in (43) where one is left with two regimes only, corresponding to a unique lightcone.
Because of the perfect decoupling, the initial quench in the two TLLs induces a quench in the massive sector only, whereas the massless sector remains at equilibrium in the ground state. In fact the massless quench induced in the general case u 1 ≠ u 2 results from an effective sudden change in the TLL parameters, but in this case they turn out to be equal. Specifically, one has Γ K + = 1 if u 1 = u 2 (as can be checked from Eq. (19)). All previously reported results match this particular limit, but for some correlations in a singular manner requiring the restating of the short-distance cutoff.

Same TLL parameters
Contrarily to the limit of equal velocities, the one of equal TLL parameters does not bring major simplifications: the lightcone structure in (43), for example, remains since it is clearly related to the presence of different velocities only.
However, we note that in this case, in the correlation function (43), there is no exponential decay in space in the "prethermal" regime (namely 2u 2 t > x > u 1 − u 2 t) for C + 2,0 . Although there is still non-trivial time dependence in (43), the lack of exponential in x can be interpreted as a prethermal "temperature" equal to zero. This last fact can be understood, at speculative level, by noting that in the prethermal regime the difference in the velocities is small compared to the the considered spacetime scale (x t ≫ u 1 −u 2 ). Therefore, as K 1 = K 2 , the total system is basically equivalent to two identical TLLs, and in that case the symmetric mode remains in its ground state [111]. In this regime the power law decay in space, coming from the massless sector, which becomes dominant, is correctly captured by our low energy approach.
We have shown that, for what concerns the large scale properties, this non-equilibrium dynamics decouples into two independent sectors, inducing an effective quench in each of them: one starting from a massive initial state and one from a massless one with an effective TLL parameter. Each of them can be studied by means of the known techniques reported in Sections 3.1 and 3.2. The equal-time correlations of the coupled system map to correlation functions at different times of two uncorrelated modes. We have also discussed that, while the leading term and the light cone structure are always well captured by our approach, this is not the case for subleading power law corrections, generated by the coupling between these two modes (i.e., γ ≠ 0 in Eq. (27)).
Moreover, a direct inspection of the correlation functions shows that while for vertex operators the leading contribution is given by the exponential decay of the massive mode and the massless mode acts with subleading multiplicative power law corrections, for derivative operators the large scale properties are determined by the power law decay of the massless mode, while the massive part constitutes an additive correction here. Therefore while vertex operators mimic a thermal like behavior, derivative operators behave as at T = 0 or more generally in a GGE [141].
A final interesting remark is that for the quench induced in the massless sector, the product of the speed of sound and the Luttinger parameter before and after the quench is equal (i.e., vΓ = u B K + in our notation). This fact suggests that the quench respects Galilean invariance [39].
This work paves the way to the study of quenches in systems consisting of more than a single TLL. We conclude by providing some future possible developments in this direction. The most natural extension of our calculation would be to apply our CFT approach to initial states in which the orthogonal modes correspond to two massive or two massless theories. To this aim, it is sufficient to change the parametrization of the initial squeezed state (27) with two eigenvalues with absolute values equal to one (massive case) or smaller than one (massless case).
Another generalization would be to consider a larger number of initially coupled TLLs. In this case, there are many different physical situations requiring different numbers of massive and massless modes. In particular, in the case where the TLLs correspond to many different tubes, it would be interesting to investigate how two-dimensional non-equilibrium physics emerges from the coupling of 1D systems, as done in a very different setup in [85].
Finally we mention possible connections with the works on conformal interfaces [142,[144][145][146][147], that we plan to investigate in the future [143]. Our framework, in fact, can be in principle reformulated in a full path integral fashion via the unfolded picture [144,145], where the initial state (living in the tensor product of two CFTs, i.e., CFT 1 × CFT 2 ) is mapped to an interface (connecting two spatial regions in a single CFT, namely CFT 1 ⋃ CFT 2 ). acknowledges support from ERC under Consolidator grant number 771536 (NEMO).

A Calculations in the massive sector
In the path integral approach in imaginary time developed in [33,35], quantities as the one in Eq. (37) are mapped to correlation functions in a strip geometry with boundary conditions corresponding to conformally invariant boundary states. In (1+1)-dimensional BCFT, those are computed exploiting the transformations of correlation functions of (primary) operators under conformal maps. Let us consider, for example, two geometries in the complex plane with a boundary (say G1 and G2) with coordinate w and z, related by the conformal map w(z). Then the correlations of primary operators φ i in the two geometries are related as with h i andh i being the holomorphic and anti-holomorphic dimensions of φ i .

A.1 Two-point function in the slab
To get the two-point function in (37), the object that we need to evaluate is a path integral over a slab (of width W ) with the operators inserted, i.e., with V α = e iαθ and θ is a bosonic field. Moreover σ i (i = 1, 2) are imaginary times, to be analitically continued to the values σ i → it i + W 2 at the very end of the calculation. Moving to complex coordinates (with points labelled by w = r+iσ), the correlation function of vertex operators V α i (w i ,w i ) on the slab geometry is first mapped by a conformal transformation to the upper-half plane (UHP) (with coordinate z s.t. Im(z) > 0). The two-point function in the UHP is related to a four-point function of chiral vertex operators V α (z i ) on the complex plane, z ∈ C [148]. The details of the calculation can be found in [35]. The final result is conveniently expressed in terms of where J = (π W ) 2(h 1 +h 2 ) denotes the Jacobian factor in (55), h i = α 2 i (8K A ) is the conformal weight of the chiral operator V α i (z), γ i = πσ i W , z ij = z i − z j , zī =z i , and R i = e π W r i (in our case R 1 = R 2 ). Plugging into this expression the actual coordinates, we get Finally, by analytically continuing σ i to real times and taking t i ≫ W , we obtain Upon specifying the values of α i (i = 1, 2), the above equation allows to access the massive component (37) of the one-point function C ± 1,0 (t) (cf. Eq. (36) in the main text). In particular, the difference between the symmetric and antisymmetric correlation boils down to the sign of α 2 . It is easy to realize that this is equivalent to consider a different sign in the corresponding time t 2 (cf. Eq (57)). Collecting all these observations, the correlation of interest read

A.2 Four-point function in the slab
A similar calculation can be carried over for the the four-point function in the A-sector in Eq. (42). In the path integral formulation, the object of interest is where all fields and variables are defined above. The calculation works in the exact same way, the only difference being the number of operator insertions The relations (58) hold with now i = {1, 2, 3, 4}, and Eq. (62) corresponds to the special case Upon analytic continuation to real times, and taking all scales much larger than W , we get to the expression For C ± 2,0 (t), everything simplifies to

B Calculations in the massless sector
For computing the contribution from the massless sector (B), we rely on the approach of [34], based on Bogoliubov transformations. In [34], the author focuses on two-point correlation functions at equal times. When the two sound velocities are different (u 1 ≠ u 2 ), we end up in correlators at different times that we provide in what follows.

B.1 Two-point function at different times
Here we derive Eqs. (39) and (40) in the main text, that enter in the B-sector contribution to Eq. (36). This is just the two-point function for a quench K 0 → K f in a TLL (below the unique sound velocity is set to 1), for which we use the notations introduced in Section 3.2. We consider t 1 ≠ t 2 and, without losing generality, we take t 1 > t 2 .
We compute the general correlation with α ∈ R and, working in the Heisenberg picture, the expectation value is on the ground state of the Luttinger liquid hamiltonian with Luttinger parameter K 0 . The correlation in the exponent in Eq. (66) can be decomposed as where each of the terms above is a two-point function of θ at equal or different times. Then, using the following decomposition in modes for the field (in terms of the post-quench ladder and taking the thermodynamic limit (L → ∞) we get where B is the Bogoliubov matrix and we further defined Finally, in Eq. (69) we denoted as [⋅] ij the elements of a given matrix. Note that (69) is in general not real, however we will only be interested in real combinations of terms like in Eq. (67). For different times (t 1 ≠ t 2 ) one finds while, at equal times (t 1 = t 2 ≡ t), it simplifies to For t = 0, Eq. (73) simpifies to K f K 0 , so that the correlations like (69) only depends on K 0 as they should. Eq. (67) then reads where we defined Note that the leading term in (74) diverges as ∼ 1 p, giving rise to a power decay in (66).
In the case of equal spatial points, Eq. (74) simplifies to (cos(2pt 1 ) + cos(2pt 2 )) ± cos(p(t 1 + t 2 )) . (76) For K 0 = K f , we are computing a correlation function at equilibrium in the ground state. Accordingly, the expression above becomes time translational invariant (only the term involving the times difference survives).

B.2 Four-point function at different times
Since the theory is quadratic, the calculation of higher point correlation functions can always be reduced to that of two-point functions. We will see it explicitly below in the case of the four-point function considered in the main text in Eq. (44). We start by noting that which follows directly from Wick theorem. Then, we proceed by splitting the exponent in the rhs of (77) in three pieces as follows This splitting is particularly convenient because each term is infared finite, so that no cutoff is needed at small p. The first two terms in (78) are of the form (74) evaluated at equal times. Performing the integral (with an UV cutoff) we get For the last term in (78), we find Putting everything together and performing trivial algebraic simplifications, we get C Calculations in the exact state: Bogoliubov approach For comparison, we briefly sketch the calculations for the same correlation functions within the Bogoliubov approach. More details can be found in Ref. [116].

C.1 One-point function of θ ±
The one-point function in Eq. (35) can be written as with where used the decomposition and Hence, the expectation value in the exponent in (82) is = cos(u i pt) cos(u j pt)⟪θ i,p θ j,−p ⟫ γ + sin(u i pt) sin(u j pt)α i,p α j,−p ⟪n i,p n j,−p ⟫ γ .
The small p expansion of (86) reads with A ij (p) and B ij (p) regular for p → 0. The leading contribution in (82) comes from the term ∝ 1 p 2 . This contribution was explicitly computed in [116] and gives an exponential decay in (82). Now, we consider the next-to-leading contribution ∝ 1 p in (87). The explicit expression for B is where we defined and Π ij = lim p→0 p α i,p α j,−p ⟪n i,p (t)n j,−p ⟫ γ = Θ a 2 ij .

C.2 Two-point function of θ ±
We can similarly derive the two-point function (41), i.e.

D Calculations in the exact state: Coherent states
Some of the calculations in the main text are more easily done in the coherent states basis (in a path integral fashion), that we now review. To begin with, we consider a simple squeezed state of the form with W p ∈ C. Let us define the coherent states z p ⟩ as follows b p z p ⟩ = z p z p ⟩, z p ⟩ = e zpb † p −z * p bp 0⟩, (99) ⟨z p w p ⟩ = e − 1 2 zp 2 + wp 2 −2z * p wp , I = p dz p dz p π z⟩⟨z (100) where b ( †) p are operators, z p , w p ∈ C, and z⟩ = ⊗ p z p ⟩. The norm of ψ⟩ ⟨ψ ψ⟩ = p dz p dz p π ⟨ψ z⟩⟨z ψ⟩, is computed as follows. Using the definitions in (99), we find ⟨ψ z⟩ = p>0 ⟨0 e W p bpb−p z p , z −p ⟩ = p>0 e W p zpz−p e − 1 2 ( zp 2 + z−p 2 ) .
Moreover, using also that ⟨z ψ⟩ = ⟨ψ z⟩ * , the norm of ψ⟩ is written as a gaussian integral, which can be computed explictly Moving to correlation functions, and taking into account the above normalization, we similarly find where we used the commutation relations to make b p act on z⟩.
We then consider the squeezed state of interest in this work, namely of the form where in w ab = w p ab (a, b ∈ {A, S}) the p-dependence is implicit, and, for simplicity, we assumed them to be real. First we (re)define the coherent states as and we start again from computing the norm of ψ⟩. Repeating the same steps above, we find and ⟨z ψ⟩ = ⟨ψ z⟩ * . Using this, the norm of the state (105) where we defined the vectorẐ p = (z A,p , z B,p , z A,−p , z B,−p , z * A,p , z * B,p , z * A,−p , z * B,−p ), I 4 is a 4 × 4 identity matrix, andM p results in a 8 × 8 symmetric matrix. Expoiting its gaussian nature, the above integral can be evaluated analytically to get Similarly, correlation functions can be evaluated making use of the property of gaussian integral For example, with the definitions above, ⟨ψ b A,−q b † A,q ψ⟩ = A B,p can be collected in the following 4 × 4 matrix where m kl = −2M −1 p,kl , and we exploited the symmetries of B p . Expectation values are understood on the normalized state ψ⟩ N . In particular, we want to consider the O(p 0 ) of B p , for the squeezed state in (105) with W p = W p (cfr. Eq. (27) in the main text), so that expectation values are given by ⟪⋅⟫ γ .
Using the following definitions (analogous to (89) and (90) and by expanding θ A B,p and n A B,p in terms of b † A B,p , one can check that cos 2ϕ − 1 cos 2ϕ + 1 (114) namely they are independent on the value of γ. This is not the case for Π ab , in which case one finds Finally, note that τ B never enters in the above expressions.