Numerical methods for design of metamaterial photonic crystals and random metamaterials

. Two-dimensional metamaterial photonic crystals (2DMPCs) composed of dispersive metamaterials in a positive-refractive-index medium were investigated by incorporating ﬁ nite-difference time-domain calculations into the auxiliary differential equation method. A distinct band gap was formed and the effects of positional and size disorder when the dispersive metamaterials are aligned in air were elucidated. In addition, using the self-consistent ﬁ nite-difference frequency-domain method, an eigenmode analysis of 2DMPCs with positional disorder was performed. Finally, a numerical methodfor the inverse design of binary randommetamaterial multilayers was proposed.


Introduction
Electromagnetic wave propagation in negative-refractiveindex materials embedded within positive-refractive-index materials is attracting significant interest because it raises fundamental questions about the physical properties of electromagnetic waves [1][2][3][4][5]. In early studies on the topic, one-dimensional multilayers consisting of alternating slabs of positive-and negative-refractive-index materials [6] and two-dimensional metamaterial photonic crystals (2DMPCs) [7][8][9] were investigated. In a previous work [7], I investigated 2DMPCs composed of dispersive left-handed materials, considering the Drude-type dispersion responses for the dielectric permittivity and magnetic permeability of the metamaterial. I found that a stop band appears in the case of transverse magnetic (TM) polarization when the dispersive left-handed metamaterials embedded in the positive-refractive-index medium are spatially isolated from each other.
During experiments, materials usually exhibit some degree of disorder. This positional disorder can have a significant impact on properties such as the photonic density of states (DOS) of photonic crystals (PhCs) as well as the performance of PhC-based devices. Numerical and experimental studies [10][11][12] on disordered photonic materials have shown that system periodicity or quasiperiodicity is not necessary for the formation of a photonic band gap. In reference [13], the DOS of two-dimensional photonic crystals with dielectric rods were numerically investigated under the condition of positional disorder. These results suggested that band tails and gap states are present in random systems.
In this study, I incorporated finite-difference timedomain calculations into the auxiliary differential equation (ADE-FDTD) method to analyze a 2DMPC composed of an array of circular rods with disorder. I performed numerical simulations to quantitatively analyze how this type of system disorder affects the characteristics of electromagnetic waves within the system. In the analysis of electromagnetic waves, finite-difference frequency-domain (FDFD) method is considered to be a promising approach [14][15][16][17]. Using the self-consistent FDFD method, an eigenmode analysis of a 2DMPC with positional disorder was also performed. Finally, a numerical method for the inverse design of binary random metamaterial multilayers is proposed. Section 2.1 describes the results of the calculations of the band gap in the disordered 2DMPC and Section 2.2 describes the results of a large-scale eigenmode analysis of the disordered 2DMPC. In Section 3, a numerical method suitable for the optimal design of random metamaterial multilayers is described. Finally, Section 4 states the conclusions of the study.
2 Two-dimensional metamaterial photonic crystals (2DMPCs) with disorder 2.1 Analysis of metamaterial composites using ADE-FDTD method In this section, the modeling of 2DMPCs consisting of a periodic array of parallel metamaterial rods with circular cross-sections in air is described. The Drude model for dispersive metamaterials, in which both the permittivity, e (v), and permeability, m (v), are frequency dependent, was employed: and where e 0 and m 0 are the permittivity and permeability, respectively, in vacuum; v pe and v pm are the electric and magnetic plasma frequencies, respectively; and g e and g m are the electric and magnetic collision frequencies, respectively [7]. I considered a 2DMPC consisting of a periodic array of circular rods with radius r. The lattice constant of the square (and triangular) lattice was a (Fig. 1). It was assumed that, in this two-dimensional system, electromagnetic waves propagate in a plane perpendicular to the rods, and two types of polarizations of the waves were considered: transverse electric (TE) polarization and TM polarization. In this section, I discuss the TM polarization case.
In the following section, I consider two different types of disorders within the system. The first is positional disorder, d, i.e., disorder in the positions of the metamaterial rods, and the other is size disorder, D, which is related to the radius of the rods (Fig. 2). In the case of the positional disorder, the random deviations (d i,x and d i,y ) of the i-th rod from the corresponding lattice position in the x-y plane exhibit the relationship ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi d 2 i;x þ d 2 i;y q < d and are uniformly distributed (Fig. 2a). With respect to the size disorder, the radius, r i , of the i-th rod exhibits a uniform distribution within the range r À D < r i < r + D (Fig. 2b).
The ADE-FDTD method was employed to treat the dispersive medium [18]. The relationship between the electric flux density,D (r, t), and the electric field, E (r, t), is given as follows: Similarly, the relationship between the magnetic flux density, B (r, t), and the magnetic field, H (r, t), is given by For a given wavevector, q, Bloch's periodic boundary conditions were considered such that where R is an arbitrary lattice vector in the unit cell. In this study, I considered the power spectrum, P (v), which corresponds to the square of the photonic DOS. At t ≈ t 0 , an impulse is imposed on the system such that where DE i 0 t ð Þ is the z-component of the electric field imposed on a randomly chosen cell, i 0 , which is introduced by the spatial decomposition by the Yee cell [14,15]. The differential operator in equation (6) is required to eliminate the "DC component" of the impulse. After the passage of a sufficient amount of time (t s >> t 0 + t p ), the time correlation function, f (t), can be calculated as follows: The power spectrum, P (v) , of the system can be calculated from Here, ⟨ ... ⟩ {q} denotes the average of the set of wavevectors, q, in the Brillouin zone. The peak position of the frequency spectrum of the time-correlation function is obtained for each wavevector, q, and the power spectrum, P (v), of the system can be obtained from the set of wavevectors in the Brillouin zone. Calculations corresponding to different wavevectors, q, can be performed independently and are suitable for parallel computations. In equation (8), a window function, w (t), is introduced to obtain a smooth power spectrum with a finite impulse response for spectral analysis. This function is given by In general, the value of the power spectrum, P (v), becomes zero in the frequency regime, where the stop band of the system exists.
In a previous study, I confirmed the existence of band gaps in 2DMPCs based on periodic arrays of cylindrical or square columns of metamaterials [7,19]. Specifically, I confirmed that the frequency at which the band gap appears increases as the size of the metamaterial rods or square pillars is increased. This suggests that the frequency at which the average refractive index of the system becomes zero increases with the metamaterial fraction.
Next, I discuss the results of the analysis of the effects of randomness on 2DMPCs. To confirm the accuracy of the calculation results, independent runs were performed by changing the source position, i 0 , of the electromagnetic wave. In the random system (d ≠ 0 or D ≠ 0), different samples were used for the positional and size disorder conditions. It was confirmed that there was no difference between the properties of the band gaps of the two independent samples. In the figures described below, I plot the average results for these samples. Figure 3a shows the frequency dependence of the power spectrum, P (v) , for TM polarization in the square-lattice 2DMPCs, for which the square lattice constant is a. Furthermore, the radius of the rods, r, was 0.17a, and the positional disorder was d.
The supercell system size was fixed at L x = L x = 5a, and the number of metamaterial rods was set to 25. The parameters for the electric and magnetic Drude models were given by v pe = v pm = 1.0 (c/a), and low loss parameters were assumed by setting g e = g m = 10 À5 (c/a). The mesh resolution for the ADE-FDTD calculations was taken to be 0.01a. Figure 3b shows the calculation results for P (v) for a triangular lattice with positional disorder d, where the system size of the supercell was defined as L x = 4.0a and L y = 3.46a, and the number of metamaterial rods was set to 16. Figure 3c shows the calculation results for P (v)for a square lattice with size disorder D, where the system size of the supercell was L x = L x = 5a, and the number of metamaterial rods was 25. The results shown are the averages of three independent samples. The CPU time required for the calculation of each sample was 400 node hours (the product of the number of nodes and the time in hours); these were performed on a FUJITSU PRIMEHPC FX1000 supercomputer. As shown in the figures, the existence of a 2D stop band was confirmed. It was also confirmed that the frequency region of the stop band was insensitive to the presence of positional disorder in the system. This result is apparently in contrast to the properties of conventional photonic crystals reported previously [13]. However, this result is in keeping with that of a previous work on acoustic metamaterial systems [20], where the appearance of a band gap in the vibrational DOS was observed for random acoustic metamaterials containing resonant particles whose effective mass was negative within a particular frequency regime. In this case, a band gap appeared even when the resonant particles were not periodically arranged but were distributed randomly within the medium. These results indicate that band gap formation in photonic and phononic crystals is governed by the Bragg reflection based on the periodicity of the system, while band gap formation in metamaterial crystals (metamaterial composites) depends on a different principle. During experiments, it is difficult to ensure that the system is completely free from imperfections. This is useful behavior for practical applications because the band gap can be made to appear simply by randomly embedding metamaterials within the medium in question without having to control the disorder in the system.

