Photonic Band Structure and Transmittance of the Superconductor Photonic Crystal

1. Introduction

The photonic crystal (PhC) is formed with a dielectric periodic structure and exhibits new electromagnetic phenomena (John, 1987). It shows some properties analog to the semiconductor, such as the photonic band structure (PBS) including photonic passing bands and photonic band gaps (PBGs), and complicated dispersion relations. In analogous to the electron transport in the semiconductor, the Bloch theorem is also applied to describe electromagnetic waves propagating in the PhC very well.

The PBS strongly depends on refracted indices of constituent materials and the geometry of the PhC. Once the materials and geometry structure of a PhC are constructed, the possible way to change its PBS is tuning the refracted indices of its constituent materials utilizing the temperature effect, the external electric field effect, or the external magnetic field effect, etc (Busch & John, 1999; Kee & Lim, 2001; Kee et al., 2000, 2001; Figotin et al., 1998; Takeda & Yoshino, 2003a, 2003b, 2003c, 2003d, 2004). For PhCs composed of ferroelectric or ferromagnetic materials, PBSs can be tuned by the external electric field effect and the external magnetic field effect (Busch & John, 1999; Figotin et al.). On the other hand, the variation on the PBS of the liquid-crystal PhC controlled by the external electric field or the temperature has also been investigated (Kee & Lim, 2001b; Takeda & Yoshino, 2003a, 2003b, 2003c, 2003d, 2004). Another potential material that can be used to tune the PBS is the superconductor by varying the temperature and the external magnetic field (Lee et al., 1995; Raymond Ooi et al., 2000; Takeda & Yoshino, 2003e).

In our previous works, we have designed a tunable PhC Mach-Zehnder interferometer composed of copper oxide high-temperature superconductors (HTSCs) utilizing the temperature modulation to reach the on and off states (Pei & Huang, 2007a). The Mach-Zehnder interferometer, whose path-length difference of two arms is fixed after designed, can be realized as an optical switching device or sensor due to the temperature effect. In the output, the signals from two arms interfere with each other, and the phases of these two signals can be modulated by HTSCs. Besides, we also discussed the superprism effect in the superconductor PhC (Pei & Huang, 2007b). The superprism effect was demonstrated experimentally by Kosaka et al. in 1998 (Kosaka et al., 1998). They found that the refracted angle of a light beam in a PhC is very sensitive to the incident angle and wavelength. The basic explanation of the superprism effect is based on the anomalous dispersion characteristics of the PBS. The propagation direction of light in the PhC is the same as the direction of the group velocity, which is determined by the equifrequency surfaces (EFS). The group velocity is normal to the EFS at a certain wave vector and is defined as$${\stackrel{\rightharpoonup }{v}}_g=gra{d}_{\stackrel{\rightharpoonup }{k}}\omega$$, where $$\stackrel{\rightharpoonup }{k}$$ and ω are the wave vector and the frequency, respectively. Notomi has published a detailed study on the superprism effect (Notomi, 2000). In our work, we not only study the transmission of light propagating through the superconductor PhC, but also pay lots of attentions on the refraction. The result shows that the refraction can be changed sensitively by the temperature of the superconductor.

In this chapter, we deduce the way to calculate the PBS of the superconductor PhC based on the plane wave expansion method first. It is not like the way to calculate the PBS of the PhC only composed of dielectric materials. Second, the finite-difference time-domain (FDTD) method for the PhC composed of dispersive materials such as superconductors are derived carefully. The time-domain auxiliary differential equations (ADEs) are introduced to represent effects of currents in dispersive materials. The ADE-FDTD algorithm can be used to calculate the transmission of the finite superconductor PhC. It has also been used in our previous works to discuss the tunability of the PhC Mach-Zehnder interferometer composed of HTSCs and the superprism effect in the superconductor PhC.

Finally, the internal-field expansion method developed by Sakoda is also introduced (Sakoda, 1995a, 1995b, 2004). This method is used to calculate the transmission of the two-dimensional PhC composed of air cylinders embedded in certain background medium. It is much like the grating theory that describes the scattering waves as Bragg waves. He successively calculated the transmission and the Bragg reflection spectra using this method, and also mentioned that the existences of the uncoupled modes (Sakoda, 1995a, 1995b). However, this method has not been yet verified on the superconductor PhC. We use this method to calculate the transmission of the finite superconductor PhC and compare the result of it with that of the ADE-FDTD method.

2. The plane wave expansion method for calculating the photonic band structure of the superconductor photonic crystal

The superconductivity of the superconductor is strongly sensitive to the temperature and the external magnetic field. We only discuss the temperature effect in this chapter. The PhC structure is composed of superconductor cylinders with triangular lattice in air as shown in Fig. 1. The two-fluid model is used to describe the electromagnetic response of a typical superconductor without an additional magnetic field (Tinkham, 2004), and it describes that the electrons occupy two states. One is the superconducting state, in which the superconducting electrons of density N_s(x, y) are paired and transport with no resistance. The definition of the superconducting state under the temperature and magnetic field effects is shown in Fig. 2. The other is the normal state, in which the normal conducting electrons of density N_n(x, y) act like electrons in general materials with a nonzero resistance. Both superconducting and normal conducting electrons coexist in the superconductor when the temperature is lower than the critical temperature. This model also characterizes the performance of high-frequency superconductive devices very well (Van Duzer & Truner, 1998).

Figure 1.

The cross-section of the two-dimensional PC formed in a triangular array.

Figure 2.

The definition of the superconducting state under the effects of the temperature and the magnetic field.

Utilizing this model, the E-polarized light with its electric field parallel to the z-axis (TM mode) is incident on a two-dimensional PhC lying in the x-y plane. In the presence of the external electric field, superconducting and normal conducting current densities J_sz and J_nz flowing along the z-axis can be expressed as the following equations (Tinkham, 2004):

$$\frac{\partial J_{sz}(x,y)}{\partial t}={\epsilon }_0{\epsilon }_1(x,y){\omega }_p^s{(x,y)}^2{E}_z(x,y)$$E1

$$\tau \frac{\partial J_{nz}(x,y)}{\partial t}+J_{nz}(x,y)=\tau {\epsilon }_0{\epsilon }_1(x,y){\omega }_p^n{(x,y)}^2{E}_z{(x,y)}_{,}$$E2

where$${\omega }_p^s(x,y)$$ and $${\omega }_p^n(x,y)$$ are the plasma frequencies of the superconducting and normal conducting electrons given by

$${\omega }_p^s(x,y)=\sqrt{N_s(x,y)e^2/{\epsilon }_0{\epsilon }_1(x,y)m}={c/\lambda (x,y)\sqrt{{\epsilon }_1(x,y)}}_{,}$$E3

$${\omega }_p^n(x,y)=\sqrt{N_n(x,y)e^2/{\epsilon }_0{\epsilon }_1(x,y)m}$$E4

λ(x, y) is the London penetration depth, ε₁(x, y) is the distribution of the dielectric constant, τ is the relaxation time, c is the wave velocity in free space, and m is the mass of the electron. Because the incident electric field is harmonic with frequency ω, the induced J_sz and J_nz also have the same oscillating period. Eqs. (1) and (2) can be further expressed as follows:

$$J_{sz}(x,y)=-i{\epsilon }_0{\epsilon }_1(x,y)\frac{{\omega }_p^s{(x,y)}^2}{\omega }{E}_z(x,y)$$E5

$$J_{nz}(x,y)={\epsilon }_0{\epsilon }_1(x,y)\frac{{\omega }_p^n{(x,y)}^2\tau }{(1+i\omega \tau )}{E}_z(x,y)$$E6

.

Substituting Eqs. (5) and (6) into the E-polarized wave equation results in the following equation (Pei & Huang, 2007a):

$$\begin{array}{l}\frac{{\partial }^2{E}_z(x,y)}{\partial x^2}+\frac{{\partial }^2{E}_z(x,y)}{\partial y^2}=-\frac{{\omega }^2{\epsilon }_1(x,y)}{c^2}{E}_z(x,y)+i{\mu }_0\omega \left\{J_{sz}(x,y)+J_{nz}(x,y)\right\}\\ =-\frac{{\omega }^2{\epsilon }_1(x,y)}{c^2}\left\{1-\frac{{\omega }_p^s{(x,y)}^2}{{\omega }^2}-\frac{{\omega }_p^n{(x,y)}^2\tau }{\omega (1+i\omega \tau )}\right\}{E}_z(x,y)\\ =-\frac{{\omega }^2}{c^2}{\epsilon }_s(x,y,\omega )E_z(x,y)\end{array}$$E7

where ε_s(x, y, ω) is the effective dielectric function given by

$${\epsilon }_s(x,y,\omega )={\epsilon }_1(x,y){\left\{1-\frac{{\omega }_p^s{(x,y)}^2}{{\omega }^2}-\frac{{\omega }_p^n{(x,y)}^2\tau }{\omega (1+i\omega \tau )}\right\}}_{.}$$E8

