Nonequilibrium current-carrying steady states in the anisotropic XY spin chain

Nonequilibrium current-carrying steady states in the anisotropic spin chain

Jarrett L. Lancaster Department of Nanoscience, Joint School of Nanoscience and Nanoengineering, The University of North Carolina at Greensboro, 2907 E. Gate City Boulevard, Greensboro, NC 27401
August 5, 2019
Abstract

Out-of-equilibrium behavior is explored in the one-dimensional anisotropic model. Initially preparing the system in the isotropic model with a linearly varying magnetic field to create a domain-wall magnetization profile, dynamics is generated by rapidly changing the exchange interaction anisotropy and external magnetic field. Relaxation to a nonequilibrium steady state is studied analytically at the critical transverse Ising point, where correlation functions may be computed in closed form. For arbitrary values of anisotropy and external field, an effective generalized Gibbs’ ensemble is shown to accurately describe observables in the long-time limit. Additionally, we find spatial oscillations in the exponentially decaying, transverse spin-spin correlation functions with wavelength set by the magnetization jump across the initial domain wall. This wavelength depends only weakly on anisotropy and magnetic field in contrast to the current, which is highly dependent on these parameters.

pacs:
05.30.-d, 71.10.Pm, 02.30.Ik, 37.10.Lk

I Introduction

Recent experiments with cold atoms have resulted in a deluge of theoretical work in the area of quantum dynamics of closed, many-body systems Polkovnikov et al. (2011); Cazalilla et al. (2011); Eisert et al. (2015). The ability to engineer faithful representations of low-dimensional models with tunable interactions has ushered in the age of quantum simulators Feynman (1982). By carefully exploiting Feshbach resonances Chin et al. (2010) and precisely tuning optical lattices one may “dial up” a large variety of model systems, such as hard-core bosons Vidmar et al. (2015a) or quantum spin chains Senko et al. (2015) and observe the isolated unitary dynamics of large numbers of strongly interacting particles.

This recent explosion of experiments probing strongly correlated dynamics leaves theorists with the task of characterizing the dynamics of observables in systems driven far from equilibrium. Definitive experimental evidence of interacting systems that do not approach thermal equilibrium within significant timescales Kinoshita et al. (2006) has provided a valuable experimental component to the theoretical investigations of integrable models with dynamics tightly constrained by the existence of a conserved quantity for each degree of freedom. Indeed, one of the fascinating results of the past decade has been the establishment of the generalized Gibbs’ ensemble (GGE) framework for obtaining long-time averages of observables in integrable systems Rigol et al. (2006, 2007, 2008); Cazalilla (2006). For weak integrability-breaking perturbations, these generalized Gibbs’ ensemble averages often correspond to long-lived “pre-thermalization plateaus” with thermalization occurring on extremely large timescales Kollar et al. (2011); Marcuzzi et al. (2013). Careful measurements of multi-point correlation functions during the relaxation of a one-dimensional Bose gas have established its approach to a state precisely captured by a statistical description for timescales on the order on milliseconds Langen et al. (2015). Additionally, a generalization of the fluctuation-dissipation allows for the formal definition of a frequency-dependent (mode-dependent) effective temperature in integrable systems far from equilibrium Bortolin and Iucci (2015).

Treatment of interactions away from equilibrium is generally challenging, especially for spatially inhomogeneous initial states. The time-dependent density matrix renormalization group (tDMRG) provides a popular numerical approach for studying the evolution of interacting systems and has been applied to the spin chain for initial states with domain-wall magnetization profiles Gobert et al. (2005). Application of low-energy, continuum methods such as bosonization to out-of-equilibrium dynamics in the spin chain has provided predictions Lancaster and Mitra (2010); Foster et al. (2010) which compare quite well to numerical treatments Foster et al. (2011); Sabetta and Misguich (2013), at least for weak interactions. Good agreement between continuum bosonization predictions and tDMRG numerical results has also been found for spatially inhomogeneous initial states in spin Heisenberg chains with next-nearest neighbor interactions Bravo et al. (2013). Understanding how the low-energy, continuum description is modified far from equilibrium is also an active area of research, with significant efforts being made to investigate the roles of operators which are irrelevant in equilibrium but strongly affect the physics away from equilibrium Mitra and Giamarchi (2011, 2012); Tavora and Mitra (2013). The majority of recent theoretical efforts have focused on the limiting cases for which the exchange couplings are equal, (, models), or in which only one coupling is nonzero (the quantum Ising model) Cazalilla et al. (2012). Experimental realizations of inhomogeneous states of weakly interacting, two-level pseudospin systems have been presented previously Weld et al. (2009), making experimental investigations of the non-equilibrium dynamics arising from spatially inhomogeneous initial states feasible.

