# Coupling of Photonic Crystal Surface Modes

Authors: Melike Behiye Yücel
Journal: Advances in Condensed Matter Physics
 (2022)

---

## Abstract

Guiding and evanescent coupling properties of surface modes bound to the interfaces of two-dimensional photonic crystals in close proximity are numerically demonstrated. Interacting photonic crystals are composed of silicon pillars in air, where their outermost layers facing each other are annular. Surface modes are identified through supercell band structure computations, while their excitation by the electromagnetic waves through a perpendicular insertion waveguide is demonstrated using finite-difference time-domain simulations. Lifting the degeneracy between the surface modes as a consequence of bringing two identical photonic crystal surfaces to a sufficient distance results in evanescent coupling in a beating manner whose beat length linearly varies between 10 and 20 periods up to a frequency at which both surface modes travel with the same group velocity. The surface mode coupling phenomenon could be employed either to enhance sensitivity or to reduce device size in bio/chemical sensor applications since the effective travelling length of surface waves increases by about 3.5 times due to evanescent coupling.


---

## Body

## 1. Introduction

Photonic crystals (PCs) are structures whose refractive indices are periodically modulated. They possess numerous interesting properties, which have been the subject of scientific and technological attention ever since they were first introduced [1]. PCs display photonic band gap (PBG) similar to their electronic counterparts in which no electromagnetic (EM) waves could be transmitted. EM wave transmission for frequencies lying in the PBG of a PC can be achieved by introducing, for instance, linear defects to form waveguides (WGs) in which EM waves are confined in and decay abruptly in the transverse direction [2–5]. Linear defect WGs have found widespread utility in applications such as wavelength division de/multiplexing [6, 7] and unidirectional light transmission [8, 9]. EM waves in two parallel WGs were shown to couple and decouple and thus enable EM power transfer in between if the dielectric region between the WGs is properly designed [10, 11]. Coupling efficiencies of such PC directional couplers have been improved [12, 13] and successfully adopted in applications such as all-optical switching [14], wavelength selective optical filters [12, 15], power splitters [16, 17], and polarizing beam splitters [18].Yet, another means for EM wave confinement and guiding for frequencies lying in the PBG are realized at the interface between a PC and a homogenous medium in which the PC is embedded [19–26]. EM waves propagating along the PC surface they bound to, namely surface waves (SWs) [19–25], have been extensively studied. SW bands within the PBG should lie below the light line of the homogenous matrix so that SWs do not penetrate either region. Surface modes excited along a corrugated PC surface at the output of a perpendicular insertion WG (IWG) were shown to focus the outgoing EM waves by suppressing diffraction [26, 27] as the modes created at the input surface enhance transmission [28]. The crucial goal of efficiently coupling PC surface and line defect WG modes for all-PC optical circuits is accomplished by modifying the WGs and interfaces [20, 29]. SWs were also employed in achieving directional emission out of a WG in a PC [20]. A 1-to-N beam splitter, which is an important device in all-PC optical circuits, was realized by adjusting the positions of point defects on a modulated surface [30].Microcavity configurations on PC surfaces are usually employed to utilize SWs in lasing and sensing applications [21–23, 31–35] since the surface modes within the semi-infinite crystal surface in which the translational symmetry along the PC boundary is broken turn into resonant states. Refractive index sensitivities ranging from 200nm/RIU to 1500nm/RIU, where RIU stands for refractive index unit, are reported for the systems employing microcavity configurations [33–35]. Moreover, dispersion properties of the SWs are highly sensitive to refractive index variations, which could take place as a result of variations in chemical composition and/or concentration, in the homogenous liquid or gas matrix in which the PC is immersed [20–23, 36]. Such properties of SWs have been in use for sensor applications in which sensitivity values of 93nm/RIU and 117nm/RIU are achieved on the surfaces of square and triangular PC lattices, respectively [37].In addition to the studies mentioned above, there are also studies in which different designs are used in sensor applications [38–40]. Wellenzohn and co-workers proposed and designed a two-dimensional photonic crystal defect waveguide biosensor based on complementary metal oxide semiconductor (CMOS)-compatible silicon-on-insulator technology operating in aqueous solutions at 1.34 μm [38]. The authors compared this operation wavelength with the 1.55 μm and reported a significantly smaller propagation loss at the wavelength they studied. On the other hand, Liu et al. proposed a configuration of D-shaped five-hole photonic crystal fiber-based surface plasmon resonance refractive index sensor with a maximum sensitivity of 20786nm/RIU and a maximum detection range of 1.30 to 1.50 [39]. In a study modeling a highly sensitive sensor based on a slot photonic crystal waveguide facilitating the enhancement of the interaction between the analyte and light for highly sensitive evanescent field absorption sensing in fluids, the symmetry of the photonic crystal waveguide modes with and without a slot has been investigated [40]. As a result, it has been shown that there is an enhancement factor of 7.6 in the evanescent field ratio compared with a slab waveguide with a similar thickness by tuning the slot width.In this article, guiding and evanescent coupling properties of surface modes bound to the interfaces of two semi-infinite PCs with homogenous host when they are brought to close proximity are numerically demonstrated and possible uses of this phenomenon in sensor applications are discussed. Surface modes are identified through supercell band structure (BS) computations, while their excitation by the waves travelling along a perpendicular IWG is demonstrated using two-dimensional (2D) finite-difference time-domain (FDTD) method [41].The study is organized as follows: Band structures of the PCs, the IWG, a single PC surface, and two interacting PC surfaces are presented in Section2. Optimization of surface modification to achieve broadband coupling, the resultant modal properties of coupling surface modes with definite parities, and spectral variation of coupling length and group velocities of coupling modes are also discussed in this section. Demonstration of the coupling between surface modes through FDTD simulations is presented and discussed in Section 3. Brief conclusions are provided inSection 4.


