%SourceDoc ./tesis.tex
\chapter{Multi-Mode Field Dynamics in Optical Resonators}
\label{chp:F-D}

Perhaps the most studied model of damping in quantum mechanics is that
of a single harmonic oscillator linearly coupled to a continuum of
harmonic oscillators \cite{Sarg74,Wall94,Gard00}. This exactly
solvable model describes the damping and noise for a \emph{single-mode} field
inside an open optical resonator due to the interaction with the
external field. In the textbooks the dynamics of the single-mode field is usually regarded only in the weak coupling
limit \cite{Sarg74,Wall94,Gard00}, yet the quantum properties of this
field for arbitrary damping strength and arbitrary heat-bath
temperature are known \cite{Ulle66a,Rise85,Haak85}. Less understood
are the damping and noise properties of \emph{multi-mode} fields in
optical resonators. Until quite recently investigations on this
subject were hampered by the lack of a quantization method that could
properly described the field leakage in cavities with overlapping
resonances. In the preceding chapter such a field quantization method
was implemented and an exact quantum representation of the
electromagnetic fields and the Hamiltonian systems was obtained. Based
on these results, we consider in this chapter the dynamics of
multi-mode fields in optical cavities.

Our starting point is the system-and-bath Hamiltonian
\eqref{eq:hlast}, which corresponds to a discrete set of harmonic
oscillators linearly coupled to the same single continuum (or even
several continua distinguished by the ``channel'' index). The
Hamiltonian \eqref{eq:hlast} is then a straight forward generalization
of the single-mode model; however the dynamics that it generates is of
a very different nature. Due to the coupling to one and the same bath,
the cavity modes does not evolve independently in time but become
coupled through both damping and noise; the emerging correlations
among the resonator modes are thus seen to be responsible for the
excess noise characteristic of multi-mode fields. Still the equations
of motion for the field operators are linear (since the Hamiltonian is
a bilinear form of the field operators) and can be integrated exactly
(up to a frequency integral that depends on the density of states of
the system at hand). Explicit expressions for the cavity and channel
field operators are given in Appendix \ref{app:2A}; they are the
generalization of the single-mode solution in
\cite{Ulle66a,Rise85,Haak85} to multi-mode fields coupled to a
continuum.

Measurements on optical resonators are typically done with detectors
located in the external region outside the cavity. The detector
therefore measures the external field, and an input--output theory is
required to relate the time evolution of the external field to the
dynamics of the system of interest. In this chapter we implement the
input--output formalism of Gardiner and Collet
\cite{Gard85,Gard00,Wall94} for multi-mode fields in resonators with
overlapping resonances. We establish the input--output relation using
the exact field equations of motion; thus the results apply even in
the case of field overdamping. We then use the input--output relation
to derive the system $S$ matrix. In particular, we will address the
case of chaotic scattering, i.e., the case when the propagation of
light in the cavity becomes chaotic due to random fluctuations of the
refractive index or due to scattering at irregular shaped mirrors. As
an application of our input--output theory we consider a linear
amplifier and show that our theory reduces to the scattering theory of
Beenakker and coworkers \cite{Been98,Been99,Patr00c} for light
propagation in random media. Additionally, the input--output relation
sets the ground for the introduction of Langevin equations to describe
the cavity field dynamics.

The main result of this chapter is the derivation of stochastic
equations for the field dynamics in resonators with overlapping
resonances. Restricting our analysis to optical fields, we can employ
the rotating wave approximation and a Markov approximation for the
system memory. Equations of motion for the cavity field are then
obtained both within the Heisenberg picture (in terms of quantum
Langevin equations) and within the Schr\"odinger picture (employing a
master equation for the reduced density operator of the cavity modes)
which include explicit expressions for all damping and noise
forces. We show that the resulting stochastic equations correctly
describe isolated as well as spectrally overlapping
resonances. Corrections would only be important if the resonance
widths were to become comparable to the resonance frequencies
themselves, a regime not encountered in optical resonators. Our
derivation clarifies the status of previously proposed stochastic
equations
\cite{Lamp99,Bard99,Bard00,Bard00a,Lee00,Lamp02,Lamp02a} which were
not derived from a rigorously quantized electromagnetic field.

In the last section of the chapter we address a representation of the
cavity electromagnetic field in terms of nonorthogonal modes. Such
modes correspond to the ``true'' resonances of the open cavity, and
arise naturally in our derivation of the stochastic equations. It is
shown that the deterministic part of the time evolution of the field
is generated by a non-Hermitian operator $\mathcal{H}$ whose complex
eigenvalues are the cavity resonances and whose eigenvectors represent
the cavity resonance modes. The field dynamics in this representation
takes a simpler form than within the orthogonal mode representation,
however one has to pay the price that the operators related to the
nonorthogonal modes are no longer standard creation and annihilation
operators.

In Sec.~\ref{sec:exactdyn} we consider the exact field dynamics in
open resonators and derive the system $S$ matrix. We then introduce
the rotating wave approximation and the Markov approximation in
Sec.~\ref{sec:LEQM} to derive the Langevin equations of motion. The
results are applied to linear amplifiers. The master equation is
derived in Sec.~\ref{sec:MEQ} and the equivalent Fokker--Plank
equation in Sec.~\ref{sec:FPEQ}. The equivalence between the Langevin
and master equation is demonstrated. We conclude with
Sec.~\ref{sec:nonorth} were the field dynamics is represented in terms
of nonorthogonal modes.

Some of the results in this chapter have been reported in
Refs.~\cite{Vivi03,Hack03}. 


\section{Exact Dynamics}
\label{sec:exactdyn}

In this section we present an exact input--output relation for
multi-mode fields in open cavities. The relation is for arbitrary
values of the outcoupling amplitudes and holds independently of the
nature of the cavity field dynamics (whether linear or
nonlinear). When the dynamics is linear, the input--output relation is
employed to derive the $S$ matrix of the system.


\subsection{Equations of Motion}
\label{sec:exacteqm}

To start with, we concentrate on the intracavity field dynamics. The
equations of motion for the cavity operators are obtained by first
formally integrating the equations of motion for the channel operators
and then substituting the solution into the equations for the resonator
operators; the time evolution of the intracavity field is thus seen to
depend on the channel field only through its initial state. Using
Hamiltonian \eqref{eq:hlast}, the Heisenberg equations of motion for
the internal and external annihilation operators $a_{\lambda}$ and
$b_{m}(\omega)$, respectively, are
\begin{align}
\label{eq:aeqm}
\dot{a}_{\lambda} (t)&= -\rmi \omega_\lambda a_{\lambda}(t) - \rmi \sum_m 
\int\limits_{\omega_{m}}^{\infty} \! \rmd \omega \: \left[ \mathcal{W}_{\lambda
  m}(\omega) b_m(\omega,t) + \mathcal{V}_{\lambda
  m}^{*}(\omega) b_m^{\dagger}(\omega,t) \right], \\
\label{eq:beqm}
\dot{b}_m(\omega,t) &= -\rmi\omega b_m(\omega) - \rmi \sum_{\lambda} \left[
  \mathcal{W}^{*}_{\lambda m}(\omega) a_{\lambda}(t) +
  \mathcal{V}^{*}_{\lambda m}(\omega) a_{\lambda}^{\dagger}(t)\right].
