# Revisiting the damped quantum harmonic oscillator

Stephen M. Barnett  School of Physics and Astronomy, University of Glasgow,
Glasgow G12 8QQ, UK<EMAIL_ADDRESS>James D. Cresser  School of
Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, UK  Department
of Physics and Astronomy, University of Exeter EX4 4QL, UK  Department of
Physics and Astronomy, Faculty of Science and Engineering, Macquarie
University, NSW 2109, Australia<EMAIL_ADDRESS>Sarah Croke
School of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, UK
<EMAIL_ADDRESS>
###### Abstract

We reanalyse the quantum damped harmonic oscillator, introducing three less
than common features. These are (i) the use of a continuum model of the
reservoir rather than an ensemble of discrete oscillators, (ii) an exact
diagonalisation of the Hamiltonian by adapting a technique pioneered by Fano,
and (iii) the use of the thermofield technique for describing a finite
temperature reservoir. We recover in this way a number of well-known and some,
perhaps, less familiar results. An example of the latter is an _ab initio_
proof that the oscillator relaxes to the mean-force Gibbs state. We find that
special care is necessary when comparing the damped oscillator with its
undamped counterpart as the former has two distinct natural frequencies, one
associated with short time evolution and the other with longer times.

###### pacs:

03.65.Yz, 42.50.Lc

## Preface: an apology

There can surely be no more intensively studied open quantum system than the
damped harmonic oscillator. This makes it all but impossible to do justice to
the vast literature on the subject, and authors seeking to write on this
system need to tread carefully and to acknowledge, freely, that much of the
relevant literature will be, unintentionally but inevitably, overlooked. The
task is made yet more hazardous by the fact that different communities in
physics have covered similar ground more or less independently. Nevertheless,
this special issue offers the opportunity to reexamine this well-studied
system, with the aim of demonstrating some techniques that should, perhaps, be
better known.

##  1 Introduction

Recent technological advances make it possible to realise simple mechanical
devices in the microscopic and nanoscopic regimes, the properties of which are
determined by quantum effects [1]. The existence of these represents a
remarkable opportunity for fundamental studies of light-matter interactions
[2] and also the potential for practical application to quantum communications
and information processing [3]. Yet they present also a challenge to existing
methods of analysis, many of which were developed to treat more rapidly
oscillating systems with weaker couplings. The strong coupling regime brings
with it some surprises, such as the possibility that quantum entanglement
might persist in the high-temperature limit [4].

The quantum theory of machines is built, to a large extent on the theory of
oscillators and strongly coupled oscillators, which are coupled to one or more
environments, each of which is at a characteristic temperature [5, 6, 7]. The
behaviour of these quantum systems is governed not simply by average
properties but also fluctuations, and has been informed by the development of
fluctuation theorems for quantum open systems [8, 9]. Such considerations
underpin developments in the rapidly advancing field of quantum thermodynamics
[10, 11].

The coupling to the environment needs to be handled with some care because of
the possibility of quantum coherence and the development of entanglement
between the system of interest and its environment, with the result that
apparently unphysical behaviours may emerge [12]. The requirement to treat the
coupling to the reservoirs with care provides the incentive to return to the
problem of a single strongly-damped oscillator and to treat this model
exactly, using the techniques identified in the abstract.

There exist at present very many mathematically and physically acceptable
methods for treating damped harmonic oscillators. Common to most, if not all,
of these developments is the treatment of the surrounding environment, or
reservoir, as an ensemble of harmonic oscillators with a broad spectrum of
oscillator frequencies. It is the dephasing brought about by this spread of
frequencies that introduces the damping, or irreversibility, in the oscillator
dynamics. The environmental harmonic oscillators may be physical oscillators
or vibration modes, as in the Caldeira-Leggett model [13, 14], or the modes of
the quantised electromagnetic field, as in many quantum optics applications
[15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29].

For weakly damped oscillators, such as those encountered regularly in quantum
optics, there are master equations and the corresponding Heisenberg-picture
operator Langevin equations. Even in this weakly damped regime, the
dissipative dynamics can be challenging with a rich structure of asymptotic
states [30, 31]. For more strongly damped systems memory effects become
important and there is a departure from Markovian evolution [32, 33, 34, 35,
36, 37, 38, 39]. Yet stronger coupling requires the inclusion of counter-
rotating interactions, which do not conserve the number of quanta. Among the
many and varied approaches adopted are the aforementioned Schrödinger-picture
master equations for the oscillator density operator [40, 41] and Heisenberg-
Langevin operator equations for oscillator observables driven by environmental
fluctuations [42, 43, 44]. Also widely used are Feynman-Vernon path integral
and related techniques [14, 45, 46, 47, 48, 49]. For finite temperatures the
environmental oscillators are considered to be prepared initially in thermal
states with Bose-Einstein statistics appropriate to the reservoir temperature.
This leads, in the Schrödinger picture, to the reservoir acting both as a
source as well as a sink of quanta. The techniques for treating this include a
product of thermal density operators for the reservoir modes [16, 18, 20, 21,
22], imaginary time methods and thermal Green functions [50, 51, 52] and also
thermofield dynamics [53, 54, 55, 56, 57, 58]. In addition to these methods,
there has also been work done on diagonalising the oscillator-reservoir
Hamiltonian in which the reservoir is formed from a collection of harmonic
oscillators [60, 61]. Our analysis develops and expands upon material in an
earlier preprint [62] (see also [63, 64, 65, 66]), which treats the
environment as a continuum. It is complementary to that adopted by Philbin who
has tackled this problem of the oscillator evolution using Green functions
[67]. We hope that the combination of his work and ours will provide a more
complete understanding.

##  2 Background

The harmonic oscillator has a special place in physics as one of the simplest
and most widely employed of physical models. The reasons for its ubiquity, no
doubt, are its simplicity and the fact that it is readily analysed. In the
quantum domain, the harmonic oscillator is barely more difficult to treat than
its classical counterpart and was one of the first dynamical systems to which
Schrödinger applied his equation [68]. Today, both the classical and quantum
forms appear in elementary courses on classical and quantum mechanics.

The damped harmonic oscillator loses energy as a result of coupling to the
surrounding environment. In the classical domain it often suffices to describe
this in terms of a simple damping coefficient, $\gamma$, and an associated
stochastic or Langevin force [69], $F(t)$, which models the effect of
environmental fluctuations on the oscillator [70, 71, 72, 73, 74]. The
dynamics is described by a simple linear differential equation of the form

$\ddot{x}+\gamma\dot{x}+\omega^{2}_{0}x=\frac{F(t)}{m},$ (1)

where $m$ is the mass of the oscillating particle. There is no requirement for
detailed knowledge of the fluctuating force, which may be considered to have a
very short correlation time with a magnitude determined by the requirements of
thermodynamic equilibrium.

The damped quantum harmonic oscillator requires that explicit account be taken
of the quantum nature of the environmental degrees of freedom [15], which are
most simply described by an ensemble of harmonic oscillators [16]. If the
damping is very weak, so that $\gamma\ll\omega_{0}$, then we can neglect
rapidly oscillating terms in the coupling between the oscillator and the
environmental oscillators by making the rotating wave approximation, which
corresponds to enforcing the conservation of the total number of vibrational
quanta, and then the Born and Markov approximations associated with weak
coupling and loss of memory in the reservoir [17]. This leads to the master
equations and Heisenberg-Langevin equations that are ubiquitous, most
especially, in quantum optics [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
29].

If the coupling is somewhat stronger then it may not be possible to make the
rotating wave approximation and we then need to retain in the Hamiltonian
terms that can create or annihilate a pair of quanta, one in the damped
oscillator and, at the same time, one in the environment. This leads to the
Caldeira-Leggett model [67, 75, 76, 77, 13, 14, 60], which we describe in the
following section, and which has been applied to study a wide variety of
quantum open systems [25, 29]. A variant on the model has been applied to the
quantum theory of light in dielectric [63, 64] and magneto-dielectric media
[78, 79, 80, 81]. An important complication that seems to be an inevitable
consequence working in this strong-coupling regime is the failure of the
Markov approximation; attempts to enforce this approximation lead to a master
equation that is unphysical in that there exist initial states for which the
dynamics leads to negative probabilities [25, 82, 83, 84, 85]. It is possible
to derive a master equation but the resulting equation is one that has within
it non-trivial time-dependent coefficients [60, 40]. This time-dependence is a
clear signature of the non-Markovian nature of the associated evolution. It
seems that this non-Markovian character is an inevitable feature of the
strongly-damped quantum harmonic oscillator.

##  3 Hamiltonian for the strongly-damped harmonic oscillator

Figure 1: Representation of a single harmonic oscillator (in red) coupled
harmonically to a bath of oscillators (in blue).

Consider a harmonic oscillator that is strongly coupled to its environment,
ultimately to be treated as a thermodynamic reservoir, modelled as large
collection of oscillators, with a range of frequencies, through their
respective positions as depicted in figure 1. We write the Hamiltonian for the
combined oscillator-reservoir system in the form [13, 60]

$\hat{H}=\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\Omega^{2}_{0}\hat{x}^{2}+\sum_{\mu}\left(\frac{\hat{p}^{2}_{\mu}}{2m_{\mu}}+\frac{1}{2}m_{\mu}\omega^{2}_{\mu}\hat{x}^{2}_{\mu}\right)-\sum_{\mu}m_{\mu}\omega_{\mu}^{2}\lambda_{\mu}\hat{x}_{\mu}\hat{x}.$
(2)

If we complete the square we can rewrite this in a minimal-coupling form to
arrive at the alternative form:

$\hat{H}=\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\omega^{2}_{0}\hat{x}^{2}+\sum_{\mu}\left[\frac{\hat{p}^{2}_{\mu}}{2m_{\mu}}+\frac{1}{2}m_{\mu}\omega^{2}_{\mu}\left(\hat{x}_{\mu}-\lambda_{\mu}\hat{x}\right)^{2}\right],$
(3)

where

$\omega_{0}^{2}=\Omega_{0}^{2}-\sum_{\mu}\frac{m_{\mu}}{m}\omega^{2}_{\mu}\lambda^{2}_{\mu}.$
(4)

Each term in our Hamiltonian (3) is strictly positive only if this quantity is
positive, corresponding to a real frequency $\omega_{0}$. If it is negative
then the second term is also negative and the Hamiltonian is not bounded from
below and hence not allowed physically. Hence the positivity of the
Hamiltonian places a physical restriction on the strength of the coupling

$\Omega_{0}^{2}>\sum_{\mu}\frac{m_{\mu}}{m}\omega^{2}_{\mu}\lambda^{2}_{\mu},$
(5)

while $\omega_{0}$ can take any positive value. 111In the literature the
change from the Hamiltonian (2) to (3) is often described as adding a counter
term [25, 42]. We find no need to speak of adding counter terms as both of the
natural frequencies $\omega_{0}$ and $\Omega_{0}$ play a role in the dynamics.

At this point it is necessary to pause and consider the fact that our damped
harmonic oscillator seems to have _two_ possible natural frequencies,
$\Omega_{0}$ and $\omega_{0}$. These appear as the potential energy in our two
forms of the Hamiltonian given in Eqs. (2) and (3). It is reasonable,
therefore, to ask which of these (if either) corresponds to the ‘true’ natural
frequency. In order to address this important point, we step aside from our
principal objective of diagonalizing the Hamiltonian to derive the Heisenberg-
Langevin equation for the oscillator position operator.

###  3.1 Heisenberg-Langevin equation: a tale of two frequencies

We find that _both_ the frequencies $\omega_{0}$ and $\Omega_{0}$ have roles
to play in the dynamics of the oscillator and that, in this sense, both fulfil
the roles of natural frequency of the damped oscillator, albeit in different
time domains. To demonstrate this we can derive form our Hamiltonian a
Heisenberg-Langevin equation of motion for the position operator of our
oscillator. The details are given in A. We find

$\ddot{\hat{x}}(t)+\int_{0}^{t}\kappa(t-t^{\prime})\dot{\hat{x}}(t^{\prime})dt^{\prime}+\left(\Omega_{0}^{2}-\kappa(0)\right)\hat{x}(t)+\kappa(t)\hat{x}(0)=\frac{\hat{F}(t)}{m}\,,$
(6)

where $\kappa(t)$ is the memory kernel:

$\kappa(t)=\sum_{\mu}\frac{m_{\mu}}{m}\omega_{\mu}^{2}\lambda_{\mu}^{2}\cos(\omega_{\mu}t)$
(7)

and $F(t)$ is the Langevin force:

$\hat{F}(t)=\omega^{2}_{\mu}\lambda_{\mu}\left(\hat{x}_{\mu}(0)\cos(\omega_{\mu}t)+\frac{\hat{p}_{\mu}(0)}{m_{\mu}\omega_{\mu}}\sin(\omega_{\mu}t)\right)\,.$
(8)

Note that this equation of motion is (essentially) exact. In order to
interpret the various terms it suffices to consider the expectation value of
this operator equation:

$\langle\ddot{\hat{x}}(t)\rangle+\int_{0}^{t}\kappa(t-t^{\prime})\langle\dot{\hat{x}}(t^{\prime})\rangle
dt^{\prime}+\left(\Omega_{0}^{2}-\kappa(0)\right)\langle\hat{x}(t)\rangle+\kappa(t)\langle\hat{x}(0)\rangle=0\,,$
(9)

and note that we can also write this in the form

$\langle\ddot{\hat{x}}(t)\rangle+\int_{0}^{t}\kappa(t-t^{\prime})\langle\dot{\hat{x}}(t^{\prime})\rangle
dt^{\prime}+\omega_{0}^{2}\langle\hat{x}(t)\rangle+\kappa(t)\langle\hat{x}(0)\rangle=0\,.$
(10)

The first of these equations is written in terms of the frequency $\Omega_{0}$
and the second in terms of $\omega_{0}$. It remains to determine the physical
role of each of these, which we can do by considering very short and longer
time scales.

####  3.1.1 Ultra-short time scales

To better appreciate what happens at very short times, we first undo the
integration by parts that led to our equation of motion and write this in the
form

$\langle\ddot{\hat{x}}(t)\rangle+\int_{0}^{t}\dot{\kappa}(t-t^{\prime})\langle\hat{x}(t^{\prime})\rangle
dt^{\prime}+\Omega_{0}^{2}\langle\hat{x}(t)\rangle=0\,.$ (11)

For a very short time, $\delta t$, this becomes

$\langle\ddot{\hat{x}}(\delta t)\rangle+\Omega_{0}^{2}\langle\hat{x}(\delta
t)\rangle+O(\delta t^{2})=0\,.$ (12)

The integral term in Eq. (11) is of order $\delta t^{2}$ by virtue of Eq. (7)
as for short times $\dot{\kappa}$ is of order $\delta t$. The combination of
Eq. (12) and $\dot{\hat{x}}=\hat{p}/m$, which is always true, leads to the
ultra-short time behaviour

$\displaystyle\langle\hat{x}(\delta t)\rangle$ $\displaystyle=$
$\displaystyle\langle\hat{x}(0)\rangle+\frac{\langle\hat{p}(0)}{m}\delta t$
$\displaystyle\langle\hat{p}(\delta t)\rangle$ $\displaystyle=$
$\displaystyle\langle\hat{p}(0)\rangle-m\Omega_{0}^{2}\langle\hat{x}(0)\rangle\delta
t\,.$ (13)

This is the short-time behaviour of a harmonic oscillator of frequency
$\Omega_{0}$. It is clear that this is the frequency of the oscillator in the
non-Markovian regime.

####  3.1.2 Longer time scales

Our first task is to firm up what we mean by longer time. To do so we note
that the function $\kappa(t)$ involves a summation of oscillating cosines of
different frequencies, one for each environment oscillator, and that these
will dephase, causing $\kappa(t)$ to decay. We define the longer-time regime
as that for which we can approximate $\kappa(t)\approx 0$. In this regime, our
our equation of motion (10) becomes

$\langle\ddot{x}(t)\rangle+\int_{0}^{t}\kappa(t-t^{\prime})\langle\dot{x}(t^{\prime})\rangle
dt^{\prime}+\omega^{2}_{0}\langle x(t)\rangle=0\,.$ (14)

