Persistence of power-law correlations in nonequilibrium steady states of gapped quantum spin chains

The existence of quasi-long range order is demonstrated in nonequilibrium steady states in isotropic $XY$ spin chains including of two types of additional terms that each generate a gap in the energy spectrum. The system is driven out of equilibrium by initializing a domain-wall magnetization profile through application of an external magnetic field and switching off the magnetic field at the same time the energy gap is activated. An energy gap is produced by either applying a staggered magnetic field in the $z$ direction or introducing a modulation to the $XY$ coupling. The magnetization, spin current, and spin-spin correlation functions are computed analytically in the thermodynamic limit at long times after the quench. For both types of systems, we find the persistence of power-law correlations despite the ground-state correlation functions exhibiting exponential decay.


I. INTRODUCTION
Many-body dynamics in closed-quantum systems has exploded in popularity as an area of intense research over the last decade [1][2][3][4][5]. The surge in popularity of low-dimensional quantum dynamics as a theoretical area of study has largely been spurred by remarkable advances in the experimental simulation of low-dimensional systems with tightly controlled and highly tunable interactions. With simulations of toy-modellike systems readily available in the form of carefully designed setups which manipulate ultracold atoms in optical traps [6][7][8][9][10][11][12], it is possible to explore questions of thermalization and relaxation experimentally using models which are simple enough to allow for a thorough analytic treatment.
Perhaps the most well-known class of such theoretically simple systems is that of quantum spin chains, which consist of a one-dimensional lattice of spin degrees of freedom. Bethe first solved [13] the Heisenberg spin chain employing a method which would become known as the "Bethe ansatz," which has since been developed into a framework suitable for investigating the exact dynamics of integrable systems away from equilibrium [14]. The simpler model was introduced and shown to map to a system of free fermions by Lieb, Shultz, and Mattis [15]. Due to the simplicity of the basic model, it has become a popular system for investigating nonequilibrium physics analytically. The spin chain with an external magnetic field leads to a surprisingly rich phase diagram, and its correlation functions were explored quite thoroughly in the work of Barouch, McCoy and collaborators [16][17][18][19].
It is noteworthy that Ref. 16 devotes a section to computing the dynamics of observables following an abrupt change in the external magnetic field-an early example of a "quantum quench" protocol [5]. In recent years, the quantum quench has become an extremely popular device for generating nonequilibrium dynamics which can be simulated experimentally through experiments with ultracold atoms with rapidly tunable external fields [20,21]. Experimentally, tuning the parameters of the initial system allows one to probe the dynamics after the quench for a wide variety of carefully chosen initial states.
From the theorist's point of view, a quantum quench is essentially a formal procedure for investigating time evolution with an arbitrary initial state. After a sudden change in system parameters, the ground state of the initial eigenstate is no longer the ground state of the Hamiltonian which generates time evolution. Very often, the initial state is no longer any eigenstate of the final Hamiltonian. While generic, nonintegrable systems coupled to an external environment would be expected to thermalize, the situation is more subtle for isolated many-body systems. The eigenstate thermalization hypothesis provides a general mechanism by which expectation values of observables in isolated, nonintegrable systems can approach thermal values despite the absence of coupling to an external environment [22][23][24].
For integrable systems, which contain a conserved quantity for every degree of freedom, the dynamics is tightly constrained. Consequently, an integrable system generally cannot relax fully to thermal equilibrium in the traditional sense and instead relaxes to a generalized Gibbs ensemble (GGE) [25][26][27]. The crossover between integrable and non-integrable systems results in a smooth transition [28]. For systems in which integrability is only weakly broken, the system often relaxes to a long-lived GGE before eventually thermalizing at very long times [29,30]. Even for systems which can be mapped to noninteracting quasiparticles, observables can exhibit nontrivial dynamics as a consequence of dephasing of the system's eigenmodes [31]. Moreover, operators which are nonlocal with respect to the quasiparticles can lead to effectively thermal behavior [32] despite the trivially integrable nature of the noninteracting system.
The majority of studies of quench dynamics focus on spatially homogeneous systems. However, spatial inhomogeneities in the initial state provide particularly simple means of generating nonequilibrium dynamics. Experimentally, domain-wall magnetization configurations can be created by application of a spatially varying magnetic field [33,34]. In the isotropic spin chain, a linearly varying magnetic field can be used to create a domain-wall magnetization profile. By suddenly switching off the magnetic field, the domain wall is observed to spread ballistically [35]. However, evidence of the system being far from equilibrium is found in the athermal relaxation of the transverse spin-spin correlation function, which acquires oscillations at the scale of the lattice [36]. The isotropic model may be recast in terms of hardcore bosons, in which these oscillations correspond to quasicondensates of the bosonic modes [25,37]. In the spin language, the oscillations are directly related to a spin current [38,39] which emerges as magnetization is being transported across the domain wall as the profile spreads. It has been shown recently that the time evolution of the domain-wall initial state is well described as an eigenstate of an emergent, effective Hamiltonian [40,41].
Domain-wall dynamics have also been investigated in the anisotropic model [42][43][44][45][46][47]. Like the isotropic model, the anisotropic model maps to free fermions so that the dynamics may be computed exactly. Aside from exact diagonalization methods [48], domain walls in the interacting model have been treated numerically using time-dependent density matrix renormalization group ( DMRG) methods [49,50] and theoretically using a bosonization approach [36]. More recently, generalized hydrodynamics [51,52] has been employed to obtain analytic results which agree well with numerical calculations.
With some notable exceptions [53,54], detailed investigations of quench dynamics in which translation-invariance is broken by the Hamiltonian generating time evolution are scarce in the literature. Breaking translation invariance in spin chain models is most easily accomplished by adding a position dependence to external magnetic fields or nearest-neighbor couplings. In the present work, we address both types of terms.
All models considered are extensions of the isotropic model as described by the Hamiltonian [15] In this paper we break translation invariance by performing each of the following two modifications to Eq. (1). First, we add a term corresponding to a "staggered" magnetic field which oscillates at the scale of the latticê where > 0 is a constant. The second type of modification is to modulate the nearest-neighbor, coupling → = 1 − (−1) for 0 < < 1. Such modulation leads to dimerization in which an energy gap appears in the spectrum [55]. The addition of a staggered magnetic field also leads to the appearance of an energy gap, and this type of perturbation has been used previously [53] as a particularly simple mechanism to generate a gap in the energy spectrum of free fermions confined to a lattice. Additionally, a self-consistent version of the staggered field perturbation has recently been used to investigate the dynamics of a Bose-Hubbard model with a particular form of global-range interactions in the thermodynamic limit, where the system maps to a free-fermion model with a self-consistency condition [56,57].
The work presented in this article focuses on the long-time behavior of one-body observables and correlation functions in the nonequilibrium steady state that forms in the center of the system as the domain wall broadens, or "melts." Previous efforts have demonstrated that the central subsystem relaxes to a GGE-like steady state when time evolution is generated by the isotropic [50] and anisotropic [46] models. In these cases, an effective momentum distribution describing the long-time limit of the central subsystem is obtained by expanding the initial momentum correlation matrix for small [50,58]. As these models map to free fermions, all observables can be computed from this effective momentum distribution. In this paper we extend the analysis to dynamics generated by the two gapped models discussed above. Due to the broken translation invariance, the effective momentum distribution is replaced by an effective Wigner function which inherits an explicit position dependence.
Despite the presence of an energy gap in the spectrum of the Hamiltonian which generates time evolution, we ultimately find the persistence of power-law correlations in the spin-spin correlation functions. The existence of power-law correlation functions in a gapped model far from equilibrium has been demonstrated previously in continuum systems [59,60] and more recently in lattice systems [57]. One noteworthy feature of our results is that the power-law correlations depend crucially on the domain-wall initial state. As we will discuss, a homogeneous quench from the ground state of the model results in at least one correlation function decaying exponentially with distance, similarly to its behavior in the ground state of the gapped model. The paper is organized as follows. Section II contains a summary of the ground state properties of the models considered to establish a baseline to which the nonequilibrium results may be compared. The main results for the long-time limit of observables in both types of gapped systems are derived in Section III. Also included in this section is a brief investigation of quenches in which the initial state does not possess a domain wall magnetization profile and a comparison of our results to domain-wall dynamics in the gapped anisotropic model. Lastly, Section IV contains a brief discussion and outlook on future work.

