Radiative heat transfer and nonequilibrium Casimir-Lifshitz forcein many-body systems with planar geometry

# Radiative heat transfer and nonequilibrium Casimir-Lifshitz force in many-body systems with planar geometry

Ivan Latella Laboratoire Charles Fabry, UMR 8501, Institut d’Optique, CNRS, Université Paris-Saclay, 2 Avenue Augustin Fresnel, 91127 Palaiseau Cedex, France    Philippe Ben-Abdallah Laboratoire Charles Fabry, UMR 8501, Institut d’Optique, CNRS, Université Paris-Saclay, 2 Avenue Augustin Fresnel, 91127 Palaiseau Cedex, France Université de Sherbrooke, Department of Mechanical Engineering, Sherbrooke, PQ J1K 2R1, Canada.    Svend-Age Biehs Institut für Physik, Carl von Ossietzky Universität, D-26111 Oldenburg, Germany    Mauro Antezza Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université de Montpellier, F- 34095 Montpellier, France Institut Universitaire de France, 1 rue Descartes, F-75231 Paris Cedex 05, France    Riccardo Messina Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université de Montpellier, F- 34095 Montpellier, France
###### Abstract

A general theory of photon-mediated energy and momentum transfer in -body planar systems out of thermal equilibrium is introduced. It is based on the combination of the scattering theory and the fluctuational-electrodynamics approach in many-body systems. By making a Landauer-like formulation of the heat transfer problem, explicit formulas for the energy transmission coefficients between two distinct slabs as well as the self-coupling coefficients are derived and expressed in terms of the reflection and transmission coefficients of the single bodies. We also show how to calculate local equilibrium temperatures in such systems. An analogous formulation is introduced to quantify momentum transfer coefficients describing Casimir-Lifshitz forces out of thermal equilibrium. Forces at thermal equilibrium are readily obtained as a particular case. As an illustration of this general theoretical framework, we show on three-body systems how the presence of a fourth slab can impact equilibrium temperatures in heat-transfer problems and equilibrium positions resulting from the forces acting on the system.

## I Introduction

Fluctuations of the electromagnetic field are responsible for momentum CasimirProcKNedAkadWet48 (); CasimirPhysRev48 (); Lifshitz55 () and heat PoldervH () exchanges. The pioneering works of Casimir CasimirProcKNedAkadWet48 () and Casimir and Polder CasimirPhysRev48 () were the first showing that such fluctuations are at the origin of an attractive force between two perfectly conducting infinite parallel planes as well as between an atom and a plate. This effect was theoretically predicted to exist even at thermal equilibrium, at zero temperature and in vacuum. A few years later, Lifshitz Lifshitz55 () and subsequently Dzyaloshinskii, Lifshitz, and Pitaevskii DzyaloshinskiiAdvPhys61 () developed a more general theory taking into account real material properties and thermal effects. Much more recently, a series of papers AntrezzaPRA04 (); AntezzaPRL05 (); AntezzaPRL06 (); AntezzaPRA08 () were focused on the effect of the absence of thermal equilibrium, showing that the presence of different temperatures in the system can not only qualitatively modify the behavior of the force (e.g. changing its power-law dependence on the distance), but also induce a repulsive force, otherwise impossible for standard geometries such as two parallel slabs. Starting from 1997, several experiments have confirmed the theoretical predictions at thermal equilibrium, both for configurations involving macroscopic bodies (mainly in the plane-plane and sphere-plane geometries) LamoreauxPRL97 (); MohideenPRL98 (); RoyPRD99 (); EderthPRA00 (); ChanPRL01 (); ChanScience01 (); BressiPRL02 (); DeccaPRL03 (); DeccaAnnPhys05 (); HarberPRA05 (); DeccaPRD07 (); KrausePRL07 (); CapassoIEEE07 (); PalasantzasAPL08 (); ChanPRL08 (); MundayPRA08 (); MundayNature09 (); JourdanEPL09 (); deManPRL09 (); ChiuPRB10 (); SushkovNatPhys11 (); ZuurbierNJP11 (); BimontePRB16 () and in the atom-plane configuration deflection (); Aspect1996 (); Shimizu2001 (); DeKieviet2003 (); Zao10 (); Pasquini1 (); Pasquini2 (). Recently, the force between a BEC and a plane has also been successfully measured out of thermal equilibrium ObrechtPRL07 (), confirming the theoretical predictions.

