*Article* **Quantitative Analysis of Photon Density of States for One-Dimensional Photonic Crystals in a Rectangular Waveguide**

#### **Ruei-Fu Jao <sup>1</sup> and Ming-Chieh Lin 2,\***


Received: 1 August 2019; Accepted: 30 October 2019; Published: 4 November 2019

**Abstract:** Light propagation in one-dimensional (1D) photonic crystals (PCs) enclosed in a rectangular waveguide is investigated in order to achieve a complete photonic band gap (PBG) while avoiding the difficulty in fabricating 3D PCs. This work complements our two previous articles (Phys. Rev. E) that quantitatively analyzed omnidirectional light propagation in 1D and 2D PCs, respectively, both showing that a complete PBG cannot exist if an evanescent wave propagation is involved. Here, we present a quantitative analysis of the transmission functions, the band structures, and the photon density of states (PDOS) for both the transverse electric (TE) and transverse magnetic (TM) polarization modes of the periodic multilayer heterostructure confined in a rectangular waveguide. The PDOS of the quasi-1D photonic crystal for both the TE and TM modes are obtained, respectively. It is demonstrated that a "complete PBG" can be obtained for some frequency ranges and categorized into three types: (1) below the cutoff frequency of the fundamental TE mode, (2) within the PBG of the fundamental TE mode but below the cutoff frequency of the next higher order mode, and (3) within an overlap of the PBGs of either TE modes, TM modes, or both. These results are of general importance and relevance to the dipole radiation or spontaneous emission by an atom in quasi-1D periodic structures and may have applications in future photonic quantum technologies.

**Keywords:** photonic crystals; photonic band gap; waveguide; complete PBG; PDOS; TE; TM

#### **1. Introduction**

Photonic crystals (PCs), also known as artificial materials, have attracted much attention in the past three decades due to the tremendous needs of gaining complete control over light propagation and emission [1,2]. PCs, according to the dimension of the periodicity, are divided into three categories: one- (1D), two- (2D), and three-dimensional (3D) crystals. Due to the periodicity, the stop band and the pass band are formed, according to the Bloch theorem [3]. Therefore, periodic dielectric materials are characterized by photonic band gaps (PBGs). As an analogy to the electronic band gaps in solid state materials [4], a PBG in PCs can prohibit the propagation of electromagnetic (EM) waves whose frequencies fall within the band gap region. The continuing success of the semiconductor industry in controlling electric properties of materials from the last century has encouraged us to manipulate the flow of light in PCs and control their optical properties. These PC-based materials are expected to have many applications in optoelectronics, optical communications, and photonic quantum technologies in the next decades [5]. In optical range, PCs have been extensively studied. It was proposed that the emission of EM radiation can be modified by the environment [6,7]. Several environments such as metallic cavities [8], dielectric cavities [9], superlattices [10–15], and 2D PCs [16] have been studied. The environmental effects have been described by the photon density of states (PDOS), which is related to the transition rate of the Fermi golden rule. In principle, a PC with a complete PBG can be best realized in a 3D system to prohibit the propagation of electromagnetic waves of any polarization traveling in any direction from any source [1]. However, the difficulty in fabricating such 3D crystals with PBGs in the optical regime prohibits the progress of many applications.

On the other hand, there has been a lot of interest in microwave and millimeter wave applications of PBG, such as the significant progress in the design of filters [17,18], microstrip antennas [19], and slow wave structures [20,21], and so on. However, the design of PBG in this frequency range is still difficult due to complexities of the modeling. There are too many parameters affecting the PBG properties, such as the number of lattice periods [22], lattice shapes [23,24], lattice spacing [25], and relative volume fraction [26–30]. Since the actual fabrication of 3D PCs remains difficult, another simpler choice is periodic dielectric or PC waveguides, which have only a one-dimensional periodic pattern [1]. The rigorous study of the PC waveguide can be traced back to the 1970s, where a more detailed review can be found [31,32]. Recently, it was demonstrated that, by considering a quantum dot spin coupled to a PC waveguide mode, the light–matter interaction can be asymmetric, leading to unidirectional emission and a deterministic entangled photon source, which might have application in future optical quantum devices [5,33]. One interesting feature of electromagnetism in dielectric media is that there is no fundamental length scale, namely the scaling properties of Maxwell's equations, i.e., the solution of problem at one length scale determines the solutions at all other length scale [1]. In a previous work [34], a multilayer dielectric window in a rectangular waveguide had been studied to achieve a wide-bandwidth transmission of a millimeter wave. A transfer matrix approach was successfully employed to discretize the dielectric function profile and the transmission functions could be calculated efficiently. In principle, the approach can be extended to study a quasi-1D PC, a PC confined in a waveguide. However, the transmission method is limited to study radiation modes in a finite-length system. In order to study the PBG phenomena such as the suppression of spontaneous emission [35] in a quasi-1D PC, the calculation of the dispersion relations or band structures (BS) and the PDOS are needed. Metallic waveguides and cavities are widely used to control microwave propagation. One of the main concerns is visible light energy is quickly dissipated within the metallic components, which makes this method of optical control almost impossible to generalize. Recently, an unconventional superconductivity in magic-angle graphene superlattices had been discovered and studied [36]. The superconductivity might help realize the metallic waveguide confinement of optics in the near future.

In this work, light propagation in 1D PCs enclosed in a rectangular waveguide or quasi-1D PCs is investigated in order to achieve a complete PBG while avoiding the difficulty in fabricating 3D PCs. This work complements two previous articles [15,16] that quantitatively analyzed omnidirectional light propagation in 1D and 2D PCs, respectively, both showing that a complete PBG cannot exist if an evanescent wave propagation is involved. The transfer matrix method is extended to study the transmittance of the quasi-1D PCs for both TE and TM polarization modes [34]. The corresponding BS are obtained by solving the eigenvalue equations with proper periodic boundary conditions following the Bloch theorem [3,4]. The formulas for evaluating the PDOS of the quasi-1D PCs for TE and TM modes are derived, respectively, for determining the PBGs. The contributions of the PDOS from each modes can be distinguished. The model is formulated in Section 2. The calculated results and discussion are presented in Section 3. The conclusions are given in Section 4.

#### **2. Formulations**

A transfer matrix approach is employed to discretize the dielectric function profile of the dielectric multilayer heterostructures and the transmission functions are calculated by matching the boundary conditions at each interfaces. In order to solve the PDOS, it is necessary to calculate the dispersion relation, and the corresponding band structures are obtained by solving the eigenvalue equations with proper periodic boundary conditions.

#### *2.1. Transfer Matrix Method*

Let us consider a waveguide with its rectangular cross section of sides *a* and *b*, and the enclosed multilayer dielectric slab with thickness, (*t*1, *t*2, *t*1, *t*2, ...) and dielectric function, (*ε*1,*ε*2,*ε*1,*ε*2, ...), as shown in Figure 1. The TE mode (H mode) and TM mode (E mode) in the rectangular waveguide are characterized by the *z* components of the magnetic field and the electric field, *Hz* and *Ez*, respectively. By definition, these components are never absent in the corresponding modes. The *z* components of the Helmholtz's equations for the inhomogeneous media are

$$\left\{\varepsilon(z)\vec{\nabla}\times\left[\frac{1}{\varepsilon(z)}\vec{\nabla}\times\vec{H}\right]\right\}\_z+\omega^2\varepsilon(z)\mu(z)H\_z=0\tag{1}$$

and

$$\left\{\mu(z)\nabla\times\left[\frac{1}{\mu(z)}\nabla\times\mathbb{E}\right]\right\}\_z + \omega^2\varepsilon(z)\mu(z)E\_z = 0. \tag{2}$$

