Thermal spin photonics in the near-field of nonreciprocal media

TL;DR

This paper explores how nonreciprocal materials can exhibit persistent thermal photon spin and heat currents at equilibrium, revealing new phenomena in near-field thermal radiation linked to spin-momentum locking and surface polaritons.

Contribution

It introduces the concepts of persistent thermal photon spin and planar heat current in nonreciprocal media, supported by theoretical analysis and a proposed imaging experiment.

Findings

Nonzero photon spin and heat flux occur in nonreciprocal slabs at thermal equilibrium.

Spin-momentum locking of evanescent waves underpins these phenomena.

Surface polaritons exhibit a spin magnetic moment in nonreciprocal photonics.

Abstract

The interplay of spin angular momentum and thermal radiation is a frontier area of interest to nanophotonics as well as topological physics. Here, we show that a thick planar slab of a nonreciprocal material, despite being at thermal equilibrium with its environment, can exhibit nonzero photon spin angular momentum and nonzero radiative heat flux in its vicinity. We identify them as the persistent thermal photon spin (PTPS) and the persistent planar heat current (PPHC) respectively. With a practical example system, we reveal that the fundamental origin of these phenomena is connected to spin-momentum locking of thermally excited evanescent waves. We also discover spin magnetic moment of surface polaritons in nonreciprocal photonics that further clarifies these features. We then propose a novel thermal photonic imaging experiment based on Brownian motion that allows one to witness these surprising features by directly looking at them using a lab microscope. We further demonstrate the universal behavior of these near-field thermal radiation phenomena through a comprehensive analysis of gyroelectric, gyromagnetic and magneto-electric nonreciprocal materials. Together, these results expose a surprisingly little explored research area of thermal spin photonics with prospects for new avenues related to non-Hermitian topological photonics and radiative heat transport.

Equations

$$\mathbf{A}({\mathbf{r}}_1)={\int }_{V_{{\mathbf{r}}_2}}\overline{\overline{\mathbf{G}}}({\mathbf{r}}_1,{\mathbf{r}}_2){\mu }_0\mathbf{J}({\mathbf{r}}_2)d^3{\mathbf{r}}_2$$

$$\mathbf{E}({\mathbf{r}}_1)=i\omega {\mu }_0{\int }_{V_{{\mathbf{r}}_2}}\overline{\overline{\mathbf{G}}}({\mathbf{r}}_1,{\mathbf{r}}_2)\mathbf{J}({\mathbf{r}}_2)d^3{\mathbf{r}}_2={\omega }^2{\mu }_0{\int }_{V_{{\mathbf{r}}_2}}\overline{\overline{\mathbf{G}}}({\mathbf{r}}_1,{\mathbf{r}}_2)\mathbf{p}({\mathbf{r}}_2)d^3{\mathbf{r}}_2$$

$${\mathbf{E}}_0({\mathbf{R}}_1,z_1)={\omega }^2{\mu }_0\int \frac{d^2{\mathbf{k}}_{\parallel }}{(2\pi {)}^2}{e}^{i{\mathbf{k}}_{\parallel }\cdot ({\mathbf{R}}_1-{\mathbf{R}}_2)}\frac{i}{2k_z}{e}^{i{k}_z(z_1-z_2)}[{\hat{\mathbf{e}}}_{s+}({\hat{\mathbf{e}}}_{s+}\cdot \mathbf{d})+{\hat{\mathbf{e}}}_{p+}({\hat{\mathbf{e}}}_{p+}\cdot \mathbf{d})]$$

$${\hat{\mathbf{e}}}_{s\pm }=\left[\begin{array}{c}\sin \varphi \\ -\cos \varphi \\ 0\end{array}\right],{\hat{\mathbf{e}}}_{p\pm }=\frac{-1}{k_0}\left[\begin{array}{c}\pm k_z\cos \varphi \\ \pm k_z\sin \varphi \\ -k_{\parallel }\end{array}\right]$$

