Spectral Green’s-function method in driven open quantum dynamics

A method based on spectral Green’s functions is presented for the simulation of driven open quantum dynamics that can be described by the Lindblad master equation in Liouville density operator space. The method extends the Hilbert space formalism and provides simple algebraic connections between the driven and nondriven dynamics in the spectral frequency domain. The formalism shows remarkable analogies to the use of Green’s functions in quantum ﬁeld theory, such as the elementary excitation energies and the Dyson self-energy equation. To demonstrate its potential, we apply the method to a coherently driven dissipative ensemble of two-level systems comprising a single “active” subsystem interacting with N “passive” subsystems—a generic model with important applications in quantum optics and dynamic nuclear polarization. The method dramatically reduces the computational cost compared with simulations based on solving the full master equation, thus making it possible to study and optimize many-body correlated states up to the physically realistic limit of an arbitrarily large N .


I. INTRODUCTION
Open quantum dynamics takes into account the environment (outer degrees of freedom) and so more accurately describes real physical phenomena compared with closed quantum dynamics based entirely on the energy operator (inner Hamiltonian). The environment tends to quench coherent states and to purge quantum information by free thermal decay. On the other hand, it can cause irreversible dynamics that makes it possible to create and keep coherent states by continuously driving the system out of its thermal equilibrium [1][2][3]. This makes driven open systems a fundamental object of quantum theory.
The mathematical description of the dynamics of open quantum systems is nonunitary and generally complex. This especially refers to correlated many-body systems that exhibit collective phenomena depending on their environment. In many cases, the Lindblad master equation approach can be used that retains the positivity of the density operator and introduces the environmental effects through Markovian quantum jump operators that enter the dissipator in a simple algebraic way [1][2][3][4][5][6]. However, the Liouville space containing the trajectories of the density operator describing the evolution of a many-body open quantum system grows exponentially with the number of constituents and the dynamics is sensitive to a large number of physical parameters. Hence efficient mathematical tools are necessary to perform adequate approximations and state space restrictions in order to gain insight into the underlying physics.
For the description of collective phenomena in large-scale closed quantum systems at thermal equilibrium, a method involving the use of Green's functions was developed that statistical physics adopted from quantum field theory [7][8][9]. This method is based on the approximate calculation of correlations between dynamic operators that leads to selfconsistent equations for the observables. Subsequently, this method was extended to nonequilibrium closed and non-Markovian open quantum systems describing various transport phenomena where, besides the standard time-domain formalism, it was transformed to an inhomogeneous spectral problem in Hilbert space [10][11][12][13][14].
Here we propose an extension of the nonequilibrium spectral approach to the important class of driven Markovian open quantum dynamics in the Liouville space of the density operator. To this end, we show that the Green's function for an inhomogeneous spectral problem can be formulated in terms of both Hamiltonian and dissipative parts of the Lindblad master equation. The steady state of the driven system is then obtained by a simple algebraic transform of the nondriven thermal equilibrium. Remarkably, there is a close analogy between our proposed spectral formalism and the use of Green's functions in quantum field theory, including the elementary excitation energies and the Dyson self-energy equation. For a demonstration, we apply the method to a coherently driven dissipative ensemble of correlated two-level systems, a basic model used in quantum optics and dynamic nuclear polarization. We show that because the computational cost of the method is significantly cheaper in comparison with the direct master equation simulation, it is possible to study and optimize the many-body correlated states in the realistically large-scale limit. This enhances possibilities in simulations of many-body driven open quantum dynamics, including the spectral response to the driving and the fast search for the optimal parameter regions.

