Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Photonic crystal topological design for polarized and polarization-independent band gaps by gradient-free topology optimization

Open Access Open Access

Abstract

Photonic crystals can be adopted to control light propagation due to their superior band gap feature. It is well known the band gap feature of photonic crystals depends significantly on the topological design of the lattices, which is rather challenging due to the highly nonlinear objective function and multiple local minima feature of such design problems. To this end, this paper proposed a new band-gap topology optimization framework for photonic crystals considering different electromagnetic wave polarization modes. Based on the material-field series-expansion (MFSE) model and the dielectric permittivity interpolation scheme, the lattice topologies are represented by using a small number of design variables. Then, a sequential Kriging-based optimization algorithm, which shows strong global search capability and requires no sensitivity information, is employed to solve the band gap design problem as a series of sub-optimization problems with adaptive-adjusting design spaces. Numerical examples demonstrated the effectiveness of the proposed gradient-free method to maximize the band gap for transverse magnetic field (TM), transverse electric field (TE), and complete modes. Compared with previously reported designs, the present results exhibit less dependency on the guess of the initial design, larger band gaps and some interesting topology configurations.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

Photonic crystals (PhCs) are optical structures constructed by the periodic arrangement of dielectric materials with different permittivities. PhCs (which are also called photonic band gap structures) prevent light with specified frequencies (i.e., a certain range of wavelengths) from propagating in certain directions [1,2], and are widely adopted in integrated optics [36]. Generally, a larger band gap means a larger design application space and better performance. So it is of great significance to design the PhCs with prescribed or maximum band gaps from a practical perspective.

Due to the periodicity of the PhCs, their optical properties can be evaluated on a unit cell (the minimum repetitive unit). For a given dielectric material, the band gaps of PhCs depend on the spatial distributions of constituents, which makes PhCs highly designable [7]. In many earlier studies, the design of PhCs is based on geometric analysis, extensive parametric studies, and even physical intuitions, which turns out to be inefficient and time-consuming [812]. Furthermore, the PhCs obtained with traditional trial and error methods may be far from optimal, and may be difficult to achieve large band gaps between two designated adjacent energy bands.

Topology optimization provides an effective way to find the material distribution with the best structural performances without prior knowledge of optimal topology [1315]. In the past two decades, the band gap designs of photonic crystals on the basis of various topology optimization methods have been reported [1627]. Moreover, the applications of topology optimization for photonics derived from photonic crystals have attracted much attention, such as photonic microcavities [28,29], photonic demultiplexers [30], metasurfaces [31,32], photonic topological insulators [33,34], waveguides [35], etc. Kao et al. [16] employed the level set method to design 2D photonic crystals with band gaps. The interface between two materials with different dielectric permittivities was moved by a generalized gradient ascent method. Men et al. [20,21] used the convex semi-definite program and the subspace method to optimize the 2D PhCs. In their study, the optimization starts with randomly generated initial designs. Although the above-mentioned methods can realize the band gap design of PhCs, the obtained solutions are highly dependent on the initial guess due to the multiple local minima feature of the band gap optimization problem.

It was well recognized that, when using gradient-based optimization algorithms, the highly nonlinear feature of objective functions in PhCs topology optimization problems may lead to low success rates and local solutions. In order to relieve this difficulty, Sigmund and Jensen [36] proposed a two-stage topology optimization strategy for PhCs, in which the first stage adopted the genetic algorithm to find an initial design on the coarse finite element mesh, and the second stage used the density-based method with a fine grid to enlarge the photonic band gaps. In the design of complete photonic band gaps, Meng et al. [37] used the bi-directional evolutionary optimization method to obtain novel topologies based on several heuristic initial designs, which were created by superposing two specific photonic structures with independent transverse magnetic field (TM) and transverse electric (TE) band gaps. However, the issue of initial solution dependency is far from being overcome due to the fact that a suitable initial design is difficult to construct.

Recently, an effective gradient-free topology optimization approach that combines the material-field series-expansion (MFSE) model and the sequential Kriging-based optimization algorithm has been proposed by Luo et al. [40]. This method has been successfully applied to microstructure designs with desired phononic and mechanical performances [41,42]. In this study, we proposed a new topology optimization framework based on the MFSE model to realize the maximum band gap design of 2D PhCs considering transverse magnetic field (TM), transverse electric field (TE), or both of them. With the MFSE method, the lattice topologies can be sufficiently represented with a small number of design variables and independent of the finite element mesh. This facilitates the adoption of a fine finite element mesh to ensure the accuracy of the finite element analysis. Then, the sequential Kriging-based optimization algorithm is employed to solve the optimization problem as a series of sub-optimization problems. Compared with the gradient-based topology optimization methods, the proposed method does not require sensitivity information. Moreover, it exhibits relatively strong global search capability, which is beneficial to design problems of photonic crystals that possess multiple local minima and have significant dependence on the initial designs.

The remainder of this study is organized as follows. Section 2 introduces the numerical simulation of PhCs with the finite element method. Section 3 proposes the general topology optimization formulation for maximizing the band gap of PhCs considering different polarization modes. Section 4 addresses the MFSE-based topology representation of the square lattice and the sequential Kriging-based optimization scheme. Section 5 gives numerical examples that demonstrate the effectiveness of the proposed method. Conclusions are drawn in Section 6.

2. Finite element evaluation of photonic crystals

Several numerical analysis methods for PhCs have been proposed, which include the plane-wave expansion method [43,44], the finite-difference time-domain method [45,46], the transfer matrix method [47,48], the finite element (FE) method [49,50], etc. In this paper, we adopt the finite element method for the performance evaluation of PhCs.

The propagation of light in a photonic crystal is governed by the macroscopic Maxwell equations coupled in time and space [1]. In the absence of source, and by assuming that electromagnetic waves are monochromatic waves, the time-harmonic Maxwell equations can be decoupled into the following system by using Fourier analysis:

$$\nabla \times \nabla \times {\textbf E}({\textbf r} )= \varepsilon ({\textbf r} ){\left( {\frac{\omega }{c}} \right)^2}{\textbf E}({\textbf r} )$$
$$\nabla \times \left( {\frac{1}{{\varepsilon ({\textbf r} )}}\nabla \times {\textbf H}({\textbf r} )} \right) = {\left( {\frac{\omega }{c}} \right)^2}{\textbf H}({\textbf r} )\textrm{ }$$
where c is the speed of light, $\varepsilon ({\bf r} )$ is the dielectric function and $\omega$ is the angular frequency of the electromagnetic wave. In general, the analysis of the 2D PhCs can be divided into TM mode and TE mode [1] according to the polarization of electromagnetic field. In the TM mode, the magnetic field is confined to the plane of wave propagation and the electric field ${\bf E}$ is perpendicular to this plane. The electric field ${\bf E}$ can be expressed as $({0,0,\textrm{E(}{\bf r}\textrm{)}} )$, where the vector ${\textbf r} = ({x,\textrm{ }y} )$ denotes the spatial coordinates. In the TE mode, the magnetic field ${\bf H}$ is perpendicular to the plane and thus can be expressed as $({0,0,\textrm{H(}{\bf r}\textrm{)}} )$. Therefore, the Maxwell equations in Eqs. (1) and (2) can be further simplified to two scalar eigenvalue equations as follows
$$- \nabla \cdot ({\nabla \textrm{E}({\textbf r} )} )= \varepsilon ({\textbf r} ){\left( {\frac{\omega }{c}} \right)^2}\textrm{E}({\textbf r} ),\;\;\;\;\;\textrm{for TM mode}$$
$$- \nabla \cdot \left( {\frac{1}{{\varepsilon ({\textbf r} )}}\nabla \textrm{H}({\textbf r} )} \right) = {\left( {\frac{\omega }{c}} \right)^2}\textrm{H}({\textbf r} ),\;\;\;\;\;\textrm{for TE mode}\textrm{ }$$

For 2D PhCs, the dielectric function satisfies $\varepsilon ({\textbf r}) = \varepsilon ({\textbf r} + {\textbf R})$, where ${\textbf R}$ is the lattice translation vector. According to the Bloch–Floquet theory [51], the electromagnetic field in periodic structures can be represented as the product of a periodic function and an exponential function as ${\textbf E}({\textbf r} )= {e^{i{\textbf k} \cdot {\textbf r}}}{{\textbf E}_{\textbf k}}({\textbf r} )\textrm{ }$ and ${\textbf H}({\textbf r} )= {e^{i{\textbf k} \cdot {\textbf r}}}{{\textbf H}_{\textbf k}}({\textbf r} )\textrm{ }$, where ${{\textbf E}_{\textbf k}}({\textbf r} )$ and ${{\textbf H}_{\textbf k}}({\textbf r} )$ are the cell periodic fields satisfying

