A Husimi Phase-Space Approach to a Driven–Dissipative Quantum Field at Finite Temperature

1. Introduction

The interaction between quantum systems and their environment is a central topic in quantum physics, since realistic systems cannot be completely isolated. Coupling to the surroundings induces dissipative dynamics, leading to decoherence and thermalization, and thus reducing quantum features such as coherence and entanglement [1,2]. Among open quantum systems, driven-dissipative quantum fields play a prominent role. In these systems, a quantized field is simultaneously subjected to coherent external driving and to dissipation arising from its coupling to a thermal reservoir [3,4,5], leading to a nontrivial interplay between coherent and incoherent dynamics [6]. The theoretical description of such systems is commonly based on the Lindblad master equation [7]. Analytical solutions for a quantized harmonic oscillator coupled to a thermal reservoir have been extensively studied in the literature and different approaches have been used, such as proposing an ansatz to solve the master equation [8], solving the differential equations that govern the evolution of the elements of the density matrix [9], using functional integrals [10], among others [11,12,13]. Superoperator formalisms have also been widely employed [14]. For example, in [15] the master equation is solved within this formalism, producing results that can be directly applied to the description of a single-mode quantized light field interacting with a squeezed thermal reservoir. In addition, the inclusion of a linear driving term in this context has also been investigated [16].

Although exact analytical solutions can be obtained for driven dissipative oscillators [3,4,5,17], they often involve cumbersome calculations or infinite series expansions, which may obscure the physical interpretation of the underlying dynamics. Alternative approaches leading to more explicit expressions have also been proposed [6]. These considerations motivate the search for formulations that provide a more transparent description of the system.

A particularly powerful approach to the study of dissipative quantum fields is provided by phase-space methods, in which the density operator is represented by a quasiprobability distribution function. For example, the dynamics of a cavity field coupled to a thermal bath has been described in terms of a Fokker–Planck equation for the Husimi Q-function [18], while in [19] the evolution of a cavity field is studied using the coherent-state expansion of the density operator, namely the Glauber–Sudarshan P-representation. More recently, [20] extended this framework to analyze the dynamics of a cavity field in a squeezed thermal environment.

Such models provide a useful framework for analyzing the statistical properties of the field [8], as well as for studying the effects of dissipation on phase-space interference [2] and for investigating the decay of quantum superpositions of displaced states, such as Schrödinger cat states and displaced thermal states, under interaction with a thermal reservoir [18].

In this work, we investigate a driven quantum field coupled to a finite-temperature reservoir. In Section 2, we introduce the physical model under study and analytically solve the corresponding master equation using superoperator techniques. Even for a coherent initial state, the solution is expressed in terms of infinite series. To overcome this difficulty, in Section 3 we adopt a phase-space description based on the Husimi Q-function, obtained by projecting the density operator onto coherent states. Within this framework, our results extend previous studies by explicitly incorporating the driving term into the analysis and providing a physically transparent phase-space description of the dynamics. In particular, we derive an analytical expression for the Husimi function of an initial superposition of coherent states and examine how its Gaussian and interference contributions evolve under driving, dissipation, and thermal fluctuations. Furthermore, we derive the corresponding Fokker–Planck equation governing the evolution of the Husimi function. This formulation allows for a comprehensive characterization of the field. In Section 4, we analyze the photon statistics of the system through the mean photon number, the standard deviation of the photon number, and the Mandel parameter, while in Section 5 we study its entropic properties by means of the Wehrl and linear entropies. Finally, the main results and concluding remarks are presented in Section 6.

2. Physical Model and Master Equation

The master equation for the density operator $\widehat{\tilde{\rho }}$ of a quantized field driven by a classical field coupled to a thermal reservoir at finite temperature is given by [4,21]

$${\displaystyle \frac{d\widehat{\tilde{\rho }}}{dt}}=-i\left[\widehat{H},\widehat{\tilde{\rho }}\right]+\gamma (n_{\mathrm{th}}+1)\mathcal{L}\left[\widehat{a}\right]\widehat{\tilde{\rho }}+\gamma n_{\mathrm{th}}\mathcal{L}\left[{\widehat{a}}^{†}\right]\widehat{\tilde{\rho }},$$

(1)

where $\widehat{H}=\omega \widehat{n}+\xi \left(\widehat{a}{e}^{i\omega t}+{\widehat{a}}^{†}{e}^{-i\omega t}\right)$ is the Hamiltonian that describes the interaction between the fields with $\widehat{n}={\widehat{a}}^{†}\widehat{a}$ the number operator, $\widehat{a}$ and ${\widehat{a}}^{†}$ the annihilation and creation operators, $\omega$ is the frequency of the quantized field and $\xi$ represents the strength of the coupling to the external source oscillating resonantly at frequency $\omega$. The dissipative part is described by Lindblad superoperators $\mathcal{L}\left[\widehat{c}\right]\widehat{\rho }=2\widehat{c}\widehat{\rho }{\widehat{c}}^{†}-({\widehat{c}}^{†}\widehat{c}\widehat{\rho }+\widehat{\rho }{\widehat{c}}^{†}\widehat{c})$, with $\widehat{c}=\widehat{a},{\widehat{a}}^{†}$, where $\gamma$ is the dissipation rate and $n_{\mathrm{th}}$ is the average number of thermal excitations in the reservoir, determined by its temperature. We solve the master equation using a superoperator approach. In the interaction picture, defined by the transformation

$$\widehat{\tilde{\rho }}=e^{-i\omega \widehat{n}t}\widehat{\rho }{e}^{i\omega \widehat{n}t},$$

(2)

the time-dependent part of the Hamiltonian $\widehat{H}$ and the free Hamiltonian of the quantized field are removed from the equation of motion, leading to

$$\begin{array}{cc}\hfill {\displaystyle \frac{d\widehat{\rho }}{dt}}& =-i\xi \left[\widehat{a}+{\widehat{a}}^{†},\widehat{\rho }\right]+\gamma (n_{\mathrm{th}}+1)\mathcal{L}\left[\widehat{a}\right]\widehat{\rho }+\gamma n_{\mathrm{th}}\mathcal{L}\left[{\widehat{a}}^{†}\right]\widehat{\rho }\hfill \\ & =\left[-i\xi {\widehat{S}}_I+\mu {\widehat{J}}_1+\beta {\widehat{J}}_2-(\mu +\beta )\widehat{L}-2\beta \right]\widehat{\rho },\hfill \end{array}$$

