{'fraction_non_alphanumeric': 0.04040510048800965, 'fraction_numerical': 0.017526368263630163, 'mean_word_length': 4.592136150234742, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Thermal Discrete Dipole Approximation (T-DDA) for near-field radiative heat transfer simulations in three-dimensional arbitrary geometries

A novel numerical method called the Thermal Discrete Dipole Approximation (T-DDA) is proposed for modeling near-field radiative heat transfer in three-dimensional arbitrary geometries. The T-DDA is conceptually similar to the Discrete Dipole Approximation, except that the incident field originates from thermal oscillations of dipoles. The T-DDA is described in details in the paper, and the method is tested against exact results of radiative conductance between two spheres separated by a sub-wavelength vacuum gap. For all cases considered, the results calculated from the T-DDA are in good agreement with those from the analytical solution. When considering frequency-independent dielectric functions, it is observed that the number of sub-volumes required for convergence increases as the sphere permittivity increases. Additionally, simulations performed for two silica spheres of 0.5 micrometer-diameter show that the resonant modes are predicted accurately via the T-DDA. For separation gaps of 0.5 micrometer and 0.2 micrometer, the relative differences between the T-DDA and the exact results are 0.35% and 6.4%, respectively, when 552 sub-volumes are used to discretize a sphere. Finally, simulations are performed for two cubes of silica separated by a sub-wavelength gap. The results revealed that faster convergence is obtained when considering cubical objects rather than curved geometries. This work suggests that the T-DDA is a robust numerical approach that can be employed for solving a wide variety of near-field thermal radiation problems in three-dimensional geometries.


Introduction
Radiation heat transfer between bodies separated by distances greater than the dominant thermal wavelength is limited by Planck's blackbody distribution. In this far-field regime, radiative heat exchange predictions in three-dimensional (3D) complex geometries are tractable using wellestablished numerical techniques such as the discrete ordinates method and Monte Carlo approaches [1,2]. In the near-field regime of thermal radiation, which refers to the case where bodies are separated by sub-wavelength gaps, heat transfer can exceed by several orders of magnitude the blackbody limit [1,[3][4][5][6][7][8][9][10][11][12][13][14]. The enhancement beyond Planck's distribution is due to the extraneous contribution to energy transport by waves evanescently confined within a distance of about a wavelength normal to the surface of a thermal source. These modes include evanescent waves generated by total internal reflection of a propagating wave at the material-gap interface as well as resonant surface waves such as surface phonon-polaritons and surface plasmon-polaritons.
With the rapid advances in nanofabrication, near-field thermal radiation is becoming an important part of heat transfer engineering. Indeed, near-field thermal radiation may find application in imaging [43], thermophotovoltaic power generation [44][45][46][47][48], nanomanufacturing [49,50], thermal management of electronic devices [51], thermal rectification through a vacuum gap [52,53] and radiative property control [54][55][56] to name only a few. Due to these numerous potential applications, there is a need for predicting near-field heat exchange in 3D complex geometries. Numerical procedures, namely the finite-difference time-domain (FDTD) method [57][58][59], the finite-difference frequency-domain (FDFD) method [60] and the boundary element method (BEM) [61] have been applied recently to near-field thermal radiation calculations. Both FDTD and FDFD approaches suffer from large computational time, while the BEM is difficult to apply when dealing with heterogeneous materials. In this work, the Thermal Discrete Dipole Approximation (T-DDA) is proposed for simulating near-field heat transfer between 3D arbitrarily-shaped objects. The Discrete Dipole Approximation (DDA), extensively used for predicting electromagnetic wave scattering by particles, is based on discretizing objects into cubical sub-volumes behaving as electric point dipoles [62][63][64][65][66]. The T-DDA follows the same general procedure as the DDA, except that the incident field is induced by thermal fluctuations of dipoles instead of being produced by an external illumination.
The objective of this paper is to formulate the T-DDA and to test the method against exact results obtained from the DGF approach. In the next section, the physical and mathematical formulation of the problem is provided. Next, the T-DDA is derived starting from the stochastic Maxwell equations and the associated solution procedure is detailed. In the fourth section, the T-DDA is verified against exact results for two spheres separated by a sub-wavelength vacuum gap; a problem involving two cubes is presented afterwards. Concluding remarks are finally provided.

Physical and mathematical formulation of the problem
The problem under consideration is shown schematically in Fig. 1 where ε 0 and µ 0 are the vacuum permittivity and permeability, respectively, E is the electric field and H denotes the magnetic field.
The first moment of the thermal current is zero (i.e., 0 ) , ( = ω r J r ), which implies that the mean radiated fields are also zero [3]. In heat transfer analysis, the quantities of interest are not the mean radiated fields, but the flux and the energy density which are proportional to the ensemble average of the spatial correlation function of currents. This correlation function is provided by the fluctuation-dissipation theorem linking the thermal current to the local temperature of the emitter [3]: , where the superscript † indicates the Hermitian operator (conjugate transpose). In Eq. (2.2), I is the unit dyadic and Θ is the mean energy of an electromagnetic state given by [68]: where  is the reduced Planck's constant (= 1.0546 × 10 -34 J·s) while k B (= 1.3807× 10 -23 J·K -1 ) is the Boltzmann constant. Note that a factor four is included in the fluctuation-dissipation theorem to account for the fact that only positive frequencies are considered when passing from the time to the frequency domain [68].
In the next section, a general formulation is proposed for modeling near-field radiative heat transfer via the T-DDA.

Volume integral equation for radiation problems
The starting point of the T-DDA formulation is to decompose the (total) electric field E as the sum of an incident field, E inc , and a scattered field, E sca : The incident field can be interpreted as the field thermally generated by point sources that is propagating in free space in the absence of scatterers. A volume integral equation for the electric field E can be determined by first taking the curl of both sides of Eq. (2.1a): The vector wave equation is determined by substituting Eq. (3.2) into Eq. (2.1b): where k ( where k 0 is the magnitude of the wavevector in vacuum. In typical DDA formulations where particles are illuminated by an external source, the incident field satisfies In radiation heat transfer, the situation is different as the incident field is generated by 8 thermally fluctuating currents. As such, the thermally generated incident field propagating in free space satisfies the following vector wave equation: The vector wave equation for the scattered field can then be obtained by subtracting Eq. (3.5) from Eq. (3.4): The scattered field in Eq. (3.6) can be interpreted as the field generated by an equivalent source , propagating in free space. Solutions for Eqs. (3.5) and (3.6) are obtained using the free space DGF denoted by G [64]: where V e is the volume of the emitting bodies, while V (= V e + V a ) is the total volume where V a is the volume of the absorbing bodies. It can be seen in Eq. (3.7a) that the integration is performed over V e only where the thermal source is non-zero, while the integration in Eq. (3.7b) is performed over the total volume V to account for the interactions between all objects. The free space DGF is given by [69]: A volume integral equation for the electric field is obtained by substituting the scattered field given by Eq. (3.7b) into Eq. (3.1): The DGF has a singularity at r r ′ = , such that the principal value method is used to circumvent this problem. In this approach, an infinitesimal volume containing the singularity point is excluded from the integral. For a spherical or cubic exclusion volume, the application of the principal value method leads to [64,70]: where P.V. stands for principal value. The core equation underlying the T-DDA method is finally obtained by substituting Eq. (3.10) into Eq. (3.9): where the incident electric field E inc is specified by Eq. (3.7a). In the next section, Eq. (3.11) is discretized in order to derive a system of linear equations.

Discretization of the volume integral equation
The first step toward the numerical solution of Eq. (3.11) is the discretization of the emitting and absorbing objects into N cubical sub-volumes on a cubical lattice. The first N e sub-volumes are located in the emitters, while the N a (= N -N e ) remaining sub-volumes are allocated to the absorbers. The discretization should be fine enough such that the dimension of each sub-volume is small compared to the radiation wavelength (more details about the discretization are provided in section 4) [63]. If this condition is satisfied, it can be assumed that the electromagnetic properties and the electric field are uniform inside each sub-volume. Equation (3.11) evaluated at the center r i of a sub-volume i can therefore be written as: where the subscripts i and j refer to sub-volumes. In Eq. (3.12), when i ≠ j, the DGF has no singularity, such that the principal value can be approximated as [64]: . For i = j, the principal value integral is given by [64]: where a i is the effective radius of sub-volume i defined as (3∆V i /4π) 1/3 . The discretized version of the volume integral equation is obtained by substituting Eqs. (3.13) and (3.14) into Eq. (3.12): Equation (3.15) is a system of N vector equations where the electric field in each sub-volume is the unknown. It is important to keep in mind that E i is stochastic since the incident field is generated by random currents. The incident field, given by Eq. (3.7a), is approximated as follows after discretization into sub-volumes: Equation (3.16) stipulates that the thermally generated incident field is nil in the emitters (i = 1, 2, …, N e ), while the incident field within the absorbers (i = N e + 1, N e + 2, …, N) is due to the k = 1, 2, …, N e emitting dipoles.
From now on, it is assumed that each sub-volume is behaving as an electric point dipole, such that Eq. (3.15) can be re-written in terms of unknown equivalent dipole moments p i instead of unknown electric fields E i using the relation [64]: where CM i α is the Clausius-Mossotti polarizability defined as: The random current density r k J in Eq. (3.16) can also be expressed in terms of a thermally fluctuating dipole moment r k p representing thermal emission [69]: As for the random current, the mean of the thermally fluctuating dipole moment is equal to zero. (each dipole i has three orthogonal components) can also be written in a compact form using the following matrix notation: where P is the 3N stochastic column vector containing the unknown dipole moments p i , inc E is the 3N stochastic column vector containing the known incident fields E inc,i , while A is the 3N by Note that r ij is the magnitude of the distance vector r ij between dipoles i and j, while ij r is the unit vector along r ij .
For the self-interaction term (i = j), the sub-matrix ii A is given by: In the next section, the solution of the stochastic system of equations is discussed.

Heat transfer calculations
The main objective in heat transfer calculations is to compute the radiative power exchanged where the linear property of mean operator has been utilized. Equation (3.32) shows that the mean incident field is zero regardless of i, since the first moment of the thermally fluctuating dipoles is zero. As a result, the mean of the dipole moment vector, P , is equal to zero. Equation (3.30) is also used to calculate the correlation matrix of P . The correlation matrix of the zero-mean dipole moment vector is defined as [72,73]: where PP R is a 3N by 3N matrix consisting of N 2 sub-matrices: where EE R is the 3N by 3N correlation matrix of the incident field consisting of N 2 sub-matrices.
Using the fact that the mean value of the incident field is zero, EE R is calculated as: where the 3 by 3 sub-matrix,

Accuracy of the T-DDA
The accuracy of the T-DDA is subjected to the same validity criteria as the DDA. According to Draine [62], there are three validity criteria associated with the DDA. The first condition stipulates that the number of sub-volumes should be large enough in order to describe the geometries of objects accurately [62]. The error introduced by this effect is known as the shape error. A general quantitative criterion providing the minimum number of sub-volumes to maintain the shape error within an acceptable range has not been established. Yurkin et al. [74] showed that shape error for a specific sub-volume size is bounded by a summation of a i and a i 2 terms (a i is the effective dipole radius associated with a cubical sub-volume), but no information has been provided regarding the coefficients of these error terms. Draine [62] suggested a minimum number of sub-volumes for a sphere in the zero-frequency limit: where a eff is the effective radius of the object defined as The last condition is concerned with the fact that magnetic dipole effects are neglected when discretizing objects into sub-volumes behaving as electric point dipoles. Indeed, even for nonmagnetic materials, magnetic dipole absorption might be comparable, or even greater, than electric dipole absorption when dealing with conductive media such as metals [62,76]. The relative importance of magnetic dipole effects reduces with decreasing the size of the subvolumes [62]. Therefore, the T-DDA can be applied to conductive media provided that the discretization is fine enough. Draine [62] suggested a criterion that combines Eq.  The T-DDA is evaluated next for different material properties, separation gaps and object shapes.

Verification of the T-DDA
The T-DDA is first tested by computing the spectral radiative conductance between two spheres and by comparing the results against those obtained from the analytical solution of Narayanaswamy and Chen [33]. The thermal conductance between two objects is defined as follows. Assuming that the first object is at temperature T while the second object is at temperature T T δ + , the spectral conductance ω G is calculated as: is the net spectral heat transfer rate between the two objects. The case of two spheres of diameter D = 200 nm separated by a distance d =10 nm is analyzed in Fig. 2; the temperature is fixed at 400 K in all cases considered. The spheres of Fig. 2(a) Table 1 for an accuracy ∆ = 0.01. It can be seen in Table 1  [Insert Table 1 Table 1. A different behavior is observed in Fig. 2(b) for a larger dielectric function. As the number of sub-volumes increases, the T-DDA oscillates around the analytical solution. Figure 2(b) shows that the best results are obtained for 32 sub-volumes. This behavior is in agreement with the observations reported by Yurkin et al. [74], where results showed that the error does not monotonically decreases as the sub-volume size decreases. Instead, the error as a function of the spatial discretization displays a local minimum. For the largest dielectric function, an oscillatory behavior can again be observed in Fig. 2(c). The amplitude of these oscillations is larger than in the previous case of Fig. 2(b). There is a nearly perfect match between the T-DDA and the analytical solution when 136 sub-volumes are used. The comments made for Fig. 2(b) are therefore applicable to the case of Fig. 2(c).
[Insert Fig. 2 about here] The T-DDA is also applied for calculating the conductance between two spheres made of silica.
The diameters of the spheres are fixed at 0.5 µm while the temperature is kept at 400 K. The frequency-dependent dielectric function of silica has been extracted from the data reported in Ref. [77]. According to the criteria given by Eqs. for developing criteria more appropriate to near-field thermal radiation calculations; this is left as a future research effort as it is beyond the scope of this paper.
[Insert Fig. 3 about here] Finally, the spectral conductance between two cubes of silica calculated with the T-DDA is shown in Fig. 4. Both cubes have side lengths of 0.5 µm and are separated by a distance of 0.5 µm; as before, the temperature is fixed at 400 K. Note that an analytical solution does not exist for this problem. Additionally, Eq. (4.1) cannot be applied to cubes for estimating the minimum number of sub-volumes required to ensure an accuracy of ∆ = 0.01. Nevertheless, a fast convergence of the T-DDA is observed in Fig. 4. Indeed, the relative difference between the total conductances for 64 and 216 sub-volumes in each cube is 1.5%, while the relative difference is 0.6% between 216 and 512 sub-volumes (result not shown). This is in agreement with Yurkin et al. [74], where it was shown that the error associated with cubically shaped objects is smaller than for curved geometries.

Conclusions
A new method called the Thermal Discrete Dipole Approximation (T-DDA) was proposed for modeling near-field thermal radiation in three-dimensional arbitrary geometries. The T-DDA is conceptually similar to the Discrete Dipole Approximation, except that the incident field originates from thermal oscillations of dipoles rather than by an external radiation source. The T-DDA was verified against exact results for two spheres separated by a sub-wavelength gap for various dielectric functions and separation distances. In all cases considered, the results obtained from the T-DDA were in good agreement with the exact results. Moreover, the resonant frequencies arising between two silica spheres were predicted accurately. For this last case, the relative errors between the analytical and the T-DDA results of total conductance were 0.35% and 6.4% for two 0.5 µm-diameter spheres separated by gaps of 0.5 µm and 0.2 µm, respectively. The T-DDA results exhibited an oscillatory behavior when the dielectric function of the materials is far from unity. It was also observed that for a fixed number of sub-volumes, the error increases as the separation gap decreases. Additionally, faster convergence was achieved when considering two cubes rather than two spheres. This paper suggests that the T-DDA is a robust, relatively simple simulation tool for predicting near-field radiative heat exchange. The convergence and the accuracy of the T-DDA need further investigation, and this is left as a future research effort.