Quantum quench, where parameters in the Hamiltonian are changed discontinuously to leave a system out of equilibrium, is an important subject in the study of nonequilibrium dynamics of isolated quantum systems [1, 2]. Though the problem is well defined, most of the known exact results essentially come from noninteracting systems, e.g. the transverse-field Ising model [3, 4], and some pairing models [5–7]. For many-body systems the situation is much more complicated, and generally completely exact treatments are far beyond reach even for relatively simple models [8–10]. In order to procure some qualitative understandings and to make decent quantitative predictions, we usually have to resort to effective theories.

Analogous to the prevalent Landau's Fermi liquid theory, a bunch of one-dimensional (1D) quantum systems can be classified into a universal category known as Luttinger liquids [11–13], where only a handful of parameters are sufficient to address the low-temperature properties of systems effectively. A bit of a surprise is that within the framework of the Luttinger liquid theory, it is possible to produce highly accurate descriptions on the nonequilibrium dynamics of the systems in their gapless (critical) phases. This observation has been justified by numerical computations on various lattice models [14, 15].

One advantage of the investigations from the effective theory side is that analytical approaches, in many cases, can be carried out at least approximately [16, 17]. In the work of Dóra et al [18], based on the analysis of an effective Luttinger Hamiltonian, an interesting observation was made, which unveiled the correspondence between the time-evolved entanglement spectrum and the nonequilibrium characteristics in the quench dynamics, including the identification of the largest momentum-space entanglement eigenvalue and the Loschmidt echo. In their work, the correspondence was made plain first via the calculations on the simple Luttinger liquid Hamiltonian, and later was ratified by the numerical results on spin-chain lattice models. In our opinion, the work highlights the effectiveness of the Luttinger liquid theory in the nonequilibrium dynamics of 1D critical systems.

Nevertheless, in the paper [18] the authors discussed primarily the noninteracting-to-interacting quenches, leaving the issue of general quenches untouched, e.g. a quench from one nonzero coupling strength to another. In this study, we would like to address the quench dynamics of the Luttinger liquid model in more general settings, by investigating the quenches that can start from an arbitrary (interacting) ground state (GS). With the assistance of su(1, 1) algebra, we procure the exact expression of the time-evolved wavefunction, and obtain the largest entanglement eigenvalue in momentum space accordingly. We note that the correspondence raised in [18] can be generalized consequently. We confirm it numerically in the XXZ-related model (see (31) for the Hamiltonian), which is a lattice realization of the Luttinger Hamiltonian (1). There we find that the agreement holds in the sense that the largest momentum-space entanglement eigenvalue coincides with the overlap of the time-evolved wavefunction with its noninteracting GS, rather than with the initial interacting GS from which the quench starts.

The rest of the paper is organized as follows. After an introduction of the Luttinger liquid Hamiltonian in terms of paired bosons in section 2, we discuss the general g₀ → g quench and calculate the aftermath time-evolved wavefunction in section 3. We note specifically that the expression of the time-evolved wavefunction depends on the choice of Bogoliubov quasiparticles. The correspondence between the wavefunction overlaps, including the Loschmidt echo, and the largest entanglement eigenvalues, is examined in section 4. Numerical results on the fermionic lattice model are presented in section 5. The discussion and conclusion is given in section 6.

As a description of a 1D Luttinger liquid, the Hamiltonian we consider here takes the following form [18]:where ${b}_{q}\ ({b}_{q}^{\dagger })$ is the annihilation (creation) operator of a bosonic particle with momentum q. For Luttinger liquids, the dispersion relation and the pairing strength are given as [12]

The above Hamiltonian can be diagonalized via the Bogoliubov transformation:where the condition of ${\left|u\right|}^{2}-{\left|v\right|}^{2}=1$ is imposed to preserve the bosonic commutations for the newly introduced a-operators. Since ω(q) and g(q) are real functions, u, v can be set to be real and parameterized by θ(q) asUnder the above Bogoliubov transformation, equation (4), (the condition of $\left|g\right|\lt \left|\omega \right|$ needs to be met in order to guarantee the existence of a solution), the pairing Hamiltonian, equation (1), can be transformed into a diagonal form in terms of the Bogoliubov operators a_q:where $\epsilon (q)=\sqrt{{\omega }^{2}(q)-{g}^{2}(q)}$. We note that for ω(q) and g(q) given by equation (2), θ is actually q independent. For simplicity, in the new Hamiltonian in equation (5) the constant term − ∑_q>0ω_q has been dropped compared to the original Hamiltonian, equation (1).

