# Tomography of near-field radiative heat exchange between
mesoscopic bodies immersed in a thermal bath

Florian Herz<EMAIL_ADDRESS>Laboratoire Charles Fabry, UMR
8501, Institut d’Optique, CNRS, Université Paris-Saclay, 2 Avenue Augustin
Fresnel, 91127 Palaiseau Cedex, France  Riccardo Messina  Laboratoire Charles
Fabry, UMR 8501, Institut d’Optique, CNRS, Université Paris-Saclay, 2 Avenue
Augustin Fresnel, 91127 Palaiseau Cedex, France  Philippe Ben-Abdallah
Laboratoire Charles Fabry, UMR 8501, Institut d’Optique, CNRS, Université
Paris-Saclay, 2 Avenue Augustin Fresnel, 91127 Palaiseau Cedex, France

(August 28, 2024)

###### Abstract

A tomographic study of near-field radiative heat exchanges between a
mesoscopic object and a substrate immersed in a thermal bath is carried out
within the theoretical framework of fluctuational electrodynamics. By using
the discrete-dipole-approximation method, we compute the power density
distribution for radiative exchanges and highlight the major role played by
many-body interactions in these transfers. Additionally, we emphasize the
close relationship between power distribution and eigenmodes within the solid
paving the way to applications for hot-spot targeting at deep sub-wavelength
scale by shape optimization.

##  I Introduction

Concentrating the heat emitted by a hot object placed at close distance to a
base material if of significant relevance for practical uses such as nano-
photolithography [1], to engrave patterning nanometer-scale structures on
material surfaces, heat-assisted magnetic recording [2, 3] for data storage in
hard-disk writing technology, or to make local temperature measurements using
the near-field scanning thermal microscopy technique [4, 5, 6]. Indeed, at
distances smaller than the thermal wavelength $\lambda_{\text{th}}=\hbar
c/k_{\rm B}T$ to the source, in which $k_{\rm B}$ is Boltzmann’s constant, the
exchanged heat flux can surpass the flux predicted by the famous Stefan-
Boltzmann’s law (blackbody limit [7]) by several orders of magnitude [8, 9,
10] due to the contribution of evanescent photons which are superimposed on
the propagative ones. In all these applications, a small section of the solid
substrate is heated due to a super-Planckian heat flux to elevate its
temperature beyond its melting point or beyond the Curie temperature of
magnetic materials to demagnetize them locally. Therefore, to understand the
underlying physics which govern the interactions of the electromagnetic field
radiated by thermal emitters with microscope tips is of crucial importance.
Despite recent progress in modeling these objects [11, 12, 13, 14, 15, 16,
17], challenges remain in fully understanding this physics. Among the open
questions is the description of many-body mechanisms [18, 19, 20] which drive
interactions and heat exchanges between mesoscopic objects and the field
radiated by surrounding thermal emitters.

In the present work, we investigate this problem in detail within the
framework of the fluctuational-electrodynamics theory [21]. We first introduce
a theoretical model to analyze the radiative heat exchanges among coupled
dipoles, a substrate, and the surrounding environment. We explore how
individual resonant modes, potential surface modes, and collective modes
contribute to these exchanges. Subsequently, employing the discrete dipole
approximation method for thermal emitters [17, 22], we compute the
distribution of the power density for radiative exchanges within a mesoscopic
object immersed in a thermal environment and positioned above a substrate at
various distances. Thus, we provide a complete tomographic image of exchanges
between each part of the object and its surrounding environment. By examining
these power exchanges and the internal power spectra across the different
regions within the solid, we reveal the diverse mechanisms driving heat
exchanges between a tip and a substrate. Beyond the practical problem of the
near-field thermal microscopy, we highlight the major role played by $N$-body
effects on the energy exchanges between a mesocopic object and its surrounding
environment. Additionally, we highlight the close connection existing between
the power distribution within a solid and its eigenmodes thereby paving the
way towards shape optimization for hot spot targeting.

The article is structured as follows: In Sec. II we introduce our discrete-
dipole-approximation method which we apply to two coupled dipoles in Sec. III
both, in vacuum (Sec. III.1) and above a substrate (Sec. III.2). Subsequently,
in Sec. IV we discuss the tomography of a 2D disk spatially (Sec. IV.1) and
spectrally (Sec. IV.2). In the spectral analysis we also investigate a ring
configuration with and without substrate (Sec. IV.2.1), special particles on a
ring (Sec. IV.2.2), and the spatial distribution for special frequencies (Sec.
IV.2.3). In Sec. V we end up with our conclusion.

##  II Theoretical framework

The system we consider is made of a collection of $N$ spherical particles of
radii $R_{\beta}$ ($\beta\in\\{1,...,N\\}$) and positions $\mathbf{r}_{\beta}$
placed in proximity to a substrate occupying the region $z<0$. Moreover, the
particles are exposed to an environmental field, which can be thought as
coming from the external boundaries, placed at large distance $z>0$ from
particles and substrate. In the context of the dipolar approximation, we are
going to describe each particle as an electric dipolar moment
$\mathbf{p}_{\beta}$. To calculate the power dissipated into each particle in
this configuration, we first need to determine the electric field and the
electric dipole moments each $\mathbf{r}_{\beta}$. The electric field
$\mathbf{E}$ at any position $\mathbf{r}$ can be decomposed into a
contribution of the environmental field $\mathbf{E}_{\text{env}}$ and the ones
which are induced by dipole moments $\mathbf{p}_{\beta}$ at the position
$\mathbf{r}_{\beta}$, yielding

$\displaystyle\mathbf{E}(\mathbf{r})$
$\displaystyle=\mathbf{E}_{\text{env}}(\mathbf{r})+\mu_{0}\omega^{2}\sum_{\beta=1}^{N}\mathds{G}_{\rm
E}(\mathbf{r},\mathbf{r}_{\beta})\mathbf{p}_{\beta}.$ (1)

Herein, $\mathds{G}_{\rm E}(\mathbf{r},\mathbf{r}_{\beta})$ describes the
electric Green’s function at observation point $\mathbf{r}$ caused by an
electric dipole moment at source point $\mathbf{r}_{\beta}$, $\omega$
represents the angular frequency and $\mu_{0}$ the vacuum permeability. Note
that since we only take into account polar materials like silicon carbide
(SiC), we can safely neglect the magnetic response of the system as it would
be necessary in the case of metals [23, 24, 25, 26, 27]. The electric dipole
moments $\mathbf{p}_{\beta}$ can be decomposed into a fluctuating part
$\mathbf{p}_{\beta,{\rm fl}}$ and an induced one $\mathbf{p}_{\beta,{\rm
ind}}$, the latter depending on the electric field in Eq. (1) at position
$\mathbf{r}_{\beta}$, so that

$\displaystyle\mathbf{p}_{\beta}$ $\displaystyle=\mathbf{p}_{\beta,{\rm
fl}}+\varepsilon_{0}\alpha_{\beta}\mathbf{E}(\mathbf{r}_{\beta})$
$\displaystyle=\mathbf{p}_{\beta,{\rm
fl}}+\varepsilon_{0}\alpha_{\beta}\mathbf{E}_{\text{env}}(\mathbf{r}_{\beta})+k_{0}^{2}\alpha_{\beta}\sum_{\gamma=1}^{N}\mathds{G}_{\rm
E}(\mathbf{r}_{\beta},\mathbf{r}_{\gamma})\mathbf{p}_{\gamma}.$ (2)

In the last expression, $\varepsilon_{0}$ denotes the vacuum permittivity,
$k_{0}=\omega/c$ the wave number in vacuum, $c$ the velocity of light in
vacuum,

$\displaystyle\alpha_{\beta}$ $\displaystyle=\frac{\alpha_{{\rm
CM},\beta}}{1-\text{i}\frac{k_{0}^{3}}{6\pi}\alpha_{{\rm CM},\beta}},$ (3)
$\displaystyle\alpha_{{\rm CM},\beta}$ $\displaystyle=4\pi
R_{\beta}^{3}\frac{\varepsilon_{\beta}-1}{\varepsilon_{\beta}+2},$ (4)

the dressed and the Clausius-Mosotti polarizabilities in the weak form of the
coupled dipole moments [19, 28, 29, 30, 31], respectively and
$\varepsilon_{\beta}$ the dielectric permittivity of particle $\beta$. In the
following, we will use SiC for both particles and substrate. Its dielectric
permittivity can be described by a Drude-Lorentz model [32]

$\displaystyle\varepsilon_{\beta}(\omega)$
$\displaystyle=\varepsilon_{\infty}\frac{\omega_{\rm LO}^{2}-\omega^{2}-{\rm
i}\omega\Gamma}{\omega_{\rm TO}^{2}-\omega^{2}-{\rm i}\omega\Gamma}$ (5)

with the following parameters: $\varepsilon_{\infty}=6.7$, $\omega_{\rm
LO}=1.827\times 10^{14}\,{\rm rad}/{\rm s}$, $\omega_{\rm TO}=1.495\times
10^{14}\,{\rm rad}/{\rm s}$, and $\Gamma=0.9\times 10^{14}\,{\rm rad}/{\rm
s}$.

We highlight that the Green’s function appearing in Eq. (2) and detailed in
Appendix A and B explicitly contains the case $\beta=\gamma$ due to the
contribution of the substrate (index s), whereas this contribution is
neglected in the vacuum contribution (index 0), so that we can write for
$\beta,\gamma=1,2,\dots,N$,

