Spectral emissivity modeling in multi-resonant systems using coupled-mode theory

: The ability to design multi-resonant thermal emitters is essential to the advancement of a wide variety of applications, including thermal management and sensing. These fields would greatly benefit from the development of more efficient tools for predicting the spectral response of coupled, multi-resonator systems. In this work, we propose a semi-analytical prediction tool based on coupled-mode theory. In our approach, a complex thermal emitter is fully described by a set of coupled-mode parameters, which can be straightforwardly calculated from simulations of unit cells containing single and double resonators. We demonstrate the accuracy of our method by predicting and optimizing spectral response in a coupled, multi-resonant system based on hBN ribbons. The approach described here can greatly reduce the computational overhead associated with spectral design tasks in coupled, multi-resonant systems.


Introduction
The ability to tailor spectral emissivity in the infrared is useful for a wide variety of applications including thermal management [1][2][3][4], energy harvesting [5,6] and sensing [7][8][9]. To this end, microstructured emitters utilizing resonant photonic modes have been investigated [10][11][12]. Localized photonic modes, such as those supported by metal-insulator-metal (MIM) structures, have been used to synthesize emissive spectra [13,14]. By combining localized resonators with different sizes, one can generate spectral peaks at multiple, desired wavelengths of interest [15]. More recent work has demonstrated the use of coupled, many-resonator systems to significantly increase the complexity of the synthesized spectrum [16]. The computational technique was based on machine learning and relied on training data generated from multi-resonator simulations [17].
Recently, an alternative approach to thermal emitter design based on temporal coupled-mode theory (CMT) has been proposed [18,19]. In this approach, a catalog is first built based on simulations of single resonators and a model of their coupling. The catalog is then used in conjunction with a semi-analytical model to calculate the electromagnetic quantities of interest. A major advantage of this approach is that once the catalog is constructed, the method can be used to maximize various objective functions with very little additional, computational cost [19]. Previous work has demonstrated the use of this method to design self-focusing thermal emitters and thermal holograms [18,19]. Here, we demonstrate that the method can be adapted for prediction of thermal emission spectra. Moreover, the use of a semi-analytic model allows for valuable insight into the relationship between structural parameters and spectral response that is not always possible with purely numerical methods.
In particular, we use a CMT approach to predict the emission spectra of hexagonal boron nitride (hBN) thermal emitters. Periodic emitters based on 2D materials such as hBN have been shown to exhibit high quality factor (Q factor) resonances in the mid-infrared [20]. Emitters with several uncoupled hBN ribbons per unit cell have been used to generate a rich, multi-resonant spectral response. Moreover, recent work has demonstrated that coupling between resonances serves as an additional degree of freedom to generate non-trivial spectral features [21]. By accounting for coupling between resonances, our model is able to accurately predict the spectral response of complex thermal emitters. While we demonstrate the utility of our approach using hBN ribbon emitters in this study, it is applicable to a wide class of coupled-resonance systems.

Methods
Our thermal emitter consists of a periodic array of hBN ribbons placed on a 1.4 µm thick MgF 2 layer and a semi-infinite Ag back reflector ( Fig. 1(a)). For generality, we consider an emitter with several hBN ribbons per unit cell generating a multi-resonant emission spectrum. Figure 1(b) shows the unit cell of an emitter consisting of N hBN ribbons oriented parallel to the x-y plane with lengths L 1 through L N . The inter-ribbon separations are denoted by d i , where i ∈ [1, N − 1] and the unit cell length is a. In this work, we assume each hBN ribbon to be composed of 6 atomic monolayers (corresponding to a total thickness of 2 nm).   1(c) presents the sequence of steps used to predict the spectrum of the emitter described in Fig. 1(b). We first begin by building a catalog of CMT parameters from simulations using the method of Rigorous Coupled Wave Analysis (RCWA) [22], which solves Maxwell's equations with no uncontrolled approximation. A periodic array containing a single hBN ribbon per unit cell gives rise to a photonic resonance. Simulations of the array with one ribbon per unit cell are used to determine the resonance frequency and decay rates as a function of the ribbon length L and a. For an array with two ribbons per unit cell, we expect the system to support two coupled photonic resonances. We perform a simulation containing two ribbons per unit cell to extract the coupling constant as a function of L 1 , L 2 , a, and d, the separation between ribbons. Using this catalog of CMT parameters, one can in principle predict the spectrum of any arbitrary N-ribbon thermal emitter.