This equation retains the possibility of non-Markovian effects in the (second)
damping term but it is clear that the natural frequency of the oscillator in
this time regime is $\omega_{0}$ and _not_ $\Omega_{0}$. Thus the question of
whether our Hamiltonian applies in the over-damped regime depends on the
relationship between the damping and $\omega_{0}$ and not $\Omega_{0}$.

We are now in a position to address the issue of whether our model Hamiltonian
can be applied in the over-damped regime. To enter this regime, we need to
reduce the natural frequency of the oscillator so that it is below the damping
rate. It is clear from the inequality (5) that there is always a lower bound
below which $\Omega_{0}$ cannot be reduced. _However_ there is no such bound
on $\omega_{0}$ and, indeed, we can set $\omega_{0}=0^{+}$ without
invalidating our model. As it is $\omega_{0}$ and not $\Omega_{0}$ that is the
natural frequency of the oscillator beyond the ultra-short time regime, it is
clear that we can use our model to describe an oscillator in the strongly
over-damped regime.

A simple example might help to illustrate the ideas presented above. We
consider the simplest case of Ohmic damping, or Ohmic friction [29]. To this
end, we go to the continuum limit in which

$\kappa(t)\rightarrow\frac{2}{\pi
m}\int_{0}^{\infty}\frac{J(\omega)}{\omega}\cos(\omega t)d\omega\,,$ (15)

with

$J(\omega)=m\gamma\omega e^{-\omega/\omega_{c}}\,.$ (16)

Note that it is essential to include a frequency cut-off in this case as,
without this, we cannot satisfy the inequality (5). We find that

$\kappa(t)=\frac{2}{\pi}\frac{\gamma\omega_{c}}{1+\omega_{c}^{2}t^{2}}\quad\Rightarrow\quad\kappa(0)=\frac{2\gamma\omega_{c}}{\pi}\,.$
(17)

This also demonstrates the necessity of a cut-off frequency.

Recall that the ultra-short time regime corresponds to times for which we
retain the term $\kappa(t)\langle x(0)\rangle$ in (9). This means times for
which $t<\omega_{c}^{-1}$, which fits with the familiar idea that the time
needs to be short compared inverse bandwidth of the reservoir. As there needs
to be, in the exact theory, a short-time non-Markovian regime, we need the
cut-off frequency in our bath coupling.

For longer times, for which $\kappa(t)\approx 0$, our equation for the
expectation value of the position reduces to (14). If the coupling to the
reservoir is sufficiently weak that the expectation value
$\langle\dot{x}\rangle$ does not vary significantly on the timescale
$\omega_{0}^{-1}$, then we can make the approximation

$\displaystyle\int_{0}^{t}\kappa(t-t^{\prime})\langle\dot{x}(t^{\prime})\rangle
dt^{\prime}$ $\displaystyle\approx$
$\displaystyle\langle\dot{x}(t)\rangle\int_{0}^{t}\kappa(t-t^{\prime})dt^{\prime}$
(18) $\displaystyle=$ $\displaystyle\gamma\langle\dot{x}(t)\rangle.$

It then follows that our equation for the expectation value of the position
becomes (in this Markovian, longer-time regime)

$\langle\ddot{x}(t)\rangle+\gamma\langle\dot{x}(t)\rangle+\omega^{2}_{0}\langle
x(t)\rangle=0\,,$ (19)

which is the familiar equation for a damped harmonic oscillator. Note however,
that it is the frequency $\omega_{0}$ and _not_ $\Omega_{0}$ that appears. We
enter the over-damped regime when $\omega_{0}<\gamma/2$. The constraints on
the natural frequencies in (4) correspond in this Ohmic damping example to

$\omega_{0}^{2}>0\quad{\rm
and}\quad\Omega_{0}^{2}>\kappa(0)=\frac{2\gamma\omega_{c}}{\pi}\,.$ (20)

We see that there is no lower bound on $\omega_{0}$ (although it must be real
and greater than zero) but there is a bound on $\Omega_{0}$: it must exceed
the geometric mean of the damping rate $\gamma$ and the cut-off frequency
$\omega_{c}$ multiplied by $2/\pi$. Thus we see again that we cannot allow the
cut-off frequency to tend to infinity. Similarly the frequency $\Omega_{0}$
will lie somewhere between $\gamma$ and the rather larger $\omega_{c}$ so we
cannot have $\Omega_{0}$ less than $\gamma$. What saves the model in the
strongly damped regime is the fact that it is $\omega_{0}$ rather than
$\Omega_{0}$ that corresponds to the natural frequency of the damped
oscillator.

As a final note in this section, we note that the presence of at least two
candidate natural frequencies is all but inevitable for a strongly damped
oscillator, as such system will, in general, always experience a non-Markovian
short-time evolution. As we shall see, the existence of these two frequencies
also complicates the question of the amount of energy associated with the
oscillator during its evolution, but especially in its steady state.

###  3.2 Continuum reservoir

Our first, perhaps, less familiar feature is to replace the discrete reservoir
of oscillators by a continuum. To proceed, we first rewrite our Hamiltonian in
terms of the familiar annihilation and creation operators:

$\displaystyle\hat{a}=\sqrt{\frac{m\Omega_{0}}{2\hbar}}\left(\hat{x}+\frac{i\hat{p}}{m\Omega_{0}}\right)$
$\displaystyle\hat{b}_{\mu}=\sqrt{\frac{m_{\mu}\omega_{\mu}}{2\hbar}}\left(\hat{x}_{\mu}+\frac{i\hat{p}_{\mu}}{m_{\mu}\omega_{\mu}}\right)$
(21)

In terms of these operators our Hamiltonian, Eq. (2), becomes

$\hat{H}=\hbar\Omega_{0}\hat{a}^{\dagger}\hat{a}+\sum_{\mu}\hbar\omega_{\mu}\hat{b}_{\mu}^{\dagger}\hat{b}_{\mu}+\sum_{\mu}\frac{\hbar}{2}V_{\mu}\left(\hat{a}+\hat{a}^{\dagger}\right)\left(\hat{b}_{\mu}+\hat{b}^{\dagger}_{\mu}\right)$
(22)

when unimportant constant shifts in the ground-state energies are removed and

$V_{\mu}=-\sqrt{\frac{m_{\mu}\omega_{\mu}}{m\Omega_{0}}}\omega_{\mu}\lambda_{\mu}.$
(23)

When written in terms of this quantity, our positivity condition (5) becomes

$\Omega_{0}>\sum_{\mu}\frac{V^{2}_{\mu}}{\omega_{\mu}}.$ (24)

At this stage we can seek to diagonalise the full Hamiltonian to find, in
effect, the normal modes of the oscillator coupled to the reservoir. This is
the approach taken by Haake and Reibold and by Ford _et al_ [60, 61]. The
dynamics are then reminiscent of the Bixon-Jortner model, with recurrences
occurring on a timescale given by the inverse of the frequency spacing of the
reservoir oscillators [17, 86, 87, 88, 89, 90], as is characteristic of
periodic and almost periodic functions [91] 222As is often the case in
science, the Bixon-Jortner model itself had an anticipation in the early work
of Fano [92]. We are grateful to Jan Petter Hansen for bringing Fano’s paper
to our attention..

We find it both simpler and also more powerful to first recast our model in
terms of a continuum description of the reservoir. To this end we introduce
continuum annihilation and creation operators, $\hat{b}(\omega)$ and
$\hat{b}^{\dagger}(\omega)$, satisfying the commutation relations

$\left[\hat{b}(\omega),\hat{b}^{\dagger}(\omega^{\prime})\right]=\delta(\omega-\omega^{\prime})$
(25)

and our Hamiltonian becomes [64]

$\displaystyle\hat{H}=\hbar\Omega_{0}\hat{a}^{\dagger}\hat{a}+\int_{0}^{\infty}d\omega\>\hbar\omega\hat{b}^{\dagger}(\omega)\hat{b}(\omega)$
$\displaystyle\qquad\qquad+\int_{0}^{\infty}d\omega\>\frac{\hbar}{2}V(\omega)\left(\hat{a}+\hat{a}^{\dagger}\right)\left[\hat{b}^{\dagger}(\omega)+\hat{b}(\omega)\right]$
(26)

and the positivity condition is then

$\Omega_{0}>\int_{0}^{\infty}d\omega\frac{V^{2}(\omega)}{\omega}\,.$ (27)

Our Hamiltonian is quadratic in the annihilation and creation operators for
the oscillator and the reservoir and hence leads to linear coupled equations
of motion for these operators. We could seek to solve these equations of
motion and this would lead to an operator Heisenberg-Langevin equation similar
to that derived above and in a number of earlier texts [42, 43, 44]. Here we
adopt the different approach of diagonalising the Hamiltonian.

###  3.3 Hamiltonian diagonalisation

Our second less familiar element is the exact diagonalisation of the
oscillator-reservoir Hamiltonian. We shall find that this can be achieved with
greater generality than is possible for the model with discrete reservoir
oscillators simply because there is greater freedom in the evaluation of
integrals than summations. Our task is to diagonalise the Hamiltonian by
finding a complete set of eigenoperators, $\hat{B}(\omega)$ and their
conjugates $\hat{B}^{\dagger}(\omega)$ that satisfy the operator equations

$\displaystyle\left[\hat{B}(\omega),\hat{H}\right]=\hbar\omega\hat{B}(\omega)$
$\displaystyle\left[\hat{B}^{\dagger}(\omega),\hat{H}\right]=-\hbar\omega\hat{B}^{\dagger}(\omega)$
(28)

for all positive frequencies $\omega$. These eigenoperators are complete if,
in addition to these they also satisfy the condition

$\left[\hat{B}(\omega),\hat{B}^{\dagger}(\omega^{\prime})\right]=\delta(\omega-\omega^{\prime}).$
(29)

These operator equations are the natural analogues of the more familiar
eigenvalue and completeness conditions for the eigenstates of a Hamiltonian
[17, 93]. In analogy with the eigenvalue problem, we expand each of the
eigenoperators as a superposition of a complete set of operators:

$\hat{B}(\omega)=\alpha(\omega)\hat{a}+\beta(\omega)\hat{a}^{\dagger}+\int_{0}^{\infty}d\omega^{\prime}\>\left[\gamma(\omega,\omega^{\prime})\hat{b}(\omega^{\prime})+\delta(\omega,\omega^{\prime})\hat{b}^{\dagger}(\omega^{\prime})\right],$
(30)

and then use the eigenoperator equations and completeness condition to
determine the coefficients in this expansion. The calculation is a little
involved but the main points are summarised in B.

We can express any of the annihilation and creation operators for the
oscillator or the reservoir in terms of our eigenoperators. To do this we
write the desired operator as a superposition of all the $\hat{B}(\omega)$ and
$\hat{B}^{\dagger}(\omega)$ operators and then use the commutation relations
to extract the coefficients in this expansion. For the oscillator operators we
find

$\displaystyle\hat{a}=\int_{0}^{\infty}d\omega\left(\alpha^{*}(\omega)\hat{B}(\omega)-\beta(\omega)\hat{B}^{\dagger}(\omega)\right)$
$\displaystyle\hat{a}^{\dagger}=\int_{0}^{\infty}d\omega\left(\alpha(\omega)\hat{B}^{\dagger}(\omega)-\beta^{*}(\omega)\hat{B}(\omega)\right).$
(31)

The requirement that these operators satisfy the familiar boson commutation
relation, $\left[\hat{a},\hat{a}^{\dagger}\right]=1$, provides a constraint on
the functions $\alpha(\omega)$ and $\beta(\omega)$ in the form

$\int_{0}^{\infty}d\omega\left[|\alpha(\omega)|^{2}-|\beta(\omega)|^{2}\right]=\int_{0}^{\infty}d\omega|\alpha(\omega)|^{2}\frac{4\Omega_{0}\omega}{(\Omega_{0}+\omega)^{2}}=1,$
(32)

where we have used (122). It is interesting to note that the correctness of
this may be verified explicitly by contour integration and that the proof
makes explicit use of the positivity condition (27) [64].

The integrand in (32) is clearly positive (or zero) for all frequencies and is
also normalized, and hence it has the mathematical form of a frequency
probability distribution:

$\pi(\omega)=|\alpha(\omega)|^{2}\frac{4\Omega_{0}\omega}{(\Omega_{0}+\omega)^{2}}.$
(33)

A number of further constraints on this quantity emerge naturally from
thinking of it as a probability distribution and from our diagonalization. We
note that the same concept has been described before: Georgievskii and Pollak
introduce an effective density of states for the diagonalized Caldeira-Leggett
model [94], while Ratchov _et al_ have expressed the properties of the damped
harmonic oscillator in terms of such a probability density [95], see also
[96].

###  3.4 Physical constraints

We leave until the next section a discussion of the physical interpretation of
the probability density $\pi(\omega)$ but consider here what can be inferred
from the fact that it has the mathematical properties of a probability
distribution. To this end let us denote the average value of a function of
frequency for this distribution by

$\langle\langle f(\omega)\rangle\rangle=\int_{0}^{\infty}d\omega
f(\omega)\pi(\omega).$ (34)

we note that $|\alpha(\omega)|^{2}$ is finite for all $\omega$ and it follows,
therefore, from (33) that the average value
$\langle\langle\omega^{-1}\rangle\rangle$ is finite.

It follows from the eigenoperator equations (3.3) that we can write the
Hamiltonian in the form

$\hat{H}=\int_{0}^{\infty}d\omega\>\hbar\omega\>\hat{B}^{\dagger}(\omega)\hat{B}(\omega)+C,$
(35)

where $C$ is an unimportant constant (even though it is formally divergent).
We can substitute into this Hamiltonian our expression for the eigenoperators
(30). We require that the coefficients of $\hat{a}^{2}$ and $\hat{a}^{\dagger
2}$ should vanish so that

$\int_{0}^{\infty}d\omega\>\omega|\alpha(\omega)|^{2}\left(\frac{\omega-\Omega_{0}}{\omega+\Omega_{0}}\right)=\int_{0}^{\infty}d\omega\frac{\pi(\omega)}{4\Omega_{0}}(\omega^{2}-\Omega_{0}^{2})=0\,,$
(36)

where we have again used (122). This implies that

$\langle\langle\omega^{2}\rangle\rangle=\Omega^{2}_{0}.$ (37)

The fact that the square of the mean value cannot exceed the mean value of the
square for any probability distribution leads us to deduce that

$\langle\langle\omega\rangle\rangle<\Omega_{0}.$ (38)

Finally we can apply the Cauchy-Schwartz inequality to provide a lower bound
on the value of $\langle\langle\omega^{-1}\rangle\rangle$:

$\displaystyle\langle\langle\omega^{-1}\rangle\rangle\langle\langle\omega\rangle\rangle\geq
1$
$\displaystyle\Rightarrow\quad\langle\langle\omega^{-1}\rangle\rangle>\frac{1}{\Omega_{0}}.$
(39)

These inequalities are useful in determining the properties of the ground
state.

##  4 Ground-state

It is immediately clear from the form of the Hamiltonian (3.2) that the ground
state of the oscillator is not that of the undamped oscillator, that is the
state annihilated by $\hat{a}$. To see this we need only note that the
interaction term $\hat{a}^{\dagger}\hat{b}^{\dagger}(\omega)$ acts on the
combined ground state of the non-interacting oscillator and reservoir, adding
a quantum to each. This suggests that the true ground state should be a
superposition of states with varying numbers of quanta in both the oscillator
and the reservoir and hence an entangled state. This situation is reminiscent
of the ground state of an atom in quantum electrodynamics, which is dressed by
virtual photons [97]. The dressing of the ground-state atom is responsible for
a number of important effects including the Casimir-Polder interaction [98,
99, 100] and the form of the polarizability of the atom [101, 102].It is
reasonable to expect that it will be similarly significant for our strongly
damped oscillator.

The true ground state, which we denote by the ket $|0\rangle$, is the zero-
eigenvalue eigenstate of all the annihilation operators $\hat{B}(\omega)$:

$\hat{B}(\omega)|0\rangle=0\qquad\qquad\forall\>\omega.$ (40)

The ground-state properties of the oscillator in this pure state are described
by a mixed state density operator obtained by tracing out the environmental
degrees of freedom:

$\hat{\rho}_{\rm Osc}={\rm Tr}_{\rm Env}\left(|0\rangle\langle 0|\right).$
(41)

The most straightforward way to determine the form of this mixed state is to
use the characteristic function [17]:

$\displaystyle\chi(\xi)$ $\displaystyle=$ $\displaystyle{\rm
Tr}\left[\hat{\rho}\exp\left(\xi\hat{a}^{\dagger}-\xi^{*}\hat{a}\right)\right]$
(42) $\displaystyle=$ $\displaystyle\langle
0|\exp\left(\xi\hat{a}^{\dagger}-\xi^{*}\hat{a}\right)|0\rangle.$

This function provides a complete description of the state and all of its
statistical properties. A brief summary of the principal properties of this
characteristic function is given in C. To evaluate the characteristic function
it we express $\hat{a}$ and $\hat{a}^{\dagger}$ in terms of the
eigenoperators333We make use of the operator ordering theorem [17]
$\displaystyle\exp(\hat{A}+\hat{B})=\exp(\hat{A})\exp(\hat{B})\exp\left(-\frac{1}{2}[\hat{A},\hat{B}]\right),$
which holds if the two operators $\hat{A}$ and $\hat{B}$ both commute with
their commutator $[\hat{A},\hat{B}]$.

$\displaystyle\chi(\xi)$ $\displaystyle=$ $\displaystyle\langle
0|\exp\left\\{\int_{0}^{\infty}d\omega\left[\left(\xi\alpha+\xi^{*}\beta\right)\hat{B}^{\dagger}(\omega)-\left(\xi^{*}\alpha^{*}+\xi\beta^{*}\right)\hat{B}(\omega)\right]\right\\}|0\rangle$
(43) $\displaystyle=$
$\displaystyle\exp\left(-\frac{1}{2}\int_{0}^{\infty}d\omega|\xi\alpha(\omega)+\xi^{*}\beta(\omega)|^{2}\right).$

This simple form, Gaussian in $\xi$, is characteristic of a squeezed thermal
state. We can rewrite this characteristic function in terms of our probability
density in the form

$\displaystyle\chi(\xi)$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2}\int_{0}^{\infty}d\omega\>\pi(\omega)\left(\frac{\omega}{\Omega_{0}}\xi_{r}^{2}+\frac{\Omega_{0}}{\omega}\xi_{i}^{2}\right)\right]$
(44) $\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2}\left(\frac{\langle\langle\omega\rangle\rangle}{\Omega_{0}}\xi_{r}^{2}+\langle\langle\omega^{-1}\rangle\rangle\Omega_{0}\>\xi_{i}^{2}\right)\right],$

where $\xi_{r}$ and $\xi_{i}$ are the real and imaginary parts of $\xi$
respectively. We note, in passing, that the quadrature operators for the
oscillator (familiar from quantum optics [17]) have unequal uncertainties:

$\displaystyle\Delta\left(\frac{\hat{a}+\hat{a}^{\dagger}}{\sqrt{2}}\right)=\sqrt{\frac{\langle\langle\omega^{-1}\rangle\rangle\Omega_{0}}{2}}$
$\displaystyle\Delta\left(\frac{-i(\hat{a}-\hat{a}^{\dagger})}{\sqrt{2}}\right)=\sqrt{\frac{\langle\langle\omega\rangle\rangle}{2\>\Omega_{0}}}.$
(45)

The product of these two exceeds $\frac{1}{2}$, as it must, by virtue of the
Cauchy-Schwartz inequality (3.4).

###  4.1 Ground-state energy

The ground-state energy of the undamped harmonic oscillator is simply
$\frac{1}{2}\hbar\omega_{0}$, where $\omega_{0}$ is the natural angular
frequency of the oscillator. For an underdamped oscillator, the oscillation
frequency is $\sqrt{\omega_{0}^{2}-\gamma^{2}/4}$, which is less than the
frequency of the undamped oscillator and this suggests, perhaps, that the
ground-state energy is similarly reduced [103]. Although plausible, this would
clearly run into difficulties in the over-damped regime in which this
characteristic frequency becomes imaginary. We shall see that the ground-state
energy of the oscillator, however we define it, _increases_ as a consequence
of the damping.

We have established that in the ground state, the reduced state of the
oscillator is Gaussian in position and momentum and, as such, is completely
determined by the first and second moments of the position and momentum. These
moments are

$\displaystyle\langle 0|\hat{x}|0\rangle=0$ $\displaystyle\langle
0|\hat{p}|0\rangle=0$ $\displaystyle\langle
0|\hat{x}^{2}|0\rangle=\frac{\hbar\langle\langle\omega^{-1}\rangle\rangle}{2m}$
$\displaystyle\langle 0|\hat{p}^{2}|0\rangle=\frac{\hbar
m\langle\langle\omega\rangle\rangle}{2}$ $\displaystyle\langle
0|\hat{p}\hat{x}+\hat{p}\hat{x}|0\rangle=0.$ (46)

Expressions of this form were derived previously by Grabert and Weiss in their
theory of the damped harmonic oscillator [42]. We note that the frequency of
the undamped oscillator, be it $\omega_{0}$ or $\Omega_{0}$, makes no explicit
appearance in these expectation values, but the averages
$\langle\langle\omega\rangle\rangle$ and
$\langle\langle\omega^{-1}\rangle\rangle$ will depend on these frequencies.

We have seen that the quantum damped harmonic oscillator is characterised by
two natural frequencies, $\omega_{0}$ and $\Omega_{0}$, and it is natural to
define the energy of the oscillator in terms of one or other of these:

$\hat{H}_{f_{0}}=\frac{\hat{p}^{2}}{2m}+\frac{1}{2}mf^{2}_{0}\hat{x}^{2}\,,$
(47)

where $f_{0}=\omega_{0}$ or $\Omega_{0}$. We shall consider a third possible
natural frequency, $\omega_{\rm diag}$, and the associated energy below. The
ground-state expectation values, Eq. (4.1), mean that we can write the
expectation value of the energy in the form

$\langle 0|\hat{H}_{f_{0}}|0\rangle=\frac{\hbar
f_{0}}{4}\left(\frac{\langle\langle\omega\rangle\rangle}{f_{0}}+\langle\langle\omega^{-1}\rangle\rangle
f_{0}\right)\,.$ (48)

Recall that the Cauchy-Schwartz inequality requires that

$\langle\langle\omega^{-1}\rangle\rangle\langle\langle\omega\rangle\rangle\geq
1\quad\Rightarrow\langle\langle\omega^{-1}\rangle\rangle\geq\frac{1}{\langle\langle\omega\rangle\rangle}$
(49)

and it follows that the ground-state energy is bounded by

$\langle 0|\hat{H}_{f_{0}}|0\rangle\geq\frac{\hbar
f_{0}}{4}\left(\frac{\langle\langle\omega\rangle\rangle}{f_{0}}+\frac{f_{0}}{\langle\langle\omega\rangle\rangle}\right)\,.$
(50)

The global minimum of this expression occurs for
$\langle\langle\omega\rangle\rangle=f_{0}$ and hence

$\langle 0|\hat{H}_{f_{0}}|0\rangle\geq\frac{1}{2}\hbar f_{0}$ (51)

for _any_ choice of frequency in the expression for the oscillator energy:

$\displaystyle\langle
0|\hat{H}_{\Omega_{0}}|0\rangle\geq\frac{1}{2}\hbar\Omega_{0}$
$\displaystyle\langle
0|\hat{H}_{\omega_{0}}|0\rangle\geq\frac{1}{2}\hbar\omega_{0}\,.$ (52)

The fact that both of these exceed the ground-state energy of the
corresponding undamped oscillator is a reflection of the fact that there is an
energy cost to be paid in order to decouple the oscillator from its
environment [12]. This increase conflicts, however, with the reduction in
ground-state energy that has previously been reported [103]. We note, further,
that the mean kinetic energy and potential energy for the oscillator alone do
not have the same value and that this is in marked contrast to the ground
state of the undamped oscillator.

###  4.2 Diagonal form of the oscillator ground state

We can diagonalise the density operator for the oscillator alone,
$\hat{\rho}_{\rm Osc}$, by means of a squeezing transformation [17] or,
equivalently, introducing a new pair of annihilation and creation operators
corresponding to a third candidate natural oscillation frequency, $\omega_{\rm
diag}$:

$\displaystyle\hat{c}=\sqrt{\frac{m\omega_{\rm
diag}}{2\hbar}}\left(\hat{x}+i\frac{\hat{p}}{m\omega_{\rm diag}}\right)$
$\displaystyle\hat{c}^{\dagger}=\sqrt{\frac{m\omega_{\rm
diag}}{2\hbar}}\left(\hat{x}-i\frac{\hat{p}}{m\omega_{\rm diag}}\right).$ (53)

To complete the diagonalisation we need only choose $\omega_{\rm diag}$ such
that the expectation values of $\hat{c}^{2}$ and $\hat{c}^{\dagger 2}$ are
zero:

$\displaystyle\langle 0|\hat{c}^{2}|0\rangle=\frac{m\omega_{\rm
diag}}{2\hbar}\langle 0|\hat{x}^{2}-\frac{\hat{p}^{2}}{m^{2}\omega_{\rm
diag}^{2}}|0\rangle=0$ $\displaystyle\Rightarrow\omega_{\rm
diag}=\sqrt{\frac{\langle\langle\omega\rangle\rangle}{\langle\langle\omega^{-1}\rangle\rangle}}.$
(54)

This frequency is the geometric mean of the two frequencies
$\langle\langle\omega\rangle\rangle$ and
$\langle\langle\omega^{-1}\rangle\rangle^{-1}$ and it is, by virtue of (38)
and (3.4), less than $\Omega_{0}$.

The mean number of $c$-quanta in the oscillator ground-state is

$\displaystyle\bar{n}_{c}=\langle
0|\hat{c}^{\dagger}\hat{c}|0\rangle=\frac{1}{2}\left(\sqrt{\langle\langle\omega\rangle\rangle\langle\langle\omega^{-1}\rangle\rangle}-1\right),$
(55)

which exceeds $0$, as it should, by virtue of (3.4). When written in terms of
the $c$-quanta, the steady-state density operator takes the form of a thermal
Bose-Einstein state, which we can write in the form

$\hat{\rho}_{\rm
Osc}=\frac{1}{\bar{n}_{c}+1}\left(\frac{\bar{n}_{c}}{\bar{n}_{c}+1}\right)^{\hat{c}^{\dagger}\hat{c}}.$
(56)

We may interpret this state as a thermal state for the oscillator at the
shifted frequency $\omega_{\rm diag}$ and at an effective “temperature”

$T_{\rm eff}=\frac{\hbar\omega_{\rm diag}}{k_{B}\ln(1+\bar{n}_{c}^{-1})}.$
(57)

We should note, however, that the true temperature in the ground-state is zero
and that this quantity and the frequency $\omega_{\rm diag}$ are at most only
parameters with which to quantify the state of the oscillator and its
entanglement with the environment. In particular, the von-Neumann entropy
associated with the steady state of the oscillator is

$S({\rm Osc})=(\bar{n}_{c}+1)\ln(\bar{n}_{c}+1)-\bar{n}_{c}\ln\bar{n}_{c}\,.$
(58)

By virtue of the Araki-Lieb inequality [104, 105, 106] and the fact that the
full state is pure, this means that this is also the total entropy of the
environment:

$S({\rm Env})=S({\rm Osc})$ (59)

and that the quantum mutual information [106], or index of correlation [107],
between the oscillator and the environment is

$S({\rm Osc}:{\rm Env})=S({\rm Env})+S({\rm Osc})-S({\rm Osc},{\rm
Env})=2S({\rm Osc}).$ (60)

##  5 Physical meaning of $\pi(\omega)$

We have seen that the physical properties of the oscillator ground state may
readily be expressed in terms of moments of the frequency given the
probability distribution $\pi(\omega)$. Here we present the case for
identifying this probability density with the contribution from the dressed
modes, associated with the eigenoperators $\hat{B}(\omega)$, to the state of
the oscillator. It is, in essence, in the spectrum of the true ground-state
(continuum) modes contributing to the oscillator state. We present three
arguments to support this interpretation.

Our first justification arises from the form of the expectation values (4.1).
We know that the ground state of a harmonic oscillator of frequency $\omega$
has

$\displaystyle\langle\hat{x}^{2}\rangle=\frac{\hbar}{2m\omega}$
$\displaystyle\langle\hat{p}^{2}\rangle=\frac{\hbar m\omega}{2}.$ (61)

If we treat the state of the oscillator as a mixture of oscillators of
different frequencies, each in its ground state, and contribution with weight
$\pi(\omega)$ then the resulting average mean-square values will be

$\displaystyle\langle\hat{x}^{2}\rangle=\int_{0}^{\infty}d\omega\>\pi(\omega)\frac{\hbar}{2m\omega}=\frac{\hbar\langle\langle\omega^{-1}\rangle\rangle}{2m}$
$\displaystyle\langle\hat{p}^{2}\rangle=\int_{0}^{\infty}d\omega\>\pi(\omega)\frac{\hbar
m\omega}{2}=\frac{\hbar m\langle\langle\omega\rangle\rangle}{2},$ (62)

which correspond to those obtained above. Note that the requirement that
$\langle\hat{x}^{2}\rangle$ must be finite imposes the condition that at zero
frequency

$\pi(0)=0.$ (63)

Our second point arises from the form of the Hamiltonian for the oscillator

$\hat{H}_{\Omega_{0}}=\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\Omega^{2}_{0}\hat{x}^{2}.$
(64)

We can, by virtue of (37), write this as a combination of potentials
corresponding to different frequencies but weighted by $\pi(\omega)$:

$\displaystyle\hat{H}_{\Omega_{0}}$ $\displaystyle=$
$\displaystyle\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\int_{0}^{\infty}d\omega\>\pi(\omega)\omega^{2}\hat{x}^{2}$
(65) $\displaystyle=$
$\displaystyle\frac{\hat{p}^{2}}{2m}+\frac{1}{2}m\langle\langle\omega^{2}\rangle\rangle\hat{x}^{2}.$

Finally we note that the effective mean energy of the oscillator, which is
associated with the diagonal form of the density operator (56) is

$\left(\bar{n}_{c}+\frac{1}{2}\right)\hbar\omega_{\rm
diag}=\frac{1}{2}\hbar\int_{0}^{\infty}d\omega\>\pi(\omega)\omega=\frac{1}{2}\hbar\langle\langle\omega\rangle\rangle,$
(66)

which neatly combines the characteristic ground state energies of the dressed
oscillators, weighted by the probability distribution $\pi(\omega)$. We note
that this mean energy is less than $\frac{1}{2}\hbar\Omega_{0}$ by virtue of
Eq (38), It is, however, always positive irrespective of whether the
oscillator evolution is under- or over-damped. The general question of whether
the ground-state energy of the damped oscillator is greater or less than that
of the undamped oscillator is difficult to answer definitively as there is no
unique form in the damped oscillator Hamiltonian for the free or undamped
oscillator frequency.

The combination of these three features (the expectation values
$\langle\hat{x}^{2}\rangle$ and $\langle\hat{p}^{2}\rangle$, the form of the
oscillator Hamiltonian and the effective mean energy) leads us to interpret
$\pi(\omega)$ as the proportion of the corresponding dressed oscillators
contributing to the properties of the damped oscillator. We emphasise that the
mathematical results obtained in the preceding section do not _require_ us to
adopt this interpretation of $\pi(\omega)$ but we find it helpful to do so.

##  6 Equilibrium state at finite temperature

We require the forms of the steady state and the dynamics for our damped
harmonic oscillator in an environment at finite temperature. In evaluating
this, we show that the equilibrium state is precisely the mean-force Gibbs
state found by tracing out the environmental degrees of freedom from the
global equilibrium Gibbs state. The analysis presents a problem as we cannot
assign a thermal state density operator in the continuum limit model of the
environment. We could take a step back and and return to a description in
terms of discrete reservoir modes, but it is more natural to adopt the
thermofield technique devised for the treatment of problems in finite
temperature quantum field theory [53, 54, 55, 56, 57, 58]. This is our third
less familiar element, and as thermofields will be a novelty for much of
intended readership, we present a brief account of the thermofield technique
in D.

It is _essential_ to realise, however, that we have not as yet established
that this equilibrium state is also the steady state of the strongly damped
oscillator. We turn to this in the following section.

We can determine the properties of the anticipated steady state by following
the same method as employed to study the ground state. The key idea is to
replace the vacuum state $|0\rangle$ of the coupled system, see Eq. (40) with
the thermal vacuum state in a doubled space, which is related to the true
vacuum state in the doubled space, $|0,\tilde{0}\rangle$, by a unitary
transformation:

