Thermalization of a Lipkin-Meshkov-Glick model coupled to a bosonic bath

TL;DR

This paper derives a Lindblad master equation for a Lipkin-Meshkov-Glick model coupled to a bosonic bath, demonstrating how the bath induces thermalization in the system and analyzing the thermalization rate for large systems.

Contribution

The work explicitly proves that a bosonic bath induces thermalization in the LMG model, which is otherwise non-thermalizing due to integrability.

Findings

Operators reach thermal equilibrium values at long times

Thermalization rate is calculated for large system sizes

Bath coupling restores thermalization in the LMG model

Abstract

We derive a Lindblad master equation that approximates the dynamics of a Lipkin-Meshkov-Glick (LMG) model weakly coupled to a bosonic bath. By studying the time evolution of operators under the adjoint master equation we prove that, for large system sizes, these operators attain their thermal equilibrium expectation values in the long-time limit, and we calculate the rate at which these values are approached. Integrability of the LMG model prevents thermalization in the absence of a bath, and our work provides an explicit proof that the bath indeed restores thermalization. Imposing thermalization on this otherwise non-thermalizing model outlines an avenue towards probing the unconventional thermodynamic properties predicted to occur in ultracold-atom-based realizations of the LMG model.

Equations

$${\mathcal{H}}_S=-\frac{\Lambda }{2S}{S}_x^2-S_z,$$

$${\mathcal{H}}_B=\sum _q{\omega }_q{a}_q^{†}{a}_q,$$

$$\mathcal{H}={\mathcal{H}}_S+{\mathcal{H}}_B+\mathcal{V},$$

$$\mathcal{V}=\sqrt{\frac{\gamma }{S}}{S}_x\otimes \mathcal{B}$$

$$\mathcal{B}=\sum _q{\gamma }_q\left(a_q^{†}+a_q\right)$$

$${\partial }_t{\rho }_{\text{S}}=\mathcal{L}{\rho }_{\text{S}}$$

$$\mathcal{L}\bullet =$$

$$-\frac{|{\kappa }_{xy}{|}^2}{{\kappa }_{xx}}[S_y,[S_y,\bullet ]]+2\mathrm{Re}({\kappa }_{xy})[S_y,[S_x,\bullet ]].$$

$${\partial }_{t}A(t)={\mathcal{L}}^{†}A(t)$$

$${\mathcal{L}}^{†}\bullet =$$

$$-\frac{|{\kappa }_{xy}{|}^2}{{\kappa }_{xx}}[S_y,[S_y,\bullet ]]+2\mathrm{Re}({\kappa }_{xy})[S_x,[S_y,\bullet ]],$$

$${\partial }_{t}A(t)={\mathcal{L}}^{†}A(t)=e^{t{\mathcal{L}}^{†}}{\mathcal{L}}^{†}A,$$

$${\mathcal{L}}^{†}(AB)\ne {\mathcal{L}}^{†}(A)B+A{\mathcal{L}}^{†}(B),$$

$$e^{t{\mathcal{L}}^{†}}(AB)\ne (e^{t{\mathcal{L}}^{†}}A)(e^{t{\mathcal{L}}^{†}}B).$$

$$[S_{\mu },S_{\nu }]=i{ϵ}_{\mu \nu \delta }{S}_{\delta }$$

$${\mathcal{L}}^{†}\bullet \approx i[\mathcal{H},\bullet ]+i\mathrm{Im}({\kappa }_{yx})\{S_y,[S_x,\bullet ]\},$$

$$\begin{gathered} {\mathcal{L}}^{†}(AB)\approx {\mathcal{L}}^{†}(A)B+A{\mathcal{L}}^{†}(B) \\ +i\mathrm{Im}({\kappa }_{yx})\left([S_y,A][S_x,B]-[S_x,A][S_y,B]\right). \end{gathered}$$

$${\partial }_{t}A(t)=i[{\mathcal{H}}_S(t),A(t)]+\frac{i\gamma }{2S}\left\{S_y(t),[S_x(t),A(t)]\right\},$$

$$\mathbf{s}=(x,y,z)\equiv (S_x,S_y,S_z)/S.$$

$$\dot{x}(t)$$

$$\dot{y}(t)$$

$$\dot{z}(t)$$

$$x^2+y^2+z^2=1,$$

$$h(t)={\mathcal{H}}_S(t)/S=-\frac{\Lambda }{2}x(t{)}^2-z(t),$$

$$\dot{h}$$

$$\ddot{x}$$

$$x(t,\gamma )$$

$$h(t,\gamma )$$

$${\partial }_{\tau }{h}_0$$

$${\partial }_{\tau }^2{x}_0$$

$${\partial }_s{h}_0$$

$$V(x)=\frac{{\Lambda }^2}{8}{x}^4+\frac{1}{2}(1+\Lambda h_0)x^2.$$

$$x_0(\tau ,s)=a\mathrm{cn}\left(\Omega \tau +u,k^2\right)$$

$$\Omega =\sqrt{a^2{\Lambda }^2/2+h_0\Lambda +1},k^2=\frac{a^2{\Lambda }^2}{4{\Omega }^2},$$

Full text

I Introduction

