Universal Theory of Light Scattering of Randomly Oriented Particles: A Fluctuational-Electrodynamics Approach for Light Transport Modeling in Disordered Nanostructures.

Disordered nanostructures are commonly encountered in many nanophotonic systems, from colloid dispersions for sensing to heterostructured photocatalysts. Randomness, however, imposes severe challenges for nanophotonics modeling, often constrained by the irregular geometry of the scatterers involved or the stochastic nature of the problem itself. In this Article, we resolve this conundrum by presenting a universal theory of averaged light scattering of randomly oriented objects. Specifically, we derive expansion-basis-independent formulas of the orientation-and-polarization-averaged absorption cross section, scattering cross section, and asymmetry parameter, for single or a collection of objects of arbitrary shape. These three parameters can be directly integrated into traditional unpolarized radiative energy transfer modeling, enabling a practical tool to predict multiple scattering and light transport in disordered nanostructured materials. Notably, the formulas of average light scattering can be derived under the principles of fluctuational electrodynamics, allowing the analogous mathematical treatment to the methods used in thermal radiation, nonequilibrium electromagnetic forces, and other associated phenomena. The proposed modeling framework is validated against optical measurements of polymer composite films with metal-oxide microcrystals. Our work may contribute to a better understanding of light-matter interactions in disordered systems, such as plasmonics for sensing and photothermal therapy, photocatalysts for water splitting and CO2 dissociation, photonic glasses for artificial structural colors, and diffuse reflectors for radiative cooling, to name just a few.

Predicting the complex optical phenomena manifesting in disordered nanomaterials represents a major challenge in the field of computational modeling. Plasmonic nanoparticle dispersions for sensing and photothermal therapy,[1,2] heterostructured photocatalysts for water splitting and CO2 dissociation,[3] diffuse reflectors for radiative cooling,[4] and porous membranes for solar water desalination[5] are a few examples where the complex inter-relation between near-field coupling and multiple scattering with a variation in particle morphology, orientation, and size impose severe limitations to theoretically predict the system optical response. Conventional modeling based on computational electromagnetics, such as the Finite Difference or Finite Elements methods, are often unsuitable to quantify the macroscopic optical properties of random media, namely, their specular and total transmittance/reflectance or the intensity distribution, all critical to assess the performance of these systems. Stochastic methods appear as the most appropriate alternative to calculate these properties, yet, with few exceptions, their applicability is limited to composite media containing subwavelength structures (effective media approximations)[6] or spherical particles (radiative transfer simulations).[7] As a result of the limitations in modeling, the majority of designs in disordered nanophotonic materials are driven by a phenomenological approach, whereby multiple samples are fabricated and tested in an iterative process that is time-consuming and expensive.

Light transport in random media is commonly addressed through the radiative transfer theory,[8,9] which describes the propagation of the light specific intensity through a composite medium containing a random distribution of independently scattering particles. It is notable that, in principle, the theory is applicable to arbitrary particle geometries and groups of particles.[9] In practice, however, radiative transfer modeling is most commonly applied to study light transport in composites with spherical particles.[7,10] This is due to the scattering properties of the spherical particles being independent of the direction and polarization of the incident light.[11] In this particular scenario, the solution of the radiative transfer equation (RTE) for unpolarized light requires three parameters: the particle’s absorption cross section, Cabs, the scattering cross section, Csca, and the asymmetry parameter, μsca.[7,12] The latter is an indicator of the scattering anisotropy and a key element to calculate the angular distribution of scattered fields.[12]

The scattering of nonspherical particles, on the other hand, varies with the incident angle and polarization, and the RTE usually becomes too complex to solve.[9] Under the independent scattering approximation, however, the correlation of the scattered field from different particles vanishes,[9] and the scattering properties of an ensemble of randomly oriented particles can be approximated by the orientation averaged from a single particle (Figure a). For unpolarized light, the RTE becomes scalar,[14] and the orientation and polarization averaged Cabs, Csca, and μsca parameter triad can be used for radiative transfer simulations of arbitrary particles, following the same methodology of spherical particles. The same principle could be applied to more complex scenarios, for instance, a medium containing denser particle distributions or even particle clusters. In this case, the averaging should now be performed over a properly chosen collection of particles for which the effects from short-range correlations (due to collective interaction and interference of scattered fields) are prevalent (Figure b).[13] Similarly, heterogeneous systems with particles of different sizes and optical properties can also be studied.