II. GROUND STATE OBSERVABLES
This section is devoted to a brief exploration of the groundstate properties of the models considered in the remainder of the paper. In the process of obtaining the low-energy physics, the basic nomenclature used throughout the paper is introduced. This section also contains the definitions of the particular observables which will be computed away from equilibrium in Sec. III.

A. Staggered field
The isotropic model with a staggered (alternating) magnetic field is described by the following Hamiltonian, We employ units in which ℏ = 1 and the lattice spacing → 1. Previous work has investigated both static [61,62] and dynamical [63][64][65] properties of the isotropic model in the presence of a staggered magnetic field. The Jordan-Wigner transformation [15,66] may be employed to write the spin operators in terms of spinless fermions, In general, Eq. (25) may be computed numerically. Often, it is possible to apply the Szegő limit theorem and related Fisher-Hartwig conjecture to extract the asymptotic exponential or power-law decay analytically for large [68][69][70]. Such a procedure only applies when ( ) is a Toeplitz matrix, having entries ( ) = − (1 ≤ , ≤ ), which is not the case here due to the term in Eq. (17) proportional to (−1) . Such a procedure does work for the model [71]. However, often one may still use the framework of the Fisher-Hartwig conjecture to understand how basic features of the correlations arise.
the asymptotic determinant of with entries ( ) = ( ) − and size can be extracted from the properties of̃ ( ) ( ) when the Toeplitz condition is satisfied. In that case,̃ ( ) can often be factored into a product of some smooth part 0 ( ), and a finite number of zeros and jump discontinuities occurring at points ,̃ The jump discontinuities are parametrized by numbers through the function . Given the form in Eq. (27), the Fisher-Hartwig conjecture states that the asymptotic form of the determinant of ( ) is given by where  is a constant independent of and 0 = ∫ − 2 ln 0 ( ).
A procedure exists for calculating , but it is quite involved whenever 0 ≠ 0 [71]. The additional position dependencẽ ( ) →̃ ( ) ( ) in the present case, wherẽ prevents one from applying the Fisher-Hartwig conjecture directly. However, formally factoring Eq. (17) into a smooth nonzero function and a product of zeros and jump discontinuities gives where cos ( ) 0 ≡ (−1) ̃ , and we takẽ < 1. Naively attempting to apply the Fisher-Hartwig conjecture, one finds Eq. (31) consists of two zeros, each with = 1 2 and two jump discontinuities, each with = 1 2 . Applying Eq. (28) predicts a vanishing power-law factor. The smooth envelope 0 ( ) results in a nonzero, exponential decay constant (32) so that one might expect where  0 (̃ ) does not depend on . Remarkably, this simple exponential decay is in agreement with a direct evaluation of Eq. (25) obtained by integrating Eq. (26) numerically using Eq. (30). A comparison of the numerical evaluation of  0 ( ) to Eq. (33) is shown in Fig. 1. It is noteworthy that Eq. (33) is obtained as a speculative prediction based on the structure of the generating function in Eq. (30). It should be emphasized that the Fisher-Hartwig conjecture is not applicable to non-Toeplitz matrices, so it is somewhat of a fortunate accident that the correct correlation length = 1∕ sinh −1̃ is obtained by its use here. Similar casual applications of the Fisher-Hartwig conjecture outside its domain of validity will not lead to correct predictions in the next sections. We also note that the correct amplitudes  0 (̃ ) do not result from this type of analysis. The values of  0 (̃ ) used in Fig. 1 are found by fitting Eq. (33) to the numerical evaluation of  0 ( ), while the decay constant is fully fixed by Eq. (33).
The exponential decay observed in  0 ( ) and 0 ( ) is an expected feature of a correlation functions in the ground state of a system with a nonzero energy gap. As̃ → 0, the energy gap vanishes and both correlation functions reduce to power laws. Note that there is no qualitative change in behavior as̃ is increased above unity, wherẽ = 1 corresponds to = . Larger values of̃ lead to difficulties in accurate numerical calculations, as the correlations decay so sharply with . For this reason, we restrict our attention tõ < 1 in the remainder of the paper, as no qualitative differences in behavior have been observed for̃ > 1. Away from equilibrium, we will find the structure of the generating function in Eq. (30) to change in such a way that asymptotic power-law decay can also be roughly observed to emerge in the scenario considered in Sec. III.

