Near-field radiative transfer between two unequal sized spheres with large size disparities

: We compute near-ﬁeld radiative transfer between two spheres of unequal radii R 1 and R 2 such that R 2 (cid:46) 40 R 1 . For R 2 = 40 R 1 , the smallest gap to which we have been able to compute radiative transfer is d = 0 . 016 R 1 . To accomplish these computations, we have had to modify existing methods for computing near-ﬁeld radiative transfer between two spheres in the following ways: (1) exact calculations of coefﬁcients of vector translation theorem are replaced by approximations valid for the limit d (cid:28) R 1 , and (2) recursion relations for a normalized form of translation coefﬁcients are derived which enable us to replace computations of spherical Bessel and Hankel functions by computations of ratios of spherical Bessel or spherical Hankel functions. The results are then compared with the predictions of the modiﬁed proximity approximation.


Introduction
Theoretical and experimental investigation of near-field radiative enhancement between bodies has been an area of much study recently due to its possible application in a range of areas including thermal photovoltaics [1][2][3][4], thermal rectification [5], nanopatterning [6], thermally assisted magnetic recording [7], thermal imaging [8] and noncontact radiative cooling [9].The theoretical analysis has relied on Rytov's theory of fluctuational electrodynamics [10] and has been used to compute near-field radiative transfer between two flat surfaces [11][12][13][14][15][16], between two dipoles [17,18], and between a dipole and a flat surface [15,19].Only recently has this been applied to other more experimentally feasible geometries such as between two mesoscopic spheres [20] and between a mesoscopic sphere and a planar surface [21,22].The adjective mesoscopic is being used here to describe objects whose sizes are of the order of the characteristic thermal wavelength of radiation λ T (λ T ≈ hc/k B T , where 2π h is the Planck's constant, c is the speed of light in vacuum, k B is the Boltzmann constant, and T is the absolute temperature).In contrast, the size of macroscopic objects is much greater than λ T .The formalism for computing radiative transfer between arbitrary bodies has also been detailed in several recent publications [22][23][24][25][26]. Recent experimental validation of this near-field enhancement for sub-micron gaps has been limited to that between a mesoscopic sphere and a plate [27][28][29][30][31] and between a scanning probe tip and a planar surface [32][33][34].
The exact computation of the near-field thermal radiative energy transfer between two adjacent bodies involves expressing the electric and magnetic dyadic Green's function (DGF) [35,36] of the vector Helmholtz equation as a sum of properly weighted vector eigenfunctions (e.g.vector spherical waves in spherical coordinates or vector planar waves in Cartesian coordinates) of the vector Helmholtz equation.The weights of the vector eigenfunctions are determined by imposing the appropriate boundary conditions on the surfaces of the bodies.The boundary conditions account for multiple reflections and scattering of the outgoing thermal radiation from the bodies.The relation between the thermally generated electric currents, which are the source for thermal radiation, and the temperature of the bodies is given by the fluctuation-dissipation theorem [37][38][39].This procedure has been adopted to compute the thermal near-field radiative transfer between two spheres [20] and between a sphere and a plane [21].
Crucial to the analysis of radiative transfer between two spheres is the computation of coefficients of the translation addition theorem to facilitate the application of boundary conditions on the surface of the adjacent objects.In particular, the vector addition theorem allows us to transform vector spherical waves in a coordinate system with origin at the center of one of the spheres into vector spherical waves in a coordinate system with origin at the center of the other sphere.While there exist exact expressions to compute the coefficients of the vector addition theorem [40,41], recursion relations have to be used to ensure their efficient computation [42].
In addition to the exact method [20] for computing radiative transfer between spherical objects, which can be computationally expensive at small gaps, there exists an approximate method to predict the radiative transfer between nonplanar surfaces by treating them as a series of parallel surfaces with varying gaps.This method, called the modified proximity approximation method (MPA) [43] differs from the proximity force approximation (PFA) [44,45] used in Casimir physics in only the way the far-field contribution to the radiative transfer (which is negligible for forces) is treated.While there has been no rigorous proof in literature in deriving such a formulation, the validity of the PFA in different limits has been investigated recently for Casimir force [46] and the near-field radiative transfer [47] between a sphere and a plane.The MPA was first developed to ensure continuity in the radiative transfer between two equal sized spheres in the far-field and the near-field limits [43].Hence it would be of interest to also study the relevance of the MPA when the radius of one of the spheres R 2 is increased such that R 2 R 1 where R 1 is the radius of the smaller sphere, since when R 2 /R 1 → ∞ it reduces to the sphere-plane configuration that has been adopted for experimental verification.
Though the theoretical framework for the analysis of near-field radiative transfer between two spheres has been established by Narayanaswamy and Chen [20], it cannot be extended as is for unequal sized spheres with large radius ratio due to computational constraints.It is our intention in this work to develop simplifications to the computational model which would enable us to extend it to the case of disproportionately sized spheres.The paper has been organized as follows.In Sec. 2 we give a brief introduction to the translation addition theorem since it is central to our work and list the recurrence relations that have been developed in literature to compute the translation coefficients.We also propose a method to normalize the translation coefficients and show that the recursion relations for the computation of these normalized translation coefficients can be expressed in terms of ratios of spherical Bessel and Hankel functions.In Sec. 3, we propose a simplification for the translation addition theorem, which we have referred to as the one-term approximation, that aids the computation of near-field conductance between spherical bodies greatly.In Sec. 4, a description of the dependence of the normalized translation coefficients on the relative size of the two spheres is given which enables further simplification of the the computational procedure.In Sec. 5 we employ the simplifications listed above to compute near-field radiative transfer between two unequal sized spheres and compare this with the predictions from the MPA.We wish to stress that the approximations/methods developed here are have relevance not only to the computation of near-field radiative transfer but also to other problems in electromagnetic scattering, especially ones where electromagnetic surface waves and/or evanescent waves play an important role.