## 2. Design and Dispersion Properties of Photonic Crystal Surfaces

The underlying 2D PC is composed of a square array of silicon pillars in air, whose refractive indices arenSi=3.46 and nAir=1.00, respectively. The PC is laid along the [10] direction (Figure 1). To facilitate surface mode coupling at wavelengths around λ0=1.55μm, lattice constant and radii of the cylinders are set to a=0.5μm and rd=0.1μm=0.2a, respectively. Such high-aspect ratio pillar-type devices can easily be fabricated through techniques, such as deep reactive ion etching [42–44].Figure 1
(Color online) geometry and organization of IWG and coupled surfaces in the square PC. Definitions of geometrical parameters are presented in the enlarged sketches. Vertical and horizontal dashed rectangles represent the supercells for calculating localized modes of respective structures. Horizontal dashed lines denoted bySU and SL represent the regions to probe mode coupling between surfaces.BS of the square PC obtained through the plane-wave expansion (PWE) method [45, 46] with 216 (64×64) plane waves implemented in RSoft Design Groups BandSOLVE software reveals a complete PBG existing along the ΓX direction for a range of wave vector k0≤k≤π/a increased by 0.001 for transverse magnetic (TM) polarized waves, i.e., the electric field being normal to the PC (xy) plane between normalized angular frequencies ωn≡ωa/2πc=0.2815 and 0.4185. Here, ω is the angular frequency and c is the speed of light in free space. Corresponding gap width and gap-over-midgap ratio values are Δωn=0.137 and 39.14%, respectively. Since it is a common practice to couple EM waves into the surface modes through perpendicular IWGs [47], a vertical WG on the left of Figure 1 denoted as IWG is introduced. As it will be discussed in detail below, its width, dIWG, is chosen such that a sufficiently broad air-like band, which covers the surface bands of interest, can be obtained. TM polarized EM wave is launched from the lower end, while the upper end is left open to prevent backscattering. Calculations of the IWG bands, as well as single/coupled surface bands, are carried out within a supercell approach [48] through the PWE method.Localized Bloch modes on PC surfaces are obtained by replacing the scatterers at the outmost layer by annuli [37] with inner and outer radii of ri and r0, respectively (Figure 1). The center-to-center distance between the annuli on the upper and lower surfaces, dSUR, is optimized such that it is not (i) too small to accommodate a single WG band rather than surface bands, or (ii) too large to prevent mode overlap, and thus evanescent coupling. The first 9 elements of the upper PC to the immediate right of the IWG, denoted as the transition region in Figure 1, are not modified such that no surface modes are supported on the upper surface within the transition region to ensure that initially surface modes on only the lower PC surface are excited by the IWG. Thus, only the TM polarized waves from the lower surface can couple to the modes of the upper surface after the transition region.Projected BSs of the IWG for two differentdIWG values, i.e., 1.5a and 4.0a, accompanied by those of a single PC surface interfacing with air host where the surface termination is annular, as well as BSs for two identical such PC surfaces in close proximity, are presented in Figure 2. These BSs are obtained by varying the respective parallel wave vector modulus, k∥, to surface. Figure 2 also presents the projection of bulk bands over the ΓX direction and the light line of air host, i.e., the ω=ck∥ line. In case of the IWG, 9 rows of scatterers on each side, accompanied by an air gap, are employed. For dIWG=1.5a and 4.0a, this corresponds to a 19.5x1 and 22x1 supercell, respectively. A single PC supercell with annular termination involves 11 rows of PC elements, 1 row of annular termination, and a 1.5a thick air layer, resulting in a 1x13.5 supercell. Supercell for coupled surfaces involves 6 rows of PC elements terminated by 1 row of annulus on each side and a dSUR wide gap, giving rise to a 1x14+dSUR/a supercell. Geometries of all supercells are shown in Figure 1, although not drawn to scale for compactness. Choice of such supercells ensures sufficient decay of the WG or surface modes into the PC bulk, sufficient decay of the surface modes into air in the single surface case, as well convergence of the obtained bands. This will be discussed in detail below in terms of mode profiles.Figure 2
(Color online) projected BSs of the IWG fordIWG=1.5a (light blue-colored dashed line with hollow circles) and dIWG=4.0a (purple-colored dashed lines with filled circles), as well as a single PC surface terminated with annular scatterers (magenta-colored dash-dotted lines) and two identical surfaces in proximity (solid lines with triangles where the red-colored filled and the pink-colored empty triangles represent the surface and WG bands, respectively). Inside the PBG, WG bands lie above the light line, while surface bands lie below the light line. The solid shaded areas are the bulk bands, whereas the graded rectangular area represents the full TM PBG. The horizontal thick dash-dotted line (black-colored) indicates the wavelength of interest, i.e., λ0=1.55μm.Figure2 shows that all IWG modes lie above the light line for either narrow (dIWG=1.5a) or wide (dIWG=4.0a) WG width. Only one band for dIWG=1.5a, which covers a small frequency range in the upper portion of the band gap much above the wavelength of interest, i.e., λ0=1.55μmωn0=0.3226, is observed. In contrast, there are three bands for dIWG=4.0a, one of which is close to the light line and varies in close resemblance with it. In fact, such a broad air-like band spanning the whole band gap can be obtained for dIWG>2.23a covering λ0ωn0, as depicted by the horizontal dash-dotted line in Figure 2. Thus, modes of this band can be utilized in exciting surface modes of the adjacent perpendicular surface in Figure 1.To optimize the single PC surface to cover a large bandwidth aroundλ0 and to lie below the light line, the inner and outer radii of the annular elements on the surface, ri and r0, respectively, are treated as optimization parameters. For ease of calculations, ri is fixed to rd, while r0 is varied in the range r0≥ri where the case ri=r0=0.20a corresponds to a hollow WG. A sufficiently broad surface band (SB), which covers λ0 and lies below the light line, can be obtained for r0=0.25a. This band can be utilized to obtain surface modes in a frequency range between the lower band gap edges, i.e., ωn=0.2815 and ωn=0.3624, with a bandwidth of 25.13% around the central frequency. Air-like band of the IWG covers the whole frequency range of the SB, thus rendering the excitation of each mode within the SB.Surface modes for adjacent PCs with annular surface elements are investigated using the supercell denoted as the coupling region in Figure1. Here, the only relevant control parameter is restricted to dSUR and it determines the degree of mode overlap between the neighboring PC surfaces. A too large value gives rise to negligible mode overlap between neighboring PC surfaces, whereas a small value results in a significant overlap and thus the WG modes instead of coupling surface modes as discussed above. Figure 2 shows that the SB of the single PC surface lying between ωn=0.2815 and 0.3624 (dash-dotted curve) splits into two SBs (solid lines accompanied by filled triangles) for a moderate dSUR of 1.5a. One of the SBs is closer to the dielectric bands of the PC, whereas the other is to the light line. Although another band lying between ωn=0.3826 and 0.4533 (solid line accompanied by empty triangles) in Figure 2 exists, it is out of scope, since it is a WG band whose modes are air-guided rather than being surface modes.The discrepancy between the two-split SBs of the neighboring surfaces is primarily determined bydSUR. When dSUR is much larger than a, SBs of the neighboring surfaces collapse into a single degenerate SB of the isolated PC surface. When dSUR, however, is sufficiently small, e.g., dSUR=1.5a, degeneracy is lifted and the two-split SBs cover a common frequency range between the top of the bulk dielectric bands, i.e., ωn=0.2815, and the top of the lower SB at ωn=0.3562 for reasons to be discussed below. The overlap region covers a frequency range of Δωn=0.0747, or equivalently a bandwidth of 23.43%, which significantly diminishes to Δωn=0.0213 as dSUR increases to 2.0a. Moreover, the corresponding WG band depicted by the solid line with hollow triangle in Figure 2 starts overlapping the SBs for dSUR>1.65a. Thus, dSUR=1.5a is the optimal value for SB coupling.Modal properties of the single PC surface band whereri=0.2a and r0=0.25a, along with those of the SBs of the adjacent PC surfaces for dSUR=1.5a at ωn0=0.3226, are presented in Figure 3. It is clearly seen that EM energy is concentrated mainly inside the annular region of the isolated PC surface (Figure 3(a)). The field strength rapidly drops toward both the air and PC bulk. Since the perpendicular electric field component, Ez, is negligible inside the PC bulk, as demonstrated by the plot of the cutline along the vertical bisector of the supercell in Figure 3(a), the adopted supercell size 1x13.5 provides a well-converged SB in Figure 2.Figure 3
(Color online) modal properties of the SB of an isolated PC surface (a) and two interacting surfaces (b) atωn0=0.3226 with annular ends. The plots on the right depict the variation of Ez along the corresponding vertical bisectors of supercells.
(a)(b)The mode profiles in Figure3(b) demonstrate that the surface mode splits into an antisymmetric mode and a symmetric mode represented by filled down and up triangles, respectively, in Figure 2. Antisymmetric mode is such that there is a node Ez=0 midway between the neighboring surfaces, whereas there is significant overlap Ez>0 for the symmetric mode (Figure 3(b). This is remarkably similar to the case of coupling two parallel PC WGs [49]. This resemblance can be utilized in obtaining unidirectional light transmission through coupling surface modes. Inspection of the surface modes at a given k∥ value reveals that antisymmetric modes are always at higher angular frequencies.An incident TM polarized wave through the IWG in Figure1 first couples to the surface mode in the transition region. Thereafter, it couples back and forth between the antisymmetric and symmetric modes until the exit on the right side of the system in Figure 1. A quantitative measure of evanescent coupling between the surface modes is the normalized spatial coupling period, denoted as lc=Lc/a, where Lc is the coupling length determined by the difference between the k values of the modes, Δk∥, as [10, 49] follows:(1)Lc=πΔk∥.Another useful quantity is the beat length, which denotes the full spatial cycle for the energy to couple back and forth once, defined aslb=2lc.Figure4(a) demonstrates that lc varies roughly between 5 and 10 for ri=0.2a, r0=0.25a and dSUR=1.5a. It first increases monotonically for ωn<0.345 and then decreases rapidly. lc=8.75lb=17.50 at ωn0=0.3226 is remarkably short and can lead to design of very compact devices out of the coupling mechanism. lc varies linearly for ωn<0.345 in which a linear fit, R2=0.9996, reveals that the rate of change with respect to ωn is 82.38, or equivalently 0.0219 per THz of frequency shift, or 0.269 per percent variation around ωn0. This is a measurable shift, and tracking the variation of lc can be used for sensing purposes.Figure 4
(Color online) variation of the normalized coupling length, accompanied by its first derivative (a) and normalized group velocities (b) of the surface modes with frequency. Inset in (a) presents a closer look forωn<0.34. Dash-dotted vertical lines indicate the frequencies at which FDTD simulations are performed, whereas the hatched rectangular areas denote the frequency ranges over which surface mode coupling is not observed.
(a)(b)A better approach for evaluating the possibility of using surface mode coupling for sensing is considering the first derivative oflc with respect to λ, i.e., dlc/λ=dλ, the result of which is also presented in Figure 4(a). Although it takes small negative values for frequencies below ωn=0.345 (equivalently λ>1.449μm), it sharply increases for higher frequencies close to the symmetric SB top, or k∥=π/a. As lc also drops sharply around the band edge, the surface mode coupling mechanism can be used in sensing devices.The coupling behavior of symmetric and antisymmetric surface modes is closely related to group velocities of the modes, which can be determined through the dispersion curves in Figure2. As the slopes of the two-split surface bands in Figure 2 are generally different, modes with different parity propagate along PC surfaces with different speeds. This, in turn, results in a phase difference between coupling symmetric and antisymmetric modes. Figure 4(b) reveals that (normalized) group velocity υg/c of the symmetric mode is higher than that of the antisymmetric mode up to a frequency ωne=0.3468 at which they become equal. For frequencies higher than ωne, the case is reversed (Figure 4(b)). No coupling is expected for ωn>0.3562 where only the antisymmetric mode is excited up to the band top. Thus, three frequencies, i.e., ωn=0.3226, 0.3468, and 0.3566, are chosen for FDTD simulations, which will be discussed in the next section to demonstrate out-of-phase, in-phase coupling, and single-mode excitation, respectively.


## 3. Demonstration of Surface Mode Coupling

Dynamic visualization of surface mode coupling is carried out through FDTD simulations [41] by RSoft Design Group’s FullWAVE software. Grid size in FDTD computations is chosen as dx=dy=a=32, while the time step Δt is such that cΔt=dx/2 to satisfy the Courant–Friedrichs–Lewy stability criterion [50]. The computational domain is surrounded by perfectly matched layer-absorbing boundaries [51]. Surface mode coupling is demonstrated through probing the electric field Ez distribution over the linear patches denoted as SU and SL in Figure 1.FDTD simulation of surface mode coupling atωn0=0.3226 is presented in Figure 5. The TM polarized wave is launched from the bottom of the IWG couples evanescently into the lower surface in the transition region. When the surface wave reaches the coupling region, it couples back and forth between the upper and lower surfaces. The coupling region is such that 2 complete beat cycles are covered (Figure 5).Figure 5
(Color online) FDTD simulation of surface mode coupling atωn0=0.3226 (top panel). The middle panel is a close-up view of the transition and coupling regions, whereas the bottom panel presents the variation of Poynting vector magnitude along the horizontal dash-dotted cutlines. The plot to right on the middle panel depicts the variation of S on the vertical dash-dotted cutline at the output. The thick hollow arrow in the top panel indicates the initially launched wave direction, whereas the shaded area in the bottom panel denotes the coupling region.The evanescent nature of SWs excited on the lower PC surface in the transition region in Figure5 reveals that decay ratios are different on the two sides of the lower PC surface, due to significantly different media in either directions, e.g., PC bulk and air. The decay length of the SW in air is considerably larger than that into the PC bulk, as evidenced by the mode profile in Figure 3(a). This, in turn, suggests that SWs on the lower PC could well excite corresponding modes on the upper PC while they travel along the x-axis.For visual purposes, the area within the dashed rectangle in the top panel of Figure5, which covers the transition and coupling regions, is enlarged and presented in a different color scale covering a narrower Ez range since the wave components propagating in the IWG have significantly larger magnitude. The middle panel in Figure 5 clearly demonstrates the coupling of surface modes between adjacent surfaces. The bottom panel in Figure 5 is a depiction of the variation of the magnitude of the Poynting vector (S) along the cutlines SU and SL in Figure 1, denoted by the horizontal arrows to the left of the middle panel of Figure 5. Only the surface modes of the lower PC are excited throughout the transition region, as S is almost zero on the upper surface, as depicted in the plot of S along the surfaces. In contrast, evanescent coupling rapidly takes place when the coupling region is encountered (Figure 5). The beat length denoted by the arrow between the two vertical dashed lines in Figure 5 is calculated to be lb=16.2, which is in very good agreement with the value (lb=16.8) predicted from Figure 4(a).The middle and bottom panels of Figure5 demonstrate that the coupling symmetric and antisymmetric modes are out of phase, as discussed in the previous section. It is more clearly visible on the plot of S that the peaks for the upper and lower PC surfaces never coincide along the surfaces, bottom panel of Figure 5. The plot of S at the output of the coupling PC surfaces, denoted by the vertical dash-dotted line in the middle panel of Figure 5, demonstrates that wave output is primarily from the lower surface at ωn0, as the coupling region length is such that as two complete beats take place over the two PC surfaces.In general, the ratio of collected power between the surfaces reveals frequency dependence as depicted in Figure6, which presents the variation of S with vertical position at the output side and frequency. The maximum frequencies of the output power tightly depend on the WG width dSUL=2y between the upper and lower PC surfaces with the lattice constant a. Maximum power peaks are obtained at y/a = 0.75 for the coupling surface modes. These modes are travelling electromagnetic waves, which are bound to the interface between PCs and uniform host (air). Figure 6 shows that most of the EM energy is output from the upper surface at ωn=0.302, while the situation converses at ωn=0.324. There is roughly an oscillatory behavior with frequency away from the band edges. For sensor applications, the device can be operated at frequencies for which wave output is almost totally from one of the surfaces such that slight variations in the refractive index of the host result in disruption of this large contrast.Figure 6
(Color online) variation of the modulus of the Poynting vector with vertical position and frequency at the output interface of the interacting PC surfaces with air.Figure7 demonstrates SW propagation at ωne=0.3468 and ωn=0.3566. Both SWs on the two surfaces in Figures 7(a) and 7(b) travel in phase with and without beating (evanescent coupling), respectively. This is clearly visible in the plots of S along the PC surfaces calculated along the horizontal dash-dotted lines in Figure 7(a) where peaks of S for the upper and lower surfaces occur at the same positions at all x values, bottom panels in Figure 7. Such in-phase operation could be employed for sensing either to reduce the sensor size or to enhance the sensitivity since the effective interaction length between the SW and the gas phase is increased by about 3.5 times.Figure 7
(Color online) simulation results forωn=0.3468 (a) and ωn=0.3566 (b) along with the variation of Poynting vector magnitude along the surfaces (bottom panels) and at the output side (to the right of FDTD plots).
(a)(b)It should be noted from Figures5 and 7(a) that SWs exit the PC surfaces unevenly for either ωn=0.3266 or ωn=0.3468. This can easily be deduced from the variation of S over the output side (vertical dash-dotted line in Figures 5 and 7(a)). However, SWs at ωn=0.3566 leave in an antisymmetric manner with a zero field strength in the midway (Figure 7(b)) as was pointed out in the discussion of Figure 4(b). There is no coupling at and above this frequency, where the propagating mode is solely the antisymmetric mode in Figure 3(b).It should be noted that peak positions of theS curves for neither upper nor the lower surface can be fitted by a smoothly varying envelope for ωn=0.3468 (Figure 7(a)). In contrast, the corresponding curves in Figure 5 are more likely to be fitted by Gaussian-like curves for ωn=0.3226. Figure 7(a) suggests that there is interplay of two beating phenomena occurring at different spatial periods. One is obviously due to evanescent coupling. The other mechanism can ascribe to beating due to the finite length of the PC surfaces where two equivalent counter-propagating modes with complementary wave numbers with respect to k∥=π/a interfere [49]. Due to this secondary beating, the beat length in Figure 7(a) is calculated as lb=14.0a, a value significantly smaller than the one expected from Figure 4(a), i.e.,lb=20.0a. The secondary beating length increases as k∥ approaches π/a, where Δk∥ of the counter-propagating modes decreases.


## 4. Conclusion

In conclusion, evanescent coupling between surface modes of parallel two-dimensional photonic crystals is numerically demonstrated. While there exists a single degenerate surface band for two identical surfaces infinitely apart, the degeneracy is lifted as the surfaces are brought together at a length scale compared with lattice periodicity. This, in turn, gives rise to the occurrence of two localized surface bands within the band gap with opposite definite parities. The discrepancy of these modes in the reciprocal space suggests that effective coupling length increases linearly with frequency up to a point where the group velocities of the two interacting modes become equal, after which a rapid fall of the coupling length is observed. Finite-difference time-domain simulations demonstrate surface mode coupling at length scales in agreement with the values predicted from the band structure calculations. Coupling modes propagate along the two surfaces in phase if the group velocities are equal. This can be utilized to achieve highly sensitive compact refractive index sensors, as the effective interaction length increases significantly due to staggering by evanescent coupling.



---

*Source: 8947410-2022-08-10.xml*