Reconstruction Approach to Quantum Dynamics of Bosonic Systems

We propose an approach to analytically solve the quantum dynamics of bosonic systems. The method is based on reconstructing the quantum state of the system from the moments of its annihilation operators, dynamics of which is solved in the Heisenberg picture. The proposed method is general in the sense that it does not assume anything on the initial conditions of the system such as separability, or the structure of the system such as linearity. It is an alternative to the standard master equation approaches, which are analytically demanding especially for large multipartite quantum systems. To demonstrate the proposed technique, we apply it to a system consisting of two coupled damped quantum harmonic oscillators.

We propose an approach to analytically solve the quantum dynamics of bosonic systems. The method is based on reconstructing the quantum state of the system from the moments of its annihilation operators, dynamics of which is solved in the Heisenberg picture. The proposed method is general in the sense that it does not assume anything on the initial conditions of the system such as separability, or the structure of the system such as linearity. It is an alternative to the standard master equation approaches, which are analytically demanding especially for large multipartite quantum systems. To demonstrate the proposed technique, we apply it to a system consisting of two coupled damped quantum harmonic oscillators.
Introduction.-One of the most intriguing problems in modern physics is understanding the dynamics of open quantum systems [1]. In general, the problem is solving the reduced dynamics of a small quantum system interacting with a large environment. Such interaction leads to seemingly irreversible processes, such as dissipation and decoherence [2]. The control of these effects is topical, for instance, in quantum information processing [3][4][5][6], where the control of dissipation [7][8][9][10][11] and routing of heat flows [12][13][14] have recently attracted great experimental interest.
The foundations for the study of dissipation in quantum systems were laid in the 1960s in terms of the influence functional formalism [15]. Subsequently, the theory of quantum dynamical semigroups [16] has led to a vast amount of theoretical work on quantum master equations [1,2,17]. Several approaches to solve master equations analytically have been presented, including algebraic methods [18][19][20][21], exact diagonalization [22,23], series expansions [24], and effective Hamiltonian approaches [25,26]. However, these techniques are technically demanding, especially for multipartite quantum systems [20].
Here, we introduce an analytical approach, alternative to the master equation techniques, to solve the complete quantum dynamics of dissipative bosonic systems. The idea is to solve the dynamics of the annihilation operators of the system in the Heisenberg picture, and to reconstruct the entire quantum state using a moment expansion of these operators. In essence, we obtain the Schrödinger picture solution while circumventing the need to neither derive nor solve a master equation for the system. The utility of this approach lies in the fact that solving the dynamics of the operators is simple, and the reconstruction step is straightforward. The method itself does not call for assumptions on the initial state or the structure of the system such as linearity.
Although the moment expansion for the quantum state of a single bosonic mode was presented as early as in 1990 by Wünsche [27], its applications have mainly been in quantum state tomography [28,29]. Here, we utilize the expansion to solve the quantum dynamics of bosonic systems.
To demonstrate the utilization of the introduced method, we consider a system of two bilinearly coupled damped quantum harmonic oscillators. Experimentally, this system can be realized for example as coupled coplanar waveguide resonators [30]. Such system is of current interest, for instance, for rapid high-fidelity measurement of superconducting qubits using Purcell filters [31,32], and for transferring heat in quantum circuits at maximal rates using exceptional points [33]. Theoretical work on the system of two coupled quantum harmonic oscillators has been presented, for example, in Refs. [34][35][36]. To the best of our knowledge, however, the analytical solution for the density operator of the composite system has not been reported.
Method.-In this section, we give a description of the reconstruction approach at a general level. The schematic process chart of the method is given in Fig. 1(a). We consider a general bosonic system consisting of N discrete modes and M continua of modes, see Fig. 1(b). In the Schrödinger picture, we assume that its Hamiltonian is of the form whereĤ is polynomial in the system operators, andâ j andB j (ω) are the annihilation operators of the discrete modes and of the continua of modes, respectively. The discrete-mode operators obey the conventional bosonic commutation relations, [â j ,â † k ] = δ j,k and [â j ,â k ] = 0, and the continuous-mode operators obey the continuousmode bosonic commutation relations, B j (ω),B † k (η) = δ j,k δ(ω − η) and B j (ω),B k (η) = 0. The continua of bosonic modes are included into the Hamiltonian to enable, for instance, first-principles modeling of dissipation in the system [37].
To solve the dynamics of the annihilation operators of the system in the Heisenberg picture, the Hamilto- (a) Schematic process chart of the proposed reconstruction approach. First, the Schrödinger picture Hamiltonian is transformed into the Heisenberg picture. Subsequently, the Heisenberg equations of motion for each annihilation operator of the system are obtained and solved. Finally, the density operator of the composite system is reconstructed by inserting the solutions for the annihilation operators, see Eqs.
(3) and (4) for details. For simplicity, the process chart describes a single-mode bosonic system. (b) Schematic of a general bosonic system consisting of N discrete modes and M continua of modes. The annihilation operatorsâj label the discrete modes, and the continuous-mode annihilation oper-atorsBj(ω) label the continua of modes. Schematic illustrations of the energy level diagrams of the discrete modes and the spectra of the continua are provided.
nian is formally transformed into the Heisenberg picture according toĤ H (t) =Û † (t)Ĥ(t)Û (t), whereÛ (t) = T e −i t 0 dτĤ(τ )/ is the temporal evolution operator, T is the time-ordering operator, and is the reduced Planck constant [1]. Since a true Hamiltonian is always Hermitian, the corresponding temporal evolution operator is unitary [38], that is,Û † (t)Û (t) =Û (t)Û † (t) =Î. Thus, the transformation into the Heisenberg picture is simple: identity operators can be inserted appropriately intoĤ H (t) such that all the system operators in the Schrödinger picture Hamiltonian can be substituted by the Heisenberg picture equivalents.
The Heisenberg equations of motion for the annihilation operators of the system reaḋ The commutators on the right sides of the above equations are straightforwardly evaluated with the help of the bosonic commutation relations. However, the existence of an explicit analytical solution to the resulting set of N + M coupled equations of motion depends on the system under interest. In the following, we assume that an explicit, but not necessarily analytical, solution to Eqs. (2a) and (2b) exists.
It is well known that the expectation value of any moment of the annihilation operators can be evaluated once the Heisenberg equations of motion for the annihilation operators are solved [18]. We utilize these expectation values to solve the dynamics of the density operator of a set of bosonic modes. In the Supplemental Material [39], we derive the following expansion for the density operator of an N -mode bosonic field,ρ(t), at any time instant t in terms of the initial normally ordered momentŝ and Ô (t) = Tr[ρ(t)Ô] = Tr[ρ(0)Ô(t)] is the expectation value of the operatorÔ. Here, we have used the fact that the expectation values of operators coincide between the pictures of the quantum mechanics. Equations (3) and (4) demonstrate that the full information on the quantum dynamics of a bosonic system is embedded in the dynamics of its annihilation operators. In the case of a single bosonic mode, N = 1, Eqs. (3) and (4) reduce to the expansion presented originally in Ref. [27].
Consequently, the insertion of the solutions for the annihilation operatorsâ j (t) into the expression for the expectation value, Eq. (4), and insertion into Eq. (3) amounts the solution for the complete quantum dynamics of the system of the discrete bosonic modes. Here, the expectation value is evaluated with the help of the initial density operator of the system. Note that the solution is analytical if and only if that ofâ j (t) is analytical; ifâ j (t) is obtained numerically, the solution is semi-analytical.
A more convenient representation of the density operator may be given in the number basis, where the elements of the density operator assume the form [39] ρ n1,m1,...,n N ,m N (t) Thus, we have reduced the problem of solving the quantum dynamics of the system into that of solving a set of coupled equations (2a) and (2b). Example: Two coupled damped harmonic oscillators.-Here, we consider a case of N = M = 2, that is, a system consisting of two discrete bosonic modes, labeled as M1 and M2, and two continua of modes, labeled as B1 and B2. Specifically, the discrete modes are damped quantum harmonic oscillators which are coupled to each other, as depicted in Fig. 2.
We model the dissipation of the discrete modes to corresponding environments using the Gardiner-Collett Hamiltonian [40]. Within the Markovian approximation, where the coupling between a mode and the corresponding environment does not depend on frequency, the Hamiltonian of the entire system readŝ where ω j are the frequencies, κ j are the energy decay rates andB j (ω) are the annihilation operators of the corresponding environments of the modes, and g is the coupling strength between the modes. The temporal evolution operator of the full system is unitary [39]. Thus, the Heisenberg picture Hamiltonian has exactly the form of Eq. (6), and the Schrödinger picture operators are replaced with the Heisenberg picture equivalents, as argued in the previous section. Consequently, the Heisenberg equations of motion for the annihilation operators are readily obtained aṡ The analytical solution to this set of equations of motion for the annihilation operator of M1 reads [39] FIG. 2. Schematic of the example system of two coupled damped quantum harmonic oscillators. The annihilation operators, the frequencies, and the decay rates of the harmonic oscillators areâj, ωj, and κj, respectively, the annihilation operators of the corresponding heat baths of the harmonic oscillators areBj(ω), and the coupling strength between the harmonic oscillators is g. where andĈ 3 is given in Ref. [39]. The detailed derivations of the results presented in this section are given in the Supplemental Material [39]. Due to the symmetry of the system, Eqs. (8)-(9b) give also the solution ofâ 2 (t) by substituting 1 → 2 and 2 → 1 in the indices of κ j , ω j , a j (0), andB j (ω, 0). As expected, Eq. (8) shows that there are two hybridized modes in the system which decay at finite dissipation rates. The excess operatorĈ 3 ensures that the bosonic equal-time commutation relations hold, [â j (t),â † k (t)] = δ j,k and [â j (t),â k (t)] = 0, and is the only term that contributes to the asymptotic behavior ofâ 1 (t) andâ 2 (t). Consequently, the asymptotic behavior depends only on the initial properties of the baths throughB 1 (ω, 0) andB 2 (ω, 0).
The complete dynamics of the modes are then reconstructed by inserting the solutions for the annihilation operators into Eqs. (4) and (5). Notably, the dynamics are obtained for arbitrary initial conditions, such as initially correlated modes and environments. For simplicity however, we consider below separable initial conditions.
We suppose that the dissipative environments and mode M2 are initially in vacuum states, that is, the initial density operator of the entire system iŝ where the density operators on the right-hand side are those of M1, M2, B1, and B2,ρ (1) (0) is an arbitrary physical density operator, and |0 is the multi-mode vacuum state. Insertion of the solutions forâ 1 (t) andâ 2 (t) given by Eq. (8) together with the initial state, Eq. (10), into Eqs. (3) and (4) yields [39] ρ n1,m1,n2,m2 (t) = f n1 where are the sums of the prefactors ofâ 1 (0) in Eq. (8) for the solutionsâ 1 (t) andâ 2 (t), respectively. Equation (11) shows that the dynamics of the density matrix elements are given as weighted sums over certain off-diagonal elements of the initial density matrix of mode M1.
Specifically, only the part of the initial state of M1 corresponding to the Hilbert space H n = span{|n , |n + 1 , . . .} where n = min [max(n 2 , n 1 + n 2 − m 2 ), max(m 1 − n 1 + m 2 , m 1 )] affects on the dynamics of the density matrix element ρ n1,m1,n2,m2 (t). Moreover, the density matrix elements have damped oscillatory behavior due to f j (t) being decaying and oscillating functions of time, and the decay rates increase with increasing n 1 + m 1 and n 2 + m 2 .
The elements of the reduced density matrices of M1 and M2 can be obtained from Eq. (11) by taking the partial trace as ρ (1) n,m = ∞ l=0 ρ n,m,l,l and ρ (2) n,m = ∞ l=0 ρ l,l,n,m , respectively, resulting in [39] This equation shows that a swap of any quantum state from M1 to M2 up to complex phase arising from the bare evolution is obtained if the modes are non-decaying and in resonance, and the interaction time is chosen such that f 1 (t) = 0 and |f 2 (t)| = 1.
For certain initial states of interest, simple solutions are obtained using Eq. (13). If the initial state of M1 is a coherent state, |α 1 (0) , the states of both of the modes remain as coherent states through the temporal evolution. The dynamics of their coherent amplitudes are given by α j (t) = f j (t)α 1 (0) [39]. If the initial state of M1 is a thermal state with the scaled inverse temper- both of the modes remain in thermal states [39]. The dynamics of their scaled inverse temperatures β j (t) = ω j /[k B T j (t)] are given by Finally, we point out that if the modes are decoupled, g = 0, Eq. (13) reduces to the result for a single damped harmonic oscillator, presented, for example, in Ref. [25].
Conclusions.-In summary, we have introduced an approach to analytically solve the quantum dynamics of bosonic systems. The essence of the method is in reconstructing the quantum state of the system under interest from the moments of its annihilation and creation operators, the dynamics of which is solved in the Heisenberg picture. The method itself does not pose requirements for the initial conditions or the structure of the system. Moreover, the proposed method is particularly practical for obtaining exact solutions for multipartite quantum systems, for which the Heisenberg equations of motion are more convenient to solve than the master equation. To demonstrate the utilization of the method, we have applied it to a system consisting of two coupled damped quantum harmonic oscillators.
In the future, the generality of the method enables it to be applied to various bosonic quantum systems. In particular, it may shed light to the effects of an initial entanglement between a system of interest and its environment, out of the reach for the conventional master equation approaches [41]. The results of the presented two-mode example may possibly be applied to the studies on launching states of a microwave resonator to as propagating waves using the so-called Schrödinger's catapult [42]. 1