Results
In the following subsection, we discuss the process of extracting CMT parameters from simulations and the underlying theoretical framework in more detail. The resulting CMT catalog is utilized in subsequent subsections to predict emission spectra of more complex thermal emitters.

Building a catalog of CMT parameters
We begin by calculating the single-resonance CMT parameters from simulations of one hBN ribbon per unit cell. Figure 2(a) shows the unit cell. The length of the hBN ribbon is L and the unit cell length is a. As an example, we compute the emission spectrum of this emitter for L = 20 nm and a = 0.2 µm using RCWA (dashed blue curve in the right panel of Fig. 2(a)). The optical constants of hBN are given by a bulk model that considers the existence of optical phonon modes in the mid-infrared [23]. Previous work has shown that using the bulk model of hBN yields qualitatively similar spectral features as the monolayer model [20]. Additionally, as a result of the atomic thickness of hBN used in this study, its dielectric function along the z direction can be neglected and optical response is accurately described by the in-plane (x-y plane) phonon modes. The dielectric function of Ag is based on a Drude model [24] while MgF 2 data is taken from literature [25]. The reflection spectrum is calculated for a normally-incident TM-polarized (H-field in the y-direction) plane wave. Details of the RCWA algorithm used in this study are provided in Ref. [22]. Due to the presence of a semi-infinite silver layer, the emitter does not have any transmission, and emissivity is given as 1 -reflectivity. The spectrum shows an emission peak close to 6.9 µm corresponding to the hBN ribbon resonance. To isolate the resonant contribution, we subtract a constant value of background emission and plot the resultant adjusted emissivity data by the green, dashed curve.
To fit the emissivity, we use temporal coupled mode theory (CMT): Here ω 0 , γ a and γ r denote the resonant frequency, absorptive decay rate and radiative decay rate of the resonance, respectively. Fitting Eq. (1) (solid black curve in the right panel of Fig. 2(a)) to the adjusted emissivity data gives γ a = 4.77 × 10 11 rad/s and γ r = 2.21 × 10 10 rad/s. The numerical values show that the thermal emitter is under-coupled: the absorptive decay rate is larger than the radiative decay rate [26]. This results in the emitter having a low peak emissivity of about 0.17.
The next step is to calculate the coupling constants between resonances. We consider two identical hBN ribbons per unit cell, separated by a distance d per unit cell (left panel of Fig. 2(b)). The length of the ribbons, L, is assumed to be 20 nm and the unit cell length is 0.2 µm. This system is described by the following CMT equations: Here a 1 and a 2 represent the mode amplitudes of the two resonances, β is the complex coupling constant, and s + and s − are the amplitudes of the incoming and outgoing plane waves respectively. For identical ribbon lengths, we assume that the resonant frequency ω 0 and decay rates γ a and γ r are equal for both resonances. Solving Eqs. (2), (3) and (4), the emission spectrum for the system is given as: Here Re(β) and Im(β) denote the real and imaginary parts of β respectively. Equation (5) represents a Lorentzian centered at ω 0 + Re(β). Even though the system supports two modes, one of them is dark as a consequence of symmetry and does not contribute to emission. Therefore, the emission spectrum is given by a single Lorentzian corresponding to the bright mode contribution. One can observe that Re(β) appears only in the frequency term of Eq. (5) and Im(β) appears in the decay rate terms. Qualitatively, one can interpret Re(β) as the frequency shift of the bright mode of the system relative to the bright mode of a system with a single ribbon per unit cell.
We simulate the emitter with two ribbons per unit cell shown in Fig. 2(b) and fit the adjusted emission spectrum using Eq. (5). The resonant frequency and decay rates are determined from the CMT analysis presented in Fig. 2(a), leaving Re(β) and Im(β) as the only free parameters. The right panel of Fig. 2(b) shows the adjusted simulation data (dashed green curve) and fitted CMT data (solid black curve) for d = 5 nm. Our model provides a good fit of the simulation data. The values of β determined from fitting are presented in Fig. 2(c) for different values of d. It can be observed that the plot is mirror-symmetric about d = 80 nm. This is due to the periodicity of the ribbon array. The unit cell of the emitter shown in Fig. 2(b) is equivalent to one in which the ribbons are separated by a distance (a -2L)d = 160 nmd. Consequently, the β value for a separation of 100 nm, for instance, is the same as that for a separation of 60 nm. Additionally, one can observe that Im(β) is two orders of magnitude smaller than Re(β). Therefore, we neglect its dependence on d and use its average value instead.

