Maxwell Eigenmode approach to the Casimir-Lifshitz Torque

TL;DR

This paper presents a systematic eigenmode-based method to derive the Casimir-Lifshitz torque between birefringent media, extending previous results and enabling analysis of various anisotropic materials.

Contribution

It introduces a transparent eigenmode decomposition approach that generalizes Barash's calculation to include magnetic birefringence and finite thickness effects.

Findings

Derived a comprehensive formula for Casimir-Lifshitz torque in birefringent media.

Extended previous results to magnetic and electric birefringence.

Provides a framework for analyzing anisotropic materials with finite thickness.

Abstract

More than forty years ago, Barash published a calculation of the full retarded Casimir-Lifshitz torque for planar birefringent media with arbitrary degrees of anisotropy. An independent theoretical confirmation has been lacking since. We report a systematic and transparent derivation of the torque between two media with both electric and magnetic birefringence. Our approach, based on an eigenmode decomposition of Maxwell's equations, generalizes Barash's result for electrically birefringent materials, and can be generalized to a wide range of anisotropic materials and finite thickness effects.

Tables

Table 1. Table 1: Comparison between the conventions of this document and those of Ref. Barash ( 1978 ) . Unlike Ref. Barash ( 1978 ) we do not include the possibility of a liquid gap. On the other hand, we allow the option of magnetic anisotropy. Hence it is assumed here that μ j ​ x = 1 = μ j ​ y subscript 𝜇 𝑗 𝑥 1 subscript 𝜇 𝑗 𝑦 \mu_{jx}=1=\mu_{jy} . The z 𝑧 z -components of the wave vector are purely imaginary since they are evaluated at imaginary frequencies.

Ref. Barash (1978)	This document
$u=r\cos \phi$	$k_x=k_{\rho }\cos \theta$
$v=r\sin \phi$	$k_y=k_{\rho }\sin \theta$
$\phi$	$\theta$
$r$	$k_{\rho }$
$\theta$	$\phi$
${\rho }_1$	$i{q}_{1o}$
${\rho }_2$	$i{q}_{2o}$
${\rho }_1^{\prime }\left(\theta \right)$	$i{q}_{1e}\left(\phi \right)$
${\rho }_2^{\prime }\left(\theta \right)$	$i{q}_{2e}\left(\phi \right)$
${\rho }_3$	$i{k}_0$
${\epsilon }_{\mathrm{I1}}$	${\epsilon }_{1x}$
${\epsilon }_{\mathrm{I2}}$	${\epsilon }_{1y}$
${\epsilon }_{\mathrm{II1}}$	${\epsilon }_{2x}$
${\epsilon }_{\mathrm{II2}}$	${\epsilon }_{2y}$
${\epsilon }_3$	1
$l$	$a$

Equations

$$\underline{\mathbf{\epsilon }}(\omega )=\left(\begin{array}{ccc}{\epsilon }_{1x}{\cos }^2\theta +{\epsilon }_{1y}{\sin }^2\theta & ({\epsilon }_{1x}-{\epsilon }_{1y})\sin \theta \cos \theta & 0\\ ({\epsilon }_{1x}-{\epsilon }_{1y})\sin \theta \cos \theta & {\epsilon }_{1x}{\sin }^2\theta +{\epsilon }_{1y}{\cos }^2\theta & 0\\ 0& 0& {\epsilon }_{1y}\end{array}\right)$$

$$\underline{\mathbf{\mu }}(\omega )=\left(\begin{array}{ccc}{\mu }_{1x}{\cos }^2\theta +{\mu }_{1y}{\sin }^2\theta & ({\mu }_{1x}-{\mu }_{1y})\sin \theta \cos \theta & 0\\ ({\mu }_{1x}-{\mu }_{1y})\sin \theta \cos \theta & {\mu }_{1x}{\sin }^2\theta +{\mu }_{1y}{\cos }^2\theta & 0\\ 0& 0& {\mu }_{1y}\end{array}\right)$$

$$\mathbf{k}\times \mathbf{E}=\frac{\omega }{c}\mathbf{B},\mathbf{k}\times \mathbf{H}=-\frac{\omega }{c}\mathbf{D}$$

$$q\left(\begin{array}{c}{E}_x\\ E_y\\ H_x\\ H_y\end{array}\right)=\left(\begin{array}{cc}\underline{0}& {\underline{\mathbf{M}}}_a\\ {\underline{\mathbf{M}}}_b& \underline{0}\end{array}\right)\left(\begin{array}{c}{E}_x\\ E_y\\ H_x\\ H_y\end{array}\right),$$

$${\underline{\mathbf{M}}}_a=$$

$$\left(\begin{smallmatrix}\tfrac{\omega}{c}\sin\theta\cos\theta(\mu_{1x}-\mu_{1y})&{-}\color[rgb]{0,0,0}k_{\rho}^{2}\frac{c}{\varepsilon_{1y}\omega}{+}\color[rgb]{0,0,0}\tfrac{\omega}{c}(\mu_{1y}\cos^{2}\theta+\mu_{1x}\sin^{2}\theta)\\ {-}\color[rgb]{0,0,0}\tfrac{\omega}{c}(\mu_{1x}\cos^{2}\theta+\mu_{1y}\sin^{2}\theta)&\tfrac{\omega}{c}\sin\theta\cos\theta(\mu_{1y}-\mu_{1x})\end{smallmatrix}\right)$$

$${\underline{\mathbf{M}}}_b=$$

$$\left(\begin{smallmatrix}\tfrac{\omega}{c}\sin\theta\cos\theta(\varepsilon_{1y}-\varepsilon_{1x})&{-}\color[rgb]{0,0,0}\tfrac{\omega}{c}(\varepsilon_{1y}\cos^{2}\theta+\varepsilon_{1x}\sin^{2}\theta){+}\color[rgb]{0,0,0}k_{\rho}^{2}\tfrac{c}{\mu_{1y}\omega}\\ \tfrac{\omega}{c}(\varepsilon_{1x}\cos^{2}\theta+\varepsilon_{1y}\sin^{2}\theta)&\tfrac{\omega}{c}\sin\theta\cos\theta(\varepsilon_{1x}-\varepsilon_{1y}).\end{smallmatrix}\right),$$

$$\displaystyle\boldsymbol{\psi}_{o}=N_{o}^{-1}\left(\begin{smallmatrix}{-}\color[rgb]{0,0,0}\mu_{1y}q_{1o}\frac{\omega}{c}\sin\theta\\ {\mu_{1y}q_{1o}\frac{\omega}{c}\cos\theta}\\ k_{1z}^{2}\cos\theta\\ -\varepsilon_{1y}\mu_{1y}\frac{\omega^{2}}{c^{2}}\sin\theta\end{smallmatrix}\right),$$

