Many-particle Sudarshan-Lindblad equation: mean-field approximation, nonlinearity and dissipation in a spin system

A system of $N$ spin-1/2 particles interacting with a thermal reservoir is used as a pedagogical example for advanced undergraduate and graduate students. We introduce and illustrate some methods, approximations, and phenomena related to dissipation and nonlinearity in many-particle physics. We start our analysis from the dynamical Sudarshan-Lindblad quantum master equation for the density operator of a system $\mathcal{S}$ interacting with a thermal reservoir $\mathcal{R}$. We derive the quantum version of the so-called Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) equations such that the master equation can be decomposed in a hierarchical set of $N-1$ equations ($N>1$). The hierarchy is broken by introducing the mean-field approximation and reducing the problem to a nonlinear single particle system. In this scenario, the Hamiltonian is nonlinear (i.e., it depends on the state of $\mathcal{S}$), although the superoperator responsible for the dissipation and decoherence of $\mathcal{S}$ remains unaffected. To provide a useful tool to students: (1) we discuss the physical approximations involved, (2) we derive the analytical solution to the mean values equations of motion resulting from the Hamiltonian, (3) we solve analytically the master equation in the stationary regime, (4) we obtain and discuss the solution of the nonlinear master equation, numerically, and finally, (5) we discuss the master equation beyond the mean-field approximation and show how to introduce higher order quantum correlations that have been previously neglected.


I. INTRODUCTION
Since the birth of quantum mechanics (QM) the investigation of irreversible processes has received special emphasis, with the aim of describing a physical system evolving in compliance with the human feeling that there exists an "arrow of time" in nature. The approaches had to go along with the observed physical systems where macroscopic phenomena evolve inexorably to an equilibrium or stationary state, so breaking the time symmetric microreversibility present in the dynamical Schrödinger equation 1 . Since its incipience several thinkers have devoted efforts to describe the irreversible processes with the care to avoid the violation of the QM principles. For a few historical reading see the papers [1][2][3][4][5][6][7][8][9][10][11][12], and for a review we recommend the books [13][14][15][16].
The standard rationale to describe formally the time irreversibility of a system of interest S is to couple it to another one, R, chosen according to the phenomenon one wants to observe -such a system is sometimes called the rest of the universe, the environment or the reservoir -and make them to interact via a convenient potential. Thence, as times go on, both systems will evolve and eventually intertwine their degrees of freedom and will share a common quantum state. However, since one does not have access (either by observation or through experiment) to the details of R, methodologically one proceeds as one commonly does in statistical theory: one gets the reduced probability function for S by integrating (or summing) over the variables representing the rest of the universe R. In quantum mechanical language one says that we calculate the trace over the system R degrees of freedom. So one remains with a dynamical equation, called master equation, for the state of the system S, represented by the density operator (or density matrix), which shall evolve irreversibly in time as originally sought. The master equation for S will contain the influence of R through the presence of extra terms depending on its parameters. This approach received a thorough attention essentially in the 1960's, as can be seen in the researches reported in Refs. [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31][32], and for some more recent applications see Refs. [33][34][35][36][37][38]. In addition, in classical statistical mechanics [39] and in the stochastic processes literature [40], a master equation is a gain-loss equation, in terms of probability transitions, for the probabilities of the possible states of a system, while in quantum mechanics, the master equation governs the density operator evolution in time [33].
The standard procedure to describe irreversibility in S, starting from equations that are time-reversible for the complete system S + R, consists in deducing a differential equation for the S state only, by calculating the trace over the degrees of freedom or washing out the detailed information about R, which leaves us with a new new modified (time-irreversible) equation for S containing the parameters that were already present in R and also news one as the absolute temperature. This approach leads to an essential equation known as Sudarshan-Lindblad (SL) equation [41][42][43], that describes irreversibility, decoherence, and dissipation in S. In addition, this approach has been applied to several recent papers related to the quantum measurement theory [44][45][46][47][48]. A recent paper was devoted to the derivation of that equation with the pedagogical purpose of being accessible for undergraduate students in physics, chemistry and mathematics, having a basic background on QM [49].
In this paper we have the concern to be quite pedagogical in order to assist advanced undergraduate and graduate students interested in many-particle physics, once it involves irreversibility, dissipation and nonlinearity for a manyparticle system. We show that for a system S constituted by N free spin-1/2 particles (or, by the mathematical isomorphism, a system of two-level atoms) interacting with a thermal reservoir R, and adopting the mean-field approximation, the single-particle Sudarshan-Lindblad equation remains the same, although the Hamiltonian acquires, additionally, nonlinear terms, i.e., they depend on the state of the system S, and the dynamical equations for the spin operators mean values become nonlinear. We obtain and digress on the solution of the nonlinear equations of the single-particle mean values operators, and how to go beyond the mean-field approximation by introducing two-particle quantum correlations.
The article is organized as follows: In Sec. II we introduce the Hamiltonian and the master equation describing the system of N spin-1/2 particles interacting with a thermal reservoir. In Sec. III, a mean-field approximation is developed and an effective single-particle Hamiltonian is obtained. In Sec. IV, we derive the analytical solution for the nonlinear mean-values operators, whose motion is generated by the nonlinear Hamiltonian without dissipative terms, and illustrate numerically the case with dissipation. In Sec V, we discuss how to go beyond the mean-field approximation by introducing two-particle correlations. Finally, in Sec. VI we present a summary and the conclusions.