Figure 1

Average scattering of randomly oriented particles. (a) A collection of independent scattering particles (in this case prolate spheroids) randomly oriented in space, is approximated by the averaged sum over all possible particle’s orientations. In the schematic at the bottom, the index m represents a particular orientation of the particle, and M is the total number of particle orientations. In both schematics, the intensity of the incident light beam (purple arrow) decays due to scattering of the particle (yellow arrows). The red dotted circle represents a characteristic domain where short-range correlations are relevant. (b) In more dense particle systems, the independent scattering approximation applies to a collection of particles within a properly chosen domain that considers the effects of short-range correlations.[13] Average scattering, thus, is calculated over a characteristic collection of particles.

The computation of orientation and polarization average scattering (average scattering, from now onward) for arbitrary nanoparticles is, however, nontrivial. Standard brute-force methods based on averaging over many plane-wave simulations at different angles of incidence and polarizations can be computationally expensive.[15] Alternatively, semianalytical solutions relying on spherical wave expansion have demonstrated considerable improvements in the efficiency of the calculations.[15−19] For example, formulas for direct computation of average scattering have been developed for axially symmetric objects, such as cylinders,[16] spheroids,[16] and clusters of spherical particles.[17] However, the restrictive use of the spherical wave basis in this approach still imposes some constraints. Such is the case for objects with no axial symmetry or with sharp edges where the expansion of the scattered fields into spherical waves is not trivial.[16] Often, in these problems, the scattered fields are more conveniently expanded through the basis relying on surface or volume discretization, such as surface currents in the Boundary Elements Method[20] (BEM) or discrete dipoles in the DDA.[21] For average scattering calculations, nonetheless, the expanded fields have to be transformed into spherical waves, and the efficiency of the method is appreciably reduced.[21]

In this Article, we present a universal theory of average light scattering from randomly oriented scatterers (single or collections of objects) of arbitrary shape and demonstrate a practical methodology for radiative transfer simulations in disordered nanostructures. The results section of the paperis organized into six sections: (i) In the first section, we derive the formulas for polarization-and-orientation-averaged Cabs, Csca, and μsca for arbitrary-shaped scatterers. The formulas are independent of the wave basis and, therefore, can be implemented by integral-equation methods for electromagnetic scattering, such as T-Matrix Method,[9] BEM,[20] or DDA.[21] Based on these results, we evaluate the accuracy limits of other expressions for average light scattering commonly used in the literature.[22,23] (ii) In the next section, we demonstrate that the formulas of average light scattering can be derived through the principles of fluctuational electrodynamics and, hence, can be computed through the mathematical methods used in studies of near-field thermal radiation,[24−27] Casimir forces,[27] and vacuum friction.[28,29] In this context, we develop a computational application to numerically compute averaged light scattering,[30] which is based on the fluctuating-surface-current BEM.[31,32] (iii) In the following section, we validate the theory and simulation code for average scattering simulations against other analytical solutions.[33] (iv) In the forth section, we analyze the advantages of the theory of average light scattering against brute-force averaging based on plane-wave simulations at different angles of incidence. (v) Next, we discuss how the three average light scattering parameters can be applied for modeling of radiative transfer in disordered nanostructures. (vi) The accuracy of the modeling framework is demonstrated in the final section, showing excellent agreement with optical measurements of polyethylene (PE) film composites with monoclinic vanadium dioxide [VO2(M)] microcrystals.

Results

Theory of Average Light Scattering of Randomly Oriented Particles

As discussed previously, light scattering from a collection of independent scatterers of arbitrary morphology and randomly oriented in space (Figure ) is equivalent to the average light scattering over all orientations and light polarizations.[14] We particularly focus on the average absorption cross section, ⟨Cabs⟩, scattering cross section, ⟨Csca⟩, and asymmetry parameter, ⟨μsca⟩, where , the P index runs over the two orthogonal polarizations, and k̂i is the direction of the incident field. The asymmetry parameter, μsca, defines the degree of anisotropy of scattering relative to k̂i:[11]where k̂s is the direction of the scattered field and psca(k̂s, k̂i) is the scattering phase function.[14] By definition of psca, . Thus, μsca > 0 (μsca < 0) represents cases of forward (backward) anisotropic scattering and μsca = 0 represents isotropic scattering.