The spin chain, first investigated by Lieb, Schultz and Mattis Lieb et al. (1961) and later generalized to include an external magnetic field Barouch et al. (1970); Barouch and McCoy (1971a); Pfeuty (1970) is a simply stated system with a rich phase diagram. Mapping to free fermions through a Jordan-Wigner transformation, interactions between quasiparticles correspond to interactions between nearest-neighboring -components of spin. Consequently, the chain and its immediate generalizations provide a fruitful playground for theoretical investigations of one-dimensional quantum dynamics far from equilibrium. In particular, magnetization dynamics arising from simple “quantum quenches,” or rapid changes in system parameters, was first investigated decades ago Barouch and McCoy (1971a) by computing observables after suddenly changing the strength of an external magnetic field. Recent efforts have concentrated on current-carrying steady states in the isotropic limit Antal et al. (1997, 1998), evolution of initial domain wall magnetization profiles Antal et al. (1999); Lancaster and Mitra (2010); Karevski (2002), and behavior of correlation functions in the Ising limit after rapid quenches of the magnetic field Rossini et al. (2010); Caneva et al. (2011). Investigations of so-called geometric quenches have also yielded exact results for the limit and numerical results for short time dynamics in the chain Alba and Heidrich-Meisner (2014). A number of previous works have also focused on quench dynamics arising from slow changes in the magnetic field Mukherjee et al. (2007) and anisotropy parameter Mukherjee et al. (2008) in the chain. Recent results have been obtained regarding sufficient conditions for the appearance of Kibble-Zurek scaling (KZS) which depend on details of the initial quasi-particle distribution function Deng et al. (2011). Field theoretic methods Calabrese and Cardy (2006) have been applied to investigate the general features of correlation functions after a quench, with some results for inhomogeneous initial states also available Sotiriadis and Cardy (2008); Calabrese et al. (2008). Analytic expressions have also been obtained for work statistics after sudden anisotropy quenches Bayocboc and Paraan (2015) for zero and nonzero temperature initial states.

As a special case of the spin chain, the one-dimensional quantum Ising model has received a significant amount of attention due to its simplicity, and there are indications that its nonequilibrium dynamics can be probed by employing circuit quantum electrodynamics (QED) Viehmann et al. (2013). Dynamical quantum phase transitions have been theoretically investigated in both the bare transverse-Ising (TI) chain Heyl et al. (2013) as well as the extended axial next-nearest neighbor Ising (ANNNI) chain Kriel et al. (2014). Additional evidence of the rich theoretical structure was provided by the discovery of underlying group structure in a perturbed transverse Ising chain Zamolodchikov (1989); Mussardo (2009) through a formal mapping to a two-dimensional classical Ising field theory. Recently, experiments have led to striking, if low-resolution, observations Coldea et al. (2010) through neutron scattering which are consistent with this emergent group structure. Ising-like models with long-range interactions have also been recently created experimentally in reduced dimensions Richerme et al. (2014); Jurcevic et al. (2015) using cold atomic gases.

Much of the theoretical work on spin-chain dynamics has involved a sudden quench in the system’s transverse field with homogeneous initial states Calabrese et al. (2011); Foini et al. (2011). In this paper, we investigate a sudden quench in the strength of anisotropy and external magnetic field while also taking the initial state to describe an inhomogeneous, domain-wall magnetization profile. While some basic features of domain walls in the Ising model have been previously considered Karevski (2002); Subrahmanyam (2003), the focus of this paper is to thoroughly investigate the time-dependent magnetization, emergent spin current and long-time limit of the two-point, equal-time correlation functions. We will present analytic results for these observables after quenching to the critical transverse-Ising limit and demonstrate the existence of an effective steady state for the central region of the system after the domain wall has effectively flattened out over a finite region of space.

The model is described by the following Hamiltonian,

(1)