$${\mathbf{E}}_{\text{ref}}({\mathbf{R}}_1,z_1)={\omega }^2{\mu }_0\int \frac{d^2{\mathbf{k}}_{\parallel }}{(2\pi {)}^2}{e}^{i{\mathbf{k}}_{\parallel }\cdot ({\mathbf{R}}_1-{\mathbf{R}}_2)}\frac{i}{2k_z}{e}^{i{k}_z(z_1+z_2)}[(r_{ss}{\hat{\mathbf{e}}}_{s+}+r_{ps}{\hat{\mathbf{e}}}_{p+})({\hat{\mathbf{e}}}_{s-}\cdot \mathbf{d})+(r_{sp}{\hat{\mathbf{e}}}_{s+}+r_{pp}{\hat{\mathbf{e}}}_{p+})({\hat{\mathbf{e}}}_{p-}\cdot \mathbf{d})]$$

$$V(\mathbf{r})=-\sum _a{f}_a(\mathbf{r})x_a(\mathbf{r})$$

$$x_a(\omega ,\mathbf{r})=\int \sum _b{\alpha }_{ab}(\omega ;\mathbf{r},{\mathbf{r}}^{\prime })f_b(\omega ,{\mathbf{r}}^{\prime })d^3{\mathbf{r}}^{\prime }$$

$$\displaystyle\langle x_{a}(\mathbf{r},\omega)x_{b}^{*}(\mathbf{r}^{\prime},\omega^{\prime})\rangle=\frac{\alpha_{ab}(\omega;\mathbf{r},\mathbf{r}^{\prime})-\alpha_{ba}^{*}(\omega;\mathbf{r}^{\prime},\mathbf{r})}{2i}\frac{1}{\omega}\underbrace{\bigg{(}\frac{\hbar\omega}{2}+\frac{\hbar\omega}{\text{exp}[\hbar\omega/k_{B}T]-1}\bigg{)}\delta(\omega-\omega^{\prime})}_{\Theta(\omega,T)}$$

$$A_j(\mathbf{r},\omega )=\int G_{jk}(\omega ;\mathbf{r},{\mathbf{r}}^{\prime }){\mu }_0{J}_k(\omega ,{\mathbf{r}}^{\prime })d^3{\mathbf{r}}^{\prime }$$

$$⟨\mathbf{A}({\mathbf{r}}_1)\otimes {\mathbf{A}}^{\ast }({\mathbf{r}}_2)⟩=\frac{\overline{\overline{\mathbf{G}}}({\mathbf{r}}_1,{\mathbf{r}}_2)-\overline{\overline{\mathbf{G}}}({\mathbf{r}}_2,{\mathbf{r}}_1{)}^{{\ast }^T}}{2i}\frac{{\mu }_0}{\omega }\Theta (\omega ,T)$$

$$⟨\mathbf{A}({\mathbf{r}}_1)\otimes {\mathbf{A}}^{\ast }({\mathbf{r}}_2)⟩=\frac{{\mu }_0\Theta (\omega ,T)}{\omega }\int \frac{k_{\parallel }d{k}_{\parallel }d\varphi }{2i(2\pi {)}^2}$$

$$+\frac{i}{2k_z}{e}^{i{k}_z(z_1+z_2)}[(r_{ss}{\hat{\mathbf{e}}}_{s+}+r_{ps}{\hat{\mathbf{e}}}_{p+}){\hat{\mathbf{e}}}_{s-}^T+(r_{sp}{\hat{\mathbf{e}}}_{s+}+r_{pp}{\hat{\mathbf{e}}}_{p+}){\hat{\mathbf{e}}}_{p-}^T]$$

$$+\frac{i}{2k_z^{\ast }}{e}^{-i{k}_z^{\ast }(z_1-z_2)}[{\hat{\mathbf{e}}}_{s-}{\hat{\mathbf{e}}}_{s-}^T+{\hat{\mathbf{e}}}_{p-}{\hat{\mathbf{e}}}_{p-}^T{]}^{{\ast }^T}$$

$$\displaystyle+\frac{i}{2k_{z}^{*}}e^{-ik_{z}^{*}(z_{1}+z_{2})}[(r_{ss}\mathbf{\hat{e}}_{s+}+r_{ps}\mathbf{\hat{e}}_{p+})\mathbf{\hat{e}}_{s-}^{T}+(r_{sp}\mathbf{\hat{e}}_{s+}+r_{pp}\mathbf{\hat{e}}_{p+})\mathbf{\hat{e}}_{p-}^{T}]^{*^{T}}\bigg{]}$$

