Revisiting the damped quantum harmonic oscillator

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.


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.

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, γ, 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 ẍ + γ ẋ + ω 2 0 x = 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 γ ω 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.

Hamiltonian for the strongly-damped harmonic oscillator
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] If we complete the square we can rewrite this in a minimal-coupling form to arrive at the alternative form: where Each term in our Hamiltonian (3) is strictly positive only if this quantity is positive, corresponding to a real frequency ω 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 while ω 0 can take any positive value.‡ ‡ In 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 ω 0 and Ω 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, Ω 0 and ω 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.

Heisenberg-Langevin equation: a tale of two frequencies
We find that both the frequencies ω 0 and Ω 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 Appendix A. We find where κ(t) is the memory kernel: and F (t) is the Langevin force: 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: and note that we can also write this in the form The first of these equations is written in terms of the frequency Ω 0 and the second in terms of ω 0 .It remains to determine the physical role of each of these, which we can do by considering very short and longer time scales.

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 For a very short time, δt, this becomes The integral term in Eq. ( 11) is of order δt 2 by virtue of Eq. ( 7) as for short times κ is of order δt.The combination of Eq. ( 12) and ẋ = p/m, which is always true, leads to the ultra-short time behaviour This is the short-time behaviour of a harmonic oscillator of frequency Ω 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 κ(t) involves a summation of oscillating cosines of different frequencies, one for each environment oscillator, and that these will dephase, causing κ(t) to decay.We define the longer-time regime as that for which we can approximate κ(t) ≈ 0. In this regime, our our equation of motion (10) becomes 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 ω 0 and not Ω 0 .Thus the question of whether our Hamiltonian applies in the over-damped regime depends on the relationship between the damping and ω 0 and not Ω 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 Ω 0 cannot be reduced.However there is no such bound on ω 0 and, indeed, we can set ω 0 = 0 + without invalidating our model.As it is ω 0 and not Ω 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 with 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 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 κ(t) x(0) in (9).This means times for which t < ω −1 c , 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 κ(t) ≈ 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 ẋ does not vary significantly on the timescale ω −1 0 , then we can make the approximation It then follows that our equation for the expectation value of the position becomes (in this Markovian, longer-time regime) which is the familiar equation for a damped harmonic oscillator.Note however, that it is the frequency ω 0 and not Ω 0 that appears.We enter the over-damped regime when ω 0 < γ/2.The constraints on the natural frequencies in (4) correspond in this Ohmic damping example to We see that there is no lower bound on ω 0 (although it must be real and greater than zero) but there is a bound on Ω 0 : it must exceed the geometric mean of the damping rate γ and the cut-off frequency ω c multiplied by 2/π.Thus we see again that we cannot allow the cut-off frequency to tend to infinity.Similarly the frequency Ω 0 will lie somewhere between γ and the rather larger ω c so we cannot have Ω 0 less than γ.What saves the model in the strongly damped regime is the fact that it is ω 0 rather than Ω 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.

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: In terms of these operators our Hamiltonian, Eq. ( 2), becomes when unimportant constant shifts in the ground-state energies are removed and When written in terms of this quantity, our positivity condition (5) becomes 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] §.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, b(ω) and b † (ω), satisfying the commutation relations b and our Hamiltonian becomes [64] and the positivity condition is then 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.§ As 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.

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, B(ω) and their conjugates B † (ω) that satisfy the operator equations for all positive frequencies ω.These eigenoperators are complete if, in addition to these they also satisfy the condition 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: 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 Appendix 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 B(ω) and B † (ω) operators and then use the commutation relations to extract the coefficients in this expansion.For the oscillator operators we find The requirement that these operators satisfy the familiar boson commutation relation, â, â † = 1, provides a constraint on the functions α(ω) and β(ω) in the form where we have used (B.10).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: 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].

Physical constraints
We leave until the next section a discussion of the physical interpretation of the probability density π(ω) 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 we note that |α(ω)| 2 is finite for all ω and it follows, therefore, from (33) that the average value ω −1 is finite.It follows from the eigenoperator equations ( 28) that we can write the Hamiltonian in the form 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 â2 and â †2 should vanish so that where we have again used (B.10).This implies that 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 Finally we can apply the Cauchy-Schwartz inequality to provide a lower bound on the value of ω −1 : These inequalities are useful in determining the properties of the ground state.

