Reference: https://arxiv.org/abs/2306.15013

There are many methods to study damped harmonic oscillators. The idea is pretty simple: one harmonic oscillator embedded in a bath of surrounding environmental harmonic oscillators each characterized by a different frequency. This spread of frequencies brings about dephasing effects in the main oscillator dynamics.

There are two regimes, as known in classical cases, of damped oscillation: weakly damped and strongly damped. The weak regime is mostly studied using the Linblad master equation, which assumed Markovian dynamics. We expect to see some oscillations before the system tends to equilibrium.

The Markov approximation in Linblad master equation typically means that, the time-scale of decay for the environment $\tau_{env}$ is much shorter than the smallest time-scale of the system dynamics $\tau_{sys}$. The Born approximation requires that the state of the environment does not significantly change as a result of the interaction with the system, and that the total state can be written as a tensor product $\rho_{tot}(t) = \rho_{sys}(t)\otimes \rho_{env}$.

The weak regime breaks down when the Markov approximation no longer holds, i.e. when the memory of the environment becomes important. This is when we enter the strong regime.

Background

The damped harmonic oscillator loses energy as a result of coupling to the surrounding environment. An undamped oscillator is described by the following equation of motion

\[ \begin{align} \ddot{x} + \omega_0^2 x = 0 \end{align} \]

To introduce damping effect, we consider a damping coefficient that is constant over time $\gamma$ and its associated external stochastic Langevin force $F(t)$.

Both $\gamma$ and $F(t)$ results from the interacting term in the Hamiltonian. Therefore there has to be some correlations between them.

The equation of motion is thus the familiar linear differential equation of the form

\[ \begin{align} \ddot{x} + \gamma\dot{x} + \omega_0^2 x = \frac{F(t)}{m} \end{align} \]

Note that knowledge regardings the external force $F(t)$ is not required.

The damped quantum harmonic oscillator requires that explicit account be taken of the quantum nature of the surrounding environment. To obtain the notorious master equation, one could make use of the rotating wave approximation that hinges on the weak damping condition, $\gamma \ll \omega_0$. The RWA breaks down as we move into strong damping regime, leading to non-Markovian dynamics.

Hamiltonian for the strongly-damped harmonic oscillator

Discrete reservoir

Let us consider a harmonic oscillator strongly coupled to its reservoir, modelled as large collection of oscillators, with a range of frequencies. The Hamiltonian describes the combined system is

\[\begin{align} H = \dfrac{p^2}{2m} + \dfrac{1}{2}m\Omega_0^2 x^2 + \sum_{\mu} \left(\dfrac{p_\mu^2}{2m_\mu}+\dfrac{1}{2}m_\mu\omega_\mu^2x_\mu^2\right) - \sum_\mu m_\mu\omega_\mu^2\lambda_\mu x_\mu x \end{align}\]

Upon completing the square, we can rewrite this in a more compact form,

\[\begin{align} H = \dfrac{p^2}{2m}+\dfrac{1}{2}m\omega_0^2 x^2 + \sum_\mu\left[\dfrac{p_\mu^2}{2m_\mu}+\dfrac{1}{2}m_\mu\omega_\mu^2(x_\mu - \lambda_\mu x)^2\right] \end{align}\]

where

\[\begin{align} \omega_0^2 = \Omega_0^2 - \sum_\mu \dfrac{m_\mu}{m}\omega_\mu^2 \lambda_\mu^2 \end{align}\]

Since the kinetic term is always positive, the Hamiltonian is bounded from below if and only if $\omega_0^2$ is positive, i.e. $\omega_0$ must be real. This imposes a physical restriction to the strength of the coupling

\[\begin{align} \Omega_0^2 > \sum_\mu\dfrac{m_\mu}{m}\omega_\mu^2\lambda_\mu^2, \end{align}\]

This is called the positivity constraint.

From this Hamiltonian a set of equations of motion can be written and ultimately linked by the Heisenberg-Langevin equation (see derivation).

\[ \begin{align} \boxed{\ddot{x}(t) + \int_0^t \kappa(t-\tau)\dot{x}(\tau)d\tau + \left[\Omega_0^2-\kappa(0)\right]x(t)+\kappa(t)x(0)= \dfrac{F(t)}{m}} \end{align} \]

where $F(t)$ is the Langevin force,

\[\begin{align} F(t) = \sum_\mu m_\mu \omega_\mu^2\lambda_{\mu}\left[x_\mu(0)\cos(\omega_\mu t) + \dfrac{p_\mu(0)}{m_\mu\omega_\mu}\sin(\omega_\mu t)\right], \end{align}\]