Another phenomenon emerging as a consequence of absence of thermal equilibrium is radiative heat transfer. Using Rytov’s theory of fluctuational electrodynamics Rytov (), Polder and van Hove PoldervH () showed that this energy-exchange mechanism, limited in the far field (for distances much larger than the thermal wavelength , close to 8 m at room temperature) by the blackbody limit predicted by Stefan-Boltzmann’s law, can overcome this value even by several orders of magnitude in the near-field regime. In particular, it was later shown that this amplification is particularly pronounced for materials supporting surface modes such as phonon polaritons JoulainSurfSciRep05 (); VolokitinRMP07 (). The radiative heat transfer has been experimentally investigated as well HargreavesPLA69 (); KittelPRL05 (); HuApplPhysLett08 (); NarayanaswamyPRB08 (); RousseauNaturePhoton09 (); ShenNanoLetters09 (); KralikRevSciInstrum11 (); OttensPRL11 (); vanZwolPRL12a (); vanZwolPRL12b (); KralikPRL12 (); KimNature15 (); StGelaisNatureNano16 (); SongNatureNano15 (); KloppstecharXiv (); WatjenAPL16 (), in a range of distances going from the nanometer region to several microns, confirming the behavior predicted theoretically.

The recent theoretical history on both topics has seen the development of a series of general theories for Casimir forces TomasPRA02 (); RaabePRA03 (); BezerraEPJC07 (); RahiPRD09 (); ReidPRL09 (); RodriguezPRA09 (); ReidPRA11 (); ReidPRA13 (), radiative heat transfer RodriguezPRL11 (); McCauleyPRB12 (); RodriguezPRB12 (); RodriguezPRB13 (); NarayanaswamyJQSRT14 (); MullerarXiv () or both in a unified approach BimontePRA09 (); KrugerPRL11 (); MessinaEurophysLett11 (); MessinaPRA11 (); KrugerPRB12 (); MessinaPRA14 (); BimontearXiv16 (). The theoretical frameworks, based on a variety of approaches (scattering matrices, Green’s functions, time-domain calculations, boundary-element method, fluctuating surface and volume currents), share the possibility of addressing bodies of arbitrary geometries and optical properties. Even if some of them can tackle the general scenario of bodies, so far only a few applications involving more than two bodies have been considered. More specifically, the heat transfer has been analyzed between three nanospheres in the dipolar approximation BenAbdallahPRL11 (); ZhuPRL16 (), between three parallel slabs ZhengNanoscale11 (); MessinaPRL12 (); MessinaPRA14 (); BenAbdallahPRL14 (); LatellaPRAppl15 (), as well as in a configuration involving one sphere between two slabs MullerarXiv (). It has to be mentioned that the radiative heat transfer in networks of more than two particles has recently received a considerable attention, but only within the dipolar approximation MessinaPRB13 (); BenAbdallahPRL13 (); LanglaisOptExpress14 (); NikbakhtJAP14 (); BenAbdallahAPL06 (); BenAbdallahPRB08 (); NikbakhtEPL15 (); BenAbdallahPRL16 (); LatellaarXiv16 (). Concerning the Casimir force, it has been discussed in the case of an atom between two slabs MessinaPRA14 (), three parallel slabs MessinaPRA14 (), and very recently the Casimir energy has been considered in the case of two and three coupled cavities when the materials undergo a phase transition from the metallic to the superconducting phase Rosa17 ().

In this paper, we focus on a system composed of planar parallel slabs made of arbitrary materials. The slabs, as well as the external environment in which they are immersed, have arbitrary temperatures. For this scenario, we derive closed-form analytical expressions for the radiative heat transfer and the Casimir force, both at and out of thermal equilibrium. To this aim, we generalize the scattering-matrix approach previously introduced for two MessinaEurophysLett11 (); MessinaPRA11 () and three MessinaPRA14 () bodies. The analytical expressions we obtain clearly highlight the nonadditive character of both momentum and energy exchange. We then consider two numerical applications, one for the Casimir force and one for heat transfer, on a system made of four slabs. For the former, we show how the equilibrium position of the central slab in a system of three slabs is modified by the introduction and lateral shift of a fourth one. For the latter, by fixing some of the temperatures in the system, we discuss the distribution of the other temperatures at local equilibrium as a function of the system parameters.