(3)

where $\mu =\gamma (n_{\mathrm{th}}+1)$ and $\beta =\gamma n_{\mathrm{th}}$, and where we have defined the following superoperators

$$\begin{array}{cccc}\hfill {\widehat{S}}_I\widehat{\rho }& =(\widehat{a}+{\widehat{a}}^{†})\widehat{\rho }-\widehat{\rho }(\widehat{a}+{\widehat{a}}^{†}),\hfill & \hfill \widehat{L}\widehat{\rho }& =\widehat{n}\widehat{\rho }+\widehat{\rho }\widehat{n},\hfill \\ \hfill {\widehat{J}}_1\widehat{\rho }& =2\widehat{a}\widehat{\rho }{\widehat{a}}^{†},\hfill & \hfill {\widehat{J}}_2\widehat{\rho }& =2{\widehat{a}}^{†}\widehat{\rho }\widehat{a}.\hfill \end{array}$$

(4)

Now, a displacement transformation can be applied to eliminate the term $(\widehat{a}+{\widehat{a}}^{†})$. Let us write

$$\widehat{\rho }=\widehat{D}\left(\eta \right)\widehat{R}{\widehat{D}}^{†}\left(\eta \right),$$

(5)

where $\widehat{D}\left(\eta \right)=exp[\eta {\widehat{a}}^{†}-{\eta }^{*}\widehat{a}]$ is the Glauber displacement operator [22], with $\eta =-i{\displaystyle \frac{\xi }{\mu -\beta }}=-i{\displaystyle \frac{\xi }{\gamma }}$. As a result, we obtain

$${\displaystyle \frac{d\widehat{R}}{dt}}=\left[\mu {\widehat{J}}_1+\beta {\widehat{J}}_2-(\mu +\beta )\widehat{L}-2\beta \right]\widehat{R}.$$

(6)

This equation can be identified as the master equation for the density matrix of a single-mode quantum field in a lossy cavity at finite temperature [14]. The solution to this equation is

$$\widehat{R}\left(t\right)=e^{-2\beta t}{e}^{\left[\mu {\widehat{J}}_1+\beta {\widehat{J}}_2-(\mu +\beta )\widehat{L}\right]t}\widehat{R}\left(0\right).$$

(7)

We note that the superoperators satisfy the following commutation relations:

$$\left[{\widehat{J}}_2,{\widehat{J}}_1\right]\widehat{\rho }=-4(\widehat{L}+1)\widehat{\rho },\left[{\widehat{J}}_1,\widehat{L}\right]\widehat{\rho }=2{\widehat{J}}_1\widehat{\rho },\left[{\widehat{J}}_2,\widehat{L}\right]\widehat{\rho }=-2{\widehat{J}}_2\widehat{\rho }.$$

(8)

In order to factorize the exponential operator in Eq. (7), we propose the following ansatz:

$$e^{\left[\mu {\widehat{J}}_1+\beta {\widehat{J}}_2-(\mu +\beta )\widehat{L}\right]t}=e^{f_3\left(t\right)}{e}^{f_2\left(t\right){\widehat{J}}_2}{e}^{f_0\left(t\right)\widehat{L}}{e}^{f_1\left(t\right){\widehat{J}}_1}.$$

(9)

Taking the time derivative of both sides of Eq. (9), and using Hadamard’s lemma [23] together with the commutation relations in Eq. (8) to evaluate similarity transformations of superoperators, we compare terms with identical superoperators and obtain a system of differential equations for functions f, with initial conditions $f_i\left(0\right)=0$ for $i=0,1,2,3$. By solving this system of equations, we obtain [14]

$$\begin{array}{cc}\hfill f_0\left(t\right)& =(\beta -\mu )t-ln\left({\displaystyle \frac{\beta e^{2t(\beta -\mu )}-\mu }{\beta -\mu }}\right)=-\gamma t-ln(1+n_{\mathrm{th}}(1-e^{-2\gamma t})),\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill f_1\left(t\right)& ={\displaystyle \frac{\mu }{2}}{\displaystyle \frac{e^{2t(\beta -\mu )}-1}{\beta e^{2t(\beta -\mu )}-\mu }}={\displaystyle \frac{n_{\mathrm{th}}+1}{2}}{\displaystyle \frac{1-e^{-2\gamma t}}{1+n_{\mathrm{th}}(1-e^{-2\gamma t})}},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill f_2\left(t\right)& ={\displaystyle \frac{\beta }{\mu }}{f}_1\left(t\right)={\displaystyle \frac{n_{\mathrm{th}}}{n_{\mathrm{th}}+1}}{f}_1\left(t\right),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill f_3\left(t\right)& =f_0\left(t\right)+\gamma (2n_{\mathrm{th}}+1)t.\hfill \end{array}$$

(13)

Therefore, in the interaction picture, the density operator will be given by

$$\widehat{\rho }\left(t\right)=\widehat{D}\left(\eta \right)\left[e^{-2\beta t+f_3\left(t\right)}{e}^{f_2\left(t\right){\widehat{J}}_2}{e}^{f_0\left(t\right)\widehat{L}}{e}^{f_1\left(t\right){\widehat{J}}_1}\left({\widehat{D}}^{†}\left(\eta \right)\widehat{\rho }\left(0\right)\widehat{D}\left(\eta \right)\right)\right]{\widehat{D}}^{†}\left(\eta \right).$$

(14)

Assuming that the initial state of the field is a coherent state with amplitude $\mathsf{\Gamma }={\mathsf{\Gamma }}_x+i{\mathsf{\Gamma }}_y$, i.e., $\widehat{\rho }\left(0\right)=|\mathsf{\Gamma }\rangle \langle \mathsf{\Gamma }|$, the action of the Glauber displacement operator and of the exponentials involving the superoperators ${\widehat{J}}_1$ and $\widehat{L}$ can be evaluated straightforwardly, yielding

