From the open Heisenberg model to the Landau-Lifshitz equation

Magnetic systems can be described by the classical Landau-Lifshitz (LL) equation or the fully quantum open Heisenberg model. Using the Lindblad master equation and the mean-field approximation, we demonstrate that the open Heisenberg model is reduced to a generalized LL equation. The open dynamic is modeled using spin-boson interactions with a common bosonic reservoir at thermal equilibrium. By tracing out the bosonic degrees of freedom, we obtain two different decoherence mechanisms: on-site dissipation and an effective spin-spin interaction mediated by bosons. Using our approach, we perform hysteresis calculations, closely connected with the Stoner-Wohlfarth theory. We compare the exact numerical master equation and the mean-field model, revealing the role of correlations originated by non-local interactions. Our work opens new horizons for the study of the LL dynamics from an open quantum formalism.


I. INTRODUCTION
Since its discovery in 1928 [1], the original Heisenberg exchange interaction between two spins JS 1 · S 2 has been extended to complex magnetic arrangements and successfully implemented in a variety of quantum systems. Nowadays, myriads of physical models describing the interaction between N spin-1/2 particles are based on particular cases of the Heisenberg Hamiltonian where α, β = x, y, z are the components of the spin operators. A suitable choice of the magnetic field B, the exchange coupling constants V ij αβ , and the topology of the system allows to understand the origin of magnetic ordering [2], phase transition [3], spin-wave excitations [4], lattice effects [5], to name a few. Furthermore, the Heisenberg dynamic is currently reproduced in different physical systems such as circuit quantum electrodynamics [6], cavity QED [7], superconducting devices [8], one-dimensional interacting spins [9], Rydberg atoms [10][11][12][13], and trapped ions [14,15]. Many of the previous setups deal with unavoidable relaxation processes induced by system-bath interactions, which is well understood in terms of the Lindblad master equation [16,17]. We use the term open Heisenberg model (OHM) to describe any system of N interacting spins with a Hamiltonian structure similar to Eq. (1) and subject to an interaction with an external bath [8,9,18].
From the classical point of view, the dynamics of spins can be described by the Landau-Lifshitz (LL) equation [19] 1 γ where M is the magnetization of the system and B eff is an effective magnetic field which includes internal and external contributions to the magnetization dynamics of the system, such as the anisotropy and the applied field.
Because of the universality of the Heisenberg model to describe magnetic properties, the following question arises. Can the OHM reproduce magnetization dynamics similar to the classical LL equation? The answer is yes, and more importantly, we shall demonstrate that the dynamics is more general under certain conditions. Here, we establish an unexplored connection between the LL equation and a particular Lindblad superoperator in a Markovian master equation. It is worth noticing that other quantum approaches have successfully addressed the microscopic derivation of the LL equation by considering quantum processes. For instance, using the spin-wave theory [44], the Fokker-Planck equation [40], the Dirac-Kohn-Sham theory [45,46], a mean-field tight-binding model [47], a non-Hermitian Hamiltonian approach [48], and the Yang-Mills equation [49]. Nevertheless, our approach is a new perspective that is useful for connecting the LL dynamics with the evolution of open quantum systems in the meanfield regime. The latter is particularly advantageous to future simulations of magnetic-like phenomena using quantum systems like trapped ions [14,15], superconducting devices [8] and cavity QED [7]. Moreover, it allows the study arXiv:2006.15658v1 [quant-ph] 28 Jun 2020 of more general environments exhibiting memory effects usually described as non-Markovian master equations [50].
The paper is organized as follows. In Sec. II, we introduce the Hamiltonian of the system and the Markovian master equation. Section III introduces the mean-field approximation leading to a generalized version of the LL equation for the magnetization dynamics. Here, we discuss the hysteresis of the system in analogy with the Stoner-Wohlfarth theory [51]. Finally, in Sec. IV we numerically solve the quantum master equation to compare our results with the mean-field model. Finally, we discuss the effect of correlations by considering both closed and open dynamics for the isotropic and anisotropic Heisenberg models.