Our derivations are based on the Lippmann–Schwinger approach, a general formalism for electromagnetic scattering phenomena.[27] In this approach, the scattered fields are given by [in this notation, ], where Ei is the incident field, is the free space Dyadic Green function, and is the scattering operator (Supporting Information, eqs S4 and S5, respectively). The mathematical form of is dictated by the expansion basis and the geometry and optical properties of the scatterer. For example, using the spherical wave basis, the operator for a spherical particle is , where freg is the spherical waves regular at the origin, Tll is the Mie scattering coefficient, and † is the conjugate transpose operator.[27] In integral-equation methods for electromagnetic scattering, such as the T-Matrix Method,[9] BEM,[20] or DDA,[21] the form of and has to be computed prior to any scattering calculation.

In this context, the formulas of Cabs and Csca for an incident plane wave of amplitude E0 are given by (details in Supporting Information, Section 1.1):[26,34]where k0 is the wavevector in free space, ⊗ is the tensor product, and is the particle’s induced potential (Supporting Information, eq S2), with . The operators and , where is the adjoin of , represent the Hermitian and anti-Hermitian part of , respectively. The trace is defined as .

The formula of μsca is, to the best of our knowledge, presented here for the first time:where the index j in k̂i and in the partial derivative ∂, represents the global coordinates of the system (e.g., j = x, y, z). The formula is derived from the Lorentz force of the scattered fields over the induced currents in an object (Supporting Information, eq S10).

Because the operators and are independent of the direction of the incident field, ⟨Cabs⟩ and ⟨Csca⟩ are uniquely determined by ⟨Ei ⊗ Ei⟩, also known as the free space self-correlator.[27] For Ei in the form of monocromatic plane waves, this term is given by the following (see derivation in Supporting Information, eq S8):

Similarly, ⟨μsca⟩ requires the following (see derivation in Supporting Information, eq S11):

Using eqs and 3, we derive the following expressions:which represent a universal recipe for average light scattering, compatible with integral-equation methods for electromagnetic scattering. As an illustrative example, in the Supporting Information, we derive the respective formulas for BEM and T-Matrix using these expressions (Supporting Information, sections 1.2 and 1.3, respectively).

The relations, eqs , 4b, and 4c, can be easily generalized for clustered particles and heterogeneous composites containing different types of particles (Supporting Information, section 1.4). The light scattering parameters for an individual particle particle n in the cluster, that is, ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩, are also obtained directly from these relations (details in Supporting Information, section 1.5).

We finalize this section by discussing a common approximation for ⟨Cabs⟩ and ⟨Csca⟩ found in the literature:[22,23]where Cabs, and Csca, (j = x, y, z) are, respectively, the absorption and scattering cross sections for a plane wave polarized in the j-direction. As demonstrated in the Supporting Information (section 1.6), the expression ⟨Cabs⟩ ≈ 1/3(Cabs, + Cabs, + Cabs,) corresponds to a particular case of eq for subwavelength objects, while the analogous expression for ⟨Csca⟩ holds only if the polarizability of the object is a diagonal tensor.

Average Light Scattering Derived from Fluctuational Electrodynamics

The trace formulas for ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ share many similarities with the relations found in studies of fluctuational electrodynamics, namely, thermal radiation and nonequilibrium electromagnetic forces.[27,31] For example, in the framework of fluctuational electrodynamics, the thermal radiation absorbed by an isolated object, Pabsth, is[27]where and , where T is the temperature of the environment, ω is the angular frequency, kB is the Boltzmann constant, and ℏ is the reduced Planck constant. On the other hand, from light scattering theory,[11]Pabsth = ∫0∞ dω4π⟨Cabs⟩Bω(T), where is the Planck distribution. This leads to the following relation:This formula can be easily adapted to compute ⟨Csca⟩ by replacing for .