General experience suggests that, after a sufficiently long time, an isolated system approaches thermal equilibrium, characterized by its initial energy but not by other details of the initial state. While occurrence of thermalization is widespread, it is also well established that certain classes of systems defy thermalization and preserve a more detailed memory of their initial conditions. Prominent examples are integrable systems, possessing an extensive number of (quasi)local conserved quantities Rigol et al. (2007); Vidmar and Rigol (2016), or many-body localized systems, where disorder effectively generates conserved quantities Nandkishore and Huse (2015); Gogolin and Eisert (2016). When a system is not isolated but coupled to an infinite bath, thermalization to a Gibbs state characterized by the temperature of the bath seems even more plausible, but a general proof of such behavior is lacking.

In the absence of a general proof we study a specific model that, due to integrability, does not thermalize when isolated from its environment, and we investigate whether a bath can indeed facilitate thermalization. The model we consider was introduced by Lipkin, Meshkov, and Glick (LMG) Lipkin et al. (1965) in the context of nuclear physics, and has subsequently found applications also in other branches of physics. The model can be written either as a bosonic Hamiltonian or, via a Schwinger mapping, in terms of pseudo-spin operators Lipkin et al. (1965). Our interest in this model is motivated by experiments with ultracold atoms for which it was shown that, under suitable circumstances, a two-component Bose–Einstein condensate can be described by an LMG model Zibold et al. (2010); Strobel et al. (2014). In the experimental realizations, the bosonic LMG Hamiltonian is restricted to a Fock space sector of fixed boson number $N$, corresponding to a fixed spin quantum number $S=N/2$ in the pseudo-spin representation (1) of the model. Because of this restriction the Hilbert space dimension of the system does not grow exponentially, but only linearly with the system size $N$ or, equivalently, the spin quantum number $S$. A number of unconventional thermodynamic properties have been shown to arise from this subexponential (linear) growth of the Hilbert space, including distribution functions that do not concentrate in the large-$S$ limit Webster and Kastner (2018).

In the abovementioned experimental realizations, the LMG model is well isolated from its environment and evolves under unitary dynamics. Since the model is solvable via Bethe Ansatz Ortiz et al. (2005), its equilibrium is in general not described by a thermal state. An experimental realization of a bath that couples to the LMG system in a controlled way, even though possibly a difficult task, may be feasible in one way or the other Webster and Kastner (2018). In this paper we theoretically model and analyze this situation by considering an LMG model coupled to an infinitely large bath of independent bosons, characterized by some spectral density. This bath, and also the form of the system–bath coupling, is not necessarily a faithful modeling of a specific future experimental realization. However, on the theoretical side, certain qualitative features of the dynamics are believed to be rather insensitive to the details of the bath. This creates the expectation that what we learn in the context of the specific bath introduced in Sec. II might be valid for a broader class of environments, hopefully including those that are relevant for potential future experimental realizations of an LMG system coupled to a bath.

In this paper, we use an LMG model (the system) weakly coupled to a bosonic bath (the bath) as a starting point for deriving, under fairly standard assumptions Breuer and Petruccione (2007), a Markovian master equation of Lindblad form Lindblad (1976) that describes the time evolution of the system’s density operator, while the bath degrees of freedom have been traced out. The main body of work that we want to report in this paper then consists of analyzing the properties of this master equation, and in particular the questions of whether a thermal state is approached in the long-time limit and, if so, at what rate. We address these questions by using the adjoint master equation, which governs the Lindblad time evolution of operators in the Heisenberg picture. By exploiting the presence of widely separated time scales, the main technical achievement of this paper is an asymptotic solution of those time-evolution equations for a relevant set of observables in the large-system limit.

Our main finding is that the presence of a bosonic bath indeed induces thermalization of the integrable LMG model for the expectation values of the system Hamiltonian density as well as of all three spin component densities $S_{x}/S$, $S_{y}/S$, and $S_{z}/S$, independently of the initial state of the system. Moreover, our results confirm, from an open-systems perspective, an unconventional aspect of the equilibrium ensemble calculation of Ref. Webster and Kastner (2018), namely the temperature independence of thermal expectation values in the large-$S$ limit. In particular, when taking the infinite-system limit at any fixed bath temperature $T$, the expectation values of the Hamiltonian and spin component densities equilibrate to their ground state values. This unconventional thermal behavior is a consequence of the subexponential growth of the Hilbert space dimension of the LMG model with $S$, which essentially prevents entropic effects from influencing the thermal properties; see Webster and Kastner (2018) for details. In addition to establishing thermalization, we calculate in Sec. IV.3 the dissipation rate at which thermal equilibrium is reached, which is found to depend on the (quantum) phase to which the equilibrium state belongs. To induce thermalization to an equilibrium state which yields expectation values different from that of the ground state, it is necessary to take the infinite-system limit while keeping the rescaled temperature $\widetilde{T}\equiv T/S$ fixed. We consider this case in Sec. V and provide evidence of thermalization on an observable-independent level by showing that a canonical Gibbs state is the unique stationary state of the Lindblad master equation.

II Model

The LMG model was proposed by Lipkin, Meshkov, and Glick in 1965 as a Hamiltonian that describes interacting bosons. In that same paper it was also recognized that, by means of a Schwinger mapping, the bosonic Hamiltonian can be expressed in terms of spin operators. Here we consider a slightly less general version of that model, which in spin language can be written as

[TABLE]

where $S_{x}$ and $S_{z}$ are the $x$- and $z$-components of a spin-$S$ vector operator, and $\Lambda$ is a coupling constant. The number of particles $N$ in the bosonic formulation of the model is related to the spin quantum number $S$ in the spin formulation via $N=2S$. In later sections of this paper we will consider the LMG model in the limit $S\to\infty$ of infinite boson number or, equivalently, the semi-classical limit of infinite spin quantum number $S$. This model exhibits a quantum phase transition at $\Lambda=1$ from the symmetric to the symmetry-broken phase. In the former the model has a non-degenerate ground state with a zero $S_{x}$ expectation value. In the limit $N\to\infty$ the ground state in the symmetry-broken phase will be two-fold degenerate with non-zero $S_{x}$ expectation values Dusuel and Vidal (2005).

We consider the prototypical environment of a bosonic bath Hamiltonian

[TABLE]

where $a_{q}$ and $a_{q}^{\dagger}$ are bosonic ladder operators and $\omega_{q}$ are the mode frequencies. The total system is described by the Hamiltonian

[TABLE]

where

[TABLE]

with

[TABLE]

couples the system to the bath. The coupling strength between system and bath is determined by $\gamma$, which we will later on consider as small. For specific parameter values, $\mathcal{H}$ reduces to models previously studied in the literature. For $S=1/2$ one obtains the spin–boson model Leggett et al. (1987), while for $\Lambda=0$ one obtains the multi-modal Dicke Hamiltonian Hepp and Lieb (1973).

The form of the coupling operator (4) is a convenient choice for the calculations. We expect, but have not shown, that the main conclusions of this paper remain valid for other coupling directions, like an $S_{y}$ operator in (4). An exception is a coupling operator (4) with $S_{z}$ in the place of $S_{x}$. In this case the total Hamiltonian conserves the parity symmetry of the spin, implying that the bath can destroy only some of the integrals of motion of the system, but not all of them, and hence thermalization cannot be expected to occur.

III Lindblad master equation

To arrive at a more tractable problem, we use standard approximations, detailed in the Appendices A and B, to derive a Born-Markov master equation of Lindblad form that approximately describes the time evolution of the density operator $\rho_{\text{S}}$ of the LMG model weakly coupled to a bath,

[TABLE]

with

[TABLE]

The entries of the $\kappa$ matrix are given in Appendix B. As part of the approximations when deriving the Lindblad master equation, the bath is at all times assumed to be in the Gibbs state $\rho_{B}\propto\exp\left(-\mathcal{H}_{B}/T\right)$ with temperature $T$ as well as uncorrelated with the system. Equation (6) describes the time evolution of the state (or density operator) of an open quantum system, analogous to the Schrödinger equation for a closed system. For closed systems it is well-known that, instead of time-evolving the state, an equivalent description of the quantum dynamics is possible in the Heisenberg picture, where operators evolve and states are constant in time. The same change of viewpoint is possible for open quantum systems, where the adjoint master equation

[TABLE]

describes the time evolution of an observable $A$ under the Lindblad dynamics of an open quantum system. Here $\mathcal{L}^{\dagger}$ denotes the adjoint Lindbladian Breuer and Petruccione (2007)

[TABLE]

which satisfies $\operatorname{{Tr}}[A(\mathcal{L}\rho_{S})]=\operatorname{{Tr}}[(\mathcal{L}^{\dagger}A)\rho_{S}]$. Equations (8) and (III) serve as the starting point for deriving the spin components’ equations of motion in the Heisenberg picture. While this derivation is a standard procedure for closed quantum systems, the non-unitary nature of the dynamics generated by the Lindbladian introduces additional technical complications here. To illustrate this, we write $A(t)$ as $A(t)=\exp(t\mathcal{L}^{\dagger})A$, in which case the adjoint master equation becomes

[TABLE]

where $\mathcal{L}^{\dagger}A$ is evaluated at $t=0$. While determining $\mathcal{L}^{\dagger}A$ is straightforward, the subsequent application of $\exp(t\mathcal{L}^{\dagger})$ will not simply transform all the operators appearing in $\mathcal{L}^{\dagger}A$ into the Heisenberg picture. This is a consequence of the Lindbladian not satisfying the Leibniz rule, i.e. in general

[TABLE]

and so

[TABLE]

This complicates the task of extracting the desired equations of motion from the adjoint master equation. In what follows we show that, in the large-$S$ limit, the Lindbladian may be approximated to restore the Leibniz rule, without neglecting the dissipative effects of the bath coupling. Within this approximation the equation of motions can then be obtained easily.

Our approximation of the adjoint master equation in the $S\rightarrow\infty$ limit is based on the scaling behavior in $S$ of the various terms appearing in (III). Within each term, every spin component contributes a factor of $S$ to the overall scaling, while every commutator contributes $S^{-1}$. This is a consequence of the fundamental commutator

[TABLE]