II. MODEL
We consider a linear spin chain composed of N spin-1/2 particles with on-site and hopping interaction terms described by the following system Hamiltonian where γ is the electronic gyromagnetic ratio and V α is the coupling between adjacent spins along the cartesian directions e α . Here, S (j) = ( /2)σ (j) is the spin operator of the j-th particle, where σ (j) α are the Pauli matrices for S = 1/2. Note that for B = B x e x and V y = V z = 0, Eq. (3) reduces to the standard transverse-field Ising model [9,12,64].
Other relevant magnetic coupling terms such as the dipole-dipole (DD) or the Dzyaloshinskii-Moriya (DM) [52,53] are crucial to model magnetic defects in solid-state systems [56] or skyrmions [54,55], respectively. However, in this work, we focus either on small systems of individual spins, or large systems of homogeneous magnetic moments. In the former, neglecting DD interaction is valid as it is several orders of magnitude smaller than the exchange coupling, and becomes relevant only as a long-range interaction [57,58]. Moreover, for large homogeneous systems, the DD interaction can be considered as an additional anisotropy [59,60], and thus not calculated explicitly. Similarly, the DM interaction is typically present in magnetic systems with broken inversion symmetry [61][62][63], which lies outside of the scope of this paper.
In general, quantum systems interact with its surrounding environment, which originates from different damping mechanisms. In our model, we assume that the spin chain is coupled to a generic boson reservoir that is in thermal equilibrium. In order to introduce dissipative effects, we focus on Markovian evolutions [16], i.e., without memory effects. Non-Markovian dynamics have been analytically solved for some particular spin systems [65,66]. However, we are interested in the dynamical properties of the system at time scales larger than the characteristic bath correlation time, i.e., where the first and second Markov approx-imations hold [50]. Hence, dissipation originates from the following system-bath interaction Hamiltonian where k runs over infinite bosonic modes of the environment and g αjk are the spin-bath coupling constants. The environment is considered as a collection of harmonic oscillators described by the bath Hamiltonian H b = k ω k a † k a k , where a k and a † k are the annihilation and creation boson operators, respectively. Spin operators S (j) α,± are defined as with S   . We remark that the interaction Hamiltonian (4) assumes that all spins couple to the same environment. As a consequence, spins will experience an effective coupling to each other in the Markovian master equation [67] [In our calculations it will appear as the last term in Eq. (9)]. The interaction Hamiltonian (4) can be written as where d k = g yjk − g xjk , e k = g xjk + ig yjk + g zjk , c k = −g xjk + ig yjk , and S (j) y . Eq. (8) reveals that the interaction Hamiltonian (4) is composed by puredephasing (d k ), amplitude damping or energy-exchange (e k ) and counter-rotating terms (c k ). For g xjk = g yjk = 0 the interaction Hamiltonian reduces to the standard amplitude damping model [68], where bosons are absorbed or emitted inducing spin flip-flop processes between the states |↓  z . The full Hamiltonian H = H s + H b + V is similar to the chain-boson model [69], however our interaction is more general since we are including pure-depashing and counter-rotating terms.
In the Markov and secular approximations, we derive the following master equation that considers effective spin interactions up to first nearest-neighbors is coupled to its nearest-neighbour via the coupling constants Vα. The common boson bath allows the on-site decay channels γα,± and the dissipative hopping terms gα,±.
(b) Diagram showing the physical interpretation of the Lindbladian Lnn(ρ) (9). Adjacent spins in the sites (j + 1) and (j) exchange energy with the boson bath and tracing out over the bath degrees of freedom we obtain the effective dissipation gα,±.
where γ α,+ and γ α,− are damping rates associated to absorption and emission processes between each site and the external boson environment, respectively. Here, j, j denotes nearest-neighbor spins by considering all terms satisfying the condition |j − j | = 1. The on-site Lindbladian L os (ρ) has already been implemented using Rydberg atoms at low temperatures [10,11], with γ z,− > 0 and γ x,η = γ y,η = γ z,+ = 0. On the other hand, we called the superoperator L nn (ρ) as the nearest-neighbours Lindbladian since it accounts for the effective energy-exchange between adjacent spins. This dissipation's source naturally appears in the master equation of multi-atomic systems coupled to light [70], and it is a pivotal result towards connecting the OHM with the LL theory, since L nn (ρ) will reproduce the damping term M × (M × B eff ) in equation (2). Further details of the microscopic derivation of the master equation (9) is presented in Appendix A, and it follows the spirit of the open dynamics for interacting qubits presented in Ref. [68]. A representation of the OHM is depicted in Fig. 1.
In the next section, we use the mean-field approximation to solve the OHM given in Eq. (9), and after introducing the non-linear single-particle dynamics, we show its connection with the LL Eq. (2).