$$\displaystyle\boldsymbol{\psi}_{e}=N_{e}^{-1}\left(\begin{smallmatrix}{-}\color[rgb]{0,0,0}k_{1z}^{2}\cos\theta\\ {-}\color[rgb]{0,0,0}\varepsilon_{1y}\mu_{1y}\frac{\omega^{2}}{c^{2}}\sin\theta\\ \varepsilon_{1y}q_{1e}\frac{\omega}{c}\sin\theta\\ -\varepsilon_{1y}q_{1e}\frac{\omega}{c}\cos\theta\end{smallmatrix}\right)$$

$$r_{1ssN}=k_0^2({\psi }_{e2}{\psi }_{o4}-{\psi }_{e4}{\psi }_{o2})+$$

$$k_0\frac{\omega }{c}({\psi }_{e1}{\psi }_{o2}-{\psi }_{e2}{\psi }_{o1}-{\psi }_{e3}{\psi }_{o4}+{\psi }_{e4}{\psi }_{o3})$$

$$+\frac{{\omega }^2}{c^2}({\psi }_{e3}{\psi }_{o1}-{\psi }_{e1}{\psi }_{o3}),$$

$$r_{1spN}=2k_0\frac{\omega }{c}({\psi }_{e1}{\psi }_{o4}-{\psi }_{e4}{\psi }_{o1}),$$

$$r_{1ppN}=k_0^2({\psi }_{e4}{\psi }_{o2}-{\psi }_{e2}{\psi }_{o4})+$$

$$k_0\frac{\omega }{c}({\psi }_{e1}{\psi }_{o2}-{\psi }_{e2}{\psi }_{o1}-{\psi }_{e3}{\psi }_{o4}+{\psi }_{e4}{\psi }_{o3})+$$

$$\frac{{\omega }^2}{c^2}({\psi }_{e1}{\psi }_{o3}-{\psi }_{e3}{\psi }_{o1}),$$

$$r_{1psN}=2k_0\frac{\omega }{c}({\psi }_{e2}{\psi }_{o3}-{\psi }_{e3}{\psi }_{o2}),$$

$$r_{1D}=k_0^2({\psi }_{e2}{\psi }_{o4}-{\psi }_{e4}{\psi }_{o2})$$

$$+k_0\frac{\omega }{c}({\psi }_{e1}{\psi }_{o2}-{\psi }_{e2}{\psi }_{o1}+{\psi }_{e3}{\psi }_{o4}-{\psi }_{e4}{\psi }_{o3})$$

$$+\frac{{\omega }^2}{c^2}({\psi }_{e1}{\psi }_{o3}-{\psi }_{e3}{\psi }_{o1}),$$

$$r_{1ssN}=\frac{{\epsilon }_{1y}{\mu }_{1y}{\omega }^2(k_0{\mu }_{1y}-q_{1e})(k_0{\epsilon }_{1y}+q_{1o}){\sin }^2\theta }{c^2}$$

$$-(k_0{q}_{1e}{\epsilon }_{1y}+k_{1z}^2)(k_{1z}^2-k_0{\mu }_{1y}{q}_{1o}){\cos }^2\theta$$

$$r_{1spN}=\frac{k_0{\epsilon }_{1y}{\mu }_{1y}\omega (k_{1z}^2-q_{1e}{q}_{1o})\sin 2\theta }{c}=r_{1psN}$$

$$r_{1ppN}=(k_{1z}^2-k_0{q}_{1e}{\epsilon }_{1y})(k_0{\mu }_{1y}{q}_{1o}+k_{1z}^2){\cos }^2\theta$$

$$+\frac{{\epsilon }_{1y}{\mu }_{1y}{\omega }^2{\sin }^2\theta (k_0{\mu }_{1y}+q_{1e})(q_{1o}-k_0{\epsilon }_{1y})}{c^2}$$

$$r_{1D}=(k_0{q}_{1e}{\epsilon }_{1y}+k_{1z}^2)(k_0{\mu }_{1y}{q}_{1o}+k_{1z}^2){\cos }^2\theta$$

$$+\frac{{\epsilon }_{1y}{\mu }_{1y}{\omega }^2{\sin }^2\theta (k_0{\mu }_{1y}+q_{1e})(k_0{\epsilon }_{1y}+q_{1o})}{c^2}.$$

$$\begin{array}{rl}\frac{E_{Cas}(a,\phi )}{A}& =\frac{k_{b}T}{4{\pi }^2}\sum _{n=0}^{\infty }(1-\frac{1}{2}{\delta }_{n,0})\times \\ & {\int }_0^{\infty }{\int }_0^{2\pi }\log (D(a,\mathbf{k},\theta ,\phi ,i{\zeta }_n))k_{\rho }\mathrm{d}{k}_{\rho }\mathrm{d}\theta \end{array}$$

$$\begin{split}&D(a,\mathbf{k},\theta,\varphi,i\zeta_{n})=\\ &\det\Big{(}\mathbf{\underline{I}}-\mathbf{\underline{r}_{1}}(\mathbf{k},\theta,i\zeta_{n})\cdot\mathbf{\underline{r}_{2}}(\mathbf{k},\theta+\varphi,i\zeta_{n})e^{-2k_{0}a}\Big{)}\end{split}$$

$$\tau (a,\phi )=-\frac{\partial E_{Cas}}{\partial \phi }.$$

$$\begin{array}{rl}& D(\mathbf{k},\phi ,a)\equiv \\ & 1+P(\mathbf{k},\phi )\exp (-2k_{0}a)+Q(\mathbf{k},\phi )\exp (-4k_{0}a).\end{array}$$

$$\begin{array}{rl}& r_{1D}{r}_{2D}P=\\ & -r_{1ppN}{r}_{2ppN}-2r_{1spN}{r}_{2spN}-r_{1ssN}{r}_{2ssN},\end{array}$$

$$\begin{array}{rl}& r_{1D}^2{r}_{2D}^{2}Q=\\ & (r_{1spN}^2-r_{1ppN}{r}_{1ssN})(r_{2spN}^2-r_{2ppN}{r}_{2ssN}).\end{array}$$

$$r_{1D}{r}_{2D}=(q_{1o}^2+k_{\rho }^2{\sin }^2\theta )(q_{2o}^2+k_{\rho }^2{\sin }^2(\theta +\phi ))\gamma$$

$$q_{jo}^2+k_{\rho }^2={\epsilon }_{jy}(k_0^2+k_{\rho }^2)j=1,2.$$

$$D_b=\frac{A_1-\frac{{\epsilon }_{2y}(i{q}_{2e}-i{q}_{2o})(B{k}_{\rho }^2{\sin }^2(\theta +\phi )+C-E(2k_{\rho }^2\sin (\theta )\cos (\phi )\sin (\theta +\phi )-k_0^2{\sin }^2(\phi )))}{-k_{\rho }^2{\sin }^2(\theta +\phi )-q_{2o}^2}}{\gamma }$$

$$A_1=((i{k}_0+i{q}_{2o})(i{k}_0+i{q}_{1o})-e^{-2a{k}_0}(i{q}_{2o}-i{k}_0)(i{q}_{1o}-i{k}_0))\times$$

