Exact solution of the master equation for interacting quantized fields at finite temperature decay

The study of the dynamics and evolution of quantum systems under environmental influence presents a significant challenge for both experimental and theoretical physics. In this context, the Lindblad master equation emerges as a fundamental tool, offering a formal framework that accurately describes the Markovian dynamics of open quantum systems [1, 2, 3, 4]. This equation not only captures the dissipation processes within the system but also incorporates thermal fluctuations, which introduce additional channels for loss and excitation. However, obtaining an exact analytical solution to the master equation remains a considerable challenge, especially in finite-temperature systems where the complexities of interactions and decoherence further complicate the pursuit of precise results [5, 6, 7, 8, 9]. Consequently, numerical methods often become the primary tool for analyzing such systems, although they typically demand extensive computational resources.

Recent developments in cavity quantum electrodynamics have led to the derivation of effective master equations for quantum systems coupled to bosonic modes, resulting in a substantial reduction of the Hilbert space. While these equations are able to describe certain system behaviors, they often lack the property of complete positivity, which can lead to inaccurate predictions regarding system dynamics [10, 11, 12]. This highlights the need for exact solutions to the Lindblad master equation, as only such solutions ensure physical consistency and accuracy in describing the behavior of open quantum systems.

On the other hand, in the analysis of the evolution of one or more electromagnetic field modes within a cavity, a technique based on superoperators has been developed that, through non-unitary transformations, reformulates the Lindblad master equation into a more tractable form. This method simplifies the analysis by eliminating complex terms, such as quantum jumps, and facilitates the application of solutions to various initial conditions, including coherent and number states [13, 14, 15, 16]. In this context, recent investigations have demonstrated that the Markovian dynamics of decaying coupled harmonic oscillators can be effectively described using a non-Hermitian Hamiltonian. This approach has garnered increasing interest due to the unique properties of these systems, particularly regarding parity-time ($\mathcal{PT}$) symmetry and exceptional points (EPs) [17, 18, 19, 5, 20, 21, 22, 23, 24, 25]. Although significant progress has been made in the zero-temperature regime, analyzing systems at finite temperatures remains a considerable challenge. In this work, we aim to extend this approach to the finite-temperature regime, exploring the possibility of obtaining analogous results under such conditions.

The structure of this article is as follows: in Section II, we introduce the theoretical model that forms the basis of our research. In Section III, we show how, by means of two non-unitary transformations, the Lindblad master equation for interacting quantized fields at finite temperature can be recast into a von Neumann-like equation with an effective non-Hermitian Hamiltonian—one of the key results of this work. Section IV presents the exact solution to the Lindblad master equation and includes an example illustrating the photon coincidence rates for two indistinguishable photons initially interacting within a cavity. Finally, in Section V, we summarize our conclusions.

We consider the interaction between two bosonic fields, denoted by $\hat{a}_{1}$ and $\hat{a}_{2}$, each subject to distinct Markovian decay processes characterized by the loss rates $\gamma_{1}$ and $\gamma_{2}$, respectively. The fields are also coupled to a thermal reservoir with an average thermal excitation number $\bar{n}{\text{th}}\geq 0$, as depicted in Fig. 1. Due to the differences in their decay rates, each field experiences its own dissipation dynamics. The Markovian evolution of the density matrix $\hat{\rho}$, in the interaction picture, is governed by the Lindblad master equation[2, 1, 3]

$$\begin{split}\frac{d\hat{\rho}}{dt}=&-\mathrm{i}g\hat{S}\hat{\rho}+\left(\bar{n}_{\text{th}}+1\right)\sum_{j=1}^{2}\gamma_{j}\left[2\hat{a}_{j}\hat{\rho}\hat{a}_{j}^{\dagger}-\{\hat{a}_{j}^{\dagger}\hat{a}_{j},\hat{\rho}\}\right]\\ &+\bar{n}_{\text{th}}\sum_{j=1}^{2}\gamma_{j}\left[2\hat{a}_{j}^{\dagger}\hat{\rho}\hat{a}_{j}-\{\hat{a}_{j}\hat{a}_{j}^{\dagger},\hat{\rho}\}\right],\end{split}$$