For HTSCs, optical characteristics show the anisotropic properties (Takeda & Yoshino, 2003e). The electric fields parallel and perpendicular to the c-axis feel different dielectric indices. However, Eq. (7) is still valid even for anisotropic materials (Lee et al., 1995). When the electric fields are parallel to the c-axis, plasma frequencies are in the microwave and far-infrared regions (Takeda & Yoshino, 2003e). In our study, the z-axis is chosen as the c-axis.

In the superconducting state, the electromagnetic wave can propagate in the range of the London penetration depth. The London penetration depth is dependent on the temperature T, which can be expressed as $$\lambda \left(T\right)={\lambda }_0/\sqrt{1-{\left(T/T_c\right)}^4}$$ (Zhou, 1999), where T_c and λ₀ are the critical temperature and the London penetration depths at the absolute zero temperature, respectively. When the temperature is above about 0.8 times the critical temperature, the London penetration depth increases rapidly and then approaches infinity as the temperature is close to T_c. Besides, $${\omega }_p^s(x,y)$$strongly depends on the London penetration depth as well as the temperature. Based on the experimental results (Shibata & T. Yamada, 1996; Matsuda, 1995) in the far-infrared region, the small contribution of the normal conducting electrons can be neglected and the plasma frequency $${\omega }_p^s(x,y)$$ can be assumed to be uniform within the rods. Then the third term on the right side of Eq. (6) can be dropped and then simplified as

$${\epsilon }_s(x,y,\omega )\approx {\epsilon }_1(x,y){\left(1-\frac{{\omega }_p^s{(x,y)}^2}{{\omega }^2}\right)}_{.}$$E9

Eq. (9) is known as the Drude model (Grosso & Parravicini, 2000) which can also be applied to the kind of PhCs constituting metallic components.

Kuzmiak et al. (Kuzmiak et al., 1994) has dealt with the two-dimensional PhC containing metallic components. We use the same method based on the plane-wave expansion to calculate the PBSs of the superconductor PhCs. In this method, the dielectric function of the PhC is directly expanded in a Fourier series. The dielectric constant of the PhC can be written explicitly in the form

$${\epsilon }_s(\stackrel{\rightharpoonup }{r},\omega )=1+\left[{\epsilon }_s(\omega )-1\right]{{\displaystyle \sum _{\stackrel{\rightharpoonup }{a}}S(\stackrel{\rightharpoonup }{r}-\stackrel{\rightharpoonup }{a})}}_{,}$$E10

where the function $$S(\stackrel{\rightharpoonup }{r})=1$$and$${\epsilon }_s(\stackrel{\rightharpoonup }{r},\omega )={\epsilon }_s(\omega )$$if $$\stackrel{\rightharpoonup }{r}$$ is inside the cylinder, and $$S(\stackrel{\rightharpoonup }{r})=1$$ and$${\epsilon }_s(\stackrel{\rightharpoonup }{r},\omega )=1$$ if $$\stackrel{\rightharpoonup }{r}$$ is outside the cross section. The expansion of $${\epsilon }_s(\stackrel{\rightharpoonup }{r},\omega )$$ in a Fourier series on reciprocal lattice vectors $$\stackrel{\rightharpoonup }{G}$$ is

$${\epsilon }_s(\stackrel{\rightharpoonup }{r},\omega )={{\displaystyle \sum _{\stackrel{\rightharpoonup }{G}}{\epsilon }_s(\stackrel{\rightharpoonup }{G},\omega )e^{i\stackrel{\rightharpoonup }{G}\cdot \stackrel{\rightharpoonup }{r}}}}_{,}$$E11

where the Fourier coefficient$${\epsilon }_s(\stackrel{\rightharpoonup }{G},\omega )$$is given by

$$\begin{array}{l}{\epsilon }_s(\stackrel{\rightharpoonup }{G},\omega )=\frac{1}{a_c}{\displaystyle {\int }_{V_0^{(2)}}{\epsilon }_s(\stackrel{\rightharpoonup }{r},\omega )e^{-i\stackrel{\rightharpoonup }{G}\cdot \stackrel{\rightharpoonup }{r}}d\stackrel{\rightharpoonup }{r}}\\ ={\delta }_{\stackrel{\rightharpoonup }{G}0}+\left[{\epsilon }_s(\omega )-1\right]\frac{1}{a_c}{\displaystyle {\int }_{V_0^{(2)}}S(\stackrel{\rightharpoonup }{r})e^{-i\stackrel{\rightharpoonup }{G}\cdot \stackrel{\rightharpoonup }{r}}d\stackrel{\rightharpoonup }{r}}\end{array}$$E12

where the integral is now over the 2D unit cell $$V_0^{(2)}$$and a_c is the area of the unit cell in the PhC. Eq. (12) can be expressed as

$${\epsilon }_s(\stackrel{\rightharpoonup }{G}=\stackrel{\rightharpoonup }{0},\omega )=1+\left[{\epsilon }_1\left(1-\frac{{\omega }_p^s{}^2}{{\omega }^2}\right)-1\right]f_{,}$$E13

$${\epsilon }_s(\stackrel{\rightharpoonup }{G}\ne \stackrel{\rightharpoonup }{0},\omega )={\epsilon }_1\left(\frac{{\epsilon }_1-1}{{\epsilon }_1}-\frac{{\omega }_p^s{}^2}{{\omega }^2}\right)\frac{2J_1(GR)}{GR}{f}_{,}$$E14

where f is the filling fraction. For the triangular superconductor PhC, the filling fraction of the superconductor rod in a unit cell is$$f=(2\pi /\sqrt{3})R^2/a^2$$. According to the Bloch theory, the electric field can be expanded in the form

$$E_z(x,y)={{\displaystyle \sum _{\stackrel{\rightharpoonup }{G}}{E}_z(\stackrel{\rightharpoonup }{k}|\stackrel{\rightharpoonup }{G})e^{i(\stackrel{\rightharpoonup }{k}+\stackrel{\rightharpoonup }{G})\cdot \stackrel{\rightharpoonup }{r}}}}_{,}$$E15

where $$\stackrel{\rightharpoonup }{k}=k_x\widehat{i}+k_y\widehat{j}$$ is the wave vector of the electromagnetic waves propagating inside the PhC. Substituting Eqs. (11) and (15) into Eq. (7), we obtain a set of equations for the coefficients$$E_z(\stackrel{\rightharpoonup }{k}|\stackrel{\rightharpoonup }{G})$$. It is a standard eigenvalue problem of a real and symmetric matrix with respect to the frequency ω. The set of equations for coefficients $$E_z(\stackrel{\rightharpoonup }{k}|\stackrel{\rightharpoonup }{G})$$ shows as follows:

$$\begin{array}{l}{\displaystyle \sum _{{\stackrel{\rightharpoonup }{G}}^{\prime }}\left[{\left(\stackrel{\rightharpoonup }{k}+\stackrel{\rightharpoonup }{G}\right)}^2{\delta }_{\stackrel{\rightharpoonup }{G},{\stackrel{\rightharpoonup }{G}}^{\prime }}+{\epsilon }_1{\left(\frac{{\omega }_p^s}{c}\right)}^2\frac{1}{a_c}{\displaystyle {\int }_{}S(\stackrel{\rightharpoonup }{r})e^{-i(\stackrel{\rightharpoonup }{G}-{\stackrel{\rightharpoonup }{G}}^{\prime })\cdot \stackrel{\rightharpoonup }{r}}d\stackrel{\rightharpoonup }{r}}\right]E_z(\stackrel{\rightharpoonup }{k}|{\stackrel{\rightharpoonup }{G}}^{\prime })}\\ ={\left(\frac{\omega }{c}\right)}^2{\displaystyle \sum _{{\stackrel{\rightharpoonup }{G}}^{\prime }}\left[\left(1-f+{\epsilon }_{1}f\right){\delta }_{\stackrel{\rightharpoonup }{G},{\stackrel{\rightharpoonup }{G}}^{\prime }}+\frac{{\epsilon }_1-1}{a_c}{\displaystyle {\int }_{}S(\stackrel{\rightharpoonup }{r})e^{-i(\stackrel{\rightharpoonup }{G}-{\stackrel{\rightharpoonup }{G}}^{\prime })\cdot \stackrel{\rightharpoonup }{r}}d\stackrel{\rightharpoonup }{r}}\right]}{E}_z(\stackrel{\rightharpoonup }{k}|{\stackrel{\rightharpoonup }{G}}^{\prime })\end{array}$$E16

Rearranging Eq. (16) that we have

