Modeling Dispersive Media in Finite-Difference Time-Domain Method for Radiative Cooling Applications

Radiative cooling offers a new strategy to address building cooling problems in an environmentally friendly way, compared with traditional cooling methods that require high-power consumption. The goal of radiative cooling is to maximize absorption at the mid-infrared (MIR) spectrum (8– $13~\mu m$ ). Thus, it is important to characterize the dispersive relationship of the radiative cooling materials (TiO2, Si3N4, SiC, SiO2) in the MIR regime. In this work, we use the multi-pole Lorentz dispersion equation to model frequency-dependent dielectric constants of radiative cooling materials, which requires particle swarm optimization to solve multidimension optimization problems. This modeling can be directly applicable to the auxiliary differential equation (ADE)-based finite-difference time-domain (FDTD) method for simulating emissivity of radiative cooling materials. Further, we limit the number of Lorentz poles to five, which minimizes the computational burden in ADE–FDTD. Then, we validate our FDTD results against an analytic prediction calculated by the transfer matrix method. Finally, we demonstrate a highly efficient radiative cooling structure based on a simple SiO2 grating. The proposed material modeling can accelerate the nano/microstructural design of the radiative cooling strategies.


I. INTRODUCTION
Radiative cooling offers a new strategy to address photovoltaic and building cooling problems in an environmentally friendly way. Traditional cooling methods require power consumption with external devices, which occupy a large portion of energy consumption in everyday life. On the other hand, passive cooling methods are broadly classified as bio-inspired cooling strategies [1]- [4] and artificial structurebased strategies [5], [6], both reflecting the solar spectrum while maximizing absorption in the mid-infrared (MIR) spectrum (8-13µm).
In this MIR spectrum, strong selective thermal emission reaches the unlimited cold heat sink of space (∼ 3K), utilizing a transparency window of the Earth's atmosphere where the radiation absorption is insignificant. A recent advance of nano-fabrication has been applied to daytime The associate editor coordinating the review of this manuscript and approving it for publication was Guido Lombardi . radiative cooling applications such as metamaterials [7]- [9], polymers [7], [10], [11], and photonic structures [12]- [14]. In order to further improve radiative cooling performance, it is essential to precisely engineer radiative cooling materials with advanced nano/microstructures [12], [15]. This brings challenges to numerical simulations of the radiative cooling materials. For example, rigorous-coupled-wave analysis (RCWA)) [16], [17], discrete dipole approximation (DDA) [18], and the transfer matrix method (TMM) [19], [20] can be used to solve photonic structures relatively fast. However, they are not preferable for studying complicated (or enormous) structures due to assumptions on their methods [21], [22]. On the other hand, solving full-wave Maxwell's equations provides accurate solutions, wherein the most representative examples of full-wave Maxwell solvers are the finite element method (FEM) [23], [24] or FDTD [11], [25]. The former utilizes unstructured mesh with steady-state analysis, which makes it powerful for analyzing accurate solutions for complex 3D geometries in discrete frequencies. However, the nature of steady-state analysis in FEM requires multiple simulations to obtain a broadband response of photonic structures. The latter solves the differential form of Maxwell equations over time and space. It employs a staircase mesh, which is known as Yee's grid. This offers potential for efficient parallelization in multicore computing. Also, the time-domain response obtained by a single simulation in FDTD can be easily converted to a broadband frequency response by Fourier transformation [22], [26]- [28] as shown in Fig. 1(b) (orange area). The main hurdle on using FDTD for a radiative cooling study is modeling dispersive materials for the MIR regime. The FDTD method requires converting frequency-dependent material parameters to an auxiliary time-domain differential equation for simulating dispersive materials [22].
In this work, we use the Lorentz dispersion model to fit the frequency-dependent dielectric constants of the radiative cooling materials via particle swarm optimization (PSO) [29]- [31]. The Lorentz dispersion equation, also known as the Drude-Lorentz oscillator model, involves modeling electrons as a damped harmonic oscillator. By superposing multiple Lorentz equations, we can efficiently model arbitrary frequency-dependent permittivities of the radiative cooling materials for FDTD simulations. A single Lorentz equation has plasma frequency, resonance frequency, and damping coefficient as its parameters, which can be fitted using the least-squares method [32] or others [33]- [37], while fitting parameters for multipole Lorentz equations is more challenging. Therefore, finding a fewer pole dispersion model that can fit the experimental complex permittivity data is significant, guaranteeing faster and memory-efficient FDTD simulations. We apply a new figure of merit to the PSO method to discover the Lorentz parameters in this work. Then, we compare our method against a conventional PSO method and genetic algorithm. Additionally, we compare dispersive FDTD simulations against analytic calculations based on the transfer matrix method (TMM) [19], [38] for broad angles and broadband to validate our result. This study may facilitate a realization of novel radiative cooling design by enabling faster and accurate simulations of the radiative cooling materials.

