Extended slow-light field enhancement in positive/negative-index heterostructures

Slow light is attracting an increasing attention in the last few years owing to its tremendous applications potential, as for example in all-optical chips krauss1 ; krauss2 ; povinell and optical memories krauss2 . Slow light has furthered the frontiers of current optics, with many works reporting on strong light-matter interactions that manifest themselves as an enhancement of optical gain krauss1 ; baba , Raman signal fan1 or non-linear optical processes krauss1 ; baba , such as third-harmonic generationeggleton . Strong light-matter interactions are associated with high-intensity electromagnetic (EM) fields. Intensity enhancing platforms have been pursuit vigorously in the recent years, which involve a spatial compression of the impinging wave into deep-subwavelength volumes. Typically, spatial light compression is facilitated by a resonant vidal ; schatz ; sfopex ; schuller interaction between light and a structured material -although a non-resonantnonres scheme has also been recently demonstrated.

Slow-light physics is opening an alternate avenue where field enhancement may occur as a result of temporal compression krauss1 ; krauss2 , thus allowing it in principle to extend over mesoscale areas of the order of the free space wavelength. This can facilitate a collective harvesting of light-matter interaction, which is important when a strong signal yield from non-linear optical processes is desired. With the exception of resonant atomic systems khurgin , all other slow-light paradigms krauss2 seem promising in terms of field enhancement, which has been indirectly evidenced by observing a dependent process fan1 ; eggleton ; marin . However, a deeper understanding is still lacking with respect to which are the key controlling attributes pertinent to the dynamics of EM-energy enhancement in slow-light platforms. Also, there seems to be a lack of a clear consensus over the role of spatial compression in slow-light systems krauss2 ; baba .

Different schemes have been employed to control the light’s dispersion and thus speed of information transfer in man-made architectures eit1 ; eit2 ; eit3 ; krauss1 ; krauss2 ; krauss3 ; krauss4 ; baba ; baba2 ; marin ; sridhar ; karalis . An interesting scheme was proposed in 2006 by C. Vandenbem et al. osapro that joins together a positive refractive index medium (PIM) with a negative refractive index medium (NIM)Veselago into a heterostructure waveguide engheta1 , shown to exhibit a flat photonic dispersionosapro . Other NIM-based waveguides had been also reported by Shadrivov et al. andrei and Tsakmakidis et al. kosmas , but with the energy being guided through the negative index material only, evanescently decaying outside. In the system of Ref. kosmas a slow light behavior leading to rainbow trapping via an adiabatic taper was predicted with frequency domain calculations. Given the on-going rapid development of photonic metamaterials soukr ; shalr , these NIM-based waveguides seem highly promising as slow-light platforms. However, the salient dispersion features relevant to an achievable ultra-high slow down factor, as well as the related unconventional propagation characteristics within a wave-optics picture have not been identified. Furthermore, the peculiar dynamics of energy accumulation in these systems, and its relation to spatial and temporal compression is certainly still outstanding. It should also be noted that efficient coupling to slow-light modes is not a trivial matter, and is the subject matter of current on-going research slcoup1 ; slcoup2 .

In this paper, we propose a particular positive-index-medium (PIM)-negative-index-medium (NIM) heterostructure as a composite bi-waveguide paradigm and discuss its extra-ordinary dispersion features in Sec. II. In Sec. III we present time-domain simulations that verify the existence of the trapped slow mode, -predicted by the frequency domain simulations of Sec. II- and also demonstrate an associated high EM enhancement covering mesoscale areas, of the order of the wavelength. We further analyze the propagation characteristics of the slow mode and report for the first time an exotic guided wave having a Poynting vector that is opposite to the direction of energy propagation in one of the participant waveguides in Sec. IV. These dynamic observations further stress on the importance of a large modal index bandwidth that we discuss in Sec. V. Furthermore, in Sec. VI we study the EM energy accumulation for the slow-mode in time-domain. We present an analytical model that underpins the observed dynamic EM energy accumulation, which reveals the participating roles of spatial and temporal compression. This model identifies the prominent features of our prototype’s band dispersion that are responsible both for observing ultra-slow light and achieving an extended high field-enhancement. Finaly we present our conclusions in Sec. VII.