$\displaystyle|0(\beta)\rangle$ $\displaystyle=$
$\displaystyle\hat{S}(\theta[\omega])|0,\tilde{0}\rangle$ (67)
$\displaystyle=$ $\displaystyle\exp\left[\int
d\omega\>\theta(\beta,\omega)\left(\hat{\tilde{B}}^{\dagger}(\omega)\hat{B}^{\dagger}(\omega)-\hat{B}(\omega)\hat{\tilde{B}}(\omega)\right)\right]|0,\tilde{0}\rangle\,.$

This state has the same single reservoir expectation values as the thermal
state and may therefore be used in its place. _If_ the coupled oscillator-
reservoir system relaxes to the thermal state of the coupled system (and we
have yet to establish this) then we can use this thermal vacuum state.

The corresponding thermal steady state of the harmonic oscillator, the mean
force Gibbs state [108, 109], will be a mixed state density operator obtained
by tracing over the environment:

$\hat{\rho}_{{\rm Osc},T}={\rm Tr}_{\rm Env}[|0(\beta)\rangle\langle
0(\beta)|]\,.$ (68)

As with the zero temperature ground state, we can determine the form of this
using the characteristic function:

$\displaystyle\chi_{T}(\xi)$ $\displaystyle=$ $\displaystyle\langle
0(\beta)|\exp\left(\xi\hat{a}^{\dagger}-\xi^{*}\hat{a}\right)|0(\beta)\rangle$
$\displaystyle=$ $\displaystyle\langle 0(\beta)|\exp\left\\{\int
d\omega\left[(\xi\alpha+\xi^{*}\beta)\hat{B}^{\dagger}(\omega)-(\xi^{*}\alpha^{*}+\xi\beta^{*})\hat{B}(\omega)\right]\right\\}|0(\beta)\rangle\,.$

We can transform this into a vacuum expectation value by applying a unitary
transformation to the annihilation and creation operators $\hat{B}(\omega)$
and $\hat{B}^{\dagger}(\omega)$:

$\displaystyle\hat{B}(\omega)$ $\displaystyle\rightarrow$
$\displaystyle\hat{B}(\omega){\rm
cosh}\theta(\beta,\omega)+\hat{\tilde{B}}^{\dagger}(\omega){\rm
sinh}\theta(\beta,\omega)$ $\displaystyle\hat{B}^{\dagger}(\omega)$
$\displaystyle\rightarrow$ $\displaystyle\hat{B}^{\dagger}(\omega){\rm
cosh}\theta(\beta,\omega)+\hat{\tilde{B}}(\omega){\rm
sinh}\theta(\beta,\omega)\,.$ (70)

Applying this transformation to our characteristic function replaces (6) by an
equivalent vacuum expectation value. Evaluating this gives

$\displaystyle\chi_{T}(\xi)$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{1}{2}\int
d\omega\>|\xi\alpha(\omega)+\xi^{*}\beta(\omega)|^{2}[{\rm
cosh}^{2}\theta(\beta,\omega)+{\rm sinh}^{2}\theta(\beta,\omega)]\right)$ (71)
$\displaystyle=$ $\displaystyle\exp\left(-\frac{1}{2}\int
d\omega\>|\xi\alpha(\omega)+\xi^{*}\beta(\omega)|^{2}{\rm
coth}(\beta\hbar\omega/2)\right)\,.$

As with the zero-temperature steady state, this is a simple Gaussian in $\xi$
and, again, is characteristic of a squeezed thermal state. When expressed in
terms of our probability density, $\pi(\omega)$, we find:

$\displaystyle\chi_{T}(\xi)$ $\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2}\int d\omega\>\pi(\omega){\rm
coth}(\beta\hbar\omega/2)\left(\frac{\omega}{\Omega_{0}}\xi_{r}^{2}+\frac{\Omega_{0}}{\omega}\xi_{i}^{2}\right)\right]$
$\displaystyle=$
$\displaystyle\exp\left[-\frac{1}{2}\left(\frac{\langle\langle\omega{\rm
coth}(\beta\hbar\omega/2)\rangle\rangle}{\Omega_{0}}\xi_{r}^{2}+\Omega_{0}\langle\langle\omega^{-1}{\rm
coth}(\beta\hbar\omega/2)\rangle\rangle\xi_{i}^{2}\right)\right]\,.$

We note that this has the same general form as the characteristic function for
the ground state, Eq (44), but with the probability density $\pi(\omega)$
replaced by a thermally-weighted density $\pi(\omega){\rm
coth}(\beta\hbar\omega/2)$. With this substitution, we can simply modify the
properties of the ground state so that, for example, the lowest moments of the
position and momentum operators in this state become

$\displaystyle\langle\hat{x}\rangle_{T}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\langle\hat{p}\rangle_{T}$ $\displaystyle=$ $\displaystyle 0$
$\displaystyle\langle\hat{x}^{2}\rangle_{T}$ $\displaystyle=$
$\displaystyle\frac{\hbar\langle\langle\omega^{-1}{\rm
coth}(\beta\hbar\omega/2)\rangle\rangle}{2m}$
$\displaystyle\langle\hat{p}^{2}\rangle_{T}$ $\displaystyle=$
$\displaystyle\frac{\hbar m\langle\langle\omega{\rm
coth}(\beta\hbar\omega/2)\rangle\rangle}{2}$
$\displaystyle\langle\hat{x}\hat{p}+\hat{p}\hat{x}\rangle_{T}$
$\displaystyle=$ $\displaystyle 0.$ (73)

We have seen that the oscillator is characterised by (at least) two different
natural frequencies, $\Omega_{0}$ and $\omega_{0}$. In terms of these, the
mean energy of the oscillator alone is

$\displaystyle\left\langle\frac{\hat{p}^{2}}{2m}+\frac{mf_{0}^{2}\hat{x}^{2}}{2}\right\rangle_{T}=\frac{\hbar
f_{0}}{4}\left(\frac{\langle\langle\omega{\rm
coth}(\beta\hbar\omega/2)\rangle\rangle}{f_{0}}+\langle\langle\omega^{-1}{\rm
coth}(\beta\hbar\omega/2)\rangle\rangle f_{0}\right),$

where $f_{0}=\omega_{0}$ or $\Omega_{0}$. This is the natural generalisation
of the ground-state energy of the oscillator given in Eq. (48)

It is interesting to pause at this point and to consider the behaviour of the
oscillator kinetic and potential energies in the high temperature limit. For a
weakly damped oscillator, we would expect both of these quantities to approach
$\frac{1}{2}k_{B}T$, the value suggested by the equipartition of energy. To
check this, we need only note that in the high temperature limit

${\rm
coth}(\beta\hbar\omega/2)\rightarrow\frac{2}{\beta\hbar\omega}=\frac{2k_{B}T}{\hbar\omega}\,.$
(75)

It follows that the high temperature limits of the kinetic and potential
energies are, respectively:

$\displaystyle\left\langle\frac{\hat{p}^{2}}{2m}\right\rangle_{T\rightarrow\infty}$
$\displaystyle=$ $\displaystyle\frac{k_{B}T}{2}$
$\displaystyle\left\langle\frac{mf_{0}^{2}\hat{x}^{2}}{2}\right\rangle_{T\rightarrow\infty}$
$\displaystyle=$
$\displaystyle\frac{k_{B}T}{2}f_{0}^{2}\langle\langle\omega^{-2}\rangle\rangle\,.$
(76)

The kinetic energy of the oscillator tends to the expected high-temperature
value, but the potential energy does not, and requires an explanation. In
pursuit of this, we note that the Cauchy-Schwartz inequality requires that
$\langle\langle\omega^{2}\rangle\rangle\langle\langle\omega^{-2}\rangle\rangle\geq
1$ and, as $\langle\langle\omega^{2}\rangle\rangle=\Omega_{0}^{2}$, it follows
that $\Omega_{0}^{2}\langle\langle\omega^{-2}\rangle\rangle\geq 1$, so that
the potential energy, when expressed in terms of $\Omega_{0}$, exceeds that
assigned by equipartition. The natural way to understand this is that the
oscillator is strongly rather than weakly coupled to its environment and the
excess thermal energy has its origin in the interaction energy with the
environment. The issue is less clear, however, if we express the potential
energy in terms of $\omega_{0}$. We shall see below that in the limit of weak
damping, when $\Omega_{0}$ and $\omega_{0}$ tend to a common value, the
probability distribution $\pi(\omega)$ becomes sharply peaked around
$\omega=\Omega_{0}$ so that the equipartition of energy for the potential
energy is restored.

##  7 Oscillator dynamics

It remains to consider the evolution of the oscillator towards equilibrium.
This will be important in practical applications of the theory but also for a
fundamental reason; we have obtained equilibrium states at zero and at finite
temperature, but have not as yet proven that the dynamics of the oscillator
causes it to evolve towards this state. Establishing this, without
approximation, is a principal aim of this section.

The exact diagonalisation of the Hamiltonian makes it straightforward to
evaluate the time-evolution of any desired property of the oscillator. All
that we need do is to express the desired observable in terms of the
eigenoperators, $\hat{B}(\omega)$ and $\hat{B}^{\dagger}(\omega)$ and then use
the time evolution of these operators, the form of which is an elementary
consequence of the fact that they are energy eigenoperators:

$\displaystyle\hat{B}(\omega,t)=\hat{B}(\omega,0)e^{-i\omega t}$
$\displaystyle\hat{B^{\dagger}}(\omega,t)=\hat{B}^{\dagger}(\omega,0)e^{i\omega
t}.$ (77)

In particular, we can determine the time-evolution of the annihilation
operator for the oscillator in this way:

$\displaystyle\hat{a}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}d\omega\left[\alpha^{*}(\omega)\hat{B}(\omega,0)e^{-i\omega
t}-\beta(\omega)\hat{B}^{\dagger}(\omega,0)e^{i\omega t}\right]$ (78)
$\displaystyle=$
$\displaystyle\left.\int_{0}^{\infty}d\omega\right\\{\alpha^{*}(\omega)e^{-i\omega
t}\left[\alpha(\omega)\hat{a}(0)+\beta(\omega)\hat{a}^{\dagger}(0)+\right.$
$\displaystyle\qquad\int_{0}^{\infty}d\omega^{\prime}\left.\>\left(\gamma(\omega,\omega^{\prime})\hat{b}(\omega^{\prime},0)+\delta(\omega,\omega^{\prime})\hat{b}^{\dagger}(\omega^{\prime},0)\right)\right]$
$\displaystyle\qquad-\beta(\omega)e^{i\omega
t}\left[\alpha^{*}(\omega)\hat{a}^{\dagger}(0)+\beta^{*}(\omega)\hat{a}(0)+\right.$
$\displaystyle\qquad\left.\int_{0}^{\infty}d\omega^{\prime}\left.\>\left(\gamma^{*}(\omega,\omega^{\prime})\hat{b}^{\dagger}(\omega^{\prime},0)+\delta^{*}(\omega,\omega^{\prime})\hat{b}(\omega^{\prime},0)\right)\right]\right\\}.$

This may be used, together with the initial state of the oscillator and the
environment, to evaluate the expectation value of any desired property. Note
that we have chosen the initial state to be one in which the oscillator and
the reservoir are uncorrelated. We allow the oscillator to be prepared in any
chosen state, but the reservoir is in a thermal state, which we describe using
a thermal vacuum state for the reservoir degrees of freedom, as given in Eq
(147). As the environment is in a stationary state, so that
$\langle\hat{b}(\omega,0)\rangle=0=\langle\hat{b}^{\dagger}(\omega,0)\rangle$,
the expectation values of the position and momentum operators take a
pleasingly simple form:

$\displaystyle\langle\hat{x}(t)\rangle=\langle\langle\cos(\omega
t)\rangle\rangle\langle\hat{x}(0)\rangle+\frac{1}{m}\langle\langle\omega^{-1}\sin(\omega
t)\rangle\rangle\langle\hat{p}(0)\rangle$
$\displaystyle\langle\hat{p}(t)\rangle=\langle\langle\cos(\omega
t)\rangle\rangle\langle\hat{p}(0)\rangle-m\langle\langle\omega\sin(\omega
t)\rangle\rangle\langle\hat{x}(0)\rangle,$ (79)

where the double angle brackets denote, as before, averages over our
probability distribution $\pi(\omega)$ as in (34). The generality of this
evolution follows simply from the linearity of the dynamics and has been noted
before, in particular by Haake and Reibold in their treatment of an oscillator
coupled to a quasi-continuum of oscillators [60]. The form of these equations
adds further support to our interpretation of $\pi(\omega)$ as a frequency
probability distribution for the damped oscillator, as they may be viewed as
the evolution of an undamped oscillator with a frequency $\omega$ averaged
using this probability distribution. The dissipation arises simply from a
dephasing amongst the different frequency components.

The evolution of the mean position and momentum, as given in (7), has the
necessary short-time form of that for an undamped oscillator

$\displaystyle\langle\hat{x}(\delta
t)\rangle=\langle\hat{x}(0)\rangle+\frac{\langle\hat{p}(0)\rangle}{m}\delta t$
$\displaystyle\langle\hat{p}(\delta
t)\rangle=\langle\hat{p}(0)\rangle-m\Omega_{0}^{2}\langle\hat{x}(0)\rangle\delta
t,$ (80)

where we have used the identity
$\langle\langle\omega^{2}\rangle\rangle=\Omega_{0}^{2}$. The effects of the
coupling to the environment enter at order $\delta t^{3}$ and this is an
indication of the essentially non-Markovian nature of the strongly-damped
oscillator. Our primary interest is in strongly damped oscillators and so we
should note that (7) includes the possibilities of both critically-damped and
over-damped evolution. The equations contain, moreover, a simple criterion for
these, which we may express in terms of our probability distribution. The
motion will be oscillatory if $\langle\langle\cos(\omega t)\rangle\rangle$ has
stationary points at times other than at $t=0$. Alternatively, we may state
that the motion is under-damped if the derivative of this quantity, that is
$\langle\langle\omega\sin(\omega t)\rangle\rangle$, is zero for any time other
than $t=0$. If it is zero only at $t=0$ then the motion is critically-damped
or over-damped.

###  7.1 Steady state

Our expression for the evolved annihilation operator (78), together with the
corresponding one for the creation operator provide a full description of the
oscillator dynamics. This is true for any initial state of the oscillator and,
moreover, for any environmental state including, of course, that associated
with a finite temperature. As an illustration let us examine the evolution of
the characteristic function foran arbitrary initial state of the oscillator
coupled to a finite-temperature environment at time $t=0$. The zero-
temperature behaviour follows, simply, in the limit $T\rightarrow 0$ or
$\beta\rightarrow\infty$. With a little effort we find (using the method of
characteristics [17])

$\chi_{T}[\xi,t]=\chi_{T}[\xi(t),0]\exp\left[-\frac{1}{2}\int_{0}^{\infty}d\omega^{\prime}\left|\int_{0}^{\infty}d\omega\>\mu(\omega,\omega^{\prime},t)\right|^{2}\right],$
(81)

where

$\displaystyle\xi(t)=\int_{0}^{\infty}d\omega\left[\eta(\omega,t)\alpha^{*}(\omega)-\eta^{*}(\omega,t)\beta(\omega)\right]$
$\displaystyle\mu(\omega,\omega^{\prime},t)=\eta(\omega,t)\gamma^{*}(\omega,\omega^{\prime})-\eta^{*}(\omega,t)\delta^{*}(\omega,\omega^{\prime})$
$\displaystyle\eta(\omega,t)=\left[\xi\alpha(\omega)+\xi^{*}\beta(\omega)\right]{\rm
coth}(\beta\hbar\omega/2)e^{i\omega t}.$ (82)

This characteristic function encodes the full dynamics and statistics of the
oscillator. As a simple illustration of this we can determine, directly, the
form of the steady state. To see this we first note $\xi(t)$ tends to zero as
$t$ tends to infinity and the different frequency components dephase so that
the prefactor in Eq (81), corresponding to the initial state of the oscillator
tends to

$\chi[\xi(\infty),0]=\chi[0,0]=1,$ (83)

which means that all memory of the initial state of the oscillator is lost.
Evaluating the long-time limit of the exponential factor in (81) requires some
care in the handling of the delta-function and principal part components. We
find