Spectral prediction for identical ribbons
Next, we test the utility of our CMT model in predicting the spectral response of more complex thermal emitters. In general, an emitter comprising of N coupled hBN ribbons per unit cell can be regarded as a set of N coupled resonances. Such a system is described by the following CMT equations [27]: Here |a⟩ = (a 1 a 2 · · · a N ) T represents the mode amplitudes of the N resonances. The matrix Ω includes the resonant frequencies and coupling constants between the resonances while Γ represents their decay rates. These matrices are written as: Here ω i , γ ai and γ ri represent the resonant frequency, absorptive decay rate and radiative decay rate of the ith resonance. We assume β ij = β ji (where i, j ∈ [1, N]). The emitter is coupled to the incoming and outgoing plane waves via the vector D = j (︂ √︁ 2γ r1 √︁ 2γ r2 · · · √︁ 2γ rN )︂ . From Eqs. (6) and (7), the emission spectrum of the system is given by: Here I is an NxN identity matrix. We use Eq. (8) with N = 3 to predict the spectrum of an emitter comprising of three hBN ribbons per unit cell ( Fig. 3(a)). The ribbons are identical with length L = 20 nm and are separated by a distance d. The unit cell length a = 0.2 µm. We assume the system can be described by three coupled, identical resonances. The resonant frequency and decay rates are obtained from the calculations with one ribbon per unit cell presented in Fig. 2(a). The coupling constants are taken from Fig. 2(c). The right panel of Fig. 3(a) presents a comparison between the adjusted emission spectrum obtained from simulations (dashed green curve) and CMT prediction (solid blue curve) for d = 1 nm and 5 nm. It can be observed that the spectra predicted from our CMT model agree well with the simulation results. The spectra show a high long-wavelength emission peak and a low short-wavelength peak corresponding to the two bright modes of the system. We refer to these as mode 1 and 2 respectively. Our model accurately predicts the mode 1 peak for both d = 1 and 5 nm, while predicting a mode 2 peak that is blue-shifted with respect to the simulation result for d = 1 nm.
In Fig. 3(b), we present the variation of the resonant frequencies of modes 1 and 2 with d. At small d, the resonators are strongly coupled and the bright modes of the system are well separated in frequency space. As d increases, the coupling between the resonators reduces resulting in decreased frequency-splitting between the modes. At large separations, the emission response of the system is governed by the bright mode of an emitter with a single ribbon per unit cell. In Fig. 3(b), we display the resonant frequency of mode 2 up to d = 10 nm; above this value the emission peak completely merges with that of mode 1 at larger separations. Our CMT model correctly predicts the resonant frequencies of modes 1 and 2 for d ≥ 5 nm. Therefore, it can be used to accurately predict the spectral response of emitters with inter-ribbon separations greater than 5 nm.

Spectral prediction for non-identical ribbons
We next use our CMT model to predict the spectral response of emitters comprised of three non-identical hBN ribbons per unit cell. Such an emitter can be described by a system of three non-identical resonances. We first determine the coupling constants for non-identical resonances. The ribbons in each unit cell have lengths L 1 and L 2 and are separated by a distance d (Fig. 4(a)).  We use Eq. (8) with N = 3 to predict the spectra of emitters comprised of three non-identical hBN ribbons per unit cell (Fig. 5(a)). The ribbons in each unit cell have lengths L 1 = 20 nm, L 2 = 30 nm and L 3 = 40 nm and are separated by distances d 1 and d 2 . The length of the unit cell is 0.2 µm. For simplicity, we only consider d 1 , d 2 = 20 and 5 nm (weak and strong coupling, respectively). Figure 5(b) shows the emission spectrum of the emitter for different values of d 1 and d 2 . The dashed black lines correspond to the peak locations of the emitter in the weak-weak coupling case. In the strong-strong coupling case, the emission peaks for the first and third resonances undergo a blue and red-shift respectively with respect to their weak coupling locations. The peak for resonance 2 undergoes a slight blue shift as a result of strong coupling to both resonances 1 and 3. In the weak-strong coupling case, the peaks for resonances 2 and 3 move apart from their weak coupling locations, while the peak location for resonance 1 remains unchanged. Similarly, for strong-weak coupling, the peak location for resonance 3 remains unchanged while the peaks for resonances 1 and 2 move apart compared to the weak-weak coupling case. Therefore, by simply changing the coupling between adjacent hBN ribbons, the spectral response of this emitter can be tuned. In all four cases, spectra predicted using our CMT model match well with the simulation results. While a very slight shift of the CMT prediction with respect to simulation can be observed in Fig. 5(b), the absolute value of the difference between the curves was less than 0.07 for all cases shown. The CMT model can thus accurately capture the effects of coupling between resonances on the spectrum.

