Something Can Come of Nothing: Surface Approaches to Quantum Fluctuations and the Casimir Force

Consider two bodies with gently curved surfaces described by (single-valued) height profiles $z=H_{1}({\bf x})$ and $z=H_{2}({\bf x})$, with respect to a reference plane $\Sigma$, where ${\bf x}\equiv(x,y)$ are Cartesian coordinates on $\Sigma$ and the $z$ axis is normal to $\Sigma$. The intervening quantum field can be a scalar or electromagnetic (EM) field. For the scalar case we can impose either Dirichlet (D) or Neumann (N) boundary conditions. The EM case, for ideal (mirror) boundaries, is obtained simply as the sum of D and N cases, corresponding to the two transverse polarizations. Extension to the case of real materials, described by a complex dielectric permittivity, is possible. The only restriction on the boundary conditions is that they should describe homogeneous and isotropic materials, so that the Casimir energy is invariant under simultaneous translations and rotations of the two profiles in the plane $\Sigma$.

The gradient expansion postulates that the Casimir energy, when generalized to two surfaces and to arbitrary fields subject to arbitrary boundary conditions, is a functional $E[H_{1},H_{2}]$ of the heights $H_{1}$ and $H_{2}$, which has a derivative expansion

	$\displaystyle E[H_{1},H_{2}]=\int_{\Sigma}d{\bf x}\,U(H)\,$		$\displaystyle\left[1+\beta_{1}(H){\bf\nabla}H_{1}\cdot{\bf\nabla}H_{1}\right.$		(4)
			$\displaystyle+\left.\beta_{2}(H){\bf\nabla}H_{2}\cdot{\bf\nabla}H_{2}+\beta_{\times}(H){\bf\nabla}H_{1}\cdot{\bf\nabla}H_{2}\right.$
			$\displaystyle+\left.\beta_{\rm-}(H)\,{\hat{\bf z}}\cdot({\bf\nabla}H_{1}{\bf\times}{\bf\nabla}H_{2})+\cdots\right],$

where $H({\bf x})\equiv H_{2}({\bf x})-H_{1}({\bf x})$ is the height difference, and dots denote higher derivative terms. Here, $U(H)$ is the energy per unit area between parallel plates at separation $H$; translation and rotation symmetries in ${\bf x}$ only permit four distinct gradient coefficients at lowest order, $\beta_{1}(H)$, $\beta_{2}(H)$, $\beta_{\times}(H)$ and $\beta_{-}(H)$ . The form of such a local expansion is motivated by the existence of well-established derivative expansions of scattering amplitudes Voronovich (1994) from which the Casimir energy can be derived Rahi et al. (2009). Arbitrariness in the choice of $\Sigma$ constrains the coefficients $\beta$ in the above expansion. Invariance of $E$ under a parallel displacement of $\Sigma$ requires that all the $\beta$’s depend only on the height difference $H$, and not on the individual heights $H_{1}$ and $H_{2}$. These coefficients are further constrained by the invariance of $E$ with respect to tilting the reference plate $\Sigma$. Under a tilt of $\Sigma$ by an infinitesimal angle $\epsilon$ in the $(x,z)$ plane, the height profiles $H_{i}$ undergo a change by $\Delta H_{i}=-\epsilon\left(x+H_{i}\frac{\partial H_{i}}{\partial x}\right)$, and the requirement that $E$ not change implies

	$\displaystyle 2\left(\beta_{1}(H)+\beta_{2}(H)\right)+2\,\beta_{\times}(H)+\,H\,\frac{d\log U}{dH}-1=0\;,$
	$\displaystyle\beta_{\rm-}(H)=0\;,$		(5)

so that the non-vanishing cross term $\beta_{\times}$ is determined by $\beta_{1}$, $\beta_{2}$ and $U$.

Equations (5) indicate that, to second order in the gradient expansion, the two-surface problem reduces to that of a single curved surface facing a plane. Therefore, $U$, $\beta_{1}$ and $\beta_{2}$ can be determined by setting $H_{1}$ or $H_{2}$ to be zero. Let us set $H_{1}=0$ and define $\beta_{2}(H)\equiv\beta(H)$. We can then determine the exact functional dependence of $\beta(H)$ on $H$ by comparing the gradient expansion Eq. (4) to a perturbative expansion of the Casimir energy around flat plates, to second order in the deformation. For this purpose, we take $\Sigma$ to be a planar surface and decompose the height of the curved surface as $H({\bf x})=d+h({\bf x})$, where $d$ is chosen to be the distance of closest separation. For small deformations $|h({\bf x})|/d\ll 1$ we can expand $E[0,d+h]$ as

$$E[0,d+h]=A\,U(d)+\,\mu(d)\tilde{h}({\bf 0})+\,\int\frac{d{\bf k}}{(2\pi)^{2}}G({k};d)|{\tilde{h}}({\bf k})|^{2}\;,$$

(6)

where $A$ is the area, ${\bf k}$ is the in-plane wave-vector and ${\tilde{h}}({\bf k})$ is the Fourier transform of $h({\bf x})$. The kernel $G({k};d)$ has been evaluated by several authors: in Ref. Emig et al. (2001) for a scalar field fulfilling D or N boundary conditions on both plates, as well for the EM field satisfying ideal metal boundary conditions on both plates; and in Ref. Neto et al. (2005) for the EM field with dielectric boundary conditions. For a deformation with small slope, the Fourier transform ${\tilde{h}}({\bf k})$ is peaked around zero. Since the kernel can be expanded at least through order $k^{2}$ about $k=0$ Voronovich (1994), we define

$$G({k};d)=\gamma(d)+\delta(d)\,k^{2}\,+\cdots\;.$$

(7)

For small $h$, the coefficients in the derivative expansion can be matched with the perturbative result. By expanding Eq. (4) in powers of $h({\bf x})$ and comparing the result with the perturbative expansion to second order in both ${\tilde{h}}({\bf k})$ and $k^{2}$, we obtain

$$U^{\prime}(d)=\mu(d)\,,\quad U^{\prime\prime}(d)=2\gamma(d)\,,\quad\beta(d)=\frac{\delta(d)}{U(d)}\;,$$

(8)

with prime denoting a derivative.

For ideal boundary conditions (D or N for a scalar field, or perfect mirrors for the EM field) the energy per unit area between parallel plates with separation $d$ is given by

$$U(d)=-\alpha\frac{\pi^{2}\hbar c}{1440d^{3}}\,.$$

(9)

The overall coefficient is $\alpha=1$ for scalar field fluctuations where both surfaces have either Dirichlet or Neumann boundary conditions, $\alpha=2$ (corresponding to the two polarizations) for electromagnetic fluctuations with perfect conductor boundary conditions, and $\alpha=-7/8$ for a scalar field with mixed Neumann/Dirichlet boundaries.