$$\begin{array}{l}{\displaystyle \sum _{{\stackrel{\rightharpoonup }{G}}^{\prime }}\left\{{\left(\stackrel{\rightharpoonup }{k}+\stackrel{\rightharpoonup }{G}\right)}^2{\delta }_{\stackrel{\rightharpoonup }{G},{\stackrel{\rightharpoonup }{G}}^{\prime }}+f{\left(\frac{{\omega }_s^p}{c}\right)}^2{\delta }_{\stackrel{\rightharpoonup }{G},{\stackrel{\rightharpoonup }{G}}^{\prime }}+f{\left(\frac{{\omega }_s^p}{c}\right)}^2\frac{2J_1\left(\left|\stackrel{\rightharpoonup }{G}-{\stackrel{\rightharpoonup }{G}}^{\prime }\right|R\right)}{\left(\left|\stackrel{\rightharpoonup }{G}-{\stackrel{\rightharpoonup }{G}}^{\prime }\right|R\right)}\right\}{E}_z(\stackrel{\rightharpoonup }{k}|{\stackrel{\rightharpoonup }{G}}^{\prime })}\\ =\frac{{\omega }^2}{c^2}{\displaystyle \sum _{{\stackrel{\rightharpoonup }{G}}^{\prime }}\left\{\left[{\epsilon }_0\left(1-f\right)+{\epsilon }_{1}f\right]{\delta }_{\stackrel{\rightharpoonup }{G},{\stackrel{\rightharpoonup }{G}}^{\prime }}-\left({\epsilon }_0-{\epsilon }_1\right)\frac{2J_1\left(\left|\stackrel{\rightharpoonup }{G}-{\stackrel{\rightharpoonup }{G}}^{\prime }\right|R\right)}{\left(\left|\stackrel{\rightharpoonup }{G}-{\stackrel{\rightharpoonup }{G}}^{\prime }\right|R\right)}\right\}{E}_z(\stackrel{\rightharpoonup }{k}|\stackrel{\rightharpoonup }{G})}\end{array}$$E17

To solve this matrix eigenvalue problem, the frequencies can be determined at a certain wave vector and the whole PBS can be obtained.

3. The E-polarized photonic band structure

In our designed device, we used high-T_c superconductor Bi_1.85Pb_0.35Sr₂Ca₂Cu_3.1O_y (Takeda & Yoshino, 2003). Previous study (Takeda & Yoshino, 2003) utilized parallel copper oxide HTSCs rods to form PhCs with square lattices repeating in two-dimensional directions (x-y plane). The authors theoretically investigated the tunability of the photonic band gap (PBG) of the two-dimensional PhC by changing temperatures of superconductors and external magnetic fields. The PhC structure we discuss here is composed of superconductor cylinders with triangular lattice in air as shown in Fig. 1. The E-polarized electromagnetic wave with the electric field parallel to the extended direction of the rod propagates in the x-y plane. Adjusting the temperature of the superconductor can control the refracted index of the superconductor as well as the PBS of the superconductor PhC. When T $$\le$$ T_c is satisfied, the dependence of the plasma frequency on the temperature is given by (Zhou, 1999)

$${\omega }_p^s(T)={\omega {}_p{}^s(0)/\sqrt{1-{\left(T/T_c\right)}^4}}_{.}$$E18

For this superconductor, the London penetration depth of the copper oxide HTSCs is λ = 23 μm at T = 5 K, the critical temperature T_c = 107 K, and the dielectric constant is ε₁ = 12 (Shibata & Yamada, 1996). When T = 5 K, we obtain$${\omega }_p^s/c\approx 1.3\times {10}^4$$cm^-1.

Figure 3.

The PBS of the PhC composed of superconductor cylinders at T = 5 K with radius of cylinders r = 0.2a and lattice constant a = 100 μm. The first Brillouin zone shows at the right lower corner. In this figure, the region between two paired horizontal lines is the PBG.

The theory discussed in Section 2 is used to calculate the PBS along the three directions ΜΓ, ΓΚ, and ΚΜ in the reduced Brillouin zone when the periodic lattice constant of the PhC is a = 100 μm, the radius of cylinders is r = 0.2a, and the overall temperature is fixed at 5 K. The PBS is shown in Fig. 3 and the first Brillouin zone at the right lower corner. The reduced Brillouin zone is denoted as the triangle ΓΚΜ. From Eq. (9), we can see that the optical response of the superconductor under the E-polarized wave is the same as that of the metal described by the Drude model. The lowest point of the first band for a metal or metal-like material is above zero frequency, which is not like a non-dispersive material whose lowest point of the first band is at zero frequency. A PBG exists from zero to a certain frequency ω_lowest, which means that the light can propagate in the PhC only in the frequency range above ω_lowest. In Fig. 3, the region between two paired horizontal lines is the PBG region. The PhC has a large second PBG, which is located in the frequency range from 0.33 to 0.47 (2πc/a). The third PBG is located in the frequency range from 0.595 to 0.605 (2πc/a).

4. The finite-difference time-domain method for the photonic crystal composed of dispersive materials

In 1966, K. S. Yee first provided the FDTD method to solve electromagnetic scattering problems (Yee, 1966). The Yee’s equations are obtained to discretize Maxwell’s equations in time and space. The fields on the nodal points of the space-time mesh can be calculated in an iteration process when the source is excited. Because the finite resource of the hardware limits the size of simulation domain, an absorbing boundary condition (ABC) needs to be set on the outer surface of the computational domain. In 1994, Berenger proposed a perfectly matched layer (PML), which is an artificial electromagnetic wave absorber with electric conductivity σ and magnetic conductivity σ^* (Berenger, 1994). The PML absorbs outgoing waves very well, so it can simulate the electromagnetic wave propagating in free space. Therefore, we apply the PML as the absorbing layer used in the FDTD method.

In the FDTD method, Maxwell's equations are solved directly in time domain via finite differences and time steps without any approximations or theoretical restrictions. The basic approach is relatively easy to understand and is an alternative to more usual frequency-domain approaches, so this method is widely used as a propagation solution technique in integrated optics. Imagine a region of space where no current flows and no isolated charge exists. Maxwell's curl equations can be written in Cartesian coordinates as six simple scalar equations. Two examples are:

$$\frac{\partial E_z}{\partial t}=-\frac{1}{\epsilon }{\left(\frac{\partial H_y}{\partial x}-\frac{\partial H_x}{\partial y}\right)}_{,}$$E19

$$\frac{\partial H_z}{\partial t}=-\frac{1}{{\mu }_0}{\left(\frac{\partial E_x}{\partial y}-\frac{\partial E_y}{\partial x}\right)}_{.}$$E20

The most common method to solve these equations is based on Yee's mesh and computes the E and H field components at points on a grid with grid points spaced Δx, Δy, and Δz apart, which are named grid sizes. The E and H field components are then interlaced in all three spatial dimensions. Furthermore, time is broken up into discrete steps of Δt. The E field components are then computed at times t = nΔt and the H at times t = (n + 1/2)Δt, where n is an integer representing the computing step. For example, the E field at a time t = nΔt is equal to the E field at t = (n - 1)Δt plus an additional term computed from the spatial variation, or curl of the H field at time t. This method results in six equations that can be used to compute the field at a given mesh point, denoted by integers i, j, k

$$E_z{|}_{i,j}^{n+1}=E_z{|}_{i,j}^n+\frac{\Delta t}{{\left({\epsilon }_0{\epsilon }_1\right)}_{i,j}}(\frac{H_y{|}_{i+1/2,j}^{n+1/2}-H_y{|}_{i-1/2,j}^{n+1/2}}{\Delta x}{-\frac{H_x{|}_{i,j+1/2}^{n+1/2}-H_x{|}_{i,j-1/2}^{n+1/2}}{\Delta y})}_{,}$$E21

$$H_z{|}_{i,j}^{n+1}=H_z{|}_{i,j}^n+\frac{\Delta t}{{\mu }_0}(\frac{E_x{|}_{i,j}^n-E_x{|}_{i,j-1}^n}{\Delta y}{-\frac{E_y{|}_{i,j}^n-E_y{|}_{i-1,j}^n}{\Delta x})}_{.}$$E22

These equations are iteratively solved in a leapfrog manner, alternating between computing the E and H fields at subsequent Δt/2 intervals. The grid sizes and time step in 2D simulations are set Δx =Δy and Δt = Δx/2.

The method for implementing FDTD models of dispersive materials utilizes ADE equations which describe the time variation of the electric current densities (Taflove & Hagness, 2005). These equations are time-stepped synchronously with Maxwell’s equations. ADE-FDTD method is a second-order accurate method.

Consider a dispersive medium whose Ampere’s Law can be expressed as

$$\nabla \times \stackrel{\rightharpoonup }{H}(t)={\epsilon }_0{\epsilon }_{\infty }\frac{\partial \stackrel{\rightharpoonup }{E}(t)}{\partial t}+\sigma \stackrel{\rightharpoonup }{E}(t)+{{\displaystyle \sum _p{\stackrel{\rightharpoonup }{J}}_p(t)}}_{,}$$E23