(1)

where $\hat{a}_{j}$ and $\hat{a}_{j}^{\dagger}$ are the annihilation and creation operators, respectively, for each mode $j=1,2$, and $g$ represents the coupling strength between the two bosonic modes [26, 27, 28, 29, 30, 31, 32, 33].

The right-hand side of Eq. (1) contains three primary contributions:

1. 1.

   The first term, $-\mathrm{i}g\hat{S}\hat{\rho}$, describes the coherent interaction between the two bosonic modes. This interaction is represented by $\hat{S}\hat{\rho}=[\hat{a}_{1}\hat{a}_{2}^{\dagger}+\hat{a}_{1}^{\dagger}\hat{a}_{2},\hat{\rho}]$, which models the exchange of excitations between the modes due to their coupling.

2. 2.

   The second term accounts for the loss of excitations due to interaction with the thermal environment. The factor $\bar{n}_{\text{th}}+1$ reflects the total number of excitations, including those induced by the thermal reservoir. The expression inside the summation models the dissipation of excitations through the annihilation operators $\hat{a}_{j}$.

3. 3.

   The third term captures the gain of thermal excitations from the environment. Here, $\bar{n}_{\text{th}}$ represents the average number of thermal excitations at equilibrium, and the term in square brackets describes the gain of excitations through the creation operators $\hat{a}_{j}^{\dagger}$.

Together, these terms provide a comprehensive description of the dynamics of the bosonic system, highlighting how both the mutual interaction and the thermal environment influence the dissipation and excitation of the modes.

In the study of open quantum systems, superoperators provide a useful framework for representing complex dynamics in a more accessible way. Unlike conventional operators in quantum mechanics, which act on state vectors, superoperators operate directly on the density matrix. This feature allows them to comprehensively capture both the coherent evolution and the dissipative effects described by the Lindblad master equation.

With this perspective, we introduce the following superoperators [13, 14, 15, 16]

$$\hat{J}_{j}^{(-)}\hat{\rho}:=2\hat{a}_{j}\hat{\rho}\hat{a}_{j}^{\dagger},\qquad\hat{J}_{j}^{(+)}\hat{\rho}:=2\hat{a}_{j}^{\dagger}\hat{\rho}\hat{a}_{j},$$

(2)

and

$$\begin{split}\hat{L}_{j}^{(-)}\hat{\rho}&:=\hat{a}_{j}^{\dagger}\hat{a}_{j}\hat{\rho}+\hat{\rho}\hat{a}_{j}^{\dagger}\hat{a}_{j},\\ \hat{L}_{j}^{(+)}\hat{\rho}&:=\hat{a}_{j}\hat{a}_{j}^{\dagger}\hat{\rho}+\hat{\rho}\hat{a}_{j}\hat{a}_{j}^{\dagger}=\left[\hat{L}_{j}^{(-)}+2\right]\hat{\rho},\end{split}$$

(3)

with $j=1,2$. From these definitions, we can derive the following relations

$$\begin{split}\left[\hat{J}_{j}^{(-)},\hat{J}_{j}^{(+)}\right]\hat{\rho}&=4\left(\hat{L}_{j}+1\right)\hat{\rho},\\ \left[\hat{L}_{j},\hat{J}_{j}^{(\pm)}\right]\hat{\rho}&=\pm 2\hat{J}_{j}^{(\pm)}\hat{\rho},\\ \left[\hat{J}_{1}^{(\pm)}+\hat{J}_{2}^{(\pm)},\hat{S}\right]\hat{\rho}&=0,\end{split}$$

(4)