and $\kappa(t)$ is the memory kernel,

\[\begin{align} \kappa(t) = \sum_\mu \dfrac{m_\mu}{m}\omega_\mu^2\lambda_\mu^2\cos(\omega_\mu t). \end{align}\]

In ultra-short time scales, this EoM can be rewritten as

\[ \begin{align} \langle\ddot{x}(t)\rangle + \Omega_0^2\langle x(t)\rangle + O(\delta t^2) = 0 \end{align} \]

which, using the Taylor expansion around $\delta t \to 0$, leads to

\[ \begin{align} \langle x(\delta t) \rangle = \langle x(0)\rangle + \dfrac{\langle p(0)\rangle}{m} \delta t \newline \langle p(\delta t)\rangle = \langle p(0)\rangle -m\Omega_0^2 \langle x(0)\rangle \delta t \end{align}\]

This is basically the behavior of a short-term undamped harmonic oscillator (since we don’t have any decay term here).

For longer time scales, the memory kernel tends to zero. The EoM becomes

\[ \begin{align} \langle\ddot{x}(t)\rangle + \int_0^t \kappa(t-\tau)\langle\dot{x}(\tau)\rangle d\tau + \omega_0^2\langle x(t)\rangle = 0. \end{align} \]

The question whether this EoM describes non-Markovian effect lies in the second integral term. However, it’s clear that on this time scale the natural frequency of the harmonic oscillator is $\omega_0$. Also, because now the natural frequency is $\omega_0$, we can actually reduce $\omega_0\to 0^+$ such that $\omega_0 < \gamma/2$ to enter the strongly over-damped regime. In the Ohmic damping case, this corresponds to

\[ \begin{align} \omega_0^2 > 0\quad\text{and}\quad \Omega_0^2 > \kappa(0) = \dfrac{2\gamma\omega_c}{\pi} \end{align} \]

In summary,

Time-scale\Regime	Underdamped $f_0 > \gamma$	Overdamped $f_0 < \gamma$
Ultra-short $\Omega_0$	Yes. Because $\Omega_0$ has no upper bound.	Maybe. Because $\Omega_0$ has a lower bound.
Longer time scales $\omega_0$	Yes. Because $\omega_0$ has no upper bound.	Yes. Because $\omega_0$ has no lower bound (except that it has to be positive).

Continuum reservoir

First we rewrite the Hamiltonian in terms of creation and annihilation operators ($\hbar=1$)

\[ \begin{align} H = \Omega_0 a^\dagger a + \sum_\mu \omega_\mu b^\dagger_\mu b_\mu + \sum_\mu \dfrac{V_\mu}{2}(a+a^\dagger)(b_\mu + b_\mu^\dagger) \end{align} \]

where $V_\mu$ being

\[\begin{align} V_\mu = -\sqrt{\dfrac{m_\mu\omega_\mu}{m\Omega_0}}\omega_\mu\lambda_\mu \end{align}\]

The positivity condition is again,

\[\begin{align} \Omega_0 > \sum_\mu \dfrac{V_\mu^2}{\omega_\mu} \end{align}\]

To make this a continuum model, we introduce a commutator

\[ \begin{align} [b(\omega), b^\dagger(\omega’)] = \delta(\omega-\omega’). \end{align} \]

The Hamiltonian is now

\[ \begin{align} H = \Omega_0 a^\dagger a + \int_0^\infty \omega b^\dagger(\omega)b(\omega)d\omega + \int_0^\infty d\omega \dfrac{V(\omega)}{2}(a+a^\dagger)[b(\omega)+b^\dagger(\omega)] \end{align} \]

Since $H$ is quadratic in the annihilation and creation opearators for the oscillator and the reservoir, it leads to linear coupled equations of motion for these operators.

Hamiltonian diagonalization

Here we briefly summarize how to diagonalize the oscillator-reservoir Hamiltonian. This can be achieved by finding a complete set of eigenoperators $B(\omega)$ and their conjugates $B^\dagger(\omega)$ that satisfy

\[ \begin{align} [B(\omega), H]=\hbar\omega B(\omega),\newline [B^\dagger(\omega), H] = -\hbar\omega B^\dagger(\omega); \end{align} \]

for all positive $\omega$. If this set of $B$ also satisfy the orthogonality condition $[B(\omega),B^\dagger(\omega’)]=\delta(\omega-\omega’)$, then this set of eigenoperators is complete.

Suppose that we know have $B$, then we can express the annihilation and creation operators of both the reservoir and the oscillator in terms of $B$ and $B^\dagger$. For example,

