Coupling of Photonic Crystal Surface Modes

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 identied through supercell band structure computations, while their excitation by the electromagnetic waves through a perpendicular insertion waveguide is demonstrated using nitedierence time-domain simulations. Lifting the degeneracy between the surface modes as a consequence of bringing two identical photonic crystal surfaces to a sucient 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. e surface mode coupling phenomenon could be employed either to enhance sensitivity or to reduce device size in bio/ chemical sensor applications since the eective travelling length of surface waves increases by about 3.5 times due to evanescent coupling.


Introduction
Photonic crystals (PCs) are structures whose refractive indices are periodically modulated. ey possess numerous interesting properties, which have been the subject of scienti c and technological attention ever since they were rst 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 con ned in and decay abruptly in the transverse direction [2][3][4][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 e ciencies of such PC directional couplers have been improved [12,13] and successfully adopted in applications such as all-optical switching [14], wavelength selective optical lters [12,15], power splitters [16,17], and polarizing beam splitters [18].
Yet, another means for EM wave con nement 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][20][21][22][23][24][25][26]. EM waves propagating along the PC surface they bound to, namely surface waves (SWs) [19][20][21][22][23][24][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 di raction [26,27] as the modes created at the input surface enhance transmission [28]. e crucial goal of e ciently 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][22][23][31][32][33][34][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 200 nm/RIU to 1500 nm/RIU, where RIU stands for refractive index unit, are reported for the systems employing microcavity configurations [33][34][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][21][22][23]36]. Such properties of SWs have been in use for sensor applications in which sensitivity values of 93 nm/RIU and 117 nm/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][39][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]. e 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 fivehole photonic crystal fiber-based surface plasmon resonance refractive index sensor with a maximum sensitivity of 20786 nm/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]. e study is organized as follows: Band structures of the PCs, the IWG, a single PC surface, and two interacting PC surfaces are presented in Section 2. 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 in Section 4.