The paper is structured as follows. In Sec. II, we present our physical system and introduce the main notation and definitions. Then, in Secs. III and IV, we formally derive the expressions of the Poynting vector and the stress tensor, respectively. In Sec. V, we give explicit expressions of the scattering coefficients, while the energy transmission coefficients and the momentum transfer coefficients are computed in Secs. VI and VII, respectively. Then, we present in Sec. VIII our numerical applications to both Casimir force and radiative heat transfer. We finally give some conclusive remarks in Sec. IX.

## Ii Many-body systems

The system we address is composed of planar slabs orthogonal to the axis and assumed to be infinite in the and directions. Each slab is at equilibrium at temperature , thermalized by means of some external source, with . Let the temperature of the left environment be and that of the right environment be , which can be seen as the equilibrium temperatures of two blackbodies at infinity to which we assign the labels and , respectively. Hereafter we refer to this configuration as an -body system. This distribution of bodies defines vacuum regions that we denote by . Moreover, for , the -th slab has thickness and is centered at on the axis, as shown in Fig. 1.

We want to compute energy and momentum fluxes across a surface within any of the vacuum regions of the system. The electric field at a point and time in region , which is created by the fluctuating currents inside the materials, can be expressed as a Fourier expansion given by

 Eγ(R,t)=∫∞−∞dω2πe−iωtEγ(R,ω), (1)

where is the frequency. We require that in order for the field to be real, where the asterisk denotes complex conjugation. In addition, the components can be decomposed using a plane-wave description MessinaPRA14 (), in such a way that a single mode of the field is specified by the frequency , the component of the wave vector on the - plane, the two polarizations , and the direction of propagation along the axis which is denoted by . The latter can take two values: indicating propagation to the right and that indicates propagation to the left. The total wave vector reads , where the component is given by

 kz=√ω2c2−k2, (2)

being the speed of light in vacuum. We note that when , the component is real and, therefore, the wave is propagative. When , is imaginary and the associated wave is evanescent. Evanescent waves are nonpropagating modes of the field: for these modes indicates the direction along which the amplitude of the wave decays. Thus, taking this plane-wave decomposition into account, the single-frequency component of the electric field can be written as

 Eγ(R,ω)=∫d2k(2π)2∑p,ϕeiKϕ⋅R^ϵϕ(ω,k,p)Eγϕ(ω,k,p), (3)

where are the components of the electric field in this decomposition and are the polarization vectors. Besides, using the Maxwell equation , the magnetic field can be obtained from the electric field and hence, the Fourier components can be expanded in terms of the plane-wave components .

Each mode of the total field in any region depends on the fields generated by all the bodies as well as on the background fields present in the left and right environments. More specifically, the mode in region results from the source field modes that we denote (see Fig. 1) with , . At this stage we observe that the fact of considering only parallel planar slabs introduces a major simplification in the calculation. As a matter of fact, the stationarity of the problem, as well as the translational invariance along the and axes safely allow us to state that the frequency , the wave vector and the polarization are conserved in any scattering process. Moreover, thanks to the cylindrical symmetry with respect to the axis, any reflection and transmission coefficient will depend only on the modulus of the wave vector. As a result, the total field in each region can be written as a linear combination of the form

 Eγϕ=∑j,ηLγϕjηEjη, (4)

in terms of the coefficients . Hereafter, if not otherwise explicitly stated, summations over indices labeling bodies run from to , including quantities associated with the environmental fields. Equation (4) will allow us in the following to relate the statistical properties of the total field in each region, needed to calculate both the Casimir force and the heat transfer, to the statistical properties of the individual source fields, directly derived from the fluctuation-dissipation theorem.

We can now move toward an explicit expression of the radiative heat transfer on each slab. To this aim, we should first observe that our assumption of infinite extension of all the slabs actually leads to a formally infinite flux. Of course, the same issue applies to Casimir force as well. Nevertheless, the translational invariance allows us to think in terms of flux and force per unit surface. As discussed in detail in Ref. MessinaPRA11 (), the transition from the calculation of the total flux (or force) to the one per unit surface basically consists of omitting an infinite Dirac delta and its coefficient in the final expressions.