$$⟨\mathbf{E}({\mathbf{r}}_1)\otimes {\mathbf{E}}^{\ast }({\mathbf{r}}_2)⟩=⟨i\omega \mathbf{A}({\mathbf{r}}_1))(i\omega \mathbf{A}({\mathbf{r}}_2){)}^{{\ast }^T}⟩={\omega }^2⟨\mathbf{A}({\mathbf{r}}_1)\otimes {\mathbf{A}}^{\ast }({\mathbf{r}}_2)⟩$$

$$⟨\mathbf{E}({\mathbf{r}}_1,\omega )\otimes {\mathbf{H}}^{\ast }({\mathbf{r}}_2,\omega )⟩=\frac{1}{{\mu }_0}⟨(i\omega \mathbf{A}({\mathbf{r}}_1,\omega ))({\nabla }_{r_2}\times \mathbf{A}({\mathbf{r}}_2,\omega ){)}^{{\ast }^T}⟩=\frac{i\omega }{{\mu }_0}⟨\mathbf{A}({\mathbf{r}}_1)\otimes {\mathbf{A}}^{\ast }({\mathbf{r}}_2)⟩({\nabla }_{r2}\times {)}^T$$

$$⟨\mathbf{H}({\mathbf{r}}_1,\omega )\otimes {\mathbf{H}}^{\ast }({\mathbf{r}}_2,\omega )⟩=\frac{1}{{\mu }_0^2}⟨({\nabla }_{r_1}\times \mathbf{A}({\mathbf{r}}_1,\omega ))({\nabla }_{r_2}\times {\mathbf{A}}^{\ast }({\mathbf{r}}_2,\omega ){)}^T⟩=\frac{1}{{\mu }_0^2}({\nabla }_{r1}\times )⟨\mathbf{A}({\mathbf{r}}_1)\otimes {\mathbf{A}}^{\ast }({\mathbf{r}}_2)⟩({\nabla }_{r2}\times {)}^T$$

$${\nabla }_{r2}\times =\left[\begin{array}{ccc}0& -{\partial }_z& {\partial }_y\\ {\partial }_z& 0& -{\partial }_x\\ -{\partial }_y& {\partial }_x& 0\end{array}\right]e^{i{k}_x(x_1-x_2)+i{k}_y(y_1-y_2)-i{k}_z^{\ast }(z_1-z_2)}=i\left[\begin{array}{ccc}0& -k_z^{\ast }& -k_y\\ k_z^{\ast }& 0& k_x\\ k_y& -k_x& 0\end{array}\right]$$

$$⟨\mathbf{E}({\mathbf{r}}_1)\otimes {\mathbf{E}}^{\ast }({\mathbf{r}}_2)⟩={\mu }_0\omega \Theta (\omega ,T)\int$$

$$+\frac{i}{2k_z}{e}^{i{k}_z(z_1+z_2)}[(r_{ss}{\hat{\mathbf{e}}}_{s+}+r_{ps}{\hat{\mathbf{e}}}_{p+}){\hat{\mathbf{e}}}_{s-}^T+(r_{sp}{\hat{\mathbf{e}}}_{s+}+r_{pp}{\hat{\mathbf{e}}}_{p+}){\hat{\mathbf{e}}}_{p-}^T]$$

$$+\frac{i}{2k_z^{\ast }}{e}^{-i{k}_z^{\ast }(z_1-z_2)}[{\hat{\mathbf{e}}}_{s-}{\hat{\mathbf{e}}}_{s-}^T+{\hat{\mathbf{e}}}_{p-}{\hat{\mathbf{e}}}_{p-}^T{]}^{{\ast }^T}$$

$$\displaystyle+\frac{i}{2k_{z}^{*}}e^{-ik_{z}^{*}(z_{1}+z_{2})}[(r_{ss}\mathbf{\hat{e}}_{s+}+r_{ps}\mathbf{\hat{e}}_{p+})\mathbf{\hat{e}}_{s-}^{T}+(r_{sp}\mathbf{\hat{e}}_{s+}+r_{pp}\mathbf{\hat{e}}_{p+})\mathbf{\hat{e}}_{p-}^{T}]^{*^{T}}\bigg{]}$$