II. EVOLUTION EQUATION FOR A N SPIN-1/2 PARTICLES INTERACTING WITH A RESERVOIR
We consider dissipation in physical systems in the framework of the System-Reservoir interaction. In this framework, the full Hamiltonian describing a generic system of interest S interacting with a reservoir R in thermal equilibrium, is given by where H S is the system of interest Hamiltonian, H R is the reservoir Hamiltonian, and H SR is the system-reservoir interaction Hamiltonian. As shown in Refs [33][34][35]49], by tracing over the reservoir degrees of freedom, considering a bilinear interaction between system and reservoir, and under the Markov approximation, the system density operator ρ S has its evolution described by the Sudarshan-Lindblad master equation The superoperator L[Â] is defined for an arbitrary system operatorÂ by where Γ is the decay rate constant, andn is the mean number of the reservoir quanta. In the right hand side (RHS) of Eq. (2), the first term describes the unitary evolution of S, while the last two terms are responsible for the non-unitary irreversible and dissipative evolution of S. By particular choices of H S and the operatorÂ, the Eq. (2) has been applied to several problems related to dissipation in quantum optics and in measurement theory [33][34][35].
Here, we are interested in describing a system consisting of N spin-1/2 particles (or two-level atoms, since their description accepts the same algebra), interacting with a reservoir, and whose constituents do not interact with each other or, because their density is quite low their direct interaction can be neglected. In our case, the Hamiltonian (1) describing N particles interacting with the reservoir (an infinite number of harmonic oscillators) is given by [14] where are the collective operators which satisfy the following commutation relations and, the single-particle spin operators for the particle i are σ 0 (i) and σ ± (i); 2 ω 0 is the transition frequency between the particles' two-states. The operators b † k and b k obey the commutation relation [b k , b † k ′ ] = δ kk ′ ; their rôle is to create and annihilate a quantum of energy, ω k , of the reservoir R; ω k is the angular frequency associated to k-th harmonic oscillator belonging to R. For the system-reservoir interaction, it is assumed that each oscillator interacts with all the particles with the same strength |g k |. So, the master equation structure for the set of particles is identical to the one corresponding to the single-particle master equation [14] where the operators σ ± , σ 0 are replaced by the collective ones S ± , S 0 . Then, identifying the operatorÂ, in Eq. (2), with the collective operator S + , the Sudarshan-Lindblad N particles master equation acquires the form where is the N -particle density operator, V is the volume of the container,n = 1/ e ωo/kB T − 1 is the mean number of the reservoir quanta at temperature T , and k B is the Boltzmann constant.