Starting from the case of heat transfer, the energy flux per unit surface in region is given by the averaged component of the Poynting vector , calculated at , where is located in zone . In Cartesian components, we have

 ⟨Sγi(R,t)⟩=ε0c2∑j,kϵijk⟨Eγj(R,t)Bγk(R,t)⟩, (5)

where indicates symmetrized statistical average, is the Levi-Civita symbol with , and is the vacuum permittivity. In order to compute this quantity, we introduce the correlation functions which are defined according to

 (6)

where the dagger denotes hermitian conjugate. Thus, using Eqs. (5) and (6), after manipulating the polarization vectors, the averaged component of the Poynting vector takes the form MessinaPRA14 ()

 Φγ≡⟨Sγz(R,t)⟩=∫∞0dω2π∫∞0dk2πk∑p,ϕ,ϕ′2ε0c2kzωϕ×[Πpwδϕϕ′+Πew(1−δϕϕ′)]Cγϕϕ′, (7)

where we have introduced the projectors on the propagating and evanescent wave sectors and , respectively, defined by

 Πpw≡θ(ω−ck),Πew≡θ(ck−ω), (8)

being the Heaviside step function. We observe that, as a consequence of Eq. (6), the energy flux is stationary and invariant under translations in the - plane. Notice that the dependence on in Eq. (7) is implicit through the correlation functions of the total field in a particular region of the system. Furthermore, we already have an expression [Eq. (4)] connecting the total field in each region to the fields emitted by the bodies and the environments. We now perform the so-called local-thermal-equilibrium approximation, i.e. we assume that each body radiates as it would do at equilibrium at its own temperature, so that the modes of the source fields corresponding to different bodies are not correlated to each other. Denoting by the correlation functions of the source fields associated to body , with , the assumption of local thermal equilibrium leads to MessinaPRA14 ()

 Cγϕϕ′=∑j,η,η′LγϕjηLγϕ′∗jη′Cjηη′. (9)

In addition, for convenience, we introduce the coefficients defined according to

 [Πpwδϕϕ′+Πew(1−δϕϕ′)]Cjηη′≡ℏω2Nj2ε0c2kzKjηη′ϕϕ′, (10)

where

 Nj(ω)=nj(ω)+12,nj(ω)=(eℏω/kBTj−1)−1, (11)

is the Boltzmann constant, and the reduced Planck constant. Taking into account the previous definitions, the energy flux (7) becomes

 Φγ=∫∞0dω2π∫∞0dk2πk∑p,jℏωNjXγ,j, (12)

where the coefficients are given by

 Xγ,j=∑ϕ,ϕ′,η,η′ϕLγϕjηLγϕ′∗jη′Kjηη′ϕϕ′. (13)

We note that the dependence on the temperature in Eq. (12) is explicit through the functions and possibly implicit in the optical properties (and thus in the scattering amplitudes) of the bodies. We show in Appendix A that the above coefficients always satisfy

 ∑jXγ,j=0. (14)

As a consequence, if all the bodies are thermalized at the same equilibrium temperature, using Eq. (14) in Eq. (12) one sees that the flux vanishes in each region .

In view of Eq. (14), we observe here that and, therefore, purely quantum contributions associated to zero-point fluctuations do not participate in the energy fluxes. Furthermore, since the net energy flux on body is given by , taking into account Eq. (14) this flux can be written as

 Φj=∫∞0dω2π∫∞0dk2πk∑p∑ℓ≠jℏωnℓ,jTℓ,j, (15)

where we have introduced and the energy transmission coefficients given by

 Tℓ,j=Xj−1,ℓ−Xj,ℓ,j=1,…,N. (16)

In this way, as can be seen from Eq. (15), energy fluxes are described with a Landauer-like formalism in many-body systems. The above energy transmission coefficients will be computed in Sec. VI, where, in particular, it can be seen that they satisfy the reciprocity relation , with .

