Reconstruction Approach to Quantum Dynamics of Bosonic Systems

TL;DR

This paper introduces a general analytical method for solving the quantum dynamics of bosonic systems by reconstructing the quantum state from operator moments, offering an alternative to traditional master equations especially for complex systems.

Contribution

The paper presents a novel reconstruction-based approach to quantum dynamics that does not rely on initial state assumptions or system linearity, simplifying analysis of large multipartite bosonic systems.

Findings

Successfully applied to two coupled damped quantum harmonic oscillators

Provides analytical solutions without initial state restrictions

Offers a scalable alternative to master equation methods

Abstract

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.

Figures

Equations

$$\hat{H}(t)=\hat{H}\Big{[}\{\hat{a}_{j},\hat{a}_{j}^{\dagger}\}_{j=1}^{N},\{\hat{B}_{j}(\omega),\hat{B}_{j}^{\dagger}(\omega)\}_{j=1}^{M};t\Big{]},$$

$${\dot{\hat{a}}}_j(t)$$

$${\dot{\hat{B}}}_j(\omega ,t)$$

$$\hat{\rho }(t)$$

$$\displaystyle\Big{\langle}\hat{c}_{k_{1},l_{1},\ldots,k_{N},l_{N}}(t)\Big{\rangle}$$

$$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,=\Bigg{\langle}\prod_{j=1}^{N}\frac{1}{k_{j}!l_{j}!}\sum_{q_{j}=-\text{min}(k_{j},l_{j})}^{\infty}\frac{(-1)^{q_{j}}(k_{j}+l_{j}+q_{j})!}{(k_{j}+q_{j})!(l_{j}+q_{j})!}$$

$$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times\left[\hat{a}_{j}^{\dagger}(t)\right]^{l_{j}+q_{j}}\hat{a}_{j}^{k_{j}+q_{j}}(t)\Bigg{\rangle},$$

$${\rho }_{n_1,m_1,\dots ,n_N,m_N}(t)$$

$$=\sum _{k_1=0}^{\text{min}(n_1,m_1)}\dots \sum _{k_N=0}^{\text{min}(n_N,m_N)}\frac{\sqrt{n_1!m_1!\dots n_N!m_N!}}{k_1!\dots k_N!}$$

$$\displaystyle\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\times\Big{\langle}\hat{c}_{n_{1}-k_{1},m_{1}-k_{1},\ldots,n_{N}-k_{N},m_{N}-k_{N}}(t)\Big{\rangle}.$$

$$\hat{H}/\hslash$$

$${\dot{\hat{a}}}_1(t)=-i{\omega }_1{\hat{a}}_1(t)-ig{\hat{a}}_2(t)-i\sqrt{\frac{{\kappa }_1}{2\pi }}\int \text{d}\omega {\hat{B}}_1(\omega ,t),$$

$${\dot{\hat{a}}}_2(t)=-i{\omega }_2{\hat{a}}_2(t)-ig{\hat{a}}_1(t)-i\sqrt{\frac{{\kappa }_2}{2\pi }}\int \text{d}\omega {\hat{B}}_2(\omega ,t),$$

$${\dot{\hat{B}}}_1(\omega ,t)=-i\omega {\hat{B}}_1(\omega ,t)-i\sqrt{\frac{{\kappa }_1}{2\pi }}{\hat{a}}_1(t),$$

$${\dot{\hat{B}}}_2(\omega ,t)=-i\omega {\hat{B}}_2(\omega ,t)-i\sqrt{\frac{{\kappa }_2}{2\pi }}{\hat{a}}_2(t).$$

$${\hat{a}}_1(t)$$

$${\lambda }_{\pm }$$

$${\hat{C}}_{1,(2)}$$

$$\hat{\rho }(0)$$

$${\rho }_{n_1,m_1,n_2,m_2}(t)$$

$$f_1(t)$$

$$f_2(t)$$

$${\rho }_{n,m}^{(j)}(t)$$

$${\hat{\rho }}^{(1)}(0)$$

$${\beta }_j(t)$$

$${\rho }_{n,m}(t)$$

Full text

††preprint: APS/123-QED

