Saturation of radiative heat transfer due to many-body thermalization

Radiative heat transfer between two bodies saturates at very short separation distances due to the nonlocal optical response of the materials. In this work, we show that the presence of radiative interactions with a third body or external bath can also induce a saturation of the heat transfer, even at separation distances for which the optical response of the materials is purely local. We demonstrate that this saturation mechanism is a direct consequence of a thermalization process resulting from many-body interactions in the system. This effect could have an important impact in the field of nanoscale thermal management of complex systems and in the interpretation of measured signals in thermal metrology at the nanoscale.

The theory of radiative heat transfer 1-6 predicts a divergence of the heat flux exchanged between two bodies kept at constant temperatures as the separation distance d between them tends to zero. During the last decade, theoretical results [7][8][9][10][11][12] have questioned this divergence and shown that it disappears when a nonlocal optical response 13 of the materials is taken into account. Recently, it has been shown that the divergence of the heat transfer can also be removed at subnanometric separation distances because of the interplay of conductive and radiative heat transfer inside the interacting bodies, which lead to the generation of temperature gradients and in turn to a saturation of the heat flux 14,15 . This effect is, however, limited to small separation distances at which new channels for heat transfer (due to phonon tunneling [16][17][18][19][20] or electron tunneling 21 ) start to play a significant role. The divergence is ultimately removed because thermal equilibrium between the bodies is established at contact 22-24 . In all these works, the interacting objects are assumed to be isolated from the environment or from other radiative sources. Here we revisit the near-field heat transfer problem between two solids when a third source of thermal radiation participates to the transfer. This situation is fundamentally different from the usual two-body description because many-body interactions are at work. Several problems in the many-body framework have recently been considered  and new thermophysical effects have been highlighted.
In this article, we investigate a heat transfer saturation mechanism due to thermalization in many-body systems under nonequilibrium conditions in the absence of nonlocal effects. While we focus on radiative heat transfer, the considered mechanism follows from the general principle of energy balance and therefore it is not unique to thermal radiation. The simplest configuration in which such a saturation mechanism can be observed is a two-body system interacting with a thermal bath. In order to describe this effect, we consider the heat transfer in the following two simple systems that may mimic many practical situations. The first one is a thin film (i.e. a membrane) that interacts with both a substrate on one side and a thermal bath on the other side, as sketched in Fig. 1(a). The second system is a small particle which also interacts with both a substrate and an external bath, see Fig. 1

