Sudden interaction quench in the quantum sine-Gordon model

We study a sudden interaction quench in the weak-coupling regime of the quantum sine-Gordon model. The real time dynamics of the bosonic mode occupation numbers is calculated using the flow equation method. While we cannot prove results for the asymptotic long time limit, we can establish the existence of an extended regime in time where the mode occupation numbers relax to twice their equilibrium values. This factor two indicates a non-equilibrium distribution and is a universal feature of weak interaction quenches. The weak-coupling quantum sine-Gordon model therefore turns out to be on the borderline between thermalization and non-thermalization.


Introduction
The last years have witnessed increasing interest in the dynamics of isolated quantum manybody systems. As has occurred in many other fields that were considered mainly academic in the past, recent advances in experiments with ultracold gases confined to optical lattices have opened up the possibility of making real tests of the long-time evolution of essentially isolated quantum systems. Subsequently, a plethora of unanswered questions have regained the attention of the scientific community.
One important issue in this context is the question of thermalization in quantum many-body systems. Coupled to environments, the dynamics of quantum systems is known to lead to equilibrium states described by the thermal ensembles of quantum statistical mechanics. But once the system is isolated and initialized in a highly non-thermal state, it remains unclear whether the final, long-time state of the system can be described with one of those ensembles, i.e. if the system has thermalized. Rigorously, it is easy to show that quantum mechanics does not allow evolution from a pure state to a thermal distribution, as unitary time evolution preserves the purity. However, to what extent thermal averages reproduce long-time quantum averages, and the conditions whereby this happens, can still be considered an open problem [1]- [14].
From this debate, new questions arose concerning the role that integrability plays in the long-time evolution of quantum systems. The dynamics of integrable systems is expected to be very constrained due to the large number of constants of motion, and hence ordinary thermalization should not be present in such systems. Experimental studies of systems near integrability point in this direction [15]. However, it has been argued that integrable systems might still relax to an ensemble described by quantum statistical mechanics if all the constants of motion are included as constraints [16]. In this respect, the thermalization debate is still relevant even regarding integrable systems [17,18].
Similar to non-equilibrium classical statistical mechanics, the lack of a general framework to study non-equilibrium quantum problems makes it necessary to focus on specific models with 3 the aim of extracting general features from them. However, the long-time evolution of far-fromequilibrium systems is a very challenging topic, and suitable techniques, both analytical and numerical, have only recently been developed. In this paper, we employ the forward-backward scheme [21,22] based on the flow equation method [23,24], which has turned out to be a reliable and powerful approach to solve Heisenberg equations of motion for operators. The main idea is to use the flow equation method to diagonalize the Hamiltonian in a controlled approximation, and then to study the time evolution problem in this diagonal basis. Hence, real time evolution in this basis becomes simple and can be extended to long times without secular terms. All the difficulties of the problem are therefore encoded in the unitary transformation.
In this paper, we apply this approach to study the real time dynamics of the sine-Gordon model after a sudden quench of the interaction. Sudden interaction quenches are very interesting for studying non-equilibrium dynamics as they provide far-from-equilibrium initial states, but simplify the theoretical calculations because these initial states are simple. In optical lattices, it is also possible to implement them experimentally [19].
The quantum sine-Gordon model is a (1 + 1)-dimensional scalar field theory with a rich phase diagram. It is an integrable model, whose exact solution can be obtained by using the Bethe ansatz [20]. However, as usual, this exact solution does not guarantee a simple calculation of the observables. Therefore, it is interesting and necessary to implement approximate schemes like the forward-backward method. The choice of the sine-Gordon model is motivated by the fact that it is a paradigm for one-dimensional translation-invariant interacting systems, with mappings that connect it to many other models. For example, it arises as the effective description of a system of interacting fermions in one dimension with backscattering, or for spin-1/2 quantum spin chains [25]. Sudden quenches have already been studied for one-dimensional fermions with density-density interactions, i.e. in the context of the Luttinger model [18]: this is an integrable model whose exact solution is quadratic once it is written in terms of bosonic excitations [25]. Likewise, interaction quenches have been studied for the exactly solvable cases of the sine-Gordon model, namely, the Luther-Emery point, where the Hamiltonian can be expressed in terms of massive free fermions, and the semiclassical limit, where the mapping is to massive free bosons [26,27]. However, in general, such simple mappings are not possible for the sine-Gordon model: despite being integrable, it cannot be expressed as a quadratic Hamiltonian, and hence one can expect a redistribution of energy between different modes.
The paper is organized as follows. In section 2, we introduce the main features of the model, including its phase diagram, and describe its solution by the flow equation method. In section 3, we implement the forward-backward scheme to study the time evolution of the mode occupation operator in the weak-coupling region of the phase diagram. In section 4, the resulting expressions are employed to analyze the effect of a sudden quench of interactions: these are the central results of this work.