\[ \begin{align} a = \int_0^\infty d\omega\left[\alpha^*(\omega)B(\omega)-\beta(\omega)B^\dagger(\omega)\right],\newline a^\dagger = \int_0^\infty d\omega \left[\alpha(\omega)B^\dagger(\omega)-\beta^* (\omega)B(\omega)\right] \end{align} \]

Combined with the requirement such that $[a, a^\dagger]=1$, this leads to

\[ \begin{align} \int_0^\infty \left[|\alpha(\omega)|^2-|\beta(\omega)|^2\right] = \int_0^\infty d\omega |\alpha(\omega)|^2\dfrac{4\Omega_0\omega}{(\Omega_0+\omega)^2}=1, \end{align} \]

where we’ve used

\[\begin{align} \beta(\omega) = \dfrac{\omega-\Omega_0}{\omega+\Omega_0}\alpha(\omega) \end{align}\]

that resulted from the commutator relations $[B,B^\dagger]$ and $[B(B^\dagger),H]$. Let us now denote

\[ \begin{align} \pi(\omega) = |\alpha(\omega)|^2\dfrac{4\Omega_0\omega}{(\Omega_0+\omega)^2} \end{align} \]

as a frequency probability distribution. This will come in handy in later analysis. Remember that this function $\pi(\omega)$ comes from the fact that we can express $a$ and $a^\dagger$ in terms of the eigenoperators in combination with bosonic commutation relation. Some physical interpretation of $\pi$ is delayed until next section. For now, we explore what can be learned from viewing this purely as a frequency distribution.

An inequality

Let us denote the average value of a function $f$ under the $\pi$ distribution,

\[ \begin{align} \langle\langle f \rangle\rangle = \int_0^\infty d\omega f(\omega)\pi(\omega), \end{align} \]

From the commutator $[B(\omega), H]=\hbar\omega B(\omega)$ and the orthogonality relation $[B(\omega), B^\dagger(\omega’)]=\delta(\omega-\omega’)$, we can write the Hamiltonian in terms of the eigenoperators

\[ \begin{align} H = \int_0^\infty d\omega \hbar \omega B^\dagger(\omega)B(\omega) + C, \end{align} \]

where $C$ is a formally divergent term. Note that both $B, B^\dagger$ and $a, a^\dagger, b, b^\dagger$ are valid bases, hence we can also write

\[\begin{align} B(\omega) = \alpha(\omega)a + \beta(\omega)a^\dagger +\int_0^\infty d\omega’\left[\gamma(\omega,\omega’)b(\omega’)+\delta(\omega,\omega’)b^\dagger(\omega’)\right] \end{align}\]

Substituting the new expression for $B$ and $B^\dagger$ into the new $H$, we find that, since $a$ and $a^{\dagger 2}$ must vanish,

\[\begin{align} \int_0^\infty d\omega\ \omega |\alpha(\omega)|^2 \left(\dfrac{\omega-\Omega_0}{\omega+\Omega_0}\right) = \int_0^\infty d\omega \dfrac{\pi(\omega)}{4\Omega_0}(\omega^2-\Omega_0^2) = 0\newline \Rightarrow \int_0^\infty \omega^2 \pi(\omega)d\omega = \Omega_0^2\newline \Rightarrow \langle\langle \omega^2 \rangle\rangle = \Omega_0^2 \end{align}\]

Note that $\langle\langle \omega \rangle\rangle^2 < \langle\langle \omega^2 \rangle\rangle$ for any distrubution, thus

\[\begin{align} \langle\langle \omega \rangle\rangle^2 &< \Omega_0^2\newline \langle\langle \omega \rangle\rangle &< \Omega_0 \end{align}\]

We arrive at a lower bound on the value of $\langle\langle\omega^{-1}\rangle\rangle$,

\[\begin{align} \langle\langle \omega^{-1}\rangle\rangle \langle\langle \omega\rangle\rangle \leq 1\newline \langle\langle \omega^{-1}\rangle\rangle > \dfrac{1}{\Omega_0}. \end{align}\]

Ground-state

In an undamped harmonic oscillator, the ground state is the state $|0\rangle$ such that

\[\begin{align} a|0\rangle = 0. \end{align}\]

The Hamiltonian derived about suggests this however is not the case for damped oscillator-reservoir system. Now, since we’ve written the Hamiltonian like

\[ H = \int_0^\infty d\omega\ \hbar\omega B^\dagger(\omega)B(\omega) \]

It’s suggestive that we choose the state $|0\rangle$ such that it is annihilated by $B(\omega)$,

\[ |0\rangle : B(\omega)|0\rangle = 0,\quad\quad\forall \omega. \]

