# Radiative Corbino effect in nonreciprocal many-body
systems

Ivan Latella, Philippe Ben-Abdallah

► To cite this version:

Ivan Latella, Philippe Ben-Abdallah. Radiative Corbino effect in nonreciprocal many-body systems. Physical Review B, 2025, 111 (12), pp.L121405. 10.1103/PhysRevB.111.L121405. hal-05065255

HAL Id: hal-05065255

https://hal.science/hal-05065255v1

Submitted on 13 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Radiative Corbino effect in nonreciprocal many-body systems

Ivan Latella¹,²,* and Philippe Ben-Abdallah³,†

¹Departament de Física de la Matèria Condensada,
Universitat de Barcelona, Martí i Franquès 1, 08028 Barcelona, Spain.

²Institut de Nanociència i Nanotecnologia (IN2UB), 08028 Barcelona, Spain.

³Laboratoire Charles Fabry, UMR 8501, Institut d'Optique, CNRS,
Université Paris-Saclay, 2 Avenue Augustin Fresnel, 91127 Palaiseau Cedex, France.

When a magnetic field is applied in the perpendicular direction to a metallic disk under the action of a radial bias voltage, a tangential electric current superimposes to the radial current due to the presence of the Lorentz force which acts on electrons. Here we introduce a thermal analog of this Corbino effect in many-body systems made of nonreciprocal bodies which interact by exchanging photons in near-field regime. In systems out of thermal equilibrium with a radial temperature gradient, we demonstrate that the Poynting field in the Corbino geometry is bent in presence of an external magnetic field, giving rise to a tangential heat flux. This thermomagnetic effect could find applications in the field of thermal management and energy conversion at nanoscale.

Radiative heat exchange in many-body systems composed of multiple thermal emitters capable of cooperative interactions is an emerging field of research, driven by the unique effects arising from their collective behavior [1]. In non-reciprocal many-body systems, the time-reversal symmetry breaking in Maxwell's equations leads to a class of peculiar phenomena such as the existence of supercurrents in systems at equilibrium [2], thermal rectification [3], a giant magneto-resistance in magneto-optical networks [4] or the existence of a photon thermal Hall effect [5]. In the present letter we demonstrate that for many-body systems made of magneto-optical objects distributed in a two dimensional circular geometry, the Poynting vector is curled in presence of a radial temperature gradient under the action of an external magnetic field applied in the direction orthogonal to the disk. This twist of the Poynting field induces a heat flux transversal to the direction of the primary temperature gradient, in a similar way to the Corbino effect [6, 7] in electronics.

To illustrate this effect, we consider the system depicted in Fig. 1, which consists of N spherical particles of radius $r_p$ made of a magneto-optical material, specifically InSb, arranged on a disk in a Corbino-like geometry. Each particle is positioned at $\mathbf{r}_j$ and has temperature $T_j$ for $j = 1, \dots, N$. A radial temperature gradient is applied to the disk in such a way that particles in the inner ring have temperature $T_I$ and particles in the outer ring have temperature $T_O$. These particles can exchange electromagnetic energy between them and with the surrounding medium which can be assimilated to a thermal bath at temperature $T_b$, which we set as $T_b = T_O$. A magnetic field of magnitude B is then applied in the direction perpendicular to the disk, so the permittivity tensor of these particles at frequency $\omega$ takes the form [8, 9]

$$ \mathcal{E}(\omega) = \begin{pmatrix} \varepsilon_1 & -i\varepsilon_2 & 0 \\ i\varepsilon_2 & \varepsilon_1 & 0 \\ 0 & 0 & \varepsilon_3 \end{pmatrix} \quad (1) $$

with

$$ \frac{\varepsilon_1(\omega)}{\varepsilon_\infty} = 1 + \frac{\omega_L^2 - \omega_T^2}{\omega_T^2 - \omega^2 - i\Gamma\omega} + \frac{\omega_p^2(\omega + i\gamma)}{\omega[\omega_c^2 - (\omega + i\gamma)^2]}, \quad (2) $$