We highlight that as a consequence of Eq. (14) and the definition (16), the energy transmission coefficients satisfy the following remarkable property: . On the one hand, using this result one immediately obtains . Since with are positive quantities, the coefficients are negative. On the other hand, consider now a situation in which all the bodies are assumed to be thermalized at , except body for which . Under these conditions and according to the previous reasoning, from Eq. (15), the net energy flux on body can be expressed as

 Φj=∫∞0dω2π∫∞0dk2πk∑pℏωnjTj,j. (17)

Notice that . The above expression allows us to interpret the coefficient as the self-emission amplitude associated to body , since it accounts for the radiation emitted by this body in the presence of the rest of the system. In this case, corresponds to the self-emission rate discussed in Ref. MullerarXiv ().

## Iv Casimir-Lifshitz forces

We now focus on the Casimir-Lifshitz forces acting on the system. The formulation we introduce here is analogous to the previous one for heat transfer, but now momentum fluxes are the relevant quantities to be considered instead of energy fluxes. Since energy and momentum fluxes are quantities with different physical nature, however, a priori one expects some differences to arise when comparing the two descriptions. For instance, the momentum transfer coefficients we introduce below are formally analogous to the energy transmission coefficients, but the former are represented by complex numbers while the latter are real. This fact is easy to understand if one bears in mind that the energy quanta are always real, whereas the momentum component becomes imaginary for photons characterizing evanescent fields. Another important difference is that nonvanishing momentum fluxes occur even at thermal equilibrium and, as is well known, purely quantum fluctuations contribute to these fluxes as well. As a consequence, reciprocity will not hold for these momentum transfer coefficients.

To describe Casimir-Lifshitz forces, let us consider the Maxwell stress tensor in a particular region of the system, whose Cartesian components read

 Tγij(R,t)=ε0[Eγi(R,t)Eγj(R,t)+c2Bγi(R,t)Bγj(R,t)]−ε02δij[|Eγ(R,t)|2+c2|Bγ(R,t)|2] (18)

with . The momentum flux in region is given by the averaged component and takes the form MessinaPRA14 ()

 Pγ≡⟨Tγzz(R,t)⟩=−∫∞0dω2π∫∞0dk2πk∑p,ϕ,ϕ′2ε0c2k2zω2×[Πpwδϕϕ′+Πew(1−δϕϕ′)]Cγϕϕ′. (19)

As for the case of heat transfer, the above expression is obtained by expanding the electric field in the plane-wave representation (3), using Maxwell’s equations to obtain the magnetic field in terms of the electric plane-wave components, manipulating the polarization vectors, and introducing the correlation functions with the help of Eq. (6). Moreover, also here we will take into account the local-equilibrium approximation to write the correlation functions of the total field in terms of the correlation functions of the source fields. According to this, using Eqs. (9) and (10), Eq. (19) can be rewritten as

 Pγ=−∫∞0dω2π∫∞0dk2πk∑p,jℏkzNjYγ,j, (20)

where the coefficients are given by

 Yγ,j=∑ϕ,ϕ′,η,η′LγϕjηLγϕ′∗jη′Kjηη′ϕϕ′. (21)

The net force per unit area acting on body can be computed as , so that using Eq. (20) we can write

 Pj=∫∞0dω2π∫∞0dk2πk∑p,ℓℏkzNℓWℓ,j, (22)

where we have introduced the previously mentioned momentum transfer coefficients defined by

 Wℓ,j=Yj−1,ℓ−Yj,ℓ,j=1,…,N. (23)

In addition, the previous expression of the net pressure , Eq. (22), can be conveniently rewritten as

 Pj=∫∞0dω2π∫∞0dk2πk∑pℏkzNj∑ℓWℓ,j+∫∞0dω2π∫∞0dk2πk∑p∑ℓ≠jℏkznℓ,jWℓ,j. (24)

Hence, in particular, we see that at thermal equilibrium at temperature , the functions in the second term of Eq. (24) vanish and, therefore, the net pressure on body reduces to

 Pjeq=∫∞0dω2π∫∞0dk2πk∑pℏkzNj∑ℓWℓ,j. (25)

The momentum transfer coefficients will be given in Sec. VII. To obtain these coefficients (and the energy transmission coefficients), in the next section we first introduce the tools needed to solve the many-body scattering problem in terms of single-body reflection and transmission coefficients.

