\section{Introduction}

In his seminal article 65 years ago, Casimir formulated a physical problem
\cite{Cas} 
which has had a tremendous influence on physics.
His pioneering analysis of the physical consequences of field
quantization under external, macroscopic boundary conditions 
is still of current relevance. Indeed, it turned out to be
one of the most prolific ideas in modern theoretical physics. 
Having a pure quantum background and being closely related to classical 
physics, the Casimir effect is a universal and wide-spread phenomenon
that can be found on all scales in the Universe.

One of the fascinating aspects of the Casimir effect is its simplicity, i.e.,
there are theoretical models that are amenable to exact solutions.
Leaving aside many important features of modern Casimir physics 
\cite{Dalvit_Book}, we consider the simple system of Ref.\cite{Cas}
with a slight generalization.
Namely, our setting consists of two half-spaces of
different materials with a plane-parallel gap in between. What
we know about the two media are their reflection coefficients (a matrix
in general) and macroscopic conditions like temperature distributions, 
macroscopic current densities, distance and relative motion parallel to 
the interfaces. We assume that these conditions are maintained 
stationary by external (generalized) forces.
We want to compute correlation
functions of the electromagnetic field in the gap. These yield for example
the average energy density of the field and the pressure exerted on the
two bodies as components of the stress (energy-momentum) tensor.
We focus in particular
on symmetrized correlations, also known as Keldysh-Green functions (KGF) 
in quantum kinetics, 
and are able to derive them under 
rather general circumstances in a stationary
Casimir geometry. We thus establish a natural non-equilibrium extension 
of Casimir's results.


Our approach starts from 
a simple observation which is a common feature among all
manifestations of the Casimir effect: the interaction between the two
interfaces disappears when their distance becomes large enough. 
On quite general grounds, we may therefore 
conclude that the electromagnetic field which we observe in
any Casimir system results from several fields that
originate in the material of the different half-spaces and in the vacuum
gap between them.
The distance dependence of the Casimir interaction is such 
that these fields become independent (uncorrelated) as the
interfaces are infinitely far away.
This fact is the cornerstone of our analysis: we use the large-distance limit
to express the KGF in a two-plate setting in terms of KGFs 
of independent half-spaces.

In this paper, we relax the assumption that the two plates share the same
temperature and state of motion. This makes it impossible to describe the
field between the plates as being in thermal equilibrium. 
The appropriate tool to calculate photon correlation functions is
then the Schwinger-Keldysh technique\cite{Schwinger_1961, Kel} of nonequilibrium
processes.
The application of the Keldysh formalism to the field of Casimir-Van der
Waals interactions has a relatively short history. The
first application, to our knowledge, is due to Janowicz, Reddig and
Holthaus \cite{Jan} in calculations of the electromagnetic heat
transfer between two bodies at different temperature. Sherkunov
used the non-equilibrium technique extensively for dispersion interactions
involving excited atoms and excited media \cite{Sher}. A general expression
for the electromagnetic force on an atom, in terms of the KGFs
of atom and radiation field, was found a few years ago
by one of the present authors \cite{vem}.

The physical processes 
behind these field-mediated interactions 
are the multiple reflections of photons between the
interfaces and their tunneling from one body to the other,
which follow from the boundary conditions for the electromagnetic field on the
interfaces.
To include these boundary conditions, we use an effective action
in the Schwinger-Keldysh technique with auxiliary fields and evaluate
generating functionals for KGFs by performing path integrals \cite{NO}.
Path integral approaches for the Casimir effect were introduced by Bordag,
Robaschik and Wieczorek \cite{Bord}. Li and Kardar considered the interaction 
between bodies mediated by a fluid with
long range correlations \cite{Meh2, Meh3}.
They applied a path integral technique to include
arbitrarily deformed bodies on which any kind of boundary condition can
be implemented. This feature makes the approach amenable to a perturbative
analysis of any deformed ideally conducting surfaces
\cite{Emig1, Emig2}. Further on, Emig
and B\"{u}scher \cite{Emig3} used the optical extinction theorem \cite{Born}
to reformulate the boundary conditions for a vacuum-dielectric interface
in  integral form that depends only on the fields on the vacuum side (i.e., in
the gap between two bodies). They then derived with path
integral techniques an effective Gaussian action for the photon gas in the
gap, providing the free energy and in particular the Casimir interaction
of dielectric bodies with arbitrary shaped surfaces. A similar approach has
been followed by Soltan et al.~\cite{Soltan} for a dispersive medium between the
bodies.
Recently, Behunin and Hu applied path integrals to problems of
the Casimir-Polder type (atom-surface interaction) in different
nonequilibrium 
situations: in Ref.\cite{Hu1} is considered the Van der Waals
interaction between an atom and a substrate in a stationary state out of
global equilibrium; Ref.\cite{Hu2} is calculating
atom-atom interactions in a quantized radiation field which is in a
nonequilibrium state.

Our theoretical approach to the non-equilibrium Casimir interaction 
where the boundary conditions involve different, locally defined
temperatures, is thus 
a synthesis between the Feynman path integral and the
Keldysh-Schwinger formalism of nonequilibrium processes \cite{KL, Alt}.

The paper is organized as follows: after preliminary notations and
introductions in Sec.\ref{s:2}, we introduce rather general expressions 
for retarded Green functions in the Casimir geometry of two planar
boundaries (Sec.\ref{s:3}). 
In Sec.\ref{s:4} is analyzed a Gaussian action in Schwinger-Keldysh space
which yields the non-equilibrium correlations for the fields (Keldysh-Green
functions). This action implements in particular the boundary conditions for
the retarded Green functions.
In Sec.\ref{s:5} we discuss
the retarded Green functions in the limit where the interfaces of the
Casimir system are infinitely removed from each other. 
Using this limiting
procedure, the Keldysh-Green function at any point in
the gap is expressed via Keldysh functions of the single interface systems
defined at the surfaces and via photon numbers in free space (Sec.\ref{s:7}). As an
application of the developed technique, we find in Sec.\ref{s:9} 
the field correlation functions in a Casimir geometry with a sliding interface.
This result is checked in
Appendix~\ref{a:sliding-reflection} where the same problem is solved using 
Rytov fluctuation electrodynamics.
We illustrate our results by analyzing the electromagnetic energy density 
between the plates in two typical non-equilibrium situations 
(Appendix~\ref{a:energy-spectrum}).
Concluding remarks are given in Sec.\ref{s:conclusions}, while technical
details are collected in the other Appendices.


\section{Preliminaries}
\label{s:2}

\subsection{Geometry of the problem}

We consider two bodies with 
parallel and homogeneous boundaries
located at $z = \pm a/2$. The boundaries are in stationary conditions: 
their temperatures are constant in time, and their relative motion 
(if any) is uniform and parallel to each other. 
We can then assume that the EM field in the cavity [$-a/2 \le z \le a/2$]
is stationary in time and homogenous in the $xy$-plane. As a consequence,
all relevant fields and correlation functions can be expanded in Fourier 
integrals with respect to frequency $\omega$ and wave vectors
${\bf q} = ( q_x, q_y )$ along the interfaces. 
We use in the following the shorthand 
\begin{equation}
	\Omega = (\omega , q_{x}, q_{y})
	\,.
	\label{eq:def-Omega}
\end{equation}
For fields like the vector potential, we get a mixed representation
\begin{equation}
	{\bf A}( \omega, q_x, q_y, z ) = {\bf A}( \Omega, z )
	\,,
	\label{eq:field-notation}