Eigenmode analysis using self-consistent FDFD calculations
While developing numerical methods for metamaterials, it is necessary to consider systems in which the permittivity and permeability exhibit a dependency on the complex frequency. Many of the proposed methods are limited to cases in which the frequency dependence of the T. Terao: EPJ Appl. Metamat. 9, 1 (2022) permittivity (permeability) follows a specific functional form such as the Debye, Lorentz, or Drude dispersion relation, and the piecewise linear recursive convolution and ADE-FDTD methods have been proposed for such dispersive media [18]. In the case of metamaterials, however, the frequency dependence of the permittivity and permeability is extremely complex, and it is difficult to apply conventional calculation methods directly. Numerical methods that are suitable for use with general dispersive media will make it possible to accurately analyze the behavior of electromagnetic waves within metamaterials. For example, the experimental data regarding the frequency dependence of the permittivity and permeability can be directly input into simulation software for analyzing the dispersion relationship and eigenmodes.
Using the FDFD method, it is possible to study the properties of electromagnetic waves in dispersive media with arbitrary frequency dependence [8,17]. An outline of the FDFD method for systems with two-dimensional periodicity is now described in brief. I begin with the Maxwell equations: where E and H denote the electric and magnetic fields, respectively. It is assumed that the electromagnetic waves propagate in a plane perpendicular to the rods. In addition, in this section, I treat the TE polarization cases, i.e., E = (E x , E y , 0) and H = (0, 0, H z ) polarization. When the current density, j, is equal to zero, the amplitude, H z (r), of the magnetic field satisfies the relationship, in the case of TE polarization, where e (r, v) and m (r, v) are the relative permittivity and relative permeability, respectively, at position r, and angular frequency v (= 2pf), and c is the speed of light in vacuum. At this stage, e (r, v) and m (r, v) are both considered to be frequency independent, as explained in the Introduction section. The numerical treatment of photonic crystals with dispersive metamaterials is discussed in Section 4. The transverse electric fields, E x (r) and E y (r), are given by where e 0 is the dielectric constant in a vacuum. In this case, the eigenvalues and matrix elements for the eigenvalue problem need to be solved self-consistently because the matrix contains the values of the permittivity and permeability, which depend on the frequency, f, of the incident wave. This calculation method is subsequently called the self-consistent FDFD method to distinguish it from the conventional FDFD method [14]. In the case of FDFD methods, it is necessary to solve the internal eigenvalue problem, where the coefficient matrix is a large sparse matrix (and in general, a non-Hermitian) matrix. However, in previous numerical studies that used the FDFD method, the system (matrices) was not large. During the numerical analysis of random metamaterial systems, as in this study, it is essential to treat larger matrices in order to make the supercell large enough.
Several methods are available for determining the eigenvalues of large sparse matrices [21][22][23]. In this study, the implicitly restarted Arnoldi method with the shift-andinvert technique was adopted [24][25][26]. I used the selfconsistent FDFD method to obtain the mode patterns of electromagnetic waves in a 2DMPC with positional disorder. In this section, I adopt the metamaterial model from the previous section. The permittivity and permeability of the surroundings (air) were e A and m A , respectively, and the permittivity and permeability of the metamaterial were e B (f) and m B (f), respectively, such that where f (GHz) is the frequency. To allow for a comparison with previous studies, the models for block B were the same as those used in previous studies [6,27]. Figure 4a shows a supercell containing metamaterial rods arranged in a square lattice with positional disorder (d = 0.25a, a = 6.0 Â 10 À2 m). The length of one side of the supercell is 3.0 Â 10 À1 m, and the radius of the metamaterial rods is 2.5 Â 10 À2 m. While using the FDFD method, the spatial partition of the system was set to 10 3 Â 10 3 , and the dimension of the matrix to be diagonalized was 10 6 .
In the eigenmode analysis based on the self-consistent FDFD method, the FDFD calculation is performed considering the permittivity and permeability values of the metamaterial at multiple frequencies f ¼ f 1 ; f 2 ; :::; f m in the frequency range in which the electromagnetic wave properties are to be investigated. The CPU time required for a single eigenvalue calculation in this process was about 40 min when a single node of the FUJITSU PRIMEHPC FX1000 system was used. If the obtained frequency f satisfies f ¼ f , it gives a solution in the original (nonlinear) eigenvalue equation. The eigenvalue analysis of the system of parallel metamaterial rods was carried out, and the calculated results were shown in Figure 4b; the obtained frequencies f were plotted on the ordinate, on the abscissa, f (GHz) was plotted, and the dotted line in the figure indicates the condition f ¼ f .
Eigenmodes with eigenfrequencies, f, of 2.97998 and 2.99896 GHz are displayed in Figures 4c and 4d, respectively. In the figures, the region in which the (real) value of H z (r) is positive is displayed in red, and the region in which it is negative is displayed in blue. Based on the calculation results, it was confirmed that the amplitude of the electromagnetic field changed rapidly in the vicinity of the metamaterial rods. Hence, through a series of calculations, it was confirmed that the implicitly restarted Arnoldi method is effective for performing calculations based on self-consistent FDFD methods.

