Kibble–Zurek scaling due to environment temperature quench in the transverse field Ising model

The Kibble–Zurek mechanism describes defect production due to non-adiabatic passage through a critical point. Here we study its variant from ramping the environment temperature to a critical point. We find that the defect density scales as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau ^{-d\nu }$$\end{document}τ-dν or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau ^{-d/z}$$\end{document}τ-d/z for thermal or quantum critical points, respectively, in terms of the usual critical exponents and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1/\tau$$\end{document}1/τ the speed of the drive. Both scalings describe reduced defect density compared to conventional Kibble–Zurek mechanism, which stems from the enhanced relaxation due to bath-system interaction. Ramping to the quantum critical point is investigated by studying the Lindblad equation for the transverse field Ising chain in the presence of thermalizing bath, with couplings to environment obeying detailed balance, confirming the predicted scaling. The von-Neumann or the system-bath entanglement entropy follows the same scaling. Our results are generalized to a large class of dissipative systems with power-law energy dependent bath spectral densities as well.

The basic idea behind Kibble-Zurek theory is that when a system is driven to [12,13] or through [6,7] the quantum critical point (QCP) by ramping some control parameter, it undergoes an adiabatic-diabatic transition [16].In the adiabatic phase, the system has enough time to adjust itself to the new thermodynamic conditions, therefore follows its equilibrium state and the defect production is negligible.On the other hand, upon entering into the diabatic regime, the relaxation time of the system is longer than the timescale associated to the drive.Therefore, the system cannot adjust itself to new equilibrium conditions and defects are inevitably produced.The density of defects depends on the rate of change of the control parameter and certain equilibrium critical exponents.
We generalize the Kibble-Zurek scaling for quantum systems containing a QCP, namely the transverse field Ising chain, and coupled to a thermalizing bath within the Lindblad equation.In this case, the relaxation is dominated by the system-bath coupling and not by the intrinsic relaxation scale of the QCP.By ramping down the environment temperature to reach the QCP, we find that the defect density obeys a universal scaling, distinct from the conventional Kibble-Zurek scenario, even when the initial temperature is relatively high.This is attributed to the enhanced relaxation due to bath-system interaction.The thermodynamic entropy of the systems also follows the same scaling.