In Fig. 1(a) we show schematics of the bi-waveguide paradigm, with its geometric characteristics designated (we assume an infinite extend in the $x-$ and $z-$ directions). We take a dielectric medium for the PIM ($\varepsilon_{2}=4.0$,   $\mu_{2}=1.0$) and consider a homogeneous left-handed medium (LHM) Veselago for the NIM, similar to the one in Refs. kirby ; Pendry ; rabia . In particular we take:

$$\varepsilon_{1}=\mu_{1}=1-\frac{\omega_{p}^{2}}{\omega^{2}},$$

(1)

with $\omega_{p}=2\pi\cdot 308\cdot 10^{12}\hskip 0.28453ptrad/s$ and $\omega$ representing the frequency of the EM wave. We investigate the $H-$polarized composite guided modes, i.e. light is guided within both the NIM and PIM layers of the bi-waveguide and has evanescent tails outside, with magnetic field along the $z-$ direction. Our design prototype has $\textrm{d}_{1}=1702.5$ nm and $\textrm{d}_{2}=2837.5$ nm.

By assuming guided waves in both NIM, and PIM regions and evanescent tails outside [see also Appendix I] we obtain the dispersion relation for the composite NIM-PIM bi-waveguide guided mode – $\omega(k_{||})$, with $k_{||}$ being the wavevector along the $x$-axis. This is given by:

$$\frac{(z_{1}+z_{2})\textrm{cos}\hskip 0.28436ptOP+(1-z_{1}z_{2})\textrm{sin}\hskip 0.28436ptOP}{(z_{1}-z_{2})\textrm{cos}\hskip 0.28436pt\Delta OP+(1+z_{1}z_{2})\textrm{sin}\hskip 0.28436pt\Delta OP}=\frac{z_{1}-z_{2}}{z_{1}+z_{2}},$$

(2)

where $OP=k_{1y}\textrm{d}_{1}+k_{2y}\textrm{d}_{2},\hskip 5.69054pt$ $\Delta OP=k_{1y}\textrm{d}_{1}-k_{2y}\textrm{d}_{2}$ and $z_{1}=\frac{k_{1y}}{\varepsilon_{1}k_{y}}$,     $z_{2}=\frac{k_{2y}}{\varepsilon_{2}k_{y}}$. $k_{y}$ represents the decaying wavevector outside the bi-waveguide while $k_{1y}$ and $k_{2y}$ represent the wavevectors along the $y$-axis inside the negative-index part and positive-index part of the bi-waveguide respectively. i.e.:

$$k_{y}=\sqrt{k_{||}^{2}-\frac{\omega^{2}}{c^{2}}},$$

(3)

$$k_{1y}=\sqrt{\varepsilon_{1}\mu_{1}\frac{\omega^{2}}{c^{2}}-k_{||}^{2}},$$

(4)

and

$$k_{2y}=\sqrt{\varepsilon_{2}\mu_{2}\frac{\omega^{2}}{c^{2}}-k_{||}^{2}},$$

(5)

with $c$ being the speed of light.

We solve Eq. (2) as a transcedental equation and show the result for the photonic dispersion in Fig. 1(b) with a solid line. Notice, the bi-waveguide mode transitions from a region of a positive band-slope to a region of a negative band-slope with increasing $k_{||}$. Essentially, this implies that the guided mode transitions from being forwards to being backwards tretyakov by passing through a regime with near-zero $d\omega/dk_{||}=\textrm{v}_{g}$ implying near-frozen light [see also dashed line in Fig. 1(c)]. Now, this by itself may not be striking and occurs also in systems where the EM energy is guided through the LHM core onlyandrei . However, what is striking is that the near-zero $\textrm{v}_{g}$ regime extends over a much larger $k_{||}$ interval, that we refer to as modal-index bandwidth from now on modalindex . This extra-ordinary feature was observed before in photonic-crystal waveguides krauss2 ; krauss3 , but in these system does not routinely come with monomodality. Cross-modal talk directs EM energy away from the slow-light channel thus being a major hurdle in the performance of photonic-crystal based platforms. On the other hand, monomodality is an inherent advantage in these metamaterial bi-waveguides engheta1 and flat-bands with a large modal-index bandwidth can be easily tailored with a proper choice of $\textrm{d}_{1}$ and $\textrm{d}_{2}$ osapro . We note that the key importance of these essential characteristics has not been discussed in other NIM-based slow-light platforms kosmas , where the light slow-down factor has not been attested in the time-domain kirby .