$\int_{0}^{\infty}d\omega^{\prime}\left|\int_{0}^{\infty}d\omega\>\mu(\omega,\omega^{\prime},\infty)\right|^{2}=\int_{0}^{\infty}d\omega\left|\xi\alpha(\omega^{\prime})+\xi^{*}\beta(\omega^{\prime})\right|^{2}{\rm
coth}(\beta\hbar\omega/2),$ (84)

so the steady-state characteristic function is

$\chi(\xi,\infty)=\exp\left(-\frac{1}{2}\int_{0}^{\infty}d\omega\left|\xi\alpha(\omega)+\xi^{*}\beta(\omega)\right|^{2}{\rm
coth}(\beta\hbar\omega/2)\right),$ (85)

which we recognise as the characteristic function for the oscillator in the
global thermal equilibrium state (71). This is a most satisfactory and exact
result.

It also means that the steady state of the oscillator is the mean-force Gibbs
state. To see this we need only note that it is given by the trace over the
environment of the full thermal equilibrium state, Eq (68). Proof of this
equivalence has also been shown in [59] by demonstrating the equality of
steady state multi-time open system correlation functions obtained by
Heisenberg-Langevin equation of motion methods to those of the closed system
thermal Gibbs state.

###  7.2 An example evolution

We have seen that both the dynamics and the steady-state of our damped
harmonic oscillator are governed by the form of the function $\pi(\omega)$.
Determining this, together with the initial conditions, provides all the
information required. We can calculate $\pi(\omega)$ directly from the
frequency-dependence of the coupling between the oscillator and its
environment or, more simply, select a form for $\pi(\omega)$ and proceed from
this. This is the approach we adopt here.

As an example we consider a $\pi(\omega)$ of the form

$\pi(\omega)=\frac{2\omega^{2}}{\pi}\cdot\frac{(\gamma_{+}+\gamma_{-})(\gamma_{-}+\Gamma)(\Gamma+\gamma)}{(\omega^{2}+\Gamma^{2})(\omega^{2}+\gamma_{+}^{2})(\omega^{2}+\gamma_{-}^{2})}\,.$
(86)

This is, perhaps, the simplest example that satisfies the necessary physical
constraints: (i) it is normalised, (ii) $\pi(0)=0$, and (iii)
$\langle\langle\omega^{2}\rangle\rangle$ is finite. We would like our example
to lead to an evolution that is close to the familiar classical behaviour. To
this end we choose $\Gamma$ to be a real decay rate and $\gamma_{\pm}$ are
either two real decay rates or complex conjugates of one another. With the
familiar classical behaviour in mind we set

$\gamma_{\pm}=\frac{\gamma}{2}\pm\sqrt{\frac{\gamma^{2}}{4}-\omega_{0}^{2}}.$
(87)

Direct evaluation of $\langle\langle\omega^{2}\rangle\rangle$ for our example
gives

$\langle\langle\omega^{2}\rangle\rangle=\Omega_{0}^{2}=\Gamma(\gamma_{+}+\gamma_{-})+\gamma_{+}\gamma_{-}\,,$
(88)

which links the parameters $\Gamma$, $\gamma_{+}$ and $\gamma_{-}$ to the
short-time natural frequency $\Omega_{0}$. When written in terms of
$\omega_{0}$ and $\gamma$, this relationship simplifies to

$\Omega_{0}^{2}=\omega_{0}^{2}+\Gamma\gamma\,.$ (89)

We can rewrite this expression in the form

$\Gamma=\frac{\Omega_{0}^{2}-\omega_{0}^{2}}{\gamma}\,,$ (90)

which emphasises the need for the decay rate $\Gamma$ as a consequence of the
fact that we require two natural frequencies, $\omega_{0}$ and $\Omega_{0}$,
to describe the quantum damped harmonic oscillator. It is straightforward to
determine from $\pi(\omega)$ the two averages
$\langle\langle\omega\rangle\rangle$ and
$\langle\langle\omega^{-1}\rangle\rangle$:

$\displaystyle\langle\langle\omega\rangle\rangle$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\cdot\frac{\gamma_{+}^{2}\gamma_{-}^{2}\ln(\gamma_{-}/\gamma_{+})+\gamma_{-}^{2}\Gamma^{2}\ln(\Gamma/\gamma_{-})+\Gamma^{2}\gamma_{+}^{2}\ln(\gamma_{+}/\Gamma)}{(\Gamma-\gamma_{+})(\gamma_{+}-\gamma_{-})(\gamma_{-}-\Gamma)}$
$\displaystyle\langle\langle\omega^{-1}\rangle\rangle$ $\displaystyle=$
$\displaystyle\frac{2}{\pi}\cdot\frac{\Gamma^{2}\ln(\gamma_{-}/\gamma_{+})+\gamma_{+}^{2}\ln(\Gamma/\gamma_{-})+\gamma_{-}^{2}\ln(\gamma_{+}/\Gamma)}{(\Gamma-\gamma_{+})(\gamma_{+}-\gamma_{-})(\gamma_{-}-\Gamma)}$
(91)

from which we can derive the expectation values of $\hat{x}^{2}$ and
$\hat{p}^{2}$ in the ground state of the damped oscillator 444When
$\gamma_{\pm}$ are complex, we can write these in the form
$\gamma_{\pm}=|\gamma_{+}|e^{\pm i\phi}$, where
$\phi=\tan^{-1}[\Im(\gamma_{+})/\Re(\gamma_{+})]$..

We can determine the evolution of the expectation values of the position and
momentum from Eq (7) by evaluating the averages $\langle\langle\cos(\omega
t)\rangle\rangle$, $\langle\langle\omega^{-1}\sin(\omega t)\rangle\rangle$ and
$\langle\langle\omega\sin(\omega t)\rangle\rangle$:

$\displaystyle\langle\langle\cos(\omega t)\rangle\rangle$ $\displaystyle=$
$\displaystyle\frac{\Gamma(\gamma_{+}+\gamma_{-})}{(\Gamma-\gamma_{+})(\gamma_{-}-\Gamma)}e^{-\Gamma
t}+\frac{\gamma_{+}(\Gamma+\gamma_{-})}{(\Gamma-\gamma_{+})(\gamma_{+}-\gamma_{-})}e^{-\gamma_{+}t}$
$\displaystyle\qquad\qquad+\frac{\gamma_{-}(\Gamma+\gamma_{+})}{(\gamma_{-}-\Gamma)(\gamma_{+}-\gamma_{-})}e^{-\gamma_{-}t}$
$\displaystyle\langle\langle\omega^{-1}\sin(\omega t)\rangle\rangle$
$\displaystyle=$ $\displaystyle\int_{0}^{t}\langle\langle\cos(\omega
t^{\prime})\rangle\rangle dt^{\prime}$
$\displaystyle\langle\langle\omega\sin(\omega t)\rangle\rangle$
$\displaystyle=$ $\displaystyle-\frac{d}{dt}\langle\langle\cos(\omega
t)\rangle\rangle\,.$ (92)

We expect to find something approximating the evolution of the classical
damped harmonic oscillator in the limit of small $\Gamma$, but departures from
this for larger values. In figure 2 we plot the evolution of
$\langle\langle\cos(\omega t)\rangle\rangle$ for a small value of $\Gamma$. As
anticipated, we find over-damped like behaviour for real values of
$\gamma_{+}$ and $\gamma_{-}$, and underdamped behaviour for complex values.
Note, however, that there is a small overshoot in the underdamped regime,
which would certainly be absent in the over-damped evolution of the classical
oscillator. This can be traced back to the $e^{-\Gamma t}$ term which,
although small, decays slowly and so has a residual influence at long times.
This is especially clear in the large $\Gamma$ regime , depicted in figure 3.
There we see that there is a significant over-shoot of the mean position in
what, classically, would be the over-damped regime. In the under-damped
regime, however, the $e^{-\Gamma t}$ term has a less dramatic effect; the
evolution for $\Gamma=0.01$ (in figure 2b) and for $\Gamma=10$ (figure 3b) are
qualitatively rather similar.

Figure 2: Evolution of $\langle\langle\cos(\omega t)\rangle\rangle$ for small
$\Gamma$. a) Overdamped regime, with parameters $\Gamma=0.01\gamma$,
$\gamma_{+}=\frac{3}{4}\gamma$, $\gamma_{-}=\frac{1}{4}\gamma$; b) Underdamped
regime with parameters $\Gamma=0.01\gamma$, $\gamma_{\pm}=\left(\frac{1}{2}\pm
5i\right)\gamma$.

Figure 3: Evolution of $\langle\langle\cos(\omega t)\rangle\rangle$ for large
$\Gamma$. a) Overdamped regime, with parameters $\Gamma=10\gamma$,
$\gamma_{+}=\frac{3}{4}\gamma$, $\gamma_{-}=\frac{1}{4}\gamma$; b) Underdamped
regime with parameters $\Gamma=10\gamma$, $\gamma_{\pm}=\left(\frac{1}{2}\pm
5i\right)\gamma$.

We have seen that the existence of two natural frequencies for the oscillator,
$\omega_{0}$ and $\Omega_{0}$, is particularly important when comparing the
short and long time behaviour of the oscillator: $\omega_{0}$ behaves as the
natural frequency, but at short times it is $\Omega_{0}$ that takes on this
role. In figure 4 we plot the short-time evolution of
$\langle\langle\cos(\omega t)\rangle\rangle$ (solid line) and compare this
with the classical evolution of for a classical damped harmonic oscillator
with natural frequency $\omega_{0}$ It is clear that the former falls off more
quickly as it must, because $\Omega_{0}>\omega_{0}$. For the parameters
chosen, $\Omega_{0}$ is just over twice the value of $\omega_{0}$. It is
interesting to note that we can match the two evolutions up to very short
times if we allow for a short-time slip so that as $t$ tends to zero,
$\langle\langle\cos(\omega t)\rangle\rangle$ becomes larger than zero [60].
This is unnecessary, however, if we take account of the existence of two
natural frequencies as we have done here.

Figure 4: Comparison of short-time behaviour in quantum and classical cases.
The exact solution for short times is plotted as a solid line, with
$\pi(\omega)$ as given in Eq. 86, and parameters $\Gamma=10\gamma$,
$\gamma_{+}=\frac{3}{4}\gamma$, $\gamma_{-}=\frac{1}{4}\gamma$; the dashed
line shows the limit $\Gamma\rightarrow 0$, in which classical damped simple
harmonic motion is recovered, with the other parameters unchanged,
$\gamma_{\pm}=\left(\frac{1}{2}\pm 5i\right)\gamma$

##  8 Weak-coupling limit

The theory developed above was designed to treat the strongly-damped harmonic
oscillator, but should also be applicable to the more familiar weakly damped
oscillator, for which the oft-employed Born and Markov approximations are
applicable and the steady state of the oscillator should be its ground state.
We show here that this is indeed the case.

We start by considering the form of the function $\alpha(\omega)$ in the weak
coupling limit. To aid our analysis we rewrite the form given in (130) as

$\alpha(\omega)=\frac{\omega+\Omega_{0}}{\Omega_{0}}\left(\frac{V(\omega)}{|V(\omega)|^{2}Y(\omega)-i\pi|V(\omega)|^{2}}\right).$
(93)

The weak damping limit corresponds to choosing the coupling to the environment
to be small or, more specifically, to $|V(\omega)|^{2}\ll\Omega_{0}$. It is
clear that in this limit, $|\alpha(\omega)|^{2}$ will be a sharply peaked
function centred around the frequency for which $Y(\omega)=0$. If the integral
part in $Y(\omega)$, as given in (B) is small555If it is not then we will need
to invoke ideas of renormalisation, a manageable complication, but one we wish
to avoid. then this frequency will be close to $\Omega_{0}$ and we can write

$|V(\omega)|^{2}Y(\omega)\approx 4(\omega-\Omega_{0})-4F(\omega),$ (94)

where

$4F(\omega)=\int_{0}^{\infty}d\omega^{\prime}\left(\frac{\mathbbmss{P}}{\omega-\omega^{\prime}}-\frac{1}{\omega+\omega^{\prime}}\right)|V(\omega^{\prime})|^{2}.$
(95)

This leads, in turn, to a corresponding approximate form for $\alpha(\omega)$:

$\alpha(\omega)\approx\frac{V(\omega)}{2}\frac{1}{\omega-\Omega_{0}-F(\omega)-i\frac{\pi}{4}|V(\omega)|^{2}},$
(96)

where we have set $\omega=\Omega_{0}$ everywhere except in the denominator. We
note that this is of the form that arises from the Fano diagonalisation of our
problem if we make the rotating wave approximation by omitting from our
original Hamiltonian all terms that are products of two creation operators or
of two annihilation operators [17, 93].

Consistency with the above approximation, which led us to set
$\omega=\Omega_{0}$ leads us to set $\beta(\omega)$ to zero:

$\beta(\omega)\approx 0,$ (97)

so that the integral over all frequency of $|\alpha(\omega)|^{2}$ is unity.
Moreover, for weak damping the thermal function ${\rm
coth}(\beta\hbar\omega/2)$ will also be slowly varying compared to the rapid
variation of $|\alpha(\omega)|^{2}$ in the vicinity of $\omega=\Omega_{0}$ and
we may replace this function by its value at $\Omega_{0}$. Hence in this limit
the steady-state characteristic function for our oscillator is

$\chi(\xi,\infty)=\exp\left(-\frac{1}{2}|\xi|^{2}{\rm
coth}(\beta\hbar\Omega_{0}/2)\right)=\exp\left(-\frac{1}{2}|\xi|^{2}(2\bar{n}+1)\right),$
(98)

where $\bar{n}$ is the mean thermal excitation number. We recognise (98) as
the symmetrically ordered characteristic function for the thermal state of the
undamped oscillator [17], as it must be. Further, we note that in this limit
our probability distribution function,
$\pi(\omega)\approx|\alpha(\omega)|^{2}$, approaches a Lorentzian centred on
$\Omega_{0}$, with some width $\gamma$ 666This is not strictly true in the
wings of the distribution, of course, as even in this limit we require
$\langle\langle\omega^{2}\rangle\rangle=\Omega_{0}^{2}$, but the corresponding
quantity for a true Lorentzian is divergent.. Thus all the complexity of of
the original problem is reduced, in the weak-damping limit to just three
parameters: a natural oscillation frequency, a damping rate and a temperature.

It was important to confirm that our more general treatment coincided, in the
right limit, with the approximate methods used for weakly damped oscillators.
We should note that even if we are working in the weakly damped regime, then
our approach offers a systematic way to treat corrections to the results
obtained using the Born and Markov approximations, which may play an important
role in modelling measurements at the limits of sensitivity.

##  9 Conclusion

We have presented an exact diagonalisation of a simple quantum model of the
damped harmonic oscillator, one that is applicable, in particular, to any
strength of the damping. As a result we have recovered the fact that much of
the behaviour of the oscillator and many of its properties can be described in
terms of a single probability function, $\pi(\omega)$, which we may interpret
as the contribution of corresponding dressed mode, at frequency $\omega$ to
the oscillator. These properties include the steady state of the oscillator at
both at zero and at finite temperature, the entanglement between the
oscillator and its environment and also its evolution, both in the familiar
under-damped regime but also in the more problematic over-damped regime.

We have applied our diagonalisation to study the properties of the true ground
state and have shown that the oscillator part of this pure entangled state
coincides with the steady-state of the oscillator in a zero-temperature
environment. The diagonalisation is not specific to any particular state of
the reservoir, however, and we have shown how it can be be applied to
environments at finite temperature. The extension to more exotic states, such
as squeezed reservoirs presents no obvious difficulties. It may be extended,
moreover, to include driving forces, coupled oscillators and multiple
reservoirs, with the latter perhaps being at different temperatures [110].
This may provide some insights into important questions of principle in the
nascent fields of quantum machines and quantum thermodynamics [10, 11, 111,
112, 113, 114, 115].

It is a pleasure to dedicate this paper to our friend and colleague, Igor Jex,
in celebration of his 60th birthday. The Hamiltonian diagonalisation upon
which much of this work is based was first calculated by Bruno Huttner, in
collaboration with SMB, 30 years ago in their work on the quantum
electrodynamics of dielectric media. We are most grateful to him and also to
Paul Radmore and Claire Gilson for helpful comments and suggestions. This work
was supported, in part, by the Royal Society through the award to SMB of a
Research Professorship, RP150122.

