Issue 
EPJ Appl. Metamat.
Volume 9, 2022
Metamaterials for Novel Wave Phenomena in Microwaves, Optics, and Mechanics



Article Number  1  
Number of page(s)  10  
DOI  https://doi.org/10.1051/epjam/2021012  
Published online  01 February 2022 
https://doi.org/10.1051/epjam/2021012
Research Article
Numerical methods for design of metamaterial photonic crystals and random metamaterials
Department of Electrical, Electronic and Computer Engineering, Gifu University, Gifu 5011193, Japan
^{*} email: terao@gifuu.ac.jp
Received:
11
October
2021
Accepted:
20
December
2021
Published online: 1 February 2022
Twodimensional metamaterial photonic crystals (2DMPCs) composed of dispersive metamaterials in a positiverefractiveindex medium were investigated by incorporating finitedifference timedomain 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 selfconsistent finitedifference frequencydomain method, an eigenmode analysis of 2DMPCs with positional disorder was performed. Finally, a numerical method for the inverse design of binary random metamaterial multilayers was proposed.
Key words: Metamaterial photonic crystals / random metamaterials / inverse design
© T. Terao, Published by EDP Sciences, 2022
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Electromagnetic wave propagation in negativerefractiveindex materials embedded within positiverefractiveindex materials is attracting significant interest because it raises fundamental questions about the physical properties of electromagnetic waves [1–5]. In early studies on the topic, onedimensional multilayers consisting of alternating slabs of positive and negativerefractiveindex materials [6] and twodimensional metamaterial photonic crystals (2DMPCs) [7–9] were investigated. In a previous work [7], I investigated 2DMPCs composed of dispersive lefthanded materials, considering the Drudetype 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 lefthanded metamaterials embedded in the positiverefractiveindex 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 PhCbased devices. Numerical and experimental studies [10–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 twodimensional 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 finitedifference timedomain calculations into the auxiliary differential equation (ADEFDTD) 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, finitedifference frequencydomain (FDFD) method is considered to be a promising approach [14–17]. Using the selfconsistent 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 largescale 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 Twodimensional metamaterial photonic crystals (2DMPCs) with disorder
2.1 Analysis of metamaterial composites using ADEFDTD method
In this section, the modeling of 2DMPCs consisting of a periodic array of parallel metamaterial rods with circular crosssections in air is described. The Drude model for dispersive metamaterials, in which both the permittivity, ϵ (ω), and permeability, μ (ω), are frequency dependent, was employed:(1)and(2)where ϵ_{0} and μ_{0} are the permittivity and permeability, respectively, in vacuum; ω_{pe} and ω_{pm} are the electric and magnetic plasma frequencies, respectively; and γ_{e} and γ_{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 twodimensional 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, δ, i.e., disorder in the positions of the metamaterial rods, and the other is size disorder, Δ, which is related to the radius of the rods (Fig. 2). In the case of the positional disorder, the random deviations (δ_{i,x} and δ_{i,y}) of the ith rod from the corresponding lattice position in the xy plane exhibit the relationship and are uniformly distributed (Fig. 2a). With respect to the size disorder, the radius, r_{i}, of the ith rod exhibits a uniform distribution within the range r − Δ < r_{i} < r + Δ (Fig. 2b).
The ADEFDTD 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:(3)
Similarly, the relationship between the magnetic flux density, B (r, t), and the magnetic field, H (r, t), is given by(4)For a given wavevector, q, Bloch's periodic boundary conditions were considered such that(5) where R is an arbitrary lattice vector in the unit cell. In this study, I considered the power spectrum, P (ω), which corresponds to the square of the photonic DOS. At t ≈ t_{0}, an impulse is imposed on the system such that(6)where is the zcomponent 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:(7)
The power spectrum, P (ω) , of the system can be calculated from(8)
Here, ⟨ ⋯ ⟩_{ {q}} denotes the average of the set of wavevectors, q, in the Brillouin zone. The peak position of the frequency spectrum of the timecorrelation function is obtained for each wavevector, q, and the power spectrum, P (ω), 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(9)
In general, the value of the power spectrum, P (ω), 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 (δ ≠ 0 or Δ ≠ 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 (ω) , for TM polarization in the squarelattice 2DMPCs, for which the square lattice constant is a. Furthermore, the radius of the rods, r, was 0.17a, and the positional disorder was δ. 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 ω_{pe} = ω_{pm} = 1.0 (c/a), and low loss parameters were assumed by setting γ_{e} = γ_{m} = 10^{−5} (c/a). The mesh resolution for the ADEFDTD calculations was taken to be 0.01a.
Figure 3b shows the calculation results for P (ω) for a triangular lattice with positional disorder δ, 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 (ω)for a square lattice with size disorder Δ, 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.
Fig. 1 Schematic illustration of 2DMPC with array of circular rods. 
Fig. 2 Disordered 2DMPC formed by supercell method: (a) positional disorder and (b) size disorder. 
Fig. 3 Power spectrum, P (ω), for 2DMPC composed of circular rods: (a) square lattice with positional disorder, (b) triangular lattice with positional disorder, and (c) triangular lattice with size disorder. 
2.2 Eigenmode analysis using selfconsistent 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 permittivity (permeability) follows a specific functional form such as the Debye, Lorentz, or Drude dispersion relation, and the piecewise linear recursive convolution and ADEFDTD 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 twodimensional periodicity is now described in brief. I begin with the Maxwell equations:(10)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,(11)in the case of TE polarization, where ϵ (r, ω) and μ (r, ω) are the relative permittivity and relative permeability, respectively, at position r, and angular frequency ω (= 2πf), and c is the speed of light in vacuum. At this stage, ϵ (r, ω) and μ (r, ω) 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(12) (13)where ϵ_{0} is the dielectric constant in a vacuum.
In this case, the eigenvalues and matrix elements for the eigenvalue problem need to be solved selfconsistently 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 selfconsistent 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 nonHermitian) 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–23]. In this study, the implicitly restarted Arnoldi method with the shiftandinvert technique was adopted [24–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 ϵ_{A} and μ_{A}, respectively, and the permittivity and permeability of the metamaterial were ϵ_{B} (f) and μ_{B} (f), respectively, such that(14) (15) (16) (17)
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 (δ = 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 selfconsistent FDFD method, the FDFD calculation is performed considering the permittivity and permeability values of the metamaterial at multiple frequencies 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 , 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, (GHz) was plotted, and the dotted line in the figure indicates the condition .
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 selfconsistent FDFD methods.
Fig. 4 (a) Supercell containing metamaterial rods arranged in square lattice with positional disorder and eigenmodes with eigenfrequencies. (b) Frequency analysis by selfconsistent FDFD method (c) f = 2.97998 GHz, and (d) f = 2.99896 GHz. 
3 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–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 positiverefractiveindex material (Ablock) and a layer of the metamaterial (Bblock), 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 Bblocks can be chosen freely, with the total number of blocks, N, being fixed. ϵ_{A} and μ_{A} in block A and ϵ_{B} (f) and μ_{B} (f) are given by equations (14), (15), (16), and (17).
In this study, I focused on the similarity between the abovedescribed problem and the Ising model in statistical physics. In these systems, the Ablocks correspond to the upward spins of the Ising model and the Bblocks 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 lowenergy 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 lowtemperature stage (T = T_{l}), the agent searches for a local minimum. At the mediumtemperature stage (T = T_{m}), the agent jumps over the local energy barriers around the current local minimum. At the hightemperature 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:(18)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 Bblocks, 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 Bblock if the corresponding block is an Ablock, 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 Ablock, it is changed to a Bblock, and vice versa. The third method uses a process called cluster flips (Fig. 7c), where the ith to the jth 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 Ablock, it is changed to a Bblock, 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 Ablocks and Bblocks 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 Ablock constituted the positiverefractiveindexmaterial layer while the metamaterial consisted of the Bblock. 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 Ablocks and Bblocks, namely, ϵ_{A} and ϵ_{B} (f), respectively, and their permeabilities, namely, μ_{A} and μ_{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.
Fig. 5 Optimized design of metamaterial multilayers. 
Fig. 6 Restart tempering method in rugged landscape. 
Fig. 7 Details of state update processes used in restart tempering method: (a) single flip, (b) multiple flips, and (c) cluster flip. 
Fig. 8 Comparison of optimized transmission rate, T (f) , and target function, T_{target} (f). 
4 Conclusions
In conclusion, I numerically analyzed the effects of randomness in metamaterial composite systems, in which the metamaterials are embedded within a positiverefractiveindex medium. Specifically, the band gap of a 2DMPC system was confirmed using the ADEFDTD 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 selfconsistent 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 ϵ (f) ≠ 0, Im μ (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.
Acknowledgments
The author acknowledges the Supercomputer Center at the University of Tokyo, for the use of their facilities. This work was supported by JSPS KAKENHI Grant Number 20K11848.
References
 V.G. Veselago, The electrodynamics of substances with simultaneously negative values of ε and μ, Sov. Phys. Usp. 10, 509 (1968) [CrossRef] [Google Scholar]
 S.A. Ramakrishnan, Physics of negative refractive index materials, Rep. Prog. Phys. 68, 449 (2005) [CrossRef] [Google Scholar]
 R.A. Shelby, D.R. Smith, S. Schultz, Experimental verification of a negative index of refraction, Science 292, 77 (2001) [CrossRef] [Google Scholar]
 V.M. Shalaev, Optical negativeindex metamaterials, Nat. Photonics 1, 41 (2007) [CrossRef] [Google Scholar]
 J.B. Pendry, Negative refraction makes a perfect lens, Phys. Rev. Lett. 85, 3966 (2000) [Google Scholar]
 J. Li, L. Zhou, C.T. Chan, P. Sheng, Photonic band gap from a stack of positive and negative index materials, Phys. Rev. Lett. 90, 083901 (2003) [CrossRef] [Google Scholar]
 T. Terao, Numerical study of twodimensional metamaterial composites: the existence of a stop band, J. Mod. Opt. 60, 1997 (2013) [CrossRef] [Google Scholar]
 T. Terao, Finitedifference frequencydomain method to study photonic crystals with dispersive metamaterials, Jpn. J. Appl. Phys. 50, 102204 (2011) [CrossRef] [Google Scholar]
 A. Soltani Vala, S. Roshan Entezar, A.A. Sedghi, Band structure of twodimensional square lattice photonic crystals of circular dispersive metamaterial rods, Eur. Phys. J. B 81, 269 (2011) [CrossRef] [Google Scholar]
 K. Edagawa, S. Kanoko, M. Notomi, Photonic amorphous diamond Structure with a 3D photonic band gap, Phys. Rev. Lett. 100, 013901 (2008) [CrossRef] [Google Scholar]
 M. Florescu, S. Torquato, P.J. Steinhardt, Designer disordered materials with large, complete photonic band gaps, Proc. Natl. Acad. Sci. U.S.A. 106, 20658 (2009) [CrossRef] [Google Scholar]
 W. Man, M. Florescu, K. Matsuyama, P. Yadak, G. Nahal, S. Hashemizad, E. Williamson, P. Steinhardt, S. Torquato, P. Chaikin, Photonic band gap in isotropic hyperuniform disordered solids with low dielectric contrast, Opt. Express 21, 19972 (2013) [CrossRef] [Google Scholar]
 R. Meisels, F. Kuchar, Densityofstates and wave propagation in twodimensional photonic crystals with positional disorder, J. Opt. A. 9, S396 (2007) [CrossRef] [Google Scholar]
 S. Guo, F. Wu, S. Albin, R. Rogowski, Photonic band gap analysis using finitedifference frequencydomain method, Opt. Express 12, 1741 (2004) [CrossRef] [Google Scholar]
 A.M. Ivinskaya, A.V. Lavrinenko, D.M. Shyroki, Modeling of nanophotonic resonators with the finitedifference frequencydomain method, IEEE Trans. Antennas Propag. 59, 4155 (2011) [CrossRef] [Google Scholar]
 D.M. Shyroki, A.V. Lavrinenko, Perfectly matched layer method in the finitedifference timedomain and frequencydomain calculations, Phys. Stat. Sol. 244, 3506 (2007) [CrossRef] [Google Scholar]
 T. Terao, Computing interior eigenvalues of nonsymmetric matrices: application to threedimensional metamaterial composites, Phys. Rev. E 82, 026702 (2010) [CrossRef] [Google Scholar]
 A. Taflove, S.C. Hagness, Computational Electrodynamics: The FiniteDifference TimeDomain Method, 3rd ed. (Artech House, Boston, 2005) [Google Scholar]
 T. Terao, Numerical analysis of metamaterial photonic crystals with positional disorder, in Proceedings of 15th international congress on artificial materials for novel wave phenomena (Metamaterials 2021), in press [Google Scholar]
 T. Terao, Manipulating sonic band gaps at will: vibrational density of states in threedimensional acoustic metamaterial composites, Waves Random Complex Media 28, 253 (2018) [CrossRef] [Google Scholar]
 G.H. Golub, C.F. van Loan, Matrix Computations, 4th ed. (John Hopkins, Baltimore, 2013) [Google Scholar]
 Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed. (SIAM, Philadelphia, 2003) [Google Scholar]
 F. Chatelin, Eigenvalues of Matrices (SIAM, Philadelphia, 2012) [CrossRef] [Google Scholar]
 R.B. Lehoucq, D.C. Sorensen, C. Yang, ARPACK Users' Guide: Solution of LargeScale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods (SIAM, Philadelphia, 1998) [CrossRef] [Google Scholar]
 U. Elsner, V. Mehrmann, F. Milde, R.A. Römer, M. Schreiber, SIAM J. Sci. Comput. 20, 2089 (1999) [CrossRef] [Google Scholar]
 T. Terao, Eigenvalue analysis of the threedimensional tightbinding model with nonHermitian disorder, Phys. Rev. B 103, 224201 (2021) [CrossRef] [Google Scholar]
 H. Daninthe, S. Foteinopoulou, C.M. Soukoulis, Omnireflectance and enhanced resonant tunneling from multilayers containing lefthanded materials, Photonics Nanostruct. Fundam. Appl. 4, 123 (2006) [CrossRef] [Google Scholar]
 A.A. Asatryan, L.C. Botten, M.A. Byrne, V.D. Freilikher, S.A. Gredeskul, I.V. Shadrivov, R.C. McPhedran, Y.S. Kivshar, Suppression of anderson localization in disordered metamaterial, Phys. Rev. Lett. 99, 193902 (2007) [CrossRef] [Google Scholar]
 J.A. Monsoriu, R.A. Depine, E.J. Silvestre, NonBragg band gaps in 1D metamaterial aperiodic multilayers, Eur. Opt. Soc. Rapid Publ. 2, 07002 (2007) [CrossRef] [Google Scholar]
 E. ReyesGómez, N. Raigoza, S.B. Cavalcanti, C.A.A. de Carvalho, L.E. Oliveira, Absorption effects on plasmon polaritons in quasiperiodic photonic superlattices containing a metamaterial, J. Phys.: Condens. Matter 22, 385901 (2010) [CrossRef] [Google Scholar]
 J.A. Monsoriu, R.A. Depine, M.L. MartínezRicci, E. Silvestre, P. Andrés, mbonacci metamaterial multilayers: location of the zeroaverage index bandgap edges, Opt. Lett. 34, 3172 (2009) [CrossRef] [Google Scholar]
 A.A. Asatryan, L.C. Botten, M.A. Bayne, V.D. Freilikher, S.A. Gredeskul, I.V. Shadrivov, R.C. McPhedran, Y.S. Kivshar, Transmission and Anderson localization in dispersive metamaterials, Phys. Rev. B 85, 045122 (2012) [CrossRef] [Google Scholar]
 D. Mogilevtsev, F.A. Pinheiro, R.R. dos Santos, S.B. Cavalcanti, S.B.L.E. Oliveira, Suppression of Anderson localization of light and Brewster anomalies in disordered superlattices containing a dispersive metamaterial, Phys. Rev. B 82, 081105 (2010) [CrossRef] [Google Scholar]
 E. Marinari, G. Parisi, Simulated tempering: a new Monte Carlo scheme, Europhys. Lett. 19, 451 (1992) [CrossRef] [Google Scholar]
 A.P. Lyubartsev, A.A. Martinovski, S.V. Shevkunov, P.N. VorontsovVelyaminov, New approach to Monte Carlo calculation of the free energy: method of expanded ensembles, J. Chem. Phys. 96, 1776 (1992) [CrossRef] [Google Scholar]
 C.J. Geyer, Markov chain Monte Carlo maximum likelihood, in E.M. Keramidas (Ed.), Computing Science and Statistics: Proceedings of the 23rd Symposium on the Interface (American Statistical Association, New York 1991), p. 156 [Google Scholar]
 K. Hukushima, K. Nemoto, Exchange Monte Carlo method and application to spin glass simulations, J. Phys. Soc. Jpn. 65, 1604 (1996) [CrossRef] [Google Scholar]
 S. Trebst, M. Troyer, U.H.E. Hansmann, Optimized parallel tempering simulations of proteins, J. Chem. Phys. 124, 174903 (2006) [CrossRef] [Google Scholar]
 M. Born, E. Wolf, with contributions by A.B. Bhatia et al., Principles of optics: electromagnetic theory of propagation, interference and diffraction of light 7th ed. (Cambridge University Press, Cambridge, 1999) [CrossRef] [Google Scholar]
 W. Salejda, A. KlauzerKruszyna, K. Tarnowski, Photonic band structure of onedimensional aperiodic superlattices composed of negative refraction metamaterials, Proc. SPIE 6581, 658112 (2007) [CrossRef] [Google Scholar]
Cite this article as: Takamichi Terao, Numerical methods for design of metamaterial photonic crystals and random metamaterials, EPJ Appl. Metamat. 9, 1 (2022)
All Figures
Fig. 1 Schematic illustration of 2DMPC with array of circular rods. 

In the text 
Fig. 2 Disordered 2DMPC formed by supercell method: (a) positional disorder and (b) size disorder. 

In the text 
Fig. 3 Power spectrum, P (ω), for 2DMPC composed of circular rods: (a) square lattice with positional disorder, (b) triangular lattice with positional disorder, and (c) triangular lattice with size disorder. 

In the text 
Fig. 4 (a) Supercell containing metamaterial rods arranged in square lattice with positional disorder and eigenmodes with eigenfrequencies. (b) Frequency analysis by selfconsistent FDFD method (c) f = 2.97998 GHz, and (d) f = 2.99896 GHz. 

In the text 
Fig. 5 Optimized design of metamaterial multilayers. 

In the text 
Fig. 6 Restart tempering method in rugged landscape. 

In the text 
Fig. 7 Details of state update processes used in restart tempering method: (a) single flip, (b) multiple flips, and (c) cluster flip. 

In the text 
Fig. 8 Comparison of optimized transmission rate, T (f) , and target function, T_{target} (f). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.