We highlight in Fig. 1 the modal-index bandwidth for near-frozen light with the dark shaded region, which transitions to the left and right to a $k_{||}$-region of forwards and backwards faster light respectively. Notice, that the composite guide mode is on the right side of the air-lightline (dashed-line) but on the left side of the lightlines in the NIM (dotted line) and in the PIM (not seen since much further to the right of figure bounds). This signifies the guided-wave nature of the mode within both the PIM and NIM parts of the bi-waveguide heterostructure. We also stress that indeed the group velocity of our system [dashed line in Fig. 1(c)] portrays the velocity of energy propagation, $\textrm{v}_{e}$ Veselago (See Appendix I). [solid line in the same figure].

In the following, we will study the dynamic behavior of light in our proposed paradigm. For this purpose we employ the Finite-Difference-Time-Domain (FDTD) method taflove ; pml – a proven excellent method to capture the time-evolution of EM waves while propagating in space even within negative index materials rabia ; ziolkowski . To excite the composite guided mode in the bi-waveguide heterostructure we employ in FDTD an Attenuated Total Reflection (ATR) set-up voulaback ; voulaatr2 in the Otto configuration otto as depicted in Fig. 2(a) voulaback . The refractive index of the prism, $\textrm{n}_{\textrm{atr}}$, varies as required to yield the desired $k_{||}$ value voulaatr1 . In this manner, we can select from the entire band of Fig. 1(b) a specific band-region centered around modal index $ck_{||}/\omega$ modalindex , with bandwidth that is inversely proportional to the beam waist of the incident Gaussian beam, depicted in the schematics of Fig. 2(a) with its magnetic field polarized out-of the prism plane.

We have confirmed numerically the band dispersion, by considering a pulsed impinging signal and recording the Fourier-transformed spectral response of properly placed detectors above the prism and laterally within the bi-waveguide. We show our results in Fig. 1 as filled circles. We highlight there three particular modes that we designate with open diamonds and name as M1, M2, and M3. These have respectively a positive, near-zero and negative band slope implying a forward, slow-light and backward propagation in each case.

With a quasi-monochromatic excitation voulaprlneg at the relevant frequency in FDTD we confirm the forwards and backwards lateral propagation as we can see in Figs. 2(b) and (d) respectively, where the time-averaged electric-field intensity is shown cropped around the waveguide area. [See also Figs. 3(a) and 3(c) respectively, where the ATR prism excitation area is included]. Now, mode M2, seen in Fig. 2(c) [or Fig. 3(b)], is mostly forwards propagating, but we see some EM energy escaping to the opposite side, due to the unavoidable non-zero $k_{||}$ bandwidth of the impigning Gaussian source. What is striking is the much higher intensity enhancement with respect to the intensity of the input beam that is associated with this mode.

We proceed in analyzing further the curious characteristics of the composite slow-light mode in the heterostructure bi-waveguide. For this purpose, we launch in the FDTD simulation a narrow-bandwidth pulsed signal mono and observe the Poynting vector within each of the joint waveguides. We monitor the $x$-component of the Poynting vector at different lateral locations, at the right-side of the prism for the forwards modes, i.e. modes M1 and M2, and at the left-side of the prism for the backwards mode, i.e. mode M3. Let, $\textrm{D}_{\textrm{det}}$ represent the distance between the lateral line detectors and the relevant side of the prism edge. For each lateral location, $\textrm{D}_{\textrm{d}et}$, one line detector is placed within the NIM part of the bi-waveguide and one line detector is placed within the PIM of the bi-waveguide. Then, for each time step, the $x$-component of the Poynting vector $\textrm{S}_{\textrm{x}}$, is integrated along the respective line detectors in the NIM and PIM parts of the bi-waveguide. We represent this integrated value as $\textrm{S}_{\textrm{x},\textrm{int}}$. In Fig. 4 we show the result for mode M2 and $\textrm{D}_{\textrm{d}et}=2.27$ microns. We show the respective results for modes M1 and M3 in Fig. 5.

We find that the pulse propagates laterally, i.e. along the x-direction, consistent with a wave-optics picture. A geometric optics picture, previously employed to explain slow-light propagation in Ref. kosmas, , clearly fails, as it implies a significant Poynting vector component along the $y-$ direction, which we did not observe here in our simulations. We observe in FDTD a truly astonishing phenomenon, where Poynting vector is antiparallel to the direction of EM energy propagation in one of the participant waveguides. Specifically, the Poynting vector consistently points towards the $-x$ direction in the NIM layer and towards the $+x$ direction in the PIM layer for all cases. Indeed, in all cases the sub-modes within each of the participant waveguides interlock together. They move jointly, as an entity, towards the same direction that is determined by the band slope of the dispersion relation, $\omega(\textrm{k}_{||})$, despite the counter-directional relation of the Poynting vector within each layer of the heterostructure waveguide.

The behavior of the guided waves within each sub-guide as an entity can be further evidenced by observing in the FDTD simulation the arrival time peatross of the pulsed signal at different locations at the right side of the prism (for modes M1 and M2 that are forwards) or the left side of the prism (for mode M3 that is backward). The pulse arrival time within the NIM and PIM parts of the bi-waveguide, at the different lateral locations, will be calculated from the detected $\textrm{S}_{\textrm{x},\textrm{int}}$ with time. Normally, in experiments or simulations the arrival time of a pulse is determined by the peak of the pulse. However, this is a rough method and can entail significant errors for systems where significant pulse broadening occurs during propagation. This is particularly true for slow light systems. J. Petross et al. peatross , have proposed an accurate measure for the pulse arrival time for non-standard systems that we employ here. Thus, we calculate the pulse arrival time, $t_{\textrm{arr}}$, at a certain lateral line detector, $\textrm{D}_{\textrm{det}}$ from the following relation:

$$t_{\textrm{arr}}=\frac{\int\limits_{0}^{t_{\textrm{sim}}}t\textrm{S}_{\textrm{x},\textrm{int}}(\textrm{D}_{\textrm{det}},t)dt}{\int\limits_{0}^{t_{\textrm{sim}}}\textrm{S}_{\textrm{x},\textrm{int}}(\textrm{D}_{\textrm{det}},t)dt},$$

(6)

with $t$ being the time-instant where $S_{\textrm{x},\textrm{int}}$ is monitored at the line detectors in the NIM and PIM part of the bi-waveguide at the $\textrm{D}_{\textrm{det}}$ lateral location. Also, $t_{\textrm{sim}}$ is the total simulation time. Note, since we have different line detectors placed in the participant waveguides of the heterostucture, arrival time is recorded separately in the NIM and PIM parts of the bi-waveguide.

We plot the results for the arrival time at the different lateral detector positions $\textrm{D}_{\textrm{det}}$, in Fig. 3 in panels (a), (b) and (c) for the modes M1, M2 and M3 respectively. The filled circles represent the arrival times calculated in the NIM part of the bi-waveguide. Conversely, the open circles represent the arrival times calculated in the PIM part of the bi-waveguide. The excellent agreement between the two ascertains further our previous observation that the composite mode is an entity, and EM energy propagates with the same speed within the entire waveguide extend, both within the NIM and PIM parts.

Now, we can take a linear fit in each case and from this we obtain the energy velocity of the composite mode. We express this in terms of the speed of light $c$ in each case and show it in the respective figures. Indeed, the observations for the energy velocity within the NIM and PIM part of the bi-waveguide are in excellent agreement, further attesting that the EM energy of the composite mode propagates with a unique speed in all parts. We have expressed the results in fractional form, so the slow down factor is immediately evident. We find a slowdown factor of $\sim 16$, for mode M1, $\sim 303$ for mode M2 and $\sim 9$ for mode M3.

To compare with the theoretical predictions from the dispersion relation we enlist from Fig. 1 the slow-down factor, that correspond to the central $\textrm{k}_{||}$ of the impinging Gaussian beam; these are $\sim 13$, for mode M1, $\sim 105$ for mode M2 and $\sim 9$ for mode M3. We observe an excellent agreement for modes M1 and M3 but we find that mode M2 is slower than expected. This is because the finite spatial extend of the impinging beam, implies a $\textrm{k}_{||}$ span within a band $\Delta\textrm{k}_{||}$. Thus in practical situation one observes the average speed from the mode contributions within this $\Delta\textrm{k}_{||}$ band. For mode M2, this ended up yielding even a higher slow-down factor than the one predicted for the central $\textrm{k}_{||}$ value.