Another relation comes from the electromagnetic friction induced by a moving photon gas on a stationary object.[28,29] To a first order approximation, the friction coefficient γf, can be expressed in terms of ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ as follows (Supporting Information, section 2.1):where ⟨Cpr⟩ = ⟨Cext⟩ – ⟨μscaCsca⟩ and ⟨Cext⟩ = ⟨Cabs⟩ + ⟨Csca⟩ are the average radiation pressure and extinction cross sections, respectively.[11]

As illustrated by eqs and 6, the formulas of average scattering can be obtained through the principles of fluctuational electrodynamics. Consequently, the vast library of analytical solutions[27,28] and numerical algorithms[31,35] developed in the context of nonequilibrium energy and momentum transfer can be used to compute ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩. For example, the thermal DDA[35] and fluctuating current BEM[31] for thermal radiation simulations, have explicit relations for Φabs that can be adapted to compute ⟨Cabs⟩ and ⟨Csca⟩. On the other hand, the fluctuating current BEM also includes routines to compute ,[32] which can be adapted for ⟨μsca⟩.

In this context, we developed a computational code for average light scattering simulations, based on the fluctuating surface-current BEM formulation for nonequilibrium energy and momentum transfer.[31,32] The code is implemented as an application of the SCUFF-EM software[36] and can be acceded here.[30] Similar to other simulation tools based on BEM,[31,32,37,38] our code supports objects of arbitrary morphology and groups of objects (Supporting Information, Figure S1), offering a convenient platform to explore the full potential of the average scattering theory presented here.

Validation of Average Light Scattering Theory against Analytical Solutions

To validate the average scattering theory and BEM simulation code, we consider the problem of light scattering by a randomly oriented sphere dimer (Figure ), which has a known solution under the T-matrix approach.[33] The dimer consists of two silver spheres of diameter, D = 200 nm, separated by a gap of (i) Δx = 2 nm and (ii) Δx = 200 nm. The light scattering parameters of a single sphere obtained from Mie Scattering Theory[11] are also plotted as a reference.

Figure 2

Average light scattering of randomly oriented silver sphere dimer. ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ of the dimer as a function of the size parameter 2Dk0 (k0 = 2π/λ). The curves ⟨Cabs⟩ and ⟨Csca⟩ are normalized to the number of spheres, Nsp = 2, and spheres volume, Vsp, for direct comparison with the absorption and scattering cross section of a single sphere. The results are computed by the BEM code for average light scattering simulations[30] and compared against the analytical solution.[33] The scattering parameters of a single sphere, that is, absorption and scattering cross section normalized to the sphere’s volume (gray areas) and asymmetry parameter (gray curve) are computed from Mie-scattering[11] and plotted as a reference. The dielectric constant of silver can be found elsewhere.[39]

The results from the average scattering theory show excellent agreement with the analytical solution by Mishchenko et al.[33] (Figure ). When Δx = 2 nm, the effects of electromagnetic coupling dominate and the average scattering curves of the dimer largely deviate from the response of a single sphere. When Δx = 200 nm, the coupling effects weaken and both ⟨Cabs⟩ and ⟨Csca⟩ approach the response of a single sphere. However, this is not the case for ⟨μsca⟩. Similar to the optical phenomenon observed in dilute particle media,[40] the scattered fields interfere constructively in the forward direction, which explains the discrepancy between the ⟨μsca⟩ curves.

As a second test, we consider the problems of average scattering from randomly oriented oblate and prolate spheroids, which has a known solution under the T-matrix approach.[33] The results are displayed in the Supporting Information (Figure S3), showing excellent agreement up to Lmax/λ ≳ 2.5, where Lmax is the longest ellipsoid axis. At shorter wavelengths, there is a discrepancy of ∼5% associated with the size of the mesh used in BEM simulations. Even with this discrepancy, the results are consistent with the predictions from the analytical solution, allowing to validate the theory presented here.

Comparison against Brute-Force Averaging Methods

Consider, for example, the computation of ⟨Csca⟩ by averaging over many plane-wave simulations at different angles of incidence (brute-force averaging). In terms of the Lippmann–Schwinger approach to scattering and using eq , this technique can be expressed aswhere m indicates a plane-wave simulation at a particular angle of incidence and polarization and M is the total number simulations.

Similar to other formulas derived in this work, eq is a general recipe that illustrates how to compute ⟨Csca⟩ using brute-force averaging. From this formula, we can deduct the steps required by integral-equation methods of electromagnetic scattering, that is:

Compute the expansion of Ei using the particular basis.

Repeat “step 1” M times.

Determine the mathematical form of and for a particular expansion basis.

Replace Ei, , and into eq .

On the other hand, as illustrated by eq , the theory of average light scattering does not require Ei. Consequently, the computation of ⟨Csca⟩ is reduced to only two steps:Because this approach requires less steps than brute-force averaging, it could potentially enable more efficient computation of ⟨Csca⟩. The same arguments applies for ⟨Cabs⟩ and ⟨μsca⟩. Note, however, that the computational cost of each technique is relative to the particular algorithm implemented. Thus, specially for problems where M is small, brute-force averaging could be performed with higher computational speed and less memory requirements.

Determine the mathematical form of and , for a particular expansion basis.

Replace and into eq .

Radiative Transfer Modeling for Random Media with Scatterers of Arbitrary Morphology

For unpolarized light and under the independent scattering approximation, the steady-state RTE for randomly oriented scatterers in a nonabsorbing host is[9]where Iλ is the specific radiative intensity (defined as the energy flux per unit solid angle), f is the volume fraction, V is the effective volume of the scatterers, k̂·∇Iλ(r, k̂) is the rate of change of Iλ at the position r and direction k̂; and ⟨psca(cos θ)⟩ is the orientation and polarization averaged scattering phase function, where cos θ = k̂·k̂′.

Commonly, solutions of eq consider approximated expressions for the phase function in terms of μsca.[12] For example, the Henyey–Greenstein model,[12]is widely used in simulations methods, such as Monte Carlo,[7] adding-doubling,[41] and discrete ordinate.[12]

As evidenced by eqs and 9, the RTE and the average light scattering parameters ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ constitute a complete set to model radiative transfer in composites with scatterers of arbitrary morphology. As demonstrated in the next section, this modeling framework enables the quantitative prediction of the macroscopic radiative properties of a composite, such as the total transmittance (Ttot), specular transmittance (Tspec), and total reflectance (Rtot).

Validation of Radiative Transport Simulation against Experiments

We demonstrate the accuracy of the previously discussed modeling framework by comparing the radiative transfer simulations against optical measurements of a composite based on VO2(M) microcrystals embedded in a polyethylene (PE) matrix (Figure a). We considered VO2(M) microcrystals given its well-defined and highly anisotropic morphology (Figure b), providing an ideal scenario to validate the theory of average scattering and radiative transfer modeling. Additionally, the refractive index[42,43] and size of microcrystals, ensures a significant contribution from both absorption and scattering in the mid infrared (IR) spectrum.[11]

Figure 3

Characterization and average scattering simulations of monoclinic vanadium dioxide [VO2(M)] microcrystals embedded into a polyethylene (PE) matrix. (a) Photograph of VO2(M)/PE composite film. (b) SEM of as-grown VO2(M) crystals, which is mainly composed of VO2(M) bars. (c) Size distribution of VO2 bars. (d) ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ of VO2(M) bars of fixed length, L = 15 μm, and variable width, W = 0.75–5.25 μm in steps of 0.5 μm. The refractive index of the host, nhost = 1.5. The average scattering parameters showed similar dependence to W for other values of L (not shown here). The legend is given by the color bar at the top of the curves. (e) An example of one of the VO2(M) crystals with flake morphology is found in the characteristic sample. (f) Computational representation of the VO2(M) flakes, which was considered for the average scattering simulations. (g) Simulated average light scattering of the VO2(M) flake. For all average scattering simulations, the refractive index of the host, nhost = 1.5, and the refractive index of the VO2(M) bars was obtained from the literature (see “film2” in Wan et al.[42]).

First, we derived the average light scattering parameters of the VO2(M) microcrystals ensemble using a characteristic sample (Supporting Information, Figure S4c). The size distribution of the bars length (L) and width (W) is shown in Figure c, which assumes bars of square cross section. We calculated ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ of single VO2(M) bars for λ = 3–8 μm, considering the range of W and L dictated by the size distribution. The spectrum λ > 8 μm is excluded in the simulations due to the large uncertainty in the refractive index of VO2(M), which is strongly conditioned by crystal orientation, growth method, strain, and partial oxidation.[42] As shown in Figure d, ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ are strongly sensitive to W. On the other hand, the three parameters are less sensitive to changes in L, with negligible variations for L > 15 μm (Supporting Information, Figure S5). In addition to the VO2(M) bars, we noted small traces of VO2(M) flakes in the sample, such as the one shown in Figure e. These VO2(M) flakes are represented by the computational model in Figure f, with the simulated average scattering parameters shown in Figure g.