where $${\stackrel{\rightharpoonup }{J}}_p(t)$$ is the polarization current. The goal of the ADE technique is to develop a simple time-stepping scheme for$${\stackrel{\rightharpoonup }{J}}_p(t)$$. In our superconductor system, $${\stackrel{\rightharpoonup }{J}}_{sz}$$and $${\stackrel{\rightharpoonup }{J}}_{nz}$$contribute to $$\sigma \stackrel{\rightharpoonup }{E}$$ and$${\stackrel{\rightharpoonup }{J}}_p$$, respectively, so Eq. (23) can be rewritten as

$$\nabla \times \stackrel{\rightharpoonup }{H}(x,y,t)={\epsilon }_0{\epsilon }_1(x,y)\frac{\partial {\stackrel{\rightharpoonup }{E}}_z(x,y,t)}{\partial t}+{\stackrel{\rightharpoonup }{J}}_{sz}(x,y,t)+{\stackrel{\rightharpoonup }{J}}_{nz}{(x,y,t)}_{.}$$E24

Another time-dependent Maxwell’s curl equation is

$$\nabla \times \stackrel{\rightharpoonup }{E}(x,y,t)=-\mu (x,y){\frac{\partial \stackrel{\rightharpoonup }{H}(x,y,t)}{\partial t}}_{,}$$E25

where μ(x,y) is the position dependent permeability of the material. Eqs. (24) and (25) can be discretized in two-dimensional space and time by the Yee-cell technique(Yee, 1966). Eqs. (1) and (2) are the required ADEs for $${\stackrel{\rightharpoonup }{J}}_{sz}$$ and$${\stackrel{\rightharpoonup }{J}}_{nz}$$. They both can be easily and accurately implemented in an FDTD code using the semi-implicit scheme where fields at time-step n+1 are created and updated by fields known at time-step n. Then, we implement Eqs. (1) and (2) in an FDTD code by finite differences, centered at time-step n+1/2:

$$\frac{J_{sz}{|}_{i,j}^{n+1}-J_{sz}{|}_{i,j}^n}{\Delta t}={\left({\epsilon }_0{\epsilon }_1{\omega }_p^s\right)}_{i,j}^2{\frac{E_z{|}_{i,j}^{n+1}+E_z{|}_{i,j}^n}{2}}_{,}$$E26

$$\frac{J_{nz}{|}_{i,j}^{n+1}-J_{nz}{|}_{i,j}^n}{\Delta t}+\frac{J_{nz}{|}_{i,j}^{n+1}+J_{nz}{|}_{i,j}^n}{2\tau }={\left({\epsilon }_0{\epsilon }_1{\omega }_p^n{}^2\right)}_{i,j}{\frac{E_z{|}_{i,j}^{n+1}+E_z{|}_{i,j}^n}{2}}_{.}$$E27

Solving Eqs. (26) and (27) for $$J_{sz}{|}_{i,j}^{n+1}$$ and$$J_{nz}{|}_{i,j}^{n+1}$$, we obtain

$$J_{sz}{|}_{i,j}^{n+1}=J_{sz}{|}_{i,j}^n+\frac{\Delta t}{2}{\left({\epsilon }_0{\epsilon }_1{\omega }_p^s\right)}_{i,j}^2\left(E_z{|}_{i,j}^{n+1}+E_z{|}_{i,j}^n\right)$$E28

$$J_{nz}{|}_{i,j}^{n+1}=\frac{\left(1-\Delta t/2\tau \right)}{\left(1+\Delta t/2\tau \right)}{J}_{nz}{|}_{i,j}^n+\frac{\Delta t}{2\left(1+\Delta t/2\tau \right)}{\left({\epsilon }_0{\epsilon }_1{\omega }_p^n{}^2\right)}_{i,j}{\left(E_z{|}_{i,j}^{n+1}+E_z{|}_{i,j}^n\right)}_{.}$$E29

Then we can evaluate Eq. (24) at time-step n + 1/2:

$$\nabla \times {\stackrel{\rightharpoonup }{H}|}_{i,j}^{n+1/2}={\left({\epsilon }_0{\epsilon }_1\right)}_{i,j}\left(\frac{{{\stackrel{\rightharpoonup }{E}}_z|}_{i,j}^{n+1}-{{\stackrel{\rightharpoonup }{E}}_z|}_{i,j}^n}{\Delta t}\right)+\frac{1}{2}\left({{\stackrel{\rightharpoonup }{J}}_{sz}|}_{i,j}^{n+1}+{{\stackrel{\rightharpoonup }{J}}_{sz}|}_{i,j}^n\right)+\frac{1}{2}{\left({{\stackrel{\rightharpoonup }{J}}_{nz}|}_{i,j}^{n+1}+{{\stackrel{\rightharpoonup }{J}}_{nz}|}_{i,j}^n\right)}_{.}$$E30

Applying Eq. (30) into the implementation of Eq. (25) in an FDTD code by finite differences, we obtain the E fields at time-step n+1:

$$E_z{|}_{i,j}^{n+1}=E_z{|}_{i,j}^n+\frac{\Delta t}{{\left({\epsilon }_0{\epsilon }_1\right)}_{i,j}}(\frac{H_y{|}_{i+1/2,j}^{n+1/2}-H_y{|}_{i-1/2,j}^{n+1/2}}{\Delta x}{-\frac{H_x{|}_{i,j+1/2}^{n+1/2}-H_x{|}_{i,j-1/2}^{n+1/2}}{\Delta y})}_{.}$$E31

Thus, the ADE-FDTD algorithm for calculating dispersive media has three processes. Starting with the assumed known values of$${\stackrel{\rightharpoonup }{E}}_z$$, $${\stackrel{\rightharpoonup }{J}}_{sz}^n$$, $${\stackrel{\rightharpoonup }{J}}_{nz}^n$$, and$${\stackrel{\rightharpoonup }{H}}^{n+1/2}$$, we first calculate the new$$E_z^{n+1}$$ components using Eq. (31). Next, we calculate the new $$J_{sz}^{n+1}$$ and $$J_{nz}^{n+1}$$ components using Eqs. (28) and (29). Finally, $$H_x^{n+3/2}$$is obtained from $$H_x^{n+1/2}$$ and $$E_z^{n+1}$$ by using Eq. (30). $$H_y^{n+3/2}$$is updated as $$H_x^{n+3/2}$$ be done.

In the end of this section, let us return to discuss the numerical stability. We choose the two-dimensional cell space steps, i.e. Δx and Δy, and the time step Δt based on the required accuracy. The space step is usually chosen less one twentieth of the smallest wavelength in order to avoid the non-physical oscillation. The time step must satisfy the well-known “Courant Condition”:

$$\Delta t\le \frac{1}{V_{\max }}{\left[\frac{1}{{(\Delta x)}^2}+\frac{1}{{(\Delta y)}^2}\right]}^{-1}{}_{,}$$E32

where V_max is the maximum wave velocity in the computational domain.

5. The transmission of the finite photonic crystal composed of the superconductor

In this section, the ADE-FDTD method is used to calculate the transmission of the finite thickness PhC from the frequency 0.01 to 1.00 (2πc/a). As we know, the PBS represents the existing mode with photon energy inside the infinite PhC; but in practice, the thickness of a PhC is always finite. So it is necessarily to calculate the transmission and further compare to the PBS in the previous section. This also verifies calculations of the PBS through the ADE-FDTD method. The triangular PhC is shown in Fig. 1 in which the interface is along the ΓΜ direction (x-direction). Light is normally incident and propagates along the ΓΚ direction (y-direction). The numbers of layers along the x- and y-directions are 40 and 30, respectively. The lattice constant along the x-direction is a₁ and that along y-direction a₂. We choose a₂ = a = 100 μm and a₁ =$$\sqrt{3}$$a₂. For simplicity, the square unit cell ΔxΔy are used in the ADE-FDTD calculations where Δx =Δy = a/30. The time increment is Δt = Δx/2c. The Gaussian wave is supposed to be incident from air on the PhC. The transmissions from 0.01 to 1.00 (2πc/a) are shown in Fig. 4. The increment of the frequency is 0.005 (2πc/a). In this simulation, we consider both currents $${\stackrel{\rightharpoonup }{J}}_{sz}$$ and$${\stackrel{\rightharpoonup }{J}}_{nz}$$. The contribution of the current $${\stackrel{\rightharpoonup }{J}}_{sz}$$ in the calculations is dominant and that of the current $${\stackrel{\rightharpoonup }{J}}_{nz}$$ is very small. It can be seen that almost zero transmission below frequency 0.16 (2πc/a) matches the prediction of the PBS. All transmissions are more than 0.80 at frequencies from 0.16 to 0.33 (2πc/a). This frequency region just corresponds to the first photonic band, in which the highest transmission is close to 1.0 at 0.28 (2πc/a). The second PBG occurs at frequencies between 0.33 and 0.47 (2πc/a). It can be seen that the transmission dramatically drops to nearly zero at 0.33 (2πc/a) and then continues almost zero until 0.47 (2πc/a). The second and third bands both occupy the frequency region from 0.47 to 0.59 (2πc/a), so the transmission becomes larger in this frequency region. From the PBS, we can predict that another sharp drop should take place in a narrow region between 0.59 and 0.61 (2πc/a), which is just the third PBG. A sharp drop after 0.59 (2πc/a) is indeed investigated and then rapid rise after 0.61 (2πc/a) from the ADE-FDTD calculations. The frequency region above 0.61 (2πc/a) and below 1.00 (2πc/a) is occupied by several bands. Hence, the most parts of this frequency region should have non-zero transmission.

Figure 4.

Transmissions calculated by the FDTD method when the Gaussian wave is incident from air into the PhC. The radius of cylinders is 0.2a and the lattice constant a is 100 μm. We set Δx = a/30 and 30 layers in the propagation direction.

Figure 5.

The transmission and reflection when the electromagnetic wave propagates through three media including two interfaces.

Another possible low transmission predicted by the PBS occurs in the vicinity of the intersection between the fifth and sixth bands. In Fig. 4, these two bands intersect at the Γ point of the first Brillouin zone when the frequency is 0.86 (2πc/a). Because the modes in the PhC only occupy k=0 states, the density of states (DOS) is very small at the intersection. On the other hand, from the viewpoint of the effectively refracted index of the PhC, they have very small effectively refracted indices in the vicinity of the intersection. Hence, the multilayer model instead of DOS can be used to explain the extremely low transmissions. In our calculations, the PhC is sandwiched between two homogeneous media. If the PhC is replaced with an effective homogeneous medium (Pei et al., 2011a, 2011b), it can define the effectively dielectric constant. ε_pc, effectively magnetic permeability μ_pc, effectively refracted index n_pc, and effective impedance η_pc in the normally incident case. The effective refracted index and effective impedance are defined as $$n_{pc}=\sqrt{{\epsilon }_{pc}{\mu }_{pc}}$$ and $${\eta }_{pc}=\sqrt{{\mu }_{pc}/{\epsilon }_{pc}}$$for the positive refraction, respectively. For the negative refraction, the refracted index and effective impedance are defined as $$n_{pc}=\sqrt{(-{\epsilon }_{pc})(-{\mu }_{pc})}$$ and $${\eta }_{pc}=\sqrt{\text{(}-{\mu }_{pc}\text{)}/\text{(}-{\epsilon }_{pc}\text{)}}$$ (Engheta & Ziolkowski, 2006). The relation between the effectively refracted index and effective impedance is

$$n_{pc}={\epsilon }_{pc}{\eta }_{pc}{}_{.}$$E33

The transmission and reflection become the problem of multiple scattering by interfaces as shown in Fig. 5 (Moreno, 2002). The total transmitted coefficient of the system consisting of one finite and two semi-infinite media with two interfaces is (Yariv & Yeh, 2002)

$$t={\frac{t_{12}{t}_{23}{e}^{-i{k}_{pc}L}}{1+r_{12}{r}_{23}{e}^{-i2k_{pc}L}}}_{,}$$E34

where k_pc=n_pcω/c and L is the length of the PhC along the y direction. The r₁₂, r₂₃, t₁₂, and t₂₃ are defined as

$$r_{12}={\frac{{\eta }_{pc}-{\eta }_i}{{\eta }_{pc}+{\eta }_i}}_{,}$$E35

$$r_{23}={\frac{{\eta }_t-{\eta }_{pc}}{{\eta }_t+{\eta }_{pc}}}_{,}$$E36

$$t_{12}={\frac{2{\eta }_{pc}}{{\eta }_{pc}+{\eta }_i}}_{,}$$E37

$$t_{23}={\frac{2{\eta }_t}{{\eta }_t+{\eta }_{pc}}}_{,}$$E38

where $${\eta }_i=\sqrt{{\mu }_i/{\epsilon }_i}$$ and $${\eta }_t=\sqrt{{\mu }_t/{\epsilon }_t}$$ are the impedances of the incident region and transmitted region, respectively. The incident and transmitted regions are both air here, so the η_i = η_t. If the effectively refracted index of the PhC is zero in Eq. (33), either ε_pc or μ_pc has to be zero. It deduces that η_pc is zero if μ_pc=0, and η_pc is infinite if ε_pc=0. By calculations, the effectively refracted index at the intersection of two bands is zero as well as the transmission. We can further check this point of view by using the ADE-FDTD method at 0.86 (2πc/a). The result shows that most of the electromagnetic wave cannot pass through the first interface between the incident region and the PhC. According to the above discussion of transmission and Eq. (37), the zero t₁₂ deduces zero effective impedance η_pc. Due to η_pc=0, we have a non-zero ε_pc and a zero μ_pc here.

Figure 6.

The determination of the refracted wave vector and refracted angle by the conservation of wave vector parallel to the interface. The outer circle represents the EFS with frequency ω in air, and the inner one the EFS with frequency ω in the PhC. The line is perpendicular to the interface along the ΓΜ direction.

From the ADE-FDTD calculation in Fig. 4, we find out that the extremely low transmission not only takes place at the intersection, but also extends to the vicinity. They locate at frequencies between 0.80 and 0.88 (2πc/a). To explain it we should calculate the effectively refracted indices in this frequency region. The effectively refracted indices can be determined from EFSs. According to the conservation rule, the incident and the refracted wave vectors are continuous for the tangential components parallel to the interface. Given the incident wave vector and an incident angle, the refracted wave vector and the refracted angle will be determined. How to determine the refracted wave vector and refracted angle by the conservation rule is shown in Fig. 6. The incident and refracted waves are on different sides of the normal line, so it is the positive refraction. By applying Snell’s law, the effectively refracted index can be further determined. But traditional Snell’s law cannot be applied if the EFSs move inward with an increasing frequency (Notomi, 2000). It needs to add minus sign on the effectively refracted index. Figs. 7–10 show the 3D EFSs of the fourth to seventh photonic bands in the first Brillouin zone. It can be seen that the frequency range of 3D EFSs matches the calculation in Fig. 3. They all form a bell shape. Some erect upward and some erect downward. The upward bell usually corresponds to the positive refraction and the downward bell usually corresponds to the negative refraction.

Crosscutting the 3D EFS at a certain frequency reduces to a two-dimensional contour. Thus, we obtain a lot of k_x and k_y at the same frequency drawn in the two-dimensional plane. Each point on the contour is the allowed propagating mode in the PhC for the chosen frequency. In the following, we further discuss the extremely low transmission at frequencies from 0.80 to 0.88 (2πc/a). EFSs of frequencies 0.81, 0.83, and 0.85 (2πc/a) for discussions are shown in Figs. 11–13. A large circle drawn at the center in each figure represents the FES of air at the same frequency. All EFSs of air in these three figures have larger radius than those in the PhC. If we use the conservation rule mentioned before, they all result in the same conclusion that the refracted angle is larger than the incident angle. It also indicates that the absolute value of the effectively refracted index is smaller than 1.0. Because each EFS shrinks with an increasing frequency, the effectively refracted index is negative. Therefore, the negative refraction takes place here. When the frequency is higher, the shape of the EFS of the fifth photonic band is closer to a circle. The circular EFS means that the PhC can be considered as a homogeneous medium at this frequency. The relations between the incident and refracted angles for these three frequencies are shown in Figs. 15(a)-(c). On the one hand, the lower curve of each figure shows the negative refraction, where the refracted angle is defined as negative for convenience. The negative angles are calculated from lines intersecting with the EFSs in the first Brillouin zone as shown in Figs. 11-13. It can be seen that the relation between the incident angle and refracted angle is much like that in a homogeneous medium. On the other hand, the upper curves for a larger incident angle in Fig. 15 (a) and (b) show the normal refraction with positive refracted angle. They are calculated from lines intersecting with the EFSs in the right repeated Brillouin zone as shown in Figs. 11-13.

Figure 7.

The EFS of the fourth photonic band in the first Brillouin zone.

Figure 8.

The EFS of the fifth photonic band in the first Brillouin zone.

Figure 9.

The EFS of the sixth photonic band in the first Brillouin zone.

Figure 10.

The EFS of the seventh photonic band in the first Brillouin zone.

