Macroscopic quantum electrodynamics in nonlocal and nonreciprocal media

The linear response of a macroscopic material to externally applied electromagnetic fields can go beyond the scope of simple descriptions via electric permittivities and magnetic permeabilities [1]. In particular, cross-susceptibilities naturally arise in chiral (meta-)materials [2], topological insulators [3] or moving media [4]. In the latter case nonlocal responses arise [5] with the additional complication that Onsager reciprocity [6] fails to hold. Onsager reciprocity, the electrodynamic manifestation of time-reversal symmetry, would also be violated in Tellegen media [7], including the recently proposed perfect electromagnetic conductor that continuously interpolates between a perfect conductor and an infinitely permeable material [8].

Chiral metamaterials with cross-susceptibilities have been constructed based on nanoscale chiral objects, such as a helix [9]. This leads to a discriminatory response of the medium to left- and right-circularly polarized light. This central feature of chiral media is important in biological systems due to the prevalence of left-handed objects in the processes crucial to life [10]. Furthermore, chiral meta-materials have been discussed as candidates for repulsive Casimir forces [11]. It should be noted that repulsive forces for magnetoelectric media were originally discussed for dielectric plates interacting with magnetic plates [12]. To implement these effects with metamaterials, the anisotropic response of the medium needs to be taken into account [13].

Topological insulators are a novel class of materials which behave as insulators in their bulk phase but allow for conduction on the surface [14, 15]. Time reversal symmetry is an important feature which ensures an extremely high stability of the surface currents. The latter make topological insulators a promising candidate for quantum computing [16]. It has recently been predicted that topological insulators [3] (or materials with a Chern–Simons interaction [17]) could be used to realize repulsive Casimir forces. A related phenomenon is the fractional quantum hall effect where the Hall current takes fractional values due to electron–electron interactions [18]. This medium can be nonlocal [19] and in contrast to topological insulators it can violate time-reversal symmetry [20] and hence Onsager reciprocity.

The impact of electric versus magnetic material properties can be studied in a systematic way by means of a duality transformation [21]. It has recently been shown that macroscopic QED [22] in isotropic magnetoelectrics obeys a discrete duality symmetry [23]. This has immediate consequences for dispersion forces in free space.

The successes in the realization of the above mentioned novel materials open the perspective on a range of new quantum phenomena related to photon-induced matter interactions, quantum dynamics and (possibly irreversible) quantum-light propagation. To make such studies possible, we will construct a quantum theory of the electromagnetic field in the most general linear absorbing media, including nonlocal, bianisotropic and Onsager reciprocity violating materials. A recent theory based on canonical quantization is a valuable step in this direction [24], which does not yet consider the most general nonlocal media (apart from the above mentioned moving media).

In addition, our theory shall answer the question under which circumstances duality can be realized as a continuous symmetry of the Maxwell equations in media; and it will shed light on the generalizations necessary to discuss moving media and quantum friction.

We begin by recalling a quantization procedure of the electromagnetic field in the presence of an absorbing medium. For alternative methods for the quantization procedure, see [25–27] and references therein. In a linearly responding medium, the effect of an external electromagnetic field on the matter can be given by Ohm's law in its most general form