III. THE NONLINEAR MASTER EQUATION
From Eq. (8) we will deduce a quantum equivalent to the Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy [50,51] of equations similar to the procedure in classical kinetic theory. The term hierarchy refers to the order of correlations between K particles, assumed as an approximation to the exact N -particle master equation (N > K). For example, an equation of evolution for the state of K particles will depend on the state of K + 1 particles, with K + 1 ≤ N . The hierarchy is said to be broken whenever the state of K + 1 particles is reduced as a sum of product of states involving less particles, i.e., introducing lower-order correlations.
Let us consider a subsystem constituted of K particles (K < N ), whose density operator is obtained by calculating the trace over the remaining K + 1, K + 2, ..., N particles in Eq. (9). Then, the reduced density operator is written as The equation of motion for ρ K is obtained by calculating trace over the master equation (8), thus obtaining explicitly, The RHS of Eq.(11) contains terms having the following structures 2 With σ ± = 1 2 (σx ± iσy), where σx, σy, and σ 0 are represented by the Pauli matrices.
with λ, λ ′ = −, 0, +. Detailing the calculations, the term (12) can be written as where the sum ..., since the trace operation is done over the particles K + 1, ..., N . The second term in the RHS of the last line in Eq. (15) is where all the N − K terms are equivalent, since the particles are assumed being identical. Therefore, using the cyclic property of the trace operation (T r(AB) = T r(BA)), it follows that Then, Eq. (15) reduces to Now, considering the trace operation (13), and breaking again the sum ..., we obtain the following four terms Evaluating explicitly each one separately we get: where the factor N − K is present because there are formally equivalent terms; In Eqs. (22) and (23) the results are zero because the trace of a commutator is null, then Eq. (19) becomes Analogously to the above procedure, for the Eq. (14) we have which results again in four terms. Looking at each one separately and following the same procedure as above, we get Therefore, Eq. (25) can be written as Thus, back to Eqs. (12)- (14), we obtain explicitly Using the results in Eqs. (31)-(33) we obtain the following master equation for the reduced K-particle system with K = 1, 2, ..., N − 1, and where h.c. stands for hermitian conjugate. It is worth noting that the master equation for the reduced density operator of K particles ρ K depends on the density operator of K + 1 particles, ρ K+1 . These N − 1 equations are similar to the BBGKY hierarchy, studied in several contexts related to the statistical mechanics of classical fluids [50,51]. For a single representative particle of the system (K = 1), the equation of motion for ρ 1 is noting that it depends on the two-particle density operator ρ 2 , and so on for the whole hierarchy. Thus, we still have a complex problem. In order to obtain a scenario easier to be handled mathematically, some kind of approximation needs to be done. A plausible physical argument to obtain some simplification for a dilute system consists in disregarding the higher order particle correlations. This is achieved when one does the factorization ρ 2 (1, 2) = ρ 1 (1)ρ 1 (2), so now a generic particle is assumed to move in a mean-field produced by all the other particles. This approximation is very similar to the well-known Hartree approximation or mean-field approximation used in atomic and condensed matter physics [52][53][54], which consist in a factorization as a product of single-particle states. Implementing this approximation, Eq. (35) reduces to a generic single-particle master equation where a single-particle effective hermitian Hamiltonian is identified as which contains nonlinear terms, obtained by grouping the second term in the RHS of Eq. (35) to the original commutator, namely T r 2 (σ ± (2)ρ 2 ) = ρ 1 (1)T r 2 (σ ± (2)ρ 1 (2)) = ρ 1 (1) σ ± , reminding that T r(ρÂ) ≡ Â is the mean value of an operatorÂ in the density operator formalism. The second term in the brackets in the Hamiltonian (37) represents a particle acted by an effective polarization field E pol = i(N − 1)Γ σ + /2V originated from the other (N − 1) particles that produce the mean-field effect, proportional to the uniform particle density (N − 1)/V . In this approximation, the structure of the Sudarshan-Lindblad dissipative terms remains the same, although the single-particle Hamiltonian acquires, additional nonlinear terms.