where $\hat{L}_{j}\hat{\rho}\equiv\hat{L}_{j}^{(-)}\hat{\rho}$. The operators defined as $\hat{A}_{0}\hat{\rho}:=\frac{1}{2}(\hat{L}_{j}+1)\hat{\rho}$ and $\hat{A}^{(\pm)}\hat{\rho}:=\frac{1}{2}\hat{J}_{j}^{(\pm)}\hat{\rho}$, satisfy the commutation relations of the $\mathfrak{su}(1,1)$ algebra, i.e. $[\hat{A}_{0},\hat{A}^{(\pm)}]\hat{\rho}=\pm\hat{A}^{(\pm)}\hat{\rho}$, $[\hat{A}^{(-)},\hat{A}^{(+)}]\hat{\rho}=2\hat{A}_{0}\hat{\rho}$. This simple discovery, yet of great significance, enables us to carry out the following transformation

$$\hat{\rho}=e^{-\eta\left[\hat{J}_{1}^{(+)}+\hat{J}_{2}^{(+)}\right]}\hat{\varrho}$$

(5)

in the Lindblad master equation (1), to obtain

$$\begin{split}\frac{d\hat{\varrho}}{dt}=&-\textrm{i}g\hat{S}\hat{\varrho}-G(\eta)\hat{\varrho}+\sum_{j=1}^{2}\gamma_{j}\left\{\hat{J}_{j}^{(-)}(1+\bar{n}_{\text{th}})-f(\eta)\hat{L}_{j}\right.\\ &+\left.\left[4\eta^{2}(1+\bar{n}_{\text{th}})+2\eta(1+2\bar{n}_{\text{th}})+\bar{n}_{\text{th}}\right]\hat{J}_{j}^{(+)}\right\}\hat{\varrho},\end{split}$$

(6)

where $G(\eta)=2(\gamma_{1}+\gamma_{2})[\bar{n}_{\text{th}}+2\eta(1+\bar{n}_{\text{th}})]$, and $f(\eta)=(4\eta+1)(1+\bar{n}_{\text{th}})+\bar{n}_{\text{th}}$. We can observe that imposing the condition

$$4\eta^{2}(1+\bar{n}_{\text{th}})+2\eta(1+2\bar{n}_{\text{th}})+\bar{n}_{\text{th}}=0,$$

(7)

we can eliminate the third term of the sum in (6). Condition (7) is fulfilled when

$$\eta=\begin{cases}\displaystyle\eta_{+}=-\frac{\bar{n}_{\text{th}}}{2(1+\bar{n}_{\text{th}})},\\[10.0pt] \displaystyle\eta_{-}=-\frac{1}{2},\end{cases}$$

(8)

such that $f(\eta_{\pm})=\pm 1$, and

$$G(\eta_{\pm})=\begin{cases}\displaystyle 0,&\mbox{if}\quad\eta=\eta_{+},\\[6.0pt] \displaystyle-2(\gamma_{1}+\gamma_{2}),&\mbox{if}\quad\eta=\eta_{-}.\end{cases}$$

(9)

Without loss of generality, we choose $\eta=\eta_{+}$, and then the Eq. (6) can be rewritten as

$$\frac{d\hat{\varrho}}{dt}=-\textrm{i}g\hat{S}\hat{\varrho}+\sum_{j=1}^{2}\gamma_{j}\left[\hat{J}_{j}^{(-)}(1+\bar{n}_{\text{th}})-\hat{L}_{j}\right]\hat{\varrho}.$$

(10)

Now, by performing the transformation

$$\hat{\varrho}=e^{-\chi\left[\hat{J}_{1}^{(-)}+\hat{J}_{2}^{(-)}\right]}\hat{\Omega},$$

(11)

we arrive at

$$\frac{d\hat{\Omega}}{dt}=-\textrm{i}g\hat{S}\hat{\Omega}+\sum_{j=1}^{2}\gamma_{j}\left[\hat{J}_{j}^{(-)}(1+\bar{n}_{\text{th}}-2\chi)-\hat{L}_{j}\right]\hat{\Omega}.$$

(12)

By setting

$$\chi=\frac{1+\bar{n}}{2},$$

(13)

the Eq. (12) reduces to a von Neumann-like equation