##  Appendix A Heisenberg equations of motion

The Heisenberg equations of motion follow directly from the Hamiltonian, and
we obtain these using the Hamiltonian in the form of (2). (We could equally
well have used the identical Hamiltonian (3).)

$\displaystyle\dot{\hat{x}}$ $\displaystyle=$ $\displaystyle\frac{\hat{p}}{m}$
$\displaystyle\dot{\hat{p}}$ $\displaystyle=$
$\displaystyle-m\Omega_{0}^{2}\hat{x}+\sum_{\mu}m_{\mu}\omega_{\mu}^{2}\hat{x}_{\mu}$
$\displaystyle\dot{\hat{x}}_{\mu}$ $\displaystyle=$
$\displaystyle\frac{\hat{p}_{\mu}}{m_{\mu}}$
$\displaystyle\dot{\hat{p}}_{\mu}$ $\displaystyle=$ $\displaystyle-
m_{\mu}\omega_{\mu}^{2}\hat{x}_{\mu}+m_{\mu}\omega^{2}_{\mu}\lambda_{\mu}\hat{x}\,.$
(99)

We seek an equation of motion for the position operator and so first eliminate
the momentum operators between the first and second and the third and fourth
equations:

$\displaystyle\ddot{\hat{x}}+\omega_{0}^{2}\hat{x}$ $\displaystyle=$
$\displaystyle\sum_{\mu}\frac{m_{\mu}}{m}\omega^{2}_{\mu}\lambda_{\mu}\hat{x}_{\mu}$
(100) $\displaystyle\ddot{\hat{x}}_{\mu}+\omega_{\mu}^{2}\hat{x}_{\mu}$
$\displaystyle=$ $\displaystyle\omega^{2}_{\mu}\lambda_{\mu}\hat{x}\,.$ (101)

The next step is to integrate the second of these equations of motion. The
complementary function is

$\displaystyle\hat{x}_{\mu}^{\rm CF}$ $\displaystyle=$
$\displaystyle\hat{x}_{\mu}(0)\cos(\omega_{\mu}t)+\frac{\dot{\hat{x}}_{\mu}(0)}{\omega_{\mu}}\sin(\omega_{\mu}t)$
(102) $\displaystyle=$
$\displaystyle\hat{x}_{\mu}(0)\cos(\omega_{\mu}t)+\frac{\hat{p}_{\mu}(0)}{m_{\mu}\omega_{\mu}}\sin(\omega_{\mu}t)\,.$

To find the particular integral we need to make a small addition to (102) by
adding a very weak damping term to give:

$\ddot{\hat{x}}_{\mu}+\varepsilon\dot{\hat{x}}_{\mu}+\omega_{\mu}^{2}\hat{x}_{\mu}=\omega^{2}_{\mu}\lambda_{\mu}\hat{x}\,,$
(103)

and work in the limit as the strictly positive quantity $\varepsilon$ tends to
zero. It is this choice of a positive (if very small) value for $\varepsilon$
that provides the irreversibility and hence the arrow of time.

Solving (103) for the particular integral requires some care and so we pause
for a moment to provide the details. Let us introduce the Fourier transform of
$\hat{x}_{\mu}$ in the form

$\bar{x}_{\mu}(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\hat{x}_{\mu}(t)e^{i\omega
t}dt\,,$ (104)

with a similar expression for the Fourier transform of $\hat{x}$. It follows
that the Fourier transform of the particular integral part of the position is
given by

$\displaystyle-\omega^{2}\bar{x}_{\mu}-i\varepsilon\bar{x}_{\mu}+\omega_{\mu}^{2}\bar{x}_{\mu}$
$\displaystyle=$ $\displaystyle\omega_{\mu}^{2}\lambda_{\mu}\bar{x}$
$\displaystyle\qquad\qquad\Rightarrow\qquad\bar{x}_{\mu}$ $\displaystyle=$
$\displaystyle-\frac{\omega_{\mu}^{2}\lambda_{\mu}\bar{x}}{\omega^{2}-\omega_{\mu}^{2}+i\varepsilon\omega}\,.$
(105)

It follows then follows that

$\displaystyle\hat{x}_{\mu}^{\rm PI}(t)$ $\displaystyle=$
$\displaystyle-\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\frac{\omega_{\mu}^{2}\lambda_{\mu}}{\omega^{2}-\omega_{\mu}^{2}+i\varepsilon\omega}\bar{x}(\omega)e^{-i\omega
t}d\omega\,,$ (106)

which is the Fourier transform of the product of two functions of $\omega$. We
can use the convolution theorem to write this in terms of the transforms. To
exploit this we write the Fourier transform of the first function as the time-
derivative of a function $K(t)$:

$\displaystyle\dot{K}_{\mu}(t)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\infty}^{\infty}\frac{\omega_{\mu}^{2}\lambda_{\mu}}{\omega^{2}-\omega_{\mu}^{2}+i\varepsilon\omega}e^{-i\omega
t}d\omega$ (107) $\displaystyle=$
$\displaystyle-\omega_{\mu}\lambda_{\mu}e^{-\varepsilon
t/2}\sin(\omega_{\mu}t)$

if $t>0$ and is zero otherwise. We can now take the limit as
$\varepsilon\rightarrow 0$ to give

$\displaystyle\dot{K}_{\mu}$ $\displaystyle=$
$\displaystyle-\omega_{\mu}\lambda_{\mu}\sin(\omega_{\mu}t)$
$\displaystyle\Rightarrow K_{\mu}$ $\displaystyle=$
$\displaystyle\lambda_{\mu}\cos(\omega_{\mu}t)\,.$ (108)

It then follows that

$\displaystyle\hat{x}_{\mu}^{\rm PI}(t)$ $\displaystyle=$
$\displaystyle-\int_{-\infty}^{\infty}\bar{\dot{K}}_{\mu}(\omega)\bar{x}(\omega)e^{-i\omega
t}d\omega$ (109) $\displaystyle=$
$\displaystyle-\int_{-\infty}^{\infty}\dot{K}_{\mu}(T)\hat{x}(t-T)dT$
$\displaystyle=$
$\displaystyle-\left[K(T)\hat{x}(t-T)\right]_{0}^{t}+\int_{0}^{t}K_{\mu}(T)\frac{d}{dT}\hat{x}(t-T)dT$
$\displaystyle=$ $\displaystyle-
K_{\mu}(t)\hat{x}(0)+K_{\mu}(0)\hat{x}(t)-\int_{0}^{t}K_{\mu}(t-t^{\prime})\dot{\hat{x}}(t^{\prime})dt^{\prime}\,,$

where we have used the facts that $K(T)=0$ for $T<0$ and we may take
$\hat{x}(\tau)=0$ for $\tau<0$.

Pulling this altogether we arrive at our desired Heisenberg-Langevin equation
for the damped harmonic oscillator. From (101) we have

$\ddot{\hat{x}}(t)+\int_{0}^{t}\kappa(t-t^{\prime})\dot{\hat{x}}(t^{\prime})dt^{\prime}+\left(\Omega_{0}^{2}-\kappa(0)\right)\hat{x}(t)+\kappa(t)\hat{x}(0)=\frac{\hat{F}(t)}{m}\,,$
(110)

where

$\displaystyle\kappa(t)$ $\displaystyle=$
$\displaystyle\sum_{\mu}\frac{m_{\mu}}{m}\Omega_{\mu}^{2}\lambda_{\mu}K(t)$
(111) $\displaystyle=$
$\displaystyle\sum_{\mu}\frac{m_{\mu}}{m}\omega_{\mu}^{2}\lambda_{\mu}^{2}\cos(\omega_{\mu}t)$

and $F(t)$ is the Langevin force:

$\hat{F}(t)=\omega^{2}_{\mu}\lambda_{\mu}\left(\hat{x}_{\mu}(0)\cos(\omega_{\mu}t)+\frac{\hat{p}_{\mu}(0)}{m_{\mu}\omega_{\mu}}\sin(\omega_{\mu}t)\right)\,.$
(112)

##  Appendix B Fano diagonalisation

We present a brief account of the exact diagonalisation of our Hamiltonian
based on methods developed by Fano for the study of configuration interactions
[116]. This idea was applied extended to weakly-coupled oscillators in a
quantum study of damped cavity modes [17, 93]. The extension to stronger
couplings, with the inclusion of counter-rotating couplings, has been given
before and applied to the quantum theory of light in dielectric media [63,
64]. We summarise here the analysis presented in [64].

Our task is to diagonalise the damped harmonic oscillator Hamiltonian777For
the Hamiltonian of interest in this paper the coupling, $V(\omega)$, is real
but treating the problem with a more general complex coupling presents no
additional difficulties.

$\displaystyle\hat{H}=\hbar\Omega_{0}\hat{a}^{\dagger}\hat{a}+\int_{0}^{\infty}d\omega\>\hbar\omega\hat{b}^{\dagger}(\omega)\hat{b}(\omega)$
$\displaystyle\qquad\qquad+\int_{0}^{\infty}d\omega\>\frac{\hbar}{2}\left(\hat{a}+\hat{a}^{\dagger}\right)\left[V(\omega)\hat{b}^{\dagger}(\omega)+V^{*}(\omega)\hat{b}(\omega)\right],$
(113)

by which we mean rewriting it in the form of a continuum of _uncoupled_ or
dressed oscillators:

$\hat{H}=\int_{0}^{\infty}d\omega\>\hbar\omega\>\hat{B}^{\dagger}(\omega)\hat{B}(\omega)+C,$
(114)

where $C$ is an unimportant constant. We proceed by writing the dressed
annihilation operators, $\hat{B}(\omega)$, as linear combinations of the bare
operators for the oscillator and bath modes:

$\hat{B}(\omega)=\alpha(\omega)\hat{a}+\beta(\omega)\hat{a}^{\dagger}+\int_{0}^{\infty}d\omega^{\prime}\>\left[\gamma(\omega,\omega^{\prime})\hat{b}(\omega^{\prime})+\delta(\omega,\omega^{\prime})\hat{b}^{\dagger}(\omega^{\prime})\right],$
(115)

where $\alpha{\omega}$, $\beta(\omega)$, $\gamma(\omega,\omega^{\prime})$ and
$\delta(\omega,\omega^{\prime})$ are to be determined.

We require the operator $\hat{B}(\omega)$ to be associated with an uncoupled
or dressed oscillator of angular frequency $\omega$. This requires us to find
its form such that the following pair of operator equations are satisfied for
every frequency, $\omega$:

$\displaystyle\left[\hat{B}(\omega),\hat{H}\right]=\hbar\omega\hat{B}(\omega)$
(116)
$\displaystyle\left[\hat{B}(\omega),\hat{B}^{\dagger}(\omega^{\prime})\right]=\delta(\omega-\omega^{\prime}).$
(117)

Substituting the ansatz (115) into (116) and comparing coefficients of the
bare creation and annihilation operators leads to the set of coupled
equations:

$\displaystyle\alpha(\omega)\Omega_{0}+\frac{1}{2}\int_{0}^{\infty}d\omega^{\prime}\left[\gamma(\omega,\omega^{\prime})V(\omega^{\prime})-\delta(\omega,\omega^{\prime})V^{*}(\omega^{\prime})\right]=\alpha(\omega)\omega$
(118)
$\displaystyle-\beta(\omega)\Omega_{0}+\frac{1}{2}\int_{0}^{\infty}d\omega^{\prime}\left[\gamma(\omega,\omega^{\prime})V(\omega^{\prime})-\delta(\omega,\omega^{\prime})V^{*}(\omega^{\prime})\right]=\beta(\omega)\omega$
(119)
$\displaystyle\frac{V^{*}(\omega^{\prime})}{2}\left[\alpha(\omega)-\beta(\omega)\right]+\gamma(\omega,\omega^{\prime})\omega^{\prime}=\gamma(\omega,\omega^{\prime})\omega$
(120)
$\displaystyle\frac{V(\omega^{\prime})}{2}\left[\alpha(\omega)-\beta(\omega)\right]-\delta(\omega,\omega^{\prime})\omega^{\prime}=\delta(\omega,\omega^{\prime})\omega.$
(121)

Our method of solution is use these to determine the functions
$\beta(\omega)$, $\gamma(\omega,\omega^{\prime})$ and
$\delta(\omega,\omega^{\prime})$ in terms of $\alpha(\omega)$ and then to
determine this remaining function by enforcing the commutation relation (117).
From (118) and (119) we see that

$\beta({\omega})=\frac{\omega-\Omega_{0}}{\omega+\Omega_{0}}\>\alpha(\omega).$
(122)

If we use this to substitute for $\beta(\omega)$ into the remaining equations
then we find

$\displaystyle
V^{*}(\omega^{\prime})\frac{\Omega_{0}}{\omega+\Omega_{0}}\alpha(\omega)=\gamma(\omega,\omega^{\prime})(\omega-\omega^{\prime})$
(123) $\displaystyle
V(\omega^{\prime})\frac{\Omega_{0}}{\omega+\Omega_{0}}\alpha(\omega)=\delta(\omega,\omega^{\prime})(\omega+\omega^{\prime}).$
(124)

Solving the second of these presents no difficulty and we find

$\delta(\omega,\omega^{\prime})=\left(\frac{1}{\omega+\omega^{\prime}}\right)V(\omega^{\prime})\frac{\Omega_{0}}{\omega+\Omega_{0}}\alpha(\omega).$
(125)

The first, however, requires careful handling because the behaviour at
$\omega=\omega^{\prime}$. Following Fano [116], we adopt the method proposed
by Dirac [117] and write

$\gamma(\omega,\omega^{\prime})=\left(\frac{\mathbbmss{P}}{\omega-\omega^{\prime}}+Y(\omega)\delta(\omega-\omega^{\prime})\right)V^{*}(\omega^{\prime})\frac{\Omega_{0}}{\omega+\Omega_{0}}\alpha(\omega),$
(126)

where $\mathbbmss{P}$ denotes that the principal part is to be taken on
integration and $Y(\omega)$ is a real function, which we determine by
substituting (126) into (118). We find

$\displaystyle
Y(\omega)=\frac{1}{|V(\omega)|^{2}}\left[\frac{2(\omega^{2}-\Omega_{0}^{2})}{\Omega_{0}}-\int_{0}^{\infty}d\omega^{\prime}\left(\frac{\mathbbmss{P}}{\omega-\omega^{\prime}}-\frac{1}{\omega+\omega^{\prime}}\right)|V(\omega^{\prime})|^{2}\right].$

If we substitute our operators, $\hat{B}(\omega)$, expressed in terms of the
function $\alpha(\omega)$ into the commutation relation (118) then we find

$\displaystyle\left[\hat{B}(\omega),\hat{B}^{\dagger}(\omega^{\prime})\right]=\alpha({\omega})\alpha^{*}(\omega^{\prime})\left\\{1-\left(\frac{\omega-\Omega_{0}}{\omega+\Omega_{0}}\right)\left(\frac{\omega^{\prime}-\Omega_{0}}{\omega^{\prime}+\Omega_{0}}\right)\right.$
$\displaystyle\quad+\int_{0}^{\infty}d\omega^{\prime\prime}\left[\left(\frac{\mathbbmss{P}}{\omega-\omega^{\prime\prime}}+Y(\omega)\delta(\omega-\omega^{\prime\prime})\right)\left(\frac{\mathbbmss{P}}{\omega^{\prime}-\omega^{\prime\prime}}+Y(\omega)\delta(\omega^{\prime}-\omega^{\prime\prime})\right)\right.$
$\displaystyle\qquad\qquad\left.\left.-\left(\frac{1}{\omega+\omega^{\prime\prime}}\right)\left(\frac{1}{\omega^{\prime}+\omega^{\prime\prime}}\right)\right]\frac{|V(\omega^{\prime\prime})|^{2}\Omega_{0}^{2}}{(\omega+\Omega_{0})(\omega^{\prime}+\Omega_{0})}\right\\}.$
(128)