The Bogoliubov Hamiltonian, equation (5) indicates that GS of the system should be the vacuum state of a-operators, i.e. a_q∣GS⟩ = 0 is satisfied for each q. It is easy to verified that ∣GS⟩ can be expressed as a pairing state of the constituent particles as:where $1/\left|u(q)\right|$ serves as the normalized factor. Throughout this paper, ∣0⟩ is reserved to specify the vacuum state of the original b-operator, i.e. b_q∣0⟩ = 0.

The issue we are going to address is the quench dynamics on the pairing strength g: starting from an interacting (entanglement) GS with respect to some (nonzero) value g₀, we switch the interaction strength abruptly from g₀ to another g value. What is the post-quench dynamics of the system?

For the Hamiltonian considered here, as mentioned before, the explicit form of the time-evolved wavefunction for the quench from the noninteracting limit (i.e. g = 0) to a finite g has been obtained, together with a discussion of the correspondence between the Loschmidt echo and the entanglement in momentum space [18, 19]. Here in this section by utilizing su(1, 1) Lie algebra, we make a further step by calculating the post-quench wavefunction for the general case of ∀ g₀ → g (under the condition of $\left|{g}_{0}\right|$ and $\left|g\right|\lt \left|\omega \right|$).

The GS for the pairing strength g₀, denoted as ∣GS₀⟩, is already given by equation (6). If we set the quench moment t = 0, thenwhere λ₀ := v₀/u₀, and the Bogoliubov coefficients u₀, v₀ are determined by equation (4) for the given g₀. Since (q, − q)-pairs with different q (q > 0) in the Hamiltonian, equation(1), are indeed decoupled from each other, for notation simplicity we would like to drop the label of q in the formalism as long as there is no misunderstanding. The summation or multiplication over q ( > 0) can be restored easily if necessary. For instance, ${b}^{\dagger }{\underline{b}}^{\dagger }$ usually stands for ${\sum }_{q\gt 0}{b}_{q}^{\dagger }{b}_{-q}^{\dagger }$, as in the case of equation (7).

For the quench g₀ → g at t = 0, the later-time wavefunction of the system is given by ∣$\Psi$(t)⟩ = U(t)∣$\Psi$(0)⟩, with the time-evolution operator U(t) = e^−iH(g)t (ℏ = 1). Note that H(g) can be expressed in a diagonal form by its corresponding Bogoliubov quasiparticles (see equation (5)), and consequently, $U(t)=\exp [-{\rm{i}}\epsilon t({a}^{\dagger }a+\underline{a}{\underline{a}}^{\dagger })]$. In this situation we have two sets of Bogoliubov quasiparticle operators, i.e. a₀ and a, corresponding to H(g₀) and H(g), respectively. It can be shown that they are related by a composition of two successive Bogoliubov transformations (see equation (3)):where $\tilde{u}=\cosh \tilde{\theta }/2$, $\tilde{v}=\sinh \tilde{\theta }/2$, and $\tilde{\theta }=\theta -{\theta }_{0}$, with θ and θ₀ the parameters for the Bogoliubov transformations, equation (4), for H(g) and H(g₀), respectively.

Equipped with the relation of a ↔ a₀ from equation (8), we can write down an alternative form of U(t) in terms of a₀ operators:where we have introduced K-operators asand the coefficients in the exponents read Note that the three Ks satisfy su(1, 1) algebra [20]:By employing the results of equation (A.3) in appendix A, the time-evolution operator U(t) in equation (9) can be factorized into a product form aswhere particularly,Then the normalized post-quench time-evolved wavefunction in terms of the a₀-Bogoliubov quasiparticles can be written asHere for simplicity we have dropped the phase factor owing to the action of ${{\rm{e}}}^{{C}_{0}{K}_{0}}$ on ∣GS₀⟩, since it is not essential in our discussions.

