text:
# Effects of surface vibrations on quadrupole response of nuclei

## I INTRODUCTION

There are two systematic features in the quadrupole excitation spectrum of heavy and medium-heavy nuclei: the giant quadrupole resonance at an energy $`\mathrm{}\omega 63A^{\frac{1}{3}}\mathrm{MeV}`$ and the lower-energy states that have often been interpreted as surface vibrations. However the exact nature of these states is still under debate since they are strongly affected by shell effects and this may be taken as an indication that they are not pure surface modes (see e. g. , p.14).

There is a vast literature on the quadrupole response of nuclei, here we refer only to the nice pedagogical introduction in  and to a recent paper  where the effect of coupling between the motion of individual nucleons and surface oscillations has been studied in a model that allows also for non-linear effects and for collisions between nucleons. The scope of our present work is more restricted than that of , we aim at studying only the effects of coupling between the motion of individual nucleons and surface oscillations of the quadrupole type. Our approach is semiclassical and is based on the solution of the linearized Vlasov kinetic equation. Solutions of this equation for finite systems have been obtained by using different boundary conditions (fixed- and moving-surface , see also ), while the fixed-surface solution can give a reasonable picture of the giant quadrupole resonance, it does not satisfactorily describe also the low-lying states, however a unified description can be achieved within this model if the moving-surface boundary conditions introduced in  are employed. Physically, this means that we are including in our model a coupling between the motion of individual nucleons and the surface vibrations. In our approach this coupling does not involve surface vibrations of different multipolarity. Thus, the present model does not include effects that are not taken into account also in the random-phase approximation.

## II Formalism

### A  Fixed-surface solution

Within the fixed-surface theory, and assuming a simplified residual interaction of separable form,

$$v(r_1,r_2)=\kappa _2r_1^2r_2^2,$$
(1)
the quadrupole response function of a spherical nucleus described as a system of $`A`$ interacting nucleons contained within a cavity of radius $`R=1.2A^{\frac{1}{3}}\mathrm{fm}`$ is given by

$$_{22}(s)=\frac{_{22}^0(s)}{1\kappa _2_{22}^0(s)},$$
(2)
where $`s=\frac{\omega R}{v_F}`$ is a convenient dimensionless variable ($`v_F`$ is the Fermi velocity). The zero-order response function $`_{22}^0(s)`$ is analogous to the single-particle response function of the quantum theory and is given explicitly by

$$_{22}^0(s)=\frac{9A}{8\pi }\frac{1}{ϵ_F}\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}\underset{N=0,\pm 2}{}C_{2N}^2_0^1𝑑xx^2s_{nN}(x)\frac{(Q_{nN}^2(x))^2}{s+i\epsilon s_{nN}(x)}.$$
(3)
Here $`ϵ_F`$ is the Fermi energy, the coefficients $`C_{2N}^2`$ are $`C_{20}^2=\frac{1}{4}`$ and $`C_{2\pm 2}^2=\frac{3}{8}`$. The functions $`s_{nN}(x)`$ are defined as

$$s_{nN}(x)=\frac{n\pi +N\mathrm{arcsin}(x)}{x}$$
(4)
and the quantity $`\epsilon `$ is a vanishingly small parameter that determines the integration path at poles. The Fourier coefficients $`Q_{nN}^2(x)`$ are the classical limit of the quantal radial matrix elements and are given by:

$$Q_{nN}^2(x)=()^nR^2\frac{2}{s_{nN}^2(x)}\left(1+N\frac{\sqrt{1x^2}}{s_{nN}(x)}\right)\mathrm{for}(n,N)(0,0)$$
(5)
and

$$Q_{00}^2(x)=R^2(1\frac{2}{3}x^2).$$
(6)
The response function (3) involves an infinite sum over $`n`$, however in practice it is sufficient to include only a few terms around $`n=0`$ in order to fulfill the energy-weighted sum rule with good accuracy. The form (6) of the coefficient $`Q_{00}^2`$ implies that the term $`(n,N)=(0,0)`$ does not contribute to the strength function, hence this term can be omitted from the sum in Eq. (3).

### B  Moving-surface solution

Within the moving-surface theory of , the collective response function (2) is replaced by

$$\stackrel{~}{}_{22}(s)=_{22}(s)+𝒮_{22}(s),$$
(7)
with $`_{22}(s)`$ still given by Eq. (2), while $`𝒮_{22}(s)`$ represents the moving-surface contribution. With the simple interaction (1) the function $`𝒮_{22}(s)`$ can be evaluated explicitly as (see Appendix)

$$𝒮_{22}(s)=\frac{R^6}{1\kappa _2_{22}^0(s)}\frac{[\chi _2^0(s)+\kappa _2\varrho _0R^2_{22}^0(s)]^2}{[C_2\chi _2(s)][1\kappa _2_{22}^0(s)]+\kappa _2R^6[\chi _2^0(s)+\varrho _0R^2]^2},$$
(8)
with $`C_2=4\sigma R^2`$ ($`\sigma 1\mathrm{M}\mathrm{e}\mathrm{V}\mathrm{fm}^2`$ is the surface tension parameter obtained from the mass formula) and $`\varrho _0=A/\frac{4\pi }{3}R^3`$ the equilibrium density.

The functions $`\chi _2^0(s)`$ and $`\chi _2(s)`$ are defined as in Refs.  and  and are given by

$$\chi _2^0(s)=\frac{9A}{4\pi }\frac{1}{R^3}\underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}\underset{N=0,\pm 2}{}C_{2N}^2_0^1𝑑xx^2s_{nN}(x)\frac{()^nQ_{nN}^2(x)}{s+i\epsilon s_{nN}(x)},$$
(9)
and

$`\chi _2(s)`$ $`=`$ $`{\displaystyle \frac{9A}{2\pi }}ϵ_F(s+i\epsilon )\{{\displaystyle \frac{1}{4}}{\displaystyle _0^1}dxx^3\mathrm{cot}[(s+i\epsilon )x]+{\displaystyle \frac{3}{8}}{\displaystyle _0^1}dxx^3(\mathrm{cot}[(s+i\epsilon )x2\mathrm{arcsin}(x)]`$ (10)
$`+`$ $`\mathrm{cot}[(s+i\epsilon )x+2\mathrm{arcsin}(x)])\}.`$ (11)
Equation (8) is the main result of the present paper, its explicit derivation is lengthy but straightforward, the main steps are outlined in the Appendix. Together with Eq. (7), Eq. (8) gives a unified expression for the quadrupole response function, including both the high-energy giant resonance and the low-energy excitations. By comparing the two response functions (2) and (7) we can appreciate the effect of the additional surface degree of freedom introduced in  and in particular the effect of coupling the motion of nucleons with surface vibrations of quadrupole type.

## III Results

In Fig.1 we display the strength function ($`E=\mathrm{}\omega `$)

$$S(E)=\frac{1}{\pi }\mathrm{Im}(E),$$
(12)
obtained for $`A=208`$ using different approximations. The dotted curve is obtained from the zero-order response function (3) and it is similar to the quantum response evaluated in the Hartree-Fock approximation. The dashed curve is obtained from the collective fixed-surface response function (2). Comparison with the dotted curve clearly shows the effects of collectivity: the main strength at about $`16\mathrm{MeV}`$ is shifted to lower excitation energy and the giant quadrupole peak becomes narrower. The strength of the interaction (1), chosen in order to reproduce the experimental value of the giant quadrupole resonance energy in $`{}_{}{}^{208}\mathrm{Pb}`$, is $`\kappa _2=\mathrm{1.\hspace{0.17em}10}^3\mathrm{MeV}\mathrm{fm}^4`$. this value is close to that suggested by the Bohr-Mottelson prescription (, p. 509)

$$\kappa _2=\frac{4\pi }{3}\frac{m\omega _0^2}{AR^2}\mathrm{0.5\hspace{0.17em}10}^3\mathrm{MeV}\mathrm{fm}^4,$$
(13)
(with $`\omega _0=41A^{\frac{1}{3}}\mathrm{MeV}`$).

We notice that the width of the giant quadrupole resonance is underestimated by the fixed-surface model, this is a well known limit of all mean-field calculations that include only Landau damping. Moreover, there is no sign of a low-energy peak in the fixed-surface response function. The solid curve instead shows the moving-surface response given by Eqs. (7)and (8). Now a broad bump appears in the low-energy part of the response and a narrower peak is situated at the giant resonance energy, thus the moving-surface solution of the Vlasov equation introduced in Ref.  accounts for both quadrupole modes, although only qualitatively. Of course the details of the low-energy excitations are determined by quantum effects, nonetheless the present semiclassical approach does reproduce the average behaviour of this systematic feature of the quadrupole response.

Another remarkable feature of the moving-surface response function is that now the giant quadrupole peak is narrower than for the fixed-surface solution. This is somewhat surprising since we could have expected that introducing a further degree of freedom would result in a smearing of the peak, however our result for the giant resonance region is very similar to that of the recent random-pase approximation (RPA) calculations of (cf. Fig.5 of ). Our model does include Landau dampig that, however, turns out to be very small in this case, clearly some additional mechanism is required in order to increase the width of the giant resonance. Two such mechanisms have been considered in Refs. , they are the coupling to surface vibrations of different multipolarity and the effect of collisions between nucleons. It would be interesting to include such effects in the present semiclassical theory, however this will be left for future work.

All the strength functions shown in Fig.1 should satisfy the following energy-weighted sum rule (EWSR) (see e.g. , p. 401):

$$_0^{\mathrm{}}𝑑EES(E)=\frac{3}{4\pi }\frac{\mathrm{}^2}{m}AR^2.$$
(14)
We have numerically checked that, when integrated up to $`E=30\mathrm{M}\mathrm{e}\mathrm{V}`$, the response function shown by the solid curve exhausts about $`98\%`$ of this sum rule. The fraction of EWSR exhausted by the dashed and dotted strength functions in the same interval is only $`80\%`$, showing that in these cases there is some more strength at higher energy.

Another interesting moment of the strength function is the inverse energy-weighted sum rule:

$$m_1=_0^{\mathrm{}}𝑑E\frac{S(E)}{E}.$$
(15)
From the three strength functions shown in Fig.1, we have three different inverse moments :

$$m_1^0=\frac{1}{2}\underset{s0}{lim}_{22}^0(s)=\frac{139}{140}\frac{1}{16\pi }\frac{AR^4}{ϵ_F}$$
(16)
(when evaluating $`lim_{s0}_{22}^0(s)`$ the term $`(n,N)=(0,0)`$ must be omitted from the sum in Eq. (3)),

$$m_1=m_1^0\frac{1}{1+2\kappa _2m_1^0}$$
(17)
and

$$\stackrel{~}{m}_1=m_1\left\{1+12\frac{140}{139}\frac{[\frac{17}{20}+2\kappa _2m_1^0]^2}{[1+2\kappa _2m_1^0][1+\frac{8}{3}\frac{4\pi \sigma R^2}{Aϵ_F}+\frac{27}{50}\frac{140}{139}\kappa _2m_1]}\right\}.$$
(18)
The inverse moment of the zero-order strength function (dotted curve in Fig.1) $`m_1^0`$ exhausts about $`97\%`$ of the sum rule (16) in the range from $`0`$ to $`30\mathrm{MeV}`$. The collective fixed-surface response (dashed curve in Fig.1) instead exhausts almost $`99\%`$ of the sum rule (17) in the same energy range while the moving-surface strength function exhusts almost $`100\%`$ of the sum rule (18), always in the same energy interval.

It can be seen from Fig.1 that, while the fixed-surface response has only one collective pole, the moving-surface quadrupole response function displays a two-pole structure. In order to get more information about the nature of the low-energy peak, we have performed calculations by putting the surface tension parameter $`\sigma `$ equal to zero. The result is shown in Fig.2, where the dotted curve corresponds to $`\sigma =0`$. As expected, the giant resonace peak is practically unaffected by the surface tension, while the low-energy peak is affected quite substantially. The surface tension increases the frequency of the low-enegy peak, which, however, is present at a non-vanishing frequency also in the absence of surface tension. In the opposite limit, if we let $`\sigma \mathrm{}`$, the fixed-surface response is obtained.

We have performed calculations of the quadrupole response functions also for other values of A corresponding to medium-heavy spherical nuclei and the results are qualitatively similar to the $`A=208`$ case, so we do not report them here.

A calculation of the isoscalar quadrupole response similar to the present one has been made in Ref. , in that case however a rather special external field has been assumed. The external force studied in  is a pressure that acts only on the surface of the nucleus, perhaps this explains why in that case very little strength was found in the region of the giant quadrupole resonance.

## IV Summary and Conclusions

We have obtained an analytical expression for the isoscalar quadrupole response function of nuclei that describes qualitatively the two main systematic features of the excitation spectrum: the giant resonance and the low-energy states. Our approach is semiclassical and is based on the solution of the linearized Vlasov kinetic equation with appropriate boundary conditions (moving surface). Comparison of our result (full curve) with the quantum response function of Ref.  shows that quantum effects modify substantially the low-energy region (a discrete state is obtained in  instead of our broad bump), while the giant resonance peak is practically the same in the quantum and semiclassical approach.

## Moving-surface response

In our semiclassical approach the nucleus is derscribed by means of the phase-space density $`f(𝐫,𝐩,t)`$. At equilibrium the density is $`f_0(𝐫,𝐩)=F(h_0(𝐫,𝐩))`$, where $`h_0`$ is the equilibrium mean-field hamiltonian. A weak external driving field $`V_{ext}(𝐫,t)=\beta (t)Q(𝐫)`$ induces small fluctuations of the equilibrium distribution $`f_0`$ and of the mean field. Since the two fluctuations are related, we have a typical self-consistency problem. For spherical system the problem can be solved by means of the partial-wave expansion

$$\delta f(𝐫,𝐩,\omega )=\underset{LMN}{}\left[\delta \stackrel{~}{f}_{MN}^{L+}(ϵ,\lambda ,r,\omega )+\delta \stackrel{~}{f}_{MN}^L(ϵ,\lambda ,r,\omega )\right]\left(𝒟_{MN}^L(\alpha ,\beta ,\gamma )\right)^{}Y_{LN}(\frac{\pi }{2},\frac{\pi }{2}),$$
(19)
where $`𝒟_{MN}^L(\alpha ,\beta ,\gamma )`$ are the Wigner rotation matrices, $`ϵ=h_0`$ is the particle energy and $`\lambda `$ its angular momentum.