$\begin{split}\mathds{G}_{\rm
E}(\mathbf{r}_{\beta},&\mathbf{r}_{\gamma})=\mathds{G}_{{\rm
E,s}}(\mathbf{r}_{\beta},\mathbf{r}_{\beta})\delta_{\beta\gamma}\\\
&+\Bigl{[}\mathds{G}_{{\rm
E},0}(\mathbf{r}_{\beta},\mathbf{r}_{\gamma})+\mathds{G}_{\rm
E,s}(\mathbf{r}_{\beta},\mathbf{r}_{\gamma})\Bigr{]}\left(1-\delta_{\beta\gamma}\right).\end{split}$
(6)

Equation (2) can also be interpreted as one line of the block matrix equation

$\displaystyle\mathbf{p}$
$\displaystyle=\mathbf{p}_{\text{fl}}+\varepsilon_{0}\boldsymbol{\alpha}\mathbf{E}_{\text{env}}+k_{0}^{2}\boldsymbol{\alpha}\mathds{G}_{\rm
E}\mathbf{p}$ (7)

in which we define the block vectors
$\mathbf{p}=(\mathbf{p}_{1},...,\mathbf{p}_{N})^{T}$ and
$\mathbf{E}=(\mathbf{E}(\mathbf{r}_{1}),...,\mathbf{E}(\mathbf{r}_{N}))^{T}$
containing the dipole moments and electric fields, respectively, of each
particle and the block matrices $\boldsymbol{\alpha}={\rm
diag}(\alpha_{1},...,\alpha_{N})$ and $\mathds{G}_{\rm E}$ for which each
block matrix component obeys $\mathds{G}_{{\rm E},\beta\gamma}=\mathds{G}_{\rm
E}(\mathbf{r}_{\beta},\mathbf{r}_{\gamma})$. With that we can recast Eq. (7)
in the following way

$\displaystyle\mathbf{p}$
$\displaystyle=\mathds{T}^{-1}\mathbf{p}_{\text{fl}}+\varepsilon_{0}\mathds{T}^{-1}\boldsymbol{\alpha}_{\text{CM}}\mathbf{E}_{\text{env}}$
(8)

with

$\displaystyle\mathds{T}_{\beta\gamma}$
$\displaystyle=\left[\mathds{1}-k_{0}^{2}\alpha_{\beta}\mathds{G}_{\text{E,s}}(\mathbf{r}_{\beta}\mathbf{r}_{\beta})\right]\delta_{\beta\gamma}$
$\displaystyle\quad-k_{0}^{2}\alpha_{\beta}\mathds{G}_{\rm
E}(\mathbf{r}_{\beta},\mathbf{r}_{\gamma})(1-\delta_{\beta\gamma}).$ (9)

In the same block-matrix notation, we can now rewrite Eq. (1) for each dipole
by inserting Eq. (8) as

$\displaystyle\mathbf{E}$
$\displaystyle=\left[\mathds{1}+k_{0}^{2}\mathds{G}_{\rm
E}\mathds{T}^{-1}\boldsymbol{\alpha}\right]\mathbf{E}_{\text{env}}+\mu_{0}\omega^{2}\mathds{G}_{\rm
E}\mathds{T}^{-1}\mathbf{p}_{\text{fl}}.$ (10)

The power $P_{\beta}$ dissipated into a given dipole $\beta=1,2,\dots,N$ can
be written as

$\displaystyle P_{\beta}$
$\displaystyle=\frac{1}{\pi}\int_{0}^{\infty}\text{d}\omega\,\omega\,\mathcal{C}_{\rm{pE},\beta}$
(11) $\displaystyle=P_{\beta,{\rm dip}\rightarrow{\rm dip}}+P_{\beta,{\rm
back}\rightarrow{\rm dip}}+P_{\beta,{\rm sub}\rightarrow{\rm dip}},$ (12)

and is proportional to the imaginary part of the coupled dipole-field
correlation function

$\mathcal{C}_{\rm{pE},\beta}=\sum_{i=x,y,z}\text{Im}\llangle
p_{\beta,i}E_{i}^{*}(\mathbf{r}_{\beta})\rrangle.$ (13)

Note that we introduced the Fourier transform $f(t)=2{\rm
Re}\left[\int_{0}^{\infty}\frac{{\rm d}\omega}{2\pi}f(\omega)e^{-{\rm i}\omega
t}\right]$. The three contributions in Eq. (12) correspond to the heat
dissipated into each single dipole due to heat transfer solely between the
particles themselves, heat dissipated into each dipole due to the background,
and heat dissipated into each dipole in presence of the substrate,
respectively. Following the method outlined in Ref. [19], each contribution
can be recast under the form

$\displaystyle P_{\zeta\rightarrow{\rm dip},\beta}$
$\displaystyle=\int_{0}^{\infty}\frac{\text{d}\omega}{2\pi}\,\hbar\omega\,\mathcal{T}_{\beta,{\zeta\rightarrow{\rm
dip}}}$ (14)

for $\zeta\in\\{{\rm dip,back,sub}\\}$, with the transmission coefficients

$\displaystyle\mathcal{T}_{\beta,{{\rm dip}\rightarrow{\rm dip}}}$
$\displaystyle=4\sum_{\gamma=1}^{N}\sum_{i,j=x,y,z}n_{\gamma\beta}\frac{\chi_{\beta}}{|\alpha_{\beta}|^{2}}\left[\mathds{T}^{-1}\right]_{\beta\gamma,ij}\chi_{\gamma}\left[\mathds{T}^{-1\dagger}\right]_{\gamma\beta,ji},$
(15) $\displaystyle\mathcal{T}_{\beta,{{\rm back}\rightarrow{\rm dip}}}$
$\displaystyle=4\sum_{\gamma,\delta=1}^{N}\sum_{i,j,l=x,y,z}n_{\text{b}\beta}\frac{\chi_{\beta}}{|\alpha_{\beta}|^{2}}\left[\mathds{T}^{-1}\right]_{\beta\gamma,ij}\alpha_{\gamma}\,\mathcal{C}_{{\rm
b},\gamma\delta,jl}\,\alpha^{*}_{\delta}\left[\mathds{T}^{-1\dagger}\right]_{\delta\beta,li},$
(16) $\displaystyle\mathcal{T}_{\beta,{{\rm sub}\rightarrow{\rm dip}}}$
$\displaystyle=4\sum_{\gamma,\delta=1}^{N}\sum_{i,j,l=x,y,z}n_{\text{s}\beta}\frac{\chi_{\beta}}{|\alpha_{\beta}|^{2}}\left[\mathds{T}^{-1}\right]_{\beta\gamma,ij}\,\alpha_{\gamma}\,\mathcal{C}_{{\rm
s},\gamma\delta,jl}\,\alpha^{*}_{\delta}\left[\mathds{T}^{-1\dagger}\right]_{\delta\beta,li},$
(17)

and introducing the Bose-Einstein occupation probabilities

$\displaystyle n_{\beta\gamma}$
$\displaystyle=\frac{1}{e^{\frac{\hbar\omega}{k_{\rm
B}T_{\beta}}}-1}-\frac{1}{e^{\frac{\hbar\omega}{k_{\rm B}T_{\gamma}}}-1},$
(18)

with the reduced Planck’s constant $\hbar$, the Boltzmann constant $k_{\rm
B}$, the temperature $T_{\gamma}$ of particle $\gamma$, and the temperatures
$T_{\text{b}}$ and $T_{\text{s}}$ of background and substrate, respectively.
The remaining correlation functions
$\mathcal{C}_{\text{b/s},\beta\gamma,ij}=\llangle E_{i,{\rm
b/s}}(\mathbf{r}_{\beta})E_{j,{\rm b/s}}^{*}(\mathbf{r}_{\gamma})\rrangle$ are
[33, 34]

$\begin{split}&\mathcal{C}_{{\rm
b},\beta\gamma,ij}=\int\frac{\text{d}^{2}k_{\perp}}{8\pi^{2}}e^{\text{i}\mathbf{k}_{\perp}(\mathbf{x}_{\beta}-\mathbf{x}_{\gamma})}\Theta_{\rm
pr}\frac{e^{-\text{i}k_{z}(z_{\beta}-z_{\gamma})}}{k_{z}}\\\
&\quad\times\sum_{n=s,p}\Bigl{(}a^{-}_{n,i}+r_{p}e^{2\text{i}k_{z}z_{\beta}}a^{+}_{n,i}\Bigr{)}\Bigl{(}a^{-}_{n,j}+r_{p}^{*}e^{-2\text{i}k_{z}z_{\gamma}}a^{+}_{n,j}\Bigr{)},\end{split}$
(19)

and

$\begin{split}&\mathcal{C}_{{\rm
s},\beta\gamma,ij}=\int\frac{\text{d}^{2}k_{\perp}}{8\pi^{2}}e^{\text{i}\mathbf{k}_{\perp}(\mathbf{x}_{\beta}-\mathbf{x}_{\gamma})}\\\
&\quad\times\sum_{n=s,p}\biggl{[}\Theta_{\rm
pr}\frac{e^{\text{i}k_{z}(z_{\beta}-z_{\gamma})}}{k_{z}}a^{+}_{n,i}a^{+}_{n,j}(1-|r_{n}|^{2})\\\
&\quad+2\Theta_{\rm
ev}\frac{e^{-|k_{z}|(z_{\beta}+z_{\gamma})}}{|k_{z}|}a^{+}_{n,i}a^{-}_{n,j}\text{Im}(r_{n})\biggr{]},\end{split}$
(20)

where we have introduced the projectors on the propagative and evanescent
sections of the spectrum

$\displaystyle\Theta_{\rm pr}$ $\displaystyle=\Theta(k_{0}-k_{\perp}),$ (21)
$\displaystyle\Theta_{\rm ev}$ $\displaystyle=\Theta(k_{\perp}-k_{0}).$ (22)