B. Dimerized hopping
An alternative way of generating an energy gap leading to qualitatively similar consequences for physical observables is to consider the dimerized, isotropic chain [55,67,72] in which the nearest-neighbor coupling oscillates in strength, → = (1 − (−1) ) so that Here 0 ≤ < 1, and the limit → 1 results in the system decoupling into isolated pairs of spins. Only values of less than unity will be considered in the remainder of this paper. Equation (35) we find a diagonal Hamiltonian with = √ cos 2 + 2 sin 2 . Diagonalization corresponds to a choice of Bogoliubov angle given by for which | | < 2 . As in the previous section, the ground state | | 0 ⟩ contains all negative-energy quasiparticle modes, Repeating the same basic steps as in the previous section, one finds the basic Majorana contractions are given by ⟨ With no external magnetic field, one obtains a vanishing magnetization in the ground state The longitudinal correlation function has been evaluated previously [67] with Interestingly, attempting to extract the correlation length for the transverse correlation function by appealing to the Fisher-Hartwig conjecture outside its domain of validity as before leads to the incorrect conclusion that → 1∕ ln(1 + ). While this does agree with the numerical evaluation of  0 ( ) for small and → ∞, it fails as → 1. On general grounds, one expects → 0 as → 1 since this limit corresponds to pairs of spins which become decoupled from the rest of the system. As the Fisher-Hartwig conjecture simply does not apply to non-Toeplitz matrices, that it does work in the case of a staggered magnetic field it is more surprising than that it does not give the correct correlation length in the present setting.
However, taking exponential decay with exactly half the correlation length of  0 ( ) does give an asymptotic form for the transverse correlation function which agrees quite well with a numerical evaluation of  0 ( ), as → ∞. The expression in Eq. (44) is shown alongside a numerical evaluation of  0 ( ) along the same lines as that discussed in the previous section for several values of in Fig. 2.
The correlation function decays exponentially with for all values of ≠ 0, and Eq. (44) correctly captures the strength of the exponential decay. Again, the amplitudes  ′ 0 are obtained by fitting Eq. (44) to the numerical evaluation of  0 ( ). Additionally, the effects of dimerization are seen in the "staircase" structure of  0 ( ) which becomes more prominent as → 1.