Translation addition theorems
Using the translation addition theorem, it is possible to transform the multipole expansion of an electromagnetic field from a coordinate system at the center of one scatterer to that of another in a multi-particle system.This becomes necessary while imposing the appropriate boundary conditions on the surface of the particles.Efficient computation of translation coefficients is essential, especially when the number of terms required for convergence in the multipole expansion is large, as is the case in near-field electromagnetic scattering problems.The simplest case of a multi-particle system is that of two spheres whose centers are translated along the zaxis as shown in Fig. 1.Emission and scattering of electromagnetic waves from such a system involves solving the vector Helmholtz equation given by: ∇ ∇ ∇ × ∇ ∇ ∇ × P P P( ( (r r r) ) ) − k 2 P P P( ( (r r r) ) ) = 0, ( where P P P( ( (r r r) ) ) is the electric or magnetic field at position vector r r r and k is the frequency dependent wave propagation constant.Following the notation in Fig. 1, k = k f = ω/c in the free-space region; k = k a = n a (ω) ω/c inside sphere 'a'; and k = k b = n b (ω) ω/c inside sphere 'b' with n a (ω) and n b (ω) being the complex refractive indices of sphere 'a' and sphere 'b' respectively.The independent divergence-free solutions of the vector Helmholtz equation in spherical coordinates are given by [20]: where M M M ( j) nm (kr r r) and N N N nm (kr r r) are vector spherical waves of order (n, m).n can take integer values from 0 to ∞.For each n, |m| ≤ n.The superscript p refers to the radial behavior of the waves.For j = 1, the M M M and N N N waves are regular waves and remain finite at the origin and z nm (θ , φ ) are vector spherical harmonics of order (n, m) and are given by: Following Ref. [42] these vector spherical waves can also be expressed in terms of the solutions to the scalar Helmholtz equation.The solutions of the scalar Helmholtz equation ∇ 2 ψ + k 2 ψ = 0 are given by: where, ψ nm are scalar spherical waves of the same order (n, m) and Y nm (θ , φ ) represent the spherical harmonics.Using these the vector spherical waves can be expressed as : and The representations for M M M and N N N waves in Eqs. ( 2), ( 3), (6), and (7) are equivalent.They differ from the corresponding representations in Ref. [42] in only a factor 1/ n(n + 1) which has nm (θ , φ ), and V V V nm (θ , φ ) are orthonormal to each other.
The scalar addition theorem for the scalar spherical waves ψ nm is given by: and the vector addition theorem for the vector spherical waves M M M ( j) nm and N N N nm waves is given by: where r r r and r r r are coordinates with respect to the coordinate systems centered at O and O respectively as shown in Fig. 1.Here, r r r = r r r + r r r .The coefficients of the vector addition theorem A νm ,nm and B νm ,nm can be expressed in terms of summation over Wigner 3j symbols [48], calculating which involve computations of a large number of factorials.For the configuration shown in Fig. 1 (with m = m due to axisymmetry condition), the vector translation coefficients are given by (using k = k f for the coordinate systems O and O positioned in free-space): where, with the symbol n ν p m µ q representing the Wigner 3j cooefficient.The above relations differ from those that are given in Ref. [41] due to the different definitions of M M M (p) nm (kr r r) and N N N (p) nm (kr r r) waves that we are using (given in Eqs. ( 6) and ( 7)).The corresponding scalar coefficient β νm,nm is given by (Eq. 25 in Ref. [41] modified to reflect axisymmetry condition for our configuration): In Eqs. ( 11), (12), and ( 14), the summation over p occurs in steps of 2, i.e., p When large number of coefficients are required (an expression for number of terms required for convergence in the summations of Eqs. ( 9) and ( 10) for radiative energy transfer calculations between two spheres has been derived in Ref. [49].We will show in Sec. 4 that n and ν ∼ 10 3 when d/R 2 ∼ 10 −3 ) computation of these Wigner 3j symbols becomes computationally tedious.Recursion relations derived for the scalar addition theorem [50] and the vector addition theorem [42] can be used, resulting in reduced computational times.The recursion relation for the scalar translation coefficient β is given by (Eq.23 and Eq.27 in Ref. [50]): and The form of the constants a + nm , a − nm and b + nm , b − nm have been retained from Ref. [50].Since the vector translation coefficients A νm,nm and B νm,nm can be related to β , they are computed directly from β as (Eq.12a and Eq.12b in Ref. [42] changed according to the axisymmetry of our configuration and definition of M M M and N N N waves): and To initialize the recursion for the scalar translation coefficient the following expression is used : where z ν (k f r ) represents the spherical Hankel function.When ν k f r , spherical Hankel and Bessel functions can be approximated by the following asymptotic forms: It follows from Eq. ( 20) and Eq. ( 21) that evaluating spherical Hankel and Bessel functions in the limit ν k f r can lead to numbers too large or too small for a given floating point format on a computer.For instance, for typical values of ν = 1000 and k f r = 200 that we encounter, z ν (k f r ) ≈ 8.3 × 10 567 and z (1) ν (k f r ) ≈ 3.05 × 10 −574 , while the maximum and minimum positive floating point numbers for double precision floating point format are 1.7976 × 10 +308 and 2.225 × 10 −308 respectively (both approximate).However, since products of such large and small terms are important for subsequent calculations (see for example, Eq. ( 25)), we need to modify the algorithm so that evaluations such large or small numbers can be minimized.This can be achieved through the "normalization" of the translation coefficients with appropriate factors, which we explain below.
A careful look at the form of the spectral radiative conductance (Eq.28 in Ref. [20]) gives us a hint about the appropriate factors.It is observed that coefficients D lM νm and D lN νm which represent the coefficients of M M M and N N N waves in the scattered field in vacuum (Eq.19 in Ref. [20]) are accompanied by factors of the form z Hence, what is important for the computation of the radiative energy transfer are not the terms D lM νm and D lN νm themselves but D lM νm z (1) only the term p = n + ν in the summation in Eq. ( 14) can be estimated to be: Using Eq. ( 20) in the limit p 1 we can show that : The error from Eq. ( 27) then reduces to For 1 % error we get n + ν ≈ 5k f r .For our computations we have used n + ν ≥ 7k f r as the criterion for employing the one-term approximation.From Fig. 2 it can be observed that when n+ν−2 (k f r ).The error (in %) in the spectral conductance at the surface phonon-polariton frequency, which is ≈ 0.061 eV for silica, (see Fig. 2 in Ref. [43]) when the one-term approximation for β νm,nm (k f r ) is used beyond n + ν = pk f r is shown in Fig. 2(b) for different values of n + ν.The errors are computed with respect to the spectral conductance value when the approximated form for translation coefficients is not employed.It can be observed that adopting n + ν = 7k f r as the criterion for employing the one-term approximation gives us an error of about 0.05 % in the spectral conductance.
n+ν−2 (k f r )) as a function of n+ν for k f r = 33 (arbitrarily chosen).The point n + ν = 7 k f r , beyond which the one-term approximation has been adopted in our computations for calculating the translation coefficients, is marked in in the figure.(b) The error in spectral conductance at the surface phonon-polariton frequency (0.061 eV) when different values of (n + ν)/k f r are chosen as the criteria for employing the one-term approximation.From the plot, a criterion n + ν = 7k f r is observed to give an error of ≈ 0.02% in the spectral conductance.The spectral conductance has been computed for two spheres of size R 1 = 13.7 µm and R 2 = 40R 1 with minimum gap d/R 1 = 0.01.
Since computation of Wigner 3j coefficient is computationally tedious, recursion relations can be derived for the approximated form for β νm,nm .Employing it is easy to show that the recursion relations turn out to be : The form of recursion given in Eq. ( 30) is computationally simpler than the corresponding recursion relation for the exact form of β νm,nm given in Eq. (15).Once the scalar translation coefficients β νm,nm are computed, the vector translation coefficients A νm,nm and B νm,nm can be got from Eqs. ( 17) and ( 18).