$$e^{f_0\left(t\right)\widehat{L}}{e}^{f_1\left(t\right){\widehat{J}}_1}\left({\widehat{D}}^{†}\left(\eta \right)\widehat{\rho }\left(0\right)\widehat{D}\left(\eta \right)\right)=e^{{-|\mathsf{\Gamma }-\eta |}^2(1-|e^{f_0}{|}^2-2f_1)}|(\mathsf{\Gamma }-\eta )e^{f_0}\rangle \langle (\mathsf{\Gamma }-\eta )e^{f_0}|,$$

(15)

where $|(\mathsf{\Gamma }-\eta )e^{f_0}\rangle$ denotes a coherent state. However, when $exp\left[f_2\left(t\right){\widehat{J}}_2\right]$ is applied to the result in Eq. (15), we obtain an infinite series, given by

$$\begin{array}{cc}\hfill e^{f_2\left(t\right){\widehat{J}}_2}{e}^{f_0\left(t\right)\widehat{L}}{e}^{f_1\left(t\right){\widehat{J}}_1}& \left({\widehat{D}}^{†}\left(\eta \right)\widehat{\rho }\left(0\right)\widehat{D}\left(\eta \right)\right)\hfill \\ & =e^{{-|\mathsf{\Gamma }-\eta |}^2(1-|e^{f_0}{|}^2-2f_1)}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{{\left(2f_2\right)}^k}{k!}}{\widehat{a}}^{†k}|(\mathsf{\Gamma }-\eta )e^{f_0}\rangle \langle (\mathsf{\Gamma }-\eta )e^{f_0}|{\widehat{a}}^k.\hfill \end{array}$$

(16)

Note that applying the creation operator to a coherent state does not lead to a simple closed-form expression for the field’s density operator; rather, it results in an infinite series that is not convenient for further analytical treatment. For this reason, we adopt a Husimi-distribution approach: by projecting the resulting expression onto coherent states, the series can be handled more efficiently, allowing us to obtain the associated Husimi function of the field. This representation yields a more tractable analytical expression, facilitating the computation of photon statistics and other relevant physical quantities.

3. Phase-Space Representation

3.1. Husimi Function

The Husimi function of the quantized field is given by the expectation value of its density matrix with respect to the coherent state $|\alpha \rangle$ [21], that is,

$$Q(\alpha ,t)={\displaystyle \frac{1}{\pi }}\langle \alpha |\widehat{\rho }|\alpha \rangle .$$

(17)

In our case, we have

$$Q(\alpha ,t)={\displaystyle \frac{1}{\pi }}\langle \alpha |\widehat{D}\left(\eta \right)\left[e^{-2\beta t+f_3\left(t\right)}{e}^{f_2\left(t\right){\widehat{J}}_2}{e}^{f_0\left(t\right)\widehat{L}}{e}^{f_1\left(t\right){\widehat{J}}_1}\left({\widehat{D}}^{†}\left(\eta \right)\widehat{\rho }\left(0\right)\widehat{D}\left(\eta \right)\right)\right]{\widehat{D}}^{†}\left(\eta \right)|\alpha \rangle .$$

(18)

Using the result in (16) and the fact that ${\widehat{D}}^{†}\left(\eta \right)|\alpha \rangle =e^{-iIm\left\{\eta {\alpha }^{*}\right\}}|\alpha -\eta \rangle$, the Husimi function (18) can be expressed as

$$Q(\alpha ,t)={\displaystyle \frac{1}{\pi }}{e}^{-2\beta t+f_3}{e}^{{-|\mathsf{\Gamma }-\eta |}^2(1-|e^{f_0}{|}^2-2f_1)}{\displaystyle \sum _{k=0}^{\infty }}{\displaystyle \frac{{\left(2f_2\right)}^k}{k!}}\langle \alpha -\eta |{\widehat{a}}^{†k}|(\mathsf{\Gamma }-\eta )e^{f_0}\rangle \langle (\mathsf{\Gamma }-\eta )e^{f_0}|{\widehat{a}}^k|\alpha -\eta \rangle .$$

(19)

Then, after carrying out the corresponding calculations, we arrive at a more compact expression for the Husimi function, namely

$$\begin{array}{cc}\hfill Q(\alpha ,t)=& {\displaystyle \frac{1}{\pi }}exp[-2\beta t+f_3]exp\left[{-|\mathsf{\Gamma }-\eta |}^2(1-2f_1)\right]\hfill \\ & \times exp\left[{-|\alpha -\eta |}^2(1-2f_2)\right]exp[2e^{f_0}\mathrm{Re}\left\{(\alpha -\eta )({\mathsf{\Gamma }}^{*}-{\eta }^{*})\}\right].\hfill \end{array}$$

(20)

Alternatively, the Husimi function can be written in compact Gaussian form, with $\alpha ={\alpha }_R+i{\alpha }_I$, and expressed in terms of the relevant physical parameters of the system, such as

$$Q(\alpha ,t)={\displaystyle \frac{\theta }{\pi }}{e}^{-\theta |\alpha -{\alpha }_0{\left(t\right)|}^2}={\displaystyle \frac{\theta }{\pi }}{e}^{-\theta {({\alpha }_R-{\alpha }_{0R}\left(t\right))}^2}{e}^{-\theta {({\alpha }_I-{\alpha }_{0I}\left(t\right))}^2},$$

(21)

where

$$\begin{array}{cc}\hfill \theta \left(t\right)=& 1-2f_2\left(t\right)={\displaystyle \frac{1}{1+n_{\mathrm{th}}(1-e^{-2\gamma t})}},\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\alpha }_0\left(t\right)=& {\alpha }_{0R}\left(t\right)+i{\alpha }_{0I}\left(t\right),\hfill \end{array}$$

(23)

with

$$\begin{array}{c}\hfill {\alpha }_{0R}\left(t\right)={\mathsf{\Gamma }}_x{e}^{-\gamma t},{\alpha }_{0I}\left(t\right)=-{\displaystyle \frac{\xi }{\gamma }}+\left({\displaystyle \frac{\xi }{\gamma }}+{\mathsf{\Gamma }}_y\right)e^{-\gamma t}.\end{array}$$

(24)