Introduction.— One of the most intriguing problems in modern physics is understanding the dynamics of open quantum systems Breuer and Petruccione (2002). 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 Weiss (2012). The control of these effects is topical, for instance, in quantum information processing Nielsen and Chuang (2011); Clarke and Wilhelm (2008); Ladd et al. (2010); Kelly et al. (2015), where the control of dissipation Geerlings et al. (2013); Tan et al. (2017); Silveri et al. (2017); Wong et al. (2019); Silveri et al. (2019) and routing of heat flows Hoi et al. (2011); Partanen et al. (2016); Pechal et al. (2016) 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 Feynman and Vernon (1963). Subsequently, the theory of quantum dynamical semigroups Lindblad (1976) has led to a vast amount of theoretical work on quantum master equations Breuer and Petruccione (2002); Alicki and Lendi (2007); Weiss (2012). Several approaches to solve master equations analytically have been presented, including algebraic methods Barnett and Radmore (1997); Chase and Geremia (2008); Bolaños and Barberis-Blostein (2015); Shammah et al. (2018), exact diagonalization Briegel and Englert (1993); Torres (2014), series expansions Lucas and Hornberger (2013), and effective Hamiltonian approaches Yi and Yu (2001); Klimov et al. (2003). However, these techniques are technically demanding, especially for multipartite quantum systems Bolaños and Barberis-Blostein (2015).

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 Wünsche (1990), its applications have mainly been in quantum state tomography Wünsche (1996); Welsch et al. (1999). 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 Pierre et al. (2019). Such system is of current interest, for instance, for rapid high-fidelity measurement of superconducting qubits using Purcell filters Jeffrey et al. (2014); Bronn et al. (2015), and for transferring heat in quantum circuits at maximal rates using exceptional points Partanen et al. (2018). Theoretical work on the system of two coupled quantum harmonic oscillators has been presented, for example, in Refs. Sandulescu et al. (1987); Chou et al. (2008); Paz and Roncaglia (2008). 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

[TABLE]

where $\hat{H}$ is polynomial in the system operators, and $\hat{a}_{j}$ and $\hat{B}_{j}(\omega)$ 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, $[\hat{a}_{j},\hat{a}_{k}^{\dagger}]=\delta_{j,k}$ and $\left[\hat{a}_{j},\hat{a}_{k}\right]=0$, and the continuous-mode operators obey the continuous-mode bosonic commutation relations, \big{[}\hat{B}_{j}(\omega),\hat{B}_{k}^{\dagger}(\eta)\big{]}=\delta_{j,k}\delta(\omega-\eta) and \big{[}\hat{B}_{j}(\omega),\hat{B}_{k}(\eta)\big{]}=0. The continua of bosonic modes are included into the Hamiltonian to enable, for instance, first-principles modeling of dissipation in the system Dutra (2005).

