%\documentclass[aps,prb,twocolumn,showpacs]{revtex4}
% \documentclass[aps,prb,showpacs]{revtex4}
\documentclass[aps,prb,preprint,showpacs]{revtex4}

\usepackage[cp866]{inputenc}
\textwidth=180mm \textheight=250mm \oddsidemargin=-3mm
\topmargin=-20mm
%\renewcommand{\baselinestretch}{1.1}
%\renewcommand{\baselinestretch}{1.3}
\usepackage{graphicx}
\usepackage{epsfig}
\usepackage{xcolor}
\usepackage{pdfpages}

\usepackage{subfigure}
\usepackage{amssymb,amsmath}


% \input{russian.sty}\relax
%\def\refname{Áèáëèîãðàôèÿ}
%{\scriptsize \color{black} \hfil \textsl{CONFIDENTIAL}}

\pagestyle{myheadings}

\begin{document}

% \markboth{}{CONFIDENTIAL}

%\def\norma#1{\vert \vert {\bf #1} \vert \vert _2}
%\def\cond#1{cond({\bf #1})}
%\def\tcond#1{\widetilde {cond}({\bf #1})}
%\def\dprime#1{ #1^{\prime\prime}}
%\def\abs#1{\vert #1 \vert}
%\def\rerror#1{\norma{\delta #1}}
%\def\ferror#1{\frac{\norma{\delta #1}}{\norma{#1}}}
%\def\d

\def\e{\begin{equation}}
\def\f{\end{equation}}
\def\l#1{\label{#1}}
\def\r#1{(\ref{#1})}
\def\=#1{\overline{\overline #1}}
\def\##1{{\bf#1\mit}}
\def\_#1{{\bf #1}}
\def\o{\omega}
\def\p{\partial}
\def\.{\cdot}
\def\x{\times}
\def\E{\epsilon}
\def\H{\Upsilon}
\def\A{\alpha}
\def\B{\beta}
\def\O{\Xi}
\def\M{\mu}
\def\D{\nabla}
\newcommand{\ds}{\displaystyle}
%\def\l#1{\label{eq:#1}}
%\def\r#1{(\ref{eq:#1})}


\title{Supplementary Material.

Nanoscale nonlocal thermal transport and thermal field emission in
high-current resonant tunnel structures
}

 \author{Michael V. Davidovich$^1$, Igor S. Nefedov$^{1,2}$,  Olga E. Glukhova$^{1,3}$, J. Miguel Rubi$^4$$^{*}$}
 \affiliation{$^1$Department of Physics, Saratov State University, Astrakhanskaya street 83, 410012 Saratov, Russia}
 \affiliation{$^2$ RUDN University, 6 Miklukho-Maklaya St, Moscow 117198, Russian Federation}