Ground-state
It is immediately clear from the form of the Hamiltonian (26) that the ground state of the oscillator is not that of the undamped oscillator, that is the state annihilated by â.
To see this we need only note that the interaction term â † b † (ω) 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 , is the zero-eigenvalue eigenstate of all the annihilation operators B(ω): 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: The most straightforward way to determine the form of this mixed state is to use the characteristic function [17]: 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 Appendix C. To evaluate the characteristic function it we express â and â † in terms of the eigenoperators This simple form, Gaussian in ξ, is characteristic of a squeezed thermal state.We can rewrite this characteristic function in terms of our probability density in the form We make use of the operator ordering theorem [17] exp which holds if the two operators Â and B both commute with their commutator [ Â, B].
where ξ r and ξ i are the real and imaginary parts of ξ respectively.We note, in passing, that the quadrature operators for the oscillator (familiar from quantum optics [17]) have unequal uncertainties: The product of these two exceeds 1 2 , as it must, by virtue of the Cauchy-Schwartz inequality (39).

Ground-state energy
The ground-state energy of the undamped harmonic oscillator is simply 1  2 hω 0 , where ω 0 is the natural angular frequency of the oscillator.For an underdamped oscillator, the oscillation frequency is ω 2 0 − γ 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 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 ω 0 or Ω 0 , makes no explicit appearance in these expectation values, but the averages ω and ω −1 will depend on these frequencies.We have seen that the quantum damped harmonic oscillator is characterised by two natural frequencies, ω 0 and Ω 0 , and it is natural to define the energy of the oscillator in terms of one or other of these: where f 0 = ω 0 or Ω 0 .We shall consider a third possible natural frequency, ω diag , and the associated energy below.The ground-state expectation values, Eq. ( 46), mean that we can write the expectation value of the energy in the form Recall that the Cauchy-Schwartz inequality requires that and it follows that the ground-state energy is bounded by The global minimum of this expression occurs for ω = f 0 and hence for any choice of frequency in the expression for the oscillator energy: 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.

Diagonal form of the oscillator ground state
We can diagonalise the density operator for the oscillator alone, ρ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, To complete the diagonalisation we need only choose ω diag such that the expectation values of ĉ2 and ĉ †2 are zero: This frequency is the geometric mean of the two frequencies ω and ω −1 −1 and it is, by virtue of ( 38) and (39), less than Ω 0 .The mean number of c-quanta in the oscillator ground-state is which exceeds 0, as it should, by virtue of (39).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 We may interpret this state as a thermal state for the oscillator at the shifted frequency ω diag and at an effective "temperature" We should note, however, that the true temperature in the ground-state is zero and that this quantity and the frequency ω 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 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: and that the quantum mutual information [106], or index of correlation [107], between the oscillator and the environment is

