Transport in electron−photon systems

Abstract

We review the description and modeling of transport phenomena among the electron systems coupled via scalar or vector photons. It consists of three parts. The first part is about scalar photons, i.e., Coulomb interactions. The second part is with transverse photons described by vector potentials. The third part is on ϕ = 0 or temporal gauge, which is a full theory of the electrodynamics. We use the nonequilibrium Green’s function (NEGF) formalism as a basic tool to study steady-state transport. Although with local equilibrium it is equivalent to the fluctuational electrodynamics (FE), the advantage of NEGF is that it can go beyond FE due to its generality. We have given a few examples in the review, such as transfer of heat between graphene sheets driven by potential bias, emission of light by a double quantum dot, and emission of energy, momentum, and angular momentum from a graphene nanoribbon. All of these calculations are based on a generalization of the Meir−Wingreen formula commonly used in electronic transport in mesoscopic systems, with materials properties represented by photon self-energy, coupled with the Keldysh equation and the solution to the Dyson equation.

Graphical abstract

Keywords

quantum transport / thermal radiation / scalar and vector photons / nonequilibrium Green's function

1 Introduction

Electron and photon interaction plays an important role in our understanding of basic physics in condensed matter, atomic physics, and nano-optics [1,2]. It is also the principle governing photosynthesis and thermal radiation, and underpins many functioning devices such as lasers and solar cells. This review focuses on one aspect of the electron-photon interaction, that is, the transfer of energy and other conserved quantities mediated by photons between metals in local thermal equilibrium and nonequilibrium states.

Radiation of energy happens for any object, due to the electrically charged nature of ions and electrons that make up materials microscopically. The fluctuation of charge and thus current is inevitable, leading to thermal radiation. At the turn of the twentieth century, Planck, in attempting to understand blackbody radiation, discovered the quantum nature of radiation and proposed his famous formula for the intensity of the radiation [3]. He is aware of the limitation of his theory, i.e., the cavity of the blackbody or the length scale of the object should be much larger than the wavelength of the photons. But it took 70 years to see the effect of geometry from so-called near-field radiative heat transfer. Heat transfer increases when the distance between objects decreases. Motivated by the earlier experimental results [4,5], Polder and van Hove worked out a theory [6], now known as fluctuational electrodynamics (FE), which in turn is based on the idea of Rytov of Maxwell’s equations with stochastic random currents [7,8]. Together with the fluctuation dissipation theorem as a corner stone [9], the theory is complete. The theory has been applied to a variety of situations and geometries [10-12] and has been reviewed extensively [13-20]. Recent experiments [21-24] mostly confirm the theory, but there remain some controversies [25-27]. An interesting recent development is the super-Planckian heat transport [28,29], e.g., two sheets side by side instead of face to face, where it is found that the enhancement of transport occurs even for the far field [30].

Another line of research, developed in parallel to the radiation problem, is the Casimir force. Casimir predicted, for ideal metal surfaces with a vacuum gap, there is an attractive force between them, due to vacuum fluctuations of the field [31]. A more complete theory is developed by Lifshitz from fluctuational electrodynamics for nonzero temperatures [32]. The effect is small, and experimental verification was only attained much later in the 1980s and 90s [33-38]. These forces can be utilized for ultra-sensitive measuring devices [39]. The recent development has been in the dynamic Casimir effect [40-42], where the objects can move and couple to the vacuum fluctuations. There is a close relationship between the two problems above. The radiation is due to thermal fluctuations, thus, a nonzero temperature is necessary. But the force can exist even at zero temperature. Both problems share some common features. Their solutions all require knowledge of the response functions of the materials and solutions to Maxwell’s equations. They differ in that energy radiation is related to the difference of the Bose functions $N$ at different temperatures, while the force is related to $N+1/2$, the Bose function with the zero-point motion contribution. It is important that we see these connections and have a unified theory.

Another transport quantity is the angular momentum related to the rotational symmetry of the problem. The radiation carries, in addition to energy and momentum, angular momentum [43-45]. The radiation field patterns are related to the orbital angular momentum [46]. This can be used for transfer of information as extra communication channels [47]. It turns out that emission of angular momentum requires a broken time reversal symmetry. More generally, nonreciprocity has been an emerging area of active research [48,49]. Both the momentum and angular momentum radiation of a single object require nonreciprocity.

In the theoretical and computational approaches for radiative heat and momentum transfers, earlier theories have been developed for the special geometry of parallel plates with solutions to the scattering problems. For arbitrary geometries, Messina et al. [50,51] and Krüger et al. [10,52] developed a formalism based on the $T$ operators and Lippmann-Schwinger equation [53]. The fluctuation of the electromagnetic field can also be viewed as coming from dipoles [54]. A general theory for many objects is developed and recently reviewed [17]. Other computational techniques are developed from the continuum engineering perspectives [55,56].

In this review, we systematically develop a nonequilibrium Green’s function approach with the aim of going beyond fluctuational electrodynamics (FE). The nonequilibrium Green’s function has been the standard approach to study mesoscopic electron transport [57,58], and also, to some extent, phononic transport in the ballistic regime [59,60]. For some reason, the use of it for photon transport is rather scarce [61,62]. The nonequilibrium Green’s function method is a powerful tool to study transport as transport, by its very nature, must deal with nonequilibrium situations. Since Maxwell’s equations are linear, the associated equations for the Green’s functions are also linear. This implies a formally exactly solvable problem, just like in the free electron or ballistic phonon case. The solutions are encapsulated as a pair of equations for two Green’s functions, namely, retarded and lesser [63]. The retarded Green’s function satisfies a Dyson equation, $D^r=$ $v+v{\mathrm{\Pi }}^r{D}^r$, where $v$ is the free Green’s function of photons when the materials are absent. This equation describes the “dynamics”. A related equation, $D^{<}=D^r{\mathrm{\Pi }}^{<}{D}^a$, known as the Keldysh equation, describes the thermal distribution. Finally, the transported physical observables can be expressed in terms of $D$ and materials properties $\mathrm{\Pi }$ known as the Meir−Wingreen formulas. Here $\mathrm{\Pi }$ is the charge-charge or current−current correlation under a random phase approximation. The advantage of this language and methodology is that it is fairly general: i) it is fully quantum-mechanical with quantum electrodynamics; ii) local equilibrium or not, where the fluctuation-dissipation theorem may or may not hold, the theory remains the same; iii) reciprocal or not, the formalism applies to both; iv) the nonequilibrium system is set up explicitly by modeling the electrons connected to fermionic baths; v) change of materials properties Π due to extreme proximity can be handled self-consistently via many-body formalism. These are not in the usual spirit of FE, which requires explicit local equilibrium assumptions and sometimes explicit reciprocity. However, if local equilibrium is assumed, we recover the standard fluctuational electrodynamics results. The global or local equilibrium implies the fluctuation-dissipation theorem [64]. In the language of NEGF, this is, $D^{<}=N(D^r-D^a)$ and ${\mathrm{\Pi }}^{<}=N({\mathrm{\Pi }}^r-{\mathrm{\Pi }}^a)$, in a concise form for the field correlation and materials properties.

Unlike FE formulated in a gauge independent manner in the electric field $\mathit{E}$ and magnetic induction $\mathit{B}$, the NEGF theory is formulated in terms of the gauge dependent scalar and vector potentials. This is necessary for a more fundamental theory as quantization requires a specification or choice of gauge. Our review consists of three parts. Part I is a theory based purely on Coulomb interaction for the electrons, which we can also formulate as a scalar field theory. In Part II, the quantization for the vector field in transverse gauge is presented. This is a standard choice in condensed matter physics and quantum optics. In the last Part III, we develop a full theory in the temporal, or also called axial gauge, which sets the scalar field to zero. The last choice is more economical, and is very closely related to the gauge independent FE. For example, the Green’s function here is the same as the dyadic Green’s function in FE up to a constant. In order not to make the paper too long, we do not attempt to review the basics of NEGF. Reader unfamiliar with the general NEGF method should consult the relevant literature cited [59,63,65,66]. We believe that the NEGF approach is not reinventing the wheel, but rather, opens up new avenues to study transport in more general settings.

Part I Scalar photons

The electromagnetic field, in the limit that the speed of light $c$ goes to infinity, reduces to a rather simple theory of the Poisson equation, $-{\mathrm{\nabla }}^2\varphi =\rho /{ϵ}_0$, relating charge density and scalar potential $\varphi$ instantaneously. We shall call such an approximation to electrodynamics the non-retardation limit. This is not a bad approximation; in fact the whole of electronic structure theories, such as the density functional theory, and the bulk of condensed matter physics, are based on such Coulomb interactions. When the typical distance is much smaller compared to $c\tau$, where $\tau$ is some characteristic time scale, the non-retardation limit is a good approximation. For heat transport mediated by electromagnetic fields, this distance scale is the thermal de Broglie wavelength for the massless photon, $\lambda \approx c\hslash /(k_{B}T)$. Here the time scale is identified with the thermal energy scale. At 300 K it is about a few micrometers. In this first part, we discuss theories for the pure Coulomb problems which is the dominant mechanism of near-field heat transport [67] at sub-micron distances.

2 Heat transfer from capacitor physics

Consider a collection of good conductors, each of which is characterized by a constant potential on the conductor, ${\varphi }_i$, and the charge on the conductor, $q_i$. Given the potentials on the conductors and at infinity, as well as the geometry of the setup, the Poisson equation determines the potential distribution in space everywhere, thus also determines the charges on the conductors by Gauss’s law. For this electrostatic equilibrium and time-independent situation, the charges are related to the potentials linearly by [68,69]

(1)$\begin{array}{l}{q}_i={\displaystyle \sum _j{C}_{ij}{\varphi }_j,}\end{array}$

where $C_{ij}$ is defined as the capacitance of conductor $i$ induced by $j$. Symbolically, we will write the above equation as $q=v^{-1}\varphi$, where $q$ and $\varphi$ denote column vectors and $v^{-1}$ is a matrix formed by the capacitance $C_{ij}$. Formally, the potential produced by the charges is just $\varphi =vq$. The total electrostatic energy is $\frac{1}{2}\sum _{i,j}{C}_{ij}{\varphi }_i{\varphi }_j$. It is clear that the matrix must be symmetric, $C_{ij}=C_{ji}$. In a simple parallel plate capacitor, we have two conductors 1 and 2. The overall charge neutrality requires $q_1+q_2=0$, and the charges should only depend on the potential difference due to gauge invariance (${\varphi }_i\to {\varphi }_i+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}$). If the total charge is not zero, it would mean that the total electrostatic energy is infinite due to a nonzero electric field outside the capacitor, which is unreasonable. These conditions imply that the four coefficients of capacitance are not independent, and should satisfy $C_{11}=C_{22}=-C_{12}=-C_{21}$, given by the usual formula of $C=C_{11}={ϵ}_{0}A/d$. Here ${ϵ}_0$ is the vacuum permittivity, $A$ is the area of the plate, and $d$ is the distance between the plates. This recovers the textbook definition of capacitance by $C=q/\mathrm{\Delta }\varphi$.

Can the capacitor system be considered as a heat transfer system mediated by the Coulomb force? In principle, yes; the charges on the conductors can fluctuate and thus transfer energy. If each conductor is maintained at a different temperature, there will be a net transfer of energy. However, if the conductors are macroscopic in size, the fluctuation will be very small. We expect substantial fluctuations only when the system is at the nanoscale. There is one constraint it must fulfill; this is the total charge conservation. We imagine that each conductor is connected to a battery maintained at a certain temperature and chemical potential. If we focus only on the conductor, the charge on the conductor is not conserved as it can go into or out of the battery from time to time. So, we need an equation to describe this process. We look towards the fluctuational electrodynamics [7], in the non-retardation limit.

On a phenomenological ground, we put down the following Langevin-like stochastic equation for the fluctuating charge capacitor system [70],

(2)$\begin{array}{l}{v}^{-1}\varphi =\delta q+\mathrm{\Pi }\varphi .\end{array}$

Here $v^{-1}$ is the capacitor matrix, $\delta q$ is a vector of the fluctuational charges, $\mathrm{\Pi }\varphi$ gives the induced charge response when the potential is changed by $\varphi$. In the above equation, we use the symbol $\varphi$ to mean the change relative to the static value due to the time-dependent change $\delta q$ in charge. The charge response is retarded, in the sense that $\mathrm{\Pi }\varphi$ depends on the history of the potential,

(3)$\begin{array}{l}\mathrm{\Pi }\varphi \to {\displaystyle {\int }^t\mathrm{\Pi }(t-t^{\prime })\varphi (t^{\prime })\mathrm{d}{t}^{\prime }.}\end{array}$

Equation (2) generalizes the static case, Eq. (1), for the fluctuational charge and field in time. It is best to think of the stochastic equation in the frequency domain, defined by

(4)$\begin{array}{l}\varphi (t)={\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\varphi (\omega ){\mathrm{e}}^{-\mathrm{i}\omega t},}\end{array}$

and similarly for $\delta q(t)$ and $\mathrm{\Pi }(t)$, then the convolution in time becomes multiplication in the frequency domain. Note that we can think of $(D^r{)}^{-1}=v^{-1}-\mathrm{\Pi }$ as a frequency-dependent capacitance matrix. It is this frequency dependence of the material response $\mathrm{\Pi }$ to the field fluctuations that gives rise to energy transfer. The Green’s function $D^r=v+v\mathrm{\Pi }{D}^r$ will be needed to express the energy transport. The solution to Eq. (2) is $\varphi =D^r\delta q$.

No matter what the detailed mechanisms are for the charge fluctuation, in thermal equilibrium, the charge fluctuations must be related to the linear response by the fluctuation-dissipation theorem [9,71],

(5)$\begin{array}{l}{\displaystyle \frac{1}{\mathrm{i}\hslash }⟨\delta q\delta q^{\mathrm{T}}{⟩}_{\omega }=\left[N(\omega )+\frac{1}{2}\right]\left[\mathrm{\Pi }(\omega )-{\mathrm{\Pi }}^{†}(\omega )\right]=\overline{\mathrm{\Pi }}.}\end{array}$

Here, the superscript “$\mathrm{T}$” is matrix transpose, and the dagger “$†$” is for Hermitian conjugate. The left-hand side is the fluctuational charge correlation, $⟨\delta q(t)\delta q(t^{\prime }{)}^{\mathrm{T}}⟩/(\mathrm{i}\hslash )$, which is a function of the time difference in steady state, Fourier transformed into $\omega$ space. $N(\omega )=$ $1/\left[\exp (\beta \hslash \omega )-1\right]$ is the Bose (or Planck) function at temperature $T=1/(k_B\beta )$. In many-body theory, the charge-charge correlation is the susceptibility $\chi =\mathrm{\Pi }+\mathrm{\Pi }v\mathrm{\Pi }+\cdots$ [72], but here it is the polarizability $\mathrm{\Pi }$ in the fluctuational electrostatics. The reason is that $\varphi$ here is the total actual field on the matter. As we will see, this is consistent with the Dyson equation and the Keldysh equation in a more rigorous treatment later. For the moment, we consider $\mathrm{\Pi }$ as a phenomenological function that is given. Subsequently, when we build a microscopic (quantum) model, its functional form can be worked out. Here with this minimum information, how can we describe the transfer of energy among the conductors when they are out of equilibrium?

In order not to deviate too much from thermal equilibrium, we assume local equilibrium. This is possible if the charge-charge correlation can be localized so that different regions are not correlated, and each region has its separate temperature and chemical potential. To be definite, let us partition the conductors as belonging to the left ($L$) or right ($R$) such that $\mathrm{\Pi }$ is block-diagonal,

(6)$\begin{array}{l}\mathrm{\Pi }=(\begin{array}{cc}{\mathrm{\Pi }}_L& 0\\ 0& {\mathrm{\Pi }}_R\end{array}).\end{array}$

The fluctuation-dissipation theorem, Eq. (5), remains the same, except that it is applied to the left sites or right sites separately with $T_L$ or $T_R$, and the charge correlation is zero between the left side and right side.

Next, we consider the average energy transfer from left to right, based on Joule heating, and the local fluctuation-dissipation theorem. On a continuum, $\mathit{j}\cdot \mathit{E}$ is the work density done by the field to the charge. Using the relation between the field and potential $\mathit{E}=-\mathrm{\nabla }\varphi$, the continuity equation $\dot{\rho }+\mathrm{\nabla }\cdot \mathit{j}=0$ for charge conservation, and an integration by parts, we also have the work as $-\dot{\rho }\varphi$ per unit volume [73]. For discrete charge, this is $-{\dot{q}}^{\mathrm{T}}\varphi$. Since Eq. (2) is linear in $\varphi$, we can consider the effect of random noises of two regions separately. Turning off $\delta q_R$, the energy transfer to the right side per unit time due to the fluctuation of charge $\delta q_L$ of the left side is

(7)$\begin{array}{l}{I}_{L\to R}=-⟨{\dot{q}}_R^{\mathrm{T}}{\varphi }_R⟩,\end{array}$

where $q_R={\mathrm{\Pi }}_R{\varphi }_R$ and ${\varphi }_R=D_{RL}^r\delta q_L$. Here we collect all the values of discrete potentials on the right side as a column vector ${\varphi }_R$. The charge $q_R$ in Eq. (7) is the induced one due to the field. These are time domain quantities, for example,

(8)$\begin{array}{l}{\varphi }_R(t)={\displaystyle \int D_{RL}^r(t-t^{\prime })\delta q_L(t^{\prime })\mathrm{d}{t}^{\prime }.}\end{array}$

This is the solution to Eq. (2) restricted to the right side for the potential. We assume that the system is in steady state and $I_{L\to R}$ is in fact independent of time. Representing all the time domain quantities by their Fourier transforms in frequency domain, after some lengthy but straightforward algebra, we find

(9)$\begin{array}{l}{I}_{L\to R}={\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\hslash \omega \mathrm{T}\mathrm{r}(D_{LR}^a{\mathrm{\Pi }}_R^a{D}_{RL}^r{\overline{\mathrm{\Pi }}}_L).}\end{array}$

Here $\mathrm{T}\mathrm{r}$ stands for matrix trace, ${\mathrm{\Pi }}_R^a=({\mathrm{\Pi }}_R{)}^{†}$, and $D_{LR}^a=(D_{RL}^r{)}^{†}$. The last factor is due to the noise correlation by the local fluctuation-dissipation theorem.

The energy pumped from $R$ to $L$ by $\delta q_R$ can be obtained similarly by swapping the index $L\leftrightarrow R$. The overall net heat current from left to right is given by the difference, $I_L=I_{L\to R}-I_{R\to L}$. The expression can be simplified using the fact that (i) $I_{L\to R}$ and $I_{R\to L}$ are real, so that we can take the Hermitian conjugate of the factors inside the matrix trace and add them, then divide by 2; (ii) we can perform cyclic permutation under trace; (iii) both $D^r$ and $\mathrm{\Pi }$ are symmetric matrices, thus, e.g., $D^a=(D^r{)}^{†}=(D^r{)}^{\ast }$. [This last condition of reciprocity is not really necessary; we can replace it by the fact that when all the baths are at the same temperature, the heat transfer must be zero, to be consistent with the second law of thermodynamics.] (iv) The integrand is even in $\omega$. With these manipulations, the expression can be simplified to a standard Caroli form,

(10)$\begin{array}{l}{I}_L={\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\hslash \omega \mathrm{T}\mathrm{r}(D_{LR}^r{\mathrm{\Gamma }}_R{D}_{RL}^a{\mathrm{\Gamma }}_L)(N_L-N_R).}\end{array}$

Here we define ${\mathrm{\Gamma }}_L=\mathrm{i}({\mathrm{\Pi }}_L-{\mathrm{\Pi }}_L^{†})$ and similarly for ${\mathrm{\Gamma }}_R$, and $N_L$ is the Bose function at temperature $T_L$, and $N_R$ is at $T_R$. We shall call this result the Landauer/Caroli formula, although the original formula was for mesoscopic electron transport [57,74,75].

Let us consider the simplest possible case of a parallel plate capacitor represented as two dots [76] with the Dyson equation in the form

(11)$\begin{array}{l}{\left(D^r\right)}^{-1}=(\begin{array}{rr}C& -C\\ -C& C\end{array})-(\begin{array}{cc}{\mathrm{\Pi }}_L& 0\\ 0& {\mathrm{\Pi }}_R\end{array}),\end{array}$

where $C={ϵ}_{0}A/d$ is the capacitance of the parallel plate capacitor. The retarded Green’s function can be solved explicitly, given

(12)$\begin{array}{l}{D}^r={\displaystyle \frac{1}{{\mathrm{\Pi }}_L{\mathrm{\Pi }}_R-({\mathrm{\Pi }}_L+{\mathrm{\Pi }}_R)C}}(\begin{array}{cc}C-{\mathrm{\Pi }}_R& C\\ C& C-{\mathrm{\Pi }}_L\end{array})\end{array}$

as a $2\times 2$ matrix. With this, the heat transfer takes the Landauer form, Eq. (10), with the transmission function given by

(13)$\begin{array}{l}\mathrm{T}\mathrm{r}(D^r{\mathrm{\Gamma }}_R{D}^a{\mathrm{\Gamma }}_L)={\displaystyle \frac{4C^2\mathrm{I}\mathrm{m}{\mathrm{\Pi }}_L\mathrm{I}\mathrm{m}{\mathrm{\Pi }}_R}{|{\mathrm{\Pi }}_L{\mathrm{\Pi }}_R-({\mathrm{\Pi }}_L+{\mathrm{\Pi }}_R)C{|}^2}}.\end{array}$

This formula tells us, at long distances, the heat transfer is proportional to $1/d^2$ and is constant for short distances, simply because the capacitance $C\propto 1/d$. The crossover distance of the two types of behaviors is controlled by the energy scale $E_{\alpha }\sim e^2/C$ [see Eq. (14) below]. The distance dependence behavior also appears in many more realistic systems, such as between two graphene sheets [77,78].

A simple model for ${\mathrm{\Pi }}_{\alpha }$ with $\alpha =L,R$ can be given based on the analytic properties over frequency $\omega$. The polarizability ${\mathrm{\Pi }}_{\alpha }$ is from a correlation of real quantities, thus the real part must be symmetric in $\omega$ and the imaginary part must be antisymmetric in $\omega$, so for low enough thermal energy $k_{B}T$, we can take

(14)$\begin{array}{l}{\mathrm{\Pi }}_{\alpha }=-{\displaystyle \frac{e^2}{E_{\alpha }}}-\mathrm{i}{\displaystyle \frac{e^2\hslash \omega }{U_{\alpha }^2}},\end{array}$

where $e$ is the magnitude of electron charge, $E_{\alpha }$ and $U_{\alpha }$ are some constants of order eV. This assumption for ${\mathrm{\Pi }}_{\alpha }$ leads to a robust Stefan-Boltzmann law of $T^4$ with a coefficient that depends on the details. Fig.1 illustrates the results of heat current as a function of distance and temperature for the simple capacitor model.

Fig.1 Heat current of the simple capacitor model, from Eqs. (10), (13), and (14). (a) Current $I_L$ vs. distance $d$ at $T_L=1000$ K and $T_R=300$ K. (b) The current $I_L$ vs. $T_L$ while fixing $d=1$ nm, $T_R=300$ K. Parameters are area $A=1$ nm2, energy scale $E_L=E_R=U_L=U_R=1$ eV.

Full size|PPT slide

Is the assumption that the left dot and right dot are uncorrelated realistic? If the capacitor is connected to the same battery, the individual fluctuations of the charges on the two plates seem to violate charge conservation. However, if the systems are more complex than two dots, the charge conservation can be fulfilled by moving the charge on the same side from one place to another. In fact, the derived Landauer/Caroli formula is generally valid for a tight-binding model of a system with $v_{ij}=1/(4\pi {ϵ}_0{r}_{ij})$ as the bare Coulomb interaction between two electrons at site $i$ and $j$ with a distance $r_{ij}$. Here a conductor becomes a hopping site. Although the validity to a two-dot capacitor can be doubtful, the validity of the theory for general electron systems with two sides not directly connected can be proved rigorously. This will be the subject of the next section.

3 Coulomb interaction model

In this section, we consider the following model of $N$ pieces of metal objects described by a free electron Hamiltonian $c^{†}Hc=\sum _{\alpha }{c}_{\alpha }^{†}{H}_{\alpha }{c}_{\alpha }$, plus a Coulomb interaction. We assume that each object labeled by $\alpha$ is not directly connected to other objects such that electric currents are absent so that we can focus on the energy (or heat) transport among the objects. This means that the matrix $H$ is block-diagonal. An example of such a system is two graphene sheets with a vacuum gap. Although the electrons cannot jump from one object to another, they can still exchange energy through a Coulomb interaction. We take the electron to be spinless, and the Coulomb term is

(15)$\begin{array}{l}{\displaystyle \frac{e^2}{2}\sum _{i,j}\frac{1}{4\pi {ϵ}_0|{\mathit{r}}_i-{\mathit{r}}_j|}{c}_i^{†}{c}_j^{†}{c}_j{c}_i.}\end{array}$

Here $i$ and $j$ run over all the sites of all the objects. This is just the standard sum of Coulomb interaction between two charges at the position ${\mathit{r}}_i$ and ${\mathit{r}}_j$ in the second quantization notation, excluding the self-interaction terms.

If our system of interest is finite, it is impossible to maintain a local thermal equilibrium and a steady-state heat transport. To realize heat transport, we connect each object to a bath so that an infinite amount of energy can be supplied. This is done in the nonequilibrium Green’s function approach (NEGF) [63] by giving a self-energy of the bath to the degrees of the particular object while the bath is maintained in thermal equilibrium at temperature $T_{\alpha }$ and chemical potential ${\mu }_{\alpha }$. Thus, the object $\alpha$ before turning on the Coulomb interaction and in local thermal equilibrium is described by two Green’s functions, the retarded one, according to

(16)$\begin{array}{l}{G}_{\alpha }^r(E)=[(E+\mathrm{i}\eta )I-H_{\alpha }-{\mathrm{\Sigma }}_{\alpha }^r{]}^{-1},\end{array}$

where $I$ is the identity matrix, $\eta >0$ is an infinitesimal quantity describing the electron relaxation, and the lesser Green’s function describing the electron occupation,

(17)$\begin{array}{l}{G}_{\alpha }^{<}(E)=-f_{\alpha }(E)[G_{\alpha }^r(E)-G_{\alpha }^a(E)],\end{array}$

where $G_{\alpha }^a={\left(G_{\alpha }^r\right)}^{†}$ is the advanced Green’s function, and $f_{\alpha }(E)=1/({\mathrm{e}}^{{\beta }_{\alpha }(E-{\mu }_{\alpha })}+1)$ is the Fermi distribution function of the bath connected to object $\alpha$, ${\mu }_{\alpha }$ is the chemical potential. If the object is not in thermal equilibrium, e.g., one object connected to two baths, we need to use the Keldysh equation $G^{<}=G^r{\mathrm{\Sigma }}^{<}{G}^a$ instead [65]. The electron Green’s functions are our starting point to characterize the materials properties, such as the polarizability $\mathrm{\Pi }$.

We now derive the Meir−Wingreen formula [79,80] for the energy current which is an exact formula if the Coulomb problem can be solved exactly. To derive the equation, we first make a notational simplification and consider only one bath and one object with the one-particle Hamiltonian

(18)$\begin{array}{l}H=H_d+V=(\begin{array}{cc}{H}^O& 0\\ 0& H^B\end{array})+(\begin{array}{cc}0& V^{OB}\\ V^{BO}& 0\end{array}),\end{array}$

here $H^O$ is the Hamiltonian of the object, and $H^B$ is the isolated bath, and $V$ is the coupling between the object and the bath. The energy transfer per unit time out of the bath is obtained by the average decrease of the bath Hamiltonian,

(19)$\begin{array}{rl}{I}_B& =-⟨\frac{\mathrm{d}{\hat{H}}^B}{\mathrm{d}t}⟩=\frac{1}{\mathrm{i}\hslash }⟨[\hat{V},{\hat{H}}^B]⟩\\ & =⟨c^{†}{V}^{OB}\dot{d}+{\dot{d}}^{†}{V}^{BO}c⟩.\end{array}$

Here we use a convention that the hatted operators are in second quantization notations with the creation operators multiplied from the left and the annihilation operators multiplied from the right to the one-particle matrix elements of the same letter without the hat, e.g., ${\hat{H}}^B=d^{†}{H}^{B}d$. We use $c,c^{†}$ for the object and $d,d^{†}$ for the bath, all of these are column or row vectors. We have used the Heisenberg equation for the rate of changes. Here the dot means evolution by ${\hat{H}}^B$ only, holding the system fixed. We can replace it by the full evolution of the full Hamiltonian, including a possible Coulomb interaction at the dot or center (but not at the bath or between dot and bath), we obtain

(20)$\begin{array}{l}\dot{d}={\displaystyle \frac{\mathrm{d}d}{\mathrm{d}t}}-{\displaystyle \frac{1}{\mathrm{i}\hslash }}{V}^{BO}c.\end{array}$

Here $\mathrm{d}/\mathrm{d}t$ represents the full evolution. Similarly, we obtain ${\dot{d}}^{†}$ by Hermitian conjugation. Putting this result into the energy current expression, we see that $c^{†}{V}^{OB}{V}^{BO}c$ terms cancel, we obtain

(21)$\begin{array}{l}{I}_B=⟨c^{†}{V}^{OB}{\displaystyle \frac{\mathrm{d}d}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{d}^{†}}{\mathrm{d}t}}{V}^{BO}c⟩.\end{array}$

This expression can be further expressed in terms of the lesser Green’s function connecting the bath with the object, $G_{kj}^{OB,<}(t,t^{\prime })=-⟨d_j^{†}(t^{\prime })c_k(t)⟩/(\mathrm{i}\hslash )$, or the reversed version $G^{BO}$. A key step is to remove the reference to the bath variables by using the Green’s function of the object only. It will have taken us too long here to repeat this argument, thus we refer to the standard NEGF textbooks [63,81] on this, given then by the so-called Langreth rule [59,82]

(22)$\begin{array}{l}{G}^{BO,<}=g_B^r{V}^{BO}{G}^{<}+g_B^{<}{V}^{BO}{G}^a,\end{array}$

here the small $g_B$ is the Green’s function of the bath when it is isolated from the object, and capital $G$ is the Green’s function of the object when it is interacting with the bath (as well as internal Coulomb interaction and other baths not in our focus). The product of two Green’s functions means a convolution in time domain and a multiplication in the energy domain. Then we have [83], using the Fourier representation of Green’s functions at time $t=0$,

(23)$\begin{array}{rl}{I}_B& =-\mathrm{i}\hslash \frac{\mathrm{\partial }}{\mathrm{\partial }t}\mathrm{T}\mathrm{r}\left(V^{OB}{G}^{BO,<}(t)-G^{OB,<}(t)V^{BO}\right){|}_{t=0}\\ & ={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\mathrm{d}E}{2\pi \hslash }E\mathrm{T}\mathrm{r}\left((G^r-G^a){\mathrm{\Sigma }}_B^{<}-G^{<}({\mathrm{\Sigma }}_B^r-{\mathrm{\Sigma }}_B^a)\right)\\ & ={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\mathrm{d}E}{2\pi \hslash }E\mathrm{T}\mathrm{r}(G^{>}{\mathrm{\Sigma }}_B^{<}-G^{<}{\mathrm{\Sigma }}_B^{>}).\end{array}$

We have used the fact that the bath self-energy is related to the isolated bath Green’s function by ${\mathrm{\Sigma }}_B=V^{OB}{g}_B{V}^{BO}$, and used the relations among Green’s functions that $G^r-G^a=G^{>}-G^{<}$, also valid for the self-energies. For a multiple-lead (or bath) setup, we simply replace $B$ by a more general index $\alpha$ for the bath or object $\alpha$. The formalism is also valid if we connect several baths to the same object, in order to establish a situation that is not in local equilibrium for the object.

As we can see from the derivation, we have not used specific properties of the objects other than the fact that the couplings between the object and bath are bilinear in the creation and annihilation operators. As a result, the Meir-Wingreen formula is valid, having the Coulomb interaction or not. However, the Coulomb interaction is one of the hardest problems in condensed matter physics, and it cannot be solved exactly. Here we use the simplest and also standard approximation for the Coulomb interaction self-energy ${\mathrm{\Sigma }}_n$ by the Fock term and for the electron screening $\mathrm{\Pi }$ by the random phase approximation (RPA) [84]. The main point below is to get rid of the electron Green’s functions, and to relate the energy current to the (scalar) photon Green’s functions and polarizability. From this, we obtain a similar Meir−Wingreen formula in terms of the photon Green’s function $D$ and polarizability ${\mathrm{\Pi }}_{\alpha }$ of object $\alpha$. Further approximation with the local thermal equilibrium gives the Landauer/Caroli formula. The Meir−Wingreen form is more general as no local equilibrium is assumed (although each object still needs to be connected to one or more baths, the bath is in thermal equilibrium).

To make progress with the Coulomb problem in energy transport, one further approximation is needed, which is the lowest order expansion approximation [85] in terms of the Coulomb interaction. Such approximation preserves energy conservation exactly. An alternative to this is the self-consistent Born approximation through iterative solutions; we will make a comment on this latter approach at the end of this section. For notational simplicity, we will here consider two-bath situation with two objects, call them $L$ and $R$. The Green’s functions for the electrons and self-energies are block diagonal since two sides are not directly connected, and the Meir-Wingreen formula needs only the block $\alpha$. We focus on the object $L$, and Green’s function $G_L^{>}$ can be expressed by the Keldysh equation as

(24)$\begin{array}{l}{G}_L^{>}=[G^r({\mathrm{\Sigma }}_L^{>}+{\mathrm{\Sigma }}_R^{>}+{\mathrm{\Sigma }}_n^{>})G^a{]}_{LL}.\end{array}$

Here ${\mathrm{\Sigma }}_{L,R}^{>}$ are the lead self-energies, and ${\mathrm{\Sigma }}_n^{>}$ is the Fock term of Coulomb interaction. For the $LL$-subblock, ${\mathrm{\Sigma }}_R$ is 0 since bath $R$ is not directly connected to object $L$. This expression requires the knowledge of the full retarded Green’s function $G^r$ with Coulomb interaction. We prefer to work with the free electron Green’s functions. An approximation we can use is the lowest order expansion,

(25)$\begin{array}{l}{G}^{>}\approx G_0^{>}+G_0^r{\mathrm{\Sigma }}_n^r{G}_0^{>}+G_0^r{\mathrm{\Sigma }}_n^{>}{G}_0^a+G_0^{>}{\mathrm{\Sigma }}_n^a{G}_0^a.\end{array}$

We obtain such terms if we expand the contour ordered Dyson equation, $G=G_0+G_0{\mathrm{\Sigma }}_{n}G\approx G_0+G_0{\mathrm{\Sigma }}_n{G}_0+\cdots$, and then take the greater component using the Langreth rule. The subscript $0$ denotes the left system that is free of Coulomb interaction, i.e., a perfect ballistic system with quadratic Hamiltonian $c^{†}Hc$ (with baths). We drop the subscript 0 from now on.

It is useful for symmetry reason we express the current by vacuum diagrams in time domain. We use the inverse Fourier transform to change the integral in energy to time, and also the Fock diagram result [83,86],

(26)$\begin{array}{l}{\mathrm{\Sigma }}_n^{>}(t,t^{\prime })=\mathrm{i}\hslash {\displaystyle \sum _{l,l^{\prime }}{M}^l{G}^{>}(t,t^{\prime })M^{l^{\prime }}{D}_{l{l}^{\prime }}^{>}(t,t^{\prime }).}\end{array}$

A matrix multiplication, $MGM$, is implied in the electron index space. Similar expressions are given for the retarded and advanced self-energies as ${\mathrm{\Sigma }}_n^r\propto G^r{D}^{<}+$$G^{>}{D}^r$, and ${\mathrm{\Sigma }}_n^a\propto G^a{D}^{<}+G^{>}{D}^a$. Here for generality, we assume that the interaction bare vertex takes the form $\sum _l{c}^{†}{M}^{l}c{\varphi }_l$, where $c$ is a column vector of electron annihilation operators and $c^{†}$ is a row vector of Hermitian conjugate, and ${\varphi }_l$ is a scalar field at site $l$, and $M^l$ is a Hermitian matrix. We will explain the scalar field $\varphi$ in more detail in the next section. Here it is sufficient to know that the electron Fock diagram is given by $GD$ (in electron many-body theory $D=v+v\mathrm{\Pi }D$ is called screened Coulomb $W$ where $v$ is the bare Coulomb interaction). We have also ignored the Hartree diagram as the Hartree term only shifts the energy levels of single-particle Hamiltonian and does not contribute to transport in our order of approximation. By plugging Eq. (25) into (23), the expansion leads to 10 terms, represented by the 10 diagrams in Fig.2. We will label these diagrams as 1 to $5$, and $1^{\prime }$ to $5^{\prime }$, as shown. The diagrammatic rule follows the usual convention with all the (real) times as dummy integration variables and space indices summed. The current is obtained by $(\mathrm{i}\hslash {)}^2/T$ times the value of the diagram. Since all the times are integration variables on an equal footing, the integral diverges due to time translational invariance, the $1/T$ factor cancels the last integral interpreted as ${\int }_{-T/2}^{T/2}\mathrm{d}t\cdots$. As an example, the graph 3 represents the contribution to current as

Fig.2 Diagrams for the heat current in the lowest order expansion. The solid lines denote the electron Green’s functions, and the dash lines for the photon Green’s functions. The boxes denote derivative of the bath self energy. The end of the arrow is associated with the first argument and beginning of the arrow with the second argument of the Green’s functions or self energies. The electron−photon vertex is associated with the matrix $M^l$.

Full size|PPT slide

(27)$\begin{array}{rl}3)& =\frac{(\mathrm{i}\hslash {)}^2}{T}\int \mathrm{d}t\mathrm{d}{t}^{\prime }\mathrm{d}{t}_1\mathrm{d}{t}_2\sum _{l,l^{\prime }}{D}_{l{l}^{\prime }}^{>}(t,t^{\prime })\\ & \times \mathrm{T}\mathrm{r}\left[M^l{G}^{>}(t,t^{\prime })M^{l^{\prime }}{G}^a(t^{\prime },t_1)\frac{\mathrm{\partial }{\mathrm{\Sigma }}^{<}(t_1,t_2)}{\mathrm{\partial }{t}_1}{G}^r(t_2,t)\right].\end{array}$

Note the partial derivative on the first argument of ${\mathrm{\Sigma }}^{<}$, which is represented by a dot in the diagram. The partial derivative can be moved around with repeated integration by parts. The trouble with this expression is that $G^a{\mathrm{\Sigma }}^{<}{G}^r$ is in the wrong order to be a $G^{<}$.

A crucial identity [57],