scaling linearly, rather than quadratically, in $S$. The general result for the commutator of arbitrary products of spin components then follows from using the Leibniz rule to reduce the calculation to the repeated application of (13). Since all the operators of interest are derived from these components, we can use this simple rule to identify the dominant terms in (III). After accounting for the $\mathcal{O}(S^{-1})$ scaling of the $\kappa$-coefficients, we find that the final three terms on the right of (III) are subdominant relative to the first two. Dropping these subdominant terms yields the approximation

[TABLE]

for which it holds that

[TABLE]

The same scaling argument as before reveals that the term in the second line of (15) is subdominant in $S$. In the large-$S$ limit we may therefore take $\mathcal{L}^{\dagger}$ to satisfy the Leibniz rule, and so $e^{t\mathcal{L}^{\dagger}}$ will satisfy $e^{t\mathcal{L}^{\dagger}}(AB)\approx(e^{t\mathcal{L}^{\dagger}}A)(e^{t\mathcal{L}^{\dagger}}B)$. Combining this with (14) and inserting $\operatorname{{Im}}(\kappa_{yx})=\gamma/(2S)$ turns the adjoint master equation in (10) into

[TABLE]

with all the operators now time-evolved in the Heisenberg picture. In Sec. IV we will analyze the dynamics that follows from this approximated form of the master equation.

Note that when identifying the dominant terms in the original Lindbladian (III) we assumed that the temperature $T$ is intensive, and does not scale with system size. This resulted in the $\mathcal{O}(S^{-1})$ scaling of the $\kappa_{xx}$ coefficient, which led to the third term on the right of (III) falling way. Allowing for the temperature to scale extensively with $N$ will impact on both the dynamics of the observables and the equilibrium state thermodynamics. We return to this discussion in Sec. V.

IV Markovian dynamics of spin observables

In this section we study the dynamics, generated by the approximated adjoint master equation in (16), of the rescaled spin observables

[TABLE]

We begin by deriving the equations of motion, which turn out to describe dynamics on two widely separated time scales. Using the two-timing method then allows us to analyze the equilibrium states and the energy dissipation rate at late times.

IV.1 Equations of motion

Inserting the three spin components into (16) and using the spin commutation relations $[S_{\mu},S_{\nu}]=i\epsilon_{\mu\nu\delta}S_{\delta}$ leads to the equations of motion

[TABLE]

In deriving these expressions we have used the fact that the commutators of the rescaled spin components are of order $\mathcal{O}(S^{-1})$, so that operator order does not matter. In the thermodynamic limit we can therefore treat $x$, $y$ and $z$ as commuting scalar variables which are constrained to the unit sphere as

[TABLE]

which is a reflection of the Casimir constraint $\bm{S}^{2}=S(S+1)$.

In the symmetric phase ($\Lambda\leq 1$) of the LMG model, the only stationary solution to (18a)–(18c) is $(x,y,z)=(0,0,\pm 1)$, while in the symmetry-broken phase ($\Lambda>1$) there are two additional stationary points at $(x,y,z)=(\pm\sqrt{1-\Lambda^{-2}},0,\Lambda^{-1})$. To analyze how the dissipative dynamics steers the system towards these points we first introduce the rescaled system energy

[TABLE]

which, together with $x(t)$, obeys the set of coupled equations

[TABLE]

From the expression for $\dot{h}$ we see that, as a result of the bath coupling, the system’s energy decays on a time scale set by $\gamma$. For weak coupling, this will be much longer than the time scale associated with the dynamics of the individual spin components. This separation of time scales suggests the use of the two-timing method Strogatz (2018), which we implement in the next section. We focus for now on the evolution of $h$ and $x$, and return to the $y$ and $z$ components at the end.

IV.2 Two-timing method

To capture the dependence of the observables on the two time scales that emerge at weak coupling, we introduce a fast time $\tau=t$ and a slow time $s=\gamma t$. These will be treated as independent variables. As an Ansatz for the solutions of (21a) and (21b) we take

[TABLE]

Inserting these forms into (21a) and (21b) and collecting powers of $\gamma$ leads to

[TABLE]

The first equation shows that $h_{0}(\tau,s)=h_{0}(s)$ is independent of $\tau$, and that the system’s energy remains constant over the short time scale associated with $\tau$. This also implies that in (23b) the $\tau$-dependence only enters through $x_{0}(\tau,s)$. The latter is the equation of motion for a particle moving in the effective quartic potential

[TABLE]

The solution for $x_{0}$ is given by Olver et al. (2010)

[TABLE]

with

[TABLE]

where $\operatorname{{cn}}$ denotes the Jacobi elliptic cosine function. Note that $a$ and $u$, which enter as integration constants in the solution for $x_{0}(\tau,s)$, are functions of $s$. To determine $a(s)$ we combine (18a), (19), and (20) into an additional constraint on $x_{0}(\tau,s)$, which reads

[TABLE]

This produces an equation for $a^{2}$ with solution

[TABLE]

which depends on $s$ through $h_{0}$. We may take $a$ to be positive, since a change in sign can be absorbed through an appropriate shift in $u$.