To solve the dynamics of the annihilation operators of the system in the Heisenberg picture, the Hamiltonian is formally transformed into the Heisenberg picture according to $\hat{H}_{\text{H}}(t)=\hat{U}^{\dagger}(t)\hat{H}(t)\hat{U}(t)$, where $\hat{U}(t)=\mathcal{T}e^{-i\int_{0}^{t}\text{d}\tau\hat{H}(\tau)/\hbar}$ is the temporal evolution operator, $\mathcal{T}$ is the time-ordering operator, and $\hbar$ is the reduced Planck constant Breuer and Petruccione (2002). Since a true Hamiltonian is always Hermitian, the corresponding temporal evolution operator is unitary Mannheim (2013), that is, $\hat{U}^{\dagger}(t)\hat{U}(t)=\hat{U}(t)\hat{U}^{\dagger}(t)=\hat{I}$. Thus, the transformation into the Heisenberg picture is simple: identity operators can be inserted appropriately into $\hat{H}_{\text{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 read

[TABLE]

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 Barnett and Radmore (1997). We utilize these expectation values to solve the dynamics of the density operator of a set of bosonic modes. In the Supplemental Material sup , we derive the following expansion for the density operator of an $N$-mode bosonic field, $\hat{\rho}(t)$, at any time instant $t$ in terms of the initial normally ordered moments

[TABLE]

where

[TABLE]

and $\langle\hat{O}(t)\rangle=\text{Tr}[\hat{\rho}(t)\hat{O}]=\text{Tr}[\hat{\rho}(0)\hat{O}(t)]$ is the expectation value of the operator $\hat{O}$. 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. Wünsche (1990).

Consequently, the insertion of the solutions for the annihilation operators $\hat{a}_{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 $\hat{a}_{j}(t)$ is analytical; if $\hat{a}_{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 sup

[TABLE]

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 Gardiner and Collett (1985). 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 reads

[TABLE]

where $\omega_{j}$ are the frequencies, $\kappa_{j}$ are the energy decay rates and $\hat{B}_{j}(\omega)$ 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 sup . 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 as

[TABLE]

The analytical solution to this set of equations of motion for the annihilation operator of M1 reads sup

[TABLE]

where

[TABLE]

and $\hat{C}_{3}$ is given in Ref. sup . The detailed derivations of the results presented in this section are given in the Supplemental Material sup . Due to the symmetry of the system, Eqs. (8)–(9b) give also the solution of $\hat{a}_{2}(t)$ by substituting $1\rightarrow 2$ and $2\rightarrow 1$ in the indices of $\kappa_{j}$, $\omega_{j}$, $\hat{a}_{j}(0)$, and $\hat{B}_{j}(\omega,0)$. As expected, Eq. (8) shows that there are two hybridized modes in the system which decay at finite dissipation rates. The excess operator $\hat{C}_{3}$ ensures that the bosonic equal-time commutation relations hold, $[\hat{a}_{j}(t),\hat{a}_{k}^{\dagger}(t)]=\delta_{j,k}$ and $\left[\hat{a}_{j}(t),\hat{a}_{k}(t)\right]=0$, and is the only term that contributes to the asymptotic behavior of $\hat{a}_{1}(t)$ and $\hat{a}_{2}(t)$. Consequently, the asymptotic behavior depends only on the initial properties of the baths through $\hat{B}_{1}(\omega,0)$ and $\hat{B}_{2}(\omega,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 is

[TABLE]

where the density operators on the right-hand side are those of M1, M2, B1, and B2, $\hat{\rho}^{(1)}(0)$ is an arbitrary physical density operator, and $|\mathbf{0}\rangle$ is the multi-mode vacuum state. Insertion of the solutions for $\hat{a}_{1}(t)$ and $\hat{a}_{2}(t)$ given by Eq. (8) together with the initial state, Eq. (10), into Eqs. (3) and (4) yields sup

[TABLE]

where

[TABLE]

are the sums of the prefactors of $\hat{a}_{1}(0)$ in Eq. (8) for the solutions $\hat{a}_{1}(t)$ and $\hat{a}_{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 $\mathcal{H}_{n}=\text{span}\{|n\rangle,|n+1\rangle,\ldots\}$ where $n=\min\left[\max(n_{2},n_{1}+n_{2}-m_{2}),\max(m_{1}-n_{1}+m_{2},m_{1})\right]$ affects on the dynamics of the density matrix element $\rho_{n_{1},m_{1},n_{2},m_{2}}(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 $\rho^{(1)}_{n,m}=\sum_{l=0}^{\infty}\rho_{n,m,l,l}$ and $\rho^{(2)}_{n,m}=\sum_{l=0}^{\infty}\rho_{l,l,n,m}$, respectively, resulting in sup

[TABLE]

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, $|\alpha_{1}(0)\rangle$, the states of both of the modes remain as coherent states through the temporal evolution. The dynamics of their coherent amplitudes are given by $\alpha_{j}(t)=f_{j}(t)\alpha_{1}(0)$ sup . If the initial state of M1 is a thermal state with the scaled inverse temperature $\beta_{1}(0)=\hbar\omega_{1}/[k_{\text{B}}T_{1}(0)]$, where $k_{\text{B}}$ is the Boltzmann constant, that is,

[TABLE]

both of the modes remain in thermal states sup . The dynamics of their scaled inverse temperatures $\beta_{j}(t)=\hbar\omega_{j}/[k_{\text{B}}T_{j}(t)]$ are given by

[TABLE]

Finally, we point out that if the modes are decoupled, $g=0$, Eq. (13) reduces to the result for a single damped harmonic oscillator,

[TABLE]

presented, for example, in Ref. Yi and Yu (2001).

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 Carmichael (2013). 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 Pfaff et al. (2017).