## V Scattering and electric field coefficients

Having introduced a formal method to compute the couplings driving energy and momentum exchanges, a procedure to relate such couplings to individual properties is required. To this aim, below we start by considering a systematic way to characterize many-body scattering processes in terms of optical properties of the single constituents of the system. Using this procedure, we subsequently determine the electric field coefficients which are necessary to express energy transmission and momentum transfer coefficients.

### v.1 Many-body scattering coefficients

The scattering operators introduced in Refs. MessinaPRA11 (); MessinaPRA14 () for two- and three-body systems, which for our geometry reduce simply to coefficients, are a useful tool that permits us to write physical quantities in a convenient way. Our aim here is to introduce such coefficients for the -body case we are concerned with.

The many-body scattering coefficients take into account the presence of different bodies at the same time and are built in terms of the single-body reflection and transmission coefficients denoted by and , respectively, where labels the associated body. For the reflection coefficients, specify the direction of propagation or decay of the outgoing field (the incoming field propagates or decays in the direction ), whereas the transmission coefficients do not depend on . On the one hand, for , these coefficients can be written as

 ρjϕ=ρje−ikz(δj+ϕ2zj),τj=τje−ikzδj, (26)

with and depending implicitly on the optical properties of the single body. For the geometry under consideration and for isotropic media, the latter take the form

 ρj=rp,j1−e2ikzjδj1−r2p,je2ikzjδj,τj=(1−r2p,j)eikzjδj1−r2p,je2ikzjδj. (27)

Here the component of the wave vector inside medium reads

 kzj=√ω2c2εjμj−k2, (28)

and the vacuum-medium Fresnel reflection coefficients (for ) are given by

 rTE,j=μjkz−kzjμjkz+kzj,rTM,j=εjkz−kzjεjkz+kzj, (29)

where and are the dielectric permittivity and magnetic permeability of body , respectively. On the other hand, for the blackbodies radiating the environmental fields ( and ) we set

 ρ0ϕ=ρN+1ϕ=0,τ0=τN+1=0. (30)

In addition, the correlation functions of the source fields, including those of the environments, are defined in terms of the previously introduced single-body reflection and transmission coefficients. These correlation functions are given by MessinaPRA11 ()

 (31)

Using this result, the coefficients introduced in Eq. (10) can be easily computed and expressed as

 (32)

In order to define the scattering coefficients for the -body case, consider a block of consecutive bodies having indexes going from to (with ), and let us denote the sequence of these bodies by . The reflection and transmission coefficients for this block, and , representing the analogues of and for a single body, are given by

 ρj→m+=ρℓ→m++(τℓ→m)2uj→ℓ−1,ℓ→mρj→ℓ−1+,ρj→m−=ρj→ℓ−1−+(τj→ℓ−1)2uj→ℓ−1,ℓ→mρℓ→m−,τj→m=τj→ℓ−1uj→ℓ−1,ℓ→mτℓ→m, (33)

where and

 uj→ℓ−1,ℓ→m=∞∑n=0(ρj→ℓ−1+ρℓ→m−)n=(1−ρj→ℓ−1+ρℓ→m−)−1. (34)

The coefficient as expressed in Eq. (33), for example, accounts for the reflection of a mode to the right due to bodies together; it has a direct contribution from the reflection produced by bodies , and a contribution that takes into account that the mode is transmitted through bodies , undergoes multiple reflections within the cavity formed by bodies and , is reflected to the right by bodies , and finally leaves the cavity by transmission through bodies . Analogously, the coefficient given in Eq. (33) represents transmission through bodies , multiple reflections between bodies and , and transmission through bodies .

According to the above expressions, the many-body scattering coefficients can be equivalently computed in several ways by choosing different allowed values of . For convenience, however, below we introduce a setup that is particularly useful for the problem at hand. We rewrite these coefficients as

 ρj→m+=^ρj→m+e−ikz(δm+2zm),ρj→m−=^ρj→m−e−ikz(δj−2zj),τj→m=^τj→mexp(−ikzm∑ℓ=jδℓ), (35)