The model and the flow equation solution
The sine-Gordon model is a ubiquitous one-dimensional model widely studied in many different areas of physics. Its classical version became very popular in the 1970s as it has non-perturbative solutions known as solitons [28]. Here we will be interested in its quantized counterpart, namely the quantum sine-Gordon model. As mentioned before, the quantum sine-Gordon model is related to one-dimensional fermions with backscattering. There are many other similar mappings such as the one-dimensional Hubbard model near half-filling, the Coulomb gas problem, the two-dimensional classical X -Y model and quantum spin chains [25]. Therefore, the quantum sine-Gordon model is a natural and important setting for understanding quench dynamics.
The Hamiltonian of the model is defined as follows: where φ(x) is a scalar bosonic field and (x) is its conjugate momentum field. In order to impose the quantum structure, they must satisfy the commutation relations The Hamiltonian contains the parameter β and the coupling constant g, which define the phase diagram of the model. In the following, we will usually use the parameter α 2 = β 2 /4π, which turns out to be directly related to the scaling dimension of the cos-interaction term. The rest of the parameters are used to regularize the theory: a is a lattice discretization parameter (its inverse 1/a plays the role of an ultraviolet cutoff) and L is the system size (its inverse is the infrared cutoff).
A schematic picture of the phase diagram of the model is shown in figure 1. At α 2 = 1 (Luther-Emory point), the Hamiltonian can be mapped to the non-interacting Thirring model [29], whose relevant degrees of freedom are fermions [31], which can be identified with quantized solitons of the sine-Gordon equation. Away from this so-called Thirring line, these fermions experience an interaction: the region α 2 > 1 corresponds to repulsive quantum solitons, whereas in the region α 2 < 1 the interaction is attractive, which leads to bound states called breathers.
Another interesting point of the phase diagram occurs near α 2 = 2. Here the system undergoes a Kosterlitz-Thouless continuous phase transition, which can be understood from the renormalization group equations for the flowing coupling constants [30]: where = 1/ √ 2πa is the ultraviolet cutoff. A graphical solution of these equations is shown in figure 2. For α 2 < 2 the coupling constant g flows to strong coupling, which signals the opening of a gap in the spectrum. This corresponds to the emergence of massive fermionic solitons as the appropriate low-energy degrees of freedom for the model. For α 2 > 2 we have the weak-coupling regime, where the coupling constant g flows to zero, and the relevant degrees of freedom are massless bosons. In this region, an approximate solution of the equations for fixed initial parameter α 0 in the weak-coupling limit |g 0 | 1 is simply In this paper, we will be mainly interested in this so-defined weak-coupling limit where the flow of α 2 can be neglected to leading order.
Next, we express the fields in terms of their bosonic modes: The various creation and annihilation operators obey the usual bosonic commutation relations: We use the Mandelstam vertex operators [31] to construct the Hamiltonian: where : O := O − 0|O|0 means normal ordering of the operator O with respect to the noninteracting ground state [32]. The Hamiltonian now reads This Hamiltonian can be studied in its entire phase diagram using the flow equation approach [24] as shown in [33,34]. The key idea behind this method is to successively apply infinitesimal unitary transformations that eventually diagonalize the Hamiltonian. The stable choice of such a transformation sequence requires energy scale separation similar to conventional renormalization approaches: as the transformation progresses, more and more energy-diagonal interaction matrix elements are eliminated. Such a scheme was proposed by Wegner [23] and independently by Glazek and Wilson [35,36]. Wegner showed that a suitable infinitesimal transformation is obtained with the following canonical generator: where H 0 is the diagonal and H int is the interaction part of the Hamiltonian. B is the flow parameter that parametrizes the diagonalizing flow. It can be related to an energy scale B = 1/ √ B in analogy to the conventional RG scheme. The key difference is that the flow equation B corresponds to an energy difference that is being eliminated, whereas in a conventional RG scheme corresponds to the UV cutoff, which is an absolute energy scale.
The flow of the Hamiltonian is given by the following differential equation: The methodological challenges come from the implementation of such a transformation and the integration of the ensuing differential equations. In most cases, approximations are required in order to get closed sets of equations. However, these approximations do not necessarily match those employed in other methods such as perturbation theory, which e.g. opens up the possibility of accessing non-perturbative regimes using flow equations. For a comprehensive review of the flow equation method and its applications, we refer the reader to [24]. Before studying the non-equilibrium dynamics of the sine-Gordon model, we first briefly review the flow equation solution of the equilibrium sine-Gordon model in order to make this paper self-contained. More details of this calculation can be found in [33,34].
It turns out to be more convenient to work with Fourier transformed vertex operators: Some relevant properties of these operators are summarized in appendix A. In this representation, the generator of the unitary transformation consists of two parts: where The actual values for these coefficients are obtained by solving the flow equations for g(B) and α(B) [33,34]: with l ≡ − 1 2 log(32B/a 2 ). This flow successfully describes the different regions of the quantum sine-Gordon phase diagram, from the weak-coupling regime close to the Kosterlitz-Thouless transition, where excitations are massless bosons, to the Thirring line at α 2 = 1 with massive solitonic excitations.
Since we are interested in the weak-coupling limit (g 0 small for fixed initial α) in this paper, we can neglect the flow of α(B) to leading order and identify α(B) with its initial value α. The effective diagonal Hamiltonian generated in the limit B = ∞ then has the following form: Here are conveniently normalized Fourier transformed vertex operators. ω k (B = ∞) is given by [33,34]