For the correlation function of the dipole moments we have used the
fluctuation-dissipation theorem [33, 19]

$\displaystyle\llangle p_{\beta,i,{\rm fl}}p_{\gamma,j,{\rm fl}}^{*}\rrangle$
$\displaystyle=\hbar\varepsilon_{0}\left(1+2n_{\beta}\right)\chi_{\beta}\delta_{\beta\gamma}\delta_{ij}$
(23)

with

$\displaystyle\chi_{\beta}$ $\displaystyle={\rm
Im}\left(\alpha_{\beta}\right)-\frac{k_{0}^{3}}{6\pi}|\alpha_{\beta}|^{2}.$
(24)

Note that this is a very general formula to describe the heat flux dissipated
into a collection of particles immersed in a vacuum background and in the
presence of a substrate. It also allows to distinguish between the heat fluxes
dissipated into each dipole meaning that we are able to discuss the spatial
distribution of heat dissipated into each single dipole. In the following we
will consider 1D and 2D systems and restrict ourselves to the simpler
configuration of identical dipoles in the thermal-equilibrium configuration
$T_{1}=...=T_{N}=T_{\rm p}$, so that $P_{{\rm dip}\rightarrow{\rm dip}}=0$
holds. At first, we will evaluate our formula in the simplest scenario of two
particles: this will allow us to highlight some physical mechanisms that we
will also encounter for more particles. Later, we will discuss in detail the
case of a 2D disk.

##  III Two particles

In order to validate our model and to get some first physical insight, let us
start with the configuration of two particles in free space and above a
substrate, shown in Fig. 1. For the materials, we choose SiC for both the
particles and the substrate, and we fix the temperatures at $T_{\rm p}=298$ K
for both particles, $T_{\rm s}=323$ K for the substrate, and $T_{\rm b}=293$ K
for the background. The particles are at edge-to-edge distance $d=200$ nm from
the substrate and $l=3R$ lateral distance to each other for an identical
radius of $R=19$ nm, ensuring the validity of the dipole approximation.

Figure 1: Sketch of two particles with temperature $T_{\rm p}=298$ K and
radius $R=19$ nm separated by a center-to-center distance $l=3R$ made of SiC
at distance $d=200$ nm above a SiC substrate at temperature $T_{\rm s}=323$ K
immersed in a vacuum background at temperature $T_{\rm b}=293$ K.

From a very fundamental point of view, this situation can be compared to two
harmonic oscillators of identical mass $m$ coupled by a spring of stiffness
$\kappa$, both bound to a wall with identical springs of stiffness $K$. This
system is described by the system of differential equations

$\displaystyle m\ddot{x}_{1/2}$ $\displaystyle=-(K+\kappa)x_{1/2}+\kappa
x_{2/1}$ (25)

in which $x_{1/2}$ are the displacements of particle 1 and 2 and
$\ddot{x}_{1/2}$ their second time derivatives. Rewriting this in matrix
notation while making the ansatz $x_{1/2}(t)=c_{1/2}e^{{\rm i}\omega t}$
yields

$\displaystyle\begin{pmatrix}m\omega^{2}-K-\kappa&\kappa\\\
\kappa&m\omega^{2}-K-\kappa\end{pmatrix}\begin{pmatrix}c_{1}\\\
c_{2}\end{pmatrix}$ $\displaystyle=\begin{pmatrix}0\\\ 0\end{pmatrix}.$ (26)

To fulfill this equation by setting the matrix determinant to zero, one
obtains two resonance frequencies, namely $\omega_{1}=\sqrt{\frac{K}{m}}$ and
$\omega_{2}=\sqrt{\frac{2\kappa+K}{m}}$. The eigenvector
$\mathbf{x}_{1}=(1,1)^{t}$ corresponding to $\omega_{1}$ describes both
oscillators moving in the same direction and $\mathbf{x}_{1}=(-1,1)^{t}$
corresponding to $\omega_{2}$ describes both oscillators moving in opposite
directions. While it is not straighforward to translate these quantities
directly to the case of two nanoparticles described as dipoles, mainly because
in the electromagnetic scenario we also have to take into account the
different polarizations, we can conclude from this simple consideration that
we should also find these two fundamental resonances in the spectrum of the
two particles.

To show where we can find these resonances regarding our system of two
dipoles, let us go back to Eqs. (1)-(2) and recast them only considering the
induced dipole moments, neglecting in other words $\mathbf{E}_{\rm env}$ and
$\mathbf{p}_{\rm fl}$. In block matrix notation, one gets

$\displaystyle\boldsymbol{0}$ $\displaystyle=\mathds{M}\mathbf{p},$ (27)
$\displaystyle\mathds{M}$
$\displaystyle=\mathds{1}-k_{0}^{2}\alpha\mathds{G}_{\rm E}$ (28)

This equation demands ${\rm det}(\mathds{M})=0$ in the non-trivial case. The
roots of this determinant will define the resonance frequencies of the system.

###  III.1 In vacuum

We start by discussing the reference scenario of two particles in vacuum, i.e.
in the absence of a substrate. In order to correctly describe the situation of
an environmental field coming from both sides of the system ($z<0$ and $z>0$)
we can take the limit in which the substrate is placed far away from the
system of particles and it is replaced by an ideal blackbody. This amounts to
neglect the evanescent contribution appearing in Eq. (20) and set $r_{i}=0$.
Under these assumptions the two correlation functions defined in Eqs.
(19)-(20) are identical, as well as the heat flux from both sides towards the
particles ($P_{{\rm back}\rightarrow{\rm dip}}=P_{{\rm sub}\rightarrow{\rm
dip}}$). Then, Eqs. (16)-(17) give the same result for $\beta=1,2$ and become

$\displaystyle\mathcal{T}_{\beta,{{\rm back/sub}\rightarrow{\rm dip}}}$
$\displaystyle=\pi k_{0}^{2}n_{\text{bp}}\chi_{\rm
p}\left(\mathcal{T}_{\parallel}+2\mathcal{T}_{\perp}\right),$ (29)

using the transmission coefficients

$\displaystyle\mathcal{T}_{\parallel/\perp}$
$\displaystyle=\frac{\bigl{[}1+k_{0}^{4}|\alpha|^{2}|G_{{\rm
E},\parallel/\perp}|^{2}\bigr{]}\frac{k_{0}}{6\pi}+2k_{0}^{2}{\rm
Re}\bigl{[}\alpha G_{{\rm E},\parallel/\perp}\bigr{]}\mathcal{C}_{{\rm
b},12,\parallel/\perp}}{|1-k_{0}^{4}\alpha^{2}G_{{\rm
E},\parallel/\perp}^{2}|^{2}},$ (30)

with

$\displaystyle G_{{\rm E},\parallel}$ $\displaystyle=\frac{e^{{\rm
i}k_{0}l}}{2\pi l}\frac{1-{\rm i}k_{0}l}{k_{0}^{2}l^{2}},$ (31) $\displaystyle
G_{{\rm E},\perp}$ $\displaystyle=\frac{e^{{\rm i}k_{0}l}}{4\pi
l}\frac{k_{0}^{2}l^{2}+{\rm i}k_{0}l-1}{k_{0}^{2}l^{2}},$ (32)

and

$\displaystyle\mathcal{C}_{{\rm b},11/22,ij}$
$\displaystyle=\frac{k_{0}}{6\pi}\delta_{ij},$ (33)
$\displaystyle\mathcal{C}_{{\rm b},12/21,\parallel}$
$\displaystyle=\frac{\sin(k_{0}l)-k_{0}l\cos(k_{0}l)}{2\pi k_{0}^{2}l^{3}},$
(34) $\displaystyle\mathcal{C}_{{\rm b},12/21,\perp}$
$\displaystyle=\frac{k_{0}l\cos(k_{0}l)-(1-(k_{0}l)^{2}\sin(k_{0}l))}{4\pi
k_{0}^{2}l^{3}}.$ (35)

This already shows that due to the symmetry of the system the two transversal
directions ($\perp$) contribute identically to the overall result while the
longitudinal direction ($\parallel$) is different.

In Fig. 2 we show the spectrum of the power absorbed by both dipoles according
to Eq. (29). There, we also depict the logarithm of the inverse of the
determinant of the matrix defined in Eq. (28) but recast into

$\displaystyle{\rm det}\left(\mathds{M}\right)$ $\displaystyle={\rm
det}\left(\mathds{M}_{\perp}\right)^{2}{\rm
det}\left(\mathds{M}_{\parallel}\right),$ (36)

with

$\displaystyle{\rm det}\left(\mathds{M}_{\perp/\parallel}\right)$
$\displaystyle=1-k_{0}^{2}\alpha G_{{\rm E},\perp/\parallel}.$ (37)

Note that the polarizability $\alpha$ is now a scalar because we consider
identical particles. Both factors $1-k_{0}^{2}\alpha G_{{\rm
E},\parallel/\perp}$ will exhibit two roots analogously to the model of two
coupled harmonic oscillators. Obviously, the two transversal directions
provide the same resonances which is mathematically described by the square.
The resonances are shown in Fig. 2 by the peaks of the inverse of the
determinant. One can clearly see that there are four peaks, as expected, of
which only two contribute to the overall spectrum, namely the lower resonance
frequency of the longitudinal modes and the higher one of the transversal
modes. That is because only these “bright” modes contribute a dipole moment
since they describe opposite movements as explained in the previous section
[35, 36]. The “dark” modes describe translations of both dipoles in the same
direction, not providing a dipole moment to the overall system and thus not
coupling to the electromagnetic field.

Figure 2: Spectral heat transfer between the two particles and the background
(blue) shown with the determinant of matrices $\mathds{M}_{\perp/\parallel}$
(red/magenta) highlighting the resonance peaks for the transversal and
longitudinal polarizations by the black dashed lines at $\omega_{1}=1.75\times
10^{14}\,{\rm rad}/{\rm s}$ and $\omega_{2}=1.759\times 10^{14}\,{\rm
rad}/{\rm s}$.

###  III.2 Above a substrate

In the presence of a substrate, the calculations become more involved because
the reflected part of the Green’s function in Eq. (45) has to be considered as
well. Additionally, the correlation functions also contain such a contribution
taking the substrate into account. Both the reflected part of the Green’s
function and the two correlation functions are given for this configuration in
Eqs. (55)-(82).

Doing the same analysis as in the previous section with respect to the
determinant of matrix $\mathds{M}$ in Eq. (28), but this time with respect to
the $xz$ components of the reflected part of the Green’s function, we do not
have a longitudinal mode as before for two particles in vacuum but a coupling
between the longitudinal mode and the transversal one in $x$ direction giving
rise to a hybridized mode and yielding the following product of two
determinants

$\displaystyle{\rm det}\left(\mathds{M}\right)$ $\displaystyle={\rm
det}\left(\mathds{M}_{\rm trans}\right){\rm det}\left(\mathds{M}_{\rm
mix}\right),$ (38)

with

$\displaystyle{\rm det}\left(\mathds{M}_{\rm trans}\right)$
$\displaystyle=\left(1-k_{0}^{2}\alpha G_{{\rm E,s},\perp}^{\rm
id}\right)^{2}-k_{0}^{2}\alpha\left(G_{{\rm E,s},y}^{\rm dif}+G_{{\rm
E},\perp}\right)^{2},$ (39)

and

$\displaystyle{\rm det}\left(\mathds{M}_{\rm mix}\right)$
$\displaystyle=\bigl{[}\alpha^{2}k_{0}^{4}G_{\rm E,s,mix}^{2}-\bigl{(}1-\alpha
k_{0}^{2}\bigl{(}G_{\rm E,s,x}^{\rm dif}+G_{{\rm
E},\parallel}\bigr{)}\bigr{)}$ $\displaystyle\quad\times\bigl{(}1+\alpha
k_{0}^{2}\bigl{(}G_{\rm E,s,z}^{\rm dif}-G_{{\rm
E},\perp}\bigr{)}\bigr{)}\bigr{]}$
$\displaystyle\quad\times\bigl{[}\alpha^{2}k_{0}^{4}G_{\rm
E,s,mix}^{2}-\bigl{(}1-\alpha k_{0}^{2}\bigl{(}G_{\rm E,s,z}^{\rm dif}+G_{{\rm
E},\perp}\bigr{)}\bigr{)}$ $\displaystyle\quad\times\bigl{(}1+\alpha
k_{0}^{2}\bigl{(}G_{\rm E,s,x}^{\rm dif}-G_{{\rm
E},\parallel}\bigr{)}\bigr{)}\bigr{]}.$ (40)

Due to the coupling, the determinant of $\mathds{M}_{\rm mix}$ has now four
roots.

The results for the heat exchange between the two particles and the background
as well as with the substrate are shown in Fig. 3. Note that, due to our
choice of temperatures, the spectral power between particles and background is
negative. Therefore, we show $-P_{{\rm back}\rightarrow{\rm dip}}$ instead.
The spectrum for the heat flux between the particles and the background (light
blue) looks the same as in Fig. 2 apart from an overall amplification. The
spectral peaks for this contribution are at the same positions as for the case
without substrate and can also be found in the heat flux between the particles
and the substrate (dark blue). As shown by the inverse determinants, these
resonance frequencies correspond to two roots of the determinant. The fact of
not seeing all of them through a contribution to the spectral flux can be
attributed to the existence of dark and bright modes as in the case without
substrate.

In addition to these two peaks for both spectral heat fluxes, there is also a
third peak in the spectral heat flux between the dipoles and the substrate. It
corresponds to the well known resonance frequency of the surface phonon
polariton (SPhP) of SiC at $\omega_{\rm SPhP}=1.787\times 10^{14}\,{\rm
rad}/{\rm s}$ at which the parallel polarized reflection coefficient in Eq.
(49) has a maximum due to $\varepsilon_{\rm sub}(\omega_{\rm SPhP})\approx-1$.
Since this mode characterizes a strong coupling between substrate and
particles, it can only be found in the heat flux between the dipoles and the
substrate and not in the one with the background. Note that both determinants
also exhibit a root at $\omega_{\rm SPhP}$ due to the inclusion of the
reflected part of the Green’s function.

Figure 3: Spectral heat transfer between the two particles and the substrate
(dark blue) or the background (light blue) shown with the determinant of
matrices $\mathds{M}_{\perp/{\rm mix}}$ (red/magenta) highlighting the
resonance peaks for the transversal and mixed polarizations by the black and
red dashed lines at $\omega_{1}=1.75\times 10^{14}\,{\rm rad}/{\rm s}$,
$\omega_{2}=1.759\times 10^{14}\,{\rm rad}/{\rm s}$, and $\omega_{\rm
SPhP}=1.787\times 10^{14}\,{\rm rad}/{\rm s}$.

Translating these validations for two particles to the upcoming case of a 2D
disk, there will be resonances stemming from the eigenmodes of the disk itself
which coincide with the case without substrate. In the spectral heat flux
between the particles and the substrate there will be a peak at $\omega_{\rm
SPhP}$ which will dominate for the near field due to its evanescent character.

##  IV 2D disk

In the next step, we want to consider a 2D disk of radius $R=500$ nm as
depicted in Fig. 4. The parameters for background and substrate remain the
same as before. The particles are identical in temperature $T_{\rm p}=298$ K,
radius $r=19$ nm, and both particles and substrate are made of SiC. In order
to fill a disk of given radius $R$ as densely as possible with nanoparticles
we consider a close-packing arrangement and keep only the nanoparticles which
are fully inside the surface of the disk. Using this approach, a disk
completely filled with particles corresponds to $N=583$ particles. The radius
$r$ in with respect to radius $R$ is chosen such that $Nr^{2}/R^{2}=0.84$. For
comparison, this ratio is $Nr^{2}/R^{2}=0.9$ for $r=3$ nm and $N=24895$ so
that the chosen radius of $r=19$ nm already represents a disk whose surface
area is almost completely filled by dipoles and, therefore, provides a
reasonable approximation of the disk geometry in the spirit of DDA.

Figure 4: Sketch of a 2D disk with radius $R=500$ nm at temperature $T_{\rm
p}=298$ K at distance $d=200$ nm above a SiC substrate at temperature $T_{\rm
s}=323$ K immersed in a vacuum background at temperature $T_{\rm b}=293$ K.
All $N=583$ particles have a radius of $r=19$ nm and are made of SiC.

In Fig. 5 we show on the total power

$P_{\zeta\rightarrow{\rm dip}}=\sum_{\beta=1}^{N}P_{\beta,\zeta\rightarrow{\rm
dip}},$ (41)

for $\zeta\in\\{{\rm back,sub}\\}$ absorbed by the system of $N$ dipoles, and
show both contributions $-P_{{\rm back}\rightarrow{\rm dip}}$ and $P_{{\rm
sub}\rightarrow{\rm dip}}$ while varying the distance between the disk and the
substrate for different amounts of layers starting with only a ring of
particles on the outside of the disk in Fig. 4 and adding more rings towards
the inside until we end up with the full disk. One can clearly see that
$P_{{\rm sub}\rightarrow{\rm dip}}$ dominates in the near field and up to
distances $d<6\,\mu$m. This distance also coincides with the coherence length
of the SPhP mode of the SiC substrate which is the dominating coupling
mechanism in the near field between substrate and particles. On the contrary,
in the far-field regime $P_{{\rm back}\rightarrow{\rm dip}}$ becomes more
important, resulting in a sign flip for the total received flux. Let us also
stress that because of the non-monotonic behavior of $P_{{\rm
back}\rightarrow{\rm dip}}$, due to interference between the radiation
directly emitted away from substrate and particles and the one reflected at
the substrate, this non-monotony will be imprinted on the overall result in
the far field as a signature of this contribution because $P_{{\rm
sub}\rightarrow{\rm dip}}$ is decaying exponentially. $P_{{\rm
back}\rightarrow{\rm dip}}$ also increases for larger distances showing an
attenuation effect due to coupling to the substrate. Note that this increase
for $P_{{\rm back}\rightarrow{\rm dip}}$ in absolute values will
asymptotically reach the blackbody limit for even larger distances.
Additionally, we want to stress that although all curves share indeed a the
same qualitative behavior, the ratio between them does not coincide with the
ratio of the number of particles involved in each configuration.

From a numerical point of view it is also interesting to note that $P_{{\rm
sub}\rightarrow{\rm dip}}$ converges faster towards the result of the full
disk, which happens already for nine layers, than $P_{{\rm
back}\rightarrow{\rm dip}}$ for which 11 layers are required, which
corresponds to all particles but the seven at the center of the disk.
Especially in the near field in which $P_{{\rm sub}\rightarrow{\rm dip}}$
dominates, one could reduce the computation time by taking less layers into
account. Nevertheless, we will see later that this does not hold for a
qualitative discussion because the different configurations will have
strikingly different spatial power distributions.

Figure 5: $-P_{{\rm back}\rightarrow{\rm dip}}$ (a) and $P_{{\rm
sub}\rightarrow{\rm dip}}$ (b) for different numbers of layers starting with
one ring on the outside and adding more towards the inside until we obtain a
full disk.

###  IV.1 Spatial power distribution

Now, we want to turn to the spatial distribution of the heat fluxes. Again, we
will look at the evolution of a ring of particles towards the full disk to see
whether certain particles or regions are more important than others. The
parameters are as before. We consider three different distances between the
(layered) disk and the substrate to cover the near field ($d=200$ nm), the far
field ($d=10\,\mu$m), and the intermediate regime ($d=1\,\mu$m). Here, we
restrict ourselves to three figures for the setup of only one outer ring, a
crown of four layers, and the full disk [37].

Figure 6: Spatial power distribution for one layer. (a)-(c): Heat flux between
particles and background. (d)-(f): Heat flux between particles and substrate.
The columns correspond to three different distances given above the upper
panel. For better reference in chapter IV.2.2 we already highlight the
particles of concern for later spectral analysis by P1-P4. Figure 7: Spatial
power distribution for four layers. Figure 8: Spatial power distribution for
the full disk.

Let us start with the description of the lower panels (d)-(f) of Figs. 6, 7
and 8, corresponding to the heat flux between particles and substrate. In all
configurations one can clearly see for $d=200$ nm that the particles close to
substrate (bottom) experience the largest flux from the substrate. This is
clearly a signature of the SPhP mode due to the strong coupling in this near-
field regime. For larger distances $d$ this signature vanishes and we are left
with bright spots in the middle. This is reminiscent of configurational modes
of each disk layout in accordance to our findings for two particles in Sec.
III, and suggests that here we are indeed observing an indirect evidence of
eigenmodes of the ring, crown, and disk, respectively. Interestingly, we can
also find bright spots at the bottom if we stay below four layers in the
intermediate regime. The bright spot at the bottom, then, smears out towards
the middle and center. This is clearly a many-body effect due to more
particles that can couple with each other and distribute the incoming heat
flux from the substrate for more particle layers. For less layers the closest
particles still obtain more heat due to evanescent waves like frustrated modes
in this intermediate regime. Due to poorer coupling, especially in the case of
a single layer where there are only two nearest neighbors for each particle,
the bottom part still obtains much more power than the upper part.
Nevertheless, the eigenmode signature dominating in the far field can be
already seen for configurations with less than four layers. This already
represents a competition between the SPhP mode and the configurational
eigenmodes of the system.

The upper panels (a)-(c) of Figs. 6, 7 and 8, corresponding to the heat flux
between particles and background, are also worth discussing because it shows
counter-intuitive behavior. First of all, the upper part of the ring, crown,
and disk emits more heat in the near-field regime for each configuration. For
less layers this is more strictly bound to the very top and smears out again
for more layers until we find brighter spots in the middle and center but
still the upper part emits more power than the lower part. This shows that the
substrate attenuates heat transfer towards the background due to coupling
which is stronger for particles closer to the substrate. Interestingly, this
brighter upper part is even more pronounced for the intermediate regime in all
configurations. In the far field we expect again the eigenmodes to dominate
the spectrum to whom we ascribe the spatial distribution of $P_{{\rm
back}\rightarrow{\rm dip}}$ for each configuration. In the intermediate regime
we would have expected something between near and far field meaning that the
upper part is still more pronounced but also a stronger signature of the
eigenmodes. This could, for example, hint towards the influence of frustrated
modes. As another interesting feature, we find in the far field that the lower
part for configurations with less than four layers is more pronounced than the
upper region.

We would like to stress the presence of the bright spots at the center which
always occur apart from the cases of only a few layers and from $P_{{\rm
sub}\rightarrow{\rm dip}}$ in the near field. This is, as we will show later
in more detail, related to the eigenmodes of the configuration which is even
an important feature if we only remove the central particle of the full disk.
Additionally, it is important to stress the exception for configurations with
only a small number of layers. For less than four layers there are many
striking differences. To show the mechanisms behind this spatial power
distribution, we will now analyze these configurations spectrally.

###  IV.2 Spectral analysis of the dissipated power

As we did before for the two-particle case, we will now discuss the spectra
for different configurations: a ring ($N=84$ particles, one layer), a crown
($N=318$, 4 layers), and the full disk ($N=583$). The spectra are shown in
Fig. 9. In the left column, panels (a), (c), and (e) show the spectral heat
flux between particles and background $P_{{\rm back}\rightarrow{\rm dip}}$ and
panels (b), (d), and (f) in the right column the one between the particles and
the substrate $P_{{\rm sub}\rightarrow{\rm dip}}$. Everything is shown for the
three regimes: near field, far field, and intermediate regime. For $P_{{\rm
back}\rightarrow{\rm dip}}$ we find either two [Fig. 9(a), (e)] or three [Fig.
9(c)] major resonance peaks which are highlighted in Fig. 9 by the arrows. We
can also find the same resonances in $P_{{\rm sub}\rightarrow{\rm dip}}$
together with an additional peak at $\omega_{\rm SPhP}$ (red dashed line)
which surpasses the other ones in the near field. Since we can see in each
figure in 9 that this peak dominates in the near field and vanishes in the far
field, this proves the strong influence of the coupling between the dipoles of
each configuration and the SPhP mode of the substrate. As expected, the
eigenmodes of the configuration dominates in the far field for $P_{{\rm
sub}\rightarrow{\rm dip}}$ and the overall spectral amplitude goes down as
shown in Fig. 5. In the case of $P_{{\rm back}\rightarrow{\rm dip}}$, the
spectrum behaves vice versa regarding that distance dependence as explained in
the previous Section due to a decreasing attenuation due to coupling with the
substrate. Interestingly, the peak at the largest frequency is always more
pronounced for the intermediate regime compared to the other two regimes. We
also relate this to our observation that the intermediate regime has a
different spatial distribution for $P_{{\rm back}\rightarrow{\rm dip}}$. For
completeness, we also show by the black dashed line the resonance frequency of
a single spherical SiC particle at which $\varepsilon_{\rm SiC}=-2$ holds. It
is clear that this frequency does not play a significant role in the overall
spectrum.

Figure 9: Spectral heat flux between the particles and the background (left,
turned positive) and the substrate (right) for three different distances to
the substrate. (a)-(b): One layer. (c)-(d): Four layers. (e)-(f): Full disk.
Additionally, we show the frequencies corresponding to $\varepsilon_{\rm
SiC}=-1$ (red, dashed) and $\varepsilon_{\rm SiC}=-2$ (black, dashed). All
parameters as before.

####  IV.2.1 Spectrum of a ring with and without substrate

In Fig. 10 we show the same as in Fig. 9(a)-(b) for $d=200$ nm in comparison
with the case without substrate. Let us first of all stress that we can
identify all resonances with roots of the determinant in Eq. (28). However,
since we have 84 particles for this configuration, we could in principle have
232 different roots in the determinant so that it is not easy to connect each
root to one resonance but rather clusters of roots. Comparing all three
graphs, one can also clearly see that, apart from the peak highlighted by the
red dashed line, all resonances stem from the eigenmodes of the ring since
they also appear in the case without substrate. At the red dashed line the
SPhP mode has a strong influence on the spectral heat flux between dipoles and
substrates.

Interestingly, $P_{{\rm back}\rightarrow{\rm dip}}$ differs significantly from
the case without substrate at the peak at $\omega=1.777\times 10^{14}\,{\rm
rad}/{\rm s}$. We conclude, therefore, that the substrate can enhance this
resonance which we also showed in Fig. 9(a) while noticing that this
enhancement is particularly important in the intermediate regime of
$d=1\,\mu$m. To go further into details we will now focus on different
particles on the ring for further spectral analysis.

Figure 10: Comparison of the spectral heat fluxes between particles and
background (blue) and substrate (red) with the case without substrate (black).
The distance between the ring and the substrate is fixed to $d=200$ nm. All
other parameters as before.

####  IV.2.2 Special particles on a ring

For the ring configuration we compare the spectra of four particles that are
highlighted in red in the legend of Fig. 11. We compare three different
distances, one for each regime, for $P_{{\rm back/sub}\rightarrow{\rm dip}}$
[(a), (c), (e) / (b), (d), (f)]. As a reference, we also put the spectrum for
the sum over all particles of the disk in red. The four particles correspond
to the one at the top (P1), two ones that are close to the middle (P2 and P3),
and one at the very bottom close to the substrate (P4), as already depicted in
Fig. 6.

For $P_{{\rm sub}\rightarrow{\rm dip}}$, the red dashed line highlights
$\omega_{\rm SPhP}$. At this frequency we clearly observe a strong resonance
for the particle close to the substrate (P4) followed by the one in the middle
that is closer to the substrate (P3) in the near-field in Fig. 11(b). The
difference in amplitude is getting less for larger distances and seems to
vanish in the far field, which is not surprising in the context of the SPhP
mode coupling between substrate and particles. The resonance at
$\omega_{1r}=1.683\times 10^{14}\,{\rm rad}/{\rm s}$ seems to favor the
particles closer to the center. For larger distances this favoring is even
more obvious. This explains why we see these particles as bright spots in Fig.
6 in the far field.

For $P_{{\rm back}\rightarrow{\rm dip}}$ we do not observe many differences in
the far field; all graphs almost overlap. In the intermediate regime and in
the near field, however, the top particles exhibits the strongest resonance at
$\omega_{1r}$ and remains dominating up until close to $\omega_{2r}$ which is
where the middle particles become more important. Note that at this frequency
for the intermediate regime $d=1\,\mu$m, the difference between the centered
particles (P2 and P3) and the top particle (P1) increases which coincides with
the inversion of graphs for different distances in Fig. 9(a), (c), (e) at this
frequency. Finally, we will show the influence of the configurational modes
and the SPhP mode spatially.

Figure 11: Comparison of the spectral heat fluxes between certain particles
(see legend and Fig. 6) and background ((a), (c), and (e)) and substrate ((b),
(d), and (f)) for different particles (blue) and the sum of all particles on
the ring divided by 15 (red). $\omega_{\rm SPhP}$ is highlighted by the red
dashed line at the bottom. All parameters as before.

####  IV.2.3 Spatial distribution for special frequencies

Let us finally look at the spatial distribution of the heat fluxes $P_{{\rm
back/sub}\rightarrow{\rm dip}}$ for different frequencies corresponding to the
major peaks highlighted in Fig. 9 for the cases of the ring ($N=84$ particles,
one layer), a crown ($N=318$, four layers), and the full disk ($N=583$). The
calculations for all figures are performed for $d=200$ nm.

For the ring we see the largest heat flux $P_{{\rm back}\rightarrow{\rm dip}}$
at $\omega_{1r}$ at the top particle which is due to the eigenmodes of the
ring being less attenuated there. For the other two frequencies we observe the
eigenmode signature attenuated for particles closer to the substrate and also
the top particles as described in the previous Section. For $\omega_{\rm
SPhP}$ and $\omega_{2r}$ the heat flux $P_{{\rm sub}\rightarrow{\rm dip}}$ is
the most pronounced at the bottom particle due to strong coupling due to the
SPhP mode and evanescent modes in general. At $\omega_{1r}$ we re-encounter
the eigenmode signature enhanced for particles closer to the substrate.

Figure 12: Spatial distribution of the spectral heat fluxes between particles
and background ((a)-(c)) and substrate ((d)-(f)) for different frequencies
indicated above the upper panel for the ring configuration (one layer). The
distance between the ring and the substrate is fixed to $d=200$ nm. The
frequencies are the ones highlighted in Fig. 9. All other parameters as
before.

For the cases of a crown and a disk we observe a similar behavior. The SPhP
mode has a clear signature by enhancing the $P_{{\rm sub}\rightarrow{\rm
dip}}$ for the bottom particles. For the other frequencies we always find the
signature of an eigenmode that is enhanced for the bottom particles due to
evanescent waves. For $P_{{\rm back}\rightarrow{\rm dip}}$ we also clearly see
the eigenmodes of the system but attenuated at the bottom part due to coupling
with the substrate. Especially for the full disk the different spatial
distribution for each eigenmode is well visible. The only exception is the
case of $\omega_{1d}$ for the full disk where we see brighter spots closer to
the substrate .

Figure 13: As in Fig. 12 but for a crown configuration (4 layers).

Conclusively, this shows the “competition” between the eigenmodes of the
configuration and the coupling with substrate which either enhances the heat
flux especially at the SPhP mode or attenuates it for the heat flux between
particles and background.

Figure 14: As in Fig. 12 but for the full disk.

##  V Conclusion

By performing a tomographic analysis of radiative heat transfers between an
object, a substrate, and its thermal environment, we have demonstrated the
crucial role played by many-body interactions in near-field regime and
highlighted the intimate connection which exists between power distribution
and eigenmodes within a solid. Our study sheds light on the fundamental
mechanisms driving heat exchanges between mesoscopic objects and their
surrounding environment. Our study paves the way to a rational design of local
hot spots at a deep sub-wavelength scale by exploiting shape optimization of
solids. This could lead to important implications in the fields of nanoscale
thermal management, heat-assisted data recording, and nanoscale thermal
imaging.

##  VI Acknowledgments

This work was supported by the French Agence Nationale de la Recherche (ANR),
under grant ANR-21-CE30-0030 (NBODHEAT). Additionally, F. H. acknowledges
financial support by the Walter Benjamin Program of the Deutsche
Forschungsgemeinschaft (eng. German Research Foundation) under project number
519479175.

##  Appendix A General Green’s functions and correlation functions

As mentioned in section II, the Green’s function for each particle consists of
a vacuum part and a reflected part due to the presence of the substrate [38].
In Eq. (6) this is described more explicitly. Throughout this work, we use the
vacuum Green’s function

$\displaystyle\mathds{G}_{{\rm E},0}(\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle=\frac{e^{{\rm
i}k_{0}\rho}}{4\pi\rho}\biggl{[}\frac{k_{0}^{2}\rho^{2}+{\rm
i}k_{0}\rho-1}{k_{0}^{2}\rho^{2}}\mathds{1}-\frac{k_{0}^{2}\rho^{2}+3{\rm
i}k_{0}\rho-3}{k_{0}^{2}\rho^{2}}\hat{\boldsymbol{\rho}}\otimes\hat{\boldsymbol{\rho}}\biggr{]}$
(42)

with

$\displaystyle\hat{\boldsymbol{\rho}}$
$\displaystyle=\frac{\mathbf{r}-\mathbf{r}^{\prime}}{|\mathbf{r}-\mathbf{r}^{\prime}|},$
(43) $\displaystyle\rho$ $\displaystyle=|\mathbf{r}-\mathbf{r}^{\prime}|.$
(44)

For the reflected contribution, we employ the following expression

$\displaystyle\mathds{G}_{\rm E,s}(\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle=\int\frac{{\rm d}^{2}k_{\perp}}{(2\pi)^{2}}e^{{\rm
i}\mathbf{k}_{\perp}\cdot(\mathbf{x}-\mathbf{x}^{\prime})}\frac{{\rm i}e^{{\rm
i}k_{z}(z+z^{\prime})}}{2k_{z}}\sum_{n=s,p}r_{i}\mathbf{a}_{i}^{+}\otimes\mathbf{a}_{i}^{-}$
(45)

with the polarization unit vectors

$\displaystyle\mathbf{a}_{s}^{\pm}$
$\displaystyle=\frac{1}{k_{\perp}}(-k_{y},k_{x},0)^{t},$ (46)
$\displaystyle\mathbf{a}_{p}^{\pm}$
$\displaystyle=\frac{1}{k_{\perp}k_{0}}(\pm k_{x}k_{z},\pm
k_{y}k_{z},-k_{\perp}^{2})^{t}$ (47)

and the Fresnel amplitude reflection coefficients

$\displaystyle r_{s}$ $\displaystyle=\frac{k_{z}-k_{z,{\rm
sub}}}{k_{z}+k_{z,{\rm sub}}},$ (48) $\displaystyle r_{p}$
$\displaystyle=\frac{\varepsilon_{\rm sub}k_{z}-k_{z,{\rm
sub}}}{\varepsilon_{\rm sub}k_{z}+k_{z,{\rm sub}}}.$ (49)

Therein, we used the definitions

$\displaystyle\mathbf{k}_{\perp}$ $\displaystyle=(k_{x},k_{y})^{t},$ (50)
$\displaystyle\mathbf{x}$ $\displaystyle=(x,y)^{t},$ (51) $\displaystyle{\rm
d}^{2}k_{\perp}$ $\displaystyle={\rm d}k_{x}{\rm d}k_{y},$ (52) $\displaystyle
k_{z}$ $\displaystyle=\sqrt{k_{0}^{2}-k_{\perp}^{2}},$ (53) $\displaystyle
k_{z,{\rm sub}}$ $\displaystyle=\sqrt{\varepsilon_{\rm
sub}k_{0}^{2}-k_{\perp}^{2}}.$ (54)

##  Appendix B Green’s function and correlation functions for two particles
above a substrate

The reflected part of the Green’s function in Eq. (45) simplifies for two
particles as in Fig. 1 to

$\displaystyle\mathds{G}_{\rm E,s}(\mathbf{r}_{1/2},\mathbf{r}_{1/2})$
$\displaystyle=G_{{\rm E,s},\perp}^{\rm
id}\left(\hat{\mathbf{x}}\otimes\hat{\mathbf{x}}+\hat{\mathbf{y}}\otimes\hat{\mathbf{y}}\right)+G_{\rm
E,s,z}^{\rm id}\hat{\mathbf{z}}\otimes\hat{\mathbf{z}}$ (55)

with

$\displaystyle G_{{\rm E,s},\perp}^{\rm id}$ $\displaystyle={\rm
i}\int_{0}^{\infty}\frac{{\rm d}k_{\perp}k_{\perp}}{8\pi k_{z}}e^{2{\rm
i}k_{z}(d+R)}\biggl{[}r_{s}-r_{p}\frac{k_{z}^{2}}{k_{0}^{2}}\biggr{]},$ (56)
$\displaystyle G_{\rm E,s,z}^{\rm id}$ $\displaystyle={\rm
i}\int_{0}^{\infty}\frac{{\rm d}k_{\perp}k_{\perp}^{3}}{4\pi
k_{0}^{2}k_{z}}e^{2{\rm i}k_{z}(d+R)}r_{p}$ (57)

if both spacial arguments of the Green’s function are identical and to

$\displaystyle\mathds{G}_{\rm E,s}(\mathbf{r}_{1/2},\mathbf{r}_{2/1})$
$\displaystyle=G_{\rm E,s,x}^{\rm
dif}\hat{\mathbf{x}}\otimes\hat{\mathbf{x}}+G_{\rm E,s,y}^{\rm
dif}\hat{\mathbf{y}}\otimes\hat{\mathbf{y}}+G_{\rm E,s,z}^{\rm
dif}\hat{\mathbf{z}}\otimes\hat{\mathbf{z}}$ $\displaystyle\quad\pm G_{\rm
E,s,mix}\left[\hat{\mathbf{x}}\otimes\hat{\mathbf{z}}-\hat{\mathbf{z}}\otimes\hat{\mathbf{x}}\right]$
(58)

with

$\displaystyle G_{\rm E,s,x/y}^{\rm dif}$ $\displaystyle={\rm
i}\int_{0}^{\infty}\frac{{\rm d}k_{\perp}k_{\perp}}{8\pi k_{z}}e^{2{\rm
i}k_{z}(d+R)}\biggl{[}r_{s}\left(J_{0}(k_{\perp}l)\pm
J_{2}(k_{\perp}l)\right)$ $\displaystyle\quad-
r_{p}\frac{k_{z}^{2}}{k_{0}^{2}}\left(J_{0}(k_{\perp}l)\mp
J_{2}(k_{\perp}l)\right)\biggr{]}$ (59)

as well as

$\displaystyle G_{\rm E,s,z}^{\rm dif}$ $\displaystyle={\rm
i}\int_{0}^{\infty}\frac{{\rm d}k_{\perp}k_{\perp}^{3}}{4\pi
k_{0}^{2}k_{z}}e^{2{\rm i}k_{z}(d+R)}r_{p}J_{0}(k_{\perp}l),$ (60)
$\displaystyle G_{\rm E,s,mix}$ $\displaystyle=\int_{0}^{\infty}\frac{{\rm
d}k_{\perp}k_{\perp}^{2}}{4\pi k_{0}^{2}}e^{2{\rm
i}k_{z}(d+R)}r_{p}J_{1}(k_{\perp}l)$ (61)

otherwise. For the correlation function of the background fields, we obtain

$\displaystyle\mathds{C}_{\text{b},11/22}$ $\displaystyle=\mathcal{C}_{{\rm
b},\perp}^{\rm
id}\left(\hat{\mathbf{x}}\otimes\hat{\mathbf{x}}+\hat{\mathbf{y}}\otimes\hat{\mathbf{y}}\right)+\mathcal{C}_{\rm
b,z}^{\rm id}\hat{\mathbf{z}}\otimes\hat{\mathbf{z}}$ (62)

with

$\displaystyle C_{{\rm b},\perp}^{\rm id}$
$\displaystyle=\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}}{8\pi
k_{z}}\biggl{[}R_{s}^{+}(k_{\perp},d)+R_{p}^{-}(k_{\perp},d)\frac{k_{z}^{2}}{k_{0}^{2}}\biggr{]},$
(63) $\displaystyle C_{\rm b,z}^{\rm id}$
$\displaystyle=\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}^{3}}{4\pi
k_{z}k_{0}^{2}}R_{p}^{+}(k_{\perp},d)$ (64)

in the case of identical particles and

$\displaystyle\mathds{C}_{\text{b},12/21}$ $\displaystyle=\mathcal{C}_{\rm
b,x}^{\rm dif}\hat{\mathbf{x}}\otimes\hat{\mathbf{x}}+\mathcal{C}_{\rm
b,y}^{12/21}\hat{\mathbf{y}}\otimes\hat{\mathbf{y}}+\mathcal{C}_{\rm b,z}^{\rm
dif}\hat{\mathbf{z}}\otimes\hat{\mathbf{z}}$
$\displaystyle\quad+\mathcal{C}_{\rm
b,xz}^{12/21}\hat{\mathbf{x}}\otimes\hat{\mathbf{z}}+\mathcal{C}_{\rm
b,zx}^{12/21}\hat{\mathbf{z}}\otimes\hat{\mathbf{x}}$ (65)

with

$\displaystyle C_{\rm b,x/y}^{\rm dif}$
$\displaystyle=\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}}{8\pi
k_{z}}\biggl{[}R_{s}^{+}(k_{\perp},d)\left[J_{0}(k_{\perp}l)\pm
J_{2}(k_{\perp}l)\right]$
$\displaystyle\quad+R_{p}^{-}(k_{\perp},d)\frac{k_{z}^{2}}{k_{0}^{2}}\left[J_{0}(k_{\perp}l)\mp
J_{2}(k_{\perp}l)\right]\biggr{]},$ (66) $\displaystyle C_{\rm b,z}^{\rm dif}$
$\displaystyle=\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}^{3}}{4\pi
k_{z}k_{0}^{2}}R_{p}^{+}(k_{\perp},d)J_{0}(k_{\perp}l)$ (67)

as well as

$\displaystyle C_{\rm b,xz}^{12/21}$
$\displaystyle=\pm\text{i}\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}^{2}}{4\pi
k_{0}^{2}}\tilde{R}_{p}^{-}(k_{\perp},d)J_{1}(k_{\perp}l),$ (68)
$\displaystyle C_{\rm b,zx}^{12/21}$
$\displaystyle=\pm\text{i}\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}^{2}}{4\pi
k_{0}^{2}}\tilde{R}_{p}^{+}(k_{\perp},d)J_{1}(k_{\perp}l)$ (69)

if one considers different particles. Here, we introduced

$\displaystyle R_{i}^{\pm}(k_{\perp},d)$ $\displaystyle=1+|r_{i}|^{2}\pm
2\text{Re}\left(e^{2\text{i}k_{z}(d+R)}r_{i}\right),$ (70)
$\displaystyle\tilde{R}_{i}^{\pm}(k_{\perp},d)$
$\displaystyle=1-|r_{i}|^{2}\pm
2\text{i}\text{Im}\left(2e^{\text{i}k_{z}(d+R)}r_{i}\right).$ (71)

In the same way we obtain for the correlation function of the substrate fields

$\displaystyle\mathds{C}_{{\rm s},11/22}$ $\displaystyle=\mathcal{C}_{{\rm
s},\perp}^{\rm
id}\left(\hat{\mathbf{x}}\otimes\hat{\mathbf{x}}+\hat{\mathbf{y}}\otimes\hat{\mathbf{y}}\right)+\mathcal{C}_{\rm
s,z}^{\rm id}\hat{\mathbf{z}}\otimes\hat{\mathbf{z}}$ (72)

with

$\displaystyle C_{{\rm s},\perp}^{\rm id}$
$\displaystyle=\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}}{8\pi
k_{z}}\biggl{[}1-|r_{s}|^{2}+\frac{k_{z}^{2}}{k_{0}^{2}}(1-|r_{p}|^{2})\biggr{]}$
$\displaystyle\quad+\int_{k_{0}}^{\infty}\frac{\text{d}k_{\perp}k_{\perp}}{4\pi|k_{z}|}e^{-2|k_{z}|(d+R)}\biggl{[}\text{Im}(r_{s})+\frac{|k_{z}|^{2}}{k_{0}^{2}}\text{Im}(r_{p})\biggr{]}$
(74)

and

$\displaystyle C_{\rm s,z}^{\rm id}$
$\displaystyle=\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}^{3}}{4\pi
k_{z}k_{0}^{2}}(1-|r_{p}|^{2})$
$\displaystyle\quad+\int_{k_{0}}^{\infty}\frac{\text{d}k_{\perp}k_{\perp}^{3}}{2\pi|k_{z}|k_{0}^{2}}e^{-2|k_{z}|d}\text{Im}(r_{p})$
(75)

for the case of two identical particles and

$\displaystyle\mathds{C}_{\text{s},12/21}$ $\displaystyle=\mathcal{C}_{\rm
s,x}^{\rm dif}\hat{\mathbf{x}}\otimes\hat{\mathbf{x}}+\mathcal{C}_{\rm
s,y}^{\rm dif}\hat{\mathbf{y}}\otimes\hat{\mathbf{y}}+\mathcal{C}_{\rm
s,z}^{\rm dif}\hat{\mathbf{z}}\otimes\hat{\mathbf{z}}$
$\displaystyle\quad+\mathcal{C}_{\rm
s,xz}^{12/21}\hat{\mathbf{x}}\otimes\hat{\mathbf{z}}+\mathcal{C}_{\rm
s,zx}^{12/21}\hat{\mathbf{z}}\otimes\hat{\mathbf{x}}$ (76)

with

$\displaystyle C_{\rm s,x/y}^{\rm dif}$
$\displaystyle=\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}}{8\pi
k_{z}}\biggl{[}\Gamma_{+}J_{0}(k_{\perp}l)\pm\Gamma_{-}J_{2}(k_{\perp}l)\biggr{]}$
$\displaystyle\quad+\int_{k_{0}}^{\infty}\frac{\text{d}k_{\perp}k_{\perp}}{4\pi|k_{z}|}e^{-2|k_{z}|(d+R)}\biggl{[}\Delta_{+}J_{0}(k_{\perp}l)\pm\Delta_{-}J_{2}(k_{\perp}l)\biggr{]}$
(77) $\displaystyle C_{\rm s,z}^{\rm dif}$
$\displaystyle=\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}^{3}}{4\pi
k_{z}k_{0}^{2}}(1-|r_{p}|^{2})J_{0}(k_{\perp}l)$
$\displaystyle\quad+\int_{k_{0}}^{\infty}\frac{\text{d}k_{\perp}k_{\perp}^{3}}{2\pi|k_{z}|k_{0}^{2}}e^{-2|k_{z}|(d+R)}\text{Im}(r_{p})J_{0}(k_{\perp}l)$
(78)

and

$\displaystyle C_{\rm s,xz}^{12/21}$
$\displaystyle=\mp\text{i}\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}^{2}}{4\pi
k_{0}^{2}}(1-|r_{p}|^{2})J_{1}(k_{\perp}l)$
$\displaystyle\quad\pm\int_{k_{0}}^{\infty}\frac{\text{d}k_{\perp}k_{\perp}^{2}}{2\pi
k_{0}^{2}}e^{-2|k_{z}|(d+R)}\text{Im}(r_{p})J_{1}(k_{\perp}l),$ (79)
$\displaystyle C_{\rm s,zx}^{12/21}$
$\displaystyle=\mp\text{i}\int_{0}^{k_{0}}\frac{\text{d}k_{\perp}k_{\perp}^{2}}{4\pi
k_{0}^{2}}(1-|r_{p}|^{2})J_{1}(k_{\perp}l)$
$\displaystyle\quad\mp\int_{k_{0}}^{\infty}\frac{\text{d}k_{\perp}k_{\perp}^{2}}{2\pi
k_{0}^{2}}e^{-2|k_{z}|(d+R)}\text{Im}(r_{p})J_{1}(k_{\perp}l)$ (80)