where from Eq. (33) we have, for ,

 ^ρj→m+=ρm+(τm)2^ρj→m−1+uj→m−1,me2ikzdm−1,^ρj→m−=ρj+(τj)2^ρj+1→m−uj,j+1→me2ikzdj,^τj→m=^τj→m−1uj→m−1,mτm, (36)

and

 uj→m−1,m=(1−^ρj→m−1+ρme2ikzdm−1)−1,uj,j+1→m=(1−ρj^ρj+1→m−e2ikzdj)−1. (37)

Here is the separation distance between the consecutive bodies and , given by

 dj=zj+1−zj−δj/2−δj+1/2, (38)

and which corresponds to the width of region . Notice that the above coefficients are to be taken as and for a single body.

### v.2 Electric field coefficients

To obtain the expressions for the energy transmission and momentum transfer coefficients, we have to determine first the coefficients relating source fields to total fields in a given region of the system. As stated by Eq. (4), the coefficients account for the contribution of the field mode , emitted by the source in direction , to the total field mode in region and direction . On this basis and using the many-body scattering coefficients, we are able to directly write down some useful relations between these coefficients that will allow us to find . For instance, recalling that and , we can write

 Lγ−jη=ργ+1→N+1−Lγ+jη,j≤γ,Lγ+jη=ρ0→γ+Lγ−jη,j>γ. (39)

The first of these relations indicates that the contribution of the source mode to the total mode is proportional to the contribution of the same source mode to the total mode . If the source is located on the left of the considered region (), the proportionality factor is the backward reflection coefficient of the block formed by all the bodies on the right of the region (see Fig. 1 for illustration). An analogous reasoning applies to the second of Eqs. (39).