The expression obtained for the Husimi function shows that the state of the system remains Gaussian throughout its evolution. The Q-function provides a convenient way to evaluate expectation values and moments of antinormal ordered operators [21,24], that is,

$$\langle {\widehat{\mathcal{O}}}_A(\widehat{a},{\widehat{a}}^{†})\rangle =\int Q(\alpha ,{\alpha }^{*}){\mathcal{O}}_A(\alpha ,{\alpha }^{*})d^2\alpha ,$$

(25)

where ${\widehat{\mathcal{O}}}_A(\widehat{a},{\widehat{a}}^{†})$ denotes an operator written in antinormal order. Using this relation, we can directly evaluate the first moment of the field. In particular, for the annihilation operator, we obtain

$$\langle \widehat{a}\rangle ={\alpha }_0\left(t\right).$$

(26)

Therefore, the evolution of the Husimi function, Eq. (21), is fully characterized by the evolution of the first moment, encoded in the displacement ${\alpha }_0\left(t\right)$, and by the width parameter $\theta$, which accounts for the effect of the thermal reservoir. In particular, ${\alpha }_0\left(t\right)$ describes the trajectory of the center of the Husimi function in phase space, which in the stationary regime approaches ${\alpha }_0=-i\xi /\gamma$, while $\theta$ reflects the competition between dissipation and thermal fluctuations, leading to a broadening of the distribution as the system approaches equilibrium. At short times, $e^{-2\gamma t}\approx 1$, so that $\theta \approx 1$, and the state remains close to its initial coherent state, i.e., with minimal broadening. As time evolves, dissipative dynamics reduces the coherent amplitude, while thermal noise from the reservoir progressively increases the fluctuations, leading to a decrease of $\theta$ and hence to a broadening of the Husimi distribution whenever $n_{\mathrm{th}}\ne 0$. In contrast, for $n_{\mathrm{th}}=0$, we have $\theta =1$ at all times, indicating that an initial coherent state remains pure during evolution. In the long-time limit, $\theta \to {(1+n_{\mathrm{th}})}^{-1}$, showing that the system reaches a stationary Gaussian state whose width is entirely determined by the reservoir temperature; that is, the field relaxes to a displaced thermal state described by

$$Q(\alpha ,t\to \infty )={\displaystyle \frac{1}{\pi }}\langle \alpha |\widehat{D}\left({\alpha }_0\right){\widehat{\rho }}_{\mathrm{th}}{\widehat{D}}^{†}\left({\alpha }_0\right)|\alpha \rangle .$$

(27)

Figure 1 illustrates this behavior by showing the evolution of the Husimi function at scaled times $\xi t=0,\kappa /4,\kappa /2,3\kappa /4,\kappa$, as well as in the asymptotic limit $\xi t\to \infty$. Here, $\kappa =\xi /\gamma$ represents the ratio between the driving strength and the dissipation rate, and the choice $\xi t=\kappa$ corresponds to a characteristic timescale at which the exponential terms in Eq. (24) decay to $e^{-1}$. The initial field state is taken to be a coherent state with amplitude $\mathsf{\Gamma }=1+i$. In this analysis, we consider $\kappa =1$ and the mean thermal excitation number $n_{\mathrm{th}}=1$.

3.2. Interference Dynamics

Having established the evolution of an initial coherent state, we now consider a nonclassical initial condition given by a superposition of two coherent states. The purpose of this analysis is to investigate how coherent driving, dissipation, and thermal fluctuations affect the interference features associated with quantum superpositions in phase space.

The initial state of the field is taken as

$$\begin{array}{cc}\hfill \left|\psi \right(0)\rangle & =C_N\left(|{\mathsf{\Gamma }}_1\rangle +|{\mathsf{\Gamma }}_2\rangle \right),\hfill \\ \hfill C_N& ={\left[2\left(1+e^{-{\displaystyle \frac{1}{2}}({\left|{\mathsf{\Gamma }}_1\right|}^2+|{\mathsf{\Gamma }}_2{|}^2)}Re\left\{e^{{\mathsf{\Gamma }}_1{\mathsf{\Gamma }}_2^{*}}\right\}\right)\right]}^{-1/2},\hfill \end{array}$$

(28)

where $|{\mathsf{\Gamma }}_j\rangle$ are coherent states with amplitudes ${\mathsf{\Gamma }}_j$, and $C_N$ is the normalization constant. Such states exhibit interference effects in phase space arising from the coherent superposition of their components.

Following the same procedure developed in the previous sections, the analytical solution of the master equation can be projected onto coherent states, yielding the corresponding Husimi function

$$\begin{array}{cc}\hfill & Q\left(\alpha \right)=Q_1\left(\alpha \right)+Q_2\left(\alpha \right)+Q_{\mathrm{int}}\left(\alpha \right)={\left|C_N\right|}^2{\displaystyle \frac{\theta }{\pi }}[e^{-\theta |\alpha -d_1{|}^2}+e^{-\theta |\alpha -d_2{|}^2}\hfill \\ \hfill & +2e^{-{\displaystyle \frac{1}{2}}\left(|{\mathsf{\Gamma }}_1{|}^2+{\left|{\mathsf{\Gamma }}_2\right|}^2\right)}Re\left\{e^{{\mathsf{\Gamma }}_1{\mathsf{\Gamma }}_2^{*}}{e}^{-\theta {\left({\alpha }_R-{\displaystyle \frac{1}{2}}({\mathsf{\Gamma }}_1+{\mathsf{\Gamma }}_2^{*})e^{-\gamma t}\right)}^2}{e}^{-\theta {\left({\alpha }_I+{\displaystyle \frac{i}{2}}({\mathsf{\Gamma }}_1-{\mathsf{\Gamma }}_2^{*})e^{-\gamma t}+{\displaystyle \frac{\xi }{\gamma }}\left(1-e^{-\gamma t}\right)\right)}^2}\right\}].\hfill \end{array}$$

(29)

The quantities $d_j\left(t\right)$ denote the centers of the Gaussian contributions,

$$\begin{array}{c}\hfill d_j\left(t\right)={\alpha }_{0Rj}\left(t\right)+i{\alpha }_{0Ij}\left(t\right),\end{array}$$