**Figure 1.** (**a**) 3D schematic of periodic multilayer heterostructure along the *z*-direction confined in a rectangular waveguide with a width *a* and a height *b* and (**b**) the corresponding dielectric function profile *ε*(*z*) in the cross-sectional view.

In these cases, the effect of losses of the medium inside the waveguide is characterized by the complex permittivity *ε*(*z*) and permeability *μ*(*z*). Thus, Equations (1) and (2) can be rearranged as:

$$\left[\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \varepsilon(z)\frac{\partial}{\partial z}\frac{1}{\varepsilon(z)}\frac{\partial}{\partial z}\right]H\_z + \omega^2\varepsilon(z)\mu(z)H\_z = 0\tag{3}$$

and

$$
\left[\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \mu(z)\frac{\partial}{\partial z}\frac{1}{\mu(z)}\frac{\partial}{\partial z}\right]E\_z + \omega^2\varepsilon(z)\mu(z)E\_z = 0. \tag{4}
$$

By symmetry, using separation of variables, Equations (3) and (4) can be split into transverse and longitudinal parts, and the problem can be simplified as solving the one-dimensional Helmholtz's equation along the *z* direction, the longitudinal parts,

$$
\varepsilon(z)\frac{\partial}{\partial z}\frac{1}{\varepsilon(z)}\frac{\partial}{\partial z}\psi(z)\_{TE} + \left[\omega^2\varepsilon(z)\mu(z) - k\_c^2\right]\psi(z)\_{TE} = 0\tag{5}
$$

and

$$
\mu(z)\frac{\partial}{\partial z}\frac{1}{\mu(z)}\frac{\partial}{\partial z}\psi(z)\_{TM} + \left[\omega^2\varepsilon(z)\mu(z) - k\_c^2\right]\psi(z)\_{TM} = 0,\tag{6}
$$

with eigenvalues *kc* which are determined by the following eigenvalue equations, the transverse parts,

$$\left(\frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + k\_\varepsilon^2\right)\phi(x, y)\_{TE} = 0\tag{7}$$

and

$$
\left(\frac{\partial^2}{\partial \mathbf{x}^2} + \frac{\partial^2}{\partial y^2} + k\_c^2\right) \phi(\mathbf{x}, y)\_{TM} = 0. \tag{8}
$$

The corresponding boundary conditions for *φ*(*x*, *y*)*TE* and *φ*(*x*, *y*)*TM* are

$$\frac{\partial}{\partial x}\phi(x,y)\_{TE}|\_{x=0,a} = 0 \quad \text{and} \quad \frac{\partial}{\partial y}\phi(x,y)\_{TE}|\_{y=0,b} = 0,\tag{9}$$

and

$$
\phi(\mathbf{x}, y)\_{TM}|\_{\mathbf{x}=0, \mathbf{a}} = 0 \quad \text{and} \quad \phi(\mathbf{x}, y)\_{TM}|\_{y=0, b} = 0. \tag{10}
$$

It then follows that

$$k\_c = \frac{2\pi}{\lambda\_c} = \sqrt{\left(\frac{m\pi}{a}\right)^2 + \left(\frac{n\pi}{b}\right)^2}.\tag{11}$$

Applying the eigenvalue equations, Equations (7) and (8), and the boundary conditions, Equations (9) and (10), to the general solutions of Equations (1) and (2), the following particular solutions can be found by separation of variables,

$$H\_z(x, y, z) = H\_0 \cos\left(\frac{m\pi x}{a}\right) \cos\left(\frac{n\pi y}{b}\right) \,\psi(z)\_{TE} \tag{12}$$

and

$$E\_z(x, y, z) = E\_0 \sin\left(\frac{m\pi x}{a}\right) \sin\left(\frac{n\pi y}{b}\right) \psi(z)\_{TM} \tag{13}$$

where *H*<sup>0</sup> and *E*<sup>0</sup> are determined by the energy of electromagnetic waves propagating inside the waveguide, and *m* and *n* are integers. The function *ψ*(*z*) is chosen for a particular solution since it represents propagating waves in the *z*-direction. A transfer matrix approach is employed to discretize the dielectric function profile of the heterostructure. For an *N*-layer dielectric-filled waveguide, *ε*(*z*) and *μ*(*z*) can be divided into *p* = 1, 2, ...,(*N* + 2) layers with a piecewise constant permittivity *ε <sup>p</sup>* and constant permeability *μp*, respectively. The discretized one-dimensional Helmholtz's equation, for the *p*th region with constant permittivity *ε <sup>p</sup>* and constant permeability *μ<sup>p</sup>* can be written as

$$\frac{d^2}{dz^2}\psi\_p(z) + k\_p^2 \psi\_p(z) = 0 \text{ for } z\_{p-1} \le z \le z\_p \tag{14}$$

with

$$k\_p = \sqrt{\omega^2 \varepsilon\_p \mu\_p - k\_c^2} = \frac{2\pi}{\lambda} \sqrt{\varepsilon\_{rp} \mu\_{rp} - \left(\frac{\lambda}{\lambda\_c}\right)^2},\tag{15}$$

where *ψp*(*z*) represents the wave function in the *p*th layer, and *kp* defines the complex wavevector in the same layer along the *z*-direction, *λ* represents the wavelength in free space at the operating angular frequency *ω*, *εrp* and *μrp* are the relative dielectric constant and permeability of the medium, respectively, and *λ<sup>c</sup>* is the cutoff wavelength. The solutions of Equations (12) and (13) can be written as a superposition of the forward and backward traveling wave functions:

$$
\psi\_p = a\_p \exp(-jk\_p z) + b\_p \exp(jk\_p z) \quad \text{for} \quad z\_{p-1} \le z \le z\_p. \tag{16}
$$

*Crystals* **2019**, *9*, 576

The boundary conditions for *ψ*(*z*) at the interface between layers *p* and (*p* + 1) at position *z* = *zp* where *p* = 1, 2, ...,(*N* + 1) are (for TE modes)

$$
\mu\_p \psi\_p(z\_p) = \mu\_{p+1} \psi\_{p+1}(z\_p) \quad \text{and} \quad \frac{d}{dz} \psi\_p(z\_p) = \frac{d}{dz} \psi\_{p+1}(z\_p), \tag{17}
$$

and (for TM modes)

$$
\varepsilon\_p \psi\_p(z\_p) = \varepsilon\_{p+1} \psi\_{p+1}(z\_p) \quad \text{and} \quad \frac{d}{dz} \psi\_p(z\_p) = \frac{d}{dz} \psi\_{p+1}(z\_p). \tag{18}
$$

By matching the boundary conditions at each discontinuity, we arrive at

$$
\sigma \begin{pmatrix} a\_{N+2} \\ b\_{N+2} \end{pmatrix} = M\_{N+1} \cdot \cdots \cdot M\_{\mathcal{P}} \cdot \cdots \cdot M\_1 \begin{pmatrix} a\_1 \\ b\_1 \end{pmatrix} = \begin{pmatrix} M\_{11} & M\_{12} \\ M\_{21} & M\_{12} \end{pmatrix} \begin{pmatrix} a\_1 \\ b\_1 \end{pmatrix},\tag{19}
$$

where for TE modes:

$$\begin{split} \boldsymbol{M}\_{p} = \frac{1}{2} \begin{bmatrix} \exp(jk\_{p+1}z\_{p}) & 0\\ 0 & \exp(-jk\_{p+1}z\_{p}) \end{bmatrix} \cdot \begin{bmatrix} \frac{\mu\_{p}}{\mu\_{p+1}} + \frac{k\_{p}}{k\_{p+1}} & \frac{\mu\_{p}}{\mu\_{p+1}} - \frac{k\_{p}}{k\_{p+1}}\\ \frac{\mu\_{p}}{\mu\_{p+1}} - \frac{k\_{p}}{k\_{p+1}} & \frac{\mu\_{p}}{\mu\_{p+1}} + \frac{k\_{p}}{k\_{p+1}} \end{bmatrix} \cdot \\ & \begin{bmatrix} \exp(-jk\_{p+1}z\_{p}) & 0\\ 0 & \exp(jk\_{p+1}z\_{p}) \end{bmatrix} \,' \end{split} \tag{20}$$

