EPJ Appl. Metamat.
Volume 9, 2022
Metamaterials for Novel Wave Phenomena in Microwaves, Optics, and Mechanics
Article Number 1
Number of page(s) 10
Published online 01 February 2022

© T. Terao, Published by EDP Sciences, 2022

Licence Creative CommonsThis is an Open Access article distributed under the terms of the Creative Commons Attribution License (, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

Electromagnetic wave propagation in negative-refractive-index materials embedded within positive-refractive-index materials is attracting significant interest because it raises fundamental questions about the physical properties of electromagnetic waves [15]. 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) [79] 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 [1012] 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 time-domain 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 [1417]. 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, ϵ (ω), 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 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, δ, 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 i-th rod from the corresponding lattice position in the x-y plane exhibit the relationship and are uniformly distributed (Fig. 2a). With respect to the size disorder, the radius, ri, of the i-th rod exhibits a uniform distribution within the range r − Δ < ri < r + Δ (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:(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 ≈ t0, an impulse is imposed on the system such that(6)where is the z-component of the electric field imposed on a randomly chosen cell, i0, 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 (ts >> t0 + tp), 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 time-correlation 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, i0, 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 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 δ. The supercell system size was fixed at Lx = Lx = 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 ADE-FDTD 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 Lx = 4.0a and Ly = 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 Lx = Lx = 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.

thumbnail Fig. 1

Schematic illustration of 2DMPC with array of circular rods.

thumbnail Fig. 2

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

thumbnail 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 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 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:(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 = (Ex, Ey, 0) and H = (0,  0,  Hz) polarization. When the current density, j, is equal to zero, the amplitude, Hz (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, Ex (r) and Ey (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 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 [2123]. In this study, the implicitly restarted Arnoldi method with the shift-and-invert technique was adopted [2426]. I used the self-consistent 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−2m). 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 103 × 103, and the dimension of the matrix to be diagonalized was 106.

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 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 Hz (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.

thumbnail Fig. 4

(a) Supercell containing metamaterial rods arranged in square lattice with positional disorder and eigenmodes with eigenfrequencies. (b) Frequency analysis by self-consistent 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,2733]. 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 positive-refractive-index material (A-block) and a layer of the metamaterial (B-block), I designed a multilayer structure that almost satisfies the required transmission spectrum, Ttarget (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. ϵ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 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 high-energy 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, Tl, Tm, and Th (Tl < Tm <Th), are considered, and the temperature of the system is changed as ThTmTlTh → ⋯, with the number of interval steps at each stage, Nt, being fixed. At the low-temperature stage (T = Tl), the agent searches for a local minimum. At the medium-temperature stage (T = Tm), the agent jumps over the local energy barriers around the current local minimum. At the high-temperature stage (T = Th), 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, Ttarget (f) is the required transmission, and fa and fb are the lower and upper frequencies, respectively (fa < fb). 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 (ji + 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 2256 ≈ 1.16 × 1077, 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 dA = 12 mm and dB = 6 mm, respectively. With respect to the expressions for the permittivities of the A-blocks and B-blocks, 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 Tl, Tm, and Th 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 Ttarget (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, Nt, at each stage was chosen to be 100, and the total number of trial steps was 2 × 105. 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 Ttarget (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.

thumbnail Fig. 5

Optimized design of metamaterial multilayers.

thumbnail Fig. 6

Restart tempering method in rugged landscape.

thumbnail Fig. 7

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

thumbnail Fig. 8

Comparison of optimized transmission rate, T (f) , and target function, Ttarget (f).

4 Conclusions

In conclusion, I numerically analyzed the effects of randomness in metamaterial composite systems, in which the metamaterials are embedded within a positive-refractive-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 ϵ (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.


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.


  1. V.G. Veselago, The electrodynamics of substances with simultaneously negative values of ε and μ, Sov. Phys. Usp. 10, 509 (1968) [CrossRef] [Google Scholar]
  2. S.A. Ramakrishnan, Physics of negative refractive index materials, Rep. Prog. Phys. 68, 449 (2005) [CrossRef] [Google Scholar]
  3. R.A. Shelby, D.R. Smith, S. Schultz, Experimental verification of a negative index of refraction, Science 292, 77 (2001) [CrossRef] [Google Scholar]
  4. V.M. Shalaev, Optical negative-index metamaterials, Nat. Photonics 1, 41 (2007) [CrossRef] [Google Scholar]
  5. J.B. Pendry, Negative refraction makes a perfect lens, Phys. Rev. Lett. 85, 3966 (2000) [CrossRef] [Google Scholar]
  6. 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]
  7. T. Terao, Numerical study of two-dimensional metamaterial composites: the existence of a stop band, J. Mod. Opt. 60, 1997 (2013) [CrossRef] [Google Scholar]
  8. T. Terao, Finite-difference frequency-domain method to study photonic crystals with dispersive metamaterials, Jpn. J. Appl. Phys. 50, 102204 (2011) [CrossRef] [Google Scholar]
  9. A. Soltani Vala, S. Roshan Entezar, A.A. Sedghi, Band structure of two-dimensional square lattice photonic crystals of circular dispersive metamaterial rods, Eur. Phys. J. B 81, 269 (2011) [CrossRef] [Google Scholar]
  10. 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]
  11. 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]
  12. 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]
  13. R. Meisels, F. Kuchar, Density-of-states and wave propagation in two-dimensional photonic crystals with positional disorder, J. Opt. A. 9, S396 (2007) [CrossRef] [Google Scholar]
  14. S. Guo, F. Wu, S. Albin, R. Rogowski, Photonic band gap analysis using finite-difference frequency-domain method, Opt. Express 12, 1741 (2004) [CrossRef] [Google Scholar]
  15. A.M. Ivinskaya, A.V. Lavrinenko, D.M. Shyroki, Modeling of nanophotonic resonators with the finite-difference frequency-domain method, IEEE Trans. Antennas Propag. 59, 4155 (2011) [CrossRef] [Google Scholar]
  16. D.M. Shyroki, A.V. Lavrinenko, Perfectly matched layer method in the finite-difference time-domain and frequency-domain calculations, Phys. Stat. Sol. 244, 3506 (2007) [CrossRef] [Google Scholar]
  17. T. Terao, Computing interior eigenvalues of nonsymmetric matrices: application to three-dimensional metamaterial composites, Phys. Rev. E 82, 026702 (2010) [CrossRef] [Google Scholar]
  18. A. Taflove, S.C. Hagness, Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed. (Artech House, Boston, 2005) [Google Scholar]
  19. 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]
  20. T. Terao, Manipulating sonic band gaps at will: vibrational density of states in three-dimensional acoustic metamaterial composites, Waves Random Complex Media 28, 253 (2018) [CrossRef] [Google Scholar]
  21. G.H. Golub, C.F. van Loan, Matrix Computations, 4th ed. (John Hopkins, Baltimore, 2013) [Google Scholar]
  22. Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd ed. (SIAM, Philadelphia, 2003) [Google Scholar]
  23. F. Chatelin, Eigenvalues of Matrices (SIAM, Philadelphia, 2012) [CrossRef] [Google Scholar]
  24. R.B. Lehoucq, D.C. Sorensen, C. Yang, ARPACK Users' Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods (SIAM, Philadelphia, 1998) [CrossRef] [Google Scholar]
  25. U. Elsner, V. Mehrmann, F. Milde, R.A. Römer, M. Schreiber, SIAM J. Sci. Comput. 20, 2089 (1999) [CrossRef] [Google Scholar]
  26. T. Terao, Eigenvalue analysis of the three-dimensional tight-binding model with non-Hermitian disorder, Phys. Rev. B 103, 224201 (2021) [CrossRef] [Google Scholar]
  27. H. Daninthe, S. Foteinopoulou, C.M. Soukoulis, Omni-reflectance and enhanced resonant tunneling from multilayers containing left-handed materials, Photonics Nanostruct. Fundam. Appl. 4, 123 (2006) [CrossRef] [Google Scholar]
  28. 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]
  29. J.A. Monsoriu, R.A. Depine, E.J. Silvestre, Non-Bragg band gaps in 1D metamaterial aperiodic multilayers, Eur. Opt. Soc. Rapid Publ. 2, 07002 (2007) [CrossRef] [Google Scholar]
  30. E. Reyes-Gó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]
  31. J.A. Monsoriu, R.A. Depine, M.L. Martínez-Ricci, E. Silvestre, P. Andrés, m-bonacci metamaterial multilayers: location of the zero-average index bandgap edges, Opt. Lett. 34, 3172 (2009) [CrossRef] [Google Scholar]
  32. 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]
  33. 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]
  34. E. Marinari, G. Parisi, Simulated tempering: a new Monte Carlo scheme, Europhys. Lett. 19, 451 (1992) [CrossRef] [Google Scholar]
  35. A.P. Lyubartsev, A.A. Martinovski, S.V. Shevkunov, P.N. Vorontsov-Velyaminov, New approach to Monte Carlo calculation of the free energy: method of expanded ensembles, J. Chem. Phys. 96, 1776 (1992) [CrossRef] [Google Scholar]
  36. 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]
  37. K. Hukushima, K. Nemoto, Exchange Monte Carlo method and application to spin glass simulations, J. Phys. Soc. Jpn. 65, 1604 (1996) [CrossRef] [Google Scholar]
  38. S. Trebst, M. Troyer, U.H.E. Hansmann, Optimized parallel tempering simulations of proteins, J. Chem. Phys. 124, 174903 (2006) [CrossRef] [Google Scholar]
  39. 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]
  40. W. Salejda, A. Klauzer-Kruszyna, K. Tarnowski, Photonic band structure of one-dimensional 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

thumbnail Fig. 1

Schematic illustration of 2DMPC with array of circular rods.

In the text
thumbnail Fig. 2

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

In the text
thumbnail 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
thumbnail Fig. 4

(a) Supercell containing metamaterial rods arranged in square lattice with positional disorder and eigenmodes with eigenfrequencies. (b) Frequency analysis by self-consistent FDFD method (c) f = 2.97998 GHz, and (d) f = 2.99896 GHz.

In the text
thumbnail Fig. 5

Optimized design of metamaterial multilayers.

In the text
thumbnail Fig. 6

Restart tempering method in rugged landscape.

In the text
thumbnail 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
thumbnail Fig. 8

Comparison of optimized transmission rate, T (f) , and target function, Ttarget (f).

In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.