\end{equation}
where the argument $\Omega$ is suppressed where no confusion 
is possible. 

We work in the Dzyaloshinskii gauge $\varphi = 0$ where 
due to the transversality condition for the electric field, 
the normal component $A_{z}$ of the vector potential can be eliminated
in favor of the tangential ones $A_{x},\, A_{y}$. Furthermore, given the plane
of incidence spanned by ${\bf q}$ and the normal to the interfaces, 
the dynamical variables of the EM field in the cavity are the following
linear combination
\begin{equation}
A_{s} = \frac{ q_{y}A_{x} - q_{x}A_{y} }{ | {\bf q}| }
\,,\qquad
A_{p} =  \frac{ {\bf q} \cdot {\bf A} }{ | {\bf q}| }
\,,
\label{eq:1.2}
\end{equation}%
which are nothing but the vector potential amplitudes of s- and p-polarized 
waves. In the following, we call the components defined in Eq.(\ref{eq:1.2}) 
the Weyl representation of the vector potential or we
say the vector potential is written in the Weyl basis \cite{Pieplow_2012}:
$\hat A = (A_s, A_p)$.

\subsection{Weyl representation of free space Green function}
\label{s:Weyl-free-space}

The Green function (GF) of the EM field plays a crucial role in the following.
Since it simply represents the vector potential due to a point current source,
it is also represented by a mixed Fourier representation 
involving a tensor
$D_{\alpha\alpha'}( \Omega ; z, z' )$ ($\alpha, \alpha' = x, y, z$).
The latter is the solution of
(given the gauge $\varphi = 0$) 
\cite{Landau_9}%
\begin{equation}
\left[ (\partial_z^{2} + q_{z}^{2})\delta_{\alpha\beta} - 
{\partial}_{\alpha } {\partial}_{\beta } \right]
D_{\beta \alpha'}( \Omega ; z, z' ) 
= 4\pi \delta_{\alpha\alpha'} \delta(z - z')
\label{eq:1.3}
\end{equation}%
where (we set $c = 1$)
\begin{equation}
q_{z}^{2} = \omega^{2} - q^{2}%
\,,\quad
{\partial}_{\alpha } = (iq_{x}, iq_{y},\partial_z )
\label{eq:1.4}
\end{equation}%
In the following, we only need the tangential part of the GF 
$D_{\sigma \sigma'}$ with $\sigma, \sigma' = x, y$. 
(See Eqs.(\ref{eq:3.15.a}--\ref{eq:3.15.c}) below for the normal components.)
Projecting
both indices into the Weyl basis according to Eq.(\ref{eq:1.2}),
we get from Eq.(\ref{eq:1.3}) the simple Helmholtz equation
\begin{equation}
(\partial_z^{2} + q_{z}^{2}) \hat{D}( \Omega ; z, z' ) 
= 4\pi \hat{g} \delta (z - z'),  
\label{eq:2.5}
\end{equation}%
where the $2\times 2$ matrix $\hat{g}$ reads 
\begin{equation}
\hat{g}=\left( 
\begin{array}{cc}
1 & 0 \\ 
0 & q_{z}^{2}/\omega^{2}%
\end{array}%
\right) 
\,.  
\label{eq:1.7}
\end{equation}

The GF in free space (no boundary conditions) is written
$\hat{ \Delta }( \Omega ; z, z' )$. As is well known, it comes
in two types: 
\textit{retarded} $\hat{\Delta}^{R}$ and \textit{%
advanced} $\hat{\Delta}^{A}$ Green functions. In the Weyl 
representation, they are given by
\begin{equation}
\hat{\Delta}^{R}( \Omega ; z, z' ) 
= \hat{\Delta}_{0} \, {\rm e}^{i q_{z} | z-z' | }  
= \hat{\Delta}^{A\ast }( \Omega ; z, z' )
\label{eq:1.8}
\end{equation}%
where%
\begin{equation}
\hat{\Delta}_{0} = \frac{2\pi }{iq_{z}}\hat{g}  
\label{eq:1.9}
\end{equation}%
and the wave vector $q_{z}$ is defined over the entire frequency axis by
\begin{equation}
q_{z} = \left\{ 
\begin{array}{ll}
\mathop{\rm sgn}\omega \sqrt{\omega^{2} - q^{2}} + i0
&\mbox{for}\quad \omega^{2} > q^{2}%
\\
i\sqrt{q^{2} - \omega^{2}}
&\mbox{for}\quad q^{2} > \omega^{2}
\end{array}%
\right.  \label{eq:1.10}
\end{equation}%
The two cases corresponding to outgoing \textit{propagating}
and to \textit{evanescent} waves, respectively.
The infinitesimal imaginary part of $q_{z}$ in ({\ref{eq:1.10}) secures the
analytical continuation of $\hat{\Delta}^R$ to complex frequencies 
in the upper half plane. 
It entails the existence of the limit%
\begin{equation}
\lim_{a \to +\infty} {\rm e}^{iq_{z}a} = 0  
\label{eq:1.11}
\end{equation}
for both propagating and evanescent waves. In addition, we have
the symmetry relation 
\begin{equation}
	q_z( - \Omega ) = - q_z^*( \Omega )
	\,.
	\label{eq:2.10}
\end{equation}
As a consequence, the retarded GF has the property
$\hat{ \Delta }^{R}( - \Omega ; z, z' ) = 
\hat{ \Delta }^{R\ast}( \Omega ; z, z' )$, 
as it should, being a response function between real-valued fields.


\section{Retarded Green functions}
\label{s:3}

\subsection{Single interface}

We include stationary and translation-invariant (in the $xy$-plane) 
boundary conditions in the retarded GF by adding reflected waves.
Let us start with a single interface at $z = - a/2$ where the GF (subscript
$-$) reads 
\begin{equation}
\hat{D}_{-}^{R}( \Omega ; z, z' ) = \hat{\Delta}^{R}( \Omega ; z, z' ) + 
\hat{R}_{-}\hat{\Delta}_{0} \,{\rm e}^{iq_{z}(z + z' + a)},  
\label{eq:3.14.a}
\end{equation}%
for $z, z' \ge - a/2$. The first term is the same as in free space.
The reflection matrix $\hat{R}_-$ at the lower
interface takes a simple diagonal form in the frame where the lower body is 
at rest
\begin{equation}
	\hat{R}_- = \left( \begin{array}{cc} 
		R^s_- & 0 \\ 
		0 & R^p_- 
	\end{array}
	\right)
	\label{eq:diagonal-R-matrix}
\end{equation}
whose matrix elements are for a
\begin{equation}
	\mbox{metal}: \quad
	R^s_- = \frac{ q_z \zeta( \omega ) - \omega 
			}{ q_z \zeta( \omega ) + \omega }
	\,,\quad
	R^p_- = \frac{ \omega \zeta( \omega ) - q_z  
			}{ \omega \zeta( \omega ) + q_z }
	\,,
	\label{eq:R-for-metal}
\end{equation}
where $\zeta( \omega )$ is the (dimensionless) impedance,
and for a 
\begin{equation}
	\mbox{dielectric}: \quad
	R^s_- = \frac{ q_z  - q_{z\varepsilon}  
			}{ q_z  + q_{z\varepsilon} }	
	\,,\quad
	R^p_- = \frac{ q_{z\varepsilon} - \varepsilon( \omega ) q_z
			}{ q_{z\varepsilon} + \varepsilon( \omega ) q_z }	
	\,,
	\label{eq:R-for-dielectric}
\end{equation}
where $\varepsilon( \omega )$ is the dielectric permittivity and
\begin{equation}
q_{z\varepsilon} = [ \varepsilon (\omega )\omega^{2} - q^{2} ]^{1/2}
\label{eq:2.4}
\end{equation}%
the wave vector in the lower medium.

If only the upper medium is present, we have a GF similar to 
Eq.(\ref{eq:3.14.a}), for $z, z' \le + a / 2$
\begin{equation}
\hat{D}_{+}^{R}( \Omega ; z, z' ) = 
\hat{\Delta}^{R}( \Omega ; z, z' ) + 
\hat{R}_{+} \hat{\Delta}_{0} \,{\rm e}^{-iq_{z} ( z + z' - a ) }
\,, 
\label{eq:3.14.b}
\end{equation}%
where $\hat{R}_+$ is the corresponding reflection matrix. 
We suppose here a generic form for $\hat{R}_{+}$.
The case of a sliding upper interface is discussed in 
Appendix~\ref{a:sliding-reflection}.

\subsection{General properties}

Let us collect a few general properties of the reflection matrices and the
GFs. It follows from Eq.(\ref{eq:d.14.a}--\ref{eq:d.14.d}) for $\hat{R}_{+}$ and
from the diagonal form
of $\hat{R}_{-}$ [Eq.(\ref{eq:diagonal-R-matrix})] that both fulfill the identity
\begin{equation}
\hat{g}\hat{R}_{\nu }^{T}=\hat{R}_{\nu} \hat{g}
\,,\quad \nu = \pm
  \label{eq:2.6}
\end{equation}%
The retarded GF~(\ref{eq:3.14.a}) therefore satisfies
\begin{equation}
\hat{D}^{R}( \Omega ; z, z' ) = [ \hat{D}^{R}( \Omega ; z', z ) ]^{T} 
\label{eq:2.8}
\end{equation}%
which is \textit{not} the reciprocity condition because the
wave vector $\mathbf{q}$ in the arguments of both sides is the same.

Taking into account Eqs.(\ref{eq:2.10}) and~(\ref{eq:d.14.a}--\ref{eq:d.14.d}),
the reflection matrices in Eq.(\ref{eq:diagonal-R-matrix}) and
Sec.\ref{a:sliding-reflection} satisfy the identity
\begin{equation}
\hat{R}_{\nu} (- \Omega ) = \hat{R}_{\nu }^{\ast} ( \Omega )  
\label{eq:2.11}
\end{equation}%
where the asterisk denotes the element-wise complex conjugation.
This entails that the symmetry relation is also valid for the retarded GF at 
a single interface
\begin{equation}
\hat{D}^{R}_\nu\left( -\Omega ;z,z' \right) =\hat{D}^{R\ast }_\nu\left(
\Omega ;z,z' \right)
\,, \quad \nu = \pm
\label{eq:2.9}
\end{equation}%
Below, we shall also deal with the advanced GF which
is defined as
\begin{equation}
	 \hat{D}^A_\nu( \Omega ; z, z' ) =
	 [ \hat{D}^R_\nu( \Omega ; z', z )]^\dag 
	 =
	 \hat{D}^{R\ast}_\nu( \Omega ; z, z' )
	\label{eq:def-advanced-GF}
\end{equation}
where the last equality follows from Eq.(\ref{eq:2.8}).

\subsection{Planar cavity}

The preceding properties carry over to the retarded and advanced GF
in the cavity formed by two interfaces. By adding up multiply reflected waves,
one finds the expression
\begin{equation}
\hat{D}^{R}( \Omega ; z, z' ) = \hat{\Delta}^{R}( \Omega ; z, z') + 
\sum\limits_{\nu ,\nu' =\pm} \hat{C}^{ \nu\nu' } 
\, {\rm e}^{iq_{z} ( \nu z + \nu' z' ) }  
\label{eq:3.1}
\end{equation}%
where
\begin{equation}
\hat{C}^{--} = \hat{U}_{+-}^{-1} \hat{R}_{+} \hat{\Delta}_{0}
\,{\rm e}^{iq_{z}a}
\,, \quad
\hat{C}^{-+} = \hat{U}_{+-}^{-1} \hat{R}_{+} \hat{R}_{-} 
\hat{\Delta}_{0}
\,{\rm e}^{2i q_{z}a}  
\label{eq:2.2}
\end{equation}%
and $\hat{C}^{++},\hat{C}^{+-}$ are defined by swapping indices
$-\leftrightarrow +$ in Eq.(\ref{eq:2.2}).
The matrix $\hat{U}_{+-}$ takes into account multiple reflections of 
photons in the cavity; it is given by the expression 
\begin{equation}
\hat{U}_{+-}=\hat{I}-\hat{R}_{+}\hat{R}_{-}\,{\rm e}^{2 i q_{z}a}
\,,
\label{eq:3.5}
\end{equation}%
where $\hat{I}$ is the unit matrix. An analogous formula gives
$\hat{U}_{-+}$. Eq.(\ref{eq:2.6}) above gives us the property%
\begin{equation}
\hat{U}_{-+}^{-1}\hat{g} = \hat{g} \hat{U}_{+-}^{-1T}  
\label{eq:2.7}
\end{equation}%
which ensures that the generalized reciprocity relation~(\ref{eq:2.8}) also
holds for the cavity GF~(\ref{eq:3.1}). Similarly, the 
relation~(\ref{eq:def-advanced-GF}) between retarded and
advanced GFs remains true as well.

\subsection{Normal components}

To conclude this Section, we give the tensor components involving
the normal direction. The following formulas can be shown from the
wave equation Eq. (\ref{eq:1.3}) and the symmetry properties above
($\lambda = s,p$ and $q =  |{\bf q}|$): 
\begin{eqnarray}
q_{z}^{2} D_{z\lambda }^{R} &=& i q\partial_z D_{p\lambda }^{R}
\,,  
\label{eq:3.15.a}
\\
q_{z}^{2}D_{\lambda z}^{R} &=& -i q\partial_{z'} D_{\lambda p}^{R}
\,, 
\label{eq:3.15.b}
\\ q_{z}^{2}D_{zz}^{R} &=& 
\frac{ q^2 }{ q_{z}^{2} } \partial_z \partial_{z'} D_{pp}^{R}
+ 4\pi \delta ( z - z' )
\,.  
\label{eq:3.15.c}
\end{eqnarray}%

\subsection{Generalized impedance matrices}

The reflection matrices $\hat{R}_\nu$ appearing in the expressions above
are the solutions to a scattering problem at the planar interfaces. We
discuss here an equivalent formulation in terms of generalized surface
impedances. These will provide the link between boundary conditions
imposed on fields and interactions with auxiliary fields restricted 
to the interfaces.

Evaluating the derivative with respect to $z$ and $z'$ of the retarded
GF~(\ref{eq:3.1}) at the interfaces $\pm a/2$, we come to the 
boundary conditions. With respect to the first coordinate $z$
(derivative $\partial_z$), we find
\begin{eqnarray}
&&
iq_{z}^{-1} \partial_z \hat{D}^{R}( -a/2, z' ) 
- \hat{Y}_{-} \hat{D}^{R}( -a/2, z' ) = 0  
\quad
\label{eq:4.1.a}
\\
&&
iq_{z}^{-1}\partial_z\hat{D}^{R}( +a/2, z')
+ \hat{Y}_{+} \hat{D}^{R}( +a/2, z' )=0  
	\label{eq:4.1.b}
\end{eqnarray}
where we defined the  matrices
\begin{equation}
\hat{Y}_{\nu} = (\hat{I} + \hat{R}_{\nu} )^{-1} (\hat{I} - \hat{R}_{\nu} )
\,,\quad
\nu =\pm  
	\label{eq:4.3}
\end{equation}
These generalize the concept of a surface admittance to a general
reflection problem. 
Indeed, from the reflection amplitudes~(\ref{eq:R-for-metal}) for a metallic
surface at rest, we get in the Weyl basis
\begin{equation}
	\mbox{metal:}\quad
	\hat{Y}_- = \frac{ 1 }{ \zeta( \omega ) }
	\left( \begin{array}{cc} 
	\omega / q_z & 0 \\ 0 & q_z / \omega
	\end{array} \right)
	\label{eq:admittance-metal}
\end{equation}
Note that despite multiple reflections, the boundary 
conditions~(\ref{eq:4.1.a}, \ref{eq:4.1.b})
are of local character: they link the fields and their normal derivatives 
at the same position with the corresponding admittance matrices.

With respect to the second coordinate $z'$ (derivative $\partial_{z'}$)
of the GF, a similar calculation yields
\begin{eqnarray}
&&
iq_{z}^{-1}\partial_{z'}\hat{D}^{R}( z, {-a/2} )
- \hat{D}^{R}( z, -a/2 ) \hat{Y}_{-}^{T} = 0  
	\label{eq:4.2.a}
\quad
\\
&&
iq_{z}^{-1}\partial_{z'}\hat{D}^{R}( z, {+a/2} ) 
+ \hat{D}^{R}( z, +a/2 ) \hat{Y}_{+}^{T} = 0  
	\label{eq:4.2.b}
\end{eqnarray}%
where the admittance matrices appear transposed. 

In the case of a single interface, we still find two boundary conditions
at $z, z' = \pm a/2$. If the ``missing'' body is the upper one, for example,
the admittance degenerates into $\hat{Y}_+ = \hat{I}$ 
from Eq.(\ref{eq:4.3}).
The boundary condition~(\ref{eq:4.1.b}) then becomes equivalent to
the Sommerfeld condition for an outgoing wave: 
$\hat{D}^R( z, z' ) \sim {\rm e}^{ i q_z z }$. This holds as long as
$z > z'$ and in both the propagating and evanescent sectors.



\subsection{Remarkable identity}

Using the boundary conditions~(\ref{eq:4.1.a}--\ref{eq:4.2.b}),
their counterparts for the advanced GF $\hat{D}^A$,
and the Green equation~(\ref{eq:2.5}), we come to the
following property of the GFs in the two-plate geometry
\begin{eqnarray}
\sum\limits_{\nu =\pm} 
	&& \hat{D}^{R}(z, \nu a/2 ) \hat{\Gamma}^{\nu } \hat{D}^{A}(\nu a/2, z') 
	\nonumber
	\\[-1ex]
	&& = 
	\hat{D}^{R}( z, z') - \hat{D}^{A}( z, z' )  
	\label{eq:4.7}
\end{eqnarray}
In Eq.(\ref{eq:4.7}) we have introduced effective source strengths
\begin{equation}
\hat{\Gamma}^{\nu } = 
- \frac12 \left(\hat{\Delta}_{0}^{-1}\hat{Y}_{\nu} - \mbox{h.c.} \right)  
	\label{eq:4.8}
\end{equation}%
which are easily shown to be antihermitian:
\(
\hat{\Gamma}^{\nu \dagger } = -\hat{\Gamma}^{\nu }  
\)
and symmetric:
\(
\hat{\Gamma}^{\nu T} = \hat{\Gamma}^{\nu }.  
\)
Therefore, they have purely imaginary matrix elements:
\(
\hat{\Gamma}^{\nu \ast }=-\hat{\Gamma}^{\nu }  
.
\)
Besides, using the parity properties~(\ref{eq:2.10}) and~(\ref{eq:2.11}) 
of $q_{z}$ and the reflection
matrices, the definition of the $\hat{\Gamma}^{\nu }$ entails
\begin{equation}
\hat{\Gamma}^{\nu }(-\Omega ) = \hat{\Gamma}^{\nu \ast}(\Omega )
= -\hat{\Gamma}^{\nu}( \Omega )
	\label{eq:4.9.d}
\end{equation}
The identity~(\ref{eq:4.7}) has a long history in macroscopic
fluctuation electrodynamics. It has been noted by Eckhardt
\cite{Eckhardt_1982}, although involving a spatial integral over volumes
where the imaginary part of the permittivity is nonzero. The
version we give here is technically somewhat simpler because only
surface sources appear. This may be related to the ``holographic principle''
stating that under certain circumstances, all relevant 
properties of a (source-free) field are encoded in a hypersurface.
This is obviously related to the classical Huyghens principle. An alternative
proof of Eq.(\ref{eq:4.7}) is given in
Appendix~\ref{a:remarkable-identity} using the Leontovich surface
impedance boundary condition.


If one of the two interfaces is missing, an identity similar to
Eq.(\ref{eq:4.7}) can be derived analogously. 
One simply has to replace the source strength
for the missing interface by
\begin{equation}
	\Gamma^\nu \mapsto \Gamma^0 = 
	- \frac12 \left(\hat{\Delta}_{0}^{-1}\hat{I} - \mbox{h.c.} \right)
	=
	- \hat{\Delta}_{0}^{-1} \Theta( \omega^2 - q^2 )
	\label{eq:4.13}
\end{equation}
where $\Theta$ is the unit step function. Note that only propagating
waves appear on the free boundary: they represent the fields incident from
infinity towards the interface.


\section{Non-equilibrium action}
\label{s:4}

We now address the key problem of this paper: evaluate correlations
for the EM field under non-equilibrium conditions. To this effect, we use
the path integral method and work with an action for the EM field.
An auxiliary field $\varphi$ is introduced to enforce the
boundary conditions at the two interfaces~\cite{NO}. This technique 
was developed in previous work \cite{Bord, Meh2, Meh3} for equilibrium
situations. 
We extend the approach to the whole Schwinger-Keldysh space and define 
a Gaussian action for two coupled Bose fields: the vector potential
$A$ and the auxiliary field $\varphi$
\begin{eqnarray}
S &=& \frac{1}{2}\int \big\{
\check{A}^{\dag} ( \Omega , z ) 
\check{\Delta}^{-1}( \Omega ; z, z' ) 
\check{A}( \Omega , z') 
+ \check{\varphi}^{\dag}( \Omega ) 
\check{F}\check{\varphi}( \Omega ) 
\nonumber
\\
&& + \check{\varphi}^{\dag}( \Omega ) 
\check{M}( \Omega, z ) \check{A}( \Omega , z ) 
+ \check{A}^{\dag}( \Omega, z ) 
\check{M}^{\dag}( \Omega, z ) \check{\varphi}( \Omega ) 
\big\}
\,.  
\nonumber
\\
\label{eq:6.1}
\end{eqnarray}
The integration is over $\Omega$ and $z$. In the action~(\ref{eq:6.1}),
the vector potential $\check{A}$ lives in Keldysh space and has two
components
that are called \emph{quantum} $\hat{A}^{\rm q}$ and 
\emph{classical} $\hat{A}^{\rm cl}$:
\begin{equation}
\check{A} = \left( 
\begin{array}{c}
\hat{A}^{\rm q} \\ 
\hat{A}^{\rm cl}%
\end{array}%
\right) 
\,,  \label{eq:6.3.a}
\end{equation}
each of which contains the familiar transverse amplitudes in the Weyl basis
\begin{equation}
\hat{A} = \left( 
\begin{array}{c}
A_{s} \\ 
A_{p}%
\end{array}%
\right) 
\,.  
\label{eq:6.3.b}
\end{equation}
The $4\times 4$ matrix $\check{\Delta}^{-1}$ is the inverse of the
Keldysh-Green (KG) matrix of the \emph{free} EM field. In the Keldysh 
basis~(\ref{eq:6.3.a}) for $\check{A}$, this KG matrix has the block
structure
\begin{equation}
\check{\Delta} = \left( 
\begin{array}{cc}
\hat{0} & \hat{\Delta}^{A} \\ 
\hat{\Delta}^{R} & \hat{\Delta}^{K}%
\end{array}%
\right)
\,,   
\label{eq:6.2}
\end{equation}%
where $\hat{\Delta}^{R,A}$ are the retarded and advanced Green functions
for free space, introduced in Sec.\ref{s:Weyl-free-space}. 
The function $\hat{\Delta}^{K}$ is the KGF for the free field, it collects
symmetrized correlations of the vector potential.
(A calculation is sketched in Sec.\ref{s:8} below.) 
Our goal is to calculate its counterpart 
in the presence of 
the two interfaces that we denote $\hat{D}^{K}$. It is given in the 
coordinate representation~\cite{X} as ($x = ( t, {\bf r} )$) 
\begin{equation}
D_{\alpha \alpha'}^{K}( x, x' ) = -i \left\langle 
		\{\hat{A}_{\alpha}( x ), \hat{A}_{\alpha'}( x' ) \}_+
	\right\rangle 
\label{eq:8.3}
\end{equation}%
where $\{ \cdot, \, \cdot \}_+$ is the anti-commutator.
We work here with the corresponding Fourier transforms with respect
to $t - t'$ and to the tangential coordinates. In the Weyl basis, this 
gives the matrix $\hat{D}^{K}( \Omega, z, z' )$.

The auxiliary field $\check{\varphi}$ in Eq.(\ref{eq:6.1}) has
eight components that are also grouped in the Keldysh structure
\begin{equation}
\check{\varphi} = \left( 
\begin{array}{c}
\hat{\varphi}^{\rm q} \\ 
\hat{\varphi}^{\rm cl}%
\end{array}%
\right) \,.  \label{eq:6.4.a}
\end{equation}%
The components $\hat{\varphi}^{\rm q, \,cl}$ themselves are
\begin{equation}
\hat{\varphi}=\left( 
\begin{array}{c}
\varphi _{-,s} \\ 
\varphi _{-,p} \\ 
\varphi _{+,s} \\ 
\varphi _{+,p}%
\end{array}%
\right) 
\,. 
\label{eq:6.4.b}
\end{equation}%
They are Weyl spinors localized in the lower 
(index $-$) and the upper ($+$) interface. 
Finally, the matrix $\check{M}( \Omega, z )$ 
in the action~(\ref{eq:6.1}) is related, as we shall see below, 
to the boundary conditions imposed at $z = \pm a/2$.

\subsection{Evaluating the path integral}

We follow the standard path integral procedure and
add source terms to the action~(\ref{eq:6.1})%
\begin{equation}
\frac{1}{2}\int (\check{J}^{\dag} \check{A} +\check{A}^{\dag} \check{J})
\end{equation}%
Then Gaussian path integrals are evaluated, first over the EM field and then 
over the auxiliary field. We get the generating functional for field correlations
whose expansion to second order in $\check{J}$ provides the following
expression for the KG matrix:
\begin{equation}
\check{D} = \check{\Delta} 
+ \check{\Delta}\check{M}^{\dag} \check{\Lambda} \check{M} \check{\Delta}  \label{eq:6.5.a}
\end{equation}
where $\check{\Lambda}$ is the solution of
\begin{equation}
\check{\Lambda}^{-1} = \check{F} - \check{M}\check{\Delta}\check{M}^{\dag}  
\label{eq:6.5.b}
\end{equation}
The KG matrix $\check{D}$ of Eq.(\ref{eq:6.5.a}) must 
have the block structure~(\ref{eq:6.2}) with off-diagonal 
blocks that are Hermitian conjugates 
[Eq.(\ref{eq:def-advanced-GF})]
and with an antihermitian Keldysh block%
\begin{equation}
[ \hat{D}^{K}(\Omega ; z, z' )]^{\dag} 
= -\hat{D}^{K}(\Omega ; z', z)
\,.  \label{eq:6.6}
\end{equation}
In addition, we impose the boundary conditions~(\ref{eq:4.1.a}--\ref{eq:4.2.b})
on the off-diagonal elements $\hat{D}^{R},\hat{D}^{A}$ of~(\ref{eq:6.5.a}).
These conditions unambiguously
define the structure of the matrices $\check{M}$ and $\check{F}$ 
in the action~(\ref{eq:6.1}) as 
\begin{equation}
\check{F}=\left( 
\begin{array}{cc}
\hat{F} & \hat{0} \\ 
\hat{0} & \hat{0}%
\end{array}%
\right)  
\label{eq:6.8}
\end{equation}%
where $\hat{F}$ is a $4\times 4$ antihermitian matrix, and 
\begin{equation}
\check{M} = \left( 
\begin{array}{cc}
\hat{0} & \hat{M} \\ 
\hat{M}^{\ast} & \hat{0}%
\end{array}%
\right)  \label{eq:6.7}
\end{equation}%
Here the matrix $\hat{M}$ collects the boundary conditions at the
two interfaces into a vector of operators in the Weyl basis 
\begin{equation}
\hat{M}(\Omega, z) = \left( 
\begin{array}{c}
\delta( z + a/2 ) [ iq_{z}^{-1}\hat{I}\partial_z - \hat{Y}_{-} ] \\ 
\delta( z - a/2 ) [ iq_{z}^{-1}\hat{I}\partial_z + \hat{Y}_{+} ]%
\end{array}%
\right)   
	\label{eq:4.5}
\end{equation}%
The boundary conditions~(\ref{eq:4.1.a}--\ref{eq:4.1.b}) then take the
integral form
(common argument $\Omega$ suppressed)
\begin{eqnarray}
0 &=& \int\! dz \, \hat{M}( z ) \hat{D}^{R}( z, z' )
	\label{eq:4.4.a}
\\
0 &=& \int\! dz' \, \hat{D}^{R}( z, z' ) \hat{M}^{T}( z' )
	\label{eq:4.4.b}
\end{eqnarray}%
with a derivative $\partial_{z'}$ acting to the left in the second line.
For the advanced GF, the boundary conditions~(\ref{eq:4.2.a}--\ref{eq:4.2.b})
can be written as integrals over 
$\hat{M}^{\ast}( z ) \hat{D}^{A}( z, z' )$
and over $\hat{D}^{A}( z, z' ) \hat{M}^{\dag}( z' )$.

With the choices~(\ref{eq:6.8}, \ref{eq:6.7}), one can show that the
Keldysh action~(\ref{eq:6.1}) acquires that so-called ``\textit{causal structure}" \cite{KL, Alt}.

Inserting the matrices $\check{M}$ and $\check{F}$ from
Eqs.(\ref{eq:6.7}, \ref{eq:6.8}) into Eq.(\ref{eq:6.5.a}), we find for the 
KG functions $\hat{D}^{R,A,K}$ the expressions
\begin{eqnarray}
\hat{D}^{R} &=& \hat{\Delta}^{R} + 
\hat{\Delta}^{R}\hat{M}^{T}\hat{\Lambda}\hat{M} \hat{\Delta}^{R}  
	\label{eq:6.9.a}
\\
\hat{D}^{A} &=& \hat{\Delta}^{A} + \hat{\Delta}^{A} \hat{M}^{\dag} 
\hat{\Lambda}^{\dag} \hat{M}^{\ast} \hat{\Delta}^{A}  
	\label{eq:6.9.b}
\\
\hat{D}^{K} &=& \big( \hat{I} + \hat{\Delta}^{R}\hat{M}^{T}
\hat{\Lambda}\hat{M} \big)
\hat{\Delta}^{K} \big( \hat{I} + \hat{M}^{\dag} \hat{\Lambda}^{\dag} 
\hat{M}^{\ast} \hat{\Delta}^{A} \big)
\nonumber
\\
&& {} + \hat{\Delta}^{R}\hat{M}^{T}\hat{\Lambda}\hat{F}
\hat{\Lambda}^{\dag} \hat{M}^{\ast} \hat{\Delta}^{A}
\label{eq:6.9.c}
\end{eqnarray}
where $\hat{\Lambda}$ is a $4\times 4$ matrix given by
\begin{equation}
\hat{\Lambda} = - (\hat{M} \hat{\Delta}^{R} \hat{M}^{T} )^{-1}  
\label{eq:6.10}
\end{equation}
The expressions~(\ref{eq:6.9.a}) and~(\ref{eq:6.9.b}) are integral forms of 
retarded and advanced GF because matrix products actually have to be read
as the concatenation of integral operators. It is trivial to check that they
satisfy the boundary conditions (\ref{eq:4.4.a}, \ref{eq:4.4.b}). 
Using the explicit form (\ref{eq:4.5}) of $\hat{M}$, we have checked in
a straightforward calculation
that Eq.(\ref{eq:6.9.a}) coincides with Eq.(\ref{eq:3.1}) for the retarded GF.

\subsection{Distribution (Keldysh-Green) functions for photons}

To analyse the expression~(\ref{eq:6.9.c}) for the KG function,
we first observe that
the first term is equal to zero because $\hat{\Delta}^{K}$ is a
solution of the homogeneous equation corresponding to Eq.(\ref{eq:2.5}). 
(See Appendix~\ref{a:homogeneous-DeltaK} for details.)
We shall argue that $\hat{F}$ that characterizes the auxiliary fields
can be taken in block-diagonal form
\begin{equation}
\hat{F}=\left( 
\begin{array}{cc}
4\hat{\Delta}_{0}\hat{P}(-)\hat{\Delta}_{0}^{\dag}  & \hat{0} 
\\ 
\hat{0} & 4\hat{\Delta}_{0}\hat{P}(+)\hat{\Delta}_{0}^{\dag} %
\end{array}%
\right)   
\label{eq:6.11}
\end{equation}%
where the quantity $\hat{P}( - )$ [$\hat{P}( + )$] correspond to the
lower [upper] boundary, respectively.
Adopting this choice, a tedious, but elementary calculation leads us to 
(common argument $\Omega$ suppressed)
\begin{equation}
\hat{D}^{K}( z, z' ) = 
\sum\limits_{\nu = \pm} 
\hat{D}^{R}( z, \nu a/2 ) \hat{P}( \nu ) \hat{D}^{A}( \nu a/2, z' )
	\,, 
\label{eq:6.12}
\end{equation}%
for the KG function. This has the same form as the Keldysh equation 
in quantum kinetics \cite{Kel}. If we would take this formal analogy serious,
then the quantities 
$\hat{P}(\Omega ;\nu )$ 
would coincide with the Keldysh polarization operators on the boundaries 
(i.e., the KG correlation of interface polarization fields). 
In the non-equilibrium theory, 
these are nonlinear functionals of the Keldysh function for
photons, so that in the semiclassical approximation, Eq.(\ref{eq:6.12}) 
would generate
a kinetic equation of Boltzmann type for the photon distribution function. 

In the two-plate geometry of the Casimir effect, however, this is not
the case because the
$\hat{P}( \nu )$ are \emph{independent} of the photon KG 
functions $\hat{D}^{R}$, $\hat{D}^{A}$, and $\hat{D}^{K}$. 
They are rather determined by given 
macroscopic states of the bodies and their interfaces and act like a 
linear driving for the EM field in the cavity. This suggests an interpretation 
of Eq.(\ref{eq:6.12}) in the spirit of Rytov's fluctuation electrodynamics:
the quantities $\hat{P}( \pm )$ encapsulate the distribution of
photon sources. They are the only piece of information needed to
produce the (non-equilibrium) distribution function
of photons in the cavity. For this reason, we call $\hat{P}( \pm )$
the \emph{photon sources} or \emph{fluctuating sources} in the following.

The above interpretation also helps to understand the 
diagonal form for the source distribution function
$\hat{F}$ in Eq.(\ref{eq:6.11}): the sources are 
located on the bodies' interfaces,
and correlations between the macroscopic states of 
the two bodies are neglected. We intend to investigate corrections to
this approximation in future work.




\section{The large distance limit}
\label{s:5}


In this Section, we consider the limit $a \to \infty$ in order to fix the strength
of photon sources by referring to the EM field in free space and near a single interface.


There are three different ways to take the limit:

({\bf A}) If we fix both points $z$, $z'$ in the cavity and go to the limit
$a \to \infty$, we recover Green functions in free space. In particular for
the retarded Green function,
\begin{equation}
      \lim\limits_{a \to \infty} \hat{D}^{R}( z, z' ) = \hat{\Delta}^{R}( z, z' )
      \label{eq:5.1a}
\end{equation}
with the free space Green function defined in Eq.(\ref{eq:1.8}).

({\bf B}) If the points stay at fixed distances $z, z' > 0$ from the lower interface,
we get 
\begin{equation}
      \lim\limits_{a \to \infty} \hat{D}^{R}( z - a/2, z' - a/2 ) 
      = \hat{D}^{R}_{-}( z, z'; a \to 0 ) 
      \equiv \hat{d}^{R}_{-}( z, z' )
      \label{eq:5.1b}
\end{equation}
where $\hat{D}^{R}_{-}$ is the Green function above a single interface. 
The notation $( \ldots ; a \to 0 )$ means that in the expression for
$\hat{D}^{R}_{-}$ [Eq.(\ref{eq:3.14.a})], 
$a$ should be set to zero.

({\bf C}) Similarly, we may take points at positions $z, z' < 0$ below the upper
interface and get the Green function below a single interface:
\begin{equation}
      \lim\limits_{a \to \infty} \hat{D}^{R}( z + a/2, z' + a/2 ) 
      = \hat{D}^{R}_{+}( z, z'; a \to 0 ) 
      \equiv \hat{d}^{R}_{+}( z, z' )
      \label{eq:5.1c}
\end{equation}
where the Green function $\hat{D}^{R}_{+}( z, z'; a \to 0 )$ is defined 
similar to case ({\bf B}) from Eq.(\ref{eq:3.14.b}).

We also need for the Keldysh-Green function Eq.(\ref{eq:6.12})
the limiting values when one position in the Green functions
recedes to infinity. Keeping $z$ fixed,
\begin{eqnarray}
	{a \to \infty}:\quad  && \hat{D}^{R}( z, \nu a/2 )
	\to
	\hat{\Delta}^{R}( z, \nu ) \, {e}^{ i q_z a/2 }
	\,,
\\
	&& 
	\hat{\Delta}^{R}( z, \nu )
	=
	\hat{\Delta}_0 \, {e}^{ - i \nu q_{z} z }
	\label{eq:5.2.a}
\end{eqnarray}
In this limit, the exponential ${\rm e}^{ i q_z a/2 }$ in the first line
restricts the support 
to propagating waves (real $q_z$);  it drops out when products with
the advanced GF $\hat{D}^{A}$ are formed because of Eq.(\ref{eq:1.8}). 
In a similar way, we
define for the single-interface limits ({\bf B}, {\bf C}) the
functions
\begin{eqnarray}
	\hat{d}_{-}^{R}( z, \nu ) &=& 
	[ \hat{I} e^{ - i \nu q_{z} z } 
		+ \hat{R}_{-} e^{ i q_{z}z} ]
		\hat{\Delta}_{0}  
\label{eq:5.2.b}
\\
	\hat{d}_{+}^{R}( z, \nu ) &=&
	[ \hat{I} e^{ - i \nu q_{z} z}
		+ \hat{R}_{+} e^{ -i q_{z}z} ]
		\hat{\Delta}_{0}
\label{eq:5.2.c}
\end{eqnarray}


\section{Keldysh--Green functions}
\label{s:7}

We now consider the KG function derived at Eq.(\ref{eq:6.12}) in the 
limit $a\rightarrow +\infty$ in order to find the photon sources 
$\hat{P}( \nu )$.
Using the notation $\hat{\Delta}^{K}$, $\hat{d}_{\pm }^{K}$ 
and $\hat{P}_{0,-,+}$ in the three limits (\textbf{A}), ({\bf B}) 
and (\textbf{C}), we get
\begin{eqnarray}
\hat{\Delta}^{K}( z, z' ) &=& \sum\limits_{\nu} 
	\hat{\Delta}^{R}(z,\nu )
	\hat{P}_{0}(\nu )
	\hat{\Delta}^{A}(\nu ,z')
	\,,  
\label{eq:7.1.a}
\\
\hat{d}_{\mp}^{K}( z, z' ) &=& \sum\limits_{\nu }
	\hat{d}_{\mp}^{R}(z,\nu )
	\hat{P}_{\mp}(\nu )
	\hat{d}_{\mp}^{A}(\nu ,z')
	\,,  
\label{eq:7.1.b}
\end{eqnarray}%
where the GFs on the right hand sides are defined in 
Eqs.(\ref{eq:5.2.a}--\ref{eq:5.2.c}) above.
Manifestly, the sources in free space $\hat{P}_{0}(\nu )$ are defined by
the quantum state of \emph{free} EM field. Similarly, the fluctuation sources 
on the physical interfaces [the lower one, $\hat{P}_{-}(-)$ in case 
({\bf B}) and the upper one, $\hat{P}_{+}(+)$ in case ({\bf C})] 
are defined by the macroscopic quantum state of the corresponding bodies.

We assume the statistical independence of the bodies and the free 
EM field \cite{Henkel_2002}
and come to the conclusion that the sources on the `free interfaces'
coincide:
\begin{equation}
\hat{P}_{-}(+) = \hat{P}_{0}(+)
\, ; \quad
\hat{P}_{+}(-) = \hat{P}_{0}(-)
\,,
\label{eq:7.2.a}
\end{equation}%
and that the sources on the macroscopic interfaces are the same for one
and for two plates:
\begin{equation}
\hat{P}(-) = \hat{P}_{-}(-)
\,; \quad
\hat{P}(+) = \hat{P}_{+}(+)
\,.  
\label{eq:7.2.b}
\end{equation}
By evaluating Eq.(\ref{eq:7.1.b}) at $z = z' = 0$,
we can also express the interface sources in terms of the KG function 
there:
\begin{equation}
\hat{P}(\nu ) = \hat{\Delta}_{0}^{-1}
	(\hat{I} + \hat{R}_{\nu })^{-1} \hat{d}_{\nu}^{K}
	(\hat{I} + \hat{R}_{\nu }^{\dagger})^{-1}
	\hat{\Delta}_{0}^{-1\dagger}
	- \hat{P}_{0}( -\nu )  
\label{eq:7.3}
\end{equation}
where $\hat{d}_{\nu }^{K}$ is the boundary value of the KG function
for a system with a single interface
\begin{equation}
\hat{d}_{\nu }^{K} = \hat{d}_{\nu }^{K}(\Omega ; 0, 0)  
\label{eq:7.4}
\end{equation}
The free space photon sources $\hat{P}_0( \nu )$
are calculated in Appendix~\ref{s:8}.
In Sec.\ref{s:9}, we work out the quantities
$\hat{P}( \nu )$ for a specific example, allowing for the upper body
to be in uniform motion relative to the lower one.

We are now ready to collect our main result. The explicit expressions for 
the GFs $\hat{D}^{R, A}$ [Eqs.(\ref{eq:3.14.a},
\ref{eq:def-advanced-GF})] and for the photon sources 
$\hat{P}$ [Eq.(\ref{eq:7.2.b})]
are inserted into the KG function~(\ref{eq:6.12}) to give
\begin{equation}
\hat{D}^{K}( z, z' ) = 
	\sum\limits_{\nu, \nu' = \pm}
	\hat{D}^{\nu\nu'}
	e^{ i (\nu q_{z}z -\nu' q_{z}^{\ast} z') }
	\,. 
\label{eq:6.5}
\end{equation}%
We find the amplitudes
\begin{eqnarray}
\hat{D}^{--} &=& \hat{T}_{+} + e^{-2 {\rm Im}\, q_{z} a}
	\hat{R}_{+}\hat{T}_{-}\hat{R}_{+}^{\dagger }
	\,,
\label{eq:7.6.a}
\\
\hat{D}^{-+} &=& \hat{R}_{+} \hat{T}_{-} e^{iq_{z}a} + 
	\hat{T}_{+}\hat{R}_{-}^{\dagger } e^{-iq_{z}^{\ast }a}
	\,,
\label{eq:7.6.b}
\\
\hat{T}_{-} &=& \hat{U}_{-+}^{-1} \hat{\gamma}_{-}
	\hat{U}_{-+}^{-1\dagger}
	\,;
\label{eq:7.7.a}
\end{eqnarray}%
where $\hat{D}^{++}$, $\hat{D}^{+-}$, $T_{+}$ are found by swapping
the subscripts $+$ and $-$, $\hat{U}_{-+}$ has been defined in
Eq.(\ref{eq:3.5}), 
and 
\begin{equation}
\hat{\gamma}_{\nu } = e^{-{\rm Im}\, q_{z} a}
	(\hat{I}+\hat{R}_{\nu })\hat{\Delta}_{0}
	\hat{P}(\nu )\hat{\Delta}_{0}^{\dagger }
	(\hat{I}+\hat{R}_{\nu }^{\dagger})  
\label{eq:7.7.b}
\end{equation}%
These expressions give the distribution function of photons in a
planar cavity (homogeneous along the $xy$-directions and
stationary in time) whatever the macroscopic state of the two
boundaries, as encoded in $\hat{P}( \nu )$. Based on the assumption of
statistical independence, the $\hat{P}( \nu )$ are given by
reference situations with a single interface [Eqs.(\ref{eq:7.3}, \ref{eq:7.4})].

In some cases (for instance in the Casimir effect), one needs to subtract
the free-space KG function to get finite results for the relevant observables:
\begin{equation}
\hat{D}_{\rm ren}^{K} = \hat{D}^{K} - \hat{\Delta}^{K}  
\label{eq:7.8}
\end{equation}%
where $\hat{\Delta}^{K}$ is defined by Eq.(\ref{eq:7.1.a}) with 
free-space photon sources $\hat{P}_{0}( \nu )$ in the problem.

For completeness, the normal ($z$-) components of the KG functions are also
given here. They are expressed in terms of tangential components, using 
the homogenous version of Eq.(\ref{eq:1.3}). We get in analogy to 
Eqs.(\ref{eq:3.15.a}--\ref{eq:3.15.c}) for the retarded GF
\begin{eqnarray}
q_{z}^{2}D_{z\lambda }^{K} &=& i q\partial_z D_{p\lambda }^{K}
\,,  
\label{eq:7.9.a}
\\
q_{z}^{2}D_{\lambda z}^{K} &=& -i q\partial_{z'} D_{\lambda p}^{K}
\,, 
\label{eq:7.9.b}
\\
q_{z}^{2}D_{zz}^{K} &=& \frac{ q^{2} }{ q_{z}^{2}} 
	\partial_z \partial_{z'} D_{pp}^{K}
\,.
\label{eq:7.9.c}
\end{eqnarray}

\section{Example: sliding interfaces}
\label{s:9}

As an application of the theory developed so far, 
let us consider the KG function in the case of two bodies in relative
motion.
In this case in the limit ($\mathbf{C}$) we have
a sliding interface (moving parallel to $x$ with velocity $v$). 
We suppose that we have equilibrium in the body's rest frame (with 
temperature $T_{+}'$),
similar to Refs.\cite{Polevoi_1990, Kyasov_2002, Dedkov_2002b}.
Using the Lorentz covariant formulation
of the fluctuation-dissipation theorem \cite{Pieplow_2013}, we get for the 
KG function in the limit ($\mathbf{C}$) the expression
($k_B = \hbar = 1$)
\begin{eqnarray}
\hat{D}_{+}^{K}( \Omega, z, z' ) &=& [\hat{d}_{+}^{R}( \Omega, z, z' )
	- \hat{d}_{+}^{A}( \Omega, z, z' )]
	\coth \frac{ \omega' }{ 2T_{+}' }
\,;
\nonumber
\\
\\
\omega' &=& \gamma (\omega - v q_{x})  
\label{eq:9.1}
\end{eqnarray}%
where $\omega'$ is the Doppler-shifted frequency.
A calculation starting from Eq.(\ref{eq:7.3}) yields the corresponding 
source located at the upper surface
\begin{equation}
\hat{P}(+) = \hat{\Gamma}^{+} \coth \frac{ \omega' }{ 2T_{+}' }
\label{eq:P-plus}
\end{equation}%
where $\hat{\Gamma}^{+}$ is given in Eq.(\ref{eq:4.8}). The
reflection matrix $\hat{R}_+$ that appears in
$\hat{d}^{R}_{+}$ [see Eq.(\ref{eq:5.2.c})] 
is calculated in 
Appendix~\ref{a:sliding-reflection}, Eqs.(\ref{eq:d.13}--\ref{eq:d.15}).
These expressions cover any relative velocities $v$.

The lower interface is in equilibrium at $T_{-}$. Therefore%
\begin{equation}
\hat{D}_{-}^{K}( \Omega, z, z' ) = [ \hat{d}_{-}^{R}( \Omega, z, z' )
	- \hat{d}_{-}^{A}( \Omega, z, z' ) ]
	\coth \frac{ \omega }{ 2T_{-} }
\label{eq:9.3}
\end{equation}%
and then using Eq.(\ref{eq:7.3}) we get for the sources there
\begin{equation}
\hat{P}(-) = \hat{\Gamma}^{-} \coth \frac{ \omega }{ 2T_{-} }
\label{eq:c.4}
\end{equation}
Inserting 
Eqs.(\ref{eq:P-plus}, \ref{eq:c.4}) 
into Eq.(\ref{eq:6.12}),
we come to the KG function for the cavity with one sliding interface
\begin{eqnarray}
&& \hat{D}^{K}( \Omega; z,z' )
\label{eq:9.5}
\\
&& =
	\hat{D}^{R}( \Omega; z,-a/2)
	\hat{\Gamma}^{-}
	\hat{D}^{A}( \Omega; -a/2, z' )
	\coth \frac{ \omega }{ 2T_{-} }
\nonumber
\\
&& \quad {} + \hat{D}^{R}( \Omega; z, +a/2)
	\hat{\Gamma}^{+}
	\hat{D}^{A}( \Omega; +a/2, z' )
	\coth \frac{ \omega' }{ 2T_{+}' }
\nonumber
\end{eqnarray}
That the frequencies differ in the two $\coth$ terms has been known
in similar contexts since the pioneering work by Frank and Ginzburg 
on the Cherenkov effect (see 
Ref.\cite{Ginzburg_1996} for an overview).
This spoils any attempt to describe the sliding geometry by a global 
equilibrium assumption even if the two bodies are of the same temperature
(as in Ref.\cite{Philbin_2009}).
The sign change between $\omega$ and $\omega'$ (anomalous Doppler
effect) has also been noted in early work on field quantization in
moving media, see, e.g., Jauch and Watson \cite{Jauch_1948}.

\section{Conclusions}
\label{s:conclusions}

We have calculated the photon distribution function between two parallel plates
(Casimir geometry) under rather general conditions: 
on the two plates, any arbitrary stationary nonequilibrium state
is allowed for. Using 
the Schwinger-Keldysh formalism of nonequilibrium field theory, 
we could express the Keldysh-Green
function of photons at any point in the gap between the plates via the Keldysh 
functions of single interface systems,
defined at the surfaces, and via photon numbers in free space. As a cross check of 
the results, we consider one plate sliding relative to the other in local equilibrium,
and we find full coincidence of the results with Rytov theory,
without any restriction in relative velocities.

Our approach is flexible enough to allow for the zero-temperature
limit $T \to 0$ to be taken. The example of the sliding plates of Sec.\ref{s:9} 
then illustrates that one does \emph{not} recover an equilibrium situation.
The resulting frictional stress that opposes the relative motion is an interesting
(and controversial \cite{Philbin_2009, Dedkov_2010a, Pendry_2010, Barton_2011}) 
manifestation of 
an unstable vacuum state.
For similar situations, we may quote the Klein paradox \cite{IV}
(electron-positron pairs created by maintaining a 
static field
configuration), and
the Schwinger-Unruh effect (thermalization of an accelerated detector in
vacuum). Indeed, one possible explanation for quantum friction involves the
creation of particle pairs in the two plates \cite{Pendry_2010, Barton_2011, Maghrebi_2013b},
as pointed out in earlier work by Polevoi within the context of Rytov theory
\cite{Polevoi_1990}.

Finally, we suggest that the developed formalism is general enough to investigate
generalizations beyond the standard fluctuation electrodynamics. The crucial
assumptions of the latter are clearly spelled out: statistical independence
of the sources localized on the macroscopic bodies. If this is relaxed, one has
to establish the photon source strengths in some other way. Concepts from
non-equilibrium kinetic theory like the Boltzmann equation or the balance
of energy and entropy exchanges are likely to be instrumental here.