III. DYNAMICS FROM DOMAIN-WALL INITIAL STATE
The main focus of this work concerns how the equilibrium correlation functions presented in the previous section are modified when the systems are far from equilibrium. To drive the system into a nonequilibrium steady state, we assume an initial state | | Ψ 0 ⟩ corresponding to a domain-wall magnetization profile, More general domain walls in which the transition region has nonzero width or the system is only partially polarized far from the central regions have also been considered in the model [35,36,46,50]. The former situation has no effect on the long-time dynamics or formation of the steady state. The latter scenario in which the system halves are only partially polarized is a straightforward extension of the main results in this paper. The basic expressions needed to compute observables are somewhat lengthy and relegated to Appendix B. Restricting attention to the state in Eq. (45), we may write | | Ψ 0 ⟩ in terms of Jordan-Wigner fermions by using the site basis and only occupying sites on the left half of the system, Our main interest for the remainder of this section is in computing the long-time limit of observables such as spin current and correlation functions with the initial state | | Ψ 0 ⟩ which is not an eigenstate of either Hamiltonian (c.f., Eqs. (7), (35)).
At long times, a non-equilibrium steady state (NESS) forms in the central region of the system in which the magnetization relaxes to zero and the spin current saturates to its asymptotic value. Local observables within this region may be computed in terms of the basic fermionic correlation function where ( ) is the effective momentum distribution of a nonequilibrium steady state which carries information about the initial domain-wall configuration. This framework has been applied previously to study relaxation of observables in the isotropic and anisotropic models [46,50]. Given ⟨ † + ⟩ NESS , the long-time limits of all local observables may be computed by application of Wick's theorem. What distinguishes the focus of this work from these previous efforts is the broken translation invariance through the staggered magnetic field or modulated hopping amplitude. These complications result in a position dependence being inherited by the effective momentum distribution, ( ) → ( ) ( ). The proper interpretation of ( ) ( ) is actually the Wigner distribution [50], though the distinction is largely unimportant in cases without explicit position dependence.
The overall strategy in this section is to exploit the quadratic nature of each Hamiltonian to obtain explicit expressions for the time-dependent fermion operators ( ) and compute the basic two-point function ⟨ † ( ) + ( ) ⟩ . As pointed out in Ref. 50, the long-time limit will depend on momentum correlations in the initial state, ⟨Ψ 0 | | † ′ | | Ψ 0 ⟩, with = + 2 , ′ = − 2 and small but nonvanishing. In this limit, The pole and residue structure encoded in Eq. (49) contain all the information needed to extract the long-time limit of the Wigner distribution ( ) ( ) in the non-equilibrium steady state. As the models considered are quadratic, Wick's theorem reduces all observables to functions of this basic distribution function. The form given in Eq. (49) is specialized to the fully-polarized semi-infinite subsystems considered here. Appendix B contains a discussion regarding generalization to domain walls of arbitrary heights.

A. Staggered magnetic field
To compute the Wigner distribution ( ) ( ) in the long-time limit, we must first obtain expressions for the time-evolved fermion operators ( ) in the Heisenberg picture. Given the Hamiltonian Eq. (7) and Bogoliubov rotation in Eq. (8), the time-evolved position-basis operators can be written as where = cos( ) − cos sin( ), = sin sin( ), with = √ cos 2 +̃ 2 . The exact expression for the basic two point function is All expectation values are with respect to the initial state | | Ψ 0 ⟩. Cross terms proportional to ⟨ † ′ + ⟩ and ⟨ † + ′ ⟩ have been dropped, as these contractions give negligible contributions (see Eq. (49)). Changing variables of integration from ( , ′ ) to ( , ) via = + 2 and ′ = − 2 , using Eq. (49) lets us perform the integration over in the long time limit → ∞ as a contour integral, yielding the NESS result where the Wigner function ( ) ( ) is given by Equation (56) underlies all of the remaining results in this section, and its form is strikingly similar to the equilibrium result (c.f. Eq. (15)). In taking the formal limit → ∞, our results apply at large but finite times for | |, | + | ≪ . Indeed, ( ) ( ) relaxes in the long-time limit to its equilibrium form but with the additional factor ( ) which is equal to ±1 and consists of four piecewise-constant segments. This combination of step functions arises in part from projecting the integral to the entire range − < < . This series of sign changes provided by ( ) has dramatic consequences for the observables and contains all of the surviving information about the initial state.

Magnetization and spin current
From Eq. (53) one may compute observables at arbitrary times. Given that̂ maps to free fermions in Eq. (10) and that the initial state in Eq. (47) can be represented by the ground state of a different quadratic Hamiltonian, it is possible to obtain the corresponding results for a finite system from a sequence of direct matrix diagonalizations, as described in Appendix A. Figure 3 shows the magnetization dynamics for various choices of with a sharp domain-wall initial state. The transient dynamics obtained numerically from the method described in Appendix A are virtually indistinguishable from the analytic formula given in Eq. (53). Using Eq. (56), the longtime limit of the magnetization in the NESS is shown to vanish, In the case of = 0, this vanishing magnetization represents relaxation to the ground state average in the isotropic spin chain. In the gapped model, the ground state possesses a staggered magnetization given by Eq. (18). Thus, while the magnetization relaxes to zero, the vanishing magnetization is quite different from the equilibrium profile, showing that the system remains in a highly excited state of the Hamiltonian̂ .
The spin current , which vanishes in the system's ground state, becomes nonzero due to the net flow of magnetization from the left side to the right side of the system. The precise form of follows from the continuity equation for the magnetization, The "light cone" swept by the nonzero current corresponds to the growing central region in which the magnetization relaxes (asymptotically) to zero. As ∕ becomes larger, the domain wall spreads more slowly and the resulting steady-state current becomes smaller.
Applying the Heisenberg equation of motion for time evolution generated by the Hamiltonian̂ in Eq. (7), one can identify which is the same as the expression for spin current in the isotropic model [46]. For a domain-wall initial state, Eq. (56) can be used in Eq. (63) to obtain ⟨ ⟩ It should be noted that the algebraic decay in Eq. (65) is qualitatively similar to that obtained from a quench from an initial state which is spatially homogeneous but supports a nonzero spin current [60].