$$\frac{d\hat{\Omega}}{dt}=-\textrm{i}\left(\hat{H}_{\text{eff}}\hat{\Omega}-\hat{\Omega}\hat{H}_{\text{eff}}^{\dagger}\right),$$

(14)

for an effective non-Hermitian Hamiltonian

$$\hat{H}_{\text{eff}}=-\textrm{i}\gamma_{1}\hat{a}_{1}^{\dagger}\hat{a}_{1}-\textrm{i}\gamma_{2}\hat{a}_{2}^{\dagger}\hat{a}_{2}+g(\hat{a}_{1}\hat{a}_{2}^{\dagger}+\hat{a}_{1}^{\dagger}\hat{a}_{2}).$$

(15)

It is proposed that

$$\hat{\Omega}=|\psi\rangle\langle\psi|,$$

(16)

with

$$\textrm{i}\frac{d}{dt}|\psi\rangle=\hat{H}_{\text{eff}}|\psi\rangle\quad\text{and}\quad-\textrm{i}\frac{d}{dt}\langle\psi|=\langle\psi|\hat{H}_{\text{eff}}^{\dagger},$$

(17)

a pair of Schrödinger-type equations. Thus, the problem has been reduced to solve the equations in (17), for the non-Hermitian Hamiltonian (15).

From the results obtained [Eqs. (14) to (17)], we emphasize that the two non-unitary transformations (5) and (11) implemented in this work allow us to rewrite the Lindblad master equation for the interaction of two fields with decay at finite temperature as a von Neumann-like equation with an effective non-Hermitian Hamiltonian. Notably, this Hamiltonian is the same as that found in the case of two quantized fields with decay at zero temperature [16]. However, the evolution of the system in this new regime experiences passive dissipation processes and thermal excitations from the environment. Despite these differences, we have demonstrated that a direct equivalence can be established between both approaches through this pair of transformations. This connection highlights the importance of non-unitary transformations in understanding thermal dynamics, representing one of the main contributions of this work.

In order to diagonalize the effective non-Hermitian Hamiltonian (15), the following operators are introduced

$$\begin{split}\hat{N}&=\hat{a}_{1}^{\dagger}\hat{a}_{1}+\hat{a}_{2}^{\dagger}\hat{a}_{2},\\ \hat{J}_{x}&=\frac{1}{2}\left(\hat{a}_{1}\hat{a}_{2}^{\dagger}+\hat{a}_{1}^{\dagger}\hat{a}_{2}\right),\\ \hat{J}_{y}&=\frac{\textrm{i}}{2}\left(\hat{a}_{1}^{\dagger}\hat{a}_{2}-\hat{a}_{1}\hat{a}_{2}^{\dagger}\right),\\ \hat{J}_{z}&=\frac{1}{2}\left(\hat{a}_{2}^{\dagger}\hat{a}_{2}-\hat{a}_{1}^{\dagger}\hat{a}_{1}\right),\end{split}$$

(18)

satisfying the corresponding commutation relations $[\hat{N},\hat{J}_{k}]$, with $k=x,y,z$, and $[\hat{J}_{k},\hat{J}_{l}]=\textrm{i}\epsilon_{klm}\hat{J}_{m}$. The Hamiltonian (15) it is thus expressed as

$$\hat{H}_{\text{eff}}=-\textrm{i}\left[\frac{\gamma}{2}\hat{N}+\Delta\hat{J}_{z}\right]+2g\hat{J}_{x},$$

(19)

where $\gamma=\gamma_{1}+\gamma_{2}$, and $\Delta=\gamma_{2}-\gamma_{1}$. Now, we make use of the non-unitary transformation $\mathcal{\hat{R}}=e^{\xi\hat{J}_{y}}$, with $\xi$ a constant (in general complex) to be determined. Thus, we obtain

$$\hat{\mathcal{H}}=\mathcal{\hat{R}}^{-1}\hat{H}_{\textrm{eff}}\mathcal{\hat{R}}=c_{N}\hat{N}+c_{x}\hat{J}_{x}+c_{z}\hat{J}_{z},$$