We now want to relate to for . On the one hand, since the (left) environment only radiates to the right, the coefficient must vanish for . In addition, this coefficient also vanishes for , because the field emitted by this source to the left is not reflected back into the system but is absorbed by the environment. On the other hand, for sources such that , one realizes that a mode emitted to the right by the source is equivalent to a mode emitted by this source to the left undergoing the following scattering process: multiple reflections in the cavity formed by bodies and body , reflection to the right by bodies , and transmission through body . Taking this into account, we thus can write

 Lγ+j−={0,j=0ρ0→j−1+u0→j−1,jτjLγ+j+,0

Notice that vanishes when , so in this case, as previously explained, vanishes as well. Following similar arguments one can also deduce that

 Lγ−j+={τjuj,j+1→N+1ρj+1→N+1−Lγ−j−,γ

We highlight at this point that in virtue of properties (39), (40), and (41), the coefficients are completely determined if and are explicitly known for and for , respectively. Again, these coefficients can be easily obtained using the many-body scattering coefficients; let us focus on for . We first note that when , the contribution of mode to the total field is simply given by the factor accounting for multiple reflections in the cavity formed by bodies and bodies , i.e. . When , in addition to the previous factor, the contribution of mode is affected by a factor that describes the scattering from region to the considered region . With this, we obtain

 Lγ+j+=u0→γ,γ+1→N+1{u0→j,j+1→γτj+1→γ,j<γ1,j=γ. (42)

Finally, by symmetry or by employing analogous arguments, one arrives at the conclusion that

 Lγ−j−=u0→γ,γ+1→N+1⎧⎨⎩1,j=γ+1uγ+1→j−1,j→N+1×τγ+1→j−1,j>γ+1. (43)

We have determined the coefficients , so we are now able to perform the calculation of energy and momentum fluxes.

## Vi Energy transmission coefficients

In this section we determine the energy transmission coefficients . To proceed, we introduce the set of coefficients defined by

 ^Tjγ≡j∑ℓ=0Xγ,ℓ, (44)

As discussed below, the energy transmission coefficients can be fully determined in terms of . In Appendix A we show that the coefficients take the form ()

 (45)

whereas, in accordance with Eq. (14),

 ^TN+1γ=∑ℓXγ,ℓ=0. (46)

We observe that these coefficients satisfy the symmetry property

 ^Tjγ=^Tγj,j,γ=0,…,N. (47)

Furthermore, the radiative heat transfer problem can be equivalently formulated in terms of the coefficients . According to the definition (44), the energy flux (12) can be rewritten as

 Φγ=∫∞0dω2π∫∞0dk2πk∑pN∑j=0ℏωnj,j+1^Tjγ. (48)

Moreover, taking into account that the net energy flux on body is given by , this flux becomes

 Φj=∫∞0dω2π∫∞0dk2πk∑pN∑ℓ=0ℏωnℓ,ℓ+1(^Tℓj−1−^Tℓj). (49)

Thus, in view of Eqs. (48) and (49), we see that this formulation is particularly useful when a sequence of consecutive bodies in the system are thermalized at the same temperature; in this case the functions cancel out and the corresponding terms do not contribute to the fluxes.

To establish the relation between and , we first can write

 Xγ,0=^T0γ,Xγ,j=^Tjγ−^Tj−1γ,j=1,…,N,Xγ,N+1=−^TNγ, (50)

where the first two relations follow directly from the definition (44), and the last one is obtained using the fact that and Eq. (46). Thus, replacing Eq. (50) in the definition of the energy transmission coefficients given by Eq. (16) leads to the desired relation:

 T0,j=^T0j−1−^T0j,Tℓ,j=^Tℓj−1−^Tℓ−1j−1−^Tℓj+^Tℓ−1j,TN+1,j=−^TNj−1+^TNj, (51)

where in the these expressions .

In addition, we highlight that vanish if body is removed from the system, which can be achieved by letting and . Accordingly and in view of Eq. (51), the associated transmission coefficient also vanishes under these conditions, as expected, since this coefficient represents the energy exchange channel between bodies and . Moreover, taking into account the property (47), from Eq. (51) one deduces that the energy transmission coefficients satisfy the reciprocity relation for .

Finally, we note that the contribution of the environmental fields to the coefficients has to be evaluated by means of Eqs. (30) separately for the cases . In the remaining coefficients , those for which , the contribution of the environmental fields can be straightforwardly evaluated since, for instance, and . Once this is done, the many-body scattering coefficients have to be expressed using Eqs. (35) to remove the dependence of the reflection coefficients on the positions and make explicit the dependence on the separation distances (the whole system is invariant under translations along the axis).

## Vii Momentum transfer coefficients

In this section we determine the momentum transfer coefficients introduced in Eq. (23). The procedure we adopt here is close to that we followed in Sec. VI to obtain the transmission coefficients .

In analogy with Eq. (44), we now introduce the set of coefficients defined as

 ^Wjγ≡j∑ℓ=0Yγ,ℓ−12∑ℓYγ,ℓ, (52)

 ^WN+1γ=12∑ℓYγ,ℓ. (53)

Moreover, in Appendix A we show that the set of coefficients take the form ()

 (54)

On the other hand, from Eq. (52) we can write

 Yγ,0=^W0γ+^WN+1γ,Yγ,j=^Wjγ−^Wj−1γ,j=1,…,N,Yγ,N+1=−^WNγ+^WN+1γ. (55)

Thus, using these relations, from Eq. (23) we get

 W0,j=^W0j−1+^WN+1j−1−^W0j−^WN+1j,Wℓ,j=^Wℓj−1−^Wℓ−1j−1−^Wℓj+^Wℓ−1j,WN+1,j=−^WNj−1+^WN+1j−1+^WNj−^WN+1j, (56)

where in these expressions . This fully determines the coefficients in terms of .

As for the case of energy, the momentum fluxes can be equivalently formulated in terms of the coefficients . Using the definition (52), the momentum flux (20) can be rewritten as

 (57)

Since the net force per unit area acting on body is given by , we get

 Pj=∫∞0dω2π∫∞0dk2πk∑pN∑ℓ=0ℏkznℓ,ℓ+1(^Wℓj−1−^Wℓj)+∫∞0dω2π∫∞0dk2πk∑pℏkz(N0+NN+1)×(^WN+1j−1−^WN+1j). (58)

Again, in view of Eq. (58), we see that this formulation is particularly useful when a sequence of consecutive bodies in the system are thermalized at the same temperature.

We note that taking into account that

 ∑ℓWℓ,j=2(^WN+1j−1−^WN+1j), (59)

at thermal equilibrium Eq. (58) becomes as given by Eq. (25). Using Eq. (54) we immediately deduce that

 N∑j=1(^WN+1j−1−^WN+1j)=^WN+10−^WN+1N=0