Dependence of normalized translation coefficients on the radius ratio of the spheres
For near-field scattering problems when the closest gap between the two spheres d R 1 , R 2 (here, d = r − R 1 − R 2 ), the normalized translation coefficients β N nm,νm show a marked dependence on the radius ratio R 2 /R 1 .This behavior can be used to improve the efficiency of computation of near-field quantities considerably.In the one-term approximation, β N νm,nm varies with n and ν mainly through the ratio: Figures 3(a), 3(b), 3(c) and 3(d) show the contour plots of of logarithm of the expression in Eq. ( 31) for two spheres of radius ratio R 2 /R 1 = 1, 3, 10, and 20, respectively, as a function of n and ν.The radius of the smaller sphere, R 1 , and the minimum gap between the spheres, d, are taken to be 10 µm and 50 nm for all the cases.When R 2 /R 1 = 1 it can be observed from Fig. 3(a) that the dominant terms in the matrix are present in a band along the diagonal with terms on either side of the band decreasing exponentially.This behavior has been utilized to simplify the computation of translation coefficients and the linear equations arising from imposing boundary conditions.Terms in the matrix with absolute value less than 10 −6 can be approximated to be zero and sparse routines can be employed to solve the linear equations.This approximation results in an error of less than 0.05 % in the final computed value of the conductance.For R 2 /R 1 > 1 the distribution of the dominant terms in the matrix changes markedly.From Fig. 3(c) and 3(d), where R 2 /R 1 1, it can be observed that only the first few rows of the matrix need be computed since the absolute values of the terms in the higher rows fall off exponentially.This substantially frees up memory resources while also reducing the computation time considerably.The errors from the two approximations (one term approximation, and making use of only the dominant (n, ν) terms in the normalized translation coefficient β N νm,nm ) are compared in Table 1.
The variation of β νm,nm with R 2 /R 1 can be determined from the asymptotic behavior of the expression in Eq. (31), where From Ref. [48], the asymptotic form for z (1) Similar asymptotic forms for z as a function of n and ν for two spheres with successive radius ratios R 2 /R 1 = (a) 1, (b) 3, (c) 10, and (d) 20 with R 1 = 10 µm, and the minimum gap maintained at 50 nm for all the cases.The dashed-lines denotes the contour line for a value of -6 which is taken as the cutoff point below which values for the normalized vector translation coefficients are approximated to zero.The line of maximum (shown as dotted lines) given by Eq. ( 35) has been superimposed on these contour plots or n+ν.When n 1/2, ν 1/2, the magnitude of the expression in Eq. ( 31) is approximately: When δ reduces to zero, employing cosh  41)) to predict the far-field contribution to the conductance demonstrate explicitly the role of translation coefficients in attaining convergence of near-field radiative transfer calculations between two spheres.For our numerical calculations we have used Eq. ( 39) as the criterion for the number of vector spherical waves since it is typically of higher value than that represented in Eq. ( 38).

