Relaxation of the Ising spin system coupled to a bosonic bath and the time dependent mean field equation

Authors: Máté Tibor Veszeli, Gábor Vattay

Abstract

The Ising model does not have strictly defined dynamics, only a spectrum. There are different ways to equip it with time dependence, e.g., the Glauber or the Kawasaki dynamics, which are both stochastic, but it means there is a master equation that can also describe their dynamics. These equations can be derived from the Redfield equation, where the spin system is weakly coupled to a bosonic bath. In this paper, we investigate the temperature dependence of the relaxation time of a Glauber-type master equation, especially in the case of the fully connected, uniform Ising model. The finite-size effects were analyzed with a reduced master equation and the thermodynamic limit with a time-dependent mean field equation.

1 Introduction

Spin models are versatile because they are simple, yet able to demonstrate fundamental phenomena, like phase transition [1–3]. The spin models originated from solid state physics, where the interaction between the electron spins can be due to direct exchange [4], indirect exchange [5], superexchange [6] or double exchange [7], but many other complex physical systems can be modelled using a simple Ising or Heisenberg model, like nuclear spins [8, 9], and even social situations [10]. It is also important in modern applied physics since one branch of adiabatic quantum computers—e.g. the D-Wave system [11]—are based on finding the global minimum of an artificial spin system [12, 13].

The Ising model is defined via its energy or, in the quantum case, where it is called Heisenberg model via its Hamiltonian operator. All the equilibrium properties can be determined from these quantities, e.g. heat capacity, magnetic moment, susceptibility etc., since they can be derived from the partition function. On the other hand to calculate inequilibrium properties like the relaxation time, we need to know the dynamical equation of the system.

The Ising model does not have a natural dynamics, and although to the Heisenberg model we can associate the Schrödinger or the Heisenberg equation, it will only generate a unitary time evolution. Therefore it cannot be responsible for a final, thermal distribution. For this we need some interaction with the surrounding environment.

Close to thermal equilibrium it is often assumed, that the dynamics is driven by the gradient of a mean field free energy, since at equilibrium it must be zero [14–17]. This phenomenological approach can account for some dynamical critical behaviour, like the dynamical slowing down, but like every mean field approximation it is valid only if the fluctuations are small. For example in one dimension a master equation gives a better description [18].

A more fundamental approach is to derive the effective dynamical equations from a system-plus-bath model, where the bath is a thermal reservoir. This method let us track how the parameters of the macroscopic equation depends on microscopic properties. To describe such a system we must use the tools of open quantum systems [19, 20] like the Redfield [21] and the Lindblad equation [22]. These equations have countless applications in quantum biology [23, 24], quantum optics [19], cold atomic gases [25], chemical physics [26] as well as it also being relevant in quantum computing [27–29].

In solid matter the interaction with the phonons is always present and the electrons, as charged particles are coupled to the photons, thus we can assume that the spin system is in a bosonic bath. Quantum dissipation and relaxation of spin systems in a bosonic bath as well as in magnetic field have been investigated by many authors. [30–33].

In this paper the interaction between an adiabatic computer and its environment is meant to be small, so the weak coupling Lindblad equation is be used, but of course, there are improved methods to describe open quantum systems, like slippage initial condition [34, 35], the Nakajima-Zwanzig equation [36, 37] or the polaron transformation [38].

The structure of this paper is the following. In section 2 we briefly summarize the derivation of a Glauber-type master equation [39] based on the Redfield equation. Using this microscopic approach we can see how the Fourier transform of the bath correlation function appears in the final master equation. In section 3 we investigate the temperature dependence of the eigenvalues of the transition matrix because they contain relevant information on the time scales of the system, e.g., the relaxation time. We also give an upper bound to the smallest nonzero eigenvalue. In section 4 the dynamics of the uniform, fully connected Ising model is investigated, and we show that the relaxation time diverges in the thermodynamic limit as the temperature approaches the critical temperature. In section 5 a time dependent mean field equation is revisited, which is used in section 6 to extend the previous results to infinite sizes.

2 Master equation of quantum Ising system