If we are instead interested in the expression of ∣$\Psi$(t)⟩ in terms of the original creation (annihilation) operator b^† (b) with respect to its vacuum ∣0⟩, we can employ the following procedure. First note that according to equation (B.8) in appendix B, the GS of H(g₀) can be also expressed by an a₀-pair generating operator acting on the vacuum ∣0⟩:where K₊ is defined in equation (10). Combined with the result in equation (15), we haveThe next step is to express the exponential term of ${K}_{+}(={a}_{0}^{\dagger }{\underline{a}}_{0}^{\dagger })$ in terms of b-operators, which is done in equation (B.6) in appendix B. By replacing τ → C₊ − λ₀, and u → u₀, v → v₀ there, we finally getwhere we have put back the normalization factor explicitly, andwith λ₀ = v₀/u₀.

Up to now we have obtained two equivalent forms for the same post-quench wavefunction: one is written in terms of an a₀-operator with respect to its vacuum ∣GS₀⟩ (see equation (15)), the other in terms of a b-operator with respect to the vacuum ∣0⟩ (see equation (18)).

A few comments on C₊(t) could be worth while. First we note that C₊(t) in equation (14) can also be written aswhere the definition of $\tilde{u}(t)$ and $\tilde{v}(t)$ can be found in equation (C.2) in appendix C, with the only difference of θ being replaced by $\tilde{\theta }$. In particular, when g₀ = 0, we have u₀ = 1, v₀ = 0, θ₀ = 0, which leads to $\tilde{\theta }=\theta $, and $\tilde{u}(t)\to u(t)$, $\tilde{v}(t)\to v(t)$. Then we return to the familiar case of 0 → g quench, which has been addressed in [18] in detail. Actually equations (15) and (20) can be regarded as a direct application of the result in 0 → g quench to the general g₀ → g quench, under the correspondence b → a₀, ∣0⟩ → ∣GS₀⟩, together with the transformation given by equation (8).

Secondly, we see that from equation (14), the quantity C₊(t) depends not only on g, but also on g₀ as well, since $\tilde{\theta }=\theta -{\theta }_{0}$. In contrast, where λ(t) in equation (19) is concerned, the dependence of g and g₀ there can be separated in a more explicit way. We can easily verify thatwhere u(t) and v(t) are given by equation (C.2). They appear as coefficients in the equation of motion for the b-operator (see appendix C), depending solely on the after-quench Hamiltonian H(g). For u₀, v₀, they are determined by g₀ via equation (4).

From the above analysis, we see that starting from a pairing state with the form of equation (7), when the system is subject to a g-quench, the aftermath time-evolved wavefunction still adapts a simple pairing form. Depending whether the a₀-pairing or the b-pairing picture is chosen, the form of the wavefunction can be either equations (15) or (18). Indeed an observation can be simply made: for any given paired state that is the null state of a given class of Bogoliubov quasiparticle annihilation operators, it can be equally expressed as a paired state in terms of another class of Bogoliubov quasiparticle operators, with respect to their own corresponding vacuum state. Different classes of Bogoliubov operators are related via SU(1,1) (Bogoliubov) transformations as in equation (3).

One of the motivations of the present work for considering the general g₀ → g quench is to examine an interesting observation made by Dóra et al in [18], which has been mentioned in section 1. In that paper it has been shown that for the 0 → g quench about the Hamiltonian, equation (1), the largest entanglement eigenvalue in the momentum space is always identical to the Loschmidt echo. The observation has been confirmed on various lattice models to which the Luttinger liquid theory can be applied, including the prestigious XXZ model in its gapless phase [18].

Here for the general g₀ → g quench, the Loschmidt echo (or the return probability) can be defined as usual as the overlap of the initial GS wavefunction and the final state wavefunction [19, 21], i.e.From the expression of ∣$\Psi$(t)⟩ in equation (15), we havesince ∣$\Psi$(0)⟩ = ∣GS₀⟩. Or more precisely, ${ \mathcal L }(t)={\prod }_{q\gt 0}\left(1-{\left|{C}_{+}(t;q)\right|}^{2}\right)$ if the q dependence is restored. On the other hand, analogous to the discussion in [18], it is straightforward to show that the largest eigenvalue of the reduced density matrix ρ_A(t) for the right-movers of a₀-Bogoliubov quasiparticles, denoted as ${P}_{\max }^{(1)}$ here, readswhere $\tilde{u}(t)$ has been introduced in equation (20). It coincides with ${ \mathcal L }(t)$ in equation (23), i.e.Thus similar to the 0 → g quench, the correspondence between the Loschmidt echo and the largest entanglement eigenvalue in the momentum space is established.