In the approach of  the functions $`\delta \stackrel{~}{f}_{MN}^{L\pm }`$ satisfy the integral equation (see Ref.  for details)

$`\delta \stackrel{~}{f}_{MN}^{L\pm }(ϵ,\lambda ,r,\omega )=F^{}(ϵ){\displaystyle \frac{e^{\pm i[\omega \tau (r)N\gamma (r)]}}{\mathrm{sin}[\omega \tau (R)N\gamma (R)]}}\omega p_r(R)\delta R_{LM}(\omega )`$ (20)
$`+\left\{{\displaystyle _{r_1}^r}𝑑r^{}\stackrel{~}{B}_{MN}^{L\pm }(ϵ,\lambda r^{})e^{[i\omega \tau (r^{})N\gamma (r^{})]}+\stackrel{~}{C}_{MN}^L(ϵ,\lambda ,\omega )\right\}e^{\pm i[\omega \tau (r)N\gamma (r)]},`$ (21)
with

$$\stackrel{~}{B}_{MN}^{L\pm }(ϵ,\lambda ,r,\omega )=F^{}(ϵ)\left(\frac{}{r}\pm \frac{iN}{v_r(ϵ,\lambda ,r)}\frac{\lambda }{mr^2}\right)\left[\beta (\omega )Q_{LM}(r)+\delta \stackrel{~}{V}_{LM}(r,\omega )\right],$$
(22)
$`\delta \stackrel{~}{V}_{LM}(r,\omega )={\displaystyle \frac{8\pi ^2}{2L+1}}{\displaystyle \underset{N=L}{\overset{L}{}}}\left|Y_{LN}({\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi }{2}})\right|^2`$ (23)
$`\times {\displaystyle }dϵ{\displaystyle }d\lambda \lambda {\displaystyle }{\displaystyle \frac{dr^{}}{v_r(r^{})}}v_L(r,r^{})[\delta \stackrel{~}{f}_{MN}^{L+}(ϵ,\lambda ,r^{},\omega )+\delta \stackrel{~}{f}_{MN}^L(ϵ,\lambda ,r^{},\omega )]`$ (24)
and

$`\stackrel{~}{C}_{MN}^L(ϵ,\lambda ,\omega )=\{e^{i2[\omega \tau (R)N\gamma (R)]}{\displaystyle _{r_1}^R}dr\stackrel{~}{B}_{MN}^{L+}(ϵ,\lambda ,r)e^{i[\omega \tau (r)N\gamma (r)]}`$ (25)
$`{\displaystyle _{r_1}^R}dr\stackrel{~}{B}_{MN}^L(ϵ,\lambda ,r)e^{i[\omega \tau (r)N\gamma (r)]}\left\}\right\{1e^{i2[\omega \tau (R)N\gamma (R)]}\}^1,`$ (26)
which is equivalent to the linearized Vlasov equation with the moving-surface boundary conditions

$$\delta \stackrel{~}{f}_{MN}^{L+}(R)\delta \stackrel{~}{f}_{MN}^L(R)=2F^{}(ϵ)i\omega p_r(R)\delta R_{LM}(\omega ).$$
(27)
The surface fluctuations $`\delta R_{LM}(\omega )`$ are related to the functions $`\delta \stackrel{~}{f}_{MN}^{L\pm }`$ by

$`\delta R_{LM}(\omega )={\displaystyle \frac{8\pi ^2}{2L+1}}{\displaystyle \frac{R^2}{C_L}}{\displaystyle \underset{N=L}{\overset{L}{}}}\left|Y_{LN}({\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi }{2}})\right|^2`$ (28)
$`{\displaystyle 𝑑ϵ𝑑\lambda \lambda p_r(R)\left[\delta \stackrel{~}{f}_{MN}^{L+}(ϵ,\lambda ,R,\omega )+\delta \stackrel{~}{f}_{MN}^L(ϵ,\lambda ,R,\omega )2F^{}(ϵ)\delta \stackrel{~}{V}_{LM}(R,\omega )\right]}.`$ (29)
For a separable interaction of the multipole-multipole type, like (1), the integral equation (20) can be reduced to an algebraic equation that, in the particular case $`L=2`$ and $`V_{ext}=\beta (t)r^2Y_{2M}(\widehat{𝐫})`$, gives Eqs.(7) and (8) for the quadrupole response function.


format: text

text:
# The production mechanism of radio jets in AGN and quasar grand unification

### The production of powerful radio jets

There is good observational evidence that radio jet power is closely linked to both the mass of the black hole and its accretion rate. Links between the mass of the black hole and its radio luminosity have been found in both radio-weak elliptical galaxies (Franceschini, Vercellone & Fabian 1998) and quasars (Laor 2000; Lacy et al. 2001; Boroson 2002). Essentially all steep-spectrum radio-selected quasars and BLLac objects have black holes more massive than $`10^8M_{}`$ (e.g. McLure & Jarvis 2002; Barth, Ho & Sargent 2002; O’Dowd, Urry & Scarpa 2002), and all powerful radio galaxies exist in massive hosts, almost invariably giant ellipticals. Recent claims that a significant fraction of flat-spectrum radio-luminous quasars have black holes $`10^{67}M_{}`$ (Oshlack, Webster & Whiting 2002), have been questioned by Jarvis & McLure (2002) on the basis of orientation biases.

Radio luminosity is roughly linearly correlated with accretion rate for quasars (e.g. Lacy et al. 2001). At low accretion rates (in the ADAF regime), however, black holes seem to be able to deliver more radio jet power than expected if the linear correlation continues to low accretion rates, i.e. low luminosity AGN tend to be radio louder than most quasars (Ho 2002).

Most models predict than black hole spin should be important for determining radio jet luminosity. Indeed, the residual scatter of radio luminosity about the best-fit combination of black hole mass and accretion rate (a radio “fundamental plane”) is 1-2 orders of magnitude, and probably not just due to measurement error, suggesting that spin, and/or some other parameter(s) need to be taken into account. Recent results from a study of the Fe K$`\alpha `$ line in MCG-6-30-15 indicates that the black hole in this radio-quiet Seyfert galaxy may be spinning rapidly, however (Fabian, these proceedings).

### Broad absorption line quasars

The fraction of highly radio-luminous BALQs is small compared to that in samples of radio-quiet quasars. If broad absorption lines are only seen in objects accreting close to the Eddington rate (Boroson 2002), then this might be explained as a selection effect. Most radio-loud quasars have high black hole masses, thus a radio-loud quasar accreting at near the Eddington rate would be highly luminous, and hence very rare.

### Acknowledgments.

I thank Aaron Barth and Mike Brotherton for helpful discussions. This work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA.

## References

Barth A.J., Ho L.C., Sargent W.L.W., 2002, ApJ, in press (astro-ph/0209562)

Boroson T.A., 2002, ApJ, 565, 78

Franceschini, A., Vercellone, S., & Fabian, A.C. 1998, MNRAS, 297, 817

Ho, L. 2002, ApJ, 564, 120

Jarvis M.J. & McLure R.J., 2002, MNRAS, submitted (astro-ph/0208390)

Lacy, M., Laurent-Muehleisen, S.A., Ridgway, S.E., Becker, R.H. & White R.L. 2001, ApJ, 551, L17

Laor, A. 2000, ApJ, 543, L111

McLure, R.J. & Jarvis M.J. 2002, MNRAS, in press (astro-ph/0204473)

O’Dowd M., Urry C.M., Scarpa R., 2002, ApJ, in press (astro-ph/0207597)

Oshlack, A.Y.K.N., Webster, R.L. & Whiting, M.T. 2002, ApJ, 576, 810


text:
# Reduced density matrices of oscillator systems

## 1 Statement of the problem

This part is devoted to the multidimensional generalizations of the harmonic oscillator model. As a rule, these are various systems of interacting oscillators, some of which may be placed in external uniform electric and magnetic fields. More precisely this class of systems can be characterized by a Hamiltonian that is assumed to be in the form of the general inhomogeneous quadratic form

$$\widehat{H}=\frac{1}{2}\widehat{𝐪}𝐁\widehat{𝐪}+𝐂\widehat{𝐪}.$$
(1.1)
Here, $`𝐪`$ is conceived as a $`2N`$-dimensional vector ($`N`$ being the number of degrees of freedom) whose components are linear combinations of Cartesian coordinates and momenta conjugated to them. The most frequent choice is $`𝐪=(p_1,\mathrm{},p_N,x_1,\mathrm{},x_N)`$, but in the presence of an external magnetic field it is more convenient to deal not with the canonical momenta $`p_j`$ but with the kinetic momenta $`\pi _j=p_jeA_j(𝐱)/c`$, where $`𝐀(𝐱)`$ is the vector potential. Other choices are also possible, taking into account the concrete physical applications. For instance, the components of the vector $`𝐪`$ may be bosonic annihilation and creation operators constructed from the coordinate and momentum operators, etc. $`𝐁`$ is a symmetric $`2N\times 2N`$ matrix, which may depend on time, as also the $`2N`$-dimensional vector $`𝐂`$.

Multidimensional systems with Hamiltonian (1.1) were the subject of investigation in numerous papers: see, e.g., , - and the references therein. In the present paper we consider the case where the multidimensional vector $`𝐪`$ may be split in two parts: $`𝐪=(𝐐,𝝃)`$, where the vector $`𝐐`$ describes a small subsystem, whereas the vector $`𝝃`$ is related to a “thermostat”. In the simplest cases the evolution of a unidimensional harmonic oscillator coupled with a “thermostat” was studied, e.g., in ,-. We investigate the most general quadratic “interaction Hamiltonians” and compare the results of this “microscopic” approach with different phenomenological models considered in reviews . Emphasis is placed on application to the problem of harmonic oscillator relaxation, including the case where a (charged) oscillator is placed in a uniform magnetic field.

## 2 Phenomenological Fokker-Planck equation

The evolution of any closed system is governed by the quantum Liouville equation for the statistical operator $`\widehat{\rho }`$,

$$i\mathrm{}\widehat{\rho }/t=\widehat{H}\widehat{\rho }\widehat{\rho }\widehat{H}.$$
(2.1)
For any quadratic Hamiltonian (1.1) this equation results in the linear equation for the average values of the $`2N`$-vector $`\widehat{𝐪}`$:

$$\dot{𝐪}=\mathrm{\Sigma }𝐁𝐪\mathrm{\Sigma }𝐂,$$
(2.2)
where $`𝐪=\text{Tr}(\widehat{\rho }\widehat{𝐪})`$, and antisymmetric nondegenerate $`2N\times 2N`$ matrix $`\mathrm{\Sigma }`$ with $`c`$-number coefficients is defined via the commutation relations between the operators $`\widehat{q}_\alpha `$,

$$[\widehat{q}_\alpha ,\widehat{q}_\beta ]=i\mathrm{}\mathrm{\Sigma }_{\alpha \beta },\mathrm{\Sigma }=\mathrm{\Sigma }_{\alpha \beta },\alpha ,\beta =1,2,\mathrm{}2N.$$
(2.3)
The operator equation (2.1) can be transformed into a partial differential equation upon the choice of some concrete representation for the statistical operator. For instance, in the coordinate representation we write $`𝐪=(𝐩,𝐱)`$, p and x being the $`N`$-vectors, and split the matrix B into $`N\times N`$ blocks:

$$𝐁=\begin{array}{cc}𝐛_1& 𝐛_2\\ 𝐛_3& 𝐛_4\end{array},𝐛_1=\stackrel{~}{𝐛}_1,𝐛_2=\stackrel{~}{𝐛}_3,𝐛_4=\stackrel{~}{𝐛}_4,𝐛_j=b_j^{mn}$$

(the tilde designates a transposed matrix). Then we obtain the following second-order equation for the density matrix $`\rho (x,x^{})`$:

$`i\mathrm{}{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2}}b_1^{mn}\left({\displaystyle \frac{^2\rho }{x_m^{}x_n^{}}}{\displaystyle \frac{^2\rho }{x_mx_n}}\right)i\mathrm{}b_2^{mn}\left(x_n{\displaystyle \frac{\rho }{x_m}}+x_n^{}{\displaystyle \frac{\rho }{x_m^{}}}\right)`$ (2.4)
$`+`$ $`{\displaystyle \frac{1}{2}}b_4^{mn}\left(x_mx_nx_m^{}x_n^{}\right)\rho i\mathrm{}\rho \text{Tr}𝐛_2`$
$``$ $`i\mathrm{}c_1^m\left({\displaystyle \frac{\rho }{x_m}}+{\displaystyle \frac{\rho }{x_m^{}}}\right)+c_2^m\left(x_mx_m^{}\right).`$
Here $`𝐜_1`$ and $`𝐜_2`$ are $`N`$-dimensional components of the vector $`𝐂=(𝐜_1,𝐜_2)`$. The disadvantage of Eq. (2.4) is the broken symmetry between the coordinates and momenta, which is inherent in the Hamiltonian. This symmetry can be restored if one proceeds from the complex density matrix to the real Wigner function

$$W(𝐩,𝐱)=\rho (𝐱+\xi /2,𝐱\xi /2)\mathrm{exp}(i𝐩\xi /\mathrm{})d\xi ,$$
(2.5)
$$\rho (𝐱,𝐱^{})=W(𝐩,\frac{1}{2}(𝐱+𝐱^{}))\mathrm{exp}\left[i𝐩(𝐱𝐱^{})/\mathrm{}\right]\text{d}𝐩/(2\pi \mathrm{})^N.$$
(2.6)
Applying transformation (2.5) to Eq. (2.4) we obtain an equivalent equation that is first order with respect to all the derivatives:

$$\frac{W}{t}=\left(p_nb_2^{nm}+x_nb_4^{nm}+c_2^m\right)\frac{W}{p_m}\left(p_nb_1^{nm}+x_nb_3^{nm}+c_1^m\right)\frac{W}{x_m}.$$
(2.7)
This equation demonstrates the distinction and the advantage of the Wigner function for the description of quadratic quantum systems (other remarkable features of the Wigner function were discussed, e.g., in ). Eq. (2.7) assumes an especially compact form in terms of the $`2N`$-vector q,