At this stage we have established the $\tau$-dependence of $x_{0}(\tau,s)$ and seen that the $s$ dependence enters via $h_{0}(s)$ and $u(s)$. We now turn to (23c), from which we extract the $s$-dependence of $h_{0}(s)$ in order to characterize the energy dissipation. First we average both sides of (23c) with respect to $\tau$ over the interval $[0,\tau_{f}]$, yielding

[TABLE]

A crucial element in the implementation of the two-timing method is the elimination of secular terms which grow in $\tau$ without bound. Here this amounts to requiring that $h_{1}(\tau_{f},s)/\tau_{f}\rightarrow 0$ as $\tau_{f}\rightarrow\infty$. Taking this limit in (29) eliminates the first term on the right-hand side and turns the second term into the long-time average of $(\partial_{\tau}x_{0})^{2}$. The latter can be replaced by the average over a single period of the periodic solution for $x_{0}(\tau,s)$ in (25). This produces

[TABLE]

where the period of $x_{0}(\tau,s)$ is

[TABLE]

and $K(k^{2})$ denotes the complete elliptic integral of the first kind Olver et al. (2010). The averaging in (30) also eliminates any dependence on $u(s)$, and the result depends on $s$ only through $h_{0}(s)$. We write (30) compactly as

[TABLE]

where the result of the integral is

[TABLE]

Note that $k$ is itself a function of $h_{0}$ through (26). Equations (32) and (33), which are the main results of this subsection, describe the slow dissipation of the system’s energy. Next we will reintroduce the operator character of $h_{0}(s)$ and study equilibration on the basis of Eq. (32).

IV.3 Approach to equilibrium

Since no other operators apart from $h_{0}(s)$ appear in (32), $h_{0}(s)$ will remain diagonal in the original $h(0)$ basis. In this sense (32) describes the evolution of the eigenvalues of $h_{0}(s)$. Any such eigenvalue $\epsilon(s)$ satisfies

[TABLE]

and the behavior of the eigenvalues at late times will reflect the system’s equilibrium state. The possible initial eigenvalues $\epsilon(0)$ are determined by the range of the spectrum of $h(0)=-\Lambda x^{2}/2-z$. In the large-$S$ limit this range can be obtained by treating the rescaled spin components as commuting variables and finding the extremal values of $h(0)$ subject to the unit sphere constraint (19). It is found that $\epsilon(0)\in[\epsilon_{g},1]$ where

[TABLE]

The same result follows from the standard variational approach based on spin coherent states Dusuel and Vidal (2005).

Figure 1 shows $A(\epsilon)$, which governs the evolution of $\epsilon(s)$, as a function of $\epsilon\in[\epsilon_{g},1]$ for values of $\Lambda$ corresponding to the two phases. We see that $-A(\epsilon)$ is negative everywhere, except at the boundaries of the spectral range, and in the symmetry-broken phase also at $\epsilon=-1$. In the symmetric phase the eigenvalues starting in $[-1,1)$ therefore evolve towards $\epsilon_{g}=-1$, which is a stable fixed point of the evolution. In the symmetry-broken phase eigenvalues in $[-1,1)$ converge to $\epsilon=-1$, while those in $[\epsilon_{g},-1)$ tend to $\epsilon_{g}$. The appearance of the fixed point at $\epsilon=-1$ for $\Lambda>1$ is an artifact of the large-$S$ limit and the expansion in $\gamma$. It corresponds to the scenario where the dissipation lowers $h_{0}(s)$ to precisely the value where $x_{0}(\tau,s)$, evolving according to (23b), gets “stuck” at the local maximum of the double-well potential $V(x)$ in (24). While this is inevitable when considering only the dynamics of $h_{0}$ and $x_{0}$, it clearly requires very specific initial conditions to occur in the original equations of motion of $h$ and $x$ in (21a) and (21b). The $\epsilon=-1$ fixed point is therefore not physically significant, and will not halt the system’s evolution at finite $S$, nor when the coupling $\gamma$ is treated non-perturbatively. The relevant fixed point in the symmetry-broken phase is therefore $\epsilon_{g}$, which serves as an attractor for all initial eigenvalues in $[\epsilon_{g},1)$. Returning to the operator level, we conclude that $\lim_{s\rightarrow\infty}h_{0}(s)=\epsilon_{0}\hat{I}$. Combining this with (20) and (27) implies the equilibrium expectation values

[TABLE]

and $\langle y_{0}\rangle_{\text{eq}}=0$. These coincide with the stable equilibria of (18a)–(18c) identified earlier.

To determine the dissipation rate close to equilibrium we consider the expectation value of $h_{0}(s)$ with respect to an arbitrary initial state. This expectation value evolves according to

[TABLE]

With all the eigenvalues of $h_{0}(s)$ tending towards $\epsilon_{0}$, we can at late times approximate $\langle-A(h_{0}(s))\rangle$ by $-A(\langle h_{0}\rangle(s))$. The behavior of $A(\epsilon)$ around $\epsilon=\epsilon_{g}$ is

[TABLE]

Close to equilibrium we therefore have

[TABLE]

where

[TABLE]

is the phase-dependent dissipation rate. Note the switch from rescaled time $s=\gamma t$ back to $t$. Solving (39) for $\langle h_{0}\rangle(t)$ yields

[TABLE]

with $C$ a constant.

V Extensive temperatures