Kibble-Zurek scaling through driving the environment temperature
We review first the conventional Kibble-Zurek scaling before generalizing it to thermal and quantum phase transition in open quantum systems.We study quenching to the critical point, which satisfies the same scaling as ramping through the critical point [12,13].The reduced temperature is T = T − T c with T c the critical temperature, and it is driven to the critical point as a function of time [16].Here, we use the conventional approach of statistical physics that the system exchanges energy with a large heat bath at temperature T but their interaction is negligible [43], e.g. the canonical ensemble.As a result, the temperature can appear in the Hamiltonian as a parameter through temperature dependent order parameter, external trapping potential etc., and the system is effectively a closed quantum system from the dynamics point of view.When the critical point is approached, the adiabatic-diabatic transition occurs when the rate, at which we drive the system through T (t), becomes comparable to the inverse of the relaxation time τ rel .This follows τ rel ∼ T −zν with z and ν the dynamical critical exponent and the exponent associated to the correlation length [44][45][46].The adiabatic-diabatic transition occurs when these two inverse timescales become comparable We consider linear cooling as T (t) = T0 (1 − t/τ ) with T0 the reduced initial temperature T 0 − T c and τ −1 the rate of change.From Eq. ( 1), the adiabatic-diabatic transition temperature is T (t tr ) ∼ τ −1/(1+zν) at the transition time t tr .After t tr , the system leaves the adiabatic time evolution and defect production takes place.This temperature governs the scaling properties during the diabatic region.The correlation length scales [45] as ξ ∼ 1/( T (t tr )) ν and in a d-dimensional system, the density of defects follows This equation applies for negligible system-heat bath interaction.Therefore, we now discuss the fate of the Kibble-Zurek scaling in the presence of non-negligible system-environment coupling, namely in a genuine open quantum system.In this case, the relaxation properties of the system are also influenced and even dominated by the interaction with the environment rather than the internal relaxation processes, namely the coupling to the environment plays a more important role than the intrinsic relaxation time of the system.Within a Lindblad description [26,47,48], the environment is characterized by an effective spectral density γ, which sets the characteristic damping rate, and possesses a given temperature through the temperature dependent environmental occupation numbers.In thermal equilibrium, the system itself exchanges energy with the bath and takes its temperature.
In the case of driving the environmental temperature, the adiabatic-diabatic transition is determined by effective spectral density of the environment γ.A more complicated case of energy dependent spectral density is discussed at the end of this section.Upon changing the environment temperature, the system temperature also changes.The rate of change of the system temperature should be compared to γ and not to the inherent relaxation time of the system, i.e.
We note that the r.h.s. of Eq. ( 3) contains in principle also the intrinsic relaxation rate of the system, i.e. γ + T zν .However, close to the critical point, the constant environmental coupling γ overwhelms the vanishing intrinsic relaxation rate T zν of the system.In other words, the system relaxes through the faster relaxation channel from the environment (if present) rather than the increasingly long intrinsic relaxation time.For linear cooling, the adiabaticdiabatic transition happens at time 1/γ before the critical point is reached.The temperature at this time instant is T (t tr ) = T0 /γτ .At the scale T (t tr ), the system crosses over from a mainly adiabatic time evolution, when the density matrix closely follows the equilibrium state, to a diabatic time evolution with significant defect production.
The correlation length scales with this temperature as ξ ∼ 1/( T (t tr )) ν and the defect density with respect to the thermal expectation value is for γτ ≫ 1.This applies to thermal phase transitions, driven by T , in the presence of a finite coupling to environment γ.In the limit of negligible coupling to environment, one has to consider the intrinsic relaxation time of the system instead, as discussed below Eq. ( 3), yielding Eq. (1).For a quantum phase transition, which occurs at T c = 0, however, the temperature itself does not drive the quantum phase transition, and the associated thermal correlation length [44,49] scales as ξ T ∼ T −1/z .Then, Eq. ( 4) is modified for a QCP as using again the temperature T (t tr ) = T 0 /γτ at the adiabatic-diabatic transition.For a given τ , the defect density in Eqs. ( 4) and ( 5) is suppressed compared to the conventional Kibble-Zurek case due to the larger exponent.The lower defect density is the consequence of the enhanced relaxation stemming from the bath-system interaction compared to the diverging relaxation time (and vanishing energy scale) for closed quantum systems.In addition, Eqs. ( 4) and ( 5) predict not only the τ , but also the T 0 and γ dependence of the defect density.
We can further generalize these scalings for an environment [26,27,48] with energy dependent effective spectral density γ(E) ∼ |E| s with s > 0 exponent.The s = 1 case corresponds to the common Ohmic bath [48].We find that while Eq. ( 4) remains unchanged, Eq. ( 5) is modified as This follows from realizing that at temperature T , the dominant contribution to damping [47] from environment comes from the E ∼ T states, therefore the r.h.s of Eq. (3) becomes T s through the energy dependent γ.Therefore, Eq. ( 4) remains unchanged since T s c is non-singular for any T c > 0. On the other hand, for the quantum case with T c = 0, we can realize that 1 T d T dt ∼ T s becomes similar to the conventional Kibble-Zurek relation in Eq. ( 1) with the zν → s and T → T replacements.As a result, Eq. ( 5) for the number of defects after driving the environment temperature to QCP is altered to Eq. ( 6) for a power-law spectral density.
We also briefly address the case of non-linear ramps, i.e., when the temperature reaches zero according to T (t) = T 0 (1−t/τ ) p .Following the same scaling arguments presented above, the exponent of Eq. ( 2) is modified to −pνd/(1+ pzν) in accordance with Refs.[50,51].In Eq. ( 4), the exponent changes to −pνd, while in Eq. ( 5) and ( 6) the exponents are modified to −pd/z and −pd/(z(1 + ps)), respectively.Further generalizations are also possible for a time dependent coupling, i.e. γ(t) as in Ref. [52], which is beyond the scope of the present investigation.