Next, we calculate the effectively refracted index varying with the incident angle only for negative refraction as shown in Fig. 16. From Figs. 15(a)-(c), the normally incident case belongs to negative refraction. It can be found out that the effectively refracted indices of three frequencies 0.81, 0.83, and 0.85 (2πc/a) at incident angle 0 are -0.31, -0.30, and -0.16, respectively. Using Eq. (33), we can calculate these three corresponding effective impedances. But we do not explicitly know the effectively dielectric constants for these three frequencies at normal incidence. However, according to the previous discussion about the normal incidence at 0.86 (2πc/a), we have the conclusion that the effective impedance η_pc is zero with a zero μ_pc and a non-zero ε_pc. Utilizing the similar explanation and a little correction, the effective impedance at 0.85 (2πc/a) should be very close to zero with non-zero μ_pc and ε_pc. The conclusion can also be applied to frequencies 0.81 and 0.83 (2πc/a). As a result, the effective impedances in the frequency range from 0.81 to 0.85 (2πc/a) are very small. Using Eqs. (34)-(38), we obtain extremely low transmissions at frequencies from 0.81 to 0.85 (2πc/a).

6. The internal-field expansion method

In this section, we introduce the internal-field expansion method (IFEM) to calculate the transmission of the finite thickness PhC (Sakoda, 1995a, 1995b, 2004). This method is based on

Figure 11.

The EFS of the fifth photonic band at 0.81 (2πc/a) in the repeated Brillouin zone. The largest circle at the center represents the FES of air with the same frequency.

Figure 12.

The EFS of the fifth photonic band at 0.83 (2πc/a) in the repeated Brillouin zone. The largest circle at the center represents the FES of air with the same frequency.

Figure 13.

The EFS of the fifth photonic band at 0.85 (2πc/a) in the repeated Brillouin zone. The largest circle at the center represents the FES of air with the same frequency.

Figure 14.

EFSs of the fifth photonic band when frequencies are 0.81, 0.83, and 0.85 (2πc/a).

Figure 15.

Refracted angles vs. incident angles calculated from Figs. 11–13 for T = 5 K at (a) 0.81, (b) 0.83, and (c) 0.85 (2πc/a).

Figure 16.

Effectively refracted indices versus incident angles at 0.81, 0.83, and 0.85 (2πc/a).

the internal fields expanded in Fourier series. We consider a two-dimensional PhC composed of a triangular array of air cylinders with radius of r in a dielectric background. The dielectric constants of the cylinders and the background are ε_a and ε_b, respectively. The infinitely extended direction of air holes is parallel to the z-axis. The PhC is infinitely extended in the x-direction and the width of the PhC in the y-direction is finite. Therefore, two dielectric-PhC interfaces exist at the left and right sides of the PhC. a₁ and a₂ are the lattice periods along the x- and y- directions, respectively. The region from the interface to the edge of the nearest cylinder is called the edge region, in which the width is d. The PhC has two edge regions at the left and right sides. The other region including all the cylinders is called the middle region. The total width L of the PhC in the y-direction is L = 2(r + d) + (N_L - 0.5)a₂, where N_L is the periodic number. So the total layers of cylinders are 2N_L. The configuration of the PhC is shown in Fig. 17. The first Brillouin zone is shown at the up-right corner. The region at the left-handed side of the PhC is called the incident region, and that at the right-handed side of the PhC is called the transmitted region. The plane wave in the incident region is incident on the left interface. After propagating through the PhC, the transmitted wave is through the right interface and into the transmitted region.

Figure 17.

The PhC structure with finite length in the y-direction and infinite length in the x-direction.

Because the two-fluid model is only suitable for the currents flowing along the cylinder direction, we only discuss an E-polarized plane wave incident upon the superconductor PhC here. Two interfaces are along the ΓΜ direction of the PhC. The 2D wave vector of the incident wave is denoted by $${\stackrel{\rightharpoonup }{k}}_i$$= (k_isinθ, k_icosθ) = (k_i,x,, k_i,y), where θ is the incident angle and$$k_i=\sqrt{{\epsilon }_i}\omega /c$$. ε_i is the dielectric constant of the incident region, ω is the angular frequency of the incident wave, and c is the light velocity in vacuum. The wave in the incident region consists of the incident plane wave and the reflected Bragg waves. In the transmitted region, the wave is composed of the transmitted Bragg waves. The reflected and transmitted Bragg waves are represented as space harmonics with the wave vector

$$k_{r,x}^n=k_{t,x}^n=k_x^n=k_{i,x}+G_n{}_{,}$$E39

where $$k_{r,x}^n$$ and $$k_{t,x}^n$$ are the wave-vector components parallel to the interface for the Bragg reflected and transmitted waves of order n, respectively, G_n = 2nπ/a₁ is the reciprocal lattice vector corresponding to the periodicity a₁, and n is an integer. Each component of the Eq. (39) is called the nth order phase matching condition. It means that the periodicity along the ΓΜ direction is like a diffraction grating. The wave-vector components of the nth order Bragg reflected and transmitted waves normal to the interface are

$$k_{r,y}^n=\{\begin{array}{l}-\sqrt{k_i^2-{\left(k_x^n\right)}^2}\stackrel{}{{\displaystyle }}if\stackrel{}{{\displaystyle }}{k}_i\ge \left|k_x^n\right|\\ -i\sqrt{{\left(k_x^n\right)}^2-k_i^2}\stackrel{}{}Otherwise\end{array}$$E40

$$k_{t,y}^n=\{\begin{array}{l}\sqrt{k_t^2-{\left(k_x^n\right)}^2}\stackrel{}{{\displaystyle }}if\stackrel{}{{\displaystyle }}{k}_i\ge \left|k_x^n\right|\\ i\sqrt{{\left(k_x^n\right)}^2-k_t^2}\stackrel{}{}Otherwise\end{array}$$E41

Here, $$k_t=\sqrt{{\epsilon }_t}\omega /c$$and ε_t are the transmitted wave vector and dielectric constant of the transmitted region, respectively. The electric fields in the incident region and the transmitted region are given by

$$E_{iz}(x,y)=E_0\exp (i{\stackrel{\leftharpoonup }{k}}_i\cdot \stackrel{\rightharpoonup }{r})+{{\displaystyle \sum _{n=-\infty }^{\infty }{R}_n\exp (i{\stackrel{\rightharpoonup }{k}}_r^n\cdot \stackrel{\rightharpoonup }{r})}}_{,}$$E42

$$E_{tz}(x,y)={{\displaystyle \sum _{n=-\infty }^{\infty }{T}_n\exp (i{\stackrel{\rightharpoonup }{k}}_t^n\cdot \stackrel{\rightharpoonup }{r})}}_{,}$$E43

where E₀, R_n, and T_n are the amplitudes of the electric field of the incident wave, the reflected Bragg wave, and the transmitted Bragg wave, respectively. The electric field inside the PhC satisfies the following equation derived from Maxwell’s equations:

$$L_E\left[E_{PC}(x,y)\right]=\left[\frac{1}{\epsilon (x,y,\omega )}\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}\right)+\frac{{\omega }^2}{c^2}\right]E_{PC}(x,y)=0_{\text{.}}$$E44

Now, we introduce a boundary value function f_E(x,y):

$$f_E(x,y)={\displaystyle \sum _{n=-\infty }^{\infty }\left\{\frac{y}{L}{T}_n+\left(1-\frac{y}{L}\right)\left({\delta }_{n0}{E}_0+R_n\right)\right\}\exp {(i{k}_x^{n}x)}_{,}}$$E45

where δ_nm is the Kronecker’s δ. The boundary value function f_E(x,y) satisfies the boundary conditions at each interface:

$$f_E\left(x,0\right)=E_{iz}\left(x,0\right)and{f}_E\left(x,L\right)=E_{tz}\left(x,L\right)$$E46

Moreover, we define

$${\psi }_E\left(x,y\right)=E_{PhC}\left(x,y\right)-f_E{\left(x,y\right)}_{.}$$E47

If we substitute Eq. (47) into Eq. (44), we have

$$L_E\left[{\psi }_E\left(x,y\right)\right]=-L_E{\left[f_E\left(x,y\right)\right]}_{,}$$E48

$${\psi }_E\left(x,0\right)=\psi \left(x,L\right)=0_{.}$$E49

The problem of unknown E_pc becomes to deal with the internal field. We have to solve Eq. (48) to obtain E_pc field in the PhC. If we expand ψ_E(x,y) and ε^-1(x,y,ω) in Fourier series, we have

$${\psi }_E(x,y)={{\displaystyle \sum _{n=-\infty }^{\infty }{\displaystyle \sum _{m=1}^{\infty }{A}_{nm}\exp (i{k}_x^{n}x)\sin \frac{m\pi }{L}y}}}_{,}$$E50

$$\frac{1}{\epsilon (x,y,\omega )}={{\displaystyle \sum _{n=-\infty }^{\infty }{\displaystyle \sum _{m=1}^{\infty }{\kappa }_{nm}\exp (i{G}_{n}x+k_{m}y})}}_{.}$$E51

Then the electric field in the PhC is expressed as