Validity of modified proximity approximation
The conductance G can be calculated once the DGF for the two sphere configuration is determined using the theory developed in Ref. [20], modified suitably to account for the changes described in Secs.2-4.The radiative conductance G (units WK −1 ) between the two spheres shown in Fig. 1 is defined as: where T A and T B are the temperatures of the two spheres, and P (T A , T B ) is the rate of heat transfer between them.As explained in the introduction, it is also possible to predict the radiative transfer between nonplanar surfaces by using the MPA, i.e., by treating the nonplanar surfaces as a series of parallel surfaces with varying gaps, as shown in Fig. 4(a).The effectiveness of MPA in predicting near-field radiative conductance for two equally sized spheres has been shown in [43].Here, we investigate the relevance of MPA for the case when R 2 R 1 .
Using the MPA, the conductance G MPA between two unequal sized spheres of radii R 1 and R 2 for a minimum separation gap d and temperature T can be written as: where R 1 is the radius of the smaller sphere, is the gap at distance r from the symmetry axis as shown in Fig. 4(a), h n f (z) is the near-field contribution [43] to the radiative heat transfer coefficient between two half-spaces at gap z and G c (d, T ) can be approximated by the conductance value from classical radiative transfer theory when diffraction effects are negligible.G c (d, T ) for two unequal spheres of equal emissivity ε and temperature T is given by [54]: where F 12 is the gap dependent view factor between the two spheres and can be approximated (to an accuracy of 1%) by [55]: The emissivity for a silica half-space has been computed and plotted as a function of frequency in eV in Fig. 4(b).The conductance values have been computed by integrating over the frequency range 0.041 eV to 0.164 eV, and the value of G c (d, T ) has been appropriately adjusted to reflect this [43, see supplemental information].For the radius of the spheres that we have considered in our study (R 1 = 2.5 µm, 13.7 µm) the value of G(d) at d ≈ 40 µm (when near-field effects are negligible) is taken to be the value G c at that gap and the effect of gap-dependence of view factor for smaller gaps is included by using: The comparison between G and G MPA is shown in Fig. 5(a) for two spheres with R 1 = 13.7 µm and 2.5 µm and R 2 = 40R 1 .The error between G and G MPA is plotted in Fig. 5(b).For gaps d/R 1 < 0.1, MPA is observed to be able to predict the exact computed values of the conductance with errors less than 1%.
Since the total radiative conductance between the two spheres has contributions from frequencies where there is resonant enhancement from surface phonon polaritons as well as contributions from non-resonant frequencies where surface phonon polaritons are not active (see Fig. 2 in [43]), it would be of interest to observe if the MPA accurately predicts the contribution from both these regions.Computed values (by DGF formalism) and MPA predictions of the spectral conductance G ω at a resonant (0.061 eV) and a non-resonant (0.081 eV) frequency are plotted in Fig. 6(a) and Fig. 7(a).The contribution from the far-field radiative conductance G c (d, T ) in Eq. (41) has to be appropriately modified to reflect spectral conductance.From proximity approximation theory, as d/R 1 → 0, the spectral conductance G ω (d) for two spheres of unequal radii R 1 and R 2 is expected to scale as [56] : This characteristic R 1 /d behavior is observed for the resonant frequency contributions shown in Fig. 6(a) for d/R 1 0.02 where a slope of -1 is indicated.However such behavior is not observed for the non-resonant frequency contributions shown in Fig. 7(a).The error between G ω and the spectral conductance predicted by MPA, G MPA ω , for the resonant and non-resonant frequency contributions are shown in Fig. 6(b) and Fig. 7(b) respectively.In the far-field region (d/R 1 2 for R 1 = 13.7 µm, and d/R 1 10 for R 1 = 2.5 µm) where the enhancement due to tunneling of waves (surface waves at the resonant frequency and evanescent waves at non-resonant frequency) is negligible, the form of MPA in Eq. ( 41) predicts that the variation in G ω with gap is primarily due to the changing view factor between the spheres  d/R 1 0.07, for R 1 = 13.7 µm) the variation of G ω with gap is dependent on both the changing view factor with gap as well as increased tunneling of waves.For gaps d/R 1 0.07 (below which the view factor increases by less than 1%) the enhanced radiative transfer with decreasing gap can be attributed entirely due to increased tunneling of waves.At such small gaps MPA is able to model the enhancement at the resonant frequency within ≈ 5% errors irrespective of the value of R 1 , whereas at the non-resonant frequency the error is observed to be greater than 10% when R 1 = 2.5 µm.Despite such high errors at the non-resonant frequencies, MPA is successful in predicting the overall conductance when R 1 = 2.5 µm with error less than 1% for d/R 1 0.1 as observed in Fig. 5.This apparent discrepancy has been analyzed in detail in Ref. [57] and has been attributed to decreased contribution from non-resonant frequencies as R 1 decreases.
We summarize the main contributions from this work: (a) The exponential behavior of the spherical Hankel function z n (k f r ) when n ≥ k f r has been utilized to propose a simplified form, referred to in this work as the one-term approximation, for the coefficients in the translation addition theorem.The one-term approximation is valid when n ≥ Ck f r where C is a constant which is chosen depending on the desired accuracy.The recursion relations for these simplified translation forms are also given.They are simpler (b) A method to normalize the translation coefficients has been proposed and the behavior of the normalized translation coefficients on the relative size of the two spheres has been exploited to simplify the computations of near-field radiative transfer between two spheres.
(c) These simplifications are utilized to compute near-field radiative transfer between a mesoscopic and a macroscopic sphere with size ratio of 40 and we show that there is agreement with error less than 1% between the predictions of the modified proximity approximation and the exact computations for d/R 1 < 0.1.Thus the modified proximity approximation, which was