Extending some recent results for the isotropic model Sabetta and Misguich (2013), we are able to extract long-time behavior of observables away from the two critical lines and describe the non-equilibrium steady state for with arbitrary . The particularly simple isotropic limit is known as the isotropic model, or model. However, for any nonzero value of and constant magnetic field , the model can be diagonalized, as the system maps to free fermions. In what follows we set the lattice spacing . Interactions may be included by adding the following term to Eq. (1)

(2)

For , the resulting interacting system described by the sum of Eqs. (1) and (2) is known as the chain, which is gapless for . In this paper we specialize to the non-interacting system .

Creating domain-wall magnetization profiles as initial states is an efficient method of generating non-equilibrium dynamics. For the isotropic model, the dynamics is quite easily understood as ballistic expansion of a gas of free fermions Rigol and Muramatsu (2004); Rigol et al. (2006); Lancaster and Mitra (2010). The more complex chain involves interactions between the fermions, but for weak interactions the basic picture appears largely unchanged. Bosonization, a low-energy approach for describing one-dimensional interacting systems, works surprisingly well in this far-from-equilibrium setting Lancaster and Mitra (2010); Sabetta and Misguich (2013). At long times, after the magnetization in a central region of the system has essentially relaxed to its equilibrium value, one finds that spatial oscillations are superimposed on the algebraic decay present in the transverse-spin correlation function,

(3)

with Lancaster and Mitra (2010) and the magnitude of magnetization far from the domain. Here the subscript “g.s.” refers to the expression one would obtain by considering the equilibrium behavior of the system at in which the system is in its ground state. In the language of hard-core bosons, this twist is related to quasi-condensation at during rapid expansion, which has recently been observed experimentally Vidmar et al. (2015a). Representing the system in terms of free fermions through the Jordan-Wigner transformation, both effects may be understood as emerging from an underlying particle current. In the isotropic model, this current emerges according to a conservation law requiring a local change in magnetization to be accompanied by a spin current,

(4)
(5)

Antal and collaborators studied observables in the chain with constant (nonzero) currents through the introduction of a chemical potential-like term involving the total current operator to the Hamiltonian Antal et al. (1997, 1998). Interestingly, the long-time limit of domain-wall states in the gapless phase of the model leads to a non-equilibrium steady state characterized by “twisted” single-particle Green’s functions Lancaster and Mitra (2010); Sabetta and Misguich (2013),

(6)

where , and is the “height” of the domain wall. In the noninteracting model, , so that the wavelength of spatial oscillations is directly related to the spin current. The presence of weak interactions has been shown to not strongly modify this picture in an essential way Lancaster and Mitra (2010); Sabetta and Misguich (2013), though strong interactions may greatly influence the details of the oscillations and decay of correlations Lancaster et al. (2010).

The main goal of this paper is to explore how anisotropy between the and components of nearest-neighbor exchange interactions affects the relationship between initial domain wall height, long-time current and spatial oscillations in . For , the system no longer conserves total magnetization,

(7)

so that the spin current in Eq. (5) no longer corresponds to a local conservation law. However, this spin current does possess nonzero overlap with the actual fermion current obtained after mapping the system to noninteracting fermions through a Bogoliubov rotation. For this reason, a piece of current survives in the long-time limit despite the actual operator not being totally conserved.

This paper is organized as follows. In Sec. II we review some relevant background information on the ground-state properties of and spin chains and provide a brief summary of domain-wall dynamics in spin chains. Section III contains analytic calculations of observables when a domain-wall spin configuration in the -model is subjected to a rapid change in anisotropy and external magnetic field so that the time evolution takes place with respect to the transverse Ising model on the critical line with . Emphasis is placed on the long-time average behavior of the magnetization, current and spin-spin correlation functions. In Sec. IV, an effective nonequilibrium steady state is shown to accurately describe the long-time limit of observables within the central region of the chain. Furthermore, it is shown that observables in this non-equilibrium steady state are easily calculated away from the critical lines and for more general domain walls with arbitrary jumps in magnetization. Results are summarized and discussed in Sec. V.

Ii Model

To formulate the spreading domain-wall problem as a single “quantum quench,” we consider the following protocol: The system is initialized as the ground state of an infinite spin chain with a spatially varying magnetic field,

(8)

where the field acts as a spatially varying chemical potential for the quasiparticles of the system (see below). At , the field is abruptly switched to a constant value at all lattice sites , and nonzero anisotropy between the nearest-neighbor spin interactions is switched on, resulting in time evolution being generated by the anisotropic model,