The parameters ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ of the VO2(M) microcrystals ensemble (Figure a), were estimated using the average scattering simulations of individual bars and the flake, together with the size distribution. We repeated this procedure for five different refractive indexes reported in the literature[42,43] in order to consider the variations in the optical properties of VO2(M). Using the parameters ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ of the VO2(M) microcrystals ensemble, together with Monte Carlo simulations of radiative transfer (see details in Methods), we estimate Ttot, Tspec, and Rtot of a VO2(M)/PE composite film (Figure b). The results are shown by filled areas, representing the upper and lower limits associated with the variations of the refractive index of VO2(M) and thickness of the film. The optical measurements show excellent agreement with the range predicted by simulations, which is further confirmed by comparing the spectral mean of Ttot, Tspec, and Rtot (table in Figure b). The accuracy of the simulation is further confirmed through a second test, which considers a composite film with a double concentration of VO2(M) microcrystals (Supporting Information, Figure S7).

Figure 4

Radiative transfer modeling of VO2(M)/PE film composite. (a) ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ of VO2(M) microcrystal ensemble calculated for five different refractive indexes of VO2(M), as reported by Wan et al. 2019[42] (labeled as “film1”, “film2”, “film3”, and “film4”), and Ramirez-Rincon et al. 2018.[43] The gray areas mark the upper and lower limit due to variations in the refractive index. For a given refractive index, the curves ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ were obtained as indicated by the schematic (left figure, inset), that is, average scattering simulations of bars weighted by the size distribution + average scattering of five flakes. Further details in the Supporting Information, section 4.3. The curves ⟨Cabs⟩ and ⟨Csca⟩ are normalized to the volume of the ensemble V (Supporting Information, eq S31). (b) Validation of radiative transfer theory, showing Ttot, Tspec, and Rtot of a VO2(M)/PE composite film, as obtained from optical measurements (solid lines) and simulations (filled areas). The optical properties of the PE film were extracted from optical measurements on a clear film (see Supporting Information, section 4.4). For the simulations, the absorption of PE is considered through the extinction coefficient κhost (see Methods). The composite is based on a 98 ± 4 μm thick PE film with 0.275% v/v of VO2(M) microcrystals. The upper (lower) limit in the filled areas are the results of variations in the refractive index of VO2(M) and thickness of the film (98 – 4 or 98 + 4 μm).

Conclusion

We presented a universal theory to predict the average light scattering from randomly oriented objects with arbitrary shape. The formulas of ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ can be implemented by any integral-equation method for electromagnetic scattering. Moreover, because these relations are exclusively defined in terms of the operators and , they could lead to a more efficient computation of average scattering than brute-force techniques. The general form of the average scattering formulas also provides a convenient landscape to explore the fundamental limits of scattering in random systems. For example, in analogy to the studies of scattering and absorption bounds,[26,34] the limits of forward or backward scattering of randomly oriented particles can be explored through the asymmetry parameter formula (eq ).

The demonstrated connection between average light scattering and fluctuational electrodynamics enables to extend the theory to other parameters of interest. For example, a formula for the average scattering of moving objects can be extracted from the theory of electromagnetic friction in objects at relative motion.[28] Alternatively, other expressions can be extracted directly through the self-correlators in eqs and 3, in a similar fashion than the relations obtained from the fluctuation–dissipation theorem.[25−27]

The parameters ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ are also practical for radiative transfer simulations for unpolarized light, enabling an accurate prediction of the optical properties of composites with scatterers of arbitrary shape, as demonstrated in the study of VO2(M)/PE composite films. The radiative transfer formula for randomly oriented particles (eq ) can be extended to consider other effects present in the light transport process. For example, the emitted thermal radiation from scatterers can be represented through the term at the right-hand side of the equation.[12] Similarly, the absorption of the host can be included through the term −2k0κhostIλ(r, k̂) at the right-hand side of eq .