(30)

with

$$\begin{array}{c}\hfill {\alpha }_{0Rj}\left(t\right)={\mathsf{\Gamma }}_{jx}{e}^{-\gamma t},{\alpha }_{0Ij}\left(t\right)=-{\displaystyle \frac{\xi }{\gamma }}+\left({\displaystyle \frac{\xi }{\gamma }}+{\mathsf{\Gamma }}_{jy}\right)e^{-\gamma t}.\end{array}$$

(31)

Equation (29) is naturally separated into three contributions. The first two terms, $Q_1\left(\alpha \right)$ and $Q_2\left(\alpha \right)$, correspond to the Gaussian distributions associated with the coherent states $|{\mathsf{\Gamma }}_1\rangle$ and $|{\mathsf{\Gamma }}_2\rangle$. Their centers evolve according to $d_j\left(t\right)$ and approach the stationary value $-i\xi /\gamma$, indicating that coherent driving displaces the phase-space distribution toward a shifted equilibrium position. At the same time, the separation between the two Gaussian peaks decreases as

$$|d_1\left(t\right)-d_2\left(t\right)|=|{\mathsf{\Gamma }}_1-{\mathsf{\Gamma }}_2|e^{-\gamma t},$$

(32)

demonstrating that dissipation progressively reduces the distinguishability of the two coherent components. In addition, thermal fluctuations broaden the distribution through the parameter $\theta$, further reducing the resolution of the phase-space structure.

The third term, $Q_{\mathrm{int}}\left(\alpha \right)$, contains the interference contribution and, therefore, encodes the coherence initially present in the superposition. From Eq. (29), it can be seen that the dependence on the coherent amplitudes ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$ appears through terms multiplied by $e^{-\gamma t}$. Consequently, the influence of the initial coherent components on the interference pattern progressively decreases with time, reflecting the gradual loss of information about the initial superposition induced by dissipation. At the same time, thermal fluctuations broaden the distribution through the parameter $\theta$, further reducing the visibility of the interference features. On the other hand, the driving force $\xi$ enters through the combination ${\displaystyle \frac{\xi }{\gamma }}(1-e^{-\gamma t})$, which corresponds to a displacement of the interference structure in phase space. Therefore, dissipation and thermal noise are responsible for washing out the interference features, whereas the driving field mainly shifts their position and does not contribute to their suppression.

It is important to emphasize that the persistence of the interference term in Eq. (29) does not imply the preservation of quantum coherence. As the system evolves, the dependence of the Husimi function on the initial coherent amplitudes is progressively erased, and the long-time dynamics becomes identical to that obtained from a single coherent-state initial condition. Consequently, all signatures of the initial superposition disappear and the system relaxes to a unique displaced thermal state independent of ${\mathsf{\Gamma }}_1$ and ${\mathsf{\Gamma }}_2$. Figure 2 provides a direct visualization of the behavior described by Eq. (29) for an initial even Schrödinger cat state $|\psi \left(0\right)\rangle =C_N\left(|{\mathsf{\Gamma }}_1\rangle +|-{\mathsf{\Gamma }}_1\rangle \right)$ with ${\mathsf{\Gamma }}_1={\mathsf{\Gamma }}_2=1+i$. The upper row shows the evolution of the complete Husimi function, where the two coherent components, $Q_1$ and $Q_2$, progressively approach each other and eventually merge into a single stationary distribution. The lower row displays the corresponding interference contribution, $Q_{\mathrm{int}}$, making explicit the gradual transformation of the interference structure.

3.3. Fokker–Planck Equation

Since we already know how the Husimi function evolves in time, it is natural to expect that it satisfies a corresponding Fokker–Planck equation. This equation can be derived starting from the Lindblad master equation in the interaction picture, Eq. (3), by taking its expectation value in coherent states, i.e., by projecting onto the coherent state basis as follows:

$$Tr\{|\alpha \rangle \langle \alpha |{\displaystyle \frac{d\widehat{\rho }}{dt}}\}=Tr\left\{\right|\alpha \rangle \langle \alpha |\left(-i\xi \left[\widehat{a}+{\widehat{a}}^{†},\widehat{\rho }\right]+\gamma (n_{\mathrm{th}}+1)\mathcal{L}\left[\widehat{a}\right]\widehat{\rho }+\gamma n_{\mathrm{th}}\mathcal{L}\left[{\widehat{a}}^{†}\right]\widehat{\rho }\right)\},$$

(33)

where $Tr\left\{\right|\alpha \rangle \langle \alpha \left|\widehat{\rho }\right\}=\langle \alpha |\widehat{\rho }|\alpha \rangle =\pi Q\left(\alpha \right)$. Taking into account the relation $|\alpha \rangle \langle \alpha |\widehat{a}=\left(\alpha +{\displaystyle \frac{\partial }{\partial {\alpha }^{*}}}\right)|\alpha \rangle \langle \alpha |$ [23], we obtain the corresponding Fokker–Planck equation

$${\displaystyle \frac{dQ\left(\alpha \right)}{dt}}=\left[-i\xi \left({\displaystyle \frac{\partial }{\partial {\alpha }^{*}}}-{\displaystyle \frac{\partial }{\partial \alpha }}\right)+\gamma \left({\displaystyle \frac{\partial }{\partial \alpha }}\alpha +{\displaystyle \frac{\partial }{\partial {\alpha }^{*}}}{\alpha }^{*}\right)+2\gamma (n_{\mathrm{th}}+1){\displaystyle \frac{{\partial }^2}{\partial {\alpha }^{*}\partial \alpha }}\right]Q\left(\alpha \right).$$

(34)

Figure 2. Time evolution of the Husimi function and its interference contribution for an initial even Schrödinger cat state, $|\psi \left(0\right)\rangle =C_N\left(|{\mathsf{\Gamma }}_1\rangle +|-{\mathsf{\Gamma }}_1\rangle \right)$, with ${\mathsf{\Gamma }}_1=1+i$. The upper row shows the Husimi function, while the lower row displays the corresponding interference contribution given by the third term in Eq. (29). The panels correspond to different values of the scaled time $\xi t$, illustrating the progressive merging of the coherent components and the evolution of the interference contribution toward its stationary form. Results are shown for $\kappa =1$ and $n_{\mathrm{th}}=1$.