$$\frac{W}{t}=\frac{}{q_\alpha }\left[(\mathrm{\Sigma }\mathrm{𝐁𝐪}+\mathrm{\Sigma }𝐂)_\alpha W\right].$$
(2.8)
Here the $`2N\times 2N`$ matrix $`\mathrm{\Sigma }`$, in accordance with Eq. (2.3), equals

$$\mathrm{\Sigma }=\begin{array}{cc}\mathrm{𝟎}& 𝐈_N\\ 𝐈_N& \mathrm{𝟎}\end{array},$$
(2.9)
$`𝐈_N`$ being the $`N\times N`$ unit matrix.

Now let us suppose that we have a given linear equation for the first-order average values

$$\dot{𝐪}=𝐀(t)𝐪+𝐊(t),$$
(2.10)
with an arbitrary matrix $`𝐀(t)`$ and an arbitrary vector $`𝐊(t)`$. The problem investigated in this section is whether it is possible to find an equation for the statistical operator or the Wigner function that would result in the given Eq. (2.10). It is trivial to check that Eq. (2.10) is the consequence of the equation

$$\frac{W}{t}=\frac{}{q_\alpha }\left[(\mathrm{𝐀𝐪}+𝐊)_\alpha W\right].$$
(2.11)
However, although Eq. (2.11) preserves the normalization of the Wigner function,

$$W(𝐪)\text{d}𝐪/(2\pi \mathrm{})^N=1,$$
(2.12)
it does not agree, in general, with the fundamental quantum mechanical principle of the positive definiteness of the statistical operator and (this is almost the same) with the uncertainty relations. Indeed, let us consider the operator $`\widehat{F}=\alpha _j\left(\widehat{q}_jq_j\right)`$ with arbitrary complex coefficients $`\alpha _j`$. For any stattistical operator, due to its nonnegative definiteness, the inequality $`\widehat{F}^{}\widehat{F}\text{Tr}\left(\widehat{\rho }\widehat{F}^{}\widehat{F}\right)0`$ must hold. Taking into account the structure of the operator $`\widehat{F}`$, we arrive at the conclusion on the nonnegative definiteness of the bilinear Hermitian form

$$\alpha ^{}\mathrm{\Phi }\alpha \alpha _j^{}\mathrm{\Phi }_{jk}\alpha _k,\mathrm{\Phi }_{jk}=\left(\widehat{q}_jq_j\right)\left(\widehat{q}_kq_k\right),$$

whose matrix $`\mathrm{\Phi }`$ is constructed from the centered second-order moments of the operators $`\widehat{q}_j`$, $`j,k=1,2,\mathrm{},2N`$. It is convenient to distinguish even and odd parts of the matrix $`\mathrm{\Phi }`$. The symmetric part consists of the symmetrized second moments (covariances)

$$_{ij}=_{ji}=\frac{1}{2}\widehat{q}_j\widehat{q}_k+\widehat{q}_k\widehat{q}_jq_jq_k.$$
(2.13)
The antisymmetric part, in accordance with Eq. (2.3), is expressed through the commutator matrix, so that

$$\mathrm{\Phi }=\frac{i\mathrm{}}{2}\mathrm{\Sigma }.$$

The elements of the matrix $``$ are calculated in terms of the Wigner function as follows:

$$_{ij}=q_iq_jW(𝐪)\text{d}𝐪/(2\pi \mathrm{})^N\overline{q}_i\overline{q}_j,$$

$$\overline{𝐪}=𝐪=𝐪W(𝐪)\text{d}𝐪/(2\pi \mathrm{})^N.$$

Due to Eq. (2.11), the variance matrix satisfies the equation

$$\dot{}=𝐀+\stackrel{~}{𝐀}.$$
(2.14)
The first-order averages do not influence the variances in the case under study.

If the coefficients $`\alpha _j`$ are chosen in such a way that the commutator $`[\widehat{F},\widehat{F}^{}]`$ is positive (the simplest example is $`\widehat{F}=\widehat{x}+i\widehat{p}`$), then the operator $`\widehat{F}`$ coincides within a constant factor with a boson annihilation operator. Choosing the initial state to be the vacuum state for this operator, at the initial instant we have

$$\text{Tr}\left(\widehat{F}^{}\widehat{F}\widehat{\rho }(0)\right)=\alpha ^{}\mathrm{\Phi }(0)\alpha =0.$$

Moreover, the matrix equalities $`\mathrm{\Phi }(0)\alpha =\alpha ^{}\mathrm{\Phi }(0)=0`$ hold as well. For the chosen initial state we have at $`t>0`$, due to Eq. (2.14),

$`\alpha ^{}\mathrm{\Phi }(t)\alpha `$ $`=`$ $`\alpha ^{}\left[(t){\displaystyle \frac{i\mathrm{}}{2}}\mathrm{\Sigma }\right]\alpha `$ (2.15)
$`=`$ $`\alpha ^{}\left[𝐀(0)(0)+(0)\stackrel{~}{𝐀}(0)\right]t\alpha +𝒪(t^2)`$
$`=`$ $`{\displaystyle \frac{1}{2}}i\mathrm{}t\alpha ^{}\left[𝐀(0)\mathrm{\Sigma }+\mathrm{\Sigma }\stackrel{~}{𝐀}(0)\right]\alpha +𝒪(t^2).`$
For the Hamiltonian systems we have $`𝐀=\mathrm{\Sigma }𝐁,\stackrel{~}{𝐀}=𝐁\mathrm{\Sigma }`$, and the linear with respect to time term disappears. But for an arbitrary matrix $`𝐀`$ the right-hand side of Eq. (2.15) can be negative. For example, in the model of a damped oscillator with the equations for the averages

$$\dot{x}=p,\dot{p}=\omega _0^2x2\gamma p,$$
(2.16)
the matrices entering Eq. (2.15) read

$$𝐀=\begin{array}{cc}2\gamma & \omega _0^2\\ 1& 0\end{array},𝐀\mathrm{\Sigma }+\mathrm{\Sigma }\stackrel{~}{𝐀}=\begin{array}{cc}0& 2\gamma \\ 2\gamma & 0\end{array}.$$
(2.17)
Then, for the initial vacuum state of operator $`\widehat{F}=\widehat{x}+i\widehat{p}`$, when $`\alpha =(i,1)`$, we obtain from Eq. (2.15)

$$\alpha ^{}\mathrm{\Phi }(t)\alpha =2\mathrm{}t\gamma +𝒪(t^2)<0.$$

Consequently, in the general case Eq. (2.11) is unacceptable. Thus we need more complicated generalizations. The simplest possibility is to add terms with second derivatives to the right-hand side of Eq. (2.11), i.e., to transform this equation into the Fokker-Planck equation:

$$\frac{W}{t}=\frac{}{q_\alpha }\left[(\mathrm{𝐀𝐪}+𝐊)_\alpha W\right]+D_{\alpha \beta }\frac{^2W}{q_\alpha q_\beta },$$
(2.18)
where the diffusion coefficients $`D_{\alpha \beta }=D_{\beta \alpha }`$, combined into a symmetric matrix $`𝐃=D_{\alpha \beta }`$, may depend on time but do not depend on the coordinates. The new “diffusion” terms do not change the equation for the average values (2.10), moreover, they preserve the normalization (2.12). But they enable one to “save” the nonnegative definiteness of the statistical operator. Indeed, considering the evolution of the bilinear form $`\alpha ^{}\mathrm{\Phi }(t)\alpha `$, we obtain, instead of Eq. (2.15), the equation

$$\alpha ^{}\mathrm{\Phi }(t)\alpha =t\alpha ^{}\left(2𝐃+\frac{i\mathrm{}}{2}\left[𝐀(0)\mathrm{\Sigma }+\mathrm{\Sigma }\stackrel{~}{𝐀}(0)\right]\right)\alpha +𝒪(t^2),$$

since Eq. (2.14) is replaced by

$$\dot{}=𝐀+\stackrel{~}{𝐀}+2𝐃.$$
(2.19)
Consequently, the necessary condition of the compatibility of Eq. (2.18) with the principles of quantum mechanics is the nonnegative definiteness of the matrix

$$𝐃_{}𝐃+\frac{i\mathrm{}}{4}\left[𝐀\mathrm{\Sigma }+\mathrm{\Sigma }\stackrel{~}{𝐀}\right]0$$
(2.20)
at any instant of time. Moreover, it can be proved  that the condition $`𝐃_{}0`$ is sufficient as well.

Vector q was defined above as $`𝐪=(𝐩,𝐱)`$, and the matrix $`\mathrm{\Sigma }`$ had the explicit form given by Eq. (2.9). Let us make the time-independent transformation of the variables

$$𝐪^{}=\mathrm{𝐓𝐪},det𝐓0,\text{Im}𝐓0.$$

Then Eqs. (2.10) and (2.18) preserve their forms, provided the matrices A and D are replaced by the matrices $`𝐀^{}=\mathrm{𝐓𝐀𝐓}^1`$ and $`𝐃^{}=\mathrm{𝐓𝐃}\stackrel{~}{𝐓}`$. The nonnegative definiteness of the matrix $`𝐃_{}`$ of Eq. (2.20) is equivalent to the nonnegative definiteness of the matrix $`𝐃_{}^{}=\mathrm{𝐓𝐃}_{}\stackrel{~}{𝐓}`$. The latter has, in turn, again the form (2.20), if one replaces the matrix $`\mathrm{\Sigma }`$ by $`\mathrm{\Sigma }^{}=𝐓\mathrm{\Sigma }\stackrel{~}{𝐓}`$, but this is just the transformation law of any matrix defined according to Eq. (2.3). This way we arrive at the important conclusion, that all the formulas containing the matrices A, D, and $`\mathrm{\Sigma }`$ are valid not only in the case where the components of the vector q coincide with the canonically conjugate momenta and Cartesian coordinates, but also in the general case where the components of the vector $`𝐪`$ are arbitrary Hermitian operators with $`c`$-number commutators, provided the matrix $`\mathrm{\Sigma }`$ is defined according to Eq. (2.3).

To transform Eq. (2.18) to an operator form that is independent on the concrete representation, one should take into account the following correspondence relations between the operators $`\widehat{𝐪}\widehat{\rho }`$, $`\widehat{\rho }\widehat{𝐪}`$, and their Weyl symbols (they result from Eqs. (2.5), (2.6)):

$$\widehat{𝐪}\widehat{\rho }\left(𝐪\frac{i\mathrm{}}{2}\mathrm{\Sigma }\frac{}{𝐪}\right)W(𝐪),\widehat{\rho }\widehat{𝐪}\left(𝐪+\frac{i\mathrm{}}{2}\mathrm{\Sigma }\frac{}{𝐪}\right)W(𝐪),$$
(2.21)
$$𝐪W(𝐪)\frac{1}{2}\left(\widehat{𝐪}\widehat{\rho }+\widehat{\rho }\widehat{𝐪}\right),\frac{W}{𝐪}\frac{i}{\mathrm{}}\mathrm{\Sigma }^1\left(\widehat{𝐪}\widehat{\rho }\widehat{\rho }\widehat{𝐪}\right).$$
(2.22)
Making transformations (2.22) in Eq. (2.18) we arrive at the equation

$`{\displaystyle \frac{\widehat{\rho }}{t}}`$ $`=`$ $`{\displaystyle \frac{i}{2\mathrm{}}}\left[\widehat{𝐪}\mathrm{\Sigma }^1𝐀\widehat{𝐪}\widehat{\rho }+\widehat{\rho }\widehat{𝐪}\stackrel{~}{𝐀}\mathrm{\Sigma }^1\widehat{𝐪}+\widehat{𝐪}\left(\mathrm{\Sigma }^1𝐀+\stackrel{~}{𝐀}\mathrm{\Sigma }^1\right)\widehat{\rho }\widehat{𝐪}\right]`$ (2.23)
$`+`$ $`{\displaystyle \frac{i}{\mathrm{}}}\left[\widehat{𝐪}\mathrm{\Sigma }^1𝐊\widehat{\rho }\widehat{\rho }\widehat{𝐪}\mathrm{\Sigma }^1𝐊\right]{\displaystyle \frac{1}{\mathrm{}^2}}\left[\widehat{𝐪}𝐒\widehat{𝐪}\widehat{\rho }+\widehat{\rho }\widehat{𝐪}𝐒\widehat{𝐪}2\widehat{𝐪}𝐒\widehat{\rho }\widehat{𝐪}\right],`$
where the symmetrical matrix

$$𝐒=\mathrm{\Sigma }^1𝐃\mathrm{\Sigma }$$
(2.24)
must satisfy a constraint equivalent to Eq. (2.20):

$$𝐒_{}𝐒\frac{i\mathrm{}}{4}\left[\stackrel{~}{𝐀}\mathrm{\Sigma }^1+\mathrm{\Sigma }^1𝐀\right]0.$$
(2.25)
Eqs. (2.21)-(2.25) hold for any vector q with a $`c`$-number commutator matrix $`\mathrm{\Sigma }`$ of (2.3).

Comparison of the elegant equation (2.18) for the Wigner function with the much more cumbersome equation (2.23) demonstrates once more the advantage of the Wigner representation for describing quantum systems with linear equations of motion for the averages.

There is an important difference between the quantum Fokker-Planck equation (2.18) and its classical counterpart. The classical Fokker-Planck equation contains, as a rule, second derivatives only with respect to momenta. However, such a simple set of the diffusion coefficients is unacceptable in the quantum case (although sometimes this incorrect equation was considered: see, e.g., ). This statement can be easily demonstrated on the example of system (2.16). Writing the matrix D in the form

$$𝐃=\begin{array}{cc}D_p& D_{px}\\ D_{px}& D_x\end{array}$$

and taking into account Eq. (2.17) we obtain the matrix

$$𝐃_{}=\begin{array}{cc}D_p& D_{px}i\mathrm{}\gamma /2\\ D_{px}+i\mathrm{}\gamma /2& D_x\end{array}.$$

The condition of its positive definiteness is given by the inequality -

$$det𝐃_{}D_pD_xD_{px}^2\mathrm{}^2\gamma ^2/40,$$
(2.26)
whose violation leads to the violation of the uncertainty relations , . For an arbitrary one-dimensional system (2.10) the condition $`𝐃_{}0`$ is equivalent to the inequality

$$det𝐃\mathrm{}^2(\text{Tr}𝐀)^2/16.$$
(2.27)
## 3 Fokker-Planck equation for a subsystem