However when it comes to a numerical simulation on a lattice model, it is the Fock space of the constituent particles that we are going to work with, whose vacuum is the empty state of these particles by definition. It remains open whether the low-lying entanglement spectrum obtained in this Fock space after tracing out, say, the q < 0 partition is identical to the one obtained from the effective Hamiltonian, equation (1). Indeed, from the results in the previous section, we note that the reduced density matrix ρ_A(t) and the subsequent entanglement spectrum can be evaluated by either the wavefunction in equation (15) or equation (18), where the vacuum is the null state of either a₀ operators or b operators.

We have addressed the situation with the wavefunction in equation (15) so far. For the other case, i.e. based on the wavefunction in equation (18), the entanglement between the left- and right-moving elemental excitations with respect to the b operators can be evaluated as before. The largest entanglement eigenvalue readswhere the expressions of u(t) and v(t) (with respect to the g value) can be found in equation (C.2), and we note that here different from ${P}_{\max }^{(1)}$ in equation (25),The relation between ${P}_{\max }^{(1)}$ and ${P}_{\max }^{(2)}$ can be made transparent by employing equation (19) to give

The above analysis shows that in the general g₀ → g quench, the correspondence between the Loschmidt echo defined in equation (22) and the entanglement eigenvalue holds, only when the entanglement of the left- and right-moving elementary excitations is measured with respect to the initial GS of H(g₀), i.e. ∣GS₀⟩. Here the excitations are nothing but the Bogoliubov a₀-quasiparticles. It is worthwhile asking which entanglement is indeed being evaluated in the numerical simulations on lattice models, and consequently to what kind of wavefunction overlap the largest entanglement eigenvalue corresponds. We will address this issue in the next section by examining a spinless interacting fermionic lattice model, which is closely related to the prestigious XXZ model.

The Hamiltonian of the S = 1/2 XXZ model, which describes the Heisenberg-type interaction between neighboring 1/2-spins along a chain, readswhere ${S}_{j}^{\alpha }{\rm{s}}$ are spin operators at site j. Employing the Jordan–Wigner transformation that maps the spin operators into a set of fermionic ones with extra nonlocal string terms, the spin Hamiltonian can be recast into a system of spinless interacting fermions:We see that the XX terms in the spin Hamiltonian are simply transformed into the hopping terms in the fermionic Hamiltonian. Meanwhile, the J_z-terms produce the scattering terms.

In the momentum space the Hamiltonian, equation (30) readswith a constant term J_zN/4 (N is the lattice size) being dropped. We may define a parameter Δ := J_z/J, which measures the strength of the uniaxial anisotropy along the z direction. In the following discussions, we set J = 1. The energy and time are measured in units of J and J⁻¹, respectively. It is well known that the XXZ model is critical (gapless) in the interval $\left|{\rm{\Delta }}\right|\leqslant 1$, with a spin-singlet GS. Correspondingly, the GS of the Hamiltonian, equation (31), in this regime takes place at the half filling, and its low-energy physics can be well described by the Luttinger liquid, e.g. see [12, 22, 23].

The way from the fermionic Hamiltonian of equation (31) (abbreviated as H_F) to the paired bosonic Hamiltonian of equation (1) (abbreviated as H_B) can be established via the bosonization technique [13]. Without going into technical details, we would like to emphasize that here the umklapp term is dropped since it is irrelevant in the gapless phase, and the g₂ term, which describes the forward scattering processes on both sides of the Fermi surface, is kept explicitly. In the Hamiltonian, equation (1), it appears as the pairings of the bosonic b operators with opposite momenta [12]. The constant quantity g in equation (2) is indeed identical to g₂. The Luttinger parameter reads $K=\sqrt{(v-g)/(v+g)}$ in H_B [18], the precise value of which in the XXZ model is given by $K=\pi /2(\pi -\arccos {\rm{\Delta }})$ via the Bethe ansatz. K ∈ [1/2, + ∞) in the gapless regime $\left|{\rm{\Delta }}\right|\leqslant 1$, and the noninteracting point K = 1 corresponds to the XY limit (Δ = 0).