$$((i{k}_0{\epsilon }_{2y}+i{q}_{2o})(i{k}_0{\epsilon }_{1y}+i{q}_{1o})-e^{-2a{k}_0}(i{q}_{2o}-i{k}_0{\epsilon }_{2y})(i{q}_{1o}-i{k}_0{\epsilon }_{1y}))$$

Full text

I Introduction

Casimir-Lifshitz forces Casimir (1948) are dispersion interactions between macroscopic bodies that arise from quantum mechanical and thermal fluctuations in the electromagnetic field. These forces, which can be considered a generalization of van der Waals forces to include finite light speed, depend on the electric and magnetic susceptibilities of the materials involved  Lifshitz (1955); *Lifshitz61 (cf. also Ref. Woods et al. (2016) for a review). Several decades after its first theoretical prediction, measuring the Casimir force directly became technologically feasible Lamoreaux (1997). Since this interaction is mediated by virtual and thermal photons, the frequency of which cannot be controlled directly, the Casimir force is a broadband phenomenon. In particular, as a consequence of the fluctuation-dissipation theorem, the frequency ranges where the susceptibilities change significantly actually provide a dominant contribution to the Casimir force.

From a fundamental viewpoint, the Casimir force plays a role in micron range gravitation experiments, and the search for deviations from Newtonian gravitation due to hypothetical new forces Palasantzas et al. (2015). This fuels the desire to come to precise comparisons between theoretical predictions and experimental data Almasi et al. (2015). More practically, Casimir interactions affect the actuation dynamics of nano- and micro-mechanical systems, such as switches, cantilevers, and actuators at a sub-micrometer length scale Serry et al. (1995); Chan et al. (2001a); *ChanScience2001; Miri and Golestanian (2008); Broer et al. (2015); Svetovoy et al. (2017); Klimchitskaya et al. (2018); Tajik et al. (2018); Ahn et al. (2018).

In order to calculate the Casimir potential, the Maxwell equations must be solved for the given geometry Lifshitz (1955); *Lifshitz61. Here we focus on the case of planar media where this can be done analytically. Anisotropy in the plane of reflection creates a dielectric contrast in the azimuthal direction, which gives rise to a Casimir torque Barash (1978) (see Fig. 1). Several experimental setups to detect the torque have been proposed Munday et al. (2005); Rodrigues et al. (2008); Guérout et al. (2015); Xu and Li (2017), but only recently has this phenomenon been observed experimentally Somers et al. (2018).

An exact analytical description for two planar birefringent half spaces was derived more than forty years ago by Yu. Barash Barash (1978). More recently, an alternative calculation has been presented Philbin and Leonhardt (2008), the result of which looks symbolically different. It has not been established that this result agrees analytically with that of Ref. Barash (1978), though the authors claim to have verified it numerically.

Recently, another formula was proposed for the Casimir-Lifshitz torque in appendix A of Ref.  Somers and Munday (2017a), where the result of of Ref. Lekner (1991) is inserted into the Lifshitz formula. Lambrecht et al. (2006) However, this does not include a supporting calculation that shows that it is equivalent to the result of Ref. Barash (1978). Up to the present, a transparent and independent calculation which analytically verifies the result of Barash has still been lacking.

The recent development of the experimental observation of the Casimir torque Somers et al. (2018) will instigate more investigations on this phenomenon, which necessitates a systematic and transparent formalism to describe it theoretically. The main obstacle is the failure of the usual decomposition into perpendicular $s$-polarized (or transverse electric) and $p$- polarized (transverse magnetic) modes. The reason for this failure is the fact that these are not the solutions of the Maxwell equations. Therefore it stands to reason to determine what these solutions actually are instead. This can be done by formulating the Maxwell equations as an eigenvalue problem. The Maxwell eigenmode formalism is a well-established method in electrical and electronic engineering Berreman (1972); Yeh (1980) that has been designed specifically to tackle the problem of scattering electromagnetic waves on anisotropic media. Note that the $s$- and $p$-mode decomposition actually does work if the plane of anisotropy is perpendicular to the plane of reflectance (cf., e.g., Ref. Grushin et al. (2011)). However, in such a case the Casimir torque vanishes.

Some limiting cases simplify the result of Ref. Barash (1978) considerably. An example is the limit of weak anisotropy, or more precisely, the limit of relatively small deviations from in-plane anisotropy. Although a treatment of the full anisotropic case exists Rosa et al. (2008), this is still a popular approximation Esquivel-Sirvent et al. (2010); Somers and Munday (2017b); Thiyam et al. (2018). This is valid for certain natural anisotropic materials such as calcite or quartz, but there is no reason to assume this must hold generally. Examples of materials that are not weakly anisotropic include, but are not limited to, to cuprate superconductors, Romanowsky and Capasso (2008) liquid crystals, Jákli (2013) anisotropic metamaterials, Poddubny et al. (2013) and multiferroic materials Fiebig et al. (2016). Another common simplifying assumption is the non-retarded limit of van der Waals forces Lu and Podgornik (2016); Lu (2018). Depending on the material(s), this approximation should work at separation distances of the order of 10 nm. However, even at such short distances this approximation can fail Somers and Munday (2017b). Here, we would like to make the case that such approximations are unnecessary by presenting an exact and transparent generalization of the calculation in Ref. Barash (1978).

The approach that we present is subject to the same assumptions underlying the original Lifshitz theory. Firstly, it relies on a continuous medium approximation, which is valid for wavelengths larger than the interatomic distance. A second assumption is that the medium exhibits linear dielectric response. Indeed, conventional Lifshitz theory does not extend to media with nonlinear dielectric response behavior. Such an extension requires a non-trivial generalization of the fluctuation-dissipation theorem,  Soo and Krüger (2016); *Soo2018. which is based on linear response theory

II Maxwell eigenmodes

In the case of birefringent media, one can distinguish the so-called ‘ordinary’ and ‘extraordinary’ waves. The former propagates as if the medium were isotropic, whereas the propagation of the latter depends on the medium’s orientation Landau and Lifshitz (1963).

Let the material slab be oriented in such a way that the anisotropic plane is facing the surface, which is defined as the $x$-$y$-plane in the laboratory’s coordinate system. Furthermore, the magnetic and electric anisotropy axes are assumed to be identical. Hence the electric permittivity and the magnetic permeability are given by the following tensors:

[TABLE]

[TABLE]

where $\theta$ denotes the angle between material’s optic axis and the $x$-axis of laboratory’s coordinate system (cf. Fig. 1). It must be stressed that the entries of both $\underline{\boldsymbol{\varepsilon}}$ and $\underline{\boldsymbol{\mu}}$ depend on frequency, but that the argument will be suppressed from now on.

The vectorial Maxwell equations in Fourier space are given by

[TABLE]

where $\mathbf{B}=\underline{\boldsymbol{\mu}}\cdot\mathbf{H}$ and $\mathbf{D}=\underline{\boldsymbol{\varepsilon}}\cdot\mathbf{E}$. Eq. 1 is a system of six linear equations with six unknowns, four of which are independent. This leads to the following 4$\times$4 matrix equation:

[TABLE]

where the non-zero quadrants are given by

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