Our analysis in Sec. IV suggests that the overall speed of the EM energy propagation is affected by the finite spatial extent of the impinging Gaussian beam implying contributions from within a $\Delta\textrm{k}_{||}$ band. This further highlights the particular advantage of having a wide modal index span over which the composite bi-waveguide mode is near-frozen, as we have seen in Fig. 1. For example, let’s consider a comparative design with $\textrm{d}_{1}$=2270 nm and $\textrm{d}_{2}$=3405 nm. We show the dispersion relation $\omega(k_{||})$ and corresponding energy velocity (in units of speed of light $c$) for the latter design in Figs. 7(a) and 7(b) respectively. In order to be able to easily compare the result with the system of Fig. 1 we also show the respective dispersion and energy velocity as dashed lines. Evidently the band for the comparative design is not as flat. A near-frozen mode exists only for a very narrow $k_{||}$ band.

We further investigate in FDTD, the implications of the lack of existence of a wide modal-index bandwidth. We set the ATR prism properties to excite in the comparative design the near-frozen mode of $ck_{||}/\omega=1.12$, as predicted from the photonic dispersion of Fig. 7. Yet what we observe (see Fig. 8) is an amphoteric propagation. This emanates from the strong contributions from the adjascent backwards and forwards modes, as the near-frozen mode has a narrow modal index bandwidth in the comparative waveguide. The latter cannot be avoided due to the finite spatial extend of any realistic excitation beam. These results suggest that near-frozen light modal-index bandwidth is the key feature of merit, when considering the design characteristics of the bi-waveguide.

Now that we have verified a large slow-down factor for the M2 mode, we try to understand its relation to the observed intensity enhancement seen in Fig. 2(c) [or 3(b)], which spans over the entire waveguide width, $d\sim 3\lambda_{0}$, with $\lambda_{0}$ being the free space wavelength. Since the waveguide extent in the $z$-direction is infinite, implying that in practice the waveguide can be arbitrary wide in this direction, our results essentially suggest an extended intensity enhancement covering a large mesoscale area, larger than $\lambda_{0}^{2}$. This to our knowledge constitutes a first report of this capability, as typically intensity enhancement is restricted within deep-subwavelength regions vidal ; schatz ; sfopex ; schuller ; nonres .

To understand this further, we excite the M2 slow mode with a quasi-monochromatic wave and monitor the time-evolution of the EM energy density, spatially averaged over a central region of the bi-waveguide lying directly below the ATR prism. We take the time-average of this quantity within each wave-period T, $<\textrm{U}>$, and normalize it with the time-averaged EM energy density of the source $\textrm{U}_{0}$. The result at each period T, shown in Fig. 9(b), implies an energy accumulation with an almost exponential time-response reaching a two-hundred-fold EM energy density enhancement at steady-state.

Let’s consider a simple crude model, depicted in Fig. 9(a), to explain the energy accumulation dynamics. EM energy, $\mathcal{E}$, is fed into the bi-waveguide from above via the ATR prism, at a faster rate than its relaxation sideways in the $+x$ direction [M2 is a forwards mode] counter , because of the ultra-low energy velocity of the M2 mode. In the following with brackets $<>$ we denote the spatial averaged quantities of the energy density, in an area $\textrm{L}_{\textrm{w}}\textrm{d}$ in the central part of the bi-waveguide below the ATR prism, where d represents the total width of the bi-waveguide [see schematics of Fig. 9(a)]. With $\textrm{L}_{\textrm{w}}$ we will denote the beam waist of the evanescent wave illumination emanating from the ATR prism, which is close to the Gaussian impinging wave’s beam waist. Moreover, time-averaged quantities, within a wave period, T are assumed for the EM energy densities, U. Note, we have translational symmetry in the $z$ direction, where the guide is assumed infinite. We focus on a part of the guide with a width $L_{\textrm{per}}$ along $z$. We will see that this chosen part can be arbitrary and does not influence the result. Then, within $\Delta T$, the EM energy that gets fed, into the waveguide would be:

$$\Delta\mathcal{E}_{\textrm{feed}}=F_{c}I_{0}\Delta T\textrm{L}_{\textrm{w}}\textrm{L}_{\textrm{per}},$$

(7)

where $F_{c}$ represents the coupling efficiency and $I_{0}$, the intensity of the impinging Gaussian Beam, thus:

$$\Delta\mathcal{E}_{\textrm{feed}}=F_{c}\hskip 0.85358ptc\hskip 0.85358pt{U_{0}}\Delta T\textrm{L}_{\textrm{w}}\textrm{L}_{\textrm{per}},$$

(8)

where $U_{0}$ is the energy density of the impinging beam. Then the energy that get’s relaxed sideways within $\Delta T$ would be:

$$\Delta\mathcal{E}_{\textrm{rel}}=<S_{x}>d\hskip 2.84526pt\textrm{L}_{\textrm{per}}\Delta T=\textrm{v}_{e}<U>d\hskip 2.84526pt\textrm{L}_{\textrm{per}}\Delta T$$

(9)

So, the increase in EM energy in $\Delta T$ would be, $\Delta\mathcal{E}_{\textrm{feed}}-\Delta\mathcal{E}_{\textrm{rel}}$.
This equals to $<\Delta U>d\hskip 2.84526pt\textrm{L}_{\textrm{w}}\textrm{L}_{\textrm{per}}$ yielding together with Eq. (8) and (9) that:

$$\frac{\Delta<U>}{\Delta T}=-\frac{\textrm{v}_{\textrm{e}}}{\textrm{L}_{\textrm{w}}}(<U>-U_{M})$$

(10)

where

$$U_{m}=F_{c}\hskip 2.84544ptU_{0}\frac{c}{\textrm{v}_{\textrm{e}}}\hskip 2.84544pt\frac{\textrm{L}_{\textrm{w}}}{\textrm{d}}.$$

(11)

Eq. (10), with Eq. (11) yields:

$$<U>=U_{0}\frac{c}{\textrm{v}_{\textrm{e}}}\hskip 2.84544pt\frac{\textrm{L}_{\textrm{w}}}{\textrm{d}}(1-e^{-\frac{\textrm{v}_{\textrm{e}}}{\textrm{L}_{\textrm{w}}}}t),$$

(12)

where time, t, should be an integer multiple of period T.

Thus, the estimation of the feeding and relaxation rates leads to a time-response for the EM energy accumulation that is reminiscent of a charging capacitor, with the characteristic time $\tau$ being equal to ${\textrm{L}_{\textrm{w}}}/\textrm{v}_{\textrm{e}}$. Eq. (12) implies that spatial compression contributes a factor of ${\textrm{L}_{\textrm{w}}}/{d}$ to the overall EM energy enhancement, which is actually quite modest (about four) for the bi-waveguide paradigm. Temporal compression contributes a large factor for slow light modes equal to ${c}/{\textrm{v}_{\textrm{e}}}$. The results of Fig. 9(b) also imply a high coupling efficiency $F_{c}$, more than 60$\%$ with our proposed simple ATR-based coupling scheme. We observe that the overall EM energy enhancement factor is not evenly spread within each sub-waveguide and is also larger than the corresponding electric field intensity enhancement (enhancement of $EE^{*}$). Specifically, we find the latter to be about a hundred in the NIM layer and about thirty in the PIM layer. Note also, that the characteristic time of the exponential EM energy growth, indirectly provides the light slow down factor which we found to be with an exponential fit about 90 for the case of Fig. 9(b)slowdiff .

We note that the comparative design yielded a much lower EM energy enhancement, as can be seen in Fig. 10. This is because additional energy relaxes from the faster, adjascent to the near-frozen mode, backward and forward modes, as any excitation beam inadvertendly includes contributions within a $\Delta k_{||}$ band. Eq. (12) simultaneously stresses the need for a wide modal-index bandwidth of the near-frozen mode but also for monomodality. Actually, it can be shown that the co-existence of any faster channel significantly hampers the maximum attainable EM energy enhancement.

For example, suppose that light couples inside the composite waveguide with a total coupling efficiency $F_{c}$, but to two different modes that co-exist with respected weights $F_{c1}/F_{c}$ and $F_{c2}/F_{c}$, with $F_{c1}+F_{c2}=F_{c}$ and energy velocities ${v_{e1}}$, ${v_{e2}}$. Then the balance of energy fed and energy released sideways will lead to:

$$<U>=U_{0}\frac{c}{\textrm{v}_{\textrm{e}}^{\textrm{eff}}}\hskip 2.84544pt\frac{\textrm{L}_{\textrm{w}}}{\textrm{d}}(1-e^{-\frac{\textrm{v}_{\textrm{e}}^{\textrm{eff}}}{\textrm{L}_{\textrm{w}}}}t),$$

(13)

with

$${\textrm{v}_{\textrm{e}}^{\textrm{eff}}}=\frac{{\textrm{v}_{\textrm{e}1}}F_{c1}+{\textrm{v}_{\textrm{e}2}}F_{c2}}{F_{c}}$$

(14)