II. SPECTRAL GREEN'S FUNCTIONS FOR THE LINDBLAD MASTER EQUATION
The dynamics of open quantum systems is described in terms of the Lindblad master equation [4,5], where ρ is the density operator, H is the Hamiltonian describing the (internal) energy of the system, and D is the dissipator that represents the effect of the (external) environment. The latter uses Markovian jumps represented by (dimensionless) operators X j and system-environment exchange rates γ j , The dissipator tends to return the system to the state that is in thermal equilibrium with the environment, while the Hamiltonian contains terms that drive the system out of thermal equilibrium, Dρ th = 0, [H, ρ th ] = 0. Equation (1) describes a driven open quantum dynamics in the Liouville space of the density operator. The master equation preserves the unit trace, Trρ = 1. The superoperator M transfers all operators to traceless operators and hence reduces the dimension of the Liouville space. As a consequence, M is degenerate with a nontrivial zero eigenspace. This eigenspace contains a nonthermal steady state that is eventually established in the driven system, In the fully dissipative case, the zero eigenspace is one dimensional; Eq. (3) uniquely defines the steady state for all initial conditions. In this case, the only traceless solution to Eq. (3) is trivial, Suppose the Hamiltonian can be represented in the form where ζ is a real scalar parameter, the Hermitian operators H 0,1 are nondriving, and the Hermitian operator P contains the driving terms of the Hamiltonian. We will assume that the operator H 1 is dimensionless and H 0 , P, and ζ are measured in frequency units. By extracting the thermal equilibrium part ρ = ρ th +ρ, Trρ = 0 and introducing the superoperators , the homogeneous Eq. (3) is rewritten as an inhomogeneous generalized spectral problem, where ζ plays the role of a spectral parameter and the solution ρ belongs to the subspace of traceless operators. Assuming the validity of Eq. (4) for any real value of ζ , the superoperator F 0 − P − ζ H 1 is nondegenerate and hence invertible. Then the unique solution to Eq. (6) isρ = G(ζ )Pρ th , where the superoperator G(ζ ) must satisfy the equation with the unit superoperator in the right-hand side. The superoperator G(ζ ) is independent of the thermal equilibrium, acts in the subspace of traceless operators, and plays the role of the Green's function of the inhomogeneous spectral problem (6). We call the superoperator G(ζ ) the driven spectral Green's function for Eqs. (1) and (5). By virtue of the previous equations, the steady state is written as Equation (7) can be rewritten in the form where the superoperator G 0 (ζ ) must satisfy the equation Indeed, multiplying both sides of Eq. (9) by the invertible superoperator F 0 − ζ H 1 , we obtain Eq. (7). Equation (10) uniquely defines the superoperator G 0 (ζ ) that is independent of the driving part P of the Hamiltonian. We will call the superoperator G 0 (ζ ) the nondriven spectral Green's function.
we obtain then for the steady state, The dual Eqs. (8) and (11) provide compact formulas for the steady state ρ of the master Eq. (1) as a linear transformation of the thermal equilibrium ρ th defined by the product of the driving superoperator P and the driven and nondriven spectral Green's functions G, G 0 determined by Eqs. (7) and (10) and connected via Eq. (9). The superoperator F 0 − P − ζ H 1 of Eq. (7) linearly depends on ζ , hence the driven Green's function G(ζ ) is rationally extendable into the complex plane of ζ , Here the poles ζ = ζ r and the residues G (r) are given by (suitably normalized) solutions to the homogeneous driven spectral problem (F 0 − P − ζ r H 1 )G (r) = 0. The superoperator F 0 − P − ζ H 1 is real and nondegenerate for real ζ , so the poles have nonzero imaginary parts and exist in complex conjugate pairs. By Eq. (8), the superoperator X (ζ ) and the steady state have rational expansions with the same poles. Similarly, the nondriven Green's function has poles and residues defined by the nondriven Eq. (10). Equations (7) and (10) can be considered to be a Liouville space extension of the spectral Green's function formalism in Hilbert space [10][11][12]. There are noteworthy analogies to quantum field theory. The real parts of the poles ζ = ζ r of the superoperator G(ζ ) in Eq. (12) play the roles of the elementary excitation energies. Equation (9) is a copy of the Dyson equation with the superoperators G 0 , G, and P playing the roles of the bare and dressed propagators and the self-energy [7][8][9].
The advantage of the method introduced in this section is that the use of Eqs. (7) and (8) is generally much less computationally costly than calculating the steady state as the dynamic limit or an element of the zero eigenspace by Eq. (3). Indeed, the former requires only an operator inversion, while the latter needs either calculation of an operator exponent or an operator diagonalization. In addition, the knowledge of the poles and residues of the rational structure (12) and Eq. (8) can be used to evaluate the steady state once for all values of the spectral parameter, thus justifying the importance of the method for spectroscopy. At the poles ζ = ζ r , the superoperator X (ζ ) becomes infinite. Hence, the real values of the spectral parameter closest to the poles, ζ ∼ Re ζ r , define the spectral peaks, i.e., the physical regions where the maximal response of the system to the driving should be expected. The imaginary parts Im ζ r give the Lorentzian widths of the spectral peaks. Furthermore, Eqs. (8) and (11) imply that the superoperators 1 + X (ζ ) and 1 − X 0 (ζ ) are inverse to each other, and so X (ζ ) and X 0 (ζ ) commute for all ζ and are diagonalized in the same basis. Equation (11) admits the formal expansion (convergent for X 0 close to nilpotent) providing the zero, linear, quadratic, etc. responses of the steady state to the driving. It follows from Eq. (13) that the poles of the nondriven Green's function are also poles of the driven Green's function. Equation (11) implies that the latter has extra poles defined by the scalar equation det [1 − X 0 (ζ )] = 0 extracting those values of ζ where the superoperator X 0 (ζ ) has a nontrivial fixed point (an eigenoperator with the unit eigenvalue) X 0 (ζ )ρ = ρ.
The spectral Green's functions we introduced in the frequency domain are closely related to the time-domain Green's functions [7][8][9]. In fact, in the case where the superoperators F 0 , H 1 commute, [F 0 , H 1 ] = 0 (that is the case most important for applications), the nondriven spectral Green's function G 0 (ζ ) can be written as a generalized Fourier transform, whereḠ 0 (t ) is the Green's function of the inhomogeneous nondriven dynamical problem for ζ = 0, Indeed, for any bounded inhomogeneity f , the bounded solution to Eq. (15) is written as We then have This implies that the magnitude −G 0 (ζ )ρ(0) describes the generalized spectrum of the free thermal decay of the traceless part of an initial state ρ(0). The superoperatorḠ 0 (t − t ) of Eq. (14) plays the role of the retarded Green's function that describes the free irreversible decay of correlations between the initial state and the thermal equilibrium. The polynomial resolution ("renormalization") of the perturbation series of Eq. (13) as well as the important links of the spectral Green's functions to the projection methods [15][16][17][18] are given in Appendix A.

