Extending third quantization with commuting observables: a dissipative spin-boson model

We consider the spectral and initial value problem for the Lindblad–Gorini–Kossakowski–Sudarshan master equation describing an open quantum system of bosons and spins, where the bosonic parts of the Hamiltonian and Lindblad jump operators are quadratic and linear respectively, while the spins couple to bosons via mutually commuting spin operators. Needless to say, the spin degrees of freedom can be replaced by any set of finite-level quantum systems. A simple, yet non-trivial example of a single open spin-boson model is worked out in some detail.


Introduction
Exact solutions of simple but nontrivial models describing characteristic physical phenomena are of paramount importance for understanding statistical physics in nonequilibrium.While there is an abundance of such exactly solvable nonequilibrium models in classical statistical mechanics (see, e.g., a review paper [1]), very few explicit analytic approaches are known in the quantum realm (e.g., review papers [2,3]).
A broad route to quantum non-equilibrium physics leads via the theory of open quantum systems, especially in the many-body realm.Considering a large (manybody) quantum system, one can often describe the dynamics of its (so-called central) parts within the framework of the so-called Markov approximation, neglecting the back information flow from the rest of the system (the so-called environment) to its central part.The differential equation describing the density matrix of the central system within the Markov approximation and the rotating-wave approximation [4] is the unique mathematical evolution law that preserves the hermeticity, positivity, and trace of the density matrix, called the Lindblad-Gorini-Kossakowski-Sudarshan [5,6] equation, or Lindblad equation for short.We note that the many-body Lindblad equation provides a perfect mathematical platform for preparing engineered quantum states or quantum phases of matter within cold atom and ion trap setups [7,8].
Some time ago, one of the authors developed canonical formalism of quantization over operator spaces for the diagonalization of the generator of the many-body Lindblad equation -the so-called Liouvillian superoperator, or completely solving the Lindbladian initial value problem, for a general quadratic Hamiltonian and a set of Lindblad jump operators which are linear in canonical creation/annihilation operators.The original proposal for fermionic systems [9,10] was later extended to bosonic systems [11] (see also Ref. [12] for a more abstract discussion of quantization over operator spaces), and further developed by other authors [13,14,15].Specifically, the latter reference extended the technique to include quadratic Hermitian jump operators, allowing the analytical treatment of nonequilibrium phase transitions [16] and crossovers between ballistic and diffusive as well as quantum and classical transport [17,18].Note, however, that even within the class of linear jump operators, one can discuss non-trivial critical phenomena in translationally invariant [19,20] and so-called boundary-driven systems [3,21,22] (in the latter, the jump processes are confined to the boundaries of the system).
Experience has shown that the kinds of systems that can be treated efficiently under the closed system (2nd quantization) formulation can also be treated efficiently under the open system (3rd quantization) formulation.However, there is one type of systems with very important applications, e.g., in quantum optics, that has somehow been left out so far, namely spin-boson systems.Some aspects of integrability and exact solvability of these systems in the closed system framework have been extensively discussed in the literature (see, e.g.[23] for the Rabi model and [24] for the Jaynes-Cummings and Dicke models).A recent case is the exactly solvable model of an electron in a driven harmonic oscillator with Rashba coupling [25,26].Coupling to a thermal bath can be treated numerically [27], but the analytical approach to such problems is precisely the goal of the present work.Here we essentially provide a small extension of the third quantization method that allows us to include additional degrees of freedom with finite-dimensional phase space, provided that these degrees of freedom enter the Hamiltonian and jump operators only through commuting operators.