where $k_{\rho}$ denotes the radial component of the wavevector.

The (extra)ordinary mode eigenvectors of Eq. 2, characterized by the respective subscripts $e$ ad $o$, are

[TABLE]

where $N_{o,e}$ denote normalisation constants, $q_{1o}=\sqrt{\varepsilon_{1y}\mu_{1x}\omega^{2}/c^{2}-(\mu_{1x}/\mu_{1y})k^{2}_{\rho}\cos^{2}\theta-k^{2}_{\rho}\sin^{2}\theta}$ and $q_{1e}=\sqrt{\varepsilon_{1x}\mu_{1y}\omega^{2}/c^{2}-(\varepsilon_{1x}/\varepsilon_{1y})k_{\rho}^{2}\cos^{2}\theta-k_{\rho}^{2}\sin^{2}\theta}$ are the respective eigenvalues,111The eigenvalues obtained in this work and in Ref. Lekner (1991) do not agree with those obtained in Ref. Rosa et al. (2008). and we use the shorthand notation k_{1z}^{2}\equiv{\mu_{1y}}\color[rgb]{0,0,0}{\varepsilon_{1y}\omega^{2}/c^{2}-k_{\rho}^{2}}. Note that only forward propagating modes are considered here, characterized by a positive real part, because only a single interface is taken into account for simplicity. The eigenvectors Eq. 3 determine the ratio between the electromagnetic field components inside the anisotropic medium. They will be required to obtain the Fresnel reflection coefficients. We emphasize that Eq. 3 should not be conflated with the $s$- and $p$-polarized modes, which can be defined as the eigenmodes of Eq. 2 for the isotropic case. These modes are characterized by the conditions $\{E_{x}\neq 0,E_{y}=0,H_{x}=0,H_{y}\neq 0\}$ and $\{E_{x}=0,E_{y}\neq 0,H_{x}\neq 0,H_{y}=0\}$, respectively. Clearly these conditions do not hold in this case. However, $s$- and $p$-mode decomposition to calculate the Casimir torque can still be found in the literature.

III Fresnel Reflection Coefficients

The anisotropic planar medium is faced with vacuum. As in any isotropic medium, in vacuum the electromagnetic waves can be written as a linear combination of $s$-polarized and $p$-polarized modes. This applies to both the incoming and the reflected waves. Inside the anisotropic medium, the fields can be written as a linear combination of the ordinary and extraordinary waves. The coefficients of these linear combinations are the Fresnel reflection and transmission coefficients. The condition of continuity of the tangential electromagnetic field components ensures that these coefficients are uniquely determined.

First, let the incoming wave be an $s$-wave, $\boldsymbol{\psi}_{i,s}$. Then this wave will be reflected as a linear combination between $s$- and $p$-modes, denoted by $\boldsymbol{\psi}_{rs}$. Hence the total wave in the isotropic media is $\boldsymbol{\psi}_{i,s}+\boldsymbol{\psi}_{rs}$. The first term is associated with an eigenvalue with a positive real part, whereas the eigenvalue of the latter mode has a negative real part.

Inside the anisotropic crystal, the electromagnetic wave can be written as a linear combination of ordinary and extraordinary waves. The ratio between the electromagnetic field components is fixed and the only degree of freedom is the proportionality constant for each mode. These constants are the transmission coefficients that couple the $s$-polarized wave to the ordinary and extraordinary waves: $\boldsymbol{\psi}_{ts}=E_{y,i}(t_{so}\boldsymbol{\psi}_{o}+t_{se}\boldsymbol{\psi}_{e})$ where the eigenvectors $\boldsymbol{\psi}_{o,e}$ are given by Eqs. (3), and the amplitude of the incoming wave $E_{y,i}>0$. Here the subscript $ts$ denotes the transmitted wave, that originates from an incoming $s$-wave, but is in itself not an $s$-wave.

The condition of continuity of the tangential electromagnetic field components at the interface implies that fields must be equal on both sides of the interface $\boldsymbol{\psi}_{i,s}+\boldsymbol{\psi}_{rs}=\boldsymbol{\psi}_{ts}$, which leads to the a system of four equations with four unknowns. The relevant solutions are denoted by $r_{ss}=\tfrac{r_{1ssN}}{r_{1D}}$ and $r_{sp}=\frac{r_{1spN}}{r_{1D}}$ and they can be obtained in terms of the eigenvector components $\psi_{ok}$ and $\psi_{ek}$ ($k=1..4$). The procedure is now repeated for an incoming $p$-polarized wave. This leads to another system of four equations with four unknowns. The relevant solutions of this system are denoted by $r_{pp}=\tfrac{r_{1ppN}}{r_{1D}}$ and $r_{ps}=\tfrac{r_{1psN}}{r_{1D}}$. This leads to the following entries of Fresnel reflection matrix:

[TABLE]

where $k_{0}=\sqrt{\omega^{2}/c^{2}-k_{\rho}^{2}}$ denotes the $z$-component of the wavevector in vacuum. Note that the reflection coefficients in terms of the Maxwell eigenvector components, Eq. 4 are not restricted to uniaxial materials, but are valid for biaxial half spaces as well. Of course, the explicit form of the eigenvector components will change in the biaxial case. In particular we note that the reflection matrix is symmetric if and only if $\psi_{e1}\psi_{o4}-\psi_{e4}\psi_{o1}=\psi_{e2}\psi_{o3}-\psi_{e3}\psi_{o2}$, which holds in the uniaxial case. Also note that the normalization constants cancel out of the reflection matrix.

Inserting the eigenvectors Eq. 3 into Eq. 4 yields the explicit reflection coefficients

[TABLE]

For the non-magnetic case, where $\mu_{1y}=1$ and $q_{1o}=k_{1z}$, the results Eq. 5 match Eqs. (34) and (42) of Ref. Lekner (1991).

The equivalent reflection coefficients for a non-identical second medium can be easily obtained by transforming Eq. 5 as follows: $q_{1e}\rightarrow{}q_{2e}$, $q_{1o}\rightarrow{}q_{2o}$, $k_{1z}\rightarrow k_{2z}$,$\varepsilon_{1y}\rightarrow\varepsilon_{2y}$, $\mu_{1y}\rightarrow\mu_{2y}$, and $\theta\rightarrow\theta+\varphi$, where $k_{2z}$, $q_{2e,o}$, $\varepsilon_{2y}$,$\mu_{2y}$ represent the equivalents of $q_{1e,o}$, $k_{1z}$, $\varepsilon_{1y}$, $\mu_{1y}$, respectively for medium 2, and $\varphi$ represents the angle between the optic axes of the media. (See Fig. 1). Next, the transformed eigenvectors must be inserted into Eq. 4. In general one must take into account the different propagation directions for each medium, but in this case, but in this case the reflection matrix is invariant under $\mathbf{k}\rightarrow-\mathbf{k}$.

IV Casimir Energy and Torque

Now we are in a position to determine the Casimir torque. The Casimir energy per unit area is given by the Lifshitz formula Lambrecht et al. (2006)

[TABLE]

with