(28)$\begin{array}{rl}{G}^r({\mathrm{\Sigma }}^{>}-{\mathrm{\Sigma }}^{<})G^a& =G^a({\mathrm{\Sigma }}^{>}-{\mathrm{\Sigma }}^{<})G^r\\ & =G^r-G^a=G^{>}-G^{<},\end{array}$

is needed to show that the 10 diagrams cancel among themselves and reduce to only two with the correct order. Here the self-energies are total lead self-energy (for object $L$ only). This identity is a simple consequence of the Dyson equation $(G^r{)}^{-1}=(g_c^r{)}^{-1}-{\mathrm{\Sigma }}^r$, where $g_c^r$ is the Green’s function of the isolated center. From the above equation we can show that

(29)$\begin{array}{l}{G}^a{\mathrm{\Sigma }}^{<}{G}^r=G^{<}+C,\end{array}$

here we define $C=G^a{\mathrm{\Sigma }}^{>,<}{G}^r-G^r{\mathrm{\Sigma }}^{>,<}{G}^a$, and is the same for greater and lesser components. $C$ is anti-Hermitian, $C^{†}=-C$. $C=0$ if matrices are actually 1 by 1, or if the system is reciprocal (i.e., $(G^r{)}^{\mathrm{T}}=G^r$, $(G^{<}{)}^{\mathrm{T}}=G^{<}$), but not so in general. From this, ignoring the proportionality constant, integration variables, and $M$ factors, we can write, symbolically,

(30)$\begin{array}{l}\mathrm{\Delta }3+\mathrm{\Delta }{3}^{\prime })=\mathrm{T}\mathrm{r}[(D^{>}{G}^{>}-D^{<}{G}^{<})C].\end{array}$

Here the notation $\mathrm{\Delta }$ means that the term when $G^a$ and $G^r$ are swapped to form $G^{>}$ or $G^{<}$ has been subtracted off. We show that Eq. (30) cancels all the other 8 diagrams. To this end, we define

(31)$\begin{array}{l}B=G^{>}{\mathrm{\Sigma }}^{<}-G^{<}{\mathrm{\Sigma }}^{>}.\end{array}$

Using the same identities, we have $B{G}^r=-C$, thus $B{G}^r-G^a{B}^{†}=-2C$, and $B{G}^r+G^a{B}^{†}=0$.

We can factor out common factors in the remaining diagrams. Using $B$, we can write

$\begin{array}{rl}& 1+1^{\prime })+2+2^{\prime })=D^{<}\mathrm{T}\mathrm{r}(G^{r}B{G}^r)+D^r\mathrm{T}\mathrm{r}(G^{>}B{G}^r),\\ & 4+4^{\prime })+5+5^{\prime })=D^{<}\mathrm{T}\mathrm{r}(G^a{G}^a{B}^{†})+D^a\mathrm{T}\mathrm{r}(G^{>}{G}^a{B}^{†}).\end{array}$

Further simplification is possible because

(32)$\begin{array}{l}{D}^{<}{G}^r+D^r{G}^{>}=D^{>}{G}^{>}-D^{<}{G}^{<}+D^{<}{G}^a+D^a{G}^{>}.\end{array}$

Now, putting all the terms together, and using the identities obtained, we see $\mathrm{\Delta }3+\mathrm{\Delta }{3}^{\prime })$ cancels all the rest as claimed.

The remaining two terms in the correct $G^r{\mathrm{\Sigma }}^{>,<}{G}^a$ order can be transformed into the desired form. First, we need to move the derivative to other places, for example, from graph $3-\mathrm{\Delta }3)$, we can write

(33)$\begin{array}{l}-D^{>}(t,t^{\prime })\mathrm{T}\mathrm{r}\left[G^{>}(t,t^{\prime }){\displaystyle \frac{\mathrm{\partial }}{\mathrm{\partial }t}}{G}^{<}(t^{\prime },t)\right].\end{array}$

The extra minus sign is due to an integration by parts. We can combine a similar term from $3^{\prime }-\mathrm{\Delta }{3}^{\prime })$ so that it becomes $\mathrm{\partial }{\mathrm{\Pi }}^{<}(t^{\prime },t)/\mathrm{\partial }t$, using integration by parts and cyclic permutation of trace whenever needed. We define the polarization or charge-charge correlation as

(34)$\begin{array}{rl}{\mathrm{\Pi }}_{l^{\prime }l}^{<}(t^{\prime },t)& =\frac{1}{\mathrm{i}\hslash }⟨q_l(t)q_{l^{\prime }}(t^{\prime })⟩\\ & =-\mathrm{i}\hslash \mathrm{T}\mathrm{r}[M^{l^{\prime }}{G}^{<}(t^{\prime },t)M^l{G}^{>}(t,t^{\prime })],\end{array}$

where the charge operator is $q_l=c^{†}{M}^{l}c$, and the second line is obtained assuming that Wick’s theorem [87,88] is valid. This is called RPA. ${\mathrm{\Pi }}^{>}$ is obtained by swapping the positions of the charge operators and by swapping $G^{<}\leftrightarrow G^{>}$. Fourier transforming the final expression to frequency domain, we obtain [89]

(35)$\begin{array}{l}{I}_L=-{\displaystyle \frac{1}{4\pi }{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{d}\omega \hslash \omega \mathrm{T}\mathrm{r}(D^{>}{\mathrm{\Pi }}_L^{<}-D^{<}{\mathrm{\Pi }}_L^{>}).}\end{array}$

The derivation under the lowest order of expansion for the electron Green’s function is rather long and complicated. In fact, if we take a self-consistent Born approximation [83] (also known as the self-consistent $GW$ method), we obtain the same result, with $\mathrm{\Pi }$ being computed by a self-consistent $G$. To this end, we note that since we do not allow the electrons to jump from one object to another, the Green’s functions and electron self-energies are block diagonal. We can treat each as a separate system, although they are connected through $D$. As a result, we have the following identity [63],

(36)$\begin{array}{l}\mathrm{T}\mathrm{r}\left(G_{\alpha }^{>}{\mathrm{\Sigma }}_{\alpha }^{<,\mathrm{t}\mathrm{o}\mathrm{t}}-G_{\alpha }^{<}{\mathrm{\Sigma }}_{\alpha }^{>,\mathrm{t}\mathrm{o}\mathrm{t}}\right)=0,\end{array}$

(37)$\begin{array}{l}{\mathrm{\Sigma }}_{\alpha }^{>,<\mathrm{t}\mathrm{o}\mathrm{t}}={\mathrm{\Sigma }}_{\alpha }^{>,<}+{\mathrm{\Sigma }}_{n,\alpha }^{>,<}.\end{array}$

Since the isolated center is non-dissipative, the equation reflects the conservation laws. The self-energy is the total, including baths and nonlinear term, for object $\alpha$. Using this identity, starting with the Meir−Wingreen formula in $G$ and bath self-energy ${\mathrm{\Sigma }}_L$, Eq. (23), we can replace the bath self-energy with the nonlinear Coulomb self-energy, keeping $G^{>}$ or $G^{<}$ as it is. The next step is to use the self-consistent Fock term to replace the nonlinear self-energy. This results in a double integral in energy $E$ and frequency $\omega$. With a change of variable, such as $E\to E-\hslash \omega$, and the symmetry of the Green’s function $(D^{>}(\omega ){)}^T=D^{<}(-\omega )$, we can rewrite the result in the form of Eq. (35).

Finally, if we can assume local thermal equilibrium for each object, i.e., ${\mathrm{\Pi }}_{\alpha }^{<}=-\mathrm{i}{N}_{\alpha }{\mathrm{\Gamma }}_{\alpha }$, ${\mathrm{\Pi }}_{\alpha }^{>}=-\mathrm{i}(N_{\alpha }+1){\mathrm{\Gamma }}_{\alpha }$, $\alpha =L,R$, together with the Keldysh equation $D^{>,<}=$ $D^r\sum _{\alpha }{\mathrm{\Pi }}_{\alpha }^{>,<}{D}^a$, the equation can be put into the Landauer/Caroli form, as [90]

(38)$\begin{array}{l}{I}_L={\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\hslash \omega \mathrm{T}\mathrm{r}(D^r{\mathrm{\Gamma }}_R{D}^a{\mathrm{\Gamma }}_L)(N_L-N_R),}\end{array}$

here we have used a fact that the integrand is an even function of the frequency $\omega$, thus, we can perform integration over the positive frequency, and times 2. A multiple lead formula of Büttiker type [91] is obtained if we replace $R$ by a summation index ${\alpha }^{\prime }$, i.e., for the current out of lead $\alpha$, it is $I_{\alpha }={\int }_0^{\mathrm{\infty }}\mathrm{d}\omega /(2\pi )\hslash \omega \cdot$ $\sum _{{\alpha }^{\prime }}\mathrm{T}\mathrm{r}(D^r{\mathrm{\Gamma }}_{{\alpha }^{\prime }}{D}^a{\mathrm{\Gamma }}_{\alpha })(N_{\alpha }-N_{{\alpha }^{\prime }})$.

If we apply the formula to a parallel plate geometry, the result is identical to the usual fluctuational electrodynamics result [15-17], taking the speed of light to infinity. Unlike the original Meir−Wingreen formula in terms of the electron Green’s function $G$ and bath self-energies ${\mathrm{\Sigma }}_B$, this final expression cannot be an exact result (even in the non-retardation limit) mainly due to approximations in $\mathrm{\Pi }$. Let us recapitulate the approximations going into the derivations. i) From the Hedin equations [92] point of view, which give a formally exact solution to the Coulomb problem, we have disregarded vertex corrections, so that the electron interaction self-energy is of the Fock form $GD$, and the polarization is the RPA form $GG$. We also ignored the Hartree term. ii) To compute the transport quantities, we have further used the lowest order of expansion in terms of the interaction strength $e$. Alternatively, we can also use the more accurate self-consistent $GW$ [93]. iii) To obtain the Landauer form, we have assumed local equilibrium for each object. iv) The vertex correction as well as the Hartree potential can be incorporated with Eq. (35) remaining valid. The most severe limitation appears to be the assumption of additivity of $\mathrm{\Pi }$. At the RPA level of approximation, $\mathrm{\Pi }$ is block-diagonal; but this is not true at higher order in charge $e$. It will be helpful to think of these approximations used here, when we assess the validity of the usual fluctuational electrodynamics approach.

4 Scalar field model

The usual approach to study the Coulomb interaction [94] in condensed matter physics is to eliminate the scalar field $\varphi$ and to focus on the instantaneous charge-charge interaction as given by expression (15). The quantization of the electromagnetic field is only for the transverse vector field $\mathit{A}$. One disadvantage in this approach is that the (scalar) photon Green’s function $D$ appears to be somewhat ad hoc as a convoluted construction out of the electron screening perspective, unrelated to the scalar field. Here in this section, we take the scalar $\varphi$ as fundamental [76,95] and define $D$ in terms of $\varphi$ in the usual way of NEGF. As a result, the electron Green’s function $G$ and photon Green’s function $D$ stand on an equal footing. This also makes the symmetry properties of the Green’ function $D$ transparent. The only technical problem with this is that $\varphi$, because it is instantaneous, does not have a free field dynamics.

To overcome the difficulty of no conjugate momentum for $\varphi$, we introduce a fictitious speed of light $\tilde{c}$. Then the usual machinery of quantum field theory applies [1,96]. We take the limit that $\tilde{c}$ goes to infinity at the end of the calculation. In this limit, the theory becomes equivalent to the instantaneous Coulomb problem. But it turns out that the free field energy must be negative definite, in order to be consistent with the Poisson equation, and removing the zero-point motion contribution to the Poynting vector (in the $c\to \mathrm{\infty }$ limit, it becomes the “Poynting” scalar), we must use an anti-normal order, instead of the usual normal order. A photon bath at infinity must have a negative temperature. All these exotic features actually disappear and have no physical consequences if, at the end of the calculation, we take the limit $\tilde{c}\to \mathrm{\infty }$ [95].

The electron and scalar photon coupled system can be described by the following Lagrangian:

(39)$\begin{array}{rl}L=& \int \mathrm{d}V\frac{{ϵ}_0}{2}\left(-\frac{{\dot{\varphi }}^2}{{\tilde{c}}^2}+(\mathrm{\nabla }\varphi {)}^2\right)\\ & +c^{†}\left(\mathrm{i}\hslash \frac{\mathrm{d}c}{\mathrm{d}t}-Hc\right)+e\sum _j{c}_j^{†}{c}_j\varphi ({\mathit{r}}_j).\end{array}$

Here the first term is for the free photons, the second term is for the free electrons, and the last term is the electron-photon interaction (electron charge is $-e$). This last term is expressed as $-\sum _l{c}^{†}{M}^{l}c{\varphi }_l$ in the previous section. This Lagrangian is equivalent to a Lorenz gauge choice with the vector field setting to zero. While the electrons are allowed to sit on a set of discrete lattice sites with a tight-binding Hamiltonian, the field exists in the whole space and the contribution to the Lagrangian is an integration over all $\mathit{r}$. We notice that the kinetic energy of the scalar photons is negative. With this Lagrangian, the variational principle, $\delta \int L\mathrm{d}t=0$, reproduces the Poisson equation for $\varphi$ and the Schrödinger equation in a field of $\varphi$ for the electrons. Here we treat $c^{†}$, $c$, $\varphi$ as independent variables.

The reason we start with the Lagrangian $L$ instead of Hamiltonian is that we need to identify the proper conjugate momenta for the fundamental dynamic variables, and to aid us for the canonical quantization. It is clear that the electrons are already in the quantum form with the anti-commutation relations. Here we mainly focus on the scalar field $\varphi$. The conjugate momenta are given by

(40)$P_{c_j}={\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }{\dot{c}}_j}}=i\hslash c_j^{†},$

(41)$P_{c_j^{†}}={\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }{\dot{c}}_j^{†}}}=0,$

(42)${\mathrm{\Pi }}_{\varphi }(\mathit{r})={\displaystyle \frac{\delta L}{\delta \dot{\varphi }(\mathit{r})}}=-{ϵ}_0{\displaystyle \frac{\dot{\varphi }}{{\tilde{c}}^2}}.$

The $\delta$ in the last equation means a functional derivative, as $L$ depends on $\varphi (\mathit{r})$ functionally. The ominous minus sign, which says that the momentum ${\mathrm{\Pi }}_{\varphi }$ is opposite to the velocity $\dot{\varphi }$, already is an indication of something unusual. From the conjugate momenta we obtain the Hamiltonian as $H=\sum p\dot{q}-L$, or

(43)$\begin{array}{rl}\hat{H}=& -\int \mathrm{d}V\frac{{ϵ}_0}{2}\left[\frac{{\dot{\varphi }}^2}{{\tilde{c}}^2}+(\mathrm{\nabla }\varphi {)}^2\right]\\ & +c^{†}Hc-e\sum _j{c}_j^{†}{c}_j\varphi ({\mathit{r}}_j).\end{array}$

The Hamiltonian for the electrons takes a familiar form with a potential applied to the diagonal elements, but the scalar photons take a strange negative-definite energy. Such a Hamiltonian of the scalar photons gives rise to the “inverted oscillators” [97] for the free field modes. When $\tilde{c}\to \mathrm{\infty }$, due to the existence of the Poisson equation, we can think of the total Coulomb energy as being purely due to charges, summed over each pair once, as in Eq. (15). Alternatively, we could attribute the energy completely to the field, with a positive energy density of $\frac{1}{2}{ϵ}_0(\mathrm{\nabla }\varphi {)}^2$. The Hamiltonian above has a third interpretation: the charge-field interaction is $\sum _j{q}_j{\varphi }_j$, but this overcounts the energy by a factor of 2, thus, we need to subtract a half in the field term. The canonical quantization means that, in addition to the usual fermionic anti-commutation relations, $c_i{c}_j^{†}+c_j^{†}{c}_i={\delta }_{ij}$, $c_i{c}_j+c_j{c}_i=0$, and $c_i^{†}{c}_j^{†}+c_j^{†}{c}_i^{†}=0$, we impose

(44)$\begin{array}{l}[\varphi (\mathit{r}),{\mathrm{\Pi }}_{\varphi }({\mathit{r}}^{\prime })]=\mathrm{i}\hslash \delta (\mathit{r}-{\mathit{r}}^{\prime }).\end{array}$

All the rest of the operators commute, i.e., between $\varphi$s, or between $\mathrm{\Pi }$s, or between the field and fermionic electron operators. These commutation relations completely define the quantum mechanical problem. If the commutation relations are postulated correctly, we should obtain the correct equations of motion in operator form with the Heisenberg equation. In our case, $\ddot{\varphi }=[\dot{\varphi },\hat{H}]/(\mathrm{i}\hslash )$ gives a wave equation with charges of the electrons as source, reducing to the Poisson equation in $\tilde{c}\to \mathrm{\infty }$, and $\mathrm{i}\hslash \dot{c}=[c,\hat{H}]$ is the Schrödinger equation of the electron in a potential of $(-e)\varphi$.

We are now in a position to define the scalar photon Green’s function,

(45)$\begin{array}{l}D(\mathit{r}\tau ,{\mathit{r}}^{\prime }{\tau }^{\prime })={\displaystyle \frac{1}{\mathrm{i}\hslash }}⟨T_{\tau }\varphi (\mathit{r},\tau )\varphi ({\mathit{r}}^{\prime },{\tau }^{\prime })⟩,\end{array}$

here $\mathit{r}$, ${\mathit{r}}^{\prime }$ are space positions, and $\tau$, ${\tau }^{\prime }$ are Keldysh contour times. $T_{\tau }$ is the contour order super-operator. The contour time is the pair $\tau =(t,\sigma )$ of real time and branch index. The average $⟨\cdot ⟩=\mathrm{T}\mathrm{r}(\rho \cdot )$ is unspecified at the moment; it could be a thermal equilibrium or in general of nonequilibrium steady state with a density matrix $\rho$. At this point, we need to pause and ask if we should define the Green’s function by deviations, $\delta \varphi =\varphi -⟨\varphi ⟩$. For noninteracting free fields, since the Hamiltonian is quadratic in $\varphi$, the average $⟨\varphi ⟩$ is zero, the subtraction is not necessary. For an interacting system with an interaction term odd in $\varphi$, without the subtraction, we do not even have a proper Dyson equation. In such case $\varphi$ should be understood to be a deviation from the average. These two definitions differ by a constant. An additive constant term in $D$ independent of the time produces a delta function in the frequency domain, thus in any way does not contribute to transport.

We briefly summarize some of the key definitions and properties of the Green’s functions. Readers unfamiliar with the concept should consult the literature on NEGF [59,63,66,98]. For notational brevity, we often omit the $\mathit{r},{\mathit{r}}^{\prime }$ arguments, and treat them as indexing a matrix. Due to the + (forward) and − (backward) branches the contour Green’s function $D(\tau ,{\tau }^{\prime })=D^{\sigma ,{\sigma }^{\prime }}(t,t^{\prime })$ gives four Green’s functions in real time, $D^{++}=D^t$ is time ordered, $D^{--}=D^{\overline{t}}$ is anti-time ordered, $D^{+-}=D^{<}$ is lesser, and $D^{-+}=D^{>}$ is greater. The four are not linearly independent and are constrained by $D^t+D^{\overline{t}}=D^{>}+D^{<}=D^K$. The symmetric correlation $D^K$ is known as the Keldysh component. The retarded Green’s function is $D^r=$ $D^t-D^{<}=\theta (D^{>}-D^{<})$, where the step function $\theta =\theta (t-t^{\prime })=1$ if $t>t^{\prime }$ and 0 otherwise. And the advanced is $D^a=D^{<}-D^{\overline{t}}=-(1-\theta )(D^{>}-D^{<})$, such that $D^{>}-D^{<}=D^r-D^a$. Letting $1=(\mathit{r}t)$ and $2=({\mathit{r}}^{\prime }{t}^{\prime })$, we have the symmetry in time domain as $D^{>}(1,2)=D^{<}(2,1)$, and $D^r(1,2)=D^a(2,1)$. The Fourier transform into frequency domain is defined by

(46)$\begin{array}{l}D(\omega )={\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{d}tD(t-t^{\prime }){\mathrm{e}}^{\mathrm{i}\omega (t-t^{\prime })}.}\end{array}$

In frequency domain, we have the Hermitian conjugate $D^r(\omega {)}^{†}=D^a(\omega )$, $D^{<}(\omega {)}^{†}=-D^{<}(\omega )$.

The reason for introducing the contour ordered Green’s function is that it facilitates systematic perturbative expansion at any temperature, as in a quantum field theory at zero temperature. The Green’s function $D$ can be determined with a Dyson equation on contour, $D=v+v\mathrm{\Pi }D$, if we know the self-energy $\mathrm{\Pi }$, which we can determine by a standard diagrammatic expansion [99,100]. To the lowest order of approximation, it is given by the RPA expression, Eq. (34). The contour version of this equation means, written out in full,

(47)$\begin{array}{rl}& D(\mathit{r}\tau ,{\mathit{r}}^{\prime }{\tau }^{\prime })=v(\mathit{r}\tau ,{\mathit{r}}^{\prime }{\tau }^{\prime })\\ & +\sum _{j,k}\int \mathrm{d}{\tau }_1\int \mathrm{d}{\tau }_{2}v(\mathit{r}\tau ,{\mathit{r}}_j{\tau }_1){\mathrm{\Pi }}_{jk}({\tau }_1,{\tau }_2)D({\mathit{r}}_k{\tau }_2,{\mathit{r}}^{\prime }{\tau }^{\prime }).\end{array}$

Here we use a mixed representation, while $D$ and $v$ are defined on the whole space, since the electrons are on a set of discrete sites ${\mathit{r}}_j$, the self-energy is defined only on the discrete sites. $v$ is the free Green’s function, i.e., the Green’s function when the electrons are absent. Since the thermal state is essentially fixed by the electrons, we only need the retarded version of $v$, which can be obtained by an equation of motion method. If we compute the derivative with respect to the first argument of time, $t$, twice, we find

(48)$\begin{array}{l}{ϵ}_0\left({\displaystyle \frac{1}{{\tilde{c}}^2}}{\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{t}^2}}-{\mathrm{\nabla }}^2\right)v(\mathit{r}t,{\mathit{r}}^{\prime }{t}^{\prime })=\delta (\mathit{r}-{\mathit{r}}^{\prime })\delta (t-t^{\prime }).\end{array}$

We see that $v$ is essentially the Green’s function of the wave equation. The contour Dyson equation implies a pair of real time equations as the retarded Dyson equation, $D^r=v^r+v^r{\mathrm{\Pi }}^r{D}^r$, and the Keldysh equation, $D^{<}=D^r{\mathrm{\Pi }}^{<}{D}^a$, so $v^{<}$ is never needed (this is because in the limit $\tilde{c}\to \mathrm{\infty }$, the bath at infinity has no effect. Or in other words, the Coulomb field $\varphi$ cannot propagate to infinity).

4.1 “Poynting scalar”

How do we describe the energy flux in a pure scalar field theory? We start from first principles. The energy density, according to our separation of the free field and electron-photon interaction, is $u=-\frac{1}{2}{ϵ}_0\left({\dot{\varphi }}^2/{\tilde{c}}^2+(\mathrm{\nabla }\varphi {)}^2\right)$. From this expression, we can obtain a conservation law in differential equation form,

(49)$\begin{array}{rl}\frac{\mathrm{\partial }u}{\mathrm{\partial }t}& =-{ϵ}_0\left(\frac{\dot{\varphi }\ddot{\varphi }}{{\tilde{c}}^2}+\mathrm{\nabla }\dot{\varphi }\cdot \mathrm{\nabla }\varphi \right)\\ & =-\mathrm{\nabla }\cdot \mathit{j}-\rho \dot{\varphi }.\end{array}$

Here we define the Poynting scalar as $\mathit{j}={ϵ}_0\dot{\varphi }\mathrm{\nabla }\varphi$. The dot denotes partial derivative with respect to time. In obtaining the second line, we note that $\varphi$ satisfies a wave equation with the charge density as a source. Here we treat the expression as a classical equation. To make a quantum equivalent, we need to symmetrize the products. The last term in the expression $\rho \dot{\varphi }$ is for Joule heating. If we are calculating the steady-state work, we can perform “integration by parts,” so the average Joule heating is also $-\dot{\rho }\varphi$, consistent with the expression given earlier in Section 2.

We can express the quantum dynamic variable, i.e., quantum field of $\varphi$, in terms of the creation and annihilation operators in the usual way,

(50)$\begin{array}{l}\varphi (\mathit{r})={\displaystyle \sum _{\mathit{k}}\sqrt{\frac{\hslash {\tilde{c}}^2}{2{ϵ}_0{\tilde{\omega }}_{\mathit{k}}V}}\left(b_{\mathit{k}}{\mathrm{e}}^{\mathrm{i}\mathit{k}\cdot \mathit{r}}+\mathrm{h}.\mathrm{c}.\right),}\end{array}$

except that the commutation relation acquires a minus sign, due to Eq. (44), namely,

(51)$\begin{array}{l}[b_{\mathit{k}},b_{{\mathit{k}}^{\prime }}^{†}]=-{\delta }_{\mathit{k},{\mathit{k}}^{\prime }},[b_{\mathit{k}},b_{{\mathit{k}}^{\prime }}]=0,[b_{\mathit{k}}^{†},b_{{\mathit{k}}^{\prime }}^{†}]=0.\end{array}$

Here the wave vector $\mathit{k}$ is box quantized in a finite volume of $V$, and the mode frequency is given by ${\tilde{\omega }}_{\mathit{k}}=\tilde{c}|\mathit{k}|$. With this transformation, the free field term is an inverted oscillator form, $\int u\mathrm{d}V=-\sum _{\mathit{k}}\hslash {\tilde{\omega }}_{\mathit{k}}(b_{\mathit{k}}{b}_{\mathit{k}}^{†}+1/2)$. Because of the minus sign in the commutation relation, and because of the negative-definite nature of the Hamiltonian, the role of creation and annihilation is reversed. The ground state is defined by $b_{\mathit{k}}^{†}|0⟩=0$. This is because the creation operator still has the meaning of increasing energy by one unit of $\hslash {\tilde{\omega }}_{\mathit{k}}$, and annihilation operator decreasing the energy by one unit, but due to the negative definiteness, one cannot increase forever. The implication of this feature is that in order to remove the zero-point contribution to the Poynting scalar, we must take an anti-normal order.

Let us explain this point more carefully. We see that the Poynting scalar is a quadratic form of $\varphi$. If we expand $\varphi$ in terms of $b$ and $b^{†}$, we get four types of terms, $bb$, $b^{†}{b}^{†}$, $b{b}^{†}$, and $b^{†}b$. The vacuum expectation value $⟨0|\cdot |0⟩$ is not zero for the last type. The anti-normal order is to swap the last case into the third case, removing the zero point contribution from the field. This anti-normal order (denoted by three vertical dots below) is also used on general states, so we calculate the heat flux using Poynting scalar as [95]

(52)$\begin{array}{rl}⟨\mathit{j}⟩& =\frac{1}{2}{ϵ}_0⟨⋮\dot{\varphi }\mathrm{\nabla }\varphi +\mathrm{\nabla }\varphi \dot{\varphi }⋮⟩\\ & ={ϵ}_0{\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{\pi }\hslash \omega \mathrm{R}\mathrm{e}{\mathrm{\nabla }}_{{\mathit{r}}^{\prime }}{D}^{>}(\omega ,\mathit{r},{\mathit{r}}^{\prime }){|}_{{\mathit{r}}^{\prime }=\mathit{r}},\end{array}$

which we can express in terms of the greater Green’s function and integrate over the positive frequencies. Here we transform the Green’s function into frequency domain, then the derivative with respect to time becomes $-i\omega$. The gradient operator is acting on the second argument of the position. After the operation, both positions are set to be equal. The use of greater instead of lesser Green’s function is due to the anti-normal order requirement.

Equation (50) is not the solution to the problem but only defines a transformation from $b$ to $\varphi$. However, for the free field which has the time dependence $b_{\mathit{k}}(t)=b_{\mathit{k}}{\mathrm{e}}^{-\mathrm{i}{\tilde{\omega }}_{\mathit{k}}t}$, we can obtain the free retarded Green’s function from the definition,

(53)$\begin{array}{rl}v(\mathit{r}t,{\mathit{r}}^{\prime }{t}^{\prime })& =\theta (t-t^{\prime })\frac{1}{i\hslash }⟨[\varphi (\mathit{r},t),\varphi ({\mathit{r}}^{\prime },t^{\prime })]⟩\\ & =\frac{{\tilde{c}}^2}{{ϵ}_0}\sum _{\mathit{k}}\frac{\theta (t-t^{\prime })\sin ({\tilde{\omega }}_{\mathit{k}}(t-t^{\prime }))}{{\tilde{\omega }}_{\mathit{k}}V}{\mathrm{e}}^{\mathrm{i}\mathit{k}\cdot (\mathit{r}-{\mathit{r}}^{\prime })}\\ & =\frac{1}{4\pi {ϵ}_0|\mathit{r}-{\mathit{r}}^{\prime }|}\delta (t-t^{\prime }),\tilde{c}\to \mathrm{\infty }.\end{array}$

If we Fourier transform the expression into $\omega$ space, and then take the limit $\tilde{c}\to \mathrm{\infty }$, we recover the Coulomb interaction in $\mathit{k}$ space, as $v(\mathit{k})=1/({ϵ}_0{k}^2)$, and the usual expression for the Coulomb potential in real space. We see also that the free retarded Green’s function is independent of the distribution, no matter what meaning we give to the average. The free retarded Green’s function is determined solely by the equal-time commutators, unrelated to the distribution controlled by $⟨b_{\mathit{k}}{b}_{\mathit{k}}^{†}⟩$, which determines $v^{<}$. We do not need the free $v^{<}$ as mentioned earlier.

4.2 A parallel plate capacitor as two quantum dots, scalar field

In this subsection, we treat the parallel plate capacitor problem again, but this time, by the scalar field approach [76]. Since in a parallel plate capacitor, the field is a function of only one variable, we call it $z$, the transverse direction can be integrated, giving the area of the capacitor $A$. As a result, the retarded Dyson equation, in differential form, $v^{-1}D=1+\mathrm{\Pi }D$, becomes

(54)$\begin{array}{rl}& -{ϵ}_{0}A\left[{\left(\frac{\omega }{\tilde{c}}+\mathrm{i}\eta \right)}^2+\frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2}\right]D^r(z,z^{\prime },\omega )\\ & =\delta (z-z^{\prime })+\sum _{\alpha =L,R}\delta (z-z_{\alpha }){\mathrm{\Pi }}_{\alpha }(\omega )D^r(z_{\alpha },z^{\prime },\omega ).\end{array}$

Here we have transformed the equation into frequency domain, and introduced an infinitesimal small damping $\eta$ so that the solution is the retarded one. This equation can be interpreted as the scalar potential generated by a unit active (external) charge located at $z^{\prime }$, together with the induced extra charges at the electron sites, $z_L=0$, $z_R=d$, due to the linear response to the internal field. Indeed, the induced charge at site $\alpha$ is given by $\delta q_{\alpha }={\mathrm{\Pi }}_{\alpha }\varphi (z_{\alpha })$, and ${\mathrm{\Pi }}_{\alpha }$ is the associated response function.

To compute the energy current density $j$ using the solution of $D$, we note that the dots are coupled to $D$ only at $z=0$ and $d$. As a result, we only need to know $D$ at one of these points. Also, $D^r(z,z^{\prime },\omega )=D^r(z^{\prime },z,\omega )$ is symmetric, thus the following solution for $0<z<d$ is sufficient:

(55)$D^r(z,0,\omega )=\left({\displaystyle \frac{{\gamma }^2-{\lambda }^2}{\gamma }}{\mathrm{\Pi }}_R+2\gamma \mathrm{\Omega }\right){\displaystyle \frac{1}{\mathcal{D}}},$

(56)$D^r(z,d,\omega )=\left[(1-{\gamma }^2){\mathrm{\Pi }}_L+2\mathrm{\Omega }\right]{\displaystyle \frac{\lambda }{\gamma \mathcal{D}}},$

where $\tilde{k}=\frac{\omega }{\tilde{c}}+\mathrm{i}\eta$, $\lambda ={\mathrm{e}}^{\mathrm{i}\tilde{k}d}$, $\gamma ={\mathrm{e}}^{\mathrm{i}\tilde{k}z}$, $\mathrm{\Omega }=\mathrm{i}{ϵ}_{0}A\tilde{k}$, and

(57)$\begin{array}{l}\mathcal{D}=({\lambda }^2-1){\mathrm{\Pi }}_L{\mathrm{\Pi }}_R-2\mathrm{\Omega }({\mathrm{\Pi }}_L+{\mathrm{\Pi }}_R)-4{\mathrm{\Omega }}^2.\end{array}$

To obtain the solution, we set $z^{\prime }=0$ or $d$, assuming backward moving and decaying wave to the left for $z<0$, $r{\mathrm{e}}^{-\mathrm{i}\tilde{k}z}$, standing wave $A{\mathrm{e}}^{\mathrm{i}\tilde{k}z}+B{\mathrm{e}}^{-\mathrm{i}\tilde{k}z}$ in the middle segment, $0<z<d$, and decaying wave to the right for $z>d$, $t{\mathrm{e}}^{\mathrm{i}\tilde{k}z}$. The wave has to be continuous at $0$ and $d$. But the first derivatives are discontinuous with the discontinuities determined by the $\delta$-functions on the right-hand side of the equation. This gives four boundary condition matching algebraic equations, uniquely determining the coefficients.

In the next step, we use the Keldysh equation

(58)$\begin{array}{l}{D}^{>}(z,z^{\prime })={\displaystyle \sum _{\alpha }{D}^r(z,z_{\alpha }){\mathrm{\Pi }}_{\alpha }^{>}{D}^r(z^{\prime },z_{\alpha }{)}^{\ast },}\end{array}$

where we have omitted $\omega$ argument for simplicity and used $D^a=(D^r{)}^{†}$. We can now plug $D^{>}$ into the expression for heat current density, Eq. (52). Since our problem is quasi-one-dimensional, ${\mathrm{\nabla }}_{{\mathit{r}}^{\prime }}$ is just $\mathrm{\partial }/\mathrm{\partial }{z}^{\prime }$. At this point after the space derivative, we can take the limit $\tilde{c}\to \mathrm{\infty }$ and $\eta \to 0^{+}$, the result becomes independent of location $z$, and [76]

(59)$\begin{array}{l}⟨j⟩={\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{\pi A}\text{Re}[\mathrm{i}\hslash \omega |D_{RL}{|}^2({\mathrm{\Pi }}_R^{>}\text{Im}{\mathrm{\Pi }}_L-{\mathrm{\Pi }}_L^{>}\text{Im}{\mathrm{\Pi }}_R)].}\end{array}$

where $D_{RL}=D_{LR}=1/[{\mathrm{\Pi }}_L{\mathrm{\Pi }}_R/C-({\mathrm{\Pi }}_L+{\mathrm{\Pi }}_R)]$, and $C=$ $A{ϵ}_0/d$ is the capacitance. If we assume local thermal equilibrium with the fluctuation-dissipation theorem for ${\mathrm{\Pi }}_{\alpha }^{>}$, we obtain the same Laudauer formula of Section 2. If we adopt a self-consistent calculation, local equilibrium is not a valid assumption, then the above equation is suitable for a self-consistent calculation. We note that the above formula agrees trivially with the photon version of the Meir−Wingreen formula, Eq. (35), and it also agrees with the electron version, Eq. (23), when the self-consistent iterations are converged.

4.3 Quantum dot model of Π

We consider each of the plates as a single quantum dot with maximum charge $Q$, not necessarily the unit charge $e$, with the Hamiltonian ${ϵ}_{\alpha }{c}_{\alpha }^{†}{c}_{\alpha }$ for each dot $\alpha =L$ or $R$. To set up a $GW$ calculation (also known as self-consistent Born approximation, see Fig.3), we start from the retarded Green’s function

Fig.3 The self-consistent Born approximation diagrams. (a) The Hartree self-energy diagram, (b) the Fock self-energy diagram, (c) the polarization bubble $\mathrm{\Pi }$. The double line with arrows signifies the full electron Green’s function $G$, the single dotted line denotes bare Coulomb $v$, while the double dotted line is the screened Coulomb $D$.

Full size|PPT slide

(60)$\begin{array}{l}{G}_{\alpha }^r(E)={\displaystyle \frac{1}{E-{ϵ}_{\alpha }-{\mathrm{\Sigma }}_{\alpha }^r(E)-{\mathrm{\Sigma }}_{n,\alpha }^r(E)}}.\end{array}$

For bath self-energies, we take a phenomenological Lorentz−Drude model [101],

(61)$\begin{array}{l}{\mathrm{\Sigma }}_{\alpha }^r(E)={\displaystyle \frac{{\mathrm{\Gamma }}_{\alpha }/2}{\mathrm{i}+E/E_{0,\alpha }}}.\end{array}$

If $E_{0,\alpha }\to \mathrm{\infty }$, it is called a wide-band model, as the coupling with the lead is independent of energy. Here we cut off to the energy scale of $E_{0,\alpha }$. This gives a better physical picture in real time $t$; the self-energy takes an exponential decay form.

If the dot is not in local thermal equilibrium, we cannot use the fluctuation-dissipation theorem for the electron Green’s function, and have to use the Keldysh equation

(62)$\begin{array}{l}{G}_{\alpha }^{<}(E)=G_{\alpha }^r(E)\left({\mathrm{\Sigma }}_{\alpha }^{<}(E)+{\mathrm{\Sigma }}_{n,\alpha }^{<}(E)\right)G_{\alpha }^a(E).\end{array}$

The lead, since it serves as a bath, by definition is in equilibrium, so we can use ${\mathrm{\Sigma }}_{\alpha }^{<}(E)=-f_{\alpha }(E)({\mathrm{\Sigma }}_{\alpha }^r(E)-$ ${\mathrm{\Sigma }}_{\alpha }^a(E))$, where $f_{\alpha }(E)=1/({\mathrm{e}}^{{\beta }_{\alpha }(E-{\mu }_{\alpha })}+1)$ is the Fermi function. The nonlinear interaction self-energy ${\mathrm{\Sigma }}_{n,\alpha }$, the Hartree and the Fock term coupling the two dots cannot be in local equilibrium, thus no fluctuation-dissipation theorem for it. The nonlinear self-energy can be calculated in real time as

(63)${\mathrm{\Sigma }}_{n,\alpha }^{<}(t)=\mathrm{i}\hslash Q^2{G}_{\alpha }^{<}(t)D^{<}(z_{\alpha },z_{\alpha },t),$

(64)${\mathrm{\Sigma }}_{n,\alpha }^{>}(t)=\mathrm{i}\hslash Q^2{G}_{\alpha }^{>}(t)D^{>}(z_{\alpha },z_{\alpha },t),$