Evaluating the integrals and setting the result equal to
$\delta(\omega-\omega^{\prime})$ gives888This requires the use of the
following formula for the product of two principal parts [17]:
$\displaystyle\frac{\mathbbmss{P}}{\omega-\omega^{\prime\prime}}\frac{\mathbbmss{P}}{\omega^{\prime}-\omega^{\prime\prime}}=\frac{\mathbbmss{P}}{\omega^{\prime}-\omega}\left(\frac{\mathbbmss{P}}{\omega-\omega^{\prime\prime}}-\frac{\mathbbmss{P}}{\omega^{\prime}-\omega^{\prime\prime}}\right)+\pi^{2}\delta(\omega-\omega^{\prime\prime})\delta(\omega^{\prime}-\omega^{\prime\prime}).$

$|\alpha(\omega)|^{2}=\frac{(\omega+\Omega_{0})^{2}}{\Omega_{0}^{2}|V(\omega)|^{2}}\left(\frac{1}{Y^{2}(\omega)+\pi^{2}}\right).$
(129)

Note that the diagonalisation does not fix the phase of the complex function
$\alpha(\omega)$ and we are free to choose this as we wish. A convenient
choice is to set

$\alpha(\omega)=\frac{\omega+\Omega_{0}}{\Omega_{0}V^{*}(\omega)}\left(\frac{1}{Y(\omega)-i\pi}\right).$
(130)

##  Appendix C The symmetrically ordered characteristic function

We have made use of the symmetrically ordered characteristic function,

$\chi(\xi)={\rm
Tr}\left[\hat{\rho}\exp\left(\xi\hat{a}^{\dagger}-\xi^{*}\hat{a}\right)\right]\,,$
(131)

to investigate both the dynamics and the steady state of the damped harmonic
oscillator. For completeness, we summarise here the main properties of this
function. Further details of this and also of related characteristic functions
can be found in [17].

The characteristic function is always defined and also well-behaved for any
oscillator state. As $\xi=0$ it reduces to the trace of $\hat{\rho}$ and it
follows that $\chi(0)=1$. More generally, it is the expectation value of the
unitary displacement operator:

$\hat{D}(\xi)=\exp\left(\xi\hat{a}^{\dagger}-\xi^{*}\hat{a}\right)\,.$ (132)

This operator, by virtue of its unitarity, has only eigenvalues of modulus 1
and it follows that

$|\chi(\xi)|\leq 1\,,$ (133)

with the maximum at $\xi=0$.

The most important property is that the density operator and the symmetrically
ordered characteristic function exist in one to one correspondence, analogous
to a Fourier transform pair. Cahill and Glauber [118] (see also [119])
exploited a theorem of Weyl [120] to show that

$\hat{\rho}=\int\frac{d^{\xi}}{\pi}{\rm
Tr}[\hat{\rho}\hat{D}(\xi)]\hat{D}(-\xi).$ (134)

We can extract from the characteristic function the expectation value of any
symmetrically ordered combination of $\hat{a}$ and $\hat{a}^{\dagger}$:

$S\langle\hat{a}^{\dagger
m}\hat{a}^{n}\rangle=\left.\left(\frac{\partial}{\partial\xi}\right)^{m}\left(-\frac{\partial}{\partial\xi^{*}}\right)^{n}\chi(\xi)\right|_{\xi=0}\,.$
(135)

By symmetrically ordered, we mean the average of all possible orderings, for
example:

$\displaystyle S\langle\hat{a}^{\dagger}\hat{a}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left(\hat{a}^{\dagger}\hat{a}+\hat{a}\hat{a}^{\dagger}\right)$
$\displaystyle S\langle\hat{a}^{\dagger 2}\hat{a}^{2}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{6}\left(\hat{a}^{\dagger
2}\hat{a}^{2}+\hat{a}^{\dagger}\hat{a}\hat{a}^{\dagger}\hat{a}+\hat{a}^{\dagger}\hat{a}^{2}\hat{a}^{\dagger}+\hat{a}\hat{a}^{\dagger
2}\hat{a}+\hat{a}\hat{a}^{\dagger}\hat{a}\hat{a}^{\dagger}+\hat{a}^{2}\hat{a}^{\dagger
2}\right)\,.$ (136)

##  Appendix D Thermofields

It is simplest to consider first an isolated discrete oscillator with
annihilation and creation operators $\hat{b}$ and $\hat{b}^{\dagger}$. For
such an oscillator in a thermal state at temperature $T$ the density operator
has the simple, diagonal form:

$\hat{\rho}_{T}=(1-e^{-\beta\hbar\omega})\sum_{n=0}^{\infty}e^{-n\beta\hbar\omega}|n\rangle\langle
n|\,,$ (137)

where $\beta=(k_{B}T)^{-1}$ is the inverse temperature. The mean number of
excitations is

$\bar{n}=\frac{1}{e^{\beta\hbar\omega}-1}$ (138)

and we can write the density operator in terms of this mean:

$\hat{\rho}_{T}=\frac{1}{\bar{n}+1}\sum_{n=0}^{\infty}\left(\frac{\bar{n}}{\bar{n}+1}\right)^{n}|n\rangle\langle
n|\,.$ (139)

The thermofield technique [53, 54, 55, 56, 57, 58] starts with the observation
that we can write a _pure_ state that has the same statistical properties as
the thermal mixed state (137). To construct this state we consider a doubled
state space in which we introduce a second oscillator with annihilation and
creation operators $\hat{\tilde{b}}$ and $\hat{\tilde{b}}^{\dagger}$. The two-
mode pure state, the _thermal vacuum_ :

$|0(\beta)\rangle=(1-e^{-\beta\hbar\omega})^{1/2}\sum_{n=0}^{\infty}e^{-\beta\hbar\omega/2}|n,\tilde{n}\rangle$
(140)

has precisely the same single-mode properties as the single-mode thermal
state:

$\langle 0(\beta)|f(\hat{b},\hat{b}^{\dagger})|0(\beta)\rangle={\rm
Tr}\left(f(\hat{b},\hat{b}^{\dagger})\right)\,.$ (141)

It is straightforward to show that a similar procedure can be applied can be
applied to express any mixed state in terms of a pure state in a doubled state
space [56]. When this procedure was rediscovered in quantum information
theory, it acquired the name purification [106].

The benefit of introducing the thermal vacuum state comes from the fact that
it is related to the two-mode vacuum state, $|0,\tilde{0}\rangle$, via a
unitary transformation:

$\displaystyle|0(\beta)\rangle$ $\displaystyle=$
$\displaystyle\hat{S}(\theta)|0,\tilde{0}\rangle$ (142) $\displaystyle=$
$\displaystyle\exp\left[\theta(\beta)\left(\hat{\tilde{b}}^{\dagger}\hat{b}^{\dagger}-\hat{b}\hat{\tilde{b}}\right)\right]|0,\tilde{0}\rangle\,.$

This transformation produces the desired state if we select $\theta(\beta)$
such that

${\rm sinh}^{2}\theta(\beta)=\bar{n}\,.$ (143)

Readers with a background in quantum optics may recognise $|0(\beta)\rangle$
as a two-mode squeezed vacuum state [17, 121]. The unitary nature of this
transformation means that we can convert, by means of the inverse
transformation, our effective thermal state into a vacuum state, accompanied
by a modified Hamiltonian. Before we can do this, however, we require a
Hamiltonian for the tilde oscillator. The natural way to introduce this is as
an inverted oscillator, so that our free oscillator Hamiltonian becomes

$\hat{H}_{0}=\hbar\Omega_{0}\left(\hat{b}^{\dagger}\hat{b}-\hat{\tilde{b}}^{\dagger}\hat{\tilde{b}}\right)\,,$
(144)

which has the advantage that it commutes with the unitary transformation. It
is essential, moreover, to avoid introducing undesired couplings between the
original and the tilde operators.

If our oscillator is coupled to another quantum system via its annihilation
and creation operators, then the required unitary transformation effects the
replacement

$\displaystyle\hat{b}$ $\displaystyle\rightarrow$
$\displaystyle\hat{S}(\theta)\hat{b}\hat{S}^{\dagger}(\theta)=\hat{b}\>{\rm
cosh}\theta(\beta)-\hat{\tilde{b}}^{\dagger}{\rm sinh}\theta(\beta)$
$\displaystyle\hat{b}^{\dagger}$ $\displaystyle\rightarrow$
$\displaystyle\hat{S}(\theta)\hat{b}^{\dagger}\hat{S}^{\dagger}(\theta)=\hat{b}^{\dagger}{\rm
cosh}\theta(\beta)-\hat{\tilde{b}}\>{\rm sinh}\theta(\beta)\,.$ (145)

As a simple illustration we note that the expectation value of
$\hat{b}^{\dagger}\hat{b}$ is

$\displaystyle\langle
0(\beta)|\hat{b}^{\dagger}\hat{b}|0(\beta)\rangle\rangle$ $\displaystyle=$
$\displaystyle\langle 0,\tilde{0}|(\hat{b}^{\dagger}{\rm
cosh}\theta(\beta)-\hat{\tilde{b}}\>{\rm sinh}\theta(\beta))$ (146)
$\displaystyle\qquad\times(\hat{b}\>{\rm
cosh}\theta(\beta)-\hat{\tilde{b}}^{\dagger}{\rm
sinh}\theta(\beta))|0,\tilde{0}\rangle$ $\displaystyle=$ $\displaystyle{\rm
sinh}^{2}\theta(\beta)\,.$

In place of a coupling to a single harmonic oscillator in a thermal state with
inverse temperature $\beta$, the transformed Hamiltonian has a coupling to a
regular oscillator in its ground state with a coupling strength increased by
${\rm cosh}\theta(\beta)$ and a coupling to a second _inverted_ oscillator in
its most highly excited state, with the original coupling multiplied by ${\rm
sinh}\theta(\beta)$. The inverted oscillator can only inject quanta (at least
initially) while the regular oscillator can only extract them.

To complete the picture we need only generalise this description to our
continuum operators. We do this by using our continuum thermal vacuum state in
the form

$\displaystyle|0(\beta)\rangle$ $\displaystyle=$
$\displaystyle\hat{S}(\theta[\omega])|0,\tilde{0}\rangle$ (147)
$\displaystyle=$ $\displaystyle\exp\left[\int
d\omega\>\theta(\beta,\omega)\left(\hat{\tilde{b}}^{\dagger}(\omega)\hat{b}^{\dagger}(\omega)-\hat{b}(\omega)\hat{\tilde{b}}(\omega)\right)\right]|0,\tilde{0}\rangle\,,$

where $|0,\tilde{0}\rangle$ now denotes the doubled continuum vacuum state and
$\hat{\tilde{b}}(\omega)$ is the annihilation operator corresponding to
_adding_ a quantum of frequency $\omega$ to the inverted, tilde reservoir. As
with the discrete oscillator, we can transform into an equivalent vacuum
picture with the free Hamiltonian for the reservoir changed to

$\hat{H}_{0}=\int
d\omega\>\hbar\omega\left(\hat{b}^{\dagger}(\omega)\hat{b}(\omega)-\hat{\tilde{b}}^{\dagger}(\omega)\hat{\tilde{b}}(\omega)\right)\,,$
(148)

and the continuum annihilation and creation operators transformed by the
inverse unitary transformation

$\displaystyle\hat{b}(\omega^{\prime})$ $\displaystyle\rightarrow$
$\displaystyle\hat{S}(\theta[\omega])\hat{b}(\omega^{\prime})\hat{S}^{\dagger}(\theta[\omega])=\hat{b}(\omega^{\prime}){\rm
cosh}\theta(\beta,\omega^{\prime})-\hat{\tilde{b}}^{\dagger}(\omega^{\prime}){\rm
sinh}\theta(\beta,\omega^{\prime})$
$\displaystyle\hat{b}^{\dagger}(\omega^{\prime})$ $\displaystyle\rightarrow$
$\displaystyle\hat{S}(\theta[\omega])\hat{b}^{\dagger}(\omega^{\prime})\hat{S}^{\dagger}(\theta[\omega])=\hat{b}^{\dagger}(\omega^{\prime}){\rm
cosh}\theta(\beta,\omega^{\prime})-\hat{\tilde{b}}(\omega^{\prime}){\rm
sinh}\theta(\beta,\omega^{\prime})\,,$

so that the expectation value of
$\hat{b}^{\dagger}(\omega)\hat{b}(\omega^{\prime})$ is

$\displaystyle\langle
0(\beta)|\hat{b}^{\dagger}(\omega)\hat{b}(\omega^{\prime})|0(\beta)\rangle$
$\displaystyle=$ $\displaystyle\langle
0,\tilde{0}|(\hat{b}^{\dagger}(\omega){\rm
cosh}\theta(\beta,\omega)-\hat{\tilde{b}}(\omega){\rm
sinh}\theta(\beta,\omega))$ (150)
$\displaystyle\qquad\times(\hat{b}(\omega^{\prime}){\rm
cosh}\theta(\beta,\omega^{\prime})-\hat{\tilde{b}}^{\dagger}(\omega^{\prime}){\rm
sinh}\theta(\beta,\omega^{\prime}))|0,\tilde{0}\rangle$ $\displaystyle=$
$\displaystyle\bar{n}(\omega)\delta(\omega-\omega^{\prime})\,.$

## References