Using the Eqs. 8, and appropriate perturbation series, it is then possible to compute the coefficient $\beta$ for the following five cases: a scalar field obeying D or N boundary conditions on both surfaces; D boundary conditions on the curved surface and N boundary conditions on the flat surface, or vice versa; and for the EM field with ideal metal boundary conditions. Since in all these cases the problem involves no other length apart from the separation $d$, $\beta$ is a pure number. $\beta_{\rm D}$ was first computed in Fosco et al. (2011), and found to be $\beta_{\rm D}=2/3$. Ref. Bimonte et al. (2012a) then found find $\beta_{\rm N}=2/3\,(1-30/\pi^{2})$, $\beta_{\rm DN}=2/3$, $\beta_{\rm ND}=2/3-80/7\pi^{2}$ (where the double subscripts denote the curved surface and flat surface boundary conditions respectively), and $\beta_{\rm EM}=2/3\,(1-15/\pi^{2})$. Upon solving Eqs. (II.1) one then finds $\beta_{\times}=2-\beta_{1}-\beta_{2}$, where $\beta_{1}$ and $\beta_{2}$ are chosen to be both equal to either $\beta_{\rm D}$, $\beta_{\rm N}$ or $\beta_{\rm EM}$, for the case of identical boundary conditions on the two surfaces, or rather $\beta_{1}=\beta_{\rm DN}$ and $\beta_{2}=\beta_{\rm ND}$ for the case of a scalar field obeying mixed ND boundary conditions.

Using the above values for $\beta$ and $\beta_{\times}$, it is possible to find the leading correction to PFA by explicit evaluation of Eq. (4) for the desired profiles. For example, for two spheres of radii $R_{1}$ and $R_{2}$, both with the same boundary conditions for simplicity, we obtain:

$$E=E_{\rm PFA}\left[1-\frac{d}{R_{1}+R_{2}}+(2\beta-1)\left(\frac{d}{R_{1}}+\frac{d}{R_{2}}\right)\right]\;,$$

(10)

where $E_{\rm PFA}=-(\alpha\pi^{3}\hbar cR_{1}R_{2}))/[{1440d^{2}(R_{1}+R_{2})}]$. The corresponding formula for the sphere-plane case can be obtained by taking one of the two radii to infinity. These results are in good agreement with analytic calculations in the sphere-plane system Teo et al. (2011); Teo (2013a).

The approach described above, converting a perturbative expansion in small deformations to a gradient expansion, can be performed for any cases where a perturbative expansion is possible. These include the case of two infinitely thick plates separated by $d$, each modeled as a homogeneous and isotropic dielectric material, with frequency dependent permittivities $\epsilon_{1}(\omega)$ and $\epsilon_{2}(\omega)$. Much like the crossover between retarded and non-retarded van der Waals interactions, the Casimir force also becomes dependent on material properties through $\epsilon(\omega)$ at short distances (typically in the range of 10-100 nm). If the dielectric response is dominated by a resonance at a single characteristic frequency $\overline{\omega}$, in the so-called near-field regime of separations $d\leq c/\overline{\omega}$, the Casimir energy density in Eq. (9) changes from $U(d)\propto\hbar c/d^{3}$ to $U(d)\propto\hbar\overline{\omega}/d^{2}$.

Since $U(d)$ still diverges as $d\to 0$ even in the near field regime, the PFA and gradient expansion should still be applicable. This calculation was implemented in Ref. Bimonte et al. (2012b), which applied this method to two infinitely thick plates composed of homogeneous and isotropic dielectric materials, with general permittivities $\epsilon_{1}(\omega)$ and $\epsilon_{2}(\omega)$. The resulting corrections to PFA are no longer pure numbers, but depend on temperature and separation (through material dependent parameters). While the details will not be reproduced here, we note that they have in fact been used for comparison to sphere/plate experiments Banishev et al. (2013); Bimonte et al. (2021); Bimonte (2018).

In principle, we may anticipate a gradient expansion for any situation when a shape dependent quantity diverges sufficiently rapidly as a function of separation. One example, described in Ref. Krüger et al. (2013), involves radiative heat transfer (RHT) between objects at different temperatures. The heat transfer between parallel plates at large distances is dominated by propagating photons, and is independent of separation $d$. However, as noted by Polder and Van Hove Polder and Van Hove (1971), RHT at short separation increases strongly upon decreasing separation due to tunneling of evanescent EM waves across the vacuum gap. The precise form of this enhancement depends on the material properties: If the response of the material can be characterized by single dominant frequency $\overline{\omega}$, in the near-field regime of $d\leq c/\overline{\omega}$, RHT diverges as $S(d)\propto\hbar\overline{\omega}^{2}/d^{2}$ Golyk et al. (2013). RHT due to evanescent waves has also attracted a lot of interest due to its connection with scanning tunneling microscopy and scanning thermal microscopy under ultra-high vacuum conditions Pendry (1999); Majumdar (1999). The enhancement of HT in the near-field regime (generally denoting separations small compared to the thermal wavelength, which is roughly 8 microns at room temperature) has only recently been verified experimentally Shen et al. (2009); Rousseau et al. (2009).

Much like PFA, the leading term for RHT between closely spaced curved objects has been computed by use of a corresponding proximity transfer approximation (PTA) Krüger et al. (2011); Sasihithlu and Narayanaswamy (2011); Otey and Fan (2011). Ref. Golyk et al. (2013) extended a gradient expansion approach to compute the first correction to PTA, which can be expressed as a spectral decomposition, with a correction factor $\beta(\omega)$ at each frequency.

While the above examples focused on two surfaces, the gradient expansion can also be applied to estimate the interaction between a polarizable particle and a gently curved surface. For a particle of polarizability $\alpha$, the Casimir-Polder force Casimir and Polder (1948) due to a flat surface scales with the separation $d$ as $\alpha\times\hbar c/d^{4}$ in the far-field, and as $\alpha\times\hbar\overline{\omega}/d^{3}$ in the near-field. Perturbative results for the interaction of the particle with a slightly modulated surface can again be converted to a gradient expansion for the force close to a curved surface. In Ref. Bimonte et al. (2014) the leading correction on approaching a surface with profile $H({\bf x})$ is shown to depend on its curvature (proportional to $d\nabla^{2}H$), while the next order scales as $[d\nabla^{2}H]^{2}$. If the polarizable particle is a molecule, Ref. Bimonte et al. (2016) demonstrates that in principle the shift of its rotational levels due to Casimir-Polder interactions can be used as a probe of the surface profile.

Scattering methods are much more tractable for objects with geometries for which the $T$-matrix is diagonal. One can extend the range of such situations by considering more unusual coordinate systems. For scalar fluctuations, these include elliptic cylinder Graham (2013), parabolic cylinder Graham et al. (2010), and spheroidal geometries Emig et al. (2009), for which limiting cases are a finite-width strip, half-plane, and disk respectively. One can also consider scattering in a different variable in order to study a wedge or cone using ordinary cylindrical and spherical coordinates, respectively Maghrebi et al. (2011a). In these cases one uses contour integration to replace the scattering theory sum over partial waves with a continuous integral. One can also consider complementary geometries using systems with planar gaps Kabat et al. (2010) and by applying the generalization of Babinet’s principle Maghrebi et al. (2011b), and for scalar theories one can apply worldline methods Gies et al. (2003).