[TABLE]

where $\mathbf{\underline{I}}$ denotes the 2$\times$2 identity matrix and $\mathbf{\underline{r}}_{j}$, $j=1,2$ represent the reflection matrices given by $\mathbf{\underline{r}}_{j}=\tfrac{1}{r_{jD}}\left(\begin{smallmatrix}r_{j,ss,N}&r_{j,sp,N}\\ r_{j,ps,N}&r_{j,pp,N}\end{smallmatrix}\right),$ the elements of which are given by Eq. 4 for $j=1$, and they should be transformed from medium 1 to 2 for $j=2$. For numerical convergence, all quantities are evaluated at the imaginary Matsubara frequencies $\zeta_{n}\equiv\frac{2\pi nk_{b}T}{\hbar}$ at finite temperature $T$. The Casimir torque is

[TABLE]

The form of Eq. 7 brings us to the apparent symbolic difference between this result and that of Ref. Philbin and Leonhardt (2008). First note that $D$ is a quadratic function of $\exp(-2k_{0}a)$ and let us introduce the notation

[TABLE]

If it is assumed that the reflection matrices are symmetric, (which they are in the birefringent case), the numerator of $P$ can be written as

[TABLE]

and that of $Q$ is

[TABLE]

The forms of Eqs. 10 and 11 are identical to those of Eqs. (55) and (56), respectively of Ref. Philbin and Leonhardt (2008). So at least in terms of the entries of the reflection matrices, there does not appear to be a symbolic difference between the result of this work and that of Ref. Philbin and Leonhardt (2008).

Next we will compare this result for $\mu_{1x}=1=\mu_{1y}$ to that of Ref. Barash (1978).

V Summary of proof of Barash’s formula

The question that arises now is: how does the result Eq. 8 for non-magnetic materials compare to Eq. (27) of Barash’s paper Barash (1978)? A detailed comparison can be found in the appendix. We will provide a summary here, omitting the algebraic details. We contend that Eq. 8 for non-magnetic materials is identical to Eq. (27) of Ref. Barash (1978). This was claimed in Ref. Somers and Munday (2017a) without proof. The basic idea of the proof is that this large, complicated problem is split into smaller, simpler and independent parts. This is illustrated by Fig. 2.

To address this problem, the first step is the realization that this comparison boils down to that between the arguments of the logarithms, i.e. $D$ of Eq. 7 to Eq. (21) of Ref. Barash (1978). Eq. (1) of Ref. Barash (1978), which defines the Casimir energy as the Helmholtz free energy, is identical to the Matsubara sum in the Lifshitz formula Eq. 6. For the sake of this comparison we introduce the subscripts $b$ and $L$, denoting results from Ref. Barash (1978) and this work respectively. Hence the comparison can be limited to that between $D_{b}$ and $D_{L}$.

The next step is the observation that both $D_{b}$ and $D_{L}$ are quadratic functions of $\exp(-2k_{0}a)$. (Cf. Eq. 9.) This reduces the comparison to that between the coefficients of $D_{b}$ and $D_{L}$.

This brings us to the constant (distance independent) term. This corresponds to the limit of large distances where $k_{0}a\rightarrow\infty$. It is expected that both $D_{b}$ and $D_{L}$ tend to unity as $k_{0}a\rightarrow\infty$, because the Casimir energy must tend to zero in this limit. From Eq. 7 it can be easily seen that $D_{L}\rightarrow 1$, but that $D_{b}\rightarrow 1$ is less obvious. The proof of the latter can be found in Section A.3 of the Appendix.

Here we continue to use the notation of Eq. 9, but now with the subscripts $b$ and $L$. It is important to realize that Barash wrote the coefficients in terms of the differences between the ordinary and extraordinary mode eigenvalues. More specifically, $P_{b}$ is a second degree polynomial in $q_{ie}-q_{io}$ with $i=1,2$ and $Q_{b}$ is a fourth degree polynomial in the same variables. Hence in order to come to a comparison, $P_{L}$ and $Q_{L}$ from Eq. 9 must be written in the same way. Let us introduce the following notation for the coefficients for both Barash’s and Lifshitz’s versions of $P$ and $Q$: $P_{b,L}=p_{0b,L}+p_{1,b,L}(q_{1e}-q_{1o})+p_{2b,L}(q_{2e}-q_{2o})+p_{3b,L}(q_{1e}-q_{1o})(q_{2e}-q_{2o})$ and $Q_{b,L}=q_{0b,L}+q_{1b,L}(q_{1e}-q_{1o})+q_{2b,L}(q_{2e}-q_{2o})+q_{3b,L}(q_{1e}-q_{1o})(q_{2e}-q_{2o})+q_{4b,L}(q_{1e}-q_{1o})^{2}+q_{5b,L}(q_{2e}-q_{2o})^{2}+q_{6b,L}(q_{1e}-q_{1o})(q_{2e}-q_{2o})^{2}+q_{7b,L}(q_{1e}-q_{1o})^{2}(q_{2e}-q_{2o})+q_{8b,L}(q_{1e}-q_{1o})^{2}(q_{2e}-q_{2o})^{2}$, where the comma denotes that one can choose between the subscripts $b$ or $L$.

Before moving on to the direct comparison between coefficients of $P_{L}$, $Q_{L}$ and $P_{b}$ and $Q_{b}$ respectively, it is worth noting that the labels 1 and 2 of the media and their respective angles $\theta$ and $\theta+\varphi$ are arbitrary and interchanging them should not change the physics. This symmetry leads to the useful relations Eqs. 17, 18, 19 and 20, which reduce the number of independent coefficients from 4 to 3 for $P_{b,L}$ and from 9 to 6 for $Q_{b,L}$. (See Fig. 2.)

Another step is now to rid Barash’s expression of the fractions within its numerator and denominator. This is done by multiplying numerator and denominator of both $P_{b}$ and $Q_{b}$ by the factor $(q_{1o}^{2}+k_{\rho}^{2}\sin^{2}\theta)(q_{2o}^{2}+k_{\rho}^{2}\sin^{2}(\theta+\varphi))$. Now, the denominators of $P_{b}$ and $P_{L}$ can be related as follows

[TABLE]

where $\gamma$, the denominator of $P_{b}$, is given by Eq. 15. Here we have taken advantage of the ordinary mode eigenvalue relations

[TABLE]

Note that Eq. 13 is simply the combination of the definitions of $q_{io}$ and $k_{0}$ with $\omega/c$ eliminated. With this elimination we follow Ref. Barash (1978).