Correlations
To obtain the general equal-time spin correlation functions, it is convenient to work in terms of the Majorana contractions, which can be written as The general time-dependent spin-spin correlation functions now take the form Our interest here is in the long-time limit after a homogeneous nonequilibrium steady state has formed in the central part of the system where  ( , ) →  NESS ( ) for = , . Using Eqs. (55)- (57), we obtain ⟨ Upon expanding Eq. (71) for large , we find with the above expressions holding in the limit → ∞. For the domain wall initial state, we observe the persistence of powerlaw correlations in contrast to the exponentially decaying correlation function in the ground state of the model. Additionally, there are oscillations at multiple wavelengths present in Eq. (73), as seen from the alternating functional forms for even and odd and the factor of (−1) 2 in Eq. (73). These intricate oscillations and power-law correlations turn out to be a common feature in the models considered. Next, we wish to evaluate the transverse correlations. Due to the vanishing of Eq. (70), the Pfaffian for the transverse correlation function simplifies to where is an antisymmetric matrix with entries given by Antisymmetry defines the entries along and below the diagonal in Eq. (75) while Eq. (76) is valid for all choices of and .
If not for the term proportional to (−1) , would be a Toeplitz matrix for which the Szegő limit theorems and Fisher-Hartwig conjecture [68][69][70][71] could be used to extract asymptotic behavior of the determinant in Eq. (74) as | − | → ∞. As in equilibrium, it is possible to formally define a generating function for the non-Toeplitz determinant, Comparing Eq. (77) to Eq. (31), the only structural difference between ( ) and ( ) is the presence of ( ). This factor of step functions serves to introduce additional jump discontinuities into the generating function. One observes that a possible representation for ( ) is where ( − ) ≡ exp π − − sgn( − ) . The generating functioñ ( ) ( ) should be periodic, and the role of − 1 2 ( − ) is to encode the jump discontinuity from = −0 + to = + 0 + → − + 0 + . Additionally, this factor cancels a lingering factor of 2 arising from 1 2 ( ). We stress that the Fisher-Hartwig conjecture strictly does not apply due to the dependence present in Eq. (77). However, the exponential decay constant in equilibrium was predicted by its naive application and we proceed here using the Fisher-Hartwig conjecture as a guide to explore how the nonequilibrium correlations might be expected to differ from equilibrium. Using Eq. (27), one would expect an additional contribution to the power law exponent of 4 − 1 2 2 = −1 in this nonequilibrium steady state. Recall that the power-law factor was neither predicted by the Fisher-Hartwig conjecture nor actually present in equilibrium (c.f., Eq. (33)). As before, there is also a nonzero exponential decay factor. Comparing this power-law prediction to numerical evaluation of Eq. (74), one finds that the exponential decay observed in equilibrium vanishes entirely, leaving only algebraic decay predicted by the discontinuities in ( ), ⟨̂ ̂ as → ∞. Note that while the Fisher-Hartwig conjecture does predict the −1 power law in Eq. (80), it also erroneously produces an exponential decay factor which is not observed in the actual correlation function. For that reason, Eq. (80) is effectively an ansatz which agrees well with the actual correlation function and contains the power-law decay predicted by the Fisher-Hartwig conjecture. The oscillating prefactors arise by comparison to numerical evaluation of  NESS ( ) discussed below. It is noteworthy that the Fisher-Hartwig conjecture is often unable to predict such oscillations even when the matrix in question is a legitimate Toeplitz matrix [46]. Numerical evaluation of the correlation function for even is shown in Fig. 4 for several values of̃ . Figure 5 shows the absolute value of  NESS ( ), which makes clear the separate branches and power-law decay. The result is that the correlations vanish for all odd , while those for divisible by 4 decay with a different -dependent amplitude than those with not divisible by 4. Interestingly, this set of two different amplitudes also appears in a closely-related model in which the effective energy gap is not suddenly switched to its final value but allowed to increase self-consistently as a staggered field is turned on [56,57]. In that work, the initial state was spatially homogeneous so that no persistent current developed at long times. Consequently, the correlation function  NESS ( ) showed no oscillations of the order cos( ∕2). However, the even and odd correlations split into separate branches in a manner very similar to Eq. (80). Interestingly, the decay was also algebraic but asymptotically the same as the equilibrium = 0 result, with  ( ) ∝ − 1 2 . The main results in this section are the persistence of powerlaw correlations (c.f., Eqs. (73), (80)) in a nonequilibrium steady state despite the system having exponentially decaying correlations in its ground state. In similar spin systems without the presence of an energy gap, the Fisher-Hartwig conjecture has been applied to extract exact asymptotics of the transverse correlations when beginning from the domain-wall state | | Ψ 0 ⟩ [46]. The presence of a staggered magnetic field leads to the effective momentum distribution in the NESS acquiring a position dependence and assuming the form of a Wigner distribution. This position dependence spoils the Toeplitz nature of the matrix whose determinant gives the transverse correlation function, and the Fisher-Hartwig conjecture can no longer be applied directly. However, by observing how the Wigner distribution ( ) ( ) is modified compared to its equilibrium form, one can predict the emergence of an additional power-law decay factor in the NESS correlation function. Interestingly, the exponential decay disappears entirely, leaving an enhanced power-law decay-a fact not easily seen from the explicit form of the Wigner distribution. We shall see below that the dimerized chain leads to extremely similar behavior.