II. MULTI-POLE LORENTZ MODEL OF RADIATIVE COOLING MATERIALS
The ADE approach can be applied to simulate Lorentz dispersion materials [39], [40] and many others [41]- [43] in FDTD [22]. The Lorentz model is written as where ε r is relative permittivity, ε ∞ is a constant offset for the real part of the dielectric function, n is the number of pole, ε n is a relative value of the plasma frequency compared to the resonance frequency, γ n is damping coefficient, and ω n is the resonance frequency of the oscillator. The second term denotes Lorentz oscillator model with n th poles. Then, frequency-domain Maxwell's equations can be written as After inserting Lorentz model (Eq. 1) into Eq. 2, the Maxwell's equations become It can be simplified to  [44], Si 3 N 4 [44], SiC [45], SiO 2 [46], [47]) of dispersive material, respectively. The real part of relative permittivities of TiO 2 (a), Si 3 N 4 (c), SiC (e), SiO 2 (g). The imaginary part of relative permittivities of  The quality of our fitting against the most widely used evolutionary algorithms (conventional genetic algorithm, conventional particle swarm optimization). We assume that the conventional PSO and the conventional GA algorithms have a figure of merit (F) given by where P L is a polarization vector for Lorentz dispersion materials. After putting Eq. 5 into Eq. 4 we get The ADE approach is applied to Eq. 6 to simulate dispersive materials in time-domain forms [39], [50]. Now, we apply PSO to search for proper Lorentz parameters defining frequency-dependent dielectric constants of radiative cooling materials. We choose four representative radiative cooling materials (TiO 2 [44], Si 3 N 4 [44], SiC [45], SiO 2 [46], [47]) in this work, although the proposed method can be applied to other dispersive materials. The PSO is an evolutionary optimization algorithm that mimics the behavior of nature, such as ants and bees [29]. It iteratively finds a potential solution with moving particles that have social and cognitive factors. The number of dimensions of the PSO algorithm is the number of poles multiplied by three (ω n , γ n , and ε n ). The velocity of the particles is determined by n is a subsequent velocity of each particle in n th dimension, which determines the next position of each particle in n th dimension (x i n ). ω is an inertia constant, and C 1 and C 2 are the cognitive and social constant, respectively, where both the cognitive and social constants are set to 2. rand() is a random number generation function, and pBest and gBest are a particle best position and a global best position found in the previous iterations. To evaluate the fitness of the particles, we define the figure of merit function (F) as follow: where ζ is a summation of error between frequency-dependent dielectric constants of the experimental data and Lorentz equation, and , indicate the real and imaginary parts of the dielectric constants, respectively. Constant ξ is an additional parameter to prevent divergence of the figure-of-merit function when the permittivity of the data crosses zero. It has zero value for Si 3 N 4 and SiC, while we gradually increase ξ for other materials until the fitting error is minimized. The PSO algorithm tries to minimize F by adjusting particles' velocities based on Eq. 7. The number of Lorentz poles is first set to one; then, we gradually increase the number of the poles until the error (F) satisfies a terminate condition. Increasing the number of Lorentz poles adds a computational burden [43] in the ADE-FDTD method [51] and escalates the complexity of the fitting problem. Therefore, we present our dielectric fitting results within five Lorentz poles, even though 10 to 15 poles demonstrate better accuracy. Table 1 shows optimized Lorentz parameters for each material with four to five Lorentz poles. Figure 2 shows dispersion complex curve fitting [52] of radiative cooling materials, optimized by the PSO algorithm. The optimized Lorentz curves show minimal errors for TiO 2 , Si 3 N 4 , and SiC while SiO 2 shows slight error on 9 µm and 10 µm in its real part. This can be improved further by employing 10 to 15 Lorentz poles; however, we limit our fitting up to five poles to make a consistent computational load. Plus, other experimental permittivity data for SiO 2 [44], [53], [54] demonstrate slightly different dielectric constants on 9 µm-10 µm; therefore, it might be unnecessary to overfitting one specific experimental data. We also compare the quality of our fitting against the most widely used evolutionary algorithms (Conventional Genetic Algorithm, Conventional Particle Swarm Optimization) as shown in Tab 2. We assume that the conventional PSO and the conventional GA algorithms have a figure of merit (F) given by F = ω .