and for TM modes:

$$M\_p = \frac{1}{2} \begin{bmatrix} \exp(jk\_{p+1}z\_p) & 0\\ 0 & \exp(-jk\_{p+1}z\_p) \end{bmatrix} \cdot \begin{bmatrix} \frac{\varepsilon\_p}{\varepsilon\_{p+1}} + \frac{k\_p}{k\_{p+1}} & \frac{\varepsilon\_p}{\varepsilon\_{p+1}} - \frac{k\_p}{k\_{p+1}}\\ \frac{\varepsilon\_p}{\varepsilon\_{p+1}} - \frac{k\_p}{k\_{p+1}} & \frac{\varepsilon\_p}{\varepsilon\_{p+1}} + \frac{k\_p}{k\_{p+1}} \end{bmatrix} \cdot \begin{bmatrix} \frac{\varepsilon\_p}{\varepsilon\_{p+1}} - \frac{k\_p}{k\_{p+1}}\\ 0 \end{bmatrix} \cdot \begin{bmatrix} \frac{\varepsilon\_p}{\varepsilon\_{p+1}} - \frac{k\_p}{k\_{p+1}} & 0\\ 0 & \exp(jk\_{p+1}z\_p) \end{bmatrix} \cdot \tag{21}$$

Using Equation (19) with *a*<sup>1</sup> = 1, *b*<sup>1</sup> = *r*, *aN*+<sup>2</sup> = *t*, and *bN*+<sup>2</sup> = 0, the reflection and transmission amplitudes, *r* and *t*, can be obtained, respectively, by

$$r = -\frac{M\_{21}}{M\_{22}}\tag{22}$$

and

$$t = \frac{M\_{11} \cdot M\_{22} - M\_{12} \cdot M\_{21}}{M\_{22}}.\tag{23}$$

The reflection and transmission coefficients, *R* and *T*, can be implicitly represented by *S*-parameters, *S*11(*dB*) and *S*12(*dB*), as a function of the operating frequency, , for TE modes

$$S\_{11}(dB) = 10 \log \left| \frac{M\_{21}}{M\_{22}} \right|^2 \quad \text{and} \quad S\_{12}(dB) = 10 \log \frac{\mu\_1 k\_{M+2}}{\mu\_{M+2} k\_1} \left| \frac{M\_{11} \cdot M\_{22} - M\_{12} \cdot M\_{21}}{M\_{22}} \right|^2,\tag{24}$$

and for TM modes

$$S\_{11}(dB) = 10 \log \left| \frac{M\_{21}}{M\_{22}} \right|^2 \quad \text{and} \quad S\_{12}(dB) = 10 \log \frac{\varepsilon\_1 k\_{N+2}}{\varepsilon\_{N+2} k\_1} \left| \frac{M\_{11} \cdot M\_{22} - M\_{12} \cdot M\_{21}}{M\_{22}} \right|^2. \tag{25}$$

#### *2.2. Band Structures*

In this model, *Ez* and *Hz* both are periodic with period Λ. According to the Bloch theorem, the electric and magnetic fields in a periodic layered medium are *Ez*(*z*) = *Ez*(*z* + Λ) and *Hz*(*z*) = *Hz*(*z* + Λ), respectively. The column vector of the Bloch wave satisfies the following eigenvalue equation for consistency

$$
\begin{pmatrix} a\_3 \\ b\_3 \end{pmatrix} = \frac{1}{4} \begin{pmatrix} A\_{TE/TM} & B\_{TE/TM} \\ C\_{TE/TM} & D\_{TE/TM} \end{pmatrix} \cdot \begin{pmatrix} a\_1 \\ b\_1 \end{pmatrix} = \epsilon^{\vec{\mathbf{x}}\_B \Lambda} \begin{pmatrix} a\_1 \\ b\_1 \end{pmatrix} . \tag{26}
$$

For the TE mode:

$$A\_{TE} = e^{j(k\_1t\_1 - k\_2t\_2)} \left(\frac{\mu\_1}{\mu\_2} - \frac{k\_1}{k\_2}\right) \left(\frac{\mu\_2}{\mu\_1} - \frac{k\_2}{k\_1}\right) + \\\\ \begin{split} \left(\frac{\mu\_1}{\mu\_1} - \frac{k\_1}{k\_2}\right) + \\ \left(\frac{\mu\_1}{\mu\_2} + \frac{k\_1}{k\_2}\right) \left(\frac{\mu\_2}{\mu\_1} + \frac{k\_2}{k\_1}\right) \end{split} \tag{27}$$

$$B\_{TE} = e^{-j[k\_1t\_1 - k\_2(2t\_1 + t\_2)]} \left(\frac{\mu\_1}{\mu\_2} + \frac{k\_1}{k\_2}\right) \left(\frac{\mu\_2}{\mu\_1} - \frac{k\_2}{k\_1}\right) + \\\\B\_{\frac{1}{\mu\_2}} = \frac{k\_1}{k\_2} \left(\frac{\mu\_1}{\mu\_1} - \frac{k\_1}{k\_2}\right) \left(\frac{\mu\_2}{\mu\_1} + \frac{k\_2}{k\_1}\right),\tag{28}$$

$$C\_{TE} = B\_{TE\prime}^{\*}\tag{29}$$

and

$$D\_{TE} = A\_{TE}^{\*}.\tag{30}$$

For the TM mode:

$$A\_{TM} = e^{j(k\_1 t\_1 - k\_2 t\_2)} \left(\frac{\varepsilon\_1}{\varepsilon\_2} - \frac{k\_1}{k\_2}\right) \left(\frac{\varepsilon\_2}{\varepsilon\_1} - \frac{k\_2}{k\_1}\right) + \\\\ \begin{split} \left(\frac{\varepsilon\_1}{\varepsilon\_1} - \frac{k\_1}{k\_2}\right) + \\ \left(\frac{\varepsilon\_1}{\varepsilon\_2} + \frac{k\_1}{k\_2}\right) \left(\frac{\varepsilon\_2}{\varepsilon\_1} + \frac{k\_2}{k\_1}\right) \end{split} \tag{31}$$

$$B\_{TM} = e^{-j[k\_1 t\_1 - k\_2(2t\_1 + t\_2)]} \left(\frac{\varepsilon\_1}{\varepsilon\_2} + \frac{k\_1}{k\_2}\right) \left(\frac{\varepsilon\_2}{\varepsilon\_1} - \frac{k\_2}{k\_1}\right) + \\\\ \left. \frac{1}{e^{j[k\_1 t\_1 + k\_2(2t\_1 + t\_2)]}} \left(\frac{\varepsilon\_1}{\varepsilon\_2} - \frac{k\_1}{k\_2}\right) \left(\frac{\varepsilon\_2}{\varepsilon\_1} + \frac{k\_2}{k\_1}\right) \right) \tag{32}$$

$$C\_{TM} = B\_{TM\prime}^\* \tag{33}$$

and

$$D\_{TM} = A\_{TM}^\*.\tag{34}$$

The phase factor *ejkB*<sup>Λ</sup> is thus the eigenvalue and satisfies the secular equation

$$
\begin{vmatrix} A\_{TE/TM} - \mathfrak{e}^{j k\_B \Lambda} & B\_{TE/TM} \\ \mathfrak{C}\_{TE/TM} & D\_{TE/TM} - \mathfrak{e}^{j k\_B \Lambda} \end{vmatrix} = 0. \tag{35}
$$