In the numerical simulation, we work directly on the interacting fermionic Hamiltonian, equation (31), at half filling in momentum space by the exact diagonalization (ED). The time-evolved wavefunction ∣$\Psi$(t)⟩ is obtained by the standard time-dependent Lanczos method [24]. The reduced density matrix ρ_A(t) with respect to ∣$\Psi$(t)⟩ and the consequent entanglement spectrum are procured accordingly, after tracing out the q < 0 part in the momentum space. The results for the quench of (Δ₀ = − 0.2) → (Δ = 0.4) on various lattice sizes (N = 18, 22, 26, 30) are presented in figure 1. In the main figure, the logarithms of both ${\left|\langle 0| {\rm{\Psi }}(t)\rangle \right|}^{2}$ and the largest entanglement eigenvalue ${P}_{\max }(t)$ are shown. The inset displays a comparison with that of ${\left|\langle {\mathrm{GS}}_{0}| {\rm{\Psi }}(t)\rangle \right|}^{2}$ for N= 30. Here following the previous conventions, ∣0⟩ is the noninteracting GS (at Δ = 0), and ∣GS₀⟩ is the GS at Δ₀. The quench takes place at t = 0, with ∣$\Psi$(0)⟩ = ∣GS₀⟩. For a given N the agreement between ${P}_{\max }(t)$ and ${\left|\langle 0| {\rm{\Psi }}(t)\rangle \right|}^{2}$ holds well during the time evolution. On the other hand, the deviation from ${\left|\langle {\mathrm{GS}}_{0}| {\rm{\Psi }}(t)\rangle \right|}^{2}$ is evident, as shown in the inset of figure 1. It indicates it is ${P}_{\max }^{(2)}$ in H_B, rather than ${P}_{\max }^{(1)}$, that is compatible with ${P}_{\max }$ which we have calculated in H_F. Recall that P⁽²⁾ is about the momentum-space entanglement in terms of b operators. Additionally, we observe that following a rapid decrease in the early stage, the later-time oscillation of these time-dependent quantities, which is due to the discreteness of the energy spectrum in finite systems, is evidently suppressed with increasing lattice size. We note that the short transient time before saturation taking place is around 4 J⁻¹ from the numerical results, consistent with the estimation from the effective field theory [18].

The above result can be understood as a natural ramification of the H_F ↔ H_B correspondence. Let us first consider the equilibrium case and take a look at the simple noninteracting case. It means Δ(:=J_z/J) = 0 in H_F, or equivalently, g = 0 in H_B (K = 1). In either situation, we have free particles with momenta as good quantum numbers. Both ground states are disentangled in momentum space, with only a single nonzero entanglement eigenvalue P = 1. As the interaction is turned on, which means Δ ≠ 0 on the fermionic side, we can expect that the entanglement between the q > 0 and q < 0 partitions will develop due to the interaction, yet it may not be easy to write down the wavefunction explicitly to figure it out. On the other hand, compared to the rather complicated situation on the fermionic side, the momentum-space entanglement can be much more easily estimated from the GS wavefunction of the bosonic counterpart (see equation (6)), where the entanglement between (q, − q)-pairs with respect to the b-operators can be readily read out.

If we agree with the effectiveness of the Luttinger liquid description of H_F [14, 15], we can expect that in the nonequilibrium case, the time evolution of the entanglement in the general Δ-quench over H_F should follow the features of its bosonic partner, as long as the Δ ↔ g relation is sustained. Thus, the arbitrariness in the entanglement calculation for the g₀ → g quench in H_B can be removed: we need to stick to the b-representation consistently and adopt the time-evolved wavefunction in the form of equation (18), in order to make reasonable a comparative study with its fermionic counterpart. Specifically in this section, we have ratified the speculation on the relation between the time-dependent largest entanglement eigenvalue and the overlap of wavefunctions for general quenches, which is drawn from the analysis on the bosonic case (see equation (27)), for the interacting fermionic lattice model.