III. MEAN-FIELD APPROXIMATION AND LANDAU-LIFSHITZ EQUATION
The mean-field approximation considers the many-body density matrix of the N -particle system as a separable ten-sor product of single-particle density matrices ρ j (t), such that The above factorization is more accurate as N increases [71] and it has been also discussed in the context of open quantum systems [16,72] or using the name of Hartree approximation [73]. Each density operator in Eq. (10) is described as a two-level system, where ρ j (t) = To shed more light on the magnetization dynamics of the system we introduce the magnetic moment of each particle through the relation µ (j) = −g s µ B S (j) / , where µ B = 9.27 × 10 −24 JT −1 is the Bohr magneton, and g s ≈ 2 is the g-factor. After averaging the effect of all magnetic moments the macroscopic magnetization reads where N is the number of spins, V = L 3 is a characteristic volume of the system, L is a characteristic length, and σ In what follows we focus on the particular case γ x,± = γ y,± = 0 in order to reduce the number of damping rates in our analysis. Thus, we can explicitly calculate the single-particle dynamics by solving dρ 2 /dt = Tr 1,3,...,N [L(ρ MF )], where Tr 1,3,...,N is the partial trace over the remaining N − 1 particles (without considering j = 2). After applying the partial trace, we obtain the following set of coupled non-linear equations: (14) where m α = M α /|M|, M = M x e x + M y e y + M z e z , g α = g α,− − g α,+ , and Γ = γ z,− + γ z,+ is the total damping when γ x,± = γ y,± = 0. Because of the boundary conditions of the linear chain, the first (j = 1) and last (j = N ) particles interact with a single neighbor spin, in contrast to the intermediate spins (j = 1, N ) that interact with two nearest-neighbors. Therefore, for particles j = 1, N Eqs. (12)-(14) must be modified by considering (14) reduces to the transverse-field open Ising model presented in Ref. [9]. As expected, the singleparticle dynamics is affected by the presence of an effective field induced by the other N − 1 particles. In fact, the first two terms on the right-hand side of Eqs. (12)-(14) are recognized as the effective magnetic field whose components are B eff,α = γB α +m α V α . The third term on the right-hand side of Eqs. (12)-(14) is crucial for the LL theory and for illustration, we write it in a compact form, By a direct calculation we get the constraint which plays an important role on the dynamics, since it leaves invariant the magnitude of the magnetization vector M. To understand this, we write the magnetization dynamics induced by L α asṀ α = L α /2|M|. We note that Additionally, the last term in Eqs. (12)-(14) accounts for on-site dissipation induced by the boson bath that directly affects the magnitude of the magnetization vector. Based on these observations, and considering a large number of spins, we obtain the following dynamical equation for the macroscopic magnetization vector where B eff = B + B an is the effective magnetic field responsible for the gyromagnetic precession of the magnetization vector. In our model, B an = m x V x e x + m y V y e y + m z V z e z is the anisotropy field caused by the local interaction between spins [40] [see Heisenberg Hamiltonian (3)]. On the other hand, we recognize D = (1/2γ)(g x e x + g y e y + g z e z ) as the magnetic field responsible for modifying the precession of the magnetization similar to LL equation (2). The last two terms in Eq. (18) are given by where R is the relaxation tensor [74] and R 0 is the noise induced by the boson environment. These terms are in agreement with the Bloch theory applied to magnetic systems [41].
In order to illustrate the scope of Eq. (18), we derive some well-known models as particular cases. First, for D = αB eff and Γ = 0, Eq. (18) exactly reduces to the LL equation (2), where α is the dimensionless damping factor [19,20]. Second, for D = 0, our model is reduced to which is known as the Bloch-Bloembergen equation [41,42], with M SS being the stationary state. Finally, for R 0 = 0 Eq. (18) reduces to the Callen's equation [75,76], which is a phenomenological equation used to describe dissipative spin systems. Hence, the generalized LL equation (18) for the magnetization vector is capable of reproducing different models. Also, it is closely connected to a Markovian master equation in the mean-field approximation and stands as one of the most relevant results in this work. Moreover, our microscopic model shows that a suitable choice of the common boson bath (10) can generate a particular Lindblad operator L nn (ρ) (9) which is intimately related to the damping of the precession of the magnetization vector. In the next section, we will explore the magnetic properties of the system by analyzing different hysteresis curves.