Eqs. (13) and (14) clearly imply that the maximum EM energy enhancement that can be attained is adversely affected by the co-existence of any faster channel. They also indicate that the respective saturation time to the maximum EM energy is smaller. This is consistent with what we observed in Fig. 10 for the energy accumulation for the comparative design of Fig. 7, where part of the impigning light couples to faster channels that are adjascent to the near-frozen mode.

In conclusion, we have presented a paradigm NIM-PIM heterostructure supporting an exotic composite guide mode, having an ultra-low energy velocity across a very broad modal-index bandwidth. We have observed in FDTD a most extra-ordinary propagation, where the sub-modes in each layer of the heterostructure interlock together and move jointly in the same direction despite the Poynting vector showing in the opposite direction in one of the layers. We verified numerically an efficient coupling to the slow mode, with speed about $c/300$, leading to exponential growth in the EM energy accumulation, reminiscent of a charging capacitor. Our findings suggest the possibility to achieve an extended electric field enhancement of the order of a 100, facilitated by the near-frozen waveguide mode of ultra-wide modal index bandwidth. We therefore believe, this study will inspire new designs for slow-light platforms for the collective harvesting of strong light-matter interactions.

The magnetic field distribution for the composite guided mode is given by:
$H_{z}$=$\begin{cases}\textrm{A}e^{k_{y}(y+\textrm{d}_{1})}\hskip 0.28453pte^{i(k_{||}x-\omega t)}\hskip 14.22636pt\textrm{for}\hskip 14.22636pty<-\textrm{d}_{1}\\ (\textrm{B}\sin{k_{1y}y}+\textrm{C}\hskip 1.42262pt\cos{k_{1y}y})\hskip 0.28453pte^{ik_{||}x}e^{-i\omega t}\hskip 14.22636pt\textrm{for}\hskip 14.22636pt-d_{1}\leq y\leq\textrm{0}\\ (\textrm{D}\sin{k_{2y}y}+\textrm{E}\cos{k_{2y}y})\hskip 0.56905pte^{ik_{||}x}e^{-i\omega t}\hskip 14.22636pt\textrm{for}\hskip 14.22636pt0\leq y\leq\textrm{d}2\\ \textrm{F}e^{-k_{y}(y-\textrm{d}_{2})}\hskip 1.42262pte^{ik_{||}x}e^{-i\omega t}\hskip 14.22636pt\textrm{for}\hskip 14.22636pty>\textrm{d}_{2}\\ \end{cases}$

where $y=0$ is taken at the PIM-NIM interface. The coefficients in the above $H_{z}$ field distribution are easily determined from the EM boundary conditions at the three interfaces (at $y=-\textrm{d}_{1}$, $y=0$ and $y=\textrm{d}_{2}$). Thus both the magnetic-field distribution and the corresponding electric-field distribution can be easily calculated for the composite guided mode. Then its energy velocity, $\textrm{v}_{\textrm{e}}$, along the $x$-axis [see schematics of Fig. 1(a)] can be given given by:

$$\textrm{v}_{\textrm{e}}=\frac{\bar{S}_{x}}{\bar{U}}=\frac{\frac{1}{8\pi}\int\limits_{-\infty}^{\infty}E_{y}H_{z}^{*}\,dy}{\frac{1}{16\pi}\int\limits_{-\infty}^{\infty}(\frac{\partial(\varepsilon\omega)}{\partial\omega}(E_{x}E_{x}^{*}+E_{y}E_{y}^{*})+\frac{\partial(\mu\omega)}{\partial\omega}H_{z}H_{z}^{*})\,dy},\\$$

(15)

where the expressions in the numerator and the denominator represent the time-averaged Poynting vector and the time-averaged energy-density in C.G.S. units, respectively, integrated along the finite dimension of the composite guide ($y$-axis) and outside to include the contributions from the evanescent tails in vacuum. Hence, the limits span from $-\infty$ to $\infty$ The general expressions for dispersive media Veselago is taken for the energy density, with $\varepsilon$ being either $\varepsilon_{1}$ or $\varepsilon_{2}$ and $\mu$ being either $\mu_{1}$ or $\mu_{2}$, depending which region of the bi-waveguide $y$ lies within, or $\varepsilon=\mu=1$ for the vacuum region. Evidently, for the nondispersive material regions $\frac{\partial(\varepsilon\omega)}{\partial\omega}=\varepsilon$ and$\frac{\partial(\mu\omega)}{\partial\omega}=\mu$.