B. Dimerized hopping
The calculations for the dimerized chain are quite similar to those for the chain with the staggered magnetic field, so we only sketch the main results and elaborate on the new features in this model. Using the Bogoliubov rotation in Eq. (36), the = sin sin( ).
Repeating the steps from the previous section, one finds the emergence of a well-defined nonequilibrium steady state, where where ( ) is defined in Eq. (57). Similarly to the previous model considered, the only difference between the nonequilibrium result and the equilibrium form is the presence of ( ), which introduces several jump discontinuities into the Wigner distribution.

Magnetization and spin current
Using Eq. (86), the NESS magnetization is calculated to vanish, which happens to be the equilibrium value for this model. One must resort to Eq. (61) to obtain the correct form for the current operator in the presence of nonzero dimerization, obtaining The modulating amplitude is compensated by oscillations within the expectation values of the fermion operators, and one finds ⟨ ⟩ The isotropic result is recovered in the limit → 0. Furthermore, the current vanishes as → 1, consistent with complete dimerization in which individual pairs become isolated and cease interacting with the rest of the system.

Correlations
As in the previous section, all observables follow from Eq. (86). The  NESS ( ) correlation function is again the most straightforward to calculate, giving formally the same result as obtained for the staggered field, For → ∞, the asymptotic result is Again we find power-law decay with oscillations of several wavelengths. Turning attention to  NESS ( ), the same factor of ( ) which appeared in the generating function when the staggered magnetic field was switched on speculatively suggests an additional power-law decay factor of −1 . Figure 6 shows the result of numerically evaluating the appropriate determinant for various values of , demonstrating that in the nonequilibrium steady state, the exponential decay in the ground state is replaced again by purely algebraic decay of the form −1 . As in the previous case, the odd correlations vanish and the even correlations split into two branches with different amplitudes for 2 even and 2 odd, To summarize the results for the dimerized chain, we find correlations which are virtually identical to those found in the chain with staggered magnetic field. The distinguishing characteristic of both models is an energy gap which opens at = ± 2 , resulting in a doubling of the unit cell and corresponding reduction in the first Brillouin zone. It appears that this reduction in size of the Brillouin zone, which is mathematically the root cause of an enhanced power-law, is intimately related to the stronger correlations appearing in the NESS compared to the ground state.

C. Quench from ground state of model
The results in the previous session warrant some further investigation into the models considered to understand how generic the power-law correlations are after a quench into gapped phases. Here we explore quenches beginning from the spatially homogeneous ground state of the isotropic model without a magnetization domain wall. Again, we find similar behavior when either a staggered magnetic field or dimerization term is switched on, but interestingly the power-law correlations do not emerge. Except where explicitly indicated, similar results are obtained for both the staggered magnetic field and dimerized model. Consequently, results are presented only for the staggered magnetic field. Let us consider the ground state of the model, which is obtained from Eq. (7) by setting = 0. In this case, the model is diagonalized by simple Fourier transform, so that the ground state | | 0 ⟩ is composed of all negativeenergy states, Beginning with the system in the state | | 0 ⟩, time evolution takes place under̂ given by Eq. (7). Starting from Eq. (53) and using we find The long-time limit can be taken directly, and the terms can be arranged into a single integral over from = − to = . Alternatively, one may begin from Eq. (B11) and take the limit + = − = 2 . In either case, we find The magnetization at long times follows from ⟨̂ ⟩ = 1 2 ⟨ ⟩ and can be computed exactly, which exhibits a staggered pattern with an amplitude that is not monotonic iñ , as shown in Fig. 7. The peak amplitude occurs for̃ =̃ * , wherẽ * ≈ 0.6627 satisfies Similar such non-monotonic behavior of observables has been found after gapless-to-gapped quenches in a variety of models [73], with a careful investigation of the evolution of correlation functions suggesting that this robust behavior arises from "freezing" of the light cone in these systems.
In the absence of a particle density imbalance provided by the domain-wall initial state, the current vanishes . However, the correlation functions are more interesting. By expanding Eq. (102) in powers of 1 , one finds the asymptotic behavior of the correlation function for large → ∞ given bỹ where ( ) is the "connected" correlation function in which the product of magnetizations is subtracted,  ( ) = ⟨̂ ⟩ ⟨̂ + ⟩ + ( ). Thus, the longitudinal correlation function decays as a power law. Interestingly, this −6 power law has been obtained previously in the continuum limit of this setup [59]. In that work, the corresponding "mass term" sine-Gordon model at the solvable Luther-Emery point was suddenly activated, resulting in similar power-law decay of the corresponding two-point correlation functions.
Despite the power-law decay in the longitudinal correlations, the transverse correlations still show exponential decay.