We have shown that at any fixed finite temperature the equilibrium expectation values of intensive quantities match their ground state values in the large-$S$ limit. This does not imply thermalization to the ground state, but rather that these quantities cannot resolve the difference between the ground state and the true mixed equilibrium state. To observe equilibration towards states for which, even in the large-$S$ limit, these expectation values deviate from those of the ground state requires, as we will see, the temperature $T$ to be considered an extensive quantity. This is equivalent to taking the large-$S$ limit at a fixed rescaled temperature $\widetilde{T}=T/S$. In the Lindbladian (III) this amounts to replacing $\kappa_{xx}$ by $2\gamma\widetilde{T}$ 111We consider $\nu_{1}$, and hence $\kappa_{xy}$ and $\kappa_{yx}$, as asymptotically independent of $S$. Since $\nu_{1}$ in (60) grows proportional to $T/\omega_{c}$, this essentially amounts to rescaling the cut-off frequency $\omega_{c}\to S\omega_{c}$ of the bath spectral density (54)., yielding

[TABLE]

Compared to Eq. (16), the last three terms on the right-hand side make the master equation (42) more complicated to deal with, and we were not able to analyse the long-time dynamics along the lines of Sec. IV. As an alternative approach, we show in this section analytically that a Gibbs thermal state is a stationary state of this master equation in the large-$S$ limit, and provide numerical evidence that no other stationary states exist. Trace preservation of the master equation then implies that any initial state evolves towards the unique stationary state, which establishes thermalization.

We consider the Gibbs thermal state

[TABLE]

with $h=\mathcal{H}_{S}/S$, where $\widetilde{\beta}=S\beta=1/\widetilde{T}$ is the rescaled inverse temperature and $Z=\operatorname{{Tr}}[\exp(-\beta{\cal H}_{S})]$ is the partition function. To establish the stationarity of the master equation (6) in the large-$S$ limit we compare the scaling in $S$ of $\dot{\varrho}_{\widetilde{\beta}}={\cal L}\varrho_{\widetilde{\beta}}$ with that of $\varrho_{\widetilde{\beta}}$ itself. Here the linearity of ${\cal L}$ allows us to replace $\varrho_{\widetilde{\beta}}$ by $e^{-\widetilde{\beta}h}$ for the sake of this comparison. We will show that, while $e^{-\widetilde{\beta}h}={\cal O}(S^{0})$, a delicate cancellation of leading order terms results in ${\cal L}e^{-\widetilde{\beta}h}$ being of order ${\cal O}(S^{-1})$. In this sense the relative rate of change of $\varrho_{\widetilde{\beta}}$ becomes negligible in the large-$S$ limit.

Inserting $e^{-\widetilde{\beta}h}$ into (42), the first commutator on the right-hand side of (42) is zero. The arguments used to identify the dominant terms in (III) for the adjoint Lindbladian, which eventually led to (14), can likewise be employed to determine the scaling behaviour of the remaining terms in (42). The final two terms are of order ${\cal O}(S^{-2})$ and ${\cal O}(S^{-1})$, respectively, while the second and third terms are individually of order ${\cal O}(S^{0})$. As we show next, when $\mathcal{L}$ acts on $e^{-\widetilde{\beta}h}$ there are cancellations between the latter two terms which result in their combined contribution being of order ${\cal O}(S^{-2})$. For each of these terms, we first combine the prefactor with the innermost (anti)commutator, which can be expanded in a power series in $1/S$ using the BCH formula. For the second term in (42) this yields

[TABLE]

where the omitted terms involve further nested commutators with $h$, each of which contributes a factor of $S^{-1}$ to the term’s overall scaling. Applying the same procedure to the third term in (42) produces

[TABLE]

which, up to order ${\cal O}(S^{-1})$, differs from the previous result only by a sign. Inserting (44) and (45) into the sum of the second and third term of (42) then yields a combined contribution of the form $[S_{x},{\cal O}(S^{-2})e^{-\widetilde{\beta}h}]$, which is itself of order ${\cal O}(S^{-2})$. Taken together, these arguments establish the ${\cal O}(S^{-1})$ scaling of ${\cal L}e^{-\widetilde{\beta}h}$ and therefore also the stationarity of $\varrho_{\widetilde{\beta}}$ in the large-$S$ limit.

The uniqueness of the stationary state $\varrho_{\tilde{\beta}}$ is encoded in the spectrum of the Lindbladian (42). Denoting the eigenvalues of $\mathcal{L}$ by $\lambda_{i}$ and the corresponding eigenstates by $\varrho_{i}$, it follows that $e^{t{\cal L}}\varrho_{i}=e^{t\lambda_{i}}\varrho_{i}$. Stationary states correspond to $\lambda_{i}=0$, while states with $\operatorname{{Re}}\lambda_{i}<0$ decay in time. We order the eigenstates by their real part such that $\operatorname{{Re}}\lambda_{i}\geq\operatorname{{Re}}\lambda_{i+1}$ for all $i$. The stationary state associated with $\lambda_{0}=0$ is unique if $\operatorname{{Re}}\lambda_{1}<0$. Using the Quantum Toolbox in Python Johansson et al. (2013) we calculated the gap $\lambda_{0}-\lambda_{1}$ of the Lindbladian for system sizes up to $S=80$. We find a nonzero gap that is essentially independent of $S$, supporting the conclusion that the gap persists in the large-$S$ limit. Accepting this conclusion as a fact, it follows that the Gibbs thermal state (43) is the unique stationary state, implying thermalization under the dynamics of the master equation (6) with Lindbladian (42). In Fig. 2 this conclusion is corroborated by numerical calculations of the time evolution of the energy expectation value $\left\langle h\right\rangle$, showing convergence to the predicted thermal equilibrium value independently of the initial state.