A. Hysteresis curves
Systems with hysteresis are relevant in nature because they can be experimentally manipulated to understand their response under an external force or action [77]. In magnetic materials, the hysteresis curve is the relation between the steady state (SS) of the magnetization as a function of the external magnetic field applied along an arbitrary direction. As a first step, we neglect the on-site dissipation (Γ = 0) and thus, the magnitude of M is constant during the dynamics. Moreover, we set D = αB eff with α < 1 in order to model the LL dynamics. The SS of the magnetization, dM SS /dt = 0, can be obtained by solving M SS × B eff + αM SS × (M SS × B eff ) = 0, and the non-trivial solution imposes that M SS must be parallel to the effective magnetic field. Under these assumptions, the dynamics is described by the LL equation (2), where the effective magnetic field is B eff = B + B an .
We numerically solve the time evolution of the system for the parametrized initial condition M(0) = |M(0)|(cos(φ 0 ) sin(θ 0 ), sin(φ 0 ) sin(θ 0 ), cos(θ 0 )). For simplicity, we choose the initial angles as φ 0 = 0 and θ 0 = π/40 to simulate a magnetic system slightly misaligned respect to the z-axis. In what follows, we introduce our natural units by setting γ = 1 and |M(0)| = 1. As a consequence, the magnetization components satisfy |M α (t)| ≤ 1 for t ≥ 0. In Fig. 2(a) we show the time evolution M α (t) for an effective magnetic field B eff = (B z + m z V z )e z , with B z = −2, V z = 0.5, and considering N = 500 spins. One can observe that M undergoes a dissipative precession leading to M SS = (0, 0, −1). As initially M z (0) ≈ 1, the anisotropy field at t = 0, B an (0) = m z (0)V z e z , points in the e z direction. Then, in the presence of a negatively increasing magnetic field B z < 0, the longitudinal component B eff,z = B z + m z V z becomes negative below the critical magnetic field B crit z = −m z V z , meaning that now B eff points in the −e z direction. Thus, as the magnetization follows B eff in order to reach the SS, all the spins suddenly rotate and the system ends in the final state M z (∞) = −1. Evidently, for a larger anisotropy field V z is necessary a larger negative component of the external magnetic field B z to generate this collective rotational effect. The situation is reversed when M z (0) ≈ −1, i.e B an (0) points in the −e z direction, and thus a positive external field B z > m z V z is necessary to induce the rotation towards the stationary state M z (∞) = 1. These observations explain the hysteresis curves illustrated in Fig. 2(b). Now, we investigate the effect of including an additional perpendicular field B ⊥ = B x e x + B y e y . We choose B x = B y = 1 but preserving V x = V y = 0. In Fig. 2(c)  the perpendicular field results in a SS with perpendicular components, i.e. M SS = (0.31, 0.31, −0.89). As a consequence, the hysteresis curves in Fig. 2(d) have a smooth dependence in terms of the external field B z , and in some cases (V = 0.5 and V = 1), the magnetic coercivity (width of the hysteresis curve) is zero because the anisotropy is weaker than the in-plane applied field (B 2 x + B 2 y ) 1/2 = √ 2. Note that the steep transitions observed in Fig. 2(b) are due to a field that is collinear with the anisotropy. When strong x, y components of the magnetic field are present, the anisotropy becomes less relevant within B eff and M follows more readily the direction of B.
Our hysteresis curves are in agreement with those studied using the Stoner-Wohlfarth theory [51,78,79], which minimizes the energy of a magnetic system E = −B·M−V z M 2 z using either a single domain description or a mean-field approximation. In the next section, we solve the spin dynamics beyond the mean-field approximation by numerically solving the master equation.