Using
= † − and = † + , one finds with defined by Eq. (11). One observes that Eqs. (106)-(107) are formally similar to the equilibrium result in Eqs. (16)-(17) aside from the additional factor of | cos | present in Eq. (107). As depicted in Fig. 8, the decay of the correlation function is exponential with some weak oscillations present in the decay. These oscillations are more easily observed by computing the ratio of the nonequilibrium correlation function to its value in the ground state of the gapped model (c.f., Eq. (7)) for a particular value of̃ , as shown in the inset. One observes that there are also weak, subleading corrections to the clean exponential decay which occurs in the ground state.
The quench from the ground state of the gapless model to the gapped model in Eq. (7) most closely resembles the setup in Ref. [56] in which a quench from the gapless phase to the gapped phase of a similar model was investigated. There it was found that the correlation function decayed algebraically  NESS ( ) ∼ − 1 2 . This stark contrast in behavior of correlations is likely attributable to several differences between the model presented here and the model investigated in Ref. [56]. First, we employ a constant "mass term" proportional to whereas Rieger et al. employ a staggered field in which the strength satisfies a self-consistency condition. Thus a quench from the gapless ( = 0) phase to the gapped phase leads to dynamics of the form → ( ) in which the instantaneous mass gap is determined by self-consistency. Such dynamics are entirely absent from the model considered here. Additionally, a constant magnetic field in the direction is also employed, leading to another adjustable parameter which can be used by Rieger et al. to explore a more complex phase diagram than is needed to describe the model in the present work. While we work exclusively with ℎ = 0, the algebraic decay of  NESS ( ) observed in Ref. [56] was found with ℎ ≠ 0 before and after the quench. While Fig. 8 does show some evidence of weak oscillations in the correlations, Ref. [56] finds that the even and odd correlations decay with different amplitudes in a manner qualitatively similar to what we have found for the domain-wall initial state in which the correlations vanish for odd and split into separate branches for 2 even or odd. In the case of the domain wall, power-law correlations do persist in the transverse correlation function but with an exponent which is double its value in the ground state of the model. If the final Hamiltonian is instead taken to be Eq. (35) corresponding to an energy gap provided by dimerized hopping instead of a staggered magnetic field, similar results are recovered. In addition to the exponential decay in the transverse correlation function, we find power-law decay in the longitudinal correlation function, In this section, we have computed the long-time limits of observables when the initial, domain-wall magnetization profile is replaced by the ground state of the isotropic model. The quench simply consists of suddenly switching on the term corresponding to an energy gap, and we find relaxation of  NESS ( ) to exponential decay. However, the longitudinal correlation function  NESS ( ) still retains power-law decay. There is some evidence in the literature of power-law decay of longitudinal correlations being more robust than in transverse correlations in domain-wall initial states in which the individual "halves" of the system are initially held at different nonzero temperatures [42,43]. The domain-wall magnetization profile investigated in this work is formally a zero-temperature state with the spatial inhomogeneity created by a spatially varying magnetic field, and both transverse and longitudinal correlation functions exhibit power-law decay.
The domain-wall magnetization profile is a highly excited initial state, so it is not surprising that power-law decays survive in the presence of an energy gap. The presence of a gap can be expected to only affect the qualitative nature of correlation decay in the ground state or other low-energy settings. However, the ground state of the model is also a highly excited state with respect to the Hamiltonian generating time evolution in both cases considered. Taking the staggered field model as an example, we may examine the quasiparticle number ≡ ⟨ † ⟩ in the initial state using Eq. (8), finds which is valid for all ∈ (− , ). The distribution described by Eq. (111) is smooth compared to the zero-temperature Fermi sea in the ground state, (0) = Θ 2 − | | . The occupation number is plotted in Fig. 9 for several values of̃ as a function of and also as a function of quasiparticle energy, with = √ ( cos ) 2 + 2 . One finds that despite the presence of a gap, high-energy modes in the upper band are already populated in the initial state except at̃ = 0 where the gap vanishes.
The examination of presented here refers explicitly to the homogeneous ground state of the model used as the initial state for dynamics generated bŷ in Eq. (3). As the opera-tor̂ = † commutes witĥ , these mode occupations are conserved and play an important role in the formation of the nonequilibrium steady state. The smoothing of the occupation function in Fig. 9 resembles qualitatively the smoothing that occurs in equilibrium in the presence of a finite temperature. Repackaging this distribution in terms of an effective temperature is rather speculative, as equating Eq. (111) to a Fermi-Dirac distribution with formal "temperature" ( ) would introduce a separate temperature for each mode . Such a distribution is quite different from that of thermal equilibrium in which a single temperature is sufficient to capture the entire distribution.
The results in this section emphasize that in addition to the properties of the Hamiltonian which generates the time evolution, the details of the initial state are extremely important to the nature of the non-equilibrium steady state which emerges at long times. The main result of this section is that the per-sistence of power-law decay in the transverse correlation function appears to be intimately connected to the inhomogeneous, fully-polarized magnetization profile in the initial state. By simply removing the spatial inhomogeneity, the transverse correlation function relaxes to exponential decay.