VI Conclusions

We derived a Lindblad master equation that approximates the dynamics of an LMG model weakly coupled to a bosonic bath of temperature $T$. The derivation mostly follows the standard procedure of performing Born and Markov approximations Breuer and Petruccione (2007), but avoids the secular approximation, which is usually the final step of the derivation. Avoiding this step results in a time evolution equation that is not completely positive. We restore complete positivity by a minimal modification of the jump rates, see Eq. (68) and below, to obtain a master equation that satisfies the requirements of a quantum dynamical semi-group.

In the absence of a bath, the LMG model is integrable and does not equilibrate in general. On the basis of our Lindblad master equation we analyzed whether the presence of a bath indeed enforces, as on would expect, thermalization, i.e., equilibration to a Gibbs state of temperature $T$ independently of the system’s initial state.

Our first main result establishes that, for an LMG model in the large-system (or large-$S$) limit, thermalization occurs in the sense that expectation values of the system’s Hamiltonian density $\mathcal{H}_{S}/S$ as well as expectation values of all three spin component densities $S_{x}/S$, $S_{y}/S$, and $S_{z}/S$ converge in the long-time limit to their thermal equilibrium values. A peculiar property of those thermal equilibrium values is that they agree, for all values of the bath temperature $T$, with the ground state expectation values of the LMG model. This peculiarity of the model was previously discussed in Webster and Kastner (2018) in the equilibrium context, and is confirmed in a dynamical context by the results of the present paper.

The unconventional thermal behavior of the LMG model is a consequence of the small Hilbert space dimension of the model. While in conventional spin models the Hilbert space dimension grows exponentially with the system size $N$, it grows only linearly with $N$ (or $S$) in the case of the LMG model. Due to the small number of available states, entropic effects are mostly irrelevant and the effect of temperature, which is conjugate to the entropy, is suppressed.

A strategy to force the system into highly excited thermal equilibrium states even in the absence of entropic effects consists in promoting temperature to become an extensive quantity. This is achieved formally by substituting $T\to\tilde{T}S$ in the master equation and considering $\tilde{T}$ to be constant while taking the large-$S$ limit. This procedure renders terms in the master equation that were previously subleading in $1/S$ into leading ones. Analyzing the resulting equation, we found in Sec. V that a Gibbs thermal state with rescaled temperature $\tilde{T}$ is stationary. Establishing thermalization requires to show that this Gibbs state is the only stationary state. We were not able to prove uniqueness of the stationary state analytically, but provided numerical evidence in favor of this claim.

Rescaling temperature with $S$ may seem as an artificial and unphysical procedure, but this is not the case. The large-system limit, or large-$S$ limit, is an idealization, and one is ultimately interested in large but finite systems. Thermalization in rescaled temperature $\tilde{T}$ then states that a large system of size $S$ and (unscaled) temperature $T$ will have a thermal equilibrium state very similar to that of a system of size $2S$ at temperature $2T$, etc., and the limit of such a sequence is given by a Gibbs state with rescaled temperature.

All results of this paper have been proved for a master equation that was derived for a bosonic bath with a Lorentz-Drude spectral density (54) and a system–bath coupling (4). It is unclear at present how the specific properties of an environment might look like in an experimental realization of the LMG model in a bath, but one would hope that, as in other known cases, details of the bath are not overly important and our results carry over to a larger class of environments, including experimentally relevant ones. This would then open up an avenue for imposing thermalization on the otherwise non-thermalizing LMG model and for probing its unconventional thermodynamic properties.

Acknowledgements.

M. K. gratefully acknowledges helpful discussions with I. de Vega, hospitality at the Kavli Institute for Theoretical Physics at Santa Barbara during the QSIM19 program, and financial support by the South African National Research Foundation through the Incentive Funding Programme and the Competitive Funding for Rated Researchers. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

Appendix A Born-Markov master equation

For the derivation of the master equation we consider the von Neumann equation in the interaction picture,

[TABLE]

where $\mathcal{V}(t)=U(t)\mathcal{V}U^{{\dagger}}(t)$ with $U(t)=\text{exp}(i\mathcal{H}_{0}t)$. Assuming system and bath to be uncorrelated at $t=0$, system–bath interactions to be weak, and the bath to evolve much faster than the system, we employ a sequence of standard approximations Breuer and Petruccione (2007) to derive the Born-Markov master equation

[TABLE]

where h.c. denotes the hermitian conjugate. Transforming (47) back into the Schrödinger picture one obtains

[TABLE]

where we have defined the counter term

[TABLE]

and the coupling operator

[TABLE]

with

[TABLE]