IV. DENSITY MATRIX FORMALISM
In this section, we compare the mean-field model (18) with the numerical solution of the master equation (9). Let's begin with a different initial condition for the normalized magnetization, say M(0) = (1, 0, 0) or equivalently all spins in the state ρ where |↑ z and |↓ z are the eigenstates of σ z . The magnetization components M α = (1/N ) j σ (j) α are immediately computed from the density matrix ρ(t) using the relation σ α ρ]. A detailed numerical method to solve the master equation is given in the Appendix (B), where the implementation of the Lindblad superoperator L on (ρ) is presented. In Fig. 3(a),(b) we show the time evolution of the magnetization components M x (t) and M z (t) for an effective magnetic field B eff with B x,y = 0.25, B z = −0.5, and V x /5 = V y = V z = 0.1. We observe that the mean-field model partially recover the dynamics of the exact density matrix approach, exhibiting a good agreement at shorter times.
The main mismatch occurs at longer times, i.e., the prediction of the stationary state of the system. We remark that the assumption of the mean-field approximation is the product decomposition given in Eq. (10). However, it is expected that local interactions between spins governed by V α (3) and the nearest-neighbor Lindbladian L nn (9) could generate correlations during the dynamics, even starting from an uncorrelated many-body state. The latter is discussed in Sec. IV A, where correlations are analyzed in more detail. In Fig. 3

(c),(d) we show the stationary states M SS
x and M SS z as a function of the applied magnetic field B z for N = 3 and N = 4 spins. The mean-field model recover the non-monotonic shape of the components M SS x and M SS z . We observe that differences between the mean-field and master equation increases as the external magnetic field depart from B z = 0. This can also be understood in terms of the propagation of correlation by the dynamics itself.
For instance, if correlations are created at time τ > 0 the density matrix immediately takes the form ρ(τ ) = ρ 1 (τ ) ⊗ ... ⊗ ρ N (τ ) + ρ corr (τ ), where ρ corr (τ ) accounts for such correlations. Consequently, the open dynamics will be strongly affected by the extra Liouvillian generator L(ρ corr ) = −i[H, ρ corr (τ )] + L on (ρ corr (τ )) + L nn (ρ corr (τ )). Hence, for a large magnetic field B z the propagation of correlations will be dominated by the Hamiltonian contribution on the Lindblad superoperator, i.e the term −i[B z N j=1 S (j) z , ρ corr (τ )]. As a consequence, the stationary state will be affected by these additional corrections neglected by the mean-field model, explaining the differences observed in Fig. 3(c),(d). In the next subsection, we discuss in more detail the effect of the spin correlations on the dynamics by considering both closed and open dynamics.