For perfect conductors, all of these cases can be extended to electromagnetism except for spheroidal scattering, since any translation-invariant geometry can be decomposed into polarizations obeying Dirichlet and Neumann boundary conditions, and scattering from a cone can also be diagonalized through analytic continuation of the standard techniques of Mie scattering. The limiting case of the perfectly conducting disk can nonetheless be solved exactly as an analytic calculation in oblate spheroidal coordinatesMeixner (1948); Emig and Graham (2016), for which the scattering matrix is no longer diagonal.

The wedge, parabolic cylinder, and elliptic cylinder results can be combined to yield a consistent model of the effects of edges Graham et al. (2011); Maghrebi and Graham (2011); Blose et al. (2015). This can be expressed as an approximate form for the Casimir energy per unit length of a perfectly conducting strip parallel to a perfectly conducting plane,

$$\frac{\cal E}{\hbar cL}=-\frac{\pi^{2}}{720}\frac{2d}{H^{3}}+\frac{2\beta}{H^{2}}+\frac{\gamma}{2dH}+\ldots$$

(17)

Here $H$ is the separation between the strip and the plane, $2d$ is the width of the strip, and $\beta=0.00092$ and $\gamma=-0.0040$ are numerical parameters capturing the effect a single edge and the interaction between two edges, respectively; the result for $\beta$ is consistent with limiting cases of the wedge and parabolic cylinder. One can also obtain the energy per unit length of a perfectly conducting half-plane perpendicular to a perfectly conducting plane as Graham et al. (2010)

$$\frac{{\cal E}}{\hbar cL}=\int_{0}^{\infty}\frac{qdq}{4\pi}\log\det\left(\mathbbm{1}_{\nu\nu^{\prime}}-(-1)^{\nu}{\rm k}_{-\nu-\nu^{\prime}-1}(2qH)\right)=-\frac{C_{\perp}}{H^{2}},$$

(18)

where $H$ is the separation distance, ${\rm k}_{\ell}(u)$ is the Bateman k-function Bateman (1931), the determinant is over $\nu,\nu^{\prime}=0,1,2,3\ldots$, and $C_{\perp}=0.0067415$ is obtained by numerical integration.

An interesting extension of Casimir force calculations is to the case of Casimir stresses on a single object. This problem is harder to define than the case of forces between rigid bodies, because deformation of a single object requires a variation of the physical configuration of an object, rather than just a variation of the position of fixed objects. As a result, the formally infinite contributions localized on the object no longer cancel, and we must apply renormalization techniques, which in turn require a model of the object’s material properties.

Debates about Casimir stresses originate in an early model of the electron as a charged shell for which attractive Casimir forces balance Coulomb repulsion Casimir (1953). Here the dynamics of the shell are assumed as a definition of the fundamental interactions. One finds that specifying an ideal conductor boundary condition leads to a positive energy and hence expansion of the shell Boyer (1968); Milton et al. (1978); Milton (2001), invalidating the model’s original premise. Of course, modern quantum electrodynamics provides a well-verified description of the electron within the broader framework of renormalized quantum field theory.

This result has often been misinterpreted, however, to imply that a conducting shell will experience a repulsive stress. This case is very different, because the dynamics of an actual piece of conducting material are determined by the underlying physics of the matter that makes it up, and so cannot be assumed to implement an ideal boundary condition at all frequencies. (The Casimir energy of an idealized boundary can be of interest for the mathematical “Weyl problem” of the relationship between eigenvalue spectra and geometry Abalo et al. (2012); Kolomeisky et al. (2010).) A repulsive stress is also in conflict with the well-established attractive Casimir force between conductors; in particular, the force between two hemispheres, as in any configuration of conductors with mirror symmetry, is always attractive Kenneth and Klich (2006). Similarly, while a boundary condition calculation yields a repulsive stress on a rectangular solid, the force on a piston is attractive, as shown in two dimensions in Ref. Cavalcanti (2004) and three dimensions in Ref. Hertzberg et al. (2005), where the latter is derived using a geometrical optics approximation Jaffe and Scardicchio (2004).

The resolution of this puzzle lies in contributions from within the material, which no longer cancel when the body is not rigid. Modeling the material as a position-dependent dielectric Graham et al. (2013) shows that the contribution from within the material is attractive, and larger in magnitude than the repulsive external contribution. In the limit where the dielectric becomes infinitely strong, the sphere becomes a conductor Bordag and Khusnutdinov (2008), but this limit does not commute with the limits involved in renormalization: the theory must be cut off at a high enough frequency that the material is transparent. Qualitatively similar results have been obtained by considering specific materials rather than a generic dielectric, such as a carbon nanostructure Barton (2001) or a “fish-eye” medium Leonhardt and Simpson (2011). These findings are also in agreement with the work of Deutsch and Candelas Deutsch and Candelas (1979); Candelas (1982), who showed that divergences in Casimir stresses arise from surface counterterms Symanzik (1981) that cannot be removed by renormalization, as well as more detailed calculations in models with a fluctuating scalar field Graham et al. (2003, 2002, 2004).

A key aspect of these works is that divergences associated with an ideal boundary arise in two different ways: the “sharp” limit of infinitesimal thickness, and the “strong” limit of strong potential. For stresses on a real material, these limits are intertwined: for example, as a spherical shell expands and increases its surface area, it becomes thinner, leading to a weaker electromagnetic response. A precise description of this situation requires detailed information about the the material’s electromagnetic properties, but one can gain a qualitative understanding by approximating this behavior via renormalization counterterms. In this process, one introduces a short-distance cutoff, which in a material represents a length scale set by the lattice spacing. The formally divergent contributions arising for a conducting material at short wavelength are then combined with cutoff-dependent counterterms, which are fixed so that the combined result matches experimental inputs. For a stress calculation, one must then fix the dependence of these inputs on the geometry. A reasonable estimate is to assume that the overall strength of the potential, integrated over volume, remains constant, roughly corresponding to holding the total number of charge carriers constant. In scalar models Graham et al. (2003, 2002, 2004) one can also adjust the strength of the potential used to model the object’s scattering response, while for electromagnetic models it is more appropriate to use a Drude model dielectric with fixed plasma wavelength and conductivity, and instead vary the thickness of the material so that its volume remains constant.

While this result does not precisely model a particular material like the examples listed above, it allows for a generic description of the discrepancy between the apparent repulsive stress on a spherical conducting shell with the generally attractive Casimir force for dieletric materials, including their perfectly conducting limit. What one finds Graham et al. (2013) is that the Casimir self-energy decomposes into a sum of two terms: The “traditional” contribution due to the shell’s external electromagnetic response and an “additional” contribution localized within the thickness of the shell, which consists of a renormalized integral over $r$. The former reproduces the Boyer result in the conducting limit, but the latter is of comparable magnitude with opposite sign, and renders the full result attractive, as shown in Fig. 1.

