## 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*ω* 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).

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):

where

*λ*(*x, y*) is the London penetration depth, *ε*_{1}(*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:

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

where *ε*_{s}(*x, y*, *ω*) is the effective dielectric function given by

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 *T*_{c} and *λ*_{0} 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,

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

where the function

where the Fourier coefficient

where the integral is now over the 2D unit cell *a*_{c} is the area of the unit cell in the PhC. Eq. (12) can be expressed as

where *f* is the filling fraction. For the triangular superconductor PhC, the filling fraction of the superconductor rod in a unit cell is

where *ω*. The set of equations for coefficients

Rearranging Eq. (16) that we have

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.85}Pb_{0.35}Sr_{2}Ca_{2}Cu_{3.1}O_{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* *T*_{c} is satisfied, the dependence of the plasma frequency on the temperature is given by (Zhou, 1999)

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 *ε*_{1} = 12 (Shibata & Yamada, 1996). When *T* = 5 K, we obtain^{-1}.

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.2*a*, 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:

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*

These equations are iteratively solved in a leapfrog manner, alternating between computing th*e* *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

where

Another time-dependent Maxwell’s curl equation is

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 *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:

Solving Eqs. (26) and (27) for

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

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:

Thus, the ADE-FDTD algorithm for calculating dispersive media has three processes. Starting with the assumed known values of

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”:

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*_{1} and that along *y-*direction *a*_{2}. We choose *a*_{2} = *a* = 100 *μ*m and *a*_{1} =*a*_{2}. 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*/2*c*. 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 *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.

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

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)

where *k*_{pc}=*n*_{pc}*ω*/*c* and *L* is the length of the PhC along the *y* direction. The *r*_{12}, *r*_{23}, *t*_{12}, and *t*_{23} are defined as

where *η*_{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*_{12} deduces zero effective impedance *η*_{pc}. Due to *η*_{pc}=0, we have a non-zero *ε*_{pc} and a zero *μ*_{pc} here.

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.

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

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*_{1} and *a*_{2} 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_{2}, 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.

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 _{i}sinθ, k_{i}cosθ) = (k_{i,x,,} k_{i,y}), where θ is the incident angle and_{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

where _{n} = 2nπ/a_{1} is the reciprocal lattice vector corresponding to the periodicity a_{1}, 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

Here, _{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

where E_{0}, 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:

Now, we introduce a boundary value function f_{E}(x,y):

where δ_{nm} is the Kronecker’s δ. The boundary value function f_{E}(x,y) satisfies the boundary conditions at each interface:

Moreover, we define

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

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

Then the electric field in the PhC is expressed as

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:

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

where ε_{a}=ε_{s}(ω) and

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

where

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

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

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:

where n and

## 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.

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_{1} and that along the y-direction is a_{2}. We choose a_{2} = a = 100 μm and a_{1} =_{2}. 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

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

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

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.