Implementation of the forward-backward scheme
The flow equation approach to the quantum sine-Gordon model can be used to study the real time dynamics after a sudden interaction quench. The general idea proposed in [21] consists of three steps. First, the observable is transformed into the diagonal basis of the Hamiltonian by carrying out the same sequence of infinitesimal unitary transformations: This is the so-called forward transformation of the observable, and in most cases it implies a very complicated structure of O(B = ∞). This observation is familiar from exact Bethe ansatz solutions.
The advantage of working in the diagonal basis comes from the actual time evolution, which is much easier for a diagonal Hamiltonian. Later we will see that in our model the term H diag from equation (25) can be neglected compared to H 0 for not too long times. Hence, time evolution translates into phase factors and O(B = ∞, t) is obtained easily.
The final step is the backward transformation, where the flow of the transformed and timeevolved operator back to the original basis is carried out. The result O(t) is then an approximate solution to the Heisenberg equations of motion for the operator and it is straightforward to work out its expectation value with respect to the initial (non-equilibrium) state.
This sequence of transformations constitutes the flow equation forward-backward scheme. It has already been used to study the real time dynamics of a Fermi liquid after a sudden interaction quench [6,7] and the real time dynamics after a sudden quench in the ferromagnetic Kondo model [8,9]. One of the main advantages of this approach is that it avoids the infamous problem of secular terms in perturbation theory: in time-dependent perturbative expansions, these can restrict the perturbative solution to time scales shorter than [coupling constant] −1 . In this sense, the forward-backward scheme is the quantum version of unitary perturbation theory in classical mechanics.
In this paper, we are mainly interested in the time evolution of the bosonic number operator after a sudden interaction quench. Hence, our first goal is the implementation of the forward-backward scheme for the creation/annihilation operators.