Formal solution
The aim of this paper is to solve the following Lindblad-Gorini-Kossakowski-Sudarshan master equation for n particles: where H and L µ are the Hamiltonian and Lindblad operators, respectively.ρ is the density operator describing the state of the n particles.The operators H and L µ are given by: where a and a † are n-dimensional vectors of canonical bosonic annihilation/creation operators, and σ is an n-dimensional vector of mutually commuting Hermitian operators, with a finite discrete spectrum.We also assume that [a j , σ j ′ ] = 0.The matrix H is Hermitian (H † = H), the matrix K is symmetric (K = K T ), and the matrix Ω is arbitrary.The vectors l µ , k µ , and w µ are n-dimensional vectors of constants.Below we follow the notation and formalism developed in Ref. [11] extending n bosonic degrees of freedom with additional n finite level quantum systems (e.g.spins where σ j can be considered as their z−projections).We introduce operators âν,j and â′ ν,j , where ν = 0, 1 and j = 1, ..., n: Here, the superscripts L and R indicate the left-and right-multiplication maps acting on a vector space of operators The operators (4) satisfy the almost-canonical commutation relations: In terms of these operators, we rewrite the Liouvillean as follows: where we define n × n matrices: To get rid of terms linear in â and â † we use a map â → â+s, â′ → â′ +s ′ where s, s ′ are constant n-vectors to be determined.Let us define b = (â, â′ ) T and σ = σ L , σ R T .The Liouvillean can be compactly written in a matrix form where In order to eliminate the linear terms, we must fulfill the following condition for the vector d: If S is nonsingular, then we can compute the vector d = − 1 2 S −1 Gσ.Further, we can calculate the last term in Eq. ( 10): By literally following the procedure introduced in [11], we successfully transform the Liouvillean into its final form: where The rapidities β r and the matrix P are obtained from the diagonalization of the matrix X, namely: and the matrix Z is a solution of the continuous Lyapunov equation If all the rapidities are non-negative, ∀β r ≥ 0, one can obtain non-equilibrium steady states |NESS s ⟩ and all the corresponding decay modes or more precisely ζr and where

Examples
Let us consider two simple examples using the above formalism.First, let us take H = 0 and L 1 = a+z 1 σ z , where z 1 ∈ C and σ z is the Pauli matrix for the spin degree of freedom.
Thus we have M = 1, W = |z 1 | 2 , and E = z1 (H = K = Ω = N = L = F = 0).We compute the vector d which equals to and from here we get Ŝ0 = tr M = 1.Combining the b and d vectors, we find that the result can be written in a simple form as b , where we have defined ÂL/R = âL/R + z 1 σ L/R , which are, in a sense, bosonic operators shifted by a z 1 σ z operator.
In this case X = ∆ is already diagonal (P = I), and the rapidities are For each sector s = (s L , s R ) ∈ {−1, 1} 2 , we can compute the NESS: where |α 1,2 ⟩ with α 1 = −z 1 s L and α 2 = −z 1 s R are coherent states, and |s L ⟩⟨s R | represents the spin degree of freedom.One can also compute all the decay modes |Θ s m ⟩ for each s and m.
In the second example, we consider the case where H = 0, L 1 = a + z 1 σ z and L 2 = z 2 σ z , with z 1 , z 2 ∈ C.This modification leads to a change in W compared to the previous example, and we have Additionally, an extra term appears in Ŝ0 , which is now given by: however, the vector d remains unchanged.The extra term introduced by L 2 (i.e., the term involving z 2 ) in Ŝ0 contributes only to a spin-dependent dephasing effect.