Finally, the dispersion relation for the Bloch wave function is

$$k\_B(k\_{z\prime}\omega) = \frac{1}{t\_1 + t\_2} \cos^{-1}[\frac{1}{2}(A\_{TE/TM} + D\_{TE/TM})].\tag{36}$$

#### *2.3. Photon Density of States*

The quasi-one-dimensional photonic crystal has been confined in the *xy*-plane, so the wave vectors *kx* and *ky* are determined according to the guiding modes. To perform the PDOS calculation, it is required to use the formal definition which is the number of available photon modes per unit frequency range. Then, we construct two frequencies, namely, *ω*(*kB*) = *ω* and *ω*(*kB*) = *ω* + Δ*ω*, where Δ*ω* is a small increment. We calculate the line therein, and divide it by the line segment occupied by a single mode. The differential line element in *K* space within the frequency range, is given by Δ*Lk* = Δ*kB*. Now, Δ*kB* is defined as

$$
\Delta k\_B = \frac{\Delta \omega}{|\nabla\_k \omega|}. \tag{37}
$$

Integrating over the frequency increment, we have that the total phase space line segment contributing to the frequency range (*ω*, *ω* + *dω*) is

$$\int\_{\omega\_k} dL\_k = \int\_{\omega\_k} \frac{1}{|\nabla\_k \omega|} d\omega\_\prime \tag{38}$$

where we take the limit of infinitesimal increments. The number of modes within the range (*ω*, *ω* + *dω*) is obtained by dividing the length calculated in Equation (38) by the line segment corresponding to one mode, 2*π*/Λ in the phase space. This yields

$$dN(\omega) = \frac{\Lambda}{2\pi} \int\_{\omega\_k} \frac{1}{|\nabla\_k \omega|} d\omega \equiv D(\omega) d\omega. \tag{39}$$

Because *ω* is a function of *k* , we can write

$$
\nabla\_k \omega = \frac{d\omega}{dk\_B} \mathbf{\hat{z}}.\tag{40}
$$

For the TE mode:

$$\nabla\_k \omega\_{TE} = -\frac{(t\_1 + t\_2)\sin\left[0.5k\_B(t\_1 + t\_2)\right]}{a\_1 + a\_2 + a\_3 + a\_4 + a\_5 + a\_6}\mathbf{\hat{z}}\tag{41}$$

and for the TM mode:

$$\nabla\_k \omega\_{TM} = -\frac{(t\_1 + t\_2)\sin[0.5k\_B(t\_1 + t\_2)]}{\beta\_1 + \beta\_2 + \beta\_3 + \beta\_4 + \beta\_5 + \beta\_6} \mathbf{\hat{z}},\tag{42}$$

where the functions, *α<sup>i</sup>* and *βi*, with *i* = 1, 2, ..., 6, are listed in the Appendix A. Finally, the formula for evaluating the PDOS can be expressed as

$$D(\omega)\_{TE/TM} = \frac{\Lambda}{2\pi} \int\_{\omega\_k} \frac{1}{\left| \nabla\_k \omega\_{TE/TM} \right|} d\omega. \tag{43}$$

#### **3. Results and Discussion**

All of macroscopic electromagnetism, including the propagation of light in a photonic crystal, is governed by the four macroscopic Maxwell's equations with no free charges or currents. One interesting feature of electromagnetism in PCs is that there is no fundamental length scale other than the assumption that the system is macroscopic [1]. Therefore, to study physical phenomena in PCs, one may scale a system from the optical frequency range to the microwave one and vice versa if suitable conditions are fulfilled. For the purposes of demonstration and easier verification by experimentalists, a WR28 (7.11 mm × 3.555 mm) rectangular waveguide, usually used for Ka-band millimeter waves, is chosen. The periodic dielectric heterostructure is arranged along the longitudinal (*z*) direction to form the quasi-1D PCs [34]. In the microwave or millimeter wave frequency ranges, the PC experiments are very popular and more affordable compared to those in the optical frequency ones. Nevertheless, in order to extend the results for general use, the data are normalized for solutions

at all length scale. Consider the quasi-1D PCs with the arrangement shown in Figure 1, 15 double-layer stacks (*n* = 30) have been included in the waveguide for investigation. The transmittance is calculated using the transfer matrix approach mentioned above. In order to demonstrate the stop band and pass band, the transmittances, expressed as the *S*-parameters, *S*12, in dB, of the lowest TE and TM modes for three quasi-1D PCs have been calculated and plotted in Figure 2. For the three cases, the dielectric constants used are the same as *ε*<sup>1</sup> = 3.8 (quartz) and *ε*<sup>2</sup> = 1.0 (air), while the thicknesses are varied as (*t*1, *t*2) = (1.00, 3.30), (1.00, 3.60), and (1.00, 3.90) mm, corresponding to the filling ratios *t*1/Λ = 23.26%, 21.74%, and 20.41% for the periods of the stacks Λ = 4.3, 4.6, and 4.9 mm, respectively. As one can see, the central frequencies of the PBGs for both the *TE*<sup>10</sup> and *TM*<sup>11</sup> modes, as shown in Figure 2a,b respectively, are shifted to lower values as the filling ratio decreases. As mentioned above, the frequency axes in the plots are normalized to the cutoff frequency of *TE*<sup>10</sup> for general use.

**Figure 2.** Calculated transmittance, *S*12(dB), of the quasi-1D PCs with 15 double-layer stacks enclosed in the rectangular waveguide for (**a**) *TE*<sup>10</sup> and (**b**) *TM*<sup>11</sup> modes of light propagation. The dielectric constants used are *ε*<sup>1</sup> = 3.8 and *ε*<sup>2</sup> = 1.0, and the corresponding filling ratios are *t*1/Λ = 23.26%, 21.74%, and 20.41% in the three cases. The frequency axes are normalized to the cutoff frequency of *TE*10.

In order to look for a complete PBG for all possible modes in the quasi-1D PCs, one may want to enlarge a PBG for a specific mode of choice. For the following four cases, the dielectric constants used are (*ε*1, *ε*2) = (2.3, 1.0), (3.8, 1.0), (4.9, 1.0), and (11.4, 1.0), while the thicknesses are varied as (*t*1, *t*2) = (2.15, 2.15), (1.00, 3.30), (0.72, 3.58), and (0.27, 4.03) mm, corresponding to the filling ratios *t*1/Λ = 50.00%, 23.26%, 16.74%, and 6.28% for the fixed period of the stacks Λ = 4.3 mm, respectively. The dielectric constants of 2.3, 3.8, 4.9 and 11.4 used in the chosen microwave frequency range correspond to the dielectric materials, polyethylene, quartz, phenolic resin, and barium sulfate, respectively. As one can see, the calculated transmittances, *S*12(dB), of the four quasi-1D PCs with different 15 double-layer stacks for the *TE*<sup>10</sup> and *TM*<sup>11</sup> modes are plotted in Figure 3a,b, respectively. The width of PBGs for both the *TE*<sup>10</sup> and *TM*<sup>11</sup> modes in the quasi-1D PCs is widened to larger sizes as decreasing the filling ratio while increasing the dielectric contrast between *ε*<sup>1</sup> and *ε*2. One may notice that the configurations of these quasi-1D PCs have been specially arranged so that the central frequencies of the PBGs of the *TE*<sup>10</sup> mode are kept the same while tuning the PBG sizes. However, those of the PBGs of the *TM*<sup>11</sup> mode are not centralized. This indicates that a tremendous computational effort is still needed to find a complete PBG in spite of the approach with a closed form developed is very efficient for the quasi-1D PCs.