Transverse field Ising chain
The paradigmatic example of a quantum phase transition is represented by the one-dimensional transverse field Ising model [6,44,[53][54][55] [56,57] .We demonstrate how the scaling behaviour in Eq. ( 5) emerges explicitly in a system whose dynamics is governed by the Lindblad equation.The model is described by the Hamiltonian Ising coupled spins in a transverse magnetic field as where j runs over the sites of the one-dimensional chain and J > 0. The number of sites is N and the length of the chain is L = N a with a the lattice constant.The dimensionless coupling g > 0 measures the strength of the transverse field.With a Jordan-Wigner transformation (see Methods), Fourier transformation to momentum space and a Bogoliubov transformation, the Hamiltonian reduces to , where E k = 2J (g − cos(ka)) 2 + sin 2 (ka) is the energy spectrum of the fermionic excitations and d km are fermionic operators.In the Hamiltonian, the sum runs over the wavenumbers k = (2n + 1)π/L with an integer n.This quantization corresponds to an antiperiodic boundary condition for the fermionic c operators which is in accordance with periodic boundary condition for the spins [6].
The density of states as a function of energy is calculated as The domain of G(E) is ∆ < E < 2J|g + 1| with ∆ = 2J|g − 1| being the gap.For g = 1, the spectrum is gapless and the system realizes a QCP with critical exponents z = ν = 1 [44] and d = 1.The QCP separates ferromagnetic (g < 1) and paramagnetic (g > 1) phases.In Ref. [6], the quantum quench between the two phases has been studied, i.e., when the system is driven along the g axis through the critical point at g = 1 and zero temperature T = 0.In the present paper, we consider approaching the QCP from another direction on the phase diagram as illustrated in Fig. 1.
Temperature quench to T = 0 within the Lindblad equation We now couple the transverse field Ising chain to a thermalizing bath via the Lindblad equation [18,20,26,58] and consider a quench in which the transverse field is kept constant while the environment temperature is driven linearly from a finite value to zero, see Fig. 1.This yields where D (L; ρ; L + ) = LρL + − {L + L, ρ}/2.In order to thermalize the system, two jump operators are considered for each k and m, which couple to the eigenstates of the Hamiltonian as L km,↑ = d + km and L km,↓ = d km , creating and annihilating a fermionic excitation with quantum numbers m and k, respectively.Thermalization is ensured by requiring the couplings to environment to obey detailed balance corresponding to a bath temperature T as γ k,↓ /γ k,↑ = e βE k with β = 1/T [27].
Since the goal is to investigate a time-dependent variation of temperature, the temperature dependence of the coupling constants γ k,↓ and γ k,↑ is essential.The condition of detailed balance determines their ratio only, while an explicit temperature dependence would follow by performing microscopic derivation of bath correlation functions [27,59].For a system of the type in Eq. ( 8), the temperature dependence can be written as [29,59] The transition rate γ is already independent from temperature.We emphasize that the Lindblad master equation makes sense only for small system-bath coupling γ [26,27,58].We note that the jump operators L km,σ are one of many possible choices to describe detailed balance, more jump processes between states of different wavenumbers could also be included.However, for the sake of simplicity and physical relevance, we focus on the most obvious dissipative processes.
In the followings, we assume that the temperature decreases from an initial temperature T 0 to zero as and describes linear cooling.Consequently, through β(t) = 1/T (t), the coupling constants γ k,↑ and γ k,↓ depend on time.We note that time-dependent coupling as in Eqs. ( 9) and (10) preserve the Markovian approximation leading to the Lindbladian dynamics, as was demonstrated in Refs.[60][61][62]: time-local Lindbladians are Markovian as long as the coupling constants in Eq. ( 9) are positive throughout the time evolution, and can be connected to and derived from a suitably interacting system-environment model.
Assuming thermal equilibrium at t = 0, the probability that the state km is filled with a fermion is calculated from the Lindblad equation as (see Methods) which depends on the wavenumber through E k only.

Defect density
The total density of defects is obtained as where G(E) is the density of states.The defect density represents the number of kinks in the FM state [6,44] as the expectation value of 1   2 ) or the transverse magnetization in the PM phase, n σ x n .For a perfectly adiabatic quench, the system would reach its ground state and no defects would be present.For finite quench duration, however, a finite number of defects is generated at the end of the quench due to the adiabatic-diabatic transition.
The final density of defects n(τ ) depends strongly on whether the system is critical or gapped, which influences the behaviour of the density of states at low energies.For the critical gapless system (g = 1), the density of states is constant, G(E) ∼ const.down to E = 0, while in the gapped phase (|g − 1| ≫ 0), the density of state diverges as G(E) ∼ E/(E − ∆) at the gap edge, E ∆.
In principle, for the gapped phase, the number of defects after the temperature ramp is expected to be exponentially suppressed on general ground, while at or very close to the QCP, a power-law dependence as in Eq. ( 5) is expected.
We start with the behaviour in the gapped phases.The final density of defects is obtained analytically in Methods.For near-adiabatic quenches in the gapped case (g far from 1) with both 1 ≪ β 0 ∆ and 1 ≪ γτ , we find that In both cases, the number of defects vanishes faster than power-law with γτ for long quenches, in accord with the gapped behaviour of the density of states.We note that the limit γτ ≫ 1 implies large values of τ but γ is kept small to be within the validity of Lindbladian description.For the gapless case, g = 1, the final density of defects follows a power-law dependence as We note that the same scaling remains valid even if we stop the time evolution at T f < T 0 before reaching T = 0.This temperature is reached at t f = τ (1 − T f /T 0 ) and using the time-dependence of the defect density n(t) from Methods, we obtain The first term describes the defect density in thermal state at T f , while the second term comes from the surplus defect density which scales again as τ −1 for long quenches.
We have also studied the full lattice version of the model by performing the energy and the temporal integrals numerically in Eqs.(11) and (12).Eqs. ( 13)-( 14) agree nicely with the numerically exact results in Fig. 2 including the γτ and T 0 dependences.Eq. ( 14) indeed shows the expected power-law decay in the adiabatic limit as Eq. ( 5) with z = d = 1.The exponent of the decay differs from the conventional Kibble-Zurek exponent from Eq. ( 2), which would predict n(τ ) ∼ τ −1/2 for the transverse field Ising model [6].The difference is explained by the bath-system interaction preventing the system from "critical slowing down" and enhancing relaxation throughout the diabatic region.We have also checked numerically (see Methods) that our scaling from Eq. ( 6) remains valid in the presence of an Ohmic bath, when the coupling to environment becomes energy dependent as γ ∼ |E| s with s = 1.In this case, the modified scaling reads as τ −1/2 , which is perfectly captured by our exact numerics in Methods.FIG. 2. Numerical results for the density of defects, n(τ ), at the end of the quench for various values of g. a) Cooling to or very close to the QCP.The solid, dashed and dash-dotted line corresponds to the initial temperature of T0/J = 0.1, 1 and 10, respectively.For each initial temperature, the black dotted line shows the τ −1 behavior of Eq. ( 14).b) Far from the QCP, the numerical results (solid lines) agree with Eq. ( 13) (dashed and dotted lines) in limiting cases, the initial temperature is T0/J = 0.1.For g < 1, similar τ -dependences are found.