\affiliation{$^3$ Laboratory of Biomedical Nanotechnology, I.M. Sechenov First Moscow State Medical University, Bolshaya
Pirogovskaya Street 2-4, 119991 Moscow, Russia } 
\affiliation{$^4$ Departament de Fisica de la Materia Condensada,
University of Barcelona, Marti i Franqu\`es 1, 08028, Barcelona, Spain}
 

 \date{\today}

 
\maketitle

 

The potential in the Schr\"{o}dinger equation is calculated by the formula \cite{26,27,28}:
\begin{equation} 
    \label{S1} 
V\left(x\right)\approx E_{Fc} +W_{c} \frac{\left(1-\alpha /d\right)\left(1+\delta /d\right)^{2} }{\left(1-\delta /d\right)^{2} \varepsilon }\times 
\left[1-\frac{\delta d}{\left(x+\delta \left(1-x/d\right)\right)\left(d-x+x\delta /d\right)} \right]-\frac{eU_{a} x}{d} .
 \end{equation}
Here the work functions (WF) of the cathode and anode are assumed to be the same: $W_c=W_a$, and are related to the parameter $\delta$:
$W_c=c^2/(16\pi\varepsilon_0\delta)$, $\alpha=\delta\left(2\ln{(2)}+1\right).$
In \eqref{S1}, at the cathode $V\left(0\right)=E_{F_c} $, and at the anode $V\left(d\right)=E_{Fc} -eU_{a} $. In the case of different WFs, it is necessary to add the term $\left(E_{Fa} -E_{Fc} \right)x/d$. 
%to \eqref{S1}, 

In addition to Joule heat generation at the cathode and anode, heat is also produced along the mean free path of phonons and electrons, resulting in increased heating of the surface layers of the electrodes.
To simplify our analysis, we assume that in the layer near the cathode of size $\lambda'_{c} =\left(\lambda _{cp} +\lambda _{ce} \right)/2$ the temperature is distributed linearly from $T' _{c} $ at its boundary to $T_{c} $ on the surface of the cathode, and similarly on the anode. The temperatures decrease linearly from $T'_{c} $ to $T_{0} $ at the cathode and from $T'_{a} $ to $T_{0} $ at the anode. 
In the regions from the electrodes to 
$l_{c,a}$, ballistic heat transfer occurs towards the electrodes, followed by diffusion.

Consider diffusive and ballistic heat transfer. In the case of a 3D conductor, the number of phonon modes is proportional to the cross-sectional area $S$: 
$M_{p} \left(\omega \right)=SM_{3D} \left(\omega \right),$ 
where $M_{3D}=D_p(\omega)v_g(\omega\pi/2=k^2(\omega)/(4\pi)$ is the number of phonon modes per unit volume in a 3D material. 

The heat flux depends on the number of phonon states $D_{p} \left(\omega \right)=3\left(\hbar \omega \right)^{2} /\left[2\pi \left(\hbar v_{ D} \right)^{3} \right]$, the phonon transmission coefficient $\tilde{T}_{p} \left(\omega ,l\right)=\lambda _{p} \left(\omega \right)/\left(\lambda _{p} \left(\omega \right)+l\right)$, and the group velocity $v_{g} \left(\omega \right)$, and is given by \cite{1,2,3,4}:
%%%%%%%%%%%%%
\begin{equation} 
    \label{G8} 
Q_{p} \left(T_{1} ,T_{2} \right)=\frac{1}{4} \int _{0}^{\infty }\hbar \omega \tilde{T }_{p} \left(\omega ,l\right)D_{p} \left(\omega \right) v_{g} \left(\omega \right)\left[n\left(\omega ,T_ {1} \right)-n\left(\omega ,T_{2} \right)\right]d\omega. 
\end{equation}
 We will use Debye's dispersion law with a constant average group velocity $v_{g} \left(\omega \right)=v_{D} $, and take the equal velocities averaged over frequency for three polarizations. The equilibrium distribution of phonons is the Bose-Einstein function $n\left(\omega, T\right)=\left(\exp \left(y\right)-1\right)^{-1} $, where $y=\hbar \omega /\left(k_{B} T\right)$. The Debye energy is $\hbar \omega _{D} =k_{B} T_{D} =\hbar v_{D} \left(6\pi ^{2} N\right)^{1/2} $, with $N$ the number of phonons per unit volume.
The heat flow in a conductor of  length \textit{l} and  cross section area \textit{S} with different temperatures at the ends is
\begin{equation} \label{G9}
Q_{p} \left(T_{1} ,T_{2} \right)=\frac{1}{2\pi \hbar } \int _{0}^{\infty }\hbar \omega \frac{ M_{p} \lambda _{p} }{\lambda _{p} +l} \left[n\left(\omega ,T_{1} \right)-n\left(\omega ,T_{2} \right)\right]d\left(\hbar \omega \right).
\end{equation}
In this expression, the number of phonon modes is proportional to the cross-section area: $M_{p} \left(\omega \right)=Sk^{2} \left(\omega \right)/\left(4\pi \right)$, with $k \left(\omega \right)$ the phonon dispersion. 

Considering the Debye model of dispersion with equal average velocities of longitudinal and transverse acoustic modes $v_{D} $, we obtain $M_{p} \left(\omega \right)=3S\left(\hbar \omega \right)^{2} /\left(4\pi \hbar ^{2} v_{D}^{2} \right)$. 
At the second end of the conductor, we have 
$$n\left(\omega ,T_{2} \right)\approx n\left(\omega ,T_{1} \right)+\Delta T\partial n\left(\omega ,T\right)/\partial T.$$
Introducing the phonon conductivity window $w_{p} \left(\omega \right)=-\left(3/\pi ^{2} \right)\omega \partial _{\omega } n\left(\omega \right )/\left(k_{B} T\right)$, considering a small temperature difference $\Delta T=T_{2} -T_{1}$, $T\left(x\right)=T_{1} +x \partial T/\partial x$, $\partial T/\partial x=\Delta T/l$, i.e. $T\left(0\right)=T_{1} $, $T\left(l\right)=T_{2} $, and linearizing \eqref{G9} as $Q_{p} \left(T_{ 1} ,T_{2} \right)=-S\kappa _{p} \Delta T/l$, we obtain the expression for the thermal conductivity of the phonons:
%%%%%%%%%%%%
\begin{equation} \label{G10}
G_{p}\left(T,l\right)=\frac{3k_{B}^{4} T^{3} }{8\pi ^{2} \hbar ^{3} v_{D}^ {2} } \int _{0}^{A } \frac{\lambda _{p} }{\lambda _{p} + l} \frac{y^{3} \exp \left(y\right)}{\left(\exp \left(y\right)-1\right)^{2} } dy,
\end{equation}
where $A=\frac{\hbar \omega _{D}}{k_{B} T}$,
and the thermal conductivity coefficient $\kappa _{p} \left(T\right)$ is given by
\begin{equation} \label{G11}
\kappa _{p} \left(T\right)=\frac{3k_{B}^{4} T^{3} }{8\pi ^{2} \hbar ^{3} v_{D}^ {2} } \int _{0}^{A}\lambda _{p} \frac{y\exp \left(y\right )}{\left(\exp \left(y\right)-1\right)^{2} } y^{2} dy.
\end{equation}
%%%%%%%%%%%%%%
Here the fraction in the integrand is proportional to the conductivity window $w_{p} \left(y\right)=\left(3k_{B} T/\pi ^{2} \right)y\exp \left(y\right)/\left( \exp \left(y\right)-1\right)^{2} $, we restricted the upper limit to the Debye energy, and the transmission coefficient used to obtain  \eqref{G11} is $\tilde{T}_{p} \left(\omega ,l\right)=\lambda _{p} \left(\omega \right)/l$, which is typical for massive structures $l>>\lambda _{p} $. We then obtain: $Q_{p} \left(T_{1} ,T_{2} \right)/S=-\Delta TG_{p} \left(T\right)=-\kappa _{p} \left( T\right)\Delta T/l$. For ballistic thermal conductivity,  we have $l<<\lambda _{p} $, and the integrals \eqref{G8}--\eqref{G10} with infinite limits yield, $G_{p}^{ball} \left(T\right)=\pi ^{2} k_{B}^{4} T^{3} /\left(10\hbar^{3} v_ {D}^{2} \right)$. The ballistic heat flux is thus
\begin{equation} \label{G12}
\frac{Q_{p}^{ball} \left(T_{1} ,T_{2} \right)}{S} =\frac{\pi ^{2} k_{B}^{4} }{ 40\hbar ^{3} v_{D}^{2} } \left(T_{1}^{4} -T_{2}^{4} \right).
\end{equation}
%%%%%%%%%%
The ballistic thermal conductivity coefficient can be determined by taking a small temperature difference $\Delta T=T_{2} -T_{1} <<T_{1} $ on the mean free path (MFP) $\lambda _{p} $, so that $ T_{1}^{4} -T_{2}^{4} \approx -4T_{1}^{3} \Delta T$. We thus obtain from \eqref{G12}: 
$$\kappa _{p}^{ball} =\pi ^{2} \lambda _{p} k_{B}^{4} T^{3} /\left(10\hbar ^{3} v_{D}^{2} \right).$$
On the other hand, assuming the MFP in \eqref{G11} to be constant and the upper limit to be infinite, we obtain $\kappa _{p} =\kappa_{p}^{ball} /4$.

Realistically, however, the Debye energy should constrain the upper limit, as higher frequency phonons do not contribute to heat transport.
To accurately determine the MFP $\lambda _{p}$, we will estimate the thermal conductivity coefficient \eqref{G11}.
 Setting $y=2x$, the integral is expressed in terms of the Bernoulli numbers $B_{2k}$ as \cite{5}:
\begin{equation} \label{G13}
I=4\int _{0}^{y_{D} /2}\frac{x^{3} }{\sinh ^{2} \left(x\right)} dx= -\frac{y_{ D}^{3} \coth \left(y_{D} /2\right)}{2} +\frac{3}{4} \sum _{k=0}^{\infty }\frac{B_ {2k} y_{D}^{p+2k-1} }{\left(p+2k-1\right)\left(2k\right)!} .
\end{equation}

The value of the integral \eqref{G13} for an infinite upper limit is $7.147632$. This value is correct for $y_{D} /3<\pi $. For beryllium $y_{D} =4.72$, it is necessary to use a very large number of terms in \eqref{G13}.
In the region of ultralow temperatures $y_{D} $ is large, i.e. one can take an infinite limit. 
In the region $T\sim T_{D} $, we have $y_{D} \approx 1$, and 
the integral can be calculated with high accuracy using a series expansion, giving 0.4719, and $\kappa _{p} =0.06k_{B}^{ 4} T_{D}^{3} \lambda _{p} /\left(\hbar^{3} v_{D}^{2} \right)$, i.e. this coefficient is approximately four times smaller than the one obtained when the limit is infinite. 
At room temperature, one has $I\approx$ 5.35321. 

The thermal conductivity of metals is the sum of electronic and lattice contributions but is largely determined by the electronic thermal conductivity $\kappa _{e} \left(T\right)=\left(\pi ^{2} /3\right)\sigma _{ e} Tk_{B}^{2} /e^{2} $, where $\sigma _{e} $ is the electronic conductivity. However, for beryllium, phonon (lattice) thermal conductivity also makes a significant contribution due to the very high Debye temperature and low electron dispersion potential. Comparing the conductivity of beryllium and copper, for which the electron MFP  equals 42 nm, we obtain for beryllium the electron MFP $\lambda _{eb} =9.8$ nm. From the Wiedemann-Franz law and the Lorentz number $L=2.47\cdot 10^{-8} $ W$\mathrm{\cdot}$Om/K${}^{2}$ for 300 K, we obtain the electronic thermal conductivity of beryllium $\kappa _{e} =168.4$ W/(m$\mathrm{\cdot}$K). Since $\kappa _{b} =\kappa _{e} +\kappa _{p} =201$ W/(m$\mathrm{\cdot}$K), the phonon thermal conductivity is $\kappa _{p } =32.6$ W/(m$\mathrm{\cdot}$K), from which we obtain the phonon MFP: $\lambda _{p} =27.3$ nm.

To determine the thermal balance, it is convenient to identify areas with dimensions $\lambda _{c} $ and $\lambda _{a} $ near the boundaries of the cathode and anode, in which heat is released or absorbed due to the Nottingham and Peltier-type effects, as well as to the absorption of ballistic phonons transferred through the RTN. 
For beryllium electrodes $\lambda _{c} =\lambda _{a} \sim \left(\lambda _{e} +\lambda _{p} \right)/2\sim 20$ nm. Since these areas are small, we set temperatures $T'_{\left(c,a\right)} $ at their boundaries with the main parts of the electrodes. 

We consider heat balances in four regions: two with dimensions $\lambda _{c} $ and $\lambda _{a} $ near the electrode surfaces, and two extensive regions bordering the thermostat. 
At the electrode boundaries, we set the temperatures $T_{c} $ and $T_{a} $, supposing $T_{c} >T_{a} $. 
This sets the direction of the ballistic flow from the cathode to the anode, which is determined by formula \eqref{G12}. Let us denote this flow per unit surface as $\tilde{Q}_{pca}^{ball} $. 
In the following,  tildes indicate flow densities. For a structure with a vacuum gap, $\tilde{Q}_{pca}^{ball} =0$, and only the heat flux due to photon transfer should be taken into account, which is usually small for low temperatures \cite{27}, and it will be further neglected.

In an inhomogeneous structure, phonon reflection occurs. Therefore, the materials of the grids and dielectric layers should have similar parameters. 
Most consistent with this are beryllium grids placed on diamond dielectric layers.
Since the Kapitza thermal resistances are formed at the boundaries of dielectric and conductive structures, we introduce the transmission coefficient of such a contact \cite{9} $\tilde{\tau }_{0} =4Z_{b} Z_{d}/\left(Z_{b} +Z_ {d} \right)^{2} $, 
where $Z_{(b,d)}=\rho_{(b,d)}v_{(b,d)}$ are acoustic impedances of beryllium and diamond.
There is no electron heat transfer through non-conductive diamond, so there is none for beryllium films either. 

Let us examine ballistic heat transfer in the RTN involving beryllium and diamond. Both materials share the same speed of sound: $v_{b} =v_{d} =12600$ m/s. The density of diamond is 3550 kg/m${}^{3}$, and that of beryllium is 1848 kg/m${}^{3}$, so their acoustic impedances differ only by 1.8 times, and $\tilde{\tau}_{0} =0.916$. 
Given that the density of CVD diamond is slightly lower, we consider  $\tilde{\tau}_{0} =0.93$. 
Let the cathode and anode be made of beryllium with a low work function of 3.92 eV and the highest Fermi energy of 14.6 eV among metals. 
It exhibits high conductivity and melting point, along with a high heat capacity of 1.8 kJ/(kg K). Its Debye temperature is $T_{D} = 1440$ K, and the corresponding Debye energy is 0.124 eV, which corresponds to a frequency of $\omega _{D} =1.87\cdot 10^{14} $ Hz.
Therefore, beryllium is a very good material for emission electronics. 
For diamond, the heat capacity is 502 J/(kg$\cdot$K), the thermal conductivity ranges from 900 to 2300 W/(m$\cdot$K), and the phonon diffusion factor is $\lambda_{pd} =494$ nm \cite{31}. 
For a beryllium layer of length $l_{b} $ the phonon transmission coefficient is $\tilde{\tau}_{b} \left(\lambda _{pb} ,l_{b} \right)=\lambda _{pb} /\left (\lambda _{pb} +l_{b} \right)$.
The same applies to a diamond layer of thickness $l_{d} $, but due to a large MFP, we can take $\tilde{\tau}_{d} =1$.
For an RTN with one quantum well, the total transmission coefficient is
$$\tilde{\tau}_{ca} =\tilde{\tau}_{0}^{4} \tilde{\tau}_{d} \left(\lambda _{pd } ,l_{d} \right)\tilde{\tau}_{b} \left(\lambda _{pb} ,l_{b} \right)\tilde{\tau}_{d} \left(\lambda _ {pd} ,l_{d} \right).$$ 
and with two quantum wells 
$$\tilde{\tau}_{ca} =\tilde{\tau}_{0}^{6} \tilde{\tau}_{a}^{3} \left(\lambda _{ pa} ,l_{a} \right)\tilde{\tau}_{b}^{2} \left(\lambda _{pb} ,l_{b} \right).$$
 For $l_{a} =l_{b} =2$ nm, $\tilde{\tau}_{a} =0.996$, and $\tilde{\tau}_{b} =0.9459$, we obtain the transmission coefficient for one well $\tilde{\tau}_{ca} =0.702$, and for two wells $\tilde{\tau}_{ca} =0.572$.
The flux density from the cathode to the anode is $\tilde{Q}_{p}^{ball} =\tilde{Q}_{pcc}^{ball} \tilde{\tau}_{ca} $, and from anode to cathode $-\tilde{Q}_{p}^{ball}$.
An amount of heat  $\tilde{Q}_{c} =\kappa _{c} \left(T_{c} -T'_{c} -T'_{0} \right)$ flows from the cathode to the thermostat. 
The diffusive thermal conductivity determines this flux, so $\kappa_{c} $ is the usual experimental thermal conductivity coefficient.
Similarly, the heat from the anode to the thermostat is $\tilde{Q}_{a} =\kappa _{a} \left(T_{a} -T'_{a} -T_{0} \right)$. We calculate the current density $J=J^{+} -J^{-} $, therefore the density of Joule heat sources at the cathode is uniform and equal to $q_{c} =\rho _{c} J^{2} $, and at the anode $q_{a} =\rho _{a} J^{2}$.
%where resistivities are indicated.

The release and absorption of heat with source densities $q_{c} $ and $q_{a} $ due to the Nottingham or 
the Peltier-type effects
occurs at electron mean free paths $\lambda _{ec} $ and $\lambda _{ea}$.
For the densities of heat sources per unit surface in areas near the electrode boundaries due to the Nottingham and Peltier effects, we have $q_{c,a} = q'_{c,a}+ q''_{c,a}$
Here
\begin{equation}\label{G18}
q'_{c,a} =\frac{m_{e} k_{B} T_{c,a} }{2\pi ^{2} \hbar ^{3} } \int _{0}^{\infty }\left(\mu_{c,a}-E\right)D^{\pm} \left(E,U_{a} \right)\ln \left(1+\exp \left(\frac{\mu_{c,a} -E} {k_{B}T_{c,a}} \right)\right)dE ,
\end{equation}
Quantities marked with a double prime correspond to electrons hitting the electrode. The corresponding densities are given by
\begin{equation} \label{G20}
q''_{c,a} =\frac{m_{e} k_{B} T_{a,c} }{2\pi ^{2} \hbar ^{3} } \int _{0}^{\infty } \left(E-\mu_{c,a} \right)D^{\mp} \left(E,U_{a} \right)\ln \left(1+\exp \left(\frac{\mu_{a,c} -E}{k_{B} T_{a,c} } \right)\right)dE.
\end{equation}

From equation \eqref{G20}, we obtain 
$$T'_{c,a} =\rho_{c,a} J^{2}\left(l_{c,a}-\lambda'_{c,a} \right )^{2} \lambda'_{c,a} /\left(l_{c,a} \kappa_{c,a} \right)+T_{c,a} \left(1-\lambda'_{c,a}/l_{c,a} \right)+T_{0}\lambda'_{c,a}/l_{c,a}. $$
Then one has $\Delta Q_{c,a} =0$, and we write the remaining conditions in the form $\Delta Q=\left(\Delta Q'_{c} l_{c} /\kappa _{c} \right)^{2} +\left(\Delta Q'_{a} l_{a} /\kappa _{a} \right)^{2} =0$. 
To find the temperatures $T_{c} $ and $T_{a} $, we minimize $\Delta Q(T_c,T_a)$. One then obtains
\begin{equation} \label{G22}
T_{c,a} =f_{c,a} \left(T_{c} ,T_{a} \right)=T_{0} +\frac{q_{c,a}l_{c,a} +\rho_{c,a}J^{2}\left[\lambda'_{c,a}l_{c,a} +\left(l_{c,a} -\lambda'_{c,a}\right)^{2} \right]-\sigma _{p} \left(T_{c}^{4} -T_{a}^{4} \right)l_{c} }{\kappa_{c,a}} ,
\end{equation}
\begin{equation} \label{G24}
\Delta Q'_{c,a} =q_{c,a} +\rho_{c,a} J^{2}\frac{\lambda'_{c,a}l_{c,a} +\left(l_{c,a} -\lambda'_{c,a}\right)^{2}}{l_{c,a} } +\kappa_{c,a} \frac{\left(T_{0}-T_{c,a} \right)}{l_{c,a}}-\sigma _{p} \left(T_{c}^{4} -T_{a}^{4} \right).
\end{equation}
%%%%%%%%%%%%%%%%
Here $\sigma _{p} =\pi ^{2} k_{B}^{4} \tilde{\tau}_{ca} /\left(40\hbar ^{3} v_{p}^{2 } \right)$ is analogous to the Stefan-Boltzmann constant (only for the average group sound velocity). 
For low currents, the sources are small and consequently, the temperatures remain low. Temperatures can be further reduced by lowering $\tau_{0} $
or by decreasing the lengths  $l_{c}$ and $l_{a} $.
Increasing the width of barriers $t_{n} $, $n=1,2,3 $,  also reduces the current, while it increases with the increase of the wells width $t_{g} $. 
In the latter case, additional resonant levels may arise. For low and approximately equal temperatures, the ballistic transport (last term on the right-hand side) can be neglected, and then we obtain \cite{28}
\begin{equation} \label{G26}
T_{c,a} =T_{0} +q_{c,a}l_{c,a} /\kappa_{c,a} +\rho_{c,a} J^{2} \left[\lambda'_{c,a} l_{c,a } +\left(l_{c,a} -\lambda'_{c,a} \right)^{2} \right]/\kappa_{c,a}.
\end{equation}
These equations apply at low temperatures when the temperature difference is small or there is a vacuum gap. 
In the latter case, the radiative heat transfer is approximately $\left(v_{p} /c\right)^{2} \sim 4\cdot 10^{-5} $ times less than heat transfer by phonons \cite{27}.
For vacuum gaps or low temperatures, these equations become accurate, and in the latter case, the right-hand sides do not depend on temperatures. 
For RT, the current density \textit{J} is small when the anode voltage \textit{U${}_{a}$} is low, while $D^{+} \approx D^{-}$, and at room temperature  $q'_{c,a} \approx -q''_{c,a} $, i.e. $q_{c} \approx q_{a} \approx 0$. 
Then the temperatures are determined only by the Joule heat and its outflow into the thermostat: 
$$T_{c,a} =T_{0} +\rho _{c,a} J^{2} \left[\lambda'_{c,a} l_{c,a} +\left(l_{c,a} -\lambda'_{c,a} \right)^{2} \right]/\kappa _{c,a}.$$ 

If temperatures differ greatly, and for example $T_{a} >T_{c} $, the approximate temperature at which the ballistic transport begins to dominate is $T_{a} >\sqrt[{4}]{\rho _{ a} l_{a} J^{2} } $. For $J\sim 3\cdot 10^{13} $A/m${}^{2}$, $\rho _{a} =3.5\cdot 10^{-8} $ Ohm$\mathrm{ \cdot}$m, and $l_{a} =3\cdot 10^{-7}$m, we obtain $T_{a} >800$ K.

From equation \eqref{G20} we obtain 
$$T'_{c,a} =\rho_{c,a} J^{2}\left(l_{c,a}-\lambda'_{c,a} \right )^{2} \lambda'_{c,a} /\left(l_{c,a} \kappa_{c,a} \right)+T_{c,a} \left(1-\lambda'_{c,a}/l_{c,a} \right)+T_{0}\lambda'_{c,a}/l_{c,a}. $$
Then one has $\Delta Q_{c,a} =0$, and we write the remaining conditions in the form $\Delta Q=\left(\Delta Q'_{c} l_{c} /\kappa _{c} \right)^{2} +\left(\Delta Q'_{a} l_{a} /\kappa _{a} \right)^{2} =0$. 
We decompose the integrals into two parts, one up to $3E_{F} $ and the other in the region $\left(3E_{F} ,\infty \right)$.

%%%%%%%%%%%%
To find the temperatures $T_{c} $ and $T_{a} $, we minimize $\Delta Q(T_c,T_a)$. One then obtains
\begin{equation} \label{Gg22}
T_{c,a} =f_{c,a} \left(T_{c} ,T_{a} \right)=T_{0} +\frac{q_{c,a}l_{c,a} +\rho_{c,a}J^{2}\left[\lambda'_{c,a}l_{c,a} +\left(l_{c,a} -\lambda'_{c,a}\right)^{2} \right]-\sigma _{p} \left(T_{c}^{4} -T_{a}^{4} \right)l_{c} }{\kappa_{c,a}} ,
\end{equation}
%%%%%%%%%%
Here $\sigma _{p} =\pi ^{2} k_{B}^{4} \tilde{\tau}_{ca} /\left(40\hbar ^{3} v_{p}^{2 } \right)$ is analogous to the Stefan-Boltzmann constant (only for the average group sound velocity). 
For low currents, the sources are small and consequently, the temperatures remain low. Temperatures can be further reduced by lowering $\tau_{0} $
or by decreasing the lengths  $l_{c}$ and $l_{a} $.
Increasing the width of barriers $t_{n} $, $n=1,2,3 $,  also reduces the current, while it increases with the increase of the wells width $t_{g} $. 
In the latter case, additional resonant levels may arise. For low and approximately equal temperatures, the ballistic transport (last term on the right-hand side) can be neglected, and then we obtain \cite{28}.
These equations apply at low temperatures when the temperature difference is small or there is a vacuum gap. 
In the latter case, the radiative heat transfer is approximately $\left(v_{p} /c\right)^{2} \sim 4\cdot 10^{-5} $ times less than heat transfer by phonons \cite{27}.
For vacuum gaps or low temperatures, these equations become accurate, and in the latter case, the right-hand sides do not depend on temperatures. 
For RT, the current density \textit{J} is small when the anode voltage \textit{U${}_{a}$} is low, while $D^{+} \approx D^{-}$, and at room temperature  $q'_{c,a} \approx -q''_{c,a} $, i.e. $q_{c} \approx q_{a} \approx 0$. 
Then the temperatures are determined only by the Joule heat and its outflow into the thermostat: 
$$T_{c,a} =T_{0} +\rho _{c,a} J^{2} \left[\lambda'_{c,a} l_{c,a} +\left(l_{c,a} -\lambda'_{c,a} \right)^{2} \right]/\kappa _{c,a}.$$ 

If temperatures differ greatly, and for example $T_{a} >T_{c} $, the approximate temperature at which the ballistic transport begins to dominate is $T_{a} >\sqrt[{4}]{\rho _{ a} l_{a} J^{2} } $. For $J\sim 3\cdot 10^{13} $A/m${}^{2}$, $\rho _{a} =3.5\cdot 10^{-8} $ Ohm$\mathrm{ \cdot}$m, and $l_{a} =3\cdot 10^{-7}$m, we obtain $T_{a} >800$ K.


%%%%%%%%%
For beryllium electrodes, we take these lengths to be about 10 nm. To simplify the model, it is convenient to assume that heat release occurs uniformly over certain lengths $\lambda '_{c} $ and $\lambda'_{a} $ of the order of MFP, i.e. determine the energy balance in these areas near the surfaces of the electrodes and in the remaining areas of the electrodes. 
Then the balance in the cathode region is: $\lambda'_{c} \left(q_{c} +q_{0c} \right)=\tilde{Q}_{p}^{ball} +\tilde{Q}_{c} $, where $\tilde{Q'}_{c} =\kappa_{c} \left(T_{c} -T'_{c} \right)$ is the flux to the cathode. In the corresponding area of the anode the balance reads as $\lambda'_{a} \left(q_{a} +q_{0a} \right)=-\tilde{Q}_{p}^{ball} +\tilde{Q' }_{a} $. 
Here $\tilde{Q'}_{a} =\kappa _{a} \left(T_{a} -T'_{a} \right)$ is the flow to the anode. Let us write similar balances on the cathode $\left(l_{c} -\lambda '_{c} \right)q_{c} =\tilde{Q}_{c} $ and anode $\left(l_{a} - \lambda '_{a} \right)q_{a} =\tilde{Q}_{a} $, $\tilde{Q}_{c} =\kappa _{A} \left(T_{c} -T_{0} \right)$, $\tilde{Q}_{a} =\kappa _{a} \left(T_{a} -T_{0} \right)$. We assume that heat generation due to electron scattering is proportional to length. Within the RTN itself, no heat is generated. Additionally, we neglect the temperature dependencies of the coefficients included in the heat balance formulas.

An alternative to diamond can be beryllium oxide (longitudinal wave speed of sound 11300 m/s, density 2800 kg/m$^3$). Alternatively, tungsten electrodes can be used instead of beryllium electrodes. The speed of sound in tungsten is 1700 m/s, and the density is 19300 kg/m$^3$, significantly reducing $\tilde T_{ca}$, but making RTH more resistant to heat.

It is important to note that the results presented provide stationary temperature and current values based on the heat balance. However, this does not imply that the system will necessarily reach these states when the anode voltage is turned on (for a given grid) or when both voltages are turned on simultaneously. These are non-stationary and nonlinear transient processes that should be rigorously analyzed using non-stationary equations.
In particular, processes such as explosive emission may occur during this evolution. The transition to the calculated states at temperatures below the melting point of the electrodes can be achieved through relatively smooth voltage changes, allowing the system to reach an equilibrium state with each adjustment. The value of the material parameters depends on temperature, which can be accounted for in the balance (6,7) (main text). By applying a smooth voltage change with discrete time $U_{a} \left(t_{n} \right)=U_{a} \left(n\Delta t\right)=n\Delta tU_{0} /N$, \textit{n=}1,2,{\dots},\textit{N}, it is possible to formulate non-stationary nonlinear heat conduction equations, with temperature-dependent coefficients and solve them iteratively.



\begin{thebibliography}{52}

\bibitem{26}  K.L. Jensen.~A Thermal-Field-Photoemission, Model and Its Application. In: Modern Developments in Vacuum Electron Sources, TAP \textbf{135}, 345--385 (2020). DOI: 10.1007/978-3-030-47291-7\_8.

\bibitem{27}  M.V. Davidovich, I.S. Nefedov, O.E. Glukhova, M.M. Slepchenkov. Toward the theory of resonant-tunneling triode and tetrode with CNT--graphene grids. J. Appl. Phys. \textbf{130}, 204301-11 (2021). DOI: 10.1063/5.0067763.

\bibitem{28}  M.V. Davidovich, I.S. Nefedov, O.E. Glukhova, M.M. Slepchenkov, J.M. Rubi, Field emission in vacuum resonant tunneling heterostructures with high current densities // Scientific Reports,  \textbf{13}, 19365 (2023). DOI: 10.1038/s41598-023-44900-2.

\bibitem{1}  C. Jeong, S. Datta, M. Lundstrom. Full Dispersion vs. Debye Model Evaluation of Lattice Thermal Conductivity with a Landauer Approach, J. Appl. Phys.   {\bf 109},  073718/8 (2011).

\bibitem{2}  D. E. Angelescu, M. C. Cross, M. L. Roukes. Heat transport in mesoscopic systems, Superlatt. Microstruct. {\bf 23}, 673 (1998).

\bibitem{3}  L. G. C. Rego, G. Kirczenow. Fractional exclusion statistics and the universal quantum of thermal conductance: A unifying approach, Phys. Rev. B {\bf 59}, 13080 (1999).

\bibitem{4}  C. Jeong, S. Datta, M. Lundstrom. Thermal conductivity of bulk and thin-film silicon: a Landauer approach, J. Appl. Phys. {\bf 111}, 093708 (2012).

\bibitem{5}	A.P. Prudnikov, Yu. A. Brychkov and O.I. Marichev. Integrals and Series. Vol. 1, Elementary Functions. New York: Gordon and Breach Science Publishers, 1986.

\bibitem{9} D.G. Cahill, W.K. Ford, K.E. Goodson, G.D. Mahan, A. Majumdar, H.J. Maris, R. Merlin, S.R. Phillpot, Nanoscale thermal transport, J. Appl. Phys. \textbf{93}, 793--818 (2003). DOI: 10.1063/1.1524305.

\bibitem{31}  K.L. Jensen, M.S. McDonald, M.K. Dhillon, D. Finkenstadt, A. Shabaev, M. Osofsky. Thermal-field emission from cones and wires editors-pick. J. Vacuum Sci. Technol. B \textbf{40}, 022801 (2022). DOI: 10.1116/6.0001656.

\end{thebibliography}



\end{document}