In this work, we first consider a paired bosonic model, which serves as an effective Hamiltonian for spinless Luttinger liquids. By employing su(1, 1) algebra, we derive the exact form of the time-evolved wavefunction under general interaction quenches on pairing strength. The method we employ and the results are general as long as the pairing takes place in the individual (q, − q) channels. The application to the fermionic case is straightforward, with only the ${\mathsf{su}}(1,1)$ algebra being replaced by ${\mathsf{su}}(2).$

Based on the analytic form of the wavefunction, we re-examine under the general quench setup, the correspondence between the entanglement spectrum and the wavefunction overlap, including the Loschmidt echo. We note that by employing different sets of Bogoliubov quasiparticles, the expression for a given time-evolved post-quench quantum state can be different. It leads to an ambiguity in measuring the entanglement of paired bosons. However, by working numerically on an interacting spinless fermionic model that is closely related to the XXZ spin chain, we find that the ambiguity can be removed if the concurrence between the fermionic and bosonic descriptions is required. More specifically, we confirm that on the fermionic lattice model, the time-dependent largest entanglement eigenvalue in the momentum space coincides with the overlap between the time-evolved wavefunction and its noninteracting GS, rather than the initial state when the quench starts.

It is well known that quantum entanglement plays a key role in understanding the nature of the quantum world, especially of many-body systems [25, 26]. The entanglement entropy and entanglement spectrum have been routinely investigated in the study. Up to the present most studies predominantly focus on real or orbital space entanglement. For many gapped systems, the congruity between the energy spectrum of the edge states and the entanglement spectrum obtained in this manner has been established [27–29]. Nevertheless for critical (gapless) systems, due to the absence of an energy gap, the intrinsic entanglement might not be evidenced by cuttings in real space, which is corroborated by the observation of the slow decay of the eigenvalue spectrum of the reduced density matrix in the real-space density-matrix renormalization group (DMRG) [30].

Consequently, the issue of the entanglement structure in momentum space in gapless systems was raised and explored, first for spin chains [31–34], then for other correlated systems, including, for example, the Hubbard model [35] and the disordered fermionic systems [36, 37]. In the framework of the Luttinger liquid theory, where the elementary excitations fall into left- and right-moving (weakly) interacting particles, it is particularly natural to consider partitioning the many-body Hilbert space by momentum. It is quite significant that the correspondence between the largest momentum-space entanglement eigenvalue and the wavefunction overlap in nonequilibrium dynamics, which is rather plain in the effective bosonic model, has been verified by the numerical observations on interacting lattice models. It suggests that the entanglement in momentum space can furnish a new perspective for the study of nonequilibrium dynamics of quantum critical many-body systems.

Furthermore, from the perspective of the development of numerical algorithms, the investigation on the momentum-space entanglement can help to improve the performance and identify the most fruit-bearing applications of the (time-dependent) momentum-space DMRG [30, 38]. For example, there have been reports on the persistence of the entanglement gap in momentum space beyond the physically critical points, and the saturation of the entanglement entropy in nonequilibrium dynamics [18, 33, 37]. The DMRG working in the momentum space might have the potential to outperform its traditional real-space counterpart in these regimes.

Before concluding this section, we would like to mention that as far as the correspondence in quench dynamics is concerned, there may be some interesting issues remaining for investigation, which are beyond the scope of the present work. For instance, we notice in the numerical simulations that the correspondence we have discussed largely holds even beyond the regime of the gapless phase where the Luttinger liquid theory works well. To be precise, for the quenches crossing through the Berezinsky–Kosterlitz–Thouless (BKT) critical point (Δ = 1), which separates the gapless phase from the Ising antiferromagnetic gapped phase, the correspondence sustains up to some critical value of Δ > 1 (e.g. around Δ ∼ 1.6 for Δ₀ = 0.5). In contrast, for the quenches moving in the opposite direction, i.e. from the critical phase to the gapped ferromagnetic phase across the first-order-like critical point (Δ = − 1), the behaviors are quite different, and the agreement quickly dissolves. To study and understand the extent of the validity of the predictions from the Luttinger liquid theory will be interesting.

H.L. acknowledges support from the National Natural Science Foundation of China (Grant Nos. 11474136, 11874187 and 12047501).