(9)

The system is generally gapped for nonzero and except along the critical line with and in the critical transverse-Ising limit (). Ground state observables and correlation functions in the presence of arbitrary anisotropy and magnetic fields were extensively studied by Barouch and collaborators Barouch et al. (1970); Barouch and McCoy (1971a, b); McCoy et al. (1971). In this paper we will investigate the non-equilibrium dynamics within the larger phase diagram, specifically in the range . The gapless cases provide instances in which the observable calculations may be performed entirely analytically.

ii.1 Ground state properties of homogeneous chains

It is useful to review the basic ground-state observables in the spin chain as a point of comparison for the non-equilibrium results to be presented. Equation (9) is diagonalized in terms of non-interacting fermions by a Jordan-Wigner transformation Lieb et al. (1961)

(10)
(11)

After applying Eqs. (10) and (11), Eq. (9) becomes

(12)

This quadratic Hamiltonian may be brought into diagonal form

(13)

with

(14)

where the quasiparticles are related to the original momentum-space operators according to a Bogoliubov transformation,

(15)

with

(16)
(17)
(18)

We are interested in the thermodynamic limit in which the sum in Eq. (18) may be converted to a continuous integral. The ground state is then occupied by all negative-energy quasiparticles

(19)

where the vacuum state obeys for all . This convention differs from original treatment Lieb et al. (1961) in which particle-hole symmetry is employed to simplify the algebra by taking all energy eigenvalues as positive and reinterpreting the ground state as a vacuum for all excitations. In the context of inhomogeneous quench dynamics, keeping track of the signs within the present formulation turns out to be a small price to pay for a physically transparent picture. Spin-spin correlation functions at equal times are obtained in terms of Pfaffians, as discussed in Appendix B. All observables can be written in terms of the basic contractions of Majorana operators

(20)
(21)

For example,

(22)
(23)
(24)
(25)

In the isotropic limit with , one finds Lieb et al. (1961); Barouch and McCoy (1971a)

(26)
(27)

In this isotropic case by rotational invariance. In the transverse-Ising limit, , one finds

(28)
(29)
(30)

where is the Glaisher-Kinkelin constant. While power-law decay is observed for and along the line in the plane, correlations generically decay exponentially and the spectrum is gapped everywhere except in the and critical transverse-Ising limits.

Our main focus in the remainder of this paper is on investigating the nonequilibrium behavior of these correlation functions with , as well as the simpler observables such as local magnetization and spin current. Expressions for the correlation functions away from equilibrium are given in Appendix B. Compared to and , the form for is rather simple,

(31)

We explicitly obtain the time evolution of these observables for in Sec. III to determine for which correlation functions the power-law correlations depicted in Eqs. (28)-(30) survive in the nonequilibrium setting. The behavior in other regions of the plane is discussed using a generalized Gibbs’ ensemble description of the central region at long times in Sec. IV.

ii.2 Creating domain-wall magnetization profiles

After applying the Jordan-Wigner transformation, the initial Hamiltonian Eq. (8) becomes

(32)

Let us specialize to the case of a linearly varying magnetic field, . This Hamiltonian maps to the single-particle Wannier-Stark problem Eisler et al. (2009) and may be diagonalized by employing a linear transformation,

(33)

with

(34)

where , and is the Bessel function of the first kind of order . The ground state is populated by all single-particle states of negative energy,

(35)

with . Far to the left (right) of the origin a very large positive (negative) magnetic field acts as a strong chemical potential creating a region of uniform positive (negative) magnetization. In terms of the Jordan-Wigner fermions, the left half of the system is filled with particles with unit -fermion occupation number, for while the right half contains no particles, . Explicitly, the magnetization profile is

(36)

which corresponds to a central spin gradient of width , separating the regions of uniform maximal polarization. The dynamics arising from this state after abruptly switching have previously been considered with the time-dependent magnetization Antal et al. (1999) and transverse correlations Lancaster and Mitra (2010) computed explicitly. One finds the long-time behavior independent of , and for simplicity we consider the limit , corresponding to the sharp domain wall constructed from two semi-infinite, uniform regions of oppositely-polarized spins. Effectively, this allows us to write the initial state as