$$ \frac{\varepsilon_2(\omega)}{\varepsilon_\infty} = \frac{\omega_p^2 \omega_c}{\omega[(\omega + i\gamma)^2 - \omega_c^2]}, \quad (3) $$

$$ \frac{\varepsilon_3(\omega)}{\varepsilon_\infty} = 1 + \frac{\omega_L^2 - \omega_T^2}{\omega_T^2 - \omega^2 - i\Gamma\omega} - \frac{\omega_p^2}{\omega(\omega + i\gamma)}. \quad (4) $$

Here, $\varepsilon_\infty = 15.7$ is the infinite-frequency dielectric constant, $\omega_L = 3.62 \times 10^{13}$ rad/s is the longitudinal optical phonon frequency, $\omega_T = 3.39 \times 10^{13}$ rad/s is the transverse optical phonon frequency, $\omega_p = (\frac{ne^2}{m^*\varepsilon_0\varepsilon_\infty})^{1/2}$ is the plasma frequency of free carriers with density $n = 1.36 \times 10^{19}$ cm$^{-3}$ and effective mass $m^* = 7.29 \times 10^{-32}$ kg, $\Gamma = 5.65 \times 10^{11}$ rad/s is the phonon damping constant, $\gamma = 10^{12}$ rad/s is the free carrier damping constant and $\omega_c(B) = eB/m^*$ is the cyclotron frequency [9, 10], e being the elementary charge and $\varepsilon_0$ the permittivity of vacuum. Thus, the anisotropic polarizability of the jth particle, including the radiative corrections, can be described by the tensor [11]

$$ \omega_j(\omega) = \left(\mathbb{I} - i \frac{k_0^3}{6\pi} \omega_{0j}\right)^{-1} \omega_{0j}, \quad (5) $$

where $\mathbb{I}$ denotes the identity, $\omega_{0j}$ is the quasistatic polarizability which for a sphere made of a magneto-optical material embedded in vacuum reads

$$ \omega_{0j}(\omega) = 4\pi r_p^3 (\mathcal{E} - \mathbb{I})(\mathcal{E} + 2\mathbb{I})^{-1}, \quad (6) $$

and $k_0 = \omega/c$, c being the speed of light in vacuum.

The net power absorbed by a particle in an arbitrary network can be computed using the Landauer formalism for N-body systems [12–14]. From this formalism, the non-equilibrium power received by the jth particle reads

$$ \mathcal{P}_j = 3 \int_0^\infty \frac{d\omega}{2\pi} \sum_{k=1}^N \hbar\omega[n(\omega, T_k) - n(\omega, T_b)] \mathcal{T}_{jk}(\omega), \quad (7) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Corbino geometry. A disk-like network of magneto-optical particles is subject to a radial temperature gradient in presence of an orthogonal magnetic field B. Particles in the inner and outer rings have fixed temperatures $T_I$ and $T_O$, respectively, while those in the middle ring are allowed to reach a stationary state characterized by an average temperature $T_M$. Background thermal radiation is considered at temperature $T_b = T_O$. The particles are spheres of radius $r_p = 100$ nm and $s = 5r_p$ specifies the interparticle distance in the inner ring whose radius is $R_I$. The middle and outer rings have radius $R_M = 2R_I$ and $R_O = 3R_I$, respectively. We consider $N_I = 8$, $N_M = 16$, and $N_O = 24$ particles in the inner, middle, and outer rings, respectively.</figcaption></figure>

where $n(\omega, T) = 1/(e^{\hbar\omega/k_B T} - 1)$ is the Bose distribution at temperature $T$ and $\mathcal{T}_{jk}(\omega)$ is the energy transmission coefficient between particles $j$ and $k$, $\hbar$ and $k_B$ being the reduced Planck constant and Boltzmann constant, respectively. In the dipolar approximation the transmission coefficients read [13–15]

