format: text

text:
# Dynamical Casimir Effect in a Leaky Cavity at Finite Temperature

## I Introduction

Since the pioneering work of Casimir casimir  the phenomena of quantum field theory under the influence of external conditions have attracted the interest of many authors, see e.g. bordag . The original prediction by Casimir, i.e., the attractive force generated between two perfectly conducting objects placed in the vacuum, has been verified in different experimental setups with relatively high precision lamoreauxmohideenbressi . However, its dynamic counterpart with non-stationary boundary conditions inducing interesting effects like the creation of particles out of the vacuum has not yet been observed rigorously in a corresponding experiment. The observation of quantum radiation could provide a substantial test of the foundations of quantum field theory and thus be of special relevance. Generally we understand the term quantum radiation to denote the conversion of virtual quantum fluctuations into real particles due to external disturbances. For the special case of the external disturbances being moving mirrors this phenomenon is known as the Dynamical Casimir effect.

These striking effects have been investigated by several authors, for an overview see e.g. bordag ; jaekelreynaud  and references therein. We will focus on the effect of particle creation within a constructed – resonantly vibrating – leaky cavity. This case is of special importance for an experimental verification of the dynamical Casimir effect since the generation of particles is enhanced drastically by resonance effects. Employing different methods and approaches it has already been shown for ideal cavities (see e.g. finitetemp ) that under resonance conditions (i.e., when one of the boundaries performs harmonic oscillations at twice the frequency of one of the eigenmodes of the cavity) the phenomenon of parametric resonance (see e.g. jijungparksoh ) will occur. In the case of an ideal cavity (i.e., one with perfectly reflecting mirrors) this is known to lead to an exponential growth of the resonance mode particle occupation numbers, cf. finitetemp ; dalvit ; dodonov1 ; dodonov2 ; dodonov3 .

In view of this prediction an experimental observation of quantum radiation using the dynamical Casimir effect appears to be rather simple – provided the cavity is vibrating at resonance for a sufficiently long period of time. However, this point of view is too naive since neither ideal cavities do exist nor is it possible to match the external frequency to the fundamental eigenfrequency of the cavity with arbitrary precision. Consequently, it is essential to include effects of leaks as well as effects of detuning, see also golestanian .

Investigations concerning effects of losses have been performed for example in jrrad  in 1+1 space-time dimensions based on conformal mapping methods as developed in davies . However, these considerations are a priori restricted to 1+1 dimensions and can not be obviously generalized to higher dimensions. In 3+1 dimensions the character of the mechanism generating quantum radiation – e.g., the resonance conditions – differs drastically from the 1+1 dimensional situation.

More realistic (3+1 dimensional) cavities were considered in leakydodonov  where the effects of losses were taken into account by virtue of a master equation ansatz. However, this master equation had not been derived starting from first principles. It has already been noted in leakydodonov  that the employed ansatz is adequate for a stationary cavity – but not necessarily for a dynamic one. In addition, most papers did not include temperature effects – which may contribute significantly in an experiment. It has been shown in Ref. finitetemp  that for an ideal cavity the effect of particle production at finite temperature is enhanced by several orders of magnitude in comparison with the pure vacuum contribution.

In this article we will adopt the canonical approach which has proven to be general, successful, and is – in addition – also capable of including temperature effects. However, the aforementioned approach still lacks a generalization for systems with losses. We are aiming at providing a remedy in this field schaller .

This paper is organized as follows: In Sec. II we present a model system and derive the effective Hamiltonian for the resonance case. In Sec. III we will calculate the number of created particles in the cavity after one of the walls has performed resonant oscillations by means of response theory. In Sec. IV we will derive and solve the associated master equation and show consistency with the results obtained in Sec. III. In Sec. V a non-perturbative approach is presented and compared with the other results. We derive a treshold condition - valid for leaky cavities - for a possible detuning from the fundamental resonance in Sec. VI . We shall close with a summary, a discussion, a conclusion and an outlook.

Throughout this paper natural units given by $`\mathrm{}=c=k_B=1`$ will be used.

## II General Formalism

### II.1 The Leaky Cavity

We want to investigate the effects of a non-ideal cavity in view of the dynamical Casimir effect. For that purpose we have to construct a suitable model system. One simple way to do that is to insert a dispersive mirror into an ideal cavity while keeping all other walls perfectly reflecting. Thereby two leaky cavities are formed. Particles in the left imperfect cavity are now able to leave into the right larger box (the reservoir). For reasons of simplicity we consider a rectangular cavity as depicted in Fig. 1.

The setup in Fig. 1 is not a new idea. A similar – but static – system has already been treated in scully ; banachloche . However, here in addition the left wall is moving with a prescribed trajectory during the time interval $`[0,T]`$. For ideal cavities this is known to lead to a squeezing of the vacuum state which causes the creation of particles inside the cavity, see e.g. finitetemp .

Note that we are assuming a finite reservoir with a discrete spectrum instead of an infinite one leading to a continuum of modes. Since, in an experimental setup, the vibrating cavity will most likely be surrounded by walls, etc. (imposing additional boundary conditions), this assumption should be justified.

Assuming a surrounding perfectly reflecting wall is a first idealization of the real situation. However, in order to minimize the error obtained by this procedure the experiment could be designed in this way, see also Fig. 8 in section IX below.

The ideal mirrors can be simulated by infinitely high potential walls inducing Dirichlet boundary conditions. For the additional dispersive mirror we use the $`\delta `$-type model potential proposed in Ref. calogeracos ; salomone

$`V(x;t)=\{\begin{array}{cc}\gamma \delta (xb)\hfill & \text{if }a(t)<x<c\hfill \\ \mathrm{}\hfill & \text{otherwise}\hfill \end{array},`$ (3)
see also Fig. 2. The parameter $`\gamma `$ represents the transmittance of the internal mirror, whose reflection and transmission amplitudes are determined as calogeracos

$`={\displaystyle \frac{i\gamma }{\omega +i\gamma }},𝒯={\displaystyle \frac{\omega }{\omega +i\gamma }}.`$ (4)
Note that the general procedure presented in this article is independent of the particular form of the potential – the aforementioned one has just been chosen for convenience. For a more realistic scenario one could apply square-well or Gaussian potentials. In a realistic experiment where one would want to create photons instead of scalar particles a dispersive mirror could be realized using a thin dielectric slab with a very high dielectric constant. Such a mirror could then be approximated by a space-dependent permittivity $`\epsilon (x)=1+\alpha \delta (x)`$. This will be addressed in Sec. X.2.

### II.2 Hamiltonian

Throughout this article we will use the notation of tremblingcav  where the particle production in an ideal vibrating cavity was calculated – for a more general treatment see e.g. blackhole . We consider a massless and neutral scalar field coupled to an external potential:

$`={\displaystyle \frac{1}{2}}(_\mu \mathrm{\Phi })(^\mu \mathrm{\Phi })V\mathrm{\Phi }^2.`$ (5)
The perfect mirrors can be incorporated by imposing the corresponding boundary condition on $`\mathrm{\Phi }`$. By expanding the field $`\mathrm{\Phi }`$

$`\mathrm{\Phi }(𝒓,t)={\displaystyle \underset{\mu }{}}Q_\mu (t)f_\mu (𝒓;t)`$ (6)
into a complete and orthonormal set of functions $`f_\mu (𝒓;t)`$ satisfying

$`{\displaystyle d^3rf_\mu ^{}(𝒓,t)f_\nu (𝒓,t)}=\delta _{\mu \nu },`$ (7)
$`{\displaystyle \underset{\mu }{}}f_\mu ^{}(𝒓;t)f_\mu (𝒓^{\mathbf{}};t)=\delta ^3(𝒓𝒓^{\mathbf{}}),`$ (8)
$`\{2V\mathrm{\Delta }\}f_\mu (𝒓;t)=\mathrm{\Omega }_\mu ^2(t)f_\mu (𝒓;t)`$ (9)
one can reach a more convenient form suitable for doing calculations. Since $`\mathrm{\Phi }`$ is a real field, we can choose the set $`f_\mu `$ to be real. Note that the time dependence of eigenfunctions and eigenfrequency is solely induced by the moving boundary. Inserting this expansion into Eq. (5) transforms the Lagrangian into tremblingcav

$`L={\displaystyle d^3r}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu }{}}\dot{Q}_\mu ^2{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu }{}}\mathrm{\Omega }_\mu ^2(t)Q_\mu ^2`$ (10)
$`+{\displaystyle \underset{\mu \nu }{}}Q_\mu M_{\mu \nu }(t)\dot{Q}_\nu `$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu \nu \kappa }{}}Q_\mu M_{\mu \kappa }(t)M_{\nu \kappa }(t)Q_\nu ,`$
where $`M_{\mu \nu }(t)`$ is an antisymmetric matrix given by

$`M_{\mu \nu }={\displaystyle d^3r\frac{f_\mu }{t}(𝒓;t)f_\nu (𝒓;t)}.`$ (11)
This matrix describes the coupling strength between two different modes. We introduce the canonical conjugated momenta

$`P_\mu ={\displaystyle \frac{L}{\dot{Q}_\mu }}=\dot{Q}_\mu +{\displaystyle \underset{\nu }{}}Q_\nu M_{\nu \mu }(t).`$ (12)
Furthermore we apply the usual Legendre transform to a Hamiltonian representation and perform the quantization. This yields

$`\widehat{H}={\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu }{}}\widehat{P}_\mu ^2+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu }{}}\mathrm{\Omega }_\mu ^2\widehat{Q}_\mu ^2+{\displaystyle \underset{\mu \nu }{}}\widehat{P}_\mu M_{\mu \nu }\widehat{Q}_\nu .`$ (13)
The above Hamiltonian can be sub-classified into

$`\widehat{H}=\widehat{H}_0+\widehat{H}_I^S+\widehat{H}_I^V,`$ (14)
where the single Hamiltonians are given by

$`\widehat{H}_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu }{}}\widehat{P}_\mu ^2+{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu }{}}(\mathrm{\Omega }_\mu ^0)^2\widehat{Q}_\mu ^2,`$ (15)
$`\widehat{H}_I^S`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\mu }{}}\mathrm{\Delta }\mathrm{\Omega }_\mu ^2(t)\widehat{Q}_\mu ^2,`$ (16)
$`\widehat{H}_I^V`$ $`=`$ $`{\displaystyle \underset{\mu \nu }{}}\widehat{P}_\mu M_{\mu \nu }(t)\widehat{Q}_\nu .`$ (17)
The deviation $`\mathrm{\Delta }\mathrm{\Omega }_\mu ^2(t)=\mathrm{\Omega }_\mu ^2(t)(\mathrm{\Omega }_\mu ^0)^2`$ denotes the difference of the (squared) time-dependent eigenfrequencies $`\mathrm{\Omega }_\mu ^2(t)`$ from the unperturbed ones $`(\mathrm{\Omega }_\mu ^0)^2`$. The first term $`\widehat{H}_0`$ is the Hamiltonian of harmonic oscillators. The remaining terms will further on be called squeezing interaction Hamiltonian and velocity interaction Hamiltonian. We want to point out that in the case of a static system (where the eigenfunctions $`f_\mu `$ and eigenfrequencies $`\mathrm{\Omega }_\mu `$ are constant in time) the complete interaction Hamiltonian $`\widehat{H}_I=\widehat{H}_I^S+\widehat{H}_I^V`$ will vanish. The derivation of the eigenfunctions $`f_\mu (𝒓,t)`$ and $`M_{\mu \nu }(t)`$ will be treated in the following subsection.

### II.3 Eigenmodes

As has already been mentioned, we want to find a set of functions satisfying $`\{2V(𝒓;t)\mathrm{\Delta }\}f_\mu (𝒓;t)=\mathrm{\Omega }_\mu ^2(t)f_\mu (𝒓;t)`$. Any time dependence can only be induced by the moving boundaries. At first we will just consider the spatial dependence i.e., the stationary problem. The differential equation can be treated using the separation ansatz $`f_\mu (𝒓)=f_\mu (𝒓_{})f_\mu (𝒓_{})`$ where $`f_\mu (𝒓_{})`$ depends only on the coordinate parallel to the wall velocity and $`f_\mu (𝒓_{})`$ is dependent on the perpendicular coordinates. For the special case of our model system this means $`f_\mu (𝒓)=f_\mu ^x(x)f_\mu ^y(y)f_\mu ^z(z)`$ leading to the trivial $`y`$ and $`z`$ dependence of the eigenfunctions

$`f_\mu ^y(y)`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{\mathrm{\Delta }y}}}\mathrm{sin}\left[{\displaystyle \frac{n_y\pi }{\mathrm{\Delta }y}}y\right],\mathrm{\Omega }_\mu ^y={\displaystyle \frac{n_y\pi }{\mathrm{\Delta }y}},`$ (18)
$`f_\mu ^z(z)`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{\mathrm{\Delta }z}}}\mathrm{sin}\left[{\displaystyle \frac{n_z\pi }{\mathrm{\Delta }z}}z\right],\mathrm{\Omega }_\mu ^z={\displaystyle \frac{n_z\pi }{\mathrm{\Delta }z}}`$ (19)
with $`\mathrm{\Delta }y`$ and $`\mathrm{\Delta }z`$ denoting the dimensions of the cavity and the frequencies relating via

$`\mathrm{\Omega }_\mu ^2=(\mathrm{\Omega }_\mu ^x)^2+(\mathrm{\Omega }_\mu ^y)^2+(\mathrm{\Omega }_\mu ^z)^2.`$ (20)
The remaining differential equation reads

$`\left\{2\gamma \delta (xb)_x^2\right\}f_\mu ^x(x)=(\mathrm{\Omega }_\mu ^x)^2f_\mu ^x(x)`$ (21)
where the Dirichlet boundary conditions coming from the perfect mirrors on either side can be satisfied by the ansatz

$`f_\mu ^x(x)=\{\begin{array}{cc}L_\mu \mathrm{sin}[\mathrm{\Omega }_\mu ^x(xa)]\hfill & \text{if }a<x<b\hfill \\ R_\mu \mathrm{sin}[\mathrm{\Omega }_\mu ^x(cx)]\hfill & \text{if }b<x<c\hfill \\ 0\hfill & \text{elsewhere}\hfill \end{array}.`$ (25)
The eigenfunctions have to obey the continuity conditions calogeracos

$`f_\mu ^x(xb)f_\mu ^x(xb)`$ $`=`$ $`0,`$ (26)
$`{\displaystyle \frac{f_\mu ^x}{x}}(xb){\displaystyle \frac{f_\mu ^x}{x}}(xb)`$ $`=`$ $`2\gamma f_\mu ^x(b),`$ (27)
where the latter can be obtained via integration. These conditions can be combined to an eigenvalue equation for $`\mathrm{\Omega }_\mu ^x`$

$`{\displaystyle \frac{2\gamma }{\mathrm{\Omega }_\mu ^x}}=\mathrm{cot}\left[\mathrm{\Omega }_\mu ^x(ba)\right]+\mathrm{cot}\left[\mathrm{\Omega }_\mu ^x(cb)\right]={\displaystyle \frac{2}{\eta _\mu }}.`$ (28)
Though there is no obvious analytical solution of this equation, a numerical solution can always be obtained for given cavity parameters $`\{a,b,c,\gamma \}`$. However, via introducing the dimensionless perturbation parameter $`\eta _\mu =\mathrm{\Omega }_\mu ^x/\gamma `$ it is also possible to obtain an approximate analytical solution. Note that this parameter is small $`\eta _\mu 1`$ in the limit of the internal mirror being nearly perfectly reflecting. Since the trigonometric functions are very sensitive to small frequency variations one can solve the equation using a series expansion in $`\eta _\mu `$. It is obvious that if the right hand side goes to $`\mathrm{}`$ one of the addends or even both can become relevant. This depends on the ratio $`(ba)/(cb)`$ and its inverse which are both assumed to be non-integer numbers in the following non-perturbative calculations implying that only one of the addends is dominating. Accordingly, expanding around the poles of one addend one yields a polynomial that can be solved for $`\mathrm{\Omega }_\mu ^x`$ as a series expansion in $`\eta _\mu 1`$. Depending on the chosen addend one obtains two sets of approximate eigenfrequencies

$`\mathrm{\Omega }_{n_x,l}^x`$ $`=`$ $`{\displaystyle \frac{n_x\pi }{ba}}{\displaystyle \frac{1}{2(ba)}}\eta _{n_x,l}`$
$`+{\displaystyle \frac{1}{4(ba)}}\mathrm{cot}\left(n_x\pi {\displaystyle \frac{cb}{ba}}\right)\eta _{n_x,l}^2`$
$`+𝒪\left(\eta _{n_x,l}^3\right),`$
$`\mathrm{\Omega }_{n_x,r}^x`$ $`=`$ $`{\displaystyle \frac{n_x\pi }{cb}}{\displaystyle \frac{1}{2(cb)}}\eta _{n_x,r}`$ (29)
$`+{\displaystyle \frac{1}{4(cb)}}\mathrm{cot}\left(n_x\pi {\displaystyle \frac{ba}{cb}}\right)\eta _{n_x,r}^2`$
$`+𝒪\left(\eta _{n_x,r}^3\right),`$
which constitute a determining polynomial for $`\mathrm{\Omega }_\mu ^x`$. Note that the index $`\mu =(n_x,l/r)`$ is a multi-index, where $`l`$ and $`r`$ stand for left-dominated and right-dominated, respectively. However, it can be shown easily that the quality of the linear (in $`\eta `$) approximation suffices already for moderate values of $`\gamma 50`$. The insertion of (II.3) into the ansatz (25) leads to two classes of eigenfunctions: left-dominated and right-dominated, respectively. The differences between those are clearly visible in Fig. 3.

In order to avoid the confusion arising from a set of perturbation parameters $`\{\eta _\mu \}`$ we will introduce the fundamental one via

$`\eta =\eta _{1l}={\displaystyle \frac{\mathrm{\Omega }_{1l}^x}{\gamma }},`$ (30)
to which all others are evidently related via $`\eta _\mu =\mathrm{\Omega }_\mu ^x/\mathrm{\Omega }_{1l}^x\eta `$. Note that this distinction between the classes of eigenfunctions is applicable only for small values of $`\eta `$.

Consequently, the eigenfunctions can be labeled by multi-indices $`\mu =(n_x,n_y,n_z,r/l)`$: 3 quantum numbers $`n_{x,y,z}_+`$ and a flag $`r/l`$ denoting the class (right- or left-dominated, respectively) of the eigenfunction.

Now we want to consider the effect of one moving boundary. It is taken into account by substituting $`aa(t)`$ everywhere in the eigenmodes and -frequencies. Thereby a time dependence of the eigenfunctions as well as of the eigenfrequencies is introduced. This induces a non-vanishing coupling matrix $`M_{\mu \nu }(t)`$ as well as the frequency deviation $`\mathrm{\Delta }\mathrm{\Omega }_\mu ^2(t)`$. For small oscillations of the boundary

$`a(t)=a_0+ϵ(ba_0)\mathrm{sin}(\omega t)`$ (31)
with a small amplitude $`ϵ1`$ it will be useful to separate the time dependence using

$`M_{\mu \nu }(t)`$ $`=`$ $`\dot{a}(t){\displaystyle d^3r\frac{f_\mu }{a}(𝒓;t)f_\nu (𝒓;t)}`$ (32)
$`=`$ $`\dot{a}(t)m_{\mu \nu }(t).`$
The geometry factor $`m_{\mu \nu }(t)`$ is approximately constant $`m_{\mu \nu }(t)=m_{\mu \nu }+𝒪\left(ϵ\right)`$ in this case. Consequently, one is lead to

$`M_{\mu \nu }(t)=m_{\mu \nu }\dot{a}(t)+𝒪\left(ϵ^2\right).`$ (33)
Since the time-dependence of the right-dominated modes is less complicated than that of the left-dominated ones, it is advantageous to exploit the antisymmetry of $`M_{\mu \nu }`$ which also implies an antisymmetry of $`m_{\mu \nu }`$. For the following calculations the coupling of the lowest left-dominated mode $`\mu =(1,1,1,l)`$ to some right-dominated one $`\nu =(n_x,n_y,n_z,r)`$ will be of special relevance. The $`y`$ and $`z`$ integrations simply generate Kronecker symbols and therefore the geometry factor results as

$`m_{\mu ,\nu }`$ $`=`$ $`\delta _{1,n_y}\delta _{1,n_z}{\displaystyle _a^c}𝑑xf_\mu ^x{\displaystyle \frac{f_\nu ^x}{a}}`$ (34)
$`=`$ $`{\displaystyle \frac{\delta _{1,n_y}\delta _{1,n_z}n_x(1)^{n_x}\sqrt{{\displaystyle \frac{ba}{cb}}}{\displaystyle \frac{\mathrm{\Omega }_{n_xr}^x}{\mathrm{\Omega }_{1l}^x}}}{(cb)\mathrm{sin}\left(n_x\pi {\displaystyle \frac{ba}{cb}}\right)\left[n_x^2\left({\displaystyle \frac{ba}{cb}}\right)^21\right]}}\eta `$
$`+𝒪\left(\eta ^2\right)`$
$`=`$ $`𝒪\left(\eta \right).`$
### II.4 Canonical Quantization

Aiming at the calculation of possible particle creation effects (expectation values of particle number operators) it is convenient to introduce the creation and annihilation operators

$`\widehat{a}_\mu (t)={\displaystyle \frac{1}{\sqrt{2\mathrm{\Omega }_\mu ^0}}}\left(\mathrm{\Omega }_\mu ^0\widehat{Q}_\mu (t)+i\widehat{P}_\mu (t)\right),`$ (35)
obeying the usual bosonic equal time commutation relations

$`[\widehat{a}_\mu (t),\widehat{a}_\nu ^{}(t)]`$ $`=`$ $`\delta _{\mu \nu },`$
$`[\widehat{a}_\mu (t),\widehat{a}_\nu (t)]`$ $`=`$ $`0,`$
$`[\widehat{a}_\mu ^{}(t),\widehat{a}_\nu ^{}(t)]`$ $`=`$ $`0.`$ (36)
These operators diagonalize the free Hamiltonian

$`\widehat{H}_0={\displaystyle \underset{\mu }{}}\mathrm{\Omega }_\mu ^0\left(\widehat{a}_\mu ^{}(t)\widehat{a}_\mu (t)+{\displaystyle \frac{1}{2}}\right).`$ (37)
The following calculations will most conveniently be done in the interaction picture where the dynamics of an observable is governed by $`\widehat{H}_0`$

$`{\displaystyle \frac{d\widehat{Y}}{dt}}`$ $`=`$ $`i[\widehat{H}_0,\widehat{Y}]+\left({\displaystyle \frac{\widehat{Y}}{t}}\right)_{\mathrm{explicit}}.`$ (38)
For reasons of generality and to include finite temperature effects we describe the state of a quantum system by a statistical operator whose dynamics is determined by the von Neumann equation

$`{\displaystyle \frac{d\widehat{\rho }}{dt}}`$ $`=`$ $`i[\widehat{H}_I,\widehat{\rho }].`$ (39)
Note that this equation without any explicit time dependence $`(\widehat{\rho }/t)_{\mathrm{explicit}}`$ leads to an unitary time evolution, see also finitetemp  and Sec. IV.

In this picture the time dependence of the creation and annihilation operators turns out to be

$`\widehat{a}_\mu (t)=\widehat{a}_\mu e^{i\mathrm{\Omega }_\mu ^0t}.`$ (40)
However, this trivial time dependence gives rise to the possibility of parametric resonance which enhances the chances to verify the effect of particle creation experimentally. Further-on we will denote the initial creation and annihilation operators by $`\widehat{a}_\mu (0)=\widehat{a}_\mu `$. Note that in this picture the particle number operator $`\widehat{N}_\mu =\widehat{a}_\mu ^{}\widehat{a}_\mu `$ is time independent for all modes.

### II.5 Rotating Wave Approximation

In the interaction picture the time-evolution operator is given by

$`\widehat{U}(T,0)=\widehat{𝒯}_t\mathrm{exp}\left[i{\displaystyle _0^T}\left(\widehat{H}_I^S(t)+\widehat{H}_I^V(t)\right)𝑑t\right]`$ (41)
where $`\widehat{𝒯}_t`$ denotes time-ordering. If the interaction Hamiltonian $`\widehat{H}_I`$ leads to small corrections, the usual procedure is to apply perturbation theory via expanding the exponential. Evidently, this would imply that $`\widehat{U}`$ is close to the identity. On the other hand, in order to make an experimental verification of quantum radiation feasible, the time-evolution operator should deviate significantly from the identity. Therefore a different approximation needs to be found. For the case of parametric resonance this may be accomplished via applying the rotating wave approximation (RWA), see e.g. finitetemp ; law ; wuchuleung . Within this scenario the left boundary performs harmonic oscillations obeying $`a(t)=a_0+ϵ(ba_0)\mathrm{sin}(\omega t)`$ with the dimensionless amplitude<sup>1</sup><sup>1</sup>1For uniqueness we restrict ourselves to $`ϵ>0`$. $`ϵ1`$ and the external vibration frequency $`\omega `$ during the time interval $`[0,T]`$. This also implies an oscillating time-dependence of the frequency deviation and coupling matrix

$`\mathrm{\Delta }\mathrm{\Omega }_\mu ^2(t)`$ $`=`$ $`2\mathrm{\Omega }_\mu ^0{\displaystyle \frac{\mathrm{\Omega }_\mu ^0}{a_0}}(ba_0)ϵ\mathrm{sin}(\omega t)+𝒪\left(ϵ^2\right),`$ (42)
$`M_{\mu \nu }(t)`$ $`=`$ $`m_{\mu \nu }(ba_0)\omega ϵ\mathrm{cos}(\omega t)+𝒪\left(ϵ^2\right),`$ (43)
see also Sec. II.3. Together with the trivial time-dependence of the ladder operators (40) in the interaction picture this enables us to perform the RWA. Expanding the time evolution operator (41) into an infinite series one can treat the time-ordering terms in the following way: In analogy to

$`\widehat{𝒯}_t\left[\widehat{H}_I(t_1)\widehat{H}_I(t_2)\right]`$ $`=`$ $`\mathrm{\Theta }(t_2t_1)[\widehat{H}_I(t_2),\widehat{H}_I(t_1)]`$ (44)
$`+\widehat{H}_I(t_1)\widehat{H}_I(t_2)`$
one can rewrite all these terms to yield a multiple product of Hamiltonians without time-ordering and terms involving commutators with Heaviside step functions. In the resonance case, i.e., when $`\omega =2\mathrm{\Omega }_\mu ^0`$ the terms with commutators can yield nothing but strongly oscillating integrands finitetemp , as can also be seen by Fourier-expanding the involved Hamiltonians. With the duration of the perturbation being sufficiently long, i.e., with $`\omega T1`$, the contribution of these terms to (41) will be comparably small. As a consequence, time-ordering can be neglected to all orders within the RWA.

The remaining integrals – without time-ordering – factorize and can be resummated to yield an effective time evolution operator

$`\widehat{U}_{\mathrm{eff}}(T,0)\stackrel{\mathrm{RWA}}{=}\mathrm{exp}\left(i\widehat{H}_{\mathrm{eff}}^IT\right),`$ (45)
where the effective interaction Hamiltonian

$`\widehat{H}_{\mathrm{eff}}^IT\stackrel{\mathrm{RWA}}{=}{\displaystyle _0^T}\left[\widehat{H}_I^S(t)+\widehat{H}_I^V(t)\right]𝑑t`$ (46)
still has to be calculated. The above time integration involves many oscillating terms. Again, with the duration of the vibration being sufficiently long $`\omega T1`$, i.e., after many oscillations, the time-integrated interaction Hamiltonian can be approximated in the following way: Since the time average of purely oscillating terms is rather small compared to that of constant contributions we may neglect the former ones. As a result, in the series expansion of the time-evolution operator $`\widehat{U}(T,0)`$ only those terms where the oscillation of the ladder operators is compensated by the external vibrations – represented by $`\mathrm{\Delta }\mathrm{\Omega }_\mu ^2(t)`$ and $`M_{\mu \nu }(t)`$, respectively – will be kept. Strictly speaking, in the above equation terms of $`𝒪\left(ϵ^I(\omega T)^J\right)`$ are neglected by the RWA if $`I>J`$ holds. The terms with $`J=K`$ – i.e., exactly the terms in which the oscillations of the creation and annihilation operators (40) are compensated by the external time dependence $`\mathrm{\Delta }\mathrm{\Omega }_\mu ^2(t)`$ and $`M_{\mu \nu }(t)`$\] – will be kept. (Note that terms with $`J<K`$ do not occur.)

The general squeezing interaction Hamiltonian reads

$`\widehat{H}_I^S`$ $`=`$ $`{\displaystyle \underset{\mu }{}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{\Omega }_\mu ^0}{a_0}}ϵ(ba_0)\mathrm{sin}(\omega t)[(\widehat{a}_\mu )^2(t)+(\widehat{a}_\mu ^{})^2(t)`$ (47)
$`+\widehat{a}_\mu (t)\widehat{a}_\mu ^{}(t)+\widehat{a}_\mu ^{}(t)\widehat{a}_\mu (t)]+𝒪(ϵ^2).`$
Accordingly, within the RWA only the terms fulfilling the squeezing resonance condition, see also e.g. leakydodonov ; finitetemp ; quantumrad ; braginsky

$`\omega =2\mathrm{\Omega }_\mu ^0`$ (48)
will be kept. In general $`\mu `$ can also be a right-dominated mode, but note that in this case the effective squeezing Hamiltonian would be of $`𝒪\left(\eta ^2\right)`$, since according to (II.3) the right-dominated eigenfrequencies do not depend on $`a(t)`$ up to $`𝒪\left(\eta \right)`$. Therefore we will restrict ourselves to left-dominated modes $`\mu `$ and among those in particular to the lowest one, i.e., as commonly done we consider the case of fundamental resonance

$`\omega =2\mathrm{\Omega }_{111l}^0=2\mathrm{\Omega }_L^0.`$ (49)
From now on this mode will be abbreviated by the index $`L=(1,1,1,l)`$ throughout this publication. Consequently, by virtue of

$`{\displaystyle \frac{1}{T}}{\displaystyle _0^T}\mathrm{sin}(\omega t)e^{\pm i\omega t}𝑑t\stackrel{\mathrm{RWA}}{=}\pm {\displaystyle \frac{i}{2}}`$ (50)
an effective squeezing Hamiltonian can be derived

$`\widehat{H}_{\mathrm{eff}}^S=i\xi \left[(\widehat{a}_L^{})^2(\widehat{a}_L)^2\right],`$ (51)
where $`\xi `$ is given by

$`\xi ={\displaystyle \frac{1}{4}}ϵ\mathrm{\Omega }_L^0\left({\displaystyle \frac{\mathrm{\Omega }_L^{x0}}{\mathrm{\Omega }_L^0}}\right)^2.`$ (52)
Obviously $`\widehat{H}_{\mathrm{eff}}^S`$ is a generator for a squeezing operator for the mode $`L`$ with $`\xi `$ being the squeezing parameter.

The same procedure can be applied for the velocity interaction Hamiltonian which reads in the vibration case

$`\widehat{H}_I^V`$ $`=`$ $`{\displaystyle \frac{i}{2}}{\displaystyle \underset{\mu \nu }{}}\sqrt{{\displaystyle \frac{\mathrm{\Omega }_\mu ^0}{\mathrm{\Omega }_\nu ^0}}}m_{\mu \nu }ϵ\omega (ba_0)\mathrm{cos}(\omega t)\times `$ (53)
$`[\widehat{a}_\mu ^{}(t)\widehat{a}_\nu ^{}(t)+\widehat{a}_\mu ^{}(t)\widehat{a}_\nu (t)`$
$`\widehat{a}_\mu (t)\widehat{a}_\nu ^{}(t)\widehat{a}_\mu (t)\widehat{a}_\nu (t)]+𝒪(ϵ^2).`$
However, the occurrence of inter-mode couplings now results in a different resonance condition (see also finitetemp ; quantumrad )

$`\omega `$ $`=`$ $`\left|\mathrm{\Omega }_\mu ^0\pm \mathrm{\Omega }_\nu ^0\right|.`$ (54)
Depending on the frequency spectrum of the cavity under consideration this resonance condition might be fulfilled by several pairs $`\mu \nu `$, but here we will assume for simplicity that only one such pair exists. Though in any case via

$`{\displaystyle \frac{1}{T}}{\displaystyle _0^T}\mathrm{cos}(\omega t)e^{\pm i\omega t}𝑑t\stackrel{\mathrm{RWA}}{=}{\displaystyle \frac{1}{2}}`$ (55)
an effective velocity Hamiltonian can be derived, two major distinctions should be made.

#### $``$ coupling

$`\omega =\mathrm{\Omega }_1^0+\mathrm{\Omega }_2^0`$. In this case one yields an effective velocity Hamiltonian given by

$`\widehat{H}_{\mathrm{eff}}^V`$ $`=`$ $`{\displaystyle \frac{i}{8}}(\sqrt{{\displaystyle \frac{\mathrm{\Omega }_1^0}{\mathrm{\Omega }_2^0}}}\sqrt{{\displaystyle \frac{\mathrm{\Omega }_2^0}{\mathrm{\Omega }_1^0}}})m_{12}ϵ\omega (ba_0)\times `$ (56)
$`\left(\widehat{a}_1^{}\widehat{a}_2^{}\widehat{a}_1\widehat{a}_2\right),`$
which is a non-diagonal multi-mode squeezing Hamiltonian. Note that if one wants to fulfill squeezing and velocity resonance conditions simultaneously ($`2\mathrm{\Omega }_L^0=\omega =\mathrm{\Omega }_1^0+\mathrm{\Omega }_2^0`$), the number of possible combinations reduces significantly, since then the velocity resonance condition can not be fulfilled by two distinct left-dominated modes. For reasons of brevity we do not consider this case here.

#### $``$ coupling

$`\omega =\mathrm{\Omega }_2^0\mathrm{\Omega }_1^0`$. Here the resulting effective velocity Hamiltonian does not resemble a squeezing but a hopping operator

$`\widehat{H}_{\mathrm{eff}}^V`$ $`=`$ $`{\displaystyle \frac{i}{8}}(\sqrt{{\displaystyle \frac{\mathrm{\Omega }_1^0}{\mathrm{\Omega }_2^0}}}+\sqrt{{\displaystyle \frac{\mathrm{\Omega }_2^0}{\mathrm{\Omega }_1^0}}})m_{12}ϵ\omega (ba_0)\times `$ (57)
$`\left(\widehat{a}_1^{}\widehat{a}_2\widehat{a}_1\widehat{a}_2^{}\right).`$
This coupling is of special interest since if one does not insist on simultaneously fulfilling both resonance conditions – parametric resonance might perhaps still be induced by lower external frequencies $`\omega =\mathrm{\Omega }_2^0\mathrm{\Omega }_1^0<2\mathrm{\Omega }_L^0`$ whose generation would be simpler in an experiment, see also Sec. V.3. In the case of simultaneously fulfilling both conditions several combinations may arise.

* The frequencies $`\mathrm{\Omega }_i^0`$ both belong to either right- or left-dominated modes. In dalvit  it has been shown that for an ideal cavity with special dimensions, e.g. a cubic one, a strong inter-mode coupling can occur. In that case a much smaller exponential particle creation rate has been found. Therefore in view of an experimental verification this case is counterproductive and not considered here. Instead we propose a cavity with transcendental ratios of the dimensions such that there is no resonant inter-mode coupling of similar-dominated modes. In any case such a coupling would certainly require large quantum numbers of the involved modes.
* The frequency $`\mathrm{\Omega }_2^0`$ represents a right-dominated mode and $`\mathrm{\Omega }_1^0`$ some left-dominated mode, respectively. The lowest possible right-dominated frequency $`\mathrm{\Omega }_2=\mathrm{\Omega }_R^0`$ would then be obtained when $`\mathrm{\Omega }_1^0=\mathrm{\Omega }_L^0`$. As an example, this case will be considered here. We want to stress that the used methods are nevertheless applicable to any possible combination of couplings.

Note that the situation would be completely different in 1+1 space-time dimensions where – due to the equidistant spectrum – the velocity term always contributes, see e.g. finitetemp ; dodonov1 ; dodonov2 ; dodonov3 ; dodonov4 ; dodonov5 ; dodonov6 . The coupling right-dominated mode fulfilling $`\mathrm{\Omega }_{n_x,n_y,n_z,r}=3\mathrm{\Omega }_L^0`$ will further-on be denoted with the index $`R=(n_x,n_y,n_z,r)`$. Accordingly, in our considerations the velocity Hamiltonian reads

$`\widehat{H}_{\mathrm{eff}}^V`$ $`=`$ $`i\chi \left(\widehat{a}_L^{}\widehat{a}_R\widehat{a}_L\widehat{a}_R^{}\right)`$ (58)
with

$`\chi `$ $`=`$ $`{\displaystyle \frac{1}{4}}ϵ\mathrm{\Omega }_L^0\left(\sqrt{{\displaystyle \frac{\mathrm{\Omega }_R^0}{\mathrm{\Omega }_L^0}}}+\sqrt{{\displaystyle \frac{\mathrm{\Omega }_L^0}{\mathrm{\Omega }_R^0}}}\right)(ba_0)m_{L,R}`$ (59)
being the velocity parameter of the system. Since $`\chi =𝒪\left(ϵ\mathrm{\Omega }_L^0m_{L,R}\right)=𝒪\left(ϵ\mathrm{\Omega }_L^0\eta \right)`$ \[see also (34)\] it follows that $`\chi /\xi =𝒪\left(\eta \right)1`$ in the limiting case of a nearly perfectly reflecting mirror<sup>2</sup><sup>2</sup>2Note that for inter-mode coupling between modes of the same region – the corresponding scenario has been assumed in dalvit  for an ideal cavity where one only has one region – the velocity parameter $`\chi `$ would be of $`𝒪\left(1\right)`$ instead of $`𝒪\left(\eta \right)`$.. This hopping operator is consistent with the visual picture of a semitransparent mirror.

### II.6 Response Theory

We assume our system to be initially in a state of thermal equilibrium that can be described by the canonical ensemble

$`\widehat{\rho }(t=0)=\widehat{\rho }_0={\displaystyle \frac{\mathrm{exp}\left(\beta \widehat{H}_0\right)}{\mathrm{Tr}\left\{\mathrm{exp}\left(\beta \widehat{H}_0\right)\right\}}}`$ (60)
with $`\beta `$ denoting the initial inverse temperature. The expectation value of an explicitly time-independent operator $`\widehat{Y}`$ at time $`t=T`$ is given by

$`Y(T)`$ $`=`$ $`\mathrm{Tr}\left\{\widehat{Y}\widehat{\rho }(T)\right\}`$ (61)
$`=`$ $`\mathrm{Tr}\{\widehat{Y}\widehat{𝒯}_t\mathrm{exp}(i{\displaystyle _0^T}\widehat{H}_I(t_1)dt_1)\times `$
$`\widehat{\rho }_0\widehat{𝒯}_t^{}\mathrm{exp}(+i{\displaystyle _0^T}\widehat{H}_I(t_2)dt_2)\}`$
where $`\widehat{𝒯}_t^{}`$ denotes the anti-chronological operator (anti-time ordering) – for a more involved discussion see e.g. finitetemp .

If the interaction Hamiltonian represented a small correction, one could expand the time-evolution operator into a perturbation series yielding a series expansion for $`P(T)`$. However, for the resonance case this procedure is not justified: In the rotating wave approximation in Sec. II.5 the interaction Hamiltonian simplified to

$`\widehat{H}_{\mathrm{eff}}=i\xi \left[(\widehat{a}_L^{})^2(\widehat{a}_L)^2\right]+i\chi \left[\widehat{a}_L^{}\widehat{a}_R\widehat{a}_L\widehat{a}_R^{}\right]`$ (62)
which implies for the time-evolution operator

$`\widehat{U}(T,0)\stackrel{\mathrm{RWA}}{=}\mathrm{exp}\left(i\widehat{H}_{\mathrm{eff}}^IT\right).`$ (63)
The whole expression for computing the expectation value of an operator now becomes much simpler

$`Y(T)`$ $`\stackrel{\mathrm{RWA}}{=}`$ $`\mathrm{Tr}\{\widehat{Y}\mathrm{exp}(i(\widehat{H}_{\mathrm{eff}}^S+\widehat{H}_{\mathrm{eff}}^V)T)\widehat{\rho }_0\times `$ (64)
$`\mathrm{exp}(+i(\widehat{H}_{\mathrm{eff}}^S+\widehat{H}_{\mathrm{eff}}^V)T)\},`$
but since the correction $`\widehat{H}_{\mathrm{eff}}^ST`$ is not small in the case of interest above expression is still not practical for applying perturbation theory. Exploiting the smallness of the velocity Hamiltonian it will prove useful to separate the two Hamiltonians. This can be achieved with the ansatz

$`\mathrm{exp}\left(i(\widehat{H}_{\mathrm{eff}}^S+\widehat{H}_{\mathrm{eff}}^V)\tau \right)=\mathrm{exp}\left(i\widehat{H}_{\mathrm{eff}}^S\tau \right)\widehat{\sigma }(\tau )`$ (65)
with $`\widehat{\sigma }`$ being an auxiliary operator. Differentiation with respect to $`\tau `$ yields a differential equation that can be solved for $`\widehat{\sigma }`$ using the initial condition $`\widehat{\sigma }(0)=\mathrm{𝟏}`$. Introducing the parameter ordering $`\widehat{𝒯}_\tau `$ in analogy to time ordering ($`\widehat{𝒯}_t`$) the solution for $`\widehat{\sigma }`$ can be cast into the form

$`\widehat{\sigma }(T)=\widehat{𝒯}_\tau \left[\mathrm{exp}\left(i{\displaystyle _0^T}\widehat{}_{\mathrm{eff}}^V(\tau )𝑑\tau \right)\right].`$ (66)
Here the squeezed effective velocity Hamiltonian has been introduced

$`\widehat{}_{\mathrm{eff}}^V(\tau )=\mathrm{exp}\left(+i\widehat{H}_{\mathrm{eff}}^S\tau \right)\widehat{H}_{\mathrm{eff}}^V\mathrm{exp}\left(i\widehat{H}_{\mathrm{eff}}^S\tau \right)`$ (67)
which is now dependent on the parameter $`\tau `$. Further-on we shall denote squeezed operators by calligraphy letters. By inserting above equations into the expectation value (64) one yields

$`Y(T)`$ $`\stackrel{\mathrm{RWA}}{=}`$ $`\mathrm{Tr}\{\widehat{𝒴}(T)\widehat{𝒯}_\tau \mathrm{exp}(i{\displaystyle _0^T}\widehat{}_{\mathrm{eff}}^V(\tau _1)d\tau _1)\times `$ (68)
$`\rho _0\widehat{𝒯}_\tau ^{}\mathrm{exp}(+i{\displaystyle _0^T}\widehat{}_{\mathrm{eff}}^V(\tau _2)d\tau _2)\}.`$
Please note that in this representation also the observables are squeezed

$`\widehat{𝒴}(T)=\mathrm{exp}\left(+i\widehat{H}_{\mathrm{eff}}^ST\right)\widehat{Y}\mathrm{exp}\left(i\widehat{H}_{\mathrm{eff}}^ST\right)`$ (69)
but here using the physical perturbation time $`T`$. We will refer to this picture as the squeezing interaction picture. Unfortunately the parameter ordering is reintroduced by this procedure but as the advantage of these manipulations we are now able to expand the expectation value $`Y(T)`$ into a perturbation series with powers of $`\widehat{}_{\mathrm{eff}}^V`$. Keeping only terms to second order one finds

$`Y(T)`$ $`=`$ $`\mathrm{Tr}\left\{\widehat{𝒴}(T)\widehat{\rho }_0\right\}`$ (70)
$`+`$ $`\mathrm{Tr}\left\{\widehat{𝒴}(T)[\widehat{\rho }_0,i{\displaystyle 𝑑\tau _1\widehat{}_{\mathrm{eff}}^V(\tau _1)}]\right\}`$
$`+`$ $`\mathrm{Tr}\left\{\widehat{𝒴}(T){\displaystyle 𝑑\tau _1\widehat{}_{\mathrm{eff}}^V(\tau _1)\widehat{\rho }_0𝑑\tau _2\widehat{}_{\mathrm{eff}}^V(\tau _2)}\right\}`$
$``$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}\left\{\widehat{𝒴}(T)\widehat{𝒯}_\tau {\displaystyle 𝑑\tau _1\widehat{}_{\mathrm{eff}}^V(\tau _1)𝑑\tau _2\widehat{}_{\mathrm{eff}}^V(\tau _2)\widehat{\rho }_0}\right\}`$
$``$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}\left\{\widehat{𝒴}(T)\widehat{\rho }_0\widehat{𝒯}_\tau ^{}{\displaystyle 𝑑\tau _1\widehat{}_{\mathrm{eff}}^V(\tau _1)𝑑\tau _2\widehat{}_{\mathrm{eff}}^V(\tau _2)}\right\}`$
$`+`$ $`𝒪\left((\widehat{}_{\mathrm{eff}}^V)^3\right)`$
which is now a practical expression for calculating expectation values.

## III The Quadratic Response

### III.1 Squeezing

According to the results of Sec. II.6 in the squeezing interaction picture both particle number operator and the effective velocity Hamiltonian have to be squeezed. The squeezing operator

$`\widehat{S}(\tau )=\mathrm{exp}\left(+i\widehat{H}_{\mathrm{eff}}^S\tau \right)=\mathrm{exp}\left(\xi \left[(\widehat{a}_L)^2(\widehat{a}_L^{})^2\right]\tau \right)`$ (71)
implies the following well-known transformation rules (see e.g. mandelwolf )

$`\widehat{b}_L(\tau )`$ $`=`$ $`\widehat{S}(\tau )\widehat{a}_L\widehat{S}^{}(\tau )`$ (72)
$`=`$ $`\widehat{a}_L\mathrm{cosh}(2\xi \tau )+\widehat{a}_L^{}\mathrm{sinh}(2\xi \tau ),`$
$`\widehat{b}_L^{}(\tau )`$ $`=`$ $`\widehat{S}(\tau )\widehat{a}_L^{}\widehat{S}^{}(\tau )`$ (73)
$`=`$ $`\widehat{a}_L^{}\mathrm{cosh}(2\xi \tau )+\widehat{a}_L\mathrm{sinh}(2\xi \tau ),`$
which can also be envisaged as a Bogoliubov transformation of the ladder operators. Due to the commutation relations (II.4) other modes than the fundamental resonance mode $`L`$ are not affected by squeezing. Inserting the above expressions into the effective velocity Hamiltonian (58) one can easily derive its squeezed counterpart

$`\widehat{}_{\mathrm{eff}}^V(\tau )`$ $`=`$ $`\widehat{S}(\tau )\widehat{H}_{\mathrm{eff}}^V\widehat{S}^{}(\tau )`$ (74)
$`=`$ $`i\chi \mathrm{cosh}(2\xi \tau )\left[\widehat{a}_L^{}\widehat{a}_R\widehat{a}_L\widehat{a}_R^{}\right]`$
$`+i\chi \mathrm{sinh}(2\xi \tau )\left[\widehat{a}_L\widehat{a}_R\widehat{a}_L^{}\widehat{a}_R^{}\right].`$
Note that the squeezed effective velocity Hamiltonian is now dependent on the parameter $`\tau `$ and that it is still of $`𝒪\left(\chi \right)=𝒪\left(\eta \right)`$ which justifies a perturbation-like treatment.

The same can be done for the particle number operators where again only the fundamental resonance mode $`\widehat{N}_L=\widehat{a}_L^{}\widehat{a}_L`$ is affected

$`\widehat{𝒩}_L(T)`$ $`=`$ $`\widehat{S}(T)\widehat{N}_L\widehat{S}^{}(T)=\widehat{S}(T)\widehat{a}_L^{}\widehat{S}^{}(T)\widehat{S}(T)\widehat{a}_L\widehat{S}^{}(T)`$ (75)
$`=`$ $`\left[1+2\mathrm{sinh}^2(2\xi T)\right]\widehat{a}_L^{}\widehat{a}_L`$
$`+{\displaystyle \frac{1}{2}}\mathrm{sinh}(4\xi T)\left[(\widehat{a}_L^{})^2+(\widehat{a}_L)^2\right]`$
$`+\mathrm{sinh}^2(2\xi T).`$
For brevity we will denote the hyperbolic functions by

$`𝒞(T)`$ $`=`$ $`\mathrm{cosh}(2\xi T),`$
$`𝒮(T)`$ $`=`$ $`\mathrm{sinh}(2\xi T)`$ (76)
throughout this paper.

### III.2 Expectation Values

Since we are mainly interested in the phenomenon of quantum radiation and thus therefore in the calculation of the cavity particle content after a time $`T`$ when the disturbance has ended. Due to the dynamical disturbance the system leaves the thermodynamic equilibrium, see also finitetemp . The quadratic response of the expectation value of the particle number operator can be calculated via substituting $`\widehat{Y}\widehat{N}`$ in equation (70). However, since the initial statistical operator $`\widehat{\rho }_0`$ involves arbitrarily high powers in $`\widehat{H}_0`$, it is practical to rewrite the expression obtained from (70). Utilizing the invariance of the trace under cyclic permutation and the property of time ordering

$`\{\widehat{A}(t),\widehat{B}(t^{})\}_+`$ $`=`$ $`\widehat{A}(t)\widehat{B}(t^{})+\widehat{B}(t^{})\widehat{A}(t)`$ (77)
$`=`$ $`\widehat{𝒯}\left[\widehat{A}(t)\widehat{B}(t^{})\right]`$
$`+\widehat{𝒯}^{}\left[\widehat{A}(t)\widehat{B}(t^{})\right]`$
the expectation value of interest can be cast into the more convenient form

$`N(T)`$ $`=`$ $`\mathrm{Tr}\left\{\widehat{𝒩}(T)\widehat{\rho }_0\right\}`$ (78)
$`+`$ $`\mathrm{Tr}\left\{[i{\displaystyle 𝑑\tau _1\widehat{}_{\mathrm{eff}}^V(\tau _1)},\widehat{𝒩}(T)]\widehat{\rho }_0\right\}`$
$`+`$ $`\mathrm{Tr}\{\widehat{\rho }_0{\displaystyle }d\tau _2\widehat{}_{\mathrm{eff}}^V(\tau _2)[\widehat{𝒩}(T),{\displaystyle }d\tau _1\widehat{}_{\mathrm{eff}}^V(\tau _1)]`$
$`+`$ $`{\displaystyle \frac{1}{2}}\widehat{\rho }_0[\widehat{𝒯}_\tau {\displaystyle }d\tau _1\widehat{}_{\mathrm{eff}}^V(\tau _1){\displaystyle }d\tau _2\widehat{}_{\mathrm{eff}}^V(\tau _2),\widehat{𝒩}(T)]\}`$
$`+`$ $`𝒪\left((\widehat{}_{\mathrm{eff}}^V)^3\right).`$
This form is now suitable for evaluating the traces since all commutators only concern a finite number of creation or annihilation operators. Thus the quadratic response can be brought in relation with the initial particle contents of the cavity (Bose-Einstein distributions)

$`N_L^0`$ $`=`$ $`\mathrm{Tr}\left\{\widehat{a}_L^{}\widehat{a}_L\widehat{\rho }_0\right\}={\displaystyle \frac{1}{e^{\beta \mathrm{\Omega }_L^0}1}},`$ (79)
$`N_R^0`$ $`=`$ $`\mathrm{Tr}\left\{\widehat{a}_R^{}\widehat{a}_R\widehat{\rho }_0\right\}={\displaystyle \frac{1}{e^{\beta \mathrm{\Omega }_R^0}1}},`$ (80)
where $`\beta `$ stands for the initial inverse temperature of the system. These mean occupation numbers incorporate the whole temperature dependence of the quantum radiation – as long as the back-reaction of the field on the moving mirror can be neglected. As we shall show later, the lowest order term ($`\eta =0`$) is in agreement with the results of an ideal cavity, as was considered for example in finitetemp ; dalvit . Also, since $`\widehat{}_{\mathrm{eff}}^V`$ contains only odd powers of creation and annihilation operators for a single mode, the linear response vanishes. Generally, every trace involving an odd power of ladder operators vanishes and $`\widehat{𝒩}`$ as well as $`\widehat{\rho }_0`$ do only contain even powers. Hence the last trace in Eq. (78) constitutes the quadratic answer. In contrast to an ideal cavity the terms with time ordering are here especially important since they will be found to produce leading order terms.

### III.3 Particle Creation

Using the squeezed operators (74) and (75) it is now straightforward to compute the commutators and the traces in the expectation value (78). As a result one finds for particles in the fundamental resonance mode $`L`$

$`N_L(T)`$ $`=`$ $`𝒮^2(T)`$ (81)
$`+\left[1+2𝒮^2(T)\right]N_L^0`$
$`+{\displaystyle \frac{\chi ^2}{4\xi ^2}}\left[3𝒞^2(T)2𝒞(T)12\xi T𝒮(2T)\right]`$
$`+{\displaystyle \frac{\chi ^2}{4\xi ^2}}\left[4𝒞^2(T)2𝒞(T)24\xi T𝒮(2T)\right]N_L^0`$
$`+{\displaystyle \frac{\chi ^2}{4\xi ^2}}\left[2𝒞^2(T)2𝒞(T)\right]N_R^0`$
$`+𝒪\left(\eta ^3\right).`$
As was anticipated, the lowest order term $`𝒮^2(T)+[1+2𝒮^2(T)]N_L^0`$ is in agreement with the results obtained in finitetemp  for an ideal cavity. The linear response (in $`\eta `$) vanishes. It might be of interest that the leading terms $`T𝒮(2T)`$ in the quadratic answer stem from the time-ordering which is therefore very important. One can see that at long disturbance times $`T`$ these leading terms show the failure of the quadratic approximation since the particle number would become negative at some point. This is due to the fact that (78) is a perturbation series in $`_0^T\widehat{H}_{\mathrm{eff}}^V(t)𝑑t`$ which will always become large at some time $`T`$. This problem can only be solved by including all orders in $`\eta `$, see also Sec. V.

Of course (78) can also be applied to the corresponding coupling right-dominated mode (whose particle number operator is invariant under squeezing) where one finds

$`N_R(T)`$ $`=`$ $`N_R^0`$ (82)
$`+{\displaystyle \frac{\chi ^2}{4\xi ^2}}\left[2𝒞^2(T)2𝒞(T)+1\right]`$
$`+{\displaystyle \frac{\chi ^2}{4\xi ^2}}\left[2𝒞^2(T)2𝒞(T)\right]N_L^0`$
$`+{\displaystyle \frac{\chi ^2}{4\xi ^2}}\left[2𝒞(T)+2\right]N_R^0`$
$`+𝒪\left(\eta ^3\right).`$
Again the linear answer is vanishing. For $`\eta =0`$ there would not be any created particles in the reservoir due to the dynamical Casimir effect corresponding to a perfect internal mirror.

It is remarkable that the coefficient of $`N_L^0`$ in $`N_R(T)`$ equals the coefficient of $`N_R^0`$ in $`N_L(T)`$. As we shall see in Sec.V.2, this feature is preserved to all orders in $`\eta `$.

## IV The Master Equation Approach

In this section it is our aim to derive the associated master equation for an effective statistical operator accounting for the left leaky sub-cavity or left-dominated modes, respectively. So far (3+1) dimensional vibrating leaky cavities have only been treated in different setups – see e.g. leakydalvit  – where the vibrating mirror is understood as a (quantized) harmonic oscillator coupled to the cavity field (the reservoir) or with master equations adequate rather for stationary systems – see e.g. leakydodonov . It was assumed in leakydodonov  that these master equations could also be applied when one of the boundaries was moving. The possibility of limitations to that procedure as well as the need for a rigorous derivation of the master equation for resonantly excited systems have already been expressed in leakydodonov . We want to derive such an equation starting from first principles. As a test we will also solve the obtained master equation and recalculate the quadratic answer for the left mode particles to compare with the previous results of Sec. III.3. To obtain a master equation we will closely follow the derivation given in mandelwolf .

### IV.1 Derivation of a Master equation

Throughout this section we will deploy the squeezing interaction picture where not only the time dependence induced by $`\widehat{H}_0`$ but also the dependence resulting from $`\widehat{H}_{\mathrm{eff}}^S`$ is determining the operator time evolution has already been proposed in Sec. II.6. In this picture the time evolution of the statistical operator is governed by a modified von Neumann equation

$`{\displaystyle \frac{\widehat{\rho }(t)}{t}}=i[\widehat{}_{\mathrm{eff}}^V(t),\widehat{\rho }(t)]=i\widehat{𝔏}(t)\widehat{\rho }(t).`$ (83)
The above equation defines the action of the Liouvillian super operator $`\widehat{𝔏}`$ on $`\widehat{\rho }`$ (see also Ref. fick ). By defining the projection super operator $`\widehat{𝔓}`$ via

$`\widehat{𝔓}\widehat{Y}=\widehat{\rho }_R(0)\mathrm{Tr}_R\left\{\widehat{Y}\right\}`$ (84)
for all observables $`\widehat{Y}`$, where $`\mathrm{Tr}_R`$ means taking the trace solely over the right dominated modes we can introduce a reduced density operator accounting for the left-dominated modes only

$`\widehat{\rho }_L(t)`$ $`=`$ $`\mathrm{Tr}_R\left\{\widehat{\rho }(t)\right\}.`$ (85)
Combining above equations it can be shown mandelwolf  that the dynamics of the full statistical operator $`\widehat{\rho }`$ is governed by the Zwanzig master equation

$`\widehat{𝔓}{\displaystyle \frac{\widehat{\rho }(t)}{t}}`$ $`=`$ $`\widehat{𝔓}\widehat{𝔏}(t){\displaystyle _0^t}\widehat{𝔘}(t,t^{})(1\widehat{𝔓})\widehat{𝔏}(t^{})\widehat{𝔓}\widehat{\rho }(t^{})𝑑t^{}`$ (86)
$`i\widehat{𝔓}\widehat{𝔏}(t)\widehat{𝔘}(t,0)(1\widehat{𝔓})\widehat{\rho }(0)`$
$`i\widehat{𝔓}\widehat{𝔏}(t)\widehat{𝔓}\widehat{\rho }(t)`$
where

$`\widehat{𝔘}(t,t^{})=\mathrm{exp}\left(i(1\widehat{𝔓}){\displaystyle _t^{}^t}\widehat{𝔏}(t^{\prime \prime })𝑑t^{\prime \prime }\right)`$ (87)
is the reduced time-evolution super operator. The Zwanzig master equation is exact to all orders in $`\eta `$ but usually too complicated to be solved. However, assuming an initial thermal equilibrium state and taking into account that initially our system and reservoir do not interact (no correlations) it can be simplified considerably:

* In analogy to the argumentation concerning the vanishing of the linear response in Sec. III.2 it follows that

$`\mathrm{Tr}_R\left\{\widehat{}_{\mathrm{eff}}^V\widehat{\rho }_R(0)\right\}=0`$ (88)
since $`\widehat{}_{\mathrm{eff}}^V`$ contains only odd and $`\widehat{\rho }_R(0)`$ only even powers of the creation and annihilation operators for the mode $`R`$. This can equivalently be written as

$`\widehat{𝔓}\widehat{𝔏}(t)\widehat{𝔓}\widehat{\rho }(t)=0.`$ (89)
* In our setup the initially stationary system (stationary walls) does not permit interactions between system and reservoir, since both $`M_{\mu \nu }(t_0)`$ and $`\mathrm{\Delta }\mathrm{\Omega }^2(t_0)`$ will vanish. Consequently, assuming thermal equilibrium, system and reservoir initially constitute independent subsystems which cannot be correlated, i.e., the initial statistical operator of the cavity modes factorizes

$`\widehat{\rho }_0=\widehat{\rho }(0)=\widehat{\rho }_L(0)\widehat{\rho }_R(0),`$ (90)
hence one finds (with $`\mathrm{Tr}_R\{\widehat{\rho }_R\}=1`$)

$`(1\widehat{𝔓})\widehat{\rho }(0)=0.`$ (91)

These assumptions yield a simplified Zwanzig master equation

$`{\displaystyle \frac{\widehat{𝔓}\widehat{\rho }(t)}{t}}=\widehat{𝔓}\widehat{𝔏}(t){\displaystyle _0^t}\widehat{𝔘}(t,t^{})\widehat{𝔏}(t^{})\widehat{𝔓}\widehat{\rho }(t^{})𝑑t^{}`$ (92)
which is exact but still too complicated to be solved.

In order to gain a solvable equation we will apply further approximations:

* Born approximation

Since $`\widehat{𝔏}=𝒪\left(\eta \right)`$ one can approximate the reduced time-evolution operator via $`\widehat{𝔘}(t,t^{})=\mathrm{𝟏}+𝒪\left(\eta \right)`$. This neglects terms of $`𝒪\left(\eta ^3\right)`$ if inserted into (92) and yields

$`{\displaystyle \frac{\widehat{𝔓}\widehat{\rho }(t)}{t}}=\widehat{𝔓}\widehat{𝔏}(t){\displaystyle _0^t}\widehat{𝔏}(t^{})\widehat{𝔓}\widehat{\rho }(t^{})𝑑t^{}+𝒪\left(\eta ^3\right).`$ (93)
By employing the reduced density operator $`\widehat{\rho }_L(t)=\mathrm{Tr}_R\left\{\widehat{\rho }(t)\right\}`$ one can equivalently write

$`{\displaystyle \frac{\widehat{\rho }_L(t)}{t}}`$ $`=`$ $`\mathrm{Tr}_R\left\{\widehat{𝔏}(t){\displaystyle _0^t}\widehat{𝔏}(t^{})\widehat{\rho }_R(0)\widehat{\rho }_L(t^{})𝑑t^{}\right\}`$ (94)
$`+𝒪\left(\eta ^3\right).`$
This equation governs the time evolution of the effective statistical operator $`\widehat{\rho }_L`$ accounting for the left cavity.
* Markov approximation

The retardation in equation (94), i.e., the occurrence of $`\widehat{\rho }_L(t^{})`$, complicates the calculation of $`\widehat{\rho }_L(t)`$. Iterative application of (94) implies that $`\widehat{\rho }_L(t^{})=\widehat{\rho }_L(t)+𝒪\left(\eta ^2\right)`$. Accordingly, we apply the Markov approximation, which is also known as short memory approximation, simply by replacing $`\widehat{\rho }_L(t^{})\widehat{\rho }_L(t)`$ on the right hand side. Since $`\widehat{𝔏}=𝒪\left(\eta \right)`$ we thereby neglect terms of $`𝒪\left(\eta ^4\right)`$ and yield the Born-Markov master equation

$`{\displaystyle \frac{\widehat{\rho }_L(t)}{t}}`$ $`=`$ $`\mathrm{Tr}_R\left\{{\displaystyle _0^t}\widehat{𝔏}(t)\widehat{𝔏}(t^{})\widehat{\rho }_R(0)\widehat{\rho }_L(t)𝑑t^{}\right\}`$ (95)
$`+𝒪\left(\eta ^3\right),`$
thus having maintained the level of accuracy.

Using the definition of the Liouville operator $`\widehat{𝔏}`$ in Eq. (83) one can equivalently write

$`{\displaystyle \frac{\widehat{\rho }_L(t)}{t}}`$ $`=`$ $`+\mathrm{Tr}_R\left\{\widehat{}_{\mathrm{eff}}^V(t)\widehat{\rho }_R(0)\widehat{\rho }_L(t){\displaystyle _0^t}\widehat{}_{\mathrm{eff}}^V(t^{})𝑑t^{}\right\}`$ (96)
$`\mathrm{Tr}_R\left\{\widehat{}_{\mathrm{eff}}^V(t){\displaystyle _0^t}\widehat{}_{\mathrm{eff}}^V(t^{})𝑑t^{}\widehat{\rho }_R(0)\widehat{\rho }_L(t)\right\}`$
$`+\mathrm{h}.\mathrm{c}.+𝒪\left(\eta ^3\right).`$
Finally, having evaluated both traces and after having performed the $`t^{}`$ integrations with the aid of (74), one yields the following master equation

$`{\displaystyle \frac{\widehat{\rho }_L(t)}{t}}`$ $`=`$ $`f_1(t)\left[2\widehat{a}_L^{}\widehat{\rho }_L(t)\widehat{a}_L\widehat{a}_L\widehat{a}_L^{}\widehat{\rho }_L(t)\widehat{\rho }_L(t)\widehat{a}_L\widehat{a}_L^{}\right]`$ (97)
$`+f_2(t)\left[2\widehat{a}_L\widehat{\rho }_L(t)\widehat{a}_L^{}\widehat{a}_L^{}\widehat{a}_L\widehat{\rho }_L(t)\widehat{\rho }_L(t)\widehat{a}_L^{}\widehat{a}_L\right]`$
$`+f_3(t)\left[\widehat{a}_L^{}\widehat{\rho }_L(t)\widehat{a}_L^{}+\widehat{a}_L\widehat{\rho }_L(t)\widehat{a}_L\right]`$
$`f_4(t)\left[(\widehat{a}_L^{})^2\widehat{\rho }_L(t)+\widehat{\rho }_L(t)(\widehat{a}_L)^2\right]`$
$`f_5(t)\left[(\widehat{a}_L)^2\widehat{\rho }_L(t)+\widehat{\rho }_L(t)(\widehat{a}_L^{})^2\right]`$
$`+𝒪\left(\eta ^3\right)`$
where the functions $`f_i(t)`$ are given by

$`f_1(t)`$ $`=`$ $`{\displaystyle \frac{\chi ^2}{2\xi }}𝒮(t)\left\{𝒞(t)(2N_R^0+1)N_R^01\right\},`$
$`f_2(t)`$ $`=`$ $`{\displaystyle \frac{\chi ^2}{2\xi }}𝒮(t)\left\{𝒞(t)(2N_R^0+1)N_R^0\right\},`$
$`f_3(t)`$ $`=`$ $`{\displaystyle \frac{\chi ^2}{2\xi }}\left\{𝒞^2(t)+𝒮^2(t)𝒞(t)\right\}(2N_R^0+1),`$
$`f_4(t)`$ $`=`$ $`{\displaystyle \frac{\chi ^2}{2\xi }}\left\{\left[𝒞^2(t)+𝒮^2(t)𝒞(t)\right]N_R^0+𝒞^2(t)𝒞(t)\right\},`$
$`f_5(t)`$ $`=`$ $`{\displaystyle \frac{\chi ^2}{2\xi }}\left\{\left[𝒞^2(t)+𝒮^2(t)𝒞(t)\right]N_R^0+𝒮^2(t)\right\}.`$ (98)
Via averaging over the degrees of freedom of the reservoir and by applying the Born-Markov approximation we have now rigorously derived a differential equation for an effective statistical operator $`\widehat{\rho }_L(t)`$ accounting for the leaky cavity. This effective statistical operator obeys a non-unitary time-evolution (changing entropy). There are several possibilities to check the obtained master equation: As the simplest tests one can verify that the time evolution preserves<sup>3</sup><sup>3</sup>3Note however, that under certain conditions (e.g. $`N_R^0=0`$) the positive definiteness of the density operator may not be preserved – if one applies the master equation (97) for large times beyond its region of validity. the hermiticity and the trace of $`\widehat{\rho }_L`$.

A better indication for a correct master equation is the fact that if one takes the limit of no squeezing, i.e., in this coupling $`\xi 0`$, the resulting equation corresponds to a harmonic oscillator coupled to a thermal bath: With

$`\underset{\xi 0}{lim}f_1(t)`$ $`=`$ $`\chi ^2tN_R^0,`$ (99)
$`\underset{\xi 0}{lim}f_2(t)`$ $`=`$ $`\chi ^2t(N_R^0+1),`$ (100)
$`\underset{\xi 0}{lim}f_{i=3,4,5}(t)`$ $`=`$ $`0`$ (101)
one arrives at a simplified equation

$`{\displaystyle \frac{\widehat{\rho }_L}{t}}`$ $`\stackrel{\xi 0}{=}`$ $`\gamma _D{\displaystyle \frac{N_R^0}{2}}\left[2\widehat{a}_L^{}\widehat{\rho }_L\widehat{a}_L\widehat{a}_L\widehat{a}_L^{}\widehat{\rho }_L\widehat{\rho }_L\widehat{a}_L\widehat{a}_L^{}\right]`$ (102)
$`+\gamma _D{\displaystyle \frac{N_R^0+1}{2}}\left[2\widehat{a}_L\widehat{\rho }_L\widehat{a}_L^{}\widehat{a}_L^{}\widehat{a}_L\widehat{\rho }_L\widehat{\rho }_L\widehat{a}_L^{}\widehat{a}_L\right]`$
$`+𝒪\left(\eta ^3\right).`$
Apart from the time dependence of the damping coefficient $`\gamma _D=2\chi ^2t`$ above equation is exactly the well-known master equation for a harmonic oscillator coupled to a thermal bath, see e.g. amohandbook .

The time-dependence of $`\gamma _D`$ is a remnant of the dynamic master equation describing the time-dependent system in the unphysical limit $`\xi 0`$. However, in order to have a stronger indication for the correctness of our ansatz we want to solve the master equation (97) explicitly.

### IV.2 Approximate solution of the master equation

So far we have neglected terms of $`𝒪\left(\eta ^3\right)`$. The functions $`f_i(t)`$ are already of $`𝒪\left(\eta ^2\right)`$ which makes it possible to maintain the level of accuracy by applying the additional approximation $`\widehat{\rho }_L(t)\widehat{\rho }_L(0)`$ on the right hand side of equation (97), which could also be envisaged as an additional Markov approximation. Accordingly, one is now able to yield a solution for $`\widehat{\rho }_L`$

$`\widehat{\rho }_L(T)`$ $`=`$ $`\widehat{\rho }_L(0)`$ (103)
$`+`$ $`F_1(T)\left[2\widehat{a}_L^{}\widehat{\rho }_L(0)\widehat{a}_L\widehat{a}_L\widehat{a}_L^{}\widehat{\rho }_L(0)\widehat{\rho }_L(0)\widehat{a}_L\widehat{a}_L^{}\right]`$
$`+`$ $`F_2(T)\left[2\widehat{a}_L\widehat{\rho }_L(0)\widehat{a}_L^{}\widehat{a}_L^{}\widehat{a}_L\widehat{\rho }_L(0)\widehat{\rho }_L(0)\widehat{a}_L^{}\widehat{a}_L\right]`$
$`+`$ $`F_3(T)\left[\widehat{a}_L^{}\widehat{\rho }_L(0)\widehat{a}_L^{}+\widehat{a}_L\widehat{\rho }_L(0)\widehat{a}_L\right]`$
$``$ $`F_4(T)\left[(\widehat{a}_L^{})^2\widehat{\rho }_L(0)+\widehat{\rho }_L(0)(\widehat{a}_L)^2\right]`$
$``$ $`F_5(T)\left[(\widehat{a}_L)^2\widehat{\rho }_L(0)+\widehat{\rho }_L(0)(\widehat{a}_L^{})^2\right]`$
$`+`$ $`𝒪\left(\eta ^3\right)`$
with $`F_i(T)=_0^Tf_i(t)𝑑t`$. Given this effective statistical operator for the leaky cavity one is now able to calculate the number of created particles in all left-dominated modes. Note that for considering right-dominated modes one would have to derive a statistical operator for the reservoir.

### IV.3 Particle Creation

Since we were working in the squeezing interaction picture – where the observables have to be squeezed – the expectation value of the particle number operator reads

$`N_L(T)=\mathrm{Tr}_L\left\{\widehat{𝒩}_L(T)\widehat{\rho }_L(T)\right\}.`$ (104)
Other left-dominated modes than the fundamental resonance mode $`L`$ are trivial to solve: Due to the commutation relations (II.4) their ladder operators commute with those of the resonance mode $`L`$. This implies (due to the invariance of the trace under cyclic permutations) that all higher order traces must vanish and one just yields the trivial result of their initial occupation numbers. Inserting the approximate reduced density operator obtained in Eq. (103) as well as $`\widehat{𝒩}_L`$ into the above equation, one can see immediately that zeroth and first order in $`\eta `$ agree with the previous results but showing this for the second order is a bit tedious. After some algebra one finally finds complete agreement which the previous result found in Eq. (81) of Sec. III.3 thus giving a strong indication for the validity of our master equation within the RWA approach.

### IV.4 Comparison with other results

In leakydodonov  the effects of losses are taken into account by a generalized version of the simple master equation ansatz

$`{\displaystyle \frac{d\widehat{\rho }}{dt}}=i[\widehat{\rho },\widehat{H}]+{\displaystyle \frac{\gamma _D}{2}}\left[2\widehat{a}\widehat{\rho }\widehat{a}^{}\widehat{a}^{}\widehat{a}\widehat{\rho }\widehat{\rho }\widehat{a}^{}\widehat{a}\right].`$ (105)
However, as we have observed in the previous calculations, this master equation does not adequately describe the leaky cavity under consideration:

* It is restricted to the case where the initial state of the reservoir is just the vacuum state and therefore does not include temperature effects. This has been taken into account in leakydodonov .
* In addition, even the master equation for an harmonic oscillator in a thermal bath amohandbook

$`{\displaystyle \frac{\widehat{\rho }}{t}}`$ $`=`$ $`i[\widehat{\rho },\widehat{H}]+{\displaystyle \frac{\overline{n}}{2}}\gamma _D\left[2\widehat{a}^{}\widehat{\rho }\widehat{a}\widehat{a}\widehat{a}^{}\widehat{\rho }\widehat{\rho }\widehat{a}\widehat{a}^{}\right]`$ (106)
$`+{\displaystyle \frac{\overline{n}+1}{2}}\gamma _D\left[2\widehat{a}\widehat{\rho }\widehat{a}^{}\widehat{a}^{}\widehat{a}\widehat{\rho }\widehat{\rho }\widehat{a}^{}\widehat{a}\right]`$
cannot be assumed to describe the system correctly. Even if one identifies the Hamiltonian $`\widehat{H}`$ in above equation with the effective squeezing Hamiltonian $`\widehat{H_{\mathrm{eff}}^S}`$ this master equation goes along with serious problems since the Markov approximation is not justified anymore. This complication reflects the inherent dynamic character of our system. As we have shown in Sec. IV.1 the complete master equation resembles above equation only in the limit of no squeezing $`\xi 0`$ – see Eq. (102) – and even then with a time-dependent damping constant $`\gamma _D`$.

Instead, the complete master equation (97) displays more similarities to one in a squeezed thermal bath where one has to replace the parameters by time-dependent functions. Accordingly, the dynamical system under consideration is described properly only by an explicitly time-dependent master equation.

Potential limitations to Eq. (105) have already been anticipated in leakydodonov .

## V The Non-Perturbative Approach

The previous results in Sec. III.3 and Sec. IV.3 have not been able to explain the behavior of the system in the limit of a long-lasting disturbance. The leading order term $`T𝒮(2T)`$ in (81) has a negative sign which would lead to negative particle numbers for large disturbance times $`T`$. This problem can only be solved by including all orders in $`\eta `$. In this section we present a non-perturbative approach within the RWA which enables a convenient calculation of expectation values using computer algebra systems schaller . As a further advantage we want to mention that it can in principle be generalized in a straightforward way to the case of more than just two coupling modes as was assumed in Sec. II.5.

### V.1 Time Evolution

Application of the RWA in Sec. II.5 yielded the effective time-evolution operator $`\widehat{U}_{\mathrm{eff}}=\mathrm{exp}(i\widehat{H}_{\mathrm{eff}}^IT)`$ with the effective interaction Hamiltonian (62). We want to calculate the expectation value of an observable $`\widehat{Y}`$

$`Y(T)`$ $`=`$ $`\mathrm{Tr}\left\{\widehat{Y}e^{i\widehat{H}_{\mathrm{eff}}^IT}\widehat{\rho }_0e^{+i\widehat{H}_{\mathrm{eff}}^IT}\right\}`$ (107)
$`=`$ $`\mathrm{Tr}\left\{e^{+i\widehat{H}_{\mathrm{eff}}^IT}\widehat{Y}e^{i\widehat{H}_{\mathrm{eff}}^IT}\widehat{\rho }_0\right\}.`$
In contrast to the previous sections here the full time dependence is shifted back on the operator $`\widehat{Y}`$. Since $`\widehat{Y}`$ can be expressed using creation and annihilation operators and due to unitarity of the time-evolution operator one just has to find a solution for the full time dependence of the ladder operators which is given by

$`\widehat{a}_\sigma (T)=e^{+i\widehat{H}_{\mathrm{eff}}^IT}\widehat{a}_\sigma e^{i\widehat{H}_{\mathrm{eff}}^IT}.`$ (108)
The above expression requires special care to evaluate, since $`\widehat{H}_{\mathrm{eff}}^I`$ is not a pure squeezing generator. In our considerations the effective interaction Hamiltonian is time-independent which does not necessarily hold in general. To preserve generality we will therefore introduce an auxiliary parameter $`\vartheta `$ while keeping the time $`T`$ fixed. This enables us to write

$`\widehat{a}_\sigma (\vartheta )=e^{+i\widehat{H}_{\mathrm{eff}}^IT\vartheta }\widehat{a}_\sigma e^{i\widehat{H}_{\mathrm{eff}}^IT\vartheta }.`$ (109)
Obviously we are interested in $`\widehat{a}_\sigma (T)=\widehat{a}_\sigma (\vartheta =1)`$. To this end we define a 4 dimensional column vector (see also dalvit )

$`\underset{¯}{\overset{^}{x}}(\vartheta )=\left(\begin{array}{c}\widehat{a}_L(\vartheta )\\ \widehat{a}_L^{}(\vartheta )\\ \widehat{a}_R(\vartheta )\\ \widehat{a}_R^{}(\vartheta )\end{array}\right).`$ (114)
Since $`\widehat{H}_{\mathrm{eff}}^I`$ does not depend on $`\vartheta `$, one finds

$`{\displaystyle \frac{d\underset{¯}{\overset{^}{x}}}{d\vartheta }}=iT[\widehat{H}_{\mathrm{eff}},\underset{¯}{\overset{^}{x}}(\vartheta )]=T\underset{¯}{A}\underset{¯}{\overset{^}{x}}(\vartheta )`$ (115)
where $`\underset{¯}{A}`$ is a number-valued 4 by 4 matrix acting on $`\underset{¯}{x}`$. This form can always be reached if the effective Hamiltonian is quadratic: The commutation relations (II.4) lead to a linear combination of creation and annihilation operators that can always be written as a number-valued matrix $`\underset{¯}{A}`$ acting on $`\underset{¯}{x}`$. Since $`\underset{¯}{A}`$ is independent of $`\vartheta `$ the solution is obtained via

$`\underset{¯}{\overset{^}{x}}(\vartheta )=\mathrm{exp}\left(\underset{¯}{A}T\vartheta \right)\underset{¯}{\overset{^}{x}}(0)`$ (116)
and hence

$`\underset{¯}{\overset{^}{x}}(1)=\mathrm{exp}\left(\underset{¯}{A}T\right)\underset{¯}{\overset{^}{x}}(0)=\underset{¯}{U}(T)\underset{¯}{\overset{^}{x}}(0).`$ (117)
Thus the whole problem reduces to a calculation of the time-evolution matrix $`\underset{¯}{U}(T)=\mathrm{exp}(\underset{¯}{A}T)`$. In the present case the structure of $`\widehat{H}_{\mathrm{eff}}`$ in Eq. (62) implies a very simple form of $`\underset{¯}{A}`$

$`\underset{¯}{A}=\left(\begin{array}{cccc}0& 2\xi & \chi & 0\\ 2\xi & 0& 0& \chi \\ \chi & 0& 0& 0\\ 0& \chi & 0& 0\end{array}\right).`$ (122)
The four eigenvalues of $`\underset{¯}{A}`$ are given by

$`\lambda _1`$ $`=`$ $`\xi +\sqrt{\xi ^2\chi ^2},`$
$`\lambda _2`$ $`=`$ $`\xi \sqrt{\xi ^2\chi ^2},`$
$`\lambda _3`$ $`=`$ $`\xi +\sqrt{\xi ^2\chi ^2},`$
$`\lambda _4`$ $`=`$ $`\xi \sqrt{\xi ^2\chi ^2}.`$ (123)
For reasons of brevity we shall omit the full listing of the matrix $`\underset{¯}{U}(T)=\mathrm{exp}(\underset{¯}{A}T)`$ – it can easily be calculated using some computer algebra system. Note that the exponential matrix $`\underset{¯}{U}(T)`$ is positive definite for all $`T`$ and thus will not exhibit the problems associated with the extrapolation of the used approximations of sections III and IV beyond their range of validity. In order to calculate expectation values one just needs the matrix elements of $`\underset{¯}{U}(T)`$. This becomes evident considering the time evolution of the new annihilation and creation operators $`\underset{¯}{\overset{^}{x}}(T)=\underset{¯}{U}(T)\underset{¯}{\overset{^}{x}}(0)`$, i.e. $`\widehat{x}_i(T)=_{j=1}^4U_{ij}(T)\widehat{x}_j(0)`$. Therefore the expectation values of particle number operators of the resonance modes $`\widehat{N}_L(T)=\widehat{x}_2(T)\widehat{x}_1(T)`$ and $`\widehat{N}_R(T)=\widehat{x}_4(T)\widehat{x}_3(T)`$ can be calculated simply by insertion of $`\underset{¯}{\overset{^}{x}}(T)`$. After evaluation of the remaining traces containing only the initial creation and annihilation operators $`\widehat{x}_i(0)`$ one finds the full response function to be a combination of matrix elements of $`\underset{¯}{U}(T)`$

$`N_L(T)`$ $`=`$ $`(U_{12}U_{21}+U_{14}U_{23})`$ (124)
$`+(U_{11}U_{22}+U_{12}U_{21})N_L^0`$
$`+(U_{13}U_{24}+U_{14}U_{23})N_R^0`$
and

$`N_R(T)`$ $`=`$ $`(U_{41}U_{32}+U_{34}U_{43})`$ (125)
$`+(U_{42}U_{31}+U_{41}U_{32})N_L^0`$
$`+(U_{33}U_{44}+U_{43}U_{34})N_R^0.`$
### V.2 Particle Creation

With the full knowledge of $`\underset{¯}{U}`$ we are now in a position to state the full response function of the particle number operator. Having inserted the matrix elements of $`\underset{¯}{U}`$ into Eq. (124) one finds after performing some simplifications

$`N_L(T)`$ $`=`$ $`{\displaystyle \frac{1}{4\left(\xi ^2\chi ^2\right)}}\times `$
$`\{\xi \mathrm{cosh}\left[2T(\xi +\sqrt{\xi ^2\chi ^2})\right][\xi +\sqrt{\xi ^2\chi ^2}]`$
$`+\xi \mathrm{cosh}\left[2T\left(\xi \sqrt{\xi ^2\chi ^2}\right)\right]\left[\xi \sqrt{\xi ^2\chi ^2}\right]`$
$`2\chi ^2\mathrm{cosh}\left[2T\xi \right]\stackrel{}{}2(\xi ^2\chi ^2)\}`$
$`+{\displaystyle \frac{N_L^0}{2(\xi ^2\chi ^2)}}\times `$
$`\{\xi \mathrm{cosh}\left[2T(\xi +\sqrt{\xi ^2\chi ^2})\right][\xi +\sqrt{\xi ^2\chi ^2}]`$
$`+\xi \mathrm{cosh}\left[2T\left(\xi \sqrt{\xi ^2\chi ^2}\right)\right]\left[\xi \sqrt{\xi ^2\chi ^2}\right]`$
$`\chi ^2\mathrm{cosh}\left[2T\xi \right][\mathrm{cosh}\left(2T\sqrt{\xi ^2\chi ^2}\right)+1]\}`$
$`+{\displaystyle \frac{N_R^0}{2(\xi ^2\chi ^2)}}\times `$
$`\left\{\chi ^2\mathrm{cosh}\left[2T\xi \right]\left[\mathrm{cosh}\left(2T\sqrt{\xi ^2\chi ^2}\right)1\right]\right\}.`$
This result is valid to all orders in $`\chi /\xi `$ or $`\eta `$, respectively. To show consistency with the results obtained in Sec. III.3 and Sec. IV.3 we expanded the above expression around $`\chi /\xi =0`$ up to second order and found complete agreement with Eq. (81)! However, even for large values of $`\chi /\xi =1/2`$ the quadratic approximation is a rather good one – provided that the duration of the disturbance $`T`$ is not extremely large – as one can see in Fig. 4.

Doing the same calculations for the corresponding right-dominated mode one finds as a result

$`N_R(T)`$ $`=`$ $`{\displaystyle \frac{1}{4(\xi ^2\chi ^2)}}\times `$
$`\{\xi \mathrm{cosh}\left[2T(\xi +\sqrt{\xi ^2\chi ^2})\right][\xi \sqrt{\xi ^2\chi ^2}]`$
$`+\xi \mathrm{cosh}\left[2T\left(\xi \sqrt{\xi ^2\chi ^2}\right)\right]\left[\xi +\sqrt{\xi ^2\chi ^2}\right]`$
$`2\chi ^2\mathrm{cosh}\left[2T\xi \right]\stackrel{}{}2(\xi ^2\chi ^2)\}`$
$`+{\displaystyle \frac{N_L^0}{2(\xi ^2\chi ^2)}}\times `$
$`\left\{\chi ^2\mathrm{cosh}\left[2T\xi \right]\left[\mathrm{cosh}\left(2T\sqrt{\xi ^2\chi ^2}\right)1\right]\right\}`$
$`+{\displaystyle \frac{N_R^0}{2(\xi ^2\chi ^2)}}\times `$
$`\{\xi \mathrm{cosh}\left[2T(\xi +\sqrt{\xi ^2\chi ^2})\right][\xi \sqrt{\xi ^2\chi ^2}]`$
$`+\xi \mathrm{cosh}\left[2T\left(\xi \sqrt{\xi ^2\chi ^2}\right)\right]\left[\xi +\sqrt{\xi ^2\chi ^2}\right]`$
$`\chi ^2\mathrm{cosh}\left[2T\xi \right][\mathrm{cosh}\left(2T\sqrt{\xi ^2\chi ^2}\right)+1]\}`$
where the remarkable agreement of coefficients of $`N_R^0`$ in $`N_L`$ and of $`N_L^0`$ in $`N_R`$ as was already noticed in Sec. III.3 is preserved for all orders in $`\eta `$. These terms fit the classical picture of particle transportation through the leaky membrane where the particle flux is proportional to the number of particles on the other side. Again, expanding around $`\chi /\xi =0`$ up to second order yields exact agreement with (82). Accordingly, also outside the leaky cavity particles are produced due to the dynamical Casimir effect, see also Fig. 5.

Note that at least the quadratic answer is necessary to treat particle creation effects outside the leaky cavity. It is still valid that finite-temperature corrections will enhance the pure vacuum phenomenon of particle production by several orders of magnitude. (For a direct comparison see Fig. 7 in Sec. VIII.)

### V.3 Further Remarks

We have derived a complete solution for the effective interaction Hamiltonian (62) which is valid to all orders in $`\chi /\xi =𝒪\left(\eta \right)`$. As an illustration we consider a case outside our initial intentions where $`\chi /\xi `$ also assumes large values, e.g. $`\chi /\xi 1`$. In this case the arguments of the hyperbolic functions in (V.2) and (V.2) will receive an imaginary part. The arising imaginary parts of $`N(T)`$ cancel as they have to because $`\widehat{N}`$ is a physical observable. Thus one finds that if the velocity parameter $`\chi `$ exceeds the critical value $`\chi \xi `$ the particle occupation number of the resonance modes versus the vibration time will exhibit oscillations. Of course, for the case of a nearly perfectly reflecting mirror inside this scenario is completely unrealistic since then $`\chi /\xi =𝒪\left(\eta \right)`$ will be relatively small. However, this case is not at all academic: If the label $`R`$ stood for a left-dominated mode – which is the case we excluded in our considerations so far and whose equivalent for ideal cavities has been considered in dalvit $`\chi /\xi `$ may very well become large, since $`m_{LR}`$ would then be of $`𝒪\left(1\right)`$.

Similar oscillations of the particle number were also found in the case of strong inter-mode coupling in an ideal cavity dalvit .

Note that with $`\chi /\xi =𝒪\left(\eta \right)`$ this also leads to an upper bound for the internal mirror transmittance $`\eta `$ above which (corresponding to a highly transparent mirror) one finds oscillations that correspond to inter-mode coupling rather than to system-reservoir coupling. From another perspective this phenomenon could also be envisaged as follows: Starting with an ideal cavity whose original dimensions do not permit inter-mode coupling one can insert a highly transparent mirror ($`\eta 1`$). This mirror in turn detunes the ideal cavity in such a way that it now permits inter-mode couplings as well.

It is remarkable that in Fig. 6 the phase of the two modes is shifted: When $`N_L(T)`$ is at its maximum, then $`N_R(T)`$ is at its minimum and vice versa. This fits nicely with the picture of mode hopping mediated by the inter-mode coupling $`\chi `$. One even observes a decrease of the particle number in the $`L`$-mode for small times. When defining an effective temperature finitetemp  this would correspond to an effective cooling of the $`L`$-mode. An extreme case of this consideration would be the limit of no squeezing, i.e. $`\xi =0`$. This would correspond to the possibility (see also Sec. II.5) of not fulfilling the squeezing but the velocity resonance condition. Performing the limit $`\xi 0`$ everywhere in equations (V.2) and (V.2) one would find pure oscillations of the particle numbers and no exponential growth at all. This case is therefore counterproductive for an experimental verification. Note however, that this is different for the case of $``$-coupling. The consistency with the earlier results leads to the conclusion that our approach was justified and the full response function should describe the rate of particle production correctly within the RWA.

Please note that the described procedure also holds for more than just two coupling modes: if one has e.g. $`n`$ modes fulfilling the resonance conditions given in Sec. II.5, the formalism still holds and one will have to define a $`2n`$ dimensional vector $`\underset{¯}{x}`$. Of course then $`n`$ creation and $`n`$ annihilation operators of these resonance modes will be contained in the Hamiltonian and therefore also $`\underset{¯}{A}`$ as well as $`\underset{¯}{U}`$ will be $`2n`$ by $`2n`$ matrices. The calculations will simply become more involved but can certainly be performed, e.g. by means of computer algebra systems.

## VI Detuning

So far we did assume an exact fulfillment of the resonance conditions, i.e. the vibration of the left cavity wall did match exactly twice the fundamental resonance frequency $`\mathrm{\Omega }_L^0`$. However, in real situations one will of course have to deal with deviations from this desired external frequency since it will not be possible to match it with arbitrary precision. In addition, the back-reaction of the created quanta might cause the external vibration frequency to change.

Consequently, we will now discuss the detuned situation – where $`\omega `$ assumes slightly off-resonant values. For a review of detuning effects see e.g. leakydodonov ; dalvit ; dodonovdetuned . It has been shown in the literature that there exist threshold values for the detuning, above which the exponential creation of particles disappears.

Unfortunately the RWA used in our previous considerations cannot simply be generalized to this situation. For a slight deviation from the resonance conditions in equations (49) and (54) the terms with time-ordering – see section II.5 – are no longer negligible in this way.

We will consider slightly off-resonant situations, where the external vibration frequency does not match the fundamental resonance exactly

$`\omega =2\mathrm{\Omega }_L^0(1+\delta ),`$ (128)
where $`\delta `$ denotes a small (dimensionless) deviation $`\delta 1`$. However, if one considers such a variance it is only consequent to include a possible discrepancy of the coupling resonance as well, cf. dodonovdetuned

$`\mathrm{\Omega }_R^0=\mathrm{\Omega }_L^0(3+\mathrm{\Delta }),`$ (129)
where $`\mathrm{\Delta }1`$ denotes the deviation of the coupling right-dominated mode from the $``$-coupling resonance condition with the fundamental resonance mode. We will adapt the multiple scale analysis (MSA) as proposed in dalvit  to our scenario of a leaky cavity, see also dodonovdetuned . For this purpose we restrict to the results, since the steps in dalvit ; dodonovdetuned  can strictly be followed – see also appendix A. The main difference in these considerations is that we use the eigenfunction system of subsection II.3 instead of those of an ideal cavity and that we assume the additional deviation (129) – see also dodonovdetuned .

In analogy to subsection V.1 one obtains a matrix $`\underset{¯}{A}^{}`$ governing the time evolution of the ladder operators

$`\underset{¯}{\overset{^}{x}}^{}=\left(\begin{array}{c}\widehat{a}_L\\ \widehat{a}_L^{}\\ \widehat{a}_R/\sqrt{3}\\ \widehat{a}_R^{}/\sqrt{3}\end{array}\right).`$ (134)
The creation of quanta will only be exponential – and thus noticeable – if at least one of the eigenvalues of this matrix

$`\underset{¯}{A}^{}\left(\begin{array}{cccc}i\mathrm{\Omega }_L^0\delta & 2\xi & \sqrt{3}\chi & 0\\ 2\xi & i\mathrm{\Omega }_L^0\delta & 0& \sqrt{3}\chi \\ \chi /\sqrt{3}& 0& i\mathrm{\Omega }_L^0(3\delta \mathrm{\Delta })& 0\\ 0& \chi /\sqrt{3}& 0& i\mathrm{\Omega }_L^0(3\delta \mathrm{\Delta })\end{array}\right),`$
does have a positive real part. Note that the slight disagreement between the above matrix and the one given in Ref. dodonovdetuned  is caused by the usage of a different phase ($`\mathrm{sin}`$ instead of $`\mathrm{cos}`$). With the abbreviations

$`𝒰=8\xi ^24\chi ^2+12(\mathrm{\Omega }_L^0)^2\delta \mathrm{\Delta }2(\mathrm{\Omega }_L^0)^2\mathrm{\Delta }^220(\mathrm{\Omega }_L^0)^2\delta ^2,`$
and

$`𝒱`$ $`=`$ $`+16\xi ^2(\xi ^2\chi ^2)`$
$`+64(\mathrm{\Omega }_L^0)^2\delta ^2[\xi ^2+\chi ^2+(\mathrm{\Omega }_L^0)^2\delta ^2]`$
$`+(\mathrm{\Omega }_L^0)^2\mathrm{\Delta }^2[8\xi ^2+4\chi ^2+52(\mathrm{\Omega }_L^0)^2\delta ^2+(\mathrm{\Omega }_L^0)^2\mathrm{\Delta }^2]`$
$`4(\mathrm{\Omega }_L^0)^2\delta \mathrm{\Delta }[12\xi ^2+8\chi ^2+24(\mathrm{\Omega }_L^0)^2\delta ^2+3(\mathrm{\Omega }_L^0)^2\mathrm{\Delta }^2],`$
the eigenvalues of the above matrix read (cf. dalvit ; dodonovdetuned )

$`\lambda =\pm {\displaystyle \frac{1}{2}}\sqrt{𝒰\pm 2\sqrt{𝒱}}.`$ (137)
As a consistency check we may set $`\delta =\mathrm{\Delta }=0`$ where the eigenvalues reduce to the ones given in Eq. (V.1). On the other hand, for $`\chi =0`$ one recovers the usual result of pure squeezing in an ideal cavity $`\lambda _+=\pm \sqrt{4\xi ^4\mathrm{\Omega }_L^2\delta ^2}`$.

Note that in contrast to dalvit ; dodonovdetuned  the inter-mode coupling and thus the parameter $`\chi `$ is very small $`\chi \xi `$. This enables us to expand the quantities $`𝒰`$ and $`𝒱`$ into powers of $`\chi `$. The condition for a real eigenvalue $`𝒰+2\sqrt{𝒱}>0`$ reads

$`\mathrm{\Omega }_L^2\delta ^2`$ $`<`$ $`4\xi ^22{\displaystyle \frac{4\xi ^2+\mathrm{\Omega }_L^2\delta \mathrm{\Delta }4\mathrm{\Omega }_L^2\delta ^2}{4\xi ^2\mathrm{\Omega }_L^2\delta ^2+\mathrm{\Omega }_L^2(3\delta \mathrm{\Delta })^2}}\chi ^2`$ (138)
$`+𝒪\left(\chi ^4\right).`$
Since $`\chi `$ is supposed to be small $`\chi \xi `$ one obtains a significant contribution only if $`3\delta \mathrm{\Delta }`$ and also in this case merely in the immediate vicinity of the critical value $`2\xi =\mathrm{\Omega }_L^0\delta `$. Consequently, the presence of an internal mirror of moderate quality do not drastically modify the threshold

$`\delta <{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\mathrm{\Omega }_L^x}{\mathrm{\Omega }_L^0}}\right)^2ϵ`$ (139)
for exponential particle creation. However, we would like to emphasize that the shifts of the eigenfrequencies of the cavity due to the partly permeable internal wall must be taken into account, see also section IX below.

## VII Summary

We have considered a massless scalar quantum field inside a leaky cavity modeled by means of a dispersive mirror. For the case of the lossy cavity vibrating at twice the fundamental resonance frequency we derived an effective Hamiltonian using the rotating wave approximation. Within the framework of response theory the magnitude of particle creation due to the dynamical Casimir effect was calculated. Furthermore we deduced the corresponding master equation via applying the Born-Markov approximation. We found a discrepancy to the master equations used so far (see leakydodonov ) to describe oscillating leaky cavities. We also applied a non-perturbative approach for the explicit calculation of the time evolution starting from the effective Hamiltonian. All these methods were found to lead to consistent results. In addition, the effects of a detuned external vibration frequency need to be taken into account.

It turned out that for the case of moderately low transmission coefficients (or sufficient quality factors) the rate of created particles is almost the same as for ideal cavities. The squeezing of the fundamental resonance mode as well as the strong enhancement of particle production due to the dynamical Casimir effect are preserved in the presence of transparent mirrors.

## VIII Discussion

In order to illustrate reasonable magnitudes let us specify the relevant parameters: A cavity with a typical size of $`\mathrm{\Lambda }=1`$ cm would have a fundamental resonance frequency of $`\mathrm{\Omega }_L^0150`$ GHz i.e., the corresponding coupling right dominated mode must have a frequency of about $`\mathrm{\Omega }_R^0=3\mathrm{\Omega }_L^0450`$ GHz. According to Ref. leakydodonov  we assume a dimensionless vibration amplitude $`ϵ10^8`$. Consequently, in order to create a significant amount of particles one would have to sustain the external oscillations over an interval of several milliseconds. At room temperature $`1/\beta 300`$ K one finds the initial particle occupation numbers to be $`N_L^0240`$ and $`N_R^080`$. Using the above values the squeezing parameter determines to $`\xi 150\mathrm{Hz}`$.

As the quality factor $`Q`$ of a resonator is defined as jackson

$`Q=2\pi {\displaystyle \frac{\mathrm{energy}\mathrm{in}\mathrm{cavity}}{\mathrm{energy}\mathrm{loss}\mathrm{per}\mathrm{period}}},`$ (140)
one finds as a classical estimate yields for our system

$`Q={\displaystyle \frac{2\pi }{\left|𝒯\right|^2}}=2\pi \left(1+\left({\displaystyle \frac{\gamma }{\mathrm{\Omega }_L^x}}\right)^2\right)=𝒪\left({\displaystyle \frac{1}{\eta ^2}}\right).`$ (141)
$`𝒯`$ denotes the transmission amplitude through the internal dispersive mirror and $`\gamma `$ was defined in Sec. II. Assuming a $`Q`$-factor of $`Q10^8`$ leakydodonov  (and references therein) this would imply for the corresponding perturbation parameter $`\eta =\mathrm{\Omega }_L^x/\gamma =𝒪\left(10^4\right)`$. With these values, a reasonable velocity parameter could be given by $`\chi 2\mathrm{mHz}`$.

The particle content of the leaky cavity is depicted in Fig. 7.

## IX Conclusion

According to the above considerations it is necessary to vibrate several milliseconds in order to produce measurable effects. As already stated, a cavity at finite temperature might even be advantageous – provided the cavity is still nearly ideal at its characteristic thermal wavelength. However, even after only one millisecond ($`10^8`$ periods) a classical estimate based on a quality factor of $`Q=10^8`$ would indicate drastic energy losses. On the other hand, our calculations based on a complete quantum treatment show that the effects of losses are almost negligible compared to the rate of particle creation as long as $`\eta 1`$. This leads to the conclusion that lower cavity quality factors than proposed in leakydodonov , e.g. $`Q=10^6`$ \[implying $`\eta =𝒪\left(10^3\right)`$\] would already completely suffice to justify our approximations schaller . Such quality factors are within the reach of the current experimental status. Of course our calculations are based on the assumption that the larger cavity – including both the reservoir and the leaky cavity – is perfectly conducting. The error made by this presumption is of $`𝒪\left(Q^2\right)`$ and therefore certainly negligible. Consequently, the experimental verification of the dynamical Casimir effect could be facilitated by a configuration where the vibrating cavity is enclosed by a slightly larger one as is demonstrated in Fig. 8.

A further important result schaller  of our investigations is the shift of the cavity eigenfrequencies (II.3) of $`𝒪\left(\eta \right)=𝒪\left(10^3\right)`$ which needs to be taken into account in order to make an experimental observation of quantum radiation feasible.

## X Outlook

### X.1 Multi-Mode Coupling

In section II.5 we assumed that only one of the right-dominated modes fulfills the resonance condition for the velocity Hamiltonian, i.e., exactly two modes are coupled. If, however, the reservoir cavity becomes larger, the spacing between different levels of its spectrum decreases so that eventually more than just one right-dominated mode begin to couple – at least within the range of detuning. In this case the effective velocity Hamiltonian would constitute a sum of single two-mode coupling Hamiltonians as the one in Eq. (58) – but accounting for different right modes $`R_1`$, $`R_2`$, etc. As we have observed in Fig. 4 in Sec. V.2, the quadratic answer is completely sufficient for reasonable values of $`\chi /\xi `$. Inserting the aforementioned sum of Hamiltonians into the quadratic answer one observes that the mixing terms vanish. Since the effective velocity Hamiltonian only contains odd powers of the creation and annihilation operators per mode, one can only obtain a non-vanishing trace if it involves two operators of the same (right) mode. Therefore the quadratic answer also decomposes into a sum of contributions each accounting for one right-dominated mode. Hence we expect the general structure of $`N_L(T)`$ to persist – just substitute $`\chi ^2_i\chi _i^2`$ in the leading contributions<sup>4</sup><sup>4</sup>4Of course, multi-mode coupling can also be taken care of by the non-perturbative approach of Sec. V. The required effort, however, will increase very fast, since the corresponding matrix grows with the number of coupling modes.. In order to ensure the applicability of the perturbative treatment the number $`n`$ of coupling right-dominated modes has to be small enough to satisfy $`_i\chi _i^2/\xi ^2n\chi ^2/\xi ^21`$.

Typically the spacing between two neighboring right-dominated modes is of $`𝒪\left(1/L\right)`$, i.e., the inverse of the characteristic length of the reservoir. On the other hand, the width of the resonance peak is of $`𝒪\left(1/T\right)`$. Consequently, if the duration of the disturbance $`T`$ (1 ms) exceeds the characteristic length of the reservoir (which is the case for $`L<10^5`$ m) then $`n`$ is certainly small enough.

### X.2 Electromagnetic Field

So far, we have considered a noninteracting, massless, and neutral scalar field. The next step could be to extend the calculations to the electromagnetic field. In this case several new difficulties arise:

1. The boundary conditions cannot just simply be described by Dirichlet (or Neumann) conditions. Especially for moving walls their form will be more complicated due to Ampere’s law (mixing of $`𝑬`$ and $`𝑩`$).
2. As the electromagnetic field is a gauge theory, one has to eliminate the unphysical degrees of freedom in order to quantize it. Again, for dynamic external conditions this requires special care, see e.g. quantumrad .
3. The different polarizations of photons need to be taken into account which are of special interest concerning the fulfillment of the resonance conditions.

According to jackson  the eigenmodes of the stationary cavity can be divided into TE and TM modes. For several cavities (rectangular, cylindrical, spherical) the eigenfrequencies are well-known. This enables one to determine the squeezing part of the interaction Hamiltonian.

In order to deduce the velocity Hamiltonian it will be necessary to find an appropriate model for the dispersive mirror. This can be achieved by using a thin dielectric slab with a high permittivity: $`\epsilon (x)=1+\gamma \delta (x)`$. As has been shown for a stationary system in scully  this leads to a similar eigenvalue equation as (28).

For the detection of the created field quanta some detecting device will have to be placed inside the cavity, e.g. an atom. However, the detector will always influence the created field as well. A simple approach for the modeling of a two-level system has been provided in dodonov3 ; dodonov5 . In addition, the non-adiabatic parametric modulation of the atomic Lamb shift – as has been considered in lambversuscasimir  – must be taken into account, since it will cause excitations of the atom as well.

Note that the induced quantum field will also excite the internal degrees of freedom of the cavity mirrors – an alternate description of losses should therefore also take the energy dissipation of the losses within the mirrors into account, see e.g. leakylaw .

Future work combining all these effects is of immense importance regarding experiments on quantum radiation using the dynamical Casimir effect.

## XI Acknowledgments

The authors are indebted to A. Calogeracos for fruitful discussions. R. S. acknowledges financial support by the Alexander von Humboldt foundation and NSERC. G. P. acknowledges financial support by BMBF, DFG, and GSI.

## Appendix A Multiple Scale Analysis

Starting with the Lagrangian (5) it is straightforward to show that the field operator fulfills a modified wave equation

$`\left\{\mathrm{}+2V(𝒓;t)\right\}\widehat{\phi }(𝒓,t)=0.`$ (142)
If one now follows dalvit  by introducing ladder operators via the expansion

$`\widehat{\phi }(𝒓,t)={\displaystyle \underset{n}{}}\widehat{a}_n^{\mathrm{in}}u_n(𝒓,t)+\mathrm{h}.\mathrm{c}.,`$ (143)
where

$`u_n(𝒓,t<0)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\mathrm{\Omega }_n^0}}}f_n(𝒓)e^{i\mathrm{\Omega }_n^0t},`$ (144)
$`u_n(𝒓,t>0)`$ $`=`$ $`{\displaystyle \underset{k}{}}Q_k^{(n)}(t)f_k(𝒓;t),`$ (145)
one can derive a time evolution equation for the coefficients $`Q_k^n(t)`$. Using the properties (7) of the eigenfunctions one obtains

$`\ddot{Q}_m^{(n)}(t)+\mathrm{\Omega }_m^2(t)Q_m^{(n)}(t)`$ $`=`$ $`2\lambda (t){\displaystyle \underset{k}{}}g_{mk}\dot{Q}_k^{(n)}(t)`$
$`+\dot{\lambda }(t){\displaystyle \underset{k}{}}g_{mk}\dot{Q}_k^{(n)}(t)`$
$`+\lambda ^2(t){\displaystyle \underset{k,l}{}}g_{lk}g_{lm}Q_k^{(n)}(t)`$
$`+\lambda ^2(t){\displaystyle \underset{k}{}}L_x{\displaystyle \frac{dg_{km}}{dL_x}}Q_k^{(n)}(t),`$
where $`\lambda (t)=\dot{L}_x(t)/L_x(t)`$ and $`L_x(t)=ca(t)`$ in our scenario. The antisymmetric coupling $`g_{mk}`$ is defined via

$`g_{mk}=L_x{\displaystyle _{\mathrm{cavity}}}{\displaystyle \frac{f_m}{L_x}}(𝒓)f_k(𝒓)d^3r`$ (147)
and is therefore related to the geometry factor $`m_{mk}`$ via $`g_{mk}=L_xm_{mk}`$. Note that compared to dalvit  the last line in (A) constitutes an additional term, since in our scenario the coupling between different modes may depend on the cavity parameters, see also equation (34). However, this difference is of minor relevance, since all these terms are accompanied by a factor of $`\lambda ^2(t)`$. If one assumes periodic oscillations of the cavity $`L_x(t)=L_x[1+ϵ\mathrm{sin}(\mathrm{\Omega }t)+ϵf(t)]`$, these terms can be neglected if the amplitude $`ϵ`$ is small. (The auxiliary function $`f(t)`$ is chosen to meet the continuity conditions on $`L_x(t)`$, see also dalvit .) Consequently, one can expand (A) in powers of $`ϵ1`$ to yield

$`\ddot{Q}_k^{(n)}(t)+(\mathrm{\Omega }_k^0)^2Q_k^{(n)}(t)`$ $`=`$ $`2\mathrm{\Omega }_k^0{\displaystyle \frac{\mathrm{\Omega }_k^0}{L_x}}L_xϵ\mathrm{sin}(\mathrm{\Omega }t)Q_k^{(n)}(t)`$ (148)
$`ϵ\mathrm{\Omega }^2\mathrm{sin}(\mathrm{\Omega }t){\displaystyle \underset{j}{}}g_{kj}Q_j^{(n)}(t)`$
$`+2ϵ\mathrm{\Omega }\mathrm{cos}(\mathrm{\Omega }t){\displaystyle \underset{j}{}}g_{kj}\dot{Q}_j^{(n)}(t)`$
$`+ϵ𝒪\left(f\right)+𝒪\left(ϵ^2\right).`$
This equation completely resembles the one found in dalvit . Note however, that we have to use the shifted eigenfrequencies and the eigenfunctions for leaky cavities. An approximate solution – for a more detailed discussion see dalvit  – can be obtained via introducing a new time scale $`\tau =ϵt`$ and inserting the formal expansion

$`Q_k^{(n)}(t)=Q_k^{(n)(0)}(t,\tau )+ϵQ_k^{(n)(1)}(t,\tau )+𝒪\left(ϵ^2\right)`$ (149)
with the unknown functions $`Q_k^{(n)(0/1)}`$ into equation (148). Finally, one has to sort in powers of $`ϵ`$. To lowest order one finds a free harmonic oscillator which can be solved by

$`Q_k^{(n)(0)}(t,\tau )=A_k^{(n)}(\tau )e^{i\mathrm{\Omega }_k^0t}+B_k^{(n)}(\tau )e^{i\mathrm{\Omega }_k^0t}.`$ (150)
The next order terms (proportional to $`ϵ`$) yield a driven harmonic oscillator equation for $`Q_k^{(n)(1)}`$ with the eigenfrequency $`\mathrm{\Omega }_k^0`$

$`_{t\tau }^2Q_k^{(n)(1)}+(\mathrm{\Omega }_k^0)^2Q_k^{(n)(1)}`$ $`=`$ $`2_{t\tau }^2Q_k^{(n)(0)}`$
$`2\mathrm{\Omega }_k^0L_x{\displaystyle \frac{\mathrm{\Omega }_k^0}{L_x}}\mathrm{sin}(\mathrm{\Omega }t)Q_k^{(n)(0)}`$
$`\mathrm{\Omega }^2\mathrm{sin}(\mathrm{\Omega }t){\displaystyle \underset{j}{}}g_{kj}Q_j^{(n)(0)}`$
$`+2\mathrm{\Omega }\mathrm{cos}(\mathrm{\Omega }t){\displaystyle \underset{j}{}}g_{kj}_tQ_j^{(n)(0)}.`$
In order to keep the expansion (149) convergent, above oscillator must not be at resonance. Consequently, all terms proportional to $`\mathrm{exp}(\pm i\mathrm{\Omega }_k^0t)`$ – with $`k`$ being the particular mode of interest – on the right hand side have to cancel. By imposing these conditions for the mode $`k=L`$ and for the coupling mode $`k=R`$ and inserting the frequency deviations

$`\mathrm{\Omega }`$ $`=`$ $`2\mathrm{\Omega }_L^0+h=2\mathrm{\Omega }_L^0+ϵ\alpha ,`$ (152)
$`\mathrm{\Omega }_R^0`$ $`=`$ $`3\mathrm{\Omega }_L^0+H=3\mathrm{\Omega }_L^0+ϵ\beta ,`$ (153)
one finds four linear and coupled evolution equations for the coefficients $`A_L^{(n)}(\tau )`$, $`B_L^{(n)}(\tau )`$, $`A_R^{(n)}(\tau )`$, and $`B_R^{(n)}(\tau )`$. These equations are – apart from the different couplings and the additional deviation $`H`$ – virtually identical with those presented in dalvit . Having applied the modified phase transformations

$`A_L^{(n)}(\tau )`$ $`=`$ $`e^{+i\alpha \tau /2}a_L^{(n)}(\tau ),`$ (154)
$`B_L^{(n)}(\tau )`$ $`=`$ $`e^{i\alpha \tau /2}b_L^{(n)}(\tau ),`$ (155)
$`A_R^{(n)}(\tau )`$ $`=`$ $`e^{+3i\alpha \tau /2}e^{i\beta \tau }a_R^{(n)}(\tau ),`$ (156)
$`B_R^{(n)}(\tau )`$ $`=`$ $`e^{3i\alpha \tau /2}e^{+i\beta \tau }b_R^{(n)}(\tau )`$ (157)
it is straightforward to rewrite these equations in matrix form.

$`{\displaystyle \frac{d}{d\tau }}\left(\begin{array}{c}a_k^{(n)}\\ b_k^{(n)}\\ a_j^{(n)}\\ b_j^{(n)}\end{array}\right)=\underset{¯}{M}\left(\begin{array}{c}a_k^{(n)}\\ b_k^{(n)}\\ a_j^{(n)}\\ b_j^{(n)}\end{array}\right),`$ (166)
where

$`\underset{¯}{M}=\left(\begin{array}{cccc}{\displaystyle \frac{i}{2}}\alpha & \gamma _1& \gamma _2& 0\\ \gamma _1& +{\displaystyle \frac{i}{2}}\alpha & 0& \gamma _2\\ \gamma _3& 0& {\displaystyle \frac{i}{2}}(3\alpha 2\beta )& 0\\ 0& \gamma _3& 0& +{\displaystyle \frac{i}{2}}(3\alpha 2\beta )\end{array}\right).`$
The quantities $`\gamma _1,\gamma _2,\gamma _3`$ have been introduced for convenience

$`\gamma _1`$ $`=`$ $`{\displaystyle \frac{L_x}{2}}{\displaystyle \frac{\mathrm{\Omega }_L^0}{L_x}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{(\mathrm{\Omega }_L^{x0})^2}{\mathrm{\Omega }_L^0}},`$ (168)
$`\gamma _2`$ $`=`$ $`{\displaystyle \frac{(2\mathrm{\Omega }_L^0+h)(2\mathrm{\Omega }_L^0+Hh/2)}{2\mathrm{\Omega }_L^0}}g_{LR}`$ (169)
$``$ $`2\mathrm{\Omega }_L^0g_{LR},`$
$`\gamma _3`$ $`=`$ $`{\displaystyle \frac{(2\mathrm{\Omega }_L^0+h)(2\mathrm{\Omega }_L^0+h/2)}{2(3\mathrm{\Omega }_L^0+H)}}g_{LR}{\displaystyle \frac{2}{3}}\mathrm{\Omega }_L^0g_{LR}`$ (170)
$``$ $`{\displaystyle \frac{\gamma _2}{3}}.`$