$$\begin{eqnarray} &&{\bf j}_{{\rm in}}({\bf r},t) =\int_{-\infty}^{\infty}{\rm d}\tau\int{\rm d}^3r'\, {\boldsymbol{\sf{Q}}}({\bf r},{\bf r}',\tau)\cdot {\bf E}({\bf r}',t-\tau) +{\bf j}_{{\rm N}}({\bf r},t). \end{eqnarray} \tag{ 1 }$$

Here, $\boldsymbol{{\sf{Q}}}(\mathbf{r},\mathbf{r}',\tau)$ is the conductivity tensor and j_N(r,t) is the random noise current required to fulfil the fluctuation–dissipation theorem (17) as given below. Causality requires that $\boldsymbol{\sf{Q}}(\mathbf{r},\mathbf{r}',\tau)= \boldsymbol{\sf{0}}$ for cτ < |r − r'|, in particular for all τ < 0 [28]. In frequency space, Ohm's law takes the simpler form

$$\begin{eqnarray} &&{\bf j}_{{\rm in}}({\bf r},\omega) =\int{\rm d}^3r'\, {\boldsymbol {\sf{Q}}}({\bf r},{\bf r}',\omega)\cdot {\bf E}({\bf r}',\omega) +{\bf j}_{{\rm N}}({\bf r},\omega) \end{eqnarray} \tag{ 2 }$$

with

$$\begin{eqnarray} &&{\boldsymbol {\sf{Q}}}({\bf r},{\bf r}',\omega) =\int_0^\infty{\rm d}\tau\,{\rm e}^{{\rm i}\,\omega\tau} {\boldsymbol {\sf{Q}}}({\bf r},{\bf r}',\tau). \end{eqnarray} \tag{ 3 }$$

As a result of the causality requirement the conductivity obeys the Schwarz reflection principle,

$$\begin{eqnarray} &&{\boldsymbol {\sf{Q}}}^\ast({\bf r},{\bf r}',\omega) ={\boldsymbol {\sf{Q}}}({\bf r},{\bf r}',-\omega^\ast)\quad \forall\;{\bf r},{\bf r}',\omega. \end{eqnarray}\noindent \tag{ 4 }$$

Quantization is achieved by specifying the commutator

$$\begin{eqnarray} &&[\hat{{\bf j}}_{{\rm N}}({\bf r},\omega), \hat{{\bf j}}{}_{{\rm N}}^\dagger ({\bf r}',\omega')] =\frac{\hbar\omega}{\pi}\,\mathcal{R}{\rm e}[{\boldsymbol {\sf{Q}}}({\bf r},{\bf r}',\omega)] \delta(\omega-\omega'), \end{eqnarray} \tag{ 5 }$$

where we have introduced generalized real and imaginary parts of a tensor field according to

$$\begin{eqnarray} &&\mathcal{R}{\rm e} {\boldsymbol{\sf{T}}}({\bf r},{\bf r}')=\textstyle\frac{1}{2}\bigl[ {\boldsymbol{\sf{T}}}({\bf r},{\bf r}') +{\boldsymbol{\sf{T}}}^\dagger({\bf r}',{\bf r})\bigr], \end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray} &&\mathcal{I}{\rm m}{\boldsymbol {\sf{T}}}({\bf r},{\bf r}')=\frac{1}{2{\rm i}}\bigl[ {\boldsymbol {\sf{T}}}({\bf r},{\bf r}') -{\boldsymbol {\sf{T}}}^\dagger({\bf r}',{\bf r})\bigr]. \end{eqnarray} \tag{ 7 }$$

They reduce to ordinary real and imaginary parts for orthogonal tensor fields that are reciprocal, i.e. ${\boldsymbol{\sf{T}}}^{\sf T}(\mathbf{r}',\mathbf{r})=\boldsymbol{\sf{T}}(\mathbf{r},\mathbf {r}')$. The other nontrivial current commutators follow alternatively via the rules $[\hat {b}, \hat {a}] = - [\hat {a}, \hat {b}]$ or $[\hat {a}^{\dagger}, \hat {b}^{\dagger}] = - [\hat {a}, \hat {b}]^{\dagger}$;

$$\begin{eqnarray} \Bigl[\hat{{\bf j}}^{\dagger}_{{\rm N}}({\bf r},\omega), \hat{{\bf j}}{}_{{\rm N}} ({\bf r}',\omega')\Bigr] & & = - \frac{\hbar\omega}{\pi}\,\mathcal{R}{\rm e}[{\boldsymbol {\sf{Q}}}^{*}({\bf r},{\bf r}',\omega)] \delta(\omega-\omega') \nonumber \\ & & = -\frac{\hbar\omega}{\pi}\,\mathcal{R}{\rm e}[{\boldsymbol {\sf{Q}}}^{\sf{T}}({\bf r}',{\bf r},\omega)] \delta(\omega-\omega'). \end{eqnarray} \tag{ 8 }$$

The fact that the right-hand side of this expression is a Hermitian tensor field guarantees the consistency of the commutation relations.

Combining Ohm's law with Maxwell's equations [ $\mathrm {i}\,\omega \hat {\rho }_{\mathrm {\,in}}(\mathbf {r}, \omega )= \mathbf {\nabla }\cdot \hat {\mathbf {j}}_{\mathrm {in}}(\mathbf {r}, \omega )$]

$$\begin{eqnarray} {\bf \nabla}\cdot \hat{{\bf E}}({\bf r}, \omega) &=&\frac{\hat{\rho}_{{\rm \,in}}({\bf r}, \omega)}{\varepsilon_0},\quad {\bf \nabla}\times \hat{{\bf E}}({\bf r}, \omega) -{\rm i}\,\omega\hat{{\bf B}}({\bf r}, \omega)={\bf 0}, \end{eqnarray}\noindent \tag{ 9 }$$

$$\begin{eqnarray} {\bf \nabla}\cdot \hat{{\bf B}}({\bf r}, \omega)&=&0,\quad {\bf \nabla}\times \hat{{\bf B}}({\bf r}, \omega) +\frac{{\rm i}\,\omega}{c^2}\,\hat{{\bf E}}({\bf r}, \omega) =\mu_0\hat{{\bf j}}_{{\rm in}}({\bf r}, \omega), \end{eqnarray} \tag{ 10 }$$

one finds that the electric field obeys a generalized inhomogeneous Helmholtz equation of the form

$$\begin{eqnarray} &&\fl \biggl[{\bf \nabla}\times {\bf \nabla}\times-\,\frac{\omega^2}{c^2}\biggr] \hat{{\bf E}}({\bf r},\omega) -{\rm i}\,\mu_0\omega\int{\rm d}^3r' {\boldsymbol {\sf{Q}}}({\bf r},{\bf r}',\omega)\cdot \hat{{\bf E}}({\bf r}',\omega) ={\rm i}\,\mu_0\omega\hat{{\bf j}}_{{\rm N}}({\bf r},\omega). \end{eqnarray} \tag{ 11 }$$

With the help of the Green function $\boldsymbol{{\sf{G}}}(\mathbf{r},\mathbf{r}',\omega)$ of the Helmholtz equation, defined by

$$\begin{eqnarray} &&\fl \biggl[{\bf \nabla}\times {\bf \nabla}\times -\,\frac{\omega^2}{c^2}\biggr] {\boldsymbol {\sf{G}}}({\bf r},{\bf r}',\omega) -{\rm i}\,\mu_0\omega \int {\rm d}^3s\, {\boldsymbol {\sf{Q}}}({\bf r},{\bf s},\omega)\cdot {\boldsymbol {\sf{G}}}({\bf s},{\bf r}',\omega) ={\boldsymbol \delta}({\bf r}-{\bf r}'), \end{eqnarray} \tag{ 12 }$$

where $\boldsymbol{{\sf{G}}}(\mathbf{r},\mathbf{r}',\omega)\to\boldsymbol{\sf{0}}$ for |r − r'| → ∞, the formal solution to the integro-differential equation (11) reads

$$\begin{equation} \hat{{\bf E}}({\bf r}, \omega) = {\rm i}\,\mu_{0}\omega[{\boldsymbol {\sf{G}}}(\omega)\star \hat{{\bf j}}_{{\rm N}}(\omega)]({\bf r}). \end{equation} \tag{ 13 }$$

Here, $[ \boldsymbol {{\sf {G}}}\star \hat {\mathbf {j}}_{\mathrm {N}}]$ is an abbreviation denoting the spatial convolution

$$\begin{equation*} [{\boldsymbol {\sf{G}}}\star \hat{{\bf j}}_{{\rm N}}]({\bf r}) \equiv \int{\rm d}^3r' \,{\boldsymbol {\sf{G}}}({\bf r},{\bf r}')\cdot \hat{{\bf j}}_{{\rm N}}({\bf r}'). \end{equation*}$$

By virtue of its definition (12), the Green tensor inherits the Schwarz reflection principle from the conductivity tensor (4),

$$\begin{eqnarray} &&{\boldsymbol {\sf{G}}}^\ast({\bf r},{\bf r}',\omega) ={\boldsymbol {\sf{G}}}({\bf r},{\bf r}',-\omega^\ast)\quad \forall\;{\bf r},{\bf r}',\omega. \end{eqnarray} \tag{ 14 }$$

However, as a major departure from previous treatments, we do not require the conductivity to obey reciprocity, i.e. the relation $\boldsymbol{{\sf{Q}}}^{\sf{T}}(\mathbf{r}',\mathbf{r},\omega)=\boldsymbol{{\sf{Q}}}(\mathbf{r},\mathbf{r}',\omega)$ does not necessarily hold. As a consequence, the Green tensor will not obey the Onsager principle, i.e. the relation

$$\begin{eqnarray} &&{\boldsymbol {\sf{G}}}^{\sf{T}}({\bf r}',{\bf r},\omega) ={\boldsymbol {\sf{G}}}({\bf r},{\bf r}',\omega) \end{eqnarray} \tag{ 15 }$$

will not hold in general. Recall that the Onsager principle, applied to electromagnetic field propagation, states a reversibility of optical paths [6]. According to (13), the Green tensor governs the relation between a source current j at r' along a direction e₂ and the generated electric field E at r along a direction e₁. If (15) holds, then the Onsager principle states that the roles of source and field can be reversed. A source current j at r along a direction e₁ would then give rise to an electric field E at r' along a direction e₂. In a configuration involving nonreciprocal media, this is not necessarily the case. Despite the extension to nonreciprocal media, it is still possible to derive the useful integral relation (appendix)

$$\begin{equation} \mu_0\omega [{\boldsymbol {\sf{G}}}(\omega)\star \mathcal{R}{\rm e}{\boldsymbol {\sf{Q}}}(\omega)\star {\boldsymbol {\sf{G}}}^{\dagger}(\omega)] ({\bf r},{\bf r}') = \mathcal{I}{\rm m}\,{\boldsymbol {\sf{G}}}({\bf r},{\bf r}',\omega). \end{equation} \tag{ 16 }$$

It generalizes the result from [29] to the case where Onsager reciprocity does not hold.

The theory thus far is analogous to classical electromagnetism in an absorbing medium under the assumption of classical fluctuating current sources, $\hat {\mathbf {j}}_{\mathrm {N}}\mapsto \mathbf {j}_{\mathrm {N}}$. Their strengths are governed by the fluctuation–dissipation theorem in the classical (high-temperature) limit [30].

Introducing the ground state |{0}〉 of the medium-field system according to ${\hat {\mathbf {j}}_{\mathrm {N}}(\mathbf {r},\omega )|\{0\}\rangle =\mathbf {0}}$, the currents satisfy the fluctuation–dissipation theorem as an immediate consequence of (5),

$$\begin{eqnarray} &&\bigl\langle\bigl\{ \Delta\hat{{\bf j}}_{{\rm N}}({\bf r},\omega), \Delta\hat{{\bf j}}^\dagger_{{\rm N}}({\bf r}',\omega')\bigr\} \bigr\rangle = \frac{\hbar}{\pi}\,\mathcal{I}{\rm m}[{\rm i}\,\omega {\boldsymbol {\sf{Q}}}({\bf r},{\bf r}',\omega)]\delta(\omega-\omega'). \end{eqnarray} \tag{ 17 }$$

Combining (5) and (13), one finds that the fluctuations of the electric field are also consistent with the fluctuation–dissipation theorem, as required:

$$\begin{eqnarray} &&\bigl\langle\bigl\{\Delta\hat{{\bf E}}({\bf r},\omega), \Delta\hat{{\bf E}}^\dagger({\bf r}',\omega') \bigr\}\bigr\rangle = \frac{\hbar}{\pi}\,\mathcal{I}{\rm m}\bigl[\mu_0\omega^2 {\boldsymbol {\sf{G}}}({\bf r},{\bf r}',\omega)\bigr] \delta(\omega-\omega'). \end{eqnarray} \tag{ 18 }$$

In order to verify the canonical equal-time commutation relations, we introduce the vector potential for the electromagnetic field in the Coulomb gauge, $\hat {\mathbf {A}}(\mathbf {r},\omega )=\hat {\mathbf {E}}^\perp (\mathbf {r}, \omega )/(\mathrm {i}\,\omega )$ (⊥: transverse part). Using (5) and (13), one finds

$$\begin{eqnarray} &&\bigl[\hat{{\bf E}}({\bf r},\omega), \hat{{\bf A}}{}^\dagger({\bf r}',\omega')\bigr] =\frac{{\rm i}\hbar\mu_0\omega}{\pi} \,\mathcal{I}{\rm m}[{\boldsymbol {\sf{G}}}^\perp({\bf r},{\bf r}',\omega)] \delta(\omega-\omega') \end{eqnarray} \tag{ 19 }$$

and hence

$$\begin{eqnarray} &&\bigl[\hat{{\bf E}}({\bf r}),\hat{{\bf A}}({\bf r}')\bigr] =\frac{\hbar\mu_0}{2\pi} \int_{-\infty}^\infty{\rm d}\omega\,\omega\, \bigl[{\boldsymbol {\sf{G}}}^\perp({\bf r},{\bf r}',\omega) +{}^\perp{\boldsymbol {\sf{G}}}^{\sf{T}}({\bf r}',{\bf r},\omega)\bigr], \end{eqnarray} \tag{ 20 }$$

where the Schwarz reflection principle (14) has been used. Use has been made of the left- and right-sided transverse projections, which are defined, respectively, as

$$\begin{equation} {}^{\perp}{\boldsymbol T}= {\boldsymbol \delta}^\perp\star {\boldsymbol T} ,\quad {\boldsymbol T}^{\perp}={\boldsymbol T}\star {\boldsymbol \delta}^\perp . \end{equation} \tag{ 21 }$$

Closing the integration contour in the upper half of the complex ω plane, where the Green's function is analytic, and using the asymptote $(\omega^2/c^2)\boldsymbol{{\sf{G}}}(\mathbf{r},\mathbf{r}',\omega)\to-{\boldsymbol \delta}(\mathbf{r}-\mathbf{r}')$ for |ω| → ∞, one finds the canonical commutation relation from free-space QED,

$$\begin{equation} \bigl[\hat{{\bf E}}({\bf r}),\hat{{\bf A}}({\bf r}')\bigr] =\frac{{\rm i}\hbar}{\varepsilon_0}\, {\boldsymbol \delta}^\perp({\bf r}-{\bf r}'), \end{equation} \tag{ 22 }$$

as required. We now introduce the bosonic creation and annihilation operators of the matter-field system, $\hat {\mathbf {f}}{}^\dagger$ and $\hat {\mathbf {f}}$ according to the prescription

$$\begin{eqnarray} &&\hat{{\bf j}}_{{\rm N}}(\omega) =\sqrt{\frac{\hbar\omega}{\pi}}\,{\boldsymbol {\sf{R}}}(\omega)\star \hat{{\bf f}}(\omega), \end{eqnarray} \tag{ 23 }$$

where $\boldsymbol{{\sf{R}}}$ is a square root of the positive definite tensor field $\mathcal {R}\mathrm {e}[\boldsymbol {{\sf {Q}}}]$,

$$\begin{eqnarray} &&{\boldsymbol {\sf{R}}}(\omega)\star {\boldsymbol {\sf{R}}}^\dagger(\omega) =\mathcal{R}{\rm e}[{\boldsymbol {\sf{Q}}}(\omega)]. \end{eqnarray} \tag{ 24 }$$

This solution is only unique up to a unitary matrix which does not affect the physical results [29]. Together with equation (5), this ensures bosonic commutation relations,

$$\begin{equation} [\hat{{\bf f}}({\bf r},\omega), \hat{{\bf f}}^\dagger({\bf r}',\omega')] = {\boldsymbol \delta}({\bf r}-{\bf r}') \delta(\omega-\omega'). \end{equation} \tag{ 25 }$$

The Hamiltonian of the medium-field system is then

$$\begin{eqnarray} &&\hat{H}_{{\rm F}} =\int{\rm d}^3r \int_0^\infty{\rm d}\omega\,\hbar\omega\, \hat{{\bf f}}^\dagger({\bf r},\omega) \cdot \hat{{\bf f}}({\bf r},\omega). \end{eqnarray} \tag{ 26 }$$

It leads to the free evolution of the dynamical variables as $\hat {\mathbf {f}}(\mathbf {r},\omega ,t)= \hat {\mathbf {f}}(\mathbf {r},\omega )\,\mathrm {e}^{-\mathrm {i}\,\omega t}$; hence Maxwell's equations for the electromagnetic-field operators in the Heisenberg picture are valid by construction.

The properties of media with spatially nonlocal or local responses can alternatively be described by their permittivity, permeability and magnetoelectric susceptibilities. To begin, it is convenient to cast the inhomogeneous Maxwell equations (10) into the forms (we drop the spatial and frequency arguments from now on)

$$\begin{eqnarray} &&{\bf \nabla}\cdot \hat{{\bf D}}=0,\quad {\bf \nabla}\times \hat{{\bf H}} +{\rm i}\,\omega\hat{{\bf D}} ={\bf 0} \end{eqnarray} \tag{ 27 }$$

with

$$\begin{eqnarray} &&\hat{{\bf D}}=\varepsilon_0\hat{{\bf E}} +\hat{{\bf P}},\quad \hat{{\bf H}}=\frac{1}{\mu_0}\,\hat{{\bf B}} -\hat{{\bf M}}. \end{eqnarray} \tag{ 28 }$$

The polarization and magnetization fields respond linearly to the electric and magnetic fields,

$$\begin{eqnarray} \hat{{\bf P}} &=&\varepsilon_0({\boldsymbol \varepsilon}-{\boldsymbol \xi}\star {\boldsymbol \mu}^{-1}\star {\boldsymbol \zeta}-{\boldsymbol {\sf{I}}})\star \hat{{\bf E}} +Z_0^{-1}{\boldsymbol \xi}\star {\boldsymbol \mu}^{-1}\star \hat{{\bf B}} +\hat{{\bf P}}_{{\rm N}}, \end{eqnarray}\noindent \tag{ 29 }$$

$$\begin{eqnarray} \hat{{\bf M}} &=&Z_0^{-1}{\boldsymbol \mu}^{-1}\star {\boldsymbol \zeta} \star \hat{{\bf E}} +\mu_0^{-1}({\boldsymbol {\sf{I}}}-{\boldsymbol \mu}^{-1})\star \hat{{\bf B}} +\hat{{\bf M}}_{{\rm N}}. \end{eqnarray}\noindent \tag{ 30 }$$

The medium is characterized by its permittivity, ε(r,r',ω), its permeability, μ(r,r',ω) and its magnetoelectric susceptibilities, ξ(r,r',ω) and ζ(r,r',ω). $\hat {\mathbf {P}}_{\mathrm {N}}(\mathbf {r},\omega )$ and $\hat {\mathbf {M}}_{\mathrm {N}}(\mathbf {r},\omega )$ denote the noise polarization and noise magnetization, respectively, and $Z_0=\sqrt {\mu _0/\varepsilon _0}\,$ is the vacuum impedance. In the case of nonlocal media the permittivity, permeability and magnetoelectric susceptibilities are functions of two independent spatial variables, whereas in a locally responding media they read ε(r,r',ω) = ε(r,ω)δ(r − r'), μ(r,r',ω) = μ(r,ω)δ(r − r'), ξ(r,r',ω) = ξ(r,ω)δ(r − r') and ζ(r,r',ω) = ζ(r,ω)δ(r − r'). By combining (28)–(30), the constitutive relations can be given in the more familiar form (see [31] for the nonconducting case)

$$\begin{eqnarray} \hat{{\bf D}} &=&\varepsilon_0{\boldsymbol \varepsilon}\star \hat{{\bf E}} +c^{-1}{\boldsymbol \xi}\star \hat{{\bf H}} +\hat{{\bf P}}_{{\rm N}} +c^{-1}{\boldsymbol \xi}\star \hat{{\bf M}}_{{\rm N}}, \end{eqnarray}\noindent \tag{ 31 }$$

$$\begin{eqnarray} \hat{{\bf B}} &=&c^{-1}{\boldsymbol \zeta}\star \hat{{\bf E}} +\mu_0{\boldsymbol \mu}\star \hat{{\bf H}} +\mu_0{\boldsymbol \mu}\star \hat{{\bf M}}_{{\rm N}}, \end{eqnarray}\noindent \tag{ 32 }$$

where the notational distinction between locally and nonlocally responding media as given above applies.

In order to distinguish reciprocal magnetoelectric susceptibilities from nonreciprocal ones, as previously discussed in the nondispersive case [32], one commonly writes ${\boldsymbol \xi}={\boldsymbol \chi}^{\sf{T}}-\mathrm{i}{\boldsymbol \kappa}^{\sf{T}}$ and ζ = χ + iκ. The chirality tensor ${\boldsymbol \kappa}=({\boldsymbol \zeta}-{\boldsymbol \xi}^{\sf{T}})/(2\mathrm{i})$ represents the reciprocal magnetoelectric response; whereas the nonreciprocal magnetoelectric tensor ${\boldsymbol \chi}=({\boldsymbol \zeta}+{\boldsymbol \xi}^{\sf{T}})/2$ vanishes for a reciprocal medium.

By combining Maxwell's equations with the constitutive relations (31) and (32), we note that the respective Green tensor is the solution to equation (12) with

$$\begin{eqnarray} &&\fl {\boldsymbol {\sf{Q}}} =({\rm i}\,\mu_0\omega)^{-1}{\bf \nabla}\times ({\boldsymbol \mu}^{-1}-{\boldsymbol {\sf{I}}}) \times \overleftarrow{{\bf \nabla}} + Z_0^{-1}{\bf \nabla}\times {\boldsymbol \mu}^{-1}\star {\boldsymbol \zeta}+ Z_0^{-1}{\boldsymbol \xi}\star {\boldsymbol \mu}^{-1} \times \overleftarrow{{\bf \nabla}} \nonumber\\&&-{\rm i}\varepsilon_0\omega ({\boldsymbol \varepsilon}-{\boldsymbol \xi}\star {\boldsymbol \mu}^{-1}\star {\boldsymbol \zeta}-{\boldsymbol {\sf{I}}}) \end{eqnarray} \tag{ 33 }$$

and

$$\begin{equation} \hat{{\bf j}}_{{\rm N}} =-{\rm i}\,\omega\hat{{\bf P}}_{{\rm N}} +{\bf \nabla}\times \hat{{\bf M}}_{{\rm N}}, \end{equation} \tag{ 34 }$$

where $[\boldsymbol {{\sf {T}}}\times \overleftarrow {\mathbf {\nabla }}]_{ij}(\mathbf {r},\mathbf {r}')= \epsilon _{jkl}\partial _l'T_{ik}(\mathbf {r},\mathbf {r}')$ denotes a derivative acting on the second argument of a tensor function. The Green tensor for the electric field (13) solves

$$\begin{eqnarray} &&\fl \biggl[{\bf \nabla}\times {\boldsymbol \mu}^{-1}\star {\bf \nabla} \times -\frac{{\rm i}\,\omega}{c}{\bf \nabla}\times {\boldsymbol \mu}^{-1}\star {\boldsymbol \zeta} +\frac{{\rm i}\,\omega}{c}{\boldsymbol \xi}\star {\boldsymbol \mu}^{-1}\star {\bf \nabla}-\,\frac{\omega^2}{c^2} ({\boldsymbol \varepsilon} - {\boldsymbol \xi}\star {\boldsymbol \mu}^{-1} \star {\boldsymbol \zeta})\biggr]\star {\boldsymbol {\sf{G}}} ={\boldsymbol \delta}. \end{eqnarray} \tag{ 35 }$$

The commutation relations for $\hat {\mathbf {P}}_{\mathrm {N}}$ and $\hat {\mathbf {M}}_{\mathrm {N}}$ can be deduced by substituting the real parts of (33) and (34) into (5),

$$\begin{eqnarray} &&\fl \left[\hat{{\bf P}}_{{\rm N}}({\bf r},\omega), \hat{{\bf P}}{}_{{\rm N}}^\dagger ({\bf r}',\omega')\right] =\frac{\varepsilon_0\hbar}{\pi}\,\mathcal{I}{\rm m} \left\{ {\boldsymbol \varepsilon}(\omega) - [{\boldsymbol \xi}(\omega)\star {\boldsymbol \mu}^{-1} (\omega) \star {\boldsymbol \zeta}(\omega)]\right\}({\bf r},{\bf r}') \delta(\omega-\omega'), \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} & &\fl \left[\hat{{\bf P}}_{{\rm N}}({\bf r},\omega), \hat{{\bf M}}{}_{{\rm N}}^\dagger ({\bf r}',\omega')\right] = \frac{\hbar}{2\pi{\rm i} Z_0}\! \left\{\![{\boldsymbol \xi}(\omega)\star {\boldsymbol \mu}^{-1} (\omega)] - [{\boldsymbol \zeta}^\dagger(\omega) \star {\boldsymbol \mu}^{-1\dagger}(\omega)] \right\}\!({\bf r}, {\bf r}') \delta(\omega - \omega'), \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} & &\fl \left[\hat{{\bf M}}_{{\rm N}}({\bf r},\omega), \hat{{\bf P}}{}_{{\rm N}}^\dagger ({\bf r}',\omega')\right] = \frac{\hbar}{2\pi{\rm i} Z_0}\! \left\{\![{\boldsymbol \mu}^{-1}(\omega)\star {\boldsymbol \zeta} (\omega)] - [{\boldsymbol \mu}^{-1\dagger}(\omega) \star {\boldsymbol \xi}^{\dagger}(\omega)] \right\}\!({\bf r}, {\bf r}') \delta(\omega - \omega'), \end{eqnarray} \tag{ 38 }$$

$$\begin{equation} \left[\hat{{\bf M}}_{{\rm N}}({\bf r},\omega), \hat{{\bf M}}{}_{{\rm N}}^\dagger ({\bf r}',\omega')\right] =-\frac{\hbar}{\pi\mu_0}\,\mathcal{I}{\rm m} [{\boldsymbol \mu}^{-1}({\bf r}, {\bf r}', \omega)] \delta(\omega - \omega'). \end{equation} \tag{ 39 }$$

We now introduce the bosonic creation and annihilation operators with commutation relations

$$\begin{equation} \left[\hat{{\bf f}}_{\lambda}({\bf r},\omega), \hat{{\bf f}}_{\lambda'}^\dagger({\bf r}',\omega')\right] =\delta_{\lambda\lambda'}{\boldsymbol \delta}({\bf r}-{\bf r}') \delta(\omega-\omega'),\quad (\lambda,\lambda'=e,m) \end{equation} \tag{ 40 }$$

according to

$$\begin{eqnarray} &&\left(\begin{array}{@{}c@{}}\hat{{\bf P}}_{{\rm N}}\\ \hat{{\bf M}}_{{\rm N}}\end{array}\right) =\sqrt{\frac{\hbar}{\pi}}\,\mathcal{R}\star \left( \begin{array}{@{}c@{}}\hat{{\bf f}}_e\\ \hat{{\bf f}}_m\end{array}\right)\!, \end{eqnarray} \tag{ 41 }$$

where the (6 × 6)-matrix $\mathcal {R}$ is a root of

$$\begin{eqnarray} \mathcal{R}\star \mathcal{R}^\dagger =\left(\begin{array}{@{}cc@{}}\varepsilon_0\mathcal{I}{\rm m} [{\boldsymbol \varepsilon} - {\boldsymbol \xi}\star {\boldsymbol \mu}^{-1}\star {\boldsymbol \zeta}] & \displaystyle\frac{{\boldsymbol \xi}\star {\boldsymbol \mu}^{-1} - {\boldsymbol \zeta}^\dagger\star {\boldsymbol \mu}^{-1\dagger}}{2{\rm i} Z_0}\\ \displaystyle\frac{{\boldsymbol \mu}^{-1}\star {\boldsymbol \zeta} - {\boldsymbol \mu}^{-1\dagger}\star {\boldsymbol \xi}^\dagger}{2{\rm i} Z_0} & \displaystyle-\frac{\mathcal{I}{\rm m}[{\boldsymbol \mu}^{-1}]}{\mu_0} \end{array}\right)\!. \end{eqnarray} \tag{ 42 }$$

The Hamiltonian of the body–field system is again quadratic and diagonal in the bosonic variables,

$$\begin{eqnarray} &&\hat{H}_{{\rm F}} =\sum_{\lambda=e,m}\int{\rm d}^3r \int_0^\infty{\rm d}\omega\,\hbar\omega\, \hat{{\bf f}}_{\lambda}^\dagger({\bf r},\omega) \cdot\hat{{\bf f}}_{\lambda}({\bf r},\omega). \end{eqnarray} \tag{ 43 }$$

Note that (28)–(30) imply a separation of the internal current density into electric and magnetic parts, $\hat {\mathbf {j}}_{\mathrm {in}} =-\mathrm {i}\,\omega \hat {\mathbf {P}} + \mathbf {\nabla }\times \hat {\mathbf {M}}$. This separation and the resulting explicit field quantization is not unique in spatially dispersive (i.e. nonlocal) media, as the magnetization field can be absorbed into the transverse part of the polarization field [28]. While a local magnetoelectric medium can always be described in terms of a (nonlocal) conductivity without reference to magnetic properties, the equivalent description in terms of a local permittivity, permeability and cross-susceptibility is much more accessible. These parameters are often known experimentally and they allow for a classification of electromagnetic responses.

An electromagnetic system separated into distinct electric and magnetic causes and effects can be subject to a duality transformation operation, that is, a global exchange of the electric and magnetic properties. A system invariant under such an operation is said to possess duality invariance as a symmetry [21]. This symmetry can be exploited in order to simplify the computation of dispersion forces involving, say magnetizable media, from known dispersion forces between polariable media [23].

By introducing dual-pair notation $(\hat {\mathbf {E}}^{\sf{T}},Z_0\hat {\mathbf {H}}^{\sf{T}})^{\sf{T}}$, $(Z_0\hat {\mathbf {D}}^{\sf{T}},\hat {\mathbf {B}}^{\sf{T}})^{\sf{T}}$, we may write the Maxwell equations (9) and (27) in the compact form

$$\begin{eqnarray} &&{\bf \nabla}\cdot \left(\begin{array}{@{}c@{}}Z_0\hat{{\bf D}}\\ \hat{{\bf B}}\end{array}\right) =\left(\begin{array}{@{}c@{}}0\\ 0\end{array}\right)\!, \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray} &&{\bf \nabla}\times \left(\begin{array}{@{}c@{}}\hat{{\bf E}}\\ Z_0\hat{{\bf H}}\end{array}\right) -{\rm i}\,\omega \left(\begin{array}{@{}cc@{}}0 & 1\\ -1 & 0\end{array}\right) \left(\begin{array}{@{}c@{}}Z_0\hat{{\bf D}}\\ \hat{{\bf B}}\end{array}\right) =\left(\begin{array}{@{}c@{}}{\bf 0}\\ {\bf 0}\end{array}\right)\!. \end{eqnarray} \tag{ 45 }$$

The constitutive relations (31) and (32) in condensed form read

$$\begin{eqnarray} &&\left(\begin{array}{@{}c@{}}Z_0\hat{{\bf D}}\\ \hat{{\bf B}}\end{array}\right) =\frac{1}{c}\left(\begin{array}{@{}cc@{}}{\boldsymbol \varepsilon} & {\boldsymbol \xi}\\ {\boldsymbol \zeta} & {\boldsymbol \mu}\end{array}\right)\star \left(\begin{array}{@{}c@{}}\hat{{\bf E}}\\ Z_0\hat{{\bf H}}\end{array}\right) + \mathcal{A}\star \left(\begin{array}{@{}c@{}}Z_0\hat{{\bf P}}_{{\rm N}}\\ \mu_0\hat{{\bf M}}_{{\rm N}}\end{array}\right) \end{eqnarray} \tag{ 46 }$$

with

$$\begin{equation} {\mathcal{A}}= \left(\begin{array}{@{}cc@{}}{\boldsymbol {\sf{I}}} & {\boldsymbol \xi}\\ 0 & {\boldsymbol \mu}\end{array} \right)\!. \end{equation} \tag{ 47 }$$

Maxwell's equations are invariant under duality transformations

$$\begin{eqnarray} &&\left(\begin{array}{@{}c@{}}{\bf x}\\ {\bf y}\end{array}\right)^{\circledast} =D(\theta) \left(\begin{array}{@{}c@{}}{\bf x}\\ {\bf y}\end{array}\right), \quad D(\theta) =\left(\begin{array}{@{}cc@{}}\cos\theta & \sin\theta\\ -\sin\theta & \cos\theta\end{array}\right)\!, \end{eqnarray} \tag{ 48 }$$

because D(θ) is a symplectic matrix. From the constitutive relations, as shown in (46), we find the transformed medium response functions

$$\begin{equation} \left(\begin{array}{@{}c@{}}{\boldsymbol \varepsilon}\\ {\boldsymbol \zeta}\\ {\boldsymbol \xi}\\ {\boldsymbol \mu}\end{array}\right)^{{\circledast}} = \mathcal{D} (\theta)\left(\begin{array}{@{}c@{}}{\boldsymbol \varepsilon}\\ {\boldsymbol \zeta}\\ {\boldsymbol \xi}\\ {\boldsymbol \mu}\end{array}\right) \end{equation} \tag{ 49 }$$

with

$$\begin{equation} \mathcal{D}(\theta)= D(\theta)\otimes D(\theta) \end{equation} \tag{ 50 }$$

as well as

$$\begin{eqnarray} &&\left(\begin{array}{@{}c@{}}Z_0\hat{{\bf P}}_{{\rm N}}\\ \mu_0\hat{{\bf M}}_{{\rm N}}\end{array}\right)^{\circledast} =\mathcal{A}^{{\circledast}-1}\star D(\theta)\mathcal{A}\star \left(\begin{array}{@{}c@{}}Z_0\hat{{\bf P}}_{{\rm N}}\\ \mu_0\hat{{\bf M}}_{{\rm N}}\end{array}\right)\!, \end{eqnarray} \tag{ 51 }$$

where

$$\begin{equation} \mathcal{A}^{-1} =\left(\begin{array}{@{}cc@{}}{\boldsymbol I} & -{\boldsymbol \xi}\cdot{\boldsymbol \mu}^{-1} \\ {\boldsymbol 0} & {\boldsymbol \mu}^{-1} \end{array}\right)\!. \end{equation} \tag{ 52 }$$

It is worth discussing a few special cases of bianisotropic media, their characteristic features and behaviour under duality transformations:

• Local media [ε(r,r',ω) = ε(r,ω)δ(r − r'), similarly for the other response functions]: Convolution operators reduce to ordinary matrix products, compatible with all other special cases below.
• Isotropic media (${\boldsymbol \varepsilon}=\varepsilon\boldsymbol{{\sf{I}}}$, ${\boldsymbol \mu}=\mu\boldsymbol{{\sf{I}}}$, ${\boldsymbol \xi}={\boldsymbol \zeta}=\boldsymbol{\sf{0}}$): Onsager reciprocity (15) holds; $\hat {\mathbf {P}}_{\mathrm {N}}$ and $\hat {\mathbf {M}}{}^\dagger _{\mathrm {N}}$ commute; generalized real and imaginary parts reduce to ordinary ones; discrete duality symmetry.
• Bi-isotropic media (${\boldsymbol \varepsilon}=\varepsilon\boldsymbol{{\sf{I}}}$, ${\boldsymbol \mu}=\mu\boldsymbol{{\sf{I}}}$, ${\boldsymbol \xi}=\xi\boldsymbol{{\sf{I}}}$, ${\boldsymbol \zeta}=\zeta\boldsymbol{{\sf{I}}}$): generalized real and imaginary parts in (36) and (39) reduce to ordinary ones; continuous duality symmetry.
• Anisotropic media (${\boldsymbol \xi}={\boldsymbol \zeta}=\boldsymbol{\sf{0}}$): $\hat {\mathbf {P}}_{\mathrm {N}}$ and $\hat {\mathbf {M}}{}^\dagger _{\mathrm {N}}$ commute; discrete duality symmetry.
• Reciprocal media (${\boldsymbol \varepsilon}^{\sf{T}}={\boldsymbol \varepsilon}$, ${\boldsymbol \xi}^{\sf{T}}=-{\boldsymbol \zeta}$, ${\boldsymbol \mu}^{\sf{T}}={\boldsymbol \mu}$): (15) holds; generalized real and imaginary parts reduce to ordinary ones; discrete duality symmetry.

Here, discrete duality symmetry means that the rotation angle is restricted to values θ = nπ/2 with $n\in \mathbb {Z}$. Note that duality is only realized as a continuous symmetry when Onsager-violation is allowed for. Notably, a reduction in reciprocity symmetry leads to an enhancement of duality symmetry.

In order to derive transformation laws for the Green tensor, we combine (9), (11), (28), (30) and (34) to write

$$\begin{eqnarray} &&\left(\begin{array}{@{}c@{}}\hat{{\bf E}}\\ Z_0\hat{{\bf H}}\end{array}\right) =-c\mathcal{B}\star \mathcal{G}\star \left(\begin{array}{@{}c@{}}Z_0\hat{{\bf P}}_{{\rm N}}\\ \mu_0\hat{{\bf M}}_{{\rm N}}\end{array}\right), \end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray} \mathcal{B} = \left(\begin{array}{@{}cc@{}}{\boldsymbol I} & {\boldsymbol 0} \\ -{\boldsymbol \mu}^{-1}\star {\boldsymbol \zeta} & {\boldsymbol \mu}^{-1} \end{array}\right)\!, \end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray} \mathcal{G}=\left(\begin{array}{@{}cc@{}}{\boldsymbol {\sf{G}}}_{ee} & {\boldsymbol {\sf{G}}}_{em}\\ {\boldsymbol {\sf{G}}}_{me} & {\boldsymbol {\sf{G}}}_{mm} +{\boldsymbol \mu} \end{array}\right)\!, \end{eqnarray} \tag{ 55 }$$

where we have introduced the shorthand notations $\boldsymbol{{\sf{G}}}_{ee}= (\mathrm{i}\,\omega/c)\boldsymbol{{\sf{G}}}(\mathrm{i}\,\omega/c)$, $\boldsymbol {{\sf {G}}}_{em}= (\mathrm {i}\,\omega /c) \boldsymbol {{\sf {G}}}\times \overleftarrow {\mathbf {\nabla }}{'}$, $\boldsymbol{{\sf{G}}}_{me}= \mathbf{\nabla}\times \boldsymbol{{\sf{G}}}(\mathrm{i}\,\omega/c)$ and $\boldsymbol {{\sf {G}}}_{mm}= \mathbf {\nabla }\times \boldsymbol {{\sf {G}}}\times \overleftarrow {\mathbf {\nabla }}{'}$. The transformed Green's tensors follow by applying duality transformations on both sides of this equation,

$$\begin{equation} \mathcal{G}^{{\circledast}} =\mathcal{B}^{{\circledast}-1}\star D(\theta) \mathcal{B}\star \mathcal{G}\star \mathcal{A}^{-1}\star D^{-1}(\theta) \mathcal{A}^{{\circledast}}, \end{equation} \tag{ 56 }$$

where

$$\begin{eqnarray} \mathcal{B}^{-1} =\left(\begin{array}{@{}cc@{}}{\boldsymbol I} & {\boldsymbol 0} \\ -{\boldsymbol \mu}^{-1}\star {\boldsymbol \zeta} & {\boldsymbol \mu}^{-1} \end{array}\right)\!. \end{eqnarray} \tag{ 57 }$$

In the special case of r and r' being in free space, we have $\mathcal {A},\mathcal {B}=\mathcal {I}$, so that $\mathcal {G}^{{\circledast }}=D(\theta )\mathcal {G}\star D^{-1}(\theta )$ and hence

$$\begin{equation} \left(\begin{array}{@{}c@{}}{\boldsymbol {\sf{G}}}_{ee}\\ {\boldsymbol {\sf{G}}}_{em} \\ {\boldsymbol {\sf{G}}}_{em}\\ {\boldsymbol {\sf{G}}}_{mm} + {\boldsymbol \delta}\end{array}\right)^{{\circledast}} = \mathcal{D} (\theta)\left(\begin{array}{@{}c@{}}{\boldsymbol {\sf{G}}}_{ee}\\ {\boldsymbol {\sf{G}}}_{em} \\ {\boldsymbol {\sf{G}}}_{em}\\ {\boldsymbol {\sf{G}}}_{mm} + {\boldsymbol \delta}\end{array}\right)\!. \end{equation} \tag{ 58 }$$

The Green tensors then transform like the medium response functions with (50).

Based on the general Ohm's law, we have quantized the electromagnetic field in the presence of nonlocal, nonreciprocal media which satisfies (i) the canonical commutation relations from free-space QED; (ii) the linear fluctuation–dissipation theorem; and (iii) the macroscopic Maxwell equations. A key feature of the scheme is the symmetrization of tensor fields via generalized real and imaginary parts, which is necessary whenever Onsager reciprocity does not hold. Their presence in the fluctuation–dissipation theorem is the key to avoiding restrictions on the allowed medium response.

For nonlocal and local bianisotropic media we have shown that quantization can alternatively be performed by the introduction of permittivity, permeability and magnetoelectric susceptibilities. When the latter do not vanish, the noise polarization and magnetization do not commute. We have explicitly determined the behaviour of the fields and response functions under duality transformations. The full continuous transformation group applies for bianisotropic and bi-isotropic media, but reduces to a discrete symmetry for isotropic, anisotropic and/or reciprocal media.

The scheme lays the foundation for exact studies of quantum phenomena such as dispersion forces, Förster energy transfer or environment-assisted molecular transition rates in the presence of motion or novel media with chiral or nonreciprocal properties. Moving media are a prime example for the occurrence of nonreciprocal material properties [28] which have to be thoroughly accounted for in order to understand, e.g. quantum friction. CP violation in atoms or molecules is manifest in their nonreciprocal cross-polarizability. The Curie principle, stating that certain interactions between two partners (atoms, molecules, bodies, etc) require them to possess similar properties, then allows for a detection of CP violation via atom–surface interactions provided the surface exhibits a corresponding nonreciprocity.

This work was supported by the UK Engineering and Physical Sciences Research Council (EPSRC).

To derive the integral relation (16) for the Green tensor, we write the Helmholtz equation (12) as

$$\begin{equation} \hat{{\boldsymbol {\sf{H}}}}\cdot\hat{{\boldsymbol {\sf{G}}}}=\hat{{\boldsymbol {\sf{I}}}}, \end{equation} \tag{ A.1 }$$

where $\langle \mathbf {r}|\hat {\boldsymbol {{\sf {G}}}}|\mathbf {r}'\rangle =\boldsymbol {{\sf {G}}}(\mathbf {r},\mathbf {r}',\omega )$ and $\langle \mathbf {r}|\hat {\boldsymbol {{\sf {H}}}}|\mathbf {r}'\rangle =[\mathbf {\nabla }\times \mathbf {\nabla }\times -\omega ^2/c^2] {\boldsymbol \delta }(\mathbf {r}-\mathbf {r}') -\mathrm {i}\,\mu _0\omega \boldsymbol {{\sf {Q}}}(\mathbf {r},\mathbf {r}',\omega )$. The Green operator is the right-inverse and, within any group of invertible operators, also the left-inverse of the Helmholtz operator,

$$\begin{equation} \hat{{\boldsymbol {\sf{G}}}}\cdot \hat{{\boldsymbol {\sf{H}}}}=\hat{{\boldsymbol {\sf{I}}}}. \end{equation} \tag{ A.2 }$$

From this relation and its Hermitian conjugate we find that

$$\begin{equation} \hat{{\boldsymbol {\sf{G}}}}\cdot \bigl(\hat{{\boldsymbol {\sf{H}}}}-\hat{{\boldsymbol {\sf{H}}}}{}^{\dagger}\bigr) \cdot \hat{{\boldsymbol {\sf{G}}}}{}^\dagger =\hat{{\boldsymbol {\sf{G}}}}{}^\dagger-\hat{{\boldsymbol {\sf{G}}}}. \end{equation} \tag{ A.3 }$$

In coordinate space this relation reads

$$\begin{eqnarray} &&\fl\mu_0\omega\int{\rm d}^3s\int{\rm d}^3s'\, {\boldsymbol {\sf{G}}}({\bf r},{\bf s},\omega)\,\cdot \, \mathcal{R}{\rm e}{\boldsymbol {\sf{Q}}}({\bf s},{\bf s}',\omega)\,\cdot \, {\boldsymbol {\sf{G}}}^{\dagger}({\bf r}',{\bf s}',\omega) =\mathcal{I}{\rm m}\ {\boldsymbol {\sf{G}}}({\bf r},{\bf r}',\omega) \end{eqnarray} \tag{ A.4 }$$

which in convolution notation takes the form of equation (16).