$$ \mathcal{T}_{jk}(\omega) = \frac{4}{3} \varepsilon_0 \mathrm{Im} \mathrm{Tr} [(\mathbb{T}^{-1})_{jk} \mathcal{K}_k (\mathbb{D}\mathbb{T}^{-1})_{jk}^\dagger], \quad (8) $$

where $(\mathbb{T}^{-1})_{jk}$ denotes the $3 \times 3$ submatrices of the $3N \times 3N$ block matrix $\mathbb{T}^{-1}$, whose inverse $\mathbb{T}$ is defined by the submatrices $\mathbb{T}_{jk} = \delta_{jk}\mathbb{I} - (1 - \delta_{jk})k_0^2 \alpha_j G_0^{jk}$ [16]. Here we introduced the particle susceptibility

$$ \mathcal{K}_j = \frac{\alpha_j - \alpha_j^\dagger}{2i} - k_0^2 \alpha_j \mathrm{Im}(G_0^{jj}) \alpha_j^\dagger, \quad (9) $$

the block matrix $\mathbb{D}$ defined by the submatrices $\mathbb{D}_{jk} = \frac{k_0^2}{\varepsilon_0} G_0^{jk}$, and the notation $G_0^{jk} \equiv G_0^{EE}(\mathbf{r}_j, \mathbf{r}_k, \omega)$ to specify the electric-electric Green tensor in vacuum evaluated at the positions of particles $j$ and $k$ [16]. Furthermore, in our Corbino configuration, particles in the inner and outer rings have fixed temperatures $T_I$ and $T_O$, respectively. The remaining particles are not thermostated and are allowed to reach stationary temperatures $T_j^{\mathrm{st}}$ for which the net power they absorb vanishes, $P_j = 0$ for $j$ in the middle ring. Due to the discrete nature of the system, the radial symmetry is only approximate and particles along the middle ring may possess different neighbor configurations. This makes that the temperatures $T_j^{\mathrm{st}}$ may be slightly different (although $P_j = 0$ for all of

them). Here we characterize the temperature $T_M$ of the middle ring by taking an average over $T_j^{\mathrm{st}}$.

To investigate how nonreciprocity and the breaking of time-reversal symmetry in Maxwell’s equations impact the energy exchange between particles, we examine the behavior of the mean Poynting vector $\langle \mathbf{S}(\mathbf{r}) \rangle$ at points $\mathbf{r}$ in the system. Since $\langle \mathbf{S}(\mathbf{r}) \rangle$ characterizes energy fluxes, the absorbed power (7) is given by $\mathcal{P}_j = -\int_{\Sigma_j} \langle \mathbf{S}(\mathbf{r}) \rangle \cdot d\Sigma_j$, where $\Sigma_j$ is a surface enclosing the $j$th particle. By definition, the statistical average of the Poynting vector can be written as $\langle \mathbf{S}(\mathbf{r}) \rangle = \int_0^\infty \frac{d\omega}{2\pi} \tilde{\mathbf{S}}(\mathbf{r}, \omega)$ with the stationary spectrum given by

$$ \tilde{\mathbf{S}}(\mathbf{r}, \omega) \equiv 2\mathrm{Re}\langle \mathbf{E}(\mathbf{r}, \omega) \times \mathbf{H}^*(\mathbf{r}, \omega) \rangle, \quad (10) $$

in which the Fourier components of the electric and magnetic fields $\mathbf{E}(\mathbf{r}, \omega)$ and $\mathbf{H}(\mathbf{r}, \omega)$, respectively, can be related to local fluctuating dipoles $\tilde{\mathbf{p}}_j$ according to

$$ \mathbf{E}(\mathbf{r}, \omega) = \mathbf{E}^b(\mathbf{r}, \omega) + \frac{k_0^2}{\varepsilon_0} \sum_j \mathbb{G}^{\mathrm{EE}}(\mathbf{r}, \mathbf{r}_j, \omega) \tilde{\mathbf{p}}_j, \quad (11) $$