$$⟨\mathbf{E}({\mathbf{r}}_1)\otimes {\mathbf{H}}^{\ast }({\mathbf{r}}_2)⟩=-k_0\Theta (\omega ,T)\int$$

$$+\frac{i}{2k_z}{e}^{i{k}_z(z_1+z_2)}[(r_{ss}{\hat{\mathbf{e}}}_{s+}+r_{ps}{\hat{\mathbf{e}}}_{p+}){\hat{\mathbf{e}}}_{p-}^T-(r_{sp}{\hat{\mathbf{e}}}_{s+}+r_{pp}{\hat{\mathbf{e}}}_{p+}){\hat{\mathbf{e}}}_{s-}^T]$$

$$+\frac{i}{2k_z^{\ast }}{e}^{-i{k}_z^{\ast }(z_1-z_2)}[{\hat{\mathbf{e}}}_{p-}{\hat{\mathbf{e}}}_{s-}^T-{\hat{\mathbf{e}}}_{s-}{\hat{\mathbf{e}}}_{p-}^T{]}^{{\ast }^T}$$

$$\displaystyle+\frac{i}{2k_{z}^{*}}e^{-ik_{z}^{*}(z_{1}+z_{2})}[(r_{ss}\mathbf{\hat{e}}_{p+}-r_{ps}\mathbf{\hat{e}}_{s+})\mathbf{\hat{e}}_{s-}^{T}+(r_{sp}\mathbf{\hat{e}}_{p+}-r_{pp}\mathbf{\hat{e}}_{s+})\mathbf{\hat{e}}_{p-}^{T}]^{*^{T}}\bigg{]}$$

$$⟨\mathbf{H}({\mathbf{r}}_1)\otimes {\mathbf{H}}^{\ast }({\mathbf{r}}_2)⟩=\frac{k_0^2\Theta (\omega ,T)}{{\mu }_0\omega }\int$$

$$+\frac{i}{2k_z}{e}^{i{k}_z(z_1+z_2)}[(r_{ss}{\hat{\mathbf{e}}}_{p+}-r_{ps}{\hat{\mathbf{e}}}_{s+}){\hat{\mathbf{e}}}_{p-}^T-(r_{sp}{\hat{\mathbf{e}}}_{p+}-r_{pp}{\hat{\mathbf{e}}}_{s+}){\hat{\mathbf{e}}}_{s-}^T]$$

$$+\frac{i}{2k_z^{\ast }}{e}^{-i{k}_z^{\ast }(z_1-z_2)}[{\hat{\mathbf{e}}}_{p-}{\hat{\mathbf{e}}}_{p-}^T+{\hat{\mathbf{e}}}_{s-}{\hat{\mathbf{e}}}_{s-}^T{]}^{{\ast }^T}$$

$$\displaystyle+\frac{i}{2k_{z}^{*}}e^{-ik_{z}^{*}(z_{1}+z_{2})}[(r_{ss}\mathbf{\hat{e}}_{p+}-r_{ps}\mathbf{\hat{e}}_{s+})\mathbf{\hat{e}}_{p-}^{T}-(r_{sp}\mathbf{\hat{e}}_{p+}-r_{pp}\mathbf{\hat{e}}_{s+})\mathbf{\hat{e}}_{s-}^{T}]^{*^{T}}\bigg{]}$$

$${ϵ}_0⟨E_x(\mathbf{r})E_y^{\ast }(\mathbf{r})⟩$$

$$\displaystyle+\frac{e^{2ik_{z}d}}{k_{z}}\bigg{(}-r_{ss}\sin\phi\cos\phi+r_{ps}\frac{k_{z}}{k_{0}}\cos^{2}\phi+r_{sp}\frac{k_{z}}{k_{0}}\sin^{2}\phi-r_{pp}\frac{k_{z}^{2}}{k_{0}^{2}}\sin\phi\cos\phi\bigg{)}$$