The common starting point of the bulk and of the surface approaches to computing Casimir forces between two or more bodies is the physical picture that the vacuum surrounding the bodies is filled with quantum and thermal fluctuations of the EM field. If the bodies are in thermal equilibrium at temperature $T$ with the environment, the correlators of the EM field can be derived from linear-response theory Agarwal (1975):

	$\displaystyle\langle{\hat{E}}_{i}({\bf x},t){\hat{E}}_{j}({\bf x^{\prime}},t^{\prime})\rangle_{\rm sym}$	$\displaystyle=$	$\displaystyle\hbar\int_{-\infty}^{\infty}\!\!\frac{d\omega}{2\pi}\coth\left(\frac{\hbar\omega}{2k_{B}T}\right){\rm Im}[{\cal G}^{(EE)}_{ij}({\bf x},{\bf x^{\prime}},\omega)]\;e^{-i\omega(t-t^{\prime})}\;,$
	$\displaystyle\langle{\hat{H}}_{i}({\bf x},t){\hat{H}}_{j}({\bf x^{\prime}},t^{\prime})\rangle_{\rm sym}$	$\displaystyle=$	$\displaystyle\hbar\int_{-\infty}^{\infty}\!\!\frac{d\omega}{2\pi}\coth\left(\frac{\hbar\omega}{2k_{B}T}\right){\rm Im}[{\cal G}^{(HH)}_{ij}({\bf x},{\bf x^{\prime}},\omega)]\;e^{-i\omega(t-t^{\prime})}\;,$
	$\displaystyle\langle{\hat{E}}_{i}({\bf x},t){\hat{H}}_{j}({\bf x^{\prime}},t^{\prime})\rangle_{\rm sym}$	$\displaystyle=$	$\displaystyle\hbar\int_{-\infty}^{\infty}\!\!\frac{d\omega}{2\pi}\coth\left(\frac{\hbar\omega}{2k_{B}T}\right)\left\{-{\rm i}\,{\rm Re}[{\cal G}^{(EH)}_{ij}({\bf x},{\bf x^{\prime}},\omega)]\right\}e^{-i\omega(t-t^{\prime})}\;,$		(19)

where the subscript sym on the average symbols denotes the symmetrized products of field operators, and the ${\cal G}_{ij}^{(\alpha\beta)}({\bf r},{\bf r^{\prime}},\omega)$ with $\alpha,\beta=E,H$ are the Fourier transforms of the classical Green’s functions for the system of bodies. The latter Green’s functions are obtained by solving the macroscopic Maxwell equations, subjected to the appropriate boundary conditions on the surfaces of the bodies. We underline that Eqs. (19) include both zero-point (i.e. quantum) and thermal fluctuation of the em field. The validity of Eqs. (19) is subjected to the condition that the length scales of the relevant fluctuations should be large compared to atomic distances, in such a way that the EM material properties of the bodies can be described by their macroscopic response functions. We shall assume for simplicity that the bodies are made out of homogeneous and isotropic magneto-dielectric materials, characterized by the respective frequency-dependent electric and magnetic permittivities $\epsilon(\omega)$ and $\mu(\omega)$.

The mechanical effects of the fluctuating EM field on the bodies are determined by the expectation value $\langle T_{ij}\rangle$ of the Maxwell stress tensor

$$\langle T_{ij}({\bf x})\rangle=\frac{1}{4\pi}\left\{\langle E_{i}({\bf x})E_{j}({\bf x})\rangle+\langle H_{i}({\bf x})H_{j}({\bf x})\rangle-\frac{1}{2}\delta_{ij}\left[\langle E_{k}({\bf x})E_{k}({\bf x})\rangle+\langle H_{k}({\bf x})H_{k}({\bf x})\rangle\right]\right\}\;,$$

(20)

where the correlators are evaluated for equal times $t=t^{\prime}$. To gain further insight, it is convenient to note that at points ${\bf x},{\bf x}^{\prime}$ both lying in the vacuum medium, the Green’s functions ${\cal G}^{(\alpha\beta)}_{ij}({\bf x},{\bf x}^{\prime};\omega)$ can be decomposed as follows:

$${\cal G}^{(\alpha\beta)}_{ij}({\bf x},{\bf x}^{\prime};\omega)={\cal G}^{(\alpha\beta;0)}_{ij}({\bf x}-{\bf x}^{\prime};\omega)+{\Gamma}^{(\alpha\beta)}_{ij}({\bf x},{\bf x}^{\prime};\omega)\;,$$

(21)

where ${\cal G}^{(\alpha\beta;0)}_{ij}({\bf x}-{\bf x}^{\prime};\omega)$ is the Green’s function of free space, while ${\Gamma}^{(\alpha\beta)}_{ij}({\bf x},{\bf x}^{\prime};\omega)$ describes the em field scattered by the bodies. When the above decomposition is used in Eq. (20), one finds that the expectation value of the stress tensor is accordingly decomposed as

$$\langle T_{ij}({\bf x})\rangle=\langle T^{(0)}_{ij}({\bf x})\rangle+{\Theta}_{ij}({\bf x})\;,$$

(22)

where $\langle T^{(0)}_{ij}({\bf x})\rangle$ is the free-space contribution, while ${\Theta}_{ij}({\bf x})$ is the scattering contribution. According to Eqs. (19), $\langle T^{(0)}_{ij}({\bf x})\rangle$ and ${\Theta}_{ij}({\bf x})$ involve frequency integrals of the imaginary parts of the respective Green’s functions ${\cal G}^{(\alpha\beta;0)}_{ij}({\bf x},{\bf x}^{\prime};\omega)$ and ${\Gamma}^{(\alpha\beta)}_{ij}({\bf x},{\bf x}^{\prime};\omega)$. Let us consider $\langle T^{(0)}_{ij}({\bf x})\rangle$ first. Since

$$\lim_{{\bf x}\rightarrow{\bf x}^{\prime}}{\rm Im}[{\bf\cal G}^{(EE;0)}_{ij}({\bf x},{\bf x}^{\prime},\omega)]=\lim_{{\bf x}\rightarrow{\bf x}^{\prime}}{\rm Im}[{\bf\cal G}^{(HH;0)}_{ij}({\bf x},{\bf x}^{\prime},\omega)]=\frac{2\,\omega^{3}}{3\,c^{3}}\delta_{ij}\;.$$

(23)