Forward transformation
The forward transformation requires the solution of the flow equation with the initial condition a i,k (B = 0) = a i,k . Here i = l, r corresponds to left and right movers. In order to obtain closed equations, we make an ansatz: and likewise for the right movers. The ansatz is parametrized by various functions that must be calculated by working out the commutators in equation (30). In the weak-coupling phase, this ansatz becomes exact in the infrared (low-energy) limit since the coupling constant g(B) flows to zero. It is convenient to decompose the transformation in two stages, given by η (1) (B) and η (2) (B). The lowest order contribution from the first part of the generator yields the following flow equation: Now we focus on the weak-coupling limit. Since the coupling constant g(B) flows to zero, we can use it as a perturbative parameter in the flow equations. It can be shown that in order to preserve the bosonic commutation relations during the flow it is consistent to assume an expansion of the form h (l) l,k (B) = 1 + O(g 2 ) and h (l) r,k = 0 + O(g 2 ) (see appendices B and C).
In the weak-coupling limit, the flow of the coupling constants can be approximated by equations (5) and (6), which simplifies the integration of the differential equations. The result of carrying out the whole flow is where This result is valid for α 2 < 4, which is the region we are mainly interested in. Now let us discuss the effect of the second part of the generator η (2) (B). Due to the structure of this generator, a different approach is possible. The complete infinitesimal transformation can be rewritten as The term in brackets corresponds to the transformation already worked out above. The advantage of this expression arises from the fact that we already know the effect of the exponentiated η (2) (B) on the bosons and vertex operators: e η (2) (B) a l,k e −η (2) (B) = a l,k cosh(ψ(B)) + a † r,−k sinh(ψ(B)), e η (2) (B) a † r,−k e −η (2) (B) = a † r,−k cosh(ψ(B)) + a l,k sinh(ψ(B)), Hence the effect of the second part of the transformation is Carrying out the whole transformation in second order of the coupling constant yields Here we have used a decomposition derived in appendix B: p . We will later see that in the present order of the calculation the contribution coming from the second part of the transformation can be neglected since ψ(∞) ∝ g 2 0 . This is consistent with our assumption that we ignore the flow of α 2 (B) in the weak-coupling limit: neglecting the flow of α 2 (B) in fact just corresponds to neglecting the generator part η (2) .

Time evolution in the diagonal basis
The second step in the forward-backward scheme is the time evolution of the observable in the diagonal basis. Here, however, an additional approximation is required in order to solve the time evolution problem: in the diagonal Hamiltonian (23) only the bosonic kinetic term H 0 is taken into account. We will later see that this approximation implies a maximum time scale up to which our calculation can be trusted. The time evolution dictated by H 0 is straightforward due to the simple transformation of the vertex operators: Therefore, the time-evolved annihilation operator in the diagonal basis reads

Backward transformation
The final step of the transformation requires one to undo the flow equation transformation for the time-evolved operator (44). This is straightforward due to the perturbative nature of the transformation. We make a general ansatz for the operator: . The flow equations resemble those for the forward transformation. Since the contributions from the other functions in the ansatz follow directly via the bosonic commutation relation, we will only explicitly write down the flow equation for the function with the solution By applying the first part of the generator, the final time-evolved operator in second order of the renormalized coupling constant then reads This is the main technical result of our paper, which can be used as a building block to study the time evolution of all the other observables.

Consistency check: ground state energy in perturbation theory
As a consistency check for the flow equation calculation, we now compare the flow equation result for the ground state energy E (2) 0 with second-order perturbation theory. In fact,

12
it is sufficient to evaluate the kinetic energy E (2) K,0 in the ground state since one can easily prove E (2) K,0 = −E (2) 0 . Our goal is therefore to calculate where |0 is the ground state of the interacting model. Within the flow equation formalism, this is most conveniently evaluated in the diagonal basis with the forward transformed operators where we have used 0 |n i,k |0 = 0|n i,k (B = ∞)|0 . Here |0 is the bosonic vacuum since this is trivially the ground state in the diagonal basis. If this calculation is carried out using our previous results from the flow equation formalism, the result actually does not coincide with perturbation theory. However, this is simply due to the fact that the flow equation calculation is a renormalized expansion, whereas conventional perturbation does not contain renormalization effects. In order to compare with perturbation theory 5 , we therefore artificially set g(B) = g 0 and find .
The same result can be obtained by working out the kinetic energy in second-order perturbation theory: | n|H I |0 | 2 e −λ k k(n l,k +n r,k ) , where |n ≡ |n l,k 1 , n l,k 2 , . . . , n r,k 1 , n r,k 2 , . . . . By using the matrix elements of vertex operators given in appendix A, we arrive at the following expression: where we insert the time-evolved operators in the Heisenberg picture (48). It is this conceptual simplicity that makes sudden interaction quenches very appealing for studying non-equilibrium problems. Notice that an identical result to (55) can be obtained for right movers, and therefore we restrict ourselves to explicit expressions for left movers only. An expression for (55) can be worked out readily from result (48) obtained within the forward-backward scheme: with the integral In order to obtain this result, we have made use of the properties of vertex operators summarized in appendix A and of the expression for (l) k ( p; ∞) given in (34). However, it is more convenient to express this result in terms of the equilibrium occupation numbers. Fortunately, they follow directly from the flow equation calculation: Here |0 denotes the interacting ground state, and This results in the following compact expression for the ratio of non-equilibrium to equilibrium occupation numbers: n l,k (t) n l,k EQ = I (ka, (t/a)) I EQ (ka) and likewise for right movers. The integrals in this expression must be computed numerically. Plots of these ratios are shown in figure 3 for different values of α 2 and fixed momentum, and in figure 4 for fixed α 2 and different momenta ka. One observes that the key phenomena after the quench are damped oscillations of the mode occupation on a time scale set by the lattice cutoff a. The asymptotic value of the mode occupation universally converges to twice its equilibrium value: This can be understood easily by noticing that I EQ (ka) and I (ka, t/a) only differ by replacing a factor 1/4 by sin 2 (xt/a) in the integrand. Clearly, the limit n l,k (t → ∞) just amounts to taking the time average over sin 2 (xt/a) in the integrand, which gives 1/2 and therefore I (ka, t/a → ∞) = 2 I EQ (ka).
Note that (60) implies a non-thermal mode distribution function for the asymptotic state of this closed quantum system: the equilibrium system with nonzero temperature cannot reproduce this expression.