\end{align}
Due to the presence of the antiresonant terms in the
system Hamiltonian the equations of motion for the creation and
annihilation operators are coupled. For an initial time $t_{0} < t$, formal
integration of the equations of motion for the channel operators
yields
\begin{equation}
\label{eq:bineq}
{b}_m(\omega,t) = \rme^{-\rmi\omega (t-t_0)} b_m(\omega,t_0) - \rmi
\sum_{\lambda} \int\limits^{t}_{t_0} \! \rmd t^{\prime} \:
\rme^{-\rmi\omega(t-t^{\prime})} \left[
  \mathcal{W}^{*}_{\lambda m}(\omega) a_{\lambda}(t') +
  \mathcal{V}^{*}_{\lambda m}(\omega) a_{\lambda}^{\dagger}(t') \right],
\end{equation}
where $b_{m}(\omega,t_0)$ denotes the operator $b_m(\omega)$ at
time $t_0$. In an analogous fashion, one can express $b_{m}(\omega,t)$
in terms of the operators at the final time $t_1>t$,
\begin{equation}
\label{eq:bouteq}
{b}_m(\omega,t) = \rme^{-\rmi\omega (t-t_1)} b_m(\omega,t_1) + \rmi
\sum_{\lambda} \int\limits^{t_1}_{t} \! \rmd t^{\prime} \:
\rme^{-\rmi\omega(t-t^{\prime})} \left[
  \mathcal{W}^{*}_{\lambda m}(\omega) a_{\lambda}(t') +
  \mathcal{V}^{*}_{\lambda m}(\omega) a_{\lambda}^{\dagger}(t') \right] .
\end{equation}
Although in this section we are only interested in the time evolution
of the cavity field for times $t> t_{0}$, Eq.~\eqref{eq:bouteq}
is necessary to establish the input--output relation for the system
below. Substituting Eq.~\eqref{eq:bineq} into Eq.~\eqref{eq:aeqm}, and
after some straightforward manipulations for which the identities
\begin{equation}
\label{eq:VtoW}
\mathcal{V}_{\lambda m}^{*}(\omega) =
\sum_{n}\int\limits_{\omega_{n}}^{\infty} \! \rmd \omega'
\mathcal{N}_{mn}^{\dagger}(\omega,\omega') \mathcal{W}_{\lambda
n}(\omega') = \sum_{\lambda'}\mathcal{N}_{\lambda \lambda'}^{\dagger}
\mathcal{W}_{\lambda' m}^{*}(\omega)
\end{equation}
are required, the equations of motion for the cavity annihilation
operators can be written as
\begin{equation}
\label{eq:aexact}
\dot{a}_{\lambda}(t) = -\rmi\omega_{\lambda} a_{\lambda}(t) -
\sum_{\lambda''} \int\limits_{t_{0}}^{t} \! \rmd t' \tilde{C}_{\lambda
\lambda'}(t-t') \biggl[ a_{\lambda'}(t') +
\sum_{\lambda'}\mathcal{N}_{\lambda' \lambda''}^{\dagger}
a_{\lambda''}^{\dagger}(t') \biggr] + \tilde{F}_{\lambda}(t).
\end{equation}
where the effects of the channel field on the dynamics of the cavity
operators are encoded in the operators $\tilde{F}(t)$ and the matrix
$\tilde{C}(t)$. The operators $\tilde{F}(t)$ are defined by
\begin{multline} 
\tilde{F}_{\lambda}(t) = \\
 - \rmi \sum_{m} \int\limits_{\omega_{m}}^{\infty} \!
\rmd \omega \, \mathcal{W}_{\lambda m}(\omega) \biggl[ \rme^{-\rmi
\omega (t -t_0)} b_{m}(\omega,t_0) +  \rme^{\rmi \omega (t -t_0)}
\sum_{n} \int\limits_{\omega_{n}}^{\infty} \! \rmd \omega'
\mathcal{N}_{mn}^{\dagger}(\omega, \omega')
b_{n}^{\dagger}(\omega',t_0) \biggr],
\end{multline}
and we notice that the channel operators enter only through their
values at the initial time $t_{0}$. For appropriated states of the
external system, $\tilde{F}(t)$ may be interpreted as noise
forces. The matrix elements
\begin{equation}
\tilde{C}_{\lambda \lambda'}(t) = -2 \rmi \sum_{m}
\int\limits_{\omega_{m}}^{\infty} \! \rmd \omega \,
\mathcal{W}_{\lambda m} (\omega) \mathcal{W}_{\lambda' m}^{*} (\omega)
\sin(\omega t), 
\end{equation}
on the other hand, are the channel response functions \cite{Haak85}
$\tilde{C}_{\lambda \lambda'}(t) = \bigl\langle
[\tilde{F}_{\lambda}(t), \tilde{F}_{\lambda'}^{\dagger}(0)]
\bigr\rangle, $ where $\langle\dots\rangle$ denotes the quantum
average over the initial channel state. The matrix $\tilde{C}$
introduces the response time of the channel field as a relevant time
scale in the dynamics of the cavity field. In order to treat the
channel field as a bath this time scale must be smaller than the
damping time for the cavity operators; this condition is always
accomplished provided the coupling is weak enough \cite{Gard00}.

The equations of motion \eqref{eq:aexact} hold for an arbitrary
strength of the coupling between the resonator and the channel field
and are suitable to describe even the case of overdamped resonator
dynamics, but they are far from being simple. Due to the time integral
on the right hand side the dynamics is not Markovian; the state of the
cavity field at one instant does not suffice to determine the future
state of the field. Simplifications may occur in the optical domain,
which is the most relevant regime for this thesis. There the
oscillation frequencies are typically much larger than the field
damping rates, and Eq.~\eqref{eq:aexact} can be simplified by means of
the \emph{rotating wave approximation}, which amounts to neglecting
the fast oscillating counter-rotating terms. This approximation
provides a description of the field in a narrow bandwidth (small
compared to the frequencies themselves) around a central frequency of
interest (see Sec.~\ref{sec:LEQM}).


\subsection{Input--Output Relation}
\label{sec:in-out-thexact}

The time evolution of the external field is related to the dynamics of
the internal field by the input--output relation
\cite{Gard00,Gard85,Wall94}. The starting point of the input--output
theory are the equations of motion \eqref{eq:aeqm} and \eqref{eq:beqm}
for the intracavity and channel field modes. Formal integration of the
equations for the channel operators, both for an initial and
a final condition, yields Eqs.~\eqref{eq:bineq} and
\eqref{eq:bouteq}. Subtracting the former from the latter and
integrating the result over frequency we arrive at the input--output
relation in the time domain
\begin{equation}
\label{eq:in-outt}
b_{m}^{\text{out}}(t) - b_{m}^{\text{in}}(t) = -\frac{\rmi}{2\pi}
\sum_{\lambda} \int\limits_{\omega_{m}}^{\infty} \! \rmd \omega \:
\int\limits_{t_0}^{t_1} \! \rmd t' \: \rme^{-\rmi\omega(t-t')} \bigl[
  \mathcal{W}_{\lambda m}^{*}(\omega) a_{\lambda}(t') +
  \mathcal{V}_{\lambda m}^{*}(\omega) a_{\lambda}^{\dagger}(t')
  \bigr].
\end{equation}
Equation \eqref{eq:in-outt} relates the cavity operators in the time
interval $t_0 < t < t_1$ to the external input and output field
operators
\begin{subequations}
\label{eq:aout}
\begin{eqnarray}
b^{\textrm{in}}_m(t) & \equiv & \frac{1}{2\pi} 
\int\limits_{\omega_{m}}^{\infty} \! \rmd\omega \:\rme^{-\rmi\omega
(t-t_0)} b_m(\omega,t_0), \\ 
b^{\textrm{out}}_m(t) & \equiv & \frac{1}{2\pi}
\int\limits_{\omega_{m}}^{\infty} \! \rmd\omega \:\rme^{-\rmi\omega
(t-t_1)} b_m(\omega,t_1). 
\end{eqnarray}
\end{subequations}
In analogy with the scattering formalism, we are eventually interested
in the limits $t_0 \to -\infty$ and $t_1 \to \infty$. Using the
identity \eqref{eq:VtoW} the input--output relation can be written as
\begin{subequations}
\begin{align}
\label{eq:in-outta}
b_{m}^{\text{out}}(t) - b_{m}^{\text{in}}(t) &= -\frac{\rmi}{2\pi}
  \sum_{\lambda} \int\limits_{\omega_{m}}^{\infty} \! \rmd \omega \:
  \int\limits_{t_0}^{t_1} \! \rmd t' \: \rme^{-\rmi\omega(t-t')}
  \mathcal{W}_{\lambda m}^{*}(\omega) \Bigl[ a_{\lambda}(t') +
  \sum_{\lambda'} \mathcal{N}_{\lambda \lambda'}^{\dagger}
  a_{\lambda'}^{\dagger}(t') \Bigr], \\
\label{eq:in-outtq}
&= -\frac{\rmi}{2\pi} \sum_{\lambda}
  \int\limits_{\omega_{m}}^{\infty} \! \rmd \omega \:
  \int\limits_{t_0}^{t_1} \! \rmd t' \: \rme^{-\rmi\omega(t-t')}
  \frac{1}{\sqrt{2\hbar\omega}} W_{\lambda m}^{*}(\omega)
  Q_{\lambda}(t'), 
\end{align}
\end{subequations}
where we used definition \eqref{eq:qbres} to obtain the last
line. That input and output annihilation operators in
Eq.~\eqref{eq:in-outta} are related through a linear combination of
cavity creation and annihilation operators is consequence of the exact
field dynamics which mixes creation and annihilation operators. More
precisely, Eq.~\eqref{eq:in-outtq} shows that the resonator field
dynamics couples the exterior bath to the position operators
$Q_{\lambda}$ defined in Eq.~\eqref{eq:aaham}. Fourier transformation of
Eq.~(\ref{eq:in-outtq}) yields the input--output relation in the
frequency domain. In the asymptotic limit $t_0 \to -\infty$ and $t_1
\to \infty$ the result takes the form
\begin{eqnarray}
\label{eq:in-outw}
b^{\textrm{out}} (\omega) - b^{\textrm{in}} (\omega) = -
\frac{\rmi}{\sqrt{2\hbar\omega}} W^\dagger
(\omega)Q_{\textrm{cav}}(\omega),
\end{eqnarray}
where we combined the input and output operators at frequency $\omega$
to $M$--component vectors.  The coupling amplitudes $W_{\lambda m}$
form an $N \times M$ coupling matrix $W$, and the cavity position
operators an $N$--component vector. The finite number $N$ of cavity
modes is artificial, and will eventually be taken to infinity. We note
that the input and output operators are simply related to the channel
mode operators at time $t_0$ and $t_1$, respectively, $b^{\textrm{in}}
(\omega) = \rme^{i\omega t_0} b (\omega,t_0)$ and $b^{\textrm{out}}
(\omega) = \rme^{i\omega t_1} b (\omega,t_1)$. 

It is a remarkable feature of the input--output relations
\eqref{eq:in-outt} and \eqref{eq:in-outw} that no assumption
whatsoever was made concerning the nature of the dynamics of the
intracavity field. The relations hold, not only for empty cavities but
also in the case when the internal field is coupled to other systems,
e.g., an ensemble of intracavity atoms, leading to a nonlinear
dynamics. However, in the simple case of a linear cavity dynamics
simplifications occur that will be addressed in the following section.

 
\subsection{$S$ Matrix}
\label{sec:Sm}

For linear systems one can eliminate the cavity modes from the
equations of motion to derive a linear relation between the input and
output field. Here we consider the empty cavity described by
Hamiltonian \eqref{eq:hlast}. Fourier transformation of the exact
result for the cavity position operator \eqref{eq:Qexact} yields
$Q_{\lambda}(\omega)$. Since the initial state has been specified in
the remote past $t_0 \to -\infty $, only the oscillatory terms are
kept. The result is
\begin{equation}
\label{eq:Qinout}
Q_{\textrm{cav}}(\omega) = \frac{\pi}{c^2}\sqrt{\frac{2\hbar}{\omega}}
\mathcal{G}_{\mathcal{QQ}}(\omega) W(\omega) b^{\textrm{in}}(\omega).
\end{equation}
Here $\mathcal{G}_{\mathcal{QQ}}(\omega)$ is the $N \times N$ matrix
representation of the projected resolvent \eqref{eq:Gqq}, the matrix
elements of its inverse (cf.\ Eq.~\eqref{eq:Gqq}) are
\begin{equation}
\label{eq:w-Leff}
\mathcal{G}_{\mathcal{QQ}}^{-1}(\omega)_{\lambda \lambda'}
= \frac{1}{c^2} \left[ \bigl( \omega^2 - \omega_{\lambda}^{2}
\bigr)\delta_{\lambda \lambda'} + \mathcal{P}
\int\limits_{\omega_{m}}^{\infty} \! d\omega' \frac{
\bigl[W(\omega')W^{\dagger}(\omega')\bigr]_{\lambda \lambda'}}{\omega^{\prime 2} -
\omega^{2}} + \frac{i\pi}{2\omega}
\bigl[W(\omega)W^{\dagger}(\omega)\bigr]_{\lambda \lambda'} \right].
\end{equation}
Combining Eqs.~\eqref{eq:in-outw} and \eqref{eq:Qinout} one can
eliminate the cavity operators from the input--output relation. This
yields a linear relation between the incoming and outgoing field,
\begin{equation}
\label{eq:slinear}
b^{\textrm{out}} (\omega) = S(\omega) b^{\textrm{in}} (\omega)
\end{equation}
where $S(\omega)$ is the $M\times M$ scattering matrix,
\begin{equation}
\label{eq:weids}
S(\omega) = \openone - \frac{\rmi \pi}{\omega c^2}
W^{\dagger}(\omega) \mathcal{G}_{\mathcal{QQ}}(\omega)W(\omega).
\end{equation}
This representation of the $S$-matrix for optical systems can be shown
to be equivalent to the representation used in nuclear and condense
matter physics \cite{Maha69,Been97,Ditt00} where the $S$ matrix is
parameterized by the energy. The later is recovered from
Eq.~\eqref{eq:weids} by replacing $(\go/c)^2 \to E$. Equation
\eqref{eq:weids} is a generalization of the usual Breit--Wigner result
for a single resonance to the scattering in the presence of $N$
resonances. Using the commutation relations $[b_n(\omega),b_m^\dagger
  (\omega^\prime)] = \delta_{nm} \, \delta (\omega - \omega^\prime)$
for both $b = b^{\textrm{in}}$ and $b = b^{\textrm{out}}$, one can
easily show that $S$ is unitary, $S S^\dagger = S^\dagger S
=\openone$.  The $S$ matrix describes scattering both in the limit of
isolated resonances and in the regime of overlapping resonances. In
the context of quantum optics, the regime of isolated resonances
corresponds to the weak-damping or almost closed cavity regime, where
all matrix elements of $WW^{\dagger}/\omega^2$ are much smaller than
the mean frequency separation of the resonator modes. The opposite
regime of overlapping resonances is realized when the damping rates
exceed the mean frequency spacing.

According to Eqs.\ (\ref{eq:slinear}) and (\ref{eq:weids}) the field
dynamics is governed by the {\em resonances} of the open cavity. The
resonances are the complex poles of the $S$ matrix which coincide with
the poles of $\mathcal{G}_{QQ}(\omega)$. They are the solutions of the
equation $\textrm{det}\bigl[\omega^2 - c^2
L_{\textrm{eff}}(\omega)\bigr] =0$; from Eq.\ (\ref{eq:w-Leff}) one
can see that they represent the complex eigenvalues of the internal
field dynamics in the presence of damping inflicted by the coupling to
the external radiation field. We note that the resonances determine
the field dynamics even though the underlying field quantization is
formulated in terms of closed--cavity eigenmodes (with boundary
conditions as discussed in Sec.~\ref{sec:eigenmodes}).

Equation (\ref{eq:weids}) has found widespread application in the
random matrix theory of scattering \cite{Been97,Ditt00}. It is the
starting point for the so-called Hamiltonian approach to chaotic
scattering. This approach assumes that the Hamiltonian of a
closed chaotic resonator can be represented by a random matrix drawn
from a Gaussian ensemble of random matrix theory. The eigenvalues of
the internal Hamiltonian show level repulsion and universal
statistical properties. The statistics of the scattering matrix is
derived from the distribution of the internal Hamiltonian using
Eq.~(\ref{eq:weids}).  An alternative approach to chaotic scattering
is the random $S$ matrix approach in which one directly models the
statistical properties of $S$ without introduction of a
Hamiltonian. Based on the latter approach, Beenakker and coworkers
\cite{Been98,Been99,Patr99,Patr00b,Mish01} recently computed the noise properties of
disordered and chaotic optical resonators. The Hamiltonian and the $S$
matrix approach to chaotic scattering are known to be equivalent.
However, the Hamiltonian approach has the advantage that one can
include the interaction with an atomic medium on a microscopic
level. An example is given in Sec.\ref{sec:atomic_decay} where we
study the spontaneous emission of an atom inside an open resonator.


\section{Langevin Equations in the Optical Domain}
\label{sec:LEQM}


We now turn to a description of the field dynamics in the optical
regime. The equations of motion for the cavity mode operators simplify
considerable under two, not completely independent, approximations
often valid in this regime: the \emph{rotating wave approximation} and
\emph{the Markov approximation}. The rotating wave approximation
describes the system dynamics in a spectral bandwidth which is small
compared with the typical oscillation frequencies but that still
contains a large number of cavity resonances. In this spectral band
one then assumes the cavity field damping rates $\gG$ to be smaller
than the field oscillation frequencies. The rotating wave
approximation amounts to neglecting the fast oscillating terms on the
scale of $1/\gG$ in the field equations of motion. Our implementation
of the rotating wave approximation covers the case of
spectrally-overlapping cavity modes, generalizing the standard
approach for a single resonator mode \cite{Gard00,Wall94} with weak
outcoupling. In a resonator with overlapping resonances the mode
spectral broadening is comparable to the mean frequency spacing
$\Delta \omega$ of the internal cavity resonances. Therefore
oscillatory terms with frequencies of the order of $\Delta \omega$
must be kept in the equations of motion leading to an effective
coupling among the modes through the damping. In the rotating wave
approximation the system-and-bath Hamiltonian takes the form
\begin{multline}
\label{eq:hrwa}
H_{\textrm{SB}} = \sum_{\lambda} \hbar \omega_{\lambda} \, a^{\dagger}_{\lambda}
a_{\lambda} + \sum_m \int \! \rmd\omega
\: \hbar\omega \, b_m^{\dagger} (\omega) b_m(\omega) \\
 + \hbar\sum_{\lambda} \sum_m \int \! \rmd\omega
\left[ \mathcal{W}_{\lambda m}(\omega) \, a^{\dagger}_{\lambda}
  b_{m}(\omega) + \text{H.c.} \right],
\end{multline}
which is obtained from the Hamiltonian \eqref{eq:hlast} by neglecting
the fast oscillating antiresonant terms $\mathcal{V}_{\lambda
  m}(\omega)a_{\lambda}b_{m}(\omega) + \mathcal{V}_{\lambda
  m}^{*}(\omega)a_{\lambda}^{\dagger}b_{m}^{\dagger}(\omega)$ and
extending the range of the frequency integrals from $-\infty$ to
$+\infty$. This extension is consistent with the rotating wave
approximation and necessary for the Markov approximation
\cite{Gard00}. Keeping the fast oscillating terms introduces
corrections to the dynamics of the order $\gG / \bar{\omega}$ and
$\Delta / \bar{\omega}$, where $\bar{\omega}$ is a typical frequency
in the bandwidth we are interested (see Appendix \ref{app:2B}). The
Markov approximation assumes a separation of time scales between the
typical life times of the cavity modes and the ``bath correlations
times'' (like the thermal time $\tau_{\textrm{bath}}^{\textrm{th}} =
\hbar/k_{B} T$) such that the former are much in excess of the later,
$1/\gG \gg \tau_{\text{bath}}$. This approximation simplifies the
equations of motion in a twofold way: First, the coupling amplitudes
$\mathcal{W}_{\lambda m}$ can be taken to be frequency independent over
the spectral bandwidth\footnote{This assumption implies
that the frequencies in the spectral range considered are far from the
threshold frequencies for the opening of further channels, see, e.g., Ref.~\cite{Fyod97}}. Second, on the same spectral range, the number of
photons in the bath may be assumed independent of the frequency. Together,
both approximations lead to a set of Langevin equations with white
noise that describe the resonator field dynamics for times $t >
\tau_{\textrm{bath}}$.


\subsection{Derivation of the Langevin Equations}
\label{sec:DLEQM}

The Heisenberg equations of motion in the rotating wave approximation
are obtained using the Hamiltonian \eqref{eq:hrwa} \footnote{The
Heisenberg equations of motions are linear in the cavity and channel
operators and can therefore be solved exactly. The results are given
in the Appendix \ref{app:2A}. We notice that the Langevin
equations obtained in this section could also be obtained by starting
from the exact solution \eqref{eq:solrwa} of the equations of motion
and then implementing the necessary approximations at that level, as
done in \cite{Ulle66a} for the single-mode coupled to a continuum
model. For multi-mode fields however such a path turns out to be much
more involved than directly dealing with the equations of motion.},
\begin{subequations}
\label{eq:eqmrwa}
\begin{align}
\label{eq:aeqmrwa}
\dot{a}_{\lambda} &= -\rmi \omega_\lambda a_{\lambda} - \rmi \sum_m 
\int \! \rmd \omega \; \mathcal{W}_{\lambda
  m}(\omega) b_m(\omega), \\ \label{eq:beqmrwa} \dot{b}_m(\omega) &=
-\rmi\omega b_m(\omega) - \rmi \sum_{\lambda} \mathcal{W}^{*}_{\lambda
m}(\omega) a_{\lambda}.  
\end{align} 
\end{subequations}
They coincide with Eqs.~\eqref{eq:aeqm} and \eqref{eq:beqm} if the
counter rotating terms are set equal to zero, $\mathcal{V}_{\lambda
m}(\omega) = 0$. By following the same reasoning than in
Sec.~\ref{sec:exacteqm} one arrives at the equations of motion for the
cavity annihilation operators
\begin{equation}
\label{eq:aRWA}
\dot{a}_{\lambda}(t) = -\rmi\omega_{\lambda} a_{\lambda}(t) -
\sum_{\lambda'} \int\limits_{t_{0}}^{t} \! \rmd t' C_{\lambda
\lambda'}(t-t') a_{\lambda'}(t') + F_{\lambda}(t).
\end{equation}
The noise force operators $F(t)$ and the bath response matrix $C(t)$ are
defined, respectively, by
\begin{align} 
F_{\lambda}(t) &= - \rmi \sum_{m} \int \!
\rmd \omega \, \mathcal{W}_{\lambda m}(\omega)  \rme^{-\rmi
\omega (t-t_0)} b_{m}(\omega,t_0),\\
C_{\lambda \lambda'}(t) &= \frac{1}{\pi} \int \! \rmd
\omega \, \gG_{\lambda \lambda'}(\omega)\rme^{-\rmi \omega t},
\end{align}
and the damping matrix $\gG_{\lambda \lambda'}(\omega)$ is given
by
\begin{equation}
\gG_{\lambda \lambda'}(\omega) = \pi
\bigl[\mathcal{W}(\omega)\mathcal{W}^{\dagger}(\omega)\bigr]_{\lambda
  \lambda'}.
\end{equation}
We now introduce a set of slowly changing variables by the
transformation
\begin{equation}
\label{eq:aslow}
a_{\lambda}(t) = \rme^{-\rmi \omega_{\lambda} t} \bar{a}_{\lambda}(t).
\end{equation}
Since it has been assumed that the mode oscillation frequencies are
much larger than the mode damping rates the dominant term in
\eqref{eq:aRWA} is $-\rmi \omega_{\lambda} a_{\lambda}(t)$; therefore
$\bar{a}_{\lambda}(t)$ must change slowly with time. In particular
$\bar{a}_{\lambda}(t)$ does not vary much on the time scale of
$\tau_{\textrm{bath}}$ in which the bath response matrix $C(t)$ decays
to zero. The equations of motion for the slowly varying variables are
obtained substituting Eq.~\eqref{eq:aslow} in Eq.~\eqref{eq:aRWA} and
replacing $\bar{a}_{\lambda}(t')$ with $\bar{a}_{\lambda}(t)$,
\begin{equation}
\dot{\bar{a}}_{\lambda}(t) = - \sum_{\lambda'} C_{\lambda
\lambda'}(\omega_{\lambda'}) \rme^{\rmi
  (\omega_{\lambda} - \omega_{\lambda'})t} \bar{a}_{\lambda'}(t) +
\rme^{\rmi \omega_{\lambda} t} F_{\lambda}(t).
\end{equation}
Here $C(\omega)$ is given by
\begin{equation}
\begin{split}
C_{\lambda \lambda'} (\omega_{\lambda'}) &= \int\limits_{-\infty}^{t}
\! \rmd t' C_{\lambda \lambda'}(t-t') \rme^{\rmi
  \omega_{\lambda'}(t-t')} \\ 
&= - \frac{\rmi}{\pi} \mathcal{P}
\int \! d\omega \frac{
  \gG_{\lambda\lambda'}(\omega)}{\omega - \omega_{\lambda'}} +
\gG_{\lambda\lambda'}(\omega_{\lambda'}),
\end{split}
\end{equation}
and we set $t_0\to - \infty$ because $C(t)$ is assumed to decay on the
time scale of the order of $\tau_{\textrm{bath}}$ while we are
interested in times $t-t_0>\tau_{\textrm{bath}}$. After moving back to
the original variables, the equations of motion in the rotating wave
approximation become
\begin{equation}
\label{eq:arwaf}
\dot{a}_{\lambda}(t) = -\rmi \omega_{\lambda} a_{\lambda}(t) -
\sum_{\lambda'} C_{\lambda \lambda'}(\omega_{\lambda'})
a_{\lambda'}(t) + F_{\lambda}(t).
\end{equation}
Equations \eqref{eq:arwaf} generalize the Langevin equations for
spectrally isolated modes \cite{Gard00,Wall94} to the case of overlapping
modes. We note that the resulting equations differ from the
independent oscillator equations of standard laser theory
\cite{Hake84} in two respects: First, the mode operators $a_{\lambda}$
are coupled due to the damping by the matrix $C(\omega)$. Second,
generically the noise force operators $F_{\lambda}$ are correlated,
$\bigl\langle F_{\lambda}(t), F_{\lambda'}^{\dagger}(t') \bigr\rangle
\neq \delta_{\lambda\gl'}$, as different modes coupled to the same
external channels ($\langle\dots\rangle$ denotes the quantum average
over the channel field at $t_0$). The mode coupling by both damping
and noise can be understood as a consequence of the
fluctuation--dissipation theorem.

The deviations from the independent oscillator dynamics may be
understood in the limiting case of weak damping.  This is the regime
where all matrix elements of $\gG$ are much smaller than the resonator
mode spacing $\Delta \omega$.  This regime is realized, e.g., in
dielectrics that strongly confine light due to a large mismatch in the
refractive index. To leading order in $\gG/ \Delta \omega$ only
diagonal elements contribute to the damping matrix, and
Eq.~(\ref{eq:arwaf}) reduces to the standard equation of motion for
independent oscillators. This shows that the independent--oscillator
dynamics is a limiting case of the true mode dynamics in the regime of
weak damping.

In obtaining the Langevin equation \eqref{eq:arwaf} no reference was
made to the statistic of the noise force $F(t)$. For arbitrary baths
Eq.~\eqref{eq:arwaf} may lead to non-Markovian dynamics for averaged
system observables. The dynamics simplifies if the Markov
approximation is employed. In that limit the frequency dependence of
the coupling amplitudes $\mathcal{W}$ can be dropped over the spectral
range in consideration, and the cavity mode operators $a_{\lambda}$
obeys the equation
\begin{equation}
\label{eq:aLEQ}
\dot{a}_{\lambda}(t) = -\rmi \omega_{\lambda} a_{\lambda}(t) -
\sum_{\lambda'} \gG_{\lambda \lambda'} a_{\lambda'}(t) + F_{\lambda}(t),
\end{equation}
where the small frequency shift in $C_{\lambda \lambda'}$, which is
rarely needed in practice, has been neglected. In the same limit one
arrives at white noise with $\la F_{\gl}(t) \ra = 0$, and the
correlation functions
\begin{subequations} 
\label{eq:f-2cf}
\begin{align}
\la F_{\gl}^{\dag}(t)F_{\gl'}(t')\ra &= 2\gG_{\gl'\gl} \, n_{\textrm
{th}} \, \delta(t-t'), \\
\la F_{\gl}(t)F_{\gl'}^{\dag}(t')\ra &= 2\gG_{\gl\gl'} (1+n_{\textrm
{th}})\,\delta(t-t') \, .
\end{align}  
\end{subequations} 
The second two-time correlation functions follows from the first and
the commutation relations for the channel operators $b_{m}$. The
remaining second moments vanish, $\la F_{\gl}(t)F_{\gl'}(t') \ra = 0 =
\la F_{\gl}^{\dag} (t) F_{\gl'}^{\dag} (t')\ra$, and higher-order
moments follow from the ones of orders 1 and 2 according to Gaussian
statistics. The thermal number of photons $n_{\textrm {th}}=
[\exp(\hbar\bar{\go}/kT)-1]^{-1}$ appearing in the second moments must
be taken as frequency independent throughout the spectral range under
consideration. Previously, Bardroff and Stenholm \cite{Bard00}
proposed a similar set of Langevin equations for the cavity field
coupled to a bath of atoms. Their equations are recovered from ours in
the limit $k T \ll \hbar
\bar{\omega}$.

We note that the deterministic part of the field dynamics in
Eq.~\eqref{eq:aLEQ} is generated by the \emph{non-Hermitian} matrix
\begin{equation}
\label{eq:Heff0} 
{\cal H}_{\gl\gl'}=\hbar \go_\gl\delta_{\gl\gl'} - \rmi \hbar
\gG_{\gl\gl'} \, . 
\end{equation}
Its Hermitian part corresponds to the dynamics of the closed system
describing the reversible dynamics, while its anti-Hermitian part
accounts for the irreversible dynamics due to the openness of the
system. For some applications it will prove helpful to rewrite the
above Langevin equation in the eigenbasis of the operator
$\mathcal{H}$. That basis is formed by the resonant modes of the
cavity and their eigenvalues correspond to the cavity resonances. As a
``penalty'' for that change of representation one has to work with
non-standard commutation relations for the operators associated with
the eigenvectors of ${\cal H}$ (see Sec.~\ref{sec:nonorth}).

At this point a remark concerning the field dynamics inside
chaotic optical resonators is at hand. According to the universality
hypothesis of chaotic scattering, the internal Hamiltonian of chaotic
resonators can be represented by a random matrix from the Gaussian
orthogonal ensemble of random-matrix theory \cite{Haak01}. The
eigenvalues $\omega_\lambda$ display level repulsion and universal
statistical properties. However, from Eqs.~\eqref{eq:aLEQ} and
\eqref{eq:Heff0}, the mode dynamics of open chaotic resonators is not
only determined by the eigenvalues of the internal Hamiltonian but
also by the coupling strength to the external radiation field: The
spectrum of such resonators is governed by a non--Hermitian
random matrix.  We thus encounter an interesting connection between
the spectral properties of chaotic optical resonators and
non--Hermitian random matrices \cite{Soko89,Fyod97,Chal97}.


\subsection{Inputs and Outputs}
\label{sec:in-outrwa}

Here we address the input--output theory within the rotating wave
approximation. Following a similar procedure as the one used in
Sec.~\ref{sec:in-out-thexact} one could re-derive the input--output
relation from the equations of motion generated by the Hamiltonian
\eqref{eq:hrwa}. The result, however, is seen to follow immediately
from Eq.~\eqref{eq:in-outta} upon neglecting the fast oscillating
terms proportional to the cavity creation operators $a_{\gl}^{\dag}$
and extending the frequency integration range from $-\infty$ to
$+\infty$. Then, for $t_{0} < t <t_{1}$, the input
and output field operators are related to the resonator operators by
\begin{equation}
\label{eq:in-outrwat}
b_{m}^{\text{out}}(t) - b_{m}^{\text{in}}(t) = -\frac{\rmi}{2\pi}
\sum_{\lambda} \int \! \rmd \omega \: \int\limits_{t_0}^{t_{1}} \! \rmd t'
\: \rme^{-\rmi\omega(t-t')} \mathcal{W}_{\lambda m}^{*}(\omega)
a_{\lambda}(t').
\end{equation}
We note that the input and output fields do not change in the rotating
wave approximation and are still given by Eq.~\eqref{eq:aout}. In the
frequency domain, the input--output relation is obtained by Fourier
transformation of Eq.~\eqref{eq:in-outrwat},
\begin{eqnarray}
\label{eq:in-outrwaw}
b^{\textrm{out}} (\omega) - b^{\textrm{in}} (\omega) = - \rmi
\mathcal{W}^\dagger (\omega)a(\omega),
\end{eqnarray}
where the asymptotic limit $t_0 \to -\infty$ and $t_1 \to \infty$ is
implied and, as in Eq.~\eqref{eq:in-outw}, we used a vector notation
for the operators $b^{\text{in/out}}(\go)$ and $a(\go)$, and the
coupling matrix $\mathcal{W}(\go)$.


\subsection{Linear Systems}
\label{sec:linear}

The $S$ matrix for a linear system simplifies in the rotating wave
approximation. Here we derive the $S$ matrix for an empty cavity and
then address the more general case of a cavity filled with a linear
absorbing and/or amplifying medium.


\subsubsection{S Matrix}

In Appendix \ref{app:2A} the equations of motion for the
electromagnetic field in the rotating wave approximation are solve
exactly. Here we are only interested in the asymptotic limit $t_{0}
\to -\infty$ and $t_{1} \to +\infty$. Therefore, for $a_{\gl}$ in
Eq.~\eqref{eq:arwaexact}, only the oscillatory terms proportional to
the input field must be kept. After Fourier transformation one finds
the relation between the cavity modes operators and the input
operators
\begin{equation}
\label{eq:aft}
a(\go) = 2 \pi D^{-1}(\go)\mathcal{W}(\go) b^{\textrm{in}}(\go),
\end{equation}
where $D$ is the $N \times N$ matrix with matrix elements 
\begin{equation}
\label{eq:D}
D_{\gl \gl'}(\go) = ( \omega - \omega_{\lambda} ) \delta_{\lambda \lambda'}
+ \mathcal{P} \int\limits \! d\omega' \frac{ \gG_{\lambda
\lambda'}(\go')}{\omega^{\prime} - \omega} + i\pi \gG_{\lambda \lambda'}(\go).
\end{equation} 
We now substitute Eq.~\eqref{eq:aft} into Eq.~\eqref{eq:in-outrwaw}
and eliminate the resonator operators from the input--output
relation. This yields a linear relation between the incoming and
outgoing field,
\begin{equation}
b^{\textrm{out}} (\omega) = S(\omega) b^{\textrm{in}} (\omega),
\end{equation}
involving the $M \times M$ scattering matrix in the rotating wave
approximation,
\begin{equation} \label{eq:SRWA}
S(\go) = \openone - 2\pi \rmi \mathcal{W}^{\dagger}(\omega)
D^{-1}(\omega) \mathcal{W}(\omega).
\end{equation}
The operator $D^{-1}(\go)$ corresponds to the rotating wave
approximation of the projected system resolvent
(cf. Eq.~\eqref{eq:weids}). We notice that if in addition the Markov
approximation is performed, and the frequency dependence of the
coupling amplitudes is neglected, Eq.~\eqref{eq:D} reduces to
\begin{equation}
\label{eq:DH}
D(\go) = \go - \mathcal{H}/\hbar
\end{equation}
where $\mathcal{H}$ is the non-Hermitian operator defined in
Eq.~\eqref{eq:Heff0}. Hence, in this approximation, $D^{-1}(\omega)$
corresponds to the resolvent of a system whose dynamics is
described by the effective Hamiltonian $\mathcal{H}$.

As we mentioned before, in open cavities the field dynamics is
governed by the complex poles of the $S$ matrix which correspond to
the cavity \emph{resonances}. In the rotating wave approximation they
coincide with the poles of $D^{-1}(\omega)$ and are obtained by
solving the equation $\det[D(\omega)]=0$. In Sec.~\ref{sec:Sm} we
showed that in the exact dynamics the system resonances coincide with
the poles of the projected resolvent
$\mathcal{G}_{\mathcal{QQ}}(\omega)$. It is then illustrative to write
down the relation between $\mathcal{G}_{\mathcal{QQ}}(\omega)$ and
$D(\omega)$. Considering Eq.~\eqref{eq:w-Leff} in the rotating wave
approximation one obtains,
\begin{equation} \label{eq:GRWA}
\mathcal{G}_{\mathcal{QQ}}^{-1}(\omega)_{\lambda \lambda'}
\simeq \frac{2\go}{c^2} D(\go)_{\gl \gl'},
\end{equation}
which shows that in this approximation the zeros of $D(\go)$ coincide
with the poles $\mathcal{G}_{\mathcal{QQ}}(\go)$. Furthermore, the
approximated $S$ matrix \eqref{eq:SRWA} follows from the exact $S$
matrix \eqref{eq:weids} upon replacing $\mathcal{G}_{\mathcal{QQ}}$ by
$(c^2/2\go)D^{-1}$. The importance of this result can be underline if
we make the Markov approximation. Then, the eigenvalues of the
non-Hermitian operator $\mathcal{H}$ correspond to the system
resonances, and hence the associated right (left) eigenvectors are
the ``true'' cavity modes, corresponding to solutions of the Maxwell's
equations with outgoing (incoming) boundary conditions. Although these
modes no longer form an orthogonal basis they can still be used to
represent the cavity field. Such expansion is discussed in
Sec.~\ref{sec:nonorth}.


\subsubsection{Linear Absorbing or Amplifying Medium}

The presence of an absorbing or amplifying medium within the cavity
leads to additional noise and modifies the input--output relation.
Phenomenologically, the interaction with linear media can be modeled
\cite{Gard00,Been98,Been99} by coupling the cavity modes to additional baths.
An absorbing medium is described by a thermal bath of harmonic
oscillators while an amplifying medium may be represented by a bath
of inverted harmonic oscillators at a negative temperature $-T$. The
total Hamiltonian is then given by
\begin{equation}
\label{eq:hmb}
H = H_{\textrm{SB}} + H_{\textrm{abs}}+ H_{\textrm{amp}},
\end{equation}
where $H_{\textrm{SB}}$ is the system--and--bath Hamiltonian
\eqref{eq:hrwa}, while $H_{\textrm{abs}}$ and $H_{\textrm{amp}}$
represent the absorbing and amplifying bath,
\begin{subequations}
\label{eq:hamblock}
\begin{align}
H_{\textrm{abs}} &=  \sum_l \int \! d\omega \:
\hbar\omega \, c^{\dagger}_{l}(\omega) c_{l}(\omega) + \hbar
\sum_{\lambda} \sum_l \int \! d\omega \left[ 
    \kappa_{\lambda l}(\omega) \, a^{\dagger}_{\lambda} c_{l}(\omega)
    + \textrm{H.c.}\right], \\ 
H_{\textrm{amp}} &= - \sum_k \int 
\! d\omega \: \hbar\omega \, d^{\dagger}_{k}(\omega) 
d_{k}(\omega) + \hbar \sum_\lambda \sum_k \int \! d\omega 
\left[ \sigma_{\lambda k}(\omega) \, a^{\dagger}_{\lambda} 
d_{k}(\omega) + \textrm{ H.c. }\right],
\end{align}
\end{subequations}
with coupling amplitudes $\kappa_{\gl l}$ and $\sigma_{\lambda k}$.
The operators $c_l$, $c_l^\dagger$ obey the canonical commutation
relations
\begin{equation}
\label{eq:cab}
[c_l(\omega),c_{l^{\prime}}^{\dagger}(\omega^\prime)] = \delta_{l
  l^\prime}\,\delta(\omega - \omega^\prime),
\end{equation}
and account for thermal emission within the absorbing medium. The
operators $d_k$ and $d_{k^\prime}^\dagger$ represent the amplifying
medium and have the commutation relations \cite{Gard00}
\begin{equation}
\label{eq:dam}
[d_k(\omega),d_{k^{\prime}}^{\dagger}(\omega^\prime)] = -\delta_{k
  k^\prime}\,\delta(\omega - \omega^\prime).
\end{equation}
As the Hamiltonian (\ref{eq:hmb}) gives rise to linear equations of
motion, we can compute the cavity output field using Fourier
transformation. The calculation proceeds along the lines of the
calculation presented in Sec.~\ref{sec:exactdyn}. The result is
\begin{equation}
\label{eq:slinm}
b^{\textrm{out}}(\omega) = S(\omega) b^{\textrm{in}}(\omega) + 
U(\omega) c^{\textrm{in}} (\omega) + 
  V(\omega) d^{\textrm{in}}(\omega),
\end{equation}
where $c^{\textrm{in}}$ and $d^{\textrm{in}}$ represent the input
noise of the absorbing and amplifying bath. Both are integrals over
bath operators at the initial time $t_0$,
\begin{subequations}
\label{eq:cin}
\begin{align}
c^{\textrm{in}}_l(t) & \equiv \frac{1}{2\pi}
\int \! d\omega \:\rme^{-\rmi\omega (t-t_0)} c_l(\omega,t_0), \\
d^{\textrm{in}}_k(t) & \equiv \frac{1}{2\pi} 
\int \! d\omega \:\rme^{-\rmi\omega (t-t_0)} d_k(\omega,t_0).
\end{align}
\end{subequations}
The $S$ matrix has the again the form \eqref{eq:SRWA} and the matrices
$U$ and $V$ are given by
\begin{subequations}
\begin{align}
U(\omega) &=  - 2\pi \rmi \mathcal{W}^{\dagger}(\omega)
    D^{-1}(\omega) \mathcal{K}(\omega),  \\
V(\omega) &=  - 2\pi \rmi \mathcal{W}^{\dagger}(\omega)
    D^{-1}(\omega) \Sigma(\omega) ,
\end{align}
\end{subequations}
where the $N\times L$ matrix $\mathcal{K}$ and the $N\times K$ matrix
$\Sigma$ comprise the coupling amplitudes $\kappa_{\lambda l}$ and
$\sigma_{\lambda k}$, respectively. In the presence of the absorbing
and amplifying baths, the elements of the $N \times N$ matrix
$D(\omega)$ have the form
\begin{equation}
\label{eq:dl}
D_{\lambda \lambda^{\prime}}(\omega)  =  (\omega - \omega_\lambda) 
\delta_{\lambda \lambda^{\prime}} + \mathcal{P} \int \!
d\omega^{\prime} \frac{ \Delta_{\lambda \lambda'} (\omega^\prime)} 
{\omega^{\prime} - \omega} + \rmi\pi \Delta_{\lambda
  \lambda^{\prime}} (\omega),
\end{equation}
where $\Delta$ is the matrix
\begin{eqnarray}
\Delta(\omega)  = \mathcal{W}(\omega)\mathcal{W}^{\dagger}
(\omega) + \mathcal{K}(\omega)\mathcal{K}^{\dagger}(\omega) - \Sigma
(\omega) \Sigma^{\dagger}(\omega) .
\end{eqnarray}
Using Eq.\ (\ref{eq:slinm}) and the commutation relations for the
output and input noise operators, one obtains the relation 
\begin{equation}
\label{eq:nonuni}
U U^{\dagger} - V V^{\dagger} = \openone - S S^{\dagger},  
\end{equation}
that was first derived by Beenakker \cite{Been98,Been99} using a
scattering approach to field quantization. We note that the matrix
$\openone -S S^\dagger$ is positive definite in an absorbing medium
($V =0$) and negative definite in an amplifying medium ($U=0$). The
relations (\ref{eq:slinm}) and (\ref{eq:nonuni}) are important as they
relate the intensity of the output field to the amplitudes of the
input field and the scattering matrix of the cavity. The statistical
properties of the scattering matrix are known from random matrix
theory. This allows \cite{Been98,Been99} to compute moments or even
the full distribution of the output field intensity from linear random
media.


\section{Master Equation}
\label{sec:MEQ} 

In this section we address the resonator dynamics within the
Schr\"odinger picture. We use the reduce density operator $\gr(\go)$
to represent the cavity field and derive a master equation to
describe its time evolution. We assume a time scale separation between
the damping time of the resonator field and the response time of the
bath, $1/\gG \ll \tau_{\textrm{bath}}$, which allow us to make the
Markov approximation to the Nakajima--Zwanzig equation \cite{Haak73},
and simplifies the derivation of the master equation.


\subsection{Derivation of the Master Equation}
\label{sec:dMeq}

We consider the cavity field dynamics in the rotating wave
approximation. Our starting point is the system-and-bath Hamiltonian
\eqref{eq:hrwa}, which we now write in terms of three contributions,
\begin{equation}
H_{\textrm{SB}} = H_{\textrm{S}} + H_{\textrm{B}} + H_{\textrm{Int}},
\end{equation}
where $H_{\text{S}}$ and $H_{\text{B}}$ are the isolated cavity and
bath Hamiltonian, respectively, and $H_{\textrm{Int}}$ describes the
interaction between the two subsystems; the explicit form of the
respective terms follow straightforwardly from Eq.~\eqref{eq:hrwa}.

The time evolution of the reduce density operator $\gr(t)$ of the
cavity field is given by the Nakajima--Zwanzig equation
\cite{Haak73}. This equation is obtained from the Liouville--von
Neumann equation $\dot{\gr}_{\text{tot}}(t) = -\rmi
[H_{\text{SB}},\gr_{\text{tot}}] / \hbar $ for the density operator of
the full electromagnetic field by projecting the dynamics onto the
cavity space, and tracing over the degrees of freedom of the outside
field using a reference bath state $\gr_{\text{B}}$. In the Markov
approximation, and taking as the reference state for the bath the
canonical density operator $\gr_{\text{B}} = \rme^{-H_{\text{B}}/kT}/
\tr_{\text{B}}\bigl\{\rme^{-H_{\text{B}}/kT}\bigr\} $, the master
equation for the reduce density operator is
\begin{equation}\label{eq:NZeq}
\dot{\gr}(t) = \Lambda \gr(t),
\end{equation}
with the time evolution generator operator
\begin{equation}
\label{eq:NZt}
\Lambda = -\frac{\rmi}{\hbar} [H_{\textrm{S}}, \;\cdot\; ] +
\int\limits_{0}^{+\infty} \! \rmd t' K(t').
\end{equation}
Here the integral kernel operator reads
\begin{equation}
K(t')\gr(t) = \frac{1}{\hbar^2} \tr_{\text{B}}\bigl\{[
H_{\text{Int}}(t),[H_{\text{Int}}(t'), \gr_{\text{B}} \otimes \gr(t)] ]
\bigr\}, 
\end{equation}
where $\tr_{\text{B}} \{\,\cdot\,\}$ stands for the trace over the
bath degrees of freedom, and the time dependent interaction Hamiltonian
$H_{\text{int}}(t)$ is given by
\begin{equation} \label{eq:Hint}
H_{\text{Int}}(t) = \exp \bigg(-\frac{\rmi}{\hbar} (H_{\textrm{S}} +
  H_{\textrm{B}}) t \bigg) H_{\text{Int}} \exp
\bigg(\frac{\rmi}{\hbar} (H_{\textrm{S}} + H_{\textrm{B}}) t \bigg).
\end{equation}

The evaluation of \eqref{eq:NZeq} using
Eqs.~\eqref{eq:NZt}--\eqref{eq:Hint} is a lengthy but simple
exercise. The final result is derived in Appendix \ref{app:2C} and
includes not only coupling amplitudes that depend on frequency but
also the so-called frequency-shift terms. A master equation equivalent
to the Langevin equation \eqref{eq:aLEQ} is recovered neglecting the
frequency dependence of the coupling amplitudes as well as the
frequency-shift terms. The master equation has then the form
\begin{multline}
\label{eq:ME1}
\dot{\gr} = -\rmi \sum_{\gl} \go_{\gl} [a_{\gl}^{\dag} a_{\gl}, \gr] +
(1+n_{\text{th}}) \sum_{\gl \gl'} \gG_{\gl \gl'} \big( [a_{\gl'},\gr
a_{\gl}^{\dag}] + [a_{\gl'}\gr,a_{\gl}^{\dag}] \big)\\ + n_{\text{th}}
\sum_{\gl \gl'} \gG_{\gl \gl'} \big( [a_{\gl}^{\dag}, \gr a_{\gl'}] +
[a_{\gl}^{\dag} \gr, a_{\gl'}] \big).
\end{multline}
This equation generalizes the familiar quantum optical master equation
for a single damped harmonic oscillator to many oscillators coupled by
the (off-diagonal elements of the) damping matrix $\gG_{\gl\gl'}$.
The latter coupling is important when the damping is strong enough to
cause spectral overlap of modes. The first double sum, proportional to
$1+n_{\text{th}}$ ($n_{\text{th}}=[\exp(\hbar\go/kT)-1]^{-1}$),
describes spontaneous and induced emission of photons towards the
outside while the second double sum, proportional only to
$n_{\text{th}}$, describes absorption from the outside; that
interpretation is easily checked by employing the Fock representation,
i.e., the representation in terms of eigenstates of the photon number
operators $a_{\gl}^{\dag} a_{\gl}$.

Systematic and stochastic forces are not as clearly separated here as
in the Langevin equation. In order to bring about such distinction
here as well, we may imagine the density operator $\rho(t)$ at time
$t$ anti-normally ordered in the annihilation and creation operators
(all $a$'s to the left of all $a^{\dag}$'s); further, we rearrange the
commutators in the right-hand side of the master equation \eqref{eq:ME1}
such that $\dot{\gr}$ becomes anti--normally ordered provided
$\gr(t)$ is,
\begin{multline}
\label{eq:ME2}
\dot{\gr} = -\rmi \sum_{\gl} \go_{\gl} \big([a_{\gl},\gr
  a_{\gl}^{\dag}] - [a_{\gl}\gr,a_{\gl}^{\dag}] \big) + \sum_{\gl
  \gl'} \gG_{\gl \gl'} \big([a_{\gl'},\gr a_{\gl}^{\dag}] +
    [a_{\gl'}\gr, a_{\gl}^{\dag}] \big) \\ 
+ 2n_{\text{th}} \sum_{\gl
      \gl'} \gG_{\gl \gl'} 
[[a_{\gl'},\gr],a_{\gl}^{\dag}]\, . 
\end{multline}
The latter form of the master equation preserves anti--normal ordering
of $\rho(t)$ at all times. Moreover, we have now separated reversible
drift terms (proportional to the frequencies $\go_{\gl}$), irreversible
drift terms ($\propto \gG_{\gl\gl'}$) not involving the thermal
number of quanta $n_{\text{th}}$, and noise generated diffusion terms
($\propto n_{\text{th}}$); the latter interpretation will become obvious
in the next section where a representation based on coherent states
will be employed.


\section{Fokker--Planck Equation}
\label{sec:FPEQ}

A powerful technique to find solutions of quantum Markov processes are
the so-called phase space methods
\cite{Gard00,Mand95,Wall94} in which the system operators are degraded
to $c$-numbers and the system density operator is represented by a
quasiprobability function. For Markov processes the $c$-number
equation of motion found for the quasiprobability function is usually
of the Fokker--Planck form and standard analytical methods may be used
to solve it.

In this section we employ the Glauber--Sudarshan $P$-function to describe
the cavity field. Starting from the master equation derived in the
previous section, we obtain a Fokker--Plank equation for the time
evolution of the $P$-function and find its steady state solution. In
addition, we establish the equivalence of this Fokker--Plank equation
with the Langevin equation \eqref{eq:aLEQ}.

\subsection{Derivation of the Fokker--Planck Equation}

If we consistently stick to antinormal ordering of $\rho(t)$ we may
write the commutators in the master equation \eqref{eq:ME2} as
differential operators \cite{Gard00} as $[a,(\cdot)]\to
(\partial/\partial a^{\dag})(\cdot)$ and $[(\cdot),a^{\dag}]\to
(\partial/\partial a)(\cdot)$. We may then just as well degrade
all creation and annihilation operators to complex c-number variables
\begin{equation}
a_{\gl}\to\ga_{\gl},\qquad a_{\gl}^{\dag}\to\ga_{\gl}^{*} ,
\end{equation}
and introduce the $P$-function to represent the density operator, 
\begin{equation}
\gr(t) \to P(\{\ga,\ga^{*}\},t).
\end{equation} 
The $P$-function has as its moments expectation values of normally
ordered observables,
\begin{equation}
\biggl\langle \prod_{\gl}(a_{\gl}^{\dag})^{m_{\gl}} \prod_{\gl'}
(a_{\gl'})^{n_{\gl'}} \biggr\rangle (t) = \int \prod_\gl
\rmd^2\ga_{\gl}\; (\ga_{\gl}^{*})^{m_{\gl}} (\ga_{\gl})^{n_{\gl}}
P(\{\ga,\ga^{*}\},t).
\end{equation} 
Here the differential volume element is to be interpreted according
to
\begin{equation} 
\rmd^2\ga_{\gl} = \rmd\text{Re}(\ga_{\gl}) \rmd\text{Im}(\ga_{\gl}) =
\rmd\ga_{\gl} \rmd\ga_{\gl}^{*} / 2\rmi\,.
\end{equation}
The master equation \eqref{eq:ME2} then becomes the Fokker--Planck
equation (with the shorthand $\partial / \partial\ga_{\gl} \equiv
\partial_{\gl}$, $\partial/\partial \ga_{\gl}^{*}\equiv
\partial_{\gl}^{*}$)
\begin{equation}\label{eq:FPE}
\dot{P}=\biggl[-\rmi \sum_{\gl} \go_{\gl} (\partial_{\gl}^{*}
\ga_{\gl}^{*} - \partial_{\gl} \ga_{\gl} ) + \sum_{\gl\gl'}
\gG_{\gl\gl'} (\partial_{\gl}^{*} \ga_{\gl'}^{*} + \partial_{\gl'}
\ga_{\gl} 
 + 2 n_{\text{th}} \partial_{\gl'} \partial_{\gl}^{*} ) \bigg] P.
\end{equation}
At this point the interpretation of the various terms in the master
equation \eqref{eq:ME2} given earlier becomes obvious.

The $c$-number Langevin equation equivalent to \eqref{eq:FPE} is
\cite{Gard00,Risk89} 
\begin{equation}\label{eq:cLE} 
\dot{\ga}_{\gl}(t) = -\frac{\rmi}{\hbar} \sum_{\gl'}
\mathcal{H}_{\gl\gl'} \ga_{\gl'}(t) + \varphi_{\gl}(t), 
\end{equation}
where the random forces $\varphi (t)$ have a Gaussian statistics and a
white spectrum according to
\begin{equation}
\la \varphi_{\gl}^{*}(t) \varphi_{\gl'}(t') \ra = 2 n_{\text{th}}
\gG_{\gl' \gl}\, \delta(t-t'),
\end{equation}
while holding on to $\la \varphi(t)\ra = \la \varphi_{\gl}(t)
\varphi_{\gl'}(t') \ra = 0$. We notice that Eq.~\eqref{eq:cLE} is
precisely the $c$-number equation that we would have obtained by
degrading the operators in the quantum Langevin equation
\eqref{eq:aLEQ} to $c$-numbers and making the correspondence $F(t) \to
\varphi(t)$ for the noise operators. Thus, the Fokker--Plank equation
\eqref{eq:FPE} is equivalent to the quantum Langevin equation
\eqref{eq:aLEQ}. 


\subsection{Stationary Solution of the Master Equation}

The linearity of the drift coefficients and the constancy of the
diffusion tensor in the Fokker-Planck equation \eqref{eq:FPE}
indicates that we are facing a stochastic process of the
Ornstein--Uhlenbeck type. Then the general time dependent solution of
\eqref{eq:FPE} can be constructed \cite{Risk89}. The Gaussian
distribution of the noise together with the linear evolution equation
\eqref{eq:cLE} imply that the stochastic variables must be Gaussian
distributed. In particular, the stationary $P$-function is immediately
checked to be 
\begin{equation}
\bar{P} (\{ \alpha, \alpha^{*}\}) = \prod_{\lambda}\frac{1}{\pi
n_{\text{th}}} \exp( -\alpha^{*}_{\lambda} \alpha_{\lambda}
/n_{\text{th}}).
\end{equation} 
The dissipative coupling of the system modes is no longer visible in
the stationary state; rather, we encounter the thermal equilibrium
state one would also find in the absence of spectral overlap.



\section{Cavity Resonances Representation Using Non\-orthogonal Modes}
\label{sec:nonorth}

So far we have described the resonator field dynamics in terms of
eigenmodes of a conveniently chosen \emph{closed} system, for which
the associated mode operators $\{ a_{\gl}, a_{\gl}^{\dag} \}$ obey
canonical commutation relations. It turns out, that the cavity
field dynamics is generated by the non-Hermitian effective Hamiltonian
$\mathcal{H}$ defined in Eq.~\eqref{eq:Heff0}. Hence, one should
expect the equations of motion for the cavity field to look simpler
provided the field is described in terms of the operator
$\mathcal{H}$. The complex eigenvalues of $\mathcal{H}$ are the
\emph{open} cavity resonances in the presence of the coupling to the bath,
while its right and left eigenvectors form a biorthogonal set of modes
representing the resonances of the open cavity. If one assumes
this set of modes to be complete, one can employ it to represent the
intracavity field. In what follows we show how such a representation
can be implemented by introducing a new set of operators associated
with the cavity resonance modes. We will see that the field dynamics
looks simpler in terms of these new operators, but that the new operators obey
noncanonical commutation rules.

Assuming the eigenvalues of $\cal{H}$ are nondegenerate, $\cal{H}$ can
be diagonalized by a similarity transformation
\begin{equation}
\label{eq:THT}
\go = \frac{1}{\hbar} T^{-1} \mathcal{H} T.
\end{equation}
The diagonal matrix $\go$ comprises the complex eigenvalues of
$\mathcal{H}/ \hbar $ with
\begin{equation}
\label{eq:cplx-f}
\go_{k} = \gO_{k}-i\gk_{k},
\end{equation}
where $\gO_{k}$ is the resonance frequency and $\gk_{k}$ is the
resonance broadening. The entries of the matrix $T$ are the
overlapping integrals
\begin{equation}
\label{eq:Tdef}
T_{\gl k} = \la \vect{u}_{\gl} | \vect{R}_{k} \ra \text{ and }
T_{k\gl}^{-1} = \la \vect{L}_{k} | \vect{u}_{\gl} \ra .
\end{equation}
Here ${\vect{u}_{\gl}}$ are the modes of the closed system while
${\vect{R}_{k}}$ and ${\vect{L}_{k}}$ are, respectively, the set of
right and left eigenmodes of $\mathcal{H}$. They satisfy $\mathcal{H}
| \vect{R}_{k} \ra = \hbar\go_{k} | \vect{R}_{k} \ra$ and
$\mathcal{H}^{\dag} | \vect{L}_{k} \ra = \hbar\go_{k}^{*} | \vect{L}_{k}
\ra$, from which it follows
\begin{equation}
\la \vect{L}_{l} | \vect{R}_{k} \ra = \delta_{lk}.
\end{equation}
We chose this particular normalization because it implies the completeness
relations 
\begin{equation}
\sum_{k} |\vect{L}_{k} \ra \la \vect{R}_{k}| = \sum_{k} |\vect{R}_{k}
\ra \la \vect{L}_{k}| = \openone. 
\end{equation}
The most prominent feature of this nonorthogonal basis is that in
general $\la \vect{R}_{l} |
\vect{R}_{k} \ra\not= 0$ and $\la \vect{L}_{l} | \vect{L}_{k} \ra\not=
0$. A measure of the nonorthogonality is given by the mode
Petermann factor $K$, defined by
\begin{equation}
\label{eq:Pete-fac}
K_{k} = \la \vect{R}_{k} | \vect{R}_{k} \ra \la \vect{L}_{k} |
\vect{L}_{k} \ra \,.
\end{equation}
Using Schwartz inequality one can easily show that $K \ge 1$ for
all modes; the equality only holds for orthogonal modes. For general
non-Hermitian matrices, eigenvalues and eigenvectors are
correlated. In particular, large values of $K$ are expected for almost
degenerated modes. The nonorthogonality of the cavity resonance modes
has physical observable consequences. The linewidth in lasers based on
nonorthogonal modes is enhanced by a factor $K$
\cite{Pete79,Sieg89a,Sieg89b} with respect to the value predicted by the
Schawlow--Townes formula. So far values up to $K \approx 500$ have been
measured for transverse modes in unstable cavities
\cite{Chen96,Eijk96}; recently, investigating the location of
resonance degeneracies, Berry \cite{Berr03} has propose that unstable
cavity configurations exist for which $K$ takes arbitrarily large
values.

Associated with the $k$th right and left eigenmodes of $\cal{H}$,
respectively, are the operators
\begin{subequations}
\label{eq:de-def}
\begin{align}
\label{eq:d-def}
d_{k} &= \sum_{\gl} T_{k\gl}^{-1} a_{\gl}, \\
\label{eq:e-def}
e_{k}^{\dag} &= \sum_{\gl} a_{\gl}^{\dag} T_{\gl k}.
\end{align}
\end{subequations}
The representation of the fields inside the cavity in terms of the new
operators and the modes of the open cavity is obtained in the rotating
wave approximation from the field expansion \eqref{eq:afin} and the
definitions \eqref{eq:de-def},
\begin{subequations}
\begin{align}
\vect{A}(\vect{r},t) &= \sum_{k}
\left(\frac{\hbar}{2\gO_{k}}\right)^{1/2} [ d_{k}
  \vect{R}_{k}(\vect{r}) + e_{k}^{\dag} \vect{L}_{k}^{*}(\vect{r}) ],
\\ 
\vect{E}(\vect{r},t) &= \rmi\sum_{k}
\left(\frac{\hbar\gO_{k}}{2}\right)^{1/2} [ d_{k}
  \vect{R}_{k}(\vect{r}) - e_{k}^{\dag} \vect{L}_{k}^{*}(\vect{r}) ],
\end{align}
\end{subequations}
where we used that $\go_{k} \sim \gO_{k}$ and $\gO_{k} \gg
\gk_{k}$. The new operators do not represent bosonic creation and
annihilation operators. Instead they satisfy the following commutation relations
\begin{equation}
\begin{split}
[ d_{l},d_{k}] &= 0, \qquad [ e_{l},e_{k}] = 0, \qquad [ d_{l},e_{k}]
= 0, \\ 
[ d_{l},d_{k}^{\dagger}] &= (T^{-1}T^{-1\dag})_{lk} = \la \vect{L}_{l}
| \vect{L}_{k} \ra, \\
[ e_{l},e_{k}^{\dagger}] &= (T^{\dag}T)_{lk} = \la \vect{R}_{l}
| \vect{R}_{k} \ra, \\
[ d_{l},e_{k}^{\dagger}] &= \delta_{lk},
\end{split}
\end{equation}
which are easily checked using the definitions \eqref{eq:de-def} and
the commutation relations \eqref{eq:a-comm} for the operators $\{
a_{\gl}, a_{\gl}^{\dag} \}$. The new operators $\{ d_{k},
e_{k}^{\dag} \}$ coincide with the set of operators introduced by
Lamprecht and Ritsch \cite{Lamp99,Lamp02} for the non-orthogonal modes
of unstable cavities. These authors did not derive dynamical equations
for the novel operators. That dynamics will be addressed below.


\subsection{Langevin Equations}

The Heisenberg equations of motion for the operators $d_{k}$ are
obtained by differentiating definition \eqref{eq:d-def} with
respect to $t$ and using the Heisenberg equations of motion
\eqref{eq:aLEQ} for the operators $a_{\gl}$. The result is 
\begin{equation}
\label{eq:d-eqmot}
\dot{d}_{k}(t) = -\rmi \go_{k} d_{k}(t) + F_{k}(t).
\end{equation}
The fact that these equations of motion decouple from each other is hardly
surprising since the operators $d_{k}$ were introduced for that
purpose. The noise forces are given by
\begin{equation}
\label{eq:nonf-in}
F_{k}(t) = \sum_{\gl} T_{k\gl}^{-1} F_{\gl}(t) = -2\pi\rmi
\sum_{m}(T^{-1}\mathcal{W})_{km} b_{m}^{\text{in}}(t),
\end{equation}
displaying their dependence on the input field. Here the coupling
amplitudes $\mathcal{W}_{\gl m}$ form an $N \times M$ matrix
$\mathcal{W}$ ($M=M(\go)$ is the number of open channels at frequency
$\go$ and the number of cavity modes $N$ is taken to infinity at the
end of the calculations) and $T$ is the $N
\times N$ linear transformation defined in Eq.~\eqref{eq:Tdef}.  The
Langevin equations of motion follow from the Heisenberg equations of
motion \eqref{eq:d-eqmot} by considering noise forces with a white
spectrum, such that $\la F_{k}(t) \ra = 0$, and the correlation
functions different from zero are given by
\begin{subequations}
\label{eq:cffs-non}
\begin{align}
\label{eq:cffn-non}
\langle F_{l}^{\dagger}(t) F_{k}(t^\prime) \rangle & = \rmi \langle
L_{l}|L_{k} \rangle (\go_{l} - \go_{k}^{*}) n_{\text{th}} \,
\delta(t-t^\prime), \\
\label{eq:cffan-non}
\langle F_{k}(t) F_{l}^{\dagger}(t^\prime) \rangle & = \rmi \la
L_{l}|L_{k} \ra (\go_{l} - \go_{k}^{*}) ( n_{\text{th}} + 1) \,
\delta(t-t^\prime),
\end{align}
\end{subequations}
with the thermal number of photons $n_{\text{th}} = [\exp(\hbar\bar\go
  / kT) -1 ]^{-1}$. Thus, although the deterministic part of the
dynamics does not couple the open cavity modes, the noise terms for
different modes are correlated.

\subsubsection{Input--Output Relation}

The input-output relation can be expressed in terms of the new set of
operator. In the time domain, one obtains the relation
\begin{equation}
\label{eq:in-outnont}
b_{m}^{\text{out}}(t) - b_{m}^{\text{in}}(t) = -\rmi \sum_{k}
(\mathcal{W}^{\dag} T)_{mk} d_{k}(t')\,.
\end{equation}
After Fourier transformation the above equations yield 
\begin{eqnarray}
\label{eq:in-outnonw}
b_{m}^{\textrm{out}} (\omega) - b_{m}^{\textrm{in}} (\omega) = - \rmi
\sum_{k}(\mathcal{W}^{\dagger}T)_{mk} d_{k}(\omega).
\end{eqnarray}
These conditions are obtained from the input--output
relations \eqref{eq:in-outrwat} and
\eqref{eq:in-outrwaw} employing the inverse of the definition
\eqref{eq:d-def}.


\subsection{Master Equation}

The operators $\{ d, e^{\dag} \}$ may be employed to express
the cavity field dynamics in the Schr\"odinger picture. In terms of
these operators the master equation \eqref{eq:ME2} takes the form
\begin{equation}
\label{eq:meq-d}
\dot{\gr} = -\rmi \sum_{k} \Bigl\{ \go_{k}\bigl( [e_{k}^{\dag},d_{k}
  \gr] - n_{\text{th}}[[d_{k},\gr],e_{k}^{\dag}] \bigr) - \go_{k}^{*}
\bigl( [\gr d_{k}^{\dag},e_{k}] -
n_{\text{th}}[[e_{k},\gr],d_{k}^{\dag}] \bigr) \Bigr\}.
\end{equation}
Here we used the relations $\go_{\gl}
\delta_{\gl \gl'} = \frac{1}{2\hbar} (\mathcal{H}_{\gl \gl'} +
\mathcal{H}_{\gl \gl'}^{\dag} )$ and $\gG_{\gl \gl'} =
\frac{i}{2\hbar} (\mathcal{H}_{\gl \gl'} - \mathcal{H}_{\gl
  \gl'}^{\dag} )$, together with the representation \eqref{eq:THT} of
$\mathcal{H}$, and the definitions \eqref{eq:de-def} of the operators
$\{ d,e^{\dag} \}$. We notice that in the limit $kT \ll
\hbar\bar{\go}$ the master equation \eqref{eq:meq-d} can be recast as
\begin{equation}
\dot{\gr} = -\frac{\rmi}{\hbar} (\mathcal{H}\gr -
\gr\mathcal{H}^{\dag}) + \rmi \sum_{kl} \la \vect{R}_{l} |
\vect{R}_{k}  \ra (\go_{k} - \go_{l}^{*}) d_{k}\gr d_{l}^{\dag},
\end{equation}
with $\mathcal{H} = \sum_{k} \hbar\go_{l} e_{k}^{\dag}d_{k}$. This is
the master equation proposed by Lamprecht and Ritsch
\cite{Lamp02,Lamp02a}. Our derivation provides a microscopic basis
for this equation and proves the equivalence to the master equations
\eqref{eq:ME1} and \eqref{eq:ME2} in the limit $kT \ll \hbar\bar{\go}$.