Dynamics of open quantum spin systems: An assessment of the quantum master equation approach
aa r X i v : . [ phy s i c s . c o m p - ph ] A ug Dynamics of open quantum spin systems: An assessment of the quantum master equation approach
P. Zhao and H. De Raedt ∗ Zernike Institute for Advanced Materials,University of Groningen, Nijenborgh 4,NL-9747AG Groningen, The Netherlands
S. Miyashita
Department of Physics, Graduate School of Science,University of Tokyo, Bunkyo-ku, Tokyo 113-0033, Japan
F. Jin
Institute for Advanced Simulation, J¨ulich Supercomputing Centre,Forschungszentrum J¨ulich, D-52425 J¨ulich, Germany
K. Michielsen
Institute for Advanced Simulation, J¨ulich Supercomputing Centre,Forschungszentrum J¨ulich, D-52425 J¨ulich, Germany andRWTH Aachen University, D-52056 Aachen, Germany (Dated: August 16, 2016)Data of the numerical solution of the time-dependent Schr¨odinger equation of a system containing one spin-1/2 particle interacting with a bath of up to 32 spin-1/2 particles is used to construct a Markovian quantum masterequation describing the dynamics of the system spin. The procedure of obtaining this quantum master equation,which takes the form of a Bloch equation with time-independent coefficients, accounts for all non-Markovianeffects in as much the general structure of the quantum master equation allows. Our simulation results showthat, with a few rather exotic exceptions, the Bloch-type equation with time-independent coefficients providesa simple and accurate description of the dynamics of a spin-1/2 particle in contact with a thermal bath. Acalculation of the coefficients that appear in the Redfield master equation in the Markovian limit shows thatthis perturbatively derived equation quantitatively differs from the numerically estimated Markovian masterequation, the results of which agree very well with the solution of the time-dependent Schr¨odinger equation.
PACS numbers: 03.65.-w, 05.30.-d, 03.65.YzKeywords: quantum theory, quantum statistical mechanics, open systems, quantum master equation ∗ [email protected]; Corresponding author I. INTRODUCTION
In general, a physical system can seldom be considered as completely isolated from its environment. Such closed systems canand should, of course, be studied in great detail. However, as they lack the ability to interact with the environment in which theyare embedded or with the apparatus that is used to perform measurements on it, such studies do not include the effects of the,usually uncontrollable, environment which may affect the dynamics of the system in a non-trivial manner. The alternative is toconsider the system of interest as an open system that is a system interacting with its environment.The central idea of theoretical treatments of open quantum systems is to derive approximate equations of motion of the systemby elimination of the environmental degrees of freedom [1–4]. In 1928, Pauli derived a master equation for the occupationprobabilities of a quantum subsystem interacting with the environment [5]. Since then, various methods have been developedto derive quantum master equations starting from the Liouville-von Neumann equation for the density matrix of the wholesystem [1–4, 6]. In order to obtain an equation of motion for the system which is tractable and readily amenable to detailedanalysis, it is customary to make the so-called Markov approximation, which in essence assumes that the correlations of the bathdegrees of freedom vanish on a short time span.Without reference to any particular model system, in 1970, Lindblad derived a quantum master equation which is Markovianand which preserves positivity (a non-negative definite density matrix) during the time evolution [4, 7]. The applicability of theLindblad master equation is restricted to baths for which the time correlation functions of the operators that couple the systemto the bath are essentially d -functions [8], an assumption that may be well justified in quantum optics [6].Using second-order perturbation theory, Redfield derived a master equation which does not require the bath correlations tobe approximately d -functions in time [1]. The Redfield master equation has found many applications to problems where thedynamics of the bath is faster than that of the system, for instance to the case of nuclear magnetic resonance in which the systemconsists of one spin coupled to other spins and/or to phonons. This approach and variations of it have been successfully appliedto study the natural linewidth of a two-level system [9–11], systems of interacting spins [12] and nonlinear spin relaxation [13].The Redfield master equation can be systematically derived from the principles of quantum theory but only holds for weakcoupling. However, the Redfield master equation may lead to density matrices that are not always positive, in particular whenthe initial conditions are such that they correspond to density matrices that close to the boundary of physically admissible densitymatrices [14, 15].Obviously, the effect of the finite correlation time of the thermal bath becomes important when the time scale of the systemis comparable to that of the thermal bath. Then the Markovian approximation may no longer be adequate and in derivingthe quantum master equation, it becomes necessary to consider the non-Markovian aspects and to treat the initial conditioncorrectly [4, 16–24].By introducing the concept of slippage in the initial conditions, it was shown that the Markovian equations of motion obtainedin the weak coupling regime are a consistent approximation to the actual reduced dynamics and that slippage captures theeffects of the non-Markovian evolution that takes place in a short transient time, of the order of the relaxation time of theisolated bath [14]. Provided that nonlocal memory effects that take place on a very short time scale are included, the Markovianapproximation that preserves the symmetry of the Hamiltonian yields an accurate description of the system dynamics [14].Following up on this idea, a general form of a slippage operator to be applied to the initial conditions of the Redfield masterequation was derived [8]. The slippage was expressed in terms of an operator describing the non-Markovian dynamics of thesystem during the time in which the bath relaxes on its own, relatively short, time scale. It was shown that the application of theslippage superoperator to the initial density matrix of the system yields a Redfield master equation that preserves positivity [8].Apparently, the difference between the non-Markovian dynamics and its Markovian approximation can be reduced significantlyby first applying the slippage operator and then letting the system evolve according to the Redfield master equation [8].The work discussed and cited earlier almost exclusively focuses on models of the environment that are described by a collec-tion of harmonic oscillators. In contrast, the focus of this paper is on the description of the time evolution of a quantum systemwith one spin-1/2 degree of freedom coupled to a larger system of similar degrees of freedom, acting as a thermal bath. Ourreasons for focusing on spin-1/2 models are twofold.First, such system-bath models are relevant for the description of relaxation processes in nuclear magnetic and electron spinresonance [1, 10, 25] but have also applications to, e.g. the field of quantum information processing, as most of the models usedin this field are formulated in terms of qubits (spin-1/2 objects) [26, 27].Secondly, the aim of the present work is to present a quantitative assessment of the quantum master equation approach bycomparing the results with those obtained by an approximation-free, numerical solution of the time-dependent Schr¨odingerequation of the system+bath. The work presented in this paper differs from earlier numerical work on dissipative quantumdynamics [28–32] by accounting for the non-trivial many-body dynamics of the bath without resorting to approximations, atthe expense of using much more computational resources. Indeed, with state-of-the-art computer hardware, e.g the IBM Blue-Gene/Q, and corresponding simulation software [33], it has become routine to solve the time-dependent Schr¨odinger equationfor systems containing up to 36 spin-1/2 objects. As we demonstrate in this paper, this allows us to mimic a large thermal bathat a specific temperature and solve for the full dynamic evolution of a spin-1/2 object coupled to the thermal bath of spin-1/2objects.From the numerically exact solution of the Schr¨odinger dynamics we compute the time-evolution of the density matrix ofthe system and by least-square fitting, obtain the “optimal” quantum master equation that approximately describes the sametime-evolution. For a system of one spin-1/2 object, this quantum master equation takes the form of a Bloch equation with time-independent coefficients. Clearly, this procedure of obtaining the quantum master equation is free of any approximation andaccounts for all non-Markovian effects in as much the general structure of the quantum master equation allows. Our simulationresults show that, with a few rather exotic exceptions, the Bloch-type equation with time-independent coefficients provides avery simple and accurate description of the dynamics of a spin-1/2 object in contact with a thermal bath.The paper is organized as follows. In section II, we give the Hamiltonians that specify the system, bath and system-bathinteraction. Section III briefly reviews the numerical techniques that we use to solve the time-dependent Schr¨odinger equation, tocompute the density matrix, and to prepare the bath in the thermal state at a given temperature. We also present simulation resultsthat demonstrate that the method of preparation yields the correct thermal averages. For completeness, Sec. IV recapitulatesthe standard derivation of the quantum master equation, writes the formal solution of the latter in a form that is suited for ournumerical work and shows that the Redfield equations have this form. We then use the simulation tool to compute the correlationsof the bath-operators that determine the system-bath interaction and discuss their relaxation behavior. Section V explains theleast-square procedure of extracting, from the solution of the time-dependent Schr¨odinger equation, the time-evolution matrixand the time-independent contribution that determine the “optimal” quantum master equation. This least-square procedure isvalidated by its application to data that originate from the Bloch equation, as explained in Appendix A. In Sec. VI, we specifythe procedure by which we fit the quantum master equation to the data obtained by solving the time-dependent Schr¨odingerequation and present results of several tests. The results of applying the fitting procedure to baths containing up to 32 spins arepresented in Sec. VII. Finally, in Sec. VIII, we discuss some exceptional cases for which the quantum master equation is notexpected to provide a good description. The paper concludes with the summary, given in Sec. IX. II. SYSTEM COUPLED TO A BATH: MODEL
The Hamiltonian of the system (S) + bath (B) takes the generic form H = H S + H B + l H SB . (1)The overall strength of the system-bath interaction is controlled by the parameter l . In this work, we limit ourselves to a systemwhich consists of one spin-1/2 object described by the Hamiltonian H S = − h x s x , (2)where sss n = ( s xn , s yn , s zn ) = ( s n , s n , s n ) denote the Pauli-spin matrices for spin-1/2 object n , and h x is a time-independentexternal field. Throughout this paper, we adopt units such that ¯ h = h x = / p / h x . We willuse the double notation with the ( x , y , z ) and ( , , ) superscripts because depending on the situation, it simplifies the writingconsiderably.The Hamiltonian for the system-bath interaction is chosen to be H SB = − N B (cid:229) n = (cid:0) J xn s xn s x + J yn s yn s y + J zn s zn s z (cid:1) = (cid:229) a = x , y , z s a B a = (cid:229) i = s i B i , (3)where N B is the number of spins in the bath, the J a n are real-valued random numbers in the range [ − J , + J ] and B x = B = − N B (cid:229) n = J xn s xn , B y = B = − N B (cid:229) n = J yn s yn , B z = B = − N B (cid:229) n = J zn s zn (4)are the bath operators which, together with the parameter l , define the system-bath interaction. As the system-bath interactionstrength is controlled by l , we may set J = / H B we take H B = − K N B (cid:229) n = (cid:0) s xn s xn + + s yn s yn + + D s zn s zn + (cid:1) − N B (cid:229) n = ( h xn s xn + h zn s zn ) . (5)The fields h xn and h zn are real-valued random numbers in the range [ − h x B , + h x B ] and [ − h z B , + h z B ] , respectively. In our simulationwork, we use periodic boundary conditions s a n = s a n + N B for a = x , y , z . Note that we could have opted equally well to useopen-end boundary conditions but for the sake of simplicity of presentation, we choose the periodic boundary conditions. For D =
1, the first term in Eq. (5) is the Hamiltonian of the one-dimensional (1D) Heisenberg model on a ring.As a second choice, we consider the 1D ring with Hamiltonian H B = − N B (cid:229) n = (cid:0) K xn s xn s xn + + K yn s yn s yn + + K zn s zn s zn + (cid:1) − N B (cid:229) n = ( h xn s xn + h zn s zn ) , (6)where the K xn ’s, K yn ’s, and K zn ’s are uniform random numbers in the range [ − K , K ] . Because of the random couplings, it is unlikelythat it is integrable (in the Bethe-ansatz sense) or has any other special features such as conserved magnetization etc.The bath Hamiltonians (5) and (6) both share the property that the distribution of nearest-neighbor energy levels is of Wigner-Dyson-type, suggesting that the correspondig classical baths exhibit chaos. Earlier work along the lines presented in this paperhas shown that spin baths with a Wigner-Dyson-type distribution are more effective as sources for fast decoherence than spinbaths with Poisson-type distribution [34]. Fast decoherence is a prerequisite for a system to exhibit fast relaxation to the thermalequilibrium state [35, 36]. Extensive simulation work on spin-baths with very different degrees of connectivity [37–40] suggestthat as long as there is randomness in the system-bath coupling and randomness in the intra-bath couplings, the simple models(5) and (6) may be considered as generic spin baths.Finally, as a third choice, we consider H B = − (cid:229) h n , n ′ i (cid:16) K xn , n ′ s xn s xn ′ + K yn , n ′ s yn s yn ′ + K zn , n ′ s zn s zn ′ (cid:17) − N B (cid:229) n = ( h xn s xn + h zn s zn ) , (7)where the K xn , n ′ ’s, K yn , n ′ ’s, and K zn , n ′ ’s are uniform random numbers in the range [ − K , K ] , and (cid:229) h n , n ′ i denotes the sum over all pairsof nearest neighbors on a three-dimensional (3D) cubic lattice. Again, because the random couplings and the 3D connectivity,it is unlikely that it is integrable or has any other special features such as conserved magnetization etc. As the solution of thetime-dependent Schr¨odinger equation (TDSE) for the 3D model Eqs. (7) takes about a factor of 2 more CPU time than in thecase of a 1D model with the same number of bath spins, in most of our simulations we will use the 1D models and only use the3D model to illustrate that the connectivity of the bath is not a relevant factor. III. QUANTUM DYNAMICS OF THE WHOLE SYSTEM
The time evolution of a closed quantum system defined by Hamiltonian (1) is governed by the TDSE i ¶¶ t | Y ( t ) i = H | Y ( t ) i . (8)The pure state | Y ( t ) i of the whole system S + B evolves in time according to | Y ( t ) i = e − itH | Y ( ) i = D S (cid:229) i = D B (cid:229) p = c ( i , p , t ) | i , p i , (9)where D S = D B = N B are the dimensions of the Hilbert space of the system and bath, respectively. The coefficients { c ( i , p , t ) } are the complex-valued amplitudes of the corresponding elements of the set {| i , p i} which denotes the complete setof the orthonormal states in up–down basis of the system and bath spins.The size of the quantum systems that can be simulated, that is the size for which Eq. (9) can actually be computed, is primarilylimited by the memory required to store the pure state. Solving the TDSE requires storage of all the numbers { c ( i , p , t ) | i = , , p = , . . . , N B } . Hence the amount of memory that is required is proportional to 2 N B + , that is it increases exponentiallywith the number of spins of the bath. As the number of arithmetic operations also increases exponentially, it is advisable to use13 - 15 digit floating-point arithmetic (corresponding to 16 = bytes for each pair of real numbers). Therefore, representing apure state of N B + / N B + bytes. For example, for N B =
23 ( N B = | Y ( t ) i . In practice we need storage for three vectors, andmemory for communication buffers, local variables and the code itself.The CPU time required to advance the pure state by one time step t is primarily determined by the number of operations tobe performed on the state vector, that is it also increases exponentially with the number of spins. The elementary operationsperformed by the computational kernel can symbolically be written as | Y i ← U | Y i where the U ’s are sparse unitary matriceswith a relatively complicated structure. A characteristic feature of the problem at hand is that for most of the U ’s, all elementsof the set { c ( i , p , t ) | i = , , p = , N B } are involved in the operation. This translates into a complicated scheme for efficientlyaccessing memory, which in turn requires a sophisticated communication scheme [33].We can exclude that the conclusions that we draw from the numerical results are affected by the algorithm used to solvethe TDSE by performing the real-time propagation by e − itH by means of the Chebyshev polynomial algorithm [41–44]. Thisalgorithm is known to yield results that are very accurate (close to machine precision), independent of the time step used [45].A disadvantage of this algorithm is that, especially when the number of spins exceeds 28, it consumes significantly more CPUand memory resources than a Suzuki-Trotter product-formula based algorithm [45]. Hence, once it has been verified that thenumerical results of the latter are, for practical purposes, as good as the numerically exact results, we use the latter for thesimulations of the large systems. A. Density matrix
According to quantum theory, observables are represented by Hermitian matrices and the correspondence with measurablequantities is through their averages defined as [46] h A i = Tr r ( t ) A , (10)where A denotes a Hermitian matrix representing the observable, r ( t ) is the density matrix of the whole system S + B at time t and Tr denotes the trace over all states of the whole system S + B. If the numerical solution of the TDSE for a pure state of N B + | F i is a pure state, picked randomly from the2 N B + -dimensional unit hypersphere, one can show in general that for Hermitian matrices A [47] Tr A = D h F | A | F i ± O ( D − / ) , (11)where D is the number of diagonal elements of the matrix A (= the dimension of the Hilbert space) and ± O ( x ) should be readas saying that the standard deviation is of order x . For the case at hand D = N B + , hence Eq. (11) indicates that for a large bath,the statistical errors resulting from approximating Tr A by h F | A | F i vanishes exponentially with the number of bath spins. Forlarge baths, this property makes the problem amenable to numerical simulation. Therefore, from now on, we replace the “ Tr ”by a matrix element of a random pure state whenever the trace operation involves a number of states that increases exponentiallywith the number of spins (in the present case, bath spins only).The state of the system S is completely described by the reduced density matrix r S ( t ) ≡ Tr B r ( t ) , (12)where r ( t ) is the density matrix of the whole system S + B at time t , Tr B denotes the trace over the degrees of freedom of thebath, and Tr S r S ( t ) = Tr r ( t ) =
1. In practice, as the dimension of the Hilbert space of the bath may be assumed to be large, wecan, using the “random-state technology”, compute the trace over the bath degrees of freedom as ( Tr B A ) i , j ≈ D B (cid:229) p = c ∗ ( i , p , t ) c ( j , p , t ) h i , p | A | j , p i . (13)In the case that the system contains only one spin, which is the case that we consider in the present work, the reduced densitymatrix can, without loss of generality, be written as r S ( t ) = (cid:229) a = x , y , z [ + r a ( t ) s a ] = (cid:229) k = h + r k ( t ) s k i , (14)where r x ( t ) = r ( t ) , r y ( t ) = r ( t ) and r z ( t ) = r ( t ) are real numbers. Making use of the “random-state technology”, it followsimmediately from Eq. (14) that r ( t ) = r x ( t ) = Tr S r S ( t ) s x = Tr r ( t ) s x ≈ h Y ( t ) | s x | Y ( t ) i r ( t ) = r y ( t ) = Tr S r S ( t ) s y = Tr r ( t ) s y ≈ h Y ( t ) | s y | Y ( t ) i r ( t ) = r z ( t ) = Tr S r S ( t ) s z = Tr r ( t ) s z ≈ h Y ( t ) | s z | Y ( t ) i . (15)Therefore, to obtain (accurate approximations to) the expectation values of the system operators we compute the expressionsthat appear in the left-hand side of Eq. (15) using the numerical solution of the TDSE in the form given by Eq. (9). B. Thermal equilibrium state
As a first check on the numerical method, it is of interest to simulate the case in which the system+bath are initially in thermalequilibrium and study the effects of the bath size N B and system-bath interaction strength l on the expectation values of thesystem spin.The procedure is as follows. First we generate a random state of the whole system, meaning that | F ( b ) i = e − b H / | F ih F | e − b H | F i / , (16)where b denotes the inverse temperature. As one can show that for any observable A ( t ) [47] h A ( t ) i = Tr e − b H A ( t ) Tr e − b H = h F ( b ) | A ( t ) | F ( b ) i ± O ( D − / ) , (17)we can use h F ( b ) | A | F ( b ) i to estimate h A ( t ) i . As e − b H commutes with e − itH , h A ( t ) i = h A ( t = ) i is time independent.Excluding the trivial case that [ H , A ( t )] = h F ( b ) | A ( t ) | F ( b ) i = h F ( b ) | e + itH A e − itH | F ( b ) i depends on time: indeed, ingeneral the random state | F ( b ) i is unlikely to be an eigenstate of H . Therefore, the simulation data obtained by solving theTDSE with | F ( b ) i as the initial state should display some time dependence. However, from Eq. (17), it follows directly that thetime dependent contributions will vanish very fast, namely as D − / . Hence this time dependence, an artifact of using “randomstate technology”, reveals itself as statistical fluctuations and can be ignored.For the system in thermal equilibrium at the inverse temperature b we have h s x i = tanh ( b h x ) , h s y i = , h s z i = . (18)In Fig. 1 we show simulation results for a bath at b = N B =
13 (left) and N B =
28 spins (right). If the system-bath interactionis sufficiently weak then, from Eq. (18), we expect that h s x i ≈ tanh b h x which for b h x = h s x i ≈ . N B =
13, it is clear that the spin averages fluctuate (due to the use of the random thermal state which is not aneigenstate of H ). As expected, for N B =
28 the fluctuations are much smaller, in concert with Eq. (17).Computing the time averages for a bath with N B =
13 and for the time interval [ , T ] with T = T Z T dt h F ( b ) | s x ( t ) | F ( b ) i = . ( . ) T Z T dt h F ( b ) | s y ( t ) | F ( b ) i = . ( . ) T Z T dt h F ( b ) | s z ( t ) | F ( b ) i = − . ( . ) , (19)where the numbers in parenthesis give the standard deviation. For N B =
28 and for the time interval [ , T ] with T =
200 we find1 T Z T dt h F ( b ) | s x ( t ) | F ( b ) i = . ( . ) T Z T dt h F ( b ) | s y ( t ) | F ( b ) i = . ( . ) T Z T dt h F ( b ) | s z ( t ) | F ( b ) i = . ( . ) , (20)indicating that for most practical purposes, a bath of N B =
28 spin may be sufficiently large to mimic an infinitely large bath.The numbers in Eq. (20) also give an indication of the statistical fluctuations that we may expect for a bath containing N B = l chosen, the second-order corrections in l are of the order of 0 .
01 and arehidden in the statistical fluctuations, suggesting that values of l ≤ . l , we have h s x i = h s x i S − bl ( h s x i S − ) h B x i B , (21)where h . i S and h . i B denote the thermal equilibrium averages with respect to the system and bath, respectively. For the sake ofargument, consider the case that K = h zn = h xn = h x B for all n = , . . . , N B (the same reasoning applies to the contributionsof second order in l ). Then, Eq. (21) becomes h s x i = tanh ( b h x ) + bl N B ( − tanh ( b h x )) tanh ( b h x B ) , (22)showing that the contribution of the “perturbation term” increases with the number of spins in the bath. In other words, it is notsufficient to consider small values of l . For the perturbation by the bath to be weak, it is necessary that l N B is small. In thisrespect the spin bath considered in this paper is not different from e.g. the standard spin-boson model [4]. In our simulationwork, we adopt a pragmatic approach: we simply compute the averages and compare them with the theoretical results of theisolated system (as we did above). The coupling l is considered to be small enough if the corrections are hidden in the statisticalfluctuations. -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s ( t ) > -3 t(a) < s (t)>< s (t)>< s (t)> -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s ( t ) > -3 t(b) < s (t)>< s (t)>< s (t)> FIG. 1. (color online) Time evolution of the average of the system spin as obtained by solving the TDSE with a random thermal state at b = K = − / D = J = / h x B = h z B = /
8. The system-bath interaction l = .
1. (a) N B =
13; (b) N B =
28. Lines connecting the data points are guide to the eye.
IV. QUANTUM MASTER EQUATION: GENERALITIES
We are interested in the dynamics of a system, the degrees of freedom of which interact with other degrees of freedom of a“bath”, “environment”, etc. The combination of system + bath forms a closed quantum system. When we consider the systemonly, we say that we are dealing with an open quantum system. The quantum state of the system + bath is represented by thedensity matrix r = r ( t ) which evolves in time according to ¶r ( t ) ¶ t = i [ r ( t ) , H ] , (23)where H is the Hamiltonian of the system + bath (recall that we adopt units such that ¯ h = P be the projector onto the “relevant” part and introduce the Liouville operator L A = i [ A , H ] . Denoting by Q = − P the projector on the “uninteresting” part, it follows that ¶ P r ( t ) ¶ t = PL P r ( t ) + PL Q r ( t ) , (24) ¶ Q r ( t ) ¶ t = QL P r ( t ) + QL Q r ( t ) . (25)Note that because H is Hermitian, i L , i PL P and i QL Q are Hermitian too. The formal solution of the matrix-valued,inhomogeneous, linear, first-order differential equation Eq. (25) reads as Q r ( t ) = e t QL Q Q r ( t = ) + Z t du e u QL Q
QL P r ( t − u ) , (26)as can be verified most easily by calculating its derivative with respect to time and using PP = P , PQ = QP = QQ = Q . Substituting Eq. (26) into Eq. (24) yields ¶ P r ( t ) ¶ t = PL P r ( t ) + PL Q e t QL Q Q r ( t = ) + Z t du PL Q e u QL Q
QL P r ( t − u ) . (27)We are primarily interested in the time evolution of the system. Therefore, we choose P such that it projects onto the systemvariables and we perform the trace over the bath degrees-of-freedom. A common choice for the projector P is [4, 8, 14, 20, 22,23] P A = r B Tr B A , (28)where r B = e − b H B Tr B e − b H B , (29)is the density matrix of the bath in thermal equilibrium. Accordingly, the density matrix of the system is given by r S ( t ) = Tr B P r ( t ) = Tr B r ( t ) , (30)consistent with Eq. (12).In the present work, we will mostly consider initial states that are represented by the direct-product ansatz r ( t = ) = r S r B , . (31)but occasionally, we also consider as an initial state, the thermal equilibrium state of the system + bath, that is r ( t = ) = e − b H / Tr e − b H , see Sec. III B. The direct-product ansatz Eq. (31) not only implies Q r ( t = ) = ¶r S ( t ) ¶ t = Tr B PL P r ( t ) + Z t du Tr B PL Q e u QL Q
QL P r ( t − u ) , (32)which is not a closed equation for r S ( t ) yet [20].Using the explicit form of the Hamiltonian Eq. (1), the first term in Eq. (32) may be written as Tr B PL P r ( t ) = L r S ( t ) where for any system operator X S , L X S ≡ − i ( [ H S , X S ( t )] + (cid:229) i = h B i i B (cid:2) s i , X S ( t ) (cid:3)) , (33)and h B i i B ≡ Tr B r B B i . Therefore, Eq. (32) may be written as ¶r S ( t ) ¶ t = L r S ( t ) + Z t du Tr B PL Q e ( t − u ) QL Q QL r B r S ( u ) . (34)Using representation Eq. (14), multiplying both sides of Eq. (34) by s j , performing the trace over the system degree offreedom, and denoting rrr ( t ) = ( r ( t ) , r ( t ) , r ( t )) , Eq. (34) can be written as ¶r rr ( t ) ¶ t = L rrr ( t ) + Z t du M ( t − u ) rrr ( u ) + Z t du K ( u ) , (35)where L jk = Tr S s j L s k M jk ( u ) = Tr s j PL Q e u QL Q QL r B s k K j ( u ) = Tr s j PL Q e u QL Q QL r B . (36)As we have only made formal manipulations, solving Eq. (35) of the system is just as difficult as solving Eq. (23) of the wholesystem. In other words, in order to make progress, it is necessary to make approximations. A common route to derive an equationwhich can actually be solved is to assume that l is sufficiently small such that perturbation theory may be used to approximatethe second term in Eq. (34) and that it is allowed to replace r S ( u ) in Eq. (34) by r S ( t ) [4].As the purpose of the present work is to scrutinize the approximations just mentioned by comparing the solution obtainedfrom the Markovian quantum master equation with the one obtained by solving the TDSE, we will not dwell on the justificationof these approximations and derivation of this equation itself, but merely state that the result of making these approximations isan equation that may be cast in the form ¶r rr ( t ) ¶ t = A rrr ( t ) + b . (37)In the following we will refer to Eq. (37) as “the” quantum master equation (QMEQ). In Sec. IV A we give a well-knownexample of a quantum master equation that is of the form Eq. (37).The formal solution of Eq. (37) reads as rrr ( t ) = e t A rrr ( ) + Z t e ( t − u ) A b du , (38)or, equivalently rrr ( t + t ) = e t A rrr ( t ) + Z t e ( t − u ) A b du = e t A rrr ( t ) + B , (39)where B = Z t e ( t − u ) A b du , (40)does not depend on time. Equation (39) directly connects to the numerical work because in practice, we solve the TDSE with afinite time step t .Generally speaking, as a result of the coupling to the bath, the system is expected to exhibit relaxation towards a stationarystate, meaning that rrr ( t ) ≈ rrr ( ¥ ) for t sufficiently large. If such a stationary state exists, it follows from Eq. (39) that rrr ( ¥ ) ≈ e t A rrr ( ¥ ) + B or that B ≈ ( − e t A ) rrr ( ¥ ) , yielding rrr ( t + t ) − rrr ( ¥ ) ≈ e t A ( rrr ( t ) − rrr ( ¥ )) . (41)Equation (41) suggests that the existence of a stationary state implies that there is no need to determine B . However, numericalexperiments with the Bloch equation model (see Appendix A) show that using Eq. (41), a least-square fit to solution of the Blochequation often fails to yield the correct e t A . Therefore, as explained in Sec. V, we will use Eq. (39) and determine both e t A and B by least-square fitting to TDSE or Bloch equation data.We can now formulate more precisely, the procedure to test whether or not a quantum master equation of the form Eq. (37)provides a good approximation to the data r k ( t ) = h s k ( t ) i obtained by solving the TDSE of the system interacting with the bathusing a time step t . To this end, we use the latter data to determine the matrix e t A and vector B such that, in a least square sense,the difference between the data obtained by solving Eq. (39) for a substantial interval of time and the corresponding TDSE datais as small as possible. If the values of rrr ( t ) computed according to Eq. (39) are in good agreement with the data r k ( t ) , onemight say that at least for the particular time interval studied, there exists a mapping of the Schr¨odinger dynamics of the systemonto the QMEQ Eq. (37). A. Markovian quantum master equation: Example
We consider the Redfield master equation [1] under the Markovian assumption [4, 8] d r S ( t ) dt = − i [ H S , r S ( t )] + l (cid:229) j = (cid:16) R j r S ( t ) s j + s j r S ( t ) R † j − s j R j r S ( t ) − r S ( t ) R † j s j (cid:17) , (42)where r S ( t ) is the density matrix of the system. The operators R j are given by [8] R j = (cid:229) k = Z ¥ dt C jk ( t ) e − itH s s k e + itH s , j = , , , (43)where C jk ( t ) = Tr B r B B j ( t ) B k ( ) are the correlations of the bath operators [8]. The specific form of C jk ( t ) is not of interest tous at this time (but also see Sec. VII). For what follows, it is important that the specific form Eq. (43) of the operators R j allowsus to write R j = (cid:229) k = r jk s k , (44)where r j = Z ¥ dt C j ( t ) r j = Z ¥ dt (cid:0) C j ( t ) cos 2 h x t + C j ( t ) sin 2 h x t (cid:1) r j = Z ¥ dt (cid:0) C j ( t ) cos 2 h x t − C j ( t ) sin 2 h x t (cid:1) , (45)0 | C ( .,. ) | t(a) |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | t(b) |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | t(c) |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | t(d) |C(1,1)||C(2,2)||C(3,3)| FIG. 2. (color online) The absolute values of three of the nine bath-operator correlations Eq. (47) as obtained by solving the TDSE for abath of N B =
32 spins with a random thermal state at b = H SB are J = / h x B = h z B = /
8. (a) the bath Hamiltonian H B is given by Eq. (5) with K = − / D = l = l = .
1; (c) the bath Hamiltonian H B is given by Eq. (6) with K = / l =
0; (d) same as (a) except that l = . do not depend on time (due to the Markov approximation).As a first step, we want to derive from Eq. (42), the corresponding equations in terms of the r k ( t ) ’s. This can be done by usingrepresentation Eq. (14), multiplying both sides of Eq. (42) with s k for k = , , r . We obtain d r dt = + l (cid:2)(cid:0) r I23 − r I32 (cid:1) − (cid:0) r R22 + r R33 (cid:1) r + r R21 r + r R31 r (cid:3) d r dt = + h x r + l (cid:2)(cid:0) r I31 − r I13 (cid:1) + r R12 r − (cid:0) r R11 + r R33 (cid:1) r + r R32 r (cid:3) d r dt = − h x r + l (cid:2) r I12 − r I21 + r R13 r + r R23 r − (cid:0) r R11 + r R22 (cid:1) r (cid:3) , (46)where we used the notation z = z R + iz I . It directly follows that Eq. (46) can be written in the form Eq. (37). It is straightforwardto show that this holds for quantum master equations of the Lindblad form as well. B. Bath correlations
A crucial assumption in deriving the QMEQ Eq. (37) from the exact equation Eq. (34) is that the correlations of the bathdecay on a short time scale, short relative to the time scale of the motion of the system spin [4]. Moreover, in the perturbative1derivation of quantum master equations, such as the Redfield master equation, it is assumed that the time evolution of the bathoperators is governed by the bath Hamiltonian only [4].Having the time evolution of the whole system at our disposal, we can compute, without additional assumptions or approxi-mations, the correlations C ( i , j , t ) = Tr r ( t = ) B i ( t ) B j ( ) , i , j = , , , (47)of the bath operators Eq. (4). Note that in general, Eq. (47) is complex-valued and that, because of the choice Eq. (31), C ( i , j , t ) = C i j ( t ) if l =
0. Of particular interest is the question whether, for the chosen value of the system-bath interaction l , the dynamicsof the system spin significantly affects the bath dynamics.In Fig. 2 we present simulation results of the correlations | C ( i , i , t ) | for a bath of N B =
32 spins, for different choices of the bathHamiltonian, and with and without system-bath interaction. The calculation of the nine correlations Eq. (47) requires solvingfour TDSEs simultaneously, using as the initial states the random thermal state | Y ( b ) i , B | Y ( b ) i , B | Y ( b ) i , and B | Y ( b ) i .As the whole system contains 33 spins, these calculations are fairly expensive in terms of CPU and memory cost. One suchcalculation needs somewhat less than 1TB memory to run and takes about 5 hours using 65536 BlueGene/Q processors which,in practice, limits the time interval that can be studied.In all four cases, the absolute values of correlations for i = j are much smaller than those for i = j and have therefore beenomitted in Fig. 2. The remaining three correlations decay rapidly but, on the time scale shown, are definitely non-zero at t = K a n in the range [ − / , / ] , we have h| K a n |i ≈ /
8. Onthe other hand, for the antiferromagnetic Heisenberg bath we have K = − / t >
10 thecorrelations of the former seem to have reached a stationary state whereas in the case of the latter, they do not. Moreover, usingthe full Hamiltonian ( l = .
1) instead of only the bath Hamiltonian to solve the TDSE, for t >
10 the changes to the correlationsare less pronounced if the bath is an antiferromagnetic Heisenberg model than if the bath has random interactions. Based onthese results, it seems advantageous to adopt the antiferromagnetic Heisenberg model ( K = − /
4) as the Hamiltonian of thebath.Qualitatively, in all cases, the correlations are either small for all t or decrease by about order of magnitude on a short timescale ( t < V. ALGORITHM TO EXTRACT e t A AND B FROM TDSE DATA
Recall that our primary objective is to determine the Markovian master equation Eq. (37) which gives the best (in the least-square sense) fit to the solution of the TDSE. Obviously, this requires taking into account the full motion of the system spin, notonly the decay envelope, over an extended period of time.The numerical solution of the TDSE of the full problem yields the data r k ( t ) = h s k ( t ) i . In this section, we consider thesedata as given and discuss the algorithm that takes as input the values of r k ( t ) and returns the optimal choice of the matrix e t A and vector B , meaning that we minimize the least-square error between the data { r k ( t ) } and the corresponding data, obtainedby solving Eq. (39).Denoting r k ( n ) ≡ r k ( n t ) , it follows that if Eq. (39) is assumed to hold, we must have r ( ) r ( ) . . . r ( N ) r ( ) r ( ) . . . r ( N ) r ( ) r ( ) . . . r ( N ) = ( e t A ) ( e t A ) ( e t A ) ( B ) ( e t A ) ( e t A ) ( e t A ) ( B ) ( e t A ) ( e t A ) ( e t A ) ( B ) r ( ) r ( ) . . . r ( N − ) r ( ) r ( ) . . . r ( N − ) r ( ) r ( ) . . . r ( N − ) , (48)where N is the number of time steps for which the solution of the TDSE is known. We may write Eq. (48) in the more compactform Z = YX , (49)where Z is a 3 × N matrix of data, Y is a 3 × X is a 4 × N matrix of data.We determine Y by solving the linear least square problem, that is we search for the solution of the problem min Y || Z − YX || .A numerically convenient way to solve this minimization problem is to compute the singular value decomposition [48, 49] of2 X = U SSS V T where U is an orthogonal 3 × SSS is the 3 × N matrix with the singular values of X on its diagonal, and V T is an orthogonal N × N matrix. In terms of these matrices we have Y = ZV SSS + U T , (50)where SSS + is the pseudo-inverse of SSS , which is formed by replacing every non-zero diagonal entry of
SSS by its reciprocal andtransposing the resulting matrix.Numerical experiments show that the procedure outlined above is not robust: it sometimes fails to reproduce the known e t A and B =
0, in particular in the case that e t A is (close to) an orthogonal matrix. Fortunately, a straightforward extension rendersthe procedure very robust. The key is to use data from three runs with different initial conditions. This also reduces the chancethat the estimates of e t A and B are good by accident. In practice, we take the initial states to be orthogonal (see Sec. VI for theprecise specification).Labeling the data for different initial states by superscripts we have (cid:0) Z ( ) Z ( ) Z ( ) (cid:1) = Y (cid:0) X ( ) X ( ) X ( ) (cid:1) , (51)but now Z = ( Z ( ) Z ( ) Z ( ) ) and X = ( X ( ) X ( ) X ( ) ) are 3 × N and 4 × N matrices of data, respectively. Using Eq. (50) wecompute Y = ( e t A ) ( e t A ) ( e t A ) ( B ) ( e t A ) ( e t A ) ( e t A ) ( B ) ( e t A ) ( e t A ) ( e t A ) ( B ) , (52)from which the matrix e t A and vector B immediately follow. In Appendix A, we discuss the method that we used to validate theextraction method. VI. FITTING A QUANTUM MASTER EQUATION TO THE SOLUTION OF THE TDSE
The procedure to test the hypothesis as to whether the QMEQ Eq. (37) provides a good approximation to the exact TDSE ofa (small) system which is weakly coupled to a (large) environment can be summarized as follows:1. Make a choice for the model parameters h x B , h z B , K , D , and the system-bath interaction l , for the number of bath spins N B ,the inverse temperature b of the bath, and the time step t ( t = | Y ( ) i x = | x i| f i , | Y ( ) i y = | y i| f i , and | Y ( ) i z = | ↑i| f i where | x i = ( | ↑i + | ↓i ) / √ | y i =( | ↑i + i | ↓i ) / √
2, and | f i denotes a pure state picked randomly from the 2 N B -dimensional unit hypersphere. For eachof the three initial states we may or may not use different realizations of | f i . If b >
0, prepare typical thermal statesby projection [47], that is set | Y ( ) i x = | x i| f ( b / ) i / h f ( b / ) | f ( b / ) i / (and similarly for the two other initial states)where | f ( b / ) i = e − b H B / | f i .3. For each of the three initial states, solve the TDSE for 0 ≤ t = n t ≤ T = N t . The case of interest is when T is largeenough for the system-bath to reach a steady state. For each of the three different initial states compute r i , j ( k ) ≡h Y ( k t ) | s i | Y ( k t ) i j , for i , j = x , y , z and store this data.4. Use the data r i , j ( k ) to construct the 3 × N matrix Z = ( Z ( ) Z ( ) Z ( ) ) and 4 × N matrix X = ( X ( ) X ( ) X ( ) ) [seeEq. (51)] and compute the 3 × Y , yielding the best (in the least-square sense) estimates of e t A and B .5. Use the estimates of e t A and B to compute the averages [denoted by e r i , j ( k ) ] of the three components of the systemspin operators sss ( t ) , according to Eq. (39) for each of the three different initial states. Quantify the difference of thereconstructed data, i.e. the solution of the “best” approximation in terms of the QMEQ, and the original data obtained bysolving the TDSE by the number e max ( t = k t ) = max i , j | r i , j ( k ) − e r i , j ( k ) | . (53)6. Check if the approximate density matrix of the system, defined by e r i , j ( k ) , is non-negative definite. In none of our simula-tion runs the approximate density matrix of the system failed this test.3 -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -4 t(a) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s z ( t ) > -4 t(b) FIG. 3. (color online) Comparison between the spin averages as obtained by solving the TDSE (solid lines) and the QMEQ (solid circles) with e t A and B extracted from the TDSE data. (a) initial state | y i| f i ; (b) initial state | ↑i| f i . The model parameters are: l = N B = b = K = − / D = h x B = h z B = /
8. For clarity, the system-spin averages are shown with a time interval of 100. The markers represent thedata obtained by least-square fitting to 15000 numbers generated by the TDSE solver.
Test of the procedure to fit Eq. (37) to TDSE data
If the system does not interact with the bath ( l = H = ( h x , , ) . Therefore, the l = y - and z -components of the system spin as obtained by solving the TDSE withinitial states | y i| f i and | ↑i| f i , respectively. Looking at the time interval shown in Fig. 3 and recalling that the spin componentsperform oscillations with a period p / h x , it is clear that Fig. 3 does not show these rapid oscillations. Instead, not to clutterthe plots too much, we only plotted the values at regular intervals, as indicated by the markers. For the initial state | x i| f i , the x -component is exactly constant (both for the TDSE and time evolution using the estimated e t A and B ) and therefore not shown.The difference between the spin averages obtained from the TDSE and from time evolution according to Eq. (39) (using theestimated e t A and B ) is rather small ( e max ( t ) < − for 0 ≤ t ≤ e max ( t ) are reflected in the excellent agreement between the TDSE and QMEQ (Eq. (37)) data shown inFig. 3. From these simulation data we conclude that for l =
0, the matrix e t A and vector B obtained by least-square fitting tothe TDSE data define a QMEQ that reproduces the correct values of the spin averages.The next step is to repeat the analysis for the case of weak system-bath interaction l = .
05 (recall that we already foundthat l = . e t A and B using the data of three different solutions of the TDSE. It does not fit data for individual spin componentsseparately nor does it fit data obtained from a TDSE solution of one particular choice of the initial state. Our procedure yieldsthe best globalestimates for e t A and B in the least-square sense.In Fig. 4 we illustrate the procedure for sampling and processing the TDSE data and for plotting these data along with thedata obtained from Eq. (39) using the estimated e t A and B . We present data for short times (top figures) and for the wholetime interval (bottom figures). The TDSE data (solid line) is being sampled, namely at times indicated by the t -values of themarkers, which in the case corresponds to a time steps of 0.2 [see Figs. 4(a)–4(c)]. The sampled data of the whole interval [ , ] are used to determine e t A and B by the least-square procedure described in Sec. V. In this particular case, the TDSEsolver supplies 15000 numbers to the least-square procedure. The estimated e t A and B thus obtained are then used to computethe time-evolution of the spin components, the data being represented by the markers.From Fig. 4(d), it is clear that although the QMEQ produces the correct qualitative behavior of the x -component of the systemspin, the difference with the TDSE data is significant (as is also clear from e max ( t ) ). In particular, the TDSE data of the x -component of the system spin do not show relaxation to the thermal equilibrium value, which is zero for b =
0. At first sight,this could be a signature that the fitting procedure breaks down because it is certainly possible to produce a much better fit to theTDSE data of the x -component ifwewouldfitacurvetothisdataonly. But, as explained above, we estimate e t A and B by fittingto the nine (three spin components × three different initial states) of such curves simultaneously. Apparently, the mismatch inthe x -component is compensated for by the close match of the y –component [see Fig. 4(e) and z -component (not shown)].Remarkably, the matrix e t A and vector B extracted from the TDSE data yield a QMEQ that does indicate that the systemspin relaxes to a state that is close to thermal equilibrium: The QMEQ yields a value of 0.04 for the expectation value of the4 -1-0.5 0 0.5 1 0 5 10 15 20 < s x ( t ) > t(a) -1-0.5 0 0.5 1 0 5 10 15 20 < s y ( t ) > t(b) -1-0.5 0 0.5 1 0 5 10 15 20 < s z ( t ) > t(c)-1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -3 t(d) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -3 t(e) 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0 0.2 0.4 0.6 0.8 1 e m a x ( t ) -3 t(f) FIG. 4. (color online) Comparison between the spin averages as obtained by solving the TDSE (solid lines) and the QMEQ (solid circles)with e t A and B extracted from the TDSE data. (a)–(c) show how the TDSE data (solid line) are being sampled, namely at times indicated bythe t -values of the markers, which in the case corresponds to a time steps of 0.2. (d)– (f): the sampled data of the whole interval [ , ] , inthis case 15000 numbers, are used to determine by the least-square procedure described in Sec. V, the parameters that enter the time-evolutionof the Markovian master equation Eq. (37). The latter is then used to compute the time-evolution of the spin components, the data beingrepresented by the markers. For clarity, in the bottom figures, the data are shown with a time interval of 10. The model parameters are: h x B = h z B = / l = . N B = b = K = − /
4, and D =
1. (a),(d) initial state | x i| f i ; (b),(e) initial state | y i| f i ; (c) initial state | z i| f i ;(f) the error e max ( t ) . x -component of the system spin and values less than 10 − for the other two components. From the general theory of the QMEQin the Markovian approximation [4], we know that if the correlations of the bath-operators Eq. (47) satisfy the Kubo-Martin-Schwinger condition, the stationary state solution of the QMEQ is exactly the same as the thermal equilibrium state of the system(ignoring corrections of O ( l ) [see Ref. 20 for a detailed discussion)].The mismatch between the QMEQ and TDSE data of the x -component can be attributed to the fact that a bath of N B =
13 spinsis too small to act as a bath in thermal equilibrium. However, the argument that leads to this conclusion is somewhat subtle. Asshown in Sec. III B, the random state approach applied to the system + bath yields the correct thermal equilibrium properties. Inparticular, in the case at hand ( b = N B = h F ( b = ) | s a ( t ) | F ( b = ) i ≈ a = x , y , z . Note that in this kind of calculation, the initial state | F ( b = ) i is a random state of the system + bath. In contrast,the data shown in Fig. 4(d) are obtained by solving the TDSE with the initial state | Y ( ) i x = | x i| f i (see Sec. VI). Therefore, theresults of Fig. 4(d) demonstrate that for N B =
13, the statement that | x i| f i −→ TDSE evolution −→ | e F i , (54)where e F denotes an (approximate) random state of the whole system, is not necessarily true. Otherwise, we would have h e F | s x ( t ) | e F i ≈ t large enough, in contradiction with the data shown in Fig. 4(d). Roughly speaking, one could say that abath of N B =
13 is not sufficiently “complex” to let the TDSE evolve certain initial states towards a random state of the wholesystem. For a discussion of the fact that in general, Eq. (54) does not necessarily hold, see Ref. 36.As a check on this argument, we repeat the simulation with a bath N B =
24 spins. The results are shown in Fig. 5. ComparingFigs. 4 and 5, it is clear that for long times the value of the x -component decreases as the number of spins in the bath increasesand that the agreement between the TDSE data and the fitted QMEQ data has improved considerably. This suggests that as thesize of the bath increases and with the bath initially in a random state, the TDSE evolution can drive the state to an (approximate)random state of the whole system, meaning that the whole system relaxes to the thermal equilibrium state. However, as discussedin Sec. IX there are exceptions [36].In general, we may expect that for short times, a Markovian QMEQ cannot represent the TDSE evolution very well [8, 14].But if we follow the evolution for times much longer than the typical correlation times of the bath-operators, the differencebetween the QMEQ and TDSE data for short times does not affect the results of fitting the data over the whole, large time-5 -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -3 t(a) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -3 t(b) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0 0.2 0.4 0.6 0.8 1 e m a x ( t ) -3 t(c) FIG. 5. (color online) Same as Fig. 4, except that the bath contains N B =
24 spins and l = .
05. The markers represent the data obtained byleast-square fitting to 15000 numbers generated by the TDSE solver. For clarity the data is shown with a time interval of 6. -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -3 t(a) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -3 t(b) 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0 0.2 0.4 0.6 0.8 1 e m a x ( t ) -3 t(c) FIG. 6. (color online) Same as Fig. 5, except that the bath is initially at b = interval in a significant manner. Hence there is no need to discard the short-time data in the fitting procedure. As a matter offact, the data shown in Fig. 4 indicate that the least-square procedure applied to the whole data set yields a Markovian masterequation that reproduces the short-time behavior quite well.Finally, we check that the conclusions reached so far for a bath at b = b >
0. In Fig. 6, we show thesimulation results for b =
1, for the same system and bath as the one used to obtain the data shown in Fig. 5. From Fig. 6 weconclude that the agreement between the TDSE and QMEQ data is quite good.
TABLE I. The parameters that appear in Eq. (55) as obtained by fitting the QMEQ to the TDSE data shown in Fig. 7. b i A i , A i , A i , b i − . × − + . × − − . × − − . × − − . × − − . × − + . − . × − − . × − − . − . × − − . × − − . × − + . × − − . × − − . × − − . × − − . × − + . − . × − − . × − − . − . × − − . × − − . × − + . × − + . × − − . × − − . × − − . × − + . − . × − − . × − − . − . × − − . × − VII. SIMULATION RESULTS: N B = , As already mentioned in Sec. III, in practice, there is a limitation on the sizes and time intervals that can be explored. Byincreasing the system-bath interaction l , we can shorten the time needed for the system to relax to equilibrium. On the other6 -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -3 t(a) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -3 t(b) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -3 t(c)-1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -3 t(d) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -3 t(e) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -3 t(f) FIG. 7. (color online) Simulation data for a bath with N B =
28 spins and system-bath interaction l = .
1. The model parameters are: h x B = h z B = / K = − /
4, and D = h s x ( t ) i as obtained by starting fromthe initial state | x i| f i , (a)–(c) corresponding to b = , ,
2, respectively. Bottom row: h s y ( t ) i as obtained by starting from the initial state | y i| f i , (d)–(f) corresponding to b = , ,
2, respectively. The TDSE simulations yield h| s x ( t = ) |i = . h| s x ( t = ) |i = . h| s x ( t = ) |i = .
756 for b = , ,
2, respectively, whereas for the system in equilibrium we have h s x i = , . , .
762 for b = , , hand, l should not be taken too large because when we leave the perturbative regime, the QMEQ of the form Eq. (37) cannotbe expected to capture the true quantum dynamics. From our exploratory simulations, we know that l = . N B =
32 spins.In Fig. 7 we present the results as obtained with a bath containing N B =
28 spins, prepared at b = , ,
2. Although Fig. 7 maysuggest otherwise, the maximum error max k e max ( t ) ≈ . , . , . b = , ,
2, respectively, indicating that the differencebetween the TDSE data and the QMEQ approximation increases with b . The results presented in Fig. 8 for a bath of N B = b = N B = ,
32 spins are sufficiently large to mimican infinite thermal bath. At any rate, in all cases, there is very good qualitative agreement between the TDSE and QMEQ data.From the TDSE data, we can, of course, also extract the values of the entries in the matrix A and vector b , see Eq. (37).Writing Eq. (37) more explicitly as ¶ h sss x ( t ) i ¶ t = A , h sss x ( t ) i + A , h sss y ( t ) i + A , h sss z ( t ) i + b ¶ h sss y ( t ) i ¶ t = A , h sss x ( t ) i + A , h sss y ( t ) i + A , h sss z ( t ) i + b ¶ h sss z ( t ) i ¶ t = A , h sss x ( t ) i + A , h sss y ( t ) i + A , h sss z ( t ) i + b , (55)and using, as an example, the data shown in Fig. 7, we obtain the values of the coefficients as given in Table I. From Table I, wereadily recognize that (i) A , ≈ − A , ≈ h x = /
2, (ii) thereis a weak coupling between the x - and ( y , z ) - components of the system spin and (iii) the three spin components have differentrelaxations times.As a final check whether l = . N B =
32 spins and system-bath interaction l = . b =
0. The simulation data are presented in Fig. 9. Clearly, there still isgood qualitative agreement between the TDSE and QMEQ data but, as expected, max k e max ( t ) has become larger (by a factor ofabout 3).In Table II we present results (first three rows) for the least-square estimates of the parameters that enter the QMEQ, asobtained from the TDSE data shown in Fig. 8. Taking into account that with each run, the random values of the model parameters7 -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -3 t(a) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -3 t(b) 0 0.002 0.004 0.006 0.008 0.01 0.012 0 0.2 0.4 0.6 0.8 1 e m a x ( t ) -3 t(c) FIG. 8. (color online) Simulation data for a bath with N B =
32 spins prepared at b = l = .
1. The modelparameters are: h x B = h z B = / K = − /
4, and D =
1. Figures (a,b) show TDSE data (solid lines) and QMEQ data (solid circles). (a) initialstate | x i| f i ; (b) initial state | y i| f i ; (c) the error e max ( t ) . The data obtained with the initial state | ↑i| f i is very similar as the data obtained withthe initial state | y i| f i and are therefore not shown. For clarity, the data are shown with a time interval of 0.4. The TDSE solver provided 3000numbers as input to the least-square procedure. -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -2 t(a) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -2 t(b) 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0 0.2 0.4 0.6 0.8 1 e m a x ( t ) -2 t(c) FIG. 9. (color online) Same as Fig. 8 except that b = l = . change, the order-of-magnitude agreement between the data for N B =
28 (Table I, rows 4–6) and the N B =
32 data is rather good.We also present results (middle and last three rows) for the parameters that enter the Redfield equation Eq. (46), as obtained fromthe TDSE data of the bath-operator correlations C ( i , j , t ) for 0 ≤ t ≤
40 [see Figs. 2(a) and 2(c) for a picture of some of thesedata]. From Table II, it is clear that there seems to be little quantitative agreement between a description based on the Redfieldquantum master equation (46) obtained by using the bath-operator correlations C ( i , j , t ) data and the parameters obtained fromthe least-square fit of Eq. (55) to the TDSE data. Simulations using the 3D bath Hamiltonian (7) support this conclusion (seeAppendix B).Although our results clearly demonstrate that QMEQ Eq. (37) quantitatively describes the true quantum dynamics of a spininteracting with a spin bath rather well, the Redfield quantum master equation Eq. (46) in the Markovian approximation, whichis also of the form Eq. (37), seems to perform rather poorly in comparison. The estimates of the diagonal matrix elements ofthe matrix A as obtained from the expressions in terms of the bath-operator correlations C ( i , j , t ) are too small by factors 3–7.This suggests that the approximations involved in the derivation of Eq. (46) are not merely of a perturbative nature but affect thedynamics in a more intricate manner [see Ref. 19 for an in-depth discussion of these aspects]. VIII. EXCEPTIONS
The simulation results presented in Secs. VI and VII strongly suggest that, disregarding some minor quantitative differences,the complicated Schr¨odinger dynamics of the system interacting with the bath can be modeled by the much simpler QMEQ ofthe form (37). But, as mentioned in Sec. IV, there are several approximations involved to justify the reduction of the Schr¨odingerdynamics to a QMEQ. In this section, we consider a few examples for which this reduction may fail.The first case that we consider is defined by the Hamiltonian H = − h x s x + l N B (cid:229) n = (cid:0) s xn s x + s yn s y + s zn s z (cid:1) + N B (cid:229) n = (cid:0) s xn s xn + + s yn s yn + + s zn s zn + (cid:1) . (56)8 TABLE II. First three data rows: coefficients that appear in Eq. (55) as obtained by fitting the QMEQ to the TDSE data shown in Fig. 8. Middlethree rows: the corresponding coefficients as obtained by numerically calculating the parameters r jk that appear in the Redfield quantum masterequation Eq. (46) according to Eq. (45), using the TDSE data of the bath-operator correlations shown in Fig. 2(a). Last three rows: same as themiddle three rows except that the used TDSE data of the bath-operator correlations are shown in Fig. 2(c). Note that the baths used in thesesimulations are very different (see Fig. 2), yet the relevant numbers (those with absolute value larger than 10 − ) are in the same ballpark. i A i , A i , A i , b i − . × − + . × − − . × − − . × − − . × − − . × − + . − . × − − . × − − . − . × − − . × − − . × − − . × − + . × − − . × − − . × − − . × − + . − . × − + . × − − . − . × − + . × − − . × − + . × − + . × − − . × − + . × − − . × − + . + . × − − . × − − . − . × − + . × − -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -2 t(a) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -2 t(b) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0 0.2 0.4 0.6 0.8 1 e m a x ( t ) -2 t(c) FIG. 10. (color online) Simulation data for a bath with N B =
32 spins prepared at b = l = .
2. The systemHamiltonian is given by Eq. (2). The system-bath interaction is given by Eq. (3) with J xn = J yn = J zn = /
4. The bath Hamiltonian is givenby Eq. (5) with K = − / D = h xn = h zn =
0. The full Hamiltonian is given by Eq. (56). Figures (a,b) show TDSE data (solid lines)and QMEQ data (solid circles). (a) initial state | x i| f i ; (b) initial state | y i| f i ; (c) the error e max ( t ) . The data obtained with the initial state | ↑i| f i is very similar as the data obtained with the initial state | y i| f i and are therefore not shown. With this choice of parameters of bath andsystem-bath Hamiltonians, the system does not relax to its thermal equilibrium state lim t → ¥ h s x ( t ) i = lim t → ¥ h s y ( t ) i = lim t → ¥ h s z ( t ) i = In other words, both the system-bath and intra-bath interactions are of the isotropic antiferromagnetic Heisenberg type and allinteraction strengths are constant. The simulation results for this case are presented in Fig. 10. From Fig. 10(a), it is immediatelyclear that the system does not relax to its thermal equilibrium state at b = t → ¥ h s x ( t ) i = H = − h x s x + l N B (cid:229) n = J zn s zn s z + N B (cid:229) n = s zn s zn + , (57)with system-bath interactions J zn chosen at random and distributed uniformly over the interval [ − , ] and the bath is modeled byan Ising Hamiltonian. The model Eq. (57) is known to exhibit quantum oscillations in the absence of quantum coherence [50].As the bath Hamiltonian commutes with all other terms of the Hamiltonian, the only non-zero bath correlation C ( , , t ) isconstant in time, hence one of the basic assumptions in deriving the QMEQ Eq. (37) does not hold.Because of the special structure of the Hamiltonian Eq. (57) it is straightforward to compute closed form expressions for theexpectation values of the system spin. For b = z h s x ( t ) i = − l ** B sin t p ( h x ) + B ( h x ) + B ++ , | Y ( t = ) i = | x i| f i , (58)9 -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -2 t (a) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -2 t (b) 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0.2 0.4 0.6 0.8 1 e m a x ( t ) -2 t (c) FIG. 11. (color online) Simulation data for a bath with N B =
32 spins prepared at b = l = .
2. The systemHamiltonian is given by Eq. (2). The system-bath interaction is given by Eq. (3) with J xn = J yn = J zn uniformly random between − / /
4, in which case the interaction of the system and bath spins is through the coupling of the z -components of the spins only. Thebath Hamiltonian is given by Eq. (6) with K xn = K yn = h xn = h zn = K zn uniformly random between − | x i| f i ; (b) initial state | y i| f i ;(c) the error e max ( t ) . The data obtained with the initial state | ↑i| f i are very similar as the data obtained with the initial state | y i| f i andare therefore not shown. With this choice of bath and system-bath Hamiltonians, the system does not relax to its thermal equilibrium statelim t → ¥ h s x ( t ) i = lim t → ¥ h s y ( t ) i = lim t → ¥ h s z ( t ) i =
0. For clarity, the data are shown with a time interval of 0.4. The TDSE solver provided3000 numbers as input to the least-square procedure. h s y ( t ) i = (cid:28)(cid:28) cos 2 t q ( h x ) + B (cid:29)(cid:29) , | Y ( t = ) i = | y i| f i , (59) h s z ( t ) i = − l ( h x ) ** sin t p ( h x ) + B ( h x ) + B ++ , | Y ( t = ) i = | ↑i| f i , (60)where B = B ( { s n } ) = (cid:229) N B n = J zn s n and hh X ii ≡ (cid:229) { s = ± } . . . (cid:229) { s N B = ± } |h s . . . s N B | f i| X ( { s n } ) , (61)denotes the average over all the bath-spin configurations.From Eq. (58) it follows immediately that if the system+bath is initially in the state | Y ( t = ) i = | x i| f i , we must have h s x ( t ) i ≥ − l . Hence the system will never relax to its thermal equilibrium state [for which lim t → ¥ h s x ( t ) i = N B , by assuming | f i to be a uniform superposition of the 2 N B different bath states and by approximating B , beinga sum of independent uniform random variables, by a Gaussian random variable. Then we have (after substituting B = h x u ) h s y ( t ) i = s √ p Z + ¥ − ¥ du e − u q / s cos (cid:16) th x p + u (cid:17) . (62)For large t , we can evaluate Eq. (62) by the stationary phase method and we find that h s y ( t ) i decays as 1 / √ t . Such a slowalgebraic decay cannot result from a time evolution described by a single matrix exponential e t A . In other words, the apparentagreement shown in Fig. (11) is due to the relatively short time interval covered. On the other hand, as already mentioned, themodel defined by Eq. (57) is rather exceptional in the sense that the bath correlations do not exhibit any dynamics. Hence it isnot a surprise that the QMEQ cannot capture the 1 / √ t dependence.Finally, in Fig. 12 we illustrate what happens if the l =
1, that is if the system-bath interaction becomes comparable to theother energy scales h x and K . Then, the perturbation expansion that is used to derive the QMEQ of the form Eq. (37) is nolonger expected to hold [4]. The data presented in Fig. 12 clearly show that even though the time it takes for the system to reachthe stationary state is rather short (because l = h s x ( t = ) i = .
264 [ |h s z ( t = ) i| ≤ − |h s z ( t = ) i| ≤ − ],whereas from statistical mechanics for the isolated system at b = h s x i = tanh ( / ) = .
462 [ h s y ( t = ) i = h s z ( t = ) i = -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s x ( t ) > -1 t (a) -1-0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 < s y ( t ) > -1 t (b) 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0.2 0.4 0.6 0.8 1 e m a x ( t ) -1 t (c) FIG. 12. (color online) Same as Fig. 8 except that b = l = IX. SUMMARY
We have addressed the question to what extent a quantum master equation of the form 37) captures the salient features of theexact Schr¨odinger equation dynamics of a single spin coupled to a bath of spins. The approach taken was to solve the time-dependent Schr¨odinger equation of the whole system and fit the data of the expectation values of the spin components to thoseof a quantum master equation of the form (37).In all cases in which the approximations used to derive a quantum master equation of the form (37) seem justified, it wasfound that the quantum master equation (37) extracted from the solutions of the time-dependent Schr¨odinger equation describesthese solutions rather well. The least-square procedure that is used to fit the quantum master equation (37) data to the time-dependent Schr¨odinger data accounts for non-Markovian effects and nonperturbative contributions. Quantitatively, we found thatdifferences between the data produced by the quantum master equation, obtained by least-square fitting to the time-dependentSchr¨odinger data, and the latter data increases with decreasing temperature.The main finding of this work is that the exact Schr¨odinger dynamics of a single spin-1/2 object interacting with a spin-1/2bath can be accurately and effectively described by Eq. (37) which, for convenience of the reader, is repeated here and reads as ¶r rr ( t ) ¶ t = A rrr ( t ) + b , (63)where the 3 × A and the three elements of the vector b are time independent. As the mathematical structure of the(Markovian) quantum master equation (63) is the same as that of the Bloch equation (A1), as a phenomenological description,the quantum master equation (63) offers no advantages over the latter. Of course, when the system contains more than onespin, the Bloch equation can no longer be used whereas the quantum master equation (63) still has the potential to describe thedynamics. We relegate the assessment of the quantum master equation approach to systems of two or more spins to a futureresearch project. ACKNOWLEDGEMENTS
We thank Dr. Takashi Mori for valuable discussions. Work of P.L. Zhao is supported by the China Scholarship Council(No.201306890009). Work of S.M. was supported by Grants-in-Aid for Scientific Research C (25400391) from MEXT ofJapan, and the Elements Strategy Initiative Center for Magnetic Materials under the outsourcing project of MEXT. The au-thors gratefully acknowledge the computing time granted by the JARA-HPC Vergabegremium and provided on the JARA-HPCPartition part of the supercomputer JUQUEEN [51] at Forschungszentrum J¨ulich.
Appendix A: Bloch equations
Whatever method we use to extract e t A and B , it is necessary to validate the method by applying it to a non-trivial problemfor which we know the answer for sure. The Bloch equations, originally introduced by Felix Bloch [52] as phenomenologicalequations to describe the equations of motion of nuclear magnetization, provide an excellent test bed for the extraction algorithmpresented in Sec. V.In matrix notation the Bloch equations read as d M ( t ) dt = b AM ( t ) + b b , (A1)1where M is the magnetization, b A = − / T h z − h y − h z − / T h x h y − h x − / T , (A2)and b b = M / T where M is the steady state magnetization. The transverse and longitudinal relaxation times T and T arestrictly larger than zero. The special but interesting case in which there is no relaxation corresponds to 1 / T = / T = rrr ((( ttt ))) = M ( t ) that is neededto test the algorithm described in Sec. V. In order that the identification rrr ((( ttt ))) = M ( t ) makes sense in the context of the quantummaster equation we have to impose the trivial condition that k M ( t = ) k ≤ k M k ≤ e t b A using the second-order product-formula [53] e t b A ≈ e t e A = (cid:16) e t A / m e t A / m e t A / m (cid:17) m , (A3)where b A = A + A and A = − / T − / T
00 0 − / T , (A4) A = h z − h y − h z h x h y − h x . (A5)The second-order product-formula approximation satisfies the bound k e t b A − e t e A k ≤ c t / m where the constant c = O ( k [ A , A ] k ) .Hence the error incurred by the approximation is known and can be reduced systematically by increasing m .It is straightforward to compute the closed form expressions of the matrix exponentials that appear in the second-order product-formula. We have e t A = e − t / T e − t / T
00 0 e − t / T e t A = W h x + ( h y + h z ) cos t W h x h y ( − cos t W ) + h z W sin t W h x h z ( − cos t W ) − h y W sin t W h x h y ( − cos t W ) − h z W sin t W h y + ( h x + h z ) cos t W h y h z ( − cos t W ) + h x W sin t W h x h z ( − cos t W ) + h y W sin t W h y h z ( − cos t W ) − h x W sin t W h z + ( h x + h y ) cos t W , (A6)where W = h x + h y + h z .Summarizing, the numerical solution of the Bloch equations Eq. (A1) is given by rrr ( t + t ) = e t e A rrr ( t ) + e B , (A7)where rrr ( t ) = M ( t ) and the trapezium rule was used to write b B = Z t e ( t − u ) b A b b du ≈ t (cid:16) + e t e A (cid:17) b b = e B . (A8)The approximate solution obtained from Eqs. (A7) and (A8) will converge to the solution of Eq. (A1) as t →
0. Clearly, Eq. (A7)has the same structure as Eq. (39) and hence we can use the solution of the Bloch equations as input data for testing the extractionalgorithm. Note that the extraction algorithm is expected to yield e t e A and e B , not e t b A and b B .
1. Validation procedure
We use the Bloch equation model to generate the data set D = { rrr ( k t ) | ≤ k ≤ N − } . The validation procedure consists ofthe following steps:1. Choose the model parameters h x , h y , h z , 1 / T , 1 / T and the steady-state magnetization M .22. Choose t and m .3. For each of the three initial states r ( ) ( ) = ( , , ) T , r ( ) ( ) = ( , , ) T , and r ( ) ( ) = ( , , ) T repeat the operation r ( j ) (( k + ) t ) ← e t e A r ( j ) ( k t ) + e B , k = , . . . , N − , j = , , , and store these data.4. Use the data { r ( j ) ( k t ) } to construct the matrices 3 × N matrix Z = ( Z ( ) Z ( ) Z ( ) ) and the 4 × N matrix X =( X ( ) X ( ) X ( ) ) . Then use the singular value decomposition of X to compute the matrix Y according to Eq. (50) andextract the matrix e t A and vector B from it, see Eq. (48). If one or more of the singular values are zero, the extractionfailed.5. Compute the relative errors e A = k e t e A − e t b A k / k e t b A k , (A9) e B = k e B − b B k / k b B k , (A10) e rrr = max k k rrr (( k + ) t ) − e t e A rrr ( k t ) − e B k / k rrr ( k t ) k . (A11)A necessary condition for the algorithm to yield reliable results is that the errors e A and e B are small, of the order of 10 − .Indeed, if one or more of the singular values are zero and the extraction has failed, e rrr may be (very) small but e A or e B is not.In the case that is of interest to us, the case in which the whole system evolves according to the TDSE, we do not know e t A nor B and a small value of e rrr is, by itself, no guarantee that the extraction process worked properly. Hence, it also is importantto check that all singular values are nonzero.
2. Numerical results
In Table III we present some representative results for the errors incurred by the extraction process. In all cases, the relativeerrors on the estimate of the time evolution operator and the constant term are for the present purpose, rather small. Therefore,the algorithm to extract the time evolution operator e t e A and constant term e B appearing in the time evolution equation Eq. (39)from the data obtained by solving the TDSE yields accurate results when the data are taken from the solution of the Blochequations. No exceptions have been found yet. TABLE III. The errors e A , e B , and e rrr as obtained fitting the matrix e t A and the constant term B , to the data of the numerical solution ofthe Bloch equation with three different initial conditions (see text). The Bloch equations are solved for N =
500 steps with the time step t .The value of the vector M = ( , , . ) T . The data of the whole time interval [ , N − ] were used for the least-square fitting procedure. Thecolumn labeled S i = h x h y h z / T / T t S i = e A < − e B < − e rrr < − . . . .
05 0 . X X X X . . . X X X X . . . .
01 0 1.0
X X X X
X X X X
Appendix B: Simulation results using the 3D bath Hamiltonian Eq. (7)
In this appendix, we present some additional results in support of the conclusions drawn from the simulations of using the 1Dbath Hamiltonians (5) and (6).Table IV summarizes the results of the analysis of TDSE data, as obtained with the 3D bath Hamiltonian Eq. (7) with randomintra-bath couplings and random h -fields for the bath spins. The model parameters that were used to compute the TDSE dataare the same as those that yield the results for the 1D bath presented in Table II. Comparing the first three rows (the parametersthat appear in the Markovian master equation (55) with the corresponding last three rows (the parameters r jk that appear in3 TABLE IV. First three data rows: coefficients that appear in Eq. (55) as obtained by fitting the QMEQ Eq. (37) to the TDSE data for h x = / l = . N B =
27, the 3D bath Hamiltonian Eq. (7) with random couplings ( K = /
4) and random h -fields ( h x B = h z B = / ) . Last three rows:the corresponding coefficients as obtained by numerically calculating the parameters r jk that appear in the Redfield quantum master equationEq. (46) according to Eq. (45) from the TDSE data of the bath-operator correlations. i A i , A i , A i , b i − . × − − . × − + . × − − . × − + . × − − . × − + . + . × − + . × − − . − . × − + . × − − . × − + . × − − . × − − . × − + . × − − . × − + . × + + . × − + . × − − . × + − . × − − . × − TABLE V. The same as Table IV except that the random couplings ( K = /
4) and random h-fields ( h x B = h z B = / ) . i A i , A i , A i , b i − . × − + . × − − . × − − . × − − . × − − . × − + . − . × − + . × − − . − . × − − . × − − . × − + . × − + . × − − . × − + . × − − . × − + . + . × − − . × − − . − . × − + . × − the Redfield quantum master equation Eq. (46)), we conclude that changing the connectivity of the bath does not significantlyimprove (compared to the data shown in Table I) the quantitative agreement between the data in the two sets of three rows.In Table V, we show the effect of increasing the energy scale of the bath spins by a factor of 10, reducing the relaxation timesof the bath-correlations by a factor of 10, i.e., closer to the regime of the Markovian limit in which Eq. (46) has been derived.The differences between the QMEQ estimates (first three rows) and the Redfield equation estimates (second three rows) valuesof A , and A , are significantly smaller than in those for the case shows in e.g. Table IV but the A , elements differ by a factorof four and the A ( , ) elements differ even much more. Although the results presented in Tables IV and V indicate that thedata extracted from the TDSE through Eq. (55) and those obtained by calculating the parameters r jk that appear in the Redfieldquantum master equation Eq. (46) in the Markovian limit will converge to each other, it becomes computationally very expensiveto approach that limit closer. The reason is simple: by increasing the energy-scale of the bath, it is necessary to reduce the timestep (or equivalently increase the number of terms in the Chebyshev polynomial expansion) in order to treat the fast oscillationsproperly. Keeping the same relaxation times roughly the same but taking a smaller time step requires more computation. Forinstance, it takes about 4 (20) h CPU time of 16384 BlueGene/Q processors to produce the TDSE data from which the numbersin Table IV (Table V) have been obtained. [1] A. Redfield, “On the theory of relaxation processes,” IBM J. Res. Develop. , 19–31 (1957).[2] S. Nakajima, “On quantum theory of transport phenomena,” Prog. Theor. Phys. , 948 – 959 (1958).[3] R. Zwanzig, “Ensemble method in the theory of irreversibility,” J. Chem. Phys. , 1338 – 1341 (1960).[4] H.-P. Breuer and F. Petruccione, Open Quantum Systems (Oxford University Press, Oxford, 2002).[5] W. Pauli,
Festschrift zum 60. Geburtstage A. Sommerfelds (Hirzel, Leipzig, 1928) p. 30.[6] M. B. Plenio and P. L. Knight, “The quantum-jump approach to dissipative dynamics in quantum optics,” Rev. Mod. Phys. , 101–144(1998).[7] G. Lindblad, “On the generators of quantum dynamical semigroups,” Commun. Math. Phys. , 119–130 (1976).[8] P. Gaspard and M. Nagaoka, “Slippage of initial conditions for the Redfield master equation,” J. Chem. Phys. , 5668 – 5675 (1999).[9] W.H. Louisell, Quantum Statistical Properties of Radiation (Wiley, New York, 1973).[10] A. Abragam,
Principles of Nuclear Resonance (Oxford Press, London, 1961).[11] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,
Atom-Photon Interactions (Wiley, New York, 1992).[12] Y. Hamano and F. Shibata, “Theory of exchange splitting in a strong magnetic field: I. General formulation,” J. Phys. C , 4843–4853(1984). [13] F. Shibata and M. Asou, “Theory of nonlinear spin relaxation. II,” J. Phys. Soc. Jpn. , 1234–1241 (1980).[14] A. Su´arez, R. Silbey, and I. Oppenheim, “Memory effects in the relaxation of quantum systems,” J. Chem. Phys. , 5101 – 5107 (1992).[15] P. Pechukas, “Reduced dynamics need not be completely positive,” Phys. Rev. Lett. , 1060–1062 (1994).[16] M. Sassetti and U. Weiss, “Correlation functions for dissipative two-state systems: Effects of the initial preparation,” Phys. Rev. A ,5383–5393 (1990).[17] U. Weiss, Quantum Dissipative Systems , 2nd ed. (World Scientific, Singapore, 1999).[18] Y. Tanimura, “Stochastic Liouville, Langevin, Fokker-Planck, and master equation approaches to quantum dissipative systems,” J. Phys.Soc. Jpn. , 082001–1–39 (2006).[19] H.-P. Breuer, J. Gemmer, and M. Michel, “Non-Markovian quantum dynamics: Correlated projection superoperators and Hilbert spaceaveraging,” Phys. Rev. E , 016139 (2006).[20] T. Mori and S. Miyashita, “Dynamics of the density matrix in contact with a thermal bath and the quantum master equation,” J. Phys.Soc. Jpn. , 124005 (2008).[21] M. Saeki, “Relaxation method and TCLE method of linear response in terms of thermo-field dynamics,” Physica A , 1827–1850(2008).[22] C. Uchiyama, M. Aihara, M. Saeki, and S. Miyashita, “Master equation approach to line shape in dissipative systems,” Phys. Rev. E ,021128 (2009).[23] T. Mori, “Natural correlation between a system and a thermal reservoir,” Phys. Rev. A , 040101(R) (2014).[24] H.-B. Chen, N. Lambert, Y.-C. Cheng, Y.-N. Chen, and F. Nori, “Using non-Markovian measures to evaluate quantum master equationsfor photosynthesis,” Sci. Rep. , 12753 (2015).[25] R. Kubo, “Statistical-mechanical theory of irreversible processes. I.” J. Phys. Soc. Jpn. , 570–586t (1957).[26] M. Nielsen and I. Chuang, Quantum Computation and Quantum Information , 10th ed. (Cambridge University Press, Cambridge, 2010).[27] M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, E. M.Chapple, C. Enderud, J. P. Hilton, K. Karimi, E. Ladizinsky, N. Ladizinsky, T. Oh, I. Perminov, C. Rich, M. C. Thom, E. Tolkacheva,C. J. S. Truncik, S. Uchaikin, J. Wang, B. Wilson, and G. Rose, “Quantum annealing with manufactured spins,” Nature , 194–198(2011).[28] H. Wang, “Basis set approach to the quantum dissipative dynamics: Application of the multiconfiguration time-dependent Hartree methodto the spin-boson problem,” J. Chem. Phys. , 9948–9956 (2000).[29] M. Nest and H.-D. Meyer, “Dissipative quantum dynamics of anharmonic oscillators with the multiconfiguration time-dependent Hartreemethod,” J. Chem. Phys. , 24–33 (2003).[30] D. Gelman and R. Kosloff, “Simulating dissipative phenomena with a random thermal phase wavefunctions, high temperature applicationof the surrogate Hamiltonian approach,” Chem. Phys. Lett. , 129–138 (2003).[31] D. Gelman, C.P. Koch, and R. Kosloff, “Dissipative quantum dynamics with the surrogate Hamiltonian approach. a comparison betweenspin and harmonic baths,” J. Chem. Phys. , 661–671 (2004).[32] G. Katz, D. Gelman, M.A. Ratner, and R. Kosloff, “Stochastic surrogate Hamiltonian,” J. Chem.Phys. , 034108 (2008).[33] K. De Raedt, K. Michielsen, H. De Raedt, B. Trieu, G. Arnold, M. Richter, Th. Lippert, H. Watanabe, and N. Ito, “Massively parallelquantum computer simulator,” Comp. Phys. Comm. , 121 – 136 (2007).[34] J. Lages, V. V. Dobrovitski, M. I. Katsnelson, H. A. De Raedt, and B. N. Harmon, “Decoherence by a chaotic many-spin bath,” Phys.Rev. E , 026225 (2005).[35] S. Yuan, M.I. Katsnelson, and H. De Raedt, J. Phys. Soc. Jpn. , 094003 (2009).[36] F. Jin, H. De Raedt, S. Yuan, M. I. Katsnelson, S. Miyashita, and K. Michielsen, “Approach to Equilibrium in Nano-scale Systems atFinite Temperature,” J. Phys. Soc. Jpn. , 124005 (2010).[37] S. Yuan, M.I. Katsnelson, and H. De Raedt, “Evolution of a quantum spin system to its ground state: Role of entanglement and interactionsymmetry,” Phys. Rev. A , 052109 (2007).[38] S. Yuan, M.I. Katsnelson, and H. De Raedt, “Decoherence by a spin thermal bath: Role of spin-spin interactions and initial state of thebath,” Phys. Rev. B , 184301 (2008).[39] F. Jin, K. Michielsen, M. A. Novotny, , S. Miyashita, S. Yuan, and H. De Raedt, “Quantum decoherence scaling with bath size: Importanceof dynamics, connectivity, and randomness,” Phys. Rev. A , 022117 (2013).[40] M. A. Novotny, F. Jin, S. Yuan, S. Miyashita, H. De Raedt, and K. Michielsen, “Quantum decoherence and thermalization at finitetemperatures within the canonical-thermal-state ensemble,” Phys. Rev. A , 032110 (2016).[41] H. Tal-Ezer and R. Kosloff, J. Chem. Phys. , 3967–3971 (1984).[42] C. Leforestier, R. H. Bisseling, C. Cerjan, M. D. Feit, R. Friesner, A. Guldberg, A. Hammerich, G. Jolicard, W. Karrlein, H.-D. Meyer,N. Lipkin, O. Roncero, and R. Kosloff, “A comparison of different propagation schemes for the time dependent Schr¨odinger equation,”J. Comput. Phys. , 59–80 (1991).[43] T. Iitaka, S. Nomura, H. Hirayama, X. Zhao, Y. Aoyagi, and T. Sugano, “Calculating the linear response functions of noninteractingelectrons with a time-dependent Schr¨odinger equation,” Phys. Rev. E , 1222–1229 (1997).[44] V. V. Dobrovitski and H. De Raedt, “Efficient scheme for numerical simulations of the spin-bath decoherence,” Phys. Rev. E , 056702(2003).[45] H. De Raedt and K. Michielsen, “Computational Methods for Simulating Quantum Computers,” in Handbook of Theoretical and Com-putational Nanotechnology , edited by M. Rieth and W. Schommers (American Scientific Publishers, Los Angeles, 2006) pp. 2 – 48.[46] L. E. Ballentine,
Quantum Mechanics: A Modern Development (World Scientific, Singapore, 2003).[47] A. Hams and H. De Raedt, “Fast algorithm for finding the eigenvalue distribution of very large matrices,” Phys. Rev. E , 4365 – 4377(2000).[48] G. H. Golub and C. F. Van Loan, Matrix Computations (John Hopkins University Press, Baltimore, 1996).[49] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling,
Numerical Recipes (Cambridge University Press, Cambridge, 2003). [50] V. V. Dobrovitski, H. De Raedt, M.I. Katsnelson, and B.N. Harmon, “Quantum oscillations without quantum coherence,” Phys. Rev.Lett. , 210401 (2003).[51] M. Stephan and J. Docter, “JUQUEEN: IBM Blue Gene/Q Supercomputer System at the Jlich Supercomputing Centre,” J. of Large-ScaleRes. Facil. A1 , 1 (2015).[52] F. Bloch, “Nuclear induction,” Phys. Rev. , 460–474 (1946).[53] M. Suzuki, “Decomposition formulas of exponential operators and lie exponentials with some applications to quantum mechanics andstatistical physics,” J. Math. Phys. , 601 – 612 (1985). | C ( .,. ) | tt
Numerical Recipes (Cambridge University Press, Cambridge, 2003). [50] V. V. Dobrovitski, H. De Raedt, M.I. Katsnelson, and B.N. Harmon, “Quantum oscillations without quantum coherence,” Phys. Rev.Lett. , 210401 (2003).[51] M. Stephan and J. Docter, “JUQUEEN: IBM Blue Gene/Q Supercomputer System at the Jlich Supercomputing Centre,” J. of Large-ScaleRes. Facil. A1 , 1 (2015).[52] F. Bloch, “Nuclear induction,” Phys. Rev. , 460–474 (1946).[53] M. Suzuki, “Decomposition formulas of exponential operators and lie exponentials with some applications to quantum mechanics andstatistical physics,” J. Math. Phys. , 601 – 612 (1985). | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt
Numerical Recipes (Cambridge University Press, Cambridge, 2003). [50] V. V. Dobrovitski, H. De Raedt, M.I. Katsnelson, and B.N. Harmon, “Quantum oscillations without quantum coherence,” Phys. Rev.Lett. , 210401 (2003).[51] M. Stephan and J. Docter, “JUQUEEN: IBM Blue Gene/Q Supercomputer System at the Jlich Supercomputing Centre,” J. of Large-ScaleRes. Facil. A1 , 1 (2015).[52] F. Bloch, “Nuclear induction,” Phys. Rev. , 460–474 (1946).[53] M. Suzuki, “Decomposition formulas of exponential operators and lie exponentials with some applications to quantum mechanics andstatistical physics,” J. Math. Phys. , 601 – 612 (1985). | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt
Numerical Recipes (Cambridge University Press, Cambridge, 2003). [50] V. V. Dobrovitski, H. De Raedt, M.I. Katsnelson, and B.N. Harmon, “Quantum oscillations without quantum coherence,” Phys. Rev.Lett. , 210401 (2003).[51] M. Stephan and J. Docter, “JUQUEEN: IBM Blue Gene/Q Supercomputer System at the Jlich Supercomputing Centre,” J. of Large-ScaleRes. Facil. A1 , 1 (2015).[52] F. Bloch, “Nuclear induction,” Phys. Rev. , 460–474 (1946).[53] M. Suzuki, “Decomposition formulas of exponential operators and lie exponentials with some applications to quantum mechanics andstatistical physics,” J. Math. Phys. , 601 – 612 (1985). | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt
Numerical Recipes (Cambridge University Press, Cambridge, 2003). [50] V. V. Dobrovitski, H. De Raedt, M.I. Katsnelson, and B.N. Harmon, “Quantum oscillations without quantum coherence,” Phys. Rev.Lett. , 210401 (2003).[51] M. Stephan and J. Docter, “JUQUEEN: IBM Blue Gene/Q Supercomputer System at the Jlich Supercomputing Centre,” J. of Large-ScaleRes. Facil. A1 , 1 (2015).[52] F. Bloch, “Nuclear induction,” Phys. Rev. , 460–474 (1946).[53] M. Suzuki, “Decomposition formulas of exponential operators and lie exponentials with some applications to quantum mechanics andstatistical physics,” J. Math. Phys. , 601 – 612 (1985). | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt
Numerical Recipes (Cambridge University Press, Cambridge, 2003). [50] V. V. Dobrovitski, H. De Raedt, M.I. Katsnelson, and B.N. Harmon, “Quantum oscillations without quantum coherence,” Phys. Rev.Lett. , 210401 (2003).[51] M. Stephan and J. Docter, “JUQUEEN: IBM Blue Gene/Q Supercomputer System at the Jlich Supercomputing Centre,” J. of Large-ScaleRes. Facil. A1 , 1 (2015).[52] F. Bloch, “Nuclear induction,” Phys. Rev. , 460–474 (1946).[53] M. Suzuki, “Decomposition formulas of exponential operators and lie exponentials with some applications to quantum mechanics andstatistical physics,” J. Math. Phys. , 601 – 612 (1985). | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt
Numerical Recipes (Cambridge University Press, Cambridge, 2003). [50] V. V. Dobrovitski, H. De Raedt, M.I. Katsnelson, and B.N. Harmon, “Quantum oscillations without quantum coherence,” Phys. Rev.Lett. , 210401 (2003).[51] M. Stephan and J. Docter, “JUQUEEN: IBM Blue Gene/Q Supercomputer System at the Jlich Supercomputing Centre,” J. of Large-ScaleRes. Facil. A1 , 1 (2015).[52] F. Bloch, “Nuclear induction,” Phys. Rev. , 460–474 (1946).[53] M. Suzuki, “Decomposition formulas of exponential operators and lie exponentials with some applications to quantum mechanics andstatistical physics,” J. Math. Phys. , 601 – 612 (1985). | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt
Numerical Recipes (Cambridge University Press, Cambridge, 2003). [50] V. V. Dobrovitski, H. De Raedt, M.I. Katsnelson, and B.N. Harmon, “Quantum oscillations without quantum coherence,” Phys. Rev.Lett. , 210401 (2003).[51] M. Stephan and J. Docter, “JUQUEEN: IBM Blue Gene/Q Supercomputer System at the Jlich Supercomputing Centre,” J. of Large-ScaleRes. Facil. A1 , 1 (2015).[52] F. Bloch, “Nuclear induction,” Phys. Rev. , 460–474 (1946).[53] M. Suzuki, “Decomposition formulas of exponential operators and lie exponentials with some applications to quantum mechanics andstatistical physics,” J. Math. Phys. , 601 – 612 (1985). | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(2,2)||C(3,3)| | C ( .,. ) | tt |C(1,1)||C(1,2)||C(1,3)||C(2,1)||C(2,2)||C(2,3)||C(3,1)||C(3,2)||C(3,3)| | C ( .,. ) | tt