$$\begin{array}{l}{E}_{PC}(x,y)={\displaystyle \sum _{n=-\infty }^{\infty }\left\{\frac{y}{L}{T}_n+\left(1-\frac{y}{L}\right)\left({\delta }_{n0}{E}_0+R_n\right)\right\}\exp (i{k}_x^{n}x)}\\ +{\displaystyle \sum _{n=-\infty }^{\infty }{\displaystyle \sum _{m=1}^{\infty }{A}_{nm}\exp (i{k}_x^{n}x)\sin \frac{m\pi }{L}y}}\end{array}$$E52

If we substitute Eqs. (45), (50), and (51) into Eq. (48) and compare the independent Fourier components, the equation about coefficients R_n, T_n, and A_nm are obtained as follows:

$$\begin{array}{l}\frac{{\omega }^2}{c^2}{A}_{nm}+{\displaystyle \sum _{n^{\prime }=-\infty }^{\infty }{\displaystyle \sum _{m^{\prime }=1}^{\infty }\left\{{\left(k_x^{n^{\prime }}\right)}^2+{\left(\frac{m^{\prime }\pi }{L}\right)}^2\right\}\left({\kappa }_{n-n^{\prime },m+m^{\prime }}-{\kappa }_{n-n^{\prime },\left|m-m^{\prime }\right|}\right)}}{A}_{n^{\prime }{m}^{\prime }}\\ =-\frac{2}{\pi }{\displaystyle \sum _{n^{\prime }=-\infty }^{\infty }{\left(k_x^{n^{\prime }}\right)}^2{\displaystyle \sum _{m^{\prime }=1}^{\infty }\left({\kappa }_{n-n^{\prime },m+m^{\prime }}-{\kappa }_{n-n^{\prime },\left|m-m^{\prime }\right|}\right)}}\times \frac{{\left(-1\right)}^{m^{\prime }-1}{T}_{n^{\prime }}+R_{n^{\prime }}+{\delta }_{n^{\prime }0}{E}_0}{m^{\prime }}\\ -\frac{2}{\pi }\frac{{\omega }^2}{c^2}\frac{{\left(-1\right)}^{m-1}{T}_n+R_n+{\delta }_{n0}{E}_0}{m}\end{array}$$E53

where κ_n,m is the Fourier coefficient of the inverse of ε(x,y,ω). Next, we calculate the Fourier coefficients of the configuration shown in Fig. 17. In our case, we have

$$\frac{1}{\epsilon (x,y,\omega )}=\frac{1}{{\epsilon }_b}+\left(\frac{1}{{\epsilon }_a}-\frac{1}{{\epsilon }_b}\right){{\displaystyle \sum _{j=1}^4{\displaystyle \sum _{l=-\infty }^{\infty }S\left(\stackrel{\rightharpoonup }{r}-\stackrel{\rightharpoonup }{u}\left(j,l,l^{\prime }\right)\right)}}}_{,}$$E54

where ε_a=ε_s(ω) and $$\stackrel{\rightharpoonup }{u}(j,l,l^{\prime })$$is the center of each cylinder, which are

$$\stackrel{\rightharpoonup }{u}(1,l,l^{\prime })={(l{a}_1,l^{\prime }{a}_2+R+d)}_{,}$$E55

$$\stackrel{\rightharpoonup }{u}(2,l,l^{\prime })={((l+1/2)a_1,(l^{\prime }+1/2)a_2+R+d)}_{,}$$E56

$$\stackrel{\rightharpoonup }{u}(3,l,l^{\prime })={(l{a}_1,-l^{\prime }{a}_2-R-d)}_{,}$$E57

$$\stackrel{\rightharpoonup }{u}(4,l,l^{\prime })={((l+1/2)a_1,-(l^{\prime }+1/2)a_2-R-d)}_{.}$$E58

The S here is the spatial function of the cylinder. The inverse of ε(x,y,ω) now is extended symmetrically to the negative y region (-L≦y≦0) to calculate the Fourier coefficients. Then, we obtain

$$\begin{array}{l}{\kappa }_{n,m}=\frac{1}{2a_{1}L}{\displaystyle {\int }_0^{a_1}dx}{\displaystyle {\int }_{-L}^L\frac{1}{\epsilon (x,y,\omega )}\exp \left[-i\left(G_{n}x+\frac{m\pi }{L}y\right)\right]}dy\\ =\frac{1}{{\epsilon }_b}{\delta }_{n0}{\delta }_{m0}+\frac{1}{2a_{1}L}\left(\frac{1}{{\epsilon }_a}-\frac{1}{{\epsilon }_b}\right){\displaystyle \sum _{j=1}^4{\displaystyle \sum _{l=-\infty }^{\infty }{\displaystyle \sum _{l^{\prime }=0}^{\infty }{\displaystyle {\int }_0^{a_1}dx}{\displaystyle {\int }_{-L}^{L}S\left(\stackrel{\rightharpoonup }{r}-\stackrel{\rightharpoonup }{u}\left(j,l,l^{\prime }\right)\right)e^{-i{\stackrel{\rightharpoonup }{G}}_{nm}\cdot \stackrel{\rightharpoonup }{r}}d{y}_{}}}}}\end{array}$$E59

where$${\stackrel{\rightharpoonup }{G}}_{nm}=(G_n,m\pi /L)$$. After calculating the integration, we obtain

$${\kappa }_{0,0}=\frac{f}{{\epsilon }_a}+{\frac{1-f}{{\epsilon }_b}}_{,}$$E60