Physical meaning of π(ω)
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 π(ω).
Here we present the case for identifying this probability density with the contribution from the dressed modes, associated with the eigenoperators B(ω), 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 (46).We know that the ground state of a harmonic oscillator of frequency ω has 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 π(ω) then the resulting average mean-square values will be which correspond to those obtained above.Note that the requirement that x2 must be finite imposes the condition that at zero frequency Our second point arises from the form of the Hamiltonian for the oscillator We can, by virtue of (37), write this as a combination of potentials corresponding to different frequencies but weighted by π(ω): Finally we note that the effective mean energy of the oscillator, which is associated with the diagonal form of the density operator (56) which neatly combines the characteristic ground state energies of the dressed oscillators, weighted by the probability distribution π(ω).We note that this mean energy is less than 1 2 hΩ 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 x2 and p2 , the form of the oscillator Hamiltonian and the effective mean energy) leads us to interpret π(ω) 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 π(ω) but we find it helpful to do so.

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 Appendix 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 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, 0 , by a unitary transformation: 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: As with the zero temperature ground state, we can determine the form of this using the characteristic function: We can transform this into a vacuum expectation value by applying a unitary transformation to the annihilation and creation operators B(ω) and B † (ω): Applying this transformation to our characteristic function replaces (69) by an equivalent vacuum expectation value.Evaluating this gives As with the zero-temperature steady state, this is a simple Gaussian in ξ and, again, is characteristic of a squeezed thermal state.When expressed in terms of our probability density, π(ω), we find: We note that this has the same general form as the characteristic function for the ground state, Eq (44), but with the probability density π(ω) replaced by a thermallyweighted density π(ω)coth(βhω/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 We have seen that the oscillator is characterised by (at least) two different natural frequencies, Ω 0 and ω 0 .In terms of these, the mean energy of the oscillator alone is where f 0 = ω 0 or Ω 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 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 It follows that the high temperature limits of the kinetic and potential energies are, respectively: 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 ω 2 ω −2 ≥ 1 and, as ω 2 = Ω 2 0 , it follows that Ω 2 0 ω −2 ≥ 1, so that the potential energy, when expressed in terms of Ω 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 ω 0 .We shall see below that in the limit of weak damping, when Ω 0 and ω 0 tend to a common value, the probability distribution π(ω) becomes sharply peaked around ω = Ω 0 so that the equipartition of energy for the potential energy is restored.

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, B(ω) and B † (ω) 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: In particular, we can determine the time-evolution of the annihilation operator for the oscillator in this way: 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 (D.11).As the environment is in a stationary state, so that b(ω, 0) = 0 = b † (ω, 0) , the expectation values of the position and momentum operators take a pleasingly simple form: where the double angle brackets denote, as before, averages over our probability distribution π(ω) 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 π(ω) as a frequency probability distribution for the damped oscillator, as they may be viewed as the evolution of an undamped oscillator with a frequency ω 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 (79), has the necessary short-time form of that for an undamped oscillator where we have used the identity ω 2 = Ω 2 0 .The effects of the coupling to the environment enter at order δ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 (79) includes the possibilities of both criticallydamped 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 cos(ωt) 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 ω sin(ωt) , 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.

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 → 0 or β → ∞.With a little effort we find (using the method of characteristics [17]) where 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 ξ(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 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 so the steady-state characteristic function is 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.

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 π(ω).Determining this, together with the initial conditions, provides all the information required.We can calculate π(ω) directly from the frequency-dependence of the coupling between the oscillator and its environment or, more simply, select a form for π(ω) and proceed from this.This is the approach we adopt here.
As an example we consider a π(ω) of the form This is, perhaps, the simplest example that satisfies the necessary physical constraints: (i) it is normalised, (ii) π(0) = 0, and (iii) ω 2 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 Γ to be a real decay rate and γ ± are either two real decay rates or complex conjugates of one another.With the familiar classical behaviour in mind we set Direct evaluation of ω 2 for our example gives which links the parameters Γ, γ + and γ − to the short-time natural frequency Ω 0 .When written in terms of ω 0 and γ, this relationship simplifies to We can rewrite this expression in the form which emphasises the need for the decay rate Γ as a consequence of the fact that we require two natural frequencies, ω 0 and Ω 0 , to describe the quantum damped harmonic oscillator.It is straightforward to determine from π(ω) the two averages ω and ω −1 : from which we can derive the expectation values of x2 and p2 in the ground state of the damped oscillator ¶.
We can determine the evolution of the expectation values of the position and momentum from Eq (79) by evaluating the averages cos(ωt) , ω −1 sin(ωt) and ω sin(ωt) : e −γ − t ¶ When γ ± are complex, we can write these in the form We expect to find something approximating the evolution of the classical damped harmonic oscillator in the limit of small Γ, but departures from this for larger values.
In figure 2 we plot the evolution of cos(ωt) for a small value of Γ.As anticipated, we find over-damped like behaviour for real values of γ + and γ − , 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 −Γt term which, although small, decays slowly and so has a residual influence at long times.This is especially clear in the large Γ 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 −Γt term has a less dramatic effect; the evolution for Γ = 0.01 (in figure 2b) and for Γ = 10 (figure 3b) are qualitatively rather similar.
We have seen that the existence of two natural frequencies for the oscillator, ω 0 and Ω 0 , is particularly important when comparing the short and long time behaviour of the oscillator: ω 0 behaves as the natural frequency, but at short times it is Ω 0 that takes on this role.In figure 4 we plot the short-time evolution of cos(ωt) (solid line) and compare this with the classical evolution of for a classical damped harmonic oscillator with natural frequency ω 0 It is clear that the former falls off more quickly as it must, because Ω 0 > ω 0 .For the parameters chosen, Ω 0 is just over twice the value of ω 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, cos(ωt) 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.

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 α(ω) in the weak coupling limit.To aid our analysis we rewrite the form given in (B.18) as The weak damping limit corresponds to choosing the coupling to the environment to be small or, more specifically, to |V (ω)| 2 Ω 0 .It is clear that in this limit, |α(ω)| 2 will be a sharply peaked function centred around the frequency for which Y (ω) = 0.If the integral part in Y (ω), as given in (B.15) is small + then this frequency will be close to Ω 0 and we can write where This leads, in turn, to a corresponding approximate form for α(ω): where we have set ω = Ω 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 + If it is not then we will need to invoke ideas of renormalisation, a manageable complication, but one we wish to avoid.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 ω = Ω 0 leads us to set β(ω) to zero: so that the integral over all frequency of |α(ω)| 2 is unity.Moreover, for weak damping the thermal function coth(βhω/2) will also be slowly varying compared to the rapid variation of |α(ω)| 2 in the vicinity of ω = Ω 0 and we may replace this function by its value at Ω 0 .Hence in this limit the steady-state characteristic function for our oscillator is where 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, π(ω) ≈ |α(ω)| 2 , approaches a Lorentzian centred on Ω 0 , with some width γ * .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.

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, π(ω), which we may interpret as the contribution of corresponding dressed mode, at frequency ω to the oscillator.These properties include the steady state of the oscillator * This is not strictly true in the wings of the distribution, of course, as even in this limit we require 0 , but the corresponding quantity for a true Lorentzian is divergent.
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].The next step is to integrate the second of these equations of motion.The complementary function is To find the particular integral we need to make a small addition to (A.4) by adding a very weak damping term to give: and work in the limit as the strictly positive quantity ε tends to zero.It is this choice of a positive (if very small) value for ε that provides the irreversibility and hence the arrow of time.Solving (A.5) 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 xµ in the form with a similar expression for the Fourier transform of x.It follows that the Fourier transform of the particular integral part of the position is given by It follows then follows that which is the Fourier transform of the product of two functions of ω.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): if t > 0 and is zero otherwise.We can now take the limit as ε → 0 to give Kµ = − ω µ λ µ sin(ω µ t) where we have used the facts that K(T ) = 0 for T < 0 and we may take x(τ ) = 0 for τ < 0. Pulling this altogether we arrive at our desired Heisenberg-Langevin equation for the damped harmonic oscillator.From (A.3) we have where and F (t) is the Langevin force:

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 counterrotating 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 Hamiltonian by which we mean rewriting it in the form of a continuum of uncoupled or dressed oscillators: where C is an unimportant constant.We proceed by writing the dressed annihilation operators, B(ω), as linear combinations of the bare operators for the oscillator and bath modes: For the Hamiltonian of interest in this paper the coupling, V (ω), is real but treating the problem with a more general complex coupling presents no additional difficulties.
We require the operator B(ω) to be associated with an uncoupled or dressed oscillator of angular frequency ω.This requires us to find its form such that the following pair of operator equations are satisfied for every frequency, ω: Substituting the ansatz (B.3) into (B.4) and comparing coefficients of the bare creation and annihilation operators leads to the set of coupled equations: Our method of solution is use these to determine the functions β(ω), γ(ω, ω ) and δ(ω, ω ) in terms of α(ω) and then to determine this remaining function by enforcing the commutation relation (B.5).From (B.6) and (B.7) we see that If we use this to substitute for β(ω) into the remaining equations then we find Solving the second of these presents no difficulty and we find The first, however, requires careful handling because the behaviour at ω = ω .Following Fano [116], we adopt the method proposed by Dirac [117] and write where P denotes that the principal part is to be taken on integration and Y (ω) is a real function, which we determine by substituting (B.14) into (B.6).We find If we substitute our operators, B(ω), expressed in terms of the function α(ω) into the commutation relation (B.6) then we find Evaluating the integrals and setting the result equal to δ (B.17) Note that the diagonalisation does not fix the phase of the complex function α(ω) and we are free to choose this as we wish.A convenient choice is to set Appendix C. The symmetrically ordered characteristic function We have made use of the symmetrically ordered characteristic function, 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 ξ = 0 it reduces to the trace of ρ and it follows that χ(0) = 1.More generally, it is the expectation value of the unitary displacement operator: This operator, by virtue of its unitarity, has only eigenvalues of modulus 1 and it follows that with the maximum at ξ = 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 † † This requires the use of the following formula for the product of two principal parts [17]: transform pair.Cahill and Glauber [118] (see also [119]) exploited a theorem of Weyl [120]  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 (D.1).To construct this state we consider a doubled state space in which we introduce a second oscillator with annihilation and creation operators b and b † .The two-mode pure state, the thermal vacuum: 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, 0 , via a unitary transformation: Readers with a background in quantum optics may recognise |0(β) 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 In place of a coupling to a single harmonic oscillator in a thermal state with inverse temperature β, the transformed Hamiltonian has a coupling to a regular oscillator in its ground state with a coupling strength increased by coshθ(β) and a coupling to a second inverted oscillator in its most highly excited state, with the original coupling multiplied by sinhθ(β).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

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

4 ) 5 )
We can extract from the characteristic function the expectation value of any symmetrically ordered combination of â and â † :By symmetrically ordered, we mean the average of all possible orderings, for example:+ â † ââ † â + â † â2 â † + ââ †2 â + ââ † ââ † + â2 â †2 .(C.6) Appendix D. Thermofields It is simplest to consider first an isolated discrete oscillator with annihilation and creation operators b and b † .For such an oscillator in a thermal state at temperature T the density operator has the simple, diagonal form: ρT = (1 − e −βhω ) ∞ n=0 e −nβhω |n n| , (D.1) where β = (k B T ) −1 is the inverse temperature.The mean number of excitations is n = 1 e βhω − 1 (D.2)and we can write the density operator in terms of this mean: