# WRAP: A Tool for Efficient Cross-Identification of Proper Motion Objects
Spanning Multiple Surveys

Hunter Brooks Department of Astronomy and Planetary Science, Northern Arizona
University, Flagstaff, AZ 86011, USA  IPAC, Mail Code 100-22, Caltech, 1200 E.
California Blvd., Pasadena, CA 91125, USA  J. Davy Kirkpatrick IPAC, Mail Code
100-22, Caltech, 1200 E. California Blvd., Pasadena, CA 91125, USA  Dan
Caselden Department of Astrophysics, American Museum of Natural History,
Central Park West at 79th Street, New York, NY 10024, USA  Backyard Worlds:
Planet 9  Adam C. Schneider United States Naval Observatory, Flagstaff
Station, 10391 West Naval Observatory Rd., Flagstaff, AZ 86005, USA  Aaron M.
Meisner NSF’s National Optical-Infrared Astronomy Research Laboratory, 950 N.
Cherry Ave., Tucson, AZ 85719, USA  Yadukrishna Raghu Backyard Worlds: Planet
9  Farid Cedeno Backyard Worlds: Planet 9  Jacqueline K. Faherty Department of
Astrophysics, American Museum of Natural History, Central Park West at 79th
Street, New York, NY 10024, USA  Federico Marocco IPAC, Mail Code 100-22,
Caltech, 1200 E. California Blvd., Pasadena, CA 91125, USA  Marc J. Kuchner
NASA Goddard Space Flight Center, Exoplanets and Stellar Astrophysics
Laboratory, Code 667, Greenbelt, MD 20771, USA  S.L.Casewell School of Physics
and Astronomy, University of Leicester, University Road, Leicester, LE1 7RH,
UK  The Backyard Worlds: Planet 9 Collaboration

###### Abstract

We introduce the Wide-field Retrieval of Astrodata Program (WRAP), a tool
created to aid astronomers in gathering photometric and astrometric data for
point sources that may confuse simple cross-matching algorithms because of
their faintness or motion. WRAP allows astronomers to correctly cross-identify
objects with proper motion across multiple surveys by wedding the catalog data
with its underlying images, thus providing visual confirmation of cross-
associations in real time. Developed within the Backyard Worlds: Planet 9
citizen science project, WRAP aims to aid in the characterization of faint,
high motion sources by this collaboration (and others).

Direct Imaging, Catalogs

††software: WiseView (Caselden et al. 2018)††facilities: NEOWISE, Gaia,
ESO:VISTA, UKIRT, CTIO:2MASS, FLWO:2MASS, PS2, Astro Data Lab, IRSA, CDS, MAST

##  1 Introduction

Gathering catalog photometry and astrometry is generally straight-forward,
although the process becomes notably challenging when dealing with faint or
high proper motion objects. There are multiple ways to gather this
information, like Aladin (Bonnarel et al. 2000), CDS Cross-Match Service (Boch
et al. 2012), etc. We built upon previous techniques to develop the Wide-field
Retrieval of Astrodata Program (WRAP) to improve the process. Other techniques
have several drawbacks, such as complex GUI interfaces, a need for
understanding multiple programming languages, or a lack of visual
verification. WRAP is designed to address these challenges. We start by
examining previous techniques and programs in §2. Following that, we leverage
the advantages of these prior methods to develop WRAP in §3.

##  2 Previous Techniques

Figure 1: A flowchart illustrating the WRAP procedure. Here we utilize 2MASSI
J1721039+334415 (Bardalez Gagliuffi et al. 2014). It is highlighted in red for
CatWISE, WFCAM, and PanSTARRS. The process begins by opening WRAP in a
terminal window. Then, the user inputs the coordinates, sets the search
radius, and chooses the catalogs for their search. WRAP executes API queries,
interacting with the chosen catalog APIs to retrieve relevant data. If the
designated region is found through the API, WRAP displays an image of the area
with catalog detections overlaid. Users interact by clicking on the image,
prompting WRAP to identify the nearest object and writing pertinent
information from the catalog into an output file.

Numerous additional tools have been developed to facilitate the gathering of
astrophysical data. However, discussing every software in detail here would be
too enormous a task for a research note. Instead, we have highlighted some of
the most popular techniques and software currently available to provide
context for understanding the improvements WRAP provides.

The most prevalent approach for collecting photometric and astrometric data
involves conducting cross-matches based on Right Ascension (R.A.) and
Declination (Decl.). One such program is the CDS X-Match Service (Boch et al.
2012; Pineau et al. 2020) allowing users to search various catalogs using an
input coordinate table. However, a drawback is the absence of visual
confirmation for matched data. This makes it susceptible to errors from
incorrect cross-matches if the object in question has large proper motion. In
contrast, other services offer visual confirmation, enhancing the accuracy of
gathering astrophysical data.

The Aladin catalog search feature (Bonnarel et al. 2000) shares similarities
with WRAP in terms of functionality. It enables users to load images and
superimpose them onto the corresponding source catalog. However, this method
can be hampered by slow performance, primarily when dealing with substantial
amounts of data simultaneously. Moreover, this approach is primarily designed
for single object searches, limiting its applicability for more extensive
cross-identification tasks.

Additional software includes the IRSA Viewer
(https://irsa.ipac.caltech.edu/irsaviewer), VizieR Photometry Viewer
(Ochsenbein et al. 2000), and AstroToolBox (also developed by the Backyward
Worlds: Planet 9 collaboration; Kiwy 2022). However, these tools have certain
drawbacks, briefly described here. The IRSA Viewer, for instance, lacks
catalogs beyond those supported at IPAC, limiting the available data. The
VizieR Photometry Viewer is comparable to the CDS X-Match Service in that it
lacks readily available imaging for visual confirmation. Lastly, AstroToolBox
offers both visual confirmation and cross-matching, but these features are
located in separate tabs.

##  3 WRAP Overview

To address these complications, we developed WRAP. Developed in Python 3.8
(the recommended version for using WRAP), WRAP is an open-source application
designed to have minimal dependencies and low computational expenses, ensuring
it remains freely accessible. After installation, the user inputs single or
multi-object files (in CSV, FITS, ASCII, or IPAC table formats): R.A. (deg),
Decl. (deg), search radius (arcsec), then selects the desired catalogs. The
catalogs included are: CatWISE2020 Point Source Catalog (Marocco et al. 2021),
AllWISE Source Catalog (WISE; Cutri et al. 2021), Gaia Data Release 3 (Gaia;
Gaia Collaboration et al. 2021), VISTA Science Archive catalogs (VSA; Irwin et
al. 2004), WFCAM Science Archive catalogs (WFCAM; Hambly et al. 2008), 2MASS
Point Source Catalog (2MASS; Cutri et al. 2003 and Skrutskie et al. 2006),
PanSTARRS Data Release 2 (PanSTARRS; Chambers et al. 2016), NOIRLab Source
Catalog Data Release 2 (NSC; Nidever et al. 2021), and the GALEX Catalog
(GALEX; Martin et al. 2005).

Once inputs are selected, a WISEView (Caselden et al. 2018) window will
display WISE coadd images to aid the user in finding their object. WRAP begins
by running through each catalog, showing the images with overlaid sources.
Each catalog overlay comes with customizable image features, such as image
stretch, photometric band selection, and circle size. Furthermore, users may
add notes to any selected catalog. After the user clicks on their chosen
object, WRAP identifies the nearest object through a proximity check within
the field, and the results are then saved to a .csv table (no images are saved
in WRAP).

In Figure 1, the flowchart demonstrates the intuitive functionality of WRAP,
emphasizing its user-friendly interface. The catalog overlay feature of WRAP
is essential for confirmation of high motion objects, such as the illustrated
object with its high proper motion of $\sim$1.95 arcsec yr-1 (Bardalez
Gagliuffi et al. 2014), that would be difficult to cross-identify through a
simple positional cross-match.

##  4 Conclusion

WRAP was developed to serve as a tool for astronomers seeking photometric and
astrometric data for faint and high proper motion objects. It incorporates the
most effective features from previous methods, resulting in a user-friendly
and streamlined process. This is achieved by offering visual verification,
consolidating popular catalogs in one location, and requiring no coding
experience. WRAP is accessible from GitHub at the following link:
https://github.com/huntbrooks85/WRAP. The DOI for WRAP is:
https://doi.org/10.5281/zenodo.10359982. Please refer to the ReadMe file on
the GitHub page for detailed installation and usage instructions.

##  5 Acknowledgements

This Backyard Worlds research was supported by NASA grant 2017-ADAP17-0067.
All PanSTARRS data used in this paper can be found in MAST:
http://dx.doi.org/10.17909/s0zg-jx37 (catalog 10.17909/s0zg-jx37)

## References

  * Bardalez Gagliuffi et al. (2014) Bardalez Gagliuffi, D. C., Burgasser, A. J., Gelino, C. R., et al. 2014, ApJ, 794, 143. doi:10.1088/0004-637X/794/2/143
  * Boch et al. (2012) Boch, T., Pineau, F., & Derriere, S. 2012, Astronomical Data Analysis Software and Systems XXI, 461, 291
  * Bonnarel et al. (2000) Bonnarel, F., Fernique, P., Bienaymé, O., et al. 2000, A&AS, 143, 33. doi:10.1051/aas:2000331
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1. doi:10.1051/0004-6361/202039657
  * Caselden et al. (2018) Caselden, D., Westin, P., Meisner, A., et al. 2018, Astrophysics Source Code Library. ascl:1806.004
  * Chambers et al. (2016) Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv:1612.05560. doi:10.48550/arXiv.1612.05560
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, “The IRSA 2MASS All-Sky Point Source Catalog, NASA/IPAC Infrared Science Archive.
  * Cutri et al. (2021) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2021, VizieR Online Data Catalog, II/328
  * Dye et al. (2018) Dye, S., Lawrence, A., Read, M. A., et al. 2018, MNRAS, 473, 5113. doi:10.1093/mnras/stx2622
  * Emerson et al. (2004a) Emerson, J. P., Sutherland, W. J., McPherson, A. M., et al. 2004, The Messenger, 117, 27
  * Emerson et al. (2004b) Emerson, J. P., Irwin, M. J., Lewis, J., et al. 2004, Proc. SPIE, 5493, 401. doi:10.1117/12.551582
  * Hambly et al. (2008) Hambly, N. C., Collins, R. S., Cross, N. J. G., et al. 2008, MNRAS, 384, 637. doi:10.1111/j.1365-2966.2007.12700.x
  * Irwin et al. (2004) Irwin, M. J., Lewis, J., Hodgkin, S., et al. 2004, Proc. SPIE, 5493, 411. doi:10.1117/12.551449
  * Kiwy (2022) Kiwy, F. 2022, Astrophysics Source Code Library. ascl:2201.002
  * Lawrence et al. (2007) Lawrence, A., Warren, S. J., Almaini, O., et al. 2007, MNRAS, 379, 1599. doi:10.1111/j.1365-2966.2007.12040.x
  * Marocco et al. (2021) Marocco, F., Eisenhardt, P. R. M., Fowler, J. W., et al. 2021, ApJS, 253, 8. doi:10.3847/1538-4365/abd805
  * Martin et al. (2005) Martin, D. C., Fanson, J., Schiminovich, D., et al. 2005, ApJ, 619, L1. doi:10.1086/426387
  * Nidever et al. (2021) Nidever, D. L., Dey, A., Fasbender, K., et al. 2021, AJ, 161, 192. doi:10.3847/1538-3881/abd6e1
  * Ochsenbein et al. (2000) Ochsenbein, F., Bauer, P., & Marcout, J. 2000, A&AS, 143, 23. doi:10.1051/aas:2000169
  * Pineau et al. (2020) Pineau, F.-X., Boch, T., Derrière, S., et al. 2020, Astronomical Data Analysis Software and Systems XXVII, 522, 125
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163. doi:10.1086/498708


# Heat transfer in the extreme near-field regime between metallic bodies

Mauricio Gómez Viloria<EMAIL_ADDRESS>Université
Paris-Saclay, Institut d’Optique Graduate School, CNRS, Laboratoire Charles
Fabry, 91127, Palaiseau, France  Yangyu Guo  Institut Lumière Matière,
Université Claude Bernard Lyon 1, CNRS, Université de Lyon, 69622
Villeurbanne, France  Samy Merabia  Institut Lumière Matière, Université
Claude Bernard Lyon 1, CNRS, Université de Lyon, 69622 Villeurbanne, France
Philippe Ben-Abdallah  Université Paris-Saclay, Institut d’Optique Graduate
School, CNRS, Laboratoire Charles Fabry, 91127, Palaiseau, France  Riccardo
Messina<EMAIL_ADDRESS>Université Paris-Saclay, Institut
d’Optique Graduate School, CNRS, Laboratoire Charles Fabry, 91127, Palaiseau,
France

###### Abstract

We analyze the heat transfer between two metals separated by a vacuum gap in
the extreme near-field regime. In this cross-over regime between conduction
and radiation heat exchanges are mediated by photon, phonon and electron
tunneling. We quantify the relative contribution of these carriers with
respect to both the separation distance between the two bodies and the applied
bias voltage. In presence of a weak bias ($V_{\mathrm{b}}<100$ mV), electrons
and phonons can contribute equally to the heat transfer near contact, while
the contribution of photons becomes negligible. On the other hand, for larger
bias voltages, electrons play a dominant role. Moreover, we demonstrate that
depending on the sign of this bias, electrons can either heat up the cold body
or pump heat from it by Nottingham effect. Our results emphasize some
inconsistencies in recent experimental results about heat exchanges in extreme
near-field regime and it sets a roadmap for future experiments.

##  I Introduction

Two solids at different temperatures which are separated by a vacuum gap
exchange heat by radiation. This exchange of thermal photons is limited, in
the far-field regime, by the famous Stefan-Boltzmann’s law which set the power
two perfect absorbers (i.e. blackbodies) at the same temperature as the
interacting solids can exchange. In near-field regime, (i.e. at separation
distances smaller than the thermal wavelengths of solids) the situation
radically changes and the power exchanged between these solids can surpass
this blackbody limit thanks to tunneling of evanescent (i.e. non-propagative)
photons [1]. In particular when the materials involved support surface
resonant modes, such as surface-polaritons or a continuum of hyperbolic modes,
the exchanged power can overcome the blackbody limit by several orders of
magnitude [2, 3, 4, 5]. This enhancement of heat exchange predicted by Polder
and van Hove in their seminal work, establishing the foundations of
fluctuational electrodynamics (FED), has been verified by numerous experiments
[6, 7, 8, 9, 10]. This result and the possibility to tune the radiative heat
flux at sub-wavelength scale have opened new possibilities for the development
of innovative technologies for nanoscale thermal management [11], solid-state
cooling [12, 13], heating-assisted data storage [14, 15], IR sensing and
spectroscopy [16, 17] and has paved the way to a new generation of energy-
conversion devices [18, 19, 20, 21, 22]. At closer separation distances, when
the objects are separated by atomic distances, further dramatic changes occur.
Indeed, in this cross-over regime between conduction and radiation also called
extreme near-field regime, some effects and new channels for heat transfer,
which are not taken into account by the theory of Polder and van Hove, appear.
More specifically, at the atomic scale the nonlocal optical response of
materials must be taken into account to properly describe the radiative
exchanges [40, 42]. Moreover, acoustic phonons and electrons participate to
the transfer through tunneling mechanisms.

Beyond its fundamental interest, the understanding of extreme near-field
regime is of prime importance for the ongoing miniaturization of thermal
management and energy technology. However, to date, this physics remains
largely elusive and very few experimental works have been reported [23, 24].
Moreover, these works lead to contradictory conclusions. On the one hand the
experiment performed in Kittel’s group [23] shows a strong deviation with
respect to Polder’s predictions and seems to demonstrate an extraordinary
large heat flux which is four orders of magnitude larger than the values
predicted by conventional theory of fluctuational electrodynamics. On the
other hand, the measurements carried out by Reddy’s group seems to perfectly
reproduce Polder’s predictions even down to the atomic scale. Due to the lack
of a general theory to describe the relative role of photons, phonons and
electrons in energy exchanges at this scale, this problem remains today in
debate. In this paper we introduce a general framework to describe all
channels of heat transport in extreme near-field regime between two metals
with emphasis on the transfer by electron tunneling which plays a major role
at atomic scale. We successively study heat carried by elastic vibrations
(tunneling of acoustic phonons), photons (near-field radiative heat transfer)
and finally free charges (electron tunneling) in the presence of an external
bias voltage. This work allows us, on the one hand, to emphasize some limits
in the experimental works carried out in Kittel and Reddy’s groups. On the
other hand, it allows to quantify the relative contribution of different
energy carriers with respect to the separation distance between two solids and
the applied bias voltage, setting so a potential road map for future
experiments.

During the last four years some attempts to model heat exchanges at atomic
scale have been performed [25, 26, 27, 28]. However, in these works the
Nottingham effect [29], that is the heating or cooling stemming from tunneling
electrons, has been totally ignored. Here we include this effect in the
definition of heat flux carried by electrons and demonstrate that it
significantly modify both quantitatively and qualitatively exchanges in the
crossover regime between conduction and radiation. Moreover, unlike the
previous works, the tunneling probability of electrons is calculated from a
rigrous approach based on the transfer-matrix method applied on a Thomas-Fermi
description of the electronic potential barrier. This approach allows us to
explore the transfer mediated by electrons with an arbitrary bias voltage
applied between the two solids in interaction.

Our paper is structured as follows. In Sec. II we address the geometrical
configuration of two parallel planes and discuss the formalism used to
describe the contribution to heat exchange associated with the three carriers,
namely phonons, photons and electrons. We conclude this section by showing
their relative contribution as a function of the separation distance and
applied potential bias. In Sec. III, we employ the Derjaguin approximation to
discuss the heat transfer in the extreme near field in the tip–plane
configuration, and compare our theoretical prediction to the recent
experimental results. Finally, in Sec. IV we draw some conclusions and discuss
some possible perspectives.

##  II Contribution of the different carriers in plane–plane configuration

In this paper, we consider two semi-infinite metallic media with parallel
planar surfaces at temperatures $T_{1}$ and $T_{2}$, respectively, separated
by a vacuum gap of thickness $d$ in the extreme near-field regime, namely
below 10 nm and down to the angstrom range as shown in Fig. 1.

Figure 1: Sketch of the system of interest. Two semi-infinite metallic slabs
separated by a vacuum gap of thickness $d$ are set at temperatures $T_{1}$ and
$T_{2}$ and have a chemical potential $\mu_{1}$ and $\mu_{2}$. A bias voltage
$V_{\mathrm{b}}$ is applied between the two bodies. Close to contact (i.e.
extreme near-field regime), electrons (el), phonons (ph) and photons (rad)
tunnel through the separation gap. The different arrows associated with the
electronic flux reflect the possible non-reciprocity of the latter (see text
for details).

Throughout this paper we identify medium 1 as the hotter body and medium 2 as
the colder body, i.e. $T_{1}\geq T_{2}$. We assume that the total heat flux
leaving media $i=1,2$ can be written as a sum of the contribution to the heat
flux associated with acoustic phonons (ph), electromagnetic radiation (rad)
and electronic tunnelling currents (el), as

$\Phi_{1}=\Phi^{\mathrm{(}ph)}+\Phi^{\mathrm{(}rad)}+\Phi^{\mathrm{(}el)}_{1}\;,$
(1a)
$\Phi_{2}=-\left(\Phi^{\mathrm{(}ph)}+\Phi^{\mathrm{(}rad)}+\Phi^{\mathrm{(}el)}_{2}\right)\;,$
(1b)

where the $\Phi^{\mathrm{(}ph)}$, $\Phi^{\mathrm{(}rad)}$,
$\Phi^{\mathrm{(el)}}_{i}$ are defined in the following subsections in Eqs.
(2), (7) and (18), respectively, under the form of a Landauer-like expression.
The assumption that the total heat flux can be separated as the sum of the
different contributions is expected to be valid when there is no strong
coupling mechanisms between carriers (small bias voltages and temperature
difference). Note that our sign convention in Eqs. (1a) and (1b) are such that
a positive (negative) heat flux on a given body is an outgoing (ingoing) flux,
i.e. it tends to decrease (increase) its temperature. We also account for the
presence of a bias voltage $V_{\mathrm{b}}$ between the bodies, affecting the
contribution of each of the three heat carriers (ph, rad, el), where
$V_{\mathrm{b}}>0$ is defined by taking body 2 as the positive electrode as
shown in Fig. 1. Note also that the electronic contribution to the total heat
flux is not necessarily reciprocal, i.e $\Phi_{1}\neq-\Phi_{2}$, a lack of
symmetry stemming from the difference in statistics between bosons (ph/rad)
and electrons in the presence of a bias voltage, as discussed in Sec. II.3. In
the following we describe in detail the theoretical framework used to
calculate the contribution of each carrier. As a reference, the numerical
values for the physical constants used in our calculations are provided in
App. A.

###  II.1 Acoustic phonons

In crystalline media, phonons are one of the key players of heat conduction in
the bulk. However, near the edges of the sample these quasiparticles only
modify its surface since lattice vibrations are only defined inside the
material. The existence of such perturbations of material–vacuum interfaces,
along with the presence of forces between the two surfaces, leads to the
possibility of a phonon to tunnel across a vacuum gap between the two bodies,
paving the way to an additional energy-transfer channel. This possibility was
first considered between metallic media in Refs. [34, 35] within a continuum
elastic-medium approach, where van der Waals forces play the role of
interaction mechanism between the two bodies. The heat flux between the metals
was reduced to an integral over the phonon energy of the transmission
probability, derived by solving the elastic waves equations. It was later
found that an equivalent approach could be obtained by using the fluctuation-
dissipation theorem on the surface displacements [36, 37], in the same spirit
of FED. Within the latter approach, the phonon heat flux can be written as

$\displaystyle\Phi^{(\mathrm{p}h)}(T_{1},T_{2},d)=$ (2) $\displaystyle
2\int_{0}^{\infty}\frac{\mathrm{d}\omega}{2\pi}\hbar\omega\Delta
n_{\mathrm{B}E}(\omega,T_{1},T_{2})\int_{0}^{\omega/c_{\mathrm{m}in}}\frac{\mathrm{d}k}{2\pi}\;k\;\mathcal{T}^{\mathrm{(}ph)}(k,\omega,d),$

where

$\Delta
n_{\mathrm{B}E}(\omega,T_{1},T_{2})=n_{\mathrm{B}E}(\omega,T_{1})-n_{\mathrm{B}E}(\omega,T_{2}),$
(3)

$n_{\mathrm{B}E}(\omega,T_{i})=1/[\exp(\hbar\omega/k_{\mathrm{B}}T_{i})-1]$
being the Bose–Einstein distribution, $k$ the parallel component of the
wavevector and $\omega$ the angular frequency of each mode. The Landauer-like
expression (2) of the heat flux has three ingredients, namely the energy of
each mode $\hbar\omega$, the statistics $\Delta
n_{\mathrm{B}E}(\omega,T_{1},T_{2})$ and the transmission probability (taking
values between 0 and 1) given by

$\displaystyle\mathcal{T}^{\mathrm{(}ph)}(\omega,k,d)=$ (4)
$\displaystyle\frac{4b^{2}\mathrm{Im}M_{1}(k,\omega)\mathrm{Im}M_{2}(k,\omega)}{|[1-aM_{1}(k,\omega)][1-aM_{2}(k,\omega)]-b^{2}M_{1}(k,\omega)M_{2}(k,\omega)|^{2}},$

where

$M_{i}(k,\omega)=\frac{\mathrm{i}}{\rho
c_{\mathrm{t},i}^{2}}\left(\frac{\omega}{c_{\mathrm{t},i}}\right)^{2}\frac{k_{\mathrm{l},i}}{(k_{\mathrm{l},i}^{2}-k^{2})^{2}+4k^{2}k_{\mathrm{t},i}k_{\mathrm{l},i}}$
(5)

is the mechanical susceptibility [38] associated with each medium, in which we
have defined a longitudinal acoustic wavevector
$k_{\mathrm{l},i}=\sqrt{(\omega/c_{\mathrm{l},i})^{2}-k^{2}}$ and a transverse
acoustic wavevector
$k_{\mathrm{t},i}=\sqrt{(\omega/c_{\mathrm{t},i})^{2}-k^{2}}$, related to the
longitudinal ($c_{\mathrm{l},i}$) and transverse ($c_{\mathrm{t},i}$) speeds
of sound in the media. The mechanical susceptibility $M_{i}(k,\omega)$ can be
derived by solving the elastic wave equations for a given mode
($k,\omega$)-mode in tensor form [38] and extracting the normal-normal
component to the surface (the only relevant component for plane–plane
symmetry). The integration bounds with respect to $k$ appearing in Eq. (2) are
limited by the function $M_{i}(k,\omega)$, which vanishes above
$\omega/c_{\mathrm{m}in}$, $c_{\mathrm{m}in}$ being the smallest transverse
speed of sound between the two media. In fact, evanescent transverse acoustic
waves do not contribute to the heat flux in the elastic model. Expression (4)
for the transmission probability stems from a fluctuational approach
analoguous to the one employed in the context of fluctuational electrodynamics
(a derivation can be found in App. C), where the currents are replaced by
surface displacements and $M_{i}(k,\omega)$ plays a similar role to that of
the dielectric permittivity in FED. Note that, contrary to the expression for
the radiative heat flux (see Sec. II.2), the transmission probability in the
case of phonons (4) cannot be separated in terms of phonon polarizations, as
longitudinal and transverse waves are coupled by the equations of motion.
Supertransverse waves (having a polarization vector parallel to the surface)
are not expected to contribute to the heat transfer [34, 35, 36, 37].

For metals, the terms $a$ and $b$ in the transmission probability (4) are
functions of $k$ and $d$, related to the van der Waals force, proportional to
the Hamaker constant $H$, and to the electrostatic forces through the bias
voltage, yielding

$a=\frac{H}{2\pi
d^{4}}+\epsilon_{0}\left(\frac{V_{\mathrm{b}}}{d}\right)^{2}k\coth(kd),$ (6a)
and $b=\frac{H}{4\pi
d^{2}}k^{2}K_{2}(kd)+\epsilon_{0}\left(\frac{V_{\mathrm{b}}}{d}\right)^{2}k\operatorname{csch}(kd),$
(6b)

where $K_{2}(x)$ is the second modified Bessel function of second kind and
$\operatorname{csch}$ represents the hyperbolic cosecant. The dependence on
$V_{\mathrm{b}}$ in Eq. (6) is independent from the van der Waals force and
both $a$ and $b$ are symmetric with respect to a change of sign of
$V_{\mathrm{b}}$. In the absence of bias voltage ($V_{\mathrm{b}}=0$) the
transmission probability in Eq. (4) is very sensitive to the difference
between $a$ and $b$.

An upper bound to the phononic contribution to the heat flux can be obtained
by setting the transmission to 1 in the expression of the flux, obtaining a
Stefan-Boltzmann-like law. This approximation yields
$\Phi^{(\mathrm{p}h)}_{\mathrm{m}ax}=(c/c_{\mathrm{min}})^{2}\sigma_{\mathrm{S}B}(T_{1}^{4}-T_{2}^{4})$
for bodies of the same material, where $c$ is the speed of light and
$\sigma_{\mathrm{S}B}=k_{\mathrm{B}}^{4}\pi^{2}/60\hbar^{3}c^{2}$ is the
Stefan–Boltzmann constant. This bound corresponds to a maximal conductance of
about $100$ GW m-2 K-1 at 300 K for gold semi-infinite media.

The phononic heat flow between gold media is shown in Fig. 2 as a function of
gap distance $d$ and bias $V_{\mathrm{b}}$.

Figure 2: Phonon heat flux between two semi-infinite gold slabs for various
values of the bias voltage $V_{\mathrm{b}}$ as a function of the separation
distance $d$. Inset: Phonon heat flux between two semi-infinite slabs for
$d=1$ nm as a function of $V_{\mathrm{b}}$ with a temperature difference of
100 K.

For $d\to 0$, the curves get closer: this results from the fact that the van
der Waals interaction becomes dominant and thus the heat flow is less impacted
by the bias voltage. Above 1 nm, the dependence of the bias voltage becomes
apparent, however the phonon heat flow decreases very rapidly for larger
distances. The slope of the curves for $d>1$ nm is $V_{\mathrm{b}}$-dependent.

We conclude this section by mentioning that the results obtained within a
continuous elastic model for phonons are significantly smaller than the ones
recently obtained by means of atomistic calculations, which include non-
equilibrium Green function methods of Refs. [26, 27] and molecular dynamics in
Ref. [28]. A comparison is given in App. D. However, the limit of validity of
the continuous approach to describe the phonon tunneling still remains today
an open problem which should be analyzed in detail in future works.

###  II.2 Photons

The theory describing near-field radiative heat transfer was originally
treated in the work of D. Polder and M. Van Hove [1]. The theoretical
framework they developed, based on FED, predicts a possibly unbounded flux,
which can overcome the blackbody limit
$\Phi_{\mathrm{B}B}=\sigma_{\mathrm{S}B}(T_{1}^{4}-T_{2}^{4})$ (corresponding
to a thermal conductance of 6 W$\cdot$m${}^{-2}\cdot$K-1 at 300 K) by several
orders of magnitude. More specifically, the first theoretical results
predicted a $d^{-2}$ divergence of the heat flux for subnanometer distances.
It was later suggested that such divergence in metals was an artefact
resulting from assuming a local dielectric function and could be corrected by
taking into account the nonlocal response by accounting for electron-electron
interactions [40, 41, 42].

The expression for the heat flux carried by electromagnetic waves is analogous
to the phonon contribution of Eq. (2). According to FED it is given by

$\displaystyle\Phi^{\mathrm{(}rad)}(T_{1},T_{2},d)=$ (7)
$\displaystyle\int_{0}^{\infty}\\!\frac{\mathrm{d}\omega}{2\pi}\hbar\omega\Delta
n_{\mathrm{B}E}(\omega,T_{1},T_{2})\int_{0}^{\infty}\\!\frac{\mathrm{d}k}{2\pi}\;k\\!\sum_{\alpha=\mathrm{s},\mathrm{p}}\mathcal{T}^{\mathrm{(}rad)}_{\alpha}(k,\omega,d),$

where the transmission probability

$\begin{split}&\mathcal{T}^{\mathrm{(}rad)}_{\alpha}(k,\omega,d)\\\
&\,=\begin{cases}\displaystyle\frac{(1-|r_{\alpha,1}|^{2})(1-|r_{\alpha,2}|^{2})}{|1-r_{\alpha,1}r_{\alpha,2}\exp(2\mathrm{i}k_{z}d)|^{2}},&k<\omega/c,\\\
\vspace{-0.3cm}\\\
\displaystyle\frac{4\,\mathrm{Im}\,r_{\alpha,1}\mathrm{Im}\,r_{\alpha,2}\exp(-2\,\mathrm{Im}\,k_{z}d)}{|1-r_{\alpha,1}r_{\alpha,2}\exp(-2\,\mathrm{Im}\,k_{z}d)|^{2}},&k\geq\omega/c,\\\
\end{cases}\end{split}$ (8)

can now be separated in terms of the two polarizations, given by the
transverse electric ($\alpha=\mathrm{s}$) and transverse magnetic
$(\alpha=\mathrm{p})$ contributions, where
$k_{z}=\sqrt{(\omega/c)^{2}-k^{2}}$. This expression takes into account both
propagative ($k<\omega/c$) and evanescent ($k>\omega/c$) waves.

Including extreme near-field effects, the reflection coefficients for metals
can be written as [42, 36]

$r_{\mathrm{s},i}(k,\omega)=\frac{\displaystyle
Z_{\mathrm{s},i}(k,\omega)-\frac{\omega}{c^{2}k_{z}}}{\displaystyle
Z_{\mathrm{s},i}(k,\omega)+\frac{\omega}{c^{2}k_{z}}},$ (9a)
$r_{\mathrm{p},i}(k,\omega)=\frac{\displaystyle\frac{k_{z}}{\omega}-Z_{\mathrm{p},i}(k,\omega)+\mathrm{i}\frac{\epsilon_{0}k^{2}}{\omega}\left(\frac{V_{\mathrm{b}}}{d}\right)^{2}M_{i}(k,\omega)}{\displaystyle\frac{k_{z}}{\omega}+Z_{\mathrm{p},i}(k,\omega)-\mathrm{i}\frac{\epsilon_{0}k^{2}}{\omega}\left(\frac{V_{\mathrm{b}}}{d}\right)^{2}M_{i}(k,\omega)},$
(9b)

where $Z_{\alpha,i}$ are the impedances [40, 42] given by

$\displaystyle
Z_{\mathrm{s},i}(k,\omega)=\frac{2\mathrm{i}}{\pi\omega}\int_{0}^{\infty}\frac{\mathrm{d}q_{z}}{\epsilon_{\mathrm{t},i}(K,\omega)-(cK/\omega)^{2}},$
(10a) $\displaystyle Z_{\mathrm{p},i}(k,\omega)$
$\displaystyle\,=\frac{2\mathrm{i}}{\pi\omega}\int_{0}^{\infty}\frac{\mathrm{d}q_{z}}{K}\left[\frac{q_{z}^{2}}{\epsilon_{\mathrm{t},i}(K,\omega)-(cK/\omega)^{2}}+\frac{k^{2}}{\epsilon_{\mathrm{l},i}(K,\omega)}\right],$
(10b)

where $K^{2}=k^{2}+q_{z}^{2}$ and $\epsilon_{\mathrm{l},i}(K,\omega)$ and
$\epsilon_{\mathrm{t},i}(K,\omega)$ are explicitly defined in App. B.

The mechanical susceptibility appears in the expression (9b) for p-polarized
waves and not in the expression (9a) for s-polarized waves. This modification
comes from the acoustic waves in the presence of a bias voltage, as the
oscillation of the displacements of the charged surfaces modulates the
radiative response [36, 37]. As in the case of phonons, the dependence on the
bias voltage of Eq. (9b) is symmetric with respect to the sign of the bias
voltage. When the dielectric function is nonlocal and anisotropic as in the
case of a gas of interacting electrons (metals), one has to distinguish
between longitudinal $\epsilon_{\mathrm{l},i}(\mathbf{K},\omega)$ and the
transverse $\epsilon_{\mathrm{t},i}(\mathbf{K},\omega)$ responses which are
functions of both the three-dimensional wavevector $\mathbf{K}$ and the
angular frequency $\omega$ of the modes. Explicit expressions are obtained
from Lindhard theory for jellium in the time-relaxation approximation and are
provided in App. B. In this work, we also consider the two common limits of
Lindhard expressions, the local case and the static case, presented in
subsection II.3, Eq. (16). Within the local isotropic case, we use the
expression given by Drude model as

$\epsilon_{\mathrm{l}oc}(\omega)=\epsilon_{\infty}-\frac{\omega^{2}_{\mathrm{p}l}}{\omega(\omega+\mathrm{i}\Gamma)},$
(11)

where $\omega_{\mathrm{p}l}$ is the plasma frequency, $\Gamma$ is the damping
coefficient and $\epsilon_{\infty}$ is a constant.

For large gap distances and in the case where there is no bias voltage
$V_{\mathrm{b}}=0$ and
$\epsilon_{\mathrm{l},i}(K,\omega),\epsilon_{\mathrm{t},i}(K,\omega)\to\epsilon_{\mathrm{l}oc}(\omega)$,
we recover the usual Fresnel reflection coefficients.

The relevance of the nonlocal contribution is shown in Fig. 3, where the
nonlocal contribution (dashed line) is compared to the local contribution as a
function of the gap distance $d$.

Figure 3: Local (dashed) and nonlocal (solid) radiative heat flux between two
semi-infinite gold slabs as a function of the separation distance $d$.

It is known that for metals, the local contribution diverges for small
distances. The nonlocal contribution fixes this issue, as the evanescent
contribution from the s-polarized waves saturates close to contact [42].
Nevertheless, there is a region below 1 nm, where the nonlocal contribution
exceeds the local contribution, as shown in figure. The difficulty to
experimentally observe this increase are due to the geometry of experiments,
where a tip–plane configuration is usually preferred. We discuss the effects
of geometry in Sec. III.

The effects of the bias voltage on the radiative contribution were studied in
Refs. [36, 37]. While the nonlocal contribution boosted the contribution of
the s-polarized waves for a range of gap distances, the bias voltage increases
the contribution of p-polarized waves for gaps smaller than 1 nm, as shown in
Fig. 4.

Figure 4: Radiative heat flux between two semi-infinite gold slabs with
respect to the separation distance $d$ for various values of the bias voltage
$V_{\mathrm{b}}$. Inset: radiative heat flux between two semi-infinite gold
slabs with respect to $V_{\mathrm{b}}$ for $d=1$ nm.

For distances above 1 nm, the standard FED results are recovered. The increase
of the radiative contribution decays slower than the phonon contribution and
it is also bias dependent.

###  II.3 Electrons

In the experiments described in Refs [23, 24], electronic currents are
measured at separation distances in the nanometer range. These currents are
related to tunneling electrons that can jump from the tip to the sample due to
both temperature and voltage biases. The standard expression for the electric
tunnelling current density between two semi-infinite media is given by [43,
29]

$\displaystyle J=-\frac{em_{\mathrm{e}}}{2\pi^{2}\hbar^{3}}$
$\displaystyle\int_{\tilde{E}}^{\infty}\mathrm{d}E_{z}\int_{0}^{\infty}\mathrm{d}E_{\perp}$
(12) $\displaystyle\times\Delta
n_{\mathrm{F}D}(E,T_{1},T_{2},\mu_{1},\mu_{2})\mathcal{T}^{\mathrm{(}el)}(E_{z},V_{\mathrm{b}})$

where $-e$ is the electron electric charge, $m_{\mathrm{e}}$ is the mass of
the elecrtron, $E=E_{\perp}+E_{z}$ is the total kinetic energy of an electron
decomposed in contributions stemming from velocities perpendicular and
parallel to the surface, $\tilde{E}=\max(0,-eV_{\mathrm{b}})$ and

$\begin{split}\Delta
n_{\mathrm{F}D}(E,T_{1},T_{2},\mu_{1},\mu_{2})&=n_{\mathrm{F}D}(E,T_{1},\mu_{1})\\\
&\quad-n_{\mathrm{F}D}(E,T_{2},\mu_{2}),\end{split}$ (13)

$n_{\mathrm{F}D}(E,T_{i},\mu_{i})=1/[\exp([E-\mu_{i}]/k_{\mathrm{B}}T_{i})+1]$
being the Fermi-Dirac distribution that depends on both temperature $T_{i}$
and chemical potential $\mu_{i}$.

Close to zero temperature, assuming a degenerate electron gas limit, the
current density leaving one body is bounded (in the case of two identical
metals) by
$J_{\mathrm{m}ax}=em_{\mathrm{e}}E_{\mathrm{F}}^{2}/4\pi^{2}\hbar^{3}$ (about
$2.4\times 10^{15}$ A/m2 for gold), where $E_{\mathrm{F}}$ is the Fermi
energy.

The transmission probability $\mathcal{T}^{(\mathrm{e}l)}(E,V_{\mathrm{b}})$
in Eq. (12) has only a few analytical solutions for specific electronic
barrier shapes. The semiclassical formula given by the
Wentzel–Kramers–Brillouin (WKB) method [44] allows us to estimate the
transmission probability through a smoothly varying electronic barrier as

$\displaystyle\mathcal{T}^{\mathrm{(}el)}_{\mathrm{W}KB}$
$\displaystyle(E_{z},V_{\mathrm{b}})=$ (14)
$\displaystyle\exp\left(-\frac{2\sqrt{2m_{\mathrm{e}}}}{\hbar}\int_{z_{1}}^{z_{2}}\mathrm{d}z\;\sqrt{U(z,V_{\mathrm{b}})-E_{z}}\right),$

where the integration is usually performed between the zeroes of the integrand
$z_{1}$ and $z_{2}$, in the region where the electronic barrier height
$U(z,V_{\mathrm{b}})$ is larger than the energy $E_{z}$. We define the
$z$-axis as shown in Fig. 1 and $z=0$ is located at the interface between
media 1 and vacuum. For the temperature of interest, we can safely neglect the
temperature dependence of the chemical potential and identify
$\mu_{1}=E_{\mathrm{F}}$. In the case of two metallic media, electrons need to
overcome a work function to leave the metal but also the forces generated by
the induced image charges in the surfaces of the two media. Using the image
method, the solution to (local) Poisson’s equation for an electron between two
identical perfect metals within a classical approach is given by [43]

$U_{\mathrm{c}l}(z)=E_{\mathrm{F}}+W+\frac{e^{2}}{16\pi\epsilon_{0}d}\left[\Psi(z/d)+\Psi(1-z/d)+2\gamma\right],$
(15)

where $W$ is the average work function of the media, $\epsilon_{0}$ is the
vacuum permittivity, $\Psi(z)$ is the digamma function and $\gamma$ is the
Euler-Mascheroni constant. The WKB method in combination with the classical
image potential has been used many times to describe tunnelling phenomena [43,
29, 49, 25, 26]. Nevertheless this methodology, when applied to a plane–plane
geometry, is troublesome in three aspects: (i) the work function has to be
determined experimentally or through ab-initio calculations, (ii) the
classical image potential is divergent for $z=0$ and $z=d$, far from the
assumptions of a smoothly varying potential, and (iii) the one-dimensional
$1/z$ divergences lead to exactly impenetrable barriers
($\mathcal{T}^{(\mathrm{e}l)}=0$) [45]. The last two difficulties are usually
avoided by translating the image planes inside the metals by a distance of a
few ansgtroms based on density functional theory (DFT) calculations.

Bardeen applied the Hartree-Fock method to self-consistently calculate the
electron density between two metals, which provided a smooth and continuous
electronic barrier [46]. In order to screen the effects of external charges,
electrons create the so called ”double layer” where a fraction of the
electrons leak into the vacuum, while leaving an effective positive charged
background near the surface of the metal.

The only way to obtain a regularized electronic barrier consists in
introducing electron-electron interactions, as in Lindhard theory of metals
that was already invoked for the nonlocal effects of the photonic
contribution. For the case of electrons, Kohn-Sham calculations show that
accounting just for exchange interactions leads to a slight overestimation of
the electronic barrier in the case of a jellium and one would need to self-
consistently add correlation functionals to recover the result of the DFT
potential [47]. In this work we limit ourselves to provide an estimation of
the electronic heat flux with a regularized potential using an approximate
solution ignoring the effects of the correlation potential. In the static long
wavelenght limit $\omega=0$, we can recover Thomas–Fermi screening theory
[32], where Lindhard dielectric functions approximate to

$\epsilon_{\mathrm{T}F}(K)=1+\frac{k_{\mathrm{T}F}^{2}}{K^{2}},$ (16)

where $K$ is the magnitude of the wavevector inside the metal,
$k_{\mathrm{T}F}=\sqrt{e^{2}m_{\mathrm{e}}k_{\mathrm{F}}/\pi^{2}\hbar^{2}\epsilon_{0}}$
is the inverse of Thomas–Fermi screening length (about 17 nm-1 for gold), and
$k_{\mathrm{F}}=\sqrt{2mE_{\mathrm{F}}/\hbar^{2}}$ is the Fermi wavevector
(about 12 nm-1 for gold).

Using a Green’s functions method, it has been shown that it is possible to
solve nonlocal Poisson’s equation in a jellium model in agreement with
Thomas–Fermi screening (TFA) [50, 51]. The full electronic potential used in
this work is given by

$U(z,V_{\mathrm{b}})=\left\\{\begin{matrix}\displaystyle U_{1}(z),&z\leq 0\\\
\vspace{-0.1cm}&\\\ \displaystyle
U_{\mathrm{g}ap}(z)-\frac{z}{d}eV_{\mathrm{b}},&0<z<d\\\ \vspace{-0.1cm}&\\\
\displaystyle U_{2}(z)-eV_{\mathrm{b}},&z\geq d\end{matrix}\right.,$ (17)

where the expression for $U_{s}(z)$ ($s=1,2,\mathrm{gap}$) at
$V_{\mathrm{b}}=0$ are provided in terms of integrals detailed in App. E. In
Eq. (17) we have added the effects of the bias voltage linearly inside the gap
region, approximation valid for small $V_{\mathrm{b}}<E_{\mathrm{F}}$ [50].
The regularized potential is depicted in Fig. 5 for $V_{\mathrm{b}}=1$ V,
together with the classical image potential of Eq. (15), for two different
separation distances.

Figure 5: Electronic potential barrier $U(x)$ normalized by the Fermi energy
$E_{\mathrm{F}}$ between two semi-infinite gold slabs with respect to the
position $x$ (in units of the gap thickness $d$) for a bias voltage
$V_{\mathrm{b}}=1$ V. The dashed lines represent the classical image
potential, the solid lines represent the regularized potential in the Thomas-
Fermi approximation. We set the work function of Eq. (15) to
$W=k_{\mathrm{T}F}e^{2}/8\pi\epsilon_{0}-E_{\mathrm{F}}$ in order to have the
same energy reference.

Through this method we obtain an electronic barrier of about $2E_{\mathrm{F}}$
for gold and large distances (as measured from the bottom of the band) which
automatically accounts for the work function. Below 1 nm the height of the
barrier reduces with respect to both the long distance value and the
prediction of the classical image potential (15). The potential reaches
different values inside the metal in accordance to the difference in chemical
potential due to the bias voltage. In Fig. 5 we also show the comparison
between two barriers of different gap distances ($d=2$ nm in black and $d=0.5$
nm in red), which shows the effect of reduction of the barrier near contact.
In the experiments this effect of a deviation from the large distances work
function was observed in Ref. [24].

In order to obtain a more accurate value of the transmission probability, we
drop the WKB approximation and calculate
$\mathcal{T}^{\mathrm{(}el)}(E,V_{\mathrm{b}})$ using a transfer matrix
method, based on Ref. [48]. We take slices of the potential and associate a
rectangular barrier to each slide with the height of the potential as a
function of the coordinate $z$. The total transmission probability is then
calculated as the product of the transmission matrices of all the slices. This
approach has the advantage of taking into account the features of the
electronic barrier inside and outside the metal, including the difference in
chemical potential due to the bias voltage.

In Fig. 6 we compare the current density between two gold semi-infinite media
obtained using the WKB method (14) and a transfer matrix method for the
classical image potential of Eq. (15) in the presence of bias voltage and the
TFA regularized potential of Eq. (17).

Figure 6: Absolute value of the electronic current $J$ (in units of the
saturation current $J_{\mathrm{m}ax}$ and in A/cm2) between two semi-infinite
gold slabs as a function of the bias voltage $V_{\mathrm{b}}$. Dashed lines
are calculated using WKB approximation, solid line employs the transfer matrix
method. Black line is calculated under the influence of the classic image
potential, red lines are calculated for the regularised potential of Eq. (17)
within the Thomas–Fermi approximation.

The dashed lines represent the current using WKB approximation for the
transmission probability of Eq. (14), while the solid lines represent the
current using the transfer matrix method. Due to the difficulties discussed
above, the transfer matrix method does not converge when using the classical
image potential (not shown). For small distances the classical potential is
significantly smaller in height compared to the TFA potential, thus the
classical potential overestimates the current (black dashed lines) for small
distances. Comparing the current density using the TFA potential, it can be
shown that the WKB calculation (red dashed) underestimates the current
compared to the transfer matrix method (red solid). Contrary to the WKB method
that is only sensitive to the shape of the higher part of the barrier, the
transfer matrix method has the feature of integrating over the whole shape of
the potential, including the difference of potential at the interior of the
metals which increases the transmission probability. The dependence of the
current on the separation distance $d$ is shown in Fig. 7 for different
applied bias voltages $V_{\mathrm{b}}$. We note, as expected, a strong
dependence on the applied bias. Moreover, the current is different from zero
even for $V_{\mathrm{b}}=0$, for which it is entirely due to the temperature
difference. The latter also explains the asymmetry between positive and
negative biases, which becomes less and less pronounced when increasing the
absolute value of $V_{\mathrm{b}}$.

Figure 7: Electronic current $J$ (in units of the saturation current
$J_{\mathrm{m}ax}$ and in A/cm2) between two semi-infinite gold slabs as a
function of the separation distance $d$ and for different bias voltages
$V_{\mathrm{b}}$. Dashed lines correspond to negative $V_{\mathrm{b}}$.

As the number density of electrons in Eq. (13) depends on the chemical
potential, we need to define the electronic heat flux carefully. Contrary to
the photons and phonons, for electrons there is a distinction between heat and
energy in the thermodynamic sense, and different definitions for the
electronic heat flux exist in literature [52, 53, 54, 26]. In this work, we
write the heat flux leaving media $i=1,2$, following the work of Ref. [29], as

$\displaystyle\Phi^{(\mathrm{el})}_{i}$
$\displaystyle(T_{1},T_{2},d)=\frac{m_{\mathrm{e}l}}{2\pi^{2}\hbar^{3}}\int_{\tilde{E}}^{\infty}\mathrm{d}E_{z}\int_{0}^{\infty}\mathrm{d}E_{\perp}$
(18) $\displaystyle\times(E-\mu_{i})\Delta
n_{\mathrm{F}D}(E,T_{1},T_{2},\mu_{1},\mu_{2})\mathcal{T}^{\mathrm{(}el)}(E_{z},V_{\mathrm{b}}),$

similar to Eqs. (2) and (7) for phonons and photons, respectively, where the
statistics is now fermionic and the energy in the integrand is replaced by the
heat $E-\mu_{i}$ associated with each reservoir. The equations for the
electronic heat flux are no longer symmetrical and one has to distinguish the
heat fluxes for each media. Note that Eq. (18) ensures that the total power
delivered to the system
$\Phi^{(\mathrm{e}l)}_{1}-\Phi^{(\mathrm{e}l)}_{2}=JV_{\mathrm{b}}$ equals the
power supplied by the external generator employed to keep the bias voltage
constant. At zero bias and equal chemical potentials $\mu_{1}=\mu_{2}$,
$\Phi^{(\mathrm{e}l)}_{2}(V_{\mathrm{b}}=0)=\Phi^{\mathrm{(}el)}_{1}(V_{\mathrm{b}}=0)$
as in the usual thermodynamic sense, where the hotter body at $T_{1}$ tends to
let electrons go, cooling down and heating the cooler body at $T_{2}$.
Assuming that the emitting body replaces the electrons of energy $E$ with new
ones with energy $E_{\mathrm{F}}$, the presence of the leads makes the
difference for $V_{\mathrm{b}}\neq 0$ as
$\Phi^{(\mathrm{e}l)}_{2}(V_{\mathrm{b}})\neq\Phi^{(\mathrm{e}l)}_{1}(V_{\mathrm{b}})$
in general, allowing for situations where the signs are opposite and both
bodies would have the tendency to heat up. This dependence on both temperature
and bias voltage leads to inversion effects. In thermoionic applications, the
Nottingham effect, where the heat flux changes sign at large temperatures
[29], has been explored for refrigeration [49]. Notice that the electronic
heat flux is bounded by
$\Phi_{\mathrm{e}l}^{(\mathrm{m}ax)}=J_{\mathrm{m}ax}E_{\mathrm{F}}/3e$. This
value corresponds to about $4\times 10^{15}$ W/m2 for Au-Au tunneling.

In Fig. 8 we compare the absolute value of the heat flux on the hotter medium
1 for different values of the bias voltage and gap distances.

Figure 8: Absolute value of the electronic heat flux on body 1 in plane–plane
configuration for various values of the bias voltage $V_{\mathrm{b}}$ as a
function of the separation distance $d$. The solid (dashed) lines indicate
emitted (receiving) heat from medium 1.

For large values of $V_{\mathrm{b}}$ the flux saturates. The heat flux for
$V_{\mathrm{b}}=0$ is in the thermodynamic regime where body 1 tends to cool
down by the emission of the electrons, and has the opposite sign of the curves
with $V_{\mathrm{b}}>0$ (solid lines) which represent heating of body 1. We
stress that our definition (18) of the electronic heat flux differs from the
one employed in Ref. [26]. A comparison of the numerical numerical results
obtained within the two approaches is given in App. D.

###  II.4 Comparison

The comparison of the different contribution to the heat flux (ph, rad, el) is
presented in Fig. 9 with a double logarithmic scale, with positive (negative)
values representing the tendency to cool down (heat up) body 1.

Figure 9: Comparison of the heat flux of the different heat carriers in
plane–plane configuration for various values of the bias voltage
$V_{\mathrm{b}}$ as a function of the separation distance $d$. Positive
(negative) values of the heat flux are related to outgoing (incoming) heat on
body 1.

The near-field radiative contribution and the phonon contribution are both
positive for any separation distance and applied bias voltage, corresponding
to heat emitted by body 1. The comparison between them shows that the phonon
contribution dominates over the near-field radiative one below 0.4 nm. The
electronic heat flux is also positive when there is no bias voltage and
changes drastically of sign above a few tens of mV. The solid lines represent
the total heat flux as the sum of the three contributions (ph, rad, el). We
observe that, while for distances larger than 1 nm photons are the only
relevant carriers, the electronic contribution rapidly dominates below 1 nm.
These results represent a clear deviation from FED. The fact that the
electronic contribution dominates the heat transfer is no coincidence as the
main contribution to the thermal conductivity of metals is electronic in
nature (phonons contribute about 10%). Near contact (distances around 1 Å),
the phononic part becomes relevant compared to the electronic contribution as
expected.

In Fig. 10 we map the dominating contributions to the total heat $\Phi_{1}$
flowing from the hotter medium 1 as a function of $d$ and $V_{\mathrm{b}}$.

Figure 10: Comparison of the absolute value of the heat flux $\Phi_{1}$
emitted by media 1 in plane–plane configuration for various values of the bias
voltage $V_{\mathrm{b}}$ and separation distance $d$ for $T_{1}=400$ K and
$T_{2}=300$ K. The different colors represent the leading heat carriers (ph in
red, rad in yellow, el in blue). Regions with two labels (el-rad in green, ph-
rad in purple, el-ph in orange) represent regions where the contribution of
the two carriers are comparable up to $10\%$ between the two carriers. The
label ”all” (in pink) represents the regions where all three carriers
contribute equally up to $10\%$ difference. The dashed lines separate the
regions where the electronic heat flux has a different sign (heating in the
external regions enclosed by the white dashed lines, cooling in the rest of
the figure).

Each label in the colored regions of the figure indicates the dominating heat
carrier (ph, rad, el). In some regions, two (el-rad, ph-rad, el-ph) or even
the three (all) carriers compete. The two regions close to the borders of the
figure enclosed by the white dashed lines indicate the regions where the heat
flux becomes negative (incoming heat). As expected, at about 0.8 nm radiation
clearly dominates and continues to do for larger distances. Below 0.5 nm, the
region of dominance of the electrons fills the map. Near the change of sign
(dashed white line) and about $\pm 40$ mV, the electronic contribution is
reduced and the phonons create thin stripes where their contribution competes
with the electrons. Near contact phonons compete with the electrons, as
expected, or even become the dominant contribution, but only as long as the
bias remains small, i.e. $|V_{\mathrm{b}}|<100$ mV. The figure is slightly
asymmetric with respect to the zero of the bias voltage due to the temperature
gradient. We do not map the heat flowing into media 2 but it is expected to be
similar aside from the change in bias voltage and no cooling regions.

##  III Tip–plane configuration

In order to compare our results with the existing experiments, we exploit the
Dejarguin or proximity force approximation (PFA), already employed in the
context of near-field radiative heat transfer by the authors of Refs. [23,
24]. For a spherical shaped tip of radius $R$, the net power exchanged between
the tip and the sample can be written as

$P_{\mathrm{t}ip}=(-1)^{i+1}(P^{\mathrm{(}ph)}+P^{\mathrm{(}rad)}+P^{\mathrm{(}el)}_{i}),$
(19a) in agreement with the net flux of Eq. (1), where each term is defined as
$P^{\mathrm{(}Q)}_{i}=2\pi\int_{0}^{R}\mathrm{d}r\;r\;\Phi^{\mathrm{(}Q)}_{i}(d+R-\sqrt{R^{2}-r^{2}}),$
(19b)

where $Q\in\\{\mathrm{rad,ph,el}\\}$ and $i=1,2$ is to be chosen depending on
which media corresponds to the tip and the subindex is dropped when discussing
phonons and radiation. Analogously, the PFA expression for the electric
current reads

$I=2\pi\int_{0}^{R}\mathrm{d}r\;r\;J(d+R-\sqrt{R^{2}-r^{2}})$ (20)

where $J$ is the tunneling current density defined in Eq.(12).

In the supplementary material of Ref. [23] it was already found that the
nonlocal effects were not observable in Au-Au experiments with tip–plane
configuration using the PFA. We provide a further evidence of this point in
Fig. 11, where we plot the power emitted by the tip as a function of the gap
by multiplying $\Phi^{\mathrm{(}rad)}$ by the cross-section of the tip $\pi
R^{2}$ (solid lines), and using the PFA (dashed) lines, in the presence of
local and nonlocal effects.

Figure 11: Comparison between local (red) and nonlocal (black) radiative power
emitted in tip–plane configuration as a function of the separation distance
$d$ for a spherical tip of radius $R=100$ nm. The dashed line is calculated
using the proximity force approximation, the solid lines represents the heat
flux multiplied by the tip cross-section.

When looking at the flux times the cross-section, we observe that the nonlocal
radiative heat flux increases with respect to the local theory in a region
below 1 nm, as we have already showed in Fig. 3. However, when looking at the
PFA curves in Fig. 11, we observe that the radiative flux is almost constant
and it is not possible to distinguish the nonlocal estimate from the local
one. The distance dependence of the radiative flux disappears by integrating
over the tip, which translates into a reduction of 1 to 3 orders of magnitude
with respect to the solid curves. The region of interest for the nonlocal
effects represent a small area at the edge of the tip and thus gives a
negligible contribution to the integral given in Eq. (19).

While the contribution of the bias voltage to the radiative heat flux can
exceed the influence of the nonlocal effects to the radiative transfer, its
effects are not striking either. For this purpose, we also compare the effects
of the bias voltage on the transmitted power in flux times cross-section and
PFA in Fig. 12, as function of gap distance and bias voltage (see Fig. 4).

Figure 12: Comparison between radiative power exchanged in a tip–plane
configuration for different values of the bias voltage $V_{\mathrm{b}}$ as a
function of the separation distance $d$ for a spherical tip of radius $R=30$
nm. The dashed lines represent the result obtained using the proximity force
approximation while the solid lines denote the heat flux multiplied by the tip
cross-section.

The quick ascent of the solid curves below 0.5 nm in the presence of the bias
with respect to zero bias translates also into almost no difference between
the PFA curves (dashed) and the difference between the curves is barely
noticeable for the largest bias shown (blue, 1 V) at vanishing gap distances.
Stronger biases could be used to probe this effect, but aside from being an
experimental challenge, the contribution of the other two carriers (ph and el)
increases as well and stronger couplings are expected.

In Figs. 13 and 14 we show the current and the total power emitted by the tip
calculated using PFA in a tip–plane scenario with a tip radius of $R=100\,$nm.
Note that in Fig. 14 the dips around 5 Å correspond to a sign flip of the
total power. We remark that both quantities are quite sensitive to the applied
bias and almost symmetric with respect to a sign change of $V_{\mathrm{b}}$.
Moreover, note that the tip is heated up for large enough $V_{\mathrm{b}}$ (in
absolute value) and for small distances.

Figure 13: Absolute value of the electric current exchanged as a function of
the separation distance in a tip–plane scenario, calculated within the PFA
approach, between a gold tip of radius $R=100\,$nm and at temperature
$T_{1}=400\,$K and a gold planar substrate at temperature $T_{2}=300\,$K. The
different curves correspond to different values of the applied bias voltage
$V_{\mathrm{b}}$. Dashed lines correspond to negative bias voltages and the
electronic current has the opposite sign compared to that of solid lines.
Figure 14: Absolute value of the total power emitted by the tip as a function
of the separation distance in a tip–plane scenario, calculated within the PFA
approach, between a gold tip of radius $R=100\,$nm and at temperature
$T_{1}=400\,$K and a gold planar substrate at temperature $T_{2}=300\,$K. The
curves at $V_{\mathrm{b}}=\pm 100\,\text{mV},\pm 1\,\text{V}$ switch to take
negative values below a separation distance below 1 nm.

Moving to theory-experiment comparison, we start noting that Kloppstech et al.
[23] modelled the tip as a sphere with a radius of about 30 nm, while Cui et
al. [24] had a tip of radius of 150 nm. Furthermore, the former experiment had
a positive temperature bias ($T_{1}=280$ K for the tip and $T_{2}=120$ K for
the sample), while in the latter set-up the bias was inverted ($T_{1}=445$ K
for the sample and $T_{2}=315$ K for the tip). In Fig. 15 we compare our
theoretical results stemming from PFA (solid and dashed curves) using Eqs.
(19) to the experimental results of Ref. [23] (symbols $\times$).

Figure 15: Comparison of the numerical results with the available experimental
data [23]. The curves indicate the power exchanged between the tip and the
plane with respect to the separartion distance $d$ for various values of the
bias voltage $V_{\mathrm{b}}$. Calculations are performed in the proximity
force approximation with a spherical tip of radius $R=30$ nm. The dips around
5 Å correspond to a sign flip of the total power.

While the phononic (black dashed line) and radiative (dash-dotted line)
contributions are shown for $V_{\mathrm{b}}=0$, the electronic ones are shown
for different bias voltages up to 1 V, including 0.6 V which is the value
applied in the experiment. The experimental data diverges from the FED
predictions already at 7 nm, much larger than the 1 nm scale where electrons
are expected to contribute. The authors reported that no electric current was
detectable above 1 nm (below 0.5 pA), indicating that, even if contamination
was present, it was not conductive as it did not enhance the electronic
contribution. However that does not exclude the possibility of an enhancement
of the bosonic contributions (ph and rad) or the activation of other heat
channels.

We now focus on the experimental results of Ref. [24], compared to our
theoretical predictions in Fig. 16.

Figure 16: Comparison of the numerical results with the available experimental
data [24]. The curves indicate the heat conductances in a tip-plane
configuration with respect to the separation distance $d$ with
$V_{\mathrm{b}}=0$. Calculations are performed in the proximity force
approximation with a spherical tip of radius $R=150$ nm. Experimental values
are indicated by the pink line while the pink background shows the
experimental error margin according to Fig. 2d of Ref. [24]

The figure shows the thermal conductance of the tip
$h_{\mathrm{t}ip}=P_{\mathrm{t}ip}/\Delta T$ under PFA for $V_{\mathrm{b}}=0$,
where $\Delta T=T_{1}-T_{2}$ is the temperature difference. The radiative
contribution (black dashed line) is in agreement with the experimental data
(dark pink solid line) and inside the error range (light pink background),
suggesting no deviation from FED. We also note that the error margin goes to
negative values in the experimental data (not shown in the figure). The
behaviour of the experimental data suggest an opposite conclusion with respect
to Kloppstech et al. [23] where there is a clear deviation from the radiative
heat transfer theory. Our theoretical predictions indicate that the electronic
contribution (blue solid line) should overcome the error range for 0.3 nm and
below in Fig. 16, which is not detected in the experiment, even if an
electronic current is measured at this range. The phonon contribution (red
dash-dotted curve), assuming $V_{\mathrm{b}}=0$, starts to be relevant at
distances where there are no experimental values. If a small bias of a few mV
was introduced in the experiment, the phonon contribution could be more
relevant at larger distances. A slight rise in the experimental curve (dark
pink solid) is seen going to smaller distances but the error margin does not
allow to draw any reliable quantitative conclusion.

##  IV Conclusion

In this work we have shown that, close to the contact, electrons and phonons
are the two main contributors to heat exchanges between two metallic solids,
provided that the bias voltage applied between them satisfies
$V_{\mathrm{b}}<100$ mV. For larger bias voltages electrons become dominant
and the Nottingham effect associated with these carriers gives rise to a rich
thermal behavior (including thermal-rectification effects). Hence, when
$V_{\mathrm{b}}>0$ the resistive heating associated to the electronic current
induces an asymmetry in the power received and emitted by the electrodes. In
this case the flux mediated by electrons flows in the same direction as the
flux mediated both by photons and phonons. Moreover, this flux becomes orders
of magnitude larger than the phononic flux. On the other hand, when
$V_{\mathrm{b}}<0$, the heat flux carried by electrons is of opposite sign
with respect to the flux carried by phonons and it tends to cool down the cold
metal.

The comparison of our theoretical predictions with the existing experimental
results show that some problems limit today their interpretation. In Reddy’s
experiment [24] we have seen that the value of heat flux is below the
experimental sensitivity. However, an increase of bias voltage or of
temperature gradient could probably make the study of heat transfer at atomic
scale possible. On the other hand, in Kittel’s experiment [23] the strong
spatial shift between the measured flux and the theoretical predictions at
distances of the order of $5\,$nm, where neither electrons nor phonons
contribute to the transfer, cannot be easily interpreted without considering
the presence of molecular contaminations or extra layers such as water layers.
Nevertheless, once again, a series of measurement made with various bias
voltages could probably help clarifying the current results.

Beside, we think that our work could be useful for future experimental studies
of Nottingham effect and to investigate its consequences on the heat exchanges
between two metallic solids close to the physical contact.

###### Acknowledgements.

This research was supported by the French Agence Nationale de la Recherche
(ANR), under grant ANR-20-CE05-0021-01 (NearHeat).

##  Appendix A Reference values for the constants

In the elastic acoustic wave model, we use $c_{\mathrm{l}}=3240$ m/s,
$c_{\mathrm{t}}=1200$ m/s and $\rho=1.9\times 10^{4}$ kg/m3 for gold. The van
der Waals force between two gold metals is proportional to the Hamaker
constant $H=34.76\times 10^{-20}$ J, taken from Ref. [39].

In Drude’s expression for the local dielectric function of Eq. (11) we use
$\omega_{\mathrm{p}l}=1.36\times 10^{16}$ rad/s as the plasma frequency for
gold, $\Gamma=1\times 10^{14}$ rad/s and $\epsilon_{\infty}=9.84$, fitted data
[30] from Ref. [31]. The Fermi energy is taken as $E_{\mathrm{F}}=5.53$ eV for
gold. For the nonlocal calculation, we use $v_{\mathrm{F}}=c/214$ from Ref.
[32].

##  Appendix B Nonlocal dielectric function

In this appendix we detail the longitudinal and transverse dielectric
functions in the rotating phase approximation (RPA), also known as Lindhard
theory for an interacting electron gas, employed in expressions of Sec. II.2
to describe the the nonlocal radiative heat transfer.

For low temperatures compared to the Fermi termperature $T_{1},T_{2}\ll
T_{\mathrm{F}}$, Lindhard theory together with Mermin’s time-relaxation
approximation [40, 42] describes the longitudinal and transverse dielectric
functions as

$\epsilon_{\mathrm{l}}(K,\omega)=\epsilon_{\infty}+\frac{3\omega_{\mathrm{p}l}^{2}}{\omega+\mathrm{i}\Gamma}\frac{\nu^{2}g_{\mathrm{l}}(\zeta,\nu)}{\omega+\mathrm{i}\Gamma\displaystyle\frac{g_{\mathrm{l}}(\zeta,\nu)}{g_{\mathrm{l}}(\zeta,0)}},$
(21a) and
$\displaystyle\epsilon_{\mathrm{t}}(K,\omega)=\epsilon_{\infty}-\frac{\omega^{2}_{\mathrm{p}l}}{\omega^{2}(\omega+\mathrm{i}\Gamma)}$
(21b)
$\displaystyle\times\\{\omega[g_{\mathrm{t}}(\zeta,\nu)-3\zeta^{2}g_{\mathrm{l}}(\zeta,\nu)]$
$\displaystyle+\mathrm{i}\Gamma[g_{\mathrm{t}}(\zeta,0)-3\zeta^{2}g_{\mathrm{l}}(\zeta,0)]\\},$

where $\zeta=K/2k_{\mathrm{F}}$ and
$\nu=(\omega+\mathrm{i}\Gamma)/Kv_{\mathrm{F}}$, $v_{\mathrm{F}}$ the Fermi
speed, and

$\displaystyle g_{\mathrm{l}}(\zeta,\nu)=\frac{1}{2}$
$\displaystyle+\frac{1-(\zeta-\nu)^{2}}{8\zeta}\log\left(\frac{\zeta-\nu+1}{\zeta-\nu-1}\right)$
(22a)
$\displaystyle+\frac{1-(\zeta+\nu)^{2}}{8\zeta}\log\left(\frac{\zeta+\nu+1}{\zeta+\nu-1}\right)$
and $\displaystyle g_{\mathrm{t}}(\zeta,\nu)=$
$\displaystyle\frac{3}{8}(\zeta^{2}+3\nu^{2}+1)$ (22b)
$\displaystyle-3\frac{[1-(\zeta-\nu)^{2}]^{2}}{32\zeta}\log\left(\frac{\zeta-\nu+1}{\zeta-\nu-1}\right)$
$\displaystyle-3\frac{[1-(z+\nu)^{2}]^{2}}{32\zeta}\log\left(\frac{\zeta+\nu+1}{\zeta+\nu-1}\right)$

are the Lindhard-Mermin functions.

##  Appendix C Fluctuational acoustodynamics

In the fluctuational electrodynamics theory, the components of the currents
due to thermal and quantum fluctuations are related to the dielectric function
through the fluctuation-dissipation theorem [36, 37]. Here we present an
acoustic analog for the elastic waves inside each media $i=1,2$. The
$z$-component of the oscillating displacements
$u_{i}(x,y,t)=u^{\mathrm{f}l}_{i}(x,y,t)+u^{\mathrm{i}nd}_{i}(x,y,t)$ at the
surface of each media are written as a sum of a fluctuating and induced
contributions. The induced displacements between the two surfaces
$u^{\mathrm{i}nd}_{i}$ originate from the van der Waals and electrostatic
forces. The fluctuation-dissipation theorem for the fluctuations of the
displacements $u_{i}^{\mathrm{f}l}$ inside each body reads [36, 37]

$\begin{split}\langle
u_{i}^{\mathrm{f}l}(\mathbf{k},\omega)[u_{j}^{\mathrm{f}l}(\mathbf{k}^{\prime},&\omega^{\prime})]^{*}\rangle=\hbar\;\mathrm{Im}M_{i}(k,\omega)\coth\left(\frac{\hbar\omega}{2k_{\mathrm{B}}T_{i}}\right)\\\
&\,\times(2\pi)^{3}\delta(\mathbf{k}-\mathbf{k}^{\prime})\delta(\omega-\omega^{\prime})\delta_{ij},\end{split}$
(23)

where $\delta(x)$ is the Dirac delta distribution, $\delta_{ij}$ is the
Kronecker delta and $M_{i}(k,\omega)$ are the mechanical susceptibilities for
each body defined in Eq. (5), and $(\mathbf{k},\omega)$-variables are related
to the real space displacements by a Fourier transform,

$u_{i}(x,y,t)=2\,\mathrm{Re}\int_{0}^{\infty}\frac{\mathrm{d}\omega}{2\pi}\int\frac{\mathrm{d}^{2}\mathbf{k}}{(2\pi)^{2}}u_{i}(\mathbf{k},\omega)\mathrm{e}^{\mathrm{i}(\mathbf{k}\cdot\mathbf{r}-\omega
t)}.$ (24)

The mechanical forces per unit surface $f_{i}$ acting on surfaces 1 and 2 are
obtained from derivatives of the interacting energy, integrating atom-wise
over the two surfaces [34, 35, 36, 37], given in terms of the displacements by

$f_{1}(\mathbf{k},\omega)=au_{1}(\mathbf{k},\omega)-bu_{2}(\mathbf{k},\omega),$
(25a)
$f_{2}(\mathbf{k},\omega)=au_{2}(\mathbf{k},\omega)-bu_{1}(\mathbf{k},\omega),$
(25b)

where $a$ and $b$ are defined in Eq. (6).

The mechanical susceptibilities $M_{i}(\omega,k)$ relate the induced component
of the displacements to the acting forces, such that ($i=1,2$)

$u_{i}(\mathbf{k},\omega)=u_{i}^{\mathrm{fl}}(\mathbf{k},\omega)+M_{i}(k,\omega)f_{i}(\mathbf{k},\omega).$
(26)

In order to obtain the heat flux carried by the phonons, we need to calculate
the mean power emitted by unit area given by

$\begin{split}\Phi_{\mathrm{p}h}&=\langle\dot{u}_{2}(x,y,t)f_{2}(x,y,t)\rangle-\langle\dot{u}_{1}(x,y,t)f_{1}(x,y,t)\rangle\\\
&=2\int_{0}^{\infty}\frac{\mathrm{d}\omega}{2\pi}\int\frac{\mathrm{d}^{2}\mathbf{k}}{(2\pi)^{2}}\int_{0}^{\infty}\frac{\mathrm{d}\omega^{\prime}}{2\pi}\int\frac{\mathrm{d}^{2}\mathbf{k}^{\prime}}{(2\pi)^{2}}\omega\\\
&\,\times\mathrm{Im}[\langle
u_{2}(\mathbf{k},\omega)f^{*}_{2}(\mathbf{k}^{\prime},\omega^{\prime})\rangle-\langle
u_{1}(\mathbf{k},\omega)f^{*}_{1}(\mathbf{k}^{\prime},\omega^{\prime})\rangle].\end{split}$
(27)

Solving eqs. (25) and (26), for $f_{i}(\mathbf{k},\omega)$ and
$u_{i}(\mathbf{k},\omega)$ in terms of
$u_{i}^{\mathrm{f}l}(\mathbf{k},\omega)$, allows to use the fluctuation-
dissipation theorem (23) together with (27) to recover the Eq. (2) of the
phonon flux found in the main text.

##  Appendix D Comparison with numerical results of Refs. [26] and [28]

In Fig. 17 we compare the phononic and electronic heat fluxs for
$T_{1}=280\,K$ and $T_{2}=120\,$K and zero applied bias obtained within our
formalism to the ones calculated in Refs. [26] (phonons and electrons) and
[28] (phonons only, results at 300 K) using atomistic calculations. The
comparison shows that our results can be various orders of magnitude smaller.
Concerning phonons, it has been argued [28] that this discrepancy is due to
the fact that the continuous model only takes into account the forces between
the elements of the surface and neglects the elements inside the material,
included in the atomistic model. However, as the molecular simulations
consider only a small number of atoms, it remains unclear if this prediction
hold for systems with larger dimensions. As for electrons, the discrepancy is
certainly due to the different definition of the heat flux, but could also be
due to simplified classical-image potential used in Ref. [26], along with the
WKB approximation.

Figure 17: Phononic and electronic thermal conductances
$h=\Phi_{1}/(T_{1}-T_{2})$ between two semi-infinite gold slabs as a function
of the separation distance at $V_{\mathrm{b}}=0$. The results for phonons
within the fluctuation-acoustodynamics approach based on continuous elastic
theory (red solid line) are compared to the atomistic approaches of Refs. [26]
(red dashed lines) and Ref. [28] (red crosses) at 300 K. The results we
obtained for electrons (blue solid line) are compared to the ones obtained
within WKB approximation in Ref. [26], using a different definition of the
electronic heat flux. For reference, the dotted black line corresponds to the
local radiative thermal conductance.

##  Appendix E Thomas–Fermi electronic barrier

In this appendix, we provide the explicit expression of the regularized
potential in Eq. (17) and provide a few steps for its derivation coming from
Refs. [50, 51].

In order to obtain the regularized electronic barrier for the plane–plane
configuration, it is necessary to solve one-dimensional nonlocal Poisson’s
equation in the presence of a polarizable media, i.e.

$\begin{split}\left(\frac{\partial^{2}}{\partial
z^{2}}-k^{2}\right)G(k;z,z^{\prime})-&\int\mathrm{d}z^{\prime\prime}\Pi(k;,z,z^{\prime})G(q;z^{\prime\prime},z^{\prime})\\\
&=\delta(z-z^{\prime}),\end{split}$ (28)

where $G(k;z,z^{\prime})$ is the Green function and $\Pi(k;z,z^{\prime})$ is
the polarization operator. In the case of three regions, allowing for specular
reflection and continuity of the potential, we can write the polarization
operator for a three layer system as

$\Pi(k;z,z^{\prime})=\begin{cases}\Pi_{1}(k;z-z^{\prime})+\Pi_{1}(k;z+z^{\prime}),\\\
\hskip 113.81102ptz,z^{\prime}\leq 0,\\\
\Pi_{2}(k;z-z^{\prime})+\Pi_{2}(k;z+z^{\prime}),\\\ \hskip
113.81102ptz,z^{\prime}\geq d,\\\
\Pi_{\mathrm{g}ap}(k;z-z^{\prime})+\Pi_{\mathrm{g}ap}(k;z+z^{\prime}),\\\
\hskip 113.81102pt0<z,z^{\prime}<d,\end{cases}$ (29a) being $\Pi_{s}(k;z\mp
z^{\prime})=\int_{-\infty}^{\infty}\frac{\mathrm{d}q_{z}}{2\pi}K^{2}[\epsilon_{s}(K)-1]\exp(\mathrm{i}q_{z}[z\mp
z^{\prime}]),$ (29b)

where $s=1,2$ corresponds to the regions inside media 1 and 2, and
$s=\mathrm{gap}$ is the region between the two media. Here we take the static
approximation $\epsilon_{s}(K)=\epsilon_{\mathrm{T}F}(K)$ for $s=1,2$ in Eq.
(16) and $\epsilon_{\mathrm{g}ap}=1$. In the case of classical metals of
vanishing screening length, i.e. $k_{\mathrm{T}F}\to\infty$, we recover the
classical image potential of Eq. (15).

By calculating the diagonal self-interacting element $G(k;z)=G(k;z,z)$, we can
recover the electronic barrier $U(x)$ in the absence of bias voltage
corresponding to $V_{\mathrm{b}}=0$ for the three different regions as

$U_{s}(x)=\frac{e^{2}}{4\pi\epsilon_{0}}\left\\{\frac{k_{\mathrm{T}F}}{2}-\int_{0}^{\infty}\mathrm{d}k\,k\left[G_{s}(k;z)+\frac{1}{2k}\right]\right\\},$
(30)

where we add the factor $e^{2}k_{\mathrm{T}F}/8\pi\epsilon_{0}$ to set the
zero energy reference inside the metal. In the case where there is no bias,
the equations are symmetric with respect to the center of the gap $x=d/2$, and

$\displaystyle G_{1}(k;z)=$
$\displaystyle\frac{a^{2}_{0}(k,z)}{B_{0}(k)}[a_{\mathrm{S}}(k,0)+a_{\mathrm{A}}(k,0)+2a_{0}(k,0)]$
(31a) $\displaystyle-\frac{a_{0}(k,0)+a_{0}(k,2z)}{2},$
$G_{2}(k;z)=G_{1}(k;d/2-z),$ (31b)
$\begin{split}G_{\mathrm{g}ap}(k;z)&=\frac{1}{2}\left[\frac{a^{2}_{\mathrm{S}}(k,z)}{a_{\mathrm{S}}(k,0)+a_{0}(k,0)}+\frac{a^{2}_{\mathrm{A}}(k,z)}{a_{\mathrm{A}}(k,0)+a_{0}(k,0)}\right]\\\
&\,-\frac{a_{\mathrm{A}}(k,0)+a_{\mathrm{A}}(k,2z)+a_{\mathrm{S}}(k,0)+a_{\mathrm{S}}(k,2z)}{4},\end{split}$
(31c)

where we have defined

$B_{0}(k)=2[a_{0}(k,0)+a_{\mathrm{A}}(k,0)][a_{0}(k,0)+a_{\mathrm{S}}(k,0)],$
(32a)
$a_{0}(k,z)=\frac{\exp(-\sqrt{k^{2}+k^{2}_{\mathrm{T}F}}|z|)}{\sqrt{k^{2}+k^{2}_{\mathrm{T}F}}},$
(32b)
$a_{\mathrm{S}}(k,z)=\frac{\cosh(kd/2-k\,\mathrm{mod}(z,d))}{k\sinh(kd/2)},$
(32c) and
$a_{\mathrm{A}}(k,z)=\operatorname{sign}\bigl{(}d-\mathrm{mod}(z,2d)\bigr{)}\frac{\sinh(kd/2-k\,\mathrm{mod}(z,d))}{k\cosh(kd/2)},$
(32d)

where $\operatorname{sign}(z)$ is the sign of $z$. Taking
$k_{\mathrm{T}F}\to\infty$ one recovers the expressions for the classical
image potential in Eq. (15).

For small biases compared to the chemical potential, we are allowed to add the
contribution of $V_{\mathrm{B}}$ linearly as in the expressions given in Eq.
(17) in the main text.

## References

  * [1] D. Polder and M. Van Hove, Theory of Radiative Heat Transfer between Closely Spaced Bodies, Phys. Rev. B 4, 3303 (1971).
  * [2] 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).
  * [3] A. I. Volokitin and B. N. J. Persson, Near-field radiative heat transfer and noncontact friction, Rev. Mod. Phys. 79, 1291 (2007).
  * [4] S.-A. Biehs, M. Tschikin, and P. Ben-Abdallah, Hyperbolic Metamaterials as an Analog of a Blackbody in the Near Field, Phys. Rev. Lett. 109, 104301 (2012).
  * [5] 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).
  * [6] L. Hu, A. Narayanaswamy, X. Chen, and G. Chen, Near-field thermal radiation between two closely spaced glass plates exceeding Planck’s blackbody radiation law, Appl. Phys. Lett. 92, 133106 (2008).
  * [7] S. Shen, A. Narayanaswamy, and G. Chen, Surface phonon polaritons mediated energy transfer between nanoscale gaps, Nano Lett. 9, 2909–2913 (2009).
  * [8] E. Rousseau, A. Siria, G. Jourdan, S. Volz, F. Comin, J. Chevrier, and J.-J. Greffet, Radiative heat transfer at the nanoscale, Nat. Photon. 3, 514 (2009).
  * [9] R. S. Ottens, V. Quetschke, Stacy Wise, A. A. Alemi, R. Lundock, G. Mueller, D. H. Reitze, D. B. Tanner, and B. F. Whiting, Near-field radiative heat transfer between macroscopic planar surfaces, Phys. Rev. Lett. 107, 014301 (2011).
  * [10] T. Kralik, P. Hanzelka, M. Zobac, V. Musilova, T. Fort, and M. Horak, Strong near-field enhancement of radiative heat transfer between metallic surfaces, Phys. Rev. Lett. 109, 224302 (2012).
  * [11] I. Latella, S.-A. Biehs, and P. Ben-Abdallah, Smart thermal management with near-field thermal radiation, Opt. Express 29, 16, 24816-24833 (2021).
  * [12] K. Chen, P. Santhanam, S. Sandhu, L. Zhu, and S. Fan, Heat-flux control and solid-state cooling by regulating chemical potential of photons in near-field electromagnetic heat transfer, Phys. Rev. B 91, 134301 (2015).
  * [13] L. Zhu, A. Fiorino, D. Thompson, R. Mittapally, E. Meyhofer, and P. Reddy , Near-field photonic cooling through control of the chemical potential of photons, Nature 566, 239 (2019).
  * [14] W. Srituravanich, N. Fang, C. Sun, Q. Luo, and X. Zhang, Plasmonic Nanolithography, Nano Lett. 4, 1085 (2004).
  * [15] P. Ben-Abdallah, Multitip Near-Field Scanning Thermal Microscopy, Phys. Rev. Lett. 123, 264301 (2019).
  * [16] Y. De 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).
  * [17] A. C. Jones, and M. B. Raschke, Thermal Infrared Near-Field Spectroscopy, Nano Lett. 12, 1475 (2012).
  * [18] R. S. DiMatteo, P. Greiff, S. L. Finberg, K. A. Young-Waithe, H. K. H. Choy, M. M. Masaki, and C. G. Fonstad, Enhanced photogeneration of carriers in a semiconductor via coupling across a nonisothermal nanoscale vacuum gap, Appl. Phys. Lett. 79, 1894 (2001).
  * [19] A. Narayanaswamy and G. Chen, Surface modes for near field thermophotovoltaics, Appl. Phys. Lett. 82, 3544 (2003).
  * [20] M. Laroche, R. Carminati, and J.-J. Greffet, Near-field thermophotovoltaic energy conversion, Appl. Phys. 100, 063704 (2006).
  * [21] K. Park, S. Basu, W. P. King, and Z. M. Zhang, Performance analysis of near-field thermophotovoltaic devices considering absorption distribution, J. Quant. Spectros. Radiat. Transfer 109, 305 (2008).
  * [22] I. Latella and P. Ben-Abdallah, Graphene-based autonomous pyroelectric system for near-field energy conversion, Sci. Rep. 11, 19489 (2021).
  * [23] K. Kloppstech, N. Könne, S.-A. Biehs, A. W. Rodriguez, L. Worbes, D. Hellmann, A. Kittel, Giant heat transfer in the crossover regime between conduction and radiation, Nat. Commun. 8, 14475 (2017).
  * [24] L. Cui, J. Womho, V. Fernández-Hurtado, J. Feist, F. J. García-Vidal, J. C. Cuevas, E. Meyhofer, P. Reddy, Study of radiative heat transfer in Ångström- and nanometre-sized gaps, Nat. Commun. 8, 14479 (2017).
  * [25] R. Messina, S.-A. Biehs, T. Ziehm, A. Kittel, and P. Ben-Abdallah, Heat transfer between two metals through subnanometric vacuum gaps, arXiv:1810.02628 (2018).
  * [26] T. Tokunaga, A. Jarzembski, T. Shiga, K. Park, and M. Francoeur, Extreme near-field heat transfer between gold surfaces, Phys. Rev. B 104, 125404 (2021).
  * [27] T. Tokunaga, M. Arai, K. Kobayashi, W. Hayami, S. Suehara, T. Shiga, K. Park, and M. Francoeur, First-principles calculations of phonon transport across a vacuum gap, Phys. Rev. B 105, 045410 (2022).
  * [28] Y. Guo, C. Adessi, M.Cobian, and S. Merabia, Atomistic simulation of phonon heat transport across metallic vacuum nanogaps, Phys. Rev. B 106 085403 (2022).
  * [29] J. B. Xu, K. Läuger, R. Möller, K. Dransfeld, and I. H. Wilson, Energy-exchange processes by tunneling electrons, App. Phys. A 7, 155 (1994).
  * [30] C. Sönnichsen, Plasmons in metal nanostructures, PhD thesis (2001).
  * [31] P. B. Johnson and R. W. Christy, Optical constants of the noble metals, Phys. Rev. B 6, 4370 (1972).
  * [32] N. W. Ashcroft and N. D. Mermin, Solid state physics, (Harcourt, Orlando, 1976).
  * [33] Y. Ezzahri and K. Joulain, Vacuum-induced phonon transfer between two solid dielectric materials: Illustrating the case of Casimir force coupling, Phys. Rev. B 90, 115433 (2014).
  * [34] J. B. Pendry, K. Sasihithlu, and R. V. Craster, Phonon-assisted heat transfer between vacuum-separated surfaces, Phys. Rev. B 94, 075414 (2016).
  * [35] K. Sasihithlu, J. B. Pendry, and R. V. Craster, Van der Waals Force Assisted Heat Transfer, Z. Naturforsch. 72, 181 (2017).
  * [36] A. I. Volokitin, Effect of Resonant Photon Emission in Radiative Heat Transfer and Generation, JETP Lett. 110, 397–404 (2019).
  * [37] A. I. Volokitin, Effect of an Electric Field in the Heat Transfer between Metals in the Extreme Near Field, J. Phys.: Condens. Matter 32 215001 (2020).
  * [38] B. N. J. Persson, Theory of rubber friction and contact mechanics, J. Chem. Phys. 115 3840 (2001).
  * [39] P. Pinchuk and K. Jiang, Size-dependent Hamaker Constants for Silver and Gold Nanoparticles, Proc. of SPIE 9549, 95491 (2015).
  * [40] G. W Ford and W. H. Weber, Electromagnetic interactions of molecules with metal surfaces, Phys. Rep. 113, 195 (1984).
  * [41] A. Kittel, W. Müller-Hirsch, J. Parisi, S.-A. Biehs, D. Reddig, and M. Holthaus, Near-Field Heat Transfer in a Scanning Thermal Microscope, Phys. Rev. Lett. 95, 224301 (2005).
  * [42] P.-O. Chapuis, S. Volz, C. Henkel, K. Joulain, and J.-J. Greffet, Effects of spatial dispersion in near-field radiative heat transfer between two parallel metallic surfaces, Phys. Rev. B 77, 035431 (2008).
  * [43] J. G. Simmons, Generalized Formula for the Electric Tunnel Effect between Similar Electrodes Separated by a Thin Insulating Film, J. Appl. Phys 34, 1793 (1963).
  * [44] M. V. Berry and K. E. Mount, Semiclassical approximations in wave mechanics, Rep. Prog. Phys. 35, 315 (1972).
  * [45] M. Andrews, Singular potentials in one dimension, Am. J. Phys. 44, 1064 (1976).
  * [46] J. Bardeen, Theory of the Work Function. II. The Surface Double Layer, Phys. Rev. 49, 653 (1936).
  * [47] A.G. Eguiluz, M. Heinrichsmeier, A. Fleszar, and W. Hanke, First-Principles Evaluation of the Surface Barrier for a Kohn-Sham Electron at a Metal Surface, Phys. Rev. Lett. 68, 1359 (1992).
  * [48] Y. Ando and I. Tomohiro, Calculation of transmission tunneling current across arbitrary potential barriers, J. App. Phys. 61, 1497 (1987).
  * [49] Y. Hishinuma, T. H. Geballe, B. Y. Moyzhes, and T. W. Kenny, Refrigeration by combined tunneling and thermionic emission in vacuum: Use of nanometer scale design, App. Phys. Lett. 78, 2572 (2001).
  * [50] A.M. Gabovich, L. G. Il’chenko, E. A. Pashitskii, Yu. A. Romanov, Electrostatic energy and screened charge interaction near the surface of metals with different Fermi surface shape, Surf. Sci. 13, 179 (1980).
  * [51] L.G. Il’chenko and T.V. Goraychuk, Role of the image forces potential in the formation of the potential barrier between closely spaced metals, Surf. Sci. 478, 169-179 (2001).
  * [52] U. Sivan and Y. Imry, Multichannel Landauer formula for thermoelectric transport with application to thermopower near the mobility edge, Phys. Rev. B 33, 551 (1986).
  * [53] R. Lake and S. Datta, Energy balance and heat exchange in mesoscopic systems, Phys. Rev. B 46, 4757 (1992).
  * [54] R. Sánchez and M. Büttiker, Optimal energy quanta to current conversion, Phys. Rev. B 83, 085428 (2011).