(37)
Figure 1: Domain-wall broadening in the isotropic model. As a domain-wall magnetization profile spreads, the central region relaxes to a non-equilibrium steady state. The presence of the polarized edges drives a current through the central region. This situation may be viewed as expansion of a quantum gas of interacting ( chain) or non-interacting ( chain) fermions in one dimension.

When time evolution is generated by the model, a long-time steady state emerges Antal et al. (1999); Lancaster and Mitra (2010); Sabetta and Misguich (2013). This steady state is characterized by a spin current, which is driven through the central region as the domain profile spreads, transferring net magnetization from the left side of the wall to the right side. The value of the current is set by the initial jump in magnetization between the far right and far left sides of the system and contains the extent of the central subsystem’s “memory” of the initial state. Denoting the magnetization on the far left (right) side of the system by (), we use the term “domain-wall height” to refer to . That the long-time limit of the nonequilibrium behavior can be entirely described using only the knowledge of this magnetization jump has been demonstrated by Sabetta and Misguich Sabetta and Misguich (2013) in their explicit construction of the nonequilibrium steady state describing the central subsystem at long times. This region where the magnetization has essentially relaxed, as depicted in Fig. 1, is accurately described by a type of generalized Gibbs’ ensemble which depends only on the initial magnetization jump . While the initial states corresponding to domain walls with are somewhat more complicated than the state in Eq. (37), we will address these more general domain walls by extending the results of Ref. Sabetta and Misguich, 2013 to time evolution in the anisotropic -model in Sec. IV.

Iii Domain wall dynamics in anisotropic -model at critical point

iii.1 -model dynamics

The time evolution of the noninteracting quasiparticles is simple, with . When combined with Eq. (15), one may obtain the time dependence of the original -fermions in terms of the initial operators,

(38)

where

(39)
(40)

To obtain time-dependent operators in position space, we perform a final Fourier transform,

(41)
(42)

where

(43)
(44)

iii.2 Analytic solution for the TI critical point

Equations (43) and (44) may be integrated numerically for arbitrary values of and . However, in addition to the isotropic line on the phase diagram where these expressions reduce to Bessel functions, the model possesses another special limit in which closed-form expressions result. Namely, for , the system is gapless and Eqs. (16) and (17) simplify to

(45)
(46)

Performing the integrations in Eqs. (43) and (44) we obtain

(47)
(48)

where is a Bessel function of the first kind of order . Using Eqs. (47)-(48) and the initial state given in Eq. (37), we may evaluate contractions of Eqs. (20)-(21), obtaining

(49)
(50)
(51)

In the above expressions we have suppressed the arguments common to every Bessel function. Our main interest is in the behavior of local observables at long enough times that a central region of essentially flat magnetization has formed, corresponding to the limit . As shown in Appendix A, Eqs. (49)-(51) reduce in this limit to

(52)
(53)
(54)

Strictly speaking, Eqs. (49) and (51) are valid for . For corresponding expressions may be derived which also reduce to Eqs. (52)-(54) in the limit . For the non-interacting problem under consideration, Wick’s theorem may then be used to reduce any local observable to expressions involving Eqs. (52)-(54).

iii.3 Observables

For the critical transverse-Ising model, we may obtain explicit expressions for the current, magnetization and correlation functions in the long-time limit. Furthermore, due to the simple structure of Eqs. (52)-(54), we may reduce the Pfaffian expression for to a determinant of a particular skew-symmetric matrix, which is no more difficult to evaluate than the equilibrium correlation function. To this end, the magnetization in the center region is seen to relax to zero,

(55)

As we shall demonstrate in the next section, the steady-state magnetization does not vanish in the central region for nonzero and arbitrary choices of domain-wall height . Interestingly, the magnetization far from the central region relaxes to half its initial value,

(56)

a result which has been previously obtained Karevski (2002). Since total -fermion number (the magnetization) is no longer conserved for , this suggests that the -fermion number operators share some nonzero overlap with the -fermion number operator which is conserved under the final Hamiltonian. After the initial magnetization relaxation, the domain wall dynamics appears similar in spirit to the isotropic model, and one might reasonably expect a nonzero current in the center. Indeed, one finds

(57)
(58)

Equation (58) is interesting for two reasons. First, despite the domain wall height decreasing to half its initial value after the initial relaxation, the current is only reduced to of its value in the isotropic spin chain, . Additionally, if one instead considers a current-carrying, “boosted-Fermi-sea,” initial state Antal et al. (1998); Sabetta and Misguich (2013) with maximum current allowed by the lattice, it follows that the initial current remains conserved in the long-time limit for any values of . That is for