A. Spin correlations
In order to explain the mismatch between the mean-field approach and master equation, we remark that the meanfield approximation considers that the density matrix of the system can be written as a separable tensor product, as given in Eq. (10). Therefore, we state that whenever the system departs from this representation, i.e. correlations between the spins show up, the mean-field theory will be deteriorate. To quantify spins correlations, we introduce the two-point correlation function [80], where {i, j} = 1, ..., N are the spin indexes, σ (i) α is the i-th Pauli operator. Typically, the expectation values σ (i) α are calculated by assuming the system at thermal equilibrium. However, in our microscopic model we are interested in the non-equilibrium properties of the system. Therefore, we calculate the expectation values as σ with ρ j (t) = Tr 1,...,N =j [ρ(t)] with ρ(t) being the solution of the master equation (9). For further comparison, we denote M MF as the mean-field solution of the magnetization vector which is obtained by solving Eqs. (12)- (14). In parallel, we compute M Exact from the exact master equation (9) using the definition given in Eq. (11). Differences between both approaches will quantified in terms of |M MF − M Exact | 2 . In the next subsections we analyse two different scenarios, Case I: Γ = g α = 0 (closed dynamics) and Case II: Γ, g α > 0 (Markovian open dynamics) 1. Case I: Γ = gα = 0 For a system described by a closed dynamics, (Γ = g α = 0), the time evolution of the spin chain is governed by the Heisenberg Hamiltonian (3). To better understand the role of correlations, it is instructive to distinguish between the isotropic and anisotropic Heisenberg models. A fully isotropic Heisenberg model occurs when V x = V y = V z = V . Now, we shall illustrate that isotropic interactions are connected with hidden Hamiltonian symmetries. Thus, we define the isotropic interaction Hamiltonian which is the second term of Hamiltonian (3) for V α = V . To take into account conserved quantities, we introduce the total spin operator for each component We note that the above operator is related to the components of the magnetization vector (11) through the relation M α = 2µ S tot α /(V N ). We get [H iso , S tot α ] = 0 for each component α = x, y, z which is known as the SU (2) symmetry. Hence, the total spin operator is a conserved quantity under the action of H iso , which means that the Hilbert space of the N particles separates into disjunct Hilbert subspaces with constant magnetization. As the Zeeman contribution on the Heisenberg Hamiltonian (3) is a local effect, we conclude that spin correlations cannot be generated for the isotropic case if the system does not interact with an external environment. These observations implies that the mean-field decomposition (10) is valid for the isotropic case if Γ = g α = 0. As a consequence, the two-point correlation functions C ij αβ = 0 for all i, j, α, β, revealing that the meanfield and exact models are identical in this particular case (which is also numerically corroborated). In other words, in the absence of a reservoir, the only non-local effect comes from the coupling V α S (j) α S (j+1) α , and spin correlations only arise when we depart from the ideal isotropic case V α = V .
Let us consider a spin chain with anisotropy such that V x = V + ∆ and V y = V z = V , where ∆ and V are the anisotropy shift and the isotropic contribution, respectively. For ∆ = 0, we recover the previous case where non-local interactions leaves invariant the magnetization vector. For ∆ > 0 the anisotropy is present along the x axis, while for ∆ < 0 the anisotropy is changed to the y, z plane. In addition, we use the same initial condition detailed in Sec. IV, where M(0) = (1, 0, 0) and we fix V = 0.1. In Fig. 4 we plot the time average of |M MF − M Exact | 2 for a time interval (0, 10 3 ) by considering B x = B y = 0.25, B z = −0.5, and N = 3. We observe that |M MF − M Exact | 2 monotonically increases with the anisotropy shift. This means that configurations without SU (2) symmetry (∆ = 0) can generate correlations between spins and M SS = M Exact in such cases. In the next subsection, we analyse the effect of including the losses induced by the bosonic thermal environment.

Case II: Γ, gα > 0
For the Markovian open dynamics, (Γ, g α > 0), the evolution of the system is ruled by the Lindblad master equation (9). For an isotropic Heisenberg Hamiltonian, the only source of non-locality is given by the nearest-neighbours Lindbladian L nn , which is by definition a non-local superoperator acting on adjacent spins. In Fig. 5(a) we plot the function |M MF − M Exact | 2 for the isotropic case by considering that all spin are aligned in the x direction at t = 0. We define γ z,+ = γ 0 n b and γ z,− = γ 0 [n b + 1] with γ 0 = Γ/(2n b + 1) and n b the mean number of phonons such that Γ = γ z,− + γ z,+ for L on and g z,η = γ z,η /10 for L nn . Along this work we set n b = 0.08 for our simulations.
At the beginning of the dynamics, the mean-field and master equation predict the same magnetization vector, see Fig. 5(a). However, after a critical time (which is shorter for larger values of Γ), we observe a deviation between both magnetization vectors, leading to a constant stationary difference at longer times. For comparison, in Fig. 5(b) we show the time evolution of the two-point correlation func-tion C 12 xx (t). For simplicity we only show C 12 xx , however, the same behaviour is numerically obtained for other components, i.e C 12 αβ with α, β = x, y, z. At shorter times we observe that C 12 xx (t) = 0 which means that the density matrix is given by the product state (10). This result supports the good agreement between the mean-field model and the master equation since the mean-field assumption is fully satisfied. As time increases, the system can no longer be described as a product state leading to C 12 xx (t) = 0, and the mean-field approximation fails. We remark that in the limit Γ → 0 the system converges to the isotropic closed dynamics, i.e M MF = M Exact for all times.
For an anisotropic Heisenberg Hamiltonian, we have two sources of non-locality, the interaction Hamiltonian (dominant) and the Lindbladian L nn (small contribution). Now, we simulate the previous case but adding anisotropy in the x direction, i.e. V y = V z = 0.1 and V x = {0.5, 1}. In Fig. 5(c), we observe that the addition of the anisotropy increases the mismatch between the mean-field and exact calculations. Furthermore, as one would expect, the coupling constant V x (local interactions between spins) contributes to increasing correlations, which is illustrated in Fig. 5(d). Moreover, we observe that the mean-field model can not exactly predict the steady-state of the system because of the correlations originated by non-local terms, which confirms our previous observations in Fig. 3.