$$\displaystyle+\frac{e^{-2ik_{z}^{*}d}}{k_{z}^{*}}\bigg{(}-r_{ss}^{*}\sin\phi\cos\phi-r_{ps}^{*}\frac{k_{z}^{*}}{k_{0}}\sin^{2}\phi-r_{sp}^{*}\frac{k_{z}^{*}}{k_{0}}\cos^{2}\phi-r_{pp}^{*}\frac{k_{z}^{*^{2}}}{k_{0}^{2}}\sin\phi\cos\phi\bigg{)}\bigg{]}$$

$${ϵ}_0⟨E_y(\mathbf{r})E_x^{\ast }(\mathbf{r})⟩$$

$$\displaystyle+\frac{e^{2ik_{z}d}}{k_{z}}\bigg{(}-r_{ss}\sin\phi\cos\phi-r_{ps}\frac{k_{z}}{k_{0}}\sin^{2}\phi-r_{sp}\frac{k_{z}}{k_{0}}\cos^{2}\phi-r_{pp}\frac{k_{z}^{2}}{k_{0}^{2}}\sin\phi\cos\phi\bigg{)}$$

$$\displaystyle+\frac{e^{-2ik_{z}^{*}d}}{k_{z}^{*}}\bigg{(}-r_{ss}^{*}\sin\phi\cos\phi+r_{ps}^{*}\frac{k_{z}^{*}}{k_{0}}\cos^{2}\phi+r_{sp}^{*}\frac{k_{z}^{*}}{k_{0}}\sin^{2}\phi-r_{pp}^{*}\frac{k_{z}^{*^{2}}}{k_{0}^{2}}\sin\phi\cos\phi\bigg{)}\bigg{]}$$

$${\mu }_0⟨H_x(\mathbf{r})H_y^{\ast }(\mathbf{r})⟩$$

$$\displaystyle+\frac{e^{2ik_{z}d}}{k_{z}}\bigg{(}-r_{ss}\frac{k_{z}^{2}}{k_{0}^{2}}\sin\phi\cos\phi-r_{ps}\frac{k_{z}}{k_{0}}\sin^{2}\phi-r_{sp}\frac{k_{z}}{k_{0}}\cos^{2}\phi-r_{pp}\sin\phi\cos\phi\bigg{)}$$

$$\displaystyle+\frac{e^{-2ik_{z}^{*}d}}{k_{z}^{*}}\bigg{(}-r_{ss}^{*}\frac{k_{z}^{*^{2}}}{k_{0}^{2}}\sin\phi\cos\phi+r_{ps}^{*}\frac{k_{z}^{*}}{k_{0}}\cos^{2}\phi+r_{sp}^{*}\frac{k_{z}^{*}}{k_{0}}\sin^{2}\phi-r_{pp}^{*}\sin\phi\cos\phi\bigg{)}\bigg{]}$$

$${\mu }_0⟨H_y(\mathbf{r})H_x^{\ast }(\mathbf{r})⟩$$

$$\displaystyle+\frac{e^{2ik_{z}d}}{k_{z}}\bigg{(}-r_{ss}\frac{k_{z}^{2}}{k_{0}^{2}}\sin\phi\cos\phi+r_{ps}\frac{k_{z}}{k_{0}}\cos^{2}\phi+r_{sp}\frac{k_{z}}{k_{0}}\sin^{2}\phi-r_{pp}\sin\phi\cos\phi\bigg{)}$$

$$\displaystyle+\frac{e^{-2ik_{z}^{*}d}}{k_{z}^{*}}\bigg{(}-r_{ss}^{*}\frac{k_{z}^{*^{2}}}{k_{0}^{2}}\sin\phi\cos\phi-r_{ps}^{*}\frac{k_{z}^{*}}{k_{0}}\sin^{2}\phi-r_{sp}^{*}\frac{k_{z}^{*}}{k_{0}}\cos^{2}\phi-r_{pp}^{*}\sin\phi\cos\phi\bigg{)}\bigg{]}$$

$$⟨E_x(\mathbf{r})H_y^{\ast }(\mathbf{r})⟩$$

$$\displaystyle+\frac{e^{2ik_{z}d}}{k_{z}}\bigg{(}r_{ss}\frac{k_{z}}{k_{0}}\sin^{2}\phi-r_{ps}\frac{k_{z}^{2}}{k_{0}^{2}}\sin\phi\cos\phi+r_{sp}\sin\phi\cos\phi-r_{pp}\frac{k_{z}}{k_{0}}\cos^{2}\phi\bigg{)}$$