( 1 )
l (kr) is the spherical Bessel function of order l.For j = 3, the M M M and N N N waves are outgoing spherical waves that are singular at the origin and z

( 3 )
n (kr) is the spherical Hankel function of the first kind of order n.The radial function ζ

Fig. 1 .
Fig. 1.The configuration for this study consisting of two spheres of unequal radii R 1 and R 2 (labeled sphere 'a' and sphere 'b' respectively) and the relation between r r r , r r r and r r r .The surface to surface gap between the two spheres d is given by d = r − R 1 − R 2

Fig. 4 .
Fig. 4. (a) Application of MPA for the two spheres where the curved surfaces are approximated by a series of flat surfaces with varying gaps z.(b) the plot of spectral emissivity for a silica half-plane as a function of frequency in eV which is used in the form of MPA (Eq.(41)) to predict the far-field contribution to the conductance

Fig. 5 .
Fig. 5. (a) Plot of computed values of the total conductance (dotted) and the MPA (solid line)as a function of the non-dimensional gap d/R 1 for two spheres with R 2 /R 1 = 40 .The study has been performed for R 1 = 13.7µm and 2.5µm (b) The % error between the computed values and values from the MPA are plotted as a function of d/R 1

Fig. 6 .
Fig. 6.(a) Comparison between the computed values of G ω at a surface phonon-polariton frequency (0.061 eV) (dotted) and G MPA ω (solid line) as a function of the non-dimensional gap d/R 1 for two spheres of radius R 1 = 2.5 µm, 13.7 µm and R 2 = 40R 1 (b) The % error between G ω and G MPA ω as a function of d/R 1

Fig. 7 .
Fig. 7. (a) Comparison between the computed values of G ω at frequency 0.0801 eV (dotted) and G MPA ω (solid line) as a function of the non-dimensional gap d/R 1 for two spheres of radius R 1 = 2.5µm, 13.7µm and R 2 = 40R 1 (b) The % error between G ω and G MPA ω as a function of d/R 1