$$ \mathbf{H}(\mathbf{r}, \omega) = \mathbf{H}^b(\mathbf{r}, \omega) + \frac{k_0^2}{\varepsilon_0} \sum_j \mathbb{G}^{\mathrm{HE}}(\mathbf{r}, \mathbf{r}_j, \omega) \tilde{\mathbf{p}}_j. \quad (12) $$

where $\mathbf{E}^b$ and $\mathbf{H}^b$ denote the electric and magnetic field contributions coming from the thermal bath, while the second terms in these expressions set the fields scattered by the dipoles, $\mathbb{G}^{\mathrm{EE}}$ and $\mathbb{G}^{\mathrm{HE}}$ being the full electric and magnetic Green tensors, respectively. These Green tensors read [17]

$$ \mathbb{G}^{\mathrm{EE}}(\mathbf{r}, \mathbf{r}_k, \omega) = \sum_j \mathbb{G}_0^{\mathrm{EE}}(\mathbf{r}, \mathbf{r}_j, \omega) (\mathbb{T}^{-1})_{jk}, \quad (13) $$

$$ \mathbb{G}^{\mathrm{HE}}(\mathbf{r}, \mathbf{r}_k, \omega) = \sum_j \mathbb{G}_0^{\mathrm{HE}}(\mathbf{r}, \mathbf{r}_j, \omega) (\mathbb{T}^{-1})_{jk}, \quad (14) $$

$\mathbb{G}_0^{\mathrm{HE}}(\mathbf{r}_j, \mathbf{r}_k, \omega)$ being the magnetic-electric Green tensor in vacuum. For simplicity, below we omit making explicit the dependence on frequency.

The dipole moment $\tilde{\mathbf{p}}_j = \mathbf{p}_j^f + \varepsilon_0 \alpha_j \mathbf{E}^b(\mathbf{r}_j)$ has a fluctuating contribution $\mathbf{p}_j^f$ arising from intrinsic fluctuations in the particle and an induced contribution $\varepsilon_0 \alpha_j \mathbf{E}^b(\mathbf{r}_j)$ coming from the fluctuating background electric field. Taking this into account, the averaged Poynting vector can be computed by means of the fluctuation-dissipation theorem [16, 18]

$$ \langle \mathbf{p}_j^f \otimes \mathbf{p}_k^{f*} \rangle = \frac{2\varepsilon_0}{\omega} \Theta_j \mathcal{K}_j \delta_{jk} \quad (15) $$

and [19, 20]

$$ \langle \mathbf{E}^b(\mathbf{r}_j) \otimes \mathbf{E}^{b*}(\mathbf{r}_k) \rangle = \frac{2k_0}{\varepsilon_0 c} \Theta_b \mathrm{Im}[\mathbb{G}_0^{\mathrm{EE}}(\mathbf{r}_j, \mathbf{r}_k)], \quad (16) $$

$$ \langle \mathbf{E}^b(\mathbf{r}_j) \otimes \mathbf{H}^{b*}(\mathbf{r}_k) \rangle = i \frac{2k_0}{\varepsilon_0 c} \Theta_b \mathrm{Re}[\mathbb{G}_0^{\mathrm{HE}}(\mathbf{r}_j, \mathbf{r}_k)], \quad (17) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Radiative Corbino effect characterized by the mean Poynting vector in the plane containing InSb particles. (a) In absence of the magnetic field, the energy flux exhibits a radial trend following the discrete distribution of the particles. (b) For an applied magnetic field, the Poynting vector bends due to the magneto-optical response of the particles. The radius of the particles is $r_p = 100$ nm, the interparticle distance in the inner ring is $s = 5r_p$, and the temperatures of the inner and outer rings are $T_I = 350$ K and $T_O = 300$ K, respectively.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Transmission coefficients between two InSb particles for the configuration considered in Fig. 2. (a) The coefficients coincide at zero magnetic field. (b) The coefficients are different due to nonreciprocity at nonzero applied field.</figcaption></figure>