**Figure 3.** Calculated transmittance, *S*12(dB), of the quasi-1D PCs with 15 double-layer stacks enclosed in the rectangular waveguide for (**a**) *TE*<sup>10</sup> and (**b**) *TM*<sup>11</sup> modes of light propagation. The dielectric constants used are *ε*<sup>1</sup> = 2.3, 3.8, 4.9, and 11.4 and *ε*<sup>2</sup> = 1, and the corresponding filling ratios are *t*1/Λ = 50.00%, 23.26%, 16.74%, and 6.28% in the four cases. The frequency axes are normalized to the cutoff frequency of *TE*10.

A better way to identify PBGs more accurately and efficiently is to perform the BS calculation and further compute the PDOS [15,16]. Figure 4 shows the band structures of (a) TE10 and (b) TM11 modes for the quasi-1D PC with the dielectric constants *ε*<sup>1</sup> = 3.8, *ε*<sup>2</sup> = 1 and the thicknesses *t*<sup>1</sup> = 1.00 mm, *t*<sup>2</sup> = 3.30 mm, corresponding to the case presented in Figures 2 and 3 in a black solid line with a filling ratio of 23.26%. The PBGs of the TE10 mode show no overlap with those of TM11 mode, all marked in gray stripes in the figure. However, from Figure 4, one may notice that there are no photon states in the frequency ranges below the cutoff frequency of the *TE*<sup>10</sup> mode nor within the first PBG of *TE*<sup>10</sup> which is under the cutoff frequency of the *TM*<sup>11</sup> mode. Note that only the lowest TE amd TM modes are plotted in Figure 4. As there might be other modes involved within the frequency range of the PBG of interest, it is easier to identify a complete PBG from the PDOS plots compared to the BS ones. Figure 5 shows the PDOS of TE10, TE01, TE11, and TM11 modes for the same quasi-1D PC. As one can see, the PDOS of the TE10 and TM11 modes are consistent with the BS calculations shown in Figure 4. The PDOS contributed from the TE01 and TE11 modes tend to fill up the first PBG of the TE10 mode. Other higher order modes are not considered as their cutoff frequencies are too high to contribute any PDOS in the frequency range of the PBG. Finally, the combined PDOS of TE10, TE01, TE11, and TM11 modes shows no photon states in some frequency ranges. The first one is below the cutoff frequency of the *TE*<sup>10</sup> mode, (0–0.77)f*c*, the second one is within the first PBG of *TE*<sup>10</sup> but under the cutoff frequency of the *TE*<sup>01</sup> mode, (1.26–1.48)f*c*, and the third one is the overlap of the PBGs of TE10, TE01, and TE11 modes, (1.79–1.87)f*c*. Therefore, a "complete PBG" can be obtained for some frequency ranges and categorized into three types: (1) below the cutoff frequency of the fundamental TE mode, (2) within the PBG of the fundamental TE mode but below the cutoff frequency of the next higher order mode, and (3) within an overlap of the PBGs of either TE modes, TM modes, or both.

**Figure 4.** Band structures of (**a**) TE10 and (**b**) TM11 modes for the quasi-1D PC with the dielectric constants *ε*<sup>1</sup> = 3.8, *ε*<sup>2</sup> = 1 and the thicknesses *t*<sup>1</sup> = 1.00 mm, *t*<sup>2</sup> = 3.30 mm, corresponding to a filling ratio of 23.26%. The PBGs of TE10 mode show no overlap with those of TM11 mode.

**Figure 5.** PDOS of TE10, TE01, TE11, and TM11 modes for the quasi-1D PC with the dielectric constants *ε*<sup>1</sup> = 3.8, *ε*<sup>2</sup> = 1 and the thicknesses *t*<sup>1</sup> = 1.00 mm, *t*<sup>2</sup> = 3.30 mm, corresponding to a filling ratio of 23.26%. The combined PDOS of TE10, TE01, TE11, and TM11 modes shows no photon states in some frequency ranges.

#### **4. Conclusions**

In summary, light propagation in quasi-1D PCs has been investigated quantitatively. The transmittances for both the TE and TM modes through a periodic multilayer heterostructure in a rectangular waveguide are calculated using the transfer matrix method. The corresponding band structures are obtained by solving the eigenvalue equations with proper periodic boundary conditions following the Bloch theorem. The formulas for determining the PDOS have been obtained to facilitate identifying the photonic band gaps for all the modes residing in the system. With our approach, the quantitative determination of PDOS can be performed very accurately and efficiently. A complete PBG can exist in these quasi-1D PCs, but the determination must be carefully conducted and verified. It is demonstrated that three types of "complete PBG" can be found and categorized as follows. The first type is the frequency range within which the TE and TM modes are both cutoff, the second type is for which the fundamental TE mode has a PBG while other higher order modes are cutoff, and the last type is an overlap of the PBGs of either TE modes, TM modes, or both. The model might be easier for an experimental validation in a millimeter wave frequency range while the optical counterpart might

be possibly pursued as well. We believe these results are of general importance and relevance to the dipole radiation or spontaneous emission by an atom in quasi-1D periodic structures and may have applications in future photonic quantum technologies.

**Author Contributions:** Conceptualization, M.-C.L.; methodology, R.-F.J. and M.-C.L.; software, R.-F.J. and M.-C.L.; validation, R.-F.J. and M.-C.L.; formal analysis, R.-F.J.; investigation, R.-F.J. and M.-C.L.; resources, R.-F.J. and M.-C.L.; data curation, R.-F.J.; writing—original draft preparation, R.-F.J. and M.-C.L.; writing—review and editing, R.-F.J. and M.-C.L.; visualization, R.-F.J.; supervision, M.-C.L.; project administration, M.-C.L.; funding acquisition, R.-F.J. and M.-C.L.

**Funding:** This work was partially supported by Guangdong Industry Polytechnic, P. R. China, under Grant Nos. KYRC2018-001 and RC2016-005, the research fund of Hanyang University (HY- 201400000002393), National Research Foundation of South Korea (2015R1D1A1A01061017), and the Alexander von Humboldt Foundation of Germany.

**Acknowledgments:** The authors would like to thank the late B. Y. Gu at the Institute of Physics, CAS and C. T. Chan at the Department of Physics, HKUST for the helpful comments, discussions, and encouragement.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **Appendix A. FORMULAS**

Here, we give the formulas of the functions employed in Equations (41) and (42).

$$\begin{split} \boldsymbol{\alpha}\_{1} = \exp\left[-j(k\_{1}t\_{1} + k\_{2}t\_{2})\right] \left[\exp\left(j2k\_{1}t\_{1}\right) + \exp\left(j2k\_{2}t\_{2}\right)\right] \\ \left[\frac{\omega}{k\_{1}}\left(\frac{\varepsilon\_{1}\mu\_{1}k\_{2}}{k\_{1}^{2}} - \frac{\varepsilon\_{2}\mu\_{2}}{k\_{2}}\right)\right] \left(-\frac{k\_{1}}{k\_{2}} + \frac{\mu\_{1}}{\mu\_{2}}\right), \end{split} \tag{A1}$$

$$n\_2 = \exp\left[-j(k\_1t\_1 + k\_2t\_2)\right] \left\{1 + \exp\left[j2(k\_1t\_1 + k\_2t\_2)\right]\right\}$$

$$\left[\frac{\omega}{k\_1}\left(-\frac{\varepsilon\_1\mu\_1k\_2}{k\_1^2} + \frac{\varepsilon\_2\mu\_2}{k\_2}\right)\right] \left(\frac{k\_1}{k\_2} + \frac{\mu\_1}{\mu\_2}\right),\tag{A2}$$