(65)$\begin{array}{rl}{\mathrm{\Sigma }}_{n,\alpha }^r(t)=& \theta (t)({\mathrm{\Sigma }}_{n,\alpha }^{>}(t)-{\mathrm{\Sigma }}_{n,\alpha }^{<}(t))\\ & +(-\mathrm{i}\hslash )Q^2\delta (t)\sum _{{\alpha }^{\prime }}v(z_{\alpha },z_{{\alpha }^{\prime }})G_{{\alpha }^{\prime }}^{<}(t=0).\end{array}$

There is no Hartree contribution to lesser and greater components of the interaction self-energy, since the Hartree term is proportional to $\delta (\tau ,{\tau }^{\prime })$ on contour. We also note that the Hartree potential should be the unscreened one of $v$. For the capacitor model, this is

(66)$\begin{array}{l}v={\displaystyle \frac{1}{2{ϵ}_{0}A\eta }}(\begin{array}{cc}1& {\mathrm{e}}^{-\eta d}\\ {\mathrm{e}}^{-\eta d}& 1\end{array}),\eta \to 0^{+},\end{array}$

which is divergent. Since the Hartree potential only renormalizes the onsite energy, and charge neutrality also requires it to be zero, we drop this term in actual calculation. Finally, the scalar photon self-energy or polarizability is calculated according to

(67)$\begin{array}{l}{\mathrm{\Pi }}_{\alpha }^{>}(t)=-\mathrm{i}\hslash Q^2{G}_{\alpha }^{>}(t)G_{\alpha }^{<}(-t),\end{array}$

(68)$\begin{array}{l}{\mathrm{\Pi }}_{\alpha }^{<}(t)=-\mathrm{i}\hslash Q^2{G}_{\alpha }^{<}(t)G_{\alpha }^{>}(-t),\end{array}$

(69)$\begin{array}{l}{\mathrm{\Pi }}_{\alpha }^r(t)=\theta (t)[{\mathrm{\Pi }}_{\alpha }^{>}(t)-{\mathrm{\Pi }}_{\alpha }^{<}(t)].\end{array}$

With this new $\mathrm{\Pi }$, we need to recalculate $D^r$, and thus new ${\mathrm{\Sigma }}_n$ and $G$, until convergence.

Since the nonlinear self-energy is easy to calculate in time domain, it is natural we adopt a fast Fourier transform method to go between energy and time domain. However, the errors are sometimes hard to control. An alternative to the Fourier transform is to perform convolution in energy space, never going into time domain. Energy domain functions are relatively smooth functions, but time domain function can be highly oscillatory, thus hard to control.

We show in Fig.4 the real and imaginary parts of the retarded $\mathrm{\Pi }$ for the single quantum dot, which is typical of the photon self-energy. Fig.5 presents the distance dependence of heat current density in a double logarithmic plot, calculated with the self-consistent Born approximation. For large distances, the heat current density decreases like $d^{-2}$. Such a scaling law arises from the capacitor property, as mentioned earlier, which can be manifested in the expression for the transmission coefficient. Our choice of the parameters should be considered typical. Fig.5 shows a large enhancement of heat transfer mediated by the scalar photons. Comparing the red solid and blue dashed-dotted lines at $10$ nm, we find that the heat current density is two thousand times larger than the blackbody limit. This result demonstrates that the heat transfer channel provided by electron-photon interaction is the dominant one for nano-capacitors at small separations.

Fig.4 The retarded ${\mathrm{\Pi }}_L$ of the left quantum dot. For parameter, see next Fig.5. Note that the real part is symmetry, while the imaginary part is anti-symmetric with respect to the frequency $\omega$. Reproduced from Ref. [95].

Full size|PPT slide

Fig.5 Energy current calculation under SCBA and the blackbody limit in log−log plot. The temperature of dot $L$ is $1000$ K and dot $R$ is $300$ K. The chemical potential of dot $L$ is $0$ eV and dot $R$ is $0.02$ eV. Charge $Q=1e$, onsite energy ${ϵ}_L={ϵ}_R=0$ eV, and plate area is $A=19.2\times 19.2$ nm2. The bath coupling is ${\mathrm{\Gamma }}_L=1$ eV, ${\mathrm{\Gamma }}_R=0.5$ eV, and $E_{0,L}=E_{0,R}=$ 50 eV. Reproduced from Ref. [95].

Full size|PPT slide

5 Heat transport without local equilibrium by current drive

Earlier in Section 3 we derived a Landauer-like formula and the Meir−Wingreen formula based on Coulomb interaction for an electron system. The Landauer formula assumes the local thermal equilibrium, while the Meir−Wingreen formula is a more general formula without the use of local thermal equilibrium. As an application of the Meir−Wingreen formula, Eq. (35), we consider two graphene layers separated by a vacuum gap, one of them is driven by a current [89]. Since the system is under a current drive, this is truly a nonequilibrium problem and we cannot use the Landauer formula. However, something very close to it does exist, which is a “Doppler” shifted version. The basic idea is the following: under the external field drive, the layer is nearly in equilibrium, satisfying the Kadanoff-Baym ansatz [102] in the sense,

(70)$\begin{array}{l}{G}^{<}=-f(G^r-G^a),\end{array}$

(71)$\begin{array}{l}{G}^{>}=(1-f)(G^r-G^a),\end{array}$

here $f$ is not the equilibrium Fermi function but a shifted one,

(72)$\begin{array}{l}f(ϵ)=f^0-{\displaystyle \frac{\mathrm{d}{f}^0}{\mathrm{d}ϵ}}\mathrm{\Phi }\approx f^0(ϵ-\mathrm{\Phi }).\end{array}$

Here $f^0$ is the usual equilibrium Fermi function, $f^0=1/({\mathrm{e}}^{\beta (ϵ-\mu )}+1)$, while the nonequilibrium distribution $f$ is obtained by deforming the dispersion relation by an amount $\mathrm{\Phi }$. In principle, we should solve the Boltzmann equation for $f$, but a single-mode relaxation approximation which is equivalent to the choice of ${\mathrm{\Phi }}_{\mathit{k}}$ below is more practical and simpler to use. In order for the “fluctuation-dissipation” theorem above to make sense, the Green’s functions are in mode space, i.e., it is for $G^{>,<}(\mathit{k},E)$ and the spectrum function $\mathrm{i}(G^r-G^a)$ is essentially a $\delta$-function of $E$ peaked at the electron band ${ϵ}_{\mathit{k}}$. Only the distributions are shifted; the band stays as it is.

To the lowest order of approximation, we can expand $\mathrm{\Phi }$ by Legendre polynomials, and keep only the most important angular dependence as [103]

(73)$\begin{array}{l}{\mathrm{\Phi }}_{\mathit{k}}=\hslash \mathit{k}\cdot {\mathit{v}}_1=\hslash v_{1}k\cos \theta .\end{array}$

Here ${\mathit{v}}_1$ is the drift velocity and $\theta$ is the angle between the wavevector $\mathit{k}$ and the drift velocity. The drift velocity is related to the electric current by multiplying it by the carrier charge density. Using this version of $G^{>,<}$, we can compute ${\mathrm{\Pi }}^{>,<}$ by the usual formula, which is, in $(\mathit{q},\omega )$ space,

(74)$\begin{array}{l}{\mathrm{\Pi }}_{\mathit{q}}^{<}(\omega )={\displaystyle \frac{e^2}{\mathrm{i}N}\sum _{\mathit{k}}\int \frac{\mathrm{d}E}{2\pi }{G}_{\mathit{k}+\mathit{q}}^{<}(E)G_{\mathit{k}}^{>}(E-\hslash \omega ).}\end{array}$

The real space Green’s function is related to the Fourier space one by $G_l^{>,<}(t)=\frac{1}{N}\sum _{\mathit{k}}\int \frac{\mathrm{d}E}{2\pi \hslash }{G}_{\mathit{k}}^{>,<}$$(E){\mathrm{e}}^{\mathrm{i}(\mathit{k}\cdot {\mathit{R}}_l-Et)/\hslash }$. Here for simplicity of notation, we assume one atom per unit cell at location of Bravais lattice point ${\mathit{R}}_l$, and the summation is over the first Brillouin zone. $N$ is the number of $\mathit{k}$ points. Using the Kadanoff−Baym ansatz and the identity $f^0(ϵ)[1-f^0({ϵ}^{\prime })]=N^0[ϵ-{ϵ}^{\prime })/\hslash ][f^0({ϵ}^{\prime })-f^0(ϵ)]$, one can show that a shifted fluctuation-dissipation theorem also exists, i.e.,

(75)$\begin{array}{l}{\mathrm{\Pi }}_{\mathit{q}}^{<}(\omega )=\tilde{N}(\omega )[{\mathrm{\Pi }}_{\mathit{q}}^r(\omega )-{\mathrm{\Pi }}_{\mathit{q}}^a(\omega )],\end{array}$

and similarly for ${\mathrm{\Pi }}^{>}$ with $\tilde{N}+1$, where we have

(76)$\begin{array}{l}\tilde{N}(\omega )={\displaystyle \frac{1}{{\mathrm{e}}^{\hslash (\omega -\mathit{q}\cdot {\mathit{v}}_1)/(k_{B}T)}-1}}.\end{array}$

This is a Doppler shifted Bose distribution [104-107]. The linearity of ${\mathrm{\Phi }}_{\mathit{k}}$ with respect to $\mathit{k}$ is important here as ${ϵ}^{\prime }-{\mathrm{\Phi }}^{\prime }-ϵ+\mathrm{\Phi }=\hslash \omega -{\mathrm{\Phi }}_{\mathit{k}+\mathit{q}}+{\mathrm{\Phi }}_{\mathit{k}}=\hslash (\omega -\mathit{q}\cdot {\mathit{v}}_1)$ depends only on $\mathit{q}$, otherwise $\tilde{N}$ cannot be factored out of the summation over $\mathit{k}$. Note that due to the replacement $f^0\to f$, the retarded ${\mathrm{\Pi }}^r$ for graphene is no longer Doppler shifting the equilibrium one, but needs to be computed afresh. Plugging in these formulas satisfying a modified version of the fluctuation-dissipation theorem, we still have a Landauer formula, just like Eq. (38), except that the Bose distribution function $N^0(\omega )$ is replaced by the Doppler shifted one of $\tilde{N}(\omega )=$ $N^0(\omega -\mathit{q}\cdot {\mathit{v}}_1)$. However, since the original Meir−Wingreen formula is an integration of the frequency $\omega$ from $-\mathrm{\infty }$ to $+\mathrm{\infty }$, and the Doppler shift breaks the symmetry of $\omega \to -\omega$, the Landauer version also needs to integrate from $-\mathrm{\infty }$ to $+\mathrm{\infty }$ just like the Meir−Wingreen version. This ends our theory for near-field heat transfer under current drive.

To make a concrete calculation, we need to simplify our Caroli expression for transmission, $\mathrm{T}\mathrm{r}(D^r{\mathrm{\Gamma }}_R{D}^a{\mathrm{\Gamma }}_L)$, using the fact that we have a lattice symmetry. Here the trace is over the collection of all electron sites. Since trace is invariant with respect to a similarity transform, it is more efficient if we perform the trace in wave-vector space for the transverse directions using lattice symmetry. Then the trace on the sites becomes a sum over the wavevectors. To simplify the treatment, we ignore possible local inhomogeneous effect and treat each site ${\mathit{R}}_l$ as on a Bravais lattice. What does this amount to is for the self-energy $\mathrm{\Pi }$ we consider all the charges in a unit cell as one single unit, by summing over all the charges in the unit cell as if it is located at the lattice site. The error introduced is negligible when the distance $d$ involved is much larger than the lattice constant. The lattice symmetry implies a translation invariance which means that correlation functions such as $\mathrm{\Pi }$ and $D^r$ are function of the relative distance between two points. We can make a discrete Fourier transform into the two-dimensional $\mathit{q}$ wavevector space, still keeping $z$ the transport direction in real space, as

(77)$\begin{array}{l}D(\mathit{q},z,z^{\prime })={\displaystyle \sum _{l}D({\mathit{R}}_l,z;\mathbf{0},z^{\prime }){\mathrm{e}}^{-\mathrm{i}\mathit{q}\cdot {\mathit{R}}_l}.}\end{array}$

Here ${\mathit{R}}_l=l_1{\mathit{a}}_1+l_2{\mathit{a}}_2$ is a two-dimensional lattice vector running through index $l=(l_1,l_2)$ on the crystal lattice sites, $\mathit{q}$ runs over the reciprocal space in the first Brillouin zone, laying in the plane of graphene. This convention of discrete Fourier transform ensures that the dimensions in real space and in $\mathit{q}$ are the same, i.e., the inverse capacitance, $[D]=[{\mathrm{\Pi }}^{-1}]=[U/e^2]$. As a result of this transformation for $D^r$ as well as for $\mathrm{\Pi }$, we can work in $\mathit{q}$ space, in which each value of $\mathit{q}$ is block-diagonal and we can focus on one particular $\mathit{q}$. The three-dimensional problem reduces to a quasi-one-dimensional problem.

We note that in our original definition of the scalar photon Green’s function, Eq. (45), it is defined on the continuum for any real space $\mathit{r}$. Fortunately, in the Caroli formula, we only need to know the values of $D$ on the lattice sites where electrons exist. As a result, we do not need to solve the Dyson equation covering the whole space $\mathit{r}$ but just over these discrete points $({\mathit{R}}_l,z)$ with $z=0$ for the left and $d$ for the right sheet of graphene. The Dyson equation still takes the form $D^r=v+v\mathrm{\Pi }{D}^r$, or $(D^r{)}^{-1}=v^{-1}-\mathrm{\Pi }$, here $\mathrm{\Pi }={\mathrm{\Pi }}^r$ is the retarded version. The free Green’s function for the transverse directions in $\mathit{q}$ space and $z$ direction in real space is obtained by solving the Poisson equation in mixed representation with a $\delta$ source, as

(78)$\begin{array}{l}\left(q^2-{\displaystyle \frac{{\mathrm{\partial }}^2}{\mathrm{\partial }{z}^2}}\right)v(\mathit{q},z,z^{\prime })={\displaystyle \frac{1}{{ϵ}_0\mathrm{\Omega }}}\delta (z-z^{\prime }).\end{array}$

Here $q=|\mathit{q}|$ is the magnitude of the wave vector, $\mathrm{\Omega }$ is the area of one unit cell. This extra $\mathrm{\Omega }$ factor in the denominator on the right-hand side has to do with the fact that our $\mathit{q}$ is not continuous, but discretized on a grid and our Fourier transform is the discrete version. Putting it in another way, $1/\mathrm{\Omega }$ is the value of the discretized two-dimensional Dirac $\delta$-function in the transverse direction, Fourier transformed in $\mathit{q}$ space. In the space of $\mathit{q}$, and $z$ taking 0 and $d$, the matrices are $2\times 2$, as

(79)$\begin{array}{l}(D^r{)}^{-1}={\left[{\displaystyle \frac{1}{2{ϵ}_{0}q\mathrm{\Omega }}}(\begin{array}{cc}1& {\mathrm{e}}^{-qd}\\ {\mathrm{e}}^{-qd}& 1\end{array})\right]}^{-1}-(\begin{array}{cc}{\mathrm{\Pi }}_L& 0\\ 0& {\mathrm{\Pi }}_R\end{array}).\end{array}$

For the first term as $v^{-1}=C$, the inverse can be worked out to give $C_{11}=C_{22}=2{ϵ}_{0}q\mathrm{\Omega }/(1-{\mathrm{e}}^{-2qd})$, and $C_{12}=C_{21}=-{\mathrm{e}}^{-qd}{C}_{11}$. We note that if we take the limit $q\to 0$, we obtain the same matrix as before with a capacitance in this case as ${ϵ}_0\mathrm{\Omega }/d$, i.e., the effective area of the capacitor is the area of a unit cell. The explicit expression for $D^r$ matrix elements can be easily worked out, e.g.,

(80)$\begin{array}{l}{D}_{LR}^r=D_{RL}^r={\displaystyle \frac{-C_{21}}{(C_{11}-{\mathrm{\Pi }}_L)(C_{22}-{\mathrm{\Pi }}_R)-C_{12}{C}_{21}}}.\end{array}$

The advanced version is obtained by Hermitian conjugate. We define the reflection coefficient as

(81)$\begin{array}{l}{r}_{\alpha }={\displaystyle \frac{v_{\alpha }{\mathrm{\Pi }}_{\alpha }}{1-v_{\alpha }{\mathrm{\Pi }}_{\alpha }}}=v_{\alpha }{\chi }_{\alpha },\alpha =L,R.\end{array}$

Here, $v_{\alpha }=1/(2{ϵ}_{0}q\mathrm{\Omega })$ is the Coulomb potential in two-dimension in $\mathit{q}$ space, ${\chi }_{\alpha }$ is charge-charge correlation or susceptibility, while ${\mathrm{\Pi }}_{\alpha }$ is the self-energy or polarizability. ${ϵ}_{\alpha }=1-v_{\alpha }{\mathrm{\Pi }}_{\alpha }$ is the dielectric function. Using the reflection coefficients with some algebra, we can simplify the transmission as [108]

(82)$\begin{array}{l}{T}_{\mathit{q}}(\omega )={\displaystyle \frac{4\mathrm{I}\mathrm{m}{r}_L\mathrm{I}\mathrm{m}{r}_R{\mathrm{e}}^{-2qd}}{{\left|1-r_L{r}_R{\mathrm{e}}^{-2qd}\right|}^2}}.\end{array}$

Of course, $r_{\alpha }$ is a function of the wave-vector $\mathit{q}$ as well as the angular frequency $\omega$ which we have suppressed. One might be curious about why $r_{\alpha }$ is called a reflection coefficient. Indeed, it does have to do with the wave reflection. In the traditional approach to near-field heat transfer of Polder and van Hove [6,109], one solves a wave scattering problem with transmission and multiple reflections between the plates. Ignoring the $s$-wave polarization, which is small at near distance, and focusing on the $p$-wave (electric field is in the plane of incidence) and at the non-retardation limit of $c\to \mathrm{\infty }$, one obtains exactly the same result as above using the fluctuational electrodynamics. Here our approach is based on NEGF, which in some way is simpler.