$${\kappa }_{2n,2m}={\{\begin{array}{l}\frac{2f}{N_L}\left(\frac{1}{{\epsilon }_a}-\frac{1}{{\epsilon }_b}\right){(}^{-}\cos \frac{a_{2}m\pi }{2L}\frac{\sin (a_{2}m{N}_L\pi /L)}{\sin (a_{2}m\pi /L)}\frac{J_1(G_{2n.2m}R)}{G_{2n.2m}R}\stackrel{}{{\displaystyle }}(m\ne 0)\\ 2f\left(\frac{1}{{\epsilon }_a}-\frac{1}{{\epsilon }_b}\right)\frac{J_1(G_{2n.2m}R)}{G_{2n.2m}R}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}\stackrel{}{{\displaystyle }}(n\ne 0,m=0)\end{array}}_{,}$$E61

$${\kappa }_{2n-1,2m-1}=\frac{2f}{N_L}\left(\frac{1}{{\epsilon }_a}-\frac{1}{{\epsilon }_b}\right){(-1)}^{m-1}\sin \frac{a_2(2m-1)\pi }{4L}\frac{\sin \left[a_2(2m-1)N_L\pi /2L\right]}{\sin \left[a_2(2m-1)\pi /2L\right]}{\frac{J_1(G_{2n-1.2m-1}R)}{G_{2n-1.2m-1}R}}_{,}$$E62

where f is the filling fraction of the superconductor rods in the calculation domain:

$$f={2N_L\pi R^2/a_{1}L}_{.}$$E63

Finally, we want to solve the unknown coefficients, A_nm, R_n, and T_n. Eq. (53) is not enough to solve all unknown coefficients because the number of equations is less than the number of all unknown coefficients. We need other boundary conditions to solve all A_nm, R_n, and T_n. The reminder boundary conditions is the continuity of the x components of the magnetic field, which leads to

$$\pi {\displaystyle \sum _{m=1}^{\infty }m{A}_{nm}}=\left(iL{k}_{r,y}^n+1\right)R_n-T_n+{\delta }_{n0}{E}_0{\left(iL{k}_{i,y}\right)}_{,}$$E64

$$\pi {\displaystyle \sum _{m=1}^{\infty }m{(-1)}^{m-1}{A}_{nm}}=R_n+\left(iL{k}_{t,y}^n-1\right)T_n+{\delta }_{n0}{E}_0{}_{.}$$E65

Follow the calculation processes and consider the boundary conditions for the E-polarized mode, we can determine the unknown coefficients, A_nm, R_n, and T_n. In practical calculation, we restrict the Fourier expansion up to n = ±N and m = M. So there are (2N + 1)M terms in the Fourier expansion. The total number of the unknown coefficients is (2N + 1)(M + 2). From Eqs. (42) and (43) and the boundary conditions, we also obtain (2N + 1)(M + 2) independent equations. Solving these independent equations can obtain these coefficients.

To discuss the reflection and transmission along the y-direction, we can sum up the y-components of the Poynting vectors of all waves and consider the energy flow conservation across these two interfaces. The y component of the wave vector with an imaginary value represents the evanescent wave in the incident region or the transmitted region, so it’s not necessary to consider this kind of wave in the summation. Then, we obtain the following relations for the E-polarized mode:

$${\displaystyle \sum _n\frac{c}{4\pi }\frac{\left|k_{r,y}^n\right|}{k_i}\sqrt{{\epsilon }_i}{\left|R_n\right|}^2+{\displaystyle \sum _{n^{\prime }}\frac{c}{4\pi }\frac{k_{t,y}^{n^{\prime }}}{k_t}\sqrt{{\epsilon }_t}{\left|T_{n^{\prime }}\right|}^2}}=\frac{c\cos \theta }{4\pi }\sqrt{{\epsilon }_i}{\left|E_0\right|}^2{}_{,}$$E66

where n and $$n^{\prime }$$ represent the summation over the Bragg waves with real wave vectors. Then, we can use the Eq. (52) to define the transmission and reflection:

$$T={{\displaystyle \sum _{n^{\prime }}\frac{k_{t,y}^{n^{\prime }}}{k_t}\sqrt{\frac{{\epsilon }_t}{{\epsilon }_i}}{\left|T_n{}^{\prime }\right|}^2}/\cos \theta {\left|E_0\right|}^2}_{,}$$E67

$$R={{\displaystyle \sum _n\frac{k_{r,y}^n}{k_i}{\left|R_n\right|}^2}/\cos \theta {\left|E_0\right|}^2}_{.}$$E68

Eqs. (67) and (68) will be used in Section 7.

7. The transmission calculated by internal-field expansion method

In previous Section, we have introduced the internal-field expansion method to calculate the finite thickness PhC. This method used to calculate the transmission of the electromagnetic wave propagating through the PhC is faster than the FDTD method if the size of the (2N + 1)(M + 2) × (2N + 1)(M + 2) matrix is not very large. In the original references (Sakoda, 1995a, 1995b, 2004), the author concludes that this method can be used for the general two-dimensional PhC. In the following, we use this method to calculate transmissions of the superconductor PhC. Obviously, the boundary conditions of the magnetic field in Eqs. (64) and (65) are no more suitable for the superconductor PhC if superconductor rods are embedded in air. It is the factor that the boundary conditions of the magnetic field in this method are dealt with at the interface between two homogeneous media but not between cylinders and a homogeneous medium. In the latter half part of this section, we try to overcome this problem by adding a virtual edge region. At the beginning, transmissions are directly calculated without adding a virtual edge region. Then we investigate the effect on transmissions after adding it.

Figure 18.

The transmission versus the M value when the frequency is 0.54 (2πc/a) without a virtual edge region.

Figure 19.

Transmissions versus frequency with N=5 and M=600 without a virtual edge region.

The same parameters as those in the previous section are used. The final results of this method are compared with those of the ADE-FDTD method. The wave is supposed to be normally incident on the superconductor PhC, and the propagation direction is along the ΓΚ direction (y-direction) perpendicular to the interface which is along the ΓΜ direction (x-direction). The number of layers along the x-direction is assumed to be infinite. The number of layers along the y-direction is still 30. The lattice constant along the x-direction is a₁ and that along the y-direction is a₂. We choose a₂ = a = 100 μm and a₁ =$$\sqrt{3}$$a₂. The radius of all superconductor cylinders is 0.2a.

First, the width of the edge region d is considered to be zero. N is fixed at 5 and M is determined at the situation when the transmission is convergent. The frequency is chosen at 0.54 (2πc/a). In Fig. 18, it is found out that M=600 is enough for calculation. Then N=5 and M=600 are used to calculate transmissions from 0.01 to 1.00 (2πc/a). The transmissions of the internal-field expansion method have some differences with those of the ADE-FDTD method shown in Fig. 4. The transmissions at frequencies below 0.17 (2πc/a) are not all close to zero. They are more than 0.1 when the frequencies are below 0.035 (2πc/a) and at 0.105 (2πc/a). These results are not coincident with the results of the PBS and the ADE-FDTD method. From the calculations of the PBS before, no propagation modes exist below 0.16 (2πc/a). The calculations of the ADE-FDTD method also confirm this conclusion even if the thickness of the PhC is finite. It means that the internal-field expansion method has some errors at the low frequency region. In the frequency region from 0.17 to 0.33 (2πc/a), transmissions of the internal-field expansion method and the ADE-FDTD method almost match each other except for the value at 0.175 (2πc/a).

The region from 0.33 to 0.47 (2πc/a) is the PBG region. It is found out that this region shifts to the right in the internal-field expansion method. The region extends to 0.53 (2πc/a) in Fig. 19. After the PBG region, the PBS displays two photonic bands existing between 0.47 and 0.595 (2πc/a), and a narrow PBG region between 0.595 and 0.605 (2πc/a). However, the transmissions in Fig. 19 show that high values exist between 0.53 and 0.64 (2πc/a) and nearly zero between 0.64 and 0.65 (2πc/a). In this region, they show a shift about 0.035 (2πc/a) forward higher frequency. From 0.65 to 0.845 (2πc/a), the trend of the transmissions in Fig. 19 is much similar to that of the ADE-FDTD method but the frequency region shifts to the right about 0.045 (2πc/a). In frequencies from 0.80 to 0.895 (2πc/a), the transmissions of the ADE-FDTD method show the third zero-transmission region. This region exists between 0.845 and 0.955 (2πc/a) in Fig. 19, which is 0.02 (2πc/a) larger than that of the ADE-FDTD method.

Next, we try to extend the boundary away from the edge of the cylinder by increasing the width d of the edge region. It is an imaginary boundary between air and the superconductor PhC because the background medium of the superconductor PhC is also air. In fact, such edge region doesn’t exist. The nonzero edge region implies that the results should have something to do with the width of it. Several values of d=0.5a, 1.0a, 1.5a, and 2.0a are calculated and all of them are shown in Figs. 5-20(a)-(d). After comparing all results, we find out that the nonzero edge region only affects transmissions below 0.17 (2πc/a), where the dielectric function in Eq. (9) is negative. The transmissions above 0.17 (2πc/a) are almost unchanged. So it explicitly reveals that this method is not suitable for negative dielectric function.

To summarize, some transmissions of the internal-field expansion method are close to those of the ADE-FDTD method, and some frequency regions have relative shifts between two methods. Roughly speaking, the shift is about 0.06 multiplying the frequency, so it is obvious that all zero-transmission regions below 1.00 (2πc/a) broaden in the internal-field expansion method. In Fig. 21, both results of the internal-field expansion method and the ADE-FDTD method are shown, in which the frequency scale of the ADE-FDTD method is

Figure 20.

Transmissions at (a) d = 0.5a, (b) d =1.0a, (c) d = 1.5a, and (d) d = 2.0a, respectively.

Figure 21.

The transmission versus the M value when the width of the edge region and frequencies are 10.6a and 0.1 (2πc/a), respectively.

multiplied by 1.06. It can be seen that most results of two methods can match each other much better after 0.17 (2πc/a). We find that the calculations of the internal-field expansion method exists some errors, which need to be overcome. It still cannot solve the problem, even the edge region is added in calculations. In order to match the results of the PBS and the ADE-FDTD method, the internal-field expansion method needs to be modified in some way.

8. Conclusion

This study focuses on the transmissions of the two-dimensional superconductor PhC. The PhC, composed of copper oxide high-temperature superconductor rods in a triangular array in air, can be tunable utilizing the temperature modulation. We use the plane wave expansion method introduced in Section 2 to calculate the PBS of it, which is much like a metallic PhC system described by the Drude’s model if the normal conducting current is ignored. The frequency of the fundamental mode in the superconductor PhC is far above zero. It is the reason that the dielectric function is positive when frequency is more than$${\omega }_p^s(x,y)$$, the plasma frequency of the superconducting electron. Because both the electric susceptibility and magnetic permeability have to be either positive or negative, light has the ability to propagate through the superconductor.

Then we use the ADE-FDTD method introduced in Section 4 to calculate the transmission when light is normally incident from air into the superconductor PhC. The results of the ADE-FDTD method are consistent with the PBS and also verify the frequency of the fundamental mode is more than$${\omega }_p^s(x,y)$$, which is 0.17 (2πc/a) in our demonstrated case. It can be seen that the extremely low transmissions correspond to the PBG regions. Some extremely low transmissions exist at the fifth and sixth bands. They can be explained by treating the superconductor PhC as an effective medium sandwiched between two air regions.

Finally, we use the internal-field expansion method developed by Sakoda to calculate the transmission when light is also normally incident from air into the superconductor PhC. It can be found out that transmissions below 0.17 (2πc/a) are not all close to zero. These non-zero transmissions can’t reach convergent values even we use large M in calculations. The results point out that this method can’t be directly applied on the negative dielectric constant media. We try to increase the width of the edge region to overcome this problem, but it is useless. Transmissions above 0.17 (2πc/a) can reach stable values as long as M is large enough. However, the frequency scale has to reduce 1.06 times in order to match the results of the ADE-FDTD method. To sum up, this method is successful to calculate the transmission of the PhC with air cylinders embedded in the homogeneous medium but not suitable very well for the superconductor PhC. One reason is that the boundary conditions between the superconductor PhC and air are not correct. So this method needs modification to obtain correct transmissions.