MOMENT EXPANSION FOR MULTI-MODE BOSONIC FIELD
In this section, we derive the normally-ordered moment expansion of the density operator of an N -mode bosonic field, and the corresponding expression for the elements of the density matrix in the number basis.
We begin by expressing the density operator of an N -mode bosonic field in the number basis aŝ where |k 1 , . . . , k N = N j=1 |k j j is the multi-mode number state with k j excitations in mode j. In Ref. [1], it is shown that the single-mode projector, |k j j j l j |, can be expanded in the normally-ordered moments as whereâ j (â † j ) is the annihilation (creation) operator of the mode j. Note that the density matrix elements in Eq. (S1) are defined as ρ k1,l1,...,k N ,l N = Tr ρ N j=1 |l j j j k j | , and recall that the trace operator is linear. Consequently, employing Eq. (S2) for both the density matrix elements and the projectors, we may rewrite Eq. (S1) aŝ

TWO COUPLED DAMPED QUANTUM HARMONIC OSCILLATORS
In this section, we consider a system consisting of two discrete modes and two continua of modes. Specifically, the system consists of two bilinearly coupled damped quantum harmonic oscillators, labeled as M1 and M2. We model the dissipation in each harmonic oscillator (mode) as coupling to a bath of harmonic oscillators, labeled as B1 and B2, using the Gardiner-Collett Hamiltonian [3]. Within the Markovian approximation, i.e., frequency-independent coupling strength between the system and the environmental modes, the Hamiltonian readŝ where is the reduced Planck constant, g is the coupling strength between the modes, and ω j , κ j ,â j andB j (ω) are the frequency, the decay rate, the annihilation operator and the annihilation operator of the corresponding bath of the mode Mj. We begin this section by solving the dynamics of the annihilation operators of the system. Then, we solve the complete quantum dynamics of the dissipative modes supposing that initially the heat baths are in zero temperature and only one of the modes is in an arbitrary state. Finally, we consider certain initial states of interest of the non-vacuum mode.