Design and Dispersion Properties of Photonic
Crystal Surfaces e underlying 2D PC is composed of a square array of silicon pillars in air, whose refractive indices are n Si � 3.46 and n Air � 1.00, respectively. e 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 r d � 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][43][44].
BS of the square PC obtained through the plane-wave expansion (PWE) method [45,46] with 2 16 (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 k (0 ≤ 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-overmidgap 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, d IWG , 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 r i and r 0 , respectively (Figure 1). e center-to-center distance between the annuli on the upper and lower surfaces, d SUR , 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. e 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. us, 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 different d IWG 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. ese 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 d IWG 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 d SUR wide gap, giving rise to a 1x(14 + d SUR /a) supercell. Geometries of all supercells are shown in Figure 1, although not drawn to scale for compactness. Choice of such supercells ensures su cient decay of the WG or surface modes into the PC bulk, su cient decay of the surface modes into air in the single surface case, as well convergence of the obtained bands.
is will be discussed in detail below in terms of mode pro les. Figure 2 shows that all IWG modes lie above the light line for either narrow (d IWG 1.5a) or wide (d IWG 4.0a) WG width. Only one band for d IWG 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 d IWG 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 d IWG > 2.23a covering λ 0 (ω n0 ), as depicted by the horizontal dash-dotted line in Figure 2. us, 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, r i and r 0 , respectively, are treated as optimization parameters. For ease of calculations, r i is xed to r d , while r 0 is varied in the range r 0 ≥ r i where the case r i r 0 0.20a corresponds to a hollow WG. A su ciently broad surface band (SB), which covers λ 0 and lies below the light line, can be obtained for r 0 0.25a. is 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 Figure 1. Here, the only relevant control parameter is restricted to d SUR 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 signi cant 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 (dashdotted curve) splits into two SBs (solid lines accompanied by lled triangles) for a moderate d SUR 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 airguided rather than being surface modes. e discrepancy between the two-split SBs of the neighboring surfaces is primarily determined by d SUR . When d SUR is much larger than a, SBs of the neighboring surfaces collapse into a single degenerate SB of the isolated PC surface. When d SUR , however, is su ciently small, e.g., d SUR 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. e overlap region covers a frequency range of Δω n 0.0747, or equivalently a bandwidth of 23.43%, which signi cantly diminishes to Δω n 0.0213 as d SUR 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 d SUR > 1.65a. us, d SUR 1.5a is the optimal value for SB coupling.
Modal properties of the single PC surface band where r i 0.2a and r 0 0.25a, along with those of the SBs of the adjacent PC surfaces for d SUR 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)). e eld strength rapidly drops toward both the air and PC bulk. Since the perpendicular electric eld component, E z , 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.
e mode pro les in Figure 3(b) demonstrate that the surface mode splits into an antisymmetric mode and a symmetric mode represented by lled down and up triangles, respectively, in Figure 2. Antisymmetric mode is such that there is a node (E z 0) midway between the neighboring surfaces, whereas there is signi cant overlap (|E z | > 0) for the symmetric mode (Figure 3(b).
is is remarkably similar to the case of coupling two parallel PC WGs [49].
is 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 Figure 1 rst couples to the surface mode in the transition region. ereafter, 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 l c L c /a, where L c is the coupling length determined by the di erence between the k values of the modes, Δk , as [10,49] follows: Another useful quantity is the beat length, which denotes the full spatial cycle for the energy to couple back and forth once, de ned as l b 2l c . Figure 4(a) demonstrates that l c varies roughly between 5 and 10 for r i 0.2a, r 0 0.25a and d SUR 1.5a. It rst increases monotonically for ω n < 0.345 and then decreases rapidly. l c 8.75 (l b 17.50) at ω n0 0.3226 is remarkably short and can lead to design of very compact devices out of the coupling mechanism. l c varies linearly for ω n < 0.345 in which a linear t, R 2 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 .
is is a measurable shift, and tracking the variation of l c can be used for sensing purposes.
A better approach for evaluating the possibility of using surface mode coupling for sensing is considering the rst derivative of l c with respect to λ, i.e., dl c /λ 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 l c also drops sharply around the band edge, the surface mode coupling mechanism can be used in sensing devices. e 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 Figure 2. As the slopes of the twosplit surface bands in Figure 2 are generally di erent, modes with di erent parity propagate along PC surfaces with di erent speeds.
is, in turn, results in a phase di erence 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. us, 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-ofphase, in-phase coupling, and single-mode excitation, respectively.

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]. e computational domain is surrounded by perfectly matched layer-absorbing boundaries [51]. Surface mode coupling is demonstrated through probing the electric eld (E z ) distribution over the linear patches denoted as S U and S L in Figure 1.
FDTD simulation of surface mode coupling at ω n0 0.3226 is presented in Figure 5. e 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. e coupling region is such that 2 complete beat cycles are covered ( Figure 5). e evanescent nature of SWs excited on the lower PC surface in the transition region in Figure 5 reveals that decay ratios are di erent on the two sides of the lower PC surface, due to signi cantly di erent media in either directions, e.g., PC bulk and air. e decay length of the SW in air is considerably larger than that into the PC bulk, as evidenced by the mode pro le in Figure 3(a). is, 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 Figure 5, which covers the transition and coupling regions, is enlarged and presented in a di erent color scale covering a narrower E z range since the wave components propagating in the IWG have signi cantly larger magnitude. e middle panel in Figure 5 clearly demonstrates the coupling of surface modes between adjacent surfaces. e 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). e beat length denoted by the arrow between the two vertical dashed lines in Figure 5 is calculated to be l b 16.2, which is in very good agreement with the value (l b 16.8) predicted from Figure 4(a). e middle and bottom panels of Figure 5 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. e 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     Advances in Condensed Matter Physics 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 7 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. is 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 e ective interaction length between the SW and the gas phase is increased by about 3.5 times.
It should be noted from Figures 5 and 7(a) that SWs exit the PC surfaces unevenly for either ω n 0.3266 or ω n 0.3468. is 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 eld strength in the midway (Figure 7(b)) as was pointed out in the discussion of Figure 4 ere 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 the |S| curves for neither upper nor the lower surface can be tted 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 tted by Gaussian-like curves for ω n 0.3226. Figure 7(a) suggests that there is interplay of two beating phenomena occurring at di erent spatial periods. One is obviously due to evanescent coupling. e other mechanism can ascribe to beating due to the nite 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 l b 14.0a, a value signi cantly smaller than the one expected from Figure 4(a), i.e.,l b 20.0a. e secondary beating length increases as k approaches π/a, where Δk of the counter-propagating modes decreases.

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 in nitely apart, the degeneracy is lifted as the surfaces are brought together at a length scale compared with lattice periodicity. is, in turn, gives rise to the occurrence of two localized surface bands within the band gap with opposite definite parities. e 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. is can be utilized to achieve highly sensitive compact refractive index sensors, as the effective interaction length increases significantly due to staggering by evanescent coupling.

Data Availability
No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Conflicts of Interest
e author declares that there are no known conflicts of financial interest or personal relationships that could have appeared to influence the work reported in this study.