where $\Theta_j = \hbar\omega[n(T_j) + \frac{1}{2}]$ and $\Theta_b = \hbar\omega[n(T_b) + \frac{1}{2}]$. Furthermore, from Eq. (10), the Poynting vector spectrum can be decomposed as [15] $\tilde{S}(\mathbf{r}) = \tilde{S}^{\text{tr}}(\mathbf{r}) + \tilde{S}^{0}(\mathbf{r})$ with the components of these vector fields written as

$$ \tilde{S}_{\alpha}^{\text{tr}}(\mathbf{r}) = 2\text{Re} \sum_{\beta, \gamma} \epsilon_{\alpha\beta\gamma} [\mathbf{M}^{\text{tr}}(\mathbf{r})]_{\beta\gamma}, \quad (18) $$

$$ \tilde{S}_{\alpha}^{0}(\mathbf{r}) = 2\text{Re} \sum_{\beta, \gamma} \epsilon_{\alpha\beta\gamma} [\mathbf{M}^{0}(\mathbf{r})]_{\beta\gamma}, \quad (19) $$

where $\epsilon_{\alpha\beta\gamma}$ is the Levi-Civita symbol for $\alpha, \beta, \gamma = x, y, z$. Here we have introduced the matrices

$$ \mathbf{M}^{\text{tr}}(\mathbf{r}) = \frac{2k_0^3}{\varepsilon_0 c} \sum_j (\Theta_j - \Theta_b) \mathbf{G}^{\text{EE}}(\mathbf{r}, \mathbf{r}_j) \otimes_j \mathbf{G}^{\text{HE}\dagger}(\mathbf{r}, \mathbf{r}_j) \quad (20) $$

and

$$ \mathbf{M}^{0}(\mathbf{r}) = \frac{k_0^3 \Theta_b}{i \varepsilon_0 c} \sum_j \left\{ \mathbf{G}^{\text{EE}}(\mathbf{r}, \mathbf{r}_j) \alpha_j \mathbf{G}_0^{\text{HE}}(\mathbf{r}, \mathbf{r}_j) - [\mathbf{G}^{\text{HE}}(\mathbf{r}, \mathbf{r}_j) \alpha_j \mathbf{G}_0^{\text{EE}}(\mathbf{r}, \mathbf{r}_j)]^{\dagger} \right\} \quad (21) $$

which are obtained from the electric-magnetic correlation matrix according to $\langle \mathbf{E}(\mathbf{r}) \otimes \mathbf{H}^*(\mathbf{r}) \rangle = \mathbf{M}^{\text{tr}}(\mathbf{r}) + \mathbf{M}^{0}(\mathbf{r})$. Notice that $\mathbf{M}^{\text{tr}}(\mathbf{r})$ vanishes when the system is in thermal equilibrium, but for general nonreciprocal systems $\mathbf{M}^{0}(\mathbf{r})$ does not; in particular, the latter depends on the temperature of the background and on zero-point fluctuations. However, as anticipated in Ref. [15] without proof, only $\tilde{\mathbf{S}}^{\text{tr}}(\mathbf{r})$ contributes to the heat transfer. In the supplemental material [17] we give a general demonstration on the fact that $\nabla \cdot \tilde{\mathbf{S}}^{0}(\mathbf{r}) = 0$, and therefore the absorbed power (7) is actually given by $\mathcal{P}_j = -\int_{\Sigma_j} \langle \mathbf{S}^{\text{tr}}(\mathbf{r}) \rangle \cdot d\mathbf{\Sigma}$, where $\langle \mathbf{S}^{\text{tr}}(\mathbf{r}) \rangle = \int_0^{\infty} \frac{d\omega}{2\pi} \tilde{\mathbf{S}}^{\text{tr}}(\mathbf{r}, \omega)$.