III. FDTD ANALYSIS OF RADIATIVE COOLING MATERIALS
This section presents numerical examples to validate our dispersion curve fitting of the radiative cooling materials.
We implement Lorentz parameters to MEEP [48], opensource FDTD solver, and simulate with a finite thickness of the dielectric slab. The FDTD time step is assigned by t = C n s/c 0 where C n is a courant number with values between 0 to 1, S is a spatial step size, which is often called a grid spacing, and c 0 is a speed of light.
We use s = λ min /50 and C n = 0.5. 5-µm and 2-µm-thick of the dielectric slabs are simulated with the FDTD method as shown in Fig. 3. Transmission of the incidence wave is calculated by accumulating the Fourier-transformed electric and magnetic fields for all points in the flux plane through the summation up to n th discrete time step:f wheref (ω) is a Fourier transformed field, t is the FDTD time step, and f (t) is a time-domain field. At the end of the simulation, these Fourier-transformed fields are used to VOLUME 10, 2022  compute the integral of the Poynting vectors to produce frequency-dependent transmission.
The red dots indicate Fourier-transformed FDTD simulation results, while the black curves indicate theoretical calculation based on TMM [19]. Analytic transmission through TMM can be calculated by where S = S 11 S 12 in Eq. 12, E + 0 is the global electric field, and E + is the transmitted electric field. In Eq. 13, S is the global scattering matrix. Note that there is a slight error on the transmission of the TiO 2 dielectric slab shown in Fig. 3(a). This error can be improved with 10 Lorentz poles.
To understand emissivity of radiative cooling materials, one may need to calculate angle-dependent absorptivity of the materials. Then, we can apply Kirchhoff's law, which states that, for a body in equilibrium, the emissivity is equal to the absorptivity at each frequency and polarization [55]- [57]. Alternative methods utilize uncorrelated white-noise sources [55], [58] or trace formulation [59] to Maxwell's equation. In this work, we employ Kirchhoff's reciprocity to calculate the emissivity of the radiative cooling materials. The incidence planewave is excited at the top of the dielectric slabs, while transmission and reflection monitors are placed on both sides. Then, the absorptivity can be calculated via A = 1 − T − R. Kirchhoff's law states that absorptivity for an incidence planewave with a specific angle equals the emissivity to that angle. Now, to analyze angle-dependent emissivity over the MIR regime, we excite a Bloch boundary planewave with different wavenumbers x + k 2 y in two dimensions. The Bloch boundary only supports a single wavenumber for each dimension (e.g., k x , k y ); thus, a single FDTD simulation with the specific Bloch boundary condition corresponds to different incidence angles at different frequencies.
Angle(θ)-dependent emissivity is shown in Fig. 4. We use TMM for an angle-dependent analytic calculation of the emissivity, as shown in Fig. 4(b,d,f,h). The simulated emissivity is calculated via oblique incident Bloch boundary planewaves (total 81 simulations) and then interpolated to generate surface plots shown in Fig. 4(a,c,e,g). The simulated angle-dependent emissivity matches the analytic calculations with less than 1% root mean square error. Note that the gray area shown in the simulated emissivity plots is not calculated because it requires a high-angle and a greater wavenumber, where both of them increase a required k x , which can exceed the wavenumber (k), meaning the nonpropagating mode. This can be simulated by adjusting the center wavelength of the planewave source.
A dielectric slab of SiO 2 shows low emissivity in the 8.0 to 9.5µm wavelength range due to high reflectivity and the low absorption coefficient. To further improve the emission characteristic of SiO 2 , we design a 2D grating structure, as shown in Fig. 5(a). The grating or texturing structure improves the light coupling of thermal emission [33], decreasing reflection while redistributing the angular emission profile. Figure 5(a) shows a simulation structure, which includes 2-µm size SiO 2 based grating. As shown in Fig. 5(b), the simulated broadband angle-dependent emission shows that 8.0 to 9.5µm wavelength emissivity is greatly improved compared with a planar SiO 2 structure shown in Fig. 4(g). The angle-dependent emissivity shown in this work provides a pathway to accurately model and engineer novel micro/nanostructures for radiative cooling studies. Thus, we believe that this result can be applied to developing advanced radiative cooling strategies. For a realization of VOLUME 10, 2022 this structure, one can use a holographic technique with a mode-locked Ti:sapphire laser [60] or utilize a commercially available SiO 2 grating product.

IV. CONCLUSION
We have modeled dispersive media used in radiative cooling applications for FDTD simulations. The PSO method is used to determine multipole Lorentz parameters for TiO 2 , Si 3 N 4 , SiC, SiO 2 materials. Four to five Lorentz poles are used for modeling radiative cooling materials. Then, the modeling results are validated against theoretical calculations based on the transfer matrix method. Then, multiple FDTD simulations predict the angle-dependent broadband response of the radiative cooling materials. Finally, we simulate the SiO 2 grating structure to improve a planar SiO 2 emissivity. The proposed dispersion model may provide a path to efficiently discover an advanced radiative cooling structure.