It is common to simplify (48) further by performing a secular or rotating-wave approximation, which amounts to averaging over rapidly oscillating terms in the master equation. This approximation guarantees that the resulting master equation is a quantum dynamical semi-group Breuer and Petruccione (2007). However, the approximation also imposes, possibly artificially, a Gibbs state to be a stationary state of the dynamics Spohn and Lebowitz (2007); Breuer and Petruccione (2007). Further this approximation fails close to a phase transition due to the divergence in time scales Koplov and Schaller . Here we continue without a secular or rotating-wave approximation, and we prove in Appendix B that, after some tweaking, our master equation nonetheless satisfies the desired properties of a quantum dynamical semi-group.

To simplify the master equation (48) it is convenient to separate the bath correlation function (51) into real and imaginary parts. Assuming the bath to be in a Gibbs state with inverse temperature $\beta$, one can calculate the so-called decoherence kernel to be Schlosshauer (2007)

[TABLE]

and the noise kernel to be

[TABLE]

Irreversible dynamics, and hence equilibration in an asymptotic sense, can occur only for an infinite bath with a continuous frequency spectrum Breuer and Petruccione (2007). Formally we promote the set of discrete frequencies $\omega_{q}$ to a continuous variable $\omega$ and replace $\sum_{q}\gamma_{q}^{2}$ by $\int d\omega J(\omega)$ in (52) and (53). We choose a spectral density $J$ of Ohmic type with a Drude cut-off Breuer and Petruccione (2007),

[TABLE]

where the cut-off frequency $\omega_{c}$ should be chosen such that the dynamics becomes essentially independent of $\omega_{c}$. After substitution one finds

[TABLE]

for the noise kernel and

[TABLE]

for the decoherence kernel Breuer and Petruccione (2007), both of which decay exponentially with $\tau$. The corresponding decay rates are $\omega_{c}$ and $2\pi/\beta$, such that the bath correlation time is

[TABLE]

To simplify the master equation (47) we want to evaluate the integral (50). The integrand contains $S_{x}(-\tau)$, and while an exact solution of this time-evolved operator seems out of reach, the first few terms of the Baker-Campbell-Hausdorff formula give the approximate solution

[TABLE]

valid for small bath correlation times $\tau_{B}$. Using the approximate solution (58) together with (55) and (56) one can show that (50) simplifies to Breuer and Petruccione (2007)

[TABLE]

where higher-order terms in $\tau$ and $\gamma$ have been neglected. Here we have defined

[TABLE]

where $\psi$ denotes the digamma function Olver et al. (2010).

Using (59) we can now write the master equation (48) in the form

[TABLE]

with a unitary part

[TABLE]

where

[TABLE]

and a dissipative part

[TABLE]

Inserting (59) into (64) one finds that all terms containing $\omega_{c}$ cancel. This justifies to use in (64) the simpler coupling operator

[TABLE]

as it leads to an identical master equation.

Appendix B Lindbladian

In this section we make adjustments to the master equation (61) that impose the properties of a quantum dynamical semi-group onto the dynamics, which guarantees that a quantum mechanical density operator remains a density operator under time evolution for all times and initial states. Lindblad Lindblad (1976) has shown that this is guaranteed if ${\cal U}$ and ${\cal D}$ in (61) can be written as

[TABLE]

and satisfy the following conditions:

1. (i)

$\mathcal{H}$ is hermitian. 2. (ii)

$L_{k}$ are traceless and pairwise orthonormal operators with respect to the Hilbert-Schmidt inner product, $\operatorname{{Tr}}\{L_{k}^{{\dagger}}L_{l}\}=\delta_{kl}$. 3. (iii)

The matrix $\kappa$ with elements $\kappa_{kl}$ is hermitian and positive semi-definite.

For our master equation, (i) is clearly satisfied by (1) and (63). To satisfy (ii) we choose $L_{1}=S_{x}/S$ and $L_{2}=S_{y}/S$. To assess condition (iii) we write the dissipator (64) with coupling operator (65) as

[TABLE]

with

[TABLE]

The matrix $\kappa$ in (68) is Hermitian, but not positive semi-definite, as is evident from

[TABLE]

To restore positivity, we follow section 3.6.2.1 of Breuer and Petruccione (2007) and determine the minimal modification of $\kappa_{yy}$ that leads to $\det{\kappa}=0$. From (69) we read off that this is achieved by setting $\kappa_{yy}=|\kappa_{xy}|^{2}/\kappa_{xx}$. Modifying $\kappa_{yy}$ “by hand”, while not strictly justified, is expected to be reasonable at least when $|\kappa_{xy}|^{2}/\kappa_{xx}=(\nu_{1}^{2}+1/4)/T$ is small. Writing $\kappa$ in its spectral decomposition we note that the zero-eigenvalue term comes with a prefactor of zero. In that basis it is therefore sufficient to consider, instead of (67), the dissipator

[TABLE]

with a single jump operator

[TABLE]

which can be left unnormalized. After some tedious manipulations that lead to cancellations between the unitary and dissipative terms of the master equation, we obtain the final Born-Markov master equation (6) with Lindbladian (III), where $\kappa_{xx}$, $\kappa_{xy}$, and $\kappa_{yx}$ are defined as in (68).