Numerical method for optimizing design of binary metamaterial multilayers
As stated in the previous section, it is of interest to elucidate the physical properties of random metamaterials. In this section, I discuss numerical methods suitable for use with random metamaterials. The transmission of waves in various metamaterial multilayers has already been studied by several groups [6,[27][28][29][30][31][32][33]. I address the problem of optimizing a multilayer structure that satisfies the required transmission spectrum of electromagnetic waves (Fig. 5). By selecting an appropriate arrangement of the positiverefractive-index material (A-block) and a layer of the metamaterial (B-block), I designed a multilayer structure that almost satisfies the required transmission spectrum, T target (f). This is an inverse design problem, and the positions and numbers of A-and B-blocks can be chosen freely, with the total number of blocks, N, being fixed. e A and m A in block A and e B (f) and m B (f) are given by equations (14), (15), (16), and (17).
In this study, I focused on the similarity between the above-described problem and the Ising model in statistical physics. In these systems, the A-blocks correspond to the upward spins of the Ising model and the B-blocks to the downward spins. Thus, one can treat this problem as the optimization of the problem of searching for the minimum of the cost function (which corresponds to the total energy in the Ising model), where the cost function, C, is defined as the deviation in the transmission spectrum from the target function. It is impossible to solve this problem with accuracy because the number of possible cases in the system is too large. For this purpose, a number of metaheuristic approaches have been proposed to obtain approximate solutions, and one example is simulated annealing. In such problems, the cost function has a complex rugged landscape with many local minima. Therefore, it is easy to get trapped by these local minima while solving the optimization problem.
In accordance with the conventions of statistical physics, I use the term temperature in the following section; however, the term does not mean the temperature of the system as defined by thermodynamics but simply refers to a dimensionless parameter that controls the optimization procedure. One strategy for avoiding getting trapped by the local minima is to allow the temperature to increase in order to jump over the energy barrier and search for another solution (tempering). Typical examples include the simulated tempering method [34,35] and the parallel tempering method [36,37]. These methods have been applied to various physical problems such as those related to Ising spin glasses and biopolymers [37,38]. Because parallel tempering in statistical physics is mainly applied to free energy analyses of the systems in question, it was necessary that the calculations spanned the low-energy regime (when the cost function is small) and the highenergy regime (when the cost function is large). For the present problem, however, only the ground state (the smallest cost function) was required, and a computational method specialized for this purpose would be desirable.
Given this perspective, I developed a method called the restart tempering method. In this method, three different temperatures, namely, T l , T m , and T h (T l < T m <T h ), are considered, and the temperature of the system is changed as T h ! T m ! T l ! T h ! ... , with the number of interval steps at each stage, N t , being fixed. At the low-temperature stage (T = T l ), the agent searches for a local minimum. At the medium-temperature stage (T = T m ), the agent jumps over the local energy barriers around the current local minimum. At the high-temperature stage (T = T h ), it transitions to a place far from the current state for further exploration. As a result, it is possible to crawl through a rugged landscape that contains many local minima and energy barriers and efficiently find a solution such that the cost function is as small as possible (Fig. 6). For the optimization of the design of a binary metamaterial multilayer, the cost function, C, is defined as follows: where T (f) is the transmission spectrum to be optimized, T target (f) is the required transmission, and f a and f b are the lower and upper frequencies, respectively (f a < f b ).
In the subsequent calculations, I aimed to determine the multilayer structure such that the value of the cost function, C, was minimized. First, I started with N randomly arranged A-and B-blocks, where the number of each block type was set to be arbitrary. Next, I attempted to update the state during each iteration based on a random choice using one of the following three different processes to generate a novel state in order to search for a much smaller cost function. The first method uses a process called single flip (Fig. 7a) and randomly selects one of the N blocks lined up in a row in the binary multilayer and changes it to a B-block if the corresponding block is an A-block, and vice versa, to set the trial state. The second method uses a process called multiple flips (Fig. 7b). From the N blocks that constitute the system, m blocks are selected randomly. These blocks do not have to be adjacent to each other. If each block is an A-block, it is changed to a B-block, and vice versa. The third method uses a process called cluster flips (Fig. 7c), where the i-th to the j-th consecutive blocks (1 i j N) from the left are considered to form a cluster in a set of N blocks. Here, i and j are chosen randomly. If each block  contained in this selected cluster is an A-block, it is changed to a B-block, and vice versa, so that (j À i + 1) blocks are changed at once, resulting in a new trial state. In this study, the total number of blocks constituting the binary metamaterial multilayer was set to N = 256. Thus, the total number of possible combinations of sequences composed of A-blocks and B-blocks would be 2 256 ≈ 1.16 Â 10 77 , and it would be impossible to examine all the systems in a strict sense. Hence, for the rest of the analysis, it was assumed that the A-block constituted the positive-refractive-index-material layer while the metamaterial consisted of the B-block. In addition, the thicknesses of these layers were d A = 12 mm and d B = 6 mm, respectively. With respect to the expressions for the permittivities of the A-blocks and B-blocks, namely, e A and e B (f), respectively, and their permeabilities, namely, m A and m B (f), respectively, the expressions used in the previous section (i.e., Eqs. (14), (15), (16), and (17)) were adopted. Finally, it was assumed that the electromagnetic waves were vertically incident. In the restart tempering method, the values of temperatures T l , T m , and T h were chosen to be 0.01, 0.05, and 0.30, respectively. The transmission rate, T (f), was calculated using the transfer matrix method [39].
The Rudin-Shapiro sequence [40] with N = 256 was employed to produce the target T target (f) of the aperiodic system to verify the effectiveness of the proposed method. Figure 8 shows an example of the calculation results. The dashed line represents the transmission spectrum of the Rudin-Shapiro system, and the solid line represents the result of the optimization performed using the restart tempering method. The interval, N t , at each stage was chosen to be 100, and the total number of trial steps was 2 Â 10 5 . In the Rudin-Shapiro system, the transmission rate becomes zero near the frequency at which the average refractive index becomes zero, and the profile of T (f) exhibits spikes in other frequency regimes. It can be seen that the optimized results for T (f) approximate the behavior of T target (f) with respect to its frequency dependence; the advantage of the proposed method is that trapping in the local minima during the search process is less likely. In addition, because it is specialized for optimization, it requires fewer computational resources than existing methods such as parallel tempering.  In conclusion, I numerically analyzed the effects of randomness in metamaterial composite systems, in which the metamaterials are embedded within a positiverefractive-index medium. Specifically, the band gap of a 2DMPC system was confirmed using the ADE-FDTD method. In comparison to the behavior of conventional photonic crystals, the band gap of the 2DMPC system was found to be less sensitive to the disorder of the system. With respect to the creation of photonic crystals, an essential question is how to form crystals with few defects. In contrast, in metamaterials, the properties of the system do not change significantly even when the metamaterials are placed randomly in the medium. To analyze the eigenmodes of such a system, that is, a 2DMPC with disorder, calculations based on the self-consistent FDFD method were performed. I confirmed the effectiveness of the implicitly restarted Arnoldi method for use in such calculations. Although this is not directly shown in Section 3, the method can also be applied in cases in which the imaginary part of permittivity and permeability of the metamaterial are nonzero (Im e (f) ≠ 0, Im m (f) ≠ 0) without any modification. I also proposed a method for the inverse design of random metamaterial multilayers. In particular, I proposed a novel method for which trapping in the local minima is less likely than for simulated annealing and that is less computationally demanding than parallel tempering. These methods should help increase the industrial applicability of random metamaterials in the future. In addition, this study provides valuable insights for the efficient design of metamaterial photonic crystals.