(20)

where $c_{N}=-\textrm{i}\frac{\gamma}{2}$, $c_{x}=2g\cosh(\xi)-\Delta\sinh(\xi)$, and $c_{z}=\textrm{i}\left[2g\sinh(\xi)-\Delta\cosh(\xi)\right]$. By setting

$$\tanh(\xi)=\frac{2g}{\Delta},\quad\xi=\zeta+\textrm{i}\frac{\tau\pi}{2},\quad\zeta\in\mathbb{R},\quad\tau=0,1,\dots,$$

(21)

the operator $\hat{\mathcal{H}}$ diagonalizes as follows

$$\hat{\mathcal{H}}=\begin{cases}\displaystyle-\textrm{i}\left[\frac{\gamma}{2}\hat{N}\pm\sqrt{\Delta^{2}-4g^{2}}\hat{J}_{z}\right],\quad\mbox{if}&|\Delta|>2g,\\[8.0pt] \displaystyle-\textrm{i}\frac{\gamma}{2}\hat{N}\mp\sqrt{4g^{2}-\Delta^{2}}\hat{J}_{z},\qquad\;\mbox{if}&|\Delta|<2g.\end{cases}$$

(22)

Finally, from equations (5), (11), (16), (17) and the first equality in (20) the exact solution of the Lindblad master equation (1), given an initial condition $\hat{\rho}(0)$, is

$$\begin{split}\hat{\rho}(t)&=e^{-\eta_{+}\left[\hat{J}_{1}^{(+)}+\hat{J}_{2}^{(+)}\right]}e^{-\chi\left[\hat{J}_{1}^{(-)}+\hat{J}_{2}^{(-)}\right]}\hat{U}(t)\\ &\quad\times\left\{e^{\chi\left[\hat{J}_{1}^{(-)}+\hat{J}_{2}^{(-)}\right]}e^{\eta_{+}\left[\hat{J}_{1}^{(+)}+\hat{J}_{2}^{(+)}\right]}\hat{\rho}(0)\right\}\hat{U}^{\dagger}(t),\end{split}$$

(23)

where $\chi$ is given in (13), and

$$\hat{U}(t)=\exp\left[\xi\hat{J}_{y}\right]\exp\left[-\mathrm{i}\hat{\mathcal{H}}t\right]\exp\left[-\xi\hat{J}_{y}\right],$$

(24)

is the non-unitary evolution operator associated with $\hat{H}_{\textrm{eff}}$ and $\xi$ is defined by (21).

This study explores the Markovian dynamics of a quantum system characterized by the interaction of two quantized fields subjected to finite temperature decay. By employing superoperator techniques, we have reformulated the Lindblad master equation into a von Neumann-like equation with an effective non-Hermitian Hamiltonian, enabling the exact calculation of the time evolution of arbitrary initial states within a fully quantum framework. Additionally, our analysis of photon coincidence rates offers significant insights into the behavior of indistinguishable photons in a lossy cavity. As the temperature increases, the system transitions from a regime dominated by quantum interference to one increasingly influenced by thermal noise. This transition is reflected in the reduced depth of the quantum interference minimum and the increased coincidence rates, highlighting how thermal conditions impact quantum coherence and the dynamics of the system.

L. Hernández-Sánchez acknowledges the Instituto Nacional de Astrofísica, Óptica y Electrónica (INAOE) for the collaboration scholarship granted and the Consejo Nacional de Humanidades, Ciencias y Tecnologías (CONAHCYT) for the SNI Level III assistantship (CVU No. 736710). In turn, I. A. B.-G. acknowledges the Spanish MCIN with funding from European Union Next Generation EU (PRTRC17.I1) and Consejería de Educación from Junta de Castilla y León through QCAYLE project, as well as grants PID2020-113406GB-I00 MTM funded by AEI/10.13039/501100011033, and RED2022-134301-T, also CONAHCyT (México) for financial support through project A1-S-24569, and to IPN (México) for supplementary economical support through the project SIP20232237.