D. Comparison to anisotropic model
The most notable result of this work is the persistence of power-law correlations in the long-time limit of the function  NESS ( ) despite the existence of an energy gap in the spectrum of the final Hamiltonian. As discussed above, the presence of an energy gap does not necessarily imply that correlation functions away from equilibrium should decay exponentially with distance. However, one does note that the power-law decay of the transverse correlation function changed to exponential decay when the domain wall magnetization profile was replaced by the ground state of the model. Here the currentcarrying NESS was replaced by a steady state with no spin current. A natural question is how generic these power-law correlations are for current-carrying nonequilibrium steady states in gapped systems.
In this brief section, we note a simple situation in which exponential decay of correlations appears in spite of the existence of a spin current. Consider the anisotropic model, In Ref. 46, time evolution of observables in this model was considered at a critical point = ℎ = 1 where the energy gap vanishes for initial states with a domain-wall magnetization profile. It was shown that at long times,  NESS ( ) ∼  cos 2 − 1 4 2 − as → ∞ for some -independent constant , despite the Hamiltonian generating time evolution being critical (i.e., having no energy gap). The correlation function in the ground state of the anisotropic model at this critical point is  0 ( ) ∼  0 − 1 4 , so that the nonequilibrium correlations decay much more rapidly with increasing distance than in the ground state. In the present work, we have the opposite situation in which power-law correlations are persisting in spite of the presence of an energy gap.
The reason for exponential decay in the model at the critical transverse-Ising point is likely due to the extreme mixing of quasiparticles and free fermion operators through a Bogoliubov angle similar to Eqs. (8) and (36) but which mixes creation and destruction operators = + † rather than only mixing creation or destruction operators. The ground state respects this mixing, containing a well-defined number of quasiparticles. However, the domain-wall state will have a well-defined number of particles but not a fixed number of particles.
To explore how this result is modified as we move away from the critical point, we can use Eq. (117) from Ref. 46, where  NESS ( ) = 1 4 det is used to obtain the correlations for general in which the system is gapped. Here is a legitimate Toeplitz matrix, so the Fisher-Hartwig conjecture should capture the → ∞ asymptotics. The gap in the anisotropic model does not induce a doubling of the unit cell, so this energy gap is of a different qualitative nature than those considered in the present work.
Applying the Fisher-Hartwig conjecture in Eq. (28), one finds that at long times where ( ) is a ( -dependent) constant. Thus, unlike the cases considered in the present work, the correlations decay exponentially but with the same oscillatory prefactor that is typically associated with the initial state | | Ψ 0 ⟩. Equation (115) is a particularly simple example of correlations which decay exponentially in a current-carrying steady state, suggesting that there is no rigorous link between the presence of a current and power-law correlations in gapped models.

IV. DISCUSSION
In this work, we computed the long-time behavior of observables in a gapped spin chain after a sudden quench which turns on the energy gap. Our main result is the existence of numerous examples of power-law correlations which exist in a gapped spin chain far from equilibrium. The ground-state correlations of the gapped spin chain generically decay exponentially with distance. However, as we have demonstrated, such power-law correlations can persist in the infinite-time limit when the system begins in a spatially inhomogeneous state which is not an eigenstate of the final Hamiltonian. The two mechanisms for generating an energy gap which have been considered are the application of a staggered magnetic field of constant amplitude and the dimerization of the coupling strength . In both cases, the spin chain Hamiltonian maps to a system of free fermions so that the dynamics may be investigated exactly.
At long times after the quench, both spin-spin correlations decay algebraically rather than exponentially. Interestingly, the domain-wall configuration appears to be quite instrumental in generating the quasi-long-range order, as the transverse correlation function decays exponentially when the ground state of the model is used as the initial state instead of the domain-wall state. The longitudinal correlation function, however, retains power-law decay at long times for both types of initial states. It is also interesting to note that the transverse correlation function exhibits exponential decay when the quench is performed with the anisotropic (gapped) model as the final Hamiltonian and the initial state possesses a domain wall magnetization profile. The energy gap in the anisotropic model arises from a global difference between couplings in the and directions. However, both the staggered magnetic field and dimerized hopping break the translation invariance of the system with periodic perturbations to the homogeneous system at the scale of the lattice, doubling the size of the effective unit cell. It is hypothesized that this doubling of the unit cell is closely related to the power-law decay in the transverse correlation function.
All systems considered map to free fermions. A natural question is how robust the results of this paper are with respect to interactions or even weak, integrability-breaking perturbations which are often present in experimental settings. Finely-tuned experiments with cold atoms have been able to verify similar power-law decay and oscillations within correlation functions [9] computed in noninteracting models. Accounting for non-integrable interactions by exact diagonalization is difficult in practice, where only modest system sizes can be handled [48]. Compounding this limitation, the energy gap leads to slower motion of the domain wall as shown in Fig. 3 so that finite-size effects can influence the results long before the non-equilibrium steady state is reached.
Many of the exact results regarding the transverse correlation function  ( ) both in and away from equilibrium rely on a careful application of the Fisher-Hartwig conjecture, which allows one to extract the asymptotic behavior of Toeplitz matrices as → ∞ provided certain conditions are met by the generating function. In this work, we have encountered systems which lead to the evaluation of determinants of matrices which are not of the Toeplitz form, so the Fisher-Hartwig conjecture does not apply. Interestingly, we have been able to obtain partial information about the transverse correlations by identifying a generalized generating function. However, the results are not consistently accurate. For example, the predicted exponential decay factor does not appear to vanish away from equilibrium despite the transverse correlation function decaying as a power law. It is hoped that the results in this work can serve as a test bed for possible generalizations of the Fisher-Hartwig conjecture to matrices with structure which are not exactly of the Toeplitz form.
One potentially promising direction for investigating quench dynamics in the presence of interactions is the application of field-theoretic methods in the continuum limit. Techniques such as bosonization, while rigorously proven to capture the low-energy physics of lattice models, are somewhat uncontrolled approximations away from equilibrium where operators irrelevant to low-energy behavior might strongly influence the dynamics [53]. The calculation of nontrivial correlations in exactly solvable models, as presented in this work, provides a benchmark which can be used to calibrate approximate techniques for handling interactions in the continuum limit.