Since Eq. 12 tells us how to relate the denominators of $P_{b}$ ad $P_{L}$, we are now ready to compare their numerators. The comparison is greatly simplified by assuming $q_{1e}=q_{1o}$ and $q_{2e}=q_{2o}$. This immediately establishes that $p_{0b}=p_{0L}$. We proceed with the assumption $q_{1e}\neq q_{1o}$ and $q_{2e}=q_{2o}$. The already established coefficient $p_{0L}$ is subtracted from the resulting expression for $P_{L}$ under this condition. This will yield and expression for $p_{1L}$, which is indeed identical to $p_{1b}$ if Eq. 13 is taken into account. The symmetry relation Eq. 17 allows us to skip $p_{2L}$, and we move on to $p_{3L}$. This requires the general case $q_{1e}\neq q_{1o}$ and $q_{2e}\neq q_{2o}$. From the full expression of $P_{L}$ the other terms that have been obtained so far, are subtracted, i.e. $P_{L}-(p_{0L}+p_{1L}(q_{1e}-q_{1o})+p_{2L}(q_{2e}-q_{2o}))$. The expression for $p_{3L}$ found in this way is equal to $p_{3b}$ with the eigenvalue relations Eq. 13. The detailed proof that $P_{b}=P_{L}$ can be found in Section A.5.

Now we can proceed with the coefficient $Q$. The denominator of $Q_{L}$ is the square of that of $P_{L}$. The numerator of $Q_{b}$ is multiplied with $(q_{1o}^{2}+k_{\rho}^{2}\sin^{2}\theta)(q_{2o}^{2}+k_{\rho}^{2}\sin^{2}(\theta+\varphi))\gamma$, so that both $Q_{b}$ and $Q_{L}$ have the same denominator. Consequently, $Q_{b}$ and $Q_{L}$ now have the number of coefficients. Most of these coefficients can be compared in the same way as those of $P_{b,L}$. However, for the mixed terms, the trick with the simplifying assumptions $q_{1e}=q_{1o}$ or $q_{2e}=q_{2o}$ no longer works. In this case the expressions for the coefficients $Q_{L}$ must be fully expanded in order to write them in the desired form. The full proof that $Q_{L}=Q_{b}$ can be found in Section A.6.

VI Conclusions

We have derived an exact and general expression for the Casimir torque between two half-spaces that exhibit both electric and magnetic birefringence, and thereby analytically verified the Casimir torque result of Ref. Barash (1978) as a limiting case of non-magnetically permeable media. This is the first complete and transparent analytic verification in the more than forty intervening years since the publication of the results of Ref. Barash (1978).

Our approach does not depend on approximations such as weak anisotropy or neglecting retardation effects, and can be generalized to other types of media, such as metamaterials Poddubny et al. (2013), multiferroics Fiebig et al. (2016), topological matter Tokura et al. (2019), biaxial materials Thiyam et al. (2018), and non-reciprocal materials Fuchs et al. (2017). The recently reported measurements of the Casimir torque Somers et al. (2018) should instigate further investigations for which this result is of considerable importance.

Acknowledgements.

WB thanks B. J. Hoenders for a useful email exchange. We acknowledge support from a start-up grant from Nanyang Technological University, under the number M4082095.110.

Appendix A Comparison to Barash’s result

A.1 Introduction

In this Appendix, we prove that the result obtained by Barash, Eq (21) in Barash (1978), is identical to the combination of the Lifshitz formula Lambrecht et al. (2006) with Lekner’s reflection coefficients Lekner (1991), i.e. Eq. 7 for the case $\mu_{x}=1=\mu_{y}$. To this end, it suffices to compare the argument of the logarithm in Eq. 7 , which Barash calls ‘dispersion equation’, denoted here by $D$. Barash (1978) After all, both results are written in the same form and the only thing that could be different is $D$.

Since the expression for $D$ can be quite complicated in both cases, it helps to split them into smaller, simpler parts. Firstly, note that it is a quadratic function of $\exp(-2k_{0}a)$: $D(\mathbf{k},\varphi,a)=1+P(\mathbf{k},\varphi)\exp(-2k_{0}a)+Q(\mathbf{k},\varphi)\exp(-4k_{0}a)$. (Part of this Appendix will be dedicated to showing that the constant term of Barash’s result equals unity indeed.) Next, the proof takes advantage of the fact that both $P$ and $Q$ are both polynomials in the two variables $q_{1e}-q_{1o}$ and $q_{2e}-q_{2o}$. The coefficients $P$ and $Q$ will be compared, by comparing the coefficients of these polynomials.

The notation of Ref.Barash (1978) is somewhat unconventional. For reference we have included a table (Table 1) to clarify how it compares to our notation.

The Appendix is organized as follows. After the introduction, both results will simply be given. Then it will be shown that Barash’s version of $D$ tends to one at sufficiently large distances. Next, the symmetry of between the coefficients will be discussed, thereby reducing the number of independent coefficients. Finally, the coefficients of both $P$ and $Q$ will be compared.

A.2 Barash vs. Lifshitz-Lekner’s result

In the notation of this document, Barash’s result is

[TABLE]

where the subscript $b$ denotes that this is a result by Barash. The coefficients are given by

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

[TABLE]

The argument of the logarithm in Lifshitz’ formula Eq. 7 in terms of the entries of the reflection matrix is

[TABLE]

[TABLE]

[TABLE]

where the subscript $L$ denotes that this is a combination of results by Lifshitz Lifshitz (1955); Dzyaloshinskii et al. (1961) and Lekner. Lekner (1991)

The entries of the reflection matrix from Eq. (7) from the main paper for $\mu_{y}=1$ and $k_{z}=q_{1o}$ are rewritten as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where we have used the relations $k_{0}=\sqrt{\omega^{2}/c^{2}-k_{\rho}^{2}}$ and $q_{1o}=\sqrt{\varepsilon_{1y}\omega^{2}/c^{2}-k_{\rho}^{2}}$, the latter of which holds only in the non-magnetic case. To obtain the analogous expressions for the second medium, the substitutions from from one medium to the other must be performed.

A.3 Barash’s result in the limit of large distances

At sufficiently large distances, the exponentially decreasing terms vanish and only the constant term remains. The coefficients $A$ to $C$ in Eq. 14 are in this limit

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

while $E_{\lim}=0$.

Now the numerator of $D_{b}$ for $k_{0}a\gg 1$ is:

[TABLE]

[TABLE]

[TABLE]

To obtain $D_{b}$ in this limit, this quantity must be divided by $\gamma$ from Eq. 15, which does not depend on $a$. It can be seen that this numerator is actually the same as $\gamma$, hence

[TABLE]

Since $D_{b}$ is the argument of the logarithm in the Lifshitz formula Eq. 7 the Casimir energy tends to zero as $a$ goes to infinity.

A.4 Symmetry relations between coefficients

We take advantage of the fact that $P_{b}$ is a second degree polynomial in $q_{1e}-q_{1o}$ and $q_{2e}-q_{2o}$, and that $Q_{b}$ is a fourth degree polynomial in the same variables. Generally a polynomial of degree $d$ with $n$ variables has $\binom{n+d}{d}=\binom{n+d}{n}$ coefficients. Hence in our case $P_{b}$ has 6 coefficients and $Q_{b}$ has 15. However, the number of nonzero coefficients are 4 and 9 respectively. Moreover, these coefficients are not independent of each other. Note that the labels 1 and 2 of the media are arbitrary, so switching them should not affect the torque. The constant term does not contribute to the torque and it can be ignored. Since medium 1 is associated with an optic axis with angle $\theta$ and medium 2 is associated with $\theta+\varphi$, these angles must be interchanged as well. Hence we define the following transformation