$$\displaystyle+\frac{e^{-2ik_{z}^{*}d}}{k_{z}^{*}}\bigg{(}-r_{ss}^{*}\frac{k_{z}^{*}}{k_{0}}\sin^{2}\phi+r_{ps}^{*}\sin\phi\cos\phi-r_{sp}^{*}\frac{k_{z}^{*^{2}}}{k_{0}^{2}}\sin\phi\cos\phi+r_{pp}^{*}\frac{k_{z}^{*}}{k_{0}}\cos^{2}\phi\bigg{)}\bigg{]}$$

Full text

I Derivation of Green’s function

The vector potential $\mathbf{A}(\mathbf{r}_{1})$ at $\mathbf{r}_{1}$ produced by source current density $\mathbf{J}(\mathbf{r}_{2})$ located at $\mathbf{r}_{2}$ can be calculated using Green’s function $\overline{\overline{\mathbf{G}}}(\mathbf{r}_{1},\mathbf{r}_{2})$ using the relation:

[TABLE]

Since we use Landau gauge $\mathbf{E}=i\omega\mathbf{A}$, this is same as the electric-type Green’s function that is commonly employed in the literature in the form of following equation:

[TABLE]

where $\mathbf{p}$ is polarization (dipole moment) density. Since the derivation of this Green’s function is quite well-known for isotropic media, we do not reproduce that derivation here but focus mainly on its extension to the case of bianisotropic half-space considered in the manuscript. We use Weyl’s angular spectrum representation of Green’s function. We write the position vectors as $\mathbf{r}_{j}=(\mathbf{R}_{j},z_{j})$ with transverse co-ordinates $\mathbf{R}_{j}=(x_{j},y_{j})$ for $j=[1,2]$ and wave-vectors $\mathbf{k}=(\mathbf{k}_{\parallel},k_{z})$ with transverse wavevector $\mathbf{k}_{\parallel}=k_{\parallel}(\cos\phi\mathbf{\hat{e}}_{x}+\sin\phi\mathbf{\hat{e}}_{y})$ where we introduce $\phi$ for simplicity of expressions below. In vacuum, the dispersion relation $k_{\parallel}^{2}+k_{z}^{2}=k_{0}^{2}=(\omega/c)^{2}$ follows from Maxwell’s equation, where $k_{\parallel}$ is real-valued and $k_{z}$ can be real (for $k_{\parallel}\leq k_{0}$) or complex-valued (for $k_{\parallel}>k_{0}$). In vacuum, the electric field at $\mathbf{r}_{1}$ produced by the dipole moment $\mathbf{p}=\mathbf{d}\delta(\mathbf{r}-\mathbf{r}_{2})$ located at $\mathbf{r}_{2}$ is written in the angular spectrum representation for $z_{1}\geq z_{2}$ as:

[TABLE]

The polarization vectors $\mathbf{\hat{e}}_{j\pm}$ for $j={s,p}$ with $\pm$ denoting waves going along $\pm\mathbf{\hat{e}}_{z}$ directions are:

[TABLE]