$$- ({\nabla + i{\textbf k}} )\cdot ({({\nabla + i{\textbf k}} ){{\textbf E}_{\textbf k}}({\textbf r} )} )= \varepsilon ({\textbf r} ){\left( {\frac{\omega }{c}} \right)^2}{{\textbf E}_{\textbf k}}({\textbf r} )\;\;\;\;\;\textrm{ for TM mode}, $$
and
$$- ({\nabla + i{\textbf k}} )\cdot \left( {\frac{1}{{\varepsilon ({\textbf r} )}}({\nabla + i{\textbf k}} ){{\textbf H}_{\textbf k}}({\textbf r} )} \right) = {\left( {\frac{\omega }{c}} \right)^2}{{\textbf H}_{\textbf k}}({\textbf r} )\;\;\;\;\;\textrm{ for TE mode}, $$
respectively. Here, ${\textbf k}$ is a wave vector in the first Brillouin zone. It should be noted that the unit cell, the first Brillouin zone and the lattice translation vector all depend on the lattice type. For a PhC with the square lattice as shown in Fig. 1, the lattice translation vector ${\textbf R} = ({a,a} )$, where a is the lattice constant. Due to the symmetry and periodicity of PhCs, we only need to consider all possible wave vectors ${\textbf k}$ on the irreducible Brillouin zone, especially its boundary [1].

 figure: Fig. 1.

Fig. 1. Left: The periodic lattice with the square lattice and lattice constant a. The white colour denotes air (${\varepsilon _1} = 1$), the red denotes GaAs (${\varepsilon _\textrm{2}} = 1\textrm{1}\textrm{.4}$), and the dashed box represents the unit cell; Right: the dashed box represents the first Brillouin zone in the reciprocal lattice. The irreducible zone is the gray triangular and its boundary is denoted by $\mathrm{\Gamma } \to \textrm{X} \to \textrm{M} \to \mathrm{\Gamma }$.

Download Full Size | PDF

For the given wave vector ${\textbf k}$ and dielectric function $\varepsilon ({\textbf r})$, one can use the finite element method to solve Eqs. (5) and (6), which can be written as a unified eigenvalue equation

$$\left[ {{\textbf K} - \frac{{{\omega^2}}}{{{c^2}}}{\textbf M}\textrm{ }} \right]{\textbf u} = 0$$
where ${\textbf u} = {{\textbf E}_{\textbf k}}({\textbf r} )$, ${\textbf K}({\textbf k} )={-} ({\nabla + i{\textbf k}} )\cdot ({\nabla + i\,{\textbf k}} ),\textrm{ }{\textbf M}(\varepsilon )= \varepsilon ({\textbf r} ){\textbf I}$ for the TM mode; and ${\textbf u} = {{\textbf H}_{\textbf k}}({\textbf r} )$, ${\textbf K}({{\textbf k},\varepsilon } )={-} ({\nabla + i{\textbf k}} )\cdot \left( {\frac{1}{{\varepsilon ({\textbf r} )}}({\nabla + i\,{\textbf k}} )} \right),\textrm{ }{\textbf M}\textrm{ = }{\textbf I}$ for the TE mode. Here ${\textbf I}$ denotes the identity operator.

3. Maximizing band gap optimization formulation of photonic crystals

For the topology optimization of PhCs under a certain electromagnetic wave polarization mode, our aim is to find the optimal topological configuration of the square lattice that maximizes the gap between two adjacent frequency bands. In this study, the relative band gap f, which refers to the gap-midgap ratio that is independent of the lattice constant, is defined as the objective function for designing PhCs under TM or TE mode

$$\max \quad \textrm{ }f = \frac{{{\omega ^{\textrm{top}}} - {\omega ^{\textrm{bot}}}}}{{{{({{\omega^{\textrm{top}}} + {\omega^{\textrm{bot}}}} )} / 2}}}, $$
where the upper bound of the band gap ${\omega ^{\textrm{top}}}$ and the lower bound ${\omega ^{\textrm{bot}}}$ have different forms:
  • (i) For the m-th TM band gap between the adjacent bands m and $m + 1$:
    $${\omega ^{\textrm{top}}} = \min ({\omega_{m + 1}^{\textrm{TM}}({\textbf k} )} )\textrm{, }{\omega ^{\textrm{bot}}} = \max ({\omega_m^{\textrm{TM}}({\textbf k} )} )$$
1. (ii) For the n-th TE band gap between the adjacent bands n and $n + 1$:
$${\omega ^{\textrm{top}}} = \min ({\omega_{n + 1}^{\textrm{TE}}({\textbf k} )} ),\textrm{ }{\omega ^{\textrm{bot}}} = \max ({\omega_n^{\textrm{TE}}({\textbf k} )} )$$

When designing a complete band gap of PhCs, which is constructed by simultaneously considering the n-th TE band gap and a certain TM band gap, the optimization objective is defined as:

$$\begin{array}{l} \max \quad \textrm{ }f = \max \left\{ {\frac{{\omega_{n,1}^{\textrm{top}} - \omega_{n,1}^{\textrm{bot}}}}{{{{({\omega_{n,1}^{\textrm{top}} + \omega_{n,1}^{\textrm{bot}}} )} / 2}}},\frac{{\omega_{n,2}^{\textrm{top}} - \omega_{n,2}^{\textrm{bot}}}}{{{{({\omega_{n,2}^{\textrm{top}} + \omega_{n,2}^{\textrm{bot}}} )} / 2}}},\ldots ,\frac{{\omega_{n,m}^{\textrm{top}} - \omega_{n,m}^{\textrm{bot}}}}{{{{({\omega_{n,m}^{\textrm{top}} + \omega_{n,m}^{\textrm{bot}}} )} / 2}}},\ldots } \right\},\\ \omega _{n,m}^{\textrm{top}} ={-} \frac{1}{\gamma }\ln \{{\textrm{exp} [{ - \gamma \cdot \min ({\omega_{n + 1}^{\textrm{TE}}({\textbf k} )} )} ]\textrm{ + }\textrm{exp} [{ - \gamma \cdot \min ({\omega_{m + 1}^{\textrm{TM}}({\textbf k} )} )} ]} \}\textrm{ + }\frac{1}{\gamma }\ln (2 ),\\ \omega _{n,m}^{\textrm{bot}} = \frac{1}{\gamma }\ln \{{\textrm{exp} [{\gamma \cdot \max ({\omega_n^{\textrm{TE}}({\textbf k} )} )} ]\textrm{ + }\textrm{exp} [{\gamma \cdot \max ({\omega_m^{\textrm{TM}}({\textbf k} )} )} ]} \}- \frac{1}{\gamma }\ln (2 ), \end{array}$$
where $\gamma > 1$ is the aggregation parameter of the Kreisselmeire-Steinhauser function.

It is assumed that the PhC is composed of two dielectric materials with permittivity ${\varepsilon _1}$ and ${\varepsilon _2}$ (where ${\varepsilon _1}$<${\varepsilon _2}$). A discrete variable $\chi ({\textbf x} )\in \{{0,1} \},\;{\textbf x} \in {\Omega _{\textrm{des}}}$ is introduced to represent the lattice topology, in which $\chi = 0$ means the dielectric material with ${\varepsilon _1}$ while $\chi = 1$ means the dielectric material with ${\varepsilon _2}$. Then, we constructed the topology optimization formulation for the PhCs design as follow:

$$\begin{array}{{ll}} {\mathop {\max }\limits_{\chi ({\textbf x} )} }&{f(\chi )}\\ {\textrm{s}\textrm{.t}.}&{\left[ {{\textbf K}({\bf k} )- {{\left( {\frac{\omega }{c}} \right)}^2}{\textbf M}\textrm{ }} \right]{\textbf u} = 0,\textrm{ }{\bf k} \in [{\mathrm{\Gamma } \to \textrm{X} \to \textrm{M} \to \mathrm{\Gamma }} ]}\\ {}&{\chi ({\textbf x} )\textrm{ = }0\textrm{ }\;or\textrm{ }1,\quad {\textbf x} \in {\Omega _{\textrm{des}}}} \end{array}$$