Figure 2. Time evolution of the Husimi function and its interference contribution for an initial even Schrödinger cat state, $|\psi \left(0\right)\rangle =C_N\left(|{\mathsf{\Gamma }}_1\rangle +|-{\mathsf{\Gamma }}_1\rangle \right)$, with ${\mathsf{\Gamma }}_1=1+i$. The upper row shows the Husimi function, while the lower row displays the corresponding interference contribution given by the third term in Eq. (29). The panels correspond to different values of the scaled time $\xi t$, illustrating the progressive merging of the coherent components and the evolution of the interference contribution toward its stationary form. Results are shown for $\kappa =1$ and $n_{\mathrm{th}}=1$.

The first term, proportional to $\xi$, corresponds to the coherent driving of the field and induces a displacement of the distribution in the phase space. The second term, proportional to $\gamma$, describes the dissipative dynamics and gives rise to a drift toward the stationary state. The last term, proportional to $\gamma (n_{\mathrm{th}}+1)$, represents diffusion induced by the thermal reservoir, leading to a broadening of the Husimi function.

4. Photon Statistics

The fact that an initially coherent state remains Gaussian throughout the evolution, as shown in Eq. (21), greatly simplifies its characterization and enables a straightforward evaluation of relevant physical quantities, such as photon statistics, directly from the Husimi function through Eq. (25), without the need to manipulate the density operator, since these quantities depend only on first- and second-order moments. Having discussed interference effects for superpositions of coherent states in the previous subsection, we now return to the coherent state initial condition and analyze the statistical properties of the field. In particular, we calculate the temporal evolution of the photon number, the photon number standard deviation, and the Mandel parameter to investigate the effects of driving and dissipation on the system.

4.1. Mean Photon Number

Using Eq. (25), the expectation value of the photon number operator is obtained as

$$\langle \widehat{n}\left(t\right)\rangle =\langle \widehat{a}{\widehat{a}}^{†}\rangle -1={\int }_{-\infty }^{\infty }Q(\alpha ,t){\left|\alpha \right|}^{2}d{\alpha }^2-1.$$

(35)

By performing the calculations, we obtain the expectation value of the photon number in terms of the physical parameters

$$\langle \widehat{n}\left(t\right)\rangle =n_{\mathrm{th}}+{\displaystyle \frac{{\xi }^2}{{\gamma }^2}}-2{\displaystyle \frac{\xi }{\gamma }}\left({\displaystyle \frac{\xi }{\gamma }}+{\mathsf{\Gamma }}_y\right)e^{-\gamma t}+\left({\left|\mathsf{\Gamma }+i{\displaystyle \frac{\xi }{\gamma }}\right|}^2-n_{\mathrm{th}}\right)e^{-2\gamma t}.$$

(36)

The first term, $n_{\mathrm{th}}$, accounts for the thermal contribution of the reservoir. The second term, ${\xi }^2/{\gamma }^2$, corresponds to the steady-state excitation induced by the driving. The terms proportional to $e^{-\gamma t}$ and $e^{-2\gamma t}$ describe transient dynamics associated with the decay of the initial coherent state and its relaxation toward the stationary regime. We note that the stationarity occurs when

$$\mathsf{\Gamma }=\pm \sqrt{n_{\mathrm{th}}}-i{\displaystyle \frac{\xi }{\gamma }}.$$

(37)

In the long-time limit, the exponential terms vanish and the system reaches a stationary value given by

$${\langle \widehat{n}\rangle }_{\mathrm{ss}}=n_{\mathrm{th}}+{\displaystyle \frac{{\xi }^2}{{\gamma }^2}},$$

(38)

which clearly shows the combined contribution of thermal fluctuations and coherent driving. Figure 4 shows the time evolution of the average photon number, as a function of the scaled time $\xi t$, for different initial coherent states of the system, including the stationary state, and for parameters $\kappa =1$ and $n_{\mathrm{th}}=1$. The colored curves correspond to the analytical results given by Eq. (36), while the dashed curve, in this and all subsequent figures, represents the numerical solution obtained by solving the master equation (3).

Figure 3. Average number of photons as a function of the scaled time $\xi t$ for different initial coherent states, with parameters $\kappa =1$ and $n_{\mathrm{th}}=1$, showing convergence toward a unique stationary value, independent of the initial condition.

Figure 3. Average number of photons as a function of the scaled time $\xi t$ for different initial coherent states, with parameters $\kappa =1$ and $n_{\mathrm{th}}=1$, showing convergence toward a unique stationary value, independent of the initial condition.

4.2. Photon-Number Standard Deviation

The fluctuations in the number of photons are characterized by the standard deviation,

$$\Delta \widehat{n}\left(t\right)=\sqrt{\langle {\widehat{n}}^2\rangle -{\langle \widehat{n}\rangle }^2},$$

(39)

where $\langle {\widehat{n}}^2\rangle$ can be calculated using equation (25). After some algebra, we obtain

$$\begin{array}{cc}\hfill \Delta \widehat{n}\left(t\right)=& [n_{\mathrm{th}}(1+n_{\mathrm{th}})+(1+2n_{\mathrm{th}}){\displaystyle \frac{{\xi }^2}{{\gamma }^2}}-2(1+2n_{\mathrm{th}}){\displaystyle \frac{\xi }{\gamma }}\left({\displaystyle \frac{\xi }{\gamma }}+{\mathsf{\Gamma }}_y\right)e^{-\gamma t}\hfill \\ \hfill & +\left(\left(|\mathsf{\Gamma }+i{\displaystyle \frac{\xi }{\gamma }}{|}^2-n_{\mathrm{th}}\right)(1+2n_{\mathrm{th}})-2n_{\mathrm{th}}{\displaystyle \frac{{\xi }^2}{{\gamma }^2}}\right)e^{-2\gamma t}\hfill \\ & +4n_{\mathrm{th}}{\displaystyle \frac{\xi }{\gamma }}\left({\displaystyle \frac{\xi }{\gamma }}+{\mathsf{\Gamma }}_y\right)e^{-3\gamma t}+n_{\mathrm{th}}\left(n_{\mathrm{th}}-2{|\mathsf{\Gamma }+i{\displaystyle \frac{\xi }{\gamma }}|}^2\right)e^{-4\gamma t}{]}^{{\displaystyle \frac{1}{2}}}.\hfill \end{array}$$