$$\begin{split} a\_{3} = \exp\left[-j(k\_{1}t\_{1} + k\_{2}t\_{2})\right] \left[\exp\left(j2k\_{1}t\_{1}\right) + \exp\left(j2k\_{2}t\_{2}\right)\right] \\ \left[\frac{\omega}{k\_{2}}\left(-\frac{\varepsilon\_{1}\mu\_{1}}{k\_{1}} + \frac{\varepsilon\_{2}\mu\_{2}k\_{1}}{k\_{2}^{2}}\right)\right] \left(-\frac{k\_{2}}{k\_{1}} + \frac{\mu\_{2}}{\mu\_{1}}\right), \end{split} \tag{A3}$$

$$\begin{aligned} \mathfrak{a}\_{4} &= -j \exp\left[-j(k\_{1}t\_{1} + k\_{2}t\_{2})\right] \left[\exp(j2k\_{1}t\_{1}) - \exp(j2k\_{2}t\_{2})\right] \\ &\left[\omega \left(-\frac{\varepsilon\_{1}\mu\_{1}t\_{1}}{k\_{1}} + \frac{\varepsilon\_{2}\mu\_{2}t\_{2}}{k\_{2}}\right)\right] \left(-\frac{k\_{1}}{k\_{2}} + \frac{\mu\_{1}}{\mu\_{2}}\right) \left(\frac{k\_{2}}{k\_{1}} + \frac{\mu\_{2}}{\mu\_{1}}\right), \end{aligned} \tag{A4}$$

$$n\_5 = \exp\left[-j(k\_1t\_1 + k\_2t\_2)\right] \left\{1 + \exp\left[j2(k\_1t\_1 + k\_2t\_2)\right]\right\}$$

$$\left[\frac{\omega}{k\_2}\left(\frac{\varepsilon\_1\mu\_1}{k\_1} - \frac{\varepsilon\_2\mu\_2k\_1}{k\_2^2}\right)\right] \left(\frac{k\_2}{k\_1} + \frac{\mu\_2}{\mu\_1}\right),\tag{A5}$$

$$\begin{split} \mu\_{6} &= -j \exp\left[-j(k\_{1}t\_{1} + k\_{2}t\_{2})\right] \{-1 + \exp\left[j2(k\_{1}t\_{1} + k\_{2}t\_{2})\right] \} \\ &\quad \left[\omega \left(-\frac{\varepsilon\_{1}\mu\_{1}t\_{1}}{k\_{1}} - \frac{\varepsilon\_{2}\mu\_{2}t\_{2}}{k\_{2}}\right)\right] \left(\frac{k\_{1}}{k\_{2}} + \frac{\mu\_{1}}{\mu\_{2}}\right) \left(\frac{k\_{2}}{k\_{1}} + \frac{\mu\_{2}}{\mu\_{1}}\right) \end{split} \tag{A6}$$

$$\beta\_1 = \exp\left[-j(k\_1t\_1 + k\_2t\_2)\right] \left[\exp\left(j2k\_1t\_1\right) + \exp\left(j2k\_2t\_2\right)\right]$$

$$\left[\frac{\omega}{k\_1}\left(\frac{\varepsilon\_1\mu\_1k\_2}{k\_1^2} - \frac{\varepsilon\_2\mu\_2}{k\_2}\right)\right] \left(-\frac{k\_1}{k\_2} + \frac{\varepsilon\_1}{\varepsilon\_2}\right),\tag{A7}$$

$$\beta\_2 = \exp\left[-j(k\_1t\_1 + k\_2t\_2)\right] \{1 + \exp[j2(k\_1t\_1 + k\_2t\_2)]\}$$

$$\left[\frac{\omega}{k\_1}\left(-\frac{\varepsilon\_1\mu\_1k\_2}{k\_1^2} + \frac{\varepsilon\_2\mu\_2}{k\_2}\right)\right] \left(\frac{k\_1}{k\_2} + \frac{\varepsilon\_1}{\varepsilon\_2}\right),\tag{A8}$$

$$\beta\_3 = \exp\left[-j(k\_1t\_1 + k\_2t\_2)\right] \left[\exp\left(j2k\_1t\_1\right) + \exp\left(j2k\_2t\_2\right)\right]$$

$$\left[\frac{\omega}{k\_2}\left(-\frac{\varepsilon\_1\mu\_1}{k\_1} + \frac{\varepsilon\_2\mu\_2k\_1}{k\_2^2}\right)\right] \left(-\frac{k\_2}{k\_1} + \frac{\varepsilon\_2}{\varepsilon\_1}\right),\tag{A9}$$

$$\begin{aligned} \beta\_4 &= -j \exp\left[-j(k\_1 t\_1 + k\_2 t\_2)\right] \left[\exp(j2k\_1 t\_1) - \exp(j2k\_2 t\_2)\right] \\ &\left[\omega \left(-\frac{\varepsilon\_1 \mu\_1 t\_1}{k\_1} + \frac{\varepsilon\_2 \mu\_2 t\_2}{k\_2}\right)\right] \left(-\frac{k\_1}{k\_2} + \frac{\varepsilon\_1}{\varepsilon\_2}\right) \left(\frac{k\_2}{k\_1} + \frac{\varepsilon\_2}{\varepsilon\_1}\right), \end{aligned} \tag{A10}$$

$$\beta\_5 = \exp\left[-j(k\_1t\_1 + k\_2t\_2)\right] \left\{1 + \exp\left[j2(k\_1t\_1 + k\_2t\_2)\right]\right\}$$

$$\left[\frac{\omega}{k\_2}\left(\frac{\varepsilon\_1\mu\_1}{k\_1} - \frac{\varepsilon\_2\mu\_2k\_1}{k\_2^2}\right)\right] \left(\frac{k\_2}{k\_1} + \frac{\varepsilon\_2}{\varepsilon\_1}\right),\tag{A11}$$

and

$$\beta\_{6} = -j \exp\left[-j(k\_{1}t\_{1} + k\_{2}t\_{2})\right] \{-1 + \exp[j2(k\_{1}t\_{1} + k\_{2}t\_{2})] \}$$

$$\left[\omega \left(-\frac{\varepsilon\_{1}\mu\_{1}t\_{1}}{k\_{1}} - \frac{\varepsilon\_{2}\mu\_{2}t\_{2}}{k\_{2}}\right)\right] \left(\frac{k\_{1}}{k\_{2}} + \frac{\varepsilon\_{1}}{\varepsilon\_{2}}\right) \left(\frac{k\_{2}}{k\_{1}} + \frac{\varepsilon\_{2}}{\varepsilon\_{1}}\right). \tag{A12}$$

#### **References**


c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Spatial Beam Filtering with Autocloned Photonic Crystals**

#### **Pei-Yu Wang, Yi-Chen Lai and Yu-Chieh Cheng \***

Department of Electro-Optical Engineering, National Taipei University of Technology, No. 1, Sec. 3, Chung-Hsiao E. Rd, Taipei 10608, Taiwan; peggy831110pp@gmail.com (P.-Y.W.); 26789jimmy@gmail.com (Y.-C.L.)

**\*** Correspondence: yu-chieh.cheng@mail.ntut.edu.tw

Received: 18 October 2019; Accepted: 6 November 2019; Published: 8 November 2019

**Abstract:** We have been numerically demonstrated the mechanism of spatial beam filtering with autocloned photonic crystals. The spatial filtering through different configurations of the multilayered structures based on a harmonically modulated substrate profile is considered. The paper demonstrates a series of parameter studies to look for the best spatial beam filtering performance. The optimization results show that a beam spectral width of 39.2◦ can be reduced to that of 5.92◦, leading to high potential applications for integrated optical microsystems.

**Keywords:** photonic crystal; beam shaping; angular filtering; autocloning; multilayered structures

#### **1. Introduction**

