Canonical quantization of electromagnetism in spatially dispersive media

In macroscopic electromagnetism the dependence of the dielectric functions on frequency, and the relation of this dependence to absorption, have a central place. The circumstances where dispersion and absorption can be neglected are a small subset of cases where macroscopic electromagnetism is successfully applied. In contrast, the dependence of the dielectric functions on wavevector (spatial dispersion, or nonlocal response) is normally regarded as of rather limited relevance (see e.g. [1]). Nevertheless, there are some situations where the inclusion of spatial dispersion is required in order to obtain accurate results. Recent work in the nano-optics of small metallic particles [2–4] has found that a nonlocal description can have a significant—and sometimes counterintuitive—effect on the localization of radiation. For example, in the system of touching nanowires investigated in [4] the inclusion of spatial dispersion turns a continuous spectrum into a discrete one. In this work we are interested in the consequences of a nonlocal description for the quantum mechanical properties of the electromagnetic field. We shall demonstrate that these are also nontrivial, and that a nonlocal description is sometimes necessary in order to obtain finite results.

Here we aim both to set up a canonical formulation of quantum electromagnetism in spatially and temporally dispersive media (in this respect we look to extend [5] and add to the work of Scheel and co-workers [6], and Suttorp [7]), and to show that such an extension to the theory can naturally cure some divergences that plague the local theory. The simplest example of such a divergence is evident in the zero temperature electric field correlation function in a homogeneous medium, written in terms of the electromagnetic Green function, G,

$$\begin{equation} \fl\langle0|\hat{\boldsymbol{E}}(\boldsymbol{r},t){\otimes}\hat{\boldsymbol{E}}(\boldsymbol{r}^{\prime},t)|0\rangle=\int_{0}^{\infty}\mathrm{d}\omega\boldsymbol{c}\,(\boldsymbol{r},\boldsymbol{r}^{\prime},\omega)=\frac{\hbar\mu_{0}}{\pi}\int_{0}^{\infty}\mathrm{d}\omega\,\omega^{2}{\rm Im}\,\left[\boldsymbol{G}(\boldsymbol{r},\boldsymbol{r}^{\prime},\omega)\right], \end{equation} \tag{ 1 }$$