The methodology used in the study of VO2(M)/PE composite films can be also applied to other composite media, with either dielectrics[7] or metal scatterers,[10] providing that the distance between particles is large enough to ignore the effects of short-range correlations. For more complex problems, such as clustered particles or more dense particle distributions,[44] the methodology can be extended using the formulation for multiple objects (Supporting Information, section 1.4). In this case, ⟨Cabs⟩, ⟨Csca⟩, and ⟨μsca⟩ must be obtained from simulations over a properly chosen collection of particles that better represent the effects from short-range correlations. The method, thus, could provide key insights to many problems in disordered nanophotonics, such as the effects of agglomeration into the optical absorption of gold nanostars or the impact of multiple scattering in the light trapping of heterostructured photocatalysts, as we will discuss in future works.

Methods

Fabrication and Characterization of VO2(M)/PE Composite Films

The composite was fabricated by dry mixing of VO2(M) microcrystals with low-density PE (LDPE; 42607, Sigma-Aldrich) and ultrahigh-molecular-weight PE (UHMWPE; 429015, Sigma-Aldrich) at a weight ratio of VO2(M)/LDPE/UHMWPE = 1:40:40. The mixture was then melt-pressed into a film at 200 °C. The VO2(M) crystals were produced by hydrothermal synthesis using our previously developed procedure.[45] The phase of the crystals was confirmed by X-ray diffraction and Raman spectroscopy (Supporting Information, Figure S4a and b, respectively). The volume fraction of the VO2(M) microcrystals was estimated from the weight ratio and the densities of VO2(M) (4.230 g/cm3),[46] LDPE (0.925 g/cm3), and UHMWPE (0.940 g/cm3). A micrometer was used to characterize the thickness of the film. The reported thickness corresponds to five measurements on different sections of the sample.

Optical Measurements

The total and specular transmittance and total reflectance of the VO2(M)/polyethylene composite film were measured with a Fourier transform infrared spectrometer (IRTracer-100, Shimadzu) and a mid-IR integrating sphere (Pike Technologies).

BEM Average Scattering Simulations

Simulations of orientation- and polarization-average light scattering were performed using the SCUFF-EM application AVESCATTER.[30] SCUFF-EM[36] is an open-source software for electromagnetic simulations based on the BEM. The meshing of the objects is based on triangular panels and was carried by GMSH.[47]

Monte Carlo Simulations of Radiative Transfer

Radiative transfer simulations were performed by our own code for Monte Carlo simulations of unpolarized light. The algorithm consists of simulating the trajectories of many individual photons as they interact with particles and interfaces until they are either absorbed by particles or exit the simulation domain. The initial condition of each photon is given by the position and direction of the light source. At each simulation step, the optical path (Λphoton) and fate of a photon are estimated by selecting the shortest path between the particle’s scattering (Λsca) and absorption (Λabs), the absorption of the host (Λhost), or diffraction (ΛFresnel), whereand ξ is a random number between 0 and 1; ΛFresnel is given by the shortest distance between the photon and an interface. In materials with more than one kind of particle, Λabs = min(Λabs) and Λsca = min(Λsca), where Λabs and Λsca are, respectively, the absorption and scattering path from the particle n.

In the case of diffraction (Λphoton = ΛFresnel), a photon is either reflected or transmitted by a random selection, with the probabilities of each event proportional to the respective energy flux defined by Fresnel laws. If the photon is absorbed by a particle (Λphoton = Λabs) or the host material (Λphoton = Λhost), the event is terminated and the simulation continues with a new photon at the initial conditions. For a scattered photon (Λphoton = Λsca), the new direction is determined by[7]where g = ⟨μsca⟩.

In all our simulations, we considered a slab with a large surface area in order to represent a 2D problem. As a criteria, we selected the smallest surface area by which no photon escapes through the edges. Two large monitors, above and below the slab, measure the total reflectance and transmittance, respectively. The specular transmittance was measured with a third small monitor (1 nm × 1 nm) at a 1 mm distance below the slab. In all the simulations, we considered 1000000 photons per wavelength. For the validation of our code, refer to Supporting Information, section 5.