[TABLE]

[TABLE]

where ‘$1\leftrightarrow 2$’ denotes that the subscript 1 needs to be replaced by 2 and vice versa. In other words, the transformations of from one medium to the other and their reverse have to be performed simultaneously. The condition that the torque must not change leads to

[TABLE]

[TABLE]

Now let

[TABLE]

Then we come to the following symmetry relations for the coefficients of $P$:

[TABLE]

[TABLE]

[TABLE]

And similarly for $Q$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

we have

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The symmetry relations should be valid for both the Barash and the Lifshitz-Lekner versions of these coefficients. In the latter case, this is immediately seen from Eqs. 21 and 30.

In the former case it is not so obvious that this symmetry holds. Therefore this needs to be checked. Switching the labels on Barash’s version of $D$ yields:

[TABLE]

where the primed quantities denote the transformed version of the unprimed quantities. In order to have $D^{\prime}=D$ we require $\gamma^{\prime}=\gamma$, which is relatively easily seen to hold.

Of course we are especially interested in symmetry relations between coefficients with different subscripts, because these could simplify some of the calculations. In particular we note that they reduce the number of independent coefficients from 4 to 3 for $P$ and from 9 to 6 for $Q$.

A.5 The coefficient $P$

Next, let us focus on the part of $D_{L}$ that is proportional to $\exp(-2k_{0}a)$, denoted by $P_{L}$:

[TABLE]

where $r_{psN1}r_{spN2}+r_{psN2}r_{spN1}$ can be simplified to $2r_{spN1}r_{spN2}$. Barash’s equivalent of this coefficient is

[TABLE]

[TABLE]

where $P_{A}$ to $P_{E}$ denote the respective parts of $A_{1}$ to $E$ proportional to $\exp(-2_{0}a)$.

First we will concentrate on the denominators of both expressions. The denominator of $P_{L}$ is

[TABLE]

[TABLE]

[TABLE]

In order to compare this to Barash’s expression, the numerator and denominator of the latter must be multiplied by $(q_{1o}^{2}+k_{\rho}^{2}\sin^{2}\theta)(q_{2o}^{2}+k_{\rho}^{2}\sin^{2}(\theta+\varphi))$. The denominator of $P_{b}$ is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which is identical to (23) due to the eigenvalues

[TABLE]

[TABLE]

The numerator of Eq. 22 is far more complicated than its denominator. However it is considerably simplified by the assumptions $q_{1e}=q_{1o}$ and $q_{2e}=q_{2o}$. First we will prove that the numerator of Eq. 21 is identical to that of Eq. 22 multiplied by $(q_{1o}^{2}+k_{\rho}^{2}\sin^{2}\theta)(q_{2o}^{2}+k_{\rho}^{2}\sin^{2}(\theta+\varphi))$ under these conditions. Later, these assumptions will be relaxed. If $q_{1e}=q_{1o}$ and $q_{2e}=q_{2o}$, the numerator of Eq. 22 will simplify to

[TABLE]

[TABLE]

[TABLE]

In this case the reflection matrices are both diagonal, hence the numerator of Eq. 21 becomes

[TABLE]

[TABLE]

which is equal to Eq. 25.

Now let $q_{1e}\neq q_{1o}$, while $q_{2o}=q_{2e}$. Now the numerator of $P_{b}$ is

[TABLE]

In this case only the term proportional to $(q_{1e}-q_{1o})$ needs to be considered, since we have already shown that the other terms are identical in the previous paragraph. Hence we are left with

[TABLE]

[TABLE]

Since $q_{2o}=q_{2e}$, $r_{spN2}=0=r_{psN2}$. Hence for $P_{L}$ the relevant term is

[TABLE]

[TABLE]

which can be simplified to

[TABLE]

[TABLE]

This expression is identical to Eq. 26 if Eq. 24 is taken into account.

Now let us make the reverse assumption, namely that $q_{1e}=q_{1o}$, but $q_{2e}\neq q_{2o}$. By the same token as before, we will focus on the term proportional to $q_{2e}-q_{2o}$, $p_{2}(q_{2e}-q_{2o})$. Because of the symmetry relation Eq. 17 it follows that Barash’s and Lifshitz-Lekner’s versions of this expression should be identical as well. Nonetheless we will check this here:

[TABLE]

[TABLE]

Here $P_{E}$ is omitted since it is proportional to $(q_{1e}-q_{1o})$. Now for $P_{L}$ we must evaluate

[TABLE]

[TABLE]

[TABLE]

which is identical to Eq. 27 under the conditions of Eq. 24. Note that the third factor in the second term of Eq. (24) of Ref. Barash (1978) should be $(r^{2}+\rho_{1}\rho_{2})$, i.e. with a relative plus sign rather than a minus sign. This has been confirmed by Ref. Munday et al. (2005), and now also here.

Finally, what remains is the term proportional to the product of the differences between the eigenvalues, $(q_{1e}-q_{1o})(q_{2e}-q_{2o})$. According to Barash, this is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The Lifshitz-Lekner equivalent of this expression is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which is identical to Eq. 28 if Eq. 24 is considered. This completes the proof that

[TABLE]

A.6 The Coefficient $Q$

Finally we direct our attention towards the terms proportional to $\exp(-4k_{0}a)$, denoted by $Q_{L}$:

[TABLE]

[TABLE]

In analogy to the previous subsection, the same coefficient according to Barash is

[TABLE]

[TABLE]

where $Q_{A}$ to $Q_{C}$ denote the respective parts of $A_{1}$ to $C$ proportional to $\exp(-4k_{0}a)$. $Q_{E}$ is omitted since $E$ is proportional to $\exp(-2k_{0}a)$ only.

The denominator of $Q_{L}$ is actually the square of that of $P_{L}$:

[TABLE]

Hence in order to compare $Q_{L}$ to $Q_{b}$, the numerator of the latter must be multiplied by a factor of $\gamma(q_{1o}^{2}+k_{\rho}^{2}\sin^{2}\theta)^{2}(q_{2o}^{2}+k_{\rho}^{2}\sin^{2}(\theta+\varphi))^{2}$. After all, the denominator of $Q_{b}$ is simply $\gamma$.

As before, we will start with the simplest case: $q_{1e}=q_{1o},q_{2e}=q_{2o}$. The numerator of $Q_{L}$ will simplify to

[TABLE]

[TABLE]

of which Barash’s equivalent is

[TABLE]

[TABLE]

which is clearly identical to the Lifshitz-Lekner expression.

Next we will assume again that $q_{2e}=q_{2o}$ and that $q_{1e}\neq q_{1o}$. However, contrary to the previous subsection, now the numerator of $Q$ is a quadratic function of $q_{1e}-q_{1o}$. In the previous paragraph, it has been shown that the constant term in this case is identical according to both Lifshitz-Lekner and Barash. The latter expression in this case is

[TABLE]

[TABLE]

which must be be multiplied by

[TABLE]

