Doubly Resonant Optical Periodic Structure

Periodic structures are well known in various branches of physics for their ability to provide a stopband. In this article, using optical periodic structures we showed that, when a second periodicity – very closed to the original periodicity is introduced, large number of states appears in the stopband corresponding to the first periodicity. In the limit where the two periods matches, we have a continuum of states, and the original stopband completely disappears. This intriguing phenomena is uncovered by noticing that, regardless of the proximities of the two periodicities, there is an array of spatial points where the dielectric functions corresponding to the two periodicities interfere destructively. These spatial points mimic photonic atoms by satisfying the standards equations of quantum harmonic oscillators, and exhibit lossless, atom-like dispersions.

Doubly resonant systems have compelling physical properties resulting from the interference effects. A three level atomic system with a Λ configuration has two resonant transitions, and by an appropriate coherent driving, we can generate steep positive (normal) [1][2][3][4] and steep negative (anomalous) 5-7 optical dispersions. Such steep dispersions can be exploited to create novel systems of slow light and subluminal light without any violation in the Einstein's causality principle [8][9] . Steep positive dispersions are usually obtained in the three level systems using an electromagnetic induced transparency (EIT) setup [1][2][3] , and there have been many proposals to mimic such configuration using plasmonic [10][11] and optical [12][13][14][15] double resonances. These mock versions work based on the coherent interference effects in the light scattering, and allow tuning of the positive dispersions via modifications in the geometrical structures. They also have been shown to possess scattering dark states 16 and superscattering states 17,18 . In this article, we demonstrate the intriguing optical properties of a new paradigm of doubly resonant systems that exploits structures with both short and long range spatial periodicities, and exhibiting two closely spaced Bragg resonances.
Typical periodic structures are single -period, structures (SPSs), and they exhibit stopbands [i.e., spectral regions for which wave propagations are forbidden]. The physical principle behind this stopband formation is the Bragg resonance of the SPS. Waves with frequency in the vicinity of Bragg resonance frequency, will experience a strong Bragg reflectivity, and therefore is unable to penetrate the bulk of the SPS.
In optics a SPS can be created by mean of a periodic variation of the dielectric constant with a fixed spatial period, a. We can again modulate the dielectric profile of this SPS, slowly and periodically, with a longer spatial period, a s [19][20][21] . This new periodic structure which exhibits rapid, short range periodicity (a) and slow, long range periodicity (a s ) is defined as a dual periodic structure (DPS). Intuitively, one can expect in the limit of a very large a s , the slow modulation vanishes, and consequently a DPS reduces to a SPS. However, a DPS in this limit does not fit into this simple intuition.
In the Fourier spectra, the dielectric function of a weakly modulated SPS, will exhibit one frequency peak at the fundamental spatial frequency G = 2π /a. However, for the DPS, due to the slow dielectric modulation, we will see a group of closely spaced peaks around the fundamental frequency. Assuming a DPS with only two of such closely spaced peaks (i.e., a structure with double Bragg resonances), the dielectric function can be casted as, where r is a number close to 1. In Eqn. 1, ε ε  1 0 , and for a simplicity, we assumed the strengths of the two closely spaced harmonics to be equal (i.e., the amplitudes of two cosine functions in Eqn. 1 are equal). The conservation of translational symmetry in DPS requires ε ε ( ) = ( + ) x x a s , and using Eqn. 1 it can be shown that this demands Scientific RepoRts | 6:20590 | DOI: 10.1038/srep20590 / a a s to be the least integer multiple of r. The least integer multiple exist, only if r is rational. Assuming a rational r, and taking its' least integer multiple as R, we have = a Ra s . In a DPS, the mixing of the two harmonics, G and G/r creates the spatial "beats" in the dielectric function at a longer spatial scale. The length of the beat ( = a Ra s ) is longer when the spacing between the two harmonics is closer. The closest allowed proximity between the two harmonics, G and G/r is one reciprocal lattice vector of the DPS, g = π/a 2 s . Any spacing lesser than g, is symmetrically forbidden, and hence will break the translational symmetry of the DPS. Assuming > r 1, and the spacing, G − G/r = g, it is easy to show that the rational form of For a SPS, r = 1, and from Eqn.
. For the DPS, a direct substitution of = r 1 for the limit r → 1 in Eqn. 1 leads to a plausible inference that the DPS should be identical to the SPS in the limit r → 1, and therefore recovers the original stopband of the SPS. However, the constraint r = R/(R − 1) for the DPS prevents the direct substitution of = r 1 for the limit r → 1 in Eqn. 1. The flawless method of analysing the limit r → 1 in DPS is by letting R to take a huge integer value. As an illustration, Fig. 1(a) shows a sketch of ε( ) x with the unit cell from − / a 2 s to / a 2 s is highlighted in blue. Figure 1(b) depicts the evolution of the unit cell as R is increased to a huge integer value. As we can see from Fig. 1(a,b), regardless of the proximity of r to 1, the dielectric function of the DPS is topologically different from the dielectric function of the SPS -which is an unmodulated cosine function. As a signature difference, in DPS, the destructive interference between the two cosine waves in Eqn. 1 creates an array of spatial points (i.e., the green dots in Fig. 1) that are shielded from the effect of the rapid dielectric modulation with the period a. As we shall illustrate, this array of spatial points mimics an array of photonic atoms by satisfying the standard equations of quantum harmonic oscillators 22 . These photonic atoms create edge states (at the edge of the DPS unit cell) that closes the stopband due to the rapid dielectric modulation, despite the limit r → 1. The dispersions in the vicinity of these photonic atoms, are strongly anomalous (i.e., a steep negative dispersion), and very much similar to the dispersions in mediums with inverted populations [23][24] and gain doublets [5][6][7] .