IV. DYNAMICS
By noting that the equations of motion for the single-particle operators mean values are obtained by inserting the RHS of Eq. (36) forρ 1 . Then, we obtain the following nonlinear equations of motion for the volumetric density operators mean values, The first term in the RHS of Eqs. (39), and the first and second terms in the RHS of Eq. (40) come from the evolution due the nonlinear Hamiltonian, while the last terms in both equations originate from the second term in the master equation (36), responsible for the dissipation mechanism. As a first analysis, we will disregard dissipation and consider the evolution under the effective Hamiltonian (37) only. This consideration is valid for processes occurring in a time interval much shorter than the relaxation time, when the system becomes stationary. In this case, the dissipative terms in Eqs. (39) and (40) are disregarded, and the equations of motion for the mean values reduce to that can be solved exactly, once recognizing that is a constant of motion. Thus, the equation of motion for σ 0 becomes a first order separable differential equation, whose solution is where A = tanh −1 ( σ 0 (0) /R). Substituting Eq. (45) into Eq. (42) and solving by integration we obtain the solution where cosh(A) = 1 − ( σ 0 (0) /R) In order to illustrate the theory, the excited state probability ρ ee = (1 + σ 0 )/2 is plotted in Fig. IV as function of time for ρ ee (0) = 1, Γ = 1, and for several values of N . Note that by having disregarded the dissipative terms there is no spontaneous decay. So, in the situation at hand the generic particle does not decay spontaneously when it is initially set in the excited state, namely ρ ee (0) = 1, which is an unstable equilibrium point, since it lies in the excited state with zero dipole moment. The curves in Fig. IV are obtained for ρ ee (0) slightly smaller than 1. The decay profile is due the mean field created by the other particles, and not by dissipative Sudarshan-Lindblad terms. We note that the mean energy of the system Ê = N ω o σ 0 (t) is not a conserved quantity, it depends on the initial state, ρ ee (0), being less than 1, when the spins tend to align along the effective polarization mean field.
In fact, the equations of motion generated by the Hamiltonian (37) are equivalent to those describing the superradiant emission by a collection of two-level atoms in the framework of a semi-classical radiation theory [55], described by a single-particle mean-field Hamiltonian for the individual emitters as pointed out in Ref. [56]. We turn now to Eqs. (39) and (40). For N = 1 the system reduces to a decoupled set of linear equations whose solutions are that contain a transient exponential decay and at long times attain the stationary values σ 0 (∞) = −1/(2n + 1) and σ ± (∞) = 0, respectively. For the system of Eqs. (39)- (40), when N > 1, the analytical solution is possible only for the stationary state, otherwise a numerical calculation is necessary. The stationary solution is obtained by setting d σ 0 /dt = d σ ± /dt = 0, and solving the resulting algebraic equations whose solutions are identical to the case N = 1. Therefore, the stationary state does not depend on the number of particles and is determined solely by the reservoir temperature through the parametern that depends on the absolute temperature T . The deviation from an exponential decay is due to the nonlinear terms and it occurs at transient times. To illustrate the interplay between nonlinear and dissipative effects, in Fig. IV we plotted ρ ee as function of time for ρ ee (0) = 1, Γ = 1, and setting the reservoir temperature at 0K (n = 0), for several values of N . For N = 1 we have a purely exponential decay (solid line) while for N > 1 the deviation from the exponential decay at transient times due the nonlinear mean-field is remarkable (dashed and doted lines). This deviation is evidenced in the inset within the Fig. IV drawn in logarithmic scale. Figure IV is the same as Fig. IV but forn = 0.5; we note that at higher temperature the decay is faster and ρ ee (t) attains an asymptotic nonzero value for the stationary state that will depend on the temperature on the value ofn.