In the preceding section we have shown that any given equation (2.10) may be considered as a consequence of some suitable Fokker-Planck equation for the Wigner function. The only problem is to select the set of diffusion coefficients satisfying the condition $`𝐃_{}0`$. It is clear that by taking sufficiently large diffusion coefficients one can always satisfy this condition. The problem of finding “minimal admissible” diffusion coefficients for some simple systems (such as a one-dimensional harmonic oscillator or a two-dimensional isotropic oscillator in a uniform magnetic field) was investigated in ,,.

Here we investigate the following problem. Suppose we have a large closed quantum system with $`N`$ degrees of freedom, described by Hamiltonian (1.1). Let us split the vector $`𝐪`$ in two parts: $`𝐪=(𝐐,\xi )`$, where the $`n`$-dimensional vector $`𝐐`$ describes a subsystem, while the vector $`\xi `$ relates to a reservoir. The question is: what kind of equation describes the evolution of the subsystem if one performs an averaging over the variables of the reservoir?

The Wigner function of the whole system is given by the relation

$$W(𝐪,t)=G(𝐪,𝐪^{},t)W(𝐪^{},0)\text{d}𝐪^{},$$
(3.1)
where the propagator $`G(𝐪,𝐪^{},t)`$ satisfies Eq. (2.8) and the initial condition $`G(𝐪,𝐪^{},0)=\delta (𝐪𝐪^{})`$. Since Eq. (2.8) is first-order with respect to all the variables, its propagator is extremely simple:

$$G(𝐪,𝐪^{},t)=\delta \left(𝐪𝐪_{}(t;𝐪^{})\right),$$
(3.2)
where the vector

$$𝐪_{}(t;𝐪^{})=𝐑(t)\left[𝐪^{}\mathrm{\Delta }(t)\right]$$
(3.3)
is the solution to the classical equation of motion (2.2), satisfying the initial condition $`𝐪_{}(0;𝐪^{})=𝐪^{}`$. Consequently, the $`2N\times 2N`$ matrix $`𝐑(t)`$ satisfies the equation

$$\dot{𝐑}=\mathrm{\Sigma }\mathrm{𝐁𝐑}𝒜𝐑,$$
(3.4)
and the initial condition $`𝐑(0)=𝐈_{2N}`$. The vector $`\mathrm{\Delta }(t)`$ equals zero at $`t=0`$. For $`t>0`$ it is determined from the equations

$$\dot{\mathrm{\Delta }}=𝐑^1\mathrm{\Sigma }𝐂\mathrm{\Sigma }\stackrel{~}{𝐑}𝐂.$$
(3.5)
These two forms are equivalent due to the identities

$$\stackrel{~}{𝐑}(t)\mathrm{\Sigma }^1𝐑(t)\mathrm{\Sigma }^1,𝐑(t)\mathrm{\Sigma }\stackrel{~}{𝐑}(t)\mathrm{\Sigma },𝐑^1\mathrm{\Sigma }\stackrel{~}{𝐑}\mathrm{\Sigma }^1,$$
(3.6)
which follow from Eq. (3.4).

It is clear that the Heisenberg equation of motion for the operator $`\widehat{𝐪}`$ coincides with Eq. (2.2), provided $`𝐪`$ is replaced by $`\widehat{𝐪}`$. The solution to this equation is given by Eq. (3.3) with carets over $`𝐪`$ and $`𝐪^{}`$. Therefore, the identities (3.6) mean nothing but the conservation of the commutation relations (2.3) in time, or, in other words, the canonicity of transformation (3.3) and the unitarity of the evolution operator.

Now let us proceed to the averaging of the total Wigner function over the reservoir variables $`\xi `$. Our first assumption is that the initial total Wigner function is factorized:

$$W(𝐪,0)=W_0(𝐪)W_1(\xi ).$$
(3.7)
The second assumption concerns the initial Wigner function of the reservoir. We assume it to be Gaussian,

$$W_1(\xi )=\mathrm{}^M(det𝐅)^{1/2}\mathrm{exp}\left[\frac{1}{2}(\xi \gamma )𝐅^1(\xi \gamma )\right],$$
(3.8)
with some symmetric positive definite $`2M\times 2M`$ matrix F and a $`2M`$-vector $`\gamma `$ ($`M`$ being the number of degrees of freedom of the reservoir). In particular, it may correspond to a mixed equilibrium state  or to a pure squeezed coherent state . For any physically admissible Gaussian Wigner function the covariance matrix F must satisfy a set of conditions expressing generalized uncertainty relations. The simplest among them is the inequality  $`det𝐅(\mathrm{}^2/4)^M`$. Moreover, the parameter

$$\mu =(\mathrm{}/2)^M(det𝐅)^{1/2}1$$
(3.9)
characterizes “the degree of quantum mechanical purity” of the Gaussian state:

$$\mu =\text{Tr}\widehat{\rho }^2=W^2(\xi )\text{d}\xi /(2\pi \mathrm{})^M.$$
(3.10)
To calculate the averaged Wigner function

$$W_\xi (𝐐,t)=W(𝐐,\xi )\text{d}\xi /(2\pi \mathrm{})^^M$$
(3.11)
we split the matrix R and the vector $`\mathrm{\Delta }`$ into rectangular blocks in accordance with the decomposition $`𝐪=(𝐐,\xi )`$:

$$𝐑=\begin{array}{cc}𝐑_{11}& 𝐑_{12}\\ 𝐑_{21}& 𝐑_{22}\end{array},\mathrm{\Delta }=\begin{array}{c}\mathrm{\Delta }_𝐐\\ \mathrm{\Delta }_\xi \end{array}.$$
(3.12)
Then Eqs. (3.1), (3.2), (3.7) lead to the integral

$`W_\xi (𝐐,t)`$ $`=`$ $`{\displaystyle \delta \left(𝐐𝐑_{11}\left[𝐐^{}\mathrm{\Delta }_𝐐\right]𝐑_{12}\left[\xi ^{}\mathrm{\Delta }_\xi \right]\right)}`$ (3.13)
$`\times \delta \left(\xi 𝐑_{21}\left[𝐐^{}\mathrm{\Delta }_𝐐\right]𝐑_{22}\left[\xi ^{}\mathrm{\Delta }_\xi \right]\right)`$
$`\times W_0(𝐐^{})W_1(\xi ^{})\text{d}𝐐^{}\text{d}\xi ^{}\text{d}\xi /(2\pi \mathrm{})^M`$
The integration over $`\text{d}\xi `$ is trivial: it simply removes the second delta function. To perform the integration over $`\text{d}\xi ^{}`$ we replace the first delta function with its integral representation (recall that the dimension of the vector Q is $`2n`$),

$$\delta (𝐱)=e^{i\mathrm{𝐤𝐱}}\text{d}^{2n}𝐤/(2\pi )^{2n}.$$
(3.14)
Thus we obtain two Gaussian integrals (the first one over $`\text{d}\xi ^{}`$ and the second over $`\text{d}𝐤`$) that can be calculated exactly due to the well-known formula

$$_{\mathrm{}}^{\mathrm{}}\text{d}𝐱\mathrm{exp}(\mathrm{𝐱𝐀𝐱}+\mathrm{𝐛𝐱})=[det(𝐀/\pi )]^{1/2}\mathrm{exp}\left(\frac{1}{4}\mathrm{𝐛𝐀}^1𝐛\right).$$
(3.15)
Finally we arrive at the same equation (3.1), but with the variables $`𝐐,𝐐^{}`$ instead of $`𝐪,𝐪^{}`$, $`W_\xi (𝐐,t)`$ instead of $`W(𝐪,t)`$, and with the averaged propagator

$`G_\xi (𝐐,𝐐^{},t)`$ $`=`$ $`(2\pi )^n\left[det_{}(t)\right]^{1/2}`$ (3.16)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\left(𝐐𝐑_{11}𝐐^{}\delta _{}\right)_{}^1\left(𝐐𝐑_{11}𝐐^{}\delta _{}\right)\right],`$
where the symmetric matrix $`_{}(t)`$ equals

$$_{}(t)=𝐑_{12}(t)𝐅\stackrel{~}{𝐑}_{12}(t).$$
(3.17)
Recall that the dimension of the rectangular matrix $`𝐑_{12}`$ is $`2n\times 2M`$. The vector $`\delta _{}(t)`$ equals

$$\delta _{}(t)=𝐑_{12}(t)\gamma 𝐑_{11}(t)\mathrm{\Delta }_𝐐(t)𝐑_{12}(t)\mathrm{\Delta }_\xi (t).$$
(3.18)
Direct inspection shows that the propagator (3.16) (and, consequently, its convolution with any initial function) satisfies the Fokker-Planck equation (2.18) with the following drift matrix $`𝐀`$ and vector $`𝐊`$,

$$𝐀=\dot{𝐑}_{11}𝐑_{11}^1,$$
(3.19)
$$𝐊=\dot{\delta }_{}𝐀\delta _{}.$$
(3.20)
The matrix of diffusion coefficients reads

$$𝐃=\frac{1}{2}\left(\dot{}_{}𝐀_{}_{}\stackrel{~}{𝐀}\right)=\text{sym}\left[\left(\dot{𝐑}_{12}\dot{𝐑}_{11}𝐑_{11}^1𝐑_{12}\right)𝐅\stackrel{~}{𝐑}_{12}\right],$$
(3.21)
where we have introduced the notation

$$\text{sym}𝐀\frac{1}{2}\left(𝐀+\stackrel{~}{𝐀}\right).$$
(3.22)
Thus we have proved that the evolution of the Wigner function of any subsystem of a closed quadratic quantum system is governed by some effective Fokker-Planck equation, provided the initial state of the “remaining part of the system” is Gaussian. However, the coefficients of the equation obtained in this way depend, as a rule, on time even for a time-independent Hamiltonian of the closed system.

## 4 Two coupled oscillators

Let us illustrate the formulas of the preceding section, applying them to a very simple special model, where both the “subsystem under study” and the “reservoir” are the harmonic oscillators with a single degree of freedom. This example admits exact solutions, in contrast to the more realistic situation of a reservoir with a very large number of degrees of freedom, where one has to make various simplifications and approximations to obtain a closed result that would be easy to analyze. Consider first the free Hamiltonian

$$H_0=\frac{p_1^2}{2m_1}+\frac{1}{2}m_1\omega _1^2x_1^2+\frac{p_2^2}{2m_2}+\frac{1}{2}m_2\omega _2^2x_2^2$$
(4.1)
and the most general quadratic interaction Hamiltonian

$$H_{\text{int}}=g_{pp}p_1p_2+g_{px}p_1x_2+g_{xp}x_1p_2+g_{xx}x_1x_2.$$
(4.2)
We assume all the coefficients to be time-independent. In the notation of the preceding section we should write

$$𝐪=(p_1,x_1,p_2,x_2),𝐐=(p_1,x_1),\xi =(p_2,x_2).$$

Then the matrix $`𝒜=\mathrm{\Sigma }𝐁`$ (see Eq. (3.4)) reads

$$𝒜=\begin{array}{cccc}0& m_1\omega _1^2& g_{xp}& g_{xx}\\ m_1^1& 0& g_{pp}& g_{px}\\ g_{px}& g_{xx}& 0& m_2\omega _2^2\\ g_{pp}& g_{xp}& m_2^1& 0\end{array}.$$
(4.3)
Its characteristic equation turns out to be biquadratic:

$$det(𝒜i\omega 𝐈)=\omega ^4\left(\omega _1^2+\omega _2^2+2\mathrm{\Delta }\right)\omega ^2+\omega _1^2\omega _2^2+\mathrm{\Delta }^2g=0,$$
(4.4)
where

$$\mathrm{\Delta }=det𝒜_{12}=g_{pp}g_{xx}g_{px}g_{xp},$$
(4.5)
$$g=\frac{g_{xx}^2}{m_1m_2}+\frac{m_1}{m_2}\omega _1^2g_{px}^2+\frac{m_2}{m_1}\omega _2^2g_{xp}^2+m_1m_2\omega _1^2\omega _2^2g_{pp}^2.$$
(4.6)
The solutions of Eq. (4.4) are as follows:

$`\omega _\pm `$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}([{\displaystyle \frac{1}{2}}(\omega _1^2+\omega _2^2)+(\omega _1^2\omega _2^2+\mathrm{\Delta }^2g)^{1/2}+\mathrm{\Delta }]^{\frac{1}{2}}`$
$`\pm [{\displaystyle \frac{1}{2}}(\omega _1^2+\omega _2^2)(\omega _1^2\omega _2^2+\mathrm{\Delta }^2g)^{1/2}+\mathrm{\Delta }]^{\frac{1}{2}})`$
(two other solutions are equal to $`\omega _\pm `$).

Thus, the evolution of two coupled harmonic oscillators can be described explicitly for quite arbitrary quadratic interaction Hamiltonians with time-independent coefficients. If both frequences given by Eq. (4) are real, then the particles perform harmonic oscillations. For certain parameters complex normal frequencies are possible. Then the motion becomes aperiodic. However, since any normal frequency $`\omega _+`$ or $`\omega _{}`$ is accompanied by the frequency with the opposite sign, it is impossible to obtain damped oscillations of either particle. The coordinate and momentum of any oscillator will increase with time.

To illustrate this statement, let us consider first the case (which seems the most natural) of the interaction via the coordinates, where the only nonzero coefficient in Eq. (4.2) is $`g_{xx}`$. Then $`\mathrm{\Delta }=0`$, and for sufficiently strong coupling, when $`g>\omega _1^2\omega _2^2`$, we have the real frequency $`\omega _+`$ and the pure imaginary frequency $`\omega _{}`$:

$$\omega _+\omega =\left\{\left[g+\frac{1}{4}\left(\omega _1^2\omega _2^2\right)^2\right]^{1/2}+\frac{1}{2}\left(\omega _1^2+\omega _2^2\right)\right\}^{\frac{1}{2}},$$
(4.7)
$$\omega _{}i\lambda ,\lambda =\left\{\left[g+\frac{1}{4}\left(\omega _1^2\omega _2^2\right)^2\right]^{1/2}\frac{1}{2}\left(\omega _1^2+\omega _2^2\right)\right\}^{\frac{1}{2}}.$$
(4.8)
Solving Eq. (3.4), we obtain the following formula for the matrix $`𝐑_{11}(t)`$:

$$𝐑_{11}=\begin{array}{cccc}\rho _+\mathrm{cos}\omega t+\rho _{}\mathrm{cosh}\lambda t& & m_1\left(\omega \rho _+\mathrm{sin}\omega t+\lambda \rho _{}\mathrm{sinh}\lambda t\right)& \\ & & & \\ \frac{1}{m_1}\left(\frac{\rho _+}{\omega }\mathrm{sin}\omega t+\frac{\rho _{}}{\lambda }\mathrm{sinh}\lambda t\right)& & \rho _+\mathrm{cos}\omega t+\rho _{}\mathrm{cosh}\lambda t& \end{array},$$
(4.9)
where

$$\rho _\pm =\frac{1}{2}\pm \frac{1}{4}\left(\omega _1^2\omega _2^2\right)\left[g+\frac{1}{4}\left(\omega _1^2\omega _2^2\right)^2\right]^{\frac{1}{2}}.$$
(4.10)
We see that the oscillations of the first particle actually increase. Its energy is derived from the nonpositive definite potential energy of the whole system, since the whole system turns out to be unstable when $`g>\omega _1^2\omega _2^2`$. Note that Eq. (4.9) still holds even if $`g_{xp}0`$. If $`g_{xx}=g_{xp}=0`$, but $`g_{pp}0`$ and $`g_{px}0`$, we obtain, instead of Eq. (4.9), the expression

$$𝐑_{11}=\begin{array}{cccc}\rho _+\mathrm{cos}\omega t+\rho _{}\mathrm{cosh}\lambda t& m_1\omega _1^2\left(\frac{\rho _+}{\omega }\mathrm{sin}\omega t+\frac{\rho _{}}{\lambda }\mathrm{sinh}\lambda t\right)& & \\ & & & \\ \frac{1}{m_1\omega _1^2}\left(\omega \rho _+\mathrm{sin}\omega t+\lambda \rho _{}\mathrm{sinh}\lambda t\right)& \rho _+\mathrm{cos}\omega t+\rho _{}\mathrm{cosh}\lambda t& & \end{array},$$
(4.11)
with the same values of $`\omega `$ and $`\lambda `$ (provided, of course, $`g>\omega _1^2\omega _2^2`$).

Both frequencies $`\omega _+`$ and $`\omega _{}`$ have nonzero imaginary parts provided the argument of the first square bracket in Eq. (4) is positive while the argument of the second square bracket is negative. Then the following inequality must hold:

$$\frac{1}{4}\left(\omega _1^2\omega _2^2\right)^^2+\mathrm{\Delta }\left(\omega _1^2+\omega _2^2\right)+g<0.$$
(4.12)
Assuming for the sake of simplicity that $`\omega _1=\omega _2=\omega _0`$, we can rewrite it as

$$\left(m_1m_2\right)^1\left[\left(m_1g_{xx}+m_2\omega _0^2g_{pp}\right)^2+\omega _0^2\left(m_1g_{px}m_2g_{xp}\right)^2\right]<0.$$
(4.13)
If $`\omega _0^2>0`$, the only possibility for fulfilling this inequality is to assume that the masses have opposite signs. Of course, such a system is unstable, and it can be considered only as an extremely simplified model. Nonetheless, sometimes the models of this sort were considered . So we discuss briefly this case as well. Suppose for simplicity that $`m_1=1`$, $`m_2=1`$. Then $`g<0`$. The square roots in Eq. (4) can be extracted if $`|g|=2\omega _0^2|\mathrm{\Delta }|`$. In this case the imaginary parts of $`\omega _\pm `$ differ from zero, provided $`\mathrm{\Delta }<0`$. Such a situation holds for the following relations between the coupling coefficients:

$$g_{xp}=g_{px},g_{xx}=\omega _0^2g_{pp}.$$
(4.14)
Then

$$\omega _\pm =\omega _0\pm i\gamma ,\gamma =|\mathrm{\Delta }|^{1/2}=\left(g_{px}^2+\omega _0^2g_{pp}^2\right)^{1/2}.$$
(4.15)
The explicit form of the matrix $`𝐑_{11}`$ reads

$$𝐑_{11}=\mathrm{cosh}\gamma t\begin{array}{ccc}\mathrm{cos}\omega _0t& & \omega _0\mathrm{sin}\omega _0t\\ \omega _0^1\mathrm{sin}\omega _0t& & \mathrm{cos}\omega _0t\end{array}.$$
(4.16)
Consequently, the amplitude of the oscillations increases without bound. The drift matrix $`𝐀`$ (3.19) depends on time as follows:

$$𝐀=\begin{array}{cc}\gamma \mathrm{tanh}\gamma t& \omega _0^2\\ 1& \gamma \mathrm{tanh}\gamma t\end{array}.$$
(4.17)
An interesting model described by the Hamiltonian

$$H=p_1p_2+\omega _0^2x_1x_2+\gamma \left(x_2p_2x_1p_1\right)$$
(4.18)
was proposed for the first time by Bateman . The Lagrangian

$$L=\dot{x}_1\dot{x}_2\left(\omega _0^2+\gamma ^2\right)x_1x_2+\gamma \left(x_1\dot{x}_2\dot{x}_1x_2\right)$$
(4.19)
was considered by Morse and Feshbach . This system was investigated, e.g., in  (see also ). The equations of motion in this case read

$$\ddot{x}_1+2\gamma \dot{x}_1+\left(\omega _0^2+\gamma ^2\right)x_1=0,\ddot{x}_2+2\gamma \dot{x}_2+\left(\omega _0^2+\gamma ^2\right)x_2=0.$$
(4.20)
Consequently, at the classical level we have damping in the first mode and amplification in the second mode. Moreover, at the classical level both particles are completely independent due to Eq. (4.20). The quantum picture is more complicated, since the quantum behaviour is governed not by the second-order equations of motion, but by the Hamiltonian (4.18), in which the dynamical variables of both particles are entangled. (The nonunique correspondence between the equations of motion and the Lagrangians or Hamiltonians leading to them, as well as related ambiguities of quantization, were investigated in detail in . It can be proved  that no Hamiltonian leading to Eq. (4.20) and coinciding with Eq. (4.1) for $`\gamma =0`$ exists.) In particular, the velocity of each particle is determined by the generalized momentum of the other  particle:

$$\dot{x}_1=p_2\gamma x_1,\dot{x}_2=p_1+\gamma x_2.$$

The two other Hamilton equations read

$$\dot{p}_1=\gamma p_1\omega _0^2x_2,\dot{p}_2=\gamma p_2\omega _0^2x_1.$$

In terms of the blocks of the matrix $`𝐑`$ these equations can be rewritten as follows:

$$\dot{𝐑}_{1k}=\mathrm{𝐚𝐑}_{1k}+\mathrm{𝐛𝐑}_{2k},\dot{𝐑}_{2k}=\mathrm{𝐛𝐑}_{1k}\mathrm{𝐚𝐑}_{2k},(k=1,2)$$
(4.21)
where a and b are $`2\times 2`$ matrices,

$$𝐚=\begin{array}{cc}\gamma & 0\\ 0& \gamma \end{array},𝐛=\begin{array}{cc}0& \omega _0^2\\ 1& 0\end{array}.$$

Eliminating the matrix $`𝐑_{2k}`$ we obtain the second-order equation

$$\ddot{𝐑}_{1k}2𝐚\dot{𝐑}_{1k}+\left(\gamma ^2+\omega _0^2\right)𝐑_{1k}=0.$$
(4.22)
Seeking its solution in the form $`𝐑_{1k}=\mathrm{exp}(\mathrm{\Lambda }t)𝐑_0`$ we obtain the characteristic equation

$$\mathrm{\Lambda }^22𝐚\mathrm{\Lambda }+\left(\gamma ^2+\omega _0^2\right)𝐈=0,$$

whose solution reads

$$\mathrm{\Lambda }=𝐚\pm i\omega _0𝐈.$$

Taking into account the initial conditions

$$𝐑_{11}(0)=𝐈,\dot{𝐑}_{11}(0)=𝐚,𝐑_{12}(0)=0,\dot{𝐑}_{12}(0)=𝐛,$$

we obtain finally the matrices

$$𝐑_{11}=\mathrm{cos}\omega _0t\begin{array}{cc}e^{\gamma t}& 0\\ 0& e^{\gamma t}\end{array},𝐑_{12}=\mathrm{sin}\omega _0t\begin{array}{cc}0& \omega _0e^{\gamma t}\\ \omega _0^1e^{\gamma t}& 0\end{array}.$$

Consequently, after averaging over the state of the second particle we have a system with a damped coordinate, but with a momentum increasing in time without bound. Moreover, the momentum will no longer be related to the velocity. The drift matrix (3.19) in the Fokker-Planck equation for the averaged Wigner function of the first particle depends on time as follows:

$$𝐀=\begin{array}{cc}\omega _0\mathrm{tan}(\omega _0t)+\gamma & 0\\ 0& \omega _0\mathrm{tan}(\omega _0t)\gamma \end{array}.$$

We see that the reduction of “Bateman’s mirror model” (4.18) does not lead to a damped quantum oscillator in the conventional meaning of this term.

## 5 Oscillator in a thermostat: weak coupling limit

We now proceed to a more realistic model, where the oscillator under study (its frequency will be denoted by $`\omega _0`$) is coupled to a large number of other oscillators with frequencies $`\omega _i`$. This model was the subject of investigations in many papers: see, e.g., -. A more comprehensive reference list can be found in the review . We assume that each oscillator is described by the standard Hamiltonian

$$H_i=\frac{1}{2}\left(p_i^2+\omega _i^2x_i^2\right)$$

with unit mass (this can easily be achieved by rescaling the coordinates), while the quadratic interaction Hamiltonian is chosen in the most general form:

$$H_{\text{int}}=\underset{i}{}\left(z_ip_ip_0+v_ip_ix_0+u_ix_ip_0+g_ix_ix_0\right).$$
(5.1)
The coupling constants $`z_i`$, $`v_i`$, $`u_i`$, $`g_i`$ are assumed to be time independent.

Equation (3.4) and its initial condition are equivalent to the following equations and initial conditions for the blocks of the matrix $`𝐑`$ of (3.12):

$$\dot{𝐑}_{11}=𝒜_{11}𝐑_{11}+𝒜_{12}𝐑_{21},𝐑_{11}(0)=𝐈,$$
(5.2)
$$\dot{𝐑}_{21}=𝒜_{21}𝐑_{11}+𝒜_{22}𝐑_{21},𝐑_{21}(0)=0,$$
(5.3)
$$\dot{𝐑}_{12}=𝒜_{11}𝐑_{12}+𝒜_{12}𝐑_{22},𝐑_{12}(0)=0,$$
(5.4)
$$\dot{𝐑}_{22}=𝒜_{21}𝐑_{12}+𝒜_{22}𝐑_{22},𝐑_{22}(0)=𝐈.$$
(5.5)
The $`2\times 2`$ matrix $`𝒜_{11}`$ and the $`2M\times 2M`$ matrix $`𝒜_{22}`$ ($`M`$ is the number of oscillators in the reservoir) read

$$𝒜_{11}=\begin{array}{cc}0& \omega _0^2\\ 1& 0\end{array},\begin{array}{cc}0& \text{diag}(\omega _1^2,\mathrm{},\omega _i^2,\mathrm{})\\ \text{diag}(1,\mathrm{},1,\mathrm{})& 0\end{array}.$$
(5.6)
The matrices $`𝒜_{12}`$ and $`𝒜_{21}`$ are rectangular with the dimensions $`2\times 2M`$ and $`2M\times 2`$, respectively:

$$𝒜_{12}=\begin{array}{cccccccc}v_1& \mathrm{}& v_i& \mathrm{}& g_1& \mathrm{}& g_i& \mathrm{}\\ z_1& \mathrm{}& z_i& \mathrm{}& u_1& \mathrm{}& u_i& \mathrm{}\end{array},$$
(5.7)
$$𝒜_{21}=\begin{array}{cc}u_1& g_1\\ \mathrm{}& \mathrm{}\\ u_i& g_i\\ \mathrm{}& \mathrm{}\\ z_1& v_1\\ \mathrm{}& \mathrm{}\\ z_i& v_i\\ \mathrm{}& \mathrm{}\end{array}.$$
(5.8)
In this section we consider the case where all the elements of the interaction matrices $`𝒜_{12}`$ and $`𝒜_{21}`$ are small. Then we may use perturbation theory. In the zeroth approximation we obtain from Eq. (5.2)

$$𝐑_{11}^{(0)}(t)=\mathrm{exp}\left(𝒜_{11}t\right).$$
(5.9)
Putting this expression into the right-hand side of Eq. (5.3) we obtain the first-order solution for the matrix $`𝐑_{21}`$,

$$𝐑_{21}^{(1)}(t)=\mathrm{exp}\left(𝒜_{22}t\right)_0^t\mathrm{exp}\left(𝒜_{22}\tau \right)𝒜_{21}\mathrm{exp}\left(𝒜_{11}\tau \right)\text{d}\tau .$$
(5.10)
Then Eqs. (3.19) and (5.2) lead to the following first-order approximation for the matrix $`𝐀`$ governing the evolution of the average values of the subsystem variables according to Eq. (2.18):

$$𝐀^{(1)}=𝒜_{11}+𝒜_{12}𝐑_{21}^{(1)}(t)\left[𝐑_{11}^{(0)}(t)\right]^1.$$
(5.11)
Taking into account Eqs. (5.9) and (5.10) and making the change of variable $`t\tau =x`$ in the integrand we arrive at the formula

$$\mu 𝐀𝒜_{11}=𝒜_{12}_0^t\mathrm{exp}\left(𝒜_{22}x\right)𝒜_{21}\mathrm{exp}\left(𝒜_{11}x\right)\text{d}x.$$
(5.12)
The solutions to Eqs. (5.4) and (5.5) in the same approximation read

$$𝐑_{22}^{(0)}(t)=\mathrm{exp}\left(𝒜_{22}t\right),$$
(5.13)
$$𝐑_{12}^{(1)}(t)=\mathrm{exp}\left(𝒜_{11}t\right)_0^t\mathrm{exp}\left(𝒜_{11}\tau \right)𝒜_{12}\mathrm{exp}\left(𝒜_{22}\tau \right)\text{d}\tau .$$
(5.14)
Then the diffusion matrix of the Fokker-Planck equation (2.18), due to Eqs. (3.21) and (5.4), reads