In Fig. 2, we show the in-plane component of the Poynting vector $\langle \mathbf{S}^{\text{tr}}(\mathbf{r}) \rangle$ inside the Corbino geometry in presence and without magnetic field. When $B = 0$ all particles are isotropic and the system is reciprocal. In this case, represented in Fig. 2(a), there is no symmetry breaking in the system so that there is no tangential flux through the network. On the contrary, when a magnetic field is applied orthogonally to the disk, particles become anisotropic and the symmetry of the system is broken, as shown in Fig. 2(b). The induced circulation of heat flux in the tangential direction is a generalization of the effect shown in [21] for a single particle.

The breaking of symmetry responsible for this effect is illustrated in Fig. 3, where we represent the energy transmission coefficients $\mathcal{T}_{jk}$ and $\mathcal{T}_{kj}$ between two particles $j$ and $k$ under the action of an external magnetic field. Moreover, as it can be seen in this figure, the magnitude of the transmission coefficients decreases when the magnetic field is applied, showing a reduction of the energy exchange between the particles. This reduction of
<figure><img src="image_5.png" /><figcaption>Figure 4: (a) Power absorbed by the outer ring in the considered Corbino geometry as a function of the applied magnetic field for different temperature differences $\Delta T = T_I - T_O$ with $T_O = T_b = 300$ K. (b) Corbino thermal magneto-resistance and (c) stationary temperature of the middle ring as a function of the magnetic field for the different values of $\Delta T$. (d) Variation of the temperature of the middle ring as a function of the magnetic field corresponding to (c).</figcaption></figure>

heat flux carried by the thermal photons is directly related to the occurrence of a thermal magneto-resistance $R(B)$ inside the system, as previously highlighted in [4]. Noting that the particles in the outer ring at $T_O$ exchange thermal radiation with those in the inner ring at $T_I$ (the environmental temperature is $T_b = T_O$), the thermal magneto-resistance can be computed as

$$R(B) = \frac{\Delta T}{P_O(B)}, \qquad (22)$$

where $\Delta T = T_I - T_O$ and $P_O(B) = \sum_{jO} P_j$, $j_O$ labeling particles in the outer ring. The total power absorbed by the outer ring $P_O(B)$ is shown in Fig. 4(a) for the previously considered geometrical arrangement by taking different values of $\Delta T$. It can be seen that $P_O(B)$ decreases as the field increases, so that $R(B)$ grows accordingly. Moreover, the ratio $R(B = 0)/R(B)$ as represented in Fig. 4(b) with respect to the magnitude of $B$ is almost constant for the different temperature differences we have considered, in agreement with the fact that $P_O(B)$ is a linear function of $\Delta T$ when $\Delta T$ remains relatively small. In addition, the magneto-resistance induces a drop in the stationary temperature $T_M$ of the middle ring as highlighted in Figs. 4(c) and 4(d). It is interesting to note that this temperature can undergo a variation of several degrees with changes of the magnetic field around 1 T. These results show that this Corbino effect could be used to modulate the inner temperature of a many-body system in order to locally control heat exchange through shuttling [22] and to convert heat into electric current by coupling these systems with pyroelectric devices [23].

We have demonstrated the existence of a thermal analog of the Corbino effect due to the interaction of thermal photons in magneto-optical many-body systems. We have shown, in these systems, that the Poynting vector is curled in presence of a radial temperature gradient under the action of an external magnetic field. This twist of the Poynting field induces a heat flux transversal to the

direction of the primary temperature gradient. In addition to its potential for thermal management and energy conversion, this effect could give rise to a tangential force on each component of the system, resulting in an overall torque. This opens up possibilities for developing thermal ratchets capable of converting local heat sources into mechanical work.