In general if a system is connected to a bath, then its Hamiltonian operator is Htot=H+HB+HI, (1) where H acts only on the system of interest, HB only on the bath, and HI is the interaction between the two subsystems, it can be written as HI = ∑α Aα ⊗ Bα, where Aα and Bα are system and bath operators respectively. The dynamics of the total system is described by the von Neumann equation. ρtot=−i[Htot,ρtot] (2) After the Born and the Markov approximation an effective equation can be derived for the density operator of the system of interest. (ρ≔TrBρtot). ρ˙(t)+i[H,ρ(t)]=∑α(Aαρ(t)Tα†−AαTαρ(t)+h.c.), (3) where Tα=∑β∫0∞dtCαβ(t)AβI(−t), AβI(t) is in interaction picture, Cαβ(t)=〈BαI(t)Bβ〉B is the bath correlation function, and h.c means hermitian conjugate. This is the Redfield equation in weak-coupling limit [21]. The Born approximation is valid if the interaction between the system and the bath is small, and we can use the Markov approximation if the relaxation time of the system of interest is much larger than the decaying time of the bath correlation function. After the so called secular or rotating wave approximation [40] one can get to the Lindblad equation [19, 22]: ρ˙+i[H+HLS,ρ]=∑αβ∑ωγαβ(ω)(Aβ(ω)ρAα†(ω)−12{Aα†(ω)Aβ(ω),ρ}), (4) where Aα(ω)=∑ij|i〉〈i|Aα|j〉〈j|δω,εj−εi and |i〉 is the eigenvector of H with eigenvalue εi. HLS is the Lamb shift Hamiltonian, and γαβ(ω) is the Fourier transform of the bath correlation function. γαβ(ω)≔∫−∞∞dteiωt⟨Bα†(t)Bβ(0)⟩B (5) There are two standard bosonic the bath of phonons and the bath of photons. For phonons γohm(ω)∼e−|ω|ωcω1−e−βω, which is called Ohmic case and for photons γsup(ω)∼ω31−e−βω, which is a super-Ohmic case. The frequency ωc is the cutoff frequency. If we assume, that ωc is large compared to the energy distances of the system, then e−ω|ωc|≈1. Fig 1 shows the main features of the two γ functions. The main difference is that γohm is strictly increasing and γohm(ω = 0) ∼ kBT, but γsup is non-monotonic, and γsup(0) = 0. In the easiest case γαβ ∝ δαβ.

The advantage of the weak-coupling limit is that a master equation can be derived for the diagonal elements of ρ. P˙i=∑jMijPj≡∑jWijPj−∑jWjiPi, (6) where Pi = ρii, Wij = ∑αβ γαβ(ωji)(Aα)ji(Aβ)ij and ωji = (εj − εi).

The system converges to the Boltzmann distribution if Wij satisfies the detailed balance condition i.e., Wij = Wji exp(−β(εi − εj)). Both the Ohmic and the super-Ohmic bath satisfy it, because γ(−ω)=e−βωγ(ω) (7)

If the system of interest is the Ising model, then a master equation can be derived [41], which contains only the diagonal elements of the density matrix. The Hamiltonian is H=−∑ij(i<j)Jijσizσjz−∑ihiσiz, (8) where σiz is the Pauli z-matrix and the corresponding eigenvectors are |S_⟩≡|S1,S2,…,SN⟩Si∈{±1} (9) with eigenenergies ES_=−∑ij(i<j)JijSiSj−∑ihiSi (10) The easiest way to couple the system to the bath is via a Pauli matrix i.e., Aα→σix. Using σ^z^ in the interaction instead of σ^x^ would not give any relevant dynamics since the system and the interaction Hamiltonians would commute. The peculiarity of this system is that the populations decouple even without the secular approximation.