Saturation mechanism for a membrane close to a substrate
Here we consider a substrate that we denote as body 1 and a membrane of thickness δ, denoted as body 2, separated by a distance d from body 1. The substrate is thermalized at a fixed temperature T 1 and the system interacts with a thermal bath of radiation at temperature < T T 3 1 , see Fig. 1(a). The thermal bath here acts as a third body. The temperature T 2 of the membrane is not fixed by a thermostat, so that this body can reach heat-transfer equilibrium at a stationary temperature T T 2 2 st = for which the net energy flux on the membrane vanishes. The radiative heat transfer originates from the electromagnetic field produced by the random thermal motion of charges inside the materials 1-6 . Expanding the electromagnetic field in plane-wave components characterized by frequency ω, parallel wave vector k, and polarization = p TM, TE, the energy flux (normal component of the Poynting vector) in the different vacuum regions of the system can be written as where γ = 1 indicates the region between bodies 1 and 2, γ = 2 labels the region on the right of body 2 [see Fig. 1 . In Fig. 2(a), we with respect to the separation distance d for several values of the thickness δ in the case in which both the slab and the membrane are made of silicon carbide (SiC). We observe a saturation of the heat flux at relatively large separation distances, where nonlocal effects are completely negligible. This saturation mechanism is directly related to the dependence of the temperature difference T T T 1 2 st Δ = − on the separation distance d. As shown in the inset of Fig. 2(b), ΔT is proportional to d 2 at short separations. Moreover, since the flux Φ approaches a constant at small d, the ratio Φ ΔT / (i.e., the heat transfer coefficient) scales as d 1/ 2 in this regime, as already outlined in the literature 53,54 for polar materials like SiC. Notice that the asymptotic value of Φ at short separations does not depend on the width δ: in the limit d 0 → , Φ correspond to the energy flux radiated to the environment by a single semi-infinite body (see below). Figure 1. Sketch of the system. (a) A membrane is placed close to a substrate at a separation distance d. The substrate is thermalized at a fixed temperature T 1 and the structure is immersed in an environmental bath of thermal radiation at temperature T 3 . The temperature T 2 of the membrane is free to reach a steady-state value = T T 2 2 st , for which the body achieves heat-transfer equilibrium. (b) A small particle is considered instead of the membrane.
www.nature.com/scientificreports www.nature.com/scientificreports/ to leading order in ΔT. Taking this into account, we can rewrite Eq. (1) as  Let us now consider how this saturation mechanism is modified for metallic materials. At short separation distances, it is well known that the heat transfer between metals is first dominated by the TE-polarization contribution and the d 1/ 2 behavior associated to TM waves is usually recovered at subnanometer separation distances. However, this divergence disappear because of nonlocal effects 10 . We show below that in metallic many-body systems a saturation of heat flux can exist at larger separation distances. To this aim, we consider a system made of a gold (Au) membrane suspended above a Au substrate.
The results for the heat flux and the steady-state temperature for a slab and a membrane made of Au are shown in Fig. 2(c,d), respectively. A saturation of the heat flux at short separations is observed also for this material. We highlight that the dependence on the thickness of the membrane is weak for δ larger than 100 nm. This is due to the fact that the electromagnetic field is completely screened for such thicknesses and the membrane becomes practically opaque. Furthermore, the behavior of ΔT in this case is shown in the inset of Fig. 2(d). Clearly, the temperature difference is not proportional to d 2 at small d because of the contribution of TE-polarized waves (TM polarization dominates well below the nanometer scale). Such a behaviour emphasizes that here the saturation mechanism is different from that for polar materials.
Asymptotic short-distance behavior. Here we analyze the asymptotic behavior of the heat flux and temperature difference at short distances. We discuss separately TM and TE polarizations, the former being the dominant contribution for the considered polar material and the latter for the metal. By neglecting the contribution of propagating waves at close separation distances, the coefficient a 1 defined in expression (3) reads 4Im ( ) Im ( ) We also assume that the thickness of the membrane is large as compared with the separation distance, which corresponds to the limit δ → ∞ in the expressions for the reflection coefficients and therefore, . For simplicity, we consider that the materials are identical and thus define ≡ = r r r Polar materials. The heat exchange between polar materials at short separations can be studied in the electrostatic limit. In this limit, only large wavevectors k k 0  contribute to the heat exchange, where k c / 0 ω = , and the normal component k z can be approximated by ik. Moreover, the Fresnel reflection coefficient for TM polarization takes the form denotes the dilogarithm function. Hence, in the limit d 0 → , the coefficients a 2 , b 1 , and b 2 in Eq. (4) remain finite since they have propagating waves contribution only, while a 1 diverges as − d 2 because of the contribution of evanescent waves in TM polarization. Thus, Δ → T 0 and b 2 Φ → as → d 0. More explicitly, in this limit the energy flux becomes which corresponds, as anticipated above, to the heat exchanged between bodies 1 and 2 together at temperature T 1 and an environment at temperature T 3 . Notice that when bodies 1 and 2 are made of the same material, so that Metals. For metals close to room temperature, the heat flux in the electrostatic limit is dominated by TM polarization at subnanometer separation distances and the traditional d 1/ 2 divergence is regularized by the presence of nonlocal effects 7,10 . However, at these separation distances, other mechanisms superimpose to the radiative transfer such as phonon [16][17][18][19][20] or electron tunneling 21 . Close to contact, these channels even dominate the heat transfer. For separation distances slightly larger (usually for  d 1 nm ), the radiative transfer in metals is entirely driven by TE-polarization states and nonlocal optical effects 7,10 do not play any role. In this case, the imaginary part of r TE decays with respect to k, so that the flux saturates 10 for a wave vector ω = k c / p max before increasing again close to contact, where ω p is the plasma frequency of the metal (see Methods). Typically, this saturation is observable between ∼ d 1 nm and separation distances similar to the skin depth of the metal evaluated at ω p (about 20 nm for Au). Nonetheless, this effect takes place at separation distances which are one order of magnitude smaller than the saturation distance induced by the thermalization process, as shown in Fig. 2(c) in our example for Au. In addition, since the transport is mediated by TE-polarized waves, the heat-transfer coefficient a 1 given by Eq. (5) remains finite at short separations. Although a 1 is finite in this regime, it is large as compared with a 2 in Eq. (4) because the latter only accounts for the contribution of propagating waves. Thus, from Eq. (4), one obtains b 2 Φ ≈ and therefore the flux is approximately given by Eq. (7), , which is small but finite in the considered limit. This behaviour is observed in the inset of Fig. 2(d) for Au.
Notice that the coefficient a 1 represents the inverse of the thermal resistance between bodies 1 and 2. Since a 1 becomes large in the limit of short separations, it implies that the thermal resistance between bodies 1 and 2 is small in this limit as compared with that between bodies 2 and 3. A similar saturation effect would be expected, with the appropriate thermal resistances, for other heat transport mechanisms different from thermal radiation.
Opaque membrane. In the example of the heat flux saturation for the metal, we have shown that the results are not sensitive to the thickness of the membrane when this is larger than 100 nm. This is due to the fact that, because of dissipation, the electromagnetic field is completely screened inside the material. In other words, the membrane becomes opaque when it is thick enough. Such a screening occurs also in polar materials, but for thicknesses typically larger than for metals: In our example of SiC, the opaque-membrane limit takes place at δ much larger than 100 μm.
Assuming that the membrane is opaque introduces a simplification in the heat-transfer problem, which then can be described through next-neighbor interactions only. This can be seen by noting that the factor e ik z2 δ in the optical reflection and transmission coefficients vanishes for large δ (see Methods), that is, when the membrane is opaque, because > k Im ( ) 0 z2 for dissipative materials. Under these conditions, we have k p r k ( , , )  . We emphasize that these energy transmission coefficients are expressed in terms of the single-interface reflection coefficients only. The energy flux Φ and stationary temperature T 2 st in the opaque membrane limit are shown in Fig. 2 as a function of the separation distance for SiC and Au. It can be seen that the value of Φ for the opaque membrane gives for all distances a lower bound on the steady-state heat flux. Furthermore, we observe that

Saturation mechanism for a particle close to a substrate
In the previous section, we have analyzed a mechanism of saturation of the heat exchange in a system with planar geometry. In this section, we extend the discussion to a situation in which a small particle is considered instead of a membrane. The particle is assumed small as compared with the thermal wave length, so that it can be modeled as a single dipole in the dipolar approximation.
The system thus consists of a substrate at temperature T 1 filling the half-space < z 0, a particle of radius R at temperature T 2 centered at the point = + x y d R r ( , , ) , and a radiative thermal bath at temperature T 3 surrounding the particle, see Fig. 1(b). The power absorbed by the particle at the point r and instant t is given by is the dipole moment of the particle and t E r ( , ) is the local electric field at the point r. Introducing the Fourier components ω f ( ) at frequency ω such that ( ) , where the first term is the free space contribution and the second term is the scattering contribution accounting for reflections on the surface of the substrate. The real part of the free space Green's function is divergent in the coincidence limit ′ → r r, but only its imaginary part contributes to the absorbed power and is given by ω = ∞ + +  Furthermore, the dipole moment of the particle can be decomposed into a fluctuating part p r ( , ) fl ω and an induced part resulting from the incident field E inc and the field produced by the dipole itself and then scattered by the surface, i.e.
where ε 0 is the vacuum permittivity and α ω ( ) is the dressed polarizability of the particle (see Methods). Hence, noting that the matrix ω  r r ( , , ) is diagonal, with simple manipulations, the components of the dipolar moment and local field can be written as (2020) 10 Taking the statistical average and using Eq. (13) leads to In what follows, for simplicity, we will omit writing down explicitly the dependence on positions of the fields, Green's functions and correlation matrices, since they will always be evaluated at the same point r in the coincidence limit. The incident field can be decomposed into a contribution coming from the substrate E ( ) 1 ω and a contribution from the bath field ω E ( ) The substrate field E ( ) 1 ω is a direct contribution to the total field at the point r propagating to the right, while the bath field ω E ( ) 3 accounts for a direct contribution propagating to the left and a reflected one propagating to the right. The correlation matrix of the incident field is given by where we have assumed that the substrate and bath fields are uncorrelated. In addition, when the substrate is in thermal equilibrium with the bath field at, for instance, temperature T 3 , the correlation matrix of the incident field can be computed from the fluctuation-dissipation theorem inc i nc 2 where the total Green's function of the system is used here because, in the absence of the dipole, the incident field is the total field on the right side of the substrate. Moreover, when evaluated at temperature T 1 , the correlation matrix of the substrate field can be written as where we have introduced the matrix ω ( )  whose explicit form is given in Methods. The correlation matrix of the bath field can thus be obtained from Eq. (15) using Eqs. (16) and (17) evaluated at temperature T 3 , which gives . In Methods we also give an explicit expression of the matrix  ω ( ). To complete the description of the problem, we need to know the correlation matrix of the fluctuating dipole moment which is given by 39 ( ) /(6 ) 3 2 3 . By using the correlation functions given above in Eq. (14), the spectral power (10) becomes which manifestly goes to zero at thermal equilibrium. In view of expression (19), the total power absorbed by the particle can be decomposed as www.nature.com/scientificreports www.nature.com/scientificreports/ is the power emitted by the particle due to its interaction with the bath in presence of the substrate. We are interested in a situation of heat-transfer equilibrium in which the total power absorbed by the particle vanishes for = T T . This situation is characterized by the stationary temperature of the particle and by the exchanged power , which here is studied as a function of the separation d. Since we describe the particle in the dipolar approximation, we restrict ourselves to separation distances larger than the radius of the particle, > d R (the distance between the substrate and the center of the particle is thus larger than R 2 ). We also emphasize that the approach developed above is appropriate for polar materials, but needs to be suitably modified for metals, introducing the magnetic contribution to the power absorbed by the dipole. With this in mind, here we consider a substrate and a particle made of SiC. In Fig. 3(a,b), we plot  and the associated Δ = − T T T 1 2 st , respectively, as a function of the separation d for = T 400 K 1 and = T 300 K 3 and for several particle radius. We observe again a saturation of the power exchanged between the substrate and the particle caused by the thermalization of the particle, whose temperature approaches that of the substrate as the separation is reduced. As shown in Fig. 3(c), there is a maximum value of this power after the transition from the far field to near-field regime, and then the exchanged power is clearly reduced as d is decreased. The corresponding particle equilibrium temperature is represented in Fig. 3(d). To highlight the influence of the thermal bath, for fixed radius (R 50 nm = ), in Fig. 3(c,d) we take the substrate temperature as , 250 K, and 300 K. In the far field, we observe a strong effect of the bath on the steady-state temperature of the particle, but the exchanged power  is similar in the different cases (recall that the total power absorbed by the particle is always zero in the considered situations).

Discussion
We have demonstrated the existence of a radiative saturation mechanism for near-field heat exchange in many-body systems. This saturation arises as a consequence of thermalization of the interacting bodies when the separation d between them is reduced. In contrast to the well-known saturation of heat transfer between two bodies, close to contact, resulting from the nonlocal response of the materials, the effect highlighted here exists even with purely local responses. For polar materials with planar geometry, a quadratic dependence of the temperature variation between neighboring elements is observed with respect to the separation distance. This dependence is counterbalanced by the d 1/ 2 scaling of the heat transfer coefficient and therefore, the energy flux reaches a constant value in the limit of small d. In metallic structures, where such a scaling does not apply, thermalization induces a saturation of the heat flux as well. In the considered example for Au, the saturation distance due to thermalization is one order of magnitude larger than the optical saturation distance 10 for the heat exchange , respectively, and the results are shown for several particle radius R. In (c,d), the substrate and bath temperatures are = T 400 K 1 and = T 0 K 3 , 250 K, and 300 K, respectively, while the particle radius is R 50 nm = .
between metals at fixed temperatures. Furthermore, we observe that the steady-state heat flux in planar geometry can be higher at intermediate separations than at short separations, and a similar result is found for the power exchanged between the dipole and the substrate. This behavior can be attributed to the non-trivial interplay between the interactions of the body reaching heat-transfer equilibrium (the membrane or the dipole) with the substrate and the bath. This mechanism of saturation due to thermalization can be observed in experimental measurements of radiative heat transfer in which the temperature of the active components is not completely fixed. This may be the case, for instance, of a membrane that is suspended by arms constituting a weak conductive channel for heat transport, or when a small object is attached to a cantilever whose internal temperature profile can be altered by the incoming radiative energy flux. The power absorbed by a particle in this more general scenario, which could represent a simplified model of a tip, can be described as abs s us env e xt = + +     , where  sus and  env account for the interaction with the substrate and the environment, respectively, and  ext is an external power that controls the state of the system. The term  env may include interactions with a bath of thermal radiation and also a conductive contribution arising from the structure supporting the particle. In the stationary state at which the absorbed power vanishes, the power ext  supplied to the system to maintain such a state can be used to infer the steady-state temperature of the particle and the radiative power  sus . As we have shown here, the induced many-body thermalization can affect both the temperature of the particle and the power exchanged with the substrate, so it can influence experimental measurements as well. Finally, the thermalization and the associated saturation effect could be relevant for thermal management in systems with several components interacting through thermal radiation. ε ω ε ω ω νω = − + with the background dielectric constant ε = 1 b , plasma frequency ω = . × 1 37 10 p 16 rad/s, and electron collision frequency ν = . × 5 32 10 13 rad/s. Furthermore, to describe the response of the particle, we assume that its nude polarizability is given by Correlation matrices of the substrate and bath fields. The matrix ω ( ) accounting for the correlations of the substrate field can be obtained by expanding the field in plane and evanescent waves and using the correlation function of the field modes 37,58 . This correlation function follows from the fluctuation-dissipation theorem (16). Then, the matrix ω ( ) describing the correlations of the thermal bath can be computed as B G S ( ) Im ( ) ( ) ω ω ω = − . A detailed derivation of these quantities is given in the Supplementary Information and here we give the final result: Im( ) Im( ) ,