text:
# Field quantization for open optical cavities

## I Introduction

In the past few years advances in microstructuring techniques have made it possible to manufacture novel mirrorless cavities known as random lasers. Feedback in these lasers is provided by the random scattering of light in a medium with spatially fluctuating refractive index. Random lasing has been demonstrated in semiconducting clusters and films Cao99  as well as in solutions of $`\text{TiO}_2`$ nanoparticles Cao00 ; Lin02  and in polymer systems Fro98 . Random lasers typically have dimensions of several $`\mu m`$ and may prove useful in microtechnology applications Wie00 . Fundamentally, interest in these systems derives from the complexity of the interplay of amplification, disorder and random scattering of light.

The theoretical investigation of disordered media has a long history She90 ; She95 , but most studies are restricted to passive systems in the regime of classical optics. Far reaching similarities have been revealed between the propagation of light through (passive) random dielectrics and the transport of electrons through disordered solids. Well-known manifestations of this similarity are the coherent backscattering of light Alb85  versus the weak–localization for electrons Imr97  and the strong localization of microwaves Gen91  or light Wie97  versus the Anderson localization of electrons And58 . Conceptually, this similarity may be understood from the nonlinear supersymmetric $`\sigma `$ model which proved to be the underlying field theory both for disordered solids Efe97  and for classical optical waves in random dielectrics Ela98 .

In this paper we address the quantum properties of the electromagnetic field in the presence of a random medium. A quantum treatment of light is required when one wants to compute the spectrum, the linewidth or the photon statistics of the output radiation. We study both dielectric media with a spatially nonuniform dielectric constant $`ϵ(𝐫)`$ and optical cavities defined by mirrors of arbitrary shape. These systems may be coupled to an arbitrary number of escape channels. A special case are cavities that, e.g. due to partly reflecting mirrors, are only weakly coupled to the external radiation field. Such cavities are used in standard lasers, and have been extensively studied in laser theory. The standard quantum theory of the laser Hak84  starts from the expansion of the electromagnetic field in terms of closed resonator modes. Escape from the resonator gives rise to a small perturbation of the closed resonator dynamics. To leading order in a perturbation theory in terms of the coupling to the outside world, the resonator modes may well be approximated by modes of an entirely closed system. This perturbation theory breaks down in cavities with a weak confinement of light. For such systems, the field quantization of standard resonator theory must be replaced by a generalization suitable for open resonators. We present and discuss such a generalized quantization technique in this paper, and demonstrate that it is well suited for analytical investigations of random media.

Our solution of the quantization problem is most relevant for “open” cavities like random lasers, but is significant for a more fundamental reason. Our technique allows for a rigorous treatment of a fundamental problem of quantum optics, concerning the damping of radiation due to leakage out of a cavity. Conventionally, this problem has been described by a system–and–bath Hamiltonian Sen59 ; Gar00 : The system is modeled by a discrete set of independent quantized harmonic oscillators associated with the normal modes of a completely isolated cavity. A different, continuous set of oscillators represents the bath while a coupling between the system and bath oscillators gives rise to damping. In spite of its intuitive appeal and its success in applications, this model of damping has been criticized. The formulations Sen59 ; Gar00  lack a microscopic justification as the coupling constants only enter as phenomenological parameters. Moreover, it has been argued Bar88 ; Dut00a  that the approach is restricted to good resonators with damping rates much smaller than the mode frequency spacing. We show below that such pessimism is inappropriate: We rigorously derive the system–and–bath Hamiltonian starting from the Maxwell equations. We employ a technique from nuclear and condensed matter physics known as the Feshbach projector formalism. Expressions are obtained for the coupling amplitudes. We find that the field Hamiltonian generally includes both resonant and non–resonant terms. When both are kept, the resulting dynamics correctly describes damping and even overdamping (when the damping rate becomes of the order of the frequency).

We note that there are alternative approaches to the field quantization in open resonators. These are either based on mode expansions of the electromagnetic field or they altogether abandon the notion of cavity modes and directly quantize the wave equation using Green functions. Approaches based on mode expansions include the modes–of–the–universe or true–mode approach (the field is expanded in terms of eigenstates of the total system comprising the resonator and the bath) Lan73 ; Gla91 ; Kno87 , expansions in terms of so–called quasi–modes Dal99  or in terms of non–orthogonal (or Fox–Li) modes You98 ; Lam99 ; Lam02 ; Dut00b ; Bro01  (a characterization of these approaches can be found in Ref. Dut00b ). Alternative approaches based on Green functions were developed by Grunder and Welsch Gru96  and by Loudon and co–workers Art97 . So far, to our knowledge, none of the first approaches has been applied to random media. In contrast, the second (Green function) approach has been successfully used by Beenakker and coworkers Bee98 ; Mis00  to establish a general relation between the emission of linear optical media and the underlying scattering matrix. Combining this relation with the statistical properties of the scattering matrix in disordered media, Beenakker derived the full photocount distribution of the radiation emitted from linear random media. Unfortunately, this random scattering approach suffers from a fundamental limitation: It is restricted to linear media and therefore cannot describe lasers above the lasing threshold. No such limitations hold for the approach presented in this paper as this approach is based on an Hamiltonian: atom–field interactions can be included in a standard way.

The outline of the paper is a follows. In Sec. II we describe the field quantization. We start from a global expansion of the field in terms of the eigenmodes of the Maxwell equations. The system–and–bath Hamiltonian is then derived using Feshbach’s projector technique. Our derivation is valid for spatially inhomogeneous dielectrics and fields of arbitrary polarization. This is a substantial generalization of our previous work Hac01  which was restricted to scalar fields and homogeneous dielectric media. In Sec. III we investigate the field dynamics. Using the input–output theory of Gardiner we calculate the output field in the presence of an arbitrary number of escape channels. We derive the equations–of–motion of the internal cavity modes, and show that the dynamics of these modes is coupled by damping and noise due to the external radiation field. As a final illustration of our technique we compute in Sec. IV the decay rate of a single two–level atom in a cavity of arbitrary shape. We conclude by discussing possible further applications of our technique, most notably the application to random lasers.

## II Field quantization

### II.1 Normal modes

We consider a three dimensional linear dielectric medium characterized by a scalar dielectric constant $`ϵ(𝐫)`$ that depends explicitly on position. We assume that the dielectric constant is real and frequency independent. The dielectric is surrounded by free space. Cavities defined by (ideal) mirrors are a special case with $`ϵ(𝐫)1`$ and appropriate boundary conditions on the mirrors. The electromagnetic field for the total system, comprising the resonator and the external radiation field, may be quantized using the exact eigenmodes of Maxwell’s equations. This so–called modes–of–the–universe approach Lan73 ; Gla91  serves as a starting point for the derivation of the system–and–bath Hamiltonian in Sec. II.2; we therefore summarize the main steps of this approach below.

It is convenient to formulate the quantization procedure in terms of the vector potential $`𝐀`$ and the scalar potential $`\varphi `$. We work in the Coulomb gauge which, in the absence of sources, corresponds to the choice $`\varphi =0`$ and the generalized transversality condition $`[ϵ(𝐫)𝐀]=0`$. The magnetic and electric field then follow from the potentials via the familiar relations

$$𝐄=\frac{1}{c}\frac{𝐀}{t},𝐁=\times 𝐀.$$
(1)
The electromagnetic Hamiltonian of the problem is given by

$`H={\displaystyle \frac{1}{2}}{\displaystyle 𝑑𝐫\left[\frac{c^2𝚷(𝐫,t)^2}{ϵ(𝐫)}+(\times 𝐀(𝐫,t))^2\right]},`$ (2)
where $`𝚷(𝐫)=ϵ(𝐫)\dot{𝐀}(𝐫)/c^2`$ is the canonical momentum field. The quantization of the fields may be achieved by imposing a suitable commutation relation between $`𝐀(𝐫,t)`$ and $`𝚷(𝐫,t)`$. An alternative but equivalent procedure is to expand the fields in a complete set of mode functions and to impose canonical commutation relations for the expansion coefficients. We follow the second procedure here, and expand the vector potential in terms of the exact eigenmodes $`𝐟_m(\omega ,𝐫)`$, defined as solutions of the wave equation

$`\times \left[\times 𝐟_m(\omega ,𝐫)\right]{\displaystyle \frac{ϵ(𝐫)\omega ^2}{c^2}}𝐟_m(\omega ,𝐫)=\mathrm{𝟎}.`$ (3)
The solutions automatically satisfy the transversality condition $`\left[ϵ(𝐫)𝐟_m(\omega ,𝐫)\right]=0`$. The eigenmodes are labeled by the continuous frequency $`\omega `$ and a discrete index $`m`$ that specifies the asymptotic boundary conditions far away from the dielectric (including the polarization). We consider asymptotic conditions corresponding to a scattering problem with incoming and outgoing waves. Then $`𝐟_m(\omega ,𝐫)`$ represents a solution with an incoming wave in channel $`m`$ and only outgoing waves in all other scattering channels. The definition of the channels depends on the problem at hand: For a dielectric coupled to free space, one may expand the asymptotic solutions in terms of angular momentum states. Then $`m`$ corresponds to an angular momentum quantum number. On the other hand, for a dielectric connected to external waveguides, $`m`$ may represent a transverse mode index. It is convenient to combine the solutions associated with the different channels to an $`M`$–component vector $`𝐟(\omega ,𝐫)`$ where $`M`$ is the total number of open channels at frequency $`\omega `$. The field expansions then take the form

$`𝐀(𝐫,t)`$ $`=`$ $`c{\displaystyle 𝑑\omega q(\omega ,t)𝐟(\omega ,𝐫)},`$ (4a)
$`𝚷(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{1}{c}}{\displaystyle 𝑑\omega 𝐟^{}(\omega ,𝐫)p(\omega ,t)},`$ (4b)
where the expansion coefficients $`q(\omega )`$ and $`p(\omega )`$ form $`M`$–component row and column vectors, respectively. The components $`q_m(\omega ,t)`$, $`p_m(\omega ,t)`$ are time–dependent variables.

The eigenmodes fulfill an orthonormality relation. This follows from the observation that the substitution

$$𝐟(\omega ,𝐫)=\frac{1}{\sqrt{ϵ(𝐫)}}\mathit{\varphi }(\omega ,𝐫)$$
(5)
transforms Eq. (3) into the eigenvalue problem

$`L\mathit{\varphi }(\omega ,𝐫){\displaystyle \frac{1}{\sqrt{ϵ(𝐫)}}}\times \left[\times {\displaystyle \frac{\mathit{\varphi }(\omega ,𝐫)}{\sqrt{ϵ(𝐫)}}}\right]={\displaystyle \frac{\omega ^2}{c^2}}\mathit{\varphi }(\omega ,𝐫),`$ (6)
for the Hermitean differential operator $`L`$. Choosing an orthonormal set of basis functions $`\mathit{\varphi }_m(\omega ,𝐫)`$, it follows that the associated mode functions $`𝐟_m(\omega )`$ satisfy the orthonormality condition

$$𝑑𝐫ϵ(𝐫)𝐟_m^{}(\omega ,𝐫)𝐟_m^{}(\omega ^{},𝐫)=\delta _{mm^{}}\delta (\omega \omega ^{}).$$
(7)
The functions $`\mathit{\varphi }_m(\omega )`$ form a complete set in the subspace of $`L^2`$ functions defined by the transversality condition

$$\left[\sqrt{ϵ(𝐫)}\mathit{\varphi }_m(𝐫)\right]=0.$$
(8)
The associated mode functions $`𝐟_m(\omega )`$ form a complete set in the space of transverse functions Gla91 .

Inasmuch as the fields are real, the vector potential and its canonical momentum fulfill the relations $`𝐀=𝐀^{}`$ and $`𝚷=𝚷^{}`$. Together with Eq. (4), this implies

$`q_m(\omega )`$ $`=`$ $`{\displaystyle \underset{m^{}}{}}{\displaystyle 𝑑\omega ^{}_{mm^{}}^{}(\omega ,\omega ^{})q_m^{}^{}(\omega ^{})},`$ (9a)
$`p_m^{}(\omega )`$ $`=`$ $`{\displaystyle \underset{m^{}}{}}{\displaystyle 𝑑\omega ^{}_{mm^{}}^{}(\omega ,\omega ^{})p_m^{}(\omega ^{})},`$ (9b)
where $``$ has the matrix elements

$`_{mm^{}}(\omega ,\omega ^{})={\displaystyle 𝑑𝐫ϵ(𝐫)𝐟_m(\omega ,𝐫)𝐟_m^{}(\omega ^{},𝐫)}.`$ (10)
We note that $``$ is unitary and symmetric foot1 . Moreover, $``$ only couples degenerate modes, as modes with different frequencies are orthogonal, $`(\omega ,\omega ^{})\delta (\omega \omega ^{})`$.

Substituting the field expansions (4) into the Hamiltonian (2), and using Eq. (9) and the properties of $``$, one obtains the Hamiltonian in terms of the variables $`q`$ and $`p`$,

$`H={\displaystyle \frac{1}{2}}{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \left[p_m^{}(\omega )p_m(\omega )+\omega ^2q_m^{}(\omega )q_m(\omega )\right]}.`$ (11)
Quantization is now achieved by promoting the variables $`q(\omega )`$ and $`p(\omega )`$ to operators. The Heisenberg equations of motion for $`q(\omega )`$ and $`p(\omega )`$ lead to Maxwell’s equations, provided we impose the equal time commutation relations

$`[q_m(\omega ),q_n(\omega ^{})]`$ $`=`$ $`[q_m(\omega ),q_n^{}(\omega ^{})]=0,`$
$`[p_m(\omega ),p_n(\omega ^{})]`$ $`=`$ $`[p_m(\omega ),p_n^{}(\omega ^{})]=0,`$ (12)
$`[q_m(\omega ),p_n(\omega ^{})]`$ $`=`$ $`i\mathrm{}\delta _{mn}\delta (\omega \omega ^{}).`$
Further use of Eq. (9) gives the remaining commutation relation

$`[q_m(\omega ),p_n^{}(\omega ^{})]`$ $`=`$ $`i\mathrm{}_{mn}^{}(\omega ,\omega ^{}).`$ (13)
Combining the field expansions (4) with the orthogonality condition (7), one can show that these commutation relations imply canonical commutation relations for the vector potential and the canonical momentum field.

The last step in the quantization procedure is to express the operators $`q(\omega )`$ and $`p(\omega )`$ in terms of creation and annihilation operators,

$`q_m(\omega )`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}}{2\omega }}\right]^{\frac{1}{2}}[A_m(\omega )`$ (14a)
$`+{\displaystyle \underset{m^{}}{}}{\displaystyle }d\omega ^{}_{mm^{}}^{}(\omega ,\omega ^{})A_m^{}^{}(\omega ^{})],`$
$`p_m(\omega )`$ $`=`$ $`i\left[{\displaystyle \frac{\mathrm{}\omega }{2}}\right]^{\frac{1}{2}}[A_m^{}(\omega )`$ (14b)
$`{\displaystyle \underset{m^{}}{}}{\displaystyle }d\omega ^{}_{mm^{}}(\omega ,\omega ^{})A_m^{}(\omega ^{})].`$
The latter obey the commutation relations

$`[A_m(\omega ),A_m^{}(\omega ^{})]`$ $`=`$ $`0,`$
$`[A_m(\omega ),A_m^{}^{}(\omega ^{})]`$ $`=`$ $`\delta _{mm^{}}\delta (\omega \omega ^{}).`$ (15)
In terms of the creation and annihilation operators the Hamiltonian takes the familiar form

$`H={\displaystyle \frac{1}{2}}{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \mathrm{}\omega \left[A_m^{}(\omega )A_m(\omega )+A_m(\omega )A_m^{}(\omega )\right]},`$ (16)
describing a set of independent harmonic oscillators. Finally, substituting the representations (14) into Eqs. (4), one obtains the field expansions

$`𝐀`$ $`=`$ $`c{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \left[\frac{\mathrm{}}{2\omega }\right]^{\frac{1}{2}}\left[A_m(\omega ,t)𝐟_m(\omega ,𝐫)+\text{h.c.}\right]},`$ (17a)
$`𝚷`$ $`=`$ $`{\displaystyle \frac{i}{c}}{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \left[\frac{\mathrm{}\omega }{2}\right]^{\frac{1}{2}}\left[A_m(\omega ,t)𝐟_m(\omega ,𝐫)\text{h.c.}\right]},`$ (17b)
In empty space, $`ϵ(𝐫)1`$, they reduce to the standard mode expansion of the free electromagnetic field.

### II.2 Feshbach–projection

The modes–of–the–universe approach yields a consistent quantization scheme for the electromagnetic field, but does not provide explicit information about the field inside the resonator. The dynamics of this field is important in various contexts: Long–lived resonator modes are responsible for scattering resonances and, in the presence of an amplifying medium, may turn into lasing modes. To introduce resonator modes and to discuss their dynamics, we now separate the electromagnetic field into two contributions, accounting, respectively, for the field inside and outside the resonator. Formally, the separation of space in two regions is achieved using the projection operators Fes62

$`𝒬`$ $`=`$ $`{\displaystyle _{𝐫I}}𝑑𝐫|𝐫𝐫|,`$ (18a)
$`𝒫`$ $`=`$ $`{\displaystyle _{𝐫I}}𝑑𝐫|𝐫𝐫|,`$ (18b)
where $`|𝐫`$ denotes a standard position eigenket and $`I`$ is the region of space occupied by the dielectric material. This choice for $`I`$ is convenient but not unique: the only requirement is that there is free propagation of light in the external region far away from the resonator (this allows us to define asymptotic boundary conditions below). The operators $`𝒬`$ and $`𝒫`$ depend on the choice of $`I`$, but all physical observables turn out to be independent of this choice. One easily shows that $`𝒫`$ and $`𝒬`$ are projection operators,

$`𝒫=𝒫^{};`$ $`𝒫^2=𝒫,`$ (19a)
$`𝒬=𝒬^{};`$ $`𝒬^2=𝒬.`$ (19b)
Moreover, they are orthogonal, $`𝒬𝒫=𝒫𝒬=0`$, and complete, $`𝒬+𝒫=1`$. Therefore, an arbitrary Hilbert space function $`\mathit{\varphi }`$ and the associated function $`𝐟=\mathit{\varphi }/\sqrt{ϵ}`$ may be decomposed into the projections onto the resonator and channel space

$`\mathit{\varphi }(𝐫)`$ $`=`$ $`\chi _{}(𝐫)𝝁(𝐫)+\chi _+(𝐫)𝝂(𝐫),`$ (20a)
$`𝐟(𝐫)`$ $`=`$ $`\chi _{}(𝐫)𝐮(𝐫)+\chi _+(𝐫)𝐯(𝐫),`$ (20b)
where $`\chi _{}`$ are the characteristic functions of the resonator and the channel region, respectively,

$`\chi _{}(𝐫)={\displaystyle _{𝐫^{}I}}𝑑𝐫^{}\delta (𝐫𝐫^{}),\chi _+(𝐫)`$ $`=`$ $`1\chi _{}(𝐫).`$ (21)
Acting on $`\mathit{\varphi }`$ with the differential operator $`L`$ defined in Eq. (6), we obtain

$`L\mathit{\varphi }`$ $`=`$ $`\chi _{}(𝐫)L𝝁(𝐫){\displaystyle _I}d^2𝐫^{}K(𝐫,𝐫^{})𝝁(𝐫^{})`$ (22)
$`+\chi _+(𝐫)L𝝂(𝐫)+{\displaystyle _I}d^2𝐫^{}K(𝐫,𝐫^{})𝝂(𝐫^{}),`$
where $`K(𝐫,𝐫^{})`$ is a singular differential operator defined at the boundary

$`K(𝐫,𝐫^{})𝝁(𝐫^{})`$ $`=`$ $`\left[{\displaystyle \frac{\delta (𝐫𝐫^{})𝐧^{}}{\sqrt{ϵ(𝐫)}}}\right]\times \left[^{}\times {\displaystyle \frac{𝝁(𝐫^{})}{\sqrt{ϵ(𝐫^{})}}}\right]`$ (23)
$`+\left[{\displaystyle \frac{\delta (𝐫𝐫^{})}{\sqrt{ϵ(𝐫)}}}\right]\times \left[𝐧^{}\times {\displaystyle \frac{𝝁(𝐫^{})}{\sqrt{ϵ(𝐫^{}})}}\right].`$
Here $`^{}`$ denotes a derivative with respect to $`𝐫^{}`$ and $`𝐧^{}`$ is a unit vector normal to the boundary. The action of $`K(𝐫,𝐫^{})`$ on $`𝝂(𝐫^{})`$ is defined in a similar fashion, with $`𝝁(𝐫^{})`$ replaced by $`𝝂(𝐫^{})`$.

The first (third) term on the right hand side of Eq. (22) contribute only inside (outside) the dielectric. The second and fourth term are boundary terms. They involve $`𝝁`$, $`𝝂`$ and their derivatives at the boundary; these functions must be evaluated in the limit where the boundary is approached from the cavity and the channel region, respectively. We note that the boundary terms generally gives rise to singular behavior. As a result, the action of $`L`$ usually goes beyond the Hilbert space. The range of $`L`$ within the Hilbert space is defined by the functions for which the singular terms vanish; this happens in particular for the eigenfunctions of $`L`$ in Hilbert space.

We now want to decompose the operator $`L`$ into a resonator, a channel, and a coupling contribution. However, it is not obvious how the decomposition can be carried out for the singular boundary terms. We therefore replace the boundary integrals by integrations along surfaces arbitrarily close to the boundary but located inside respectively outside the resonator region. Equation (22) then becomes

$$L\mathit{\varphi }=L_{𝒬𝒬}𝝁+L_{𝒬𝒫}𝝂+L_{𝒫𝒬}𝝁+L_{𝒫𝒫}𝝂,$$
(24)
where $`L_{𝒬𝒬}`$ and $`L_{𝒫𝒫}`$ are the projections of $`L`$ onto the resonator and channel space,

$`L_{𝒬𝒬}𝝁=\chi _{}(𝐫)L𝝁(𝐫)`$
$`{\displaystyle _I}d^2𝐫_{}^{}\left[{\displaystyle \frac{\delta (𝐫𝐫_{}^{})𝐧^{}}{\sqrt{ϵ(𝐫)}}}\right]\times \left[^{}\times {\displaystyle \frac{𝝁(𝐫_{}^{})}{\sqrt{ϵ(𝐫_{}^{})}}}\right],`$ (25a)
$`L_{𝒫𝒫}𝝂=\chi _+(𝐫)L𝝂(𝐫)`$ (25b)
$`+{\displaystyle _I}d^2𝐫_+^{}\left[{\displaystyle \frac{\delta (𝐫𝐫_+^{})}{\sqrt{ϵ(𝐫)}}}\right]\times \left[𝐧^{}\times {\displaystyle \frac{𝝂(𝐫_+^{})}{\sqrt{ϵ(𝐫_+^{})}}}\right],`$
and $`L_{𝒬𝒫}`$ and $`L_{𝒫𝒬}`$ the coupling terms

$`L_{𝒬𝒫}𝝂`$ $`=`$ $`+{\displaystyle _I}d^2𝐫_{}^{}[{\displaystyle \frac{\delta (𝐫𝐫_{}^{})𝐧^{}}{\sqrt{ϵ(𝐫)}}}]\times \left[^{}\times {\displaystyle \frac{𝝂(𝐫_+^{})}{\sqrt{ϵ(𝐫_+^{})}}}\right],`$ (26a)
$`L_{𝒫𝒬}𝝁`$ $`=`$ $`{\displaystyle _I}d^2𝐫_+^{}[{\displaystyle \frac{\delta (𝐫𝐫_+^{})}{\sqrt{ϵ(𝐫)}}}]\times \left[𝐧^{}\times {\displaystyle \frac{𝝁(𝐫_{}^{})}{\sqrt{ϵ(𝐫_{}^{})}}}\right].`$ (26b)
The shorthands $`𝐫_{}^{}`$ indicate that the integrals are to be evaluated in the limit where $`𝐫_{}^{}`$ and $`𝐫_+^{}`$ approach the boundary from inside respectively outside the resonator. One can easily show that the decomposition (24) preserves the Hermiticity of $`L`$. Moreover, $`L_{𝒬𝒬}`$ and $`L_{𝒫𝒫}`$ define Hermitean operators on the Hilbert space functions of the resonator and the channel region, respectively.

Substitution of Eq. (24) into the eigenmode equation (6) yields two coupled equations for the projections of the eigenfunctions onto the resonator and channel space:

$$\left(\begin{array}{cc}L_{𝒬𝒬}& L_{𝒬𝒫}\\ L_{𝒫𝒬}& L_{𝒫𝒫}\end{array}\right)\left(\begin{array}{c}𝝁(\omega )\\ 𝝂(\omega )\end{array}\right)=\frac{\omega ^2}{c^2}\left(\begin{array}{c}𝝁(\omega )\\ 𝝂(\omega )\end{array}\right).$$
(27)
The condition that the singular terms on the left hand side vanish yields the two matching conditions

$`𝐧\times \left[𝐮(\omega )𝐯(\omega )\right]`$ $`=`$ $`\mathrm{𝟎},`$ (28a)
$`𝐧\times \left[\times 𝐮(\omega )\times 𝐯(\omega )\right]`$ $`=`$ $`\mathrm{𝟎},`$ (28b)
for all points along the boundary of the resonator region. The gauge condition $`[ϵ𝐟(\omega )]=0`$ and the requirement $`(\times 𝐟(\omega ))=0`$ along the boundary give the further matching conditions

$`𝐧\left[ϵ𝐮(\omega )ϵ𝐯(\omega )\right]`$ $`=`$ $`0,`$ (29a)
$`𝐧\left[\times 𝐮(\omega )\times 𝐯(\omega )\right]`$ $`=`$ $`0.`$ (29b)
The four matching conditions (28), (29) together with Eq. (1) realize the well–known Jac75  boundary conditions for the electromagnetic field at an interface in the absence of surface currents or surface charges.

### II.3 Eigenmodes of the resonator and channel region

The operators $`L_{𝒬𝒬}`$ and $`L_{𝒫𝒫}`$ are self–adjoint operators in the Hilbert space of the resonator and the channel functions, respectively. The eigenfunctions of $`L_{𝒬𝒬}`$ satisfy the equation

$$\times \left[\times 𝐮_\lambda (𝐫)\right]=\frac{ϵ(𝐫)\omega _\lambda ^2}{c^2}𝐮_\lambda (𝐫).$$
(30)
From the condition that the singular term in Eq. (25) vanishes one obtains the boundary condition

$$𝐧\times \left[\times 𝐮_\lambda \right]|_I=\mathrm{𝟎}.$$
(31)
Hence, the tangential component of $`\times 𝐮_\lambda `$ vanishes at the boundary. No boundary condition for the normal component of $`\times 𝐮_\lambda `$ is required as the three components of this vector are connected through the gauge condition. We note that the eigenmodes of the resonator form a discrete set.

In a similar fashion the eigenmodes of the channel region are found from Eq. (25b). They satisfy the equation

$$\times \left[\times 𝐯_m(\omega ,𝐫)\right]=\frac{ϵ(𝐫)\omega ^2}{c^2}𝐯_m(\omega ,𝐫),$$
(32)
and the condition that the tangential component must vanish at the boundary,

$$𝐧\times 𝐯_m(\omega )|_I=\mathrm{𝟎}.$$
(33)
The channel modes form a continuum, labeled by the frequency $`\omega `$ and the index $`m`$ that specifies the asymptotic conditions at infinity. We note that the resonator modes $`𝐮_\lambda `$ have support only within the resonator and vanish in the channel region; vice versa the channel functions $`𝐯_m(\omega )`$ vanish inside the resonator and take nonzero values only within the channel region. The resonator and channel modes form complete and orthonormal basis sets for the resonator and channel region, respectively. As a result, the projectors $`𝒫`$ and $`𝒬`$ can be represented in terms of these modes,

$`𝒬`$ $`=`$ $`{\displaystyle \underset{\lambda }{}}|𝝁_\lambda 𝝁_\lambda |,`$ (34a)
$`𝒫`$ $`=`$ $`{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega |𝝂_m(\omega )𝝂_m(\omega )|}.`$ (34b)
Together with the eigenmode equation (27) this reduces the eigenmode problem to the well–known problem Fan61 ; Dit00  of a discrete number of states coupled to a continuum.

We note that the boundary conditions (31), (33) on the resonator and channel modes are a consequence of our separation of the singular terms in Eq. (22). This separation is by no means unique: Different separations are possible and generally give rise to different boundary conditions. For example, the substitution of $`\delta (𝐫𝐫^{})`$ by $`\delta (𝐫𝐫_{+}^{}{}_{}{}^{})`$ and of $`\delta (𝐫𝐫^{})`$ by $`\delta (𝐫𝐫_{}^{}{}_{}{}^{})`$ in Eq. (23), leads to a new set of boundary conditions for which the conditions on the internal and external eigenmodes are just interchanged. Further choices are possible subject to the condition that the decomposition of $`L`$ is self–adjoint. The freedom in choosing the boundary conditions is characteristic for the projector technique Fes62 ; Som02 .

It is worth emphasizing that neither the modes $`𝝁_\lambda `$ of the closed resonator nor the channel modes $`𝝂_m(\omega )`$ represent eigenmodes of the total system. The latter modes satisfy the matching conditions derived earlier but, in general, neither of the boundary conditions (31), (33). Still the eigenmodes of the total system may be expanded in terms of the resonator and channel modes as these modes form complete basis sets in the respective regions. This fact is a consequence of the convergence in Hilbert space which does not imply pointwise convergence (at the boundary). Consequently, the matching conditions must not be imposed directly at the boundary but hold, as usual Jac75 , in a limiting sense infinitesimally close to the boundary. A further discussion of the physical meaning of the resonator and channel modes, as well as explicit applications of the projector technique to the potential scattering of matter waves can be found in Ref. Som02 .

The eigenmode equation (27) in full space may now be solved by standard methods Fan61 ; Dit00 ; Coh92 . For the projection onto the channel space one has the Lippmann–Schwinger type solution

$$𝒫|\mathit{\varphi }(\omega )=|𝝂(\omega )+\frac{1}{\frac{\omega ^2}{c^2}L_{𝒫𝒫}+iϵ}L_{𝒫𝒬}|\mathit{\varphi }(\omega ),$$
(35)
where the limit $`ϵ0^+`$ is implied. Substitution into the equation for the projection onto the resonator space yields

$$𝒬|\mathit{\varphi }(\omega )=\frac{1}{\frac{\omega ^2}{c^2}L_{\text{eff}}(\omega )}L_{𝒬𝒫}|𝝂(\omega ),$$
(36)
where $`L_{\text{eff}}`$ is the non–Hermitean operator

$$L_{\text{eff}}(\omega )L_{𝒬𝒬}+L_{𝒬𝒫}\frac{1}{\frac{\omega ^2}{c^2}L_{𝒫𝒫}+iϵ}L_{𝒫𝒬}.$$
(37)
To simplify notation, we introduce the Green function of the resonator in the presence of the coupling to the channels

$`G_{𝒬𝒬}\left({\displaystyle \frac{\omega ^2}{c^2}}\right)`$ $`=`$ $`{\displaystyle \frac{1}{\frac{\omega ^2}{c^2}L_{\text{eff}}(\omega )}}.`$ (38)
Combining Eqs. (35), (36) we arrive at an expression for the eigenstates $`|\mathit{\varphi }(\omega )`$,

$`|\mathit{\varphi }(\omega )=G_{𝒬𝒬}L_{𝒬𝒫}|𝝂(\omega )`$
$`+\left[1+{\displaystyle \frac{1}{\frac{\omega ^2}{c^2}L_{𝒫𝒫}+iϵ}}L_{𝒫𝒬}G_{𝒬𝒬}L_{𝒬𝒫}\right]|𝝂(\omega ).`$ (39)
Using the expansion (34), this yields an exact representation of the eigenstates in terms of the resonator and channel modes

$$|\mathit{\varphi }(\omega )=\underset{\lambda }{}\alpha _\lambda (\omega )|𝝁_\lambda +𝑑\omega ^{}\beta (\omega ,\omega ^{})|𝝂(\omega ^{}),$$
(40)
where the expansion coefficients are given by

$`\alpha _\lambda (\omega )=𝝁_\lambda |G_{𝒬𝒬}L_{𝒬𝒫}|𝝂(\omega ),`$ (41a)
$`\beta (\omega ,\omega ^{})=𝝂(\omega ^{})|[1+`$
$`{\displaystyle \frac{1}{\frac{\omega ^2}{c^2}+iϵL_{𝒫𝒫}}}L_{𝒫𝒬}G_{𝒬𝒬}L_{𝒬𝒫}]|𝝂(\omega ).`$ (41b)
The modes $`𝐟(\omega )`$ of the electromagnetic field are recovered from $`\mathit{\varphi }(\omega )`$ using Eqs. (5), (20):

$$𝐟(\omega ,𝐫)=\underset{\lambda }{}\alpha _\lambda (\omega )𝐮_\lambda (𝐫)+𝑑\omega ^{}\beta (\omega ,\omega ^{})𝐯(\omega ^{},𝐫).$$
(42)
### II.4 Field expansions and Hamiltonian

The decomposition of the electromagnetic field modes into a resonator and a channel contribution suggests a quantization scheme different from the modes–of–the–universe approach discussed in Sec. II.1. In this section we carry out the field quantization on the basis of the resonator and channel modes. Our starting point is the expansion of the vector potential and the canonical momentum in terms of these modes; combining Eqs. (4), (42) this expansion takes the form

$`𝐀(𝐫,t)`$ $`=`$ $`c{\displaystyle \underset{\lambda }{}}Q_\lambda 𝐮_\lambda (𝐫)+c{\displaystyle 𝑑\omega Q(\omega )𝐯(\omega ,𝐫)},`$ (43a)
$`𝚷(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{1}{c}}{\displaystyle \underset{\lambda }{}}𝐮_\lambda ^{}(𝐫)P_\lambda +{\displaystyle \frac{1}{c}}{\displaystyle 𝑑\omega 𝐯^{}(\omega ,𝐫)P(\omega )},`$ (43b)
where we have defined the position operators

$`Q_\lambda `$ $`=`$ $`{\displaystyle 𝑑\omega q(\omega )\alpha _\lambda (\omega )},`$ (44a)
$`Q(\omega )`$ $`=`$ $`{\displaystyle 𝑑\omega ^{}q(\omega ^{})\beta (\omega ,\omega ^{})},`$ (44b)
and the momentum operators

$`P_\lambda `$ $`=`$ $`{\displaystyle 𝑑\omega \alpha _\lambda ^{}(\omega )p(\omega )},`$ (45a)
$`P(\omega )`$ $`=`$ $`{\displaystyle 𝑑\omega ^{}\beta ^{}(\omega ,\omega ^{})p(\omega ^{})}.`$ (45b)
The $`Q_\lambda `$ and $`P_\lambda `$ are time–dependent operators that represent complex amplitudes associated with the resonator field. Likewise, the operators $`Q(\omega )`$ and $`P(\omega )`$ are amplitudes describing the channel field. The (equal–time) commutation relations of the various operators are discussed in Appendix A. The calculation shows that operators associated with different subsystems commute. Moreover, for each subsystem $`Q`$ and $`P`$ behave like the fundamental operators for position and momentum, respectively.

To discuss the dynamical evolution of the resonator and channel operators, we must express the field Hamiltonian in terms of these operators. We use Eqs. (44), (45) and the completeness relation $`𝒬+𝒫=1`$ to invert the relation between the operators for the total system and the operator for the two subsystems

$`q(\omega )`$ $`=`$ $`{\displaystyle \underset{\lambda }{}}\alpha _\lambda ^{}(\omega )Q_\lambda +{\displaystyle 𝑑\omega ^{}\beta ^{}(\omega ,\omega ^{})Q(\omega ^{})},`$ (46a)
$`p(\omega )`$ $`=`$ $`{\displaystyle \underset{\lambda }{}}P_\lambda \alpha _\lambda (\omega )+{\displaystyle 𝑑\omega ^{}P(\omega ^{})\beta (\omega ^{},\omega )}.`$ (46b)
Substitution into Eq. (11) yields the desired expression for the field Hamiltonian. Using relations between the expansion coefficients $`\alpha `$ and $`\beta `$ that follow from the completeness and orthogonality of the modes functions (see Appendix B) we can write the result in the form

$`H`$ $`=`$ $`{\displaystyle \underset{\lambda }{}}\left[P_\lambda ^{}P_\lambda +\omega _\lambda ^2Q_\lambda ^{}Q_\lambda \right]`$ (47)
$`+{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \left[P_m^{}(\omega )P_m(\omega )+\omega ^2Q_m^{}(\omega )Q_m(\omega )\right]}`$
$`+{\displaystyle \underset{\lambda }{}}{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \left[W_{\lambda m}(\omega )Q_\lambda ^{}Q_m(\omega )+\text{ h.c.}\right]},`$
with $`W_{\lambda m}(\omega )=(c^2/2)𝝁_\lambda |L|𝝂_m(\omega )`$. This shows that the operators of the subsystems do not simply oscillate, as they would if the subsystems were completely isolated from each other. The origin for this is the third term on the right hand side of Eq. (47) which couples the motion of the resonator and channel operators. The coupling reflects the fact that the boundary of the dielectric is not completely reflecting; thus radiation may leak through the boundary to the external radiation field.

The operators $`Q`$ and $`P`$ have a standard representation in terms of creation and annihilation operators,

$`Q_\lambda `$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}}{2\omega _\lambda }}\right]^{\frac{1}{2}}\left[a_\lambda +{\displaystyle \underset{\lambda ^{}}{}}𝒩_{\lambda \lambda ^{}}^{}a_\lambda ^{}^{}\right],`$ (48a)
$`P_\lambda `$ $`=`$ $`i\left[{\displaystyle \frac{\omega _\lambda }{2\mathrm{}}}\right]^{\frac{1}{2}}\left[a_\lambda ^{}{\displaystyle \underset{\lambda ^{}}{}}𝒩_{\lambda \lambda ^{}}a_\lambda ^{}\right],`$ (48b)
where the matrix element $`𝒩_{\lambda \lambda ^{}}`$ is the overlap integral

$`𝒩_{\lambda \lambda ^{}}={\displaystyle 𝑑𝐫ϵ(𝐫)𝐮_\lambda (𝐫)𝐮_\lambda ^{}(𝐫)}.`$ (49)
The operators $`a_\lambda `$ and $`a_\lambda ^{}`$ obey the canonical commutation relations

$`[a_\lambda ,a_\lambda ^{}^{}]=\delta _{\lambda \lambda ^{}},[a_\lambda ,a_\lambda ^{}]=0.`$ (50)
In a similar fashion one derives the representation of the channel operators $`Q_m(\omega )`$, $`P_m(\omega )`$ in terms of a continuous set of creation and annihilation operators $`b_m^{}(\omega )`$, $`b_m(\omega )`$. Substituting these representations into the Hamiltonian (47) and using the symmetry and unitarity of the overlap matrices, one arrives at the Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{\lambda }{}}\mathrm{}\omega _\lambda a_\lambda ^{}a_\lambda +{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \mathrm{}\omega b_m^{}(\omega )b_m(\omega )}`$ (51)
$`+\mathrm{}{\displaystyle \underset{\lambda }{}}{\displaystyle \underset{m}{}}{\displaystyle }d\omega [𝒲_{\lambda m}(\omega )a_\lambda ^{}b_m(\omega )`$
$`+𝒱_{\lambda m}(\omega )a_\lambda b_m(\omega )+\text{h.c.}],`$
where we have omitted an irrelevant zero point contribution. The coupling matrix elements are given by

$`𝒲_{\lambda m}(\omega )`$ $`=`$ $`{\displaystyle \frac{c^2}{2\mathrm{}\sqrt{\omega _\lambda \omega }}}𝝁_\lambda |L_{𝒬𝒫}|𝝂_m(\omega ),`$ (52a)
$`𝒱_{\lambda m}(\omega )`$ $`=`$ $`{\displaystyle \frac{c^2}{2\mathrm{}\sqrt{\omega _\lambda \omega }}}𝝁_\lambda ^{}|L_{𝒬𝒫}|𝝂_m(\omega ).`$ (52b)
The notation $`𝝁_\lambda ^{}|`$ means $`𝝁_\lambda ^{}|𝐫𝝁(𝐫)`$. Finally, substituting the representation (48) into Eq. (43) we find the expansion of the intracavity field

$`𝐀(𝐫,t)`$ $`=`$ $`c{\displaystyle \underset{\lambda }{}}\left[{\displaystyle \frac{\mathrm{}}{2\omega _\lambda }}\right]^{1/2}[a_\lambda 𝐮_\lambda (𝐫)+a_\lambda ^{}𝐮_\lambda ^{}(𝐫)],`$ (53a)
$`𝚷(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{i}{c}}{\displaystyle \underset{\lambda }{}}\left[{\displaystyle \frac{\mathrm{}\omega _\lambda }{2}}\right]^{1/2}[a_\lambda 𝐮_\lambda (𝐫)a_\lambda ^{}𝐮_\lambda ^{}(𝐫)].`$ (53b)
The Hamiltonian (51) and the field expansions (53) are the key results of the quantization procedure. The field expansions of the open resonator reduce precisely to the standard expressions known from closed resonators. However, the field dynamics is fundamentally different as shown below. We note that the resonator modes are coupled to the external radiation field via both resonant ($`a^{}b`$, $`b^{}a`$) and non–resonant ($`ab`$, $`a^{}b^{}`$) terms. The non-resonant terms become important in the case of overdamping (when the mode widths are comparable to the optical frequencies) Haa85 . In most cases of interest, the widths are much smaller than the relevant frequencies; then the rotating–wave approximation can be made, which amounts to keeping only the resonant terms in the Hamiltonian. In this approximation, the Hamiltonian (51) reduces to the well–known system–and–bath Hamiltonian Sen59 ; Gar00  of quantum optics. It has been argued, that this Hamiltonian is valid only for good cavities with spectrally well–separated modes. Our derivation shows that such pessimism is inappropriate: the system–and–bath Hamiltonian does describe the dynamics of overlapping modes, provided the broadening of these modes is much smaller than their frequency (so that non–resonant terms can be neglected).

## III Field dynamics

Measurements on optical cavities are typically done with detectors located in the external region outside the cavity. The detectors therefore measure the external field, and an input–output theory is required to relate the evolution of the external field to the dynamics of the system of interest. In this section we apply the input–output formalism of Gardiner and Collet Gar00 ; Wal94  to the dynamics generated by the system–and–bath Hamiltonian. We show that there is a linear relation $`b^{\text{out}}(\omega )=S(\omega )b^{\text{in}}(\omega )`$ between the cavity input and output field, involving the scattering matrix $`S(\omega )`$ of the cavity. Our derivation of $`S(\omega )`$ differs from most applications of the input–output formalism in two respects. First, we do not impose the Markov approximation and do not require that the coupling amplitudes between cavity and channel modes are frequency independent. Second, we do not restrict ourselves to essentially one–dimensional scattering, but consider the more general situation with multiple input and output channels and a non–trivial cavity dynamics. 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 irregularly shaped mirrors.

### III.1 Input–output relation

The starting point of the input–output theory Gar00 ; Wal94  are the equations of motion for the annihilation operators of the intracavity and channel field modes. These equations take a particular simple form when the damping rates are much smaller than the frequencies of interest. This permits us to use the rotating–wave approximation, in which the Hamiltonian (51) reduces to the system–and–bath Hamiltonian

$`H_{\text{SB}}`$ $`=`$ $`{\displaystyle \underset{\lambda }{}}\mathrm{}\omega _\lambda a_\lambda ^{}a_\lambda +{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \mathrm{}\omega b_m^{}(\omega )b_m(\omega )}`$ (54)
$`+\mathrm{}{\displaystyle \underset{\lambda }{}}{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \left[𝒲_{\lambda m}(\omega )a_\lambda ^{}b_m(\omega )+\text{ h.c. }\right]}.`$
Here we have extended the range of the frequency integrals from $`\mathrm{}`$ to $`\mathrm{}`$ consistent with the rotating–wave approximation Gar00 . The Heisenberg equations of motion for the internal operators $`a_\lambda `$ and the channel operators $`b_m(\omega )`$ are given by

$`\dot{a}_\lambda `$ $`=`$ $`i\omega _\lambda a_\lambda i{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega 𝒲_{\lambda m}(\omega )b_m(\omega )},`$ (55)
$`\dot{b}_m(\omega )`$ $`=`$ $`i\omega b_m(\omega )i{\displaystyle \underset{\lambda }{}}𝒲_{\lambda m}^{}(\omega )a_\lambda .`$ (56)
Integration of Eq. (56) starting from some initial time $`t_0<t`$ yields

$`b_m(\omega ,t)`$ $`=`$ $`e^{i\omega (tt_0)}b_m(\omega ,t_0)`$ (57)
$`i{\displaystyle \underset{\lambda }{}}𝒲_{\lambda m}^{}(\omega ){\displaystyle \underset{t_0}{\overset{t}{}}}𝑑t^{}e^{i\omega (tt^{})}a_\lambda (t^{}),`$
where $`b_m(\omega ,t_0)`$ denotes the channel operator $`b_m(\omega )`$ at time $`t_0`$. In an analogous fashion, one can express $`b_m(\omega ,t)`$ in terms of the channel operators at the final time $`t_1>t`$,

$`b_m(\omega ,t)`$ $`=`$ $`e^{i\omega (tt_1)}b_m(\omega ,t_1)`$ (58)
$`+i{\displaystyle \underset{\lambda }{}}𝒲_{\lambda m}^{}(\omega ){\displaystyle \underset{t}{\overset{t_1}{}}}𝑑t^{}e^{i\omega (tt^{})}a_\lambda (t^{}).`$
Eventually, we are interested in the limits $`t_0\mathrm{}`$ and $`t_1\mathrm{}`$. Substracting Eq. (57) from Eq. (58) and integrating the result over frequency, we obtain the input–output relation in the time domain,

$`b_m^{\text{out}}(t)b_m^{\text{in}}(t)=`$
$`{\displaystyle \frac{i}{2\pi }}{\displaystyle \underset{\lambda }{}}{\displaystyle 𝑑\omega 𝒲_{\lambda m}^{}(\omega )\underset{t_0}{\overset{t_1}{}}𝑑t^{}e^{i\omega (tt^{})}a_\lambda (t^{})}.`$ (59)
Equation (III.1) relates the cavity operators in the time interval $`t_0<t<t_1`$ to the input and output field operators

$`b_m^{\text{out}}(t)`$ $``$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑\omega e^{i\omega (tt_1)}b_m(\omega ,t_1)},`$ (60a)
$`b_m^{\text{in}}(t)`$ $``$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑\omega e^{i\omega (tt_0)}b_m(\omega ,t_0)}.`$ (60b)
Fourier transformation of Eq. (III.1) yields the input–output relation in the frequency domain. In the asymptotic limit $`t_0\mathrm{}`$, $`t_1\mathrm{}`$ the result takes the form

$`b^{\text{out}}(\omega )b^{\text{in}}(\omega )=i𝒲^{}(\omega )a(\omega ),`$ (61)
where we combined the input and output operators at frequency $`\omega `$ to $`M`$–component vectors. The coupling amplitudes $`𝒲_{\lambda m}`$ form an $`N\times M`$ coupling matrix $`𝒲`$, and the cavity mode annihilation 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^{\text{in}}(\omega )=e^{i\omega t_0}b(\omega ,t_0)`$ and $`b^{\text{out}}(\omega )=e^{i\omega t_1}b(\omega ,t_1)`$.

### III.2 S–matrix

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. Substitution of Eq. (57) into the equations of motion (55) for the cavity modes yields a set of linear differential equations which is readily solved by Fourier transformation. The result is

$$a(\omega )=2\pi D^1(\omega )𝒲(\omega )b^{\text{in}}(\omega ),$$
(62)
where $`D`$ is the $`N\times N`$ matrix with the matrix elements

$`D_{\lambda \lambda ^{}}(\omega )=(\omega \omega _\lambda )\delta _{\lambda \lambda ^{}}`$
$`+{\displaystyle \underset{\mathrm{}}{\overset{t}{}}}𝑑t^{}{\displaystyle 𝑑\omega ^{}e^{i(tt^{})(\omega ^{}\omega iϵ)}[𝒲(\omega )𝒲^{}(\omega )]_{\lambda \lambda ^{}}}.`$ (63)
We have taken the limit $`t_0\mathrm{}`$ and introduced the positive infinitesimal $`ϵ`$ to regularize the integral. Note that the positive sign of $`ϵ`$ is a consequence of the fact that $`a(\omega )`$ is expressed in terms of the channel operators in the remote past $`t_0\mathrm{}`$. The integral can be evaluated using the identity

$$\underset{ϵ0^+}{lim}\frac{1}{(\omega ^{}\omega )iϵ}=𝒫\left(\frac{1}{\omega ^{}\omega }\right)+i\pi \delta (\omega ^{}\omega ).$$
(64)
The result is

$`D_{\lambda \lambda ^{}}(\omega )`$ $`=`$ $`(\omega \omega _\lambda )\delta _{\lambda \lambda ^{}}`$ (65)
$`+\mathrm{\Delta }_{\lambda \lambda ^{}}(\omega )+i\pi \left[𝒲(\omega )𝒲^{}(\omega )\right]_{\lambda \lambda ^{}},`$
where $`\mathrm{\Delta }_{\lambda \lambda ^{}}`$ is the principal value integral

$$\mathrm{\Delta }_{\lambda \lambda ^{}}(\omega )=\underset{m}{}𝒫𝑑\omega ^{}\frac{𝒲_{\lambda m}(\omega ^{})𝒲_{m\lambda }^{}(\omega ^{})}{\omega ^{}\omega }.$$
(66)
Combining Eqs. (61), (62) we can eliminate the internal operators from the input–output relation. This yields a linear relation between the incoming and the outgoing field,

$$b^{\text{out}}(\omega )=S(\omega )b^{\text{in}}(\omega ),$$
(67)
where $`S`$ is the $`M\times M`$ scattering matrix,

$$S(\omega )=𝟙\mathrm{𝟚}\pi 𝕚𝒲^{}(\omega )𝔻^\mathrm{𝟙}(\omega )𝒲(\omega ).$$
(68)
This representation of the scattering matrix is well–known in nuclear and condensed matter physics Mah69 ; Bee97 ; Dit00 . It 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^{}(\omega ^{})]=\delta _{nm}\delta (\omega \omega ^{})`$ for both $`b=b^{\text{in}}`$ and $`b=b^{\text{out}}`$, one can easily show that $`S`$ is unitary, $`SS^{}=S^{}S=𝟙`$. 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 regime, where all matrix elements of $`𝒲𝒲^{}`$ are much smaller than the mean frequency spacing of the resonator modes. The opposite regime of overlapping resonances is realized when the damping rates exceed the mean frequency spacing.

According to Eqs. (67), (68) the field dynamics is governed by the resonances of the open cavity. The resonances are the complex poles of the S–matrix. They are the solutions of the equation $`\text{det}D(\omega )=0`$; from Eq. (65) 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. II.3).

Equation (68) has found widespread application in the random matrix theory of scattering Bee97 ; Dit00 . 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. (68). 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 Bee98 ; Mis00  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. IV.

### III.3 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 Gar00 ; Bee98  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

$`H`$ $`=`$ $`H_{\text{SB}}+H_{\text{abs}}+H_{\text{amp}},`$ (69)
where $`H_{\text{SB}}`$ is the system–and–bath Hamiltonian (54) while $`H_{\text{abs}}`$ and $`H_{\text{amp}}`$ represent the absorbing and amplifying bath,

$`H_{\text{abs}}={\displaystyle \underset{l}{}}{\displaystyle 𝑑\omega \mathrm{}\omega c_l^{}(\omega )c_l(\omega )}`$
$`+\mathrm{}{\displaystyle \underset{\lambda }{}}{\displaystyle \underset{l}{}}{\displaystyle 𝑑\omega \left[\kappa _{\lambda l}(\omega )a_\lambda ^{}c_l(\omega )+\text{h.c.}\right]},`$ (70a)
$`H_{\text{amp}}={\displaystyle \underset{k}{}}{\displaystyle 𝑑\omega \mathrm{}\omega d_k^{}(\omega )d_k(\omega )}`$
$`+\mathrm{}{\displaystyle \underset{\lambda }{}}{\displaystyle \underset{k}{}}{\displaystyle 𝑑\omega \left[\gamma _{\lambda k}(\omega )a_\lambda ^{}d_k(\omega )+\text{ h.c. }\right]}.`$ (70b)
The operators $`c_l`$, $`c_l^{}`$ obey the canonical commutation relations

$$[c_l(\omega ),c_l^{}^{}(\omega ^{})]=\delta _{ll^{}}\delta (\omega \omega ^{}),$$
(71)
and account for thermal emission within the absorbing medium. The operators $`d_k`$ and $`d_k^{}^{}`$ represent the amplifying medium and have the commutation relations Gar00

$$[d_k(\omega ),d_k^{}^{}(\omega ^{})]=\delta _{kk^{}}\delta (\omega \omega ^{}).$$
(72)
As the Hamiltonian (69) 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. III.2. The result is

$$b^{\text{out}}(\omega )=S(\omega )b^{\text{in}}(\omega )+U(\omega )c^{\text{in}}(\omega )+V(\omega )d^{\text{in}}(\omega ),$$
(73)
where $`c^{\text{in}}`$ and $`d^{\text{in}}`$ represent the input noise of the absorbing and amplifying bath. Both are integrals over bath operators at the initial time $`t_0`$,

$`c_l^{\text{in}}(t)`$ $``$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑\omega e^{i\omega (tt_0)}c_l(\omega ,t_0)},`$ (74a)
$`d_k^{\text{in}}(t)`$ $``$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle 𝑑\omega e^{i\omega (tt_0)}d_k(\omega ,t_0)}.`$ (74b)
The matrices $`U`$ and $`V`$ are given by

$`U(\omega )`$ $`=`$ $`2\pi i𝒲^{}(\omega )D^1(\omega )𝒦(\omega ),`$ (75a)
$`V(\omega )`$ $`=`$ $`2\pi i𝒲^{}(\omega )D^1(\omega )\mathrm{\Gamma }(\omega ),`$ (75b)
where the $`N\times L`$ matrix $`𝒦`$ and the $`N\times K`$ matrix $`\mathrm{\Gamma }`$ comprise the coupling amplitudes $`\kappa _{\lambda l}`$ and $`\gamma _{\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

$`D_{\lambda \lambda ^{}}(\omega )`$ $`=`$ $`(\omega \omega _\lambda )\delta _{\lambda \lambda ^{}}+\mathrm{\Delta }_{\lambda \lambda ^{}}(\omega )+i\pi \mathrm{\Sigma }_{\lambda \lambda ^{}}(\omega ),`$ (76a)
$`\mathrm{\Delta }_{\lambda \lambda ^{}}(\omega )`$ $`=`$ $`𝒫{\displaystyle 𝑑\omega ^{}\frac{\mathrm{\Sigma }_{\lambda \lambda ^{}}(\omega ^{})}{\omega ^{}\omega }},`$ (76b)
where $`\mathrm{\Sigma }`$ is the matrix

$`\mathrm{\Sigma }(\omega )=𝒲(\omega )𝒲^{}(\omega )+𝒦(\omega )𝒦^{}(\omega )\mathrm{\Gamma }(\omega )\mathrm{\Gamma }^{}(\omega ).`$ (77)
Using Eq. (73) and the commutation relations for the output and input noise operators, one obtains the relation

$$UU^{}VV^{}=𝟙𝕊𝕊^{},$$
(78)
that was first derived by Beenakker Bee98  using a scattering approach to field quantization. We note that the matrix $`𝟙𝕊𝕊^{}`$ is positive definite in an absorbing medium ($`V=0`$) and negative definite in an amplifying medium ($`U=0`$). The relations (73) and (78) 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 Bee98  to compute moments or even the full distribution of the output field intensity from linear random media.

### III.4 Langevin equations for the internal modes

It is frequently impractical or impossible to eliminate the cavity modes from the equations of motion. This happens, for example, when the cavity field is coupled to strongly pumped atoms. Then the dynamics of the total system comprising the field and the atoms becomes nonlinear, and the S–matrix approach of the previous section cannot be applied. A standard method for tackling interaction problems is to solve the equations of motion for the internal modes. These equations are quantum Langevin equations in which the coupling to the external radiation field gives rise to damping and noise.

We derive the Langevin equations in the Markov approximation. In particular, we will assume that the coupling amplitudes $`𝒲_{\lambda m}(\omega )`$ are independent of frequency Fyo97  over a sufficiently large frequency band centered around the frequency $`\omega _0`$ of interest (in a laser $`\omega _0`$ is the atomic transition frequency). Substituting Eq. (57) into the equation of motion (55), and performing the frequency integral, we obtain the Langevin equations

$`\dot{a}_\lambda (t)=i\omega _\lambda a_\lambda (t)\pi {\displaystyle \underset{\lambda ^{}}{}}[𝒲𝒲^{}]_{\lambda \lambda ^{}}a_\lambda ^{}(t)+F_\lambda (t),`$ (79)
where $`F_\lambda (t)`$ is the noise operator

$`F_\lambda (t)=i{\displaystyle 𝑑\omega e^{i\omega (tt_0)}\underset{m}{}𝒲_{\lambda m}b_m(\omega ,t_0)}.`$ (80)
The Eqs. (79) generalize the Langevin equation for a single cavity mode Wal94  to the case of many modesBar99 . We note that the resulting equations differ from the independent–oscillator equations of standard laser theory Hak84  in two respects: First, the mode operators $`a_\lambda `$ are coupled by the damping matrix $`𝒲𝒲^{}`$; second, the noise operators $`F_\lambda `$ are correlated, $`F_\lambda ^{}F_\lambda ^{}\delta _{\lambda \lambda ^{}}`$, as different modes couple to the same external channels (the expectation value is defined with respect to the channel oscillators at time $`t_0`$). The mode coupling by both damping and noise can be understood as a consequence of the fluctuation-dissipation theorem.

The origin of 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 $`𝒲𝒲^{}`$ are much smaller than the resonator mode spacing $`\mathrm{\Delta }\omega `$. This regime is realized in dielectrics that strongly confine light due to a large mismatch in the refractive index. To leading order in $`𝒲𝒲^{}/\mathrm{\Delta }\omega `$ only diagonal elements contribute to the damping matrix, and Eq. (79) 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. Coupled equations of motion are found when the damping rate becomes of the order or larger than the mean frequency spacing of the internal modes.

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 Bibel . The eigenvalues $`\omega _\lambda `$ display level repulsion and universal statistical properties. From Eq. (79), 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. Therefore, the spectrum of such resonators is governed by a non–Hermitean random matrix. We thus encounter an interesting connection between the spectral properties of chaotic optical resonators and non–Hermitean random matrices Sok89 ; Fyo97 ; Cha97 .

## IV Spontaneous emission

In the preceding sections have been concerned with the quantum properties of the electromagnetic field in open resonators. We now address the interaction of the radiation field with atoms. As a simple but nontrivial problem we consider the spontaneous emission of a two–level atom inside a cavity. This problem has attracted considerable interest Coh92 ; and it was found that the cavity may drastically modify the rate of spontaneous emission from its value in free space. The reason for the effect is the modification due to the cavity of the local density of modes at the position of the atom. Most investigations of the spontaneous emission rate assumed cavities of regular shape, but recently Mis97 ; Lam99  also unstable and chaotic cavities were addressed. We show below that our system–and–bath Hamiltonian reproduces the standard result for the atomic decay rate within the Wigner–Weisskopf approximation. We express the result in terms of left and right eigenmodes of a non–Hermitean matrix and demonstrate that for chaotic resonators a statistical analysis of the decay rate is possible using random matrix theory.

We consider a single two-level atom with transition frequency $`\omega _0`$ located at the position $`𝐫_0`$ inside an open cavity. The cavity is empty and defined by external mirrors of arbitrary shape. The coupling between the atom and the field is described in the dipole approximation; the dipole strength of the atomic transition is given by $`𝐝=0|e𝐫|1`$, where $`e`$ is the elementary charge and $`|1`$, $`|0`$ the excited state and the ground state of the atom. It is convenient to introduce the lowering and raising operators $`\sigma =|01|`$ and $`\sigma ^{}=|10|`$, respectively.

In the rotating wave approximation the total Hamiltonian for the field and the atom has the form

$$H=H_{\text{SB}}+\mathrm{}\omega _0\sigma ^{}\sigma +\underset{\lambda }{}\left[g_\lambda a_\lambda \sigma ^{}+\text{h.c.}\right].$$
(81)
The first term on the right hand side is the system–and–bath Hamiltonian (54), the second term represents the free Hamiltonian of the atom and the last terms account for the coupling between the atom and the cavity modes with the coupling amplitudes

$$g_\lambda =i\left(\frac{\mathrm{}\omega _\lambda }{2}\right)^{\frac{1}{2}}𝐝𝐮_\lambda (𝐫_0).$$
(82)
We note that the rotating wave approximation enters the Hamiltonian (81) in a twofold way: First we neglected rapidly oscillating non–resonant terms in the atom–field interaction. Second, we omitted the non–resonant terms in the field Hamiltonian, assuming that the decay rate of the cavity modes is much smaller than the frequencies of interest. This second approximation is convenient but not essential for the following calculation. The modifications that arise when the second approximation is dropped are summarized in Appendix C.

To compute the spontaneous emission rate we follow the standard Wigner–Weisskopf procedure. We assume that initially the atom is in the excited state while there is no photon in the radiation field. Hence, the state of the total system at time $`t=0`$ is given by $`|1,\text{vac}`$, where vac represents the vacuum state of the electromagnetic field. Since the Hamiltonian (81) conserves the total number of atom and field excitations, the time-dependent solution of the Schrödinger equation can be written in the form

$`|\mathrm{\Phi }(t)`$ $`=`$ $`c(t)|1,\text{vac}+{\displaystyle \underset{\lambda }{}}c_\lambda (t)|0,1_\lambda `$ (83)
$`+{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega c_m(\omega ,t)|0,1_m(\omega )},`$
where $`c_\lambda (t)`$ and $`c_m(\omega ,t)`$ are, respectively, the probability amplitude to find a single photon in the cavity mode $`\lambda `$ and in channel $`m`$ with frequency $`\omega `$. The time evolution of the amplitudes $`c(t)`$, $`c_\lambda (t)`$ and $`c_m(\omega ,t)`$ follows from the Schrödinger equation; an exact solution can be obtained using Laplace transformation. However, within the framework of the Wigner–Weisskopf approximation $`c(t)`$ decays exponentially

$$c(t)=\mathrm{exp}\left[i(\omega _0+\delta \omega _0)t\frac{\gamma }{2}t\right]c(0),$$
(84)
where $`\delta \omega _0`$ is a frequency shift and $`\gamma `$ the decay rate of the intensity $`|c(t)|^2`$. We note that an exponential decay is only found if the local density of modes is smooth on the scale of the atomic decay rate. The decay rate is given by

$$\gamma =\underset{ϵ0}{lim}\text{Re}\left[\frac{2}{\mathrm{}^2}\underset{ij}{}d_id_j^{}C_{ij}(\omega _0+iϵ)\right],$$
(85)
where $`i`$, $`j`$ label the components of the dipole matrix element. Here $`C_{ij}(\omega _0+iϵ)`$ is the Fourier transform of the two time correlation function of the electric field

$$C_{ij}(tt^{})\mathrm{\Theta }(tt^{})E_i^+(𝐫_0,t)E_j^{}(𝐫_0,t^{})_{\text{vac}},$$
(86)
$`E^\pm `$ denote the positive and negative frequency part of the electric field, $`\mathrm{\Theta }(tt^{})`$ is the step function, and the average $`\mathrm{}_{\text{vac}}`$ is the quantum average over the initial state of the field.

The electric field is connected with the canonical momentum field through the relation $`𝐄(𝐫,t)=c𝚷(𝐫,t)`$. Inside the cavity both can be expanded in terms of the cavity modes using Eq. (53). Substitution into Eq. (86) reduces the field correlation function to a sum over the Green functions of the cavity modes,

$`C_{ij}(\tau )`$ $`=`$ $`i{\displaystyle \underset{\lambda \lambda ^{}}{}}{\displaystyle \frac{\mathrm{}\sqrt{\omega _\lambda \omega _\lambda ^{}}}{2}}u_{\lambda i}(𝐫_0)u_{\lambda ^{}j}^{}(𝐫_0)G_{\lambda \lambda ^{}}(\tau ),`$ (87a)
$`G_{\lambda \lambda ^{}}(\tau )`$ $``$ $`i\mathrm{\Theta }(\tau )a_\lambda (\tau )a_\lambda ^{}^{}(0)_{\text{vac}}.`$ (87b)
To compute the Green functions, we differentiate Eq. (87b) with respect to $`\tau `$ and use the equations of motion (55) of the cavity operators $`a_\lambda `$. This yields the equations of motion of the Green functions

$`\dot{G}_{\lambda \lambda ^{}}(\tau )=\delta (\tau )G_{\lambda \lambda ^{}}(0)i\omega _\lambda G_{\lambda \lambda ^{}}(\tau )`$
$`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑\omega ^{}{\displaystyle \underset{0}{\overset{\tau }{}}}𝑑t^{}e^{i\omega ^{}(\tau t^{})}\left[𝒲(\omega ^{})𝒲^{}(\omega ^{})G(t^{})\right]_{\lambda \lambda ^{}}.`$ (88)
There is no contribution from the noise term in Eq. (55) as $`b_m(\omega ,0)a_\lambda ^{}(0)_{\text{vac}}=0`$. The initial condition at $`\tau =0`$ is $`G_{\lambda \lambda ^{}}(0)=i\delta _{\lambda \lambda ^{}}`$. Equation (IV) is readily solved by Fourier transformation. The result is

$$G(\omega )=D^1(\omega ),$$
(89)
where the non–Hermitean matrix $`D`$ was defined in Eq. (65). Substituting the result into the Fourier transform of Eq. (87a), we obtain the field correlation function in the frequency domain,

$$C_{ij}(\omega )=\frac{i\mathrm{}\omega _0}{2}\underset{\lambda \lambda ^{}}{}u_{\lambda i}(𝐫_0)u_{\lambda ^{}j}^{}(𝐫_0)\left[D^1(\omega )\right]_{\lambda \lambda ^{}},$$
(90)
where we again made use of the rotating wave approximation to replace $`\sqrt{\omega _\lambda \omega _\lambda ^{}}\omega _0`$. The decay rate follows upon substitution of Eq. (90) into Eq. (85),

$$\gamma =\frac{\omega _0}{\mathrm{}}\text{Im}\left[\underset{ij}{}d_id_j^{}\underset{\lambda \lambda ^{}}{}u_{\lambda i}(𝐫_0)\left[D^1(\omega _0)\right]_{\lambda \lambda ^{}}u_{\lambda ^{}j}^{}(𝐫_0)\right].$$
(91)
The sum over modes may be simplified in the eigenbasis of the non–Hermitean matrix $`D`$. In this basis, the double sum over the mode functions $`𝐮_\lambda `$ reduces to a summation over the left and right eigenmodes of the wave equation of the open cavity,

$`\gamma `$ $`=`$ $`{\displaystyle \frac{\pi \omega _0d^2}{\mathrm{}}}\rho (𝐫_0,\omega _0),`$ (92a)
$`\rho (𝐫_0,\omega _0)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}\text{Im}\left[{\displaystyle \underset{k}{}}{\displaystyle \frac{L_k^{}(𝐫_0,\omega _0)R_k(𝐫_0,\omega _0)}{\omega _k\omega _0i\frac{\mathrm{\Gamma }_k}{2}}}\right].`$ (92b)
Here, $`\rho (𝐫_0,\omega _0)`$ is the local density of modes at the position of the atom, $`L_k`$, $`R_k`$ denote the component along $`𝐝`$ in the left and right mode $`k`$, and $`\omega _k`$ and $`\mathrm{\Gamma }_k`$ are the mode frequency and the mode broadening.

Equations (92) are the final result for the decay rate. They describe spontaneous emission not only in cavities with quasi–discrete modes but also in unstable resonators with strongly overlapping modes. In the latter case, the left eigenfunctions of the cavity may differ strongly from the corresponding right eigenfunctions. Our result agrees with the decay rate derived in Refs. Lam02 ; Bro01  using a field expansion in terms of non–orthogonal modes. Equations (92) have recently Mis97  been used to calculate the distribution $`P(\gamma )`$ of decay rates for a two–level atom inside a chaotic cavity. The local density of modes also determines the photodissociation rate of small molecules with chaotic internal dynamics Fyo98 . Our derivation of the decay rate was based on the rotating wave approximation for the system–and–bath Hamiltonian. This approximation is valid for the (typical) case in which the broadening of the resonator modes is much smaller than the atomic transition frequency. When the mode broadening is of the order of the transition frequency, one can still compute the decay rate provided the coupling between the field and the atom is sufficiently small. The calculation is done in Appendix C.

## V Conclusion

We have presented an approach to the field quantization in optical resonators. Our quantization scheme applies to fields of arbitrary polarization and holds in the presence of an arbitrary number of escape channels.

An attractive feature of our approach is that it is based on a field expansion in terms of a set of orthogonal resonator and channel modes. The creation and annihilation operators associated with these modes obey canonical commutation relations. This is in contrast to the non–standard commutation relations found in alternative procedures that directly quantize the non–orthogonal modes of a non–Hermitean eigenvalue problem. Non–Hermitean operators enter our approach only through the mode dynamics.

In the case of weak damping our approach reduces to the well–known field quantization of standard laser theory. Then the resonator modes can be approximated by eigenmodes of a closed resonator. Each mode is damped due to escape from the cavity. Deviations from this simple dynamics show up when the damping rates are comparable with the frequency spacing of the resonator modes. Then the Langevin dynamics of the internal modes is coupled by damping and noise inflicted by the external radiation field. Our field quantization provides a unified description of both the regime of weak damping and the regime of strong damping and spectrally overlapping modes.

In disordered dielectrics light scatters chaotically due to spatial fluctuations of the dielectric constant. Chaotic scattering can be included in our approach by assuming that the cavity dynamics is described by a random matrix. The matrix is non–Hermitean since the cavity is coupled to the external radiation field. Averages over an appropriate ensemble of random matrix theory yield statistical information about the physical observables. In the present paper we demonstrated the connection with random matrix theory for linear optical media. Future application of the quantization technique to random lasers may allow to extend the connection with random matrix theory to non–linear optical media.

###### Acknowledgements.

We are indebted to H.–J. Sommers for important discussions concerning the Feshbach projector technique. We also enjoyed helpful discussions with F. Haake and D. V. Savin. This work has been supported by the SFB 237 “Unordnung und große Fluktuationen” der Deutschen Forschungsgemeinschaft.

## Appendix A Commutation relations for cavity and channels operators

In this appendix we compute the (equal–time) commutation relations for the cavity and channels position and momentum operators. The reality condition on $`𝐀`$ and $`𝚷`$ implies the following relations between the cavity operators and their adjoints

$`Q_\lambda `$ $`=`$ $`{\displaystyle \underset{\lambda ^{}}{}}𝒩_{\lambda \lambda ^{}}^{}Q_\lambda ^{}^{}`$ (93)
$`P_\lambda ^{}`$ $`=`$ $`{\displaystyle \underset{\lambda ^{}}{}}𝒩_{\lambda \lambda ^{}}^{}P_\lambda ^{}.`$ (94)
Likewise, the channel operators are connected with their adjoints via the relations

$`Q_m(\omega )`$ $`=`$ $`{\displaystyle \underset{m^{}}{}}{\displaystyle 𝑑\omega ^{}𝒩_{mm^{}}^{}(\omega ,\omega ^{})Q_m^{}(\omega ^{})},`$ (95)
$`P_m^{}(\omega )`$ $`=`$ $`{\displaystyle \underset{m^{}}{}}{\displaystyle 𝑑\omega ^{}𝒩_{mm^{}}^{}(\omega ,\omega ^{})P_m^{}(\omega ^{})}.`$ (96)
The matrix elements

$`𝒩_{\lambda \lambda ^{}}`$ $`=`$ $`{\displaystyle 𝑑𝐫𝝁_\lambda (𝐫)𝝁_\lambda ^{}(𝐫)},`$
$`𝒩_{mm^{}}(\omega ,\omega ^{})`$ $`=`$ $`{\displaystyle 𝑑𝐫𝝂_m(\omega ,𝐫)𝝂_m^{}(\omega ^{},𝐫)}.`$ (97)
are the expansion coefficients of the mode functions $`𝝁_\lambda `$ ($`𝝂_m(\omega )`$) in terms of the complex conjugate functions $`𝝁_\lambda ^{}`$ ($`𝝂_m^{}(\omega )`$). One can easily show that the matrices $`𝒩`$ are unitary, symmetric and that they only couple degenerate modes.

The equal–time commutation relations follow easily from the commutation relations (II.1) for the operators of the total system and the completeness of the modes $`|\mathit{\varphi }_m(\omega )`$. As an example we show that $`[Q_\lambda ,P_m(\omega )]=0`$. Using the definitions (44), (45) of $`Q_\lambda `$ and $`P_m(\omega )`$, we obtain

$`[Q_\lambda ,P_m(\omega )]`$ $`=`$ $`{\displaystyle \underset{m^{}m^{\prime \prime }}{}}{\displaystyle 𝑑\omega ^{}𝑑\omega ^{\prime \prime }\alpha _{m^{}\lambda }(\omega ^{})\beta _{m^{\prime \prime }m}^{}(\omega ^{\prime \prime },\omega )}`$ (98)
$`\times [q_m^{}(\omega ^{}),p_{m^{\prime \prime }}(\omega ^{\prime \prime })]`$
$`=`$ $`i\mathrm{}{\displaystyle \underset{m^{}}{}}{\displaystyle 𝑑\omega ^{}\alpha _{m^{}\lambda }(\omega ^{})\beta _{m^{}m}^{}(\omega ^{},\omega )},`$
where we used the commutation relation (II.1). According to Eq. (40), the coefficients $`\alpha `$ and $`\beta `$ can be written as

$`\alpha _{m^{}\lambda }(\omega ^{})`$ $`=`$ $`𝝁_\lambda |\mathit{\varphi }_m^{}(\omega ^{}),`$
$`\beta _{m^{}m}(\omega ^{},\omega )`$ $`=`$ $`𝝂_m(\omega )|\mathit{\varphi }_m^{}(\omega ^{}).`$ (99)
Substitution into the right hand side of Eq. (98) yields

$`{\displaystyle \underset{m^{}}{}}{\displaystyle 𝑑\omega ^{}\alpha _{m^{}\lambda }(\omega ^{})\beta _{m^{}m}^{}(\omega ^{},\omega )}`$
$`={\displaystyle \underset{m^{}}{}}{\displaystyle 𝑑\omega ^{}𝝁_\lambda |\mathit{\varphi }_m^{}(\omega ^{})\mathit{\varphi }_m^{}(\omega ^{})|𝝂_m(\omega )}`$
$`=0.`$ (100)
The calculation of all remaining commutators reduces to Eq. (A) or to one of the expressions

$`{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \alpha _{m\lambda }(\omega )\alpha _{m\lambda ^{}}^{}(\omega )}`$ $`=`$ $`\delta _{\lambda \lambda ^{}},`$ (101)
$`{\displaystyle \underset{m^{\prime \prime }}{}}{\displaystyle 𝑑\omega ^{\prime \prime }\beta _{mm^{\prime \prime }}(\omega ,\omega ^{\prime \prime })\beta _{m^{}m^{\prime \prime }}^{}(\omega ^{},\omega ^{\prime \prime })}`$ $`=`$ $`\delta _{mm^{}}\delta (\omega \omega ^{}).`$
One finds that the cavity operators have the commutation relations

$`\begin{array}{ccccc}[Q_\lambda ,Q_\lambda ^{}]\hfill & =& [Q_\lambda ,Q_\lambda ^{}^{}]\hfill & =& 0,\hfill \\ [P_\lambda ,P_\lambda ^{}]\hfill & =& [P_\lambda ,P_\lambda ^{}^{}]\hfill & =& 0,\hfill \\ [Q_\lambda ,P_\lambda ^{}]\hfill & =& i\mathrm{}\delta _{\lambda \lambda ^{}},\hfill & & \\ [Q_\lambda ,P_\lambda ^{}^{}]\hfill & =& i\mathrm{}𝒩_{\lambda \lambda ^{}}^{}.\hfill & & \end{array}`$ (107)
The channel operators have the commutation relations

$`\begin{array}{ccccc}[Q_m(\omega ),Q_n(\omega ^{})]\hfill & =& [Q_m(\omega ),Q_n^{}(\omega ^{})]\hfill & =& 0,\hfill \\ [P_m(\omega ),P_n(\omega ^{})]\hfill & =& [P_m(\omega ),P_n^{}(\omega ^{})]\hfill & =& 0,\hfill \\ [Q_m(\omega ),P_n(\omega ^{})]\hfill & =& i\mathrm{}\delta _{mn}\delta (\omega \omega ^{}),\hfill & & \\ [Q_m(\omega ),P_n^{}(\omega ^{})]\hfill & =& i\mathrm{}𝒩_{nm}^{}(\omega ^{},\omega ),\hfill & & \end{array}`$ (112)
and the cavity operators commute with all channel operators. This shows that for each subsystem the operators $`Q`$ and $`P`$ behave like the basic operators of position an momentum, respectively.

## Appendix B The Hamiltonian

We show how the Hamiltonian (47) is derived from the Hamiltonian (11) that involves operators associated with the eigenmodes of the total system. We separately treat the two contributions to the Hamiltonian (11) involving integrals over momentum and position operators, respectively. We start with the contribution

$$T=\frac{1}{2}\underset{m}{}𝑑\omega p_m^{}(\omega )p_m(\omega ).$$
(113)
Substitution of the representation (46) for $`p_m(\omega )`$ reduces the right hand side to three integrals which can be done using Eqs. (A)–(A). The result has the form

$$T=\frac{1}{2}\underset{\lambda }{}P_\lambda ^{}P_\lambda +\frac{1}{2}\underset{m}{}𝑑\omega P_m^{}(\omega )P_m(\omega ).$$
(114)
The second contribution

$`V={\displaystyle \frac{1}{2}}{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \omega ^2q_m^{}(\omega )q_m(\omega )}`$ (115)
is more difficult to compute due to the presence of the term $`\omega ^2`$ in the integral over frequency. Substituting the representation (46) for $`q_m(\omega )`$ into Eq. (115), we obtain

$`V={\displaystyle \frac{1}{2}}{\displaystyle \underset{\lambda \lambda ^{}}{}}Q_\lambda ^{}V_{\lambda \lambda ^{}}^{(1)}Q_\lambda ^{}`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{m^{}m^{\prime \prime }}{}}{\displaystyle 𝑑\omega ^{}𝑑\omega ^{\prime \prime }Q_m^{}^{}(\omega ^{})V_{m^{}m^{\prime \prime }}^{(2)}(\omega ^{},\omega ^{\prime \prime })Q_{m^{\prime \prime }}(\omega ^{\prime \prime })}`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{m^{}}{}}{\displaystyle \underset{\lambda }{}}{\displaystyle 𝑑\omega ^{}\left(Q_m^{}^{}(\omega ^{})V_{m^{}\lambda }^{(3)}(\omega ^{})Q_\lambda +\text{h.c.}\right)},`$ (116)
where $`V^{(1)}`$, $`V^{(2)}`$, and $`V^{(3)}`$ are integrals over the coefficients $`\alpha `$ and $`\beta `$,

$`V_{\lambda \lambda ^{}}^{(1)}`$ $`=`$ $`{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \omega ^2\alpha _{m\lambda }(\omega )\alpha _{m\lambda ^{}}^{}(\omega )},`$ (117)
$`V_{m^{}m^{\prime \prime }}^{(2)}`$ $`=`$ $`{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \omega ^2\beta _{m^{}m}(\omega ^{},\omega )\beta _{m^{\prime \prime }m}^{}(\omega ^{\prime \prime },\omega )},`$ (118)
$`V_{m^{}\lambda }^{(3)}(\omega ^{})`$ $`=`$ $`{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \omega ^2\beta _{m^{}m}(\omega ^{},\omega )\alpha _{m\lambda }^{}(\omega )}.`$ (119)
We compute these integrals using relations that follow from the eigenmode equation (27). Projecting this equation onto the cavity states $`𝝁_\lambda |`$ and the channel states $`𝝂_m(\omega )|`$, we obtain the set of coupled equations

$`\left[{\displaystyle \frac{\omega ^2}{c^2}}{\displaystyle \frac{\omega _\lambda ^2}{c^2}}\right]\alpha _{m\lambda }(\omega )={\displaystyle \frac{2}{c^2}}{\displaystyle \underset{m^{}}{}}{\displaystyle 𝑑\omega ^{}W_{\lambda n}(\omega ^{})\beta _{nm}(\omega ^{},\omega )},`$
(120)
$`\left[{\displaystyle \frac{\omega ^2}{c^2}}{\displaystyle \frac{\omega ^2}{c^2}}\right]\beta _{nm}(\omega ^{},\omega )={\displaystyle \frac{2}{c^2}}{\displaystyle \underset{\lambda }{}}W_{\lambda m}^{}(\omega ^{})\alpha _{\lambda m}(\omega ).`$
We used the definition $`W_{\lambda m}(\omega )=(c^2/2)𝝁_\lambda |L|𝝂_m(\omega )`$ and the definitions (A) of $`\alpha `$ and $`\beta `$. We multiply both sides of Eq. (B) by $`\alpha _{\lambda ^{}m}^{}(\omega )`$ and then sum over $`m`$ and integrate over $`\omega `$. This yields

$`V_{\lambda \lambda ^{}}^{(1)}`$ $`=`$ $`\omega _\lambda ^2{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \alpha _{\lambda ^{}m}^{}(\omega )\alpha _{\lambda m}(\omega )}`$ (122)
$`+2{\displaystyle \underset{m^{}}{}}{\displaystyle 𝑑\omega ^{}W_{\lambda m^{}}(\omega ^{})}`$
$`\times \left[{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \alpha _{\lambda ^{}m}^{}(\omega )\beta _{m^{}m}(\omega ^{},\omega )}\right].`$
The term in the square brackets vanishes according to Eq. (A). The remaining term on the right hand side can be simplified using Eq. (A), and Eq. (122) reduces to

$`V_{\lambda \lambda ^{}}^{(1)}`$ $`=`$ $`\omega _\lambda ^2{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega 𝝁_\lambda |\mathit{\varphi }_m(\omega )\mathit{\varphi }_m(\omega )|𝝁_\lambda ^{}}`$ (123)
$`=`$ $`\omega _\lambda ^2\delta _{\lambda \lambda ^{}}.`$
The expressions $`V^{(2)}`$ and $`V^{(3)}`$ can be computed in a similar fashion. Combining results the second contribution to the Hamiltonian takes the form

$`V`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\lambda }{}}\omega _\lambda ^2Q_\lambda ^{}Q_\lambda +{\displaystyle \frac{1}{2}}{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \omega ^2Q_m^{}(\omega )Q_m(\omega )}`$ (124)
$`+{\displaystyle \underset{m}{}}{\displaystyle \underset{\lambda }{}}{\displaystyle 𝑑\omega \left[W_{\lambda m}(\omega )Q_\lambda ^{}Q_m(\omega )+\text{ h.c.}\right]}.`$
The sum of the contributions (114) and (124) yields the Hamiltonian (47).

## Appendix C Spontaneous decay

We show here how to derive the atomic decay rate $`\gamma `$ without using the rotating wave approximation for the field Hamiltonian. Our starting point is Eq. (86), which written in terms of the exact modes of the total system takes the form Gla91

$`C_{ij}(tt^{})`$ $`=`$ $`\mathrm{\Theta }(tt^{}){\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \frac{\mathrm{}\omega }{2}}`$ (125)
$`\times f_{mi}(\omega ,𝐫_0)f_{mj}^{}(\omega ,𝐫_0)e^{i\omega (tt^{})}.`$
The Fourier transform of this equation is readily evaluated,

$$C_{ij}(\omega _0+iϵ)=i\underset{m}{}𝑑\omega \frac{\mathrm{}\omega }{2}\frac{f_{mi}(\omega ,𝐫_0)f_{mj}^{}(\omega ,𝐫_0)}{(\omega _0\omega )+iϵ}.$$
(126)
Substitution into Eq. (85) yields the golden rule result

$`\gamma `$ $`=`$ $`{\displaystyle \frac{\pi }{\mathrm{}^2}}{\displaystyle \underset{ij}{}}{\displaystyle \underset{m}{}}{\displaystyle 𝑑\omega \mathrm{}\omega }`$ (127)
$`\times \left[d_id_j^{}f_{mi}(\omega ,𝐫_0)f_{mj}^{}(\omega ,𝐫_0)\delta (\omega \omega _0)\right].`$
The atom is located inside the cavity. Therefore, we can use the mode expansion (42) to replace the modes of the total system by the cavity modes,

$`\gamma ={\displaystyle \frac{\pi }{\mathrm{}^2}}{\displaystyle \underset{ij}{}}{\displaystyle 𝑑\omega \mathrm{}\omega d_id_j^{}}`$
$`\times \left[{\displaystyle \underset{\lambda \lambda ^{}}{}}{\displaystyle \underset{m}{}}u_{\lambda i}(𝐫_0)u_{\lambda ^{}j}^{}(𝐫_0)\alpha _{\lambda m}(\omega )\alpha _{\lambda ^{}m}^{}(\omega )\right]\delta (\omega \omega _0).`$
The quantity in square brackets is proportional to the local density of states $`\rho (𝐫_0,\omega )`$,

$`{\displaystyle \frac{i\pi c^2}{\omega }}{\displaystyle \underset{m}{}}{\displaystyle \underset{\lambda \lambda ^{}}{}}𝐮_\lambda (𝐫_0)\alpha _{\lambda m}(\omega )\alpha _{\lambda ^{}m}^{}(\omega )𝐮_\lambda ^{}(𝐫_0)`$
$`=2i\text{Im}𝐫_0|G_{𝒬𝒬}\left({\displaystyle \frac{\omega ^2}{c^2}}\right)|𝐫_0.`$ (129)
Therefore, the atomic decay rate has the same form as in Eq. (92a) but with a modified local density of modes

$`\rho (𝐫_0,\omega _0)`$ $`=`$ $`{\displaystyle \frac{2\omega _0}{\pi c^2}}\text{Im}\left[{\displaystyle \underset{k}{}}{\displaystyle \frac{l_k^{}(𝐫_0,\omega _0)r_k(𝐫_0,\omega _0)}{\sigma _k(\omega _0)\left(\frac{\omega _0}{c}\right)^2}}\right],`$ (130)
where $`l_k`$ and $`r_k`$ are the left and right eigenmodes and $`\sigma _k`$ the eigenvalues of the non-Hermitian operator $`L_{\text{eff}}(\omega _0)`$. Equation (92) is recovered in the rotating wave approximation. Then the eigenmodes of $`L_{\text{eff}}`$ are simply related to the eigenmodes of the non–Hermitean matrix $`D^1`$, $`l_kL_k`$ and $`r_kR_k`$, and the eigenvalues of $`L_{\text{eff}}(\omega _0)`$ can be approximated by

$$\sigma _k\left(\frac{\omega _k}{c}\right)^2\frac{i\omega _0}{c^2}\mathrm{\Gamma }_k,$$
(131)
where $`\omega _k`$ and $`\mathrm{\Gamma }_k`$ are the mode frequency and the mode broadening.


format: text

text:
# Stability of gold nanowires at large Au-Au separations

## Abstract

The unusual structural stability of gold nanowires at large separations of gold atoms is explained from first-principles quantum mechanical calculations. We show that undetected light atoms, in particular hydrogen, stabilize the experimentally observed structures, which would be unstable in pure gold wires. The enhanced cohesion is due to the partial charge transfer from gold to the light atoms. This finding should resolve a long-standing controversy between theoretical predictions and experimental observations.

Gold is known as the noblest metal of all due to its inert behavior as a bulk material and low reactivity of its perfect surfaces. This property made gold probably the most valuable and desirable metal in modern human history. At the same time, it still challenges scientists with some puzzling features. Among them there are striking structural and electrical properties of gold nanocontacts. Such nanocontacts are obtained by the mechanically controllable break-junction (MCB) technique or by the scanning tunneling microscopy (STM), when the tip is first driven into the gold sample, and then retracted forming a nanowire . On stretching the nanowire undergoes a number of structural rearrangements finally producing a monoatomic gold bridge suspended between the tip and sample. In most experiments such nanowires are found to be amazingly stable, reaching the interatomic distances of about 3.6 - 4.0 $`\AA `$  i.e. significantly larger than the nearest-neighbor distance in bulk gold (2.88 $`\AA `$). The conductance of the nanocontacts is quantized in steps of $`2e^2/h`$ approaching the theoretical one-dimensional limit of $`2e^2/h`$ when the suspended wire becomes monoatomic . A different approach to produce and study gold nanowires is based on the transmission electron microscopy (TEM). In this technique gold nanowires are fabricated by perforating holes in a gold thin film using a focused electron beam with current density $``$ 100 $`A/cm^2`$. The structural reorganization of the nanowire is further studied by the electron beam with current density either reduced to $``$ 30 $`A/cm^2`$  or kept the same as in the preparation process . The results of the former experiments  generally agree with the results of the STM and MCB based experiments, finding long gold interatomic distances of 3.4 - 3.6 $`\AA `$ . On the contrary, the latter experiment  does not confirm the existence of long Au-Au bonds, reporting the bulk-like spacing of $``$ 2.9 $`\AA `$.

The issue of the novel gold nanowires has already been addressed in a large body of theoretical and experimental works  which describe many important features of gold nanocontacts, including their formation, structural rearrangements and conducting properties. However, a convincing explanation of the puzzling structural stability of gold monoatomic wires at large Au - Au separations, seen in most experiments, is still lacking. In the present report we show from first-principles quantum mechanical calculations that the unusual structural stability of gold nanowires at large separations of gold atoms can be explained by the presence of undetected light atoms able to stabilize the experimentally observed structures, which would be unstable in pure gold wires. In particular, we find hydrogen to be a primary suspect.

To understand the origin of the phenomena we have performed first-principles calculations in the framework of the density functional theory (DFT) . Two sets of data have been calculated at T = 0 K, first, for the gold nanowires supported by gold tips consisted of 28 atoms, with monoatomic part consisted of 2 to 6 atoms, and second, for infinite monoatomic chains with the similar number of gold atoms treated as inequivalent in the periodically reproduced unit cells. In agreement with previous findings , the results obtained for both sets (with and without supporting tips) are very similar, indicating that the stability of a nanocontact is essentially determined by its monoatomic part. No magnetization is found for any studied configuration. The binding energy curve for an infinite gold monoatomic wire is given in Fig. 1(a) together with the equilibrium atomic configurations and the ranges of their stability. Zigzag is an energetically favorable atomic configuration. Upon stretching the wire becomes linear at the Au-Au distance of about 2.7 $`\AA `$ remaining stable only in a narrow interval up to about 2.9 $`\AA `$, when the Peierls distortion triggers the dimerization of the wire accompanied by a metal-insulator transition  (Fig. 1(a)). Decreasing the Au-Au distance below the range shown in Fig. 1(a) leads to a series of deepening local energy minima and corresponds to the chain folding into the bulk-like structure.

These results suggest the wire to break (and obviously stop being conductive) already at interatomic distances of 2.9 $`\AA `$, in agreement with previous calculations  and the experiment by Takai et. al . They, however, by no means explain the stability of gold monoatomic wires with the spacing of 3.4 - 4.0 $`\AA `$ observed in most of other experiments . Therefore the existing first-principles theory does not provide a complete description of the phenomenon. It is especially remarkable regarding the fact that very different gold systems (bulk, surface, dimers, atomic clusters), and in particular their geometry, are known to be accurately described in the framework of the DFT . It seems to be a reasonable question to ask if there could be anything missing in the previous theoretical considerations of gold nanowires. For instance, it is tempting to check whether small amounts of light atoms, which are hard to detect in standard experiments, may influence the stability of gold nanocontacts.

We start with an observation that very similar results for the Au-Au spacing (3.4 - 4.0 $`\AA `$) have been obtained in different experiments, with and without (ultra) high vacuum conditions, in situ sample preparation, and controlled low temperature . Especially in the case of the ultra-high vacuum the probability to encounter most of the light atoms, such as, for instance, carbon and oxygen, is expected to be low. In particular it has been shown by Ugarte  and Rodrigues and Ugarte  that in the TEM-based experiments intense electron irradiation is helpful against carbon contamination. Nevertheless, such experiments  still show large Au-Au distances in gold nanowires. There is, however, at least one element, hydrogen, which is present even at high vacuum conditions . It is also known that gold becomes unrecognizably active when deposited as small particles on a substrate . This enhancement of gold reactivity is caused by a high concentration of steps and surface defects in small particles, as well as by the substrate-induced strain . Note that such low-dimensional gold structures are well presented on the surfaces of the STM tips. Gold surfaces interact with atomic hydrogen, which is able to form a stable layer on the surface . Dissolving into the bulk, hydrogen tends to segregate at the surface and subsurface layers . Gold has also been reported to play a role of a hydrogen trap in gold alloys that may be considered as an indication of a strong interaction between gold and hydrogen atoms .

We have extensively studied this interaction from first-principles in the systems of different scales: from nanostructures to the fcc bulk gold . Our results show that hydrogen molecules do dissociate when embedded into the gold matrix. Moreover, being dissolved in small quantities in the fcc bulk gold, hydrogen binds strongly to Au atoms with characteristic Au-H distances of about 1.75 $`\AA `$. Moderate stress and deformations of the gold host lattice do not break these bonds making hydrogen follow Au atoms.

The most profound influence of the hydrogen presence has been observed in gold nanowires. The binding energy curve for a one-dimensional gold-hydrogen wire is presented in Fig. 1(b). There is a sequence of structural transformations similar to that of a pure gold wire (see Fig. 1(a)). A zigzag arrangement of atoms is most stable for Au-Au distances smaller than 2.7 $`\AA `$, followed upon stretching by the linear configuration of the gold atoms with zigzag-like positions of the hydrogen atoms. This configuration is stable up to about 3.4 $`\AA `$ of Au-Au distance. On further stretching the chain becomes linear with hydrogen placed in between the gold atoms reaching the Au-Au spacing of about 3.8 $`\AA `$ and Au - H spacing of 1.85 $`\AA `$. Further stretching results in dimerization eventually followed by the complete break of the wire. Thus monoatomic gold-hydrogen wire shows an enhanced stability at large Au-Au distances that well matches the experimentally observed spacing of 3.6 - 4.0 $`\AA `$ . The cohesion is twice as strong as the one of a pure gold chain (see Fig. 1). As experiments on gold nanowires are often performed at room temperature, we have also studied the stability of monoatomic Au-H wires at finite temperatures by first-principles molecular dynamics (MD) simulations . According to our MD simulations the structures of the nanowires obtained at T = 0 K remain intact at least up to T = 300 K.

The nature of bonding in gold nanocontacts containing hydrogen can be understood by analyzing the charge density as well as the difference between the actual charge density and the atomic densities of gold and hydrogen constituting the system (Fig. 2). Charge density presented in Fig. 2(a) reveals pronounced electron density bridges caused by hydrogen situated in between the gold atoms. This resembles the hydrogen-type bonding where the charge transfer from Au to H, which can be seen in the charge density difference plot (Fig. 2(b)), is a key component. Therefore, an enhanced effective binding between the gold atoms via Au-to-H charge transfer mechanism results in the striking stability of the gold-hydrogen nanowires.

Further, in experiments  the conductance of gold nanocontacts has been measured at a low d.c. bias (about 10 mV). We have also studied the stability and ability of the gold-hydrogen chains to conduct at a similar bias, which is simulated by shifting the chemical potential at one end of the wire with respect to the other in the way described in Ref. . The band structure of such a wire (Fig. 3) clearly shows the existence of one conduction channel, in agreement with theory and most experimental findings . Additional forces appearing on atoms due to the applied bias are of about 0.001 eV/$`\AA `$ that is by no means sufficient for breaking the bonds.

In general one cannot rule out the possibility of the presence of other light elements. First, small amount of light atoms can be already present in the sample used to fabricate a nanowire. Second, not all experiments on gold nanowires employ high vacuum conditions. We have considered carbon and oxygen as possible contaminants (see Ref. for calculational details) and found that both of them stabilize monoatomic wires up to the Au-Au distances of $``$ 4.5 $`\AA `$. The enhanced cohesion is due to the bonding mechanism similar to the one operating in Au-H wires. The equilibrium Au-Au distances in Au-O and Au-C wires are, however, somewhat larger ($``$ 3.6 $`\AA `$ in both systems). Thus we notice that due to the larger equilibrium Au-Au separation in such systems, they should be in a "compressed" state in the region of experimentally observed Au-Au distances, contrary to the Au - H system. Therefore, if in a particular experiment oxygen or carbon atoms are responsible for the increased stability of gold nanowires, one would at least expect a different dependence of the force on the tips displacement from what is actually seen in the STM-based experiments .

Further, we comment on the experimental results by Takai et al.  reporting the bulk-like spacing in gold monoatomic wires and, therefore, confirming the theoretical results for pure gold wires. We suggest that this experiment was most probably free of contamination due to a permanent use of intense electron irradiation (current density $``$ 100 $`A/cm^2`$) during the observations. As a result, gold nanowires were stable at a bulk-like atomic separation of 2.9 $`\AA `$ within a very short time interval. It is in a sharp contrast with other experiments , where an electron beam with current density of 100 $`A/cm^2`$ was only used to perforate neighboring holes in gold films while the observations were performed at essentially lower current density ($``$ 30 $`A/cm^2`$). As we already mentioned, these experiments resulted in long-living monoatomic wires with large interatomic distances.

To summarize, we have shown that the striking stability of the suspended gold nanowires may be attributed to the hydrogen atoms, which mediate gold-gold interactions, dramatically increasing binding. Hydrogen, being always present but rarely controlled in experiments, alters the properties of gold nanowires, the crucial components of nanometer-scale devices developed for microelectronics and biophysics applications . In addition we note that other light atoms, like carbon and oxygen, are also able to stabilize gold nanowires. Our calculations show that Au-O and Au-C monoatomic wires are stable up to the Au-Au distances of 4.5 $`\AA `$. We suggest that further experiments, incorporating hydrogen and other light atoms detection, may lead to important implications for the modern atomic-size technologies.

We thank M. Springborg, K. Hermannson, Y. Andersson, K. Forsgren, G. Källén, S. Kubatkin and D. Erts for discussions. Support from the Swedish Foundation for Strategic Research and Swedish Research Council is gratefully acknowledged.