The σix operator acting on |S〉 only flips the ith spin, so the WSS′ matrix element is WS_S_′=∑iγ(ωS_′S_)(σix)S_′S_(σix)S_S_′={γ(ωS_′S_)|iftheHammingdistancebetweenS_andS_′is10|otherwise. (11) With Eqs (6) and (11) we have a dynamics for the Ising model. P˙S_=∑S_′MS_S_′PS_′, (12) where MS_S_′=WS_S_′−δS_S_′∑S_″WS_″S_′ is the transition matrix. This matrix is temperature dependent, and it has at least one zero eigenvalue, which is the eigenvalue of the equilibrium PS_eq=e−βES_Z (13) For constant temperature the general solution of (12) is PS_(t)=∑S_′∑μe−λμtPμ,S_RPμ,S_′LPS_(t=0), (14) where PμRs are the right, and PμLs are the left eigenvectors of M with −λ~μ~ eigenvalues. All the λ~μs are nonnegative. If the system is ergodic, then there is only one zero eigenvalue, and the other λs are positive. Let the smallest positive be λmin~ and the largest be λmax. The relaxation time is tr = 1/λmin. This is the time scale in which all but the equilibrium mode dies out. The other relevant time scale is 1/λmax, which is the characteristic time of the fastest mode. If, for example, this spin system is a quantum computer, then the fastest mode is more important because if the computation is slower than this time scale, then the environment is not negligible. In other words, λmin is important if we want the system to relax thermally, and λmax is important if we want to avoid any thermal influence.

3 Temperature dependence of the eigenvalues

Both the smallest and the largest eigenvalues carry relevant information, and since M(β) is temperature dependent λmin(β) and λmax(β) are too.

At high temperatures we can determine the temperature dependence of all λs by simply Taylor expanding γ(ω; β) for small β. γ(ω;β)=ηωα1−e−βω≈ηωα−1β, (15) where α = 1 in the Ohmic, and α = 3 in the super-Ohmic case. The transition matrix inherits this temperature MSS′(β) ∼ β^−1^, and hence λ ∼ β^−1^.

Despite the high-temperature limit, where the elements of the dynamical matrix M diverges, in the low-temperature limit they converge. limβ→∞γ(ω;β)={0|ω≤0ηωα|ω>0. (16) It means all the eigenvalues also converge. As a consequence, we cannot slow down arbitrary all the modes by reducing the temperature. We have an upper limit in time for quantum computing. Of course, this calculation is valid only for a time-independent system, but the main features apply to more general cases.

Without an external magnetic field (h = 0) at zero temperature, the equilibrium Boltzmann distribution prefers only the two spin configurations with the lowest PS_eq(T=0)=12(δS_,S_g+δS_,−S_g), (17) where Sg and −Sg are the ground states. At zero temperature, there is one more eigenvector with a zero Pmin,S_R=12(δS_,S_g−δS_,−S_g) (18) The question is how λmin(β) behaves at low temperature. We can give an upper bound. First let us introduce the following symmetric M˜S_S_′=MS_S_′PS_′eqPS_eq≡MS_S_′e−βES_′−ES_2 (19) This transformation doesn’t affect the eigenvalues, and the eigenvectors transform like P˜μS_=Pμ,S_RPS_eq. (20) Since M˜ is symmetric its right and left eigenvectors are the same, and now the variational method applies to λmin≤−∑S_,S_′Π˜S_M˜S_S_′Π˜S_′, (21) where Π˜ is an arbitrary vector with ∑S_Π˜S_2=1, and it must be perpendicular to the equilibrium vector (P˜S_eq≡PS_eq), because λmin is the second smallest eigenvalue of M˜. Let Π˜=12(δS_S_g−δS_−S_g). Then according to (21) λmin≤−12(M˜Sg,Sg−M˜Sg,−Sg−M˜−Sg,Sg+M˜−Sg,−Sg)=−M˜Sg,Sg=−MSg,Sg=∑SWS,Sg=∑Sd(S,Sg)=1γ(ωSg,S;β), (22) where d(S, Sg) is the Hamming distance, and the S ↦ −S symmetry was used. In the bosonic bath λmin(β)≤∑S_d(S_,S_g)=1η(ΔES_)α1eβΔES_−1, (23) where ΔES_≔ES_−ES_g>0. At low temperatures this is the sum of some e−βΔES_ functions, so λmin(β) can be estimated from above with an exponential function.

Fig 2 shows λmin(β), λmax(β) and −MS_g,S_g (the upper bound) for a 4 × 4, ferromagnetic, 2D Ising model with Ohmic bath. The dashed vertical line marks the critical temperature (βcJ=ln(1+2)2≈0.44). The left figure is in log-log scale, where we can see, that at low temperature the eigenvalues has a β^−1^ temperature dependence, and λmax(β) converges, and the right figure with lin-log scale shows, that λmin(β) goes to zero exponentially.

4 Eigenvalues of the uniform Ising model

The M matrices are 2^N^ × 2^N^ large; therefore, we cannot see how the eigenvalues behave at the thermodynamic limit. However, the uniform, fully connected Ising model is so symmetric that an effective equation can be derived with the same relaxation time as the original equation.

The energy of the model is ES_=−JN∑i,j=1(i>j)NSiSj. (24) The 1/N factor is to keep the energy extensive and J > 0. Given an S microstate, it consists of N↑ spins with Si = 1 and N↓ spins with Si = −1. The number of spins is constant, i.e., N↑ + N↓ = N = fix. The energy of such a configuration is ES_=−JN[N↑(N↑−1)2+N↓(N↓−1)2−N↑N↓], (25) If N is fixed, then the energy is the function of only N↑. The symmetry of the system is that we can permute the spins in any way, the energy and the M matrix remains the same. If in the dynamics the initial condition also has this symmetry, then the PS will inherit this property. The slowest mode propagates between the two deepest valleys of the energy landscape, which are the ↑ ↑ …↑ and ↓↓ …↓. Assume that initially P↓↓…↓(t = 0) = 1, and we want to determine relaxation time, where P↓↓…↓(tr) ≈ P↑ ↑…↑(tr). Since both the equations and the initial condition have the permutation symmetry all the probabilities, which have the same up spin have the same value, e.g. for 3 spins P↑↓↓(t) = P↓↑↓(t) = P↓↓↑(t) ∀t. The probability can only flow between spin configurations if the Hamming distance between them is 1. Let us introduce the following PN↑=∑S_′PS_=(NN↑)P↑…↑︸N↑↓…↓︸N−N↑, (26) where the prime denotes that only such configurations count where there are N↑ up spin. We can give a closed set of differential equations which only contain these new PN↑ probabilities. P˙N↑=(N−N↑+1)WN↑,N↑−1PN↑−1+(N↑+1)WN↑,N↑+1PN↑+1−(N↑WN↑−1,N↑+(N−N↑)WN↑+1,N↑)PN↑, (27) where WN↑,N↑+1=γ(EN↑+1−EN↑)=γ(−2N(2N↑−N+1)). (28) This master equation has only N + 1 variables instead of 2^N^, thus is easy to simulate for large systems. A comparison between the quantum and the thermal simulated annealing of the fully connected Ising model was investigated by Wauters et al. using a similarly reduced master equation [42]. Eq (27) has the form P˙N↑=∑N↑′=0NMN↑,N↑′redPN↑′, (29) and we want to determine the lowest (nonzero) eigenvalue of M^red^, which is the same as the lowest (nonzero) eigenvalue of M. The matrix M^red^ is sparse, because it is a tridiagonal matrix, i.e., only the main diagonal, the first diagonal below and above the main diagonal are nonzero. Fig 3a shows the temperature dependence of λmin for different system sizes. As N increases we can see, that around the critical temperature (which is βcJ = 1) the behaviour of the system changes. At Fig 3b we can see better that above the critical temperature (T > Tc) for large N values λmin converges, meaning for every system size there is a finite relaxation time. At the critical temperature (T = Tc), it follows a power law (λmin ∝ N^−0.5^). Below the critical temperature (T < Tc) λmin goes to 0 for large N, but doesn’t follow a power law. This behaviour is the famous critical slowing down phenomenon.

From the N → ∞ thermodynamic limit we can determine the dynamical critical exponent. Fig 4 shows λmin(T, N → ∞) as a function of the reduced temperature (T−TcTc). This follows an easy power law, because λmin ∝ T − Tc. In the next section we will see that this result can be obtained from the mean field approximation.

5 Time dependent mean field equation

Since the primary interest is the magnetization (mi ≔ 〈Si〉), we would like to derive a differential equation for it. Penrose showed that with some approximations, this is possible [43], and as the number of neighbors increases, this becomes more and more accurate.

Using the definition of mi and the master equation we get m˙i=∑S_P˙S_Si=∑S_S_′WS_S_′PS_′Si−∑S_S_′WS_′S_PS_Si′=∑S_′S_WS_′S_PS_(Si′−Si) (30) The WS′S matrix component is nonzero if the Hamming distance between S′ and S is one. Introducing Λi(S_,n)={Si|i≠n−Si|i=n (31) we can rewrite the double sum in (30). m˙i=∑S_∑n=1NWΛ_(S_,n),S_PS_(Λi(S_,n)−Si)=∑S_WΛ_(S_,i),S_PS_(−2Si)=−2⟨WΛ_(S_,i),S_Si⟩ (32) In the second step the (Λ~i(S, n) − Si) = 2Si~ δin identity was used. The nonzero elements of W are a function of the energy WΛ_(S_,i),S_=γ(ES_−EΛ_(S_,i))≡γ(−2h˜iSi), (33) where h˜i=∑jJijSj+hi, so this is still the function of the S random variable, but because Jii = 0 it is not a function of Si. Since Si can be only 1 or −1 the γ(−2h˜iSi) as a function of Si must have the γ(−2h˜iSi)≡γ(−2h˜i)+γ(2h˜i)2+γ(−2h˜i)−γ(2h˜i)2Si (34) form. Using (7) yields γ(−2h˜iSi)=γ(2h˜i)[e−2βh˜i+12+e−2βh˜i−12Si]=γ(2h˜i)e−2βh˜i+12[1−tanh(βh˜i)Si], (35) then substituting back to (32) gives m˙i=−〈γ(2h˜i)(e−2βh˜i+1)(Si−tanh(βh˜i))〉. (36) Eq (36) is similar to the Callen equation 44, 45, where the averaging is outside the hyperbolic function. In order to get a closed equation for the expected values, the average must move inside, and instead of the Si random variables, their mi expected values must be written. m˙i=−γ(2(ΣjJijmj+hi))(1+e−2β(∑jJijmj+hi))×(mi−tanh(β(ΣjJijmj+hi))) (37) This mean field approximation is valid as long as the fluctuations are small, which holds if every spin interacts with many other spins. This can be due to long range interaction or in high spatial dimensions [46]. The right-hand side of Eq (37) contains the self-consistent equation from the equilibrium statistical physics; hence if the equation of state is satisfied, then m˙i=0.

Eq (37) contains both the real-time and the temperature of the bath. The temperature can also be time-dependent, and in that case, we could get thermal annealing, but if the temperature is constant, we can determine the relaxation time and the dynamical critical exponent. If m(t) = m^eq^ + δm(t), where m^eq^ is the equilibrium solution and δm(t) is small, then the linearized equation of (37) is δm˙i=−bi(m_eq)∑j{(δij−βJijcosh2(β∑kJikmkeq+hi))δmj}, (38) where bi(m_eq)=γ(2(ΣjJijmjeq+hi))(1+e−2β(∑jJijmjeq+hi)). (39) Using the 1cosh2(x)≡1−tanh2(x) identity, and the equation of state we get to δm˙i=−bi(m_eq)∑j{(δij−βJij(1−(mjeq)2))δmj}. (40) Eq (40) contains the inverse susceptibility of the mean field Ising model. χij−1≔∂2FMFA∂mi∂mj=−Jij+Tδij1−mi2, (41) where FMFA(m_,h_,T)=−12∑ijJijmimj−∑ihimj+T∑i[1+mi2ln(1+mi2)+1−mi2ln(1−mi2)]. (42) Substituting the inverse susceptibility back into (40) yields δm˙i=−bi(m_eq)β(1−(mieq)2)∑jχij−1δmj≡−Γi∑jχij−1δmj. (43) This is a well-known equation in the theory of dynamical critical phenomena [47], but it is usually derived from the m_˙=−Γ∂m_FMFA phenomenological equation. Now we can see, how it is related to a master equation and the spin-boson model. If the system is symmetric in a sense, that all the spins behave the same, then (43) simplifies to δm˙=−Γχ−1δm, (44) where χ−1=∑jχij−1.

6 Time dependent mean field equation for the uniform Ising model

As before in section 4 the uniform Ising model will be studied because in the equilibrium case in the thermodynamic limit, it gives back the exact results. According to (42) the mean field free energy is FMFA(m,h,T)N=−12Jm2−hm+T(1+m2ln(1+m2)+1−m2ln(1−m2)), (45) and the time dependent mean field equation is m˙=−γ(2(Jm+h))(1+e−2β(Jm+h))×(m−tanh(β(Jm+h))). (46) If h = 0 the critical temperature is Tc = J, and above this temperature the equilibrium solution is m^eq^ = 0. The inverse susceptibility is χ−1=−J+T≡T−Tc, (47) therefore λmin=Γ(T−Tc), (48) where Γ = 2γohm(0; β)β = 2η in the Ohmic bath. In the super-Ohmic bath, this is zero because γsup(0) = 0, i.e., there is no transition between states with the same energy. In a real physical system it means that besides the Glauber dynamics some other effects are not negligible, e.g. 2 spin flips. Eq (48) is the same result that we have already seen in Fig 4. As in equilibrium statistical physics, the mean field approximation gives back the exact result for the uniform model in the thermodynamic limit.

At the critical temperature the inverse susceptibility is zero, the linear term vanishes, and we need the higher order terms. Taylor expanding (46) at T = Tc around m = m^eq^ ≡ 0 up to the third order gives. δm˙=−23ηJδm3 (49) which has the δm(t)∝t−12 (50) solution for large ts, which means there isn’t a characteristic time.

Below the critical temperature, the linearized equation is good again; only m^eq^ changes. On the other hand Eq (46) can give a different solution from the (29) master equation. If we want to compare these two equations, the initial condition must also be the same, which gives a restriction to the initial condition of (29). In the mean field approximation, the probability is a product of the one-particle PS_MFA=∏i1+miSi2, (51) which means the initial probability of (29) is PN↑(t=0;m)=(NN↑)(1+m2)N↑(1−m2)N↓ (52) If h = 0, then the master equation must converge to the m = 0 solution, but the time-dependent mean field equation finds one of the minima of the free energy if initially m ≠ 0. If h is finite and m is in the valley of the global minimum (point A in Fig 5), then in the N → ∞ limit the solution of the master equation converges to the mean field solution (Fig 6).

On the other hand, if initially, m is in the valley of the local minimum (point B in Fig 5), then the solution of the mean field equation converges to this local minimum, but the solutions of the master equation are different. The probabilities converge to the equilibrium PN↑eq(h,N)∝(NN↑)e−βEN↑(h,N) (53) distribution, so they tend to approach m^eq^ for large N values in the t → ∞ limit, but as N grows, so does the relaxation time. In the thermodynamic limit the relaxation time diverges as in Fig 3. In the N → ∞ limit the solution of the master equation converges to the solution of the mean field equation and none of them will approach the global minimum, because the relaxation time will be infinite.

7 Conclusion

In this work, we have studied a Glauber-type master equation derived from the spin-boson model. The most relevant dynamical properties are encoded in the eigenvalues of the transition matrix of the master equation. They are temperature dependent and behave significantly differently below, above, and at the critical temperature as a function of the system size. In the case of the uniform, fully connected Ising model, in the thermodynamic limit, above the critical temperature the relaxation time follows a power tr ∝ (T − Tc)^−1^.

Using the time-dependent mean field equation, we could also investigate the thermodynamic limit. Its dynamics differ significantly from the finite size master equation if the initial condition is close to a local minimum of the free energy, which means that the relaxation time diverges. If it is close to the global minimum, then the solution of the finite size master equation converges to the solution of the time-dependent mean field equation.