it is clear that $\langle T^{(0)}_{ij}({\bf x})\rangle$ is expressed by a formally divergent frequency integral. One notes, however, that according to Eq. (23), $\langle T^{(0)}_{ij}({\bf x})\rangle$ represents a homogeneous and isotropic pressure acting the surfaces of the bodies, which is independent of the material properties of the bodies. Since this pressure cannot give rise to an overall force on the body, just a stress as discussed above, we can neglect it here. Let us turn now to the scattering contribution ${\Theta}_{ij}({\bf x})$. By performing a Wick rotation to the imaginary frequency axis, one finds that ${\Theta}_{ij}({\bf x})$ has the expression

	$\displaystyle{\Theta}_{ij}({\bf x})$	$\displaystyle=$	$\displaystyle\frac{k_{B}T}{2\pi}\left.\sum_{n=0}^{\infty}\right.\!^{\prime}\left[{\Gamma}^{(EE)}_{ij}({\bf x},{\bf x};{\rm i}\,\xi_{n})+{\Gamma}^{(HH)}_{ij}({\bf x},{\bf x};{\rm i}\,\xi_{n})\right.$		(25)
			$\displaystyle\left.-\frac{1}{2}\delta_{ij}\left({\Gamma}^{(EE)}_{kk}({\bf x},{\bf x};{\rm i}\,\xi_{n})+{\Gamma}^{(HH)}_{kk}({\bf x},{\bf x};{\rm i}\,\xi_{n})\right)\right]\;.$

where $\xi_{n}=2\pi nk_{B}T/\hbar$ are the Matsubara imaginary frequencies, and the prime in the summations means that the $n=0$ term is taken with a weight of one-half. The sum on the right-hand side of the above equation is finite, because the scattering parts of the Green’s functions ${\Gamma}^{(\alpha\beta)}_{ij}({\bf x},{\bf x}^{\prime};\omega)$ remain finite in the coincidence limit ${\bf x}^{\prime}\rightarrow{\bf x}$.

The Casimir force on body $r$ can be now obtained by integrating ${\Theta}_{ij}({\bf x})$ on any surface $S_{r}$ drawn in the vacuum that surrounds that body and excludes all other bodies, yielding

$$F_{i}^{(r)}=\oint_{{S}_{r}}d^{2}\sigma\,\hat{n}_{j}({\bf x})\;{\Theta}_{ji}({\bf x})\;.$$

(26)

The surface integral on the right-hand side of this equation can be recast in a remarkably simple form, by taking advantage of the peculiar structure of the scattering Green’s functions ${\Gamma}^{(\alpha\beta)}_{ij}({\bf x},{\bf x}^{\prime})$ (for brevity, from now on we shall not display the dependence of the Green’s functions on the Matsubara frequencies $\xi_{n}$). This structure becomes manifest once ${\Gamma}^{(\alpha\beta)}_{ij}({\bf x},{\bf x}^{\prime})$ is expressed in terms of the so-called $T$-operator Waterman (1965, 1971). We recall that the $T$-operator $T_{kl}^{(\rho\sigma)}({\bf y},{\bf y}^{\prime})$ is defined to give the currents induced at the point ${\bf y}$ in the interior of the bodies by an external em field at point ${\bf y}^{\prime}$, Using the $T$-operator, the scattering parts of the Green’s functions can be expressed as

$${\Gamma}^{(\alpha\beta)}_{ij}({\bf x},{\bf x}^{\prime})=\sum_{r,r^{\prime}=1}^{N}\int_{V_{r}}d^{3}{\bf y}\int_{V_{r^{\prime}}}d^{3}{\bf y}^{\prime}{\cal G}^{(\alpha\rho;0)}_{ik}({\bf x}-{\bf y})T_{kl}^{(\rho\sigma)}({\bf y},{\bf y}^{\prime}){\cal G}_{lj}^{(\sigma\beta;0)}({\bf y}^{\prime}-{\bf x}^{\prime})\;.$$

(27)

Despite its simple meaning, the $T$-operator is difficult to compute in general, and its expression is known only for bodies with sufficiently symmetric geometries, mostly in cases where the $T$-operator is diagonal. A more powerful representation of the scattering Green’s functions can be derived from the equivalence principle of classical electromagnetic theory Harrington (2001). According to this principle, the induced currents existing in the interior of the bodies can be replaced by fictitious surface currents, which produce the same scattered field as the induced currents. By following this equivalence principle, one arrives at an alternative representation of the scattering Green’s functions consisting of a double surface integral extended onto the surfaces of the bodies Reid et al. (2013); Bimonte and Emig (2021):

	$\displaystyle{\Gamma}^{(\alpha\beta)}_{ij}({\bf x},{\bf x}^{\prime})$	$\displaystyle=$	$\displaystyle-\sum_{r,r^{\prime}=1}^{N}\int_{V_{r}}d^{3}{\bf y}\int_{V_{r^{\prime}}}d^{3}{\bf y}^{\prime}\,\delta(F_{r}({\bf y}))\,\delta(F_{r^{\prime}}({\bf y}^{\prime}))\;$		(28)
		$\displaystyle\times$	$\displaystyle{\cal G}^{(\alpha\rho;0)}_{ik}({\bf x}-{\bf y}))\left(M^{-1}\right)_{kl}^{(\rho\sigma)}({\bf y},{\bf y}^{\prime}){\cal G}_{lj}^{(\sigma\beta;0)}({\bf y}^{\prime}-{\bf x}^{\prime})\;,$

where $\delta(x)$ is the Dirac delta function and $F_{r}({\bf y})=0$ is the equation of the surface $\Sigma_{r}$ of the $r$-th body. The surface operator $M_{ij}^{(\alpha\beta)}({\bf y},{\bf y}^{\prime})$, whose inverse appears in the integral on the right-hand side, has a remarkably simple expression Reid et al. (2013); Bimonte and Emig (2021)