Entropy
While closed quantum systems are typically described by a wavefunction, open quantum systems possess a density matrix.This allows us to calculate the thermodynamic entropy after the temperature ramp, which also quantifies the entanglement between the system of interest and its environment.The entropy change represents a useful measure of the adiabaticity of the quench.For the transverse field Ising model, after a perfectly adiabatic quench, the system is expected to reach the pure ground state with vanishing entropy.For a quench of finite duration, however, some entropy is unavoidably generated.The irreversible entropy production from Kibble-Zurek theory was touched upon in Ref. [63].Based on the occupation probabilities p ks (t) in Eq. ( 11), the entropy is The final entropy, S(τ ), depends on the final occupation probabilities p(E k , τ ).FIG. 3. Numerical results for the entropy, S(τ ), at the end of the quench for various values of g. a) Cooling to or very close to the QCP.The solid, dashed and dash-dotted line corresponds to the initial temperature of T0/J = 0.1, 1 and 10, respectively.For each initial temperature, the black dotted line shows the τ −1 behavior of Eq. ( 17).b) Numerical results far from the QCP with the initial temperature T0/J = 0.1, the black dashed lines denote 2N n(τ ).
In the gapped phase with 1 ≪ β 0 ∆ ≪ γτ , the final entropy is obtained as S(τ ) ≈ 2N n(τ ), displayed in Fig. 3. Similarly to the defect density, the most interesting case involves quench to the gapless, critical system with g = 1.
For long quenches, we get The final entropy is plotted in Fig. 3 and is in good agreement with the numerical result for long quenches.For quenches terminating away from the QCP, the entropy gets exponentially suppressed in τ , similarly to the defect density.We expect that the residual entropy should scale as S(T = T 0 /γτ ) in general, where S(T ) is the equilibrium thermal entropy of the system.The temperature T 0 /γτ gets imprinted into the dynamics of the system at the adiabatic-diabatic transition.For the transverse field Ising chain, the thermal entropy scales [64] as S(T ) ∼ T , which explains the scaling in Eq. ( 17).