$`𝐃`$ $`=`$ $`\text{sym}\left(𝒜_{12}𝐑_{22}^{(0)}(t)𝐅\stackrel{~}{𝐑}_{12}^{(1)}(t)\right)`$ (5.15)
$`=`$ $`\text{sym}\left\{𝒜_{12}\mathrm{exp}\left(𝒜_{22}t\right)𝐅{\displaystyle _0^t}\mathrm{exp}\left(\stackrel{~}{𝒜}_{22}\tau \right)\stackrel{~}{𝒜}_{12}\mathrm{exp}\left(\stackrel{~}{𝒜}_{11}\tau \right)\text{d}\tau \mathrm{exp}\left(\stackrel{~}{𝒜}_{11}t\right)\right\}.`$
The explicit forms of the matrices $`𝐑_{11}^{(0)}(t)`$ and $`𝐑_{22}^{(0)}(t)`$ are as follows:

$$\mathrm{exp}\left(𝒜_{11}t\right)=\begin{array}{cc}\mathrm{cos}\omega _0t& \omega _0\mathrm{sin}\omega _0t\\ \omega _0^1\mathrm{sin}\omega _0t& \mathrm{cos}\omega _0t\end{array},$$
(5.16)
$$\mathrm{exp}\left(𝒜_{22}t\right)=\begin{array}{cc}\text{diag}(\mathrm{cos}\omega _it)& \text{diag}(\omega _i\mathrm{sin}\omega _it)\\ \text{ diag}(\omega _i^1\mathrm{sin}\omega _it)& \text{diag}(\mathrm{cos}\omega _it)\end{array},$$
(5.17)
If we choose the thermostat variance matrix to be

$$𝐅=\begin{array}{cc}\text{diag}\left(\omega _i^2f_i\right)& 0\\ 0& \text{diag}\left(f_i\right)\end{array},$$
(5.18)
(in particular, $`𝐅`$ may be an equilibrium variance matrix for the thermostat variables), then it describes the steady-state solution of Eq. (2.14) in the absense of the interaction:

$$𝐅(t)=\mathrm{exp}\left(𝒜_{22}t\right)𝐅\mathrm{exp}\left(\stackrel{~}{𝒜}_{22}t\right)=𝐅.$$
(5.19)
Consequently, one may interchange the matrices $`\mathrm{exp}\left(𝒜_{22}t\right)`$ and $`𝐅`$ in Eq. (5.15) in accordance with the formula

$$\mathrm{exp}\left(𝒜_{22}t\right)𝐅=𝐅\mathrm{exp}\left(\stackrel{~}{𝒜}_{22}t\right).$$
(5.20)
Then Eq. (5.15) assumes the form ($`x=t\tau `$)

$$𝐃=\text{sym}\left\{_0^t\mathrm{exp}\left(𝒜_{11}x\right)𝒜_{12}\mathrm{exp}\left(𝒜_{22}x\right)\text{d}x𝐅\stackrel{~}{𝒜}_{12}\right\}.$$
(5.21)
The explicit expressions for the matrix elements of $`2\times 2`$ matrix $`\mu `$ of (5.12) are as follows:

$$\mu _{11}=\frac{1}{2}\underset{i}{}\left[\mathrm{\Delta }_iS_i^{(+)}\omega _0^1G_iS_i^{()}+\kappa _iC_i^{(+)}\right],$$
(5.22)
$$\mu _{12}=\frac{1}{2}\underset{i}{}\left[\omega _0\kappa _iS_i^{()}+\omega _0\mathrm{\Delta }_iC_i^{()}+G_iC_i^{(+)}\right],$$
(5.23)
$$\mu _{21}=\frac{1}{2}\underset{i}{}\left[\omega _0^1\kappa _iS_i^{()}\omega _0^1\mathrm{\Delta }_iC_i^{()}Z_iC_i^{(+)}\right],$$
(5.24)
$$\mu _{22}=\frac{1}{2}\underset{i}{}\left[\mathrm{\Delta }_iS_i^{(+)}\omega _0Z_iS_i^{()}\kappa _iC_i^{(+)}\right].$$
(5.25)
We have introduced the notation

$$S_i^{(\pm )}=\frac{\mathrm{sin}(\omega _i\omega _0)t}{\omega _i\omega _0}\pm \frac{\mathrm{sin}(\omega _i+\omega _0)t}{\omega _i+\omega _0},$$
(5.26)
$$C_i^{(\pm )}=\frac{1\mathrm{cos}(\omega _i\omega _0)t}{\omega _i\omega _0}\pm \frac{1\mathrm{cos}(\omega _i+\omega _0)t}{\omega _i+\omega _0}.$$
(5.27)
Other parameters are the bilinear combinations of the coupling constants:

$$\begin{array}{ccc}\mathrm{\Delta }_i=g_iz_iu_iv_i,& & \kappa _i=\omega _iz_iv_i+g_iu_i/\omega _i,\\ G_i=\omega _iv_i^2+g_i^2/\omega _i,& & Z_i=\omega _iz_i^2+u_i^2/\omega _i.\end{array}$$
(5.28)
Eq. (5.21) results in the following expressions for the elements of the diffusion matrix $`𝐃`$:

$$D_{11}=\frac{1}{2}\underset{i}{}\omega _if_i\left[G_iS_i^{(+)}+\omega _0\mathrm{\Delta }_iS_i^{()}\omega _0\kappa _iC_i^{()}\right],$$
(5.29)
$$D_{22}=\frac{1}{2}\underset{i}{}\omega _if_i\left[Z_iS_i^{(+)}+\omega _0^1\mathrm{\Delta }_iS_i^{()}+\omega _0^1\kappa _iC_i^{()}\right],$$
(5.30)
$$D_{12}=\frac{1}{2}\underset{i}{}\omega _if_i\left[\kappa _iS_i^{(+)}+(2\omega _0)^1\left(\omega _0^2Z_iG_i\right)C_i^{()}\right].$$
(5.31)
We see that in the general case both the drift matrix and the diffusion matrix have rather complicated time dependences. However, under certain conditions the formulas can be simplified, if we proceed to the continuum limit. This means that we assume the number of oscillators in the reservoir to be very large and the frequencies $`\omega _i`$ to be so close to each other that we may replace the sums by integrals over $`\text{d}\omega _i\text{d}\omega `$. Then we need to calculate integrals of the following type:

$$\sigma _\pm =\phi (\omega )\frac{\mathrm{sin}(\omega \pm \omega _0)t}{\omega \pm \omega _0}\text{d}\omega ,$$
(5.32)
$$\sigma _\pm =\phi (\omega )\frac{1\mathrm{cos}(\omega \pm \omega _0)t}{\omega \pm \omega _0}\text{d}\omega .$$
(5.33)
From physical considerations it is clear that the elements of the matrices $`\mu `$ and $`𝐃`$ are determined mainly by the terms in Eqs. (5.22)-(5.25) and (5.29)-(5.31) that correspond to frequencies near $`\omega _0`$, since the most effective interactions between the oscillator under study and the thermostat oscillators take place under resonance condition. Indeed, for a sufficiently smooth function $`\phi (\omega )`$ and for $`t\omega _0^1`$, only points belonging to the domain $`|\omega \omega _0|t^1`$ make a significant contribution to the integrals (5.32) and (5.33), due to the rapid oscillations of the trigonometric functions outside this domain. Thus, assuming that $`\omega =\omega _0`$ in all the functions except $`\mathrm{sin}(\omega +\omega _0)t`$, we may evaluate the integral $`\sigma _+`$ as follows:

$$\sigma _+\frac{\phi (\omega _0)}{2\omega _0}_{\omega _1}^{\omega _2}\mathrm{sin}(\omega +\omega _0)t\text{d}\omega =\frac{\phi (\omega _0)}{2\omega _0t}\left[\mathrm{cos}(\omega _1+\omega _0)t\mathrm{cos}(\omega _2+\omega _0)t\right].$$

Consequently, we may neglect the value of $`\sigma _+`$ at $`\omega _0t1`$. The same is true for the integral

$$\phi (\omega )\frac{\mathrm{cos}(\omega +\omega _0)t}{\omega +\omega _0}\text{d}\omega 𝒪\left(\frac{1}{\omega _0t}\right).$$

As concerns the integral $`\sigma _{}`$, its value does not depend on time for $`t\omega _0^1`$, since making the substitutions $`x=\omega \omega _0`$, $`y=xt`$, we obtain

$$\sigma _{}=\phi (\omega _0)\frac{\mathrm{sin}xt}{x}\text{d}x=\phi (\omega _0)\frac{\mathrm{sin}y}{y}\text{d}y=\pi \phi (\omega _0).$$
(5.34)
To evaluate the integral $`\phi (\omega )\frac{\mathrm{cos}(\omega \omega _0)t}{\omega \omega _0}\text{d}\omega `$ we use the Taylor expansion $`\phi (\omega )=\phi (\omega _0)+\phi ^{}(\omega _0)(\omega \omega _0)+\mathrm{}`$, make the substitution $`x=\omega \omega _0`$, and expand the limits of integration from $`\mathrm{}`$ to $`\mathrm{}`$. Then the first integral vanishes, because the function $`\mathrm{cos}(xt)/x`$ is odd. The second integral decreases at least as $`1/t`$ when $`t\mathrm{}`$. Therefore, we may assume that the integrals $`\delta _\pm `$ do not depend on time at $`\omega _0t1`$:

$$\delta _\pm =\frac{\phi (\omega )}{\omega \omega _0}\text{d}\omega .$$
(5.35)
In the case of the “minus” sign the principal value of the integral is implied: it is designated with the symbol $`^{}`$. We arrive at the following expressions for the matrix elements of the matrices $`\mu `$ and $`𝐃`$ in the continuum limit:

$$\mu _{11}=\frac{1}{2}\pi \nu (\omega _0)\left[\mathrm{\Delta }(\omega _0)+\omega _0^1G(\omega _0)\right]+^{}\frac{\omega \kappa (\omega )\nu (\omega )}{\omega ^2\omega _0^2}\text{d}\omega ,$$
(5.36)
$$\mu _{12}=\frac{1}{2}\pi \nu (\omega _0)\omega _0\kappa (\omega _0)+^{}\frac{\nu (\omega )}{\omega ^2\omega _0^2}\left[\omega _0^2\mathrm{\Delta }(\omega )+\omega G(\omega )\right]\text{d}\omega ,$$
(5.37)
$$\mu _{21}=\frac{1}{2}\pi \nu (\omega _0)\omega _0^1\kappa (\omega _0)^{}\frac{\nu (\omega )}{\omega ^2\omega _0^2}\left[\mathrm{\Delta }(\omega )+\omega Z(\omega )\right]\text{d}\omega ,$$
(5.38)
$$\mu _{22}=\frac{1}{2}\pi \nu (\omega _0)\left[\mathrm{\Delta }(\omega _0)+\omega _0Z(\omega _0)\right]^{}\frac{\omega \kappa (\omega )\nu (\omega )}{\omega ^2\omega _0^2}\text{d}\omega ,$$
(5.39)
$$D_{11}=\frac{1}{2}\pi \nu (\omega _0)f(\omega _0)\omega _0\left[\omega _0\mathrm{\Delta }(\omega _0)+G(\omega _0)\right]\omega _0^2^{}\frac{\omega f(\omega )\kappa (\omega )\nu (\omega )}{\omega ^2\omega _0^2}\text{d}\omega ,$$
(5.40)
$$D_{22}=\frac{1}{2}\pi \nu (\omega _0)f(\omega _0)\left[\mathrm{\Delta }(\omega _0)+\omega _0Z(\omega _0)\right]+^{}\frac{\omega f(\omega )\kappa (\omega )\nu (\omega )}{\omega ^2\omega _0^2}\text{d}\omega ,$$
(5.41)
$$D_{12}=\frac{1}{2}\pi \nu (\omega _0)\omega _0f(\omega _0)\kappa (\omega _0)+^{}\frac{\omega f(\omega )\nu (\omega )}{\omega ^2\omega _0^2}\left[\omega _0^2Z(\omega )G(\omega )\right]\text{d}\omega .$$
(5.42)
Here $`\nu (\omega )`$ is the density of states function, whereas the functions $`Z(\omega )`$, $`G(\omega )`$, etc. are obvious generalizations of functions defined in Eq. (5.28).

In principle, the frequency dependences of the coupling constants can be chosen in such a way that the integrals in Eqs. (5.36)-(5.42) vanish. For instance, this is possible provided the corresponding combinations of functions $`\kappa `$, $`\nu `$, $`f`$, $`G`$, $`Z`$, $`\mathrm{\Delta }/\omega `$, understood as functions of the argument $`x=\omega ^2`$, do not change their values under the reflection in the point $`x_0=\omega _0^2`$, and these functions (or at least the density of states) decrease sufficiently rapidly with distance from the point $`x_0`$ in both directions. In such a case all the coefficients $`\mu _{ik}`$ and $`D_{ik}`$ are determined by the values of the aforementioned functions at the point $`\omega _0`$. Moreover, the diffusion coefficients are proportional to the corresponding elements of the matrix $`\mu `$:

$$D_{11}=\omega _0^2f_0\mu _{11},D_{22}=f_0\mu _{22},D_{12}=f_0\mu _{12}=\omega _0^2f_0\mu _{21}.$$
(5.43)
We see that in the continuous weak coupling limit the reduced Wigner function of the oscillator obeys (under certain conditions) the Fokker-Planck equation (2.18) with time-independent coefficients at times $`t\omega _0^1`$. Let us check, however, whether condition (2.27) is fulfilled. Since $`\text{Tr}𝐀=\text{Tr}\mu `$, then due to Eq. (5.43) we must check the inequality

$$\omega _0^2f_0^2\mu _{11}\mu _{22}f_0^2\mu _{12}^2\mathrm{}^2\left(\mu _{11}+\mu _{22}\right)^2/16.$$
(5.44)
Taking into account Eqs. (5.36)-(5.39) and Eq. (5.28), we arrive at the inequality (all the functions are taken at the point $`\omega =\omega _0`$)