[TABLE]

[TABLE]

Here we are concerned only with the terms proportional to $(q_{1e}-q_{1o})$ and $(q_{1e}-q_{1o})^{2}$ of the product between these expressions. The former is according to Barash:

[TABLE]

[TABLE]

which can be simplified to

[TABLE]

[TABLE]

and the latter is

[TABLE]

[TABLE]

The equivalent expressions according to Lifsitz-Lekner can be obtained by simplifying the factors $r_{ppN1}r_{ppN2}$, $r_{ssN1}r_{ssN2}$, $r_{spN1}r_{ssN2}$, and $-r_{ppN2}r_{psN1}$ first. From this we gather the terms proportional to $q_{1e}-q_{1o}$:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which is identical to Eq. 32 with the eigenvalues Eq. 24. So

[TABLE]

The terms proportional to $(q_{1e}-q_{1o})^{2}$ are

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which is identical to Eq. 33 with the eigenvalues Eq. 24. Hence

[TABLE]

We will now assume again that $q_{1e}=q_{1o}$ and that $q_{2e}\neq q_{2o}$, and focus on the terms proportional to $q_{2e}-q_{2o}$ and $(q_{2e}-q_{2o})^{2}$. The term proportional to $(q_{2e}-q_{2o})^{2}$, according to Barash, can be simplified to

[TABLE]

[TABLE]

According to Lekner-Lifshitz this term is

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which is the same as Eq. 36 with the eigenvalue relations Eq. 24. In other words

[TABLE]

The symmetry relation Eq. 18 shows that the Lekner-Lifshitz version of the term proportional to $q_{2o}-q_{2e}$ is identical to its Barash’s version:

[TABLE]

Finally, we arrive at the mixed terms, i.e. those proportional to $(q_{1e}-q_{1o})^{j}(q_{2o}-q_{2e})^{k}$ with $j,k\in\{1,2\}$. There are in total four such terms, three of which are independent. First we will establish the Barash (denoted by subscript $b$ and Lifshitz-Lekner (subscript $L$) variants of these coefficients independently.

Let us start with the Lifshitz variant. This has the advantage that each coefficient can be written as a product of terms associated with medium 1 and those associated with medium 2. This makes it possible to calculate only the terms associated with the first medium, and then multiply that with a similar expression corresponding to the other medium. The part of $Q_{L}$ associated with medium 1, denoted by $Q_{1L}$, can be written as follows:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

which has to multiplied with the same expression transformed to medium 2 to obtain the total $Q_{L}$. From now on we will limit ourselves to the part of $Q_{L}$ corresponding to medium 1 only. The goal is now to write this expression as a second order polynomial in $q_{1e}-q_{1o}$. Then the total $Q_{L}$ is a product of this polynomial with the corresponding polynomial for medium 2. For this, we gather the coefficients $q_{0L}$ to $q_{8L}$. The first step is to expand $Q_{1L}$ completely as a polynomial in $q_{1e}$. This leads to the rather lengthy expression, so it is more insightful to write this implicitly in terms of the dummy variable coefficients $q^{\prime}_{1i}$, (where the first subscript labels the medium, and the second one the order) i.e.

[TABLE]

which can be transformed into a polynomial in $q_{1e}-q_{1o}$ as

[TABLE]

which leads to the following relations between the dummy coefficients:

[TABLE]

[TABLE]

[TABLE]

Now the same procedure is repeated for medium 2, leading to the analogous coefficients $q_{20}$, $q_{21}$, and $q_{22}$. Then the total $Q_{L}$ is written as

[TABLE]

which leads to the following relevant Lifshitz coefficients, now written explicitly:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The reason why we consider $q_{7L}$ and not $q_{6L}$ will become apparent later.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Next we will move on to the Barash variant of these coefficients, i.e. $q_{3b}$, $q_{7b}$, and $q_{8b}$. Unfortunately these are not so straightforwardly determined, since Barash’s result is not separated into a product of terms associated with medium 1 and medium 2. Rather, $Q_{b}=Q_{1b}+Q_{2b}$ is split into two terms as in Eq. 31, where both terms contain a mixture of expressions depending on both media. The first term of Eq. 31 is a quadratic polynomial in $q_{1e}-q_{1o}$ and $q_{2e}-q_{2o}$:

[TABLE]

[TABLE]

for some dummy coefficients $d_{i}$. So this term contributes only to $q_{3b}$ and $q_{6b}$, and the other mixed terms must be extracted from the second term

[TABLE]

[TABLE]

[TABLE]

Now it becomes apparent why we considered $q_{7L}$ in Eq. 40: it is more convenient since $q_{7b}$ is contained in one term, unlike $q_{6b}$. The same applies to $q_{8b}$, which is the most convenient to start with, as it is the highest order term. Simply multiplying the terms $~{}(q_{1e}-q_{1o})$ and $~{}(q_{2e}-q_{2o})$ within the parentheses yields

[TABLE]

[TABLE]

[TABLE]

which is identical to Eq. 41 under the condition of Eq. 24. Hence

[TABLE]

Next we will direct our attention to the coefficient $q_{7b}$. This one is a bit harder to determine. It is useful to remember that $Q_{2b}$ is essentially a polynomial of the form

[TABLE]

where $x=q_{1e}-q_{1o}$ and $y=q_{2e}-q_{2o}$ and the coefficients $a_{i}$ are complicated expressions. Here $q_{7b}$ would be the coefficient of this polynomial proportional to $xy^{2}$, which is

[TABLE]

More explicitly, in our case this coefficient is

[TABLE]

[TABLE]

which identical to Eq. 40 under the conditions of Eq. 24. Hence

[TABLE]

Due to the symmetry relation Eq. 20 the same equality can be claimed for the sixth coefficients:

[TABLE]

Finally, we consider the most difficult coefficient to obtain, $q_{3b}$. In accordance with the previous analogy, we can write the total $Q_{b}$ as

[TABLE]

Here $q_{3b}$ would be the coefficient multiplied with $xy$, which is

[TABLE]

Explicitly $q_{3b}$ is then, after simplification,

[TABLE]

[TABLE]

which is identical to Eq. 39 under conditions of Eq. 24. Hence

[TABLE]

In this subsection we have taken advantage of the symmetry relations Eq. 18 to obtain Eq. 38 and Eq. 20 to obtain Eq. 46. So the natural next step is to prove that Eqs. 18 and 20 hold indeed. This is simply a matter of determining the relevant coefficients, starting with

[TABLE]

[TABLE]

which is the transformed version of $q_{1b}$ from Eq. 32. Now finally we move on to the remaining coefficient

[TABLE]

[TABLE]

[TABLE]

which is Eq. 44 with the media swapped.

This completes the proof that

[TABLE]

Since $P_{b}=P_{L}$ (see Section A.5) and the constant term of $D_{b}$ equals unity (see Section A.3), this also proves that the result Eq. 7 for $\mu_{x}=1=\mu_{y}$ is identical to Eq. (21) of Ref. Barash (1978). To our knowledge this is the first independent proof of this result.