Angular/spatial filtering devices based on photonic crystals (PhCs) [1,2] provide diffraction of the angular components of an incident beam. The effect of a PhCs-based spatial/angular filtering device that works on a spatial frequency spectrum relies on an angular band-gap [3–7]. For spatial filtering, a range of angular components of a beam can be removed due to the angular band-gaps, that is, the waves can be reflected in a backward direction [3–5] or deflected at large angles in a forward one [6,7].

Furthermore, double-periodic photonic structures enable manipulation of the zero diffraction order of a transmitted beam [8]. For example, some angular components of an incident light source diffract from the zero diffraction order to the other orders at resonance conditions. On the other hand, some angular components, out of resonance, directly propagate through the PhCs. In this way, low-angle-pass or high-angle-pass filtering devices are achievable through a proper interplay among the grating characteristics.

In particular, the PhCs filtering has been already implemented for intracavity angular filtering in an integrated platform such as microchip lasers [9]. Such a PhCs-based confocal filtering device presents an alternative method for replacing conventional filtering devices [10], but has a critical disadvantage that is the presence of an optical axis [11]. Therefore, the transmitted axisymmetric concentric ring structures result in the limitation of angular filter merely for on-axis incident light.

It is noted that there have been other approaches to spatial filtering such as passive [12] or light-induced [13] Bragg gratings and pulse-induced population density gratings [14]. However, these alternative methods require not only more sophisticated schemes but additional optical components, leading to limited applications in the compact micro-systems. Therefore, more compact PhCs-based angular filters are desirable solutions for the use in the microlasers, e.g., autocloned PhCs, which preserves the initial modulation during the deposition of multilayers [15]. For example, a multilayered photonic microstructure based on a sinusoidal or braze profile was demonstrated experimentally as one of the most compact PhCs-based angular filters [16]. However, this proposed filter presents the weak filtering performance, so further investigation is required for practical use. For example, a compact filter with a transverse invariance performs both the narrow angular bandwidth and the high efficiency. The filtering for the application of such a compact filter toward a Gaussian beam has not been well studied in the prior study [16]. Although the fabrication of the autocloned PhCs has been demonstrated experimentally [17,18], it is still unavoidable that the variation of the amplitude (Amp) of the harmonic modulation of the autoclaved PhCs increases with the number of layers [19]. Further studies for feasible parameters in fabrication such as less number of layers are a concern.

In this paper, we provide a numerical study of an angular/spatial filtering based on multilayered PhCs gratings with a harmonically modulated substrate profile. The multilayered gratings are all-dielectric and periodic, where the variation of the periods or wavelength results in a low- or high-angle-pass filter. A spectral width (SW), defined as the full spectral width at half maximum (FWHM), is calculated by using the finite-difference time-domain (FDTD) method [20]. To enhance the spatial filtering, we focus on narrowing filtering angular distributions by the design of the low-angle-pass filter and the results are compared with the study in [16].

#### **2. Numerical Far-Field Simulations for Autocloned PhCs**

The configuration of the proposed multilayered PhCs grating for spatial filtering is schematically shown in Figure 1. The multilayers consisting of alternating high-refractive-index (NH) and low-refractive-index (NL) layers with the number of layers (N) on top of a grating substrate (Nsub) that results in a sinusoidal profile with a peak-to-peak Amp. Such a multilayered grating with a sinusoidal profile provides a transversal and a longitudinal modulations of the period (D*x*) and an alternating thicknesses (D*z*) of the multilayers. The thicknesses of the high- and low-refractive-index layers are equal to D*z*/2. The wavelength of the incident light is equal to λ = 582 nm and a transverse electric (TE) polarization is considered. A more detailed calculation procedure is further described in the following section, Method.

**Figure 1.** Illustration of the proposed autocloning photonic crystals (PhCs) with a harmonic modulation. The inset shows the parameters of the wavy structure. Dx and Dz represent the horizontal and longitudinal periods, respectively. N is the number of layers. Low- and high-refractive-index materials are shown in red and blue colors, respectively. The parameter Amp indicates the peak-to-peak value of the multilayer structure, also regarded as the amplitude of the harmonic modulation. The beam width (BW) represents the full spectral width at half maximum (FWHM) of a launched beam. The steady-state near-field plane is defined as the output plane of the simulation domain.

We first analyze diffraction patterns of the periodic multilayered structure by using a rigorous coupled-wave analysis (RCWA) method [21] and its 0th-order transmission is as plotted in Figure 2a. The parameters of the studied structure are identical as those in [16]: N = 33 layers, Dx = 1.67 μm, D*z*=0.24 μm, NSub =1.5, NH =1.42, and NL =1.3. Since the parameter of Amp was not mentioned in [16],

we found the narrowest SW by scanning Amp up to 0.5 μm. The narrowest SW of 9.21◦ is obtained for Amp = 0.27 μm. Furthermore, the angular dependence of the 0th-order diffraction efficiency (DEt0) (Appendix A Method) presents a low-pass filtering design where the angular components at around 10◦ are removed or coupled out to the others, as shown in Figure 2b. Thus, the far-field distribution of a Gaussian beam, regarded as different plane-wave components at different angles of incidence, can be narrowed down. Figure 2c,d show that the filtered far-field distributions for two incident beams with different beam widths (BW1= 1 λ and BW2 = 2 λ). Their filtered SWs, SW1' (9.21◦) and SW2' (9.42◦), are narrower than SW1 (39.2◦) and SW2 (20.99◦), respectively.

**Figure 2.** Numerical far-field simulation of autocloning-mode-based PhCs: (**a**) the diffraction map of 0th order versus angles of incidence; (**b**) diffraction efficiency of 0th-order transmission (DEt0) with respect to angles of incidence; far-field distributions for two different incident Gaussian beam widths: (**c**) BW1 = λ; (**d**) BW2 = 2λ. The parameters of the structure in (**a**) and (**b**) are as following: N = 33 layers, Amp = 0.27 μm, D*z* = 0.24 μm, and D*x* = 1.67 μm. These two spectral width SW1 and SW2 mean the FWHM of the spectral width of the Gaussian beam for BW1 = λ and BW2 = 2λ, respectively. SW1 and SW1' in (**c**) represent the normalized far-field distributions of the Gaussian beam passing without PhCs in blue and with the PhCs in red for BW1 = λ. SW2 and SW2' in (**d**) represent the normalized far-field distributions of the Gaussian beam passing without PhCs in blue and with the PhCs in red for BW2 = 2λ.

#### **3. Minimum SWs for Di**ff**erent Configurations of Autocloned PhCs**

In this section, a series of numerical calculations is executed for the best filtering effect. We focus on scanning structural parameters such as Amp and Dz for obtaining a narrow FWHM of SW. First, three configurations of the wavy structures are considered with different layer numbers (N = 40, 30, and 20). Three related maps of their SWs are calculated by scanning two parameters Amp from 0.005 to 0.5 μm and Dz from 0.2 to 0.4 μm, as shown in Figure 3a. It is noted that their transversal periods are identical, referred to [15].

**Figure 3.** (**a**) The SWs of autocloning-mode-based PhCs structures with different amplitudes of harmonic modulation and longitudinal periods Dz for different layer structures. The smallest SW can be obtained at Dz indicated by the dashed white lines. (**b**) Variation of the SW as a function of Amp at Dz indicated by the dashed white lines in (**a**).

The smallest SWs for the 40-, 30-, and 20-layer configurations are found at the longitudinal periods Dz of 0.22, 0.24, and 0.3 μm, respectively. The corresponding variances of the SWs at the specific Dz, highlighted by white dashed lines in Figure 3a, are shown in Figure 3b. The optimal Amps of these three harmonic structures with N = 40, 30, and 20 are 0.19, 0.275, and 0.435 μm, respectively. As a result, in the three configurations with the different layer numbers, the 40-layer structure presents the smallest SW of 9.32◦ where BW of 1 λ is considered.