$$(4f\omega /\mathrm{})^2\left[2\mathrm{\Delta }^2+\mathrm{\Delta }(\omega Z+G/\omega )\right]\left[2\mathrm{\Delta }+(\omega Z+G/\omega )\right]^2.$$
(5.45)
It cannot be satisfied for an arbitrary choice of coupling constants. For instance, it is violated if $`\mathrm{\Delta }=0`$, or when any three of four coefficients $`g`$, $`z`$, $`u`$, $`v`$ vanish.

This result seems paradoxical. Indeed, we started from the exact equation of motion for the density matrix of a closed system, found the exact solution to this equation, and after this we performed averaging over the thermostat degrees of freedom. Since the laws of quantum mechanics were not violated at any step, the reduced density matrix must be positive definite at any time for quite arbitrary coupling constants. On the other hand, if, for instance, $`z_i=u_i=v_i=0`$, then inequality (5.45) does not hold. Hence, following the reasonings given in Sec. 2, we could obtain a nonpositive definite density matrix in the process of evolution!

This apparent contradiction is resolved in the following way. Inequality (2.27) is a necessary and sufficient condition ensuring that any density matrix that was positive definite at any instant of time will remain positive definite at all subsequent moments. But in the problem under study we have the selected instant $`t=0`$: this is just the moment when the interaction with the thermostat was turned on. Since in the presence of the interaction with the environment the evolution of the oscillator density matrix is nonunitary, the set of density matrices $`\rho (t)`$ arising from all initially admissible density matrices $`\rho (0)`$ does not coincide with the set of all admissible density matrices. In particular, correct initial density matrices cannot  turn into the specific ones that could become nonpositive definite at some instant of time in the case of violation of inequality (5.45). Therefore, there is no need to check conditions like (2.27) or (5.44), (5.45) when the density matrix of the subsystem is obtained by reduction of the exact density matrix of the closed system: the reduced density matrix (calculated with the proper accuracy) turns out to be positive definite automatically.

However, if the goal is the derivation of a self-consistent Fokker-Planck equation on the basis of a “microscopic” model of the oscillator interacting with a large reservoir, then we have to recognize that the underlying “microscopic” model cannot be quite arbitrary: its parameters must satisfy rather strong restriction (5.45), in order to prevent the appearance of the unphysical solutions when this equation is applied to arbitrary initial states.

Returning to the analysis of Eqs. (5.36)-(5.43), we notice that the matrix

$$𝐅_0=\begin{array}{cc}\omega _0^2f_0& 0\\ 0& f_0\end{array}$$
(5.46)
satisfies, due to Eq. (5.43), the relation

$$\mathrm{𝐀𝐅}_0+𝐅_0\stackrel{~}{𝐀}+2𝐃=0.$$
(5.47)
This means that $`𝐅_0`$ is the steady state solution to Eq. (2.19) for the oscillator variance matrix, independently of the concrete values of the coefficients of the drift matrix $`\mu _{ik}`$. In particular, if

$$f_0=\frac{\mathrm{}}{2\omega _0}\mathrm{coth}\left(\frac{\mathrm{}\omega _0}{2kT}\right),$$
(5.48)
then the matrix (5.46) coincides with the equilibrium variance matrix of the oscillator. Consequently, the steady state solution of the Fokker-Planck equation with the coefficients (5.36)-(5.42) is the equilibrium Wigner function of the oscillator, if the thermostat itself is in the equilibrium state. Moreover, we may imagine a situation with

$$f_0=\frac{\mathrm{}}{2\omega _0}\underset{i=0}{\overset{n}{}}\alpha _i\mathrm{coth}\left(\frac{\mathrm{}\omega _0}{2kT_i}\right),\underset{i=1}{\overset{n}{}}\alpha _i=1,\alpha _i0,$$
(5.49)
Then, the oscillator under study exhibits relaxation to the Gaussian steady state with the variance matrix given by Eqs. (5.46) and (5.49). Such a situation may be realized when the thermostat consists of several large independent subsystems possessing their own temperatures (see, e.g., ). In this case $`n`$ is the number of subsystems, and $`\alpha _i`$ is the “weight” of each subsystem.