$$M_{ij}^{(\alpha\beta)}({\bf y},{\bf y}^{\prime})=\left\{\begin{array}[]{ll}\!\Pi^{(r)}_{ik}({\bf y})\left[{\cal G}^{(\alpha\beta;r)}_{kl}({\bf y}-{\bf y}^{\prime})+{\cal G}^{(\alpha\beta;0)}_{kl}({\bf y}-{\bf y}^{\prime})\right]\Pi^{(r)}_{lj}({\bf y}^{\prime})&{\rm if}\;\;{\bf y},{\bf y}^{\prime}\in\Sigma_{r}\\ \Pi_{ik}^{(r)}({\bf y})\;{\cal G}^{(\alpha\beta;0)}_{kl}({\bf y}-{\bf y}^{\prime})\;\Pi_{lj}^{(s)}({\bf y}^{\prime})&{\rm if}\;\,{\bf y}\in\Sigma_{r},{\bf y}^{\prime}\in\Sigma_{s},\\ &\;\;\;\;r\neq s\end{array}\right.\,.$$

(29)

where $\Pi^{(r)}_{ik}({\bf y})$ is the projector onto the plane tangent at $\Sigma_{r}$ at ${\bf y}$, and ${\cal G}^{(\alpha\beta;r)}_{ij}({\bf y}-{\bf y}^{\prime})$ are the Green’s functions for an infinite homogeneous and isotropic magneto-dielectric medium identical to that of body $r$. By inspection of Eqs. 27 and 28, it is apparent that both representations of the scattering Green’s functions have a common ${\hat{G}}^{(0)}{\hat{\cal K}}\,{\hat{G}}^{(0)}$ structure, where the ${\hat{\cal K}}$ operator coincides either with Waterman’s $T$-operator, or with (minus) the inverse of the surface operator ${\hat{M}}$. Importantly, both the $T$-operator and the surface operator ${\hat{M}}$ satisfy the reciprocity relations Harrington (2001) enjoyed by the Green’s functions. In Ref. Bimonte and Emig (2021) it is shown that whenever the scattering Green’s function has the ${\hat{G}}^{(0)}{\hat{\cal K}}\,{\hat{G}}^{(0)}$ structure, with a ${\hat{\cal K}}$ satisfying reciprocity, the force formula in Eq. (26) can be recast in the remarkably simple form

$${\bf F}^{(r)}=k_{B}T\left.\sum_{n=0}^{\infty}\right.\!\!^{\prime}\,{\rm Tr}\left[{\hat{\cal K}}({\rm i}\,\xi_{n})\frac{\partial}{\partial{\bf x}_{r}}\hat{\cal G}^{(0)}({\rm i}\,\xi_{n})\right]\;,$$

(30)

where the trace operation ${\rm Tr}$ denotes an integral over the volumes occupied by the bodies and a sum over both the spatial indices $i,j$ and internal indices $\alpha,\beta$, while the symbol ${\partial}/{\partial{\bf x}_{r}}$ denotes a derivative with respect to a rigid translation of the $r$-th body. When ${\hat{\cal K}}$ is identified with the $T$-operator, Eq. (30) reproduces the “TGTG” formula derived in Ref. Kenneth and Klich (2006). If ${\hat{\cal K}}$ is instead identified with (minus) the surface operator $(\hat{M})^{-1}$, then Eq. (30) reproduces the surface force formula derived in Ref. Reid et al. (2013). We thus see that Eq. (30) encompasses in a single compact formula the bulk and the surface formulations of the Casimir force.

Starting from Eq. (30), it is possible to compute the Casimir free-energy of the system of bodies. Direct integration of Eq. (30) leads to a formally divergent result ${\cal F}_{\rm bare}$. However, a finite expression is easily recovered by subtracting from ${\cal F}_{\rm bare}$ the self-energies of the individual bodies. The details of the subtraction procedure differ slightly within the bulk and surface formulations, but the final formula for the renormalized Casimir free energy ${\cal F}$ has an identical form in both approaches Bimonte and Emig (2021). In the simple case of two bodies, it reads

$${\cal F}=k_{B}T\left.\sum_{n=0}^{\infty}\right.\!^{\prime}{\rm Tr}\,\log\left[1-\hat{\cal{K}}_{1}\,\hat{\cal G}^{(0)}\hat{\cal{K}}_{2}\hat{\cal G}^{(0)}\right]\;.$$

(31)

Within the bulk approach, the operator $\hat{{\cal K}}_{r}$ coincides with the $T$-operator $T_{ij}^{(\alpha\beta;r)}({\bf y},{\bf y}^{\prime})$ of body $r$ in isolation, while in the surface approach $\hat{{\cal K}}_{r}$ coincides with (minus) the inverse of the surface operator $M_{ij}^{(\alpha\beta;r)}({\bf y},{\bf y}^{\prime})$ defined by the first line of Eq. 29.

We underline that Eqs. (30) and (31) are valid for any shapes and relative dispositions of the bodies. In particular, both Equations hold for a system of interleaved bodies, possibly nested one inside the other. A case deserving special consideration is that of a “separable” configuration of two bodies, namely a configuration in which the bodies can be separated from one another by a plane drawn in the vacuum between them. Three typical examples of separable configurations are those of a sphere opposite a plate, a system of two spheres, or a sphere and a cylinder, which represent the configurations adopted in the vast majority of experiments. In Bimonte and Emig (2021) it is shown that in a separable configuration, the general formula Eq. 31 reproduces the famous scattering formula Emig et al. (2007); Kenneth and Klich (2008); Lambrecht et al. (2006); Rahi et al. (2009) in which the operators $\hat{{\cal K}}_{r}$ are replaced by the the scattering matrices of the bodies and $\hat{\cal G}^{(0)}$ is replaced by the so-called translation matrices.

As in the previous subsection, we consider again $N$ dielectric bodies with properties described above. In the Euclidean path integral quantization of the electromagnetic field, the Casimir free energy at finite temperature $T$ can be obtained as

$${\mathcal{F}}=-k_{B}T\left.\sum_{n=0}^{\infty}\right.\!^{\prime}\log\frac{\mathcal{Z}(\xi_{n})}{\mathcal{Z}_{\infty}(\xi_{n})}\,,$$

(32)

where again the sum runs over the Matsubara frequencies $\xi_{n}=2\pi nk_{B}T/\hbar$, with a weight of $1/2$ for $n=0$. The partition function $\mathcal{Z}$ is given by a path integral that we shall derive now. The partition function $\mathcal{Z}_{\infty}$ describes the configuration of infinitely separated bodies and subtracts the self-energies of the bodies from the bare free energy.

We express the action of the electromagnetic field in the absence of free sources in terms of the gauge field ${\bf A}$. We chose the transverse gauge with $A_{0}=0$. The functional integral will then run over ${\bf A}$ only. The electric field is given by ${\bf E}=ik{\bf A}\to-\kappa{\bf A}$ and the magnetic field by ${\bf B}=\nabla\times{\bf A}$. Then the action in terms of the induced sources at fixed frequency $\kappa$ is given by

$\displaystyle\hat{S}[{\bf A}]$

$\displaystyle=$

$\displaystyle-\frac{1}{2}\int_{\mathbb{R}^{3}}d^{3}{\bf x}\,\left[{\bf A}^{2}\epsilon_{\bf x}\kappa^{2}+\frac{1}{\mu_{\bf x}}(\nabla\times{\bf A})^{2}\right]-\kappa\sum_{r=1}^{N}\int_{V_{r}}d^{3}{\bf x}\,{\bf A}\cdot{\bf P}_{r}\,,$

(33)

for fluctuations ${\bf A}$ of the gauge field, and induced bulk currents ${\bf P}_{r}$ inside the objects. The inverse of the kernel of the quadratic part of this action is given by the Green’s tensor which for spatially constant $\epsilon_{r}$ and $\mu_{r}$ of body $r$ is given by

$${\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(AA;r)}({\bf x},{\bf x}^{\prime})=\mu_{r}\left({\bf 1}-\frac{1}{\epsilon_{r}\mu_{r}\kappa^{2}}\nabla\otimes\nabla\right)\frac{e^{-\sqrt{\epsilon_{r}\mu_{r}}\kappa|{\bf x}-{\bf x}^{\prime}|}}{|{\bf x}-{\bf x}^{\prime}|}\,,$$

(34)

From the relation between the gauge field ${\bf A}$ and the electric field ${\bf E}$ the relation $-\kappa^{2}{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(AA;r)}({\bf x},{\bf x}^{\prime})={\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(EE;r)}({\bf x},{\bf x}^{\prime})$ follows, which allows us to compare the results below to those of the stress-tensor-based derivation. Next, we define the classical solutions $\bm{\mathcal{A}}_{r}$ of the vector wave equation in each region $V_{r}$, obeying $\nabla\times\nabla\times\bm{\mathcal{A}}_{r}+\epsilon_{r}\mu_{r}\kappa^{2}\bm{\mathcal{A}}_{r}=-\kappa\mu_{r}{\bf P}_{r}$. Then the source terms of Eq. (33) can be written as an integral over the surface of body,

$$-\kappa\int_{V_{r}}d^{3}{\bf x}{\bf A}\cdot{\bf P}_{r}=\frac{1}{\mu_{r}}\int_{\Sigma_{r}}d^{3}{\bf x}\left[{\bf A}_{-}\cdot({\bf n}_{r}\times(\nabla\times\bm{\mathcal{A}}_{r}))+(\nabla\times{\bf A})_{-}\cdot({\bf n}_{r}\times\bm{\mathcal{A}}_{r})\right]\,.$$

(35)

The values of the gauge field ${\bf A}$ and its curl $\nabla\times{\bf A}$ appearing in this expression are those when the surface is approached from the inside, denoted by ${\bf A}_{-}$ and $(\nabla\times{\bf A})_{-}$. It is important to note that in the above surface integral, ${\bf A}$ and $\nabla\times{\bf A}$ multiply vectors that are tangential to the surface, and hence only the tangential components of ${\bf A}$ and $\nabla\times{\bf A}$ contribute to the integral. Hence, we can use the continuity conditions of the tangential components of ${\bf E}$ and ${\bf H}$,

$${\bf n}_{r}\times{\bf E}_{-}={\bf n}_{r}\times{\bf E}_{+}\,,\quad\frac{1}{\mu_{r}}{\bf n}_{r}\times(\nabla\times{\bf E})_{-}=\frac{1}{\mu_{0}}{\bf n}_{r}\times(\nabla\times{\bf E})_{+}\,,$$

(36)

to write the source terms also as

$$-\kappa\int_{V_{r}}d^{3}{\bf x}{\bf A}\cdot{\bf P}_{r}=\int_{\Sigma_{r}}d^{3}{\bf x}\left[\frac{1}{\mu_{r}}{\bf A}_{+}\cdot({\bf n}_{r}\times(\nabla\times\bm{\mathcal{A}}_{r}))+\frac{1}{\mu_{0}}(\nabla\times{\bf A})_{+}\cdot({\bf n}_{r}\times\bm{\mathcal{A}}_{r})\right]\,,$$

(37)

where the values of the gauge field ${\bf A}$ and its curl $\nabla\times{\bf A}$ appearing in this expression are now those when the surface is approached from the outside, denoted by ${\bf A}_{+}$ and $(\nabla\times{\bf A})_{+}$.

Now we shall see the advantage of having expressed the source integrals in terms of the values of ${\bf A}$ and $\nabla\times{\bf A}$ when the surfaces are approached from either the outside or the inside of the objects. In the region $V_{0}$, the field ${\bf A}\equiv{\bf A}_{0}$ is fully determined by its values on the surfaces $\Sigma_{r}$ and $\epsilon_{0}$, $\mu_{0}$, which are constant across $V_{0}$. When integrating out ${\bf A}_{0}$, one computes the two-point correlation function of ${\bf A}_{+}$ and $(\nabla\times{\bf A})_{+}$ on the surfaces $\Sigma_{r}$, and hence the behavior of ${\bf A}_{0}$ inside the regions $V_{r}$ with $r>0$ is irrelevant. Following the same arguments for ${\bf A}\equiv{\bf A}_{r}$ inside the objects, the behavior of ${\bf A}_{r}$ outside of region $V_{r}$ is irrelevant for computing the correlations of ${\bf A}_{-}$ and $(\nabla\times{\bf A})_{-}$ on the surfaces $\Sigma_{r}$. Hence, we can replace in the action the spatially dependent $\epsilon_{\bf x}$ by $\epsilon_{0}$ when the coupling of ${\bf A}_{0}$ to the surface fields $\bm{\mathcal{A}}_{r}$ is represented by Eq. (37), and similarly replace $\epsilon_{\bf x}$ by $\epsilon_{r}$ when the coupling of ${\bf A}_{r}$ to the surface fields $\bm{\mathcal{A}}_{r}$ is represented by Eq. (35).

That this is justified can also be understood as follows. The field ${\bf A}_{0}$ in region $V_{0}$ can be expanded in a basis of functions that obey the wave equation with $\epsilon_{0}$. The same can be done for ${\bf A}_{r}$ in the interior of each object, i.e., ${\bf A}_{r}$ can be expanded in a basis of functions that obey the wave equation with $\epsilon_{r}$ in $V_{r}$. For each given set of expansion coefficients in $V_{0}$ there are corresponding coefficients within each region $V_{r}$ that are determined by the continuity conditions at the surfaces $\Sigma_{r}$. The functional integral over ${\bf A}$ then corresponds to integrating over consistent sets of expansion coefficients that are related by the continuity conditions. The two-point correlations of ${\bf A}_{+}$ and $(\nabla\times{\bf A})_{+}$ on the surfaces $\Sigma_{r}$ are then fully determined by the integral over the expansion coefficients of ${\bf A}_{0}$ in $V_{0}$ only, and the interior expansion coefficients play no role. Equivalently, the two-point correlations of ${\bf A}_{-}$ and $(\nabla\times{\bf A})_{-}$ on the surfaces $\Sigma_{r}$ are then fully determined by the integral over the expansion coefficients of ${\bf A}_{r}$ in $V_{r}$ only, and now the exterior expansion coefficients are irrelevant. Hence, in the functional integral, the integration of ${\bf A}$ can be replaced by $N+1$ integrations over the fields ${\bf A}_{r}$, $r=0,\ldots,N$, where each ${\bf A}_{r}$ is allowed to extend over unbounded space with the action for a free field in a homogeneous space with $\epsilon_{r}$, $\mu_{r}$. The multiple counting of degrees of freedom that results from $N+1$ functional integrations poses no problem since the (formally infinite) factor in the partition function cancels when the Casimir energy is computed from Eq. (32).

With this representation, we can write the partition function as a functional integral over the fluctuations ${\bf A}_{r}$, separately in each region $V_{r}$, and the surface fields $\bm{\mathcal{A}}_{r}$ on body $r$, with the action

	$\displaystyle\hat{S}[\{{\bf A}_{r}\},\{\bm{\mathcal{A}}_{r}\}]$	$\displaystyle=$	$\displaystyle-\frac{1}{2}\sum_{r=0}^{N}\int_{\mathbb{R}^{3}}d^{3}{\bf x}\left[{\bf A}_{r}^{2}\epsilon_{r}\kappa^{2}+\frac{1}{\mu_{r}}(\nabla\times{\bf A}_{r})^{2}\right]$
		$\displaystyle+$	$\displaystyle\sum_{r=1}^{N}\int_{\Sigma_{r}}d^{3}{\bf x}\,\left[\frac{1}{\mu_{r}}{\bf A}_{0}({\bf n}_{r}\times(\nabla\times\bm{\mathcal{A}}_{r}))+\frac{1}{\mu_{0}}(\nabla\times{\bf A}_{0})({\bf n}_{r}\times\bm{\mathcal{A}}_{r})\right]$
		$\displaystyle+$	$\displaystyle\sum_{r=1}^{N}\int_{\Sigma_{r}}d^{3}{\bf x}\,\left[\frac{1}{\mu_{r}}{\bf A}_{r}({\bf n}_{r}\times(\nabla\times\bm{\mathcal{A}}_{r}))+\frac{1}{\mu_{r}}(\nabla\times{\bf A}_{r})({\bf n}_{r}\times\bm{\mathcal{A}}_{r})\right],.$

Now, the fluctuations ${\bf A}_{r}$ can be integrated out easily, noting that the two-point correlation function $\langle{\bf A}_{r}({\bf x}){\bf A}_{r^{\prime}}({\bf x}^{\prime})\rangle=0$ for all $r$, ${r^{\prime}}=0,\ldots,N$ with $r\neq{r^{\prime}}$. This yields the partition function

	$\displaystyle\mathcal{Z}(\xi)$	$\displaystyle=$	$\displaystyle\prod_{r=1}^{N}\int{\mathcal{D}}\bm{\mathcal{A}}_{r}\exp\left[-\frac{\beta}{2}\left(\sum_{r=1}^{N}\int_{\Sigma_{r}}\!\!\!d^{3}{\bf x}\int_{\Sigma_{r}}\!\!\!d^{3}{\bf x}^{\prime}\bm{\mathcal{A}}_{r}({\bf x})L_{r}({\bf x},{\bf x}^{\prime})\bm{\mathcal{A}}_{r}({\bf x}^{\prime})\right.\right.$		(39)
		$\displaystyle+$	$\displaystyle\left.\left.\sum_{r,{r^{\prime}}=1}^{N}\int_{\Sigma_{r}}\!\!\!d^{3}{\bf x}\int_{\Sigma_{{r^{\prime}}}}\!\!\!d^{3}{\bf x}^{\prime}\bm{\mathcal{A}}_{r}({\bf x})M_{r{r^{\prime}}}({\bf x},{\bf x}^{\prime})\bm{\mathcal{A}}_{r^{\prime}}({\bf x}^{\prime})\right)\right]$

with the kernels

	$\displaystyle L_{r}({\bf x},{\bf x}^{\prime})=$	$\displaystyle\frac{1}{\mu_{r}^{2}}\left[\,\nabla\times\nabla\times{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(AA;r)}({\bf x},{\bf x}^{\prime})({\bf n}_{r}\times\reflectbox{$\vec{\reflectbox{$\cdot$}}$}\,)({\bf n}^{\prime}_{r}\times\vec{\cdot}\,)\right.$
		$\displaystyle+\left.\nabla\times{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(AA;r)}({\bf x},{\bf x}^{\prime})({\bf n}_{r}\times(\nabla\times\reflectbox{$\vec{\reflectbox{$\cdot$}}$}\,))({\bf n}^{\prime}_{r}\times\vec{\cdot}\,)\right.$
		$\displaystyle+\left.\nabla\times{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(AA;r)}({\bf x},{\bf x}^{\prime})({\bf n}_{r}\times\reflectbox{$\vec{\reflectbox{$\cdot$}}$}\,)({\bf n}^{\prime}_{r}\times(\nabla^{\prime}\times\vec{\cdot}\,))\right.$
		$\displaystyle+\left.{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(AA;r)}({\bf x},{\bf x}^{\prime})({\bf n}_{r}\times(\nabla\times\reflectbox{$\vec{\reflectbox{$\cdot$}}$}\,))({\bf n}^{\prime}_{r}\times(\nabla^{\prime}\times\vec{\cdot}\,))\right]$
	$\displaystyle M_{r{r^{\prime}}}({\bf x},{\bf x}^{\prime})=$	$\displaystyle\frac{1}{\mu_{0}^{2}}\,\nabla\times\nabla\times{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(AA;0)}({\bf x},{\bf x}^{\prime})({\bf n}_{r}\times\reflectbox{$\vec{\reflectbox{$\cdot$}}$}\,)({\bf n}^{\prime}_{r^{\prime}}\times\vec{\cdot}\,)$
		$\displaystyle+\frac{1}{\mu_{0}\mu_{r}}\,\nabla\times{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(AA;0)}({\bf x},{\bf x}^{\prime})({\bf n}_{r}\times(\nabla\times\reflectbox{$\vec{\reflectbox{$\cdot$}}$}\,))({\bf n}^{\prime}_{r^{\prime}}\times\vec{\cdot}\,)$
		$\displaystyle+\frac{1}{\mu_{0}\mu_{r^{\prime}}}\,\nabla\times{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(AA;0)}({\bf x},{\bf x}^{\prime})({\bf n}_{r}\times\reflectbox{$\vec{\reflectbox{$\cdot$}}$}\,)({\bf n}^{\prime}_{r^{\prime}}\times(\nabla^{\prime}\times\vec{\cdot}\,))$
		$\displaystyle+\frac{1}{\mu_{r}\mu_{r^{\prime}}}\,{\stackrel{{\scriptstyle\leftrightarrow}}{{\cal G}}}^{(AA;0)}({\bf x},{\bf x}^{\prime})({\bf n}_{r}\times(\nabla\times\reflectbox{$\vec{\reflectbox{$\cdot$}}$}\,))({\bf n}^{\prime}_{r^{\prime}}\times(\nabla^{\prime}\times\vec{\cdot}\,))\,,$		(40)

where the arrow over the placeholder $\cdot$ indicates to which side of the kernel it acts. This notation implies that the derivatives are taken before the kernel is evaluated with ${\bf x}$ and ${\bf x}^{\prime}$ on the surfaces $\Sigma_{r}$, i.e., information about the behavior of basis functions is required in an infinitesimal vicinity of the surfaces. There is an important simplification of this representation: The bilinear form described by the kernel $L_{r}$ is degenerate on the space of functions over which the functional integral runs, i.e., $\int_{\Sigma_{r}}\!\!\!d^{3}{\bf x}\int_{\Sigma_{r}}\!\!\!d^{3}{\bf x}^{\prime}\bm{\mathcal{A}}_{r}({\bf x})L_{r}({\bf x},{\bf x}^{\prime})\bm{\mathcal{A}}_{r}({\bf x}^{\prime})=0$ for all $\bm{\mathcal{A}}_{r}({\bf x})$ that are regular solutions of the vector wave equation $\nabla\times\nabla\times\bm{\mathcal{A}}_{r}+\epsilon_{r}\mu_{r}\kappa^{2}\bm{\mathcal{A}}_{r}=0$ inside region $V_{r}$. This implies that the kernel $L_{r}$ can be ignored in the above functional integral over regular waves $\bm{\mathcal{A}}_{r}$ inside the objects, and the partition function in Eq. 32 is given by the functional determinant of the kernel $M_{rr^{\prime}}$ alone, which can be computed in some basis for regular waves inside the bodies, evaluated at the surfaces only.