DISCUSSION
We have studied the effect of cooling the environment temperature to a thermal or quantum critical point, and its influence on Kibble-Zurek scaling.We find that the diverging relaxation time, associated to critical points, gets replaced by the inverse coupling to environment, which results in a suppressed, but universal scaling of the defect density.By investigating the dissipative version of the transverse field Ising chain, we verify this prediction and also find that by ramping down the temperature to a gapped ground state, the defect density follows an exponential scaling with the ramp time.The system-bath entanglement entropy follows the same universal scaling and should be accessible experimentally, together with the defect density [8].
For a power-law energy dependent bath spectral function with exponent s, the obtained Kibble-Zurek scaling in Eq. ( 6) is identical to the conventional one in Eq. ( 2) when zν = 1/s.Typically these numbers are of order one, therefore one can easily get identical scaling for a temperature quench to the conventional Kibble-Zurek scenario.In particular, as we demonstrated, an Ohmic bath with s = 1 for the transverse field Ising chain produces the conventional behaviour for the defect density with exponent 1/2.Experimentally, our results can be tested similarly to Ref. [11].
The universal scaling of the defect density in terms of the quench duration for ramping to a quantum critical point applies to a large variety of open quantum systems, ranging from energy independent through subohmic and ohmic to super ohmic bath spectral densities.Not only are these results relevant in highlighting universal features during near-adiabatic cooling processes but can also be beneficial for quantum thermodynamics [65] for efficient heat pumps or quantum refrigerators [66].Moreover, understanding defect production through temperature variations close to quantum critical points promises to be important for smart design of adiabatic quantum computation protocols [67] in open quantum systems.

