Near-Field Radiative Heat Transfer Eigenmodes
Stephen Sanders, Lauren Zundel, Wilton J. M. Kort-Kamp, Diego A. R. Dalvit, Alejandro Manjavacas
NNear-Field Radiative Heat Transfer Eigenmodes
Stephen Sanders, Lauren Zundel, Wilton J. M. Kort-Kamp, Diego A. R. Dalvit, and Alejandro Manjavacas
3, 1, ∗ Department of Physics and Astronomy, University of New Mexico, Albuquerque, New Mexico 87106, United States Theoretical Division, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, United States Instituto de ´Optica (IO-CSIC), Consejo Superior de Investigaciones Cient´ıficas, 28006 Madrid, Spain (Dated: February 12, 2021)The near-field electromagnetic interaction between nanoscale objects produces enhanced radiativeheat transfer that can greatly surpass the limits established by far-field black-body radiation. Here,we present a theoretical framework to describe the temporal dynamics of the radiative heat transferin ensembles of nanostructures, which is based on the use of an eigenmode expansion of the equationsthat govern this process. Using this formalism, we identify the fundamental principles that determinethe thermalization of collections of nanostructures, revealing general but often unintuitive dynamics.Our results provide an elegant and precise approach to efficiently analyze the temporal dynamics ofthe near-field radiative heat transfer in systems containing a large number of nanoparticles.
The thermal radiation exchanged between macroscopicbodies separated by macroscopic distances is accuratelydescribed by Planck’s law [1]. However, this descrip-tion breaks down when the distance between the objectsor their size becomes significantly smaller than the so-called thermal wavelength, which, for a temperature T ,is λ T = 2 π ¯ hc/ ( k B T ). In this limit, the contribution ofnear-field components of the electromagnetic field [2–10],together with the strong responses provided by the elec-tromagnetic resonances of nanostructures [11–18], resultsin an enhanced radiative heat transfer (RHT), which cansurpass the black-body limit by several orders of magni-tude [19–23].Near-field RHT is usually described within the frame-work of fluctuational electrodynamics [23, 24]. In par-ticular, when considering collections of finite nanostruc-tures, a dipole approximation, in which each nanoparti-cle is modeled as a fluctuating dipole, can be exploited[16, 23, 25–27]. By doing so, it is possible to calculate thepower transferred between the different constituents fora particular fixed distribution of temperatures [16, 28–30]. However, if one is interested in understanding thetemporal evolution of the temperatures of the particlesin the ensemble, this approach presents several disadvan-tages. Specifically, since the power transferred betweenthe particles in the ensemble depends on their tempera-tures, which change over time, it is necessary to performa new calculation at each step in the temporal evolu-tion [31–34]. As a result, this approach does not providedirect insight into the fundamental principles that deter-mine the thermalization dynamics and, in addition, canbe very computationally demanding when the number ofparticles is sufficiently large.Here, we present a different approach to describe thethermalization dynamics of an ensemble of nanoparti-cles. Our approach is based on the linearization of theequations that govern the power transferred between thenanoparticles, which allows us to convert them into aneigenvalue problem. By doing so, we are able to finda set of RHT eigenmodes for the ensemble, which com- pletely describe the temporal evolution of the system un-der any possible initial temperature distribution. Usingthis approach, we are able to identify the general prin-ciples that control the thermalization process mediatedby near-field RHT, which often give rise to unintuitivebehaviors. This insight leads us to explore exotic scenar-ios, including dynamics in which the temperature of aparticle oscillates around the equilibrium temperature asit thermalizes. The simplicity of this formalism makes itan elegant and efficient method to describe the dynamicsof the near-field RHT in ensembles ranging from only afew nanoparticles up to thousands.We consider an ensemble of N nanospheres with radii R i and temperatures T i , placed at positions r i and sur-rounded by an environment at T , which we fix to 300 Kfor the remainder of this letter. We assume that, forall particles, R i (cid:28) λ T and all interparticle distances d ij = | r i − r j | ≥ R i , R j ), but significantly smallerthan λ T . Therefore, we can model the nanoparticles asfluctuating dipoles with electric polarizabilities α i . Fol-lowing previous works [25–27, 31], the power absorbedby particle i is given by (see the Appendix for details) P i = N (cid:88) j =1 (cid:90) ∞ dωf ij ( ω ) [ n ( ω, T j ) − n ( ω, T )] , (1)where n ( ω, T ) = [exp(¯ hω/k B T ) − − is the Bose-Einstein distribution at temperature T and f ij ( ω ) =(2¯ hω/π )Tr (cid:2) Im (cid:8) A ij Im { χ j } C + ij (cid:9)(cid:3) . In this expression,“+” represents the conjugate transpose, the trace is takenover Cartesian components, and the different matrices,which have dimensions 3 N × N , are defined as: A =[ I − α G ] − , C = (cid:0) G + G (cid:1) A , and χ = α − G α + α ,with I being the identity matrix, α a matrix with theparticle polarizabilities, G the dipole-dipole interactiontensor, and G = iω c I . This model can be readily gen-eralized to particles with magnetic response by includinga magnetic polarizability [16, 30].The temporal evolution of the temperatures of thenanoparticles is determined by the ratio between the a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b Eigenmode λ μ ( s - ) Chain (c) λ = 2.3 s -1 λ = 2884.5 s -1 λ = 7833.6 s -1 λ = 11163.1 s -1 λ = 2.3 s -1 λ = λ = 8470.8 s -1 λ = 12596.6 s -1 (b) R = 25 nm, d = 4 R R d d Rd λ μ (s -1 ) R = 25 nm d = 4 Rd = 12 R max-max Δ (a) Square
FIG. 1. (a) Schematics of the systems under study, consistingof four SiC nanospheres arranged in either a chain (left) ora square (right). (b) Decay rates of the RHT eigenmodes ofthe two systems assuming R = 25 nm and d = 4 R . The insetsdisplay the components of each of the RHT eigenmodes andthe explicit value of the associated decay rate. (c) Decay ratesfor both systems for different values of d . power they absorb P i and their heat capacities γ i . Byexpanding n ( ω, T j ) around T , as n ( ω, T j ) ≈ n ( ω, T ) +∆ T j ∂n ( ω, T ) /∂T | T = T , with ∆ T j = T j − T , we can lin-earize Eq. (1) to obtain the differential equation that gov-erns the temporal evolution of the nanoparticle temper-atures, ddt ∆ T ( t ) = − H ∆ T ( t ) . (2)Here, the matrix H = Γ − F is the product of the inverseof a diagonal matrix Γ containing the heat capacities ofthe nanoparticles γ i and a symmetric matrix F , whosecomponents are defined as F ij = − (cid:90) ∞ dωf ij ( ω ) ∂n ( ω, T ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) T = T . As shown in the Appendix, the particular structure of H ensures its diagonalizability. This allows us to write thegeneral solution of Eq. (2) using its eigenvalues λ µ andeigenvectors ∆ T ( µ ) as∆ T ( t ) = N (cid:88) µ =1 a µ e − λ µ t ∆ T ( µ ) , (3) where the coefficients a µ are obtained from the weightedinner product between ∆ T ( µ ) and the vector containingthe initial temperatures of the nanoparticles ∆ T (0) as a µ = (cid:80) Ni =1 γ i ∆ T i (0)∆ T ( µ ) i , with the eigenvectors satisfy-ing (cid:80) Ni =1 γ i ∆ T ( µ ) i ∆ T ( ν ) i = δ µν . Therefore, we concludefrom Eq. (3) that the dynamics of the near-field RHT ofan ensemble of nanoparticles can be completely under-stood by analyzing its RHT eigenmodes and the corre-sponding decay rates given, respectively, by the eigenvec-tors and eigenvalues of H . Importantly, H is a positivedefinite matrix and therefore all of the decay rates arepositive, which ensures that the ensemble thermalizes as t → ∞ .This approach assumes that the temperature depen-dence of the material properties of the nanoparticles canbe neglected. Furthermore, its accuracy improves as theratios max( | ∆ T j | /T ) and ¯ hω / ( k B T ) decrease. Here, ω represents the characteristic frequency of the electro-magnetic response of the nanoparticles. For the systemsunder consideration, the results of the eigenmode ap-proach have very good agreement with the non-linearizedfull calculation up to max( | ∆ T j | /T ) ≈ /
3, as shown inFig. S1.To illustrate the developed framework, we first con-sider a simple example, although the conclusions we draware general to any ensemble of nanoparticles. In partic-ular, we analyze the two systems depicted in Fig. 1(a),consisting of N = 4 identical SiC spherical nanoparticlesarranged in either a chain or a square (see Fig. S2 fora similar analysis of a system with N = 2197). We ob-tain the polarizability of the particles from the dipolarMie coefficient [35] using the dielectric function ε ( ω ) = ε ∞ (cid:2) ω − ω ) / ( ω − ω − iωτ − ) (cid:3) , with ε ∞ = 6 . hω T = 98 . hω L = 120 meV, and ¯ hτ − = 0 .
59 meV[36]. Figure 1(b) analyzes the RHT eigenmodes of thechain (black) and the square (gray) assuming that theparticles have a radius R = 25 nm and are separated by d = 4 R . The chain has four distinct eigenmodes, whilethe larger symmetry of the square results in two of itsmodes being degenerate. Since particles with the sametemperature do not exchange heat with one another, ev-ery ensemble, including the two analyzed here, must al-ways have an eigenmode with equal amplitude in all par-ticles. This eigenmode, which we label as µ = 1, rep-resents a net transfer of heat between the ensemble andthe environment and, as we explain below, always has theslowest decay rate. The orthogonality of the eigenmodesforces the rest of them to satisfy (cid:80) Ni =1 γ i ∆ T ( µ> i = 0,which physically means that they represent processes inwhich the heat stored in the ensemble remains constant.Therefore, every eigenmode with µ > µ increases, the length scale over whichthe sign of the components alternates, and hence the -3 -2 -1 Δ T ( K ) N = 49R = 25 nm (a) -2 -1 -3 Δ T ( K ) -4 -2 -1 -3 -5 t (s) N = 490 R = 25 nm d min = 4 R (b) d = 4 Rd = 12 R FIG. 2. (a) Thermalization dynamics for a square array of N = 49 SiC nanoparticles with R = 25 nm under differentinitial conditions. The red, blue, and green curves display theevolution of the temperature of the particle of that color in theinset schematics, when such particle is initially at ∆ T = 49 Kand all of the others at ∆ T = 0 K. The gray curves representthe case in which all of the particles of the ensemble are ini-tially at ∆ T = 49 /N K. In all cases, solid and dashed curvescorrespond to d = 4 R and d = 12 R , respectively. (b) Sameas (a), but for an ensemble of N = 490 SiC nanoparticleswith R = 25 nm, randomly distributed inside a spherical vol-ume of radius 600 nm with a minimum interparticle distance d min = 4 R (see schematics). near-field RHT occurs, decreases. This is consistent withthe increase of the associated decay rate, whose valueis dominated by terms proportional to ( R i R j ) /d ij . Onthe contrary, λ describes the net radiation exchange be-tween the ensemble and the environment, which scales as( R i /λ T ) . Therefore, for near-field RHT ( i.e. , d ij (cid:28) λ T ), λ always has the smallest value among all of the decayrates, although, as shown in Fig. 1(c), the difference be-tween λ and the rest of the decay rates is reduced byincreasing the distance between the particles.We know from Eq. (3) that the thermalization of anensemble of particles is initially dominated by the eigen-modes with largest decay rates. However, for sufficientlylong time, this process is controlled by the first eigen-mode, which, as discussed above, has equal amplitude inall particles and, consequently, its decay rate is the small-est. Therefore, in the limit t → ∞ , the thermalizationdynamics of a given ensemble of particles depends exclu-sively on a ∝ (cid:80) Ni =1 γ i ∆ T i (0), or, in other words, the to- tal heat initially stored in the ensemble. This gives rise tointeresting behaviors, as we illustrate in Fig 2(a). There,we analyze the thermalization dynamics of a square arrayof N = 49 identical SiC particles with R = 25 nm and d = 4 R (solid curves). We consider different initial tem-perature distributions, all of them corresponding to thesame value of a . Specifically, the gray curve displaysthe time evolution of the temperature of the nanopar-ticles when all of them begin at ∆ T = 1 K. On theother hand, the colored curves represent different scenar-ios where only one particle, indicated in the schematicsusing the same color, is initially hot at ∆ T = 49 K. Onemight anticipate that when all of the particles in the ar-ray begin at ∆ T = 1 K, the system would thermalizemost quickly to the environment. However, as seen inFig 2(a), this is not the case. Instead, in all of the sce-narios under consideration, all of the particles approachthe equilibrium identically as e − λ t .Interestingly, for the scenarios in which only one par-ticle is initially hot, the thermalization process happensover two steps: first, all of the particles in the array con-verge to ∆ T = 1 K and, second, the whole array thermal-izes to the environment. This behavior is the result ofthe large difference between λ and the rest of the decayrates, as shown in Fig. S3. Therefore, if such differenceis decreased by, for instance, increasing the interparticledistance to d = 12 R , the two-step behavior fades away,as shown by the dashed curves.Although, so far, we have only considered ordered dis-tributions of particles, our conclusions apply to any arbi-trary ensemble of particles. For example, in Fig 2(b), weconsider an ensemble of N = 490 identical SiC nanopar-ticles with R = 25 nm randomly arranged within a spher-ical volume of radius 600 nm, as shown in the inset. Asin panel (a), we compare the thermalization process forfour different initial conditions; in three of them, oneparticle, marked in the schematics with the same coloras its corresponding curve, begins at ∆ T = 49 K, while,in the fourth (gray curve), all of the particles begin at∆ T = 0 . a takes the same valuefor all of the cases, they all approach the thermalizationto the environment identically, despite their very differ-ent initial temperature distributions.Another interesting scenario to consider is when theinitial distribution of temperatures is orthogonal to thefirst RHT eigenmode and hence a = 0. Physically, thismeans that, although the system is not thermalized, thetotal amount of heat initially stored in the ensemble iszero. In this case, the thermalization process is governedentirely by the eigenmodes describing the near-field RHTbetween the particles, since a net transfer of heat tothe environment (which is described by the first RHTeigenmode) is forbidden. To illustrate this, in Fig. 3(a),we study the thermalization dynamics of the array ofFig. 2(a) with d = 4 R , for the initial temperature distri-butions depicted in the insets. In both of them, one par- (a) -40-2002040 Δ T ( K ) -4 -2 -1 -3 -5 -6 t (s) N = 49 R = 25 nm d = 4 R (b) -2.0-1.00.01.02.0 Δ T ( K ) FIG. 3. (a) Thermalization dynamics for a hot (red curves)and a cold (blue curves) particle in a square array of N = 49SiC nanoparticles with R = 25 nm and d = 4 R . The red andblue particles are initially at ∆ T = 49 K and ∆ T = −
49 K,respectively, while the rest of the particles in the array are at∆ T = 0 K. We consider the two different cases depicted in theinset schematics, which are displayed with solid and dashedcurves, respectively. For comparison, the gray curve repre-sents the thermalization dynamics when all of the particles ofthe ensemble are initially at ∆ T = 49 /N K. (b) Zoom of thedata displayed in panel (a) around ∆ T = 0 K. ticle begins at ∆ T = 49 K and another at ∆ T = −
49 K,while the rest of the array is at ∆ T = 0 K, so a = 0.The corresponding results are displayed using solid anddashed curves, as indicated by the legend, with red andblue colors describing, respectively, the evolution of thetemperature of the hot and cold particles. As expected,in both cases, the thermalization of the array occurs ona time scale ∼ λ − ≈ − s. This is much faster thanthe thermalization when all of the nanoparticles begin at∆ T = 1 K (gray curve), even though, in that case, theparticles have to undergo a temperature change of only1 K (see panel (b) for a zoom around ∆ T = 0 K). Thereason is, again, the large difference between λ µ> and λ .Interestingly, the closer look provided in Fig. 3(b) re-veals an unituitive behavior: when the hot and cold par-ticles are placed next to each other (dashed curves), thetemperature of the initially cold particle rises beyond∆ T = 0 K and subsequently approaches it from above.We attribute this behavior to the difference in the lo-cal environment of the two nanoparticles; while the hot -5 -3 -2 -4 -6 -7 t (s) -50-2502550 Δ T ( K ) d d dd R R = 10 nm d = 4 R (a)(b) λ = 2.3 s -1 λ = 1.6×10 s -1 λ = 1.9×10 s -1 λ = 1.3×10 s -1 λ = 5.0 s -1 max-max Δ FIG. 4. (a) Thermalization dynamics for a chain with N = 5SiC nanoparticles arranged as shown in the upper schematics.We assume that R = 10 nm, d = 4 R , and the nanoparticlesare initially at (from left to right): ∆ T = 50, 50, 25, − −
50 K. (b) RHT eigenmodes of the chain and value ofthe associated decay rate. one lies on the corner of the array, the cold one is situ-ated in the interior and is therefore surrounded by moreparticles. This creates an imbalance in the cooling andheating rates of the two particles.We can use the RHT eigenmode framework to gainmore insight into this oscillatory behavior. To that end,we analyze a simpler system that exhibits similar oscil-latory dynamics but in a more pronounced way. In par-ticular, we consider the chain of N = 5 SiC nanoparti-cles with R = 10 nm and d = 4 R , shown in the upperschematics of Fig. 4. From left to right, the particlesare initially at ∆ T = 50, 50, 25, −
50, and −
50 K. Thedifferent curves in panel (a) show the temporal evolu-tion of the temperature of the particle with a matchingcolor. As the particles thermalize, their temperaturesoscillate around ∆ T = 0 K, with the center one (yel-low) crossing this value four times over the course of thethermalization process. The origin of this exotic behav-ior becomes clear by considering the RHT eigenmodes ofthe system, which are shown, with their correspondingdecay rates, in Fig. 4(b). Specifically, the initial stage ofthe thermalization is dominated by the eigenmode withthe largest decay rate, which corresponds to a near-fieldRHT process happening almost exclusively between thecenter nanoparticle and its nearest neighbor. After that,the contribution of the next fastest eigenmode drives thethermalization of both of those particles with their next-nearest neighbor. This pattern repeats with each suc-cessive eigenmode, resulting in the observed oscillatorybehavior of ∆ T .In conclusion, we have presented a theoretical frame-work to characterize the temporal dynamics of the near-field RHT in arbitrary ensembles of nanoparticles. Ourapproach is based on an eigenmode expansion of theequations that govern the RHT, obtained upon their lin-earization. The resulting set of eigenmodes completelycharacterize the RHT between the constituents of theensemble and their environment and therefore allow usto express, in a simple closed form, the temporal evolu-tion of the temperatures of the ensemble for any initialcondition. Exploiting this formalism, we have identifiedgeneral characteristics of the dynamics of RHT, whichoften present themselves in unintuitive ways. In par-ticular, we have shown that an ensemble of nanoparti-cles beginning with a fixed amount of stored heat alwaysapproaches thermalization identically, regardless of howthat heat is initially distributed. Similarly, when the to-tal initial heat stored in an ensemble is zero, the systemreaches thermal equilibrium faster than the case wherethere is any initially stored heat, no matter how large of a temperature change is required in the former casecompared to the latter. We have also predicted and ex-plained an exotic behavior in which the temperature ofnanoparticles oscillates around the equilibrium value asthey thermalize. Our results provide an insightful andcomputationally efficient approach to study the thermal-ization dynamics mediated by the near-field RHT, whichwill facilitate the systematic investigation of the impactthat novel phenomena, such as topology [37] and non-reciprocity [38, 39], have on this process. Furthermore,this framework can be exploited to analyze the combinedtransfer of energy and momentum mediated by the fluc-tuations of the electromagnetic field [40].This work has been sponsored by the U.S. NationalScience Foundation (Grant No. DMR-1941680) and theMinisterio de Ciencia, Innovaci´on y Universidades ofSpain (Grant TEM-FLU PID2019-109502GA-I00). L.Z.acknowledges support from the Department of EnergyComputational Science Graduate Fellowship (Grant No.DE-SC0020347). D.D. and W.K.K. acknowledge finan-cial support from the Laboratory Directed Research andDevelopment program of Los Alamos National Labora-tory under project LDRD 20210327ER. AppendixDerivation of Equation (1)
Here, we follow the approach from Refs. [25–27, 31] to derive Eq. (1). Within the dipolar approximation, the powerabsorbed by particle i in an ensemble with N elements can be written as P i = (cid:28) E i ( t ) · ∂ p i ( t ) ∂t (cid:29) , where E i is the electric field at the position of the particle, p i is its dipole moment, and (cid:104)(cid:105) stands for theaverage over thermal fluctuations. Shifting to the frequency domain through the Fourier transform defined as p i ( t ) = (cid:82) ∞−∞ dω π p i ( ω ) e − iωt for the dipole moment, and similarly for the field, P i can be rewritten as P i = − (cid:90) ∞−∞ dωdω (cid:48) (2 π ) e − i ( ω − ω (cid:48) ) t iω (cid:104) E ∗ i ( ω (cid:48) ) · p i ( ω ) (cid:105) , (4)where ∗ represents the complex conjugate. The electric field and the dipole moment appearing in this expression arethe self-consistent solutions of the many-body scattering problem for the ensemble with sources p fl and E fl . Thesesources are the fluctuating dipole and fields arising from the finite temperature of the particles and their environment.By solving this scattering problem, we can write E i and p i as p i = N (cid:88) j =1 (cid:2) A ij p fl j + B ij E fl j (cid:3) , E i = N (cid:88) j =1 (cid:2) C ij p fl j + D ij E fl j (cid:3) , (5)in terms of the following matrices with dimensions 3 N × N : A = [ I − α G ] − , B = A α , C = (cid:0) G + G (cid:1) A ,and D = I + C α . Here, I represents the identity matrix, α is a matrix that contains the polarizabilities of thenanoparticles, G = i k I , and G is the dipole-dipole interaction tensor. The components of G are zero for i = j and G ij = e ikd ij d ij (cid:2) ( kd ij ) + ikd ij − (cid:3) I × − e ikd ij d ij (cid:2) ( kd ij ) + 3 ikd ij − (cid:3) d ij d + ij d ij , for i (cid:54) = j , where d ij = r i − r j is the vector describing the distance between particles i and j , I × is the 3 × k = ω/c , and “+” represents the conjugate transpose.Substituting the solutions given in Eq. (5) into the expression of the power absorbed by dipole i shown in Eq. (4),we obtain P i = − (cid:90) ∞−∞ dωdω (cid:48) (2 π ) e − i ( ω − ω (cid:48) ) t iω N (cid:88) j,k =1 (cid:68)(cid:2) C ij p fl j + D ij E fl j (cid:3) + (cid:2) A ik p fl k + B ik E fl k (cid:3)(cid:69) . In order to perform the average over fluctuations, we use the fluctuation-dissipation theorem [41, 42] (FDT), whichtakes the form (cid:104) p fl i ( ω ) p fl+ j ( ω (cid:48) ) (cid:105) = 4 π ¯ hδ ( ω − ω (cid:48) )Im { χ i } δ ij (cid:20) n ( ω, T i ) + 12 (cid:21) for the dipole fluctuations and (cid:104) E fl i ( ω ) E fl+ j ( ω (cid:48) ) (cid:105) = 4 π ¯ hδ ( ω − ω (cid:48) )Im { G ij + G ij } (cid:20) n ( ω, T ) + 12 (cid:21) , for the electric field fluctuations. In these expressions, n ( ω, T i ) represents the Bose-Einstein distribution for temper-ature T i , with T being the temperature of the environment, and χ = α − G α + α . Then, using these expressionsand noting that any cross terms involving dipole and field fluctuations vanish, since they are uncorrelated, we obtain P i = (cid:90) ∞ dω N (cid:88) j =1 (cid:2) f ij n ( ω, T j ) + f ij n ( ω, T ) (cid:3) , where f ij = (2¯ hω/π )Tr (cid:2) Im { A ij Im { χ j } C + ij } (cid:3) and f ij = (2¯ hω/π )Tr (cid:104)(cid:80) Nj (cid:48) =1 Im { B ij Im { G jj (cid:48) + G jj (cid:48) } D + ij (cid:48) } (cid:105) , with thetrace taken over Cartesian components. Finally, since the power absorbed by particle i must vanish when thetemperatures of all particles are equal to T , regardless of the actual value of T , we have that f ij = − f ij , whichyields Eq. (1). Diagonalizability of the matrix H
As explained in the main text, the thermalization dynamics of an ensemble of nanoparticles, induced by the near-field radiative heat transfer (RHT), can be characterized by analyzing the matrix H . This matrix is the productof a positive definite diagonal matrix Γ − , whose entries are the inverse of the heat capacities γ i of the differentnanoparticles, and a symmetric matrix F . When all of the particles in the ensemble are identical, the matrix H is realand symmetric and therefore diagonalizable by the spectral theorem [43]. The situation is more complicated when theparticles in the array have different heat capacities. In this case, despite still being real, H is not symmetric because Γ − does not commute with F . However, H is still diagonalizable, as we show in the following. First, it is worthnoting that √ Γ exists and is symmetric because Γ is a diagonal positive definite matrix. Then, let us consider thesimilarity transformation √ ΓH √ Γ − = H S , (6)where we have introduced the matrix H S = √ Γ − F √ Γ − . Clearly, H S is symmetric because it is equal to itstranspose. Therefore, H S is diagonalizable and, consequently, has a complete set of eigenvalues λ µ and correspondingeigenvectors ∆ τ ( µ ) that satisfy H s ∆ τ ( µ ) = λ µ ∆ τ ( µ ) . Furthermore, from Eq. (6), it is clear that H has the same eigenvalues as H S and eigenvectors given by ∆ T ( µ ) = √ Γ − ∆ τ ( µ ) . Although, in general, these vectors are not orthogonal under the usual inner product, they are orthogonalusing an inner product weighted by the heat capacities N (cid:88) i =1 γ i ∆ T ( µ ) i ∆ T ( ν ) i = N (cid:88) i =1 ∆ τ ( µ ) i ∆ τ ( ν ) i = δ µν , -4 -2 -1 -3 -5 -6 t (s) -4 -2 -1 -3 -5 -6 t (s) -4 -2 -1 -3 -5 -6 t (s) -4 -2 -1 -3 -5 -6 t (s) Δ T ( t ) / Δ T ( ) Δ T ( t ) / Δ T ( ) Δ T ( t ) / Δ T ( ) R dd Rd N = 4 R = 25 nm d = 4 R ChainSquare(a)(b) (c)(d) (e)ChainChain ChainSquare Square
RHT Eigenmodes Δ T (0) = 10 K Δ T (0) = 50 K Δ T (0) = 100 K FIG. S1. (a-c) Comparison of the thermalization dynamics obtained using the near-field RHT eigenmode formalism (dashedcurves) with equivalent results obtained using the non-linearized full calculation (solid curves) for a chain of N = 4 identicalSiC nanoparticles with R = 25 nm and d = 4 R , as shown in the schematics. The different colored curves correspond todifferent values of ∆ T (0) for the particles, as indicated by the legend, while the different panels correspond to different initialtemperature distributions. Specifically, in panel (a) all of the particles begin at the same ∆ T (0), while in panels (b) and (c)only one of the particles (see inset) is initially hot and the rest are thermalized to the environment. In each case, we plot theevolution of the temperature of the particle(s) that is (are) initially hot. (d-e) Same as (a-c), but for a square arrangementof the particles. In all cases, the agreement between the RHT eigenmode formalism and the non-linearized full calculation isexcellent. where we have used the orthonormality of ∆ τ ( µ ) . Then, the solution of Eq. (2) is given by∆ T ( t ) = N (cid:88) µ =1 a µ e − λ µ t ∆ T ( µ ) , where the coefficients a µ are defined as a µ = N (cid:88) i =1 γ i ∆ T i (0)∆ T ( µ ) i . (a) Eigenmode λ μ ( s - ) max-max Δ (b)
900 1500 1800 2100120060030001015200525 λ = 2.19 s -1 λ = λ = λ = 3.53×10 s -1 λ = 2.49×10 s -1 λ = λ = λ = 2.47×10 s -1 N = 2197 R = 25 nm d = 4 R FIG. S2. (a) Decay rates of the RHT eigenmodes of a cubic array of N = 13 = 2197 identical SiC nanoparticles with R = 25 nmand d = 4 R . (b) First and last four RHT eigenmodes of the system and their corresponding decay rates. Eigenmodes 2 − − (a) max-max λ = 2.32 s -1 (b) λ = λ = 1.14×10 s -1 λ = 1.80×10 s -1 λ = λ = 1.74×10 s -1 Eigenmode λ μ ( s - )
15 25 30 35201050 40 45 5010152005 N = 49 R = 25 nm d = 4 R Δ FIG. S3. (a) Decay rates of the RHT eigenmodes of the system considered in Figs. 2 and 3, consisting of a square array of N = 49 identical SiC nanoparticles with R = 25 nm and d = 4 R . (b) First and last three RHT eigenmodes of the system andtheir corresponding decay rates. Eigenmodes 2 and 3, as well as 47 and 48, have a twofold degeneracy due to the symmetryof the ensemble. As in Fig. S2, the first eigenmode has equal amplitude across all particles, while, in the last one, nearestneighbors alternate signs and the amplitude decreases closer to the edges of the array. ∗ Corresponding author: [email protected][1] F. Reif,
Fundamentals of Statistical and Thermal Physics (McGraw-Hill, New York, 1965).[2] A. Narayanaswamy, S. Shen, L. Hu, X. Chen, and G. Chen, Appl. Phys. A , 357 (2009).[3] E. Rousseau, A. Siria, G. Jourdan, S. Volz, F. Comin, J. Chevrier, and J. J. Greffet, Nat. Photon. , 514 (2009).[4] R. S. Ottens, V. Quetschke, S. Wise, A. A. Alemi, R. Lundock, G. Mueller, D. H. Reitze, D. B. Tanner, and B. F. Whiting,Phys. Rev. Lett. , 014301 (2011).[5] R. St-Gelais, B. Guha, L. Zhu, S. Fan, and M. Lipson, Nano Lett. , 6971 (2014).[6] H. Chalabi, E. Hasman, and M. L. Brongersma, Phys. Rev. B , 014302 (2015).[7] K. Kim, B. Song, V. Fern´andez-Hurtado, W. Lee, W. Jeong, L. Cui, D. Thompson, J. Feist, M. T. H. Reid, F. J. Garc´ıa-Vidal, et al., Nature , 387 (2015).[8] B. Song, D. Thompson, A. Fiorino, Y. Ganjeh, P. Reddy, and E. Meyhofer, Nat. Nanotechnol. , 509 (2016).[9] R. St-Gelais, L. Zhu, S. Fan, and M. Lipson, Nat. Nanotechnol. , 515 (2016).[10] K. Shi, Y. Sun, Z. Chen, N. He, F. Bao, J. Evans, and S. He, Nano Lett. , 8082 (2019).[11] G. Domingues, S. Volz, K. Joulain, and J. J. Greffet, Phys. Rev. Lett. , 085901 (2005).[12] A. I. Volokitin and B. N. J. Persson, Rev. Mod. Phys. , 1291 (2007).[13] P. Ben-Abdallah, K. Joulain, J. Drevillon, and C. Le Goff, Phys. Rev. B , 075417 (2008).[14] A. Narayanaswamy and G. Chen, Phys. Rev. B , 075125 (2008).[15] G. V. Dedkov and A. A. Kyasov, J. Comput. Theor. Nanosci. , 2019 (2010).[16] A. Manjavacas and F. J. Garc´ıa de Abajo, Phys. Rev. B , 075466 (2012).[17] A. Manjavacas, S. Thongrattanasiri, J. J. Greffet, and F. J. Garc´ıa de Abajo, Appl. Phys. Lett. , 211102 (2014).[18] F. V. Ramirez, S. Shen, and A. J. H. McGaughey, Phys. Rev. B , 165427 (2017).[19] M. P. Bernardi, D. Milovich, and M. Francoeur, Nat. Commun. , 12900 (2016).[20] R. Yu, A. Manjavacas, and F. J. Garc´ıa de Abajo, Nat. Commun. , 2 (2017).[21] A. Fiorino, D. Thompson, L. Zhu, R. Mittapally, S.-A. Biehs, O. Bezencenet, N. El-Bondry, S. Bansropun, P. Ben-Abdallah,E. Meyhofer, et al., ACS Nano , 5774 (2018).[22] J. C. Cuevas and F. J. Garca-Vidal, ACS Photonics , 3896 (2018).[23] S.-A. Biehs, R. Messina, P. S. Venkataram, A. W. Rodriguez, J. C. Cuevas, and P. Ben-Abdallah, 0 , arXiv:2007.05604v1(2020).[24] D. Polder and M. Van Hove, Phys. Rev. B , 3303 (1971).[25] P. Ben-Abdallah, S. A. Biehs, and K. Joulain, Phys. Rev. Lett. , 114301 (2011).[26] M. Nikbakht, J. Appl. Phys. , 094307 (2014).[27] M. Nikbakht, EPL (Europhysics Letters) , 14004 (2015).[28] P. Ben-Abdallah, A. Belarouci, L. Frechette, and S.-A. Biehs, App. Phys. Lett. , 053109 (2015).[29] J. Dong, J. Zhao, and L. Liu, J. Quant. Spectrosc. Radiat. Transfer , 114 (2017).[30] J. Dong, J. Zhao, and L. Liu, Phys. Rev. B , 125411 (2017).[31] R. Messina, M. Tschikin, S.-A. Biehs, and P. Ben-Abdallah, Phys. Rev. B , 104307 (2013).[32] Y. Wang and J. Wu, AIP Adv. , 025104 (2016).[33] J. Song, L. Lu, B. Li, B. Zhang, R. Hu, X. Zhou, and Q. Cheng, Int. J. Heat Mass Transf. , 119346 (2020).[34] L. Zundel and A. Manjavacas, Phys. Rev. Applied , 054054 (2020).[35] V. Myroshnychenko, J. Rodr´ıguez-Fern´andez, I. Pastoriza-Santos, A. M. Funston, C. Novo, P. Mulvaney, L. M. Liz-Marz´an,and F. J. Garc´ıa de Abajo, Chem. Soc. Rev. , 1792 (2008).[36] E. D. Palik, Handbook of Optical Constants of Solids (Academic Press, San Diego, 1985).[37] A. Ott and S.-A. Biehs, Phys. Rev. B , 115417 (2020).[38] L. Zhu and S. Fan, Phys. Rev. Lett. , 134303 (2016).[39] A. Ott, R. Messina, P. Ben-Abdallah, and S.-A. Biehs, Appl. Phys. Lett. , 163105 (2019).[40] S. Sanders, W. J. M. Kort-Kamp, D. A. R. Dalvit, and A. Manjavacas, Commun. Phys. , 71 (2019).[41] S. M. Rytov, Theory of Electric Fluctuations and Thermal Radiation (Air Force Cambridge Research Center, Bedford,MA, 1959).[42] A. Manjavacas and F. J. Garc´ıa de Abajo, Phys. Rev. Lett. , 113601 (2010).[43] S. Friedberg, A. Insel, and L. Spence,