The above objective function is highly nonlinear and non-differentiable. Although some aggregated techniques [25,52,53] that construct a differentiable global function approximately can be adopted, the success of the gradient-based optimization process usually cannot be guaranteed due to the initial design dependence and the multiple local solutions of the objective functions. In this study, a general gradient-free strategy that combines the MFSE model and the sequential Kriging-based optimization algorithm is proposed to achieve the maximal band gap of 2D PhCs.

4. Gradient-free topology optimization framework

4.1 Material-field topological description of the square lattice

To represent the structural topologies of the square lattices of 2D PhCs, the material-field series expansion (MFSE) method [54] is adopted, which can clearly describe the topologies boundaries with a small number of design variables. In the topology optimization problem of maximizing the photonic crystal band gaps, we only take 1/4 of the unit cell as the diamond symmetric design domain ${\Omega _{\textrm{des}}}$ (see Fig. 2) due to the periodicity and symmetry of the PhC unit cell. Introducing a material-field function $\varphi ({\textbf x} )\in [ - 1,\;\;1]$, ${\textbf x} \in {\Omega _{\textrm{des}}}$ with certain spatial correlation, the distribution of dielectric materials in the design domain ${\Omega _{\textrm{des}}}$ is determined by the value of material-field function as follows.

$$\chi ({\textbf x} )\textrm{ = }\left\{ {\begin{array}{{l}} {\textrm{0,}\quad \textrm{if } - \textrm{1} \le \varphi ({\textbf x} )\le 0,\;\;\;\textrm{ occupied by material}\;\textrm{with }{\varepsilon_1}}\\ {\textrm{1,}\quad \textrm{if }\;0 < \varphi ({\textbf x} )\le 1,\textrm{ }\;\;\;\;\;\;\textrm{occupied by material}\;\textrm{with }{\varepsilon_2}} \end{array}} \right.$$

Here, the values of material-field function $\varphi ({\textbf x} )$ at any two positions (such as ${{\textbf x}_i},{{\textbf x}_j}$) are spatially correlative. Considering the diamond symmetry of the dielectric material distribution in the design domain ${\Omega _{\textrm{des}}}$, the spatial correlation $C({{{\textbf x}_i},{{\textbf x}_j}} )$ can be defined by a distance-related function that has an exponential form as

$$C({{{\textbf x}_i},{{\textbf x}_j}} )= \textrm{exp}\left( {\frac{{ - \min ({{{{\big \Vert}{{{\textbf x}_i} - {{\textbf x}_j}} {\big \Vert}}^2},{{{\big \Vert}{{{\textbf x}_i} - {{\tilde{{\textbf x}}}_j}} {\big \Vert}}^2}} )}}{{l_c^2}}} \right)$$
where ${\tilde{{\textbf x}}_j}$ is the symmetrical position of ${{\textbf x}_j}$ along the diagonal line of the design domain, $l_c^{}$ is the correlation length of the material field which controls the degree of topology complexity and it is typically 10–30% of the dimensions of the design domain. According to the detailed discussions on the impact of the correlation length $l_c^{}$ in [54], as the value of $l_c^{}$ decreases more detailed microstructural geometry features will be obtained in the optimized design, which may be beneficial to high-order frequency band gap design.

 figure: Fig. 2.

Fig. 2. 1/4 symmetric topology design domain of the PhC unit cell and the spatial correlation definition of the diamond symmetric material-field, where the red box indicates the design domain, the variables ${d_1}$ and ${d_2}$ denotes the Euclidean distance between point ${{\textbf x}_i}$ and points ${{\textbf x}_j}$ and ${\tilde{{\textbf x}}_j}$, respectively.

Download Full Size | PDF

According to the bounded field theory [55], the material field $\varphi ({\textbf x} )$ can be linearly expanded on the basis of the eigenpairs of the material-field’s correlation function. In the numerical implementation, the whole design domain is discretized into ${N_\textrm{P}}$ observation points ${{\textbf x}_i}\;\;({i = 1,2,\ldots ,{N_\textrm{P}}} )$, which are uniformly distributed in ${\Omega _{\textrm{des}}}$. Therefore, the material-field function $\varphi ({\textbf x} )$ is defined as

$$\varphi ({\textbf x},\eta ) = \sum\limits_{k = 1}^{{N_\textrm{P}}} \, \left( {{\eta_k}\frac{1}{{\sqrt {{\lambda_k}} }}{\mathbf \psi }_k^\textrm{T}{{\textbf C}_d}({\textbf x})} \right),\;\;\;\textrm{ }{\textbf x} \in {\Omega _{des}}$$
where the vector ${{\textbf C}_d}({\textbf x}) = {\{{C({{\textbf x},{{\textbf x}_1}} )\;\;C({{\textbf x},{{\textbf x}_2}} )\;\; \cdots \;\;C({{\textbf x},{{\textbf x}_{{N_\textrm{P}}}}} )} \}^\textrm{T}}$ is calculated by Eq. (14), ${\eta _k}\;\;({k = 1,2,\ldots {N_\textrm{p}}} )$ are the coefficients to be determined, ${\lambda _k}$ and ${{\mathbf \psi }_k}$ are the k-th eigenvalue and eigenvector of the correlation matrix ${\textbf C}$, respectively. Note that the correlation matrix ${\textbf C}$ is a ${N_\textrm{p}} \times {N_\textrm{p}}$ dimensional, symmetric positive-definite real matrix composed of ${N_\textrm{p}}$ vectors ${{\textbf C}_d}({\textbf x})$. Due to the contribution of eigenvectors with very small eigenvalues to the material field description is negligible, the series expansion function in Eq. (15) can be truncated with the M largest eigenvalues. For a given truncation error $\varepsilon > 0$, the value of M is determined by
$$\sum\limits_{k = 1}^M \, {\lambda _k} \approx ({1 - \varepsilon } )\cdot \left( {\sum\limits_{k = 1}^{{N_P}} \, {\lambda_k}} \right)$$
where the value of truncation error is chosen to be 0.005 in this study, and the corresponding $M \approx \textrm{50}$. Such a setting is enough for ensuring the accurate approximation of the material field of the PhC unit cell (see Eq. (15)). For the details of choosing $\varepsilon$, the readers can refer to [40]

Consequently, the reduced series expansion of the material field can be expressed as

$$\varphi ({\textbf x}) \approx \sum\limits_{k = 1}^M {{\eta _k}\frac{{{\mathbf \psi }_k^\textrm{T}{{\textbf C}_d}({\textbf x})}}{{\sqrt {{\lambda _k}} }}\,} = {{\mathbf \eta }^\textrm{T}}{{\mathbf \Lambda }^{ - 1/2}}{{\mathbf \Phi }^\textrm{T}}{{\textbf C}_d}({\textbf x})$$
where ${\mathbf \eta } = {\{{{\eta_1},{\eta_2},\; \cdots \;,\;{\eta_M}} \}^\textrm{T}}$ is the vector of undetermined coefficients or design variables, ${\mathbf \Lambda } = diag({{\lambda_1},{\lambda_2}, \cdots ,{\lambda_M}} )$ and ${\mathbf \Phi } = \{{{{\mathbf \psi }_1},\;{{\mathbf \psi }_2},\; \cdots \;,{{\mathbf \psi }_M}} \}$ are respectively the diagonal matrix composed of the first M largest eigenvalues and corresponding eigenvectors.

Meanwhile, in order to ensure the boundness of the material field, the material-field function $\varphi ({\textbf x} )$ at every observation point ${{\textbf x}_i}\textrm{ }(i = 1,2,\ldots ,{N_\textrm{p}})$ is constrained as $\varphi ({{{\textbf x}_i}} )\in [ - 1,\;1]$. The constraint form is as follow

$${\boldsymbol{\mathrm{\eta}}^\textrm{T}}{{\textbf W}_i}\boldsymbol{\mathrm{\eta}} \le 1,\textrm{ }({i = 1,2,\ldots ,{N_\textrm{p}}} )$$
where the characteristic matrices ${{\textbf W}_{\textbf i}} = {{\mathbf \Lambda }^{ - 1/2}}{{\mathbf \Phi }^\textrm{T}}{{\textbf C}_d}({{\textbf x}_i}){{\textbf C}_d}{({{\textbf x}_i})^\textrm{T}}{\mathbf \Phi }{{\mathbf \Lambda }^{ - 1/2}}$.

4.2 MFSE-based topology optimization formulation

To satisfy the topology representation in Eq. (13), a Heaviside function projection is used to map the value of the material field onto the corresponding lattice topology. The formulation is as follows

$$\tilde{\varphi }({\textbf x} )= \left\{ {\begin{array}{{c}} {\varphi ({\textbf x}) \cdot {e^{ - \beta }} + 1 - {e^{ - \beta \varphi ({\textbf x})}},\;\;\;\;\;\textrm{if }\varphi ({\textbf x}) \ge 0}\\ {\varphi ({\textbf x}) \cdot {e^{ - \beta }} - 1 + {e^{\beta \varphi ({\textbf x})}},\;\;\;\;\;\textrm{if }\varphi ({\textbf x}) < 0} \end{array}} \right.$$
where $\beta $ is a parameter controlling the smoothness of the mapping. If $\beta $ tends to infinity, the mapped material field $\tilde{\varphi }({\textbf x} )$ is equal to a discrete function with the value of −1 or 1. In the practical optimization process, $\beta $ will be gradually increased from 0 to 200 with a continuation approach.

Supposing that the design domain is discretized into ${N_{\textrm{ele}}}$ finite elements, we select the value $\tilde{\varphi }({{\textbf x}_j})$ at the center ${{\textbf x}_j}$ of each element as the elemental average value of material field. Note that $\tilde{\varphi }({{\textbf x}_j})$ can be easily calculated by Eqs. (17) and (19). Then, the relative density ${\rho _j}$ of the j-th element is related to $\tilde{\varphi }({{{\textbf x}_j}} )$ as

$${\rho _j} = \frac{{1 + \tilde{\varphi }({{\textbf x}_j})}}{2},\textrm{ (}j = 1,2,\ldots ,{N_{\textrm{ele}}}\textrm{),}$$
and the spatial permittivity distribution can be described with the following interpolations
$$\varepsilon ({{\rho_j}} )= {\varepsilon _1} + ({{\varepsilon_2} - {\varepsilon_1}} ){\rho _j},\;\;\;\;\;\textrm{ for TM mode}$$
$$\frac{1}{{\varepsilon ({{\rho_j}} )}} = \frac{1}{{{\varepsilon _1}}} + \left( {\frac{1}{{{\varepsilon_2}}} - \frac{1}{{{\varepsilon_1}}}} \right){\rho _j},\;\;\;\;\;\textrm{for TE mode}$$

Therefore, the MFSE-based topology optimization formulation for the 2D PhCs design is expressed as

$$\begin{array}{{ll}} {\textrm{Find}}&{\boldsymbol{\mathrm{\eta}} = {{\{{{\eta_1},{\eta_2},\ldots ,{\eta_M}} \}}^\textrm{T}}}\\ {\mathop {\max }\limits_{\boldsymbol{\mathrm{\eta}}} }&{f({\boldsymbol{\mathrm{\eta}}} )}\\ {\textrm{s}\textrm{.t}\textrm{.}}&{\left[ {{\bf K}({\bf k} )- {{\left( {\frac{\omega }{c}} \right)}^2}{\bf M}\textrm{ }} \right]{\textbf u} = 0,\textrm{ }{\bf k} \in [{\mathrm{\Gamma } \to \textrm{X} \to \textrm{M} \to \mathrm{\Gamma }} ]}\\ {}&{{{\boldsymbol{\mathrm{\eta}}}^\textrm{T}}{{\bf W}_i}{\boldsymbol{\mathrm{\eta}}} \le 1,\textrm{ }({i = 1,2,\ldots ,{N_\textrm{p}}} )} \end{array}$$
where the coefficients ${\boldsymbol{\mathrm{\eta}}} = {\{{{\eta_1},{\eta_2},\; \cdots \;,\;{\eta_M}} \}^\textrm{T}}$ are taken as the design variables. With the MFSE method, the number of design variables in the topology optimization problem has been successfully reduced to M. In general, M is much smaller than the given number of observation points ${N_\textrm{p}}$ and the number of finite elements ${N_{\textrm{ele}}}$. It should be noted that the finite element mesh is decoupled from the lattice topology represented by the material field, which means that a coarse-to-fine mesh strategy of FE simulations can be easily implemented for improving the optimization efficiency. In this study, the entire unit cell is first discretized into 40×40 uniform square finite elements and then refined to 80×80 mesh during the optimization procedure.

4.3 Sequential Kriging-based optimization algorithm

To search the optimum of the band gap topology optimization problem in Eq. (23) without the sensitivity information, the sequential Kriging-based optimization algorithm [40] is used, which is based on the Kriging surrogate model and adopts an adaptive design domain adjusting technique. Thus, the topology optimization model in Eq. (23) is converted to the unconstrained and minimization formulation as follow

$$\mathop {\min \textrm{ }}\limits_{\boldsymbol{\mathrm{\eta}}} \textrm{ }{f_{obj}}({\boldsymbol{\mathrm{\eta}}}) ={-} f({\boldsymbol{\mathrm{\eta}}} )+ {p_0} \cdot \mathop {\max \textrm{ }}\limits_i ({({{{\boldsymbol{\mathrm{\eta}}}^{\textbf T}}{{\textbf W}_i}{\boldsymbol{\mathrm{\eta}}} - 1} ),0} )\textrm{, }$$
where ${p_0}$ is the barrier multiplier.

Figure 3 shows the flow chart of PhCs optimization design based on the proposed gradient-free optimization strategy. In the sequential Kriging-based optimization algorithm, the unconstrained problem in Eq. (24) is reformulated into several sub-optimization problems. For the first sub-optimization problem within the initial design space ${\boldsymbol{\mathrm{\eta}}} \in {\Omega _0}$, the Kriging surrogate model is constructed based on the Latin hypercube sampling and updated by using the sampling criteria such as the expectation improvement (EI) and the minimization prediction (MP) techniques. Then, the k-th sub-optimization problem is performed sequentially by Kriging-based optimization algorithm after adaptively shrinking the sub-design space ${\boldsymbol{\mathrm{\eta}}} \in {\Omega _k}$, which centers at the optimum of the previous sub-optimization problem. The convergence criterion for each sub-optimization problem is that the relative change of the objective functions among ten successive adding samples should be less than 0.5%. For a detailed introduction of the sequential Kriging surrogate model algorithm, please reference the literature [40].

 figure: Fig. 3.

Fig. 3. Flow chart of the gradient-free topology optimization strategy for photonic crystal band gap design.

Download Full Size | PDF

5. Numerical examples and discussion

In this section, we consider the topological design of 2D PhCs with the square lattice, which is composed of air (${\varepsilon _1} = 1$) and gallium arsenide (i.e. GaAs, ${\varepsilon _2} = 11.4$). The lattice constant is $a\textrm{ } = \textrm{ }1\textrm{um}$ and the correlation length of the material field is set to be ${l_\textrm{c}} = 0.075\textrm{um}$.

In this study, we set the sample centre of the initial sub-optimization problem at ${\boldsymbol{\mathrm{\eta}}}{{\kern 1pt} _c} = 0$, and the corresponding design space is denoted by ${\Omega _0}\textrm{ = }\{{{{|{{\boldsymbol{\mathrm{\eta}}} - {{\boldsymbol{\mathrm{\eta}}}_{\boldsymbol c}}} |}_\infty } \le 0.4} \}$.

5.1 Band gap design of photonic crystals under TM mode

The PhCs design problem under the TM mode is first considered. Figure 3 gives the optimized PhCs geometries with the maximum band gaps that from the first to the tenth band, together with the corresponding optimized band diagrams. These optimized PhCs are represented by 3×3 unit cells, in which the red and white regions respectively represent GaAs and air, and the unit cell is outlined with the dashed box.

As can be seen from Fig. 4, isolated GaAs with circular or more complex shapes is scattered over the square lattice to achieve the maximal band gap of PhCs under the TM mode. The optimized designs obtained with the present method have similar topologies to those reported in [20,36], though minor differences also exist. In general, the size of isolated GaAs structures decreases as the increasing order of band gap. However, for high-order frequency bands (e.g., the seventh to tenth band gaps), the size and positions of the isolated GaAs structures in the unit cell are almost unchanged, but their shapes become different. From the band diagrams of optimized designs in Fig. 4, it is seen that large band gaps are successfully obtained for all the specified frequency bands. The maximum relative band gap under TM mode is 46.20%, which locates between the seventh and the eighth bands.

 figure: Fig. 4.

Fig. 4. Optimized 3 × 3 unit cells and their band diagrams under the TM mode. Figures (a) to (j) correspond to the first to the tenth band, where the red and white regions respectively represent GaAs and air, and the dashed box represents the unit cell. In the band diagram, the vertical axis is the normalized frequency, and the horizontal axis is the wave vector on the boundary of the irreducible Brillouin zone. The sweep path is $\mathrm{\Gamma } \to \textrm{X} \to \textrm{M} \to \mathrm{\Gamma }$. Label TM mode denotes the TM polarization of the electromagnetic field. The dispersion curves of different bands are represented by different colours. The cyan region represents the band gap, and the value in it is the relative band gap.

Download Full Size | PDF

The iteration history of the optimization problem that maximizing the seventh band gap is given in Fig. 5. As mentioned above, the proposed method uses ${\Omega _0}$ as the initial sub-design space. The objective function values of the initial samples in the first sub-optimization loop are greater than zero, which indicates there is no band gap. With the sequential sampling and updating of 15 sub-optimization problems (represented by different colours), the optimization procedure converged after 2715 finite element calculations. Compared with conventional gradient-based optimization methods which may converge in a few hundreds of iterations, the proposed method does require more FE evaluations due to the Latin hypercube sampling (i.e. FE evaluations). However, these samples contribute the present algorithm’s features of relatively strong global search capability and less dependence on initial designs.

 figure: Fig. 5.

Fig. 5. Optimization process of maximizing the seventh band gap for PhCs under TM mode. Different colours represent different sub-optimization problems.

Download Full Size | PDF

The objective function gradually decreases from a positive value to a negative value, which means the band gap occurred and increased, and finally converged to the maximum value of 46.20%. Figure 6 gives the evolution history of the unit cell topology. The clear 0–1 topology is obtained in the final loop. Due to the spatial correlation of the material-field function, the MFSE-based optimization model is inherently capable of avoiding the checkerboard patterns and mesh dependency.

 figure: Fig. 6.

Fig. 6. Evolution history of the unit cell topology for maximizing the seventh band gap of PhCs under TM mode. Obtained topologies after (a) the first, (b) the second, (c) the fourth, (d) the sixth, (e) the eighth sub-optimization problems, and (f) the final optimum.

Download Full Size | PDF

Figure 7 shows the optimization results of maximizing the seventh band gap under the TM mode starting with different initial sample space ${\Omega _0}$. We first consider two different centers ${\boldsymbol{\mathrm{\eta}}}{{\kern 1pt} _\textrm{0}}$ of the initial sample space to produce material fields with uniform density distribution ${\rho _j} = 0.4\textrm{ }({j = 1,2,\ldots ,{N_{\textrm{ele}}}} )$ and ${\rho _j} = 0.6\textrm{ }({j = 1,2,\ldots ,{N_{\textrm{ele}}}} )$, respectively. The optimized results are shown in Figs. 7(a) and 7(b), and their corresponding seventh relative band gaps are 46.04% and 45.32%. From the periodicity point of view, the PhC topology configurations in Figs. 7(a) and 7(b) are very similar.

 figure: Fig. 7.

Fig. 7. Optimized 3 × 3 unit cell designs from different initial sample spaces. (a) Optimized design 1 (optimized from initial Configuration 1 with uniform density distribution ${\rho _j} = 0.4({j = 1,2,\ldots ,{N_{\textrm{ele}}}} )$). (b) Optimized design 2 (optimized from initial Configuration 1 with uniform density distribution ${\rho _j} = 0.6({j = 1,2,\ldots ,{N_{\textrm{ele}}}} )$). (c) Initial Configuration 3 (central crystal column, radius $r/a = 0.25$.). (d) Optimized design 3 corresponding to initial Configuration 3. (e) Initial Configuration 4 (central hole, radius $r/a = 0.25$). (f) Optimized design 4 corresponding to initial Configuration 4.

Download Full Size | PDF

Further, we defined two opposite initial designs (i.e., Configuration 3 and Configuration 4) as shown in the Figs. 7(c) and 7(e): the circular region at the center of the unit cell represents crystal column or hole and its radius is 0.25a, while the surrounding region represents the opposite material. We started the proposed optimization method with the initial sample space centered at Configuration 3 and Configuration 4. The obtained results are shown in Figs. 7(d) and (f), and their corresponding seventh relative band gaps are 45.96% and 46.03%. Obviously, these optimized results are all surprisingly similar to the solution in Fig. 4(g). This reveals that the proposed gradient-free topology optimization method is able to yield meaningful solutions without initial solution dependence for this problem. In addition, Table 1 gives a comparison between the present optimized results and the previously reported results. It is seen that the relative band gaps with the proposed gradient-free method are comparable or larger than the solutions obtained with gradient-based optimization methods [20,25]. Therefore, the proposed method is shown to be effective for the considered band gap topology design of PhCs in the sense that that it has a relatively strong global search capability.

Tables Icon

Table 1. Comparison of the relative band gap solutions of PhCs under TM mode.

5.2 Band gap design of photonic crystals under TE mode

In this section, the band gap topology optimization of PhCs under the TE mode is considered. The optimized PhCs and their corresponding band diagrams are given in Fig. 8. It is observed from the figure that the unit cell topologies under the TE mode are distinctly different from the optimized results under the TM mode, which is in agreement with the conclusion drawn in Ref. [1]: ‘the TM band-gaps are favoured in a lattice of isolated high- $\varepsilon $ regions, and TE band-gaps are favoured in a connected lattice’. Compared with previous reports [20,25] (listed in Table 2), the obtained designs for high-order bands present novel geometric configurations with larger band gaps, which illustrates the photonic band gap optimization problem for the TE mode is more complicated and has many local minima. By using the proposed method, the maximum relative band gap of PhCs under the TE mode is 45.27%, which is the gap between the seventh and eighth bands.

Figure 9 shows the percentage of the GaAs material consumptions of optimized unit cells for each band to the entire unit cell. Obviously, under the TE mode, the optimized unit cells require larger amounts of GaAs material than the ones under the TM mode. Under the TM mode, the optimal topology for maximal band gap appears as small isolated GaAs columns with regular or irregular shapes; while under the TE mode, it requires more GaAs materials to appear as closed-walled structures.

 figure: Fig. 8.

Fig. 8. Optimized 3 × 3 unit cells and their band diagrams under the TE mode. Figures (a) to (j) correspond to the first to the tenth band. Label TE mode denotes the TE polarization of the electromagnetic field.

Download Full Size | PDF

 figure: Fig. 9.

Fig. 9. Volume fraction of GaAs vs. band number in the optimized designs. Material consumption under TM mode is expressed with blue solid line and solid circle, and is expressed with red dashed line and solid circle under TE mode.

Download Full Size | PDF

Tables Icon

Table 2. Comparison of the relative band gap solutions of PhCs under TE mode

5.3 Complete band gap design of photonic crystals

In general, one may overlap the band gap under TM mode with the band gap under TE mode to generate a narrow complete band gap that can prohibit both TM and TE waves. However, due to the different formation mechanisms of PhCs for the TM and TE modes, it is almost impossible to obtain a desired, large complete band gap of PhCs through experience. In previous researches, the genetic algorithm and gradient-based optimization algorithms have been employed to design PhCs with complete band gaps [3739,5661]. These researches showed that the maximization of complete band gap of PhCs via the existing topology optimization method relies heavily on the carefully constructed initial designs.

In this section, we consider several complete band gap designs (i.e., the second to the ninth band gaps for TE mode and a certain band gap for TM mode) to illustrate the effectiveness of the proposed method on the topological design of PhCs. The corresponding solutions and band diagrams are shown in Fig. 10. Here, the band diagram under TE mode is represented by red dashed lines, and the band diagram under TM mode is represented by the blue solid lines. The optimization results show the second-to-ninth order band gaps for TE modes and the third, fourth, eighth, ninth, ninth, twelfth, twelfth, fifteenth order for TM band gaps respectively construct a large complete band gap (the corresponding relative band gap values are 19.31%, 10.53%, 15.93%, 20.98%, 12.75%, 20.46%, 10.18% and 13.38%). It can be seen from the figure that these optimized cell topologies are mostly composed of subpartitioning closed-wall structures and irregular crystal columns. Compared with the PhCs under a single polarization mode, the optimal topology configuration of complete band gap is more complicated and therefore is difficult to predict in advance. The proposed optimization model and gradient-free method can effectively yield a desired solution without requiring a carefully choice of the initial guesses.

 figure: Fig. 10.

Fig. 10. Optimized 3 × 3 unit cells photonic crystals with a complete band gap containing: (a) the 2th order TE and 3th order TM; (b) the 3th order TE and 4th order TM; (c) the 4th order TE and 8th order TM; (d) the 5th order TE and 9th order TM; (e) the 6th order TE and 9th order TM; (f) the 7th order TE and 12th order TM; (g) the 8th order TE and 12th order TM; (h) the 9th order TE and 12th order TM.

Download Full Size | PDF

5.4 Application in photonic crystal slabs

In this example, we explore the performances of photonic crystals with two-dimensional periodicity and a finite thickness (the so-called photonic crystal slabs [1]). The band diagrams of two-dimensional photonic crystals and the corresponding photonic crystal slabs yield some similarities. Although the finite thickness introduces new behavior (such as the light cone), some of the two-dimensional photonic crystals can be directly used to construct the photonic crystal slabs [62].

The performances of slabs depend on factors such as light cone, thickness, substrates, and mirror symmetries. If the photonic crystal slabs have mirror symmetry with respect to the z=0 plane (where z axis is the out-of-plane direction), the electromagnetic wave modes can be classified as odd mode (TM-like) and even mode (TE-like). Conventionally, rod slabs are adopted to achieve the band gaps under the TM-like mode, and hole slabs are used to achieve the band gaps under the TE-like mode. The chosen thicknesses of slabs is related to the lattice constants. As is known that the hole slabs can achieve band gaps with lower thicknesses than rod slabs [1], we extrude the two optimized two-dimensional configurations (i.e., the optimized designs in Fig. 4(a) and Fig. 8(a) for the 1st order TM/TE mode) to construct photonic crystal slabs with thicknesses of 2a and 0.6a, respectively. These two slabs are both surrounded by air, and their structures and corresponding band gap diagrams are shown in Fig. 11.

 figure: Fig. 11.

Fig. 11. Two photonic crystal slabs obtained by extruding the optimized configurations in Fig. 4(a) and Fig. 8(a), respectively. (a) and (b) are their corresponding band gap diagrams. The gray line is the light line, and the blue/red dashed lines represent the band diagrams under TM/TE-like modes, respectively. The gaps of the rod/hole slabs are represented by the cyan region.

Download Full Size | PDF

It can be seen that these two photonic crystal slabs both have relatively large band gaps below the light line. Under the TM-like mode the relative band gaps of the rod slab is 19.81%, and that of the hole slab under TE-like mode is 16.82%. Some narrow band gaps and even complete band gaps below the light line can also be observed in Figs. 11(a) and (b), which is not observed with the corresponding two-dimensional photonic crystals. Though it is interesting to explore the band gap feature of photonic crystal slabs, the focus of this study is the design of maximizing of two-dimensional photonic crystals.

6. Conclusions

In this study, a new gradient-free topology optimization strategy is proposed to obtain the maximum band gap design of 2D PhCs, considering the TM mode, the TE mode or the complete mode. With the MFSE-based topology representation of PhCs unit cell, the band gap design problem is converted into a no-constrained form with a few design variables and is efficiently solved with a surrogate model optimization algorithm.

Numerical results indicate that the proposed method can effectively yield PhCs with relative large band gaps under different polarization modes. For the 2D PhCs 5with square lattice, the maximum relative band gaps we obtained are 46.20% for the TM mode and 45.27% for the TE mode. To our best knowledge, both of them are the largest band gaps that have been reported so far. Furthermore, the present method does not require a careful choice of the experience-based guess of the PhCs unit cell, especially for the design problem of complete band gaps.

The proposed algorithm exhibits relatively strong global search capability, and is easy to implement without the sensitivity information. In further studies, it will be extended to deal with other topology optimization problems related to photonic crystals, such as waveguide structures, photonic crystal slabs, specified frequency band gaps or complex types of PhCs lattices (e.g. hexagonal lattice) where the multiple local optimum feature and dependency on initial designs are significant.

Funding

the LiaoNing Revitalization Talents Program of China (XLYC1807187); National Key Research and Development Program of China (2017YFB0203604); National Natural Science Foundation of China (11772077).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, Photonic crystals: molding the flow of light (Princeton University Press, 2008).

2. J. G. Fleming, S. Y. Lin, I. El-Kady, R. Biswas, and K. M. Ho, “All-metallic three-dimensional photonic crystals with a large infrared bandgap,” Nature 417(6884), 52–55 (2002). [CrossRef]  

3. S. Fan, P. R. Villeneuve, J. D. Joannopoulos, and H. A. Haus, “Channel drop filters in photonic crystals,” Opt. Express 3(1), 4–11 (1998). [CrossRef]  

4. S. Fan, J. N. Winn, A. Devenyi, J. C. Chen, R. D. Meade, and J. D. Joannopoulos, “Guided and defect modes in periodic dielectric waveguides,” J. Opt. Soc. Am. B 12(7), 1267–1272 (1995). [CrossRef]  

5. M. Soljačić, M. Ibanescu, S. G. Johnson, Y. Fink, and J. D. Joannopoulos, “Optimal bistable switching in nonlinear photonic crystals,” Phys.Rev.E 66(5), 055601 (2002). [CrossRef]  

6. M. F. Yanik and S. Fan, “Stopping and storing light coherently,” Phys.Rev.A 71(1), 013803 (2005). [CrossRef]  

7. E. E. Hart, A. Sóbester, K. Djidjeli, M. Molinari, K. S. Thomas, and S. J. Cox, “A geometry optimization framework for photonic crystal design,” Photonics Nanostruct. Fundam. Appl. 10(1), 25–35 (2012). [CrossRef]  

8. M. Doosje, B. J. Hoenders, and J. Knoester, “Photonic bandgap optimization in inverted fcc photonic crystals,” J. Opt. Soc. Am. B 17(4), 600–606 (2000). [CrossRef]  

9. J. H. Wu, A. Q. Liu, and L. K. Ang, “Band gap optimization of finite photonic structures using apodization method,” J. Appl. Phys. 100(8), 084309 (2006). [CrossRef]  

10. M. M. Hossain, G. Chen, B. Jia, X. H. Wang, and M. Gu, “Optimization of enhanced absorption in 3D-woodpile metallic photonic crystals,” Opt. Express 18(9), 9048–9054 (2010). [CrossRef]  

11. C. Mee, P. Contu, and P. Pintus, “One-dimensional photonic crystal design,” J. Quant. Spectros. Radiat. Transfer 111(1), 214–225 (2010). [CrossRef]  

12. B. Thubthimthong and F. Chollet, “Design and simulation of a tunable photonic band gap filter,” Microelectr. Eng. 85(5-6), 1421–1424 (2008). [CrossRef]  

13. M. Bendsøe and N. Kikuchi, “Generating optimal topologies in structural design using a homogenization method,” Comput. Methods Appl. Mech. Eng. 71(2), 197–224 (1988). [CrossRef]  

14. O. Sigmund and K. Maute, “Topology optimization approaches,” Struct. Multidisc. Opt. 48(6), 1031–1055 (2013). [CrossRef]  

15. J. D. Deaton and R. V. Grandhi, “A survey of structural and multidisciplinary continuum topology optimization: post 2000,” Struct. Multidisc. Opt. 49(1), 1–38 (2014). [CrossRef]  

16. C. Y. Kao, S. Osher, and E. Yablonovitch, “Maximizing band gaps in two-dimensional photonic crystals by using level set methods,” Appl. Phys. B 81(2-3), 235–244 (2005). [CrossRef]  

17. P. I. Borel, A. Harpøth, L. H. Frandsen, M. Kristensen, P. Shi, J. S. Jensen, and O. Sigmund, “Topology optimization and fabrication of photonic crystal structures,” Opt. Express 12(9), 1996–2001 (2004). [CrossRef]  

18. Q. Liang, D. Li, Y. Gai, and H. Han, “Ultra-wide bandgap of gradient dielectric constant photonic crystal,” Mater. Lett. 79, 48–50 (2012). [CrossRef]  

19. A. Takezawa and M. Kitamura, “Phase field method to optimize dielectric devices for electromagnetic wave propagation,” J. Comput. Phys. 257, 216–240 (2014). [CrossRef]  

20. H. Men, N. C. Nguyen, R. M. Freund, P. A. Parrilo, and J. Peraire, “Bandgap Optimization of Two-Dimensional Photonic Crystals Using Semidefinite Programming and Subspace Methods,” J. Comput. Phys. 229(10), 3706–3725 (2010). [CrossRef]  

21. H. Men, N. C. Nguyen, R. M. Freund, K. M. Lim, P. A. Parrilo, and J. Peraire, “Design of photonic crystals with multiple and combined band gaps,” Phys. Rev. E 83(4), 046703 (2011). [CrossRef]  

22. H. F. Zhang, G. W. Ding, H. M. Li, and S. B. Liu, “Complete photonic band gaps and tunable self-collimation in the two-dimensional plasma photonic crystals with a new structure,” Phys. Plasmas 22(2), 022105 (2015). [CrossRef]  

23. H. W. Dong, Y. S. Wang, T. X. Ma, and X. X. Su, “Topology optimization of simultaneous photonic and phononic bandgaps and highly effective phoxonic cavity,” J. Opt. Soc. Am. B 31(2014).

24. H. Dong, Y. Wang, and C. Zhang, “Topology Optimization of Chiral Phoxonic Crystals With Simultaneously Large Phononic and Photonic Bandgaps,” IEEE Photonics J. 9(2), 1–16 (2017). [CrossRef]  

25. F. Meng, X. Huang, and B. Jia, “Bi-directional evolutionary optimization for photonic band gap structures,” J. Comput. Phys. 302, 393–404 (2015). [CrossRef]  

26. F. Meng, Y. Li, S. Li, H. Lin, B. Jia, and X. Huang, “Achieving Large Band Gaps in 2D Symmetric and Asymmetric Photonic Crystals,” J. Lightwave Technol. 35(9), 1670–1676 (2017). [CrossRef]  

27. F. Meng, B. Jia, and X. Huang, “Topology-Optimized 3D Photonic Structures with Maximal Omnidirectional Bandgaps,” Adv. Theory Simul. 1(12), 1800122 (2018). [CrossRef]  

28. X. Liang and S. G. Johnson, “Formulation for scalable optimization of microcavities via the frequency-averaged local density of states,” Opt. Express 21(25), 30812 (2013). [CrossRef]  

29. F. Wang, R. E. Christiansen, Y. Yi, J. Mrk, and O. Sigmund, “Maximizing the quality factor to mode volume ratio for ultra-small photonic crystal cavities,” Appl. Phys. Lett. 113(24), 241101 (2018). [CrossRef]  

30. A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vukovi, “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nat. Photonics 9(2015).

31. Z. Lin, V. Liu, R. Pestourie, and S. G. Johnson, “Topology optimization of freeform large-area metasurfaces,” Opt. Express 27(11), 15765–15775 (2019). [CrossRef]  

32. H. Chung and O. D. Miller, “High-NA, Achromatic, Visible-Frequency Metalenses by Inverse Design,” in Frontiers in Optics + Laser Science APS/DLS, OSA Technical Digest (Optical Society of America, 2019), JW4A.67.

33. R. E. Christiansen, F. Wang, O. Sigmund, and S. Stobbe, “Designing Photonic Topological Insulators with Quantum-Spin-Hall Edge States using Topology Optimization,” Nanophotonics (2019).

34. B. Xia, H. Fan, and T. Liu, “Topologically Protected Edge States of Phoxonic Crystals,” Int. J. Mech. Sci. 155, 197–205 (2019). [CrossRef]  

35. J. S. Jensen and O. Sigmund, “Topology optimization for nano-photonics,” Laser Photonics Rev. 5(2), 308–321 (2011). [CrossRef]  

36. O. Sigmund and K. Hougaard, “Geometric Properties of Optimal Photonic Crystals,” Phys. Rev. Lett. 100(15), 153904 (2008). [CrossRef]  

37. F. Meng, S. Li, Y. F. Li, B. Jia, and X. Huang, “Microstructural design for 2D photonic crystals with large polarization-independent band gaps,” Mater. Lett. 207, 176–178 (2017). [CrossRef]  

38. S. Li, H. Lin, F. Meng, D. Moss, X. Huang, and B. Jia, “On-Demand Design of Tunable Complete Photonic Band Gaps based on Bloch Mode Analysis,” Sci. Rep. 8(1), 14283 (2018). [CrossRef]  

39. Y. Chen, F. Meng, G. Li, and X. Huang, “Designing photonic materials with complete band gaps by topology optimization,” Smart Mater. Struct. 28(2018).

40. Y. Luo, J. Xing, and Z. Kang, “Topology optimization using material-field series expansion and Kriging-based algorithm: An effective non-gradient method,” Comput. Methods Appl. Mech. Eng. 364, 112966 (2020). [CrossRef]  

41. X. Zhang, J. Xing, P. Liu, Y. Luo, and Z. Kang, “Realization of full and directional band gap design by non-gradient topology optimization in acoustic metamaterials,” Extreme Mech. Lett. 42, 101126 (2021). [CrossRef]  

42. P. Liu, Y. Yan, X. Zhang, Y. Luo, and Z. Kang, “Topological design of microstructures using periodic material-field series-expansion and gradient-free optimization algorithm,” Mater. Des. 199, 109437 (2021). [CrossRef]  

43. P. Griffin, P. Nagel, and R. D. Koshel, “The plane-wave expansion method,” J. Math. Phys. 15(11), 1913–1917 (1974). [CrossRef]  

44. V. F. D. Poggetto and A. L. Serpa, “Elastic wave band gaps in a three-dimensional periodic metamaterial using the plane wave expansion method,” Int. J. Mech. Sci. 184, 105841 (2020). [CrossRef]  

45. A. Taflove and S. Hagness, “Computational Electrodynamics: The Finite-Difference Time-Domain Method,” Artech House 5, 629–670 (2005). [CrossRef]  

46. M. Qiu, “Analysis of guided modes in photonic crystal fibers using the finite-difference time-domain method,” Microw. Opt. Technol. Lett. 30(5), 327–330 (2001). [CrossRef]  

47. I. V. Soboleva, V. V. Moskalenko, and A. A. Fedyanin, “Giant Goos-Hnchen Effect and Fano Resonance at Photonic Crystal Surfaces,” Phys. Rev. Lett. 108(12), 123901 (2012). [CrossRef]  

48. C. Luo and Y. Xie, “The Design of Optical Switch for One-Dimension Photonic Crystal Constructed with TiO2 and Liquid Crystal,” Adv. Mat. Res. 721, 90–94 (2013). [CrossRef]  

49. J. Jin, “The Finite Element Method in Electromagnetics,” John Wiley & Sons, New York (2002).

50. I. Andonegui and A. J. Garcia-Adeva, “The finite element method applied to the study of two-dimensional photonic crystals and resonant cavities,” Opt. Express 21(4), 4072–4092 (2013). [CrossRef]  

51. C. Kittel, “Introduction to Solid State Physics (A. Klemm),” Phys. Today 61, 59–60 (2008).

52. G. Kreisselmeier and R. Steinhauser, “Systematic Control Design by Optimizing a Vector Performance Index,” in Computer Aided Design of Control Systems, M. A. Cuenod, ed. (Pergamon, 1980), pp. 113–117.

53. J. He and K. Zhan, “Achieving directional propagation of elastic waves via topology optimization,” Ultrasonics 82, 1–10 (2018). [CrossRef]  

54. Y. Luo and J. Bao, “A material-field series-expansion method for topology optimization of continuum structures,” Comput. Struct. 225, 106122 (2019). [CrossRef]  

55. Y. Luo, J. Zhan, J. Xing, and Z. Kang, “Non-probabilistic uncertainty quantification and response analysis of structures with a bounded field model,” Comput. Methods Appl. Mech. Eng. 347, 663–678 (2019). [CrossRef]  

56. A. Cerjan and S. Fan, “Complete photonic band gaps in supercell photonic crystals,” Phys. Rev. A 96(5), 051802 (2017). [CrossRef]  

57. L. Shen, Z. Ye, and S. He, “Design of two-dimensional photonic crystals with large absolute band gaps using a genetic algorithm,” Phys. Rev. B 68(3), 035109 (2003). [CrossRef]  

58. H. Li, L. Jiang, W. Jia, H. Qiang, and X. Li, “Genetic optimization of two-dimensional photonic crystals for large absolute band-gap,” Opt. Commun. 282(14), 3012–3017 (2009). [CrossRef]  

59. W. Jia, J. Deng, A. Danner, and X. Li, “Optimization and fabrication of the photonic crystals for large absolute band gaps,” in Proceedings of 2011 International Conference on Electronics and Optoelectronics, (2011), V4-152–V154-154.

60. D. Wang, Z. Yu, Y. Liu, P. Lu, L. Han, H. Feng, X. Guo, and H. Ye, “The optimal structure of two dimensional photonic crystals with the large absolute band gap,” Opt. Express 19(20), 19346–19353 (2011). [CrossRef]  

61. P. Shi, K. Huang, and Y. Li, “Photonic crystal with complex unit cell for large complete band gap,” Opt. Commun. 285(13-14), 3128–3132 (2012). [CrossRef]  

62. A. Cerjan and S. Fan, “Complete photonic bandgaps in supercell photonic crystals,” in Conference on Lasers and Electro-Optics, OSA Technical Digest (online) (Optical Society of America, 2018), STh3A.7.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (11)

Fig. 1.
Fig. 1. Left: The periodic lattice with the square lattice and lattice constant a. The white colour denotes air ( ${\varepsilon _1} = 1$ ), the red denotes GaAs ( ${\varepsilon _\textrm{2}} = 1\textrm{1}\textrm{.4}$ ), and the dashed box represents the unit cell; Right: the dashed box represents the first Brillouin zone in the reciprocal lattice. The irreducible zone is the gray triangular and its boundary is denoted by $\mathrm{\Gamma } \to \textrm{X} \to \textrm{M} \to \mathrm{\Gamma }$ .
Fig. 2.
Fig. 2. 1/4 symmetric topology design domain of the PhC unit cell and the spatial correlation definition of the diamond symmetric material-field, where the red box indicates the design domain, the variables ${d_1}$ and ${d_2}$ denotes the Euclidean distance between point ${{\textbf x}_i}$ and points ${{\textbf x}_j}$ and ${\tilde{{\textbf x}}_j}$ , respectively.
Fig. 3.
Fig. 3. Flow chart of the gradient-free topology optimization strategy for photonic crystal band gap design.
Fig. 4.
Fig. 4. Optimized 3 × 3 unit cells and their band diagrams under the TM mode. Figures (a) to (j) correspond to the first to the tenth band, where the red and white regions respectively represent GaAs and air, and the dashed box represents the unit cell. In the band diagram, the vertical axis is the normalized frequency, and the horizontal axis is the wave vector on the boundary of the irreducible Brillouin zone. The sweep path is $\mathrm{\Gamma } \to \textrm{X} \to \textrm{M} \to \mathrm{\Gamma }$ . Label TM mode denotes the TM polarization of the electromagnetic field. The dispersion curves of different bands are represented by different colours. The cyan region represents the band gap, and the value in it is the relative band gap.
Fig. 5.
Fig. 5. Optimization process of maximizing the seventh band gap for PhCs under TM mode. Different colours represent different sub-optimization problems.
Fig. 6.
Fig. 6. Evolution history of the unit cell topology for maximizing the seventh band gap of PhCs under TM mode. Obtained topologies after (a) the first, (b) the second, (c) the fourth, (d) the sixth, (e) the eighth sub-optimization problems, and (f) the final optimum.
Fig. 7.
Fig. 7. Optimized 3 × 3 unit cell designs from different initial sample spaces. (a) Optimized design 1 (optimized from initial Configuration 1 with uniform density distribution ${\rho _j} = 0.4({j = 1,2,\ldots ,{N_{\textrm{ele}}}} )$ ). (b) Optimized design 2 (optimized from initial Configuration 1 with uniform density distribution ${\rho _j} = 0.6({j = 1,2,\ldots ,{N_{\textrm{ele}}}} )$ ). (c) Initial Configuration 3 (central crystal column, radius $r/a = 0.25$ .). (d) Optimized design 3 corresponding to initial Configuration 3. (e) Initial Configuration 4 (central hole, radius $r/a = 0.25$ ). (f) Optimized design 4 corresponding to initial Configuration 4.
Fig. 8.
Fig. 8. Optimized 3 × 3 unit cells and their band diagrams under the TE mode. Figures (a) to (j) correspond to the first to the tenth band. Label TE mode denotes the TE polarization of the electromagnetic field.
Fig. 9.
Fig. 9. Volume fraction of GaAs vs. band number in the optimized designs. Material consumption under TM mode is expressed with blue solid line and solid circle, and is expressed with red dashed line and solid circle under TE mode.
Fig. 10.
Fig. 10. Optimized 3 × 3 unit cells photonic crystals with a complete band gap containing: (a) the 2th order TE and 3th order TM; (b) the 3th order TE and 4th order TM; (c) the 4th order TE and 8th order TM; (d) the 5th order TE and 9th order TM; (e) the 6th order TE and 9th order TM; (f) the 7th order TE and 12th order TM; (g) the 8th order TE and 12th order TM; (h) the 9th order TE and 12th order TM.
Fig. 11.
Fig. 11. Two photonic crystal slabs obtained by extruding the optimized configurations in Fig. 4(a) and Fig. 8(a), respectively. (a) and (b) are their corresponding band gap diagrams. The gray line is the light line, and the blue/red dashed lines represent the band diagrams under TM/TE-like modes, respectively. The gaps of the rod/hole slabs are represented by the cyan region.

Tables (2)

Tables Icon

Table 1. Comparison of the relative band gap solutions of PhCs under TM mode.

Tables Icon

Table 2. Comparison of the relative band gap solutions of PhCs under TE mode

Equations (24)

Equations on this page are rendered with MathJax. Learn more.

× × E ( r ) = ε ( r ) ( ω c ) 2 E ( r )
× ( 1 ε ( r ) × H ( r ) ) = ( ω c ) 2 H ( r )  
( E ( r ) ) = ε ( r ) ( ω c ) 2 E ( r ) , for TM mode
( 1 ε ( r ) H ( r ) ) = ( ω c ) 2 H ( r ) , for TE mode  
( + i k ) ( ( + i k ) E k ( r ) ) = ε ( r ) ( ω c ) 2 E k ( r )  for TM mode ,
( + i k ) ( 1 ε ( r ) ( + i k ) H k ( r ) ) = ( ω c ) 2 H k ( r )  for TE mode ,
[ K ω 2 c 2 M   ] u = 0
max   f = ω top ω bot ( ω top + ω bot ) / 2 ,
ω top = min ( ω m + 1 TM ( k ) ) ω bot = max ( ω m TM ( k ) )
ω top = min ( ω n + 1 TE ( k ) ) ,   ω bot = max ( ω n TE ( k ) )
max   f = max { ω n , 1 top ω n , 1 bot ( ω n , 1 top + ω n , 1 bot ) / 2 , ω n , 2 top ω n , 2 bot ( ω n , 2 top + ω n , 2 bot ) / 2 , , ω n , m top ω n , m bot ( ω n , m top + ω n , m bot ) / 2 , } , ω n , m top = 1 γ ln { exp [ γ min ( ω n + 1 TE ( k ) ) ]  +  exp [ γ min ( ω m + 1 TM ( k ) ) ] }  +  1 γ ln ( 2 ) , ω n , m bot = 1 γ ln { exp [ γ max ( ω n TE ( k ) ) ]  +  exp [ γ max ( ω m TM ( k ) ) ] } 1 γ ln ( 2 ) ,
max χ ( x ) f ( χ ) s .t . [ K ( k ) ( ω c ) 2 M   ] u = 0 ,   k [ Γ X M Γ ] χ ( x )  =  0   o r   1 , x Ω des
χ ( x )  =  { 0, if  1 φ ( x ) 0 ,  occupied by material with  ε 1 1, if  0 < φ ( x ) 1 ,   occupied by material with  ε 2
C ( x i , x j ) = exp ( min ( x i x j 2 , x i x ~ j 2 ) l c 2 )
φ ( x , η ) = k = 1 N P ( η k 1 λ k ψ k T C d ( x ) ) ,   x Ω d e s
k = 1 M λ k ( 1 ε ) ( k = 1 N P λ k )
φ ( x ) k = 1 M η k ψ k T C d ( x ) λ k = η T Λ 1 / 2 Φ T C d ( x )
η T W i η 1 ,   ( i = 1 , 2 , , N p )
φ ~ ( x ) = { φ ( x ) e β + 1 e β φ ( x ) , if  φ ( x ) 0 φ ( x ) e β 1 + e β φ ( x ) , if  φ ( x ) < 0
ρ j = 1 + φ ~ ( x j ) 2 ,  ( j = 1 , 2 , , N ele ),
ε ( ρ j ) = ε 1 + ( ε 2 ε 1 ) ρ j ,  for TM mode
1 ε ( ρ j ) = 1 ε 1 + ( 1 ε 2 1 ε 1 ) ρ j , for TE mode
Find η = { η 1 , η 2 , , η M } T max η f ( η ) s .t . [ K ( k ) ( ω c ) 2 M   ] u = 0 ,   k [ Γ X M Γ ] η T W i η 1 ,   ( i = 1 , 2 , , N p )
min   η   f o b j ( η ) = f ( η ) + p 0 max   i ( ( η T W i η 1 ) , 0 )
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.