The ground-state of the interested oscillator in this purestate $|0\rangle$ is a mixed state obtained by tracing out the environment,

\[ \begin{align} \rho_O = \text{tr}_E\left(|0\rangle\langle 0|\right) \end{align} \]

Using such density operator $\rho_O$ to calculate stuff is not very wise. Instead, we employ the characteristic fuction $\chi(\xi)$,

\[ \begin{align} \chi(\xi) = \text{tr}\left[\rho\exp\left(\xi a^\dagger - \xi^* a\right)\right] \end{align} \]

where $\rho=|0\rangle\langle 0|$. The characteristics function can be used to evaluate any operator function of $a$ and $a^\dagger$. In fact, information regardings the characteristic function is enough to reconstruct the density matrix $\rho$. In our context, with the probability distribution function $\pi(\omega)$, the characteristic function reads

\[ \begin{align} \chi(\xi) = \exp\left[-\dfrac{1}{2}\left(\dfrac{\langle\langle\omega\rangle\rangle}{\Omega_0}\Re{(\xi)}^2+\langle\langle\omega^{-1}\rangle\rangle\Omega_0\Im{(\xi)}^2\right)\right] \end{align} \]

Ground-state energy

The ground-state energy of the undamped harmonic oscillator is simply $1/2\hbar\omega_0$. For an undamped oscillator, the oscillation frequency is

\[ \begin{align} \sqrt{\omega_0^2-\gamma^2/4} < \omega_0, \end{align} \]

implying that the ground-state energy is perhaps reduced. However, things might get tricky when it comes to the overdamped regime, where $\omega_0 < \gamma/2$, making the oscillation frequency imaginary. It turns out that the ground-state energy of a damped oscillator, whether $\omega_0$ or $\Omega_0$ being the natural frequency, increases in comparision to the ground-state energy of an undamped one, which is quite nice to think about.

We can proceed to investigate the reduced state of the oscillator, which is completely characterized by the average and the variance of position and momentum. These moments are, using the characteristic function

\[ \begin{align} \text{tr}\left[\rho x\right] &= 0,\newline \text{tr}\left[\rho p\right] &= 0,\newline \text{tr}\left[\rho x^2\right] &= \dfrac{\hbar}{2}\dfrac{\langle\langle\omega^{-1}\rangle\rangle}{m},\newline \text{tr}\left[\rho p^2\right] &= \dfrac{\hbar}{2}m\langle\langle\omega\rangle\rangle,\newline \end{align} \]

where we have made use of the fact that $x\propto (a+a^\dagger)$ and $p\propto (a-a^\dagger)$. Note that the value of $\langle\langle \omega\rangle\rangle$ is dependent on $\pi(\omega)$, which ultimately a function of $\Omega_0$ (also affecting $\omega_0$). Hence these moments are dependent on the natural frequency of the oscillator, be it $\omega_0$ or $\Omega_0$.

It is then instructive to define a new frequency $f_0=\{\omega_0, \Omega_0\}$. From second moments of $x$ and $p$, the energy determined from the reduced Hamiltonian of the oscillator is

\[\begin{align} \text{tr}\left[\rho H_{f_0}\right] = \dfrac{\hbar f_0}{4}\left(\dfrac{\langle\langle \omega\rangle\rangle}{f_0}+\langle\langle \omega^{-1}\rangle\rangle f_0\right) \end{align}\]

Note that $\langle\langle \omega^{-1}\rangle\rangle \leq \langle\langle\omega\rangle\rangle^{-1}$, the energy is thus bounded below by

\[ \begin{align} \text{tr}\left[\rho H_{f_0}\right] = \dfrac{\hbar f_0}{4}\left(\dfrac{\langle\langle \omega\rangle\rangle}{f_0}+\langle\langle \omega\rangle\rangle^{-1} f_0\right) \end{align} \]

where

\[ \text{min}\left[\dfrac{\hbar f_0}{4}\left(\dfrac{\langle\langle \omega\rangle\rangle}{f_0}+\langle\langle \omega\rangle\rangle^{-1} f_0\right)\right] = \dfrac{1}{2}\hbar f_0. \]

Thus for any choice of frequency $f_0$, the oscillator energy can never be smaller than the vacuum state energy of an undamped oscillator.

This reflects the fact that there is an energy cost to be paid in order to decouple the oscillator from its environment (when tracing out the environmental degrees of freedom).

Physical meaning of $\pi(\omega)$

tbd.

Equilibrium state at finite temperature

tbd.

Oscillator dynamics

Now that we’ve obtained the equilibrium state, which is also the steady state of the strongly damped oscillator, what remains is to consider the evolution of the oscillator towards equilibrium.

tbd.