It is well known that equation (1) diverges as r → r', even in free space due to the integration over the zero-point fields at all frequencies. However the integrand, c(r,r',ω), should certainly not diverge in this limit for it represents the intensity of the electric field at a fixed frequency, which determines single-frequency phenomena such as spontaneous emission rates. In vacuum c (r,r',ω) is finite, but inside a dielectric it diverges. Given that the electromagnetic Green function satisfies

$$\begin{eqnarray*} &&\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{G}(\boldsymbol{r},\boldsymbol{r}^{\prime},\omega)-\frac{\omega^{2}\varepsilon(\omega)}{c^{2}}\boldsymbol{G}(\boldsymbol{r},\boldsymbol{r}^{\prime},\omega)={\mathbbm{1}}_{3}\delta^{(3)}(\boldsymbol{r}-\boldsymbol{r}^{\prime}) , \end{eqnarray*}$$

it has a longitudinal part

$$\begin{equation} \fl\boldsymbol{G}_{\parallel}(\boldsymbol{r},\boldsymbol{r}^{\prime},\omega)=-\frac{c^{2}}{\omega^{2}\varepsilon}\boldsymbol{\delta}_{\parallel}(\boldsymbol{r}-\boldsymbol{r}^{\prime})\to{\rm Im\,}[\boldsymbol{G}_{\parallel}(\boldsymbol{r},\boldsymbol{r}^{\prime},\omega)]=\frac{c^2\,{\rm Im\,}[\varepsilon]}{\omega^{2}|\varepsilon|^{2}}\boldsymbol{\delta}_{\parallel}(\boldsymbol{r}-\boldsymbol{r}^{\prime}), \end{equation} \tag{ 2 }$$

which contributes an infinite amount to c (r,r,ω) when Im [ε(ω)] ≠ 0. This divergence stems from an assumption that the current flowing within the dielectric can have an arbitrarily rapid variation in space, as can be seen from the equal time correlation function [5],

$$\begin{equation} \langle0|\hat{\boldsymbol{\jmath}}(\boldsymbol{r},t){\otimes}\hat{\boldsymbol{\jmath}}(\boldsymbol{r}^{\prime},t)|0\rangle=\hbar\pi{\mathbbm{1}}_{3}\delta^{(3)}(\boldsymbol{r}-\boldsymbol{r}^{\prime})\int_{0}^{\infty}\mathrm{d}\omega\,\omega^{2}{{\rm \,Im\,}}[\varepsilon(\omega)], \end{equation} \tag{ 3 }$$

where the delta function on the right hand side of (3) shows that this current is made up from arbitrarily large component wavevectors. The divergence would therefore be removed if we assumed that electromagnetic energy is always absorbed over some nonzero volume of space, replacing $\boldsymbol {\mathbbm{1}}_{3}\delta ^{(3)}(\boldsymbol {r}-\boldsymbol {r}^{\prime }){{\mathrm { Im\,}}}[\varepsilon (\omega )]$ in (3) with Im [ε(r,r',ω)]. Once this description is adopted then it is possible to obtain finite results for c(r,r,ω).

The link between such divergences encountered in macroscopic quantum electrodynamics (QED) and the failure to account for spatial dispersion has been previously alluded to (e.g. [8, 9]), but to the authors' knowledge there has been no attempt to incorporate spatial dispersion into macroscopic QED from the outset to eliminate these divergences. Perhaps part of the reason for this is that one requires a knowledge of the entire wavevector dependence of the permittivity of material samples, and such a characterization is currently unavailable. Both the theoretical and experimental understanding of spatial dispersion are quite underdeveloped compared to other aspects of macroscopic electromagnetism, and there is much scope in future years for important developments in this topic.

The structure of this paper is as follows: the standard description of spatial dispersion in a homogeneous medium is outlined in section 2. In section 3 we give a Lagrangian whose field equations are the macroscopic Maxwell equations in a homogeneous medium with nonlocal response. This allows the classical theory to be canonically quantized (section 4) and the Hamiltonian is diagonalized in section 5. We apply the resulting quantum description to thermal and zero-point fields in a homogeneous medium and show that the resulting electromagnetic field intensities are no longer infinite when we assume a particular model of spatial dispersion (section 6). In section 7 we describe a range of basic questions concerning the Casimir effect (zero-point and thermal electromagnetic fields) that cannot be tackled, even approximately, if spatial dispersion is not included.

In this section we briefly recall the form of the classical macroscopic Maxwell equations when spatial dispersion is included. A detailed account can be found in [1, 9, 10]. The electromagnetic response of all real materials is nonlocal in space as well as in time, and the D and H fields must take the form

$$\begin{equation} \boldsymbol{D}(\boldsymbol{r},t)=\varepsilon_0\int_{-\infty}^t{\rm d} t' \int{\rm d}^3\boldsymbol{r}' \, \boldsymbol{\varepsilon}(\boldsymbol{r},\boldsymbol{r}',t-t')\,{\boldsymbol{\cdot}}\,\boldsymbol{E}(\boldsymbol{r}',t'), \end{equation} \tag{ 4 }$$

$$\begin{equation} \boldsymbol{B}(\boldsymbol{r},t)=\mu_0\int_{-\infty}^t{\rm d} t' \int{\rm d}^3\boldsymbol{r}' \, {\boldsymbol{\mu}}(\boldsymbol{r},\boldsymbol{r}',t-t')\,{\boldsymbol{\cdot}}\,\boldsymbol{H}(\boldsymbol{r}',t'), \end{equation} \tag{ 5 }$$

where ε(r,r',t − t') is the relative permittivity and $\boldsymbol{\mu}(\boldsymbol{r}, \boldsymbol{r}^{\prime }, t -t^{\prime})$ is the relative permeability. If the system as a whole exhibits time-reversal symmetry (i.e. the medium is neither in motion, nor is there an external magnetic field), then the dielectric functions obey the symmetry relation [1]

$$\begin{equation} \fl \varepsilon_{ij}(\boldsymbol{r},\boldsymbol{r}',t-t')=\varepsilon_{ji}(\boldsymbol{r}',\boldsymbol{r},t-t'), \quad \mu_{ij}(\boldsymbol{r},\boldsymbol{r}',t-t')=\mu_{ji}(\boldsymbol{r}',\boldsymbol{r},t-t'). \end{equation} \tag{ 6 }$$

Taking the simplest case first, we specialize to a homogeneous medium, and assume that the dielectric functions in (4) and (5) depend only on the difference between the coordinates, r − r'. In frequency and wavevector space the dielectric functions then depend on a single wavevector, k and are complex functions obeying [1, 9]

$$\begin{equation} \varepsilon_{ij}^*(\boldsymbol{k},\omega)=\varepsilon_{ij}(-\boldsymbol{k},-\omega), \quad \mu_{ij}^*(\boldsymbol{k},\omega)=\mu_{ij}(-\boldsymbol{k},-\omega). \end{equation} \tag{ 7 }$$

In terms of the wavevector, the symmetry relation (6) can also be written as

$$\begin{equation} \varepsilon_{ij}(\boldsymbol{k},\omega)=\varepsilon_{ji}(-\boldsymbol{k},\omega), \quad \mu_{ij}(\boldsymbol{k},\omega)=\mu_{ji}(-\boldsymbol{k},\omega). \end{equation} \tag{ 8 }$$

To further simplify the situation, we assume an isotopic dielectric with a centre of symmetry, where the permittivity takes the form [1, 9]

$$\begin{equation} \boldsymbol{\varepsilon}(\boldsymbol{k},\omega)=\varepsilon_\perp(k,\omega)\left(\mathbbm{1}-\frac{\boldsymbol{k}\otimes\boldsymbol{k}}{k^2}\right)+\varepsilon_\parallel(k,\omega)\frac{\boldsymbol{k}\otimes\boldsymbol{k}}{k^2}, \end{equation} \tag{ 9 }$$

where ε_⊥(k,ω) and ε_∥(k,ω) are the, respectively, referred to as the transverse and longitudinal permittivities, which depend only on the magnitude of the wavevector, k = |k|. Note that the Kramers–Kronig relations still hold for ε_⊥(k,ω) and ε_∥(k,ω) at each value of k, as a consequence of the delayed time response in (4).

The electromagnetic response captured by the tensor permittivity (9) can also be implemented by an (ω,k)-dependent scalar permittivity and permeability [9]. This is seen by confirming that the macroscopic Maxwell equations with permittivity (9) (μ = 1) are identical to those with scalar ε(k,ω) and μ(k,ω) given by [9]

$$\begin{equation} \varepsilon(k,\omega)=\varepsilon_\parallel(k,\omega), \quad 1-\mu^{-1}(k,\omega)=\frac{\omega^2}{c^2k^2}\left[\varepsilon_\perp(k,\omega)-\varepsilon_\parallel(k,\omega)\right]. \end{equation} \tag{ 10 }$$

The electromagnetic energy dissipated by the medium naturally splits into absorption of transverse fields and longitudinal fields [9]; this means that the transverse and longitudinal permittivities independently contribute to the absorption properties of the medium. It then follows that in a dissipative medium we must have [9]

$$\begin{equation} \mathrm{Im}\,\varepsilon_\perp(k,\omega)>0, \quad \mathrm{Im}\,\varepsilon_\parallel(k,\omega)>0. \end{equation} \tag{ 11 }$$

Interestingly, the same dissipative medium, when described by ε(k,ω) and μ(k,ω) in (10), does not in general have Im μ(k,ω) > 0 [9]. This is clear from imposing (11) on (10), which gives Im ε(k,ω) > 0 but does not constrain the sign of Im μ(k,ω). The imaginary part of μ(k,ω) can be positive or negative in a dissipative medium because ε(k,ω) and μ(k,ω) do not govern the absorption of linearly independent components of the electromagnetic fields, in contrast to ε_⊥(k,ω) and ε_∥(k,ω).

In what follows we present the canonical quantization of a homogeneous, isotropic medium described by the permittivity (9). This is the simplest case, and the theory can be readily generalized to more complicated nonlocal dielectric functions.

Following a similar procedure to [5, 19], the quantum theory of electromagnetism in a spatially dispersive medium is constructed through initially finding an action, S, whose equations of motion are the macroscopic Maxwell equations with the constitutive relations (4) and (5).

The required expression for S is a simple generalization of the action for the macroscopic Maxwell equations in the absence of spatial dispersion [5]: S remains a functional of the electromagnetic potentials ϕ (r,t) and A (r,t), and a continuum of reservoir fields X_ω(r,t). The reservoir fields—although unobservable—are a necessary part of the description, for they allow the absorption of electromagnetic radiation to be described within a closed system. Nevertheless, the reservoir ought not to be confused with the real microscopic degrees of freedom of the medium, and we emphasize that the classical theory developed in this section will contain no more information than we can glean from the macroscopic Maxwell equations.

To incorporate spatial dispersion into the picture, we use a nonlocal coupling between the electromagnetic field and the reservoir fields, finding the required action to be of the form

$$\begin{equation} S\,[\phi,\boldsymbol{A},\boldsymbol{X}_\omega]=S_{\mathrm{em}}\,[\phi,\boldsymbol{A}]+S_{\mathrm {X}}\,[\boldsymbol{X}_\omega]+S_{\mathrm{int}}\,[\phi,\boldsymbol{A},\boldsymbol{X}_\omega], \end{equation} \tag{ 12 }$$

where S_em is the free electromagnetic action

$$\begin{equation} S_{\mathrm{em}}[\phi,\boldsymbol{A}]=\frac{\varepsilon_{0}}{2}\int{\rm d}^4 x\left(\boldsymbol{E}^{2}-c^{2}\boldsymbol{B}^{2}\right), \end{equation} \tag{ 13 }$$

where the fields are related to the potentials by E = −∇ϕ − ∂_tA and B = ∇ × A. The quantity S_X is the action for the free reservoir oscillators,

$$\begin{equation} S_{\mathrm {X}}\,[\boldsymbol{X}_\omega]=\frac{1}{2}\int{\rm d}^4 x\int_0^\infty\mathrm{d}\omega\left[(\partial_t\boldsymbol{X}_\omega)^{2}-\omega^2\boldsymbol{X}_\omega^{2}\right], \end{equation} \tag{ 14 }$$

and S_int is the interaction term, coupling the electromagnetic fields to the reservoir

$$\begin{equation} \fl S_{\mathrm{int}}\,[\phi,\boldsymbol{A},\boldsymbol{X}_\omega]=\int_{-\infty}^\infty{\rm d} t \int{\rm d}^3 \boldsymbol{r} \int{\rm d}^3 \boldsymbol{r}' \int_0^\infty\mathrm{d}\omega\, \boldsymbol{X}_\omega(\boldsymbol{r},t)\,{\boldsymbol{\cdot}}\,\boldsymbol{F}(\boldsymbol{r}-\boldsymbol{r}',\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{E}(\boldsymbol{r}',t). \end{equation} \tag{ 15 }$$

The interaction (15) features a spatially nonlocal coupling bi-tensor F(r − r',ω) that depends on two scalar functions α₁ (|r − r'|,ω) and α₂(|r − r'|,ω),

$$\begin{equation} \boldsymbol{F}(\boldsymbol{r}-\boldsymbol{r}',\omega)=\alpha_1(|\boldsymbol{r}-\boldsymbol{r}'|,\omega)\mathbbm{1} +\boldsymbol{\nabla}\alpha_2(|\boldsymbol{r}-\boldsymbol{r}'|,\omega)\otimes\stackrel{\leftarrow}{\boldsymbol{\nabla'}}. \end{equation} \tag{ 16 }$$

The scalar functions appearing in (16) are defined in k-space in terms of the imaginary parts of the transverse and longitudinal permittivities of an arbitrary homogeneous medium

$$\begin{equation} \alpha_1(k,\omega)=\left[\frac{2\varepsilon_0}{\pi}\omega\,\mathrm{Im}\,\varepsilon_\perp(k,\omega)\right]^{1/2}, \end{equation} \tag{ 17 }$$

$$\begin{equation} \alpha_2(k,\omega)=\frac{1}{k^2}\sqrt{\frac{2\varepsilon_0\omega}{\pi}}\left[\sqrt{\mathrm{Im}\,\varepsilon_\parallel(k,\omega)}-\sqrt{\mathrm{Im}\,\varepsilon_\perp(k,\omega)}\right]. \end{equation} \tag{ 18 }$$

The real parts of the permittivities appearing in (17) and (18) are determined by the imaginary parts through the Kramers–Kronig relation [1, 17]

$$\begin{equation} \mathrm{Re}\,\varepsilon_{\{\perp,\parallel\}}(k,\omega')-1=\frac{2}{\pi}\mathrm{P}\int_0^\infty\mathrm{d}\omega\frac{\omega\,\mathrm{Im}\,\varepsilon_{\{\perp,\parallel\}}(k,\omega)}{\omega^2-\omega'^2}. \end{equation} \tag{ 19 }$$

Therefore the appearance of only the imaginary parts of the dielectric function within (12) should not be surprising, for we do not have the freedom to choose the real and imaginary parts of the susceptibilities independently. As the permittivities depend on the magnitude of the distance (or wavevector), it is evident from (16) that the nonlocal coupling tensor, F possesses the symmetry, F_ij(r − r',ω) = F_ji(r' − r,ω). In the limiting case of a homogeneous medium lacking spatial dispersion, ε_⊥(k,ω) = ε_∥(k,ω) = ε(ω), the action (12)–(18) reduces to that in [5], for the particular case where μ = 1.

We now show that the classical equations of motion derived from (12) are the macroscopic Maxwell equations in a spatially dispersive medium. Variation of the action with respect to the potentials ϕ and A leads to

$$\begin{equation} \varepsilon_0\boldsymbol{\nabla}\,{\boldsymbol{\cdot}}\,\boldsymbol{E}+\int_0^\infty\mathrm{d}\omega\int{\rm d}^3 \boldsymbol{r}'\,\boldsymbol{\nabla}\,{\boldsymbol{\cdot}}\,\boldsymbol{F}(\boldsymbol{r}-\boldsymbol{r}',\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{X}_\omega(\boldsymbol{r}',t)=0, \end{equation} \tag{ 20 }$$

$$\begin{equation} -\frac{1}{\mu_0}\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{B}+\varepsilon_0\partial_t\boldsymbol{E}+\int_0^\infty\mathrm{d}\omega\int{\rm d}^3 \boldsymbol{r}'\,\boldsymbol{F}(\boldsymbol{r}-\boldsymbol{r}',\omega)\,{\boldsymbol{\cdot}}\,\partial_t\boldsymbol{X}_\omega(\boldsymbol{r}',t)=0, \end{equation} \tag{ 21 }$$

while variation with respect to X_ω gives

$$\begin{equation} -\partial_t^2\boldsymbol{X}_\omega-\omega^2\boldsymbol{X}_\omega+\int{\rm d}^3 \boldsymbol{r}'\,\boldsymbol{F}(\boldsymbol{r}-\boldsymbol{r}',\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{E}(\boldsymbol{r}',t)=0. \end{equation} \tag{ 22 }$$

The equation of motion (22) is that of a field of driven harmonic oscillators, which has the general solution

$$\begin{equation} \boldsymbol{X}_\omega\,(\boldsymbol{k},\omega')=\frac{\boldsymbol{F}\,(\boldsymbol{k},\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{E\,}(\boldsymbol{k},\omega')}{\omega^2-(\omega'+{\rm i} 0^+)^2}+\delta(\omega-\omega')\,\boldsymbol{Z}_\omega\,(\boldsymbol{k})+\delta(\omega+\omega')\,\boldsymbol{Z}^*_\omega\,(\boldsymbol{k}). \end{equation} \tag{ 23 }$$

The infinitesimal positive number 0⁺ in (23) amounts to a retarded boundary condition on the motion of the oscillator (i.e. the oscillator responds to the behaviour of the electric field in the past), while Z_ω(k) is an arbitrary function that, in the classical case, can be chosen to specify the configuration of the reservoir at some initial time, t = t₀.

The equations satisfied by the electromagnetic fields (20) and (21) can be put in the desired form through applying the result (23). This is most easily carried out in the Fourier domain. The Fourier-transformed versions of (20) and (21) both contain the quantity $\int _0^\infty \mathrm {d}\omega \boldsymbol {F}(\boldsymbol {k},\omega ){\boldsymbol {\cdot }}\boldsymbol {X}_\omega (\boldsymbol {k},\omega ')$, which is evaluated using (23) as follows. Firstly, it is clear from (16)–(18) that the square of the nonlocal coupling tensor is proportional to the dissipation in the medium

$$\begin{equation} \boldsymbol{F}(\boldsymbol{k},\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{F}(\boldsymbol{k},\omega)=\frac{2\varepsilon_0\omega}{\pi}\mathrm{\,Im}\boldsymbol{\varepsilon}\,(\boldsymbol{k},\omega), \end{equation} \tag{ 24 }$$

where ε(k,ω) is the (tensor) permittivity (9). Rewriting the first term in (23) using the identity

$$\begin{equation} \frac{1}{\omega^2-(\omega'+{\rm i} 0^+)^2}={P}\left(\frac{1}{\omega^2-\omega'^2}\right)+\frac{{\rm i}\pi}{2\omega}[\delta(\omega-\omega')-\delta(\omega+\omega')], \end{equation} \tag{ 25 }$$

and employing (24) and the Kramers–Kronig relation (19), we find the required expression

$$\begin{equation} \fl\int_0^\infty\mathrm{d}\omega\,\boldsymbol{F}(\boldsymbol{k},\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{X}_\omega(\boldsymbol{k},\omega')=\varepsilon_0\left[\boldsymbol{\varepsilon}(\boldsymbol{k},\omega')-\mathbbm{1}\right]\,{\boldsymbol{\cdot}}\,\boldsymbol{E}(\boldsymbol{k},\omega')+\boldsymbol{F}(\boldsymbol{k},\omega')\,{\boldsymbol{\cdot}}\,\boldsymbol{Z}_{\omega}'(\boldsymbol{k}). \end{equation} \tag{ 26 }$$

Substituting this result (26) in the Fourier-transformed versions of (20) and (21), we obtain the following form for the electromagnetic field equations:

$$\begin{equation} {\rm i} \boldsymbol{k}\,{\boldsymbol{\cdot}}\,\boldsymbol{D}(\boldsymbol{k},\omega)=\sigma(\boldsymbol{k},\omega), \quad {\rm i} \boldsymbol{k}\boldsymbol{\times}\boldsymbol{H}(\boldsymbol{k},\omega)=-{\rm i}\,\omega\boldsymbol{D}(\boldsymbol{k},\omega)+\boldsymbol{j}(\boldsymbol{k},\omega), \end{equation} \tag{ 27 }$$

where the D and H fields are given by

$$\begin{equation} \boldsymbol{D}(\boldsymbol{k},\omega)=\varepsilon_0\boldsymbol{\varepsilon}(\boldsymbol{k},\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{E}(\boldsymbol{k},\omega), \quad \boldsymbol{H}(\boldsymbol{k},\omega)=\mu_0^{-1}\boldsymbol{B}(\boldsymbol{k},\omega), \end{equation} \tag{ 28 }$$

and the charge and current densities are related to the initial configuration of the reservoir

$$\begin{equation} \sigma(\boldsymbol{k},\omega)=-{\rm i} \boldsymbol{k}\,{\boldsymbol{\cdot}}\,\boldsymbol{F}(\boldsymbol{k},\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{Z}_{\omega}(\boldsymbol{k}), \quad \boldsymbol{j}(\boldsymbol{k},\omega)=-{\rm i}\,\omega \boldsymbol{F}(\boldsymbol{k},\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{Z}_{\omega}(\boldsymbol{k}). \end{equation} \tag{ 29 }$$

As anticipated, equations (27) and (28) are the macroscopic Maxwell equations in a homogeneous, spatially dispersive medium with permittivity (9) (the other two Maxwell equations are identities in terms of the electromagnetic potentials). It is evident from (29) that the charge conservation law

$$\begin{equation} \partial_t\sigma(\boldsymbol{r},t)+ \boldsymbol{\nabla}\,{\boldsymbol{\cdot}}\, \boldsymbol{j}(\boldsymbol{r},t)=0 \end{equation} \tag{ 30 }$$

is identically true.

Having shown that the action (12) reproduces the classical theory of electromagnetism in a spatially dispersive medium, we are now in a position to develop the quantum theory. Before doing so we give the expression for the electromagnetic Green function in the medium, which will be seen to determine the properties of the field operators in the quantum theory (see section 5). The electric field within the medium satisfies an inhomogeneous wave equation, which can be derived through substituting B(k,ω) = k × E(k,ω)/ω, into the second of (27)

$$\begin{equation} \fl \left[\left( \frac{\omega^2}{c^2}\varepsilon_\perp(k,\omega)-k^2\right)\left(\mathbbm{1}-\frac{\boldsymbol{k}\otimes\boldsymbol{k}}{k^2}\right)+ \frac{\omega^2}{c^2k^2}\varepsilon_\parallel(k,\omega)\boldsymbol{k}\otimes\boldsymbol{k} \right]\,{\boldsymbol{\cdot}}\,\boldsymbol{E}(\boldsymbol{k},\omega)=-{\rm i}\mu_0\omega\boldsymbol{j}(\boldsymbol{k},\omega). \end{equation} \tag{ 31 }$$

This has the solution

$$\begin{equation} \boldsymbol{E}(\boldsymbol{k},\omega)={\rm i}\mu_0\omega\boldsymbol{G}(\boldsymbol{k},\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{j}(\boldsymbol{k},\omega), \end{equation} \tag{ 32 }$$

where G(k,ω) is the inverse of the matrix on the left hand side of (31), known as the Green bi-tensor

$$\begin{equation} \fl \boldsymbol{G}(\boldsymbol{k},\omega)=-\left[ \frac{\omega^2}{c^2}\varepsilon_\perp(k,\omega)-k^2\right]^{-1}\left(\mathbbm{1}-\frac{\boldsymbol{k}\otimes\boldsymbol{k}}{k^2}\right)- \frac{c^2}{\omega^2k^2\varepsilon_\parallel(k,\omega)}\boldsymbol{k}\otimes\boldsymbol{k}. \end{equation} \tag{ 33 }$$

For reference, we note that the Green bi-tensor (33) satisfies

$$\begin{equation} \frac{\omega^2}{c^2}\boldsymbol{G}^*(\boldsymbol{k},\omega) \,{\boldsymbol{\cdot}}\, \mathrm{Im}\,{\boldsymbol{\varepsilon}}(\boldsymbol{k},\omega) \,{\boldsymbol{\cdot}}\, \boldsymbol{G}(\boldsymbol{k},\omega)= \mathrm{Im}\,\boldsymbol{G}(\boldsymbol{k},\omega) . \end{equation} \tag{ 34 }$$

The quantization of the field theory given in (12)–(18) does not fundamentally differ from [5]. The distinction between this theory and [5] is that here we have no coupling between the reservoir and the magnetic field, and the coupling of the reservoir to the electric field is spatially nonlocal. The procedure is nevertheless the same: impose the canonical commutation relations between the fields and their associated canonical momenta, and construct the Hamiltonian operator from the action, (12).

The canonical momenta of the dynamical variables ϕ, A and X_ω are, respectively,

$$\begin{equation} \fl {\Pi}_\phi=0, \quad \boldsymbol{\Pi}_{\mathrm{A}}=-\varepsilon_0\boldsymbol{E}-\int_0^\infty \mathrm{d}\omega\int{\rm d}^3 \boldsymbol{r}'\,\boldsymbol{F}(\boldsymbol{r}-\boldsymbol{r}',\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{X}_\omega(\boldsymbol{r}',t)\,\ \mathrm{and} \quad \boldsymbol{\Pi}_{X_\omega}=\partial_t\boldsymbol{X}_\omega. \end{equation} \tag{ 35 }$$

Following the standard Coulomb-gauge formulation of QED [18], we eliminate ϕ as a dynamical variable and impose the constraints

$$\begin{equation} \boldsymbol{\nabla}\,{\boldsymbol{\cdot}}\, \boldsymbol{A}=0\,\ \mathrm{and} \quad \boldsymbol{\nabla}\,{\boldsymbol{\cdot}}\,\boldsymbol{\Pi}_ A=0. \end{equation} \tag{ 36 }$$

With these constraints the equal–time commutation relation between $\hat {\boldsymbol {A}}$ and $\hat {\boldsymbol {\boldsymbol {\Pi }}}_{\!A}$ take the usual form [10, 18]

$$\begin{equation} \left[ \hat{\boldsymbol{A}}(\boldsymbol{r},t),\hat{\boldsymbol{\boldsymbol{\Pi}}}_{\!A} (\boldsymbol{r}',t) \right]={\rm i}\hbar \, \boldsymbol{\delta}_{\mathrm{T}}(\mathbf{r}-\mathbf{r}'), \end{equation} \tag{ 37 }$$

where δ_T(r − r') is the transverse delta function

$$\begin{equation} \boldsymbol{\delta}_{\mathrm{T}}(\boldsymbol{r}-\boldsymbol{r}')=\mathbbm{1}\delta(\boldsymbol{r}-\boldsymbol{r}')+\boldsymbol{\nabla} \otimes \boldsymbol{\nabla} \left(\frac{1}{4\pi|\boldsymbol{r}-\boldsymbol{r}'|}\right). \end{equation} \tag{ 38 }$$

No such constraints are imposed on the reservoir, and we have the usual commutation relations

$$\begin{equation} \left[\hat{\boldsymbol{X}}_{\omega}(\boldsymbol{r},t),\hat{\boldsymbol{\boldsymbol{\Pi}}}_{\!X_{\omega}'} (\boldsymbol{r}',t)\right]={\rm i}\hbar\mathbbm{1}\delta\,(\omega-\omega')\,\delta\,(\boldsymbol{r}-\boldsymbol{r}'). \end{equation} \tag{ 39 }$$

The constraints (36) give the following division of the electric-field operator into transverse and longitudinal parts, $\hat {\boldsymbol {E}}=\hat {\boldsymbol {E}}_{\mathrm {T}}+\hat {\boldsymbol {E}}_{\mathrm {L}}$:

$$\begin{eqnarray} \hat{\boldsymbol{E}}_{\mathrm{T}}&=&-\partial_t\hat{\boldsymbol{A}},\nonumber\\ \hat{\boldsymbol{E}}_{\mathrm{L}}&=&-\boldsymbol{\nabla}\hat{\phi}=-\frac{1}{\varepsilon_0}\left[\int_0^\infty \mathrm{d}\omega\int{\rm d}^3 \boldsymbol{r}'\,\boldsymbol{F}(\boldsymbol{r}-\boldsymbol{r}',\omega)\,{\boldsymbol{\cdot}}\,\hat{\boldsymbol{X}}_\omega(\boldsymbol{r}',t)\right]_{\mathrm{L}}. \end{eqnarray} \tag{ 40 }$$

The construction of the Hamiltonian follows the steps explained in detail within [5] for a very similar action, with the result

$$\begin{eqnarray} &&\fl \hat{H}=\int{\rm d}^3\boldsymbol{r}\Bigg\{\frac{1}{2\varepsilon_0}\left[\hat{\boldsymbol{\boldsymbol{\Pi}}}_{\!A}+\int_0^\infty\! \mathrm{d}\omega\int\! {\rm d}^3 \boldsymbol{r}'\,\boldsymbol{F}(\boldsymbol{r}-\boldsymbol{r}',\omega)\,{\boldsymbol{\cdot}}\,\hat{\boldsymbol{X}}_\omega(\boldsymbol{r}',t)\right]^{2}+\frac{1}{2\mu_0}(\boldsymbol{\nabla}\boldsymbol{\times}\hat{\boldsymbol{A}})^2 \nonumber \\ &&+\frac{1}{2}\int_0^\infty\mathrm{d}\omega\left[\hat{\boldsymbol{\boldsymbol{\Pi}}}_{\!X_{\omega}}^2+\omega^2\hat{\boldsymbol{X}}_\omega^2\right] \Bigg\}. \end{eqnarray} \tag{ 41 }$$

In the absence of any magnetic properties the Hamiltonian (41) is of the same form as that in [5] with the replacement $\alpha (\boldsymbol {r},\omega )\,\hat {\boldsymbol {X}}_{\omega }\,(\boldsymbol {r})\to \int \mathrm {d}\boldsymbol {r}^{\prime }\boldsymbol {F}\,(\boldsymbol {r}-\boldsymbol {r}^{\prime },\omega )\,\hat {\boldsymbol {X}}_{\omega }(\boldsymbol {r}^{\prime })$. As one might expect, with the constraints (36), the equations of motion for the operators are formally identical to the classical equations of motion (20)–(22). Due to the similarity of the form of the Hamiltonian, the demonstration of this is almost identical to that given in [5].

The Hamiltonian (41) is a quadratic combination of the field operators, and can be brought into the diagonal form

$$\begin{equation} \hat{H}=\int{\rm d}^3 \boldsymbol{r}\int_0^\infty\mathrm{d}\omega\,\hbar\omega\hat{\boldsymbol{C}}^\dagger(\boldsymbol{r},\omega)\ \,{\boldsymbol{\cdot}}\,\ \hat{\boldsymbol{C}}(\boldsymbol{r},\omega). \end{equation} \tag{ 42 }$$

This form of Hamiltonian amounts to the statement that the total energy of the system is equal to the integral over frequency of the number of quanta in the normal modes, weighted by the energy ℏω. We may think of an excitation of each normal mode as something like a 'polariton', which in this context would refer to a point-like excitation within the medium, and an associated electromagnetic field. These point-like excitations may seem strange considering the discussion given in the introduction. However the inclusion of spatial dispersion has the consequence that these excitations are now nonlocally related to the current and charge density within the medium. For instance, the state $\hat {\boldsymbol {C}}^\dagger _{\omega }\,(\boldsymbol {r},t)|0\rangle$ now corresponds to an excitation of a current within the medium that is smeared out over a region of space centred at r. These 'polaritons' are bosons, and the vector creation $\hat {\boldsymbol {C}}^\dagger (\boldsymbol {r},\omega )$ and annihilation $\hat {\boldsymbol {C}}(\boldsymbol {r},\omega )$ operators therefore obey the commutation relations

$$\begin{equation} \left[\hat{\boldsymbol{C}}(\boldsymbol{r},\omega),\hat{\boldsymbol{C}}^\dagger(\boldsymbol{r}',\omega')\right]=\mathbbm{1} \delta\,(\omega-\omega')\, \delta\,(\boldsymbol{r}-\boldsymbol{r}'). \end{equation} \tag{ 43 }$$

The time-dependence of these operators is given by

$$\begin{equation} \hat{\boldsymbol{C}}(\boldsymbol{r},t,\omega)=\hat{\boldsymbol{C}}(\boldsymbol{r},\omega)\mathrm{e}^{-{\rm i}\,\omega t}, \end{equation} \tag{ 44 }$$

which is consistent with (42) and (43). The diagonalization procedure leading to (42) is described in detail in [5] (more general versions of the procedure, required for microscopic models of dielectrics, are given in [19, 20]). Given the similarity to [5], we confine ourselves to a sketch of the argument.

The transformation of (41) into (42) is achieved through a canonical transformation of the fields, $\hat {\boldsymbol {A}}$, $\hat {\boldsymbol {\Pi }}_{\boldsymbol {A}}$, $\hat {\boldsymbol {X}}_{\omega }$ and $\hat {\boldsymbol {\Pi }}_{\boldsymbol {X}_{\omega }}$. In general this transformation can be written in the form,

$$\begin{equation} \hat{\boldsymbol{A}}(\boldsymbol{r},t)=\int{\rm d}^3 \boldsymbol{r}'\int_0^\infty\mathrm{d}\omega\, [\boldsymbol{f}_A(\boldsymbol{r}-\boldsymbol{r}',\omega)\ \,{\boldsymbol{\cdot}}\,\ \hat{\boldsymbol{C}}(\boldsymbol{r}',t,\omega)+{\rm h.c.}], \end{equation} \tag{ 45 }$$

$$\begin{equation} \hat{\boldsymbol{\boldsymbol{\Pi}}}_{\!A}(\boldsymbol{r},t)=\int{\rm d}^3 \boldsymbol{r}'\int_0^\infty\mathrm{d}\omega\, [\boldsymbol{f}_{{\Pi}_A}(\boldsymbol{r}-\boldsymbol{r}',\omega)\,{\boldsymbol{\cdot}}\,\hat{\boldsymbol{C}}(\boldsymbol{r}',t,\omega)+{\rm h.c.}], \end{equation} \tag{ 46 }$$

$$\begin{equation} \hat{\boldsymbol{X}}_\omega(\boldsymbol{r},t)=\int{\rm d}^3 \boldsymbol{r}'\int_0^\infty\mathrm{d}\omega^{\prime}\, [\boldsymbol{f}_{X}(\boldsymbol{r}-\boldsymbol{r}',\omega,\omega')\,{\boldsymbol{\cdot}}\,\hat{\boldsymbol{C}}(\boldsymbol{r}',t,\omega')+{\rm h.c.}], \end{equation} \tag{ 47 }$$

$$\begin{equation} \hat{\boldsymbol{\boldsymbol{\Pi}}}_{\!X_{\omega}}(\boldsymbol{r},t)=\int{\rm d}^3 \boldsymbol{r}'\int_0^\infty\mathrm{d}\omega^{\prime}\, [\boldsymbol{f}_{{\Pi}_X}(\boldsymbol{r}-\boldsymbol{r}',\omega,\omega')\,{\boldsymbol{\cdot}}\,\hat{\boldsymbol{C}}(\boldsymbol{r}',t,\omega')+{\rm h.c.}], \end{equation} \tag{ 48 }$$

where the f_{{A,ΠA,X,ΠX}} are c-number bi-tensors. It is also useful to define the bi-tensor f_E that gives the expansion of the electric-field operator, $\hat {\boldsymbol {E}}$ in terms of the $\hat {\boldsymbol {C}}^\dagger (\boldsymbol {r},\omega )$ and $\hat {\boldsymbol {C}}(\boldsymbol {r},\omega )$ operators. Inserting (46) and (47) into our earlier expression for the electric field operator (35), we find this to be

$$\begin{eqnarray} &&\fl \boldsymbol{f}_{E}(\boldsymbol{r}-\boldsymbol{r}',\omega)=-\frac{1}{\varepsilon_0}\boldsymbol{f}_{{\Pi}_{A}}(\boldsymbol{r}-\boldsymbol{r}',\omega)-\frac{1}{\varepsilon_0}\int_0^\infty \! \mathrm{d}\omega' \! \! \int \! {\rm d}^3 \boldsymbol{r}''\,\boldsymbol{F}(\boldsymbol{r}-\boldsymbol{r}'',\omega')\,{\boldsymbol{\cdot}}\,\boldsymbol{f}_{X}(\boldsymbol{r}''-\boldsymbol{r}',\omega',\omega).\nonumber\\ \end{eqnarray} \tag{ 49 }$$

The commutation relations satisfied by the various operators imply that the inverse of the transformation (45)–(48) is given by

$$\begin{eqnarray} &&\fl \hat{\boldsymbol{C}}(\boldsymbol{r},t,\omega)=-\frac{{\rm i}}{\hbar}\int{\rm d}^3 \boldsymbol{r}'\left\{ \vphantom{+\int_0^\infty \mathrm{d}\omega'\left[\hat{\boldsymbol{X}}_{\omega'}(\boldsymbol{r}',t)\,{\boldsymbol{\cdot}}\,\boldsymbol{f}^*_{{\Pi}_X}(\boldsymbol{r}'-\boldsymbol{r},\omega',\omega)-\hat{\boldsymbol{\boldsymbol{\Pi}}}_{\!X_{\omega'}}(\boldsymbol{r}',t)\,{\boldsymbol{\cdot}}\,\boldsymbol{f}^*_{X}(\boldsymbol{r}'-\boldsymbol{r},\omega',\omega)\right]}\hat{\boldsymbol{A}}(\boldsymbol{r}',t)\,{\boldsymbol{\cdot}}\,\boldsymbol{f}^*_{{\Pi}_{A}}(\boldsymbol{r}'-\boldsymbol{r},\omega)-\hat{\boldsymbol{\boldsymbol{\Pi}}}_{\!A}(\boldsymbol{r}',t)\,{\boldsymbol{\cdot}}\,\boldsymbol{f}^*_{A}(\boldsymbol{r}'-\boldsymbol{r},\omega) \right. \nonumber \\ &&\left. +\int_0^\infty \mathrm{d}\omega'\left[\hat{\boldsymbol{X}}_{\omega'}(\boldsymbol{r}',t)\,{\boldsymbol{\cdot}}\,\boldsymbol{f}^*_{{\Pi}_X}(\boldsymbol{r}'-\boldsymbol{r},\omega',\omega)-\hat{\boldsymbol{\boldsymbol{\Pi}}}_{\!X_{\omega'}}(\boldsymbol{r}',t)\,{\boldsymbol{\cdot}}\,\boldsymbol{f}^*_{X}(\boldsymbol{r}'-\boldsymbol{r},\omega',\omega)\right] \right\} . \nonumber \\ \end{eqnarray} \tag{ 50 }$$

It must now be demonstrated that it is possible to find a transformation (45)–(48) that will turn (41) into (42), without disturbing the canonical commutation relations. This transformation should also be such that (50) is consistent with (45)–(48). We do this briefly in words—a more rigorous argument can be found in [5].

Given that the field operators (45)–(48) formally satisfy the classical equations of motion (20)–(22), and the $\hat {\boldsymbol {C}}_{\omega }$ have an assumed time dependence given by (44), the transformation can only work if there is a direct relation between the solutions to the classical equations of motion in the frequency domain and the f_{{A,ΠA,X,ΠX}}. For example, the frequency-domain expression for the classical electric field is given by the Fourier transform of (32) with respect to k. A comparison between the operator expression for $\hat {\boldsymbol {E}}$ and its classical counterpart then shows that the $\hat {\boldsymbol {C}}_{\omega }$ operator is analogous to the classical quantity Z_ωe^−iωt introduced in (23). One can therefore replace Z_ωe^−iωt with $\zeta (\omega )\,\hat {\boldsymbol {C}}_{\omega }$ in the solution to the classical problem and thereby obtain (45)–(48): having done this, the operators will satisfy the classical equations of motion, and this will be automatically consistent with (42). The only quantity left to determine is ζ (ω), which is a function of frequency chosen so that the commutation relations (37), (39) and (43) are all simultaneously satisfied. The necessary analysis can be found in [21], where it is determined that $\zeta \,(\omega )=2\pi \sqrt {\hbar /2\omega }$. The quantity f_E is then equal to

$$\begin{equation} \boldsymbol{f}_E(\boldsymbol{k},\omega)=\mu_0\omega\sqrt{\frac{\hbar\omega}{2}}\,\boldsymbol{G}(\boldsymbol{k},\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{F}(\boldsymbol{k},\omega). \end{equation} \tag{ 51 }$$

Similarly the bi-tensorial coefficients associated with the electromagnetic operators are

$$\begin{equation} \boldsymbol{f}_A(\boldsymbol{k},\omega)=-\frac{{\rm i}}{\omega}\left[\boldsymbol{f}_E(\boldsymbol{k},\omega)\right]_{\mathrm{T}}, \end{equation} \tag{ 52 }$$

$$\begin{equation} \boldsymbol{f}_{{\Pi}_{A}}(\boldsymbol{k},\omega)= -\varepsilon_0\boldsymbol{\varepsilon}(\boldsymbol{k},\omega)\,{\boldsymbol{\cdot}}\,\boldsymbol{f}_E(\boldsymbol{k},\omega)-\sqrt{\frac{\hbar}{2\omega}}\,\boldsymbol{F}(\boldsymbol{k},\omega) \end{equation} \tag{ 53 }$$

and those of the reservoir are

$$\begin{equation} \fl\boldsymbol{f}_{X}(\boldsymbol{k},\omega',\omega)= \frac{{\rm i}}{\omega}\boldsymbol{f}_{{\Pi}_X}(\boldsymbol{k},\omega',\omega), \end{equation} \tag{ 54 }$$

$$\begin{eqnarray} &&\fl \boldsymbol{f}_{{\Pi}_X}(\boldsymbol{k},\omega',\omega)= -\frac{{\rm i}\,\omega}{2\omega'}\left[ \frac{1}{\omega'-\omega-{\rm i} 0^+} + \frac{1}{\omega'+\omega} \right]\boldsymbol{F}(\boldsymbol{k},\omega')\,{\boldsymbol{\cdot}}\,\boldsymbol{f}_E(\boldsymbol{k},\omega) -{\rm i}\sqrt{\frac{\hbar\omega}{2}}\,\mathbbm{1}\delta(\omega-\omega').\nonumber \\ \end{eqnarray} \tag{ 55 }$$

The external charge and current densities (29) appearing in the classical problem are now operators

$$\begin{equation} \hat{\sigma}(\boldsymbol{k},\omega)=-2\pi{\rm i}\left[\frac{\hbar\varepsilon_0}{\pi}\,\mathrm{Im}\,\varepsilon_\parallel(k,\omega)\right]^{1/2}\boldsymbol{k}\,{\boldsymbol{\cdot}}\,\hat{\boldsymbol{C}}(\boldsymbol{k},\omega),\end{equation} \tag{ 56 }$$

$$\begin{equation} \hat{\boldsymbol{\jmath}}(\boldsymbol{k},\omega)=-2\pi{\rm i}\sqrt{\frac{\hbar\omega}{2}}\,\boldsymbol{F}(\boldsymbol{k},\omega) \,{\boldsymbol{\cdot}}\,\hat{\boldsymbol{C}}(\boldsymbol{k},\omega). \end{equation} \tag{ 57 }$$

From (43), the k-space creation and annihilation operators have the commutation relations

$$\begin{equation} \left[\hat{\boldsymbol{C}}(\boldsymbol{k},\omega),\hat{\boldsymbol{C}}^\dagger(\boldsymbol{k}',\omega')\right]=(2\pi)^3\mathbbm{1} \delta\,(\omega-\omega')\,\ \delta\,(\boldsymbol{k}-\boldsymbol{k}'), \end{equation} \tag{ 58 }$$

so that the charge and current operators (56) satisfy

$$\begin{equation} \left[\hat{\sigma}(\boldsymbol{k},\omega),\hat{\sigma}^\dagger(\boldsymbol{k}',\omega')\right]=(2\pi)^5\frac{\hbar k^{2}\varepsilon_0}{\pi}\mathrm{\,Im}\,\varepsilon_{\parallel}(k,\omega)\delta(\omega-\omega')\, \delta\,(\boldsymbol{k}-\boldsymbol{k}'),\end{equation} \tag{ 59 }$$

$$\begin{equation} \left[\hat{\boldsymbol{\jmath}}(\boldsymbol{k},\omega),\hat{\boldsymbol{\jmath}}^\dagger(\boldsymbol{k}',\omega')\right]=(2\pi)^5\frac{\hbar\omega^2\varepsilon_0}{\pi}\mathrm{Im}\,\boldsymbol{\varepsilon}(\boldsymbol{k},\omega)\delta(\omega-\omega')\, \delta\,(\boldsymbol{k}-\boldsymbol{k}'), \end{equation} \tag{ 60 }$$

where we have applied (24). It is interesting that, without any dependence of ε_∥ on k, (59) implies that the correlation between the charge density within the medium grows with k. When Im [ε] depends on k such that it tends to zero as k → ∞, then the commutation relation (60) will modify the correlation function for the current density from the delta function given in (3) to a correlation that is spread out over a finite region of space. This is the anticipated effect of spatial dispersion discussed in the introduction.

We are now in a position to return to the problem raised in the introduction. It was shown that the arbitrarily localized current densities appearing in the local theory of macroscopic quantum electromagnetism give rise to divergent results. Here we explicitly demonstrate that the introduction of spatial dispersion allows us to obtain finite results. Our attention is confined to the equilibrium properties of the electromagnetic field within a homogeneous medium.

In thermal equilibrium, the 'polaritons' identified through the diagonalization of the Hamiltonian in section 5, have an occupation number given by the Planck distribution [22]

$$\begin{equation} \left\langle\hat{\boldsymbol{C}}^\dagger(\boldsymbol{r},\omega)\otimes\hat{\boldsymbol{C}}(\boldsymbol{r}',\omega')\right\rangle={\mathcal{N}}\,(\omega)\,\mathbbm{1}\delta\,(\omega-\omega')\, \delta\,(\boldsymbol{r}-\boldsymbol{r}'), \end{equation} \tag{ 61 }$$

where

$$\begin{equation} \mathcal{N}(\omega):=\frac{1}{\mathrm{e}^{\hbar\omega/k_{\mathrm{B}}T}-1}, \end{equation} \tag{ 62 }$$

so that the current operator (56) has the thermal correlation

$$\begin{equation} \fl \left\langle\hat{\boldsymbol{\jmath}}^\dagger(\boldsymbol{k},\omega)\otimes\hat{\boldsymbol{\jmath}}(\boldsymbol{k}',\omega')\right\rangle=(2\pi)^5\mathcal{N}(\omega)\frac{\hbar\omega^2\varepsilon_0}{\pi}\mathrm{Im}\,\boldsymbol{\varepsilon\,}(\boldsymbol{k},\omega)\,\delta\,(\omega-\omega')\, \delta\,(\boldsymbol{k}-\boldsymbol{k}), \end{equation} \tag{ 63 }$$

where we have again applied (24). It is now straightforward to use (61) and (51) to compute the equal-time field correlation functions in thermal equilibrium. With use of the Green-function relation (34) the results take the same form as in the local theory

$$\begin{equation} \fl \left\langle\hat{\boldsymbol{E}} (\boldsymbol{r}, t)\otimes\hat{\boldsymbol{E}}(\boldsymbol{r}', t)\right \rangle=\frac{\hbar\mu_0}{\pi}\int_0^\infty\mathrm{d}\omega\,\omega^2\mathrm{coth}\left(\frac{\hbar\omega}{2k_{\mathrm{B}} T}\right)\mathrm{Im}\,\boldsymbol{G}(\boldsymbol{r}-\boldsymbol{r}',\omega),\end{equation} \tag{ 64 }$$

$$\begin{equation} \fl \left\langle\hat{\boldsymbol{B}} (\boldsymbol{r}, t)\otimes\hat{\boldsymbol{B}}(\boldsymbol{r}', t)\right \rangle=\frac{\hbar\mu_0}{\pi}\int_0^\infty\mathrm{d}\omega\,\mathrm{coth}\left(\frac{\hbar\omega}{2k_{\mathrm{B}} T}\right)\boldsymbol{\nabla}\boldsymbol{\times}\mathrm{Im}\,\boldsymbol{G}(\boldsymbol{r}-\boldsymbol{r}',\omega)\boldsymbol{\times}\stackrel{\leftarrow}{\boldsymbol{\nabla}'}, \end{equation} \tag{ 65 }$$

where the notation ${\times}\stackrel{\leftarrow}{\boldsymbol{\nabla}'}$ denotes a curl with respect to the right-hand index, i.e. $\boldsymbol{V}(\boldsymbol{r}){\times}\stackrel{\leftarrow}{\boldsymbol{\nabla}'}= \boldsymbol{\nabla} {\times}\boldsymbol{V}(\boldsymbol{r})$ for a vector field V(r).

In this nonlocal theory we again consider the intensity of the electric field at a fixed frequency, c(r,r,ω), which at T > 0 K is, from (64)

$$\begin{equation} \boldsymbol{c}(\boldsymbol{r},\boldsymbol{r},\omega)=\frac{\hbar\mu_{0}}{\pi}\omega^{2}\,\mathrm{coth}\left(\frac{\hbar\omega}{2 k_{\mathrm{B}}T}\right)\lim_{\boldsymbol{r}\to\boldsymbol{r}^{\prime}}{\rm Im}[\boldsymbol{G}(\boldsymbol{r}-\boldsymbol{r}^{\prime},\omega)] . \end{equation} \tag{ 66 }$$

For want of established expressions for the nonlocal permittivity of a real medium, we examine the behaviour of (66) in a medium described by the permittivity

$$\begin{equation} \varepsilon_{\perp}(\omega)=1-\frac{A}{\omega(\omega+\mathrm{i}\,\gamma)}, \end{equation} \tag{ 67 }$$

$$\begin{equation} \varepsilon_{\parallel}(k,\omega)=1-\frac{A}{\omega(\omega+\mathrm{i}\,\gamma)-\beta^{2}k^{2}}, \end{equation} \tag{ 68 }$$

where A and β are arbitrary constants, and γ > 0 governs the absorption into the medium. The nonlocal permittivity defined by (67) and (68) is of the same form as the hydrodynamic Drude model (see e.g. [2]), and can also be considered a limiting case of the Hopfield model [24]. In general both ε_⊥ and ε_∥ would depend on k, but in these models of homogeneous media the nonlocality is confined to ε_∥. The imaginary part of the Green function in the limit as r → r' is related to its spatial Fourier transform by

$$\begin{equation} \lim_{\boldsymbol{r}\to\boldsymbol{r}^{\prime}}{\rm Im\,}[\boldsymbol{G}(\boldsymbol{r}-\boldsymbol{r}^{\prime},\omega)]=\int\frac{\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3}}{\rm Im\,}[\boldsymbol{G}(\boldsymbol{k},\omega)] . \end{equation} \tag{ 69 }$$

Inserting (33) into (69) then gives

$$\begin{equation} \fl \lim_{\boldsymbol{r}\to\boldsymbol{r}^{\prime}}{{\rm Im\,}}[\boldsymbol{G}(\boldsymbol{r}-\boldsymbol{r}^{\prime},\omega)]=\int\frac{\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3}}\Bigg\{\frac{\left(\mathbbm{1}-\frac{\boldsymbol{k}{\otimes}\boldsymbol{k}}{k^{2}}\right)\frac{\omega^{2}}{c^{2}}{\rm Im\,}[\varepsilon_{\perp}(\omega)]}{\left|\frac{\omega^{2}}{c^{2}}\varepsilon_{\perp}(\omega)-k^{2}\right|^{2}}+\frac{\boldsymbol{k}{\otimes}\boldsymbol{k} c^{2}{\rm Im\,}[\varepsilon_{\parallel}(k,\omega)]}{\omega^{2}k^{2}|\varepsilon_{\parallel}(k,\omega)|^{2}}\Bigg\} . \end{equation} \tag{ 70 }$$

In a local description, where ε_∥ is only frequency dependent, the integral over the longitudinal part of the integrand in (70) diverges. However,

$$\begin{equation} {\rm Im\,}[\varepsilon_{\parallel}(k,\omega)]=\frac{A\omega\gamma}{(\omega^{2}-\beta^{2}k^{2})^{2}+\omega^{2}\gamma^{2}} \end{equation} \tag{ 71 }$$

tends to zero as k⁻⁴, which is sufficient to make (70) converge. We now explicitly evaluate (70), first performing the integral over the angle of the k vector,

$$\begin{equation} \int_{0}^{2\pi}\frac{\mathrm{d}\phi_{k}}{2\pi}\int_{0}^{\pi}\frac{\sin(\theta_{k})\,\mathrm{d}\theta_{k}}{2\pi}\boldsymbol{k}{\otimes}\boldsymbol{k}=\frac{k^{2}}{3\pi}\mathbbm{1} , \end{equation} \tag{ 72 }$$

where $\boldsymbol {k}=k\,[\sin (\theta _{k})\cos (\phi _{k})\hat {\boldsymbol {x}}+\sin (\theta _{k})\sin (\phi _{k})\hat {\boldsymbol {y}}+\cos (\theta _{k})\hat {\boldsymbol {z}}]$. The integrals over the magnitude of the wavevector k may be performed through extending the range of integration from [0,∞) to (− ∞,∞), closing the contour of integration in either the upper or lower half k-plane, and applying the residue theorem. This gives

$$\begin{eqnarray*} &&\int_{-\infty}^{\infty}\frac{k^{2}{\rm d}k}{4\pi^{2}}\frac{1}{\left|k^{2}-\frac{\omega^{2}}{c^{2}}\varepsilon_{\perp}(\omega)\right|^{2}}=\frac{1}{4\pi}\frac{{\mathrm{Re}}\,[\sqrt{\varepsilon_{\perp}(\omega)}]}{\frac{\omega}{c}{\rm Im\,}[\varepsilon_{\perp}(\omega)]} , \\ &&\int_{-\infty}^{\infty}\frac{k^{2}\mathrm{d}k}{4\pi^{2}}\frac{{\rm Im\,}[\varepsilon_{\parallel}(k,\omega)]}{|\varepsilon_{\parallel}(k,\omega)|^{2}}=\frac{A}{4\pi\beta^{3}}{\mathrm{Re}}\,[\sqrt{\omega^{2}-A+{\rm i}\,\omega\gamma}], \end{eqnarray*}$$

which finally leads to the following expression for the intensity of the electric field at a fixed frequency:

$$\begin{equation} \fl {\boldsymbol{c}}(\boldsymbol{r},\boldsymbol{r},\omega)=\frac{\hbar}{\pi\varepsilon_{0}}\,\mathrm{coth}\left(\frac{\hbar\omega}{2 k_{\mathrm{B}}T}\right)\mathbbm{1}\Bigg\{\frac{\omega^{3}}{6\pi c^{3}}{\mathrm{Re}}\,[\sqrt{\varepsilon_{\perp}(\omega)}]+\frac{A}{12\pi\beta^{3}}{\mathrm{Re}}\,[\sqrt{\omega^{2}-A+{\rm i}\,\omega\gamma}]\Bigg\} . \end{equation} \tag{ 73 }$$

In the limit T → 0 K, the first bracketed term is proportional to the transverse spontaneous emission rate, already found in previous studies of quantum electromagnetism in absorbing media (e.g. [23]). Meanwhile the second term—proportional to the longitudinal spontaneous emission rate—arises from the nonlocality of the longitudinal permittivity, and diverges in the local limit β → 0 as shown in (2). In the hydrodynamic Drude model, the polarizability of the medium is due to the motion of the conduction electrons, which are treated as a charged fluid. The quantity β is related to a pressure term in the equations of motion for this fluid that serves to smooth out rapid variations in the density [2]. Equation (73) shows that when the electromagnetic field is coupled to such a model medium and quantized, then this pressure term naturally regulates the divergence of the electric field intensity at a fixed frequency. Furthermore, if we make the replacement, coth(ℏω/2k_BT) → coth (ℏω/2k_BT) − 1 in (73), then a finite result may be obtained for the purely thermal contribution to the total (rather than single frequency) equilibrium field intensity in a homogeneous medium (64). It would be interesting if processes such as spontaneous emission act as a probe of the nonlocal properties of a medium through the longitudinal emission rate.

An extension of our results to inhomogeneous materials, including the case of piece-wise homogeneous configurations, is of course required to treat most problems of interest. In particular, extension of the thermal and zero-point results to materials with boundaries will give a full treatment of the Casimir effect that incorporates spatial dispersion. In this section we describe some fundamental problems in Casimir theory that can be tackled with a full account of spatial dispersion in real materials. Our purpose is to point out the deep significance of spatial dispersion for the Casimir effect, as this significance is not widely recognized in the literature.

The Lifshitz formula [25] for the Casimir force between parallel half-spaces does not take account of spatial dispersion. The change in the force due to the nonlocal response of conduction electrons in a metal has been calculated in [26, 27] and found to be small for current experimental regimes. Although the contribution of spatial dispersion is minor for Casimir forces between separate objects (except at very small separations), it is dominant for some components of the Casimir stress-energy tensor at the boundaries of those objects. In particular, if spatial dispersion is neglected, both the Casimir energy density and the lateral Casimir stress diverge at planar boundaries. This fact was noted in [12, 13] for zero-point fields, but the divergence is present also for the purely thermal stress-energy even if zero-point fields are dropped. The divergent result for the thermal energy density at material boundaries, when only the frequency dependence of the permittivity is included, is well known in surface science (see the review [14], and in particular equation (48) therein, which diverges at the surface). The behaviour of purely thermal radiation shows that the divergences in the Casimir stress-energy at material boundaries are not due to zero-point radiation specifically and therefore they are not related to the regularization of zero-point fields. The divergences are the result of integrating over evanescent fields on the boundary with arbitrarily large lateral wavevectors; the divergence occurs if the permittivity is taken to be independent of the wavevector, which is clearly an incorrect assumption when wavevectors up to infinite values are crucial in the calculation. A full treatment of the problem will require detailed knowledge of the wavevector dependence of the permittivities of real material samples.

Divergences in the Casimir stress-energy tensor also occur for curved material boundaries when spatial dispersion is ignored. There is a distorting pressure on curved boundaries due to thermal and zero-point fields and a calculation of this pressure requires the correct (finite) component of the Casimir stress tensor perpendicular to the boundary. The perpendicular Casimir stress vanishes for flat boundaries but it diverges for curved boundaries if the permittivity is taken to depend only on frequency. This diverging distortion force, or self-force, is most familiar in the case of a dielectric ball and a spherical dielectric shell [15]. The result of Boyer [28] that the zero-point self-force on a perfectly conducting, infinitely thin spherical shell is directed outwards, has no direct physical relevance since such an object does not exist. Even in the case of Boyer's shell, however, the standard, experimentally supported, Lifshitz theory shows that the self-force is in fact inwardly directed and infinite, rather than outwardly directed and finite. All these self-force divergences are again present for purely thermal radiation so they are not zero-point phenomena per se. The divergences can again be traced to the contribution of fields with arbitrarily large wavevectors along the boundary, together with the assumption that the material properties do not change with wavevector. Using a cut-off at some large value of the lateral wavevector removes the divergence [16] but would not be expected to give an accurate estimate of Casimir self-forces; instead the full details of the spatial dispersion of the object in question are needed.

Finally, the Casimir self-forces on inhomogeneous materials, where the permittivity changes continuously with position, also diverge when spatial dispersion is ignored [11]. This fact shows that the divergences at sharp boundaries, mentioned above, are not removed if the boundary is smoothed so that the permittivity decreases continuously to one [11]. In summary, a better understanding of the nonlocal response of real materials is required to address some basic problems in the theory of the Casimir effect, and these problems can also be posed for purely thermal radiation.

We have shown that it is possible to formulate consistently a canonical theory of quantum electromagnetism in nonlocal dielectric media. Although we specialized to the case of homogeneous and isotropic media, it is straightforward to extend (12) to treat the general case.

Our treatment of spatial dispersion was motivated by a basic demonstration that the local theory of macroscopic quantum electromagnetism can fail to yield finite predictions when applied within bulk media. When a model for spatial dispersion is adopted such as the hydrodynamic Drude model, or the Hopfield model, then finite results for field intensities can be obtained (section 6). In particular, it is interesting that in the local theory the purely thermal contribution to the mean field intensity is predicted to be infinite within a homogeneous medium, a result that becomes finite in our treatment. Despite these encouraging results, further work is required to determine which model is appropriate for describing the nonlocal behaviour of a given material sample.

We have pointed out that some well-known divergences that arise in Casimir theory are due to the neglect of nonlocal response (section 7). A realistic treatment of nonlocality in materials is thus essential for an accurate estimate of distortion-forces on dielectric bodies due to thermal and zero-point fields.

We acknowledge financial support from EPSRC under Program Grant EP/I034548/1.