Quantum correlations
To complement the previous analysis based on correlations we note that one remaining question is whether the observed correlations are quantum in nature or just a statistical mixing of the density matrix, i.e. a mixed state. Quantum correlations (QC) are one of the most fundamental concepts in Quantum Information Theory [81]. Even more, QC provides a useful resource to speed up several tasks in quantum computing [82]. In particular, the concurrence defined by Wootters [83] in the context of a twoqubit system, is a widely used measure that account for QC based on the separability of the system. For a system with N spins we trace over N − 2 spins and we obtain the two-qubit density matrix ρ ij (t) = Tr 1,...,N =i,j [ρ(t)], where particles i and j are two different arbitrary particles of the spin chain. Then, we calculate the spin-flipped density ma- y ), and the Concurrence is given by where the α 1 , α 2 , α 3 , α 4 are the square root of the eigenvalues of ρ ijρij in decreasing order.
The concurrence is a bounded function, 0 ≤ C ij ≤ 1, where C ij = 0 means zero entanglement between particles i and j. In Fig. 6 we plot the concurrence C 12 for N = 3 by considering the open anisotropy case illustrated in Fig. 5(c),(d). First, we observe that QC reaches a nonnegligible maximum value, i.e max[C 12 ] ≈ 0.26 for V x = 1. However, for the parameters used in the simulation, we note that QC appears in a short time window, which is related to the temporal region where both master equation and mean-field models are different. Second, QC increases with the local coupling term V x , which is expected since the interaction between adjacent spins creates bi-partite entangled states. More generally, one could find QC by performing a local measurement on the spins instead of tracing out them. This allows us to define a localizable entanglement [84], which has been shown that is closely related to the two-point correlation function [84]. Nevertheless, local measurements will involve additional resources that we are not considering here, and thus, it is out of our scope.