We thank fruitful discussions with M. Antezza and S.-A. Biehs. I.L. acknowledges financial support from the Spanish Government through Grant No. PID2021-126570NBI00 (MICIU/FEDER, UE). P.B.-A. acknowledges financial support from the French Agence Nationale de la Recherche (ANR), under grant ANR-21-CE30-0030 (NBODHEAT).

<footnote>* Electronic address: ilatella@ub.edu</footnote>

<footnote>† Electronic address: pba@institutoptique.fr</footnote>

<footnote>[1] S.-A. Biehs, R. Messina, P. S. Venkataram, A. W. Rodriguez, J. C. Cuevas, and P. Ben-Abdallah, Rev. Mod. Phys. 93, 025009 (2021).</footnote>

<footnote>[2] L. Zhu and S. Fan, Phys. Rev. Lett. 117, 134303 (2016).</footnote>

<footnote>[3] A. Ott and S.-A. Biehs, Phys. Rev. B 101, 155428 (2020).</footnote>

<footnote>[4] I. Latella and P. Ben-Abdallah, Phys. Rev. Lett. 118, 173902 (2017).</footnote>

<footnote>[5] P. Ben-Abdallah, Phys. Rev. Lett. 116, 084301 (2016).</footnote>

<footnote>[6] O. M. Corbino, Il Nuovo Cimento VI, 397-420 (1911).</footnote>

<footnote>[7] K. M. Koch, Phys. Rev. 82, 460 (1951).</footnote>

<footnote>[8] E. Moncada-Villa, V. Fernández-Hurtado, F.J. García-Vidal, A. García-Martín and J.C. Cuevas, Phys. Rev. B, 92, 125418, (2015).</footnote>

<footnote>[9] E. D. Palik, R. Kaplan, R. W. Gammon, H. Kaplan, R. F. Wallis and J. J. Quinn, Phys. Rev. B 13, 2497 (1976).</footnote>

<footnote>[10] S. Law, R. Liu, and D. Wasserman, J. Vac. Sci. Technol. B 32, 052601 (2014).</footnote>

<footnote>[11] S. Albaladejo, R. Gómez-Medina, L. S. Froufe-Pérez, H. Marincho, R. Carminati, J. F. Torrado, G. Armelles, A. García-Martín, and J. J. Sáenz, Opt. Express 18, 3556 (2010).</footnote>

<footnote>[12] P. Ben-Abdallah, S.-A. Biehs, and K. Joulain, Phys. Rev.</footnote>
Lett. **107**, 114301 (2011).

[13] R. M. Abraham Ekeroth, A. García-Martín, and J. C. Cuevas, Phys. Rev. B **95**, 235428 (2017).

[14] R. M. Abraham Ekeroth, P. Ben-Abdallah, J. C. Cuevas and A. García-Martín, ACS Photonics **5**, 705 (2018).

[15] A. Ott and S.-A. Biehs, Phys. Rev. B, **101**, 155428 (2020).

[16] R. Messina, M. Tschikin, S.-A. Biehs, and P. Ben-Abdallah, Phys. Rev. B **88**, 104307 (2013).

[17] See Supplemental Material at [link] for a derivation of

the full Green tensors and for a proof of $\nabla \cdot \tilde{S}^0(r) = 0$.

[18] H. B. Callen and T. A. Welton, Phys. Rev. **83**, 34 (1951).

[19] G. S. Agarwal, Phys. Rev. A **11**, 230 (1975).

[20] W. Eckhardt, Phys. Rev. A **29**, 1991 (1984).

[21] A. Ott, P. Ben-Abdallah, and S.-A. Biehs, Phys. Rev. B **97**, 205414 (2018).

[22] I. Latella, R. Messina, J. M. Rubi, and P. Ben-Abdallah, Phys. Rev. Lett. **121**, 023903 (2018).

[23] I. Latella and P. Ben-Abdallah, Sci Rep 11, 19489 (2021).