DPS as a Metamaterial Cavity
For the purpose of the numerical illustration, throughout this article, we use ε = .
2 56 0 , and ε = . 0 16 1 . An optical structure with such dielectric constants, and the dielectric profile as in Eqn. 1 with a large R, can be realized in many different ways. Some of the techniques include, fabrication of porous silicon via electrochemical anodization with varying current density 25 , deposition of a dual periodic multilayer using a logical combination technique 19 , holographic interferometry, that make uses laser beam interferences on photosensitive materials 26 , and the deposition of silicon oxynitride with varying stoichiometry of oxygen and nitrogen [27][28] .
Using a plane wave expansion method 29 , we can solve the dispersion relation, ω λ = / a versus k, where ω, λ and k are the normalized frequency, freespace wavelength and wavevector, respectively. Firstly, consider the limiting case, when r → 1. As r = R/(R − 1), the dispersion curve of the DPS in this limit can be obtained asymptotically by increasing R to a huge integer value, using an extended zone scheme 30 . For a very large R, the dispersion of the DPS converges to a continuous curve, Ω( ) k , shown in Fig. 2(a) [blue curve]. In the same diagram, we have also plotted the dispersion curves of the SPS, and a homogenous medium with dielectric constant ε 0 . When r → 1, from the direct substitution of r = 1 in Eqn. 1, one can expect the dispersion of the DPS to be identical to the dispersion of the SPS. However, this is only true for wavevectors far from G/2. For wavevectors far from G/2, both SPS and DPS behave as linear homogenous materials of dielectric constants ε 0 . Near k = G/2, the dispersion relation of the DPS is remarkably different from the dispersion relation of the SPS despite r → 1 [ Fig. 2(a)]. In the SPS, the . The anomalous dispersion, and the resulting superluminal group index of the DPS in the limit r → 1, is very much similar to anomalous dispersions in the nonlinear mediums with population inversions [23][24] , and gain doublets [5][6][7] . Besides these active nonlinear structures, anomalous dispersion has been also shown as a result of scattering in passive structures that facilitates tunnelling of light [8][9][31][32][33][34] . However, such system is very lossy since the anomalous dispersion is for the evanescent solution of the system. On the other hand, in the case of DPS (with the limit r → 1), the anomalous dispersion is for the real solutions (i.e., propagating waves) of the system.
The refractive indices ( n p and n g ) obtained for the limit r → 1 is indeed useful to describe and infer the behaviour of the DPSs with finite values of R. For a finite R, the DPS behaves like a multimode optical cavity made of an artificial material with the dispersive refractive index ω ( ) n p . For an illustration, in Fig. 3(a), we show the dispersion curve for R = 500. As we can see from this figure, the dispersion curve for R = 500, which is in similar shape as Ω( ) k [ Fig. 2(a)], is discontinuous at each half of the Brillouin zone (BZ), forming discrete bands in the vicinity of ω c . Imporatanltly, these discrete bands are flat, signifying the dispersions of slow or localized optical modes 35 . The dispersion curve for R = 500 in the reduced zone scheme is shown in Fig. 3(b). Note that, the dispersion curve in the reduced zone scheme is better known in the name of photonic band structure 29,36 . In Fig. 3(c-f), we show the photonic band structures for R = 750, 1000, 2000, and 5000. The frequency positions of the flat bands, to a very good approximation is given by where m is a integer that indexes the flat band in the photonic band structure. The indexing scheme is shown in Fig. 3(b). The frequency spacing between the flat bands, ω ∆ , can be expressed using the group index of the DPS in the limit r → 1 [ ω ( ) n g ] as This frequency spacing for the DPS with the finite R is the same as the frequency spacing in a Fabry Perot cavity made with a dispersive dielectric material of refractive index ω ( ) n p , with a length Ra. From Fig. 3(b-f), we can see that the density of the flat bands increases as a function of R. This is because the discretization step (i.e., the length of half ) decreases as R increases. Further, these figures also indicate that the density of flat bands are not uniform across the frequencies of the stopband. For each R in Fig. 2(b-f), the density of the flat bands is maximum near SPS band edge frequencies, and the density is minimum near the SPS stopband centre. These density variations are due the dispersive nature of ω ( ) n p , and can be easily understood from the frequency spacing [Eqn. 3] which is inversely proportional to ω ( ) n g . The values of ω ( ) n g are maximum and minimum for SPS stopband edges and stopband centre, respectively [see the n g plot in Fig. 2(c)]. Therefore, the densities of flat bands are maximum and minimum for frequencies near SPS stopband edges and stopband centre, respectively.