(59)

with , Eq. (24) gives an initial current . Performing time evolution with this initial state and using Eqs. (38)-(40) one finds . One can imagine carrying out the quench in a two-step procedure, first relaxing the spatially varying magnetic field and leaving , . The domain wall spreads ballistically, giving rise to a central steady-state current . After a long time has passed, a sudden quench in anisotropy and magnetic field , will result in some transient dynamics, but the current at long times will average to . Alternatively, if the spatially varying magnetic field is switched off at the same time the anisotropy and homogeneous field are turned on, the current will approach at long times. In this way, the precise details of the quench protocol affect the long-time steady state even in the absence of interactions.

The long-time limit of the correlation functions in the central region may also be investigated. Using Eq. (31) and Eqs. (52)-(54) we obtain a power-law decay for ,

(60)

Turning attention to the transverse correlations, Eq. (149) for may be organized as Barouch and McCoy (1971a)

(61)

where

(62)
(63)
(64)

Here and . Furthermore, according to Eq. (52) in the central region in the long-time limit. Equations (53) and (54) indicate that and depend on and only through the combination . Using the skew-symmetric nature of Q, Eq. (61) may be reduced to

(65)
(66)

where . The matrix in Eq. (66) is an example of a Toeplitz matrix, satisfying . Due to the antisymmetry of Q through , only even values of result in a non-vanishing . The phase may be factored out and is purely real, with for even . Our interest is in extracting the asymptotic behavior of this determinant for large , which is facilitated by defining a generating function, , satisfying

(67)

According to the Fisher-Hartwig conjecture Deift et al. (2011); Its and Korepin (2009); Basor and Morrison (1994), if has the form

(68)

with jump discontinuities occurring at isolated points characterized by exponents and zeros characterized by exponents , then the asymptotic form of the determinant determinant for large is given by

(69)

where is a numerical constant independent of and

(70)

A Fourier transform of Eq. (53) yields

(72)

for . Taking and accounting for the jump discontinuity at , we may represent

(73)

with . This representation makes explicit the two Fisher-Hartwig singularities arising at (, ) and (, ). The numerical prefactor leads to exponential decay and,

(74)

Employing the straightforward methods in Ref. Ovchinnikov, 2007, we may calculate the constant

(75)

where is the Glaisher-Kinkelin constant. Collecting all the pieces using Eq. (69) yields

(76)

The results from under the condition that is even, and . Previously a similar result was obtained Sengupta et al. (2004) by quenching the transverse-Ising model from the ground state of to the critical point . In that case, it was found that exactly. Here we find a power-law decay factor in addition to the exponential decay. Curiously, the same exponent appears in the ground state correlation function, which is given by purely power-law decay [Eq. (29)]. A graphical depiction of at long times is shown in Fig. 2. The Pfaffian expression for may be evaluated numerically or reorganized into the basic structure of Eq. (61). From either approach, it follows that

(77)

despite time evolution being generated by the anisotropic model at the Ising critical point. This equality in transverse correlations is reflective of the rotational symmetry present in the initial state.

Figure 2: (Color online) Long-time limit of transverse correlations for after a quench to the transverse-Ising limit with an initial magnetization domain wall of height . Oscillations are superimposed on an overall exponential decay (inset). The system retains memory of its initial state through both the short wave-length of oscillations and the equivalence between and (not shown).

iii.4 Homogeneous quench

As a specific point of comparison, we may also calculate in the long-time limit after a rapid anisotropy quench from the ground state of the model to the same critical transverse-Ising Hamiltonian. The motivation behind this brief digression is to understand the effects of the anisotropy quench separately from the effects due to time evolution of a spatially inhomogeneous initial state. That is, in this section we briefly consider the simpler quench in which the system begins in the ground state of

(78)

This ground state is given by

(79)

and time evolution takes place with respect to the critical transverse-Ising model, Eq. (9) with , as in the previous section. This situation is described by the same time-dependent operators, which may be written in terms of the Majorana operators,

(80)
(81)

with the initial occupation numbers

(82)

We have suppressed the time-dependence in , for brevity. Taking the long-time limit, one finds

(83)
(84)

Despite being far from equilibrium, the structure of Eq. (61) reduces to a Toeplitz determinant as in equilibrium,