In this appendix, we show that the operator U(τ) with the formcan be expressed aswhere K₀ and K_± satisfy su(1, 1) Lie algebra (see equation (12)) [39], and x, y, and z are constants. The coefficients in the exponents in equation (A.2), which satisfy the initial condition of C₊(0) = C₋(0) = C₀(0) = 0, are given bywhere r_1,2 are the roots of a quadratic equation of t: yt² + zt + x = 0. Note that the solution, though quite general, is not applicable in all situations, e.g. it fails when r₁ = r₂. We will discuss a specific example for this case in appendix B.

Following the procedure detailed in [20], the derivation of the above result is quite straightforward. First note that the operator U(τ) satisfies a differential equation asSubstituting the second form of U(τ) from equation (A.2) into the equation, and using the commutation relations of Ks, we can obtain the differential equations for C-coefficients asThe above equations, although nonlinear, can be integrated analytically. For instance, the equation for C₊ can be written aswhich can be solved by separation of variables: ${\rm{d}}{C}_{+}/(x+{{yC}}_{+}^{2}+{{zC}}_{+})={\rm{d}}\tau $. In the case of the quadratic equation yt² + zt + x = 0 having distinct roots, we get the solution in equation (A.3).

In this section as a supplement to the discussion in appendix A, let us consider an exponential form of the pair creation operators of the Bogoliubov quasiparticles, $U(\tau ):= \exp (\tau {a}^{\dagger }{\underline{a}}^{\dagger })$. Note that ${a}^{\dagger }={{ub}}^{\dagger }+v\underline{b}$ (see equation (3), and we also confine ourselves to real u, v), we havewithThe K operators, as a realization of su(1, 1) Lie algebra, are given byThough we cannot apply directly the result in equation (A.3) since here r₁ = r₂ = − u/v, the coefficients inequation (A.2) can still be solved as follows. From equation (A.6), we haveThenandWe obtain consequently

As an application of equation (B.6), let us consider the GS of the Hamiltonian, equation (1), which, in our simplified notation, readsWe say that the GS can also be expressed asThe conclusion can be verified in a straight way by means of replacing τ in equation (B.6) with − v/u, and using the fact that ∣0⟩ is the vacuum state of b, i.e. b∣0⟩ = 0.

At the end of this section, an extra comment might be interesting. Note that equation (B.8) can be rewritten asHere we see that there is some kind of duality between equations (B.7) and (B.9). For the vacuum of the constituent particles, i.e. ∣0⟩, we can equally regard it as a kind of pairing state of Bogoliubov quasiparticles as well, with respect to its own vacuum state.

Consider the operator b in the Heisenberg picture b(t) = e^iHt b e^−iHt, with H given by equation (1). ThenwhereIn the above derivation, the Hamiltonian in another form of equation (5) and the Bogoliubov transformation, equation (3) (and equation (4)) have been used.

We note that with u(t), v(t) defined above, the time-evolved wavefunction for the 0 → g quench can be written as ${\rm{\Psi }}(t)\sim \frac{1}{\left|{u}^{* }(t)\right|}{{\rm{e}}}^{\tfrac{{v}^{* }(t)}{{u}^{* }(t)}{b}^{\dagger }{\underline{b}}^{\dagger }}| 0 \rangle $ [18]. It suggests that in order to obtain ∣$\Psi$(t)⟩, there could be some convenient way by employing the time-dependent features of b(t). This is true actually. In the following we show that by assuming $| {\rm{\Psi }}(t) \rangle =\gamma (t){{\rm{e}}}^{\lambda (t){b}^{\dagger }{\underline{b}}^{\dagger }}| 0 \rangle $, we can quickly obtain the coefficient λ(t). First from equation (C.1), we haveOn the other hand, from the definition of b(t), we haveThe comparison of these two equations producesActually by employing the same strategy with the same assumption about ∣$\Psi$(t)⟩, we can show that for the general g₀ → g quench, the coefficient λ(t) in the exponent iswhere λ₀ := v₀/u₀, as for equation (7). It is identical to the result of equation (21) in the main text, which has been derived in a somewhat intricate but rigorous way by employing su(1, 1) algebra.