For $z_{1}<z_{2}$, the integrand is modified and contains the term $e^{-ik_{z}(z_{1}-z_{2})}[\mathbf{\hat{e}}_{s-}(\mathbf{\hat{e}}_{s-}\cdot\mathbf{d})+\mathbf{\hat{e}}_{p-}(\mathbf{\hat{e}}_{p-}\cdot\mathbf{d})]$. For consistency, we will stick with $z_{1}\geq z_{2}$ in the following discussion. The scattered/reflected field is calculated by considering the reflection of the incident field at the interface ($z=0$). Since the waves propagate in $-\mathbf{\hat{e}}_{z}$ direction to reach the interface, the incident field will be of the form $e^{ik_{z}z_{2}}[\mathbf{\hat{e}}_{s-}(\mathbf{\hat{e}}_{s-}\cdot\mathbf{d})+\mathbf{\hat{e}}_{p-}(\mathbf{\hat{e}}_{p-}\cdot\mathbf{d})]$. It undergoes reflection at the interface where polarization vectors change to $\mathbf{\hat{e}}_{s-}\rightarrow r_{ss}\mathbf{\hat{e}}_{s+}+r_{ps}\mathbf{\hat{e}}_{p+}$ and $\mathbf{\hat{e}}_{p-}\rightarrow r_{sp}\mathbf{\hat{e}}_{s+}+r_{pp}\mathbf{\hat{e}}_{p+}$. The Fresnel reflection coefficient $r_{jk}$ for $j,k=[s,p]$ describes the amplitude of $\mathbf{\hat{e}}_{j}$-polarized reflected light due to unit amplitude $\mathbf{\hat{e}}_{k}$-polarized incident light. For isotropic media, cross-polarization Fresnel coefficients $r_{sp},r_{ps}$ are zero which simplifies the calculation. But they are not necessarily zero for general bianisotropic media. The reflected field then acquires an additional phase of $e^{ik_{z}z_{1}}$ upon reaching the position $\mathbf{r}_{1}$ along with the overall transverse phase accrual same as $e^{ik_{\parallel}\cdot(\mathbf{R}_{1}-\mathbf{R}_{2})}$. This results in the scattered/reflected field at $\mathbf{r}_{1}$ given below:

[TABLE]

By writing the total field $\mathbf{E}(\mathbf{r}_{1})=\mathbf{E}_{0}(\mathbf{r}_{1})+\mathbf{E}_{\text{ref}}(\mathbf{r}_{2})$ and using Eq 2, the Green’s function $\overline{\overline{\mathbf{G}}}(\mathbf{r}_{1},\mathbf{r}_{2})$ is derived for the geometry considered in the manuscript.

II Derivation of Fluctuation Dissipation Theorem (FDT)

We follow Landau’s discussion in Ref. Pitaevskii et al. (2014) (Statistical Physics Part 2, Chapter 8) to obtain the vector potential correlations in the vacuum half-space when both the vacuum and the material half-spaces are at the same thermodynamic temperature $T$ (FDT of first kind). Let’s consider the linear response theory developed by Kubo. In this theory, we consider a discrete set of quantities denoted by $x_{a}$ for $(a=1,2,...)$ which describe the behavior of the system under certain external interactions. These interactions are described by external forces $f_{a}$ such that interaction energy has the form:

[TABLE]

The quantities $x_{a}$ are further related to to the forces $f_{a}$ through linear generalized susceptibilities $\alpha_{ab}(\mathbf{r},\mathbf{r}^{\prime})$ (linear response). In the Fourier domain, they can be written as:

[TABLE]

The spectral distribution of the fluctuating quantities $x_{a}(\omega,r)$ is related to the generalized susceptibilities by Kubo’s fluctuation dissipation relation given by:

[TABLE]

For electromagnetic fields, $x_{a}(\omega,\mathbf{r})\rightarrow A_{j}(\omega,\mathbf{r})$ ($j=x,y,z$ component of vector potential). The interaction with the externally induced current is given by $V=-\mathbf{j}\cdot\mathbf{A}$ where $\mathbf{j}(\omega,\mathbf{r})$ is the generalized force. Since vector potential and current density are related by the Green’s function:

[TABLE]

the generalized susceptibility becomes $\alpha_{ab}(\omega;\mathbf{r},\mathbf{r}^{\prime})=\mu_{0}G_{jk}(\omega;\mathbf{r},\mathbf{r}^{\prime})$. Making these substitutions in Kubo’s linear FDT given by Eq. 8, one retrieves FDT for vector potential components written in the matrix form as:

[TABLE]

$[..]^{T}$ denotes the matrix transpose and $[..]^{*}$ denotes complex conjugation. The vector quantities are written as column vectors such that $\mathbf{A}=[A_{x},A_{y},A_{z}]^{T}$. Substituting the Green’s function obtained in section I, the vector potential correlations are:

[TABLE]

