{'fraction_non_alphanumeric': 0.03231997165605409, 'fraction_numerical': 0.00902479111102363, 'mean_word_length': 4.926854876341578, 'pattern_counts': {'":': 0, '<': 2, '<?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': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fluctuational Electrodynamics in Atomic and Macroscopic Systems: van der Waals Interactions and Radiative Heat Transfer

We present an approach to describing fluctuational electrodynamic (FED) interactions, particularly van der Waals (vdW) interactions as well as radiative heat transfer (RHT), between material bodies of vastly different length scales, allowing for going between atomistic and continuum treatments of the response of each of these bodies as desired. Any local continuum description of electromagnetic (EM) response is compatible with our approach, while atomistic descriptions in our approach are based on effective electronic and nuclear oscillator degrees of freedom, encapsulating dissipation, short-range electronic correlations, and collective nuclear vibrations (phonons). While our previous works using this approach have focused on presenting novel results, this work focuses on the derivations underlying these methods. First, we show how the distinction between"atomic"and"macroscopic"bodies is ultimately somewhat arbitrary, as formulas for vdW free energies and RHT look very similar regardless of how the distinction is drawn. Next, we demonstrate that the atomistic description of material response in our approach yields EM interaction matrix elements which are expressed in terms of analytical formulas for compact bodies or semianalytical formulas based on Ewald summation for periodic media; we use this to compute vdW interaction free energies as well as RHT powers among small biological molecules in the presence of a metallic plate as well as between parallel graphene sheets in vacuum, showing strong deviations from conventional macroscopic theories due to the confluence of geometry, phonons, and EM retardation effects. Finally, we propose formulas for efficient computation of FED interactions among material bodies in which those that are treated atomistically as well as those treated through continuum methods may have arbitrary shapes, extending previous surface-integral techniques.


I. INTRODUCTION
Quantum and thermal fluctuations in electromagnetic (EM) fields are modified in the presence of polarizable objects. In thermal equilibrium, these fluctuating fields can transfer momentum, effecting van der Waals (vdW) or (synonymously) Casimir interactions, while out of thermal equilibrium, they can transfer energy, effecting thermal radiation and heat transfer between bodies. vdW interactions are of particular importance to molecular and low-dimensional structures both large and small, determining binding energies, stable conformations of polymorphic noncovalent crystals, and mechanical properties [1][2][3][4][5]. Recent studies of vdW interactions in molecular materials have illustrated the importance of modeling vdW interactions beyond the regime of pairwise additivity [6][7][8][9], which is valid only for isolated atoms/small molecules or (equivalently) sufficiently dilute bulk media, though these works have only considered distance regimes where the EM field response may be taken as the Coulomb potential in the electrostatic limit. On the other hand, theoretical studies of Casimir interactions among macroscopic bodies [10][11][12][13][14] have demonstrated nonmonotonic and repulsive forces among conducting objects even in vacuum at much larger distance scales where the speed of light (EM retardation) matters, but such continuum treatments are generally restricted to size and distance regimes large enough that continuous local empiricallyfitted dispersive susceptibilities accurately model the polarization response, so they are unable to accurately capture the atomistic nature and nonlocality (spatial dispersion) of the re-sponse of smaller molecular systems. Meanwhile, theoretical descriptions of radiative heat transfer (RHT) have been largely restricted to macroscopic bodies modeled with continuum local susceptibilities [15][16][17][18], demonstrating large enhancements as well as suppression factors arising from the tunneling of surface waves at short body separations, as compared to the predictions of the Planck blackbody law; in contrast to the case of vdW interactions, only a handful of investigations of RHT have focused on atomistic structures [19][20][21][22]. Comparatively more work has been pursued in the context of conductive heat transport by electrons and phonons [23][24][25][26], but even there, existing models tend to be fully atomistic and therefore restricted to treatments of small molecules or simple geometries with a high degree of (e.g. translational or rotational) symmetry. Since phonons and plasmons arise from and are influenced by EM interactions, respectively, fundamental questions remain surrounding the transition between radiative and conductive heat transfer at subnanometric gaps [20,[27][28][29].
Recently, we proposed a theoretical approach that conjoins atomistic treatments of molecular and low-dimensional structures with continuum treatments of macroscopic bodies in the context of fluctuational electrodynamics (FED) to enable description of EM fluctuation effects over a wide range of distance and geometric scales (from atom-to micron-scale gaps and from molecular to macroscopic media), including situations in which continuum approximations fail for a subset of the interacting bodies but not for others. We have called this the retarded many-body (RMB) framework of mesoscale FED. Our work has illustrated the importance of retardation effects in small molecular systems (where they are typically assumed to be negligible) and of geometry in determining the impact of collective, long-range EM fluctuations (i.e. polaritons) that cannot be appropriately captured by pairwiseadditive approximations [30]. Furthermore, we have shown that phonons in molecular structures can delocalize the polarization response of large molecules, leading for instance to nontrivial corrections to vdW interactions at room temperature relative to purely quantum fluctuations [31]. Similar consideration of phonons, nonlocal response, and longrange EM effects play a critical role in describing heat transfer among proximate molecules, with ab-initio atomistic modeling of the molecular response of materials enabling accurate descriptions of the transition from radiative to conductive heat transfer within the same unified theoretical framework [32]. A sample of results from these papers is in Fig. 1.
Early experiments on vdW interactions and RHT focused on measuring these phenomena in simple planar or spherical geometries validating predictions from continuum models of material response [33][34][35][36][37][38][39][40][41][42][43][44][45][46], while later experiments within the continuum domain have gone beyond such simple geometries [47][48][49][50][51][52]. More recent experiments in the atomistic and continuum domains of FED have begun to emerge, probing the edges of the regimes of validity of prior theoretical treatments, suggesting the need for new theoretical frameworks to better treat such multiscale problems. These experiments include measurements of vdW forces between organic molecules, macromolecular arrays, or single-layer sheets, and planar metallic or dielectric substrates without retardation [53][54][55], as well as between nanoparticles and nanotubes [56,57], that explore situations beyond the pairwise additive regime, measurements of Casimir-Polder forces on ground-state and Rydberg atoms, molecules, and Bose-Einstein condensates near planar substrates and gratings [58][59][60][61][62] where EM retardation is relevant, measurements of nearfield RHT between metallic tips and substrates at nanometric gaps [20,[27][28][29], and observations of thermal conductances in single-molecule junctions [63,64]. Such experiments are relevant to the engineering and operational understanding of molecular-scale devices [63,[65][66][67][68], heat management in electronic and thermophotovoltaic devices [52,[69][70][71], and manipulation of living cells and nanoparticles used in nanomedicine [1,[72][73][74][75], among other applications. Accurate explanation of all of these experimental results at many different length scales will require consideration of the interplay of phonons, retarded EM response, and complex geometric effects at the mesoscale, and suggests that our RMB framework may be well-suited to answer such questions. This paper accompanies a computational code which has been published as open source for others to use and extend, so the main goal of this paper is to provide rigorous derivations of the formulas underlying this code as well as our previous works [30][31][32]. In particular, we give detailed derivations of the most general formulation of mesoscopic FED, and show how our RMB approach combines sophisticated scat-tering and electronic calculation techniques in regimes where accurate atomic-scale descriptions of response are required for molecules while continuum permittivity models suffice for larger bulk objects can be exploited. Computational efficiency demands fast calculation of system matrices representing scattering among microscopic degrees of freedom (expanded in a basis of Gaussian functions): these matrices are found to greatly simplify into semi-analytical formulas involving Gaussian integrals in general, as well as Ewald summations in periodic media, thereby speeding up matrix assembly, so we provide detailed derivations of those formulas without assuming the absence of retardation or the validity of point dipolar approximations. The generality of our method allows for easy extension to other mesoscale EM phenomena of interest, including deterministic phenomena like absorption or scattering [76,77] and fluctuational phenomena like fluorescence [78], and while we do not focus on such phenomena in this paper nor implement such computational routines in our code, the open source nature of our code lowers barriers to pursuing this line of work in the future. Additionally, our method is general enough to consider macroscopic environments of arbitrary geometries and material properties, but our code and prior works [30][31][32] have almost exclusively focused on idealized perfect electrically conducting (PEC) planes as the archetypal macroscopic body or else have made further approximations involving compact molecules in the presence of more realistic macroscopic bodies. Thus, in this paper, we provide rigorous derivations of the extension of our method to treat arbitrary compact molecular and macroscopic bodies in conjunction with each other, relaxing those assumptions about the macroscopic body being a PEC plane. These derivations are based on the surface-integral formulation of Maxwell's equations [79][80][81], and while we do not computationally implement these formulas, our existing code as well as the code required for macroscopic computations [82] are both free and open source, making such a conjunction more feasible for future work.
Prior treatments of vdW interaction and RHT phenomena in macroscopic bodies have generally been related under the rubric of FED, allowing for exploitation of state-of-the-art classical computational EM techniques. These include finitedifference [17,18,[83][84][85][86], spectral [15,22,87,88] and Toperator [11,13,[89][90][91][92], surface integral [79][80][81], and volume integral [78] methods. All of these methods depend on primarily local empirical models for macroscopic susceptibilities and typically treat macroscopic objects as having hard boundaries, while accounting for EM retardation and scattering to all orders without making approximations about the smoothness of the object surfaces nor the diluteness of the media involved. Among these, only finite-difference time domain methods can handle spatially dispersive and potentially nonlinear polarizability response functions, and can handle any object geometry equally well; however, this comes at the cost of needing to discretize all of space and needing to step through time to a sufficient extent to obtain converged results, making this method computationally inefficient in most cases. The other methods are frequency domain methods, which precludes consideration of material nonlinearity, but each has its pros and cons beyond that; it is worth noting that all of the other aforementioned methods besides finitedifference have the advantage of discretizing only the degrees of freedom (DOFs) associated with each object, without needing to discretize the space in between. Spectral and T-operator methods converge most quickly for systems with continuous translational or rotational symmetries, but for arbitrary geometries, this convergence is drastically diminished. Surface and volume integral approaches can be implemented using localized rather than spectral basis functions, allowing for more efficient treatment of arbitrary geometries; volume integral methods require discretization of the volumes of every object, which is beneficial for objects with spatially varying susceptibilities or temperature gradients, while surface integral approaches are typically formulated to work only with homogeneous materials in each object. All of these methods can in principle handle linear spatially dispersive materials, but the susceptibilities describing those spatially dispersive materials tend to originate from phenomenological descriptions such as the hydrodynamic model, which cannot easily be applied outside of the simplest situations of spheres or planar substrates.
Atomistic treatments of heat transfer in molecular and larger systems have come in a few different flavors. Molecular dynamics [93][94][95][96] essentially solves Newton's equations of motion by treating each atom or molecule as a particle interacting with other such particles in a potential; this allows for conceptual simplicity, and the use of the time domain allows for treating anharmonic effects, but the treatment of noncovalent interactions among atoms and molecules is typically via crude pairwise additive approximations. Atomistic Green's function methods [19,20,[23][24][25][26] in the frequency domain allow for computation of heat transport by electrons or phonons, with harmonic oscillator models of electronic and phononic coupling often obtained from ab-initio calculations (while having the pitfall of being unable to treat anharmonicity); however, treatment of radiative heat transfer (i.e. via photons) is often neglected or heuristically approximated in a pairwise fashion, and apart from that, these methods tend to only be applicable to either small molecules or large bulk media with no ability to effectively bridge between the two disparate length scales. In the context of vdW interactions, recent methods [1-9, 74, 97, 98] have been successful by modeling electrons in molecules as effective oscillators coupled by long-range EM interactions, such that even though the underlying electronic model does not capture the inherent delocalization of electrons in metallic systems, long-range collective EM effects (i.e. polaritons) can be properly captured at zero temperature. We have been able to extend these methods [30][31][32] to include the effects of EM retardation, complex macroscopic geometries in the vicinity of molecules, and phonons in molecular systems as treated in an ab-initio manner, showing that especially at finite temperature, both vdW interactions and heat transfer are strongly influenced by the nonlocal response brought about in molecules by phonons, and that these interactions can be significantly modified by the presence of other large macroscopic bodies even when idealized as PEC planes.
This paper is organized as follows. After introducing matters of notation & convention and reviewing Maxwell's equations in Section II, we detail the general FED formulas for mesoscopic systems involving molecular and macroscopic bodies in Section III, starting with fully general formulas irrespective of the continuum approximation, and then showing how molecular and macroscopic DOFs can be partitioned. Following that, in Section IV, we briefly discuss how certain formulas simplify for compact molecules interacting in a fixed macroscopic environment, and the pitfalls therein, and demonstrate new results for interactions among small biological molecules in the presence of a PEC plane. We then use Section V to provide rigorous detailed derivations of the extensions of our general formulas to systems of infinite extent with Bloch periodicity, and discuss vdW interactions and RHT between two parallel graphene sheets in vacuum as an example of the versatility of our method. Finally, we describe how to extend our formalism to treat FED involving arbitrary compact molecules and macroscopic bodies in the surface integral equation formulation of Maxwell's equations in Section VI, and conclude in Section VII.

II. PRELIMINARIES
A vector field u i (x) will be denoted as |u ; we stress that the notation |u is a shorthand for a semiclassical field, not a second-quantized state. The standard conjugated inner product is defined as u|v = . In terms of this, the Hermitian and anti-Hermitian parts of a square operator (whose domain and range are the same), both of which are themselves Hermitian operators, are respectively defined as sym(A) = (A + A † )/2 and asym(A) = ( ; in terms of these, for square operators, Re(A) = (A + A ⋆ )/2 and x independent of basis. Once again, these are all in the position space representation; other bases can be used as convenient. Additionally, all quantities will be evaluated in the frequency domain, so the dependence on the frequency ω will generally be notationally suppressed for brevity and will only appear explicitly as needed for clarity.
FIG. 1. Diverse sample of past results in the RMB framework. Top-left, adapted from [30]: schematic of general systems that may be investigated in the RMB framework, involving atom-scale and macroscopic bodies. Bottom-left, adapted from [30]: vdW interaction energy (at T = 0) of a C500-fullerene or 250 atom-long carbyne wire in the parallel or perpendicular orientations with respect to a gold conical tip at distance z, relative to the interaction energy with respect to a semi-infinite planar gold slab at the same z. Middle, adapted from [32]: exponent (power law) of RHT power with respect to distance d separating two parallel 500 atom-long carbyne wires, each at distance z from a PEC plane. Middle inset: RHT power at each z, normalized to 4 times the equivalent blackbody emission power, as functions of d. Right, adapted from [31]: exponent (power law) of vdW interaction free energy at zero or room temperatures between an infinite graphene sheet and a parallel PEC plane with respect to the separation z, comparing RMB predictions (solid) to continuum RPA predictions without doping (dashed). Right inset: ratios of vdW interaction free energies at room to zero temperatures in RMB or RPA frameworks.
Maxwell's equations may be written in the frequency domain as describing the propagation of electric fields |E due to free polarization sources |P (0) in the presence of a set of polarizable bodies, which could each be low-dimensional atomistic structures or bulk media, with collective susceptibility V. The susceptibility is related to the permittivity via , and relates the total polarization density |P to the total electric field |E via |P = |P (0) + V|E . We focus solely on reciprocal media, in which the relation V = V ⊤ holds in position space ). Treating inherently nonreciprocal materials like topological insulators [99][100][101][102] (which could break reciprocity in the presence of an applied magnetic field) or other intrinsic nonreciprocal magneto-optic media [103,104] would require an extension of this formalism beyond the scope of this paper, but metamaterials that exhibit emergent nonreciprocal magneto-optic responses in an effective medium framework [105,106] can be treated using our framework if the underlying materials obey reciprocity. As these definitions are common in continuum EM theory but may be less familiar in the context of quantum chemistry, we point out three items of note. The first point is that the susceptibility in continuum EM theory is often denoted χ, especially for homogeneous, local, and isotropic material responses, and may more generally be written as χ ij (x, x ′ ). This contrasts with conventions in quantum chemistry, which define the charge density response as χ(x, , the latter of which is identical to our susceptibility V ij (x, x ′ ) [98]; strictly speaking, the charge density response that corresponds to the susceptibility is typically taken as the noninteracting charge density response. To avoid confusion, we exclusively use the notation V (and its position space representation) for the susceptibility. Readers may refer to Appendix B for more details. The second point, related to the first, is that in quantum chemical settings where retardation effects are unimportant, the induced polarization density ∆P(x) is of less interest than the bound charge density ρ B (x) = −∇ · ∆P(x), so the same physical bound charge density can be reproduced by shifting ∆P(x) by ∇ × Q(x) for an arbitrary gauge field Q(x): this is a reflection of the fact that in the absence of EM retardation, all electric fields are longitudinal and irrotational, so the addition of incompressible (solenoidal) vector fields cannot change the electrostatic properties of the system. However, such gauge freedom in ∆P(x) is lost when retardation is important, as can be seen by rewriting Maxwell's equations as ∇ × (∇×) − ω 2 c 2 I |E = ω 2 c 2 |P where |P = |P (0) + |∆P : the existence of transverse radiative (electrodynamic) fields destroys such electrostatic gauge invariance.
The third point, related to both of the prior two, is that it is more common in quantum chemical treatments of vdW interactions to speak of the free charge density ρ (0) (x) than the free polarization field P (0) (x) (and likewise the charge density response instead of the susceptibility). In general, the two are related in the frequency domain by − i ω∇ · P (0) (x) − i ωρ (0) (x) = 0. It is true that the free charge density becomes independent of the free polarization field exactly at ω = 0 (i.e. the static regime). However, our formulations of vdW interactions and thermal radiation depend on integrals over frequency in which the contribution at exactly ω = 0 is infinitesimal (and vanishes in the specific case of thermal radiation). For this reason and also to fully account for finite frequency effects (i.e. EM retardation) as well as anisotropy, we consistently use the free polarization field |P (0) instead of ρ (0) (x). For the same reason, we use the electric field |E and vacuum Maxwell Green's function G (0) , which are generalizations of the static potential φ(x) and Coulomb kernel v(x, x ′ ) = 1/(4π|x−x ′ |) common in quantum chemical treatments of vdW interactions, as the former two include far-field EM retardation effects. We again refer readers to Appendix B for more details.
Maxwell's equations (1) may be formally inverted to yield |E = G|P (0) , where we define the total Maxwell Green's function as the operator solving Maxwell's equations in the presence of all susceptibilities: We point out that the assumption of reciprocal media implies that G is reciprocal,

III. GENERAL SCATTERING AMONG MOLECULAR AND MACROSCOPIC STRUCTURES
In this section, we start with the most general formulation of EM scattering among molecular and macroscopic bodies in order to derive expressions for the vdW interactions and thermal radiation among collections of such bodies, which we specifically do in Section III A. These formulas depend only on the T-operators describing the EM scattering properties and response of individual bodies in isolation to all orders of scattering, and the vacuum Green's function propagating fields between pairs of bodies in a manner that only depends on the relative separations and orientations of the bodies. While these formulas are not new [11,92], they underscore the fact that molecular and macroscopic bodies can be treated together, on the same footing, in a unified formalism. In anticipation of our exposition of the computational details of the description of molecular and macroscopic DOFs, we then describe in Section III B how to equivalently rewrite the formulas for vdW interactions and thermal radiation by partitioning the total response of the system into molecular and macroscopic components. Finally, we give details about the basis representations of molecular and macroscopic response quantities in Section III C and Section III D, respectively, with further derivations of the expression of G (0) in the molecular basis in Section III E. We emphasize that although we focus in this paper on vdW interactions and thermal radiation, the EM scattering formalism is fully general, and the basis representation of molecular response can be applied to problems including those involving deterministic absorption or scattering, local density of states [76,77], or fluorescent emission [78], among others.
After this section, the following three sections each deal with a special case of the general formulas we present for vdW interactions and thermal radiation. The first case is when the macroscopic bodies do not change in separation or orientation relative to each other and when consideration of heat transfer may be restricted just between molecules. In the context of vdW interactions involving molecules, there might be only one macroscopic body present, like a thick metallic substrate or an atomic force microscopy (AFM) tip, in which case the question of relative displacements or orientations among multiple macroscopic bodies is moot. In the context of thermal radiation, consideration of energy exchange may be restricted to molecules if again only one macroscopic body is present, like an AFM tip, and it is in thermal equilibrium with its environment, while the molecules are maintained at a hotter temperature; this could be the case for measurements done at room temperature on biological molecules in hotter samples of living tissues or organisms. This special case allows for exploiting the EM field response (Green's function) of the collection macroscopic bodies without the molecules, which can be computed using a larger variety of methods that do not make reference to T-operators. The second case is an extension of the first case for extended molecular structures that obey spatial periodic boundary conditions. In that situation, we further derive analytical expressions for the vdW interaction energy and thermal radiation among molecular bodies in the presence of macroscopic bodies of commensurate periodicity, as well as analytical formulas for the expression of G (0) in the set of periodic molecular basis functions in a manner closely related to Ewald summation. Such a situation could arise, for example, when computing vdW interactions or thermal radiation for extended organic molecular crystals like aspirin in the vicinity of planar or periodically nanostructured metallic substrates. As we make clear in those sections, however, our code only implements these classes of computations for molecular bodies in vacuum or in the presence of a single PEC plane for computational simplicity. The third case, which is much more general, is for compact molecular and macroscopic bodies when the macroscopic bodies are characterized by spatially piecewise-constant permittivities, for which we may reformulate our method to exploit the surface integral equation (SIE) formulation of Maxwell's equations, as that yields significant computational benefits in arbitrary macroscopic geometries over more typical formulations, like spectral T-operator or volume integral equation (VIE) formulations. Such a situation could arise, for example, when considering vdW interactions or thermal radiation among a collection of proteins, polynucleotides, compact low-dimensional carbon allotropes, larger metallic nanoparticles, and an AFM tip, as may be relevant in more complex novel biomedical settings. That said, while we give the mathematical details of the method in this paper, we have not yet implemented this functionality in our code and leave that to future work.
A. Scattering, vdW interactions, and thermal radiation among general polarizable bodies To start, we consider a collection of N polarizable bodies labeled n ∈ {1, . . . , N } with susceptibilities V n . As we clarify later, we assume that the electronic structures and shortrange interaction properties of each polarizable body are unaffected by the presence of other bodies and that the bodies are otherwise spatially disjoint, so the the total susceptibility V = N n=1 V n may be written as a direct sum over the disjoint constituent susceptibilities, and V n = P n VP n is written in terms of the projection operators P n onto the polarizable material DOFs of body n; this means V is block-diagonal in the space of polarizable bodies. With this in mind, Maxwell's equations may be written in integral form as where |P (0) are free polarization sources in the polarizable bodies, while |E (0) refers to incident fields produced by sources outside of the system of polarizable bodies under consideration (so it does not include the lowest-order radiated fields G (0) |P (0) , which are already accounted in |P ). These equations can be self-consistently solved to yield where we define the T-operator of the total system as , describing scattering to all orders within and between all polarizable bodies; application of V −1 to |P (0) is allowed as the susceptibilities are nonsingular in the spaces spanned by the DOFs of the polarizable bodies. We note is exactly the fully interacting charge density response in the nonretarded approximation (under the random phase approximation), just as i, is the noninteracting charge density response. While we emphasize that these polarization, scattering, and radiation operators can be applied to a broad class of deterministic as well as stochastic EM problems, in this paper we particularly consider vdW interactions and thermal radiation. Both vdW interactions and thermal radiation arise from quantum and thermal fluctuations in the polarizations of material bodies. If the free polarization sources |P (0) and external incident fields |E (0) are taken to arise from quantum and thermal fluctuations, then their correlations are given through the fluctuation-dissipation theorem [59,107] (restoring explicit dependence on frequency for the sake of clarity) which relates fluctuations in free polarizations or ambient vacuum fields to dissipation quantities, respectively material absorption or free-space far-field radiation; these are defined in terms of the Planck factor Θ(ω, [108]. We note also that the fluctuating free polarization sources are uncorrelated from the ambient vacuum fields: We point out that this is a generalization of the fluctuation-dissipation theorem for the free charge density commonly used in quantum chemical treat- relates fluctuations in the free charge density ρ (0) (ω, x) to the dissipation given by the charge density response χ(ω, x, x ′ ) [109,110] (where it is worth noting that this is to be distinguished from V, which is often denoted as the susceptibility The total vdW free energy in a system of polarizable bodies at thermal equilibrium may be written in a Hellmann-Feynman form as the interaction between the total polarizations and fields [98,111]: where λ is the Hellmann-Feynman adiabatic connection parameter which linearly rescales V and G (0) , and where the expectation value . . . is taken over time, or equivalently over ensembles by ergodicity. By writing the polarizations and fields in the frequency domain (where in a slight abuse of notation, the same symbol is used for time and frequency domain quantities), we may then rewrite the total vdW energy as in a basis-independent manner. Using the results of (5) (where rescaling by λ is implicit for now), algebraic manipulations yield this into the formula for F tot allows for reduction of the integration to a single frequency variable instead of two, so henceforth frequency dependence will again be implicit in the notation. Additionally, restoring the factors of the adiabatic coupling coefficient λ means asym(G (0) T) = asym ∞ n=1 λ 2n (G (0) V) n , so the integration over λ can be done as´1 0 λ 2n−1 dλ = 1 2n , while Tr[asym(A)] = Im(Tr[A]) for any operator A means the imaginary part operation can be applied to the whole integral; it can then be seen that − This therefore allows for writing the total vdW free energy as where all quantities in the log-determinant expression depend on ω; we note that this formula reduces to the adiabatic connection formula for the vdW interaction energy in the nonretarded regime, performing an integration by parts using the definitions of the charge density response and G (0) in terms of V and the Coulomb kernel, respectively. Hence, there are two final steps needed to reach the desired expressions for the vdW interaction free energy. The first is that the interaction energy is the difference between the total energy in two different geometric configurations. If we maintain the assumption that the material properties of each polarizable body do not change with respect to geometric configuration, so that the susceptibilities V n only change trivially by virtue of rigid geometric transformations, then we can identify I − G (0) V = T −1 V, define T ∞ as the T-operator corresponding to each object in isolation (which by assumption only affects long-range EM scattering among the various bodies), and then rewrite the differences in the integrands between the desired configuration and the reference configuration of each object in isolation as ln(det[T ∞ T −1 ]). The second is that while the above frequency integral may be evaluated directly, it is analytically and numerically more desirable to perform a Wick rotation to positive imaginary frequency, where Θ(ω, T ) has simple poles whose residues may be evaluated easily [59,83]. Thus, we derive the vdW interaction free energy as where all quantities are evaluated at the Matsubara frequencies ω l = i ξ l for ξ l = 2πkBT l , and where the prime indicates a weight of 1/2 at l = 0 to avoid double-counting (as this was originally an integral over the entire real frequency axis). We point out that this Matsubara summation procedure is mathematically like a Riemann sum, and as T → 0, this sum converges to an integral. Indeed, as T → 0, the entropic contributions to the free energy vanish, and we recover the familiar expression for the zero-temperature vdW interaction though through the rest of this paper, we will use the notation F (T ) to denote the vdW interaction free energy of a given system at temperature T (so interactions at T = 0 will be denoted F (0)).
We now turn to thermal emission and RHT among polarizable bodies. Each of these phenomena can be described as the net work (with respect to relevant temperature differences) done by fields on the polarizations of one body labeled n, where those fields have been radiated by fluctuating sources in another body labeled m (which may be the same as n). Consequently, the ambient fluctuating vacuum fields |E (0) in (5) are irrelevant and may therefore be neglected, while the fluctuating sources are |P (0) m = P m |P (0) . To start, the total power may be written as via Poynting's theorem, where J(t, x) = ∂ ∂t P(t, x), and where the expectation value . . . may again be considered over time or equivalently over ensembles through ergodicity. Using the prior expressions for the Fourier transforms as well as the projection operators allows for rewriting in a basis-independent manner.
Following similar steps as with the vdW derivation, using the fact that V is block-diagonal and is therefore invertible in the space of material DOFs, algebraic manipulations yield where all quantities depend on ω. Plugging this in and again using the causality properties of the relevant response quantities to reduce the integral over frequencies to the positive axis allows for writing the power as We note that this only depends on the temperature of one of the bodies in question, and not that of another body or the ambient environment. In general, we can still define the dimensionless radiation spectrum from body m to body n at each frequency as and then, in terms of that, define general frequency integrated power quantities P =´∞ 0 W (ω) dω 2π , where W is defined as for thermal emission of body m into an environment of ambient temperature T env in terms of the sign function s nm = 1 − 2δ nm , or as for RHT between bodies m and n. In the context of thermal emission and RHT, as the Planck function Θ only appears in the form of differences at different temperatures, the zeropoint contribution ω/2 drops out, so it is helpful to redefine Θ(ω, T ) = ω/(exp( ω/(k B T )) − 1) without the zero-point term.

B. Partitioning molecular and macroscopic DOFs
None of the formulas in the prior subsection made particular reference to whether the polarizable bodies were atomistic or continuous, nor to any particular basis set, but as will become clear shortly, it is useful for the purposes of physical interpretation and computational convenience to introduce that distinction. That said, before proceeding, we clarify that the terms "molecular" and "macroscopic" are not absolute descriptors, but depend on the details of the configuration of polarizable bodies. As a general rule of thumb, bodies that are smaller than about 5 nm in at least one dimension or in feature size must be treated in an ab-initio manner incorporating atom-scale effects. We call such bodies "molecular" and use the label "mol" as a superscript or subscript associated with relevant response quantities, because prior work has focused on compact molecules and finite-size low-dimensional atomistic systems. As even certain structures of infinite extent in multiple dimensions must generally be treated atomistically, we refer to such structures as "molecular" for semantic consistency; as an example, graphene can arguably be visualized as a polycyclic aromatic hydrocarbon of infinite extent with no termination points where hydrogen atoms may lie. If none of the above conditions hold, then bodies may be well-described by coarse-grained continuum local bulk models of material susceptibility; even at such small length scales (e.g. metallic spherical nanoparticles), we term such bodies "macroscopic" due to the accuracy of bulk material modeling, and associate the label "mac" as a superscript or subscript associated with relevant response quantities. Even this categorization is not complete, because as the distances between proximate bodies fall below O(1 nm), each body would need to be treated atomistically; this would apply even to bulk metal substrates, at least with respect to the atoms closest to the other body. Thus, we assume the validity of continuum models for macroscopic bodies that are at least 5 nm along each dimension and feature and at least 1 nm away from any other body; if these conditions are violated, the "macroscopic" body would need to be treated atomistically, but this can anyway be done in our formalism.
Formally, we separate V = V mol + V mac as a sum of susceptibilities for disjoint collections of objects, where V mol = N mol s=1 V s and V mac = Nmac a=1 V a are each block-diagonal in their respective sets of DOFs. As described above, the assumption of disjointness will hold for any pair of bodies that are sufficiently separated and of appropriate dimensionality that short-range electronic exchange and correlation effects may be neglected for material DOFs between the two bodies. With this, we further define block 2 × 2 matrices via this separation of molecular (top block row and left block column) from macroscopic (bottom block row and right block column) DOFs: encodes scattering properties of the collection of molecules (macroscopic bodies) in a particular geometric configuration relative to each other in vacuum in the absence of macroscopic bodies (molecules), while the off-diagonal blocks G (0) propagate EM fields in vacuum between molecular and macroscopic DOFs. These formulas can simplify physical interpretation and computational implementation of various deterministic as well as stochastic EM phenomena involving molecules in conjunction with macroscopic bodies, though we specifically focus on vdW interactions as well as thermal emission and thermal radiation.
In general, the vdW interaction free energy or force may be desired in situations where one or more molecular or macroscopic bodies are taken together as a compound object; for example, if the force on an AFM tip in proximity with a graphene sheet adsorbed at a particular small separation to a metallic surface is desired, then the reference configuration would be the tip isolated from the graphene and metal surfaces, but the graphene sheet would remain at the same small separation from the adsorbent metallic surface. In such a case, the relevant reference configuration would not correspond to every body being isolated from each other in vacuum. In analogy to (18), we may define the block matrix T ∞ as the reference configuration of molecular and macroscopic bodies via where T mol(mac)∞ encode the positions, displacements, and orientations among only molecular or macroscopic bodies in the given reference configuration; for example, if the reference configuration is of all bodies infinitely separated, the representation of the off-diagonal blocks of G (0) will vanish. Ex-panding T ∞ and the determinant in the vdW summand blockwise leads to an expression for the summand where conceptually, the first term is the vdW interaction energy purely among macroscopic bodies in vacuum (in the absence of molecules) relative to their reference configuration, while the second term is the vdW interaction energy among molecular bodies in a scattering background created by the macroscopic bodies, relative to their reference configuration accounting for the change in the macroscopic bodies' positions and orientations from the corresponding reference configuration too. This formula has the additional benefit of making explicit the full interchangeability of molecular and macroscopic DOFs, as a fully mathematically equivalent formula arises simply by exchanging the labels mol ↔ mac and associated basis functions, showing how our formulation really does treat molecular and macroscopic DOFs on an equal footing; in particular, performing this exchange allows for writing the vdW interaction energy as the sum of that purely between molecular bodies in vacuum (in the absence of macroscopic bodies) relative to their reference configuration and the vdW interaction energy among macroscopic bodies in a scattering background created by the molecular bodies, relative to their reference configuration accounting for the change in the molecular bodies' positions and orientations too.
The formula for thermal emission and RHT in (15) holds for general molecular or macroscopic bodies, treating both sorts of bodies on the same footing and featuring the same benefits and pitfalls as (10). For the same reason, it may be more beneficial to explicitly separate molecular from macroscopic DOFs as in (18). In fact, the possibility of energy exchange between molecules and macroscopic bodies, going beyond energy exchange between molecules or macroscopic bodies alone, allows for such a separation to more clearly illustrate the richness of the mathematical formalism and computational & physical implications. Physically, heat transfer among macroscopic bodies in the presence of molecules can be realized via molecular junctions, heat transfer between a macroscopic body and a molecule could be realized via a metallic probe or nanoparticle locally heating a cancerous protein, and heat transfer between molecular bodies in the presence of macroscopic bodies could be seen in energy exchange between a hot graphene sheet and a cooler fullerene in the vicinity of a thick metallic substrate. Therefore, it behooves us to more fully draw out the formulas for thermal radiation in each such case. The radiation spectrum between just molecules m and n, in the presence of other molecules and macroscopic bodies, may be written as after using (18) and performing further operator manipulations, showing that the macroscopic bodies merely form a scattering background for energy exchange among molecular bodies. Likewise, the radiation spectrum between just macroscopic bodies m and n may be written exactly as (21) after exchanging the labels mol ↔ mac, showing again that the molecules merely form a scattering background for energy exchange among the macroscopic bodies. Finally, if m is a molecular body while n is a macroscopic body, the heat transfer may be written (with the reverse again obtained under the substitution mol ↔ mac) as: after manipulating operators and using the fact that the susceptibility operators are block-diagonal, so P n V = P n V n P n ; this expression clearly shows symmetry in the equation when the molecular and macroscopic bodies are interchanged. Note that while (6) and past T-operator and VIE formulations of thermal radiation make use of V, we choose to write our expressions in terms of V −1 as much as possible, because as we will shortly make clear, the molecular basis expansion we use directly gives V −1 mol without need for further inversion. The only thing remaining to describe scattering among molecular and macroscopic bodies is to represent T mol(mac) in appropriate basis sets and the off-diagonal blocks G (0) in the basis functions connecting molecular and macroscopic bodies. Such a representation will make the practical computational aspects and physical interpretations of formulas for vdW interactions and thermal radiation more clear.

C. Basis expansions of molecular DOFs
We write the molecular susceptibility as V mol = pi,qj α pi,qj |f pi f qj |. In general, the molecular susceptibility must account for the contributions of electrons, phonons, and other collective modes to the response; especially in metallic systems, this typically requires delocalized basis functions |f pi . However, for insulating or weakly conducting molecular systems, we may model the molecule as being made of nuclei that are harmonically coupled to nearest neighbors within each molecule and effective valence electronic harmonic oscillators associated 1-to-1 with a corresponding nucleus; this accurately captures the features of molecular response salient to fluctuational EM phenomena, like vdW interactions and thermal radiation, at ultraviolet frequencies via the effective electronic oscillators and at infrared frequencies via phonons arising from the coupled nuclear oscillators, and is valid for low temperatures where the harmonic approximation holds. In particular, the molecules together have N atoms labeled p located at positions r p , each of which has an effective electronic oscillator of charge q ep and mass m ep (which might not be equal to the fundamental electron charge or mass), coupled to its corresponding nucleus via an isotropic harmonic spring of constant k ep and damped isotropically with coefficient b ep = m ep γ ep (written in terms of a damping rate γ ep ), and a nucleus of mass m Ip coupled (in addition to its own electronic oscillator) to its nearest neighbors within each molecule via anisotropic spring constants K pq and damped isotropically with coefficient b Ip = m Ip γ Ip . The quantities q ep , m ep , k ep , K pq , and atomic coordinates r p are computed for each molecular body (or cluster, if a set of molecular bodies exhibits a more strongly correlated electronic structure even for nuclear separations beyond a few bond lengths) separately via density functional theory (DFT) calculations in conjunction with Hirshfeld partitioning [3][4][5]98], while m Ip is given from elemental data, and the damping rates γ ep and γ Ip are taken from empirical data. As the molecular DOFs are all damped coupled harmonic oscillators, with only the effective electronic oscillators directly coupling to electric fields (neglecting the nonlinear magnetic contribution to the Lorentz force, as may be done at typical operating temperatures as the relevant speeds of the material DOFs are nonrelativistic), the frequency domain equations of motion are simply where (Q e , M e , M I , B e , B I , K e , K I ) collect the parameters q ep , m ep , m Ip , b ep , b Ip , k ep , and K pq respectively into 3N × 3N matrices. These equations of motion determine the nuclear displacements x I and electronic dipole moments p e = Q e x e in response to an electric field e e obtained by evaluating |E at the atomic positions r p (leading to a 3Ndimensional vector); note that in this model, only the electronic oscillators directly couple to the electric field. Solving for p e = αe e gives the susceptibility matrix entering the basis expansion of V mol . The distinction between the ultraviolet contributions primarily from the electronic oscillators and the infrared contributions primarily from phonons arises due to M e and M I differing by 4 orders of magnitude, in contrast to the comparable magnitude of K e to K I . Additionally, as we model the electrons and nuclei as harmonic oscillators, we use Gaussian basis functions where the widths σ p , rather than being phenomenological, microscopically capture the nonlocal response of each molecule at each frequency by virtue of the definition σ p (ω) = (α p (ω)/3) 1/3 /(2 √ π) [30][31][32]112] in terms of α p (ω) = | q,j α pj,qj (ω)|/3. This choice of contracting the molecular susceptibility and averaging over the Cartesian tensor components to yield isotropic atomic fragment polarizabilities is consistent with previous expressions for isotropic local molecular susceptibilities used to construct Gaussian basis functions in the absence of phonons [4,8,30,98], and is also consistent with similar expressions deriving atomic polarizabilities from screened molecular susceptibilities [4,8,98] (Toperators, though those by definition include long-range EM interactions, unlike our bare expressions for V). Physically, this definition accounts not only for the change in the response at any given atom due to nonlocal internuclear couplings K I , but also for the full spatial extent of the nonlocality by summing over contributions from other atoms as well, though it does not explicitly preserve the anisotropy of the response in the Gaussian widths; the latter point, which could become especially important for low-dimensional materials like carbyne or graphene, is not further addressed in this work, but will be the subject of future work. Mathematically, we have found that while at imaginary frequency (relevant to vdW interactions) the polarizabilities α p (i ξ) will always be positive, at real frequency (relevant to thermal radiation and other EM scattering phenomena), the absolute value is necessary to ensure real positive Gaussian widths when constructing the basis functions, because the polarizability matrix α will in general be complex-valued and will have some diagonal or offdiagonal elements that have negative real parts at frequencies above electronic or phononic resonances. Additionally, the summation over other atoms q (as opposed to an alternative like | j α pj,pj (ω)|/3 which only accounts for the response at a given atom) is necessary to ensure positive-definiteness of asym(T) at real ω (or of T at ω = i ξ) in the molecular basis, though we have not been able to conclusively prove this statement. We point out that for most compact molecules as well as extended low-dimensional structures, the Gaussian widths will be less than O(5 nm).
We also require computation of the matrix elements f pi |G (0) f qj for atoms within and between molecules, in order to represent T mol . The use of Gaussian basis functions fortunately leads to analytical expressions for these matrix elements, which we first state and qualitatively discuss here, de-riving these expressions shortly afterwards. The expression is written analytically in terms of the dimensionless quantities q ≡ (ω/c) 2(σ 2 p + σ 2 q ) and ρ ≡ |r p − r q |/ 2(σ 2 p + σ 2 q ), thereby obviating the need for time-consuming numerical cubature over the volumes of the basis functions and in turn speeding up evaluation of the basis representation of T mol . We note that while our choice of basis functions |f pi is effectively a Galerkin discretization reminiscent of VIE formulations of Maxwell's equations, the number of basis functions is determined directly by the number of atoms as opposed to being chosen arbitrarily for numerical convergence. For nonzero Gaussian widths σ 2 p + σ 2 q , these matrix elements are finite even in the coincidence limit |r p − r q | → 0, unlike those of G (0) (ω, x, x ′ ), though the latter can be attained in the limit σ 2 p + σ 2 q → 0; this approach to a finite value captures the screening of long-range EM interactions due to short-range electronic response.

D. Basis expansions of macroscopic DOFs
For macroscopic bodies treated using continuous dielectric functions, while the expression T −1 mac = V −1 mac −G (0) can technically be used, it is not necessarily the most efficient way to obtain a basis representation {|b β } for T mac . In particular, while VIE methods do use this expression in conjunction with localized voxel or Schaubert-Wilton-Glisson basis functions {|b β } to represent T mac , other methods like scattering methods in planar or spherical waves, or finite-difference methods, may represent T mac in a way that is mathematically equivalent, but less obviously so, to the above expression. In any case, the choice of macroscopic basis |b β will also affect the computation and convergence properties of the representation of G (0) as b β |G (0) f qj connecting molecular and macroscopic DOFs; the term G (0) |f qj is evaluated analytically in position space in the same way as f pi |G (0) f qj but in the limit σ p → 0 and with r p replaced by a generic x, as those limits applied to |f pi yield a Dirac delta function, while the convergent, smooth, analytic properties of G (0) |f qj facilitate analytical or numerical evaluation of the matrix elements b β |G (0) f qj . Thus, the choice of macroscopic basis |b β should in practice account for the convergence properties of both T mac and b β |G (0) f qj .

E. Expression of G (0) in the molecular basis
The following is a brief digression deriving (26), which will be beneficial to demonstrate how analytical expres-sions for (26) exist without needing high-dimensional numerical cubature, and to later extend similar formulas in the particular case of periodic molecular structures. The derivations below assume ω = i ξ, so evaluation at real ω can be obtained by substituting ξ = − i ω in the results at the end. We define G 0ij (i ξ, x, This means the inner product may be written as´´f Performing integrations by parts given vanishing surface terms to put the derivatives on f p , noting the form of f p as dependent only on |x − x p | allows for writing ∂ j f p (x) = −∂ rpj f p (x), and bringing the derivatives with respect to the Gaussian basis function centers r p outside of the integrals over x and x ′ allows for rewriting the inner product as To simplify this calculation, it is necessary to write g 0 using the following integral representation: in which case the inner product simply turns into a set of Gaussian integrals. This can be seen in the product, f p (x)g 0 (iξ, x, x ′ )f q (x ′ ), whose spatial dependence comes only in the exponential term, the exponent of which can be written as − 1 2 X ⊤ AX + J ⊤ X given the 6-dimensional vectors and tensor, Exploiting the well-known identity of multivariate Gaussian integration,ˆe and accounting for the remaining factors in the multiplication f p g 0 f q , one obtains in terms of an integral over the auxiliary variable u. This integral may be evaluated directly through use of a few variable substitutions. The first is to transform to dimension- , and . This transforms the integral into such that all dimensional terms are prefactors of the integral, which itself is dimensionless. The second is to transform to w = v √ 1+v 2 so that the semi-infinite integration range is mapped to the finite interval [0, 1], and (1+v 2 ) −3/2 dv = dw, yielding the integral after direct evaluation. It is this expression that is finally used to obtain (26) at real frequency (after substituting ξ = − i ω).

IV. COMPACT MOLECULES IN A FIXED MACROSCOPIC ENVIRONMENT
For the case of vdW interactions, we start by considering (20) for the case where the macroscopic bodies (if there are multiple) are fixed relative to each other, so that T mac = T mac∞ . Typically, the energy differences we choose to measure are set to be relative to a configuration where the molecules are infinitely separated from the macroscopic bodies, so the off-diagonal blocks between molecules and macroscopic bodies in the reference configuration satisfy G 0∞ → 0. This allows for simplification of the integrand to ln(det(T mol∞ ) det(T −1 mol − G (0) T mac G (0) )). At this point, we may define as the EM field response in the presence of only the fixed macroscopic bodies, and as an effective T-operator encoding the scattering properties among the molecules in a modified EM environment due to the presence of fixed macroscopic bodies in the background. This allows for rewriting the integrand as ln(det(T mol∞ T ′−1 mol )). For the case of thermal radiation, the definitions in (33) and (34) can only be used to simplify (21) in the case where the labels m and n are only for molecules; it turns out that a fuller consideration of macroscopic DOFs is required when at least one of m or n is a macroscopic body, so (33) and (34) are insufficient in those cases. Thus, if we focus on the case of heat exchange only among molecules, the macroscopic bodies again form a fixed background that only act to modify the field response experienced by the molecules, so (21) is changed to yield as the thermal energy exchange among molecules.
Our definitions (33) and (34) are useful because T mac is stipulated to be fixed given that the macroscopic bodies will never change in separation or orientation relative to each other, so the field response G mac may be computed using a much broader range of computational methods, such as finitedifference or multipole methods in addition to spectral or VIE T-operator methods, in which T mac by itself may be practically more difficult to extract. However, this benefit can also be seen as a pitfall in itself: in almost every situation where G mac can be computed analytically or numerically, the matrix elements f pi |G mac f qj will generally require slow 6-dimensional numerical cubature, as the Gaussian widths defining the basis functions are not guaranteed to be small enough at any given frequency and geometric configuration to be approximated as point dipoles. The only exception is if G mac can be analytically written in terms of G (0) , as is true, for example, in the case of a PEC plane thanks to image theory, which can be used as a good approximation for a thick planar metallic substrate at frequencies below the ultraviolet; this would allow for using the analytical formulas of (26), avoiding the need for costly numerical cubature. (An approximate exception, seen in Fig. 1 (bottom-left), comes from [30] where we ignored the effective nuclear DOFs: this yielded much smaller Gaussian basis function widths, so we approximated G mac − G (0) above a gold plate or a gold cone using numerical techniques by approximating the Gaussian basis functions as Dirac delta functions. However, this approximation is not valid when nuclear DOFs are considered.) Thus, our code can currently only treat molecular bodies interacting either in vacuum or in the presence of a single PEC plane. This is the implementation we have used for our past works [30][31][32] which primarily consider compact molecules above a single PEC plane, and we point readers to those works for more detailed discussions of specific example systems. We also point out that the Gaussian basis functions |f pi describing the molecular DOFs do not have compact support, which means those basis functions will nontrivially overlap with a PEC plane (or equivalently with their images on the other side of the plane) if present, potentially leading to unphysical results; this is not problematic for Gaussian widths smaller than 1 nm as we never consider smaller separations anyway, while for molecular bodies with large Gaussian widths especially at smaller separations, even if the results are not rigorously justifiable, we retain them as a useful approximation to the short-range EM interaction effects of the molecular body near a metallic surface. As a simple example of interactions involving small com-pact biological molecules, in Fig. 2, we consider (a) the vdW interaction free energy F (T ), or (b) the RHT coefficient dP dT , between the nucleotides guanine and cytosine as functions of the orientation of cytosine, given by the clockwise rotation angle ϕ about the z-axis through the center of mass of cytosine, in which the two molecules are displaced from each other horizontally by a distance d = 1 nm, and both are displaced vertically by the same distance z above a PEC plane; all calculations are done at the normal human body temperature T = 310 K. We compare F (T ) and dP dT for different values of z: each quantity changes by much less than 1% when z is increased beyond 100 nm, so we only consider z ∈ {1 nm, 100 nm}. The vdW interaction free energy shows clear differences at each z, indicating that there is a significant contribution from the vertical force by the PEC plane to the overall interaction for z ≤ 100 nm. However, this is completely independent of ϕ, because if the two functions of ϕ are overlaid upon each other to have the same value at ϕ = 0, they consistently remain well within 1% of each other for all ϕ. This means that the vdW torque − ∂F ∂ϕ at each ϕ is essentially independent of z; at any z, for that value of d and initial orientation of molecules, there are two stable and two unstable equilibria for the vdW torque. Meanwhile, the RHT coefficient is likewise essentially independent of z, as is clear from the figure. The apparent independence of these quantities from the distance to a PEC plane is due in both cases to consideration of the interactions between two small, compact, chemically heterogeneous molecules of complicated shapes, yielding weaker polarization responses, as opposed to the interactions between a low-dimensional compact or extended low-dimensional carbon allotropes of simple high-symmetry shapes, which would yield stronger polarization responses. This has previously been observed in comparisons between the interactions with a metal plate of low-dimensional carbon allotropes versus complicated proteins [30].
As a more complex illustrative example leading to nontrivial interaction behaviors, in Fig. 3, we consider (a) the vdW interaction torque, or (b) the RHT coefficient dP dT , of a complicated four-body system involving guanine & cytosine along with a 250 atom-long carbyne wire, all of which lie above a PEC plane. The two small molecules' centers of mass, as well as the bottom of the carbyne wire, are fixed at a distance z above the PEC plane, while the angle ϕ that the wire makes with respect to the horizontal plane is varied. At ϕ = 0, the wire comes exactly in between the two small molecules, such that the centers of mass of these three molecules lie along a line perpendicular to the wire, with each small molecule's center of mass a distance 2d = 2 nm away from the other. Effectively, the wire can be thought of as a "switch" to be lifted, as we consider variations in the torque and heat transfer coefficient with respect to the orientation ϕ for two different values of z, namely z = 1 nm or z = 100 nm. All of these quantities are computed at the normal human body temperature T = 310 K.
The behavior of the torque is dominated by the interactions of the wire with the PEC plane, and guanine and cytosine, In both calculations, the bottom atom of the wire and the centers of mass of guanine and cytosine are at a distance z, namely 1 nm (blue) or 100 nm (red), above a PEC plane, and the wire may be oriented with respect to the horizontal axis at an angle ϕ. Furthermore, at ϕ = 0 (corresponding to the wire lying parallel to the PEC plane), the centers of mass of the wire, guanine, and cytosine lie in a line perpendicular to the wire axis and parallel to the PEC plane, with guanine and cytosine each lying at d = 1 nm on opposite sides of the projection of the wire onto the horizontal plane. Both calculations are at T = 310 K. Top inset: schematic of the carbyne wire at an angle ϕ above the PEC plane. Bottom inset: plan (top) view of the wire, guanine, and cytosine. In both schematics, the wire is not to scale, but everything else is to scale. being small molecules with no obvious symmetries & weak polarizabilities, only have a marginal impact on the interaction torque. In particular, the configuration at z = 100 nm is far enough from the PEC plane that the torque is always very close to zero, as an isolated carbyne wire in free space should not exhibit any torque at all; for ϕ ≥ 15 degrees (i.e. even when the wire is somewhat close to parallel to the PEC plane), the wire is far enough from the small molecules that the torque is effectively negligible, while it only has a magnitude of 6.2 × 10 −22 N · m (with the negative sign indicating attraction to the PEC plane and other molecules) even at the parallel orientation ϕ = 0. By contrast, at z = 1 nm, the proximate PEC plane produces a much stronger attraction even at relatively large ϕ (i.e. when the wire is closer to perpendicular to the PEC plane), and the attractive torque of 6×10 −20 N · m at ϕ = 0 is nearly 2 orders of magnitude larger than the corresponding torque at z = 100 nm.
Qualitatively subtler effects emerge when considering the heat transfer coefficient specifically between guanine and cytosine, because in this situation, the wire and the PEC plane both act only to modify the environmental EM scattering properties, with the wire orientation angle ϕ further varying this; thus, aggregate many-body effects must be considered. At z = 100 nm, the wire being parallel to the PEC plane and lying between the two small molecules facilitates heat transfer between the two small molecules compared to when it is perpendicular to the PEC plane (such that the small molecules are exchanging energy effectively in vacuum), by virtue of modifying the EM scattering properties: the heat transfer coefficient decreases from 8.96 × 10 −13 W/K at ϕ = 0 to 3.73 × 10 −13 W/K at ϕ = 90 degrees. Qualitatively similar behavior is observed at z = 1 nm, but screening from the proximate PEC plane (i.e. from the image dipoles of each molecule) consistently decreases the heat transfer coefficients, and this decrease is not uniform with respect to ϕ, as the ratio of the heat transfer coefficient at z = 1 nm to its counterpart at z = 100 nm decreases nonmonotonically from 0.88 at ϕ = 0 to 0.8 at ϕ = 90 degrees.
While our past works [30][31][32] have focused on power laws or ratios of vdW interaction free energies or RHT powers with respect to distances, and have focused on high-symmetry lowdimensional carbon allotropes present individually or in pairs in vacuum or near a PEC plane, here we show the greater generality of our framework in treating compact molecules, modeling FED phenomena in more complex many-body systems of carbon allotropes and biological molecules in the presence of a PEC plane, and considering rotational dependence in addition to distance dependence. We expect that our framework, as given, may have fruitful applications to predictions of highly nontrivial vdW interaction energy and RHT power behaviors in even more complex systems of compact molecules near metallic substrates as experiments become more sensitive.

V. EXTENDED MOLECULAR STRUCTURES IN AN ARBITRARY ENVIRONMENT
In this section, we extend expressions from the previous sections for the vdW interaction energy and heat exchange among a collection of molecular and macroscopic bodies to consider spatially extended geometries with commensurate spatial periodicity. The imposition of Bloch periodicity leads to nontrivial expressions of the fluctuation-dissipation theorem as well as novel formulas extending (26) to spatially extended structures. Given this, we first review the definitions and relations among EM fields, polarizations, and response functions in periodic geometries, and derive relevant FED formulas taking care with the nontrivial changes imposed by Bloch periodicity in Section V A. Then, we exploit Ewald summations and integral techniques to derive fast-converging semianalytical expressions for the matrix elements that describe scattering among periodic molecular basis functions in vacuum in Section V B. Finally, we consider examples of relevant atomistic systems exhibiting periodic boundary conditions to which this formalism may be applied, particularly vdW interactions and RHT between two parallel graphene sheets in vacuum, in Section V C as a demonstration of the versatility of our method.

A. Scattering, vdW interactions, and thermal radiation among polarizable bodies with periodic boundaries
Consider a collection of polarizable bodies labeled n ∈ {1, . . . , N } that all obey periodic boundary conditions given by lattice vectors R; these polarizable bodies will first be treated in a fully general manner (without regard to whether they should be treated atomistically or continuously), and are assumed to be disjoint, such that V = N n=1 V n . The periodicity of this system allows for writing this more explicitly in terms of projections at each unit cell. Denoting P R as the projection operator onto the subspace spanned by the DOFs in the unit cell at R, the resolution of the identity can be written as I = R P R . Further defining |P R = P R |P (and likewise |E R = P R |E ) allows for writing |P R = At this point, it is useful to expand the real-space representations of these quantities in terms of Bloch-periodic functions within each unit cell. This means expressing and its inverse, |P R = Vuc (2π) d´B Z |P k e i k·R d d k, (with similar expressions for |E ) in terms of the Bloch wavevector k, which is assumed to lie within the first Brillouin zone (BZ) which has volume (2π) d /V uc given in terms of the unit cell volume V uc . These come from the completeness rela- The susceptibility V R,R ′ , by virtue of representing a periodic system, has translational symmetry across unit cells, so Since the susceptibility V R,0 is the polarization response within the unit cell centered at R to an electric field applied to atoms in unit cell 0, namely δ R,0 I, its reciprocal space representation is given by: where the choice of unit cell 0 is arbitrary due to the discrete translational symmetry underlying this system; namely, changing the summand V R,0 to V R,R ′ changes V k to e − i k·R ′ V k , reflecting the Bloch periodicity of the response. Additionally, reciprocity V = V ⊤ in position space implies that V R,0 = (V 0,R ) ⊤ , from which it follows that (V k ) ⊤ = V −k in Bloch space. Hence, the relationship between the polarization and electric field in reciprocal space is and the transformation to reciprocal space partially diagonalizes the problem, reducing it to one that can be solved within the unit cell. All of the above relations also hold when considering the vacuum electromagnetic field Green's function G (0) relating |E = G (0) |P and solving Maxwell's equations in vacuum under the same periodicity as the susceptibility V. Thus, the Green's function in reciprocal space can be written as for k in the BZ, and |E k = G for this system. These equations become easier to manipulate in reciprocal space. In particular, (41) becomes for k in the BZ, so these can be formally solved to yield just as in (5), If the free polarization sources and incident fields arise from quantum and thermal fluctuations, they satisfy the fluctuationdissipation theorem after exploiting reciprocity to equate asym(V) = Im(V) and asym(G (0) ) = Im(G (0) ) in position space. At this point, it becomes necessary to transform the fluctuation-dissipation theorems into Bloch space. As the fluctuating fields and polarizations are correlated only with themselves and not with each other, then |E R ′ (ω ′ )| = 0 as before. We start with from using the real space fluctuation-dissipation theorem. If both sides are integrated over k ′ in the BZ, then this yields

so this finally yields the integrated reciprocal space fluctuation-dissipation theorem
which in turn yields the Bloch space fluctuation-dissipation theorem: As the same reciprocity properties of V k hold for G (0) k , then it also follows that (at thermal equilibrium), Having thus derived the fluctuation-dissipation theorems for systems with Bloch periodicity, we may now derive the vdW free energy at equilibrium temperature T . This once again requires evaluation of the quantity P(λ), E(λ)) = Tr(|E(λ) P(λ)|) . Taking λ = 1 for now (restoring explicit factors of λ later), Tr[|E(ω) P(ω ′ )|] = R Tr[|E R (ω) P R (ω ′ )|] . Using the fact that, and that, in analogy to the case of compact molecules, after using the fluctuation-dissipation theorems for periodic structures and the properties of Dirac delta functions, the rest of the derivation follows essentially identically to the case of compact polarizable bodies. The integration of the Dirac delta function leaves a sum over R of a quantity independent of R; this physically reflects the invariance of this periodic problem with respect to discrete translations, and the fact that periodic structures are infinite, albeit with the interaction free energy per unit cell remaining finite. Ultimately, the vdW interaction free energy among a collection of extended polarizable bodies, per unit cell, is given by: where the prime again implies a half weight on the l = 0 term in the sum.
One can also derive a compact formula for the radiation spectrum per unit cell between extended polarizable bodies by following the same steps as in the case of compact molecules. In particular, the radiation spectrum between polarizable bodies m and n (which may be the same) at a given ω can be written as where dependence on (ω, k) is made implicit. The formulas for the far-field emission W (m) and the heat transfer W m→n remain the same in terms of Φ (m) n , though they are now thermal emission or RHT spectra per unit cell.

B. Fast molecular scattering matrix elements evaluations via Ewald summation
The above formulas do not make explicit reference to molecular or macroscopic bodies, but just as for general polarizable bodies, the corresponding DOFs may be separated to yield formulas that yield greater physical insight and ease of implementation. In particular, the molecular basis functions are the same as in the compact case, and just as for general polarizable bodies, if the configuration of Bloch periodic macroscopic bodies is fixed, they form a scattering background with a modified Green's function in which vdW interactions and radiative energy exchange may be computed among periodic molecular structures. In practice, just as in the compact case, for most macroscopic geometries, the matrix elements f pi |G mac f qj would need to be computed using costly 6-dimensional numerical cubature; Bloch periodic boundary conditions adds another cost in the form of summing over lattice vectors too. For this reason, our code only implements computations where G mac can be expressed analytically in terms of G (0) , namely when either no macroscopic body is present (i.e. vacuum) or a single PEC plane is present (which can again be computed via image theory). This may not be such a severe practical limitation though, as materials like graphene, which have become of great recent scientific interest, can be treated atomistically in our model, so interactions among graphene sheets and molecular crystals in vacuum may be considered without significant issues. Additionally, the same caveats as for compact molecular structures apply with respect to the overlaps of the Gaussian basis functions with a PEC plane. With this in mind, we now turn to deriving the expressions for the vacuum Green's function matrix elements in the molecular basis in periodic geometries. The expression of the vacuum Green's function in terms of (32) and the facility in analytically performing the resulting spatial integrals over Gaussian basis functions ensures that the formulas we obtain are analytical and fast converging over the real and reciprocal lattice summations; the expressions bear many similarities with Ewald summation, while the nonzero Gaussian widths ensure that certain divergences are mitigated, just as for isolated (non-periodic) basis functions. We perform the following derivations at ω = i ξ, and notationally suppress the functional dependence on ω for brevity; formulas valid for real ω can be obtained by substituting ξ = − i ω at the end results.
Our use of Gaussian basis functions of relatively large widths (especially so when one considers phonons [31,32]), ensures that in periodic geometries, the field responses G (0) k can no longer be treated from the perspective of simple point dipoles. Instead, one must directly compute the matrix elements using the definitions of the basis functions |f pi ≡ |f p e i , where the widths of the Gaussian basis functions f p depend on (i ξ, k) via the susceptibility matrix α k ; notationally, |f p+R,i refers to the periodic image of |f pi at lattice vector R, and is represented in position space as f p (x − R)e i . Performing this summation over the real lattice yields slow conditional convergence, so the goal is to transform this summation into equivalent fast and absolutely convergent sums, accounting for the nontrivial Gaussian screening widths. In particular, this involves rewriting and then splitting the integral in eq. (32) over u, from [0, ∞) to the ranges [0, κ) and [κ, ∞), where κ is a user-specified Ewald splitting parameter that controls the speed of convergence [113,114]. Explicitly, this involves writing such that G  [113,114], but with appropriate changes accounting for the molecular basis functions having a finite Gaussian spread rather than corresponding to point dipoles.
The term G (0)SR kpi,qj is evaluated over the real lattice, giving expressions independent of periodic dimensionality. In particular, making the same variable substitutions v = 2(σ 2 , and θ = 2(σ 2 p + σ 2 q )ξ/c, then carrying out the integration with respect to w over the range [ν, 1) yields for any periodic lattice. The term G (0)LR kpi,qj is evaluated over the reciprocal lattice, leading to different expressions for different periodic dimensionalities. For a 1D-periodic system, the lattice vectors lie along a single direction with R = na (where a = |a|), and the reciprocal lattice vectors are likewise g = nb, where b = 2πa/a 2 . Defining r p −r q = ∆r a/a+∆r ⊥ where ∆r is the component of the displacement between the two atoms along the periodic axis a and ∆r ⊥ is the orthogonal projection, then in the integrand, |r p + R − r q | 2 = (∆r + na) 2 + ∆r 2 ⊥ . The real lattice sum is expressed as: Defining the function, allows for use of the Poisson summation formula is the Fourier transform with respect to the coordinates along the periodic axis. Using the facts that b = 2πa/a and g = nb, and that (k − 2πn/a)∆r = (k − g) · (r p − r q ) by definition, it follows that the integral over w, can be written in terms of the reciprocal lattice sum, having defined η 2 = θ 2 + 2(σ 2 p + σ 2 q )|k − g| 2 and ρ ⊥ = |∆r ⊥ | √ 2(σ 2 p +σ 2 q ) . One further variable substitution y = ν 2 w 2 and an expansion of the exponential term involving ρ ⊥ in terms of its Taylor series finally yields the long-range contribution for a 1D-periodic system along an arbitrary axis of periodicity, given by: in terms of the exponential integral functions E 1 (x) = ∞ x t −1 e −t dt and E s+1 (x) = s −1 (e −x − xE s (x)), which are closely related but not identical to incomplete gamma and hypergeometric functions [115].
For a 2D-periodic system, the lattice vectors and reciprocal lattice vectors can be defined to lie in a plane with orthonormal vectors e 1 & e 2 , with vector e ⊥ lying normal to the plane. This allows for writing R = R 1 e 1 + R 2 e 2 , g = g 1 e 1 + g 2 e 2 , and k = k 1 e 1 + k 2 e 2 , regardless of lattice geometry, and r p − r q = ∆r 1 e 1 + ∆r 2 e 2 + ∆r ⊥ e ⊥ . This means k · R = k 1 R 1 + k 2 R 2 and |r p + R − r q | 2 = (∆r 1 + R 1 ) 2 + (∆r 2 + R 2 ) 2 + (∆r ⊥ ) 2 . Once again, from the integral over w, the function f (l 1 e 1 + l 2 e 2 ) = exp(θ 2 /4 − i(k 1 l 1 + k 2 l 2 )) can be used in the Poisson summation formula where the Fourier transform, is written in terms of η as above and ρ ⊥ = . Performing the integration over w finally yields for a 2D-periodic geometry parallel to an arbitrary plane. In principle, the infinite set of real lattice vectors R and reciprocal lattice vectors g must be used for the above summations. In practice, however, these sums are fast-converging allowing for truncation after a relatively small number of vectors R and g, provided an appropriate choice of the Ewald parameter κ. The optimal value of this parameter strongly depends on the separation |r p − r q |, frequency ω, and effective Gaussian width 2(σ 2 p + σ 2 q ), and the last among those in particular depends heavily on the material properties of the body in addition to the geometry; a full convergence analysis is beyond the scope of this work.

C. vdW interactions and RHT between parallel graphene sheets
In our past work [31], we have considered only graphene in the RMB framework as it has proved to numerically work between two parallel graphene sheets at separation d at zero (blue) and room (red) temperatures. Inset: ratio of the free energy F(T ) to the zero-temperature PEC planar interaction energy EPEC.
well with the aforementioned Ewald summation procedure; by contrast, hexagonal boron nitride does not seem to produce such good convergence properties in practice, and we have not tried other extended periodic media using our atomistic description in RMB. We also note, as we have discussed in detail in our previous work [31], that the atomistic treatment of graphene in RMB ignores the electromagnetic effects of the interplay between delocalized electrons and phonons, so our use of graphene in this paper is meant merely to qualitatively illustrate salient behaviors in fluctuational EM interactions and to show the convergence and power of the RMB framework, not to provide high-precision quantitative results to compare with other theories. As we have already considered the interactions between a graphene sheet and a parallel PEC plane [31], we now consider the interactions between two parallel graphene sheets in vacuum separated by distance d.
For this system, we consider the vdW interaction free energies in Fig. 4 at zero temperature T = 0 and room temperature T = 300 K; numerical difficulties in this system preclude consideration of separations outside of the range d ∈ [1 nm, 50 nm]. The power laws ∂ ln |F (T )| ∂ ln(d) at both temperatures show significant deviations from the conventional pairwise prediction of −4 as well as the predictions in the nonretarded random phase approximation of −3 [116]. In particular, both increase from values more negative than −3 at small d to those less negative than −3 at larger d, and the room temperature power law in particular exhibits more sensitivity to d at for d < 5 nm. These behaviors are because the static Gaus- between two parallel graphene sheets at separation d at low (blue) and room (red) temperatures. Inset: ratio of the exchanged power P (T ) to the blackbody emission power PBB(T ) at the same temperature. sian widths for graphene near the center of the Brillouin zone, which is where the integrand is dominant, are around 3 nm, so the Gaussian basis functions overlap at such small d; the vdW free energy integrand is more sensitive to static effects at room temperature than at zero temperature, so the power law is more complicated at room temperature. Meanwhile, at larger d, electromagnetic retardation interplays with the more complicated material dispersion due to phonons, leading in the case of two graphene sheets in vacuum to an initial increase rather than a decrease in the power law. All of these effects are qualitatively very similar to those seen for two parallel long carbyne wires in vacuum [31], with quantitative differences arising in the behaviors of these two carbon allotropes due to greater dimensionality and the lack of finite size effects in the case of graphene. Additionally, as exemplified in the ratio of F (T ) to the corresponding zero-temperature interaction energy of two PEC planes E PEC = π 2 cA 720d 3 , the free energies themselves are smooth and monotonic functions of d, though the nonmonotonic behavior of the ratio is exemplified in the behavior of the corresponding power laws.
We also consider the RHT powers in Fig. 5 at low temperature T = 100 K and room temperature T = 300 K for one of the sheets, where the other sheet is consistently assumed to be maintained at zero temperature. The power laws ∂ ln |P (T )| ∂ ln(d) deviate significantly from the prediction of −4 by a pairwise summation of near-field RHT for two sheets, in both cases behaving nonmonotonically and remaining less negative than −3 in the range of separations of interest due to the confluence of factors involving the overlap of Gaussian basis functions particularly for d < 5 nm and the complicated interplay of geometry, material dispersion, and electromagnetic retardation at all separations. That said, the room temperature power law remains more negative than the low temperature power law due to the existence of resonances in Φ at higher frequency that depend more strongly on separation and are not exponentially suppressed as they would be at low temperature. The room temperature power law in particular behaves qualitatively very similarly to that for radiative heat transfer between two fullerenes in vacuum [32], with quantitative differences again arising in the behaviors of these two carbon allotropes due to greater dimensionality and the lack of finite size effects in the case of graphene; for the system of two graphene sheets, numerical difficulties again preclude consideration of separations outside of the range d ∈ [1 nm, 50 nm], but we expect based on the very similar results of fullerene that as d drops below 1 nm, the power law for decreasing d would continue to drop toward a local minimum and then sharply increase and essentially saturate near zero, corresponding to a saturation rather than a divergence of the RHT power itself with decreasing d due to the strongly nonlocal material response of graphene as captured in the atomistic model used in the RMB framework. Furthermore, the RHT powers themselves monotonically decay with increasing separation and are significantly larger than the corresponding blackbody emission powers at each temperature, though the normalized power is larger at low temperature largely because the corresponding blackbody emission power is so much less there than at room temperature.

VI. SIE FORMULATION OF INTERACTIONS AMONG MOLECULES AND MACROSCOPIC BODIES
We turn to the SIE formulation of Maxwell's equations for the special case of macroscopic bodies defined by sharp boundaries between regions where the permittivity is local and homogeneous. This allows for writing scattering quantities involving macroscopic bodies in terms of surface DOFs and the homogeneous Green's functions on each side of a surface, rather than volumetric DOFs and associated susceptibilities. In brief, rather than solving a discretization of Maxwell's equations in differential or integral form in the full volume of a body, we instead assign fictitious electric and magnetic currents to boundaries between permittivity regions and solve for them by enforcing continuity of the tangential electric and magnetic fields across each boundary, so that the fields radiated by the fictitious currents are the scattered fields accounting for multiple scattering within and between bodies. Note that magnetic surface currents are needed even for bodies with vanishing magnetic susceptibility, as the fictitious surface currents are simply the tangential components of the total fields. The use of surface DOFs already provides a drastic reduction in computational complexity over methods that use volumetric DOFs; while the macroscopic surface basis functions {|b β } may be spectral or other arbitrary functions, particular computational gains can be realized via localized basis functions, such as Rao-Wilton-Glisson (RWG) basis functions, where the ability to heterogeneously mesh a surface allows for treatment of general macroscopic surface shapes with arbitrary features. Below, we define the SIE operators and provide formulas for the vdW interactions and RHT among molecules and macroscopic objects in this framework; we do not demonstrate any particular computational implementation of these formulas, leaving that for future work.
Application of the general formulas using the SIE method requires appropriate modifications and operator substitutions. In particular, given a collection of macroscopic objects labeled by the index n, the DOFs are defined on their corresponding surfaces, with interactions mediated by the exterior vacuum Green's functions G (0) (by our assumption, though the exterior medium could in principle be a different nontrivial permittivity) and within the macroscopic body interiors G (0,n) ; having assumed that the macroscopic bodies are made of homogeneous, local, isotropic susceptibilities, we clarify that G (0,n) is the homogeneous Maxwell Green's function corresponding to the bulk material constituting macroscopic body n, as if its boundaries didn't exist. We further assume for the purposes of these derivations that the macroscopic bodies have distinct surfaces and are not embedded in each other, though the SIE formulation is general enough to allow for relaxation of those assumptions [80]. General scattering problems then obtained via a SIE scattering operator whose inverse given by −W −1 mac = G (0) + n S n G (0,n) S n [80,81], such that the scattering Green's function outside of the collection of macroscopic bodies is G (0) W mac G (0) ; here, S n is a projection operator onto the surface (rather than volumetric) DOFs of macroscopic body n. For the purposes of vdW interactions as well as thermal emission or heat transfer only among molecules, the macroscopic bodies only affect the EM field scattering properties, so the replacements T mac → W mac (and analogously T mac∞ → W mac∞ for vdW interactions) are sufficient when evaluating (10) and (21) in conjunction with (18).
In situations where one seeks to compute energy exchange between a collection of molecules and a macroscopic body, as may be useful for localized heating of a molecule by a AFM tip [64], it is incumbent to perform additional simplifications (beyond the substitution T mac → W mac ). This is because the macroscopic DOFs are only defined at their surfaces, without any reference to volumetric degrees of freedom, so the SIE formulation leads more naturally to a definition of heat transfer in terms of the Poynting flux through the surface of a given macroscopic body n due to fluctuating volumetric polarization sources in molecule m. It is useful to start with the result of first performing the aforementioned substitution along with P n → S n into (22): This expression can be further rewritten to obtain a formula that is conceptually and technically similar to previously derived formulas for heat transfer between macroscopic bodies [80]. First, the combination W −1 mac + G (0) is block diagonal in the space of macroscopic bodies, such that S n (W −1 mac + G (0) ) = −S n G (0,n) S n . Hence, the above expression can be rewritten as: Next, if the blockwise inversion to evaluate T from (18) is performed accounting for the identity Finally, we may define a modified SIE operator − W ′−1 mac = G (0) + G (0) T mol G (0) + n S n G (0,n) S n (61) in analogy to (34), as an effective SIE operator where the exterior medium is no longer vacuum but encodes the scattering properties (to infinite order) of the molecules as a background medium. This allows for writing the radiative energy exchange between a molecule m and a macroscopic body n can be written as: S n asym(G (0,n) )S n × W ′ mac G (0) T mol P m in a more compact way. Conceptually, this formula describes the energy transfer as a Poynting flux from volumetric sources in molecule m, whose correlations are proportional to asym(V −1 † m ), through the surface of macroscopic body n via S n G (0,n) S n , where scattering between all of the molecules and macroscopic bodies is accounted to all orders via the combination of W ′ mac and T mol . It is exactly this substitution in (61) that further allows for computing the heat transfer among macroscopic bodies in the presence of molecules in the SIE framework. In such a case, the energy flow from the surface of one macroscopic body (due to fluctuations in its interior) through the surface of another is desired, with the molecules simply modifying the scattering properties of the medium exterior to the macroscopic objects. This makes the heat transfer between macroscopic bodies m & n Φ (m) n = 4 Tr asym(G (0,m) )S m W ′ † mac S n × asym(G (0,n) )S n W ′ mac S m the same as that in [80], using (61) in the presence of the molecular bodies. In all of these formulas, the molecular DOFs are expressed in terms of the Gaussian basis functions |f pi as usual, while the macroscopic DOFs are expressed in terms of basis functions denoted |b β : the latter may in principle be either spectral or localized basis functions, but localized RWG basis functions are preferred for convergence in arbitrary macroscopic geometries that do not have a high degree of translational or rotational symmetry. This means that the expression of T mol in terms of |f pi , along with the matrix elements f pi |G (0) f qj , b β |G (0) b β ′ , and b β |G (0,n) b β ′ (for a macroscopic body labeled n of a given homogeneous susceptibility) are needed: routines to compute these matrix elements have already been implemented, the former two in our new code and the latter two in the SCUFF-EM boundary element solver [82]. However, on top of this, the matrix elements f pi |G (0) b β need to be computed as well: this has not yet been implemented, but may be done through appropriate conjunction of the SCUFF-EM code with our code as both are open source software. We do note that just as for molecular bodies above a PEC plane, these derivations assume that the molecular basis functions can be associated purely with the space external to the macroscopic bodies, which is not exactly true given that Gaussian basis functions do not have compact support, and this assumption becomes somewhat more questionable when the center of a basis function is less than one Gaussian width away from the boundary of a macroscopic body. That said, this approach should still qualitatively capture the effects of screening on interactions between molecules and macroscopic bodies even at such short separations, and is an improvement on our previous approximation of molecular basis functions as point dipoles in their interactions with macroscopic bodies [30] (which was only justifiable in the absence of phonons so that the Gaussian widths were much smaller than the considered separations between molecules and macroscopic bodies).

VII. CONCLUDING REMARKS
The RMB formulation of fluctuational electrodynamics makes clear that "molecular" and "macroscopic" bodies can be treated on the same footing, given appropriate atomistic or continuum descriptions of each. It allows atomistic descriptions of material bodies based on coupled effective electronic and nuclear oscillators, accounting for short-range electronic correlations and phonons whose properties are obtained from ab-initio density functional theory calculations, and is in principle compatible with arbitrary continuum descriptions of material response as well. It can in principle be extended to account for material bodies treated with continuum response theories when such bodies have arbitrary shapes (beyond simple planar structures). Furthermore, the power of this formulation lies in the analytical formulas for the electromagnetic interaction matrix elements of material bodies treated atomistically using the aforementioned oscillator model, sidestepping questions of convergence common to finite-volume or discrete-dipole computational techniques by assigning Gaussian basis functions created from material response properties obtained within the RMB framework itself.
There are several shortcomings and open questions that require further attention. Chief among them is that the atomistic oscillator model is physically accurate only for insulating or weakly conducting system, and is less appropriate for strongly metallic or semimetallic systems where electron delocalization effects are more visible in conjunction with phonons and long-range electromagnetic interactions. This has been discussed in detail in our prior work [31], particularly concerning how the RMB framework can capture the salient geometric and phononic properties of graphene and related atomically-thin materials like hexagonal BN, which will be similar, but cannot capture the inherent electron delocalization in graphene that is absent from polar dielectric media like hexagonal BN. Related to this, in extended media where the effects of phonon and electron delocalization would be most relevant if present, the Ewald summation procedure applied to Gaussian basis functions constructed from the susceptibility within each unit cell is not guaranteed in practice to yield numerically well-behaved results: for example, proper convergence is obtained for infinite sheets of graphene, but not for infinite sheets of hexagonal BN.
Even for compact molecules, the widths of the Gaussian basis functions encode information about the anisotropy of the molecule as a whole but are not themselves anisotropic for each atomic basis function. It remains to be seen for a broader variety of molecules interacting at the mesoscale the extent to which such a change in the basis functions may make a difference, but that is beyond the scope of this work. Furthermore, for compact molecules and extended atom-scale structures, DFT calculations may yield effective internuclear coupling matrices K I that go far beyond nearest neighbors, but numerical convergence of such long-range couplings is not always guaranteed in practice; therefore, some care must be Quantum chemistry Continuum FED Relationship Tij(ω, x, x ′ ): χRPA(ω, x, x ′ ) = χRPA = χ (0) + χ (0) vχRPA T = V + VG (0) T i,j ∂i∂jT nonret ij (ω, x, x ′ ) TABLE I. Glossary: note that T nonret is computed at each ω such that V is evaluated at that ω but G (0) is evaluated at zero frequency.
taken in plugging those matrices into code built on the RMB framework, and it may be necessary to restrict couplings to nearest or next-nearest neighbors.
The extension to include continuum bodies of arbitrary geometries has yet to be computationally implemented in practice. That is beyond the scope of this work, but we imagine this to be relatively straightforward. However, the extent to which the possible overlap of relatively large atomistic Gaussian basis functions with hard boundaries of continuum bodies may degrade the accuracy of predictions in the RMB framework remains to be tested through direct comparisons with relevant past [53][54][55] and future experiments; these would be the ultimate tests of the reliability and raison d'être of our RMB framework.
Finally, we note that while the measurements of nearfield RHT between metallic tips and substrates at nanometric gaps [20,[27][28][29] can be modeled using the RMB framework if electrons and phonons are localized within each body to preclude the possibility of conduction between bodies, the current RMB framework is unable to model total heat transfer when both radiative and conductive processes between bodies are relevant, particularly in the extreme near-field. Extending the RMB framework and associated code to handle such situations is the subject of an upcoming manuscript.
Further testing in diverse combinations of molecular and macroscopic bodies will doubtless yield more questions, so addressing all of these issues will be the subject of many future works. We anticipate that other researchers may be able to make use of our code both to model mesoscale FED phenomena in specific systems and to answer some of these broader outstanding questions.