V. BEYOND THE MEAN-FIELD APPROXIMATION
Here we present a brief discussion on how to go beyond the single-particle mean-field approximation, which consists in the factorization ρ 2 (1, 2) = ρ 1 (1)ρ 1 (2), that disregards the two-particle correlations that exist, either in the initial conditions, or being generated by the interaction with the environment, if the gas is quite diluted, as the system of particles that evolve dynamically in time. Under stricter conditions, the correlations between the particles must be considered in order to describe the dynamics more faithfully, and this means that we have to add, at least, a two-particle operator χ(1, 2), such that having the properties T r 2χ (1, 2) = T r 1χ (1, 2) = 0, where T r k means calculating the trace over the degrees of freedom of particle k . The operatorχ(1, 2) contains all kind of two-particle correlations, it can be introduced within the initial conditions for ρ 2 (1, 2), and its evolution in time will be ruled by the master equation. However, according to the BBGKY equation (34), if we set in its left hand side (LHS) the time derivative of a two-particle density operator, the RHS will contain a dependence on a three-particle density operator, ρ 3 (1, 2, 3), and in order to solve the differential equation one should go a step further, namely set the time derivative on the LHS for a three-body operator and we shall get on the RHS of the equation a dependence on a four-particle operator. This procedure goes on until we set the exact equation for the N -particle problem. As this procedure increases in complexity, we have adopted the quite simple single-particle mean-field approximation. To go beyond this approximation, and still to extend the BBGKY equation that involves two-particle correlations, we can adopt the following ansatz for the three-particle density operator, which satisfies the requirements T r 3 ρ 3 (1, 2, 3) = ρ 2 (1, 2), T r 2 ρ 3 (1, 2, 3) = ρ 2 (1, 3) and T r 1 ρ 3 (1, 2, 3) = ρ 2 (2, 3). By introducing the correlation operator (49) in Eq. (50) we get Thus, the BBGKY equation reduces to a closed system of coupled differential equations involving the one-body density operator and the two-body correlation operator. The set of the dynamical equations for ρ 1 (k) andχ(j, k) is still not easy to solve analytically, one could adopt a further simplification, that is known as the linear response approximation, where one looks for the correlations nearby a stationary or equilibrium state of the system, represented by the one-particle operator obtained under the mean-field approximation, namely, ρ st 1 (k). This ansatz reduces the density operator (51) to and the BBGKY becomes a dynamical equation for the correlation operatorχ(j, k) only.

VI. SUMMARY AND CONCLUSIONS
The aim of this paper is to be pedagogical for advanced undergraduate and graduate students by introducing some techniques and approximations that arise in the physics that involves dissipation and nonlinearity in many-particle systems. Here we have considered a system S of N independent spin-1/2 particles interacting with a thermal reservoir R. We started with the dynamical quantum Sudarshan-Lindblad master equation for the N -particle density operator. From this point on we derived the quantum version of the BBGKY hierarchy of K equations for 1 ≤ K < N . Then, by breaking the hierarchy we neglected interparticle correlations, obtaining thus an equation for a representative singleparticle. This is equivalent to introduce the mean-field approximation in the equation for an N particle system. This scenario leads to a quite similar Sudarshan-Lindblad equation for the single-particle density operator; however, a main difference arises: the Hamiltonian is now nonlinear, although the Sudarshan-Lindblad superoperator, responsible for the dissipation and decoherence in S, remains unaffected. In order to provide a useful tool to students interested in many-particle problems, we have discussed the physical approximations involved, we derived the analytical solution for the nonlinear mean-values operators, whose motion is generated by the nonlinear Hamiltonian without dissipative terms, and illustrate numerically the case with dissipation. We also included a discussion on how to go beyond the mean-field approximation by introducing two-particle correlations.