Note that even though we eventually compute the correlations for $\mathbf{r}_{1}=\mathbf{r}_{2}$, we still need to expand all the terms since $\overline{\overline{\mathbf{G}}}(\mathbf{r}_{2},\mathbf{r}_{1})^{*^{T}}\neq\overline{\overline{\mathbf{G}}}(\mathbf{r}_{1},\mathbf{r}_{2})^{*^{T}}$. Furthermore, following calculations involve curl operators that act differently on different terms due to different phase factors making it necessary to calculate each term carefully.

III Derivation of Electromagnetic field correlations

The calculation of field correlations is straightforward in Landau gauge since $\mathbf{E}=i\omega\mathbf{A},\mathbf{B}=\nabla\times\mathbf{A}$:

[TABLE]

The curl operator $\nabla\times$ acts only on the exponential phase factor and not the $s,p$ polarization vectors ($\mathbf{\hat{e}}_{s\pm},\mathbf{\hat{e}}_{p\pm}$) and therefore leads to the above simplified form. The matrix form of the curl operator therefore depends on the exponential phase in each term. For instance, for a term in the correlations $\langle\mathbf{A}(\mathbf{r}_{1})\otimes\mathbf{A}^{*}(\mathbf{r}_{2})\rangle$ that has the phase factor $e^{ik_{x}(x_{1}-x_{2})+ik_{y}(y_{1}-y_{2})-ik_{z}^{*}(z_{1}-z_{2})}$, the curl operator $\nabla_{r2}\times$ is:

[TABLE]

Here $k_{x}=k_{\parallel}\cos\phi$, $k_{y}=k_{\parallel}\sin\phi$ and $k_{z}^{2}=k_{0}^{2}-k_{\parallel}^{2}$. By performing these matrix operations one can obtain the above electromagnetic field correlations. The electromagnetic field correlations are:

[TABLE]

In the following, we look at the spin angular momentum density and Poynting flux along certain $\mathbf{\hat{e}}_{z}$ direction (perpendicular to surface) and along $\mathbf{\hat{e}}_{x}$ direction (parallel to surface). We evaluate these quantities at spatial point $\mathbf{r}_{1}=\mathbf{r}_{2}=\mathbf{r}=(0,0,d)$.

Calculation of spin density and Poynting flux perpendicular to surface

[TABLE]

The heat flux density along $\mathbf{\hat{e}}_{z}$ direction is given by the Poynting flux $P_{z}=\langle\mathrm{Re}\{E_{x}^{*}H_{y}-E_{y}^{*}H_{x}\}\rangle$. From the above expressions, we get:

[TABLE]

The heat flux along $\mathbf{\hat{e}}_{z}$ direction is always zero irrespective of the material type. Similarly, we obtain electric and magnetic contributions to the spin angular momentum density along $\mathbf{\hat{e}}_{z}$ direction.

[TABLE]

It follows that the total spin angular momentum density along $\mathbf{\hat{e}}_{z}$ is always zero. Note that the individual contributions above can be nonzero in presence for $r_{sp}+r_{ps}\neq 0$ which is true for nonreciprocal materials. The above results show that heat and total angular momentum flux rates perpendicular to the surface are always zero at thermal equilibrium. Since this is a thermodynamic requirement, these results prove the consistency of fluctuational electrodynamic theory with thermodynamics.

Calculation of spin density and Poynting flux parallel to the surface

[TABLE]

Upon integration over angle $\phi$ for constant terms and cancellation of various other terms, the simplified final expressions for spin densities and flux rates are:

[TABLE]

For isotropic media, $r_{sp}=r_{ps}=0$ and $r_{ss},r_{pp}$ do not depend on the angle $\phi$. Upon integrating over $\phi$, all the Poynting flux rates and spin densities in the vicinity of isotropic media are zero. For generic biansiotropic media, one needs to compute these quantities by integrating over $\phi$. We find that for reciprocal bianisotropic media, the flux rates and spin-densities are again zero. While this requires numerical validation by integrating over angle $\phi$, it can also be proved based on time-reversal symmetry arguments as discussed in the manuscript. For nonreciprocal media (broken time-reversal symmetry), nonzero heat current and spin angular momentum are expected at thermal equilibrium. We identify them as persistent planar heat current (PPHC) and persistent thermal photon spin (PTPS) in the manuscript. The presence of PTPS and PPHC parallel to the surface does not lead to any thermodynamic contradiction.