Furthermore, the diffraction efficiency for the center of the filtered beam is also considered as depicted in Figure 4. The variations of the SW and DEt0 with respect to Dz are studied for the three different configurations of the wavy structures. The SWs of 9.32◦, 10.18◦, and 12.68◦ and DEt0 values of 42%, 60%, and 58% are achievable for the 40-, 30-, and 20-layer structures, respectively. Although the narrowest SW of 9.32◦ of a filtering beam is achievable by the configuration with 40 layers, it brings in the lowest diffraction efficiency.

**Figure 4.** The smallest SWs and the 0th-order diffraction efficiency (DEt0) varied with the longitudinal period Dz for three autocloning-mode-based PhCs with N = 40 layers (**a**), 30 layers (**b**), and 20 layers (**c**). The minimum of the SW of 9.32◦ can be obtained for the 40-layer structure at Dz = 0.22 μm and Amp = 0.19 μm. The minimum of the SW of 10.18◦ can be obtained for the 30-layer structure at Dz = 0.24 μm and Amp = 0.275 μm. The minimum of the SW of 12.68◦ can be obtained for the 20-layer structure at Dz = 0.3 μm and Amp = 0.435 μm. The DEt0 values of the 40-, 30-, and 20-layered structure are 42%, 60%, and 58%, respectively. The transversal period D*x* is constant, i.e., 1.67 μm.

For all the above simulations, the transversal period Dx of 1.67 μm is considered for the low-pass filtering design. As generally known, the transversal period of PhCs plays an important role in diffraction patterns. As the transversal period Dx is changed, the far-field distribution changes dramatically, leading to the low- and high-pass filtering effects. Hence, Figure 5 shows the variation of far-field distributions by varying Dx for three configurations with 40-, 30-, and 20-layer structures. The best low-angle-pass filtering performance (SW = 8.89◦ and DEt0 = 60%) is found at the 30-layer structure with Amp = 0.275 μm, Dz = 0.24 μm, and Dx = 1.7 μm.

**Figure 5.** 2-D maps of normalized far-field distributions with respect to the angle of incidence and transversal period D*x* for the multilayered structures with 40 layers (**a**), 30 layers (**b**), and 20 layers (**c**). In the 40-layer structure, the minimum FWHM is 9.23◦ at Dx = 1.69 μm with the DEt0 of 39%. In the 30-layer structure, the minimum FWHM is 8.89◦ at Dx = 1.7 μm with the DEt0 of 60%. In the 20-layer structure, the minimum FWHM is 10.07◦ at Dx = 2.1 μm with the DEt0 of 49%.

A more optimal parameter Dx of 1.7 μm is found after scanning the far-field distribution with respect to Dz, as shown in Figure 4. Even though the longitudinal period Dx of 1.67 μm is considered for the previously studied structure, it has shown that the procedure of the search of optimal parameters such as Amp, layers, and Dz for the low-pass angular filtering is beneficial. For example, as shown in Figures 3 and 4, we obtain SWs of 12.68◦ and 9.32◦ and DEt0 values of 58% and 42% for the two configurations with 20 layers and 40 layers, respectively. Although the 40-layer structure presents stronger spatially filtering, it requires more layers to achieve the narrower SW; however, the DEt0 is lower.

Furthermore, the number of layers is also a concern for the fabrication of a wavy-like multilayered structure because the experimental modulation of the wavy structure may be reduced with the increased number of layers [19]. Even if the fabrication of the 20-layer structure may be controllable, a weaker filtering effect is obtained. Obviously, the fabricating feasibility concerning the number of layers should be also considered at the stage in the design process. Therefore, a series of the analysis of structural parameters is helpful to define the optimization range and consider fabrication feasibility.

In our work, the optimization of this wavy structure by using the simplex method has been demonstrated, as shown in Figure 6. All parameters of the multilayered structure are optimized simultaneously and the target of the optimization is to realize a low-pass filter design that achieves the narrowest SW for a Gaussian beam. To ensure fabrication feasibility, the optimizing range of the number of the layer should not be more than 40 layers. The optimizing ranges of Amp and the longitudinal period Dz should be less than 0.5 μm and 0.4 μm, respectively, as referred to the previous study. As a result, the final optimal parameters of the structure are obtained as follows: N = 27, Amp = 0.36 μm, Dz = 0.26 μm, and Dx = 2 μm. The far-field distributions in Figure 6b,c show the narrowest SW1' = 5.92◦ and SW2' = 6.58◦ for two Gaussian beams with BW1 = λ and BW2 = 2λ, respectively. The DEt0 of the optimal result is 64%. The narrower SW and the higher diffraction efficiency are achieved by the optimized structure with the less number of layer, compared with that of the previous study. Therefore, the optimizing results not only present the best spatial filtering to narrow the SW, but also demonstrate the use of the structures with fewer layers is a more feasible approach to a manufacturing process.

**Figure 6.** The optimization results of the wavy structure by using the simplex method. (**a**) The illustration of the optimal structures with parameters: N = 27, Amp = 0.36 μm, Dz = 0.26 μm, and Dx = 2 μm. Far-field distributions for two different incident Gaussian beam widths BW1 = λ (**b**) and BW2 = 2λ (**c**), where λ is the operating wavelength. The narrowest FWHM of a Gaussian beam with BW = λ is 5.92◦.

#### **4. Conclusions**

This paper has demonstrated the PhCs filtering effect for a light beam. The autocloned PhCs present transversal and longitudinal periods, the key element to modulate spatial spectra. The paper has studied different configurations of the multilayered structures based on the harmonic modulation. The narrowest SW of 8.89◦ and the diffraction efficiency of 60% are obtained for a 30-layer structure after a series of the scanning procedure in the structural parameters.

Considering feasible fabrication, the optimization has been further studied for practical use by using the simplex method. The number of layers N for our optimal structure is reduced to 27 and a stringer spatial filtering performance provides an SW of 5.92◦ and the diffraction efficiency of 0th transmission is 64%. For an autocloned PhC with several layers, another approach, such as an autocloned blazed modulation, could be more applicable for fabrication, although the fabrication of a tilted blazed profile is not as simple as that of the sinusoidal modulation.

**Author Contributions:** Conceptualization, Y.-C.C.; methodology, software and validation, Y.-C.C., Y.-C.L. and P.-Y.W.; investigation, writing—original draft preparation, writing—review and writing—editing, Y.-C.C. and P.-Y.W.; visualization, Y.-C.C., P.-Y.W. and Y.-C.L.; supervision, Y.-C.C.; funding acquisition, Y.-C.C.

**Funding:** This work was financially supported by the Young Scholar Fellowship Program by the Ministry of Science and Technology (MOST) in Taiwan, under grant number MOST108-2636-M-027-001.

**Acknowledgments:** The authors would like to thank the Headquarters of University Advancement at National Cheng Kung University (NCKU) for the funding support of the Higher Education Sprout Project of Ministry of Education (MOE).

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Appendix A. Method**

The diffraction efficiency of zero-order transmission is calculated by using the RCWA method, which solves full vectorial versions of Maxwell's equations. It covers a wide range of scattering problems on the wavy structures with horizontally periodic boundary conditions and vertically perfectly matched layer (PML) boundary. The incident and transmitted fields are defined at the simulation domain along the z-direction. A plane wave incidence is used in the RCWA simulation.

The far-field distribution is obtained by using two methods. First, the FDTD method provides steady-state near-field distributions at an output plane while an incident Gaussian beam with a BW considered. Second, we use the discrete Fourier transform (DFT) to convert the steady-state near-field distribution to the far-field distribution. The SW is defined as the FWHM of the far-field distribution. The Gaussian beam is launched at the position of 1 μm above the wavy structure. The horizontal

simulation domain is equal to 50 periods of the proposed wavy structure and the vertical boundary condition is the PML.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