III. APPLICATION TO ENSEMBLE OF TWO-LEVEL SYSTEMS
We now illustrate the method of spectral Green's functions by its application to a driven dissipative ensemble of correlated two-level quantum systems-the generic model system to study collective phenomena in quantum optics, magnetism, and quantum information [19][20][21].
The model Hamiltonian that we initially consider is built of one irradiated subsystem (called "active") described by the spin-1/2 angular momentum S and N nonirradiated subsystems (that we call "passive") characterized by spin-1/2 angular momenta I (k) , 1 k N. In the rotating wave approximation, we have (in frequency units) It is assumed that the active and passive level separation frequencies satisfy the condition ω S ω I and the effective irradiation acts along the x axis orthogonal to the quantization z axis and has the strength ω 1 and frequency ω 0 ω I . Then, = ω S − ω 0 characterizes the offset of the irradiation frequency from the level separation frequency of the active subsystem, while the level separation frequency ω I of the passive subsystems remains unchanged. The term H IS describes the interactions of the passive subsystems with the active subsystem that take into account single-quantum passive spin coherences I (k) y and the active-passive interaction strengths A k . This term is the only coherent term of the dipole-dipole interactions that commutes with S z and is preserved in the rotating wave approximation. By a suitable rotation of the transverse spin components, I (k) ± → I (k) ± e ±iφ k , we can always achieve that the interaction strengths A k coincide with their absolute values. We can then write where A is the average absolute value of the active-passive interaction strengths and the factors a k 0 characterize the contributions of the passive subsystems to the active-passive interactions.
To describe the dissipation, we assume that the spin interaction strengths |A k | ω I,S are much smaller than the level separation frequencies. Then the thermal equilibrium is well described by the Boltzmann distribution of the energies along the quantization axis, In the case where the active subsystem is "cold" and the passive subsystems are "hot" with respect to the thermal energy, 012224-3 hω I kT hω S , the thermal equilibrium is approximated as ρ th = 2 −N (1/2 − S z ) where the active subsystem is in the ground state p S ∼ 1, while the passive subsystems have all equally populated levels p I ∼ 0. The typical Lindblad dissipator preserves the thermal equilibrium and has the form Here, 1 , 2 , γ 1 , γ 2 > 0 are the effective active and passive longitudinal and transverse relaxation rates, and V ± are dimensionless jump operators given by Eq. (17). We also take into account the passive transverse relaxation in the simplest collective average form.
The chosen model has two important applications. In quantum optics, it describes an optically irradiated system of unlike two-level atoms, where the active subsystem represents the pumped (solid atomic or molecular) gain medium, while the passive ensemble plays the role of the population inverted amplifier [22,23]. In high-field solid-state dynamic nuclear polarization, it describes a microwave irradiated electronnuclear paramagnetic system, where the active subsystem is formed by a microwave irradiated unpaired electron spin (of a free radical or paramagnetic ion), while the passive subsystems belong to nuclear spins in the proximity of the electron [24,25]. The driving is caused by a time-periodic (optical or microwave) excitation and the rotating wave approximation is applied. It is assumed that the passive dissipation is dominated by the collective relaxation mechanisms [23,[26][27][28]. It is important for our study that in both cases, the passive longitudinal relaxation is relatively slow, so that the following condition is well satisfied: In optics, this is because ω I ω S and so γ 1 ∼ (ω I /ω S ) 3 1 , as follows from the spontaneous emission theory [22,23]. In dynamic nuclear polarization, Eq. (19) is satisfied in the highfield low-temperature limit where γ 1 ∼ (A/ω I ) 2 (1 − p 2 S ) 1 , in accordance with the theory of nuclear relaxation by paramagnetic impurities [27,28]. Hence, condition (19) holds independently of the transverse relaxation rates γ 2 , 2 .
Our next step is to consider the "solid effect" resonance, where the active frequency offset is comparable to the passive frequency, ∼ ω I . In this case, the active spin flips are "synchronized" with the passive spin flops. Using the adiabatic elimination method [15,17], the Hamiltonian (17) is transformed to a two-spin flip-flop Hamiltonian, Here, ζ is the resonance offset ζ = − ω I . The dissipator (18) remains unchanged. Similarly, the case ∼ −ω I leads to an effective two-spin flip-flip Hamiltonian [15,17]. The master equation with the Hamiltonian and dissipative parts defined by Eqs. (20) and (18) preserves the subspace 0 of zero-quantum coherences, [I z + S z , 0 ] = 0. Since ρ th ∈ 0 , the driven dynamics and the steady state are closed in 0 . The dim 0 = [2(N + 1)]!/[(N + 1)!] 2 exponentially grows with N. The typical volume of computer memory limits the feasibly fast spectral simulation within the full master equation to N < 10, which is far from a physically realistic assumption. Remarkably, the method we introduced in the previous section enables one to extend the feasible number of the passive subsystems up to the physically realistic limit of an arbitrarily large N.
In the notations of the previous section, The density operator admits the decomposition with ρ 0,z,± containing only the passive spin components. We have ρ (0,1) ∈ (0,1) , ρ th ∈ (0) where the subspaces (0,1) built of zero-quantum and single-quantum coherences of the active subsystem satisfy the conditions of the projection method described in Appendix A. Equations (A3) and (A4) imply that the projections ρ (0,1) of the steady state are found independently from the equations Applying, to both sides of the first of Eqs. (23), the superoperator F 0 − ζ H 1 , and using Eq. (10) we see that the first of Eqs. (23) is equivalent to the equation For any operator ρ (0) ∈ (0) , we obtain with ρ ± containing only passive spin components and D (1) denoting the longitudinal part of the dissipator in Eq. (18). Indeed, ρ (0) commutes with both H 0 , H 1 and S ± S z = ∓S ± /2. For any operators ρ ± , we have where D (1) I is the longitudinal part of the passive dissipator in Eq. (18) and is defined in Eq. (19). The latter implies that the passive longitudinal relaxation makes a negligible contribution to the dynamics in the subspace (1) . As a result, by virtue of Eq. (10) that defines the nondriven Green's function, Since [P ± , ρ ∓ S ± ] = 0, we obtain

012224-4
Equation (24) can be rewritten then as It is seen that the right-hand sides of Eqs. (25) are fully determined by the dimensionless magnitudes a k that participate in the expressions for V ± and P ± in Eqs. (17) and (21) and by three dimensionless physical parameters: the effective spectral parameter f , the effective irradiation strength ω, and the ratio γ = 1 /γ 1 of the longitudinal relaxation rates of the passive and active subsystems. Note, also, that L(P + + P − ) = L(P + ) + L(P − ). It follows from Eq. (25) that in the limit of a small driving ω → 0 or large values of the spectral parameter f → ±∞, the effect of the driving is quenched, so the thermal equilibrium is preserved. In terms of the decomposition defined by Eqs. (22), we have and no polarization of the passive ensemble is created. In the physically reasonable limit γ = 1 /γ 1 → ∞, large driving values ω → ∞, and small values of the spectral parameter f → 0, the steady state is approximated by the operator that annihilates the active longitudinal dissipation and commutes with both P ± . In terms of Eqs. (22), this gives which corresponds to the fully polarized (population inverted) state of the passive ensemble, The intermediate values of the spectral parameter generate an "absorption line" I z (ζ ) that is zero at ζ → ∞ and has a peak of the maximal polarization (26) at ζ = 0. The shape and the width of the absorption line can be estimated in the "mean-field" approximation obtained by setting the magnitudes a k in Eq. (17) to all be equal, a k = 1. In Eqs. (21) and (25), we then obtain V ± = I ± , the passive subsystems become identical, and the dynamics is fully defined by the components of the total passive spin. This case is simplified by representation of the passive ensemble by a single angular momentum I with the spin quantum number I = N/2 similar to the Dicke model [26,[29][30][31][32]. The corresponding occupation numbers are defined as n = (n + − n − )/2 = −I, −I + 1, . . . , I, where n ± are the numbers of the passive subsystems respectively in the excited and ground state. In this ansatz [33], The operators P 0 and ρ 0 , ρ z are diagonal in the basis generated by the occupation numbers, so the contribution of the Hamiltonian part of Eq. (25) becomes zero. For any diagonal operatorρ = nρ n |n n|, we obtain, by virtue of Eq. (27), In terms of Eq. (22) and the new notations, Eq. (25) is equivalent to the system From the point of view of applications, the limit of large active or small passive longitudinal relaxation γ → ∞ is particularly important. In this limit, Eq. (28) enables the analytical solution valid for any N (in the basis generated by the occupation numbers, denotingη = 1 + η), The second of Eqs. (23) then gives The spectral character of the steady state is determined by the poles of the driven Green's function that annihilate the denominators in Eqs. (29) and (30), This gives N + 1 pairs of poles that are exactly calculated as (m = 1, . . . , N) The first pair ζ 0 is the poles of the nondriven Green's function. Equation (29) can be used to estimate two important steady-state characteristics of the driven open quantum system, i.e., the polarization and self-correlation of the total z component of the passive spin, I z = Tr (ρ 0 I z ), I 2 z = Tr ρ 0 I 2 z . By virtue of Eq. (29), proceeding to continuous integral approximations, we obtain the expressions valid for arbitrarily large values of N,  Fig. 1(a). It is evident that at ζ ∼ 0, the passive ensemble is almost fully polarized (population inverted) with I z ∼ −N/2 that is accompanied with the creation of correlations between the passive subsystems. According to Eqs. (31), the poles of the driven Green's function are distributed on the complex plane symmetrically and densely around the origin ζ = 0 [see Fig. 1(b)], leading to a single "absorption line" composed of 2N Lorentzian peaks at ζ = Re ζ m with widths |Im ζ m |. It follows from this analysis that for large values of η 0 , the spectral width grows linearly with , exactly as for a single passive subsystem N = 1, while it grows nonlinearly as ∼ √ N with the number of passive subsystems. Indeed, for large N, we have max |ζ m | = |ζ 1 | ∼ √ η 0 (N + 1)/2π . To illustrate the applicability of the method to optimization problems, consider now an ensemble of many active subsystems, each "serving" N passive subsystems. Physically, the active transverse relaxation rate 2 that influences the passive polarization is caused by active spin-spin interactions and grows quadratically with the spatial concentration c of active subsystems. We can write 2 = 0 2 ξ 2 , where ξ = c 0 /c is the dimensionless relative concentration with respect to some reference concentration c 0 , and 0 2 is the rate for c = c 0 . For 2 γ 2 + 1 /2, the simulation by the first of Eqs. (32) implies that the total peak polarization of the passive ensemble ξ I z at ζ = 0 has an active concentration optimum whose location and peak value both increase with N; see Fig. 1(c). Equations (28) enable us to also analyze the effect of the relative longitudinal relaxation described by the parameter γ . It follows from these equations that the active and passive polarizations are connected as S z + I z /γ + 1/2 = 0. Restricting to the limit η 0 → ∞, Eqs. (28) are resolved by a simple recurrence in the basis generated by the occupation numbers. As a result, for large N, the peak ζ = 0 dependence of the active and passive polarizations on the relative relaxation parameter γ is well described as For N 1, the value γ /N = 1 can be treated then as the critical value for the second-order (continuous) phase transition between regimes dependent on and independent of the longitudinal relaxation: for γ < N, the active subsystem is fully saturated, while for γ > N, the passive ensemble is fully polarized; see Fig. 1(d). This links our model to phase transitions predicted in the Dicke model [31,32].
The details of the derivation of Eqs. (32) and (33) are given in Appendix B. There we also briefly discuss the general situation where the relaxation is dominated by individual spin mechanisms giving links to the kinetic Monte Carlo algorithm [1,17,18].

IV. CONCLUSION
We have proposed a method of simulation of driven Markovian open quantum dynamics based on Green's functions in a spectral frequency domain. We demonstrated that the method is computationally highly efficient and opens up further ways in the simulation, spectroscopy, and optimization of manybody quantum dynamics in the realistically large-scale limit.

ACKNOWLEDGMENT
This work was funded by the British Engineering and Physical Science Research Council (EPSRC) through Grant No. EP/N03404X/1 (W.K.).
We have X 0 = G 0 P, where P is proportional to the commutation superoperator with the driving P. Hence, all operators commuting with P belong to the zero subspace V 0 of the superoperator X 0 . It means that the latter is degenerate, having the zero eigenvalue of multiplicity m = dimV 0 N 0 , where N 0 is the dimension of the Hilbert space of the quantum problem. Indeed, all operators diagonal in the basis where P is diagonal belong to V 0 . The characteristic polynomial of X 0 then has the form where π (x) is a polynomial with nonzero roots. Equation (A1), then, implies whereπ (x) is a polynomial. It means that Eq. (A1) describes a renormalization of the major coefficients of the generally divergent infinite perturbation series of Eq. (13) of the main text in such a way that the series is truncated to an always convergent polynomial expression. Suppose (that is typically the case) that the Liouville space = (0) + (1) is decomposed into two components that are invariant in the nondriven system and coupled by the driving, with ρ th ∈ (0) . We obtain, for the dynamics of the density operator projections ρ (0,1) ∈ (0,1) , where A (0,1) are the restrictions of the nondriven superoperator to the subspaces (0,1) , If the nondriven dynamics in the subspace (1) is much faster than its exchange with the subspace (0) , then the subspace (1) can be adiabatically eliminated. The dynamics in the subspace (1) is well approximated by the quasiequilibrium, ρ (1) = (A (1) ) −1 Pρ (0) = G 0 (ζ )Pρ (0) .

APPENDIX B: MATHEMATICS OF MODEL EXAMPLE
Equations (32) of the main text describing the polarization and self-correlation of the total z component of the passive spin are obtain by proceeding from the discrete set of the occupation numbers to the continuous interval x = n/I ∈ [−1, 1] and replacing discrete summations over n by integrals with respect to x using the smallness of the discrete step 1/I in the interval [−1, 1] [26]. For example, Note, finally, that Eqs. (25) are of the Lindblad form. They can be treated by unraveling in Hilbert space using the kinetic Monte Carlo method [1]. Here only four jump operators V ± , P ± are involved in the computation scheme. In the case where the transverse relaxation of the passive ensemble is dominated by the individual dephasing mechanism 2γ 2 k L(I (k) z ) with the strong rate γ 2 γ 1 , the dynamics remains closed in the subspace spanned by the z components I (k) z of the passive spins. In this subspace the collective Lindblad terms L(V ± ), L(P ± ) are split into sums of individual terms k L(I (k) ± ), k L(I (k) ∓ S ± ). In the latter case, the kinetic Monte Carlo scheme is reduced to sign permutations in a subsequence of N + 1 symbols. This extends the feasible number of passive subsystems from N ∼ 10 to N ∼ 10 3 ; see Refs. [17,18] for details.