(40)

The standard deviation exhibits nontrivial transient dynamics characterized by multiple exponential decay rates, reflecting the interplay between the initial coherent state, coherent driving, dissipation, and thermal fluctuations. In the long-time limit, all transient contributions vanish, and the system approaches a stationary value given by

$$\Delta {\widehat{n}}_{\mathrm{ss}}={\left[n_{\mathrm{th}}(1+n_{\mathrm{th}})+(1+2n_{\mathrm{th}}){\displaystyle \frac{{\xi }^2}{{\gamma }^2}}\right]}^{1/2},$$

(41)

which shows that the steady-state fluctuations arise from both thermal noise and coherent driving. This behavior is illustrated in Figure 4, where the time evolution of the standard deviation of the photon number is displayed for different initial coherent states.

Figure 4. Standard deviation of the photon number as a function of the scaled time $\xi t$ for different initial coherent states, with $\kappa =1$ and $n_{\mathrm{th}}=1$.

Figure 4. Standard deviation of the photon number as a function of the scaled time $\xi t$ for different initial coherent states, with $\kappa =1$ and $n_{\mathrm{th}}=1$.

4.3. Mandel Parameter

To further characterize the statistical properties of the field, we consider the Mandel parameter Q (hereafter denoted by $Q_M$ to avoid confusion with the Husimi Q-function), which provides a measure of the deviation of the photon statistics from a Poissonian distribution. It is defined in terms of the first and second moments of the photon-number operator as

$$Q_M={\displaystyle \frac{\langle {\left(\Delta \widehat{n}\right)}^2\rangle -\langle \widehat{n}\rangle }{\langle \widehat{n}\rangle }}.$$

(42)

Negative values of $Q_M$ indicate sub-Poissonian statistics (nonclassical light), $Q_M=0$ corresponds to Poissonian statistics as in a coherent state, and $Q_M>0$ describes super-Poissonian behavior associated with enhanced fluctuations. Using the results obtained for the mean photon number and its variance, we find

$$\begin{array}{cc}\hfill Q_M\left(t\right)=& {\displaystyle \frac{n_{\mathrm{th}}}{\langle \widehat{n}\left(t\right)\rangle }}\left[\left(1-e^{-2\gamma t}\right)\left(n_{\mathrm{th}}(1-e^{-2\gamma t})+2\left({\left[{\displaystyle \frac{\xi }{\gamma }}-\left({\displaystyle \frac{\xi }{\gamma }}+{\mathsf{\Gamma }}_y\right)e^{-\gamma t}\right]}^2+{\mathsf{\Gamma }}_x^2{e}^{-2\gamma t}\right)\right)\right]\hfill \\ \hfill =& n_{\mathrm{th}}\left(1-e^{-2\gamma t}\right)\left(2-{\displaystyle \frac{n_{\mathrm{th}}\left(1-e^{-2\gamma t}\right)}{\langle \widehat{n}\left(t\right)\rangle }}\right).\hfill \end{array}$$

(43)

Several conclusions follow from the above expression. First, since $\langle \widehat{n}\rangle \ge 0$ and $(1-e^{-2\gamma t})\ge 0$ for $t\ge 0$, one has $Q_M\left(t\right)\ge 0$ at all times, indicating that the field exhibits only Poissonian or super-Poissonian statistics during evolution. In particular, for $t=0$ one obtains $Q_M\left(0\right)=0$, as expected for an initial coherent state. Moreover, when $n_{\mathrm{th}}=0$, the Mandel parameter remains zero at all times, showing that in the absence of thermal excitations, the field preserves Poissonian statistics. Therefore, the emergence of super-Poissonian behavior is entirely induced by the thermal reservoir. Figure 5 shows the time evolution of the Mandel parameter for different initial coherent amplitudes of the field. As expected for coherent states, all curves start from $Q_M=0$, corresponding to Poissonian photon statistics. As time evolves, the interaction with the thermal reservoir induces super-Poissonian behavior, leading to positive values of $Q_M$. Eventually, all trajectories converge toward the same stationary value.

Figure 5. Mandel parameter $Q_M$ as a function of the scaled time $\xi t$ for different initial coherent amplitudes of the field, with parameters $\kappa =1$ and $n_{\mathrm{th}}=1$.

Figure 5. Mandel parameter $Q_M$ as a function of the scaled time $\xi t$ for different initial coherent amplitudes of the field, with parameters $\kappa =1$ and $n_{\mathrm{th}}=1$.

We note that in the steady state, the Mandel parameter is given by

$$Q_{Mss}=n_{\mathrm{th}}{\displaystyle \frac{n_{\mathrm{th}}+2{\left(\frac{\xi }{\gamma }\right)}^2}{n_{\mathrm{th}}+{\left(\frac{\xi }{\gamma }\right)}^2}}.$$

(44)

Figure 6 illustrates the dependence of the steady-state Mandel parameter on the ratio $\kappa =\xi /\gamma$ and on the average thermal occupation number $n_{\mathrm{th}}$.

Figure 6. Steady-state Mandel parameter $Q_M$ as a function of $\kappa$ and $n_{\mathrm{th}}$.

Figure 6. Steady-state Mandel parameter $Q_M$ as a function of $\kappa$ and $n_{\mathrm{th}}$.

5. Entropic Measures

Entropy measures provide complementary information about the evolution of the field state beyond photon statistics. Due to the Gaussian form of the state, the known Wehrl entropy [25], an entropy of a quantum state in terms of coherent states, can be directly evaluated from the Husimi function, reinforcing the usefulness of the phase-space approach adopted here. We also calculate the linear entropy, which quantifies the loss of purity caused by the interaction with the thermal reservoir.

5.1. Wehrl Entropy

The Wehrl entropy [25] is an entropy defined over the Husimi function Q that captures phase-space uncertainties of a quantum state. Beyond its foundational role, it has proven useful in characterizing the degree of squeeze of squeezed states [26], and has been proposed as a measure of the nonclassicality of bosonic states [27,28]. It is defined as