We see that the final steady state of the oscillator weakly interacting with a thermostat does not depend on the concrete values of the coupling constants, provided the thermostat was initially in an equilibrium or quasi-equilibrium (described by Eq. (5.49) state. Now let us analyze possible forms of the drift matrix $`𝐀`$. Due to the property $`\mu _{12}=\omega _0^2\mu _{21}`$ its characteristic equation $`det(𝐀\lambda 𝐈)=0`$ reads

$$\lambda ^2\left(\mu _{11}+\mu _{22}\right)\lambda +\omega _0^2+\mu _{11}\mu _{22}\mu _{12}\mu _{21}=0.$$
(5.50)
The solutions to this equations can be written as $`\lambda _{1,2}=\gamma \pm i\omega _{}`$, with

$$\omega _{}=\left[\omega _0^2\frac{1}{4}\left(\mu _{11}\mu _{22}\right)^2\mu _{12}\mu _{21}\right],$$
(5.51)
$$\gamma =\frac{1}{2}\text{Tr}\mu =\frac{1}{2}\left(\mu _{11}+\mu _{22}\right)=\frac{1}{4}\pi \nu (\omega _0)\left[2\mathrm{\Delta }_0+\omega _0Z_0+\omega _0^1G_0\right].$$
(5.52)
The last formula holds even when the integrals in Eqs. (5.36) and (5.39) are not equal to zero. Assuming the damping coefficient $`\gamma `$ to be positive and comparing Eqs. (5.52) and (5.45), we may rewrite inequality (5.45) as follows (recall that all functions are taken at $`\omega =\omega _0`$):

$$(4f\omega /\mathrm{})^2\mathrm{\Delta }2\mathrm{\Delta }+(\omega Z+G/\omega )=4\gamma /[\pi \nu (\omega )].$$
(5.53)
In particular, for an equilibrium thermostat at zero temperature (when $`f=\mathrm{}/2\omega `$) we obtain the relation

$$2\mathrm{\Delta }\omega Z+G/\omega ,$$
(5.54)
which due to Eq. (5.28) turns into the inequality

$$(u+v)^2+(\omega zg/\omega )^20.$$
(5.55)
Thus we arrive at a striking conclusion: the microscopic model of the interaction between the oscillator and the equilibrium reservoir at zero temperature, based on the general quadratic Hamiltonian (5.1), leads to the Fokker-Planck equation with time independent coefficients, valid for any physically admissible initial states of the oscillator, if and only if the coupling constants at the main oscillator frequency satisfy the restriction

$$u_0=v_0,z_0=g_0/\omega _0^2.$$
(5.56)
In this case Eqs. (5.36)-(5.39) yield

$$\mu _{11}=\mu _{22},\mu _{12}=\mu _{21}=0,$$
(5.57)
so we have the unique drift matrix

$$\mu =\begin{array}{cc}\gamma & 0\\ 0& \gamma \end{array}$$
(5.58)
with the damping coefficient

$$\gamma =\pi \nu (\omega _0)\mathrm{\Delta }_0,\mathrm{\Delta }_0=u_0^2+g_0^2/\omega _0^2.$$
(5.59)
Due to Eq. (5.57) the frequency $`\omega _{}`$ in formula (5.51) exactly equals the oscillator eigenfrequency $`\omega `$<sub>0</sub>.

The equations of motion for the average values of the coordinate and momentum read

$$\dot{p}=\gamma p\omega _0^2x,$$
(5.60)
$$\dot{x}=p\gamma x.$$
(5.61)
One can see that consistent quantum mechanical consideration do not result in the conventional classical equations (2.16). To understand the origin of Eqs. (5.60) and (5.61), let us introduce the annihilation and creation operators

$$\widehat{a}=\frac{\omega _0\widehat{x}+i\widehat{p}}{\sqrt{2\mathrm{}\omega _0}},\widehat{a}^{}=\frac{\omega _0\widehat{x}i\widehat{p}}{\sqrt{2\mathrm{}\omega _0}},$$
(5.62)
which are the most natural for the description of a quantum oscillator. It turns out that precisely Eqs. (5.60), (5.61) lead to uncoupled equations for $`\widehat{a}`$ and $`\widehat{a}^{}`$:

$$\dot{\widehat{a}}=i\omega _0\widehat{a}\gamma \widehat{a},\dot{\widehat{a}}^{}=i\omega _0\widehat{a}^{}\gamma \widehat{a}^{}.$$
(5.63)
Denoting the annihilation and creation operators for the thermostat oscillators by $`\widehat{c}_i`$ and $`\widehat{c}_i^{}`$ we may rewrite the interaction Hamiltonian (5.1) (for thermostat oscillators whose frequencies coincide with $`\omega _0`$) as follows:

$$\widehat{H}_{\text{int}}(\omega _0)=\mathrm{}\underset{\omega _i=\omega _0}{}\left(\delta _0\widehat{a}\widehat{c}_i^{}+\delta _0^{}\widehat{a}^{}\widehat{c}_i\right),$$
(5.64)
$$\delta _0=g_0/\omega _0+iu_0,\mathrm{\Delta }_0=|\delta _0|^2.$$

Hamiltonian (5.64) was considered in almost every paper devoted to the models of a quantum damped oscillator (see, e.g., review ). But frequently it was chosen only because it is the simplest one. We have shown in fact that it is the only possible quadratic interaction Hamiltonian ensuring the validity of the resulting Fokker-Planck equation for any initial states taken at any initial times. This is probably related to the quantum nature of the interactions between the systems: each act of interaction must consist in the annihilation of a quantum in one system and its creation in another system. Precisely the Hamiltonian (5.64) expresses this property in the most distinct form.

Nonetheless we have no sufficient grounds for claiming that other interaction (quadratic) Hamiltonians should be excluded. But they will result in the Fokker-Planck equation with time-dependent drift and diffusion matrices describing nonexponential relaxation.

Recall that we assumed the integrals in Eqs. (5.36)-(5.42) to equal zero. What will happen if we abandon this assumption, but assume instead that the constraints (5.56) hold for all the coupling constants? In this case we have $`\kappa (\omega )0`$ at all frequencies. Furthermore, $`G(\omega )\omega ^2Z(\omega )\omega \mathrm{\Delta }(\omega )`$. Consequently, the integral terms “survive” only in the off-diagonal elements of both matrices $`\mu `$ and $`𝐃`$:

$$\mu _{11}=\mu _{22}=\gamma ,$$

$$\mu _{12}=^{}\frac{\nu (\omega )\mathrm{\Delta }(\omega )}{\omega ^2\omega _0^2}\left[\omega _0^2+\omega ^2\right]\text{d}\omega ,$$

$$\mu _{21}=^{}\frac{\nu (\omega )\mathrm{\Delta }(\omega )}{\omega ^2\omega _0^2}\text{d}\omega ,$$

$$D_{11}=\omega _0^2D_{22},D_{22}=\gamma f(\omega _0),D_{12}=\frac{1}{\pi }f(\omega )\gamma (\omega )\text{d}\omega $$
(5.65)
(the function $`\gamma (\omega )`$ is defined by Eq. (5.59) with $`\omega _0`$ replaced by $`\omega `$).

The presence of nonzero coefficients $`\mu _{12}`$ and $`\mu _{21}`$ may be interpreted as some kind of renormalization of mass and eigenfrequency of the main oscillator due to the interaction with the environment. However, inequality (2.26) is obviously violated at zero temperature, when $`f=\mathrm{}/2\omega `$, since due to Eq. (5.65) coefficient $`D_{12}`$ is strictly negative for all temperatures (evidently, both functions $`\nu (\omega )`$ and $`f(\omega )`$ are positive). This example shows once more that the self-consistent Fokker-Planck equations with time-independent coefficients can be derived from microscopic models only in exceptional cases.

## 6 Oscillator in a magnetic field. Weak coupling with a thermostat

We now consider, within the framework of the same scheme, the case where the subsystem under study is a two-dimensional isotropic oscillator with eigenfrequency $`\omega _0`$ and mass $`m`$ placed in a uniform magnetic field $``$ characterized by the cyclotron frequency

$$\omega =e/mc.$$
(6.1)
The Hamiltonian of this subsystem reads

$$\widehat{H}_0=\frac{1}{2m}\left(\widehat{\pi }_x^2+\widehat{\pi }_y^2\right)+\frac{1}{2}m\omega _0^2\left(x^2+y^2\right),$$
(6.2)
where $`\widehat{\pi }_x`$ and $`\widehat{\pi }_y`$ are the operators of the kinetic momentum projections, related to the canonical momentum $`𝐩`$ and the vector potential $`𝐀`$ in the usual way:

$$\pi =𝐩(e/c)𝐀(x,y),[\widehat{\pi }_x,\widehat{\pi }_y]=i\mathrm{}m\omega .$$
(6.3)
The problem of constructing the Fokker-Planck equation for this subsystem is reduced to that solved in the previous section, because the Hamiltonian (6.2) can be expressed as a sum of two oscillator Hamiltonians:

$$\widehat{H}_0=\mathrm{}\omega _+\left(\widehat{a}^{}\widehat{a}+\frac{1}{2}\right)+\mathrm{}\omega _{}\left(\widehat{b}^{}\widehat{b}+\frac{1}{2}\right).$$
(6.4)
The annihilation operators

$$\widehat{a}=(2m\mathrm{}\mathrm{\Omega })^{1/2}\left[\widehat{\pi }_x+i\widehat{\pi }_y+m\omega _{}\left(\widehat{y}i\widehat{x}\right)\right],$$
(6.5)
$$\widehat{b}=(2m\mathrm{}\mathrm{\Omega })^{1/2}\left[\widehat{\pi }_xi\widehat{\pi }_ym\omega _+\left(\widehat{y}+i\widehat{x}\right)\right]$$
(6.6)
satisfy the commutation relations

$$[\widehat{a},\widehat{a}^{}]=[\widehat{b},\widehat{b}^{}]=1,[\widehat{a},\widehat{b}]=[\widehat{a},\widehat{b}^{}]=0.$$
(6.7)
The frequencies are defined as follows:

$$\omega _\pm =\frac{1}{2}(\mathrm{\Omega }\pm \omega ),\mathrm{\Omega }=\left(\omega ^2+4\omega _0^2\right)^{1/2},\omega _+\omega _{}=\omega _0^2.$$
(6.8)
In the continuous weak coupling limit the self-consistent equations of motion for the first-order average values of the operators $`\widehat{a}`$ and $`\widehat{b}`$ are given by Eq. (5.63), provided one replaces $`\omega _0`$ with $`\omega _+`$ and $`\omega _{}`$. Furthermore, two different damping coefficients are possible: they are determined by the density of states and coupling constants at the frequencies $`\omega _\pm `$ (see Eqs. (5.59) and (5.64)):

$$\gamma _\pm =\pi \nu (\omega _\pm )|\delta (\omega _\pm )|^2.$$
(6.9)
The relations inverse to Eqs. (6.5) and (6.6) read

$$\widehat{\pi }_x=(m\mathrm{}/2\mathrm{\Omega })^{1/2}\left[\omega _+\left(\widehat{a}+\widehat{a}^{}\right)+\omega _{}\left(\widehat{b}+\widehat{b}^{}\right)\right],$$
(6.10)
$$\widehat{\pi }_y=i(m\mathrm{}/2\mathrm{\Omega })^{1/2}\left[\omega _+\left(\widehat{a}^{}\widehat{a}\right)+\omega _{}\left(\widehat{b}\widehat{b}^{}\right)\right],$$
(6.11)
$$\widehat{x}=i(\mathrm{}/2m\mathrm{\Omega })^{1/2}\left[\widehat{a}\widehat{a}^{}+\widehat{b}\widehat{b}^{}\right],$$
(6.12)
$$\widehat{y}=(\mathrm{}/2m\mathrm{\Omega })^{1/2}\left[\widehat{a}+\widehat{a}^{}\widehat{b}\widehat{b}^{}\right].$$
(6.13)
The average values of the coordinates and the kinetic momenta obey equations resulting from equations of the form of (5.63):

$$\dot{\pi }_x=\alpha \pi _x+\omega \pi _ym\omega _0^2x+m\omega _0^2ϵy,$$
(6.14)
$$\dot{\pi }_y=\omega \pi _x\alpha \pi _ym\omega _0^2ϵxm\omega _0^2y,$$
(6.15)
$$\dot{x}=m^1\pi _xm^1ϵ\pi _y\eta x,$$
(6.16)
$$\dot{y}=m^1ϵ\pi _x+m^1\pi _y\eta y.$$
(6.17)
We have introduced the notation

$$\alpha =\left(\gamma _+\omega _++\gamma _{}\omega _{}\right)/\mathrm{\Omega },\eta =\left(\gamma _+\omega _{}+\gamma _{}\omega _+\right)/\mathrm{\Omega },ϵ=\left(\gamma _{}\gamma _+\right)/\mathrm{\Omega }.$$
(6.18)
The second-order equations of motion read

$$\ddot{x}+\left(\gamma _{}+\gamma _+\right)\dot{x}\omega \dot{y}+\left(\omega _0^2+\gamma _{}\gamma _+\right)x\left(\gamma _{}\omega _+\gamma _+\omega _{}\right)y=0,$$
(6.19)
$$\ddot{y}+\left(\gamma _{}+\gamma _+\right)\dot{y}+\omega \dot{x}+\left(\omega _0^2+\gamma _{}\gamma _+\right)y+\left(\gamma _{}\omega _+\gamma _+\omega _{}\right)x=0.$$
(6.20)
We see that “one-photon” interaction with a thermostat of the form of (5.64) results in coordinate-dependent forces perpendicular to the vector $`𝐫=(x,y)`$ and proportional to the damping coefficients. The necessity of introducing such forces was shown earlier in  within the framework of a phenomenological approach. As was demonstrated in these papers, if the second-order equations of motion (for a charged particle or oscillator placed in a uniform magnetic field) contain the term $`\gamma \dot{𝐫}`$, then it is impossible to satisfy simultaneously Eq. (2.19) with the equilibrium matrix $`(T)`$ for all temperatures (including $`T=0`$) and condition (2.20), unless a force of the form $`𝐟=[𝐡\times 𝐫]`$ is introduced. Now we have arrived at the same result on the base of the microscopic model. Moreover, the relation between the velocity and the kinetic momentum becomes much more complicated than in the conservative case (see Eqs. (6.16) and (6.17).

It is clear from the preceding section that the steady state Wigner function at $`t\mathrm{}`$ coincides with the equilibrium distribution that was found in refs. . Since this distribution is Gaussian, it is completely determined (see Eq. (3.8)) by the equilibrium variance matrix

$$^{(\text{eq})}=\begin{array}{cccc}_\pi & 0& 0& _a\\ 0& _\pi & _a& 0\\ 0& _a& _\rho & 0\\ _a& 0& 0& _\rho \end{array},$$
(6.21)
$$_\pi =\frac{m\mathrm{}\mathrm{\Omega }}{4Q}\left[\left(1+\frac{\omega ^2}{\mathrm{\Omega }^2}\right)\mathrm{sinh}\stackrel{~}{\mathrm{\Omega }}2\frac{\omega }{\mathrm{\Omega }}\mathrm{sinh}\stackrel{~}{\omega }\right],$$
(6.22)
$$_\rho =\frac{\mathrm{}\mathrm{sinh}\stackrel{~}{\mathrm{\Omega }}}{m\mathrm{\Omega }Q},$$
(6.23)
$$_a=\frac{\mathrm{}}{2Q}\left[\frac{\omega }{\mathrm{\Omega }}\mathrm{sinh}\stackrel{~}{\mathrm{\Omega }}\mathrm{sinh}\stackrel{~}{\omega }\right],$$
(6.24)
where

$$Q(\beta )=\mathrm{cosh}\stackrel{~}{\mathrm{\Omega }}\mathrm{cosh}\stackrel{~}{\omega },\stackrel{~}{\mathrm{\Omega }}=\frac{1}{2}\beta \mathrm{}\mathrm{\Omega },\stackrel{~}{\omega }=\frac{1}{2}\beta \mathrm{}\omega ,$$
(6.25)
and $`\beta =1/kT`$ is the inverse temperature of the thermostat.

The drift matrix $`𝐀`$ corresponding to Eqs. (6.14)-(6.17) reads

$$𝐀=\begin{array}{cccc}\alpha & \omega & m\omega _0^2& m\omega _0^2ϵ\\ \omega & \alpha & m\omega _0^2ϵ& m\omega _0^2\\ m^1& m^1ϵ& \eta & 0\\ m^1ϵ& m^1& 0& \eta \end{array}.$$
(6.26)
Putting this matrix into Eq. (5.47) with the matrix $`^{(\text{eq})}`$ instead of $`𝐅_0`$ we obtain the diffusion matrix

$$𝐃=\begin{array}{cccc}D_\pi & 0& 0& D_a\\ 0& D_\pi & D_a& 0\\ 0& D_a& D_\rho & 0\\ D_a& 0& 0& D_\rho \end{array}$$
(6.27)
with the following coefficients,

$$D_\pi =\frac{m\mathrm{}}{2\mathrm{\Omega }Q}\left[\left(\gamma _+\omega _+^2+\gamma _{}\omega _{}^2\right)\mathrm{sinh}\stackrel{~}{\mathrm{\Omega }}\left(\gamma _+\omega _+^2\gamma _{}\omega _{}^2\right)\mathrm{sinh}\stackrel{~}{\omega }\right],$$
(6.28)
$$D_a=\frac{\mathrm{}}{2\mathrm{\Omega }Q}\left[\left(\gamma _+\omega _+\gamma _{}\omega _{}\right)\mathrm{sinh}\stackrel{~}{\mathrm{\Omega }}\left(\gamma _+\omega _++\gamma _{}\omega _{}\right)\mathrm{sinh}\stackrel{~}{\omega }\right],$$
(6.29)
$$D_\rho =\frac{\mathrm{}}{2m\mathrm{\Omega }Q}\left[\left(\gamma _++\gamma _{}\right)\mathrm{sinh}\stackrel{~}{\mathrm{\Omega }}\left(\gamma _+\gamma _{}\right)\mathrm{sinh}\stackrel{~}{\omega }\right].$$
(6.30)
In particular, at zero temperature $`(\beta =\mathrm{})`$ we obtain

$`D_\pi ^{(\text{low})}`$ $`=`$ $`{\displaystyle \frac{m\mathrm{}}{2\mathrm{\Omega }}}\left(\gamma _+\omega _+^2+\gamma _{}\omega _{}^2\right),`$
$`D_a^{(\text{low})}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2\mathrm{\Omega }}}\left(\gamma _+\omega _+\gamma _{}\omega _{}\right),`$
$`D_\rho ^{(\text{low})}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2m\mathrm{\Omega }}}\left(\gamma _++\gamma _{}\right).`$ (6.31)
In the opposite, high-temperature, case $`(\beta 0)`$ we have

$$D_\pi ^{(\text{high})}=mkT\alpha ,D_\rho ^{(\text{high})}=\frac{kT\eta }{m\omega _0^2},D_a^{(\text{high})}=kTϵ.$$
(6.32)
It is noteworthy that all three diffusion coefficients remain nonzero even in the high temperature limit, which is usually identified with the quasiclassical limit. Recall that in classical statistical mechanics it is usually implied that the only nonzero diffusion coefficient is $`D_\pi `$.

Various sets of the diffusion coefficients compatible with inequality (2.20) and leading to an equilibrium steady state (with the variance matrix (6.21)) in the limit of infinitely small damping were constructed within the framework of the phenomenological approach in . However, none of them had a structure similar to that given by Eqs. (6.28)-(6.30).

For example, only the coefficient $`D_\pi `$ was proportional to the temperature in the high-temperature limit, whereas the other diffusion coefficients decreased as $`1/kT`$, in contrast to Eq. (6.32). This difference is due to at least two causes. First, it was assumed in  that the elements $`A_{41}`$ and $`A_{32}`$ of the drift matrix must be zero, although the elements $`A_{14}`$ and $`A_{23}`$ could be nonzero. Eq. (6.26) shows that within the framework of the microscopic approach such a choice is impossible, since all these coefficients are proportional to the parameter $`ϵ`$. Furthermore, in the aforementioned papers we admitted the possibility that some coefficients of the drift matrix (those related to the damping) could depend on temperature. In principle, such a possibility (i.e., the time dependence of the coupling constants in the interaction Hamiltonian) is not excluded within the framework of the microscopic approach as well. Then the high temperature limit of the diffusion matrix coefficients could be quite different from Eq. (6.32).

The expressions (6.26)-(6.30) are simplified in the special case of the equal damping coefficients, when $`\gamma _+=\gamma _{}=\gamma _0`$:

$$𝐀^{(0)}=\begin{array}{cccc}\gamma _0& \omega & m\omega _0^2& 0\\ \omega & \gamma _0& 0& m\omega _0^2\\ m^1& 0& \gamma _0& 0\\ 0& m^1& 0& \gamma _0\end{array},$$
(6.33)
$$D_\pi ^{(0)}=\frac{m\mathrm{}\gamma _0}{2\mathrm{\Omega }Q}\left[\left(\omega ^2+2\omega _0^2\right)\mathrm{sinh}\stackrel{~}{\mathrm{\Omega }}\omega \mathrm{\Omega }\mathrm{sinh}\stackrel{~}{\omega }\right],$$
(6.34)
$$D_a^{(0)}=\frac{\mathrm{}\gamma _0}{2\mathrm{\Omega }Q}\left[\omega \mathrm{sinh}\stackrel{~}{\mathrm{\Omega }}\mathrm{\Omega }\mathrm{sinh}\stackrel{~}{\omega }\right],$$
(6.35)
$$D_\rho ^{(0)}=\frac{\mathrm{}\gamma _0}{m\mathrm{\Omega }Q}\mathrm{sinh}\stackrel{~}{\mathrm{\Omega }}.$$
(6.36)
If coefficient $`\gamma _{}`$ tends to zero sufficiently rapidly as $`\omega _{}0`$, then the set of diffusion coefficients (6.28)-(6.30) possesses the finite limit for a free particle in a magnetic field, when $`\omega _0=\omega _{}=0`$, $`\mathrm{\Omega }=\omega `$, $`\gamma _+(\omega _+\omega )=\gamma `$:

$$D_\pi =\frac{1}{2}\gamma m\mathrm{}\omega \mathrm{coth}\stackrel{~}{\omega },$$
(6.37)
$$D_a=\frac{1}{2}\gamma \mathrm{}\mathrm{coth}\stackrel{~}{\omega },$$
(6.38)
$$D_\rho =\frac{\gamma \mathrm{}}{2m\omega }\mathrm{coth}\stackrel{~}{\omega },$$
(6.39)
$$\alpha =\gamma ,\eta =0,ϵ=\gamma /\omega .$$
(6.40)
In this case the operators $`\widehat{b}`$ and $`\widehat{b}^{}`$ become the integrals of motion, whose real and imaginary parts are connected with the center-of-orbit operators in a uniform magnetic field -.

One should remember, however, that the results obtained in this section can be justified only under rather strong limitations imposed on the interaction Hamiltonian. First, it must be written in the specific form (5.64) at the resonant frequencies. Secondly, the off-resonance terms must ensure the disappearance of the integral terms in Eqs. (5.36)-(5.42). In particular, the density of states must decrease sufficiently rapidly as $`\omega \mathrm{}`$.

All these conditions are violated, for example, in the case where the role of a reservoir is played by a quantized electromagnetic field coupled to the oscillator by means of the standard interaction Hamiltonian in the dipole approximation:

$$\widehat{H}_{\text{int}}=\frac{e}{m}(2\pi \mathrm{})^{1/2}\underset{j,\sigma }{}\frac{\left(\widehat{\pi }\tau _{j,\sigma }\right)}{\sqrt{\omega _j}}\left[\widehat{c}(𝐤_j,\sigma )+\widehat{c}^{}(𝐤_j,\sigma )\right],$$
(6.41)
where $`\widehat{c}(𝐤_j,\sigma )`$ is the operator of annihilation of a photon with wave vector $`𝐤_j`$, frequency $`\omega _j=c|𝐤_j|`$, and polarization $`\sigma `$; $`\tau _{j,\sigma }`$ is the unit polarization vector perpendicular to the vector $`𝐤_j`$. The density of states is proportional to $`\omega ^2`$ in this case, and the integrals in Eqs. (5.36)-(5.42) diverge. Consequently, in this case the radiation damping leads to nonexponential relaxation.