{'fraction_non_alphanumeric': 0.03455613014436544, 'fraction_numerical': 0.008780435251023486, 'mean_word_length': 4.703379416282642, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Long-range Kitaev chain in a thermal bath: Analytic techniques for time-dependent systems and environments

We construct and solve a"minimal model"with which nonequilibrium phenomena in many-body open quantum systems can be studied analytically under time-dependent parameter changes in the system and/or the bath. Coupling a suitable configuration of baths to a Kitaev chain, we self-consistently derive a Lindblad master equation which, at least in the absence of explicit time dependencies, leads to thermalization. Using the method of Third Quantization we derive time-evolution equations for the correlation matrix, which we relate to the occupation of the system's quasiparticle modes. These results permit analytic and efficient numeric descriptions of the nonequilibrium dynamics of open Kitaev chains under a wide range of driving protocols, which in turn facilitate the investigation of the interplay between bath-induced dissipation and the generation of coherent excitations by nonadiabatic driving. We advertise this minimal model of maximum simplicity for the study of finite-temperature generalizations of Kibble-Zurek ramps, Floquet physics, and many other nonequilibrium protocols of quantum many-body systems driven by time-varying parameters and/or temperatures.


I. INTRODUCTION
Experimental advances continue to pave the way for probing the nonequilibrium dynamics of many-body quantum systems and for the realization of quantum-based technologies. While intriguing quantum effects are usually most strikingly on display in isolated systems, laboratory realizations of quantum systems will always be affected by some coupling to the environment. Hence, on the theory side, there is a need to complement the experimental developments with analytic and numeric techniques that provide insight into the nontrivial interplay between unitary and dissipative dynamics in open quantum systems at finite temperature.
The goal of this paper is to demonstrate how a combination of established techniques, augmented by model-specific adjustments, can give rise to faithful and exactly-solvable opensystem versions of quadratic fermionic lattice models at finite temperature. To be specific, and for its wealth of physical phenomena to explore, we choose as a setting a long-range Kitaev chain [1,2]. In isolation, this model of a one-dimensional topological superconductor has served as a setting for the investigation of topological phase transitions and topological order. Our aim is to derive for this model a Lindblad master equation that leads to thermalization at late times, so that temperature can be used as a control-knob of the system. Specifically, to explore a broad range of nonequilibrium phenomena, we want to permit time-dependent system Hamiltonians as well as timedependent bath temperatures. This will allow us to study not only finite-temperature effects in Floquet or Kibble-Zurek protocols and the concomitant effects on topological properties of the Kitaev chain, but additionally opens up possibilities to study nonequilibrium quantum dynamics under temperature quenches, temperature ramps, or temperature oscillations.
We start from a microscopic picture in which each site of the Kitaev chain is weakly coupled to its own independent thermal bath, with all the baths characterized by the same temperature. We employ a procedure by D'Abbruzzo and Rossini [3], tailored to quadratic Hamiltonians, which allows us to derive the corresponding Lindblad master equation in a self-consistent manner. In the absence of explicit time-dependencies in the Hamiltonian and the baths, this construction ensures thermalization of the Lindblad equation at late times. This is in contrast to previous studies of nonequilibrium dynamics of open Kitaev chains [4][5][6], where Lindblad equations with local jump operators were postulated ad hoc, which may be reasonable choices for the modeling of dissipative effects at shorter times, but fail to thermalize at late times.
The Lindblad master equation we derive, albeit nonlocal in the jump operators, is quadratic in the fermionic operators. Lindbladians of this type can be solved analytically by a diagonalization technique termed Third Quantization [7]. Going beyond this general technique, we can exploit specific properties of the Kitaev chain and the thermal environment to further simplify the solution, which further increase the accessible system sizes when evaluating the resulting expressions on a computer. In addition to providing eigenstates and eigenvalues of the Lindbladian, the formalism of Third Quantization also provides a convenient setting for studying the dynamics of the correlations between fermionic operators [8]. Here, exploiting again specific properties of the Lindbladian, most notably its translational invariance, the time evolution of the correlation matrix can be broken down into 4 × 4 blocks. Via Wick's theorem, the resulting matrix differential equations give access to the dynamics of arbitrary correlation functions and mode occupation numbers of the open Kitaev chain.
As a first illustration of the usefulness of our results, we calculate the time evolution of the mode occupation numbers of a Kitaev chain with nearest-neighbor hopping and pairing which is in contact with a thermal environment at positive temperature and for which a parameter in the Hamiltonian is ramped (increased linearly in time). Via numerical integration of the time-evolution equations of the mode occupation numbers we obtain data for system sizes of several thousand lattice sites. Plots of these data show the pronounced creation of excitations when ramping the parameter across the critical parameter value of the quantum phase transition, which can be interpreted as a finite-temperature manifestation of the quantum Kibble-Zurek effect [9]. The simplicity of our model not only allows for the efficient numeric treatments of large systems, but also, in certain limits, for the derivation of a rate equation for mode occupation numbers, similar to the one postulated for a spin model in Ref. [9]. Based on this rate equation, scaling laws governing the excitation density at the end of the ramp can be derived analytically. The full extent of the power and versatility of our results is on display in the companion paper "Universal cooling dynamics towards a quantum critical point," in which, based on the analytic results for the time-evolution of mode occupation numbers, Kibble-Zurek-type scaling laws for protocols with time-dependent temperatures are derived [10].
The paper is organized as follows. In Sec. II we introduce the long-range Kitaev chain, outline the procedure through which it can be diagonalized, and provide an overview of the model's zero-temperature phase structure. Section III deals with the coupling of the chain to a thermal environment, starting from an interacting system-bath Hamiltonian to eventually obtain a Lindblad master equation. This equation is treated using the method of Third Quantization in Sec. IV. Section V is dedicated to deriving the equations of motion for fermionic correlations, and relating these to the occupation of the system's quasiparticle modes. In Sec. VI we present results for a particular parameter ramp across the system's critical point, performed at a fixed finite temperature. Here we highlight the competing processes which emerge both in the numeric results and in suitably adapted versions of the dynamic equations. A summary and conclusions follow in Sec. VII.

II. LONG-RANGE KITAEV CHAIN
We introduce a version of the Kitaev chain containing longrange hopping and pairing terms, and outline the procedure through which this model's Hamiltonian can be diagonalized. Aided by these results, we also give a brief overview of the model's zero temperature phase structure. While the results of this section are fairly standard and can be found in the literature, the diagonalized Hamiltonian and the corresponding linear transformations will be useful when dealing with the open-system version of the Kitaev chain in Sec. IV B.

A. Model definition
The long-range Kitaev chain is defined by the Hamiltonian [2,11] where the operators † and create, respectively destroy, a spinless fermion at site on a chain with sites. We enforce periodic boundary conditions, which generate a translationally invariant closed-ring geometry. The parameters and Δ appearing in S are the hopping and pairing strengths, and is the on-site chemical potential. The parameters and govern the power-law decay of the hopping and pairing terms [11] with the distance The factor of 2 1/ in /2 only survives in the hopping terms, and prevents an overcounting of terms with a hopping range of = /2. In the limit , → ∞ we recover Kitaev's original model with nearest-neighbor hopping and pairing [1]. We restrict the discussion to systems where , > 1 in order to avoid complications related to nonextensivity and large-system limits.

B. Diagonalization
The Kitaev chain is an exactly solvable many-body system and can be diagonalized through a combination of Fourier and Bogoliubov transformations; see Refs. [2,3,11,12] for more comprehensive accounts. We introduce Fourier transforms of the fermionic site operators, which allow the Hamiltonian (1) to be recast as with † = † − . Here the lattice momenta ∈ {2 / : ∈ Z } run over the first Brillouin zone. The 2 × 2 Bogoliubov-de Gennes (BdG) Hamiltonian in Eq. (4) reads where we have introduced with the Bogoliubov angle set by In terms of the fermions defined in Eq. (7) the Hamiltonian takes the diagonal form with dispersion relation In the derivation of the Lindblad master equation in Sec. III it will be convenient to combine the Fourier and Bogoliubov transformations into a single canonical transformation which relates the site fermions directly the fermions of the quasiparticle modes. By organizing the site operators into column vectors = [ 1 , . . . , ] T and doing the same for the operators of the quasiparticle modes, the combined transformation can be written as [3] Here, V and S are × matrices with entries In general, for the T transformation to preserve the fermionic anticommutation relations it is required that It is straightforward to verify that V and S indeed satisfy these conditions.

C. Phase structure
The Kitaev chain exhibits two second-order quantum phase transitions at zero temperature. These transitions are signaled by the closure of the system's excitation gap, which occurs when the mode energy , ( ) in Eq. (10) vanishes for a particular momentum . For this to happen, have to hold simultaneously. To determine the phase boundaries, we note that >1 ( ) vanishes only at = 0, . Disregarding the null-pairing case, this reveals the existence of two critical lines in the ( , / ) or ( , / ) planes: one corresponding to (0) = − / and the other to ( ) = − / , as shown in Fig. 1. In this figure as well as in the discussion that follows, either the hopping or pairing are taken to be nearest-neighbor processes, i.e., either = ∞ or = ∞. We have split the phase diagram in Fig. 1 into two regions, shortrange and weakly long-range, with short-range referring to hopping and pairing exponents > 2 and > 2, respectively. Weakly long-range systems are those for which ∈ (1, 2) or ∈ (1, 2). For nearest-neighbor hopping ( = ∞) and arbitrary range pairing ( > 1) we have ∞ ( ) = cos( ), and so the phase boundaries are always located at / = ∓1. For nearest-neighbor pairing ( = ∞) and arbitrary range hopping ( > 1) the more complicated behavior of ( ) results in the phase boundaries deviating from / = ∓1 at finite .

III. KITAEV CHAIN COUPLED TO A BATH
Starting from system and bath Hamiltonians, we derive a Lindblad master equation describing the open-system dynamics of the Kitaev chain. We start from a microscopic systembath configuration where each site on the chain is coupled to an individual bosonic bath. These baths all have the same temperatures and are also otherwise identical. This is depicted schematically in Fig. 2. Together, the chain and baths form an isolated quantum system which evolves unitarily under the Hamiltonian where S is the system Hamiltonian (9) and I S and I E are identity operators on the Hilbert spaces of the system and the environment, respectively. The environment, consisting of identical thermal baths with mode frequencies , is defined by Here , is the bosonic operator associated with mode of the th bath. For the system-environment interaction we choose a linear coupling between the bath modes and sites on the chain, Schematic diagram of the system-environment interaction, where the circles represent sites on the chain and the squares are the independent bosonic baths. Each lattice site is coupled to its own bath, and all the baths together compose the environment.
The coupling strength between bath mode and the th lattice site is denoted .

A. Markovian master equation
The density operator of the system-baths composite evolves unitarily under the Hamiltonian (15) according to the von Neumann equation. Assuming that the baths are prepared in thermal equilibrium states at a temperature , and that their coupling to the system has a negligible effect on their dynamics, it is possible to trace out the environment's degrees of freedom. Following a standard series of approximations, the chain's dynamics is described by the master equation [13] where is the density operator of the chain and LS is the Lamb-shift Hamiltonian. The dissipation superoperatorD on the right-hand side of the master equation reads [3,13] where and are site labels and ranges over all energy differences between S eigenvalues. Here, is the Fourier transform of an environment correlation function. The remaining elements on the right-hand side of Eq.
(20) are the operators obtained from decomposing the operators = ( + † ) in Eq. (18) into eigenoperators of the system Hamiltonian, where denotes an eigenvalue of S and Π( ) is the projection operator onto the corresponding eigenspace [13].

B. Calculation of the dissipator and jump operators
Next we follow the procedure described in sections IV.B and IV.C of Ref. [3] to simplify the dissipator (20) and to bring the master equation (19) into its final form. Starting with the environment correlation functions (21) and following step by step the procedure of Sec. IV.C of Ref. [3], we find that is the density of states of a single bath. This result forΓ ( ) brings the dissipator (20) into the form (25) The eigenoperators ( ) in Eq. (23) can be found by relating the site operators to those of the quasiparticle modes via the linear transformation in Eq. (11). This leads to [3] where = V + S * , with V and S as defined in Eq. (12). To use Eq. (26) in order to further simplify the dissipator (25) and eventually extract the jump operators, we need to keep track of the degeneracies between the energies of the and − quasiparticle modes. Assuming for the moment that there are no zero-energy modes, we can insert Eq. (26) into (25) to producê with Φ = Φ * = * . Further simplifications arise through constraints on the coefficients Φ , which are inherited from constraints on the matrices S and V in Eq. (12) and enter via = V + S * . In addition to the conditions (13a) and (13b) from the Bogoliubov transformation, the particle-hole and time-reversal symmetries of the Kitaev Hamiltonian (1) impose the condition See Appendix B for details. Overall these conditions ensure that is unitary, It now follows that which reduces the dissipator in Eq. (27) tô To finally extract the jump operators, we combine the master equation (19) and the simplified dissipator (31) to obtain the master equation for the Kitaev chain in Lindblad form as Here we have introduced , previously contained in the coefficients in Eq. (18), as a parameter for controlling the system-bath coupling strength. The jump operators in Eq. (32) are read off as with bath coupling constants The first term of Eq. (32) contains the Lamb-shift correction LS , which amounts to a bath-induced modification of the system's unitary dynamics [3,13]. In Appendix A we show that LS is diagonal in the fermions and hence does not introduce a coupling between these modes, but only a shift proportional to in the mode energies in the unitary part of the master equation. Since the system-bath coupling strength is weak, we will neglect the contribution of LS in what follows. In the derivation of Eqs. (32)-(34b) we assumed that the energies of the quasiparticle modes are all strictly positive. As discussed in Sec. II C, this is not the case at one of the critical points of the Kitaev chain, where the energy of either the = 0 or = mode vanishes. If, for example, 0 = 0 then, through a similar series of steps as before, we would obtain a modified = 0 contribution to the dissipator in Eq. (31) which is of the form (35) This suggests that there is now a single jump operator 0 = √︁ 2Γ(0) ( 0 + † 0 ) associated with this mode. Despite this, it can be shown that the different structures of the = 0 terms in Eqs. (31) and (35) do not impact on the dynamics of the mode occupation † 0 0 , which is our primary quantity of interest in later applications. For our purposes, the Lindblad master equation (32) will therefore suffice to treat all the quasiparticle modes, including the zero-energy mode at a critical point.

IV. DIAGONALIZING THE LINDBLADIAN IN LIOUVILLE SPACE
Prosen developed in Ref. [7] a general technique for solving Lindblad master equations in which the Hamiltonian is quadratic and the jump operators are linear in the fermionic operators. This approach makes use of a canonical quantization procedure in Liouville space and goes under the name of Third Quantization. We establish some groundwork for this method in Sec. IV A and then apply it to the master equation (32) of the Kitaev chain in Sec. IV B. We also prove that the master equation leads to thermalization at late times.

A. Third Quantization
We consider a general Lindblad master equation of the form where the superoperatorL is called the Liouvillian. The Hamiltonian is assumed to be quadratic in some set of fermionic modes with = 1, . . . , , and the jump operators are assumed to be linear. The formalism of Third Quantization then allowsL to be diagonalized, yielding an exact solution of the master equation. We outline this procedure for the case whereL is time-independent, see Ref. [14] for details on the time-dependent case. The first step is to introduce 2 Hermitian Majorana operators , written in terms of the (Dirac) fermions as  Note the use of hats to indicate superoperators acting on K. In turn, this allows for the introduction of 4 Majorana superoperator mapŝ which satisfyˆ=ˆ † and the anticommutation relation {ˆ,ˆ} = . Using these Majorana maps, the master equation (36) governing the time evolution of = | can be expressed in the Liouville-Fock picture as If we restrict our study to observables which are products of an even number of Majorana operators, then the antisymmetric structure matrix A of the Liouvillian has entries [7] Assuming that the structure matrix A defined through Eqs. (44a)-(44d) can be diagonalized, there exist 4 eigenvectors , = 1, . . . , 4 with corresponding (generally complex-valued) eigenvalues 1 , − 1 , 2 , − 2 , . . . , 2 , − 2 where Re( ) ≥ 0. To be clear, 2 −1 has eigenvalue , while the eigenvalue of 2 is − . The ordering is such that Re( 1 ) ≥ · · · ≥ Re( 2 ). Hereafter the eigenvalues will be referred to as rapidities. It can be shown that these rapidities coincide with the eigenvalues of the 2 × 2 matrix [8] where H is defined in Eq. (38) and M = Re(M) is the real part of M in Eq. (45). Using the eigenvectors as rows, we construct the 4 × 4 matrix V and impose the normalization condition It follows that V diagonalizes A such that V −T AV T = D with D = diag{ 1 , − 1 , 2 , − 2 , . . . , 2 , − 2 }. Using the V matrix one can introduce a new set of superoperators in terms of which the LiouvillianL is diagonal, and also calculate the steady state expectation value of any quadratic observable . For the latter application, we first introduce the notion of a nonequilibrium steady state |NESS , which satisfieŝ L |NESS = 0 and hence is a stationary solution of the Lindblad equation (43). The nonequilibrium steady state is unique if, and only if, the rapidity spectrum { } does not contain zero. In this case the steady state expectation value of is [7] where |I represents the identity operator acting on the fermionic Fock space F . The expectation values of more complicated observables containing even products of the fermionic operators can be determined using Wick contractions, see Ref. [15] for details.

B. Application to the open Kitaev chain
The master equation (32) derived for the Kitaev chain in Sec. III shares the form of Eq. (36) and is therefore solvable through the Third Quantization procedure. To this end, the generic Dirac fermions will be identified with those of the chain's quasiparticle modes, in which case the Majorana fermions in Eq. (37) become such that the diagonalized Kitaev chain Hamiltonian (9) becomes . Similarly, the jump operators (33) may be written as with the bath coupling constants Γ ,± defined in Eqs. (34a) and (34b). The structures of S and ,± , in which 2 only couples to 2 −1 , ensure that the various quasiparticle modes evolve independently of each other, allowing the problem to be factorized into copies of a single fermion coupled to a bath. This is reflected in the structure matrix A, which turns out to be block-diagonal, with a 4 × 4 block associated to each quasiparticle mode. In the notation of Ref. [7] the block corresponding to the mode is and Γ ,1 = Γ ,+ + Γ ,− and Γ ,2 = Γ ,+ − Γ ,− . The eigenvalues of A then yield the two rapidities The V matrix, constructed using the eigenvectors of A as rows, exhibits a similar block-diagonal structure. The four eigenvectors of A , paired according to their eigenvalues and normalized as per Eq. (47), yield the associated submatrix of V as

C. Thermalization of the open Kitaev chain
Except for a single mode precisely at one of the system's critical points, the rapidities in Eq. (54) are nonzero, which implies that the steady state of the Lindblad master equation for the Kitaev chain is unique. We now show that this steady state is a thermal state with a temperature matching that of the environment. In fact, the set of two-point correlations fully characterize this state [7]. For the present example, these correlations are readily found by combining the general result in Eq. (48) with the expression for V in Eq. (55). The only nonzero off-diagonal correlations are found to be where, from Eqs. (34a) and (34b), one obtains Upon rewriting Eq. (56) in terms of the operators for the quasiparticle modes via Eq. (49), we find, as expected, that the steady-state occupation of these modes is given by a Fermi-Dirac distribution, † = 1 2 1+ This confirms that the coupling of each lattice site to an individual bath, as set out in Sec. III, indeed thermalizes the system.

V. DYNAMICS OF TWO-POINT CORRELATIONS AND MODE OCCUPATION NUMBERS
The Third Quantization formalism described in Sec. IV A also provides a convenient setting for studying the dynamics of the correlations between fermionic variables [8]. In the following two subsections we outline this approach and show how it can be applied to the Kitaev chain to study the time evolution of mode occupation numbers during various types of time-dependent protocols.

A. Correlation dynamics
We combine all two-point correlations of the Majorana fermions into the correlation matrix C with elements The dynamics of C can be derived from the Liouville-Fock master equation (43) and is governed by = modes. In particular, the Lindblad dynamics captured by Eq. (60) cannot modify the correlations between Fourier modes which are not time-reversed pairs. If only correlations between {+ , − } mode pairs are present initially, as would be the case if (0) is a thermal state, then this will remain the case at all later times.
We first consider the correlations associated with a particular {+ , − } pair of Fourier modes, keeping in mind that we are ultimately interested in the occupations of the two associated quasiparticle modes. The relevant Majorana operators are The Using the Bogoliubov transformation (7), the solutions to these equations also provide access to the mode occupation numbers †

VI. EXAMPLE: PARAMETER RAMP AT NONZERO TEMPERATURE
We now turn to applications that illustrate the usefulness of our results for probing the nonequilibrium dynamics of the open Kitaev chain. In particular, we study the excitation dynamics resulting from ramping protocols in which a parameter in the Kitaev chain Hamiltonian is driven towards, or across, its quantum critical value. This offers the opportunity to explore the emergence of Kibble-Zurek physics and universal scaling laws in the context of an open quantum system. In a first scenario, we fix the parameters in the system Hamiltonian at their critical values and then ramp the bath temperature from a positive initial value down to zero. We explore this protocol for cooling the system to criticality, and the scaling laws emerging from it, in our companion paper [10]. Here we focus instead on ramps of a parameter in the Hamiltonian, performed at a fixed positive temperature, where the parameter crosses its critical value. The zero-temperature version of this protocol has been studied previously for the Kitaev chain [11], and also for a range of other systems exhibiting quantum phase transitions [16]. Central to these studies are the scaling laws obtained for the excitation density at the end of the ramp, which can be understood within the framework of the Kibble-Zurek mechanism as a consequence of the breakdown of adiabatic evolution at the critical point. Although quantum phase transitions are strictly zero temperature phenomena, it is expected that a system at small nonzero temperatures will be sensitive to the proximity of the critical point, and that in certain parameter regimes universal scaling behavior might still emerge. This question was investigated using Keldysh techniques in Refs. [9,17] for a spin model, but with a different form of system-bath coupling. More recently, the authors of Ref. [6] considered a class of open systems governed by Lindblad master equations with local jump operators and argued that Kibble-Zurek scaling can persist in the presence of weak dissipation. Here, too, the Kitaev chain was used for numeric benchmarking. In contrast, our open-system model has nonlocal jump operators, which are essential for enforcing thermal-ization at late times, which in turn allows for a well-defined temperature that can be associated with the environment.
The specific example we consider is a Kitaev chain initialized in a thermal state with temperature , which matches the temperature of the baths. The bath spectral density J , which appears in the bath coupling constants (34a) and (34b), is arbitrary at this point. In our numeric implementation we will specialize to an ohmic density J ( ) = exp(− / ) with dimensionless parameter and a high bath cutoff frequency . The chemical potential in the Kitaev chain Hamiltonian (1) is then ramped across its critical value while keeping the bath temperature fixed (see Fig. 1 of Ref. [10]). We consider the = 0 critical point, as identified in Sec. II C. The ramp is performed linearly, with ( ) increasing from to at a constant velocity , such that Solving for the correlations from Eqs. (65a) and (65b) numerically, and then using Eq. (66), allows us to monitor the system's excitation density throughout the ramp. There are two physical mechanisms which impact on the dynamics of E ( ). The first are the coherent excitations, which result from the breakdown of the adiabaticity of the system's unitary dynamics close to the critical point. This is an inevitable consequence of the critical slowing down of the system dynamics due to the closing of the excitation gap at the critical point. The second mechanism is that of incoherent thermal excitations, and possible dissipation, resulting from the system-bath coupling. As emphasized in Refs. [9,17] these processes cannot generally be disentangled to write E ( ) as a sum of two independent contributions, except in limiting cases. Instead, there is a complicated interplay between the two mechanisms. Despite this, it is possible to cast the compact equations (65a) and (65b) into a form with a somewhat more transparent structure. To this end, we invert the Bogoliubov transformation in Eq. (7) and reformulate these evolution equations as ones describing the dynamics of correlations. Here it is crucial to account for the time-dependence that the ramping protocol introduces into the Bogoliubov angle ( ) via the latter's dependence on ( ). The resulting equations are Let us unpack how the coherent and incoherent excitation mechanisms are contained in these equations. Setting = 0 would eliminate the coupling to the bath and isolate the coherent excitations. The resulting dynamics is then equivalent to that of the standard two-level Landau-Zener problem. In this case the rate at which excitations are generated, i.e., the right hand side of Eq. (69a), is proportional to , which is itself a measure of how rapidly the chain's Hamiltonian varies during the ramp. As suggested by the adiabatic theorem, larger will lead to greater deviations from adiabatic dynamics, and a more rapid generation of excitations. In contrast, retaining only the term proportional to in Eq. (69a) isolates the bathinduced dynamics. In this case the expression for Γ ,2 /Γ ,1 in (57) allows us to write with ( / ) = (exp( ( / )) + 1) −1 the Fermi-Dirac distribution. This rate equation describes the relaxation of the mode occupation to its thermal equilibrium value. The associated relaxation rate of the mode is the prefactor = 2 Γ ,1 . Note that varying the bath temperature while keeping the system Hamiltonian fixed would imply that = 0, in which case Eq. (70) provides a full account of the excitation dynamics. This rate equation will serve as a starting point for analyzing the scaling laws resulting from temperature ramps at the critical point in our companion paper [10]. For the present example, we turn to numeric results in order to illustrate the interplay between the two mechanisms highlighted above.
We present two sets of numeric results for the excitation density (68), obtained through an exact numeric treatment of the coupled differential equations (65a) and (65b) using the Bulirsch-Stoer method with a semi-implicit extrapolation modification due to Bader and Deuflhard [18]. This numeric approach was necessitated by the high degree of stiffness exhibited by these equations. For the nearest-neighbor Kitaev chain, Fig. 3 shows the time evolution of the excitation density during the ramp, and includes for comparison results for the closed chain as well. For all but the slowest ramp velocity we observe a rapid increase in the excitation density as adiabatic evolution breaks down in the vicinity of the critical point at = = −1. However, in contrast to the case of an isolated chain with only coherent excitations (dashed curves in Fig. 3), we observe (i) the generation of thermal excitations close to , particularly for slow velocities, and (ii) bath-induced relaxation after the critical point is crossed. The former suggests an interesting interplay between the two mechanisms of excitation, with the ramp velocity being a determining factor in which of the two is dominant. The second observation is a consequence of the system leaving the quantum critical region and crossing into the region where < Δ [17]. To better illustrate observation (i), we show in Fig. 4 the final excitation density E ( ) as a function of the ramp velocity for a ramp terminating at = . Here it is possible to identify a quantitative crossover velocity cr which separates the < cr region where incoherent bath-induced excitations are dominant from the > cr region where coherent excitations will dominate. The range of ramp velocities for which each mechanism is dominant is crucial for understanding why scaling behavior was absent in previous studies [6,9,17] for certain regimes of parameter space.

VII. SUMMARY AND CONCLUSIONS
The objective of the present paper was to construct and solve a "minimal model" with which nonequilibrium phenomena in many-body open quantum systems can be studied analytically under time-dependent parameter changes in the system and/or the bath. In particular, the following properties were asked for: (i) The master equation modeling the open-system dynamics is exactly solvable, such that large systems can be dealt with. (ii) The open system thermalizes at late times under the time evolution induced by the master equation, (iii) The corresponding closed system undergoes an equilibrium quantum phase transition, ideally accompanied by further interesting physical properties, like topological order, whose fate in nonequilibrium conditions and under the influence of a thermal environment can then be studied. Our model of choice is a Kitaev chain to which identical baths are coupled individually at each lattice site. We demonstrate the working of a selection of methods which are well suited to this model and these types of studies, and which open up avenues for a range of applications to be explored. The quadratic nature of the Kitaev chain, in combination with linear systembath couplings, allows for a self-consistent derivation of the corresponding Lindblad master equation by a method due to D'Abbruzzo and Rossini [3], resulting in a particularly reliable and faithful description of the open-system dynamics. The resulting Lindbladian (32)-(34b) is quadratic in the fermionic operators and can hence be diagonalized in Liouville space by the method of Third Quantization [7]. Suitable model-specific adaptations of this method result in an analytic solution of a relatively simple form. Specifically, we obtain a set of decoupled 4 × 4 matrix differential equations (63), which in turn give access to the dynamics of arbitrary correlation functions and mode occupation numbers of the open Kitaev chain. For nonequilibrium protocols with time-dependent Hamiltonian or bath parameters, these differential equations have timedependent coefficients, but numerical integration can nonetheless give access to systems consisting of thousands of lattice sites.
A broad variety of nonequilibrium situations and protocols of quantum many-body systems at finite temperature can be studied with the methods proposed and discussed in this paper. This includes protocols in which system parameters are quenched or ramped at fixed finite temperature, like in the example discussed in Sec. VI, or ramps of the bath temperature, a situation that is explored in detail in the companion paper [10]. Further generalizations include protocols in which a Hamiltonian parameter and the bath temperature are ramped simultaneously, as well as nonlinear ramps, all of which we are planning to explore in future work. Monitoring excitation densities during such ramps can further the understanding of the interplay between the coherent excitations which result from a breakdown in adiabaticity, and the bath-induced heating and cooling of the system. A similar interplay of coherent and dissipative effects, but in a rather different setting, may be studied by considering oscillating Hamiltonian and/or bath parameters, allowing for the investigation of Floquet physics at finite temperature.
The derivations in this paper, and also the applications just mentioned, lead to expressions of moderate complexity, partly because of translation invariance of system and bath (or, in an open-system description, translation invariance of the corresponding master equation). Several interesting extensions of our results become possible once translation invariance is abandoned. This could be done by introducing disorder into the hopping or pairing strengths of the Kitaev chain, resulting in a setting in which the fate of many-body localization in the presence of a thermal bath can be investigated. Another way of breaking translation invariance is by introducing a non-uniform set of thermal baths, with different lattice sites coupled to baths at different temperatures. Even in the absence of explicitly time-dependent parameters or temperatures, such a setting gives rise to nonequilibrium behavior, characterized by nonequilibrium steady states. Augmenting this setting with a time-dependent parameter in the Hamiltonian, universal scaling features may be discovered when the Hamiltonian is parametrically driven across its critical value. To the best of our knowledge, the scaling behavior of excitations on top of nonequilibrium steady states is an unexplored field at present, which may greatly benefit from the concepts and methods put forward and advertised in this paper.