$$S_W\left(\widehat{\rho }\right)=-\int d^2\alpha Q\left(\alpha \right)ln\left(\pi Q\right(\alpha \left)\right).$$

(45)

Substituting the result in Eq. (21) into the definition above, we obtain a closed-form expression for the Wehrl entropy

$$S_W\left(\widehat{\rho }\right)=1+ln\left[1+n_{\mathrm{th}}(1-e^{-2\gamma t})\right].$$

(46)

Several physically meaningful conclusions can be drawn from this expression. First, at $t=0$ the Wehrl entropy reduces to $S_W=1$, which is the known minimum value attained by a coherent state [29], consistent with the initial condition of the system. Second, in the zero-temperature limit ($n_{\mathrm{th}}=0$), the entropy remains equal to unity for all times, reflecting the fact that the reservoir does not introduce thermal fluctuations and the state retains its coherent character. Third, for $n_{\mathrm{th}}>0$ the entropy grows monotonically in time, quantifying the progressive increase of the phase-space uncertainty due to thermal decoherence. Finally, as $t\to \infty$ the entropy saturates at the asymptotic value $S_{Wss}=1+ln(1+n_{\mathrm{th}})$, which depends solely on the reservoir temperature and corresponds to the thermal equilibrium state of the field, independently of the driving strength $\xi$. Figure 7 shows the time evolution of the Wehrl entropy for different values of $\kappa$. It can be observed that increasing $\gamma$, that is, reducing $\kappa$, causes the entropy to approach its stationary value more rapidly.

5.2. Linear Entropy

As a complementary measure of mixedness that directly quantifies how driving and dissipation affect the purity of the evolved state, we compute the linear entropy, given by

$$S=1-Tr\left\{{\widehat{\rho }}^2\right\}=1-{\displaystyle \frac{\theta }{2-\theta }}=1-{\displaystyle \frac{1}{1+2n_{\mathrm{th}}(1-e^{-2\gamma t})}},$$

(47)

where $Tr\left\{{\widehat{\rho }}^2\right\}$ is the purity of the state. Unlike the Wehrl entropy, $S_L$ vanishes exactly for a pure state and approaches unity for a maximally mixed state, providing a more direct measure of the degree of mixing. Linear entropy displays a qualitatively similar behavior to $S_W$: it grows monotonically from $S_L\left(0\right)=0$, saturates at the asymptotic value $S_{ss}=1-1/(1+2n_{\mathrm{th}})$, and the approach to equilibrium is governed by the same timescale ${\left(2\gamma \right)}^{-1}$. This is illustrated in Figure 8, where the effect of increasing $\gamma$ is again to accelerate thermalization.

It is important to note that the linear entropy is determined exclusively by the reservoir parameters, namely the damping rate $\gamma$ and the mean thermal excitation number $n_{\mathrm{th}}$, and is independent of the driving strength $\xi$. This follows from the fact that the classical driving induces only a unitary displacement in phase space, whereas purity, or the linear entropy, is invariant under unitary transformations. It is worth noting the particular case $n_{\mathrm{th}}=0$, which corresponds to a reservoir at zero temperature. In this limit, the linear entropy vanishes identically $S=0$, at all times. This result has a clear physical interpretation: in the absence of thermal fluctuations, the only effect of the reservoir is to damp the field amplitude, while the classical driving continuously displaces the state in phase space. Since coherent states remain coherent under both amplitude damping at zero temperature and coherent displacement, the system preserves its purity throughout evolution, regardless of the value of $\gamma$. Therefore, the loss of purity originates entirely from thermal fluctuations in the reservoir.

6. Conclusions

We have presented a comprehensive analytical treatment of a driven quantum field coupled to a finite-temperature reservoir using a phase-space approach based on the Husimi Q-function. Starting from the Lindblad master equation, by applying a rotating-frame transformation followed by a displacement operation, the explicit time dependence introduced by the coherent driving was eliminated, reducing the problem to the well-known master equation for a lossy cavity field at finite temperature, solvable via superoperator techniques.

A central result is the derivation of a closed Gaussian form for the Husimi function of an initially coherent state. The Gaussian structure is preserved throughout the evolution, whereas the stationary state corresponds to a displaced thermal state. This representation shows that coherent driving determines the trajectory of the phase-space distribution through a displacement of its center, while the thermal reservoir controls its broadening and stationary width. In particular, for a zero-temperature reservoir, the Gaussian width remains unchanged, indicating that an initial coherent state preserves its purity during the evolution. Furthermore, the Husimi representation enables an analytical description of the superposition of coherent states, showing that the driving field determines the position of the interference structure in phase space, whereas dissipation and thermal fluctuations progressively erase the signatures of the initial superposition. Consequently, the long-time dynamics becomes independent of the initial coherent amplitudes and converges to a unique stationary displaced thermal state. Furthermore, we derive the corresponding Fokker–Planck equation, which confirms that the Husimi function encodes the complete dynamics of the system while identifying the roles of coherent driving, dissipation, and thermal fluctuations in shaping the phase-space distribution.

All relevant photon statistics were obtained in closed form: the mean photon number and standard deviation converge to stationary values reflecting the combined effect of thermal excitations and driving, while the Mandel parameter remains non-negative at all times, confirming that sub-Poissonian statistics are entirely absent and that super-Poissonian behavior is exclusively induced by the thermal reservoir.

The entropic analysis revealed that both the Wehrl and linear entropies grow monotonically from their initial values and saturate at equilibrium values determined solely by the reservoir temperature, independently of the driving strength. This independence reflects the unitary nature of coherent displacement, which cannot alter the purity of the state. In the zero-temperature limit, purity is preserved throughout the entire evolution.

In summary, the Husimi-function approach provides an elegant and analytically tractable framework for characterizing driven-dissipative quantum fields, yielding closed-form expressions for all relevant quantities without requiring explicit manipulation of the density operator. The results obtained constitute a self-consistent and complete characterization of the field’s phase-space, statistical, and entropic properties and generalize previous studies of lossy cavity fields by explicitly incorporating the effect of a resonant classical driving term.