in the case of two different particles. Here, we introduced

$\displaystyle\Gamma_{\pm}$
$\displaystyle=1-|r_{s}|^{2}\pm\frac{k_{z}^{2}}{k_{0}^{2}}(1-|r_{p}|^{2}),$
(81) $\displaystyle\Delta_{\pm}$
$\displaystyle=\text{Im}(r_{s})\pm\frac{|k_{z}|^{2}}{k_{0}^{2}}\text{Im}(r_{p}).$
(82)

## References

  * Srituravanich  _et al._ [2004] W. Srituravanich, N. Fang, C. Sun, Q. Luo, and X. Zhang, Plasmonic nanolithography, Nano Lett. 4, 1085 (2004).
  * Challener  _et al._ [2009] W. Challener, C. Peng, A. Itagi, D. Karns, W. Peng, Y. Peng, X. Yang, X. Zhu, N. Gokemeijer, Y.-T. Hsia, G. Ju, R. Rottmayer, M. Seigler, and E. Gage, Heat-assisted magnetic recording by a near-field transducer with efficient optical energy transfer, Nat. Photonics 3, 220 (2009).
  * Stipe  _et al._ [2010] B. Stipe, T. Strand, C. Poon, H. Balamane, T. Boone, J. Katine, J.-L. Li, V. Rawat, H. Nemoto, A. Hirotsune, O. Hellwig, R. Ruiz, E. Dobisz, D. Kercher, N. Robertson, T. Albrecht, and B. Terris, Magnetic recording at 1.5 pb m-2 using an integrated plasmonic antenna, Nat. Photonics 4, 484 (2010).
  * Wilde  _et al._ [2006] Y. D. Wilde, F. Formanek, R. Carminati, B. Gralak, P.-A. Lemoine, K. Joulain, J.-P. Mulet, Y. Chen, and J.-J. Greffet, Thermal radiation scanning tunnelling microscopy, Nature 444, 740 (2006).
  * Kittel  _et al._ [2008] A. Kittel, U. Wischnath, J. Welker, O. Huth, F.Rüting, and S.-A. Biehs, Near-field thermal imaging of nanostructured surfaces, Appl. Phys. Lett. 93, 193109 (2008).
  * Jones and Raschke [2012] A. Jones and M. Raschke, Thermal infrared near-field spectroscopy, Nano Lett. 12, 1475 (2012).
  * Planck [1914] M. Planck,  _The Theory of Heat Radiation_ (P. Blakiston’s Son & Co., Philadelphia, 1914).
  * Polder and Hove [1971] D. Polder and M. V. Hove, Theory of radiative heat transfer between closely spaced bodies, Phys. Rev. B 4, 3303 (1971).
  * Volokitin and Persson [2007] A. Volokitin and B. Persson, Near-field radiative heat transfer and noncontact friction, Rev. Mod. Phys. 79, 1291 (2007).
  * Joulain  _et al._ [2005] K. Joulain, J.-P. Mulet, F. Marquier, R. Carminati, and J.-J. Greffet, Surface electromagnetic waves thermally excited: Radiative heat transfer, coherence properties and Casimir forces revisited in the near field, Surf. Sci. Rep. 57, 59 (2005).
  * Nguyen  _et al._ [2017] K. L. Nguyen, O. Merchiers, and P.-O. Chapuis, Near-field radiative heat transfer in scanning thermal microscopy computed with the boundary element method, J. Quant. Spectrosc. Radiat. Transf. 202, 154 (2017).
  * Narayanaswamy  _et al._ [2008] A. Narayanaswamy, S. Shen, and G. Chen, Near-field radiative heat transfer between a sphere and a substrate, Phys. Rev. B 78, 115303 (2008).
  * McCauley  _et al._ [2012] A. P. McCauley, M. T. H. Reid, M. Krüger, and S. G. Johnson, Modeling near-field radiative heat transfer from sharp objects using a general three-dimensional numerical scattering technique, Phys. Rev. B 85, 165104 (2012).
  * Edalatpour and Francoeur [2016] S. Edalatpour and M. Francoeur, Near-field radiative heat transfer between arbitrarily shaped objects and a surface, Phys. Rev. B 94, 045406 (2016).
  * Herz  _et al._ [2018] F. Herz, Z. An, S. Komiyama, and S.-A. Biehs, Revisiting the dipole model for a thermal infrared near-field spectroscope, Phys. Rev. Appl. 10, 044051 (2018).
  * Herz and Biehs [2022a] F. Herz and S.-A. Biehs, Generalized coupled dipole method for thermal far-field radiation, Phys. Rev. B 105, 205422 (2022a).
  * Ben-Abdallah [2019] P. Ben-Abdallah, Multitip Near-Field Scanning Thermal Microscopy, Phys. Rev. Lett. 123, 264301 (2019).
  * Ben-Abdallah  _et al._ [2011] P. Ben-Abdallah, S.-A. Biehs, and K. Joulain, Many-body radiative heat transfer theory, Phys. Rev. Lett. 107, 114301 (2011).
  * Messina  _et al._ [2013] R. Messina, M. Tschikin, S.-A. Biehs, and P. Ben-Abdallah, Fluctuation-electrodynamic theory and dynamics of heat transfer in systems of multiple dipoles, Phys. Rev. B 88, 104307 (2013).
  * Biehs  _et al._ [2021] S.-A. Biehs, R. Messina, P. S. Venkataram, A. W. Rodriguez, J. C. Cuevas, and P. Ben-Abdallah, Near-field radiative heat transfer in many-body systems, Rev. Mod. Phys. 93, 025009 (2021).
  * Rytov  _et al._ [1989] S. Rytov, Y. Kravtsov, and V. Tatarskii,  _Principles of Statistical Radiophysics, Vol. 3: Elements of Random Fields_ (Springer, Berlin, 1989).
  * Draine [1988] B. Draine, The discrete-dipole approximation and its application to interstellar graphite grains, Astrophys. J. 333, 848 (1988).
  * Dong  _et al._ [2017] J. Dong, J. Zhao, and L. Liu, Radiative heat transfer in many-body systems: Coupled electric and magnetic dipole approach, Phys. Rev. B 95, 125411 (2017).
  * Tomchuk and Grigorchuk [2006] P. M. Tomchuk and N. I. Grigorchuk, Shape and size effects on the energy absorption by small metallic particles, Phys. Rev. B 73, 155423 (2006).
  * Dedkov and Kyasov [2007] G. Dedkov and A. A. Kyasov, Thermal radiation of nanoparticles occurring at a heated flat surface in vacuum, Tech. Phys. Lett. 33, 305 (2007).
  * Chapuis  _et al._ [2008a] P.-O. Chapuis, M. Laroche, S. Volz, and J.-J. Greffet, Radiative heat transfer between metallic nanoparticles, Appl. Phys. Lett. 92, 201906 (2008a).
  * Chapuis  _et al._ [2008b] P.-O. Chapuis, M. Laroche, S. Volz, and J.-J. Greffet, Near-field induction heating of metallic nanoparticles due to infrared magnetic dipole contribution, Phys. Rev. B 77, 125402 (2008b).
  * Ekeroth  _et al._ [2017] R. A. Ekeroth, A. García-Martín, and J. Cuevas, Thermal discrete dipole approximation for the description of thermal emission and radiative heat transfer of magneto-optical systems, Phys. Rev. B 95, 235428 (2017).
  * Lakhtakia [1992] A. Lakhtakia, Strong and weak forms of the method of moments and the coupled dipole method for scattering of time-harmonic electromagnetic fields, Int. J. Mod. Phys. C 3, 583 (1992).
  * Albaladejo  _et al._ [2010] S. Albaladejo, R. Gómez-Medina, L. Froufe-Pérez, H. Marinchio, R. Carminati, J. Torrado, G. Armelles, A. García-Martín, and J. Sáenz, Radiative corrections to the polarizability tensor of an electrically small anisotropic dielectric particle, Opt. Express 18, 3556 (2010).
  * Carminati  _et al._ [2006] R. Carminati, J.-J. Greffet, C. Henkel, and J. Vigoureux, Radiative and non-radiative decay of a single molecule close to a metallic nanoparticle, Opt. Commun. 261, 368 (2006).
  * Palik [1985] E. Palik,  _Handbook of Optical Constants of Solids_ (Academic Press, San Diego, 1985).
  * Eckhardt [1984] W. Eckhardt, Macroscopic theory of electromagnetic fluctuations and stationary radiative heat transfer, Phys. Rev. A 29, 1991 (1984).
  * Messina and Antezza [2011] R. Messina and M. Antezza, Scattering-matrix approach to Casimir-Lifshitz force and heat transfer out of thermal equilibrium between arbitrary bodies, Phys. Rev. A 84, 042102 (2011).
  * Herz and Biehs [2022b] F. Herz and S.-A. Biehs, Thermal radiation and near-field thermal imaging of a plasmonic Su-Schrieffer-Heeger chain, Appl. Phys. Lett. 121, 181701 (2022b).
  * Downing  _et al._ [2017] C. Downing, E. Mariani, and G. Weick, Retardation effects on the dispersion and propagation of plasmons in metallic nanoparticle chains, J. Phys. Condens. Matter 30, 025301 (2017).
  * [37] See Supplemental Material at [link to be inserted] for three videos containing the power distribution for the remaining 10 layouts – one for each distance.
  * Sipe [1987] J. Sipe, New Green-function formalism for surface optics, J. Opt. Soc. Am. B 4, 481 (1987).