Diagonalization of the Hamiltonian in the transverse field Ising model
With the Jordan-Wigner transformation σ z j = c j + c + j m<j e iπ nm and σ x j = 1−2n j with nj = c + j c j as introduced in Ref. [6], the Hamiltonian reads By applying Fourier transform to momentum space, we obtain where A k = g − cos(ka) and B k = sin(ka).The Hamiltonian is diagonalized by the Bogoliubov transformation where After Bogoliubov transformation the Hamiltonian reads where E k = 2J (g − cos(ka)) 2 + sin 2 (ka) is the energy spectrum of the fermionic excitations.The energy dependent density of states G(E) is calculated based on the fact that the density of states should be preserves both in momentum and energy space, expressed as G(E)dE = 2 L 2π dk where the factor 2 stems from the m-degeneracy.Substituting the spectrum and expressing the wavenumber with the energy leads to Derivation of the number of defects after temperature quench in the transverse field Ising model In this section, the derivation of the density of defects is presented.The number of defects is defined as N S (t) = km p k (t) where p k (t) = d + km d km (t) is the occupation probability of the fermionic state corresponding to the quantum numbers k and m.
In order to determine the dynamics of p k (t), let us recall the Lindblad equation where Time evolution of p(t) FIG. 4. The time evolution of occupation probability for different quench duration.The solid line is the numerical solution of Eq. ( 27) while the dashed line corresponds to the system in thermal equilibrium at the instantaneous temperature.and γ k,↓ (t) = γ − γ k,↑ with the time-dependent temperature T (t) = T 0 (1 − t/τ ).Note that both the unitary and the dissipative terms are diagonal in k and m.Hence, for each km, the dynamics is restricted to a two-dimensional Hilbert space spanned by the states that km is empty or occupied.Using the empty and occupied states as basis, the density matrix can be represented by with the real-valued functions p k (t) = d + km d km and complex-valued functions q k (t).Based on the Lindblad equation, the following equations are derived.
In our model, we assume that the initial condition is the thermal equilibrium state corresponding to the initial temperature T 0 , i.e., p(0) = (1 + e β0E ) −1 with β 0 = T −1 0 and q(0) = 0.The inhomogeneous differential equations in Eqs.(26) are solved by and q k (t) = 0 for the time interval 0 < t < τ .After the temperature quench, i.e., when the cooling has already ended and the temperature is constant zero, we obtain p(t > τ ) = p(τ )e −γ(t−τ ) .The time dependence of p(t) is evaluated numerically and is shown in Fig. 4 for several quench durations.In the figure, dashed lines show the probability for the system in thermal equilibrium at the instantaneous temperature.It can be observed that for long quenches, the time evolution follows closely the equilibrium values but for short quenches, they differ significantly.
The final number of defects is obtained by summing up Eq. ( 27) leading to N S (t) = F (β 0 )e −γt + γ t 0 e −γ(t−t ′ ) F (β(t ′ ))dt ′ where is the expectation value of the total number of fermionic excitations in the system where ∆ = 2J|g − 1| is the gap.We are mostly interested in the low-temperature behavior, i.e., when the temperature is much lower than the bandwidth during the whole quench, T 0 ≪ 2J.In this situation, only low-energy states are occupied for which the density of states is approximated as and the function F (β) is computed as 4πgβJ e −β∆ if g is far from 1 ln 2 2πβJ if g = 1 (31) where the upper limit of the integral in Eq. ( 29) has been set to infinity.Numerical investigations show that the approximate functions in Eq. ( 31) are in good agreement with the numerically evaluated Eq. ( 29) at low temperatures.Substituting Eqs.(31) if g is far from 1.In the formula, Φ(x) is the error function defined as Φ(x) = 2 √ π x 0 e −y 2 dy.If g = 1, At the end of the quench, t = τ , the density of defects is calculated as for g being far from 1 and for g = 1.
Defect density in the transverse field Ising model coupled to an ohmic thermal bath In this section, we assume that the transverse field Ising chain is coupled to an environment with the coupling constants γ k,↓ (t) = γ(E k ) e β(t)E k 1 + e β(t)E k (37) with γ(E) = γ 0

E
2J and E k is the positive-valued spectrum of fermionic excitations.In contrast to the previous model, the coupling constants are characterized by an effective spectral density proportional to the energy as typical for Ohmic environment [48].The normalization with 2J has been introduced to preserve the dimension of γ 0 .

FIG. 1 .
FIG. 1. Illustration of the phase diagram and the linear cooling protocol, denoted by vertical arrows for the transverse field Ising chain at fixed transverse field.The QCP at g = 1 separates ferro-and paramagnetic phases.