Initial condition and expectation values of spin
Here, we present a calculation of the time evolution of the expectation value of spin, ⟨σ µ (t)⟩ = Tr (σ µ ρ(t)) (µ = x, y, z), given an initial condition.For this calculation, we will focus on the first example, where L 1 = a + z 1 σ z .We start by writing a general solution for the density matrix: where c s m are coefficients related to the initial condition |ρ 0 ⟩ = s |ρ s 0 ⟩ ⊗ |s L ⟩⟨s R |, with |ρ s 0 ⟩ being the bosonic part.To compute the coefficients c s m , we need to solve the following system of linear equations for c s m obtained by combining Eqs. ( 23), (25), and (27): We consider a scenario where the system is prepared in the initial state |ρ 0 ⟩ = |0⟩⟨0|⊗ 1 2 (I +σ x ).In Fig. 1, we illustrate the numerical convergence of the coefficients c s m for different truncations of the linear system Eq.( 28), and two choices of the parameter z 1 .
Once the coefficients c s m are calculated, one can compute the expectation values ⟨σ µ (t)⟩.For the example discussed above, the expectation values at t = 0 are ⟨σ x (0)⟩ = 1, ⟨σ y (0)⟩ = ⟨σ z (0)⟩ = 0.It is evident that only the non-zero components will exhibit time-dependence, which, in this example, is ⟨σ x (t)⟩.To test our results, we also computed ⟨σ x (t)⟩ using the Python library QuTiP [28,29], which allows for the computation of the time evolution of the density matrix following the Lindblad master equation.Numerical calculations for the time dependence of ⟨σ x (t)⟩ (Fig. 2a) are within numerical precision fully consistent with the conjectured analytical expression: In Fig. 2b, we illustrate the discrepancies between ⟨σ x ⟩ computed using the coefficients c s m obtained from solving the linear system (28) for different truncations (trunc) and the analytical expression (29).As expected, the discrepancy decreases as we increase trunc.However, beyond trunc ≈ 18, the solution does not improve further due to limitations in numerical precision when inverting the linear system (28) to obtain the coefficients.
We can analytically justify the expression in the limit |z 1 | → 0 (29).To obtain the coefficients c s m , we perform a Taylor expansion of the expression (28) up to quadratic terms in z 1 .This expansion results in a 4 × 4 system of equations to solve for each s = (s L , s R ): where the solution for the coefficients is given by: We also need the expectation values of σ x for the slowest decay modes to obtain the desired result.Following the notation in [11], i.e., (A|ρ⟩ = Tr(Aρ), we have: (σ where δ is the Kronecker delta.The final result in the limit |z 1 | → 0 for the slowest decay modes (λ (0,0) = 0, λ (0,1) = λ (1,0) = −1) is then: If we take the logarithm and approximate it for |z 1 | → 0, we arrive at the analytical expression (29).Additionally, from Fig. 2a, we observe that Eq. (33) approximates the analytical expression (29) and the QuTiP solution well even for z 1 = 0.2 and t ∼ 1.
It is essential to emphasize that the expressions ( 29) and (33) are valid only for the specific initial condition |ρ s 0 ⟩ ∼ |0⟩⟨0|.For different initial conditions, such as |ρ s 0 ⟩ ∼ |n⟩⟨n| with n ̸ = 0, one must include higher orders of m in the computation of the coefficients c s m to obtain the correct result.For completeness, let us provide the solution for the second example with L 1 = a + z 1 σ z and L 2 = z 2 σ z .As mentioned before, the additional term |z 2 | 2 σ L − σ R 2 from Eq. ( 26) just adds trivially to the spin-dependent dephasing.Thus, combined with the conjectured expression (29), we obtain: In Fig. 3, we present the time evolution (34) for various values of z 2 , expressed in terms of z 1 .Evidently, for large t, the solutions exhibit exponential decay with a decay rate proportional to |z 2 | 2 .

Limiting cases
Lastly, let us consider two limiting cases.In the first case suppose again |z 1 | → 0. Then we have ⟨σ In the second case, suppose slightly modified L 1 = wa + σ z = w(a + σ z /w) and |w| → 0 such that |w| 2 t ≪ 1 for some value of t.Now we get and we reconstructed the exponential decoherence of the Lindbladian L = σ z with the rate γ = 2 (note the factor 2 in the definition of the Liouvillean (1)).

Conclusion
In conclusion, we have studied the dynamics of a quantum system -which is quadratic in bosonic degrees of freedom and linear in additional commuting (classical) degrees of freedom -in the framework of the formalism of the third quantization.In particular, we succeeded in diagonalizing the Liouville superoperator and obtaining nonequilibrium states and the corresponding decay modes for various combinations of eigenvalues of the mutually commuting Hermitian operators.A special case is the harmonic oscillator with spin degrees of freedom.Analyzing two specific examples solved exactly, we observed intriguing behavior with exponential decay rates of the spin expectation values.
For the first example, we analyzed the time evolution of the expectation value of spin ⟨σ x (t)⟩ for L 1 = a + z 1 σ z .We have verified numerically the dependence on the parameter z 1 , but the analytic proof of full generality remains an open question.In the second example, we considered a modified system with the additional term L 2 = z 2 σ z in the Lindbladian.The analysis of this case showed an exponential decay with a decay rate proportional to |z 2 | 2 for large t.The results emphasize the influence of the parameter z 2 on the long-term relaxation behavior of the system.Moreover, we have investigated two limiting cases and recovered the known time dependence for the Lindbladians L = a and L = σ z .
The results presented provide valuable contributions to understanding the dynamics of open quantum systems and provide a basis for further studies of related physical phenomena in which spin degrees of freedom are coupled with bosonic ones.For example, these findings could be particularly important for understanding decoherence in spin systems coupled to a bosonic thermal bath and confined in quadratic potential traps [27].Such studies could provide important insights into the behavior of quantum systems under realistic experimental conditions and pave the way for advances in quantum technology and quantum information processing.