V. CONCLUSIONS
In summary, we established an unexplored microscopic connection between the open Heisenberg model and the Landau-Lifshitz equation. Starting from a generic spinboson interaction Hamiltonian, we derived a Markovian master equation, and applying the mean-field approximation, we found a generalized LL equation. Consequently, we recognized the microscopic origin of anisotropy effects, dissipative magnetic fields, and relaxation processes induced by the Heisenberg Hamiltonian and external boson bath.
First, we focused on the hysteresis curves for longitudinal and transverse magnetic fields, reaching a good agreement with the Stoner-Wohlfarth theory. We also solved the non-equilibrium dynamics by numerical calculations of the master equation for a small number of spins. We compared the mean-field and master equation, revealing that a mean-field model phenomenologically describes the main magnetic features such as temporal behavior (oscillations and decay) and stationary states as a function of external fields, although with some deviations on their exact behavior.
Using the two-point correlation function and the concurrence, we showed that these deviations are due to a correlation originated from magnetic anisotropy and the non-local Lindbladian between spins, making the central assumption of the mean-field approximation invalid. We expect these deviations to be negligible as the number of spins increases. Finally, our model can be used to connect the dynamics of open quantum systems with magnetic-like systems.  −βH b )) is the partition function, β = (k B T ) −1 is the inverse temperature, and H b = k ω k a † k a k is the bath Hamiltonian (harmonic oscillators). In addition, we employ the born approximation [50], where is assumed that at any time the full density matrix can be decomposed as an uncorrelated product state, i.e.ρ s+b (t) =ρ s (t) ⊗ ρ b . The latter is valid in the weak-coupling limit, which is fulfilled when g αjk max(γ|B|, V α ). Under these assumptions, we have Tr b (a k ρ b ) = Tr b (a † k ρ b ) = 0, and therefore the first term of the right-hand of Eq. (A5) vanishes. As a consequence, we derive the following convolution dynamics where A k αj;α j = g * αjk g α j k n(ω k ), B k αj;α j = g αjk g * α j k [n(ω k ) + 1] are the coupling terms associated with absorption and emission processes, respectively. Here, n(ω k ) = [exp( ω k /k B T ) − 1] −1 is the mean number of bosons at thermal equilibrium. Now, we applied the first and second Markov approximations [50] and we assume that ρ s (t ) ≈ρ s (t) and that the integral contribution can be evaluated at larger times, i.e. for t → ∞. Now, we introduce the spectral decomposition [16,85] where δ(ω ba − ω) is a Kronecker function, i.e. δ(x) = 1 for x = 0, and δ(x) = 0 otherwise. The quantum states |a , |b are eigenstates of the Heisenberg Hamiltonian (3), with V ij αβ = δ αβ δ i,j−1 V α and α = β. In the interaction picture, the following relations are satisfied in the frequency domain α ,+ (t ) into Eq. (A6) using the spectral decomposition (A9) and (A10) we obtain the oscillating functions exp(±i(ω − ω)t). In the secular approximation, we neglect the terms ω = ω due to the condition τ b T s , where τ b ∼ g −1 αjk and T s ∼ 1/ max(γ|B|, V α ) ∼ T s are the bath and system characteristic times, respectively. Therefore, in the secular and Markov approximations, we obtain the following Lindblad master equation in the Schrödinger picture Finally, we make the last approximations to derive the phenomenological master equation presented in Eq. (9). First, we consider a nearest-neighbor interaction model to transfer energy between adjacent spins, which implies that we only consider contributions satisfying the condition |j − j | = 1. Second, we assume that the anisotropy induced by the common reservoir has the same form as the Heisenberg model presented in Sec. II, then α = α . Finally, following the Einstein model's spirit for the heat capacity in solid-state physics, we introduce a phenomenological average resonant frequency ω ∼ ω 0 for all spins and therefore ω = ω 0 . Under these assumptions the Lindblad master equation (A11) reduces to Eq. (9).

Appendix B: Solving the master equation
To numerically solve the master equation we adopt the following general solution [9,86] where R k and L k are the right and left eigenmatrices given by the equations L(R k ) = λ k R k and L † (L k ) = λ k L k , respectively. The Lindblad generator L is defined from the structure of the Markovian master equationρ = L(ρ). The matrices R k and L k must to satisfy the orthonormality condition Tr(R k L k ) = δ kk , c k = Tr(ρ(0)L k ) are coefficients with ρ(0) being the initial state, and λ k the corresponding eigenvalues of the right eigenmatrices R k . For numerical purposes is convenient to sort the eigenvalues λ k = λ R k +iλ I k by choosing 0 = λ R 1 ≤ λ R 2 ... ≤ λ N d 2 , where N d = 2 2N is the number of eigenvalues of the system. The zero eigenvalue λ 1 = 0 is related to the stationary state since others eigenvalues with k > 1 satisfy λ R k < 0 [17] leading to dissipative terms ∝ e −|λ R k |t t→∞ −−−→ 0. To compute the matrices R k and L k we employ the formalism presented in Ref. [87], where the strategy is to rewrite the effect of the Lindblad generator L on a more involved vector space. To this end, the many-body density matrix is mapped to a new vector space as follow where |k, l = |k ⊗|l is the new vector basis constructed by the initial vector basis |k and C = ( k,l |ρ kl | 2 ) 1/2 is a normalization factor. In this new vector space spawned by the basis |k, l the master equation can be rewritten as [87] d|ρ dt =L|ρ The operatorL is a complex matrix with 2 2N × 2 2N elements, and is the dissipative term of the Markovian master equation (boson reservoir). The algorithm to solve the open dynamics is quite simple. First one compute the eigenvalues and eigenvectors ofL andL † (in the new basis), then using the map (B2) we rewrite the right and left eigenmatrices in the initial Hilbert space, and finally, we employ the general solution given in Eq. (B1).