## References

  * [1] Chan J, Alegre T P M, Safavi-Naeini A H, Hill J T, Gröblacher S, Aspelmeyer M and Painter O 2011 Laser Cooling of a Nanomechanical Oscillator into its Quantum Ground State _Nature_ 478 89–92
  * [2] Gigan S, Böhm H R, Paternostro M, Blaser F, Langer G, Hertzberg J B, Schwab K C, Bäuerle D, Aspelmeyer M and Zeilinger A 2006 Self-Cooling of a Micromirror by Radiation Pressure _Nature_ 444 67–70
  * [3] Stannigel K, Rabl P, Sørensen AS, Zoller P and Lukin M D 2010 Optomechanical Transducers for Long-Distance Quantum Communication _Phys. Rev. Lett._ 105 220501
  * [4] Galve F, Pachón L A and Zueco P 2010 Bringing Entanglement to the High Temperature Limit _Phys. Rev. Lett._ 105 180501
  * [5] Zhang K, Bariani F and Meystre P 2014 Quantum Optomechanical Heat Engine _Phys. Rev. Lett._ 112 150602
  * [6] Joshi C, Öhberg P, Cresser J D and Andersson E 2014 Markovian Evolution of Strongly Coupled Harmonic Oscillators _Phys. Rev._ A 90 063815
  * [7] Brunelli M, Xuereb A, Ferraro A, De Chiara G, Kiesel N and Paternostro M 2015 Out-of-Equilibrium Thermodynamics of Quantum Optomechanical Systems _New J. Phys._ 17 035016
  * [8] Campisi M, Talkner P and Hänggi P 2009 Fluctuation Theorem for Arbitrary Open Quantum Systems _Phys. Rev. Lett._ 102 210401
  * [9] Talkner P, Campisi M and Hänggi P 2009 Fluctuation Theorems in Driven Open Quantum Systems _J. Stat. Mech._ P02025
  * [10] Gemmer J, Michel M and Mahler G 2004 _Quantum Thermodynamics_ (Berlin: Springer)
  * [11] Vinjanampathy S and Anders J 2016 Quantum Thermodynamics _Contemp. Phys._ 57 545–579
  * [12] Allahverdyan AE and Nieuwenhuizen TM 2000 Extraction of Work from a Single Thermal Bath in the Quantum Regime _Phys. Rev. Lett._ 85 1799–1802
  * [13] Caldeira A O and Leggett A J 1983 Path Integral Approach to Quantum Brownian Motion _Physica_ 121 A 587–616
  * [14] Caldeira A O and Leggett A J 1983 Quantum Tunnelling in a Dissipative System _Ann. Phys. N.Y._ 149 374–456
  * [15] Senitzky I R 1960 Dissipation in Quantum Mechanics. The Harmonic Oscillator _Phys. Rev._ 119 670–679
  * [16] Louisell W H 1964 _Radiation and Noise in Quantum Electronics_ (New York: McGraw-Hill)
  * [17] Barnett S M and Radmore P M 1997 _Methods in Theoretical Quantum Optics_ (Oxford: Oxford University Press)
  * [18] Louisell W H 1973 _Quantum Statistical Properties of Radiation_ (New York: Wiley)
  * [19] Grabert H 1982 _Projection Operator Techniques in Nonequilibrium Statistical Mechanics_ (Berlin: Springer-Verlag)
  * [20] Peřina J 1991 _Quantum Statistics of Linear and Nonlinear Optical Phenomena_ 2nd ed. (Dordrecht: Kluwer)
  * [21] Meystre P and Sargent III M 1991 _Elements of Quantum Optics_ 2nd ed. (Berlin: Springer-Verlag)
  * [22] Carmichael H 1993 _An Open Systems Approach to Quantum Optics_ (Berlin: Springer-Verlag)
  * [23] Walls D F and Milburn G J 1994 _Quantum Optics_ (Berlin: Springer-Verlag)
  * [24] Carmichael H J 1999 _Statistical Methods in Quantum Optics 1_ (Berlin: Springer)
  * [25] Breuer H-P and Petruccione F 2002 _The Theory of Open Quantum Systems_ (Oxford; Oxford University Press)
  * [26] Gardiner C W and Zoller P 2004 _Quantum Noise_ 2nd ed. (Berlin: Springer)
  * [27] Ficek Z and Swain S 2004 _Quantum Interference and Coherence_ (New York: Springer)
  * [28] Wiseman H M and Milburn G J 2010 _Quantum Measurement and Control_ (Cambridge: Cambridge University Press)
  * [29] Weiss U 2012 _Quantum Dissipative Systems_ 4th ed. (New Jersey: World Scientific)
  * [30] Navotný J, Alber G and Jex I 2012 Asymptotic Properties of Quantum Markov Chains _J. Phys. A: Math. Theo._ 45 485301
  * [31] Navotný J, Maryška J and Jex I 2018 Quantum Markov Processes: from Attractors to Explicit Forms of Asymptotic States _Eur. Phys. J. Plus_ 133 310
  * [32] Barnett S M and Stenholm S 2001 Hazards of Reservoir Memory _Phys. Rev._ A 64 033808
  * [33] Maniscalco S 2007 Complete Positivity of a Spin-1/2 Master Equation with Memory _Phys. Rev._ A 75 062103
  * [34] Kaikka P, Cresser J D and Maniscalco S 2011 Comparing Different Non-Markovianity Measures in a Driven Qubit System _Phys. Rev._ A 83 012112
  * [35] Breuer H-P 2012 Foundations and Measures of Non-Markovianity _J. Phys. B: At. Mol. Opt. Phys._ 45 154001
  * [36] Pernice A, Helm J and Strunz W T 2012 System-Environment Correlations and Non-Markovian Dynamics _J. Phys. B: At. Mol. Opt. Phys._ 45 154007
  * [37] Vacchini B 2012 A Classical Appraisal of Quantum Definitions of Non-Markovian Dynamics _J. Phys. B: At. Mol. Opt. Phys._ 45 154007
  * [38] Chraśinński D and Maniscalco S 2014 Degree of Non-Markovianity of Quantum Evolution _Phys. Rev. Lett._ 112 120404
  * [39] Hall M J W, Cresser J D, Li L and Andersson E 2014 Canonical Form of Master Equations and Characterization of Non-Markovianity _Phys. Rev._ A 89 042120
  * [40] Hu B L, Paz J P and Zhang Y 1992 Quantum Brownian Motion in a General Environment: Exact Master Equation with Nonlocal Dissipation and Colored Noise _Phys. Rev._ D 45 2843–2861
  * [41] Diosi L 1993 Caldeira-Leggett Master Equation and Medium Temperatures _Physica_ 199 517–526
  * [42] Grabert H and Weiss U 1984 Quantum Theory of the Damped Harmonic Oscillator _Z. Phys. B - Cond. Matt._ 55 87–94
  * [43] Ford G W, Lewis J T and O’Connell R F 1985 Quantum Oscillator in a Blackbody Radiation Field _Phys. Rev. Lett._ 55 2273–2276
  * [44] Ford G W, Lewis J T and O’Connell R F 1985 Quantum Langevin Equation _Phys. Rev._ A 37 4419–4428
  * [45] Feynman R and Hibbs A R 2010 _Quantum Mechanics and Path Integrals_ (Emended ed.) (New York: Dover)
  * [46] Zinn-Justin J 2002 _Quantum Field Theory and Critical Phenomena_ (Oxford: Oxford University Press)
  * [47] Kleinert H 2009 _Path Integrals in Quantum Mechanics, Statistics, Polymer physics and Financial markets_ (fifth ed.) (New Jersey: World Scientific)
  * [48] Smith C M and Caldeira A O 1987 Generalized Feynman-Vernon Approach to Dissipative Quantum Systems _Phys. Rev._ A 36 3509–3511
  * [49] Grabert H, Schramm P and Ingold G-L 1988 Quantum Brownian Motion: the Functional Integral Approach _Phys. Rep._ 168 115–207
  * [50] Abrikosov A A, Gorkov L P and Dzaloshinski I E 1975 _Methods of Quantum Field Theory in Statistical Physics_ (New York: Dover)
  * [51] Lifshitz E M and Pitaevskii L P 1980 _Statistical Physics part 2_ (Oxford: Butterworth Heinemann)
  * [52] Le Bellac M 1996 _Thermal Field Theory_ (Cambridge: Cambridge University Press)
  * [53] Takahashi Y and Umezawa H 1975 Thermo Field Dynamics _Collect. Phenom._ 2 55–80. Reprinted in 1996 _Int. J. Theo. Phys._ 10 1755–1805
  * [54] Umezawa H, Matsumoto H and Tachiki M 1982 _Thermo Field Dynamics and Condensed States_ (Amsterdam: North-Holland)
  * [55] Barnett S M and Knight P L 1985 Thermofield Analysis of Squeezing and Statistical Mixtures in Quantum Optics _J. Opt. Soc. Am._ B 2 467–479
  * [56] Barnett S M and Dalton B J 1987 Liouville Space Description of Thermofields and their Generalisations _J. Phys. A: Math. Gen._ 20 411–418
  * [57] Umezawa H 1993 _Advanced Field Theory, Micro, Macro and Thermal Physics_ (New York: American Institute of Physics)
  * [58] Dalton B J, Jeffers J and Barnett S M 2015 _Phase Space Methods for Degenerate Quantum Gases_ (Oxford: Oxford University Press)
  * [59] Subaşı Y, Fleming C H , Taylor J M, and Hu B L 2012 Equilibrium states of open quantum systems in the strong coupling regime _Phys. Rev._ E 86 061132
  * [60] Haake F and Reibold R 1985 Strong Damping and Low-Temperature Anomalies for the Harmonic Oscillator _Phys. Rev._ A 32 2462–2475
  * [61] Ford G W, Lewis J T and O’Connell R F 1988 Independent Oscillator Model of a Heat Bath: Exact Diagonalization of the Hamiltonian _J. Stat. Phys._ 53 439–455
  * [62] Barnett S M, Cresser J D and Croke S 2015 Theory of the Strongly-Damped Quantum Harmonic Oscillator _arXiv_ 1508:02442 [quant-ph]
  * [63] Huttner B and Barnett S M 1992 Dispersion and Loss in a Hopfield Dielectric _Europhys. Lett._ 18 487–492
  * [64] Huttner B and Barnett S M 1992 Quantization of the Electromagnetic Field in Dielectrics _Phys. Rev._ A 46 4306–4322
  * [65] Dutra S M 2005 Cavity Quantum Electrodynamics (Hoboken New Jersey: Wiley)
  * [66] Philbin T G and Anders J 2016 Thermal Energies and Quantum Damped Oscillators Coupled to Reservoirs _J. Phys. A: Math. Theo._ 49 215303
  * [67] Philbin T G 2012 Quantum Dynamics of the Damped Harmonic Oscillator _New J. Phys._ 14 083043
  * [68] Schrödinger E 1926 Quantisierung als Eigenwertproblem _Ann. Phys._ 79, 489-527; reprinted in translation as “Quantisation as a Problem of Proper Values (Part II)” in Schrödinger E 1982 _Collected Papers on Wave Mechanics_ 3rd ed. (Providence, Rhode Island: AMS Chelsea Publishing)
  * [69] Langevin P 1908 Sur la Théorie du Mouvement Brownien _Comptes Rend. Acad. Sci. (Paris)_ 146 530–533; reprinted in a translation by Gythiel A 1997 “On the Theory of Brownian Motion” _Am. J. Phys._ 65 1079–1081
  * [70] Uhlenbeck G E and Ornstein L S 1930 On the Theory of Brownian motion _Phys. Rev._ 36, 823–841
  * [71] Wang M C and Uhlenbeck G E 1945 On the Theory of Brownian motion II _Rev. Mod. Phys._ 17 323–342
  * [72] Risken H 1989 _The Fokker-Planck Equation_ 2nd ed. (Berlin: Springer)
  * [73] Lemons D S 2002 _An Introduction to Stochastic Processes in Physics_ (Baltimore: Johns Hopkins University Press)
  * [74] Mazo R M 2002 _Brownian Motion_ (Oxford: Oxford University Press)
  * [75] Agarwal G S 1971 Brownian Motion of a Quantum Oscillator _Phys. Rev._ A 4 739–747
  * [76] Ullersma P 1966 An Exactly Solvable Model for Brownian Motion I. Derivation of the Langevin Equation _Physica_ 32 27-55
  * [77] Ullersma P 1966 An Exactly Solvable Model for Brownian Motion II. Derivation of the Fokker-Planck equation and the Master Equation _Physica_ 32 56–73
  * [78] Kheirandish F and Amooshahi 2006 Electromagnetic Field Quantisation in a Linear Polarizable and Magnetizable Medium _Phys. Rev._ A 74 042102
  * [79] Kheirandish F and Soltani M 2008 Extension of the Huttner-Barnett Model to a Magnetodielectric Medium _Phys. Rev._ A 78 012102
  * [80] Kheirandish F, Amooshahi M and Soltani M 2009 A Canonical Approach to Electromagnetic Field Quantization in a Nonhomogeneous and Anistropic Magnetodielectric Medium _J. Phys. B: At. Mol. Opt._ 42 075504
  * [81] Amooshahi M 2009 Canonical Quantisation of the Electromagnetic Field in an Anisotropic Polarizable and Magnetizable Medium _J. Math. Phys._ 50 062301
  * [82] Munro W J and Gardiner C W 1996 Non-Rotating-Wave Master Equation _Phys. Rev._ A 53 2633–2640
  * [83] Stenholm S 1996 Simple Quantum Dynamics in _Quantum Dynamics of Simple Systems_ G. L. Oppo, S. M. Barnett, E. Riis and M. Wilkinson eds. (Bristol: Scottish Universities Summer Schools in Physics and Institute of Physics) pp. 267–313
  * [84] Barnett S M and Cresser J D 2005 Quantum Theory of Friction _Phys. Rev._ A 72 022107
  * [85] Barnett S M, Jeffers J and Cresser J D 2006 From Measurements to Quantum Friction _J. Phys: Condens. Matter_ 18 S401–S410
  * [86] Bixon M and Jortner J 1968 Intramolecular Radiationless Transitions _J. Chem. Phys._ 48 715–725
  * [87] Kyröla E 1984 Photoexcitation of a Quasi-Continuum: Connections to Few-Level Dynamics _J. Opt. Soc. Am._ B 1 737–74
  * [88] Milonni P W, Ackerhalt J R, Galbraith H W and Shil M-L 1988 Exponential Decay and Quantum Mechanical Spreading in a Quasicontinuum Model _Phys. Rev._ A 28 32–39
  * [89] Eberly J H, Yeh J J and Bowden C M 1982 Interupted Coarse-Grained Theory of Quasi-Continuum Photoexcitation _Chem. Phys. Lett._ 86 76–80
  * [90] Radmore P M and Tarzi S 1987 Photo-Excitation of a Structured Continuum: a Soluble Model _J. Mod. Opt._ 34 1409–1420
  * [91] Bohr H 2018 _Almost Periodic Functions_ (New York: Dover)
  * [92] Fano U 1935 Sullo Spettro di Assortbimento dei Gas Nobili Presso il Limite dello Spettro d’Arco _Il Nuovo Cimento_ 12 154–161
  * [93] Barnett SM and Radmore PM 1988 Quantum Theory of Cavity Quasimodes _Opt. Commun._ 68 364–368
  * [94] Georgievskii Y and Pollak E 1995 Activated Rate Processes: Anharmonic Corrections to the Quantum Rate _J. Chem. Phys._ 103 8910–8920
  * [95] Ratchov A, Faure F and Hekking F WJ 2005 Loss of Quantum Coherence in a System Coupled to a Zero-Temperature Environment _Eur. J. Phys._ B 46 519–528
  * [96] Ghosh A, Bandyopadhyay M, Dattagupta S and Gupta S 2023 Quantum Brownian Motion: A Review arXiv:2306.02665
  * [97] Compagno C, Passante R and Persico F 1995 _Atom-Field Interactions and Dressed Atoms_ (Cambridge: Cambridge University Press)
  * [98] Casimir H B G and Polder D 1948 The Influence of Retardation on the London-van der Waals Forces _Phys. Rev._ 73 360–372
  * [99] Power E A 1964 _Introductory Quantum Electrodynamics_ (London: Longmans)
  * [100] Craig D P and Thirunamachandran T 1984 _Molecular Quantum Electrodynamics_ (London: Academic Press)
  * [101] Loudon R and Barnett S M 2006 Theory of the Linear Polarizability of a Two-Level Atom _J. Phys. B: At. Mol. Opt. Phys._ 39 S555–S563
  * [102] Milonni PW, Loudon R, Berman PR and Barnett SM 2008 Linear Polarizabilities of Two- and Three- Level Atoms _Phys. Rev._ A 77 043835
  * [103] Philbin T G and Horsley S A R 2013 Damping the Zero-Point Energy of a Harmonic Oscillator _arXiv_ : 1304.0977v2 [quant-ph]
  * [104] Araki H and Lieb E H 1970 Entropy iIequalities _Commun. Math. Phys._ 18 160–170
  * [105] Wehrl A 1978 General Properties of Entropy _Rev. Mod. Phys._ 50 221–260
  * [106] Barnett S M 2009 _Quantum Information_ (Oxford: Oxford University Press)
  * [107] Barnett S M and Phoenix S J D 1989 Entropy as a Measure of Quantum Optical Correlation _Phys. Rev._ A 40 2404–2409
  * [108] Kirkwood J G 1935, Statistical mechanics of fluid mixtures, J. Chem. Phys. 3, 300
  * [109] Trushechkin A S, Merkli M, Cresser J D, and Anders J 2022 AVS Quantum Science 4, 012301
  * [110] Martinez E A and Paz J P 2013 Dynamics and Thermodynamics of Linear Quantum Open Systems _Phys. Rev. Lett._ 110 130406
  * [111] Vaccaro J A and Barnett S M 2011 Information Erasure without an Energy Cost _Proc. R. Soc._ A 467 1770–1778
  * [112] Barnett S M and Vaccaro J A 2013 Beyond Landauer Erasure _Entropy_ 15 4956–4968
  * [113] Roncaglia A J, Cerisola F and Paz J P 2014 Work as a Generalized Measurement _Phys. Rev. Lett._ 113 0250601
  * [114] Binder F, Vinjanampathy S, Modi K and Goold J 2015 Operational Thermodynamics of Open Quantum Systems _Phys. Rev._ E 91 032119
  * [115] Brandao F, Horodecki M, Ng K, Oppenheim J and Wehner S 2015 The Second Laws of Quantum Thermodynamics _Proc. Natl. Acad. Sci. USA_ 112 3275–3279
  * [116] Fano U 1956 Effects of Configuration Interaction on Intensities and Phase Shifts _Phys. Rev._ 124 1866–1878
  * [117] Dirac P A M 1927 Über die Quantenmechanik der Stoßvorgänge _Z. Phys._ 44 585–595; reprinted in translation as “On the Quantum Mechanics of Collisions” in R. H. Dalitz (ed.) 1995 _The Collected Works of P. A. M. Dirac 1924–1948_ (Cambridge: Cambridge University Press)
  * [118] Cahill K E and Glauber R J 1969 Ordered Expansions of Boson Amplitude Operators _Phys. Rev._ 177 1857–1881
  * [119] Barnett S M and Dalton B J 2023 Glauber-Sudarshan P-Representations for Fermions _Phys. Scr._ 98 044006
  * [120] Weyl H 1950 _The Theory of Groups and Quantum Mechanics_ (New York: Dover) p. 274
  * [121] Milburn G J 1984 Multimode Minimum Uncertainty Squeezed States _J. Phys. A: Math. Gen._ 17 737–745