Conclusions
We have demonstrated that an interaction quench in the weak-coupling phase of the sine-Gordon model leads to interesting dynamics that is quite different from a quench of the forward scattering only. Quenching the forward scattering does not induce any dynamics for the bosonic occupation numbers, which are in fact constants of motion in this case [18]. On the other hand, we have found that a quench of the backscattering term in the weak-coupling phase leads to a real time dynamics where the excitation energy of the quench is converted into bosonic mode occupations that oscillate on a time scale set by the ultraviolet cutoff and eventually reach twice their equilibrium values (60). This factor 2 is universal for weak interaction quenches and has previously been seen for the momentum distribution function in the non-equilibrium Hubbard model [6,7] and for the magnetization of the non-equilibrium ferromagnetic Kondo model [8,9]. It occurs for quenches in quantum systems where (i) second-order perturbation theory is valid (at least up to a certain time scale) and (ii) for observables such as the mode occupation number operator that commute with H 0 (for a proof and more details on the conditions, see [7]).
Since such general results are important for developing a better understanding of nonequilibrium quantum many-body systems in general, we need to critically re-examine the approximations in our calculation. Due to the weak-coupling behavior of the running coupling constant, the second-order calculation presented here becomes more and more reliable in the infrared limit. Therefore, higher order corrections to the universal factor 2 will vanish in the low-energy limit. However, we had to make the additional approximation to neglect the time evolution generated by H diag (∞) from (25). Therefore, result (60) can only be trusted up to the time scale Whether our central result really holds beyond the time scale τ k or only in a time window a t τ k cannot be answered based on our calculation. The interaction quench in the weak-coupling phase of the sine-Gordon model therefore occupies an interesting place between the ferromagnetic Kondo model, where the factor 2 is in fact asymptotically exact [8,9], and the non-equilibrium Hubbard model in d 2 dimensions, where the factor 2 describes the prethermalization regime [6,7] before the system eventually thermalizes. The integrability of the sine-Gordon model could make one expect that the non-equilibrium distribution function (62) remains stable for all times and does not approach a thermal limit form. However, the conserved quantities in the sine-Gordon model do not impose any obvious constraints on the dynamics of the momentum distribution function, unlike in the case of quenching the forward scattering [18].
Further studies of this question would be very worthwhile, either numerical or based on an exact solution. We have seen that the weak-coupling quench in the sine-Gordon model is on the borderline between thermalization and non-thermalization seen through the eyes of the mode distribution function. This is similar to the role played by the celebrated Fermi-Pasta-Ulam problem for classical many-body systems [37,38]. A better understanding of the weak-coupling quench in the quantum sine-Gordon model could be an important step in elucidating the fundamental question of thermalization in the quantum world.

Appendix A. Properties of vertex operators
We summarize some important properties of vertex operators (k, k > 0). More details can be found in [33]. Using some straightforward algebra, one can verify that this is indeed fulfilled for the flow equations derived in this paper.