Finally, to obtain the total heat current, one sums over all the modes $\mathit{q}$ in the first Brillouin zone with $N$ sampling points, and integrates over the frequency, i.e., $I_L={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\mathrm{d}\omega /(4\pi )\hslash \omega \sum _{\mathit{q}}{T}_{\mathit{q}}(\omega )({\tilde{N}}_L-{\tilde{N}}_R)$. The number of $\mathit{q}$ points is related to the actual area of the plates by $A=N\mathrm{\Omega }$. To obtain the heat transfer per unit area, we divide by the area, i.e., $(1/A)\sum _{\mathit{q}}\cdots =\int {\mathrm{d}}^2\mathit{q}/(2\pi {)}^2\cdots$. It is important that the points of $\mathit{q}$ are in the first Brillouin zone, i.e., the Wigner−Seitz cell with the $\mathrm{\Gamma }$-point at the cell center. This is because although the materials have lattice periodicity with ${\mathrm{\Pi }}_{\mathit{q}}$ a periodic function in the reciprocal space, the photon free Green’s function $v$ does not has such periodicity, our forcing $D^r$ running on lattice sites is an approximation.

After this long preparation, in Fig.6, we illustrate the effect of drift velocity to heat transfer in the two-layer graphene setup. We emphasize that the current drive not only changes the Bose distribution by a Doppler shift, but also, the self-energy $\mathrm{\Pi }$ or in turn the reflection coefficient $r$ needs to be recalculated anew. The changes in ${\mathrm{\Pi }}_{\alpha }$ or $r_{\alpha }$ are necessary so that the integrand for the heat transfer is not divergent or at least still integrable at the point when $\omega -\mathit{q}\cdot {\mathit{v}}_1=0$. Details are in [89]. Here we drive with velocity $v_1$ in $x$ direction for the left layer (called 1 in the original paper), and compute the heat transfer out of the left layer. Although the drive breaks symmetry in frequency for each given $\mathit{q}$, after integration over all $\mathit{q}$, the even symmetry in $\omega$ is restored. Fig.6(a) plots the $\mathit{q}$ integrated result $f(\omega )$ for fixed $\omega$. Further integration over frequency gives the total heat transfer. We notice originally without drift, heat is flowing from right to left (negative value on the plot) since the right is hotter. But as the drift velocity increases, the heat transfer reverses sign, going from cold to hot. This is understandable as the left layer of graphene is no longer in local thermal equilibrium, and it is not a broken down of the second law of thermodynamics. Fig.6(b) demonstrates the effect of drift to the distribution of the total integrand over $q_x$ integrated over $q_y$ and $\omega$. Without drift, the distribution in $q_x$ is symmetric with respect to $q_x=0$, while drift breaks this symmetry, and it turns into having both positive and negative contributions when $v_1$ is large, causing a cancellation effect for the total heat transfer. In Fig.7, the reversal of the heat transfer direction is clearly shown when the temperature of the right-side sheet is varied. There is a particular balance temperature where heat transfer is zero even though the temperatures of the two sheets are different. Above this temperature, we have a net cooling effect for the right side due to the current driven on the left side. We comment that driving a conductor with current will produce Joule heat as well as electro-luminescence, the effects of electron−phonon and electron-transverse-photon interactions. This extra heat has not been taken into account in our theory.

Fig.6 (a) Integrated spectral transfer function $f(\omega )$ as a function of frequency and (b) $g(q_x)$ as a function of wave vector in the driven direction, with different drift velocities: no drift (blue dash-dot line), total heat current density $I_L/A=-0.84$ MW/m2; $v_1=5.0\times {10}^5\text{m/s}$ (red dash line), $-0.30$ MW/m2; and $v_1=9.0\times {10}^5\text{m/s}$ (black solid line), $+0.09$ MW/m2. The temperatures are $T_L=300$ K and $T_R=320$ K. The chemical potential of graphene $\mu$ is set as 0.1 eV. Gap distance $d$ is set as 10 nm. The damping parameter is $\eta =9$ meV. Reproduced from Ref. [89].

Full size|PPT slide

Fig.7 Heat current density from left to right layer as a function of $T_R$ (temperature of the right layer) with drift velocity $v_1=9.0\times {10}^5\text{m/s}$ (red solid line). The dotted line is a reference line for zero current density. The green circle indicates the point for “off temperature”. The chemical potential of graphene $\mu$ is set as 0.1 eV. Temperature of left layer of graphene $T_L$ is set as 300 K. Vacuum gap distance $d$ is set as 10 nm. Reproduced from Ref. [89].

Full size|PPT slide

6 Full counting statistics for energy transfer

At the nanoscale due to thermal agitation, the measured results themselves are fluctuational quantities [110,111]. The full counting statistics [112,113] here means that we compute not only energy but also high order moments in a transport setup. To compute the higher moments, it is convenient we compute the total heat $Q$ of a fixed duration in a two-time measurement [114] of the left bath ${\hat{H}}_L$. The energy of the left bath is measured at an earlier time $t_0$ and then measured again at time $t$, the decrease in energy is the transfer of the heat to the right out of the left bath. According to the standard measurement interpretation of quantum mechanics, the result of a measurement is an eigenvalue of ${\hat{H}}_L$. We will be interested in the long time $(t-t_0)\to \mathrm{\infty }$ which gives a simpler result. Other protocols of measurements are also possible, but the two-time measurement results in a simpler mathematics, although it is not clear how this measurement of energy of the bath can be done experimentally when the bath is supposed to be infinitely large. For the long-time result, the details of measurement protocols should not matter.

We consider the system as two blocks of metal with a Hamiltonian ${\hat{H}}_L+{\hat{H}}_R+{\hat{H}}_{\gamma }+\hat{V}={\hat{H}}_0+\hat{V}$. ${\hat{H}}_{\gamma }$ is the negative-definite free scalar photon Hamiltonian. The last term $\hat{V}$ is the Coulomb interaction in the scalar field form, $\hat{V}=\sum _{jkl}{M}_{jk}^l{c}_j^{†}{c}_k{\varphi }_l$. Only parts of the sites near the interface between the left and right blocks have the Coulomb interaction. Deep into the baths, we set the Coulomb interaction to zero (due to screening). We prepare the system to be in the decoupled product initial state given by the density matrix

(83)$\begin{array}{l}\hat{\rho }(t_0)\propto {\mathrm{e}}^{-{\beta }_L({\hat{H}}_L-{\mu }_L{\hat{N}}_L)}{\mathrm{e}}^{-{\beta }_R({\hat{H}}_R-{\mu }_R{\hat{N}}_R)}{\mathrm{e}}^{-{\beta }_{\gamma }{\hat{H}}_{\gamma }},\end{array}$

here the left and right baths are in the grand-canonical ensembles and the last factor is for the scalar photons in canonical ensemble. We can set ${\beta }_{\gamma }=0$ without loss of generality.

We can prove a very general formula by defining a “partition function” or moment generating function as

(84)$\begin{array}{l}Z(\xi )=⟨T_{\tau }{\mathrm{e}}^{-\frac{\mathrm{i}}{\hslash }\int \mathrm{d}\tau {\hat{V}}^x(\tau )}{⟩}_{{\hat{H}}_0}.\end{array}$

Here, the exponential is contour ordered, and the integral is over the contour from $t_0$ on the upper branch to $t$ and then back to $t_0$ from the lower branch. ${\hat{V}}^x$ is the interaction term, but Heisenberg evolved by $\hslash x$ in the interaction picture by the left bath, i.e.,

(85)$\begin{array}{l}{\hat{V}}^x={\mathrm{e}}^{\mathrm{i}x{\hat{H}}_L}\hat{V}{\mathrm{e}}^{-\mathrm{i}x{\hat{H}}_L}.\end{array}$

The amount of $x$ is contour dependent, it is $-\xi /2$ on the upper $(+)$ branch and $\xi /2$ on the lower $(-)$ branch. With this definition of $Z$, the $n$-th moment of $Q$ is obtained by $n$-th derivative of $Z$, $⟨Q^n⟩={\mathrm{\partial }}^{n}Z(\xi )/\mathrm{\partial }(\mathrm{i}\xi {)}^n$ evaluated at $\xi =0$ and the $n$-th order cumulant is

(86)$\begin{array}{l}⟨⟨Q^n⟩⟩={\displaystyle \frac{{\mathrm{\partial }}^n\ln Z(\xi )}{\mathrm{\partial }(\mathrm{i}\xi {)}^n}}{|}_{\xi =0}.\end{array}$

The advantage of working with the cumulants instead of moments is that they are all linear in $t-t_0$ at long time. The first order moment and cumulant are the same, $⟨Q⟩=⟨⟨Q⟩⟩=\mathrm{\partial }\ln Z/\mathrm{\partial }(\mathrm{i}\xi )=(t-t_0)I_L$, $I_L$ is the current. The second cumulant is just the variance or the fluctuations of the current, $⟨⟨Q^2⟩⟩=⟨Q^2⟩-⟨Q{⟩}^2=(t-$ $t_0){\sigma }^2$.

To show the validity of Eq. (84), let ${\hat{H}}_L|{\phi }_a⟩=a|{\phi }_a⟩$, here $a$ is the eigenvalue of the isolated left side with the eigenstate $|{\phi }_a⟩$. If a measurement is performed at time $t_0$ obtaining the eigenvalue $a$, and then measured again at time $t$ obtaining the eigenvalue $b$, the generating function is

(87)$\begin{array}{rl}Z(\xi )& =\sum _{a,b}{\mathrm{e}}^{\mathrm{i}(a-b)\xi }P(b,a)\\ & =\sum _{a,b}{\mathrm{e}}^{\mathrm{i}(a-b)\xi }\mathrm{T}\mathrm{r}[\hat{\rho }(t_0){\mathrm{\Pi }}_{a}U(t_0,t){\mathrm{\Pi }}_{b}U(t,t_0)].\end{array}$

Here $P(b,a)$ is the probability of being in state $b$ at time $t$ given that it is in state $a$ at an initial time $t_0$, which can be expressed by the evolution operator $U$ of the full Hamiltonian and the projectors of the respective states,

(88)$\begin{array}{l}{\mathrm{\Pi }}_a=|{\phi }_a⟩⟨{\phi }_a|.\end{array}$

We see that $Z(0)=1$ due to the probability normalization. Taking the derivative with respect to $(i\xi )$ $n$ times, and then set $\xi$ to 0, we obtain the expectation value of $Q^n=(a-b{)}^n$ over the probability $P(b,a)$. we have $[{\mathrm{\Pi }}_a,{\hat{H}}_L]=0$, and because of the choice of the product initial state, we also have $[{\mathrm{\Pi }}_a,\hat{\rho }(t_0)]=0$. Since

(89)$\begin{array}{l}{\displaystyle \sum _a{\mathrm{e}}^{\mathrm{i}a\xi }{\mathrm{\Pi }}_a={\mathrm{e}}^{\mathrm{i}\xi {\hat{H}}_L},{\displaystyle \sum _b{\mathrm{e}}^{-\mathrm{i}b\xi }{\mathrm{\Pi }}_b={\mathrm{e}}^{-\mathrm{i}\xi {\hat{H}}_L},}}\end{array}$

we can express $Z(\xi )$ as

(90)$\begin{array}{rl}Z(\xi )& =\mathrm{T}\mathrm{r}[\hat{\rho }(t_0){\mathrm{e}}^{\mathrm{i}\xi {\hat{H}}_L}U(t_0,t){\mathrm{e}}^{-\mathrm{i}\xi {\hat{H}}_L}U(t,t_0)]\\ & =\mathrm{T}\mathrm{r}[\hat{\rho }(t_0)U^{\frac{\xi }{2}}(t_0,t)U^{-\frac{\xi }{2}}(t,t_0)].\end{array}$

Here we have split the exponential factors into two halves and used a cyclic permutation of the trace. The superscript on $U$ denotes an extra Heisenberg evolution with ${\hat{H}}_L$, i.e., $U^x={\mathrm{e}}^{\mathrm{i}x{\hat{H}}_L}U{\mathrm{e}}^{-\mathrm{i}x{\hat{H}}_L}$. This extra $x$ dependence can be transferred from $U$ into the Hamiltonian, to give ${\hat{H}}^x={\hat{H}}_0+{\hat{V}}^x$. Here the noninteracting Hamiltonian ${\hat{H}}_0$ is unaffected as ${\hat{H}}_L$ commutes with the three free terms, ${\hat{H}}_L$, ${\hat{H}}_R$, and ${\hat{H}}_{\gamma }$. At this point, we transform the expression into the interaction picture, given Eq. (84) with the average evaluated with respect to ${\hat{H}}_0$. In the interaction picture, the effect of the Heisenberg evolution is to shift the time argument of the left side as $c_L(\tau )=c_L(t-\hslash \xi /2)$ on the upper branch, and $c_L(t+\hslash \xi /2)$ on the lower branch. This in turn means a corresponding shifts of time argument for the self-energy ${\mathrm{\Pi }}_L$.

The rest of the steps follow a standard diagrammatic expansion. We work at the level of RPA for the Coulomb interaction. Since the photon Hamiltonian is quadratic in the scalar field $\varphi$, an odd number of $\hat{V}$ evaluates to zero. Thus the lowest order of nonzero term, when the exponential is expanded, is $1/[2(\mathrm{i}\hslash {)}^2]\cdot$$\int \mathrm{d}{\tau }_1\int \mathrm{d}{\tau }_2⟨T_{\tau }\hat{V}({\tau }_1)\hat{V}({\tau }_2){⟩}_0$. Applying Wick’s theorem, we can write the result as $\frac{1}{2}{\mathrm{T}\mathrm{r}}_{\tau j}(v{\mathrm{\Pi }}^{(0)})$, here $v$ is the bare Coulomb in $j$ and $\tau$ space which contains a $\delta$ function in the $\tau$ variable as Coulomb interaction is instantaneous, and ${\mathrm{\Pi }}^{(0)}$ is the RPA bubble diagram result. We will drop the superscript $(0)$ for notational simplicity. If we continue the expansion to higher orders, and collect only these bubble diagrams that form a ring (see Fig.8), we can sum the diagrams as a logarithm due to the $1/n$ factor at the $n$-th order. These diagrams are the same as the RPA result for the grand potential in a Matsubara Green’s function formulation [94]. Finally, we have the RPA expression for the generating function as

Fig.8 The first few terms of diagrammatic expansion of $\mathrm{l}\mathrm{n}Z$.

Full size|PPT slide

(91)$\begin{array}{l}\ln Z(\xi )=-{\displaystyle \frac{1}{2}}\mathrm{T}\mathrm{r}\ln [I-v({\mathrm{\Pi }}_L^x+{\mathrm{\Pi }}_R)].\end{array}$

Here the trace is performed in the combined contour time $\tau$ and site $j$ space, that is, $\mathrm{T}\mathrm{r}(\cdot )=\sum _j\int \mathrm{d}\tau (\cdot )$. $I$ is the identity is this space. To aid the Fourier transform in the full time domain, it is convenient to take the limit $t_0\to -\mathrm{\infty }$. We eliminate the bare Coulomb $v$ in favor of the screened Coulomb $D$ by introducing ${\mathrm{\Pi }}_L^A={\mathrm{\Pi }}_L^x-{\mathrm{\Pi }}_L$, then $I-v({\mathrm{\Pi }}_L^x+{\mathrm{\Pi }}_R)=v(v^{-1}-{\mathrm{\Pi }}_L-{\mathrm{\Pi }}_R)-v{\mathrm{\Pi }}_L^A=v(D^{-1}-$ ${\mathrm{\Pi }}_L^A)=v{D}^{-1}(I-D{\mathrm{\Pi }}_L^A)$. Since the factor $v{D}^{-1}$ does not depend on $\xi$, it is an additive constant to $\ln Z$ and will not contribute to the derivatives, so we drop it and redefine the generating function as

(92)$\begin{array}{l}\ln Z(\xi )=-{\displaystyle \frac{1}{2}}\mathrm{T}\mathrm{r}\ln \left(I-D{\mathrm{\Pi }}_L^A\right).\end{array}$

If we Taylor expand in $\xi$, the linear term gives the current $I_L$. After simplifying from the contour time to real time, and then Fourier transform to frequency domain for the Green’s functions and self-energies, we recover the Meir−Wingreen formula, Eq. (35). If we expand up to second order in $\xi$, we obtain the variance. Explicit formula for the variance of heat transfer has been derived by Herz et al. [110] within a scattering operator formalism. If we use local equilibrium, we can express the generating function in terms of a matrix $D^r{\mathrm{\Gamma }}_R{D}^a{\mathrm{\Gamma }}_L$ which is identical to the usual Levitov−Lesovik formula [113,115,116].

7 Density functional theory calculation based on Coulomb interaction

So far, our electron models have been in the tight-binding form where the electrons are allowed to sit only at a discrete set of sites. This has the computational efficiency advantage of dealing with finite-dimensional Hermitian matrices, enough to capture the solid-state band structure for a lattice, for example. In fact, our systems can be any structure unrestricted by lattice periodicity, e.g., tips of two triangles made of graphene [117]. Assuming that electrons sit only at certain discrete places is an approximation. The density functional theory (DFT) [118] offers a parameter-free approach to real complex materials. In DFT approach the electrons are treated as distributed in the whole space continuously. Another feature of the plane-wave based DFT is that we must work on a periodic unit cell or super-cell. As a result, we must consider the fluctuations of the electron density within a cell. The theory developed earlier based on the Green’s function $D$ and the self-energy $\mathrm{\Pi }$ needs some revisions, but it is only of a technical nature and no new conceptual difficulty remains.

The definition $D$ remains the same as given by Eq. (45), whether the electrons are treated as discrete or continuous degrees of freedom. For a continuum of electron density, $\rho$, it is convenient to define the contour version $\mathrm{\Pi }$ as

(93)$\begin{array}{l}\mathrm{\Pi }(\mathit{r}\tau ;{\mathit{r}}^{\prime }{\tau }^{\prime })={\displaystyle \frac{1}{\mathrm{i}\hslash }}⟨T_{\tau }\rho (\mathit{r},\tau )\rho ({\mathit{r}}^{\prime },{\tau }^{\prime }){⟩}_{\mathrm{i}\mathrm{r}}.\end{array}$

Here the subscript $\mathrm{i}\mathrm{r}$ means that we take only the irreducible diagrams in a Feynman-diagrammatic expansion with the Coulomb interaction. The irreducible ones are those that cannot be cut into two disconnected pieces by a single bare scalar photon line $v$. Without it, the charge-charge correlation is $\chi =\mathrm{\Pi }{ϵ}^{-1}$; $ϵ=1-v\mathrm{\Pi }$ is the longitudinal dielectric function. Under the random phase approximation, we just take the lowest order bubble diagram. We can lump the charges in a cell to a point. Assuming all the cells having the same volume $\mathrm{\Omega }$, the relation between a continuous description and the earlier discrete description for $\mathrm{\Pi }$ is ${\mathrm{\Pi }}_{ij}={\mathrm{\Omega }}^2\mathrm{\Pi }({\mathit{r}}_i,{\mathit{r}}_j)$. Here ${\mathit{r}}_i$ is a point in the cell $i$. Exactly where it is in the cell does not matter provided that the function varies with the position smoothly enough. From this relation, we can see that the Dyson equation in a continuous charge description is obtained by replacing the summations by integrals over space, i.e., in Eq. (47), we replace $\sum _{j,k}\int \mathrm{d}{\tau }_1\int \mathrm{d}{\tau }_{2}v(\mathit{r}\tau ,{\mathit{r}}_j{\tau }_1){\mathrm{\Pi }}_{jk}({\tau }_1,{\tau }_2)\cdot$ $D({\mathit{r}}_k{\tau }_2,{\mathit{r}}^{\prime }{\tau }^{\prime })$ by the integral

$\int {\mathrm{d}}^3{\mathit{r}}_1\int \mathrm{d}{\tau }_1\int {\mathrm{d}}^3{\mathit{r}}_2\int \mathrm{d}{\tau }_{2}v(\mathit{r}\tau ,{\mathit{r}}_1{\tau }_1)\mathrm{\Pi }({\mathit{r}}_1{\tau }_1,{\mathit{r}}_2{\tau }_2)D({\mathit{r}}_2{\tau }_2,{\mathit{r}}^{\prime }{\tau }^{\prime }).$

The quantities like the above expressed in real space $\mathit{r}$ are not convenient for actual computation. The translation invariance in a crystal means we should evoke the convolution theorem in Fourier transform so that it becomes multiplication in $\mathit{q}$ space. However, if we also take into account the local inhomogeneity inside a cell, the Fourier transform is slightly more complicated for which we elaborate below.

For a periodic system, the single-electron wave function satisfies the Bloch theorem [119], $\varphi (\mathit{r})={\mathrm{e}}^{\mathrm{i}\mathit{k}\cdot \mathit{r}}u(\mathit{r})$; that is, a plane wave modulated by a periodic function, $u(\mathit{r})=u(\mathit{r}+\mathit{R})$. Here $\mathit{R}=l_1{\mathit{a}}_1+l_2{\mathit{a}}_2+l_3{\mathit{a}}_3$ specifies the Bravais lattice sites by the unit cell vectors ${\mathit{a}}_i$ and integers $l_i$. This fact implies that the average electron density is a cell-periodic function,

(94)$\begin{array}{l}⟨\rho (\mathit{r})⟩=⟨\rho (\mathit{r}+\mathit{R})⟩,\end{array}$

since under the Kohn-Sham DFT framework, the electron density is simply the sum of $|\varphi (\mathit{r}){|}^2$ of the occupied bands. The operator $\rho$ itself before the thermodynamic average does not have this cell periodicity. A periodic function can be Fourier expanded as

(95)$\begin{array}{l}⟨\rho (\mathit{r})⟩={\displaystyle \sum _{\mathit{G}}{\rho }_{\mathit{G}}{\mathrm{e}}^{\mathrm{i}\mathit{G}\cdot \mathit{r}},}\end{array}$

here the capital letter $\mathit{G}$ is the reciprocal lattice vector, running over all the values $\mathit{G}=n_1{\mathit{b}}_1+n_2{\mathit{b}}_2+n_3{\mathit{b}}_3$, $n_i$ are integers and ${\mathit{a}}_i\cdot {\mathit{b}}_j=2\pi {\delta }_{ij}$, $i,j=1,2,3$. The Fourier coefficients are obtained by ${\rho }_{\mathit{G}}=(1/\mathrm{\Omega }){\int }_{\mathrm{\Omega }}⟨\rho (\mathit{r})⟩{\mathrm{e}}^{-\mathrm{i}\mathit{G}\cdot \mathit{r}}{\mathrm{d}}^3\mathit{r}$, integrating over one unit cell with cell volume $\mathrm{\Omega }$.

The charge-charge correlation involves two positions, possibly in different unit cells. For simplicity, let us consider the static correlation $f(\mathit{r},{\mathit{r}}^{\prime })=⟨\rho (\mathit{r})\rho ({\mathit{r}}^{\prime })⟩$. The lattice cell translation means that, if we shift both positions by a common Bravais lattice vector, the correlation between the new pair should be the same, i.e., we have

(96)$\begin{array}{l}f(\mathit{r}+\mathit{R},{\mathit{r}}^{\prime }+\mathit{R})=f(\mathit{r},{\mathit{r}}^{\prime }).\end{array}$

Note that this is different from a continuous translation symmetry of $f$ being a function of the difference $\mathit{r}-{\mathit{r}}^{\prime }$ only. We will show that the cell translation symmetry can be Fourier expanded with double series $\mathit{G}$ and ${\mathit{G}}^{\prime }$ and a Fourier integral in the first Brillouin zone, as

(97)$\begin{array}{l}f(\mathit{r},{\mathit{r}}^{\prime })={\displaystyle \sum _{\mathit{G},{\mathit{G}}^{\prime }}{\int }_{1\mathrm{B}\mathrm{Z}}\frac{{\mathrm{d}}^3\mathit{q}}{(2\pi {)}^3}{\tilde{F}}_{\mathit{G}{\mathit{G}}^{\prime }}(\mathit{q}){\mathrm{e}}^{\mathrm{i}(\mathit{G}+\mathit{q})\cdot \mathit{r}-\mathrm{i}({\mathit{G}}^{\prime }+\mathit{q})\cdot {\mathit{r}}^{\prime }}.}\end{array}$

In the above, if we keep only the $\mathit{G}={\mathit{G}}^{\prime }=\mathbf{0}$ term, we obtain a continuous translation symmetry, $f(\mathit{r},{\mathit{r}}^{\prime })=$ $f(\mathit{r}-{\mathit{r}}^{\prime })$; this is a long-wave approximation. The extra nonzero $\mathit{G}$-vector terms reflect the local inhomogeneity within a cell. To proof this result, first we can choose a Bravais vector such that the first argument is in the first cell, and define $F(\mathit{r},\mathbf{\Delta })=f(\mathit{r},\mathit{r}+\mathbf{\Delta })$. Here $\mathbf{\Delta }$ is still arbitrary running over the whole space. If we consider $F$ as a function of $\mathit{r}$ fixing the difference $\mathbf{\Delta }={\mathit{r}}^{\prime }-\mathit{r}$, $F$ is a periodic function of $\mathit{r}$. So, we can write

(98)$\begin{array}{l}F(\mathit{r},\mathbf{\Delta })={\displaystyle \sum _{\mathit{G}}{c}_{\mathit{G}}(\mathbf{\Delta }){\mathrm{e}}^{\mathrm{i}\mathit{G}\cdot \mathit{r}}.}\end{array}$

The Fourier coefficients $c_{\mathit{G}}$ are still a function of $\mathbf{\Delta }$ which varies continuously in the full space. For the variable $\mathbf{\Delta }$ we can make a Fourier integral transform into $\mathit{q}$. This is

(99)$\begin{array}{l}{c}_{\mathit{G}}(\mathbf{\Delta })={\displaystyle \int \frac{{\mathrm{d}}^3\mathit{q}}{(2\pi {)}^3}{\tilde{f}}_{\mathit{G}}(\mathit{q}){\mathrm{e}}^{-\mathrm{i}\mathit{q}.\mathbf{\Delta }}.}\end{array}$

The integral over the full Fourier space $\mathit{q}$ can be split into pieces of reciprocal space cell ${\mathit{G}}^{\prime }$ with a change of variable to each cell by $\mathit{q}\to {\mathit{G}}^{\prime }+\mathit{q}$. Then the new $\mathit{q}$ variable varies in the first Brillouin zone only. After some regrouping and simplification, noting that $\mathbf{\Delta }={\mathit{r}}^{\prime }-\mathit{r}$, and defining ${\tilde{F}}_{\mathit{G}{\mathit{G}}^{\prime }}(\mathit{q})={\tilde{f}}_{\mathit{G}-{\mathit{G}}^{\prime }}({\mathit{G}}^{\prime }+\mathit{q})$, we obtain the desired result, Eq. (97). Tracing back the steps, we can compute the Fourier expansion coefficient as

(100)$\begin{array}{rl}& {\tilde{F}}_{\mathit{G}{\mathit{G}}^{\prime }}(\mathit{q})=\sum _{\mathit{R}}\frac{1}{\mathrm{\Omega }}{\int }_{\mathrm{\Omega }}{\mathrm{d}}^3\mathit{r}{\int }_{\mathrm{\Omega }}{\mathrm{d}}^3{\mathit{r}}^{\prime }f(\mathit{r}+\mathit{R},{\mathit{r}}^{\prime }){\mathrm{e}}^{-\mathrm{i}\phi },\\ & \phi =(\mathit{G}+\mathit{q})\cdot \mathit{r}-({\mathit{G}}^{\prime }+\mathit{q})\cdot {\mathit{r}}^{\prime }+\mathit{q}\cdot \mathit{R}.\end{array}$

Here both of the integral variables $\mathit{r}$ and ${\mathit{r}}^{\prime }$ are in the first unit cell. With the $\mathit{G}$, ${\mathit{G}}^{\prime }$, and $\mathit{q}$ variables for the correlation functions, we also have a convolution theorem. That is, the expression of type $v\mathrm{\Pi }D$ in the Dyson equation can be written as a matrix multiplication indexed by $\mathit{G},{\mathit{G}}^{\prime }$ for each given $\mathit{q}$. Similarly, a trace by integration over position $\mathit{r}$ can now be expressed as a trace in $\mathit{G}$ as a matrix trace and integration of $\mathit{q}$ in the first Brillouin zone.

7.1 Adler−Wiser formula

We now give a formula for the retarded ${\mathrm{\Pi }}^r$ expressed by the Kohn−Sham or independent single-particle orbitals known as the Adler−Wiser formula [120,121]. The retarded formula can be computed according to ${\mathrm{\Pi }}^r(t)=\theta (t)({\mathrm{\Pi }}^{>}(t)-{\mathrm{\Pi }}^{<}(t))$, and then Fourier transformed into frequency domain. The lesser and greater components can be read-off from the contour expression, Eq. (93), as ${\mathrm{\Pi }}^{>}(t)=⟨\rho (t)\rho (0)⟩/(\mathrm{i}\hslash )$ and ${\mathrm{\Pi }}^{<}(t)=⟨\rho (0)\rho (t)⟩/(\mathrm{i}\hslash )$. We remind the reader here that $\rho$ is a quantum operator, which can be expressed in the quantum field as $\rho =(-e){\mathrm{\Psi }}^{†}\mathrm{\Psi }$, where $\mathrm{\Psi }$ is space and time dependent which we expand in the mode space,

(101)$\begin{array}{l}\mathrm{\Psi }(\mathit{r}t)={\displaystyle \sum _{nk\sigma }{c}_{n\mathit{k}\sigma }(t){\varphi }_{n\mathit{k}\sigma }(\mathit{r}).}\end{array}$

Here the eigenmodes are labeled by the band index $n$, the wave vector $\mathit{k}$, and spin $\sigma =\uparrow ,\downarrow$. However, for a spin-independent problem, the net effect of spin degeneracy is simply multiplying the final expression of the polarizability by a factor of 2. In the following, we will treat our electrons as spinless and then keep a factor of 2 for ${\mathrm{\Pi }}^r$. The Kohn-Sham wave function ${\varphi }_{n\mathit{k}}$ must be normalized to 1 in the whole system of volume $V=N\mathrm{\Omega }$ in order to give the correct electron density. In mode space, the Hamiltonian is diagonal, thus the time-dependence for the annihilation operator is simply the free evolution, $c_{n\mathit{k}}(t)=c_{n\mathit{k}}{\mathrm{e}}^{-\mathrm{i}{ϵ}_{n\mathit{k}}t/\hslash }$, here ${ϵ}_{n\mathit{k}}$ is the electron band energy.

In evaluating the density-density correlation, we encounter terms of the form $⟨c^{†}c{c}^{†}c⟩$ which we apply Wick’s theorem to factor into product of two $c$’s. Noting that $⟨cc⟩$ or $⟨c^{†}{c}^{†}⟩$ is 0, the remaining terms are related to the Fermi function, i.e., $⟨c_j^{†}{c}_l⟩={\delta }_{jl}{f}_j$, and $⟨c_j{c}_l^{†}⟩=$ ${\delta }_{jl}(1-f_j)$. Here we have used a short-hand notation $j\equiv (n\mathit{k})$, $l\equiv (n^{\prime }{\mathit{k}}^{\prime })$. With some algebra, we can express the retarded scalar photon self-energy in the frequency domain as

(102)$\begin{array}{l}{\mathrm{\Pi }}^r(\mathit{r},{\mathit{r}}^{\prime },\omega )=2e^2{\displaystyle \sum _{jl}\frac{(f_j-f_l){\varphi }_j(\mathit{r}){\varphi }_j^{\ast }({\mathit{r}}^{\prime }){\varphi }_l({\mathit{r}}^{\prime }){\varphi }_l^{\ast }(\mathit{r})}{{ϵ}_j-{ϵ}_l-\hslash \omega -\mathrm{i}\eta }.}\end{array}$

Here the extra small damping $\mathrm{i}\eta$ in the denominator is necessary, so that the poles in a complex frequency plane is below the real axis, and the inverse transform has the property ${\mathrm{\Pi }}^r(t)=0$ when $t<0$. This $\eta$ parameter also has a physical meaning. We can interpret it as the inverse of electron lifetime of the quasi particle. Finally, using our general transformation formula for the $(\mathit{r},{\mathit{r}}^{\prime })$ correlation to the $\mathit{G},{\mathit{G}}^{\prime },\mathit{q}$ space, the self-energy becomes

(103)$\begin{array}{rl}& {\mathrm{\Pi }}_{\mathit{G}{\mathit{G}}^{\prime }}^r(\mathit{q},\omega )\\ & =\frac{2e^2}{N\mathrm{\Omega }}\sum _{j,l}\frac{⟨l|{\mathrm{e}}^{-\mathrm{i}(\mathit{G}+\mathit{q})\cdot \hat{\mathit{r}}}|j⟩⟨j|{\mathrm{e}}^{\mathrm{i}({\mathit{G}}^{\prime }+\mathit{q})\cdot \hat{\mathit{r}}}|l⟩(f_j-f_l)}{{ϵ}_j-{ϵ}_l-\hslash \omega -\mathrm{i}\eta }.\end{array}$

Here we define the matrix element $⟨l|\hat{A}|j⟩\equiv$ ${\int }_{N\mathrm{\Omega }}{\mathrm{d}}^3\mathit{r}{\varphi }_l^{\ast }\hat{A}{\varphi }_j$; the integration is over the full volume with $N$ unit cells. This is the Adler−Wiser formula. Note that the single particle operator $\hat{A}$ is a shift of the momentum, so the matrix elements are zero unless the momentum $\mathit{k}$ in $j$ is related to that in $l$ by $\mathit{k}={\mathit{k}}^{\prime }+\mathit{q}$ modulo $\mathit{G}$. Efficient evaluation of this expression with massively parallel algorithms has been implemented in the well-known BerkeleyGW package [122,123], but only for zero temperature.

7.2 Solving the Dyson equation

Here we consider the parallel plate geometry for near-field heat transfer. Since the first-principles codes require periodic supercell while the heat transfer problem is intrinsically a non-periodic problem in the transport direction, there is a fundamental conflict. We can imagine putting two slabs of materials into the simulation cell with sufficient vacuum gaps. Although the transverse directions are intrinsically periodic, the transfer direction, if we still use periodic boundary condition, may requires large vacuum gap to void artificial interactions. As a result, it is best we treat each slab separately by DFT and combine the results [124]. We also assume the periodicity in the transverse directions for the two slabs is the same, otherwise, how to combine them is a problem. Finally, if the slab is thick, we work in a mixed representation of ${\mathit{G}}_{\perp }$ and $z$, and consider the self-energy in the form ${\mathrm{\Pi }}_{{\mathit{G}}_{\perp }{\mathit{G}}_{\perp }^{\prime }}^r({\mathit{q}}_{\perp },z,z^{\prime },\omega )$. Here the $z$ dependence is obtained by Fourier transforming $G_z,G_z^{\prime },q_z$ into real space, and ${\mathit{G}}_{\perp }=(G_x,G_y)$ lays in the plane.

Very likely the above approach is still too computationally intensive; so far no one was able to do a calculation for near-field heat transfer. In the following, we make further approximation, this is to treat the slab as infinitely thin, in such a way that the density of the electrons is confined strictly in 2D at location $z=0$ and $d$. In this way, the $z$ variables become discrete. We can define surface charge density by integrating over $z$ of the volume density, $\sigma =\int \rho \mathrm{d}z$, and define the surface version of $\mathrm{\Pi }$. A careful analysis shows that this surface version can be obtained from the volume version in $\mathit{q}$ space just by multiplying the supercell length $c$ in $z$ direction, i.e., [125]

(104)$\begin{array}{l}{\mathrm{\Pi }}_{{\mathit{G}}_{\perp }{\mathit{G}}_{\perp }^{\prime }}^{2\mathrm{D}}({\mathit{q}}_{\perp },\omega )=c{\mathrm{\Pi }}_{({\mathit{G}}_{\perp },0),({\mathit{G}}_{\perp }^{\prime },0)}^{3\mathrm{D}}(({\mathit{q}}_{\perp },0),\omega ).\end{array}$

Here in 3D the wave vector $\mathit{q}=({\mathit{q}}_{\perp },q_z)$ and $\mathit{G}=({\mathit{G}}_{\perp },G_z)$, and for both of them the $z$-components are set to 0 on the right-hand side of the equation. This is a convenient formula to use as existing DFT codes are for 3D problems.

The Dyson equation needs to be transferred into a mixed representation, that is, ${\mathit{r}}_{\perp }=(x,y)$ variables are transformed to ${\mathit{G}}_{\perp }$ and ${\mathit{q}}_{\perp }$ and $z$ variable stays. Using our general transformation, Eqs. (97) and (100), specialized to 2D, we obtain

(105)$\begin{array}{rl}D({\mathit{q}}_{\perp },\omega )=& v({\mathit{q}}_{\perp })+\\ & v({\mathit{q}}_{\perp })(\begin{array}{cc}{\mathrm{\Pi }}_0({\mathit{q}}_{\perp },\omega )& 0\\ 0& {\mathrm{\Pi }}_1({\mathit{q}}_{\perp },\omega )\end{array})D({\mathit{q}}_{\perp },\omega ).\end{array}$

Here $v$ and $D$ are matrices in the combined $z,{\mathit{G}}_{\perp }$ space, $z=0$ and $d$ only. ${\mathit{G}}_{\perp }$ runs over an energy cut-off choice. If we take only ${\mathit{G}}_{\perp }={\mathit{G}}_{\perp }^{\prime }=\mathbf{0}$, our matrices $v$ and $D$ will be $2\times 2$ which gives a result where local inhomogeneity is ignored. ${\mathrm{\Pi }}_0$ and ${\mathrm{\Pi }}_1$ are the 2D polarizability matrix located at $z=0$ and $d$, respectively. This equation for the retarded $D$ is easily solved, in the form $(I-v\mathrm{\Pi })D=v$, by calling numerical packages such as LAPACK [126].

Before closing this long theory session, we also need to transform the bare Coulomb Green’s function into mixed representation. In real space $\mathit{r}$ and time domain $t$, it is given by Eq. (53). We note that the bare Coulomb potential is a function of the difference $\mathit{r}-{\mathit{r}}^{\prime }$ only. This means in the $\mathit{G}$ representation, it is a diagonal matrix. Fourier transforming into full 3D space, we get $v(\mathit{q})=1/({ϵ}_0{q}^2)$. A two-dimensional expression is obtained if we perform an inverse Fourier transform for $q_z$ back to real space $z$ using the Cauchy residue theorem,

(106)$\begin{array}{l}{v}_{{\mathit{G}}_{\perp }{\mathit{G}}_{\perp }^{\prime }}(z,z^{\prime },{\mathit{q}}_{\perp })={\delta }_{{\mathit{G}}_{\perp },{\mathit{G}}_{\perp }^{\prime }}{\displaystyle \frac{{\mathrm{e}}^{-|{\mathit{q}}_{\perp }+{\mathit{G}}_{\perp }||z-z^{\prime }|}}{2{ϵ}_0|{\mathit{q}}_{\perp }+{\mathit{G}}_{\perp }|}}.\end{array}$

Finally, the Caroli formula remains the same with sum over $\mathit{q}$ in the first Brillouin zone and trace over ${\mathit{G}}_{\perp }$ as matrix index.

7.3 Example calculation of multiple-layer graphene

Using the first-principles method introduced above, the near-field heat flux between the monolayer graphene has been performed in Ref. [78]. As the first-principles method can be easily applied to different systems without further model treatment [127,128], we present here calculations of the heat flux between two parallel graphene sheets with finite layers. We start from the ground state calculations by using DFT as implemented in QUANTUM ESPRESSO [129,130]. The norm-conserving pseudopotential generated by the Martins-Troullier method [131] with the Perdew−Burke−Ernzerhof exchange-correlation functional [132] in the generalized gradient approximation has been adopted. The Kohn-Sham wave functions are expanded using the plane-wave basis set with a 60 Ry energy cut-off. The Fermi−Dirac smearing with a 0.002 Ry smearing width is employed to treat the partial occupancies. The in-plane lattice constants are $a=b=2.46$ Å and the inter-layer distance for multiple-layer graphene is $3.40$ Å. To avoid interactions from the neighboring lattice in the $z$ direction, a large lattice constant of $c=24.6$ Å is set in the $z$ direction of the unit cell.

The scalar photon self-energy $\mathrm{\Pi }$ of each side is calculated on top of the ground state band structure by using the BerkeleyGW package [122,123]. A $90\times 90\times 1$ Monkhorst−Pack [133] grid is used to sample the first Brillouin zone for the nonlocal polarizability, while the long-wave ($q\to 0$) polarizability is obtained from a much finer $300\times 300\times 1$ grid. To avoid divergence of the Coulomb potential, we use a small value of $q={10}^{-5}$ a.u. in the calculation of contributions from the long-wave polarizability. The energy broadening factor $\eta$ is set to $0.05$ eV, which corresponds to an electron relaxation time of ${10}^{-14}$ s [108]. We neglect the local field effects that are important only for systems with inhomogeneous geometry [134]. Then we solve the Dyson equation, Eq. (105), and calculate the transmission coefficient from the Caroli formula. Lastly, we integrate over frequencies to get the heat flux.

In Fig.9, we show the calculated heat flux of two parallel single and multiple-layer graphene sheets as a function of gap sizes. The vertical coordinate is the ratio of the calculated near-field heat flux to the black-body radiative heat flux given by the Stefan−Boltzmann law $J_{bb}=\sigma (T_0^4-T_1^4)$, with $\sigma \approx 5.67\times {10}^{-8}$ W/(m²·K⁴). As shown, the near-field heat flux is remarkably larger than that of the black-body radiation for all three systems. For monolayer graphene, a converged ratio around $5\times {10}^4$ is shown which agrees well with a previous report that used a tight-binding method to calculate the density response function of graphene [77]. The saturation of heat flux in the extreme near field originates from the nonlocal effect of wave vectors, which is a typical behavior of thermal radiation mediated by $p$-polarized evanescent waves [135]. Without spatial dispersion, the heat flux calculated from a local response function shows a $1/d$ dependence at short separation, which agrees with the previous report [136]. At extreme small distances, the heat flux between bilayer graphene approximately doubles the value of the monolayer graphene. However, heat flux between trilayer graphene sheets even becomes slightly smaller than that of the bilayer graphene. This may be due to the fact that we treat each side as infinitely thin in Eqs. (105) and (106) for simplicity. With further increases of the sample layers, this treatment is not good in the extreme near field as we assume the gap size should be larger than the inter-layer distance.

Fig.9 Distance dependence of the near-field heat flux ratio between two parallel single and multiple-layer graphene sheets. The temperatures of two sides are fixed at $T_L=$ 1000 K and $T_R=$ 300 K.

Full size|PPT slide

With an increase of the gap size, the heat flux decreases monotonically for all three systems. Nevertheless, the magnitude of the heat flux is smaller for the multiple-layer graphene. We suspect that increasing the layer number decreasing the energy transfer is due to a screening between the layers. At the distance between 7 Å to 3 nm, the influence from the finite layer is not significant and all three systems show similar results. Due to the exponential factor that appears in the 2D Coulomb potential in Eq. (106), the long-wave ($q\to 0$) contribution becomes dominant at large distances. When $d>100$nm, the heat fluxes for all three systems show an asymptotic dependence of $1/d^2$, which is consistent with the result of near-field heat flux between parallel plate capacitors [76].

Part II Vector photon and Coulomb gauge

So far in Part I, we have focused on the scalar potential and ignored the vector potential in the electrodynamics. The picture of the Coulomb interaction is a valid approximation when the length scale is shorter than some typical length scale of order micrometers or less at room temperature, but is certainly not correct for long distances. We know that the electromagnetic waves can propagate to infinity, but only from the transverse part of the field. Due to charge neutrality, the Coulomb interaction decays much faster with distance and cannot have any effect at infinity. In this part, we treat the energy transport taking into account the full electromagnetic field contributions. We study the thermal radiation from a cluster of objects modeled as a collection of tight-binding electrons. This is more than just the ideal blackbody radiation, which is independent of the details of the materials. Here again, we characterize the systems with a version of $\mathrm{\Pi }$, but it is now a tensor associated with the vector field $\mathit{A}$, or the current−current correlation at the RPA level of approximation.

8 General formulation with transverse vector field

8.1 Lagrangian and Hamiltonian, gauge invariance

To add the contribution from the transverse vector field $\mathit{A}$, we start from the Lagrangian of the scalar field version, Eq. (39), by the Peierls substitution [137,138] of the tight-binding Hamiltonian $H$, and an extra transverse field piece, obtaining

(107)$\begin{array}{rl}L=& \int \mathrm{d}V\frac{{ϵ}_0}{2}\left(-\frac{{\dot{\varphi }}^2}{{\tilde{c}}^2}+(\mathrm{\nabla }\varphi {)}^2\right)\\ & +\frac{1}{2}\int \mathrm{d}V\left[{ϵ}_0{\left(\frac{\mathrm{\partial }\mathit{A}}{\mathrm{\partial }t}\right)}^2-\frac{1}{{\mu }_0}{\left(\mathrm{\nabla }\times \mathit{A}\right)}^2\right]\\ & +\mathrm{i}\hslash c^{†}\frac{\mathrm{d}c}{\mathrm{d}t}-\sum _{j,l}{c}_j^{†}{H}_{jl}{c}_l\exp \left(-\mathrm{i}\frac{e}{\hslash }{\int }_l^j\mathit{A}\cdot \mathrm{d}\mathit{r}\right)\\ & +e\sum _j{c}_j^{†}{c}_j\varphi ({\mathit{r}}_j).\end{array}$

Here the second line is from the “kinetic” and “potential” energy of the free transverse field. The word “transverse” means that the vector field satisfies $\mathrm{\nabla }\cdot \mathit{A}=0$. The meaning is clearer if this equation is Fourier transformed into the wave-vector space, which is $\mathrm{i}\mathit{q}\cdot {\mathit{A}}_{\mathit{q}}=0$, which says that the direction of $\mathit{A}$ is perpendicular to the direction of the wavevector, thus transverse. From the Lagrangian above, we recover the Hamiltonian as (taking the limit of $\tilde{c}\to \mathrm{\infty }$)

(108)$\begin{array}{rl}\hat{H}=& \frac{{ϵ}_0}{2}\int \mathrm{d}V\left[-(\mathrm{\nabla }\varphi {)}^2+{\left(\frac{\mathrm{\partial }\mathit{A}}{\mathrm{\partial }t}\right)}^2+c^2{\left(\mathrm{\nabla }\times \mathit{A}\right)}^2\right]\\ & +\sum _{j,l}{c}_j^{†}{H}_{jl}{c}_l\exp \left(-\mathrm{i}\frac{e}{\hslash }{\int }_l^j\mathit{A}\cdot \mathrm{d}\mathit{r}\right)-e\sum _j{c}_j^{†}{c}_j\varphi ({\mathit{r}}_j).\end{array}$

Here the integral on the exponential is a line integral from site $l$ to site $j$ in a straight path. We check that the Hamiltonian is gauge invariant in the sense, that if we make the replacement,

(109)$\begin{array}{l}\varphi \to \varphi +{\displaystyle \frac{\mathrm{\partial }\chi }{\mathrm{\partial }t}},\end{array}$

(110)$\begin{array}{l}\mathit{A}\to \mathit{A}-\mathrm{\nabla }\chi ,\end{array}$

(111)$\begin{array}{l}{c}_j\to c_j\exp \left(\mathrm{i}{\displaystyle \frac{e}{\hslash }}{\chi }_j\right),\end{array}$

the result will remain the same independent of $\chi$, where ${\chi }_j=\chi ({\mathit{r}}_j,t)$ is an arbitrary smooth function of space and time. In a sense, the requirement of the gauge invariance uniquely fixes the form of the electron-photon interaction. It is noted that the Peierls substitution form of a tight-binding model has a fundamental limitation [139] as it cannot describe transitions among electronic states at the same location. We use it for its simplicity, and it is a good starting point to describe metals.

8.2 Quantization, current operator, and Green’s functions

We now discuss the quantization of the electromagnetic field. Because of the transverseness of the vector field, we see that the scalar and vector fields are not coupled, thus we can quantize $\varphi$ as in Part I. For the vector field, from the Lagrangian, the conjugate momentum is

(112)$\begin{array}{l}\mathbf{\Pi }={\displaystyle \frac{\delta L}{\delta \dot{\mathit{A}}}}={ϵ}_0\dot{\mathit{A}}.\end{array}$

However, due to the transverse nature of the field $\mathit{A}$, the three components, $A_x$, $A_y$, and $A_z$ cannot be treated as independent quantities, thus we cannot postulate commutation relation in the usual way. The true degrees of freedom are demonstrated more clearly in the Fourier space after the transformation. This is the standard approach in Coulomb gauge [2]. To make a long story short, we just give the commutation relation as

(113)$\begin{array}{l}[A_{\mu }(\mathit{r}),{\mathrm{\Pi }}_{\nu }({\mathit{r}}^{\prime })]=\mathrm{i}\hslash {\delta }_{\mu \nu }^{\mathrm{\perp }}(\mathit{r}-{\mathit{r}}^{\prime }),\end{array}$

where the right-hand side is the transverse $\delta$-function defined by an inverse Fourier transform

(114)$\begin{array}{l}{\delta }_{\mu \nu }^{\mathrm{\perp }}(\mathit{r})\equiv {\displaystyle \int \frac{{\mathrm{d}}^3\mathit{k}}{(2\pi {)}^3}\left({\delta }_{\mu \nu }-\frac{k_{\mu }{k}_{\nu }}{|\mathit{k}{|}^2}\right){\mathrm{e}}^{\mathrm{i}\mathit{k}\cdot \mathit{r}}.}\end{array}$

Here $\mu$, $\nu$ take the Cartesian directions $x$, $y$, or $z$, and ${\delta }_{\mu \nu }=1$ if $\mu =\nu$ and 0 otherwise, which is the usual Kronecker delta. It is worth noting that the transverse $\delta$-function is nonlocal and decays in space as $1/r^3$ at long distance.

Together with the earlier commutation relations for $\varphi$ and $c_j,c_j^{†}$, the problem is completely specified. We can apply the Heisenberg equations of motion for $c_j$, $\varphi$ and $\mathit{A}$, obtaining a Schrödinger equation for $c_j$ with a Peierls substituted $H$ and extra external potential due to $\varphi$, and a Poisson equation for $\varphi$ (after taking the limit $\tilde{c}\to \mathrm{\infty }$) with the usual charge density $\rho =(-e)\sum _j{c}_j^{†}{c}_j\delta (\mathit{r}-{\mathit{r}}_j)$ as the source. Finally, the equation for $\mathit{A}$ is

(115)$\begin{array}{l}{\displaystyle \frac{1}{c^2}\frac{{\mathrm{\partial }}^2\mathit{A}}{\mathrm{\partial }{t}^2}-{\mathrm{\nabla }}^2\mathit{A}={\mu }_0{\mathit{j}}_{\mathrm{\perp }}(\mathit{r}).}\end{array}$

Here the transverse current is

(116)$\begin{array}{l}{\mathit{j}}_{\mathrm{\perp }}(\mathit{r})={\displaystyle \frac{1}{\mathrm{i}\hslash }\left[\mathbf{\Pi }(\mathit{r}),\sum _{j,l}{c}_j^{†}{H}_{jl}{c}_l\exp \left(-\mathrm{i}\frac{e}{\hslash }{\int }_l^j\mathit{A}\cdot \mathrm{d}\mathit{r}\right)\right].}\end{array}$

These are consistent with Maxwell’s equations.

Due to the presence of the vector potential on the exponential, the commutator is hard to compute explicitly. But the integral is proportional to the lattice spacing $a$, which is small. In the continuum limit, we take $a$ to zero, as a result, we need to keep only to the second order in the expansion. The third and higher orders vanish in the continuum limit. But again, the transverse $\delta$-function causes some complication. Formally, we can write

(117)$\begin{array}{l}{\mathit{j}}_{\mathrm{\perp }}(\mathit{r})=P_{\mathrm{\perp }}\mathit{j}={\displaystyle \int \mathrm{d}{\mathit{r}}^{\prime }{\mathit{\delta }}^{\mathrm{\perp }}(\mathit{r}-{\mathit{r}}^{\prime })\cdot \mathit{j}({\mathit{r}}^{\prime }).}\end{array}$

Here ${\mathit{\delta }}^{\mathrm{\perp }}$ is the $3\times 3$ tensor. There are two terms to the current $\mathit{j}$, a paramagnetic term independent of the vector potential, and a diamagnetic term proportional to $\mathit{A}$, just like the electron-photon interaction in a first quantization formulation. The explicit form depends on how one approximates the line integral. Here we adopt the trapezoidal rule for the integral,

(118)$\begin{array}{l}{\displaystyle {\int }_l^j\mathit{A}\cdot \mathrm{d}\mathit{r}\approx \frac{1}{2}({\mathit{A}}_j+{\mathit{A}}_l)\cdot ({\mathit{R}}_j-{\mathit{R}}_l),}\end{array}$

where ${\mathit{R}}_l$ and ${\mathit{R}}_j$ are the respective locations of the two sites. The field is evaluated at these sites. Using this approximation, we can give an explicit formula for the paramagnetic term as

(119)$\begin{array}{l}\mathit{j}(\mathit{r})={\displaystyle \frac{1}{2}\sum _{j,l}{\mathit{I}}_{jl}(\delta (\mathit{r}-{\mathit{R}}_j)+\delta (\mathit{r}-{\mathit{R}}_l)),}\end{array}$

(120)$\begin{array}{l}{\mathit{I}}_{jl}={\displaystyle \frac{(-e)}{\mathrm{i}\hslash }{c}_j^{†}{H}_{jl}{c}_l({\mathit{R}}_j-{\mathit{R}}_l)=(-e){\hat{\mathit{V}}}_{jl}.}\end{array}$

When an electron hops from the site $l$ to $j$, it is not really clear where the current is located. It could be attributed to the middle of the sites, or one of the sites. Here we take an average of current being associated with site $j$ or $l$. The local total current obtained by volume integrating the current density around the two sites, ${\mathit{I}}_{jl}$, has a more useful interpretation; it is the velocity of the electron when hopping from site $l$ to $j$, times the charge of electron, $(-e)$.

By Taylor expanding the Peierls substituted Hamiltonian, and using the same trapezoidal approximation for the line integral, we can write the interaction part of the Hamiltonian as

(121)$\begin{array}{l}{\hat{H}}_{\mathrm{i}\mathrm{n}\mathrm{t}}={\displaystyle \sum _{jkl\mu }{c}_j^{†}{M}_{jk}^{l\mu }{c}_k{A}_{l\mu }=-\sum _l{\mathit{I}}_l\cdot {\mathit{A}}_l,}\end{array}$

where we can express the tensor $M$ in terms of the velocity matrix as $M_{jk}^{l\mu }=(e/2)({\delta }_{jl}+{\delta }_{kl})V_{jk}^{\mu }$. The index $\mu$ labels the Cartesian directions, $x$, $y$, or $z$. We also introduce the volume integrated current around the site $l$, as defined by the above equation, which we can write compactly as $I_{l\mu }=-c^{†}{M}^{l\mu }c$, where $c$ is a column vector of the annihilation operators, $c^{†}$ is a row vector of creation operators, and $M^{l\mu }$ is a Hermitian matrix in the electron site space. The current associated with the site is useful to define the current−current correlation on a discrete lattice. Finally, we can also obtain the current density, Eq. (119), by the functional derivative of the interaction Hamiltonian with respect to the vector field by $\mathit{j}(\mathit{r})=-\delta {\hat{H}}_{\mathrm{i}\mathrm{n}\mathrm{t}}/\delta \mathit{A}(\mathit{r})$.

The diamagnetic term is usually not important as it is higher order, also, it contributes only a purely real, diagonal $\mathrm{\Pi }$, so it is not dissipative. In the continuum limit, we have a simple expression for the current ${\mathit{j}}_A=-e^2(n/m)\mathit{A}$, which is the London equation [140]. Here $n$ is the electron density and $m$ is its mass. On lattice, if we use the trapezoidal approximation, we get a complicated form of the type $H_{\mathrm{i}\mathrm{n}\mathrm{t}}^{A^2}=\frac{1}{2}\sum _{l,l^{\prime },\mu ,\nu ,j,k}{c}_j^{†}{N}_{jk}^{l\mu ,l^{\prime }\nu }{c}_k{A}_{l\mu }{A}_{l^{\prime }\nu }$.

Finally, based on the vector field $\mathit{A}$ and the current ${\mathit{I}}_l$ we can define two similar contour ordered Green’s functions, call them $D$ and $\mathrm{\Pi }$ as they play very similar role as in the scalar photon theory, except that $D$ and $\mathrm{\Pi }$ will have additional directional indices as $3\times 3$ tensors or dyadic.

(122)$\begin{array}{l}{D}_{\mu \nu }(\mathit{r}\tau ;{\mathit{r}}^{\prime }{\tau }^{\prime })={\displaystyle \frac{1}{\mathrm{i}\hslash }⟨T_{\tau }{A}_{\mu }(\mathit{r},\tau )A_{\nu }({\mathit{r}}^{\prime },{\tau }^{\prime })⟩,}\end{array}$

(123)$\begin{array}{l}{\mathrm{\Pi }}_{l\mu ,l^{\prime }\nu }(\tau ;{\tau }^{\prime })={\displaystyle \frac{1}{\mathrm{i}\hslash }⟨T_{\tau }{I}_{l\mu }(\tau )I_{l^{\prime }\nu }({\tau }^{\prime }){⟩}_{\mathrm{i}\mathrm{r}}.}\end{array}$

Here $D$ is defined in the full space of $\mathit{r}$, while $\mathrm{\Pi }$ is only on the discrete sites. In computing $\mathrm{\Pi }$ under the random phase approximation, we take only the lowest order in such an expansion, ignoring the $A^2$ terms in the Hamiltonian. With the electron-photon interaction of the form $-\sum \mathit{I}\cdot \mathit{A}$, we still have the standard Dyson equation of the form, $D=d+d\mathrm{\Pi }D$, except that the sizes of the matrices are 3 times larger.

8.3 Free Green’s function

We need to determine the free photon Green’s function $d$ associated with the transverse vector field $\mathit{A}$ (earlier for the scalar field $\varphi$ it was called $v$). This can be done with the standard procedure of quantization of the field. We can express the vector field with the annihilation and creation operator as [94,141]

(124)$\begin{array}{l}\mathit{A}(\mathit{r})={\displaystyle \sum _{\mathit{q}\lambda }\sqrt{\frac{\hslash }{2{ϵ}_0{\omega }_{\mathit{q}}V}}{\mathit{e}}_{\mathit{q}\lambda }{a}_{\mathit{q}\lambda }{\mathrm{e}}^{\mathrm{i}\mathit{q}\cdot \mathit{r}}+\mathrm{h}.\mathrm{c}.,}\end{array}$

where ${\omega }_{\mathit{q}}=c|\mathit{q}|$ is the free photon dispersion relation, ${\mathit{e}}_{\mathit{q}\lambda }$ are the two mutually perpendicular unit polarization vectors indicated by the index $\lambda$, which are also orthogonal to $\mathit{q}$. This condition gives a transverse field for $\mathit{A}$, i.e., $\mathbf{\nabla }\cdot \mathit{A}=0$. $a_{\mathit{q}\lambda }$ is the associated annihilation operator for the mode $\mathit{q},\lambda$. These are the standard bosonic operators satisfying $[a_{\mathit{q}\lambda },a_{{\mathit{q}}^{\prime }{\lambda }^{\prime }}]=0$, $[a_{\mathit{q}\lambda },a_{{\mathit{q}}^{\prime }{\lambda }^{\prime }}^{†}]={\delta }_{\mathit{q}{\mathit{q}}^{\prime }}{\delta }_{\lambda {\lambda }^{\prime }}$. The field exists in a finite volume $V$ with periodic boundary conditions, thus the wave vectors $\mathit{q}$ are discrete. This equation together with a corresponding equation for the conjugate field $\mathbf{\Pi }(\mathit{r})$ is viewed as a variable transformation between $a_{\mathit{q}\lambda }$ and $\mathit{A}(\mathit{r})$. This is the correct transformation if the commutation relation between $\mathit{A}$ and $\mathbf{\Pi }$, Eq. (113), is fulfilled and the Hamiltonian is diagonal, ${\hat{H}}_{\gamma }=$ $\sum _{\mathit{q}\lambda }\hslash {\omega }_{\mathit{q}\lambda }(a_{\mathit{q}\lambda }^{†}{a}_{\mathit{q}\lambda }+1/2)$. Indeed, these claims can be verified.

For the free field, the time dependence for the annihilation operator is trivially $a_{\mathit{q}\lambda }(t)=a_{\mathit{q}\lambda }{\mathrm{e}}^{-\mathrm{i}{\omega }_{\mathit{q}}t}$. Using this result, plugging into the definition of the retarded Green’s function, we obtain

(125)$\begin{array}{l}{d}_{\mu \nu }^r(\mathit{r},t)={\displaystyle \frac{1}{{ϵ}_{0}V}\sum _{\mathit{q}}\left({\delta }_{\mu \nu }-\frac{q_{\mu }{q}_{\nu }}{q^2}\right){\mathrm{e}}^{\mathrm{i}\mathit{q}\cdot \mathit{r}}{d}_{\mathit{q}}(t),}\end{array}$

(126)$\begin{array}{l}{d}_{\mathit{q}}(t)=-\theta (t){\displaystyle \frac{\sin ({\omega }_{\mathit{q}}t)}{{\omega }_{\mathit{q}}}{\mathrm{e}}^{-\eta t}.}\end{array}$

The factor in the brackets takes care of the transverseness of the Green’s function, which is the transverse projector in $\mathit{q}$ space. The transverse projector appears because two of the polarization vectors and the unit vector $\mathit{q}/q$ form an orthonormal basis. We can use the completeness relation to eliminate the reference to the polarization vectors, ${\mathit{e}}_{\mathit{q}\lambda }$. The second line defines the retarded Green’s function for a single mode in time domain. In frequency domain, it is $1/[(\omega +\mathrm{i}\eta {)}^2-{\omega }_{\mathit{q}}^2]$. If we take the volume to infinity, the summation can be turned into an integral in $\mathit{q}$. The final expression in real space and frequency domain can be obtained with the help of the residue theorem with the contour integral [142] on the complex plane of $q$, given as [143]

(127)$\begin{array}{rl}\overleftrightarrow{d^r}(\mathit{r},\omega )=& -\frac{1}{4\pi {ϵ}_0{c}^{2}r}\{{\mathrm{e}}^{\mathrm{i}\frac{\omega }{c}r}(\overleftrightarrow{\mathit{U}}-\hat{\mathit{R}}\hat{\mathit{R}})\\ & +[-\frac{{\mathrm{e}}^{\mathrm{i}\frac{\omega }{c}r}}{\mathrm{i}\frac{\omega }{c}r}+\frac{{\mathrm{e}}^{\mathrm{i}\frac{\omega }{c}r}-1}{{\left(\mathrm{i}\frac{\omega }{c}r\right)}^2}](\overleftrightarrow{\mathit{U}}-3\hat{\mathit{R}}\hat{\mathit{R}})\}.\end{array}$

Here, we express the $3\times 3$ tensor in the Cartesian directions in the dyadic notation, $r=|\mathit{r}|$ is the magnitude of the vector, and $\hat{\mathit{R}}=\mathit{r}/r$ is a unit vector in the radial direction. $\overleftrightarrow{\mathit{U}}$ is the unit or identity matrix.

Another equivalent way to obtain the retarded free Green’s function is from the equation it must satisfy. We observe that the transverse vector field $\mathit{A}$ satisfies a wave equation with the transverse current as a source, Eq. (115). The Green’s function maps the current to the vector field, $\mathit{A}=-\overleftrightarrow{d^r}\cdot \mathit{j}$. The Green’s function satisfies the same equation as the vector potential, but with a transverse $\delta$ function as a source. When the wave equation and the transverse $\delta$ function is expressed in $\omega$ and $\mathit{q}$ space, we have a simple result for the Green’s function as

(128)$\begin{array}{l}\overleftrightarrow{d^r}(\mathit{q},\omega )={\displaystyle \frac{\overleftrightarrow{\mathit{U}}-\hat{\mathit{q}}\hat{\mathit{q}}}{{ϵ}_0[(\omega +\mathrm{i}\eta {)}^2-c^2{q}^2]},}\end{array}$

where $\hat{\mathit{q}}=\mathit{q}/|\mathit{q}|$ is a unit vector. Fourier integral transforming into real space $\mathit{r}$, we obtain the explicit formula of Eq. (127).

In computing transport quantities in a planar geometry, such as between graphene sheets or surface of a lattice, we need the free Green’s function in mixed representation, i.e., the transverse direction is in wave vector space, ${\mathit{q}}_{\perp }$, but the transport direction, say $z$, is in real space. The retarded Green’s function in the mixed representation can be obtained by inverse Fourier transform $q_z$ back to real space,

(129)$\begin{array}{rl}& d^{r,\mu \nu }({\mathit{q}}_{\perp },\omega ,z,z^{\prime })\\ & ={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\mathrm{d}{q}_z}{2\pi }\frac{\left({\delta }_{\mu \nu }-\frac{q_{\mu }{q}_{\nu }}{q^2}\right){\mathrm{e}}^{\mathrm{i}{q}_z(z-z^{\prime })}}{a^2{ϵ}_0[(\omega +\mathrm{i}{0}^{+}{)}^2-c^2{q}_{\perp }^2-c^2{q}_z^2]},\end{array}$

where $q^2=|\mathit{q}{|}^2=q_{\perp }^2+q_z^2$. Here $a^2$ is the area of a unit cell, since we will consider a discrete set of ${\mathit{q}}_{\perp }$. This integral can be performed using the residue theorem. After somewhat lengthy and tedious calculation, we obtain [144]

(130)$\begin{array}{l}{d}^{r,\alpha \beta }={\delta }_{\alpha \beta }d-q_{\alpha }{q}_{\beta }F,\alpha ,\beta =x,y,\end{array}$

(131)$\begin{array}{l}{d}^{r,\alpha z}=d^{r,z\alpha }=\mathrm{s}\mathrm{g}\mathrm{n}(z-z^{\prime })q_{\alpha }(B-A)/C,\end{array}$

(132)$\begin{array}{l}{d}^{r,zz}=q_{\perp }^{2}F.\end{array}$

We have introduced the shorthand notations $A=$ ${\mathrm{e}}^{\mathrm{i}{\tilde{q}}_z|z-z^{\prime }|}$, $B={\mathrm{e}}^{-q_{\perp }|z-z^{\prime }|}$, $d=A/(a^2{ϵ}_{0}2\mathrm{i}{c}^2{\tilde{q}}_z)$, $F=(A/{\tilde{q}}_z+$ $\mathrm{i}B/q_{\perp })/C$, $C=a^2{ϵ}_{0}2\mathrm{i}{\omega }^2$, and ${\tilde{q}}_z=\pm \sqrt{[(\omega +\mathrm{i}{0}^{+})/c{]}^2-q_{\perp }^2}$, where the sign is chosen such that $\mathrm{I}\mathrm{m}{\tilde{q}}_z>0$. We note that the free Green’s function is transverse in the sense $\mathrm{i}{q}_x{d}^{r,x\beta }+\mathrm{i}{q}_y{d}^{r,y\beta }+\frac{\mathrm{\partial }}{\mathrm{\partial }z}{d}^{r,z\beta }=0$ for $\beta =x,y,z$.

9 Poynting vector and energy transport

At a place in vacuum without material presence, the energy can be transferred through the electromagnetic waves. This is described by the Poynting vector, $\mathit{S}=\mathit{E}\times \mathit{B}/{\mu }_0$, in electrodynamics. The meaning of the Poynting vector is made clear through Poynting’s theorem which reflects the energy conservation [145],

(133)$\begin{array}{l}{\displaystyle \frac{\mathrm{\partial }u}{\mathrm{\partial }t}}+\mathrm{\nabla }\cdot \mathit{S}=-\mathit{E}\cdot \mathit{j},\end{array}$

here $u=1/2({ϵ}_0{E}^2+B^2/{\mu }_0)$ is the field energy density. The right-hand side is the Joule heating contribution where $\mathit{j}$ is the electric current density. If we consider a volume $V$ and integrate over the volume, using Gauss’s theorem, the divergence over the volume becomes a surface integral, so

(134)$\begin{array}{l}I={\displaystyle \int \mathrm{d}\mathbf{\Sigma }\cdot \mathit{S}}\end{array}$

is the energy flux going out of the volume where $\mathrm{d}\mathbf{\Sigma }$ is the surface element with an outward norm.

In order to use this classical expression for the energy current, we need to consider several additional features in a quantum theory. We need to perform an ensemble average by expressing the Poynting vector in terms of the Green’s function $D$. We need to worry about the operator order as $\mathit{E}$ and $\mathit{B}$ in general are not commuting quantities. We also need to remove the zero-point motion contribution to the Poynting vector, as intuitively we do not expect that the zero-point motion of the electromagnetic waves transfers energy when objects are stationary (see, however Ref. [42]). Lastly, we need to worry about our choice of gauge. Since we have used the Coulomb gauge here, we have $\mathit{E}=-\mathrm{\nabla }\varphi -\mathrm{\partial }\mathit{A}/\mathrm{\partial }t$, here the vector field is transverse. This split of $\mathit{E}$ into a scalar potential term and a transverse vector potential term means we can also split the Poynting vector into two corresponding terms, $\mathit{S}={\mathit{S}}_{//}+{\mathit{S}}_{\perp }$, where ${\mathit{S}}_{//}=-\mathrm{\nabla }\varphi \times \mathit{B}/{\mu }_0$ and ${\mathit{S}}_{\perp }=-\mathrm{\partial }\mathit{A}/\mathrm{\partial }t\times (\mathrm{\nabla }\times \mathit{A})/{\mu }_0$. The longitudinal contribution ${\mathit{S}}_{//}$ can be transformed back to a volume integral by the divergence theorem and the curl of $\mathit{B}$ is then related to the electric current. In steady state, the average of a time derivative is zero. Using this property and Maxwell’s equations, we can show that ${\mathit{S}}_{//}$ outside matter is the same as the “Poynting scalar” discussed earlier, plus a cross correlation term between the longitudinal and transverse fields, $-{ϵ}_0\varphi {\mathrm{\partial }}^2\mathit{A}/\mathrm{\partial }{t}^2$. In this section, we shall focus on the transverse contribution ${\mathit{S}}_{\perp }$.

As for the remaining issues, in general for a product $AB$ of two Hermitian operators, $A$ and $B$, the result is not Hermitian, and its expectation is not guaranteed to be real. Thus, we will replace the product by the symmetrized version, $(AB+BA)/2$. To remove the zero-point motion contribution, we need to impose a normal order [61,141], so the final form is $:(AB+BA):/2$. We elaborate this point in the next subsection in some detail.

9.1 Operator order

Given a Schrödinger picture Hermitian operator $A$, the Heisenberg operator $A(t)$ is assumed to be defined for all $t$, by $A(t)={\mathrm{e}}^{\mathrm{i}\hat{H}t/\hslash }A{\mathrm{e}}^{-\mathrm{i}\hat{H}t/\hslash }$. We Fourier decompose the operator as

(135)$\begin{array}{rl}A(t)& =A^{+}(t)+A^{-}(t)\\ & ={\int }_0^{+\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\tilde{A}(\omega ){\mathrm{e}}^{-\mathrm{i}\omega t}+{\int }_{-\mathrm{\infty }}^0\frac{\mathrm{d}\omega }{2\pi }\tilde{A}(\omega ){\mathrm{e}}^{-\mathrm{i}\omega t}.\end{array}$

$A^{+}(t)$ is the positive frequency part of $A(t)$ by integrating over the positive frequencies in the Fourier space, and $A^{-}(t)$ is called the negative frequency part of $A(t)$. The positive part of the frequency is associated with annihilation operators and negative one with creation operators [146]. The normal order or anti-normal order is defined in terms of $A^{\pm }(t)$. The fundamental assumptions for the operators are

(136)$\begin{array}{l}{A}^{+}(t)|0⟩=0,⟨0|A^{-}(t)=0,\end{array}$

where $|0⟩$ is the vacuum state. What we have in mind is the free photon field, but we assume it is generally valid. Also, since $A(t)$ is Hermitian, $A^{+}(t)=[A^{-}(t){]}^{†}$.

Consider a steady-state Green’s function formed by two operators, $A$ and $B$, and defined as

(137)$\begin{array}{l}{D}_{AB}^{>}(t)={\displaystyle \frac{1}{\mathrm{i}\hslash }}⟨A(t)B(0)⟩,B(0)=B.\end{array}$

The decomposition of the positive and negative frequency parts of $A(t)$ naturally leads to positive and negative frequency parts of $D_{AB}^{>}(t)$, thus we must have

(138)$\begin{array}{l}{\displaystyle \frac{1}{\mathrm{i}\hslash }⟨A^{+}(t)B⟩={\int }_0^{+\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }{\tilde{D}}_{AB}^{>}(\omega ){\mathrm{e}}^{-\mathrm{i}\omega t},}\end{array}$

(139)$\begin{array}{l}{\displaystyle \frac{1}{\mathrm{i}\hslash }⟨A^{-}(t)B⟩={\int }_{-\mathrm{\infty }}^0\frac{\mathrm{d}\omega }{2\pi }{\tilde{D}}_{AB}^{>}(\omega ){\mathrm{e}}^{-\mathrm{i}\omega t}.}\end{array}$

That is, the positive frequency part $A^{+}(t)$ contributes to positive $\omega$ in $D^{>}$ only. Similarly, $D_{AB}^{<}(t)=⟨BA(t)⟩/(\mathrm{i}\hslash )$ can be decomposed analogously. The normal order is defined by

(140)$\begin{array}{l}:A^{-}{B}^{+}:=A^{-}{B}^{+},\end{array}$

(141)$\begin{array}{l}:B^{+}{A}^{-}:=A^{-}{B}^{+},\end{array}$

(142)$\begin{array}{l}:A^{+}{B}^{+}:=A^{+}{B}^{+},\end{array}$

(143)$\begin{array}{l}:A^{-}{B}^{-}:=A^{-}{B}^{-},\end{array}$

i.e., the right-most operator should be the annihilation operator if it is not already so. The normal order has no effect if it is already normal ordered. We assume that the normal order has no effect if both are +, or −. The normal ordering is a linear operation. So we can express the equal-time correlation as

(144)$\begin{array}{l}⟨:AB:⟩=⟨[A^{+},B^{+}]⟩+⟨A^{-}B+B{A}^{+}⟩,\end{array}$

(145)$\begin{array}{l}⟨:BA:⟩=⟨[B^{-},A^{-}]⟩+⟨A^{-}B+B{A}^{+}⟩.\end{array}$

Our basic assumption is that the positive (negative) frequency part contains only annihilation (creation) operators, and they commute. So the first term is zero for whatever meaning of the average. As a result, we have, consisting with [147],

(146)$\begin{array}{rl}⟨:AB:⟩& =⟨:BA:⟩=⟨A^{-}B+B{A}^{+}⟩\\ & =\mathrm{R}\mathrm{e}\mathrm{i}\hslash {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{\pi }{\tilde{D}}_{AB}^{<}(\omega ).\end{array}$

Here we use the Fourier decomposition at time $t=0$ and use a symmetry of the Green’s function, ${\tilde{D}}_{AB}^{>}(\omega )=$ ${\tilde{D}}_{BA}^{<}(-\omega )$. We see that the two terms, $A^{-}B$ and $B{A}^{+}$, are Hermitian conjugate of each other, so the result is explicitly real. The last formula, Eq. (146), has the effect of removing the zero-point motion contribution by taking integral only for the positive frequencies for the lesser Green’s function, remembering that due to the fluctuation-dissipation theorem, $D^{<}$ is proportional to the Bose function $N(\omega )$ in thermal equilibrium which decays to 0 at high frequency exponentially. If this were $D^{>}$, we have $N(\omega )+1$, which may run into problem of divergence when integrated over the positive frequencies. We note that a symmetrization of the operator product is not necessary as the normal order automatically makes the result symmetric with respective to the product order.

9.2 Average transverse Poynting vector

With the preparation in the above subsection of a formula to express the equal-time correlation of a normal ordered product as a positive frequency integral of the lesser Green’s function in the frequency domain, we can work out an expression of the average of the Poynting vector,

(147)$\begin{array}{rl}⟨S_{\perp }^i⟩& =\frac{1}{{\mu }_0}⟨:({\mathit{E}}_{\perp }\times \mathit{B}{)}_i:⟩\\ & =\frac{1}{{\mu }_0}⟨[:(-\dot{\mathit{A}})\times (\mathrm{\nabla }\times \mathit{A}):{]}_i⟩\\ & =-\frac{1}{{\mu }_0}\sum _{ijklm}{ϵ}_{ijk}{ϵ}_{klm}{\left[\frac{\mathrm{\partial }}{\mathrm{\partial }t}\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_l^{\prime }}⟨:A_j(\mathit{r},t)A_m({\mathit{r}}^{\prime },t^{\prime }):⟩\right]}_{\stackrel{t^{\prime }=t}{{\mathit{r}}^{\prime }=\mathit{r}}}\\ & =\mathrm{R}\mathrm{e}\sum _{ijklm}{ϵ}_{ijk}{ϵ}_{klm}{\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{{\mu }_0\pi }\hslash \omega {\left[\frac{\mathrm{\partial }}{\mathrm{\partial }{x}_l^{\prime }}{D}_{mj}^{<}({\mathit{r}}^{\prime },\mathit{r},\omega )\right]}_{{\mathit{r}}^{\prime }=\mathit{r}}.\end{array}$

Here for the Poynting vector, we look for the $i$-th Cartesian component, expressed in terms of the transverse vector field $\mathit{A}$. The average Poynting vector depends on explicitly the location $\mathit{r}$ which we have suppressed in notation. Since we are interested in steady-state average, the Poynting vector does not depend on time. The vector cross products are written out explicitly with the Levi-Civita symbol, ${ϵ}_{ijk}$, which is ${ϵ}_{xyz}=1$ and antisymmetric for each permutation of any two indices, e.g., ${ϵ}_{ijk}=-{ϵ}_{jik}$. In order to express the final result in terms of an $AA$ correlation, i.e., the $D^{<}$ photon Green’s function, we need to pull the time and space derivatives outside the average by changing the variable $\mathit{r}$ to ${\mathit{r}}^{\prime }$ and $t$ to $t^{\prime }$, and changing them back to $\mathit{r}$ and $t$ after the derivatives are performed. We use Eq. (146) in the last step, identifying $A_m({\mathit{r}}^{\prime },t^{\prime })$ as the operator $A$ and $A_j(\mathit{r},t)$ as $B$. Notice that time derivative in time domain becomes $(+\mathrm{i}\omega )$ in the frequency domain as $t$ is the second argument in $D^{<}$.

To compute the total energy current, we need to make a dot product with the surface norm $\mathit{n}$ and to integrate over the surface. We can simplify the formula a bit by summing over the index $k$ in the product of Levi−Civita symbols,

(148)$\begin{array}{l}{\displaystyle \sum _k{ϵ}_{ijk}{ϵ}_{klm}={\delta }_{il}{\delta }_{jm}-{\delta }_{im}{\delta }_{jl}.}\end{array}$

Using this identity, we can write the surface integral for the energy transfer as

(149)$\begin{array}{l}I={\displaystyle \int \mathrm{d}\mathrm{\Sigma }{\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{{\mu }_0\pi }\hslash \omega \mathrm{R}\mathrm{e}[(\mathit{n}\cdot {\mathbf{\nabla }}^{\prime })\mathrm{T}\mathrm{r}(D^{<})-\mathrm{T}\mathrm{r}({\mathbf{\nabla }}^{\prime }\mathit{n}\cdot D^{<})].}\end{array}$

Here ${\mathbf{\nabla }}^{\prime }$ is the gradient operator acting on the first argument of $D^{<}$, the combination ${\mathbf{\nabla }}^{\prime }\mathit{n}$ is a dyadic, i.e., as a $3\times 3$ matrix with component $\mathrm{\partial }/\mathrm{\partial }{x}_i^{\prime }{n}_j$ as a differential operator acting on $D^{<}$, $\mathit{n}$ is the unit norm of the surface element, $\mathrm{d}\mathrm{\Sigma }$ is the magnitude of the surface element, the dot · is a scalar product in the first term and matrix multiplication in the second term. The trace is over the matrix indices.

9.3 Radiation at far field

We consider a cluster of some materials of finite sizes and compute the radiation at far field. At sufficiently far distances from the matter, the longitudinal part from ${\mathit{S}}_{//}$ decays to zero and only the transverse field in ${\mathit{S}}_{\perp }$ propagates to infinity. So for the energy radiation to infinity, we only need to compute the contribution from the above formula. In the far field, we have another simplification which we can use, that is, the electromagnetic waves at far field are spherical waves of the form ${\mathrm{e}}^{\mathrm{i}\frac{\omega }{c}R}$, where $R=|\mathit{R}|$ is the distance to the coordinate origin. Since the photon Green’s function satisfies the same equation as the field, the Green’s function also takes this form. As a result, the gradient operator ${\mathbf{\nabla }}^{\prime }$ can be replaced by the vector $\mathrm{i}\frac{\omega }{c}\hat{\mathit{R}}$, here $\hat{\mathit{R}}=\mathit{R}/R$ is the unit vector from the origin to the observation point $\mathit{R}$. We can choose a large sphere of radius $R$ to perform the surface integral as a solid angle integration, $\mathrm{d}\mathrm{\Sigma }=\mathrm{d}\mathrm{\Omega }{R}^2$. Using this observation and replacing the general surface norm $\mathit{n}$ by $\hat{\mathit{R}}$, we obtain

(150)$\begin{array}{l}I=\mathrm{R}\mathrm{e}{\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{c{\mu }_0\pi }(\mathrm{i}\hslash {\omega }^2)\int \mathrm{d}\mathrm{\Omega }{R}^2\mathrm{T}\mathrm{r}[(\overleftrightarrow{\mathit{U}}-\hat{\mathit{R}}\hat{\mathit{R}})D^{<}].}\end{array}$

Here $\overleftrightarrow{\mathit{U}}$ is the unit matrix, $P=\overleftrightarrow{\mathit{U}}-\hat{\mathit{R}}\hat{\mathit{R}}$ is the transverse projector, having the property $P^2=P$. Two matrices $P$ and $D^{<}$ are multiplied and then trace is taken.

We can express $D^{<}$ by the Keldysh equation as $D^{<}=D^r{\mathrm{\Pi }}^{<}{D}^a$, where $D^r$ is the retarded Green’s function, and $D^a=(D^r{)}^{†}$ is the advanced version, and ${\mathrm{\Pi }}^{<}$ is the material property we call self-energy (with respect to the photons). Here the three matrices are multiplied, which implies a sum over the sites as well as directions, i.e.,

(151)$\begin{array}{l}{D}_{ij}^{<}(\mathit{r},{\mathit{r}}^{\prime },\omega )={\displaystyle \sum _{l{l}^{\prime }\alpha \beta }{D}_{i\alpha }^r(\mathit{r},{\mathit{r}}_l,\omega ){\mathrm{\Pi }}_{l\alpha ,l^{\prime }\beta }^{<}(\omega )D_{\beta j}^a({\mathit{r}}_{l^{\prime }},{\mathit{r}}^{\prime },\omega ).}\end{array}$

We still have to solve the Dyson equation for the retarded Green’s function, $D^r=d^r+d^r{\mathrm{\Pi }}^r{D}^r$. However, for the far field problem, the corrections to the free Green’s function are rather small since ${\mathrm{\Pi }}^r{D}^r\sim (v/c{)}^2$ is of the order of the ratio of electron velocity to the speed of light squared. As a result, it is sufficient and is a good approximation just using the free Green’s function $d^r$. Also, since $\mathit{r}=\mathit{R}\to \mathrm{\infty }$, it does not matter where the atoms are located with respect to $\mathit{R}$, so we also set all the positions ${\mathit{r}}_l$ in the argument of the Green’s function at the origin $\mathbf{0}$. We call this a monopole approximation. Then our approximate Green’s function is

(152)$\begin{array}{l}{D}^r(\mathit{R},{\mathit{r}}_l,\omega )\approx -{\displaystyle \frac{1}{4\pi {ϵ}_0{c}^{2}R}}{\mathrm{e}}^{\mathrm{i}\frac{\omega }{c}R}(\overleftrightarrow{\mathit{U}}-\hat{\mathit{R}}\hat{\mathit{R}}).\end{array}$

Putting this result into Eq. (150), using the cyclic property of trace, and the fact that $P$ is a projector, $P^3=P$, we find

(153)$\begin{array}{l}I=\mathrm{R}\mathrm{e}{\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}\omega \frac{\mathrm{i}\hslash {\omega }^2}{16{\pi }^3{ϵ}_0{c}^3}\int \mathrm{d}\mathrm{\Omega }\mathrm{T}\mathrm{r}[(\overleftrightarrow{\mathit{U}}-\hat{\mathit{R}}\hat{\mathit{R}}){\mathrm{\Pi }}^{<}(\omega )].}\end{array}$

Here the total ${\mathrm{\Pi }}^{<}(\omega )=\sum _{l,l^{\prime }}{\mathrm{\Pi }}_{l{l}^{\prime }}^{<}$ is the sum over all the sites. Since the only angular dependence is in $\hat{\mathit{R}}=(\sin \theta \cos \varphi ,\sin \theta \sin \varphi ,\cos \theta )$, the integration over the solid angle divided by $4\pi$ is equivalent to an average over the projector. We can verify by a direct integration that

(154)$\begin{array}{l}\overline{\hat{\mathit{R}}\hat{\mathit{R}}}={\displaystyle \frac{1}{4\pi }\int \mathrm{d}\mathrm{\Omega }\hat{\mathit{R}}\hat{\mathit{R}}=\frac{1}{3}\overleftrightarrow{\mathit{U}}.}\end{array}$

Using this simple result, we obtain the final expression for energy radiation as [148]

(155)$\begin{array}{l}I={\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}\omega \frac{-\hslash {\omega }^2}{6{\pi }^2{ϵ}_0{c}^3}\mathrm{I}\mathrm{m}\mathrm{T}\mathrm{r}({\mathrm{\Pi }}^{<}(\omega )).}\end{array}$

As a simple application of this formula, we reproduce the textbook result of dipole radiation. Consider a charge moving according to $z(t)=a\cos ({\omega }_{0}t)=$ $\frac{1}{2}a({\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}+{\mathrm{e}}^{+\mathrm{i}{\omega }_{0}t})$ in $z$ direction. $a$ is the amplitude of the oscillation with frequency ${\omega }_0$. The velocity is $v(t)=\mathrm{d}z(t)/\mathrm{d}t=\frac{1}{2}(-\mathrm{i}{\omega }_{0}a)({\mathrm{e}}^{-\mathrm{i}{\omega }_{0}t}-{\mathrm{e}}^{+\mathrm{i}{\omega }_{0}t})$. We need the total current (current density integrated over the volume), which is $J(t)=qv(t)$, here $q$ is the charge of the particle. The self-energy $\mathrm{\Pi }$ is the current-current correlation. As we are dealing with a classical charge, and athermal, we do not distinguish ${\mathrm{\Pi }}^{<}$ with ${\mathrm{\Pi }}^{>}$ and just call it $\mathrm{\Pi }$. $\mathrm{\Pi }$ is a 3 by 3 matrix, but since the particle is moving in $z$ direction, we only have the ${\mathrm{\Pi }}_{zz}$ component nonzero, which is, in real time,

(156)$\begin{array}{l}\mathrm{\Pi }(t,t^{\prime }{)}_{zz}={\displaystyle \frac{q^2}{\mathrm{i}\hslash }⟨v(t)v(t^{\prime })⟩.}\end{array}$

We see if we plug in the formula for the velocity, we do not have a time-translationally invariant result. This is because we assume that the oscillator has no damping. If we introduce a damping to the oscillator, then the $t-t^{\prime }$ dependence remains but $t+t^{\prime }$ will be damped out at long time. So performing the average has the effect of getting rid of the $t+t^{\prime }$ dependence. Fourier transform the result into frequency domain, we find

(157)$\begin{array}{l}{\mathrm{\Pi }}_{zz}(\omega )={\displaystyle \frac{2\pi ({\omega }_{0}p{)}^2}{\mathrm{i}4\hslash }[\delta (\omega +{\omega }_0)+\delta (\omega -{\omega }_0)],}\end{array}$

here $p=qa$ is the dipole moment. The correlation is $\delta$ peaked at $\pm {\omega }_0$. Putting this expression in the radiation power formula, we find

(158)$\begin{array}{l}I={\displaystyle \frac{{\omega }_0^4{p}^2}{12\pi {ϵ}_0{c}^3}.}\end{array}$

This is the result in electrodynamics for dipole radiation [68].

10 A two-dot model

The simple example given above does not involve thermal equilibrium (or nonequilibrium). In this section, we give a different example where the material system is modeled explicitly. Obviously, an isolated electric dipole cannot oscillate forever. In order to have a steady state established, we must supply continuously energy. This energy is pumped in through the effect of the baths. Here we consider a two-dot model consisting of two electron sites, see Fig.10. An electron can hop from the left bath into the left dot, and then it may hop again to the right dot, and eventually hop to the right bath. This is a toy model for electroluminescence — the emission of light due to electric current [149,150]. The hopping between the dots generates electric current, which couples the electron with photons in space, generating radiation. A special situation is that when the temperatures and chemical potentials of the two baths are equal, then it is a thermal radiation problem. The role of the baths is to supply and to dump the electrons, and we assume that they do not couple to the field.

Fig.10 (a) Two-dot model connected to the left and right baths with strength $\eta$ electron hops between two sites with hopping parameter $t$. (b) Energy levels of the two-dot model and relation to the bath distribution functions. When electron jumps from the excited state to the ground state, it releases energy $2t=\hslash \omega$.

Full size|PPT slide

The double-dot Hamiltonian before coupling to the baths is given by

(159)$\begin{array}{l}{\hat{H}}_C=-t(c_1^{†}{c}_2+c_2^{†}{c}_1)=(c_1^{†}{c}_2^{†})H_C(\begin{array}{c}{c}_1\\ c_2\end{array}).\end{array}$

Here $c_1^{†}$ creates an electron at site 1, and $c_2^{†}$ at site 2. The isolated center has four many-body eigenstates with no particles, energy $E_1=0$, and the vacuum state ${\mathrm{\Psi }}_1=|0⟩$, one particle at the ground state of the one-particle state with energy $E_2=-t$ and ${\mathrm{\Psi }}_2=$ $\frac{1}{\sqrt{2}}(c_1^{†}+c_2^{†})|0⟩$ in a symmetric combination, and one particle at excited state with $E_3=+t$ with an antisymmetric combination, ${\mathrm{\Psi }}_3=\frac{1}{\sqrt{2}}(c_1^{†}-c_2^{†})|0⟩$, and finally, $E_4=0$ and ${\mathrm{\Psi }}_4=c_1^{†}{c}_2^{†}|0⟩$ with both sites occupied. The energy levels of the isolated center give us a clear picture of when light will be emitted based on energy conservation. If an electron flows from left to right without hopping between the one-particle states, it cannot emit photon. The process that emits photon is the one that an electron comes at the excited state and stays there for sufficiently long time, then it spontaneously jumps to the ground state and loses energy to the space at infinity.

We can consider the two sites as part of a 1D chain exposed to the coupling to photon field. In such a strong system-bath coupling regime, the emission turns out rather weak as most of the time, the electrons simply travel through the wire without emission. So a weak coupling of the double-dot to the baths is preferred. We use the simplest possible bath coupling as the wide-band model, where the coupling is constant independent of the energy. With the effect of the baths, the retarded Green’s function of the electron is given by

(160)$\begin{array}{l}(E-H_C-{\mathrm{\Sigma }}^r)G^r=I,\end{array}$

where ${\mathrm{\Sigma }}^r={\mathrm{\Sigma }}_L^r+{\mathrm{\Sigma }}_R^r$ is the self-energy due to the baths. In the wide-band approximation, we take

(161)$\begin{array}{l}{\mathrm{\Sigma }}_L^r=(\begin{array}{cc}-\mathrm{i}\eta & 0\\ 0& 0\end{array}),{\mathrm{\Sigma }}_R^r=(\begin{array}{cc}0& 0\\ 0& -\mathrm{i}\eta \end{array}),\end{array}$

such that the left dot (1) is coupled to the left bath, and right dot (2) to the right bath. We note that the damping parameter is related to the usual notation ${\mathrm{\Gamma }}_L\equiv \mathrm{i}({\mathrm{\Sigma }}_L^r-{\mathrm{\Sigma }}_L^a)$ by $\eta =({\mathrm{\Gamma }}_L{)}_{11}/2$, and similar for ${\mathrm{\Gamma }}_R$. With this treatment of the baths, it is completely equivalent to the usual replacement of $E\to E+\mathrm{i}\eta$ to the free electron Green’s function. The 2 by 2 matrix elements for the Green’s function are obtained from a matrix inverse $(E+\mathrm{i}\eta -H_C{)}^{-1}$ as

(162)$\begin{array}{l}{G}_{11}^r=G_{22}^r={\displaystyle \frac{E+\mathrm{i}\eta }{(E+\mathrm{i}\eta {)}^2-t^2}},\end{array}$

(163)$\begin{array}{l}{G}_{12}^r=G_{21}^r={\displaystyle \frac{-t}{(E+\mathrm{i}\eta {)}^2-t^2}}.\end{array}$

Having obtained the retarded Green’s function and clarified the bath self-energies, our next step is to use them to obtain the photon self-energy ${\mathrm{\Pi }}^{<}$ as it enters explicitly for the radiation power. Initially, the self-energy as defined by Eq. (123) is location dependent. But in calculating the total radiation, a long-wave approximation was used as the typical thermal or even optical wave lengths are much larger comparing to the atomic distances. This leads to a total ${\mathrm{\Pi }}^{<}$ which is only directional dependent as a 3 by 3 matrix. While the local current is expressed by the coupling matrix $M^{l\mu }$, summing over $l$ gives the velocity matrix times the electron charge, i.e.,

(164)$\begin{array}{l}{I}_{\mu }={\displaystyle \sum _l{I}_{l\mu }=(-e)c^{†}{V}^{\mu }c,\mu =x,y,z.}\end{array}$

Under the random phase approximation (RPA), we treat the electrons as free electrons not coupled to the field (but can be coupled to the baths through bath self-energy). Since they are free electrons, we can apply Wick’s theorem, taking care that the electron operators are anti-commuting under the contour order sign. A general term in $\mathrm{\Pi }(\tau ,{\tau }^{\prime })$ takes the form $⟨T_{\tau }{c}_1^{†}{c}_2{c}_3^{†}{c}_4⟩$. Equal-time decomposition gives a constant in time due to time translational invariance. This constant can be absorbed into a redefinition of $D$, thus it does not appear as a self-energy term. Alternatively, the constant results in a $\delta$ function in frequency domain which does not contribute to transport. For a normal metal, the particle non-conserving terms, $⟨T_{\tau }{c}_1^{†}{c}_3^{†}⟩$ or $⟨T_{\tau }{c}_2{c}_4⟩$, is zero. This left with only one possible combination, $c_2{c}_3^{†}$ and $c_1^{†}{c}_4$. After swapping into the correct order for the Green’s function, which produces a minus sign, we find that the contour ordered $\mathrm{\Pi }$ can be expressed in terms of the electron Green’s function as

(165)$\begin{array}{l}{\mathrm{\Pi }}_{\mu \nu }(\tau ,{\tau }^{\prime })=-i\hslash e^2\mathrm{T}\mathrm{r}[V^{\mu }G(\tau ,{\tau }^{\prime })V^{\nu }G({\tau }^{\prime },\tau )].\end{array}$

Here the trace is over the electron site space, and the contour ordered electron Green’s function is $G_{jk}(\tau ,{\tau }^{\prime })=⟨T_{\tau }{c}_j(\tau )c_k^{†}({\tau }^{\prime })⟩/(\mathrm{i}\hslash )$. The lesser component is then given by ${\mathrm{\Pi }}_{\mu \nu }^{<}(t)=-\mathrm{i}\hslash e^2\mathrm{T}\mathrm{r}[V^{\mu }{G}^{<}(t)V^{\nu }{G}^{>}(-t)]$, here time-translational invariance is assumed. Since the radiation power formula requires the frequency domain one, we Fourier transform into $\omega$, obtaining

(166)$\begin{array}{l}{\mathrm{\Pi }}_{\mu \nu }^{<}(\omega )=-\mathrm{i}{e}^2{\displaystyle {\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\mathrm{d}E}{2\pi }\mathrm{T}\mathrm{r}[V^{\mu }{G}^{<}(E)V^{\nu }{G}^{>}(E-\hslash \omega )].}\end{array}$

The retarded component ${\mathrm{\Pi }}^r$ can be computed similarly with the replacement $\cdot G^{<}\cdot G^{>}\to \cdot G^r\cdot G^{<}+\cdot G^{<}\cdot G^a$.

Given the Hamiltonian matrix element as $H_{jk}$, the velocity matrix elements can be constructed from it as

(167)$\begin{array}{l}{V}_{jk}^{\mu }={\displaystyle \frac{1}{\mathrm{i}\hslash }{H}_{jk}(R_j^{\mu }-R_k^{\mu }).}\end{array}$

On a continuum, the velocity is the rate of change of position, $\mathit{v}=\mathrm{d}\mathit{r}/\mathrm{d}t=[\mathit{r},H]/(\mathrm{i}\hslash )$. Taking the matrix element with the states where the position operator is diagonal, we obtain the above in a tight-binding representation. For the two-dot model, the electron can move only in one direction with a spacing $a$, call it $x$, then

(168)$\begin{array}{l}{V}^x=\left(\begin{array}{cc}0& -\mathrm{i}\frac{at}{\hslash }\\ \mathrm{i}\frac{at}{\hslash }& 0\end{array}\right).\end{array}$

Plugging this result into the general formula, Eq. (166), we have

(169)$\begin{array}{rl}{\mathrm{\Pi }}_{xx}^{<}(\omega )& =\mathrm{i}{e}^2{\left(\frac{at}{\hslash }\right)}^2{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\frac{\mathrm{d}E}{2\pi }[G_{12}^{<}(E)G_{12}^{>}(E-\hslash \omega )\\ & +G_{21}^{<}(E)G_{21}^{>}(E-\hslash \omega )-G_{11}^{<}(E)G_{22}^{>}(E-\hslash \omega )\\ & -G_{22}^{<}(E)G_{11}^{>}(E-\hslash \omega )].\end{array}$

The integral can be performed approximately if we use the fact that $\eta$ is small. In this limit, the spectrum of the system is two sharp peaks at $\pm t$. The lesser Green’s function is obtained from the Keldysh equation by

(170)$\begin{array}{l}{G}^{<}(E)=G^r(E)(\begin{array}{cc}2\mathrm{i}{f}_L\eta & 0\\ 0& 2\mathrm{i}{f}_R\eta \end{array})G^a(E).\end{array}$

Here $f_L$ and $f_R$ are the Fermi functions associated with the left and right bath. $G^{>}$ is obtained by replacing $f$ by $f-1$. The calculation becomes somewhat tedious, but in the limit $\eta \to 0$, it simplifies greatly. Skipping some details, the final result for power after integrating over $E$ and $\omega$ is

(171)$\begin{array}{l}I={\displaystyle \frac{1}{3\pi }}{\displaystyle \frac{{\mathrm{e}}^2{v}_0^2{t}^2}{{\hslash }^2{ϵ}_0{c}^3}}(f_L+f_R)(2-f_L^{\prime }-f_R^{\prime }),\end{array}$

here $v_0=at/\hslash$, $f_L\equiv f_L(t)=1/({\mathrm{e}}^{{\beta }_L(t-{\mu }_L)}+1)$ is the electron occupation at the excited state from the left bath, and $f_L^{\prime }\equiv f_L(-t)$ is the electron occupation at the ground state from the left bath. $f_R$ and $f_R^{\prime }$ are similarly defined. At high bias with ${\mu }_L\to +\mathrm{\infty }$ and ${\mu }_R\to -\mathrm{\infty }$, we have $f_L=f_L^{\prime }=1$ and $f_R=f_R^{\prime }=0$. In this limit, we see that the power obtained is identical to the dipole oscillator result if we identify the frequency by $\hslash \omega =2t$ and a dipole moment with $p=ea/2$. This value of dipole moment is consistent with the matrix element $⟨{\mathrm{\Psi }}_{+}|(-e)\hat{x}|{\mathrm{\Psi }}_{-}⟩$ between the excited and ground state of the position operator. The radiation power at high bias can also be written as $\hslash \omega /\tau$ where $1/\tau$ is the spontaneous decay rate as in the Weisskopf−Wigner theory [151].

At the beginning of this section, we mentioned that the energy emitted is supplied by the electron baths. This point can be made quantitatively. We denote the energy out of the left bath as $I_L$, and that of right bath as $I_R$. We can also think of the infinity as a bath, but it only absorbs energy, $I_{\mathrm{\infty }}=-I\le 0$. Conservation of energy means $I_L+I_R+I_{\mathrm{\infty }}=0$. This conservation can be checked explicitly. We can calculate the energies of the left and right baths from the Meir−Wingreen formula, Eq. (23). Here we must consider the coupling of electrons with the field, thus the electron Green’s function should be the one with the nonlinear self-energy due to photons, i.e. $G_n=G+G{\mathrm{\Sigma }}_n{G}_n$, where $G$ is the free electron one. The electron−photon couplings are weak, so we use again the lowest order expansion approximation, i.e., $G_n^{<}\approx G^{<}+G^r{\mathrm{\Sigma }}_n^{<}{G}^a$. We also note $\mathrm{T}\mathrm{r}(G^{>}{\mathrm{\Sigma }}^{<}-G^{<}{\mathrm{\Sigma }}^{>})=0$, where ${\mathrm{\Sigma }}^{>,<}$ is the sum of the left and right bath self-energy. We are not interested in the energies that come out of the left bath and go into the right bath. If we compute the total, $I_L+I_R$, the first term from the free electron $G^{<}$ vanishes, and we shall focus on the second term only. The nonlinear self-energy ${\mathrm{\Sigma }}_n^{<}$ is similar to the scalar photon theory case and is due to the Fock diagram. The Fock diagram needs $D^{<}$ which is obtained from the Keldysh equation. Omitting some calculation details, the energy that is lost to space from left bath is

(172)$\begin{array}{rl}{I}_L=& \frac{e^2{v}_0^2{t}^2}{6\pi {\hslash }^2{ϵ}_0{c}^3}[f_L(2-f_L^{\prime }-f_R^{\prime })\\ & +(1-f_L^{\prime })(f_L+f_R)].\end{array}$

$I_R$ is obtained by swapping $L\leftrightarrow R$, and the sum is equal to $I$ obtained earlier.

Part III Full theory in temporal gauge

In Part II, we discussed the electromagnetic interactions with electrons in the transverse gauge by the vector potential, and in Part I, we take care of the Coulomb interaction with the scalar potential. In principle, taking both together we have a full theory of the electrodynamics. The transverse gauge, also known as the Coulomb gauge, is a standard approach in condensed matter physics for ease of quantization. However, it is not the most economical, as we have to calculate the charge-charge correlation, current-current correlation, or possibly a cross correlation between charge and current as a 4 by 4 matrix for $\mathrm{\Pi }$. Since charge and current are related by the continuity equation, these correlation functions are related. In the standard fluctuational electrodynamics (FE) of Polder and von Hove type [6,7], it is usually formulated by the $\mathit{E}$ and $\mathit{B}$ fields which are gauge independent quantities. It is then possible to reformulate our vector-field based theory in a different gauge, known as the $\varphi =0$ gauge or temporal gauge [152,153], which is more directly related to the gauge independent FE theory. Since we have banished the scalar field, we only need to consider a general vector field $\mathit{A}$ without the transverse requirement, and we only need to compute the current-current correlation. The drawback of the $\varphi =0$ gauge is that we need to impose additional condition on the quantum states so that Gauss’s law $\mathrm{\nabla }\cdot \mathit{E}=\rho /{ϵ}_0$ is satisfied. This extra complication seems not to hinder our formulation as we usually never consider the states explicitly in an NEGF formulation.

11 Hamiltonian, etc.

In the $\varphi =0$ gauge, the Lagrangian is the same as in Eq. (107), except that the scalar field $\varphi$ is set to 0, and the transverse condition on $\mathit{A}$ is not imposed. The Hamiltonian, by a similar step, is the same with $\varphi =0$, i.e.,

(173)$\begin{array}{rl}\hat{H}=& \frac{{ϵ}_0}{2}\int \mathrm{d}V\left[{\left(\frac{\mathrm{\partial }\mathit{A}}{\mathrm{\partial }t}\right)}^2+c^2{\left(\mathrm{\nabla }\times \mathit{A}\right)}^2\right]\\ & +\sum _{j,l}{c}_j^{†}{H}_{jl}{c}_l\exp \left(-\mathrm{i}\frac{e}{\hslash }{\int }_l^j\mathit{A}\cdot \mathrm{d}\mathit{r}\right).\end{array}$

Note that the Hamiltonian, Eq. (108), is gauge invariant. Here we commit a gauge choice by setting $\varphi =0$. The gauge is not completely fixed as we can still make transformation $\mathit{A}\to \mathit{A}-\mathrm{\nabla }\chi$, $c_j\to \exp (\mathrm{i}e{\chi }_j/\hslash )c_j$ with a time-independent $\chi$. Since the vector field here is the full one, the canonical commutation relation uses the normal Dirac $\delta$ function instead of the transverse delta function,

(174)$\begin{array}{l}[A_{\mu }(\mathit{r}),{\mathrm{\Pi }}_{\nu }({\mathit{r}}^{\prime })]=\mathrm{i}\hslash {\delta }_{\mu \nu }\delta (\mathit{r}-{\mathit{r}}^{\prime }),\end{array}$

where the conjugate momentum ${\mathrm{\Pi }}_{\nu }={ϵ}_0{\dot{A}}_{\nu }=-{ϵ}_0{E}_{\nu }$. The Heisenberg equation of motion for the vector field is

(175)$\begin{array}{l}{v}^{-1}\mathit{A}=-{ϵ}_0\left[{\displaystyle \frac{{\mathrm{\partial }}^2\mathit{A}}{\mathrm{\partial }{t}^2}}+c^2\mathrm{\nabla }\times (\mathrm{\nabla }\times \mathit{A})\right]=-\mathit{j}.\end{array}$

Here we define $v^{-1}$ as the $3\times 3$ matrix in the directional index and as a differential operator defined by the middle term acting on $\mathit{A}$. It turns out that its inverse is the free Green’s function, i.e., $\mathit{A}=-v\mathit{j}$. The awkward minus sign is needed so that $v$ is consistent with the NEGF definition of the retarded Green’s function by the commutator, $\theta (t)⟨[A_{\mu }(\mathit{r},t),A_{\nu }(\mathbf{0},0)]⟩/(\mathrm{i}\hslash )$. Formally, the current density on the right-hand side is given by a functional derivative with respect to $\mathit{A}$ or a commutator of the conjugate variable of $\mathit{A}$ to the Peierls substitution term,

(176)$\begin{array}{rl}\mathit{j}& =-\frac{\delta }{\delta \mathit{A}}\sum _{jl}{c}_j^{†}{H}_{jl}{c}_l\exp \left(-\frac{\mathrm{i}e}{\hslash }{\int }_l^j\mathit{A}\cdot \mathrm{d}\mathit{r}\right)\\ & =\frac{{ϵ}_0}{\mathrm{i}\hslash }\left[\dot{\mathit{A}}(\mathit{r}),\sum _{jl}{c}_j^{†}{H}_{jl}{c}_l\exp \left(-\frac{\mathrm{i}e}{\hslash }{\int }_l^j\mathit{A}\cdot \mathrm{d}\mathit{r}\right)\right].\end{array}$

This all looks easy in the $\varphi =0$ gauge. However, the catch is that Gauss’s law is missing from the Lagrangian, so we have to impose it as an extra condition. The question is when and where. It can be shown that Gauss’s law, if it is satisfied at one particular time, is also satisfied at all times, i.e., $[\mathrm{\nabla }\cdot \mathit{E}-\rho /{ϵ}_0,\hat{H}]=0$. This means that we can impose the requirement as an initial condition, i.e., initial states. The physical states $\mathrm{\Psi }$ must be selected such that $(\mathrm{\nabla }\cdot \mathit{E}-\rho /{ϵ}_0)\mathrm{\Psi }=0$ [154]. Alternatively, we can impose the condition on the operator $\mathit{A}$. The problem with $\mathit{A}$ is that the Fock state space created by it is larger than physically allowed. In the following, in determining the free Green’s function below, we will apply Gauss’s law when solving for $\mathit{A}$.

11.1 Free field Green’s function

The solution to the Green’s function in the absence of matter depends on the choice of boundary conditions. In full space, the easiest way to obtain the retarded Green’s function in the temporal gauge is to solve the equation of motion, $v^{-1}\mathit{A}=-\mathit{j}$, in Fourier space. To incorporate Gauss’s law, we note $\mathrm{\nabla }\times (\mathrm{\nabla }\times \mathit{A})=\mathrm{\nabla }(\mathrm{\nabla }\cdot \mathit{A})-{\mathrm{\nabla }}^2\mathit{A}$. Using $\mathrm{\nabla }\cdot \mathit{E}=\rho /{ϵ}_0$, $\mathit{E}=-\mathrm{\partial }\mathit{A}/\mathrm{\partial }t$, together with charge conservation, $\mathrm{\partial }\rho /\mathrm{\partial }t+\mathrm{\nabla }\cdot \mathit{j}=0$, in frequency domain, the divergence of $\mathit{A}$ can be expressed as proportional to the divergence of the current. We move this term to the right-hand side, thus it becomes a source term of a wave equation. We obtain the frequency and wave-vector space equation with the replacement $\mathrm{\partial }/\mathrm{\partial }t\to -\mathrm{i}\omega$, $\mathrm{\nabla }\to \mathrm{i}\mathit{q}$, as

(177)$\begin{array}{l}\mathit{A}(\omega ,\mathit{q})=-{\displaystyle \frac{\overleftrightarrow{\mathit{U}}-\mathit{q}\mathit{q}/(\omega /c{)}^2}{{ϵ}_0({\omega }^2-c^2{q}^2)}}\cdot \mathit{j}(\omega ,\mathit{q}).\end{array}$

The coefficient in front of $\mathit{j}$ is the Green’s function in $(\omega ,\mathit{q})$ space. This is nearly the same as Eq. (128) for the transverse gauge except that the numerator is not a projector. To be qualified as retarded, we also need a displacement for the frequency, $\omega \to \omega +\mathrm{i}\eta$, with an infinitesimal positive $\eta$. The Green’s function in real space is obtained by inverse Fourier transform, as [143,155]

(178)$\begin{array}{rl}{v}^r(\mathit{r},\omega )=& -\frac{{\mathrm{e}}^{\mathrm{i}\frac{\omega }{c}r}}{4\pi {ϵ}_0{c}^{2}r}\{(\overleftrightarrow{\mathit{U}}-\hat{\mathit{R}}\hat{\mathit{R}})\\ & +[-\frac{1}{\mathrm{i}\frac{\omega }{c}r}+\frac{1}{{\left(\mathrm{i}\frac{\omega }{c}r\right)}^2}](\overleftrightarrow{\mathit{U}}-3\hat{\mathit{R}}\hat{\mathit{R}})\}.\end{array}$

Here $\overleftrightarrow{\mathit{U}}$ is the identity dyadic, $\hat{\mathit{R}}=\mathit{r}/r$ is the radial direction unit vector. This is almost the same as Eq. (127) for the Coulomb gauge, except that they differ at the last term at $1/r^3$. Our notation $v^r$ is related to the usual dyadic Green’s function by $v^r=-{\mu }_0\overleftrightarrow{\mathit{G}}$. The mixed representation is useful for planar geometry, for which we Fourier transform the variable in $z$ direction back to real space, with the result

(179)$\begin{array}{rl}& v^r(\omega ,{\mathit{q}}_{\perp },z,z^{\prime })=\frac{{\mathrm{e}}^{\mathrm{i}{q}_z|z-z^{\prime }|}}{2\mathrm{i}{q}_z{ϵ}_0{c}^2{a}^2}\\ & \times \left(\begin{array}{ccc}1-\frac{q_x^2}{(\omega /c{)}^2}& -\frac{q_x{q}_y}{(\omega /c{)}^2}& -s\frac{q_x{q}_z}{(\omega /c{)}^2}\\ -\frac{q_y{q}_x}{(\omega /c{)}^2}& 1-\frac{q_y^2}{(\omega /c{)}^2}& -s\frac{q_y{q}_z}{(\omega /c{)}^2}\\ -s\frac{q_z{q}_x}{(\omega /c{)}^2}& -s\frac{q_z{q}_y}{(\omega /c{)}^2}& 1-\frac{q_z^2}{(\omega /c{)}^2}\end{array}\right).\end{array}$

Here $s$ is the sign of $(z-z^{\prime })$, and $q_z=\sqrt{(\omega /c{)}^2-q_{\perp }^2}$ for the propagating mode when $\omega /c>q_{\perp }$ and $q_z=$ $\mathrm{i}\sqrt{q_{\perp }^2-(\omega /c{)}^2}$ for the evanescent mode when $\omega /c<q_{\perp }$, $q_{\perp }=|{\mathit{q}}_{\perp }|$. We assume ${\mathit{q}}_{\perp }=(q_x,q_y)$ takes a discrete set of values in the first Brillouin zone of a square lattice with lattice constant $a$.

12 A unified theory for energy, momentum and angular momentum transfer

12.1 Conservations of energy, momentum, and angular momentum

The three conservation laws, energy, momentum, and angular momentum, are due to the symmetries of time translation, space displacement and rotation. We have already discussed the equation describing the energy conservation through the Poynting theorem, namely, $\mathrm{\partial }u/\mathrm{\partial }t+\mathrm{\nabla }\cdot \mathit{S}=-\mathit{E}\cdot \mathit{j}$. We now consider the other two, momentum and angular momentum conservations. The momentum conservation is related to the momentum density of the field. The Poynting vector $\mathit{S}$ is the energy flux, i.e., the energy current per unit cross-section area. So, the magnitude of the Poynting vector is $S=uc$ where $c$ is the speed of light, since photons move with the speed of light. The relation between energy of a photon and momentum is $cp=ϵ$ since photon is a massless relativistic particle. This means that the momentum density is given by $u/c$, or in terms of the Poynting vector, is $S/c^2$, or in vector form, $\mathit{S}/c^2={ϵ}_0\mathit{E}\times \mathit{B}$. We compute the rate of change of the momentum density. With the help of Maxwell’s equations, we find [68]

(180)$\begin{array}{l}-{\displaystyle \frac{1}{c^2}}\frac{\mathrm{\partial }\mathit{S}}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \overleftrightarrow{\mathit{T}}=\rho \mathit{E}+\mathit{j}\times \mathit{B},\end{array}$

here $\overleftrightarrow{\mathit{T}}={ϵ}_0\mathit{E}\mathit{E}+\frac{1}{{\mu }_0}\mathit{B}\mathit{B}-u\overleftrightarrow{\mathit{U}}$ is Maxwell’s stress tensor, and the right-hand side is the Lorentz force applied to object per unit volume, $\mathit{f}$. The conservation equation associated with the angular momentum is similarly obtained by noting that $\mathit{l}=\mathit{r}\times \mathit{S}/c^2={ϵ}_0\mathit{r}\times (\mathit{E}\times \mathit{B})$ is the angular momentum density. Taking again the rate of change of $\mathit{l}$ and using Maxwell's equations, we have [143,156]

(181)$\begin{array}{l}-{\displaystyle \frac{\mathrm{\partial }\mathit{l}}{\mathrm{\partial }t}}-\mathrm{\nabla }\cdot (\overleftrightarrow{\mathit{T}}\times \mathit{r})=\mathit{r}\times \mathit{f},\end{array}$

where $\mathit{f}=\rho \mathit{E}+\mathit{j}\times \mathit{B}$. We can obtain the angular momentum conservation equation (181) from the momentum conservation, Eq. (180), by multiplying it by $\mathit{r}\times$, and using the fact that $\overleftrightarrow{\mathit{T}}$ is symmetric. As a result, we can “pull” the divergence operator $\mathrm{\nabla }\cdot$ out, i.e., $\mathit{r}\times (\mathrm{\nabla }\cdot \overleftrightarrow{\mathit{T}})=-\mathrm{\nabla }\cdot (\overleftrightarrow{\mathit{T}}\times \mathit{r})$. Note that when a tensor $\overleftrightarrow{\mathit{T}}$ is cross-multiplied by $\mathit{r}$, the result is still a tensor, with component $(\overleftrightarrow{\mathit{T}}\times \mathit{r}{)}_{ij}=\sum _{kl}{T}_{ik}{x}_l{ϵ}_{klj}$. When a tensor is dotted with another vector, the result is a vector, e.g., $(\mathrm{\nabla }\cdot \overleftrightarrow{\mathit{T}}{)}_i=\sum _j\mathrm{\partial }/\mathrm{\partial }{x}_j{T}_{ji}$. For nonsymmetric tensor, such as $\overleftrightarrow{\mathit{T}}\times \mathit{r}$, dot from the left is different from dot from the right. We have a symmetry $-\mathbf{\Sigma }\cdot (\overleftrightarrow{\mathit{T}}\times \mathit{r})=(\mathit{r}\times \overleftrightarrow{\mathit{T}})\cdot \mathbf{\Sigma }$.

From the conservation equations, we have the physical interpretation that the Poynting vector $\mathit{S}$ gives the energy flux, and the Maxwell stress tensor $\overleftrightarrow{\mathit{T}}$ integrated over an enclosing surface with outward norm is the force applied to the body, and $-\overleftrightarrow{\mathit{T}}\times \mathit{r}$ the torque applied to the body. Since we are mainly interested in steady-state average, the average of a rate of change of a finite quantity is zero. Since averaging and taking partial derivative commute, we find $⟨\mathrm{\partial }a/\mathrm{\partial }t⟩=\mathrm{\partial }⟨a⟩/\mathrm{\partial }t=0$. Using this result, together with the conservation equations, we can compute the energy, momentum, and angular momentum transfer in two ways — by integrating over the surface enclosing some body or by a volume integral of the expressions on the right-hand side of the equations.

These classical expressions are changed to quantum operators in a quantum theory, with symmetrized product if necessary, and also earlier in part II we argue for a normal order in order to remove the zero-point motion contribution. However, here we have a second thought for the normal ordering. The reason is that it is precisely the zero-point motion that gives rise to the Casimir force. If we continue using normal order, we would not be able to predict a Casimir force. In fact, we need to use a symmetric order, $\frac{1}{2}⟨(AB+BA)⟩$. This means that the Green’s function will be the Keldysh one of $D^K=D^{>}+D^{<}$ that enters the expressions for physical observables of the transport quantities. In order to show that this is not ad hoc and arbitrary, we will say that the symmetric order is fundamental, and show explicitly, for energy transport and also for non-confined geometries such as radiation to infinity, the zero-point motion contribution cancels by itself. This will be demonstrated in Section 13.1. With these considerations, the power emitted, the force applied, and torque applied to a body $\alpha$ enclosed with a surface are

(182)$\begin{array}{l}{I}_{\alpha }={\displaystyle {\propto }_{\mathbf{\Sigma }}\mathit{S}\cdot \mathrm{d}\mathbf{\Sigma }={\int }_V(-\mathit{E}\cdot \mathit{j})\mathrm{d}V,}\end{array}$

(183)$\begin{array}{l}{\mathit{F}}_{\alpha }={\displaystyle {\propto }_{\mathbf{\Sigma }}\overleftrightarrow{\mathit{T}}\cdot \mathrm{d}\mathbf{\Sigma }={\int }_V\mathit{f}\mathrm{d}V,}\end{array}$

(184)$\begin{array}{l}{\mathit{N}}_{\alpha }={\displaystyle {\propto }_{\mathbf{\Sigma }}\mathit{r}\times \overleftrightarrow{\mathit{T}}\cdot \mathrm{d}\mathbf{\Sigma }={\int }_V\mathit{r}\times \mathit{f}\mathrm{d}V.}\end{array}$

Here $V$ is the volume enclosed by a simply connected surface $\mathbf{\Sigma }$. We have used Gauss’s theorem to change the divergence over the volume into a surface integral with outward norm.

12.2 Expressing transport quantities by Green’s functions

We comment on the surface integral expressions first. We have already given a general expression for heat transfer in terms of the Green’s function $D$ by Eq. (149), except that now $D$ is defined by the full $\mathit{A}$ without imposing the transverse requirement. Also, since we have decided to use the symmetric operator order, $D^{<}$ there should be replaced by the Keldysh version $D^K/2=(D^{<}+D^{>})/2$. The rest remain. For the force and torque with Maxwell’s stress tensor, we can write down similar expressions, but it is messy and not particularly illuminating. However, the expressions simplify if we take the long distance far field limit by integrating over a sphere of radius $R$. Instead of trying to derive these formulas explicitly here, we will give the results as a special case of the more general Meir−Wingreen formula, after we have elaborated on the concept of bath at infinity.

For the volume integrals, we remove the explicit charge $\rho$ dependence in favor of the current $\mathit{j}$. Due to charge conservation, we can compute the charge from the current. The Lorentz force per unit volume is $\mathit{f}=\rho \mathit{E}+\mathit{j}\times \mathit{B}=-\rho \mathrm{\partial }\mathit{A}/\mathrm{\partial }t+\mathit{j}\times (\mathbf{\nabla }\times \mathit{A})$. In steady state, $⟨\mathrm{\partial }(\rho \mathit{A})/\mathrm{\partial }t⟩=⟨(\mathrm{\partial }\rho /\mathrm{\partial }t)\mathit{A}⟩+⟨\rho \mathrm{\partial }\mathit{A}/\mathrm{\partial }t⟩=0$. We can move the time derivative from the vector potential to charge with a minus sign. Using the continuity equation, $\mathrm{\partial }\rho /\mathrm{\partial }t=$ $-\mathbf{\nabla }\cdot \mathit{j}$, we can write the force as (valid after taking average)

(185)$\begin{array}{rl}\mathit{f}& =-(\mathbf{\nabla }\cdot \mathit{j})\mathit{A}+\mathit{j}\times (\mathbf{\nabla }\times \mathit{A})\\ & =-\sum _{\mu }{\mathrm{\partial }}_{\mu }(j_{\mu }\mathit{A})+\sum _{\nu }{j}_{\nu }\mathbf{\nabla }{A}_{\nu }.\end{array}$

Here we have used the triple cross-product formula and combined the charge term with one of the cross-product terms. The index $\mu$ or $\nu$ takes the Cartesian directions, $x$, $y$, or $z$. Note that the first term is a divergence. If we integrate over a volume large enough to enclose the object where at the surface there is no current, then the first term is zero. We can then use for the total force calculation inside the integral as

(186)$\begin{array}{l}\mathit{f}={\displaystyle \sum _{\nu }{j}_{\nu }\mathbf{\nabla }{A}_{\nu }.}\end{array}$

However, the first term in Eq. (185) is no longer a divergence in the torque calculation in $\mathit{r}\times \mathit{f}$. By an integration by parts in the space index $\mu$, we can move the derivative to $\mathit{r}$. The same argument that at the surface there is no electric current can be used to eliminate the surface contribution. This gives the integrand for the total torque as

(187)$\begin{array}{l}\mathit{r}\times \mathit{f}=\mathit{j}\times \mathit{A}+\sum _{\nu }{j}_{\nu }(\mathit{r}\times \mathbf{\nabla })A_{\nu }.\end{array}$

Here the first term is interpreted to be the spin part of the contribution to angular momentum transfer, as it is independent of a choice of the coordinate origin, while the second term is interpreted as due to orbital angular momentum [157].

Now we have transformed each of the integrands for the three transport quantities in terms of the current density $\mathit{j}$ and vector potential $\mathit{A}$. The rest of the steps can be dealt with in a unified way. First, we define a new Green’s function as an intermediate quality by

(188)$\begin{array}{l}{F}_{\mu \nu }(\mathit{r}\tau ;{\mathit{r}}^{\prime }{\tau }^{\prime })={\displaystyle \frac{1}{i\hslash }}⟨T_{\tau }{A}_{\mu }(\mathit{r},\tau )j_{\nu }({\mathit{r}}^{\prime },{\tau }^{\prime })⟩.\end{array}$

Using this Green’s function, which reflects the interaction between the field and matter, we can write

(189)$\begin{array}{l}{I}_{\alpha }=\mathrm{R}\mathrm{e}{\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\hslash \omega \mathrm{T}\mathrm{r}\left[F^K(\omega )\right],}\end{array}$

(190)$\begin{array}{l}{\mathit{F}}_{\alpha }=\mathrm{R}\mathrm{e}{\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\mathrm{i}\hslash \mathrm{T}\mathrm{r}\left[{\mathrm{\nabla }}_{\mathit{r}}{F}^K(\omega )\right],}\end{array}$

(191)$\begin{array}{l}{\mathit{N}}_{\alpha }=\mathrm{R}\mathrm{e}{\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\mathrm{T}\mathrm{r}[\mathrm{i}\hslash \mathit{r}\times {\mathrm{\nabla }}_{\mathit{r}}{F}^K(\omega )-\mathit{S}{F}^K(\omega )].}\end{array}$

Note that with the symmetric order of any two Hermitian operators, we have $\frac{1}{2}⟨AB+BA⟩=$ $\mathrm{R}\mathrm{e}{\int }_0^{\mathrm{\infty }}\mathrm{i}\hslash D_{AB}^K(\omega )\mathrm{d}\omega /(2\pi )$, as the correlation defined in time domain by $\mathrm{i}\hslash D_{AB}^K(t)=⟨A(t)B+BA(t)⟩$ is real. As a result, in frequency domain the real part is symmetric, and imaginary part antisymmetric. By integrating from $-\mathrm{\infty }$ to $+\mathrm{\infty }$ only the real part of $D_{AB}^K(\omega )$ survives. We have used the symmetry to write the integrals with the positive frequency only. In the above expressions, we should view the Keldysh component $F^K$ as a matrix indexed by both position $\mathit{r}$ as well as direction $\nu$. The trace is in the combined space, i.e., $\mathrm{T}\mathrm{r}[\cdots ]={\int }_{V_{\alpha }}\mathrm{d}V\sum _{\nu }\cdots$. The volume integral covers the object $\alpha$ only. The gradient operator ${\mathbf{\nabla }}_{\mathit{r}}$ acts on the first argument of $F^K=F^K(\mathit{r},{\mathit{r}}^{\prime },\omega )$ which is associated with the argument of the vector field $\mathit{A}$. The factor in the energy current, $\hslash \omega$, is due to the time derivative with respect to $\mathit{A}$, Fourier transformed to the frequency domain. The first term of the torque expression, $\hat{\mathit{L}}=\mathit{r}\times (\frac{\hslash }{\mathrm{i}}{\mathbf{\nabla }}_{\mathit{r}})$, is the single particle orbital angular momentum operator in the position space, while $S_{\nu \gamma }^{\mu }=(-\mathrm{i}\hslash ){ϵ}_{\mu \nu \gamma }$ is the spin operator in the Cartesian direction acting on $F^K$.

We need to connect our $F$ back to our earlier Green’s function of the field $D$ and the materials properties $\mathrm{\Pi }$. In fact, such a relation does exist, it is

(192)$\begin{array}{l}F=-D{\mathrm{\Pi }}_{\alpha }.\end{array}$

This equation should be best viewed as defined on the Keldysh contour, and is a convolution in space as well as Keldysh contour. In the next subsection, we will prove this result and point out an additivity assumption needed for its validity. Here, if this expression is assumed, then the Keldysh component is obtained by the Langreth rule, $-F^K=D^r{\mathrm{\Pi }}_{\alpha }^K+D^K{\mathrm{\Pi }}_{\alpha }^a$. With this, we obtain the Meir−Wingreen formulas for the transport quantities as [158]

(193)$\begin{array}{l}{I}_{\alpha }={\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }(-\hslash \omega )\mathrm{R}\mathrm{e}\mathrm{T}\mathrm{r}(D^r{\mathrm{\Pi }}_{\alpha }^K+D^K{\mathrm{\Pi }}_{\alpha }^a),}\end{array}$

(194)$\begin{array}{l}{\mathit{F}}_{\alpha }={\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\mathrm{R}\mathrm{e}\mathrm{T}\mathrm{r}[\hat{\mathit{p}}(D^r{\mathrm{\Pi }}_{\alpha }^K+D^K{\mathrm{\Pi }}_{\alpha }^a)],}\end{array}$

(195)$\begin{array}{l}{\mathit{N}}_{\alpha }={\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\mathrm{R}\mathrm{e}\mathrm{T}\mathrm{r}[\hat{\mathit{J}}(D^r{\mathrm{\Pi }}_{\alpha }^K+D^K{\mathrm{\Pi }}_{\alpha }^a)].}\end{array}$

Here $\hat{\mathit{p}}=\frac{\hslash }{\mathrm{i}}\mathbf{\nabla }$ is the momentum operator, and $\hat{\mathit{J}}=\mathit{r}\times \hat{\mathit{p}}+\mathit{S}$ is the total angular momentum operator.

12.3 Properties of the Meir−Wingreen formulas

In this subsection, we investigate some of the properties implied by the Meir−Wingreen formulas for energy, momentum, and angular momentum transfers. We first discuss the conservation laws. We assume $N$ well-separated objects with $\alpha =1$, 2, ···, $N$ localized in a bounded region, see Fig.11. In order to account for the loss of the transported quantities to infinity, it is necessary to introduce one more object $(N+1)$ as the “object” at infinity. The conservations of the physical observables are then

Fig.11 $N$ objects of arbitrary shapes, each experiencing energy radiation, force and torque exerted to. There is one more special object called $\mathrm{\infty }$ which is the space outside the sphere of radius $R$. The norm vector $\mathit{n}$ is a unit vector pointing outwards from the enclosing surface.

Full size|PPT slide

(196)$\begin{array}{l}{\displaystyle \sum _{\alpha =1}^{N+1}{I}_{\alpha }=0,\sum _{\alpha =1}^{N+1}{\mathit{F}}_{\alpha }=\mathbf{0},\sum _{\alpha =1}^{N+1}{\mathit{N}}_{\alpha }=\mathbf{0}.}\end{array}$

These equations are obviously true, as the volume integrals can also be obtained by surface integrals enclosing the objects. Each surface separating the objects is used twice with opposite sense of the “outward” norm. Since the $\alpha$ dependence is only through the self-energy ${\mathrm{\Pi }}_{\alpha }$, it is necessary to introduce also a self-energy for the infinity. With the total self-energy $\mathrm{\Pi }=\sum _{\alpha }^{N+1}{\mathrm{\Pi }}_{\alpha }$, the conservations are satisfied for the three quantities if

(197)$\begin{array}{l}{D}^r{\mathrm{\Pi }}^K+D^K{\mathrm{\Pi }}^a=D^r{\mathrm{\Pi }}^K(I+D^a{\mathrm{\Pi }}^a)=0.\end{array}$

Here to get the second expression we have used the Keldysh equation $D^K=D^r{\mathrm{\Pi }}^K{D}^a$. Since the factor $D^r{\mathrm{\Pi }}^K$ cannot be zero in general, we require that $I+D^a{\mathrm{\Pi }}^a=0$. This is indeed valid if we recall the Dyson equation is

(198)$\begin{array}{l}{D}^r=v^r+v^r{\displaystyle \sum _{\alpha }^N{\mathrm{\Pi }}_{\alpha }^r{D}^r,}\end{array}$

which we can also write as $(v^r{)}^{-1}{D}^r=I+\sum _{\alpha }^N{\mathrm{\Pi }}_{\alpha }^r{D}^r$. If we identify the self-energy at infinity (or of the environment) as the differential operator [52],

(199)$\begin{array}{l}{\mathrm{\Pi }}_{\mathrm{\infty }}^r=-(v^r{)}^{-1},\end{array}$

and move the $-{\mathrm{\Pi }}_{\mathrm{\infty }}^r{D}^r$ term to the right side, we find $I+{\mathrm{\Pi }}^r{D}^r=0$. Here ${\mathrm{\Pi }}^r$ is a sum from object 1 to $N+1$ with ${\mathrm{\Pi }}_{N+1}^r\equiv {\mathrm{\Pi }}_{\mathrm{\infty }}^r$. Taking the Hermitian conjugate of this equation, we obtain the needed identity. We note that in the Dyson equation, the object at infinity is only an absorbing boundary condition as quantities transported to infinity cannot come back. The self-energies in the Dyson equation on the right-hand side in Eq. (198) does not include an “object” at infinity if the equation is solved in an unbounded domain. It is useful to think of each of the objects, including the object at infinity as a bath, supplying energy, momentum, and angular momentum for the photon fields. This derivation presents to us an explicit expression for the bath at infinity as a differential operator. In a later section, we give an algebraic expression defined on the surface of a sphere of radius $R\to \mathrm{\infty }$.

We can give an alternative argument for the self-energy of the environment (bath at infinity) as $-(v^r{)}^{-1}$. We recall that the contour ordered Green’s function $D=v+v\mathrm{\Pi }D$ implies a pair of equations in real time, one of them is the Dyson equation, Eq. (198), the other is the Keldysh equation for which we can write in alternative forms,

(200)$\begin{array}{rl}{D}^K& =D^r{\mathrm{\Pi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}^K{D}^a+(I+D^r{\mathrm{\Pi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}^r)v^K(I+{\mathrm{\Pi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}^a{D}^a)\\ & =D^r{\mathrm{\Pi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}^K{D}^a+D^r(v^r{)}^{-1}{v}^K(v^a{)}^{-1}{D}^a\\ & =D^r\{{\mathrm{\Pi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}^K+(2N_{\mathrm{\infty }}+1)[(-v^r{)}^{-1}+(v^a{)}^{-1}]\}{D}^a.\end{array}$

Here we define the self-energy of the objects to be ${\mathrm{\Pi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}=\sum _{\alpha =1}^N{\mathrm{\Pi }}_{\alpha }$, excluding the bath at infinity. The first line above is a general form of the Keldysh equation [63] without evoking a special property of $v^K$. For isolated systems, the second term on the right is zero because in this case, an isolated system is nondissipative, satisfying $v^{-1}{v}^K=0$ where $v^{-1}$ is interpreted as the differential operator. We can use the Dyson equation to obtain the second line. It is clear if the objects are absent, ${\mathrm{\Pi }}_{\mathrm{o}\mathrm{b}\mathrm{j}}=0$, we get $D^K=v^K$, which is the Green’s function for the free field. Assuming the free field in the absence of the objects is thermal, i.e., satisfying the fluctuation-dissipation theorem at temperature $T_{\mathrm{\infty }}$, that is, $v^K=(2N_{\mathrm{\infty }}+1)(v^r-v^a)$, we obtain the last line, from which we see that $(-v^r{)}^{-1}$ serves as the retarded self-energy for the environment.

The Meir−Wingreen formulas are the most general ones where each of the objects could be in some arbitrary nonequilibrium state. For the energy formula, we have a more symmetric form by adding the Hermitian conjugate inside the trace and divided by 2, and then using the general relations, $D^r-D^a=D^{>}-D^{<}$, $D^K=D^{>}+D^{<}$, and similarly for the self-energies $\mathrm{\Pi }$, as $I_{\alpha }=$ $-{\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\hslash \omega \mathrm{T}\mathrm{r}(D^{>}{\mathrm{\Pi }}_{\alpha }^{<}-D^{<}{\mathrm{\Pi }}_{\alpha }^{>})$. This is the same as for the scalar photon version, Eq. (35), except that here $D$ and ${\mathrm{\Pi }}_{\alpha }$ are tensors with directional indices.

The next consequence of the Meir−Wingreen formula we discuss is a derivation of the Landauer−Bütticker formula for the energy transport when local equilibrium is valid. By local equilibrium we mean that each object has a version of the fluctuation-dissipation theorem for the self-energies,

(201)$\begin{array}{l}{\mathrm{\Pi }}_{\alpha }^K=(2N_{\alpha }+1)({\mathrm{\Pi }}_{\alpha }^r-{\mathrm{\Pi }}_{\alpha }^a),\alpha =1,2,\cdots ,N,N+1.\end{array}$

Here $N_{\alpha }=1/({\mathrm{e}}^{{\beta }_{\alpha }\hslash \omega }-1)$ is the Bose function for object $\alpha$ at a local temperature of $T_{\alpha }=1/(k_B{\beta }_{\alpha })$. Using this expression for ${\mathrm{\Pi }}_{\alpha }^K$ and the Keldysh equation $D^K=D^r(\sum _{\alpha }^{N+1}{\mathrm{\Pi }}_{\alpha }^K)D^a$, we can completely eliminate the Keldysh components in favor of the retarded Green’s functions and self-energies. Introducing ${\mathrm{\Gamma }}_{\alpha }=\mathrm{i}({\mathrm{\Pi }}_{\alpha }^r-{\mathrm{\Pi }}_{\alpha }^a)$ as the spectrum of the bath, adding the Hermitian conjugate of $D^r{\mathrm{\Pi }}_{\alpha }^K+D^K{\mathrm{\Pi }}_{\alpha }^a$ and divided by 2 to take care of the real part, using the symmetry $(D^K{)}^{†}=-D^K$ and similarly for the self-energies, and finally the identity $\mathrm{i}(D^r-D^a)=D^r(\sum _{\beta }^{N+1}{\mathrm{\Gamma }}_{\beta })D^a$ [see Eq. (28)], we obtain, after some straightforward algebra,

(202)$\begin{array}{l}{I}_{\alpha }={\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{2\pi }\hslash \omega \sum _{\beta =1}^{N+1}(N_{\alpha }-N_{\beta })\mathrm{T}\mathrm{r}(D^r{\mathrm{\Gamma }}_{\beta }{D}^a{\mathrm{\Gamma }}_{\alpha }).}\end{array}$

We remind the reader that here the trace means integrating over the volume and summing over the Cartesian directions, and the Green’s functions have the arguments and indices, e.g., $D_{\mu \nu }^r(\mathit{r},{\mathit{r}}^{\prime },\omega )$. We note originally the Bose function enters as $2N_{\alpha }+1=\coth ({\beta }_{\alpha }\hslash \omega /2)$ where the $1$ is the contribution from zero-point motion. However, in the final formula, the Bose function enters only as a difference, the 1s have been canceled out. Thus, for energy transport, with local thermal equilibrium, the zero-point motion never contributes to energy transport. The same cannot be said about force and angular momentum or if the system is not in local equilibrium. Due to the presence of the extra differential operator $\hat{\mathit{p}}$ or $\hat{\mathit{J}}$, a similar derivation fails to go through, thus there is no equivalent Landauer formula for the force and angular momentum unless the operator $\hat{O}$ in front of $D$ commutes with ${\mathrm{\Pi }}_{\alpha }$. If there were such a formula in the sense that it has a factor $N_{\alpha }-N_{\beta }$, we would not have Casimir forces.

The Caroli-like expressions for the transmission coefficients have been derived for energy transfer between objects modelled as dipoles [54,159], and also for fluctuational-surface-current formulation in Ref. [56]. If we define the multiple-bath transmission coefficients as $T_{\beta \alpha }=\mathrm{T}\mathrm{r}(D^r{\mathrm{\Gamma }}_{\beta }{D}^a{\mathrm{\Gamma }}_{\alpha })$, in general it is not symmetric with respect to a swap of the two baths unless there are in total only two baths. For systems that are reciprocal, i.e., ${\mathrm{\Pi }}_{\alpha }^{\mathrm{T}}={\mathrm{\Pi }}_{\alpha }$, we can show that we do have $T_{\beta \alpha }=T_{\alpha \beta }$. Unsymmetric transmission implies an energy current known as super-current between objects even at thermal equilibrium [160], but the total $I_{\alpha }$ out of object $\alpha$ is still zero. Symmetric or not, there is a sum rule [57,161], $\sum _{\beta }{T}_{\beta \alpha }=\sum _{\beta }{T}_{\alpha \beta }$, which is just a consequence of the total current conservation, $\sum _{\alpha }{I}_{\alpha }=0$.

12.4 Prove F = −DΠ_α

In proving the Meir−Wingreen formulas, we have used this relation, which means, in full form in space, contour time, and direction,

(203)$\begin{array}{rl}& F_{\mu \nu }(\mathit{r}\tau ;{\mathit{r}}^{\prime }{\tau }^{\prime })\\ & =-\int {\mathrm{d}}^3{\mathit{r}}^{″}\int \mathrm{d}{\tau }^{″}\sum _{\lambda }{D}_{\mu \lambda }(\mathit{r}\tau ;{\mathit{r}}^{″}{\tau }^{″}){\mathrm{\Pi }}_{\alpha ,\lambda \nu }({\mathit{r}}^{″}{\tau }^{″};{\mathit{r}}^{\prime }{\tau }^{\prime }),\end{array}$

which is a matrix multiplication in the direction indices, and convolution in space and contour time. Note that here the self-energy in continuum form is for the object $\alpha$ so that $F$ function is associated with a particular object, which we have suppressed the $\alpha$ subscript for notational simplicity. We can argue about this relation not so rigorously with linear response. Within the lowest order of random phase approximation, the electrons and photons are not coupled directly so we can take the averages in each space separately. Focusing on the electrons first, then the response of the electrons due to the internal vector field is the current $\mathit{j}=-{\mathrm{\Pi }}_{\alpha }\mathit{A}$, here we have not yet evaluated the average on the photon space, so $\mathit{j}$ and $\mathit{A}$ are still quantum operators. Again, this equation means convolution in contour time and space, and ${\mathrm{\Pi }}_{\alpha }$ is a $3\times 3$ tensor in directional index space. Let us be slightly more precise by writing $j(2)=$ $-\int d(3){\mathrm{\Pi }}_{\alpha }(2,3)A(3)$, here we use the abbreviation $(n)\equiv ({\mathit{r}}_n,{\tau }_n,{\mu }_n)$, and $\int d(3)\cdot$ means integrations over space, contour time, and summing over the directions. Since we still treat $j(2)$ as a quantum operator, we can put the linear response into the definition of the Green’s function $F(1,2)=⟨T_{\tau }A(1)j(2)⟩/(\mathrm{i}\hslash )=-\int d(3)⟨T_{\tau }A(1)\cdot$ ${\mathrm{\Pi }}_{\alpha }(2,3)A(3)⟩/(\mathrm{i}\hslash )$. We can pull ${\mathrm{\Pi }}_{\alpha }$ out from the average as it is just a number, so we have, using the definition of $D$, $F(1,2)=-\int d(3)D(1,3){\mathrm{\Pi }}_{\alpha }(2,3)$. As the (bosonic) contour ordered function is symmetric with respect to the permutation of their arguments by definition, ${\mathrm{\Pi }}_{\alpha }(2,3)={\mathrm{\Pi }}_{\alpha }(3,2)$, we obtain $F=-D{\mathrm{\Pi }}_{\alpha }$.

But there is a more rigorous proof of this relation. It is convenient to start with the total $F$, that is, the sum of the contributions of all the objects (not counting infinity as one) generated by the total current, $\mathit{j}=\sum _{\alpha }{\mathit{j}}_{\alpha }$. As a consequence of the Dyson equation, we can show that $F=-D\mathrm{\Pi }$ is an exact result. To demonstrate this, we note that the contour ordered Dyson equation can be written alternatively as $D=v+v\mathrm{\Pi }D=v+D\mathrm{\Pi }v$. Multiplying $v^{-1}$ from the right, we find $D{v}^{-1}=I+D\mathrm{\Pi }$. Here the contour time differentiation is with respect to the second argument ${\tau }^{\prime }$. We can compute the left-hand side explicitly and connect to $F$. We first note that the vector field $\mathit{A}$ and the total current $\mathit{j}$ are related at the Heisenberg operator level as $v^{-1}\mathit{A}=-\mathit{j}$ with $v^{-1}=$ $-{ϵ}_0({\mathrm{\partial }}^2/\mathrm{\partial }{\tau }^2+c^2\mathrm{\nabla }\times \mathrm{\nabla }\times )$ acting on $\mathit{A}$. Here we have generalized the real time $t$ to the Keldysh contour time $\tau$. Acting $v^{-1}$ on $D$ from right is the same as acting on the second $\mathit{A}$ in $D$. Without the contour ordering, we can immediately replace $\mathit{A}{v}^{-1}$ by $-\mathit{j}$, which gives $-F$. Note that $v^{-1}$ is symmetric in space indices so acting from right is the same as acting from left. The contour order operator introduces extra terms. We can express $D$ in terms of the step functions without the contour operator $T_{\tau }$ as $D(\tau ,{\tau }^{\prime })=\theta (\tau ,{\tau }^{\prime })⟨A(\tau )A({\tau }^{\prime })⟩/(\mathrm{i}\hslash )+\theta ({\tau }^{\prime },\tau )⟨A({\tau }^{\prime })A(\tau )⟩/(\mathrm{i}\hslash )$, here we have suppressed the space and direction arguments since they are irrelevant for the reasoning. The space differentiation operation by $\mathrm{\nabla }$ does not contribute extra terms as it straightly goes inside the contour order $T_{\tau }$. Taking the first derivative $\mathrm{\partial }/\mathrm{\partial }{\tau }^{\prime }$ to $D$ leads to an extra term of the form $-\frac{1}{\mathrm{i}\hslash }\delta (\tau ,{\tau }^{\prime })\cdot$ $⟨[A(\tau ),A({\tau }^{\prime })]⟩$. The delta function appears because $\mathrm{\partial }\theta (\tau ,{\tau }^{\prime })/\mathrm{\partial }{\tau }^{\prime }=-\delta (\tau ,{\tau }^{\prime })$, and $\mathrm{\partial }\theta ({\tau }^{\prime },\tau )/\mathrm{\partial }{\tau }^{\prime }=\delta (\tau ,{\tau }^{\prime })$. The commutator is now at equal time because of the delta-function factor, and thus is zero. If we continue to the second derivative ${\mathrm{\partial }}^2/\mathrm{\partial }{\tau }^{\prime 2}$ that appears in $v^{-1}$, we find by similar steps a new term $\delta (\tau ,{\tau }^{\prime })⟨[\dot{A},A]⟩$. Using the canonical commutation relation, Eq. (174), we find that this is precisely the identity $I$ in space, contour time, and directional index. Thus, we have $-F=D\mathrm{\Pi }$.

This is not yet the equation we are supposed to prove, which is $F_{\alpha }=-D{\mathrm{\Pi }}_{\alpha }$ for object $\alpha$. In order for $F$ to have additivity over each object, we require that $\mathrm{\Pi }$ is additive, i.e.,

(204)$\begin{array}{l}\mathrm{\Pi }={\displaystyle \sum _{\alpha }^N{\mathrm{\Pi }}_{\alpha }.}\end{array}$

This additivity can be made true by selecting those Feynman diagrams of $\mathrm{\Pi }$ with the right-most vertex associated with the second argument of $\mathrm{\Pi }$ with object $\alpha$ for ${\mathrm{\Pi }}_{\alpha }$, while the first argument has no restriction. The additivity helps us to separate the contributions for the transported quantities unambiguously from each object. If we insist that ${\mathrm{\Pi }}_{\alpha }$ means strictly from object $\alpha$, then additivity is certainly true at the RPA level of approximation, but it fails at higher order of approximation for $\mathrm{\Pi }$. A good example is the diagram shown in Fig.12, known as the Aslamazov−Larkin diagram [162]. The additivity of the self-energy breaks down at $O(e^6)$, where $e$ is the magnitude of electron charge. We note that the additivity is false at the $e^4$ order for $\chi$, while is at $e^6$ for $\mathrm{\Pi }$. This is the main reason that we prefer to work with the irreducible diagrams from $\mathrm{\Pi }$, instead of working with the current-current correlation $\chi$.

Fig.12 One of the higher order contributions to the self-energy $\mathrm{\Pi }$ in a diagrammatic expansion known as Aslamazov-Larkin diagram. It consists of 6 electron Green’s functions (solid line with arrow) and 2 bare interactions $v$ (dotted line). The electron line of a closed loop must belong to the same object, say, $\alpha$ and $\beta$ for the two loops, respectively. Due to the existence of the $v$ lines which introduce mutual interaction, it is impossible to write the self-energy as a sum over each object individually.

Full size|PPT slide

13 Baths at infinity

Since the electromagnetic field can propagate to infinity, any collection of finite objects will have some energy transmitted to infinity. In order to account for the conservation laws, we must also count the “object” at infinity. An exception is the Polder and von Hove geometry of two materials with a gap; here in this problem, there is no need to consider the “infinity” and the space are all occupied by the objects. We can treat the empty space for $|\mathit{r}|>R$ outside a big sphere of radius $R$ as the infinity. It has the property that any energy, or momentum, or angular momentum sent to infinity is absorbed, and never reflected back. Mathematically, when we solve the retarded Dyson equation, we can treat the surface of the sphere at $R$ as an absorbing boundary condition. Or if the problem is solved in the full space including that outside the sphere, we must seek for a decaying solution that goes to zero at $R\to \mathrm{\infty }$.

Any object in our formalism is represented by the self-energy $\mathrm{\Pi }$, including the object at infinity. Usually, we set the temperature at infinity to zero, but we can also ascribe a temperature at infinity. In this case, the finite objects are enclosed in a black-body cavity at temperature $T_{\mathrm{\infty }}$ from the environment. The self-energy of the environment can be expressed as a differential operator defined in the whole space, Eq. (199), ${\mathrm{\Pi }}_{\mathrm{\infty }}^r=-v^{-1}$. The question is, does it have an imaginary part? Since an infinite domain must be dissipative, it does. The effect of this differential operator is realized if it acts on the actual solution of the problem. Very often it is transformed so that it does not appear explicitly in the end. This is conceptually simple, but computationally not very useful. A more useful point of view is to think of the bath (the object) at infinity as defined precisely at the sphere $|\mathit{r}|=R$. In this way, the degrees of freedom of the bath are solely specified by the solid angle $\mathrm{\Omega }$.

Here in this section, we derive an expression for the bath at infinity defined on the sphere locally as [95,158]

(205)$\begin{array}{l}{\mathrm{\Pi }}_{\mathrm{\infty }}^r=-\mathrm{i}{ϵ}_{0}c\omega (\overleftrightarrow{\mathit{U}}-\hat{\mathit{R}}\hat{\mathit{R}}),\end{array}$

and discuss its consequence. To obtain this result, we consider a dust model [64] of the bath at infinity, by saying that the space outside $R$ is dissipative. The dissipation is obtained phenomenologically by setting $\omega \to \omega +\mathrm{i}\eta$ in the free Green’s function, with some small positive $\eta$. This is equivalent to give the space $|\mathit{r}|>R$ as a dielectric medium with a local dielectric function $ϵ\approx 1+2\mathrm{i}\eta /\omega$, or a constant, infinitesimally small conductivity $2\eta {ϵ}_0$. The free Green’s function is a good approximation if one of the space coordinate approaches infinity, even if there is matter present at a finite region. We assume the asymptotic form of spherical wave solution $D^r\sim {\mathrm{e}}^{\mathrm{i}(\omega +\mathrm{i}\eta )|\mathit{r}|/c}$. The effect of the dusts is to introduce a damping over the purely oscillatory solution. To leading order in large $r$, we have $v^{-1}{D}^r\approx -2\mathrm{i}\eta {ϵ}_0\omega D^r$. Comparing with the Dyson equation, $v^{-1}{D}^r=I+{\mathrm{\Pi }}^r{D}^r$, we can identify the prefactor as the dust self-energy in a local form. Note that the retarded and advanced self-energies or Green’s functions are related by Hermitian conjugate. In the dust picture for the bath at infinity, since the dissipative bath has been modeled explicitly, we can put all the objects, including the dusts, in a large finite box as an overall isolated system. In this case $-v^{-1}$ will no longer have an imaginary part and is not interpreted as the self-energy of the bath anymore. Thus, the meaning of the differential operator depends on the boundary conditions choosing.

The bath self-energy appears in the Meir-Wingreen formula as $\mathrm{T}\mathrm{r}[\hat{O}(D^r{\mathrm{\Pi }}_{\mathrm{\infty }}^K+D^K{\mathrm{\Pi }}_{\mathrm{\infty }}^a)]$ where $\hat{O}$ is an operator acting on the first argument of $D^r$ and $D^K$. If we ascribe a temperature to the bath at infinity, the Keldysh bath self-energy is related to the retarded/advanced version by the fluctuation-dissipation theorem, ${\mathrm{\Pi }}_{\mathrm{\infty }}^K=$ $(2N_{\mathrm{\infty }}+1)({\mathrm{\Pi }}_{\mathrm{\infty }}^r-{\mathrm{\Pi }}_{\mathrm{\infty }}^a)$. The trace here will be interpreted in two ways. Originally, the trace is supposed to be the volume integral over all space with ${\mathrm{\Pi }}_{\mathrm{\infty }}^r$ the differential operator $-v^{-1}$. Now we have an explicit formula for the self-energy over the volume outside the sphere due to the dusts. The field intensity decay produces a finite volume integral. Due to the spherical symmetry, we can work in polar coordinates. The integration over the volume is a solid angle integral and radial ${\int }_R^{\mathrm{\infty }}\mathrm{d}r{r}^2\cdots$. The exponential decay in intensity, $D^r{D}^a\sim {\mathrm{e}}^{-2\eta |\mathit{r}|/c}$, gives us a length $c/(2\eta )$ to multiply the volume expression of the self-energy. In the limit $\eta \to 0^{+}$, we obtain a finite answer for the solid angle expression independent of $\eta$, which can be interpreted as a new self-energy defined on the sphere of radius $R$ as given by Eq. (205). The extra transverse projector which has angular dependence takes care that electrodynamic field in the far field is transverse.

13.1 No zero-point motion contribution at infinity?

In this subsection, we give an argument based on the Meir−Wingreen formula that there is no zero-point motion contribution to transport at infinity. If this were not true, we might run into problem of a divergent contribution to energy transport, for example. This is equivalent to say, that we can replace the Keldysh Green’s functions and self-energies by the lesser components, $D^K\to 2D^{<}$ and ${\mathrm{\Pi }}_{\mathrm{\infty }}^K\to 2{\mathrm{\Pi }}_{\mathrm{\infty }}^{<}$ in the Meir−Wingreen formulas, Eqs. (193) to (195), when applied to the quantities transmitted to infinity. In order to show that this is true, we need an assumption that each object has local thermal equilibrium (perhaps this assumption can be relaxed). Applying the fluctuation-dissipation theorem, then the terms $2N_{\alpha }+1$ in the Keldysh Green’s functions or self-energies can be split as twice the lesser components plus extra terms. We need to check that the extra terms are zero,

(206)$\begin{array}{l}\mathrm{R}\mathrm{e}\mathrm{T}\mathrm{r}[\hat{O}(D^r({\mathrm{\Pi }}_{\mathrm{\infty }}^r-{\mathrm{\Pi }}_{\mathrm{\infty }}^a)+D^r({\mathrm{\Pi }}^r-{\mathrm{\Pi }}^a)D^a{\mathrm{\Pi }}_{\mathrm{\infty }}^a)]=0?\end{array}$

Here ${\mathrm{\Pi }}^r=({\mathrm{\Pi }}^a{)}^{†}$ without a subscript is the total self-energy (including the bath at infinity). Using the identity $D^r-D^a=D^r({\mathrm{\Pi }}^r-{\mathrm{\Pi }}^a)D^a$, the expression can be simplified as to check if the real part of $\mathrm{T}\mathrm{r}[\hat{O}(D^r{\mathrm{\Pi }}_{\mathrm{\infty }}^r-D^a{\mathrm{\Pi }}_{\mathrm{\infty }}^a)]$ is zero. For the case of $\hat{O}=-\hslash \omega$, the trace is imaginary since the two terms are related by complex conjugate. But the power $I_{\mathrm{\infty }}$ must be real, so the result must be zero. For the general situations when $\hat{O}=-\hslash \omega$, $\hat{\mathit{p}}$, or $\hat{\mathit{J}}$, we note that the three factors $\hat{O}$, $D^r$, and ${\mathrm{\Pi }}_{\mathrm{\infty }}^r=(-v{)}^{-1}$ are operators in the space $\mathit{r}$ and direction $\mu$. The operator $\hat{O}$ is Hermitian. Importantly, $\hat{O}$ commutes with ${\mathrm{\Pi }}_{\mathrm{\infty }}^r$, $[\hat{O},{\mathrm{\Pi }}_{\mathrm{\infty }}^r]=0$; this is a consequence of time translation, space displacement, and rotation symmetries of the free field, which we can check explicitly using the operator representation of ${\mathrm{\Pi }}_{\mathrm{\infty }}^r$. Using these properties of the operators and the cyclic permutation property of trace, we find again that the two terms are related by complex conjugate and the expression (206) has no real part. With these arguments, we can compute the physical observables at infinity as

(207)$\begin{array}{l}{I}_{\mathrm{\infty }}(\hat{O})=\mathrm{R}\mathrm{e}{\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{\pi }\mathrm{T}\mathrm{r}[\hat{O}(D^r{\mathrm{\Pi }}_{\mathrm{\infty }}^{<}+D^{<}{\mathrm{\Pi }}_{\mathrm{\infty }}^a)].}\end{array}$

Here $\hat{O}=-\hslash \omega$, $\hat{\mathit{p}}$, and $\hat{\mathit{J}}$ for energy, momentum, and angular momentum transfer, respectively. Due to our sign convention, $I_{\mathrm{\infty }}=I_{\mathrm{\infty }}(-\hslash \omega )$ is the energy out of infinity, i.e., $-I_{\mathrm{\infty }}$ is the energy absorbed at infinity, ${\mathit{F}}_{\mathrm{\infty }}$ is the force applied to infinity (momentum transferred to infinity), and similarly ${\mathit{N}}_{\mathrm{\infty }}$ is the torque applied to infinity. If the temperature at infinity is set to 0, ${\mathrm{\Pi }}_{\mathrm{\infty }}^{<}=0$, only the second $D^{<}{\mathrm{\Pi }}_{\mathrm{\infty }}^a$ term contributes. We obtain for the energy transferred to infinity per second (power) by the objects from a finite region as,

(208)$\begin{array}{l}{I}_{\mathrm{\infty }}={\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{\pi }\mathrm{R}\mathrm{e}\mathrm{T}\mathrm{r}[(-\hslash \omega )D^r\sum _{\alpha =1}^N{\mathrm{\Pi }}_{\alpha }^{<}{D}^a{\mathrm{\Pi }}_{\mathrm{\infty }}^a].}\end{array}$

We have already derived the same formula using the Poynting vector expectation value and normal order on a sphere of radius $R\to \mathrm{\infty }$, see Eq. (150). Here $\mathrm{T}\mathrm{r}[\cdots ]=\int \mathrm{d}\mathrm{\Omega }{R}^2\sum _{\mu }\cdots$, i.e., integrating over the surface of the sphere and summing over the direction index. We obtain the formulas for force and torque similarly, just by replacing the operator $(-\hslash \omega )$ by the momentum or angular momentum operator, which can be shown to agree with a direct calculation using Maxwell’s stress tensor on the surface of a large sphere.

14 A consistency check with Krüger et al. theory

In the literature, the susceptibility $\chi$ is used as a primary material property instead of the self-energy (or polarizability) $\mathrm{\Pi }$. Here $\chi$ is the charge-charge correlation for the scalar field theory, and is the current-current correlation for the full theory. The self-energy $\mathrm{\Pi }$ is obtained from $\chi$ by keeping only the irreducible diagrams with respect to the interaction $v$ in a diagrammatic expansion. The two are related by $\chi =\mathrm{\Pi }+\mathrm{\Pi }v\mathrm{\Pi }+\cdots =\mathrm{\Pi }+v\mathrm{\Pi }\chi$. This gives rise to alternative expressions for the same physical quantities. For the scalar field result of the Landauer formula, Eq. (10), if we assume that the system is reciprocal [48], i.e., both $\mathrm{\Pi }$ and $\chi$ are symmetric, we can transform the Caroli formula into a form given by Yu et al. [73]. The reciprocity allows us to use ${\mathrm{\Pi }}^r-{\mathrm{\Pi }}^a=2\mathrm{i}\mathrm{I}\mathrm{m}{\mathrm{\Pi }}^r$ needed in the transformation [70]. Note that the assumed reciprocal relation, $\mathrm{\Pi }={\mathrm{\Pi }}^{\mathrm{T}}$, does not hold when the electrons are in a magnetic field.

Similarly, in pioneering works, Krüger et al. [10,18,52] give a series of formulas for various cases of energy transfer and (Casimir) forces for arbitrary objects. Again, if we assume reciprocity, the Meir−Wingreen formulas can also be transformed into the Krüger et al. form. It seems the assumed reciprocity, ${\mathrm{\Pi }}^{\mathrm{T}}=\mathrm{\Pi }$, is important for this equivalence to hold. The symmetry of $\mathrm{\Pi }$ implies the symmetry of $\chi$ as the free Green’s function is symmetric. For the case of current−current correlation, the transpose here is in the combined space of coordinate and directional index, i.e., $\mathrm{\Pi }={\mathrm{\Pi }}^{\mathrm{T}}$ means

(209)$\begin{array}{l}{\mathrm{\Pi }}_{\mu \nu }(\mathit{r},{\mathit{r}}^{\prime },\omega )={\mathrm{\Pi }}_{\nu \mu }({\mathit{r}}^{\prime },\mathit{r},\omega ).\end{array}$

Since the retarded $\mathrm{\Pi }$ is related to the dielectric function by ${\mathrm{\Pi }}^r=-{ϵ}_0{\omega }^2(ϵ-1)$, in a local theory, reciprocity means ${ϵ}^T=ϵ$ as a 3 by 3 matrix, which is the original meaning of “reciprocity” used by Lorentz [163].

As an illustration, let us consider a simple case of one object at a temperature $T$ emitting energy to infinity which is at zero temperature. Our result is given by Eq. (208). To transform into Krüger et al. form of Eq. (37), the key term in our formula is $\mathrm{T}\mathrm{r}(D^{<}{\mathrm{\Pi }}_{\mathrm{\infty }}^a)$. The trace here is interpreted as a volume integral over the whole space and sum over the directional index. The retarded Green’s function is $D=v+v\chi v$, here all of them are retarded version. We omitted the superscript $r$ for simplicity. The lesser component of $D$ can be expressed in terms of the retarded Green’s function by the Keldysh equation $D^{<}=D^r{\mathrm{\Pi }}^{<}{D}^a=ND(\mathrm{\Pi }-{\mathrm{\Pi }}^{†})D^{†}$. We have used the fluctuation-dissipation theorem for the object, where $N$ is the Bose function at temperature $T$; bath at infinity is at zero temperature. We now need a relation between the imaginary part of $\mathrm{\Pi }$ and the imaginary part of $\chi$. Using the relation $\chi =\mathrm{\Pi }+\mathrm{\Pi }v\mathrm{\Pi }+\cdots$, we can check that

(210)$\begin{array}{l}(1+\chi v)(\mathrm{\Pi }-{\mathrm{\Pi }}^{†})(1+v^{†}{\chi }^{†})=\chi -{\chi }^{†}+\chi (v^{†}-v){\chi }^{†}.\end{array}$

With this identity, we have $D^{<}{\mathrm{\Pi }}_{\mathrm{\infty }}^a=Nv(1+\chi v)\cdot$ $(\mathrm{\Pi }-{\mathrm{\Pi }}^{†})(1+v^{†}{\chi }^{†})v^{†}(-v^{-1}{)}^{†}=-Nv((\chi -{\chi }^{†})+\chi (v^{†}-v){\chi }^{†})$. We have used the representation for the bath at infinity by ${\mathrm{\Pi }}_{\mathrm{\infty }}^a=(-v^{-1}{)}^{†}$ which gets canceled by a factor $v^{†}$ from the $D^{†}$ term on the right side of the Keldysh equation. Finally, putting the integral over frequency back, using the reciprocity so that $\chi -{\chi }^{†}=2\mathrm{i}\mathrm{I}\mathrm{m}\chi$, and similarly for $v$, we get

(211)$\begin{array}{l}-I_{\mathrm{\infty }}={\displaystyle {\int }_0^{\mathrm{\infty }}\frac{\mathrm{d}\omega }{\pi }2\hslash \omega N\mathrm{I}\mathrm{m}\mathrm{T}\mathrm{r}[v\mathrm{I}\mathrm{m}\chi -v\chi (\mathrm{I}\mathrm{m}v){\chi }^{\ast }].}\end{array}$

This is Eq. (37) in Krüger et al. if we note a slight notational change, $\chi =-\mathcal{T}/{\mu }_0$, where $\mathcal{T}$ is the scattering operator, and $v=-{\mu }_0{\mathcal{G}}_0$ the free Green’s function, with the product $v\chi ={\mathcal{G}}_0\mathcal{T}$ remaining the same. Here ${\mu }_0$ is the vacuum permeability.

To show the equivalence of the Landauer/Caroli formula with Krüger et al. formulation with heat transfer or force between two objects, 1 and 2, we need to solve the Dyson equation in a block matrix form as we need only part of the matrix elements of $D^r$ connecting the two objects, as ${\mathrm{\Pi }}^r$ is block diagonal. By doing this, we can express $D_{12}^r$ in terms of ${\chi }_1$ and ${\chi }_2$ of individual objects together with the free field $v_{12}$ connecting the two objects. Similar steps show indeed they are equivalent, e.g., Eq. (57) in Ref. [52] and also is equivalent to Yu et al. form when $\mathrm{I}\mathrm{m}v=0$. A recent result for torque by Strekha et al. [164] is believed to agree with our formula also.

15 Breaking reciprocity by current drive – far field results

In this last section, we apply the above developed general theory for the transport of energy, momentum, and angular momentum due to a driven current in a nanoscale piece of metal or semiconductor, more specifically, a segment of nanoribbon of graphene as an example, see Fig.13. The driven current is realized by applying a bias voltage to the two leads connecting the piece, typical of mesoscopic ballistic transport setup.

Fig.13 Illustration of bias-induced emissions of energy, momentum, and angular momentum from the conducting channel of a two-terminal transport device of a graphene strip. The length of the central channel is $l$. The semi-infinite left and right leads are extensions of the pristine graphene strip.

Full size|PPT slide

In applying the radiation-to-infinity formulas, Eq. (208), for practical calculations, we need to make two more approximations. First, we replace the full Green’s function, $D$, by the free field one, $v$. This means that we ignore the multiple reflections inside matter, and take only the first term in the Dyson expansion, $D=v+v\mathrm{\Pi }v+\cdots$. The correction to $v$ is a kind of screening effect given by a factor $(1+v\chi )={ϵ}^{-1}$, so that $D={ϵ}^{-1}v$. For the far field, the longitudinal component (corresponding to the scalar field $\varphi$ in a transverse gauge) decays to zero quickly, thus does not contribute to far-field radiation. The screening effect for the transverse components is much smaller. We can make an order of magnitude estimate of the omitted terms based on a dimensional analysis. The self-energy is $\mathrm{\Pi }=$ $-{\omega }^2{ϵ}_0{a}^3(ϵ-1)$ on a lattice in a discrete representation, where $a$ is lattice constant and $ϵ$ is a dimensionless dielectric function. The free Green’s function is $v\sim 1/(4\pi {ϵ}_0{c}^{2}R)$; we take the distance $R\sim a$. Multiplying the two, we find $v\mathrm{\Pi }\approx v\chi \sim (\omega a/c{)}^2$, which is a ratio of two velocities. If we take the velocity $\omega a$ to be the electron speed $at/\hslash$ where $t$ is a hopping parameter of order eV, we find $\mathrm{\Pi }v\sim {10}^{-4}$. If we take $\hslash \omega$ to be the thermal energy of order $k_{B}T$, the result is even smaller. Using the approximation $D\approx v$ works for small molecules or graphene sheet where most of the atoms are exposed to space, the emission and then re-absorption are negligible. But inside the bulk of a material, even though the omitted single scattering term is small, the cumulative effect is large, so $D\approx v$ is not a good approximation. In such a case, we need to solve the Dyson equation $D=v+v\mathrm{\Pi }D$ by some means.

The second approximation is to use multipole expansion for the Green’s function. Since the object is finite, and the observation is at $R\to \mathrm{\infty }$, multipole expansion is a good approximation when the wavelength of the field is much larger than the relevant size of the object. Multipole expansion simply means we Taylor expand the second argument in

(212)$\begin{array}{l}{D}_{\mu \nu }(\mathit{R},\mathit{r})=D_{\mu \nu }(\mathit{R},\mathbf{0})+\mathit{r}\cdot ({\mathrm{\nabla }}_{{\mathit{r}}^{\prime }}{D}_{\mu \nu }(\mathit{R},{\mathit{r}}^{\prime })){|}_{{\mathit{r}}^{\prime }=\mathbf{0}}+\cdots .\end{array}$

We have omitted the frequency $\omega$ argument for simplicity. After replacing $D$ by $v$, the function is translationally invariant, which depends only on the difference $|\mathit{R}-\mathit{r}|$. So, we can also write as $v(\mathit{R},\mathit{r})=v(\mathit{R})-$ $\mathit{r}\cdot {\mathrm{\nabla }}_{\mathit{R}}v(\mathit{R})+\cdots$. It is sufficient to keep to linear order in $\mathit{r}$ for our purpose.

We now discuss the evaluation of the Keldysh term $D^r{\mathrm{\Pi }}^{<}{D}^a$ at far field and the solid angle integration. For the energy, force, and the spin part of the torque contribution, the second $(\overleftrightarrow{\mathit{U}}-3\hat{\mathit{R}}\hat{\mathit{R}})$ term in Eq. (178) can be omitted as it decays to zero faster than $1/R$. Only the first $1/R$ transverse term contributes to the far field. The $1/R^2$ factor from the product of $D^r$ and $D^a$ is canceled by the surface area element $R^2\mathrm{d}\mathrm{\Omega }$, picking up a finite result in the limit $R\to \mathrm{\infty }$. In this limit, we can replace the gradient operation ${\mathrm{\nabla }}_{\mathit{R}}v$ by $\mathrm{i}(\omega \hat{\mathit{R}}/c)v$. We have already given the monopole term after the solid angle integration by Eq. (155). The linear terms in $\mathit{r}$ do not contribute to the total energy radiation, because it contains an odd number of $\hat{\mathit{R}}$. The next non-vanishing terms appear in quadruple form, $\mathit{r}\mathit{r}{\mathrm{\Pi }}^{<}$. We can also give an order of magnitude estimate of the extra correction terms. Since each time we take a gradient, we obtain a dimensionless factor $\mathit{r}\cdot \hat{\mathit{R}}\omega /c$. So, the correction terms are smaller by a factor $(a\omega /c{)}^2$ than the monopole term. Here we agree that the typical distance is the lattice constant. This is at a similar order of magnitude as due to the screening corrections. Since $\omega /c\sim 1/\lambda$ is the inverse wavelength, the higher order terms are smaller by a factor $(a/\lambda {)}^2$.

For the force, if we keep only the monopole term, the result is zero as it contains an average of an odd number of the unit vectors $\hat{\mathit{R}}$. In this case, the dipole term is essential in order to have a nonzero value. We need to compute the average of four $\hat{\mathit{R}}$’s for the solid angle integration. We can check explicitly

(213)$\begin{array}{l}\overline{{\hat{R}}_{\alpha }{\hat{R}}_{\beta }{\hat{R}}_{\gamma }{\hat{R}}_{\mu }}={\displaystyle \frac{1}{15}}({\delta }_{\alpha \beta }{\delta }_{\gamma \mu }+{\delta }_{\alpha \gamma }{\delta }_{\beta \mu }+{\delta }_{\alpha \mu }{\delta }_{\beta \gamma }).\end{array}$

Here the overline means average over the solid angles, $\frac{1}{4\pi }\int \mathrm{d}\mathrm{\Omega }\cdots$. The Greek subscripts indicate the Cartesian components of the unit vectors. We also recall $\overline{{\hat{R}}_{\alpha }{\hat{R}}_{\beta }}=\frac{1}{3}{\delta }_{\alpha \beta }$. After some algebra, the force is [158]

(214)$\begin{array}{rl}{F}_{\mathrm{\infty }}^{\mu }=& {\int }_0^{\mathrm{\infty }}\mathrm{d}\omega \frac{\hslash {\omega }^3}{60{\pi }^2{ϵ}_0{c}^5}\sum _{\alpha ,l,l^{\prime }}[4{\mathrm{\Pi }}_{l\alpha ,l^{\prime }\alpha }^{<}(r_{\mu }^l-r_{\mu }^{l^{\prime }})\\ & -(r_{\alpha }^l-r_{\alpha }^{l^{\prime }}){\mathrm{\Pi }}_{l\alpha ,l^{\prime }\mu }^{<}-{\mathrm{\Pi }}_{l\mu ,l^{\prime }\alpha }^{<}(r_{\alpha }^l-r_{\alpha }^{l^{\prime }})].\end{array}$

For the monopole expressions, like the total energy emission, we only need the sum total of ${\mathrm{\Pi }}_{l\mu ,l^{\prime }\nu }$ over the sites $l$ and $l^{\prime }$. For the force, it is the first moment respect to the distance ${\mathit{r}}^l-{\mathit{r}}^{l^{\prime }}$ that is needed. For systems with reciprocity, $({\mathrm{\Pi }}^{<}{)}^T={\mathrm{\Pi }}^{<}$, this expression is zero. Thus, we need to break reciprocity in order to have a nonvanishing force.

For the orbital angular momentum contribution to the torque, since $\mathit{r}\times \hat{\mathit{p}}$ term is proportional to $R$ at large $\mathit{r}$, the leading $1/R^2$ term in the product of the Green’s functions is canceled to zero. We need to work harder to pick up the $1/R^3$ terms from ${\mathrm{\nabla }}_{\mathit{R}}{D}^r{D}^a$. In this case, we need to keep track of the next order term in the gradient operation. It turns out that the contribution from the spin part is equal to the contribution from the orbital angular momentum part (different from the conclusion in [165]), given a total that is in agreement with a direct calculation with Maxwell’s stress tensor result [148],

(215)$\begin{array}{l}{N}_{\mathrm{\infty }}^{\mu }={\displaystyle {\int }_0^{\mathrm{\infty }}\mathrm{d}\omega \frac{\hslash \omega }{6{\pi }^2{ϵ}_0{c}^3}\sum _{\alpha \beta }{ϵ}_{\mu \alpha \beta }{\mathrm{\Pi }}_{\beta \alpha }^{<}.}\end{array}$

Here ${ϵ}_{\mu \alpha \beta }$ is the anti-symmetric Levi−Civita symbol, ${\mathrm{\Pi }}_{\alpha \beta }^{<}$ without the site indices $(l,l^{\prime })$ means the sum total. Since $({\mathrm{\Pi }}^{<}{)}^{†}=-{\mathrm{\Pi }}^{<}$, there is no need to take the real part; the expression is real.

15.1 Graphene strip calculation

We consider the radiation of energy, as well as force and torque experienced by a finite piece of a zigzag graphene nano ribbon, as shown in Fig.13. The nonequilibrium state is modeled in the usual way of mesoscopic transport by connecting the piece with baths from the left and right. The two leads are in their respective equilibrium states, where the electrons follow the Fermi distribution $f_{L(R)}=1/\mathrm{e}\mathrm{x}\mathrm{p}[(E-{\mu }_{L(R)})/(k_B{T}_{L(R)})+1]$. Here, ${\mu }_{L(R)}$ and $T_{L(R)}$ are chemical potential and temperature of the left (right) lead, respectively. The C−C bond length is $a=1.4$ Å. In the numerical calculation, we choose the hopping parameter $t=2.7$ eV, the length of the central channel $l=5\sqrt{3}a$, and with $W=9$ (number of zigzag lines). We set the temperatures of the two leads to be equal, with $T_L=T_R=300$ K.

We show in Fig.14 the results of the force, power, and torque as a function of the chemical potentials of the two leads. Both the power and torque are symmetric with respect to ${\mu }_L=\pm {\mu }_R$, while the force is antisymmetric with ${\mu }_L=\pm {\mu }_R$. Fig.14(a) only shows the $x$ component of the force, as the $y$ component is very small, and the $z$ component is zero due to symmetry. The force on the strip due to the light emission is a nonequilibrium effect induced by the driving current in the $x$ direction, the direction of electron transport. The torque perpendicular to the surface in Fig.14(c) is due to the asymmetric structure in the central channel with width of $W=9$, while we have checked that the driving current cannot induce nonzero torque for mirror symmetric strip when the width is even, for example, $W=8$. As we can see from the figures that the force and torque are rather small even with large drive. It is a challenge to find realistic systems that give signals that experimentalist can detect.

Fig.14 Results of (a) force component in the $x$ direction, (b) power, and (c) torque component in the $z$ direction as a function of the chemical potentials ${\mu }_L/t$ and ${\mu }_R/t$ for the two-terminal graphene strip system.

Full size|PPT slide

Radiation of angular momentum from a benzene molecule driven out of equilibrium by two leads unsymmetrically attached to the molecule is presented in Ref. [148], while the emission of angular momentum from a two-dimensional Haldane model is calculated in Ref. [166]. In Ref. [158], the emissions of energy, momentum, and angular momentum from a semi-infinite graphene edge are calculated using the above formalism, taking into account translational invariance in the driven direction by going into wave-vector space. We note here, for the angular momentum it is truly an edge effect as the total is not proportional to the area but only length of the graphene. Also, it is truly a nonequilibrium effect, as in thermal equilibrium, both the momentum and angular momentum emissions (or the negative of force and torque applied to graphene) are zero. The nonequilibrium electron states are set up in a cheaper way by setting the Fermi functions based on the sign of the group velocity due to the left and right baths far away, which is valid in the ballistic regime. Such problem is beyond the ability of the usual fluctuational electrodynamics.

We setup a periodic boundary condition in the electron transport direction (open boundary condition in the width direction) [158]. With periodic boundary condition, the self-energy ${\mathrm{\Pi }}^{<}$ can be calculated in an eigenmode representation of the electrons. In the limit $\eta \to 0^{+}$, the lesser Green’s function can be represented in the eigenmode as

(216)$\begin{array}{l}{G}^{<}(E)=2\pi \mathrm{i}{\displaystyle \sum _n{f}_n|n⟩⟨n|\delta (E-E_n),}\end{array}$

where $H|n⟩=E_n|n⟩$ solves the single electron eigenvalue problem, and $f_n$ is the Fermi function in state $n$ with energy $E_n$. $G^{>}$ is obtained by a replacement of $f_n$ to $f_n-1$. Then ${\mathrm{\Pi }}^{<}(\omega )$ can be computed according to Eq. (166). This gives a more efficient method as the self-energy is a sum of delta functions, from which, the frequency integration can be performed analytically. In such a setup, the different baths are given implicitly based on the sign of the group velocity of electrons to the Fermi distribution.

The power emitted in this mode-space approximation is a Fermi-golden rule result [167],

(217)$\begin{array}{rl}-I_{\mathrm{\infty }}=& \frac{4\alpha }{3\hslash c^2}\sum _{\mu ,n,n^{\prime }}(E_n-E_{n^{\prime }}{)}^2\theta (E_n-E_{n^{\prime }})\\ & \times |⟨n|V^{\mu }|n^{\prime }⟩{|}^2{f}_n(1-f_{n^{\prime }}),\end{array}$

where $\alpha \approx 1/137$ is the fine-structure constant, $\theta$ is the step function, and $V^{\mu }$ is the velocity matrix in $\mu$ direction. We observe energy emission when an electron jumps from an occupied state to an empty state. Similar formulas are obtained for the torque and force with the help of the velocity matrices [158].

16 Conclusion

In this review, we presented the NEGF formalism for photon transport in the framework of Coulomb interactions among electrons in the non-retardation limit, as well as scalar and transverse vector potential formulations. In the near field, Coulomb interaction is the most important contribution, equivalent to keeping only $p$ polarization (TM mode) in a parallel plate geometry. For far-field radiation, we must also compute the vector potential contribution due to a finite speed of light, since only this term contributes to thermal radiation at infinity. To keep track of both the scalar and vector potentials in transverse gauge is a bit clumsy, thus, we also offer another choice of gauge which is the more economical temporal gauge. Not only the transmission of energy, but also of momentum and angular momentum, can be computed in a unified fashion with a Meir−Wingreen-like formula.

The solution to electrodynamics is formulated to solve the retarded Dyson equation; this is completely equivalent to solving the Lippmann−Schwinger equation in a scattering approach. The distribution or correlation function of the field is given by the Keldysh equation. This is equivalent to applying the fluctuation−dissipation theorem in fluctuational electrodynamics (FE). In FE, it is applied in the last step, but here we in some sense incorporate it upfront. The NEGF approach offers the option of not applying the local thermal equilibrium used in FE, thus opening the door for more general settings of nonequilibrium steady states. We have given several examples of this sort. The steady states are established through the applications of multiple electron baths at different temperatures or chemical potentials.

The Meir−Wingreen formula, or if local equilibrium is used for energy, the Landauer formula with the Caroli form for the transmission, is simple. But the complexity hides in the solution to the Dyson equation. In the Krüger et al. approach, one focuses on a single object at a time, and then combine the results to get the full solution, incorporating multiple reflections. This gives computational efficiency, but with a more involved formula for the transport quantity. Another advantage of treating each object separately as a scattering problem is that one can handle moving objects, by Lorentz boost. It is not clear how such problems can be handled in the NEGF framework. An individual treatment of each object by computing its scattering operator will fail to work if the property of one object is strongly influenced by nearby objects due to extreme proximity.

Our point of view here is to couple the calculation of the materials properties with the electrodynamic calculation closely so that the materials properties are calculated ab initio through the random phase approximation (RPA) for the electrons. One can also use the machinery of many-body physics to go beyond RPA, such as the self-consistent GW calculation or the contribution to the self-energy $\mathrm{\Pi }$ from the Aslamazov−Larkin diagram. Nonlinear effects in the field can be handled by including high order diagrams. These mutual correlations, not in the standard FE, may be important at extreme near field. Unlike the frequency-dependent local dielectric function, the photon self-energy $\mathrm{\Pi }$ is nonlocal to start with, which is a better physical quantity when the structure is described atomistically. In the examples shown, we only considered electrons as our materials, but phononic contributions can be, in principle, incorporated. In fact, there is a formula for the longitudinal inverse dielectric function due to phonons in a crystal given in Ref. [168] as

(218)$\begin{array}{rl}& {ϵ}^{-1}(\mathit{q},\omega )=\frac{1}{{ϵ}_{\mathrm{\infty }}}\\ & +\frac{e^2}{\mathrm{\Omega }{ϵ}_0{ϵ}_{\mathrm{\infty }}}\sum _{n,\kappa ,{\kappa }^{\prime }}\frac{Z_{\kappa }{Z}_{{\kappa }^{\prime }}}{\sqrt{M_{\kappa }{M}_{{\kappa }^{\prime }}}}\frac{\hat{\mathit{q}}\cdot {\mathit{e}}_{\kappa }(\mathit{q}n)\hat{\mathit{q}}\cdot {\mathit{e}}_{{\kappa }^{\prime }}^{\ast }(\mathit{q}n)}{(\omega +i\eta {)}^2-{\omega }_n^2(\mathit{q})},\end{array}$

where ${ϵ}_{\mathrm{\infty }}$ is the high-frequency dielectric constant from electrons, $Z_{\kappa }$ is the Born effective charge for the ion $\kappa$ in a unit cell, $\mathrm{\Omega }$ is unit cell volume, ${\mathit{e}}_{\kappa }(\mathit{q}n)$ is the phonon polarization vector of wavevector $\mathit{q}$ and phonon branch $n$, and ${\omega }_n(\mathit{q})$ is the phonon dispersion relation (note the relation ${ϵ}^{-1}=1+v\chi$). But the real challenge is to treat electrons, phonons, photons, and their mutual interactions together consistently.

References

[1] C.TannoudjiJ. Dupont-RocG.Grynberg, Photons and Atoms: Introduction to Quantum Electrodynamics, Wiley, 1989

[2] J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, A. Rubio. Strongly correlated electron–photon systems. Nature, 2022, 606(7912): 41

[3] M.Planck, The Theory of Heat Radiation, 2nd Ed., P. Blakiston’s Son & Co., Philadelphia, 1914

[4] C. M. Hargreaves. Anomalous radiative transfer between closely-spaced bodies. Phys. Lett. A, 1969, 30(9): 491

[5] G. A. Domoto, R. F. Boehm, C. L. Tien. Experimental investigation of radiative transfer between metallic surfaces at cryogenic temperatures. J. Heat Transfer, 1970, 92(3): 412

[6] D. Polder, M. van Hove. Theory of radiative heat transfer between closely spaced bodies. Phys. Rev. B, 1971, 4(10): 3303

[7] S.M. Rytov, Theory of Electric Fluctuations and Thermal Radiation, Air Force Cambridge Research Center, Bedford, MA, 1953

[8] S.M. RytovY. A. KravtsovV.I. Tatarskii, Principles of Statistical Radiophysics 3, Springer, Berlin, 1989

[9] H. B. Callen, T. A. Welton. Irreversibility and generalized noise. Phys. Rev., 1951, 83(1): 34

[10] M. Krüger, T. Emig, M. Kardar. Nonequilibrium Electromagnetic Fluctuations: Heat transfer and interactions. Phys. Rev. Lett., 2011, 106(21): 210404

[11] C. R. Otey, L. Zhu, S. Sandhu, S. Fan. Fluctuational electrodynamics calculations of near-field heat transfer in non-planar geometries: A brief overview. J. Quant. Spectrosc. Radiat. Transf., 2014, 132: 3

[12] G. Tang, L. Zhang, Y. Zhang, J. Chen, C. T. Chan. Near-field energy transfer between graphene and magneto−optic media. Phys. Rev. Lett., 2021, 127(24): 247401

[13] K. Joulain, J. P. Mulet, F. Marquier, R. Carminati, J. J. Greffet. Surface electromagnetic waves thermally excited: Radiative heat transfer, coherence properties and Casimir forces revisited in the near field. Surf. Sci. Rep., 2005, 57(3−4): 59

[14] S. Basu, Z. M. Zhang, C. J. Fu. Review of near-field thermal radiation and its application to energy conversion. Int. J. Energy Res., 2009, 33(13): 1203

[15] B. Song, A. Fiorino, E. Meyhofer, P. Reddy. Near-field radiative thermal transport: From theory to experiment. AIP Adv., 2015, 5(5): 053503

[16] A. I. Volokitin, B. N. J. Persson. Near-field radiative heat transfer and noncontact friction. Rev. Mod. Phys., 2007, 79(4): 1291

[17] S. A. Biehs, R. Messina, P. S. Venkataram, A. W. Rodriguez, J. C. Cuevas, P. Ben-Abdallah. Near-field radiative heat transfer in many-body systems. Rev. Mod. Phys., 2021, 93(2): 025009

[18] G. Bimonte, T. Emig, M. Kardar, M. Krüger. Nonequilibrium fluctuational quantum electrodynamics: Heat radiation, heat transfer, and force. Annu. Rev. Condens. Matter Phys., 2017, 8(1): 119

[19] C. Henkel. Nanoscale thermal transfer – An invitation to fluctuation electrodynamics. Zeitschrift für Naturforshchung A, 2017, 72(2): 99

[20] M.PascaleM. GiteauG.T. Papadakis, Perspective on near-field radiative heat transfer, arXiv: 2210.00929 (2022)

[21] A. Kittel, W. Müller-Hirsch, J. Parisi, S. A. Biehs, D. Reddig, M. Holthaus. Near-field heat transfer in a scanning thermal microscope. Phys. Rev. Lett., 2005, 95(22): 224301

[22] S. Shen, A. Narayanaswamy, G. Chen. Surface phonon polaritons mediated energy transfer between nanoscale gaps. Nano Lett., 2009, 9(8): 2909

[23] R. S. Ottens, V. Quetschke, S. Wise, A. A. Alemi, R. Lundock, G. Mueller, D. H. Reitze, D. B. Tanner, B. F. Whiting. Near-field radiative heat transfer between macroscopic planar surfaces. Phys. Rev. Lett., 2011, 107(1): 014301

[24] K. Kim, B. Song, V. Fernández-Hurtado, W. Lee, W. Jeong, L. Cui, D. Thompson, J. Feist, M. T. H. Reid, F. J. García-Vidal, J. C. Cuevas, E. Meyhofer, P. Reddy. Radiative heat transfer in the extreme near field. Nature, 2015, 528(7582): 387

[25] L. Cui, W. Jeong, 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., 2017, 8(1): 14479

[26] 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., 2017, 8(1): 14475

[27] T. Tokunaga, A. Jarzembski, T. Shiga, K. Park, M. Francoeur. Extreme near-field heat transfer between gold surfaces. Phys. Rev. B, 2021, 104(12): 125404

[28] V. Fernández-Hurtado, A. I. Fernández-Domínguez, J. Feist, F. J. García-Vidal, J. C. Cuevas. Super-Planckian far-field radiative heat transfer. Phys. Rev. B, 2018, 97(4): 045408

[29] J. C. Cuevas. Thermal radiation from subwavelength objects and the violation of Planck’s law. Nat. Commun., 2019, 10(1): 3342

[30] D. Thompson, L. Zhu, R. Mittapally, S. Sadat, Z. Xing, P. McArdle, M. Qazilbash, P. Reddy, E. Meyhofer. Hundred-fold enhancement in far-field radiative heat transfer over the blackbody limit. Nature, 2018, 561(7722): 216

[31] H. B. G. Casimir. On the attraction between two perfectly conducting plates. Proc. K. Ned. Akad. Wet., 1948, 51: 793

[32] E. M. Lifshitz. The theory of molecular attractive forces between solids. Sov. Phys. JETP, 1956, 2: 73

[33] P.H. G. M. van BloklandJ.T. G. Overbeek, van der Waals forces between objects covered with a chromium layer, J. Chem. Soc. Faraday Trans. I 74(0), 2637 (1978)

[34] S. K. Lamoreaux. Demonstration of the Casimir force in the 0.6 to 6 μm range. Phys. Rev. Lett., 1997, 78(1): 5

[35] U. Mohideen, A. Roy. Precision measurement of the Casimir force from 0.1 to 0.9 μm. Phys. Rev. Lett., 1998, 81(21): 4549

[36] J. M. Obrecht, R. J. Wild, M. Antezza, L. P. Pitaevskii, S. Stringari, E. A. Cornell. Measurement of the temperature dependence of the Casimir−Polder force. Phys. Rev. Lett., 2007, 98(6): 063201

[37] G. L. Klimchitskaya, U. Mohideen, V. M. Mostepanenko. The Casimir force between real materials: Experiment and theory. Rev. Mod. Phys., 2009, 81(4): 1827

[38] J. L. Garrett, D. A. T. Somers, J. N. Munday. Measurement of the Casimir force between two spheres. Phys. Rev. Lett., 2018, 120(4): 040401

[39] A. Stange, D. K. Campbell, D. J. Bishop. Science and technology of the Casimir effect. Phys. Today, 2021, 74(1): 42

[40] C. M. Wilson, G. Johansson, A. Pourkabirian, M. Simoen, J. R. Johansson, T. Duty, F. Nori, P. Delsing. Observation of the dynamical Casimir effect in a superconducting circuit. Nature, 2011, 479(7373): 376

[41] S. Vezzoli, A. Mussot, N. Westerberg, A. Kudlinski, H. Dinparasti Saleh, A. Prain, F. Biancalana, E. Lantz, D. Faccio. Optical analogue of the dynamical Casimir effect in a dispersion-oscillating fibre. Commun. Phys., 2019, 2(1): 84

[42] K. Y. Fong, H. K. Li, R. Zhao, S. Yang, Y. Wang, X. Zhang. Phonon heat transfer across a vacuum through quantum fluctuations. Nature, 2019, 576(7786): 243

[43] M. F. Maghrebi, A. V. Gorshkov, J. D. Sau. Fluctuation-induced torque on a topological insulator out of thermal equilibrium. Phys. Rev. Lett., 2019, 123(5): 055901

[44] M. Katoh, M. Fujimoto, H. Kawaguchi, K. Tsuchiya, K. Ohmi, T. Kaneyasu, Y. Taira, M. Hosaka, A. Mochihashi, Y. Takashima. Angular momentum of twisted radiation from an electron in spiral motion. Phys. Rev. Lett., 2017, 118(9): 094801

[45] X. Gao, C. Khandekar, Z. Jacob, T. Li. Thermal equilibrium spin torque: Near-field radiative angular momentum transfer in magneto−optical media. Phys. Rev. B, 2021, 103(12): 125424

[46] M. L. N. Chen, L. J. Jiang, W. E. I. Sha. Orbital angular momentum generation and detection by geometric-phase based metasurfaces. Appl. Sci. (Basel), 2018, 8(3): 362

[47] E. Nagali, F. Sciarrino, F. De Martini, L. Marrucci, B. Piccirillo, E. Karimi, E. Santamato. Quantum information transfer from spin to orbital angular momentum of photons. Phys. Rev. Lett., 2009, 103(1): 013601

[48] V. S. Asadchy, M. S. Mirmoosa, A. Dìaz-Rubio, S. Fan, S. A. Tretyakov. Tutorial on electromagnetic nonreciprocity and its origins. Proc. IEEE, 2020, 108(10): 1684

[49] C. Khandekar, S. Buddhiraju, P. R. Wilkinson, J. K. Gimzewski, A. W. Rodriguez, C. Chase, S. Fan. Nonequilibrium lateral force and torque by thermally excited nonreciprocal surface electromagnetic waves. Phys. Rev. B, 2021, 104(24): 245433

[50] R. Messina, M. Antezza. Casimir−Lifshitz force out of thermal equilibrium and heat transfer between arbitrary bodies. Europhys. Lett., 2011, 95(6): 61002

[51] R. Messina, M. Antezza. Scattering-matrix approach to Casimir−Lifshitz force and heat transfer out of thermal equilibrium between arbitrary bodies. Phys. Rev. A, 2011, 84(4): 042102

[52] M. Krüger, G. Bimonte, T. Emig, M. Kardar. Trace formulas for nonequilibrium Casimir interactions, heat radiation, and heat transfer for arbitrary objects. Phys. Rev. B, 2012, 86(11): 115423

[53] B. A. Lippmann, J. Schwinger. Variational principles for scattering processes I. Phys. Rev., 1950, 79(3): 469

[54] P. Ben-Abdallah, S. A. Biehs, K. Joulain. Many-body radiative heat transfer theory. Phys. Rev. Lett., 2011, 107(11): 114301

[55] A. W. Rodriguez, O. Ilic, P. Bermel, I. Celanovic, J. D. Joannopoulos, M. Soljačić, S. G. Johnson. Frequency-selective near-field radiative heat transfer between photonic crystal slabs: A computational approach for arbitrary geometries and materials. Phys. Rev. Lett., 2011, 107(11): 114302

[56] A.W. RodriguezM.T. H. ReidS.G. Johnson, Fluctuating-surface-current formulation of radiative heat transfer for arbitrary geometries, Phys. Rev. B 86, 220302(R) (2012)

[57] S.Datta, Electronic Transport in Mesoscopic Systems, Cambridge Univ. Press, 1995

[58] M.D. Ventra, Electrical Transport in Nanoscale Systems, Cambridge Univ. Press, 2008

[59] J. S. Wang, J. Wang, J. T. Lü. Quantum thermal transport in nanostructures. Eur. Phys. J. B, 2008, 62(4): 381

[60] Z.Z. YuG. H. XiongL.F. Zhang, A brief review of thermal transport in mesoscopic systems from nonequilibrium Green’s function approach, Front. Phys. 16(4), 43201 (2021)

[61] M. Janowicz, D. Reddig, M. Holthaus. Quantum approach to electromagnetic energy transfer between two dielectric bodies. Phys. Rev. A, 2003, 68(4): 043823

[62] U. Aeberhard. Theory and simulation of quantum photovoltaic devices based on the non-equilibrium Green’s function formalism. J. Comput. Electron., 2011, 10(4): 394

[63] H.HaugA. P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors, 2nd Ed., Springer-Verlag, 2008

[64] W. Eckhardt. Macroscopic theory of electromagnetic fluctuations and stationary radiative heat transfer. Phys. Rev. A, 1984, 29(4): 1991

[65] L. V. Keldysh. Diagram technique for nonequilibrium processes. Sov. Phys. JETP, 1965, 20: 1018

[66] J. S. Wang, B. K. Agarwalla, H. Li, J. Thingna. Nonequilibrium Green’s function method for quantum thermal transport. Front. Phys., 2014, 9(6): 673

[67] G. D. Mahan. Tunneling of heat between metals. Phys. Rev. B, 2017, 95(11): 115427

[68] J.D. Jackson, Classical Electrodynamics, 3rd Ed., John Wiley & Sons, 1999

[69] I. Smolić, B. Klajn. Capacitance matrix revisited. Prog. Electromagn. Res. B Pier B, 2021, 92: 1

[70] J. S. Wang, Z. Q. Zhang, J. T. Lü. Coulomb-force-mediated heat transfer in the near field: Geometric effect. Phys. Rev. E, 2018, 98(1): 012118

[71] R.KuboM. TodaN.Hashitsume, Statistical Physics II — Nonequilibrium Statistical Mechanics, 2nd Ed., Springer, 1991

[72] G.F. GiulianiG.Vignale, Quantum Theory of the Electron Liquid, Cambridge Univ. Press, 2005

[73] R. Yu, A. Manjavacas, F. J. García de Abajo. Ultrafast radiative heat transfer. Nat. Commun., 2017, 8(1): 2

[74] R. Landauer. Spatial Variation of Currents and Fields Due to Localized Scatterers in Metallic Conduction. IBM J. Res. Develop., 1957, 1(3): 223

[75] C. Caroli, R. Combescot, P. Nozieres, D. Saint-James. Direct calculation of the tunneling current. J. Phys. C, 1971, 4(8): 916

[76] J. S. Wang, J. Peng. Capacitor physics in ultra-near-field heat transfer. Europhys. Lett., 2017, 118(2): 24001

[77] J. H. Jiang, J. S. Wang. Caroli formalism in near-field heat transfer between parallel graphene sheets. Phys. Rev. B, 2017, 96(15): 155437

[78] T. Zhu, J. S. Wang. Generalized first-principles method to study near-field heat transfer mediated by Coulomb interaction. Phys. Rev. B, 2021, 104(12): L121409

[79] Y. Meir, N. S. Wingreen. Landauer formula for the current through an interacting electron region. Phys. Rev. Lett., 1992, 68(16): 2512

[80] A. P. Jauho, N. S. Wingreen, Y. Meir. Time-dependent transport in interacting and noninteracting resonant-tunneling systems. Phys. Rev. B, 1994, 50(8): 5528

[81] G.StefanucciR.van Leeuwen, Nonequilibrium Many-Body Theory of Quantum Systems, Cambridge Univ. Press, 2013

[82] D.C. Langreth, in: Linear and Nonlinear Electron Transport in Solids, NATO Advanced Study Institute Series, Vol. 17, edited by J. T. Devreese and V. E. van Doren, Springer, Boston, MA, 1976, p. 3

[83] J. T. Lü, J. S. Wang. Coupled electron and phonon transport in one-dimensional atomic junctions. Phys. Rev. B, 2007, 76(16): 165418

[84] D.BohmD. Pines, A collective description of electron interactions (III): Coulomb interactions in a degenerate electron gas, Phys. Rev. 92(3), 609 (1953)

[85] M.PaulssonT. FrederiksenM.Brandbyge, Modeling inelastic phonon scattering in atomic- and molecular-wire junctions, Phys. Rev. B 72, 201101(R) (2005)

[86] L. K. Dash, H. Ness, R. W. Godby. Nonequilibrium electronic structure of interacting single-molecule nanojunctions: Vertex corrections and polarization effects for the electron−vibron coupling. J. Chem. Phys., 2010, 132(10): 104113

[87] A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle Systems, McGraw-Hill, 1971

[88] G. W. Ford, M. Kac, P. Mazur. Statistical mechanics of assemblies of coupled oscillators. J. Math. Phys., 1965, 6(4): 504

[89] J. Peng and J. S. Wang, Current-induced heat transfer in double-layer graphene, arXiv: 1805.09493 (2019)

[90] Z. Q. Zhang, J. T. Lü, J. S. Wang. Energy transfer between two vacuum-gapped metal plates: Coulomb fluctuations and electron tunneling. Phys. Rev. B, 2018, 97(19): 195450

[91] M. Büttiker. Symmetry of electrical conduction. IBM J. Res. Develop., 1988, 32(3): 317

[92] L. Hedin. New method for calculating the one-particle Green’s function with application to the electron-gas problem. Phys. Rev., 1965, 139(3A): A796

[93] R. M. Martin, L. Reining, and D. M. Ceperley, Interacting Electrons, Cambridge Univ. Press, 2016

[94] G. D. Mahan, Many-Particle Physics, 3rd Ed., Kluwer Academic, 2000

[95] J. Peng, H. H. Yap, G. Zhang, and J. S. Wang, A scalar photon theory for near-field radiative heat transfer, arXiv: 1703.07113 (2017)

[96] S. Weinberg, The Quantum Theory of Fields, Volume 1: Foundations, Cambridge Univ. Press, 2005

[97] R. J. Glauber. Attenuators, and Schrödinger’s Cat. Ann. N. Y. Acad. Sci., 1986, 480(1): 336

[98] R. A. Jishi, Feynman Diagram Techniques in Condensed Matter Physics, Cambridge Univ. Press, 2013

[99] J. Rammer, Quantum Field Theory of Non-equilibrium States, Cambridge Univ. Press, 2007

[100] H. Bruus and K. Flensberg, Many-Body Quantum Theory in Condensed Matter Physics: An Introduction, Oxford Univ. Press, 2004

[101] S. Datta. Nanoscale device modeling: The Green’s function method. Superlattices Microstruct., 2000, 28(4): 253

[102] L. P. Kadanoff and G. Baym, Quantum Statistical Mechanics, W. A. Benjamin, Inc, 1962

[103] J. M. Ziman, Electrons and Phonons, Clarendon Press, Oxford, 1960

[104] B. van Duppen, A. Tomadin, A. N. Grigorenko, M. Polini. Current-induced birefringent absorption and non-reciprocal plasmons in grapheme. 2D Mater., 2016, 3: 015011

[105] D. Svintsov, V. Ryzhii. Comment on “Negative Landau damping in bilayer graphene”. Phys. Rev. Lett., 2019, 123(21): 219401

[106] T. A. Morgado, M. G. Silveirinha. Negative Landau damping in bilayer graphene. Phys. Rev. Lett., 2017, 119(13): 133901

[107] B. Shapiro. Fluctuation-induced forces in the presence of mobile carrier drift. Phys. Rev. B, 2017, 96(7): 075407

[108] O. Ilic, M. Jablan, J. D. Joannopoulos, I. Celanovic, H. Buljan, M. Soljačić. Near-field thermal radiation transfer controlled by plasmons in graphene. Phys. Rev. B, 2012, 85(15): 155422

[109] J. B. Pendry. Radiative exchange of heat between nanostructures. J. Phys.: Condens. Matter, 1999, 11(35): 6621

[110] F. Herz, C. Kathmann, S. A. Biehs. General trace formula for heat flux fluctuations. Europhys. Lett., 2020, 130(4): 44003

[111] J. L. Wise, N. Roubinowitz, W. Belzig, D. M. Basko. Signature of resonant modes in radiative heat current noise spectrum. Phys. Rev. B, 2022, 106(16): 165407

[112] J. S. Wang, B. K. Agarwalla, H. Li. Transient behavior of full counting statistics in thermal transport. Phys. Rev. B, 2011, 84(15): 153412

[113] G. Tang, J. S. Wang. Heat transfer statistics in extreme-near-field radiation. Phys. Rev. B, 2018, 98(12): 125401

[114] M. Campisi, P. Hänggi, P. Talkner. Colloquium: Quantum fluctuation relations: Foundations and applications. Rev. Mod. Phys., 2011, 83(3): 771

[115] B. K. Agarwalla, B. Li, J. S. Wang. Full-counting statistics of heat transport in harmonic junctions: Transient, steady states, and fluctuation theorems. Phys. Rev. E, 2012, 85(5): 051142

[116] L. S. Levitov, G. B. Lesovik. Charge distribution in quantum shot noise. JETP Lett., 1993, 58(3): 230

[117] G. Tang, H. H. Yap, J. Ren, J. S. Wang. Anomalous near-field heat transfer in carbon-based nanostructures with edge states. Phys. Rev. Appl., 2019, 11(3): 031004

[118] R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford Univ. Press, 1989

[119] N. W. Ashcroft and N. D. Mermin, Solid State Physics, Saunders College Publishing, 1976

[120] S. L. Adler. Quantum theory of the dielectric constant in real solids. Phys. Rev., 1962, 126(2): 413

[121] N. Wiser. Dielectric constant with local field effects included. Phys. Rev., 1963, 129(1): 62

[122] M. S. Hybertsen, S. G. Louie. Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies. Phys. Rev. B, 1986, 34(8): 5390

[123] J. Deslippe, G. Samsonidze, D. A. Strubbe, M. Jain, M. L. Cohen, S. G. Louie. BerkeleyGW: A massively parallel computer package for the calculation of the quasiparticle and optical properties of materials and nanostructures. Comput. Phys. Commun., 2012, 183(6): 1269

[124] F. Xuan, Y. Chen, S. Y. Quek. Quasiparticle levels at large interface systems from many-body perturbation theory: The XAF-GW method. J. Chem. Theory Comput., 2019, 15(6): 3824

[125] F. A. Rasmussen, First Principles Calculations of Electronic Excitations in 2D Materials, Ph. D. thesis, Technical University of Denmark, 2016

[126] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen, LAPACK Users’ Guide, 3rd Ed., Society for Industrial and Applied Mathematics, Philadelphia, PA, 1999

[127] T. Zhu, Z. Q. Zhang, Z. Gao, J. S. Wang. First-principles method to study near-field radiative heat transfer. Phys. Rev. Appl., 2020, 14(2): 024080

[128] T. Zhu, M. Antezza, J. S. Wang. Dynamical polarizability of graphene with spatial dispersion. Phys. Rev. B, 2021, 103(12): 125421

[129] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car. . Quantum ESPRESSO: A modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter, 2009, 21(39): 395502

[130] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli. . Advanced capabilities for materials modelling with Quantum ESPRESSO. J. Phys.: Condens. Matter, 2017, 29(46): 465901

[131] N. Troullier, J. L. Martins. Efficient pseudopotentials for plane-wave calculations. Phys. Rev. B, 1991, 43(3): 1993

[132] J. P. Perdew, K. Burke, M. Ernzerhof. Generalized gradient approximation made simple. Phys. Rev. Lett., 1996, 77(18): 3865

[133] H. J. Monkhorst, J. D. Pack. Special points for Brillouin-zone integrations. Phys. Rev. B, 1976, 13(12): 5188

[134] T. Zhu, P. E. Trevisanutto, T. C. Asmara, L. Xu, Y. P. Feng, A. Rusydi. Generation of multiple plasmons in strontium niobates mediated by local field effects. Phys. Rev. B, 2018, 98(23): 235115

[135] P. O. Chapuis, S. Volz, C. Henkel, K. Joulain, J. J. Greffet. Effects of spatial dispersion in near-field radiative heat transfer between two parallel metallic surfaces. Phys. Rev. B, 2008, 77(3): 035431

[136] P. Rodriguez-López, W.-K. Tse, D. A. R. Dalvit. Radiative heat transfer in 2D dirac materials. J. Phys. :Condens. Matter, 2015, 27: 214019

[137] R. Peierls. Zur Theorie des Diamagnetismus von Leitungselektronen. Eur. Phys. J. A, 1933, 80(11−12): 763

[138] M. Graf, P. Vogl. Electromagnetic fields and dielectric response in empirical tight-binding theory. Phys. Rev. B, 1995, 51(8): 4940

[139] J. Li, D. Golez, G. Mazza, A. J. Millis, A. Georges, M. Eckstein. Electromagnetic coupling in tight-binding models for strongly correlated light and matter. Phys. Rev. B, 2020, 101(20): 205140

[140] P. G. de Gennes, Superconductivity of Metals and Alloys, CRC Press, 1999

[141] R. Loudon, The Quantum Theory of Light, 3rd Ed., Oxford Univ. Press, 2000

[142] G. B. Arfken, H. J. Weber, and F. E. Harris, Mathematical Methods for Physicists, 7th Ed., Academic Press, 2013

[143] O. Keller, Quantum Theory of Near-Field Electrodynamics, Springer, Berlin, 2011

[144] J. S. Wang and J. Peng, A microscopic theory for ultra-near-field radiation, arXiv: 1607.02840 (2016)

[145] D. J. Griffiths, Introduction to Electrodynamics, 4th Ed., Cambridge Univ. Press, 2017

[146] N. N. Bogoliubov and D. V. Shirkov, Quantum Fields, Addison-Wesley, 1982

[147] G. S. Agarwal. Quantum electrodynamics in the presence of dielectrics and conductors. I. Electromagnetic-field response functions and black-body fluctuations in finite geometries. Phys. Rev. A Gen. Phys., 1975, 11(1): 230

[148] Z.-Q.ZhangJ.-T. LüJ.-S.Wang, Angular momentum radiation from current-carrying molecular junctions, Phys. Rev. B 101, 161406(R) (2020)

[149] K. Kuhnke, C. Große, P. Merino, K. Kern. Atomic-scale imaging and spectroscopy of electroluminescence at molecular interfaces. Chem. Rev., 2017, 117(7): 5174

[150] Z. Q. Zhang, J. S. Wang. Electroluminescence and thermal radiation from metallic armchair carbon nanotubes with defects. Phys. Rev. B, 2021, 104(8): 085422

[151] V. Weisskopf, E. Wigner. Berechnung der natürlichen Linienbreite auf Grund der Diracschen Lichttheorie. Eur. Phys. J. A, 1930, 63(1−2): 54

[152] W. Heisenberg, W. Pauli. Zur Quantentheorie der Wellenfelder II. Eur. Phys. J. A, 1930, 59(3−4): 168

[153] M. Creutz. Quantum electrodynamics in the temporal gauge. Ann. Phys., 1979, 117(2): 471

[154] E. Fradkin, Quantum Field Theory: An Integrated Approach, Princeton Univ. Press, 2021

[155] L. Novotny and B. Hecht, Principles of Nano-Optics, 2nd Ed., Cambridge Univ. Press, 2012

[156] S. M. Barnett. Optical angular-momentum flux. J. Opt. B, 2002, 4(2): S7

[157] S. M. Barnett, L. Allen, R. P. Cameron, C. R. Gilson, M. J. Padgett, F. C. Speirits, A. M. Yao. On the natures of the spin and orbital parts of optical angular momentum. J. Opt., 2016, 18(6): 064004

[158] Y. M. Zhang, T. Zhu, Z. Q. Zhang, J. S. Wang. Microscopic theory of photon-induced energy, momentum, and angular momentum transport in the nonequilibrium regime. Phys. Rev. B, 2022, 105(20): 205421

[159] R. M. Abraham Ekeroth, A. García-Martín, J. C. Cuevas. Thermal discrete dipole approximation for the description of thermal emission and radiative heat transfer of magneto-optical systems. Phys. Rev. B, 2017, 95(23): 235428

[160] L. Zhu, S. Fan. Persistent directional current at equilibrium in nonreciprocal many-body near field electromagnetic heat transfer. Phys. Rev. Lett., 2016, 117(13): 134303

[161] I. Latella, P. Ben-Abdallah. Giant thermal magnetoresistance in plasmonic structures. Phys. Rev. Lett., 2017, 118(17): 173902

[162] L. G. Aslamazov, A. I. Larkin. Effect of fluctuations on the properties of a superconductor above the critical temperature. Sov. Phys. Solid State., 1968, 10: 875

[163] H. A. Lorentz, Het theorema van Poynting over de energie in het electromagnetisch veld en een paar algemeene stellingen over de voortplanting van het licht, Verslagen der Afdeeling Natuurkunde van de Koninklijke Akademie van Wetenschappen 4, 176 (1895)

[164] B. Strekha, S. Molesky, P. Chao, M. Krüger, A. W. Rodriguez. Trace expressions and associated limits for nonequilibrium Casimir torque. Phys. Rev. A, 2022, 106(4): 042222

[165] R. Khrapko. Unknown spin radiation. J. Phys. Conf. Ser., 2019, 1172(1): 012055

[166] Y. M. Zhang, J. S. Wang. Far-field heat and angular momentum radiation of the Haldane model. J. Phys.: Condens. Matter, 2021, 33(5): 055301

[167] O. V. Kibis, M. R. da Costa, M. E. Portnoi. Generation of terahertz radiation by hot electrons in carbon nanotubes. Nano Lett., 2007, 7(11): 3414

[168] O. V. Dolgov, E. G. Maksimov. The dielectric function of crystalline systems. Modern Problems in Condensed Matter Sciences, 1989, 24: 221

Acknowledgements

We thank Gaomin Tang, Jingtao Lü, and Mauro Antezza for collaborations. J.-S. W. thanks Mehran Kardar for hosting a visit at MIT. He also thanks Shanhui Fan, Philippe Ben Abdallah, Matthias Krüger, and Zhuomin Zhang for discussion. This work is supported by NSFC under grant No. 12204346, MOE tier 2 grant R-144-000-411-112, and MOE Academic Research Tier 1 Fund A-8000990-00-00. Parts of the manuscript were written while visiting Kavli Institute for Theoretical Physics, University of California Santa Barbara, supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

RIGHTS & PERMISSIONS

2023 The Authors