Computational efficiency of the CMT model
The CMT approach allows for fast optimization of multiple-ribbon systems. To illustrate the computational speedup provided by CMT, we consider maximization of emissivity at a single wavelength of interest. For simplicity, we consider an emitter composed of N identical hBN ribbons per unit cell ( Fig. 6(a)) and vary the ribbon length L to maximize emissivity at 7µm. We fix the inter-ribbon separation d = 5 nm and choose the length of the unit cell a = 0.2 µm.  Figure 6(b) shows the computational time required to calculate a single emissivity spectrum, given ribbon length L and number of resonators per unit cell N. For the RCWA simulation, the required time was approximately 25s on the computer used. For the CMT simulation, the calculation was more than 1000 faster on the same computer (note that the blue curve has been multiplied by 1000x to fit on the same graph). In calculating this time, we assume that the resonance parameters and coupling constants were already known, due to pre-computation of an appropriate catalog. The pronounced speedup provided by CMT thus arises because evaluation of Eq. (8) requires far fewer calculations than the full-field simulations of the electromagnetic fields performed in RCWA.
Using the CMT method, we next optimized L for each value of N. Given a pre-computed catalog of L-dependent resonance parameters and coupling constants, we used the MATLAB patternsearch function to optimize L with respect to the target function ε(λ = 7µm), evaluated via Eq. (8). The results for optimal L and maximized emissivity at the target wavelength are shown in Fig. 6(c) and 6(d), respectively. To verify the prediction, we compared the spectrum for a structure optimized by CMT to direct simulation. Figure 6(e) plots the adjusted emission spectrum of the optimized emitter with N = 8 obtained from CMT simulations (solid blue curve). The spectrum has a tall emission peak close to 7 µm and a very short peak at about 6.95 µm. This result shows good agreement with the spectrum given by direct, RCWA simulations (dashed green curve).
A similar method could be used to optimize other target functions of interest. For example, one could maximize or minimize emissivity at another wavelength or define a target spectral function ε(λ) for the multi-ribbon system (e.g. top hat, step function, etc.). After an initial, up-front computational investment in building an appropriate catalog (e.g. containing multiple ribbon lengths and separations), the catalog can then be reused to perform a variety of optimization tasks.

Conclusion
We proposed a CMT-based model for predicting the emission spectra of complex thermal emitters. In this study, we investigated emitters consisting of a periodic array of hBN ribbons on top of an MgF 2 layer and a semi-infinite Ag back reflector. We began by building a catalog of CMT parameters by simulating thermal emitters with single and double hBN ribbons per unit cell. This was subsequently used to predict emission spectra of structures consisting of three hBN ribbons per unit cell. We showed that our model can accurately predict the emission response of such emitters for inter-ribbon separations greater than 5 nm. To demonstrate the computational efficiency of our CMT model, we considered the task of maximizing emissivity at a desired wavelength. We found that the CMT method reduces the computation time by three orders of magnitude relative to RCWA, provided that a catalog has already been constructed.
These results thus suggest a number of interesting directions for future work. While we have focused on hBN ribbon emitters in this study, we expect our approach to be applicable to a broad class of metallic and dielectric microstructures supporting photonic resonances. Moreover, our approach is suited for a variety of spectral optimization tasks. The most significant computational cost comes from building the catalog. Once built, the same catalog can be reused to optimize a given objective function at minimal additional cost. For example, one could maximize or minimize emissivity at single or multiple wavelengths of interest, or even determine the collection of emitters that gives a "best fit" to a specified emissive spectrum. We thus expect that the orders-of-magnitude reduction in computation time offered by the CMT model can potentially allow optimization within a previously inaccessible design space.

Disclosures. The authors declare no conflicts of interest.
Data Availability. Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.