DPS as a Photonic Harmonic Oscillator
In order to perceive the intriguing dispersion of the DPS, let us examine the dielectric function in the vicinity of green dots in Fig. 1. For the sake of discussion let us pick x = a s /2. Using the approximation θ θ ≈ sin for a small angle θ, it can be easily shown that near x = a s /2, and for a large R Eqn. 1 becomes ε ε ε As can be seen from this equation, in the proximity of x = a s /2, the strength of the rapid dielectric modulation-the modulation with the period a (i.e., the amplitude of sin When x is exactly a 2 s , the strength of the rapid dielectric modulation is zero, and we have ε( ) x = ε 0 . This means, at this position the DPS is completely shielded from the effect of the SPS (i.e., the rapid dielectric modulation). So, any light with a frequency in the vicinity of the SPS stopband centre (i.e., ω ω = c ) tends to concentrate at = x a 2 s . A slight deviation from = x a 2 s , causes the light to face the rapid dielectric modulation of the SPS in a linearly increasing strength, ε ε ε , and as a consequence the light will be reflected back towards = x a 2 s . As the stopband resulting from the rapid dielectric modulation, is proportional to the strength of the s , the reflection will be stronger as the deviation from = x a 2 s increases. As we shall prove in the following, this scenario is analogous a lossless harmonic oscillator with linearly increasing restoring force 30 .
In order to analyse the localized modes at = x a 2 s , let us first move the origin of the x-axis from 0 to a 2 s . In the new coordinate system, Eqn. 1 becomes Fig. 4(a), we illustrate ε( ) x , and its slowly varying amplitude functions in the new coordinate system. As we can readily see from this figure, ε gx 1 provides a very good approximation to the slowly varying amplitude over half of the unit cell [i.e., from − a s /4 to where ( ) E x is the electric field. In order to solve Eqn. 4, assume the DPS to be a SPS with a linearly perturbed dielectric function in the slow spatial scale. Consequently, ( ) E x can be expressed as a linear combination of the SPS's modes,  is the center of SPS bandgap. Changing the variable, x, to a dimensionless position variable, α = z x , and eliminaing q and p in Eqns 6 and 7, respectively, we have two independent second order differential equations, . In Fig. 4(b) we plot ( ) S x n using Eqn. 14 for n = 0 to 5, and R = 750. In the same figure, we have also plotted the similar quantity obtained from the exact numerical calculation (based on the plane wave expansion method 29 . Figure 4(c) shows a similar plot to Fig. 4(b), however for R = 2000. As we can see from these figures, both analytical and exact calculations are in very good agreement. The degree of agreement reduces when the slowly varying electric field amplitude moves far from the centre of the DPS's unit cell.
In order to determine the continuous dispersion relation in the limit r → 1, let us write the wavevector in the extended zone scheme as = ± k n G g 2 2 . Using this k, and the frequency expression, ω ω = ± In Fig. 4(d), we compare the dispersion curve obtained from Eqn. 15, with respect to the exact numerical calculation (i.e., the blue curve in Fig. 2). As we can see from this figure, we have a good agreement between these two curves for frequencies in the stopband of the SPS, and in the vicinity of ω c , the agreement is perfect. The discrepancies in the wavevector values of the analytically obtained curve with respect to the numerically evaluated curve  Fig. 4(a); green dot in Fig. 1(a)] is essentially a tunnelling process, and therefore superluminal in nature 8-9,31-34 .

Application and Optical Performances of DPS
DPSs can be used for designing high quality broadband, and multichannel slow light devices. The harmonic modes of the DPS exhibit Gaussian evanescent tails, which decay smoothly. Therefore it naturally generates resonant peaks of high quality factors 39 , despite of a geometrical structure with a low refractive index contrast. This is favourable for many applications such as high temperature realization of Bose-Einstein condensation of exciton polaritons 40 , and realization low threshold nonlinear optical devices, using the abundant low refractive index optical materials.
In the metamaterial cavity section, we showed that for a finite R, the DPS exhibits many flat bands (for example see Fig. 3). In the transmission spectrum, these flat bands will appear as sharp resonant peaks. Figure 5(a) shows the schematic of a single unit cell DPS (N = 1). This single unit cell DPS acts as a metamaterial cavity. The dielectric constant of the ambience is taken as ε ε + 2 0 1 to match with the dielectric constant at the edge of the unit cell [ Fig. 5(a)]. The device with the schematic as shown in Fig. 5(a) can be easily fabricated using the holographic interferometric techniques 26 . Figure 5 Each sharp transmission peak in Fig. 5(b), can be labelled with an integer m using the frequency indexing scheme used in Fig. 3(b). We showed the frequency labelling for the transmission peaks of R = 400 in Fig. 5(b). This labelling will assist us to compare the locations of simulated transmission peaks with the theory developed in this paper. Figure 5(c) compares the positions of the transmission peaks with the resonant frequencies obtained from Eqn. 2 [i.e., from the continuous dispersion curve: blue curve in Fig. 2(a)] and Eqn. 12 [i.e., from the theory of photonic harmonic oscillators]. As the figure suggests, both of these frequency expressions serve as good approximations to the frequencies of the transmission peaks. The most profound transmission peak of the DPS is the peak with m = 0. This peak has the highest quality factor, and the peak is visible even for small values of R. In Fig. 5(d), we plotted the frequency of the m = 0 resonant transmission peak (obtained with the full numerical simulation; see the Methods section) as a function of R. As R increases the frequency of the m = 0 peak move towards to the centre of the SPS bandgap, ω c . The position of the m = 0 transmission peak is in perfect agreement with the prediction of the harmonic oscillator theory which gives . In order to assess the quality of transmission peaks with m = 0, in Fig. 5(e), we plot the logarithmic values of the their quality factors as a function of R, for ε 1 values of . 0 16 and . 0 08. As we can see from this figure, the quality factor increases exponentially as R increases, and the quality factor is comparatively high when the dielectric modulation is of the DPS is high.
For N > 1 the finite DPS is indeed a coupled system of harmonic oscillators [ Fig. 1]. An essential nature of any coupled oscillators is the splitting of the resonant peaks [41][42] . The DPS with N > 1, and ε 1 = 0.16 [the schematic is shown in Fig. 1] is numerically simulated using an ambience dielectric constant of ε 0 [see the Methods section for the details]. Figure 6(a) shows the resulting transmission spectrum near the frequency ω r 0 for N = 2, 3, 4 and 5. When N = 2, there is one transmission peak with frequency of ω r 0 . This peaks corresponds to the localized mode of m = 0, in the vicinity of one green dot in Fig. 1. When N = 3, there are two of such modes couples together, and as a result the original peak at N = 2 splits in two peaks [see Fig. 6(a)]. In general, for the DPS with N unit cells, the coupling of modes in adjacent oscillators results in N-1 closely spaced peaks. The splitting of the peaks is systematically depicted in Fig. 6(b). The frequency span of these closely spaced peaks equal to the bandwidth of the m = 0 flat band of the infinite system ( → ∞ N ), and the span can be obtained from the DPS band structure calculations. For R = 50, and ε 1 = 0.16 the frequency span obtained from the band structure calculation is from ω = 0.309 to 0.3098. If we choose to work around the telecommunication wavelength, then, for a = 480 nm, this normalized frequencies translates to wavelengths from 1549 to 1553 nms. This frequency span is adjustable. If we would like to have a narrow frequency span, then R has to be increased to generate a flatter band.
The closely spaced peaks for any given N in Fig. 6(a) display non-uniform quality factors. The outermost peak [see Fig. 6(a,b)] exhibits the largest quality factor. In Fig. 6(c), we show the quality factor of the outermost peak

Constructing a high dielectric contrast DPS
As we have mentioned in the metamaterial cavity section, the DPS with the dielectric function as in Eqn. 1 can be implemented in many different ways [25][26][27][28] . Most of these methods generates continuous dielectric profile, and therefore they have two important limitations: 1) generating structures with large dielectric modulations (i.e., large ε 1 ); 2) generating structures that are amenable for mass production via lithographical techniques. Therefore, in this section, we would like to introduce the method of creating dielectric profiles as in Eqn. 1, however, with a large dielectric contrast. The method can be implemented either using a multilayer deposition or standard lithographical techniques.
Recall that Eqn. 1 is actually a cosine series. Therefore, in general, we will have a DPS as long as the Fourier series of any periodic dielectric function, at least in an approximation, takes the form of Eqn. 1. Thus, what is really needed to form a DPS, is a dielectric function with two closely spaced frequency peaks (at frequencies G and G/r) in its' spatial Fourier spectra.
Consider a dielectric profile, ε ( ) x A as in Fig. 7(a) which has a period a. This is a binary profile with alternating dielectric constants of ε a and ε b . The fundamental harmonic of ε ( ) x A occurs at the frequency G = 2π /a. In Fig. 7(b), we have a similar dielectric function ε ( ) x B , but with a period ra. The fundamental harmonic ε B (x) is therefore, at the frequency 2π /ra = G/r. Now if we linearly combine ε ( ) [see Fig. 7(c)] then by the linearity of the Fourier transform, the new function ε( ) x will have two Fourier peaks at frequencies G and G/r. Therefore to a good approximation this results in dielectric profile similar to Eqn. 1. Note that in Fig. 7(a,b), the dielectric functions are binary valued, however the dielectric function in Fig. 7(c) is not a binary profile. In order to implement this dielectric profile, we need three materials with dielectric constants ε a , ε b , and ε ε . ( + ) 0 5 a b . Although the linear combination looks simple in its' operation, the real implementation requires a third material with the dielectric constant ε ε . ( + ) 0 5 a b . This condition can be relaxed, if we use a logical combination, instead of the linear combination. The output of a logical combination is always binary, and therefore if we combined ε ( ) x A and ε ( ) x B , using a logical operation at every x, then we will obtain a dielectric profile with the binaries ε a and ε b . Further, if we propely choose the duty cycle of ε ( ) x A and ε ( ) x B , this will also give two strong harmonics at frequencies G and G/r. For an example in Fig. 7(a,b), assume the length of the ε a portion within each period is 0.2a. If we treat ε a and ε b to be equivalent to the binaries 1 and 0, and logically combine, , using a logical OR combination, then the result is the dielectric profile as shown in Fig. 7(d). The Fourier transform of this [ε a and ε b are treated as binaries 1 and 0, respectively] (e) Fourier transform of the periodic dielectric function (period = a s ) in "(d)". Here, we assumed ε = .
3 4 a (silicon) and ε = . 1 45 b (silicon dioxide). As we can clearly see, the Fourier transform exhibits two profound peaks at the frequencies G and G/r. Therefore, using the method of linear and logical combinations, we can generate non-continuous, high dielectric contrast DPS structures. The dielectric profiles in Fig. 7(c,d) can be easily fabricated either via multilayer deposition, or lithographical techniques. One important point to note when handing high dielectric contrast logically or linearly combined structures is that the their Fourier transforms will also consists the higher order harmonics. Nevertheless, the continuous dispersion curve for the limit → r 1 is still obtainnable via exact numerical calculations, and the photonic harmonic oscilator theory will serve as a qualitative model that gives a good physical perspective.

Conclusion
In conclusion, we have presented the unique dispersion properties of a DPS with two closely spaced harmonics. Our discussion confirms that the anomalous dispersion of the DPS in the vicinity of ω 0 , is due to the linear perturbation in the dielectric function of the SPS. As we have shown, this linear perturbation is analogous to a presence of a harmonic oscillator, which pulls the light back towards the equilibrium position [i.e., x = 0 in Fig. 4(a), new coordinate system].
One of the key signatures of the DPS is the large density of flat bands with modes of Gaussian tails. This can be used designing high quality broadband, and multichannel slow light devices. The anomalous dispersion of the DPS also can be engineered to generate new class of passive superluminal, and dispersion controlling devices. Although we have presented the DPS in one dimension for an optical wave, the idea can be easily extended, to generate harmonic oscillators, and therefore lossless effective dispersive metamaterials, in higher dimensions, and other physical wave systems.

Methods
The continuous dispersion curve for the limit → ∞ R , and the photonic band structures for finite values of R are obtained using the plane wave expansion method 29 . The continuous dispersion curve is obtained by increasing R to a huge value, until the results are converged. Specifically, in this paper we used R = 10000 to obtain the dispersion curve [Figs 2(a) and 4(d)]. The transmission spectrums for the DPS are obtained using transfer matrix method (TMM) 43 . In TMM simulations, we slice the continuous dielectric profile [Eqn. 1], into large number of spatial steps with uniform dielectric constants. We verified the TMM results independently, using the finitedifference time domain (FDTD) 44 simulation. For FDTD simulation we used the freely available software, MEEP 45 .