Modelling electromagnetically induced transparency media using the finite-difference time-domain method

Time-domain electromagnetic modelling of complex structures which include both non-dispersive media and media exhibiting electromagnetically induced transparency (EIT) require a general, robust calculation method such as finite-difference time-domain (FDTD). We propose a complex-valued, exact two-pole representation of the permittivity of a three-level system which is suitable for integration into the FDTD algorithm via the auxiliary differential equation method. Our calculation model confirmed reported results which were calculated with an approximate representation of the EIT permittivity. Additionally, propagation calculations which mimic slow light experiments were performed. A major advantage of our representation is the ease with which changes in the control field Rabi frequency can be implemented by using a time-dependent permittivity.


Introduction
Many interesting optical phenomena, such as electromagnetically induced transparency (EIT) and ultra-slow light propagation have been observed in highly dispersive media [1,2]. Accordingly, much effort has been invested in developing numerical techniques to model such media. Two such finite-difference time-domain (FDTD) techniques which have risen to prominence in recent years are the piecewise linear recursive convolution (PLRC) method and the auxiliary differential equation (ADE) method [3,4]. Both methods have identical accuracy and memory requirements in the case of Lorentz media [5]. However, the latter makes no assumptions about the linearity of the media in the calculations. Typically, when using the ADE method to model dispersive media, the permittivity is represented by a multi-pole Lorentz media. One can then model EIT behaviour by using a three-pole expansion where a gain line is sandwiched between two absorption lines [6]. However, the gain medium can introduce instabilities into the FDTD algorithm, and the resulting permittivity lacks a simple connection to experimental parameters, e.g., the control field Rabi frequency. With these physical links, such an implementation of the FDTD method would complement EIT experimental studies, while allowing arbitrary dielectric structures to occupy the same computational domain as the EIT media. Such a mixed-media domain is not easily realizable in other calculation methods such as the Maxwell-Bloch algorithm [7]. In this work, we describe the successful implementation of the FDTD method to model electromagnetic behaviour in media exhibiting EIT using a permittivity representation which incorporates experimental parameters.
Modelling EIT in the time-domain poses certain challenges because of the extraordinarily high dispersion which arises from a resonance condition established by a quantum interference effect known as coherent population trapping [8]. Figure 1(a) is a schematic of the atomic system considered in the current study, a three-level scheme, while figures 1(b) and (c) are plots of the real and imaginary parts of the permittivity ε in the vicinity of the ab-transition frequency. (However, the procedure can be applied to any other EIT scheme, such as ladder or V.) The system consists of three atomic energy levels, a, b and c. If a strong field, the control field, is applied to the medium and is tuned to the resonance frequency of the ac-transition, a transparency window opens for a probe pulse with frequency ω when the two-photon resonance condition is fulfilled. In the current study, we assume that the probe pulse is weak in comparison to the control field. When two-photon resonance is satisfied, the single absorption peak in Im{ε} (solid line in figure 1(c)) splits into two peaks (dashed line) located equidistant from the ω ab resonance frequency and separated spectrally by a magnitude equal to the control field Rabi frequency, c = p ac E c /h, where p ac is the dipole moment of the ac-transition and E c is the electric field strength. Essentially, the probability amplitudes of the ab-and ac-transitions destructively interfere, creating a situation where the probe beam does not lose energy to the ab-atomic transition. Therefore, the probe beam is not absorbed as it propagates in the EIT medium.

Formulation
In the FDTD method, Maxwell's curl equations are solved in an time-iterative manner over a spatial grid defined on the computational domain. Partial derivatives are represented by finite differences, resulting in expressions for E and H field equations at their respective locations on the grid that are fully explicit, second order accurate, and dependent only on the value of neighbouring grid points at the current, n-th, and previous, n − 1, time steps. To incorporate dispersive media into this time-domain method, we used the ADE method, in which an additional set of differential equations is used to link a frequency dependent permittivity ε(ω) to the fieldupdate equations. For more comprehensive treatment of the FDTD and ADE methods, the reader is directed to [3].
The ADE method offers the flexibility of modelling dispersion using either a Debye or Lorentz response. At optical frequencies, a multiterm Lorentz dispersive medium is typically used to represent ε of the dispersive media in the domain of the computational method. Thus we propose that an EIT medium can be modelled using the ADE method if a suitable form of ε can be derived from the susceptibility. The susceptibility of a three-level EIT medium can be expressed as [9,10] where = ω − ω ab is the probe beam detuning from the resonance frequency ω ab of the ab-transition, γ ab is the loss due to the ab-transition, γ bc is the decoherence rate, N is the density of excitation centres, and the p ab is the dipole moment of the ab-transition. A partial fraction expansion of (1) follows if we consider the poles 1 Assuming χ 1 and rearrange terms in (2), we can then write the permittivity ε(ω) = ε 0 (1 + χ) compactly in an exact, complex two-pole representation: where ε b is the background permittivity and A l and B l are complex coefficients given by and Figures 1(b) and (c) are representative plots of (3) assuming ε b = 1.0, γ ab /ω ab = 0.006, N|p ab |/(ε 0h ) = 0.01 s −1 and γ bc = 0. For convenience, we normalize frequencies by ω ab since the FDTD calculations are performed with a unit length that is normalized by the wavelength of interest. As shown there is a single absorption peak in Im{ε} when the control field is off (solid line), and there are two absorption peaks and the highly dispersive region between them in Re{ε} when the control field has a value of c /ω ab = 0.143 (dashed line).
The standard FDTD implementation is then performed in combination with the ADE method to derive the field-update equations at the n + 1 time step. The frequency dependence of ε is linked to the electric flux density D via the polarization P l of the l-th pole pair, which in turn, introduces a phasor polarization current J l (t) = ∂ P l /∂t. Using the semi-implicit scheme, we can write an expression for the updated J l J n+1 where t is the time increment. The next step is to solve Ampere's Law for the electric field at the next temporal point, E n+1 , however this requires J n+1/2 l , which up to this point is unknown. Again, we use a semi-implicit scheme along with (6) to derive an expression for J n+1/2 l , which we then substitute into Ampere's Law. Finally, the modified update equation for the E field can be derived. For example, the updated electric field component in the z-direction is given by where σ is the conductivity and Similar expressions can be derived for E x and E y . Since a current term J n+1 l is introduced into the equations, the update algorithm involves one extra step. Thus starting from the current time step components E n z , J n l and H n+1/2 x,y , the E n+1 z components are first calculated using equation (7). Next, the J n+1 l terms are calculated from equation (6), followed by the H n+3/2 x,y components in the usual manner, and the cycle is repeated.
In all, the update algorithm involves only one extra step since a current term was introduced into the equations to link the frequency dependence of ε into the algorithm. Thus the performance of the FDTD calculations does not suffer severely from the inclusion of dispersive media.

Strip waveguide microcavity calculations
As mentioned in the introduction, earlier reports attempted to model the EIT dispersive behaviour using a three-pole Lorentz representation for ε. We performed calculations on the same corrugated strip waveguide microcavity system reported in reference [6] to compare our exact representation of the EIT permittivity with that of the approximation. The geometry of the microcavity system is shown in figure 2. It consists of a silicon strip waveguide with three air holes (radius r = 0.35a) on each side of a defect region. The defect region was formed by increasing the spacing between the 3rd and 4th holes to 1.4a, where a is the lattice spacing, i.e., the nearest-neighbour spacing between holes. The computational domain is terminated on all sides by a perfectly-matched-layer (PML) absorbing boundary condition (BC) [11]. In figure 2(a), the location of the source pulse (TE polarization, e.g., E-field in the plane of the page) is indicated by the black ellipse on the far left of the waveguide, while the detector locations are shown as the vertical lines on either side of the microcavity structure. Comparisons in the transmission spectrum of the microcavity were made between a dispersion-free system, i.e., a homogeneous material in the waveguide, and a system in which only the material in the defect region of the waveguide shows EIT, as indicated by the darker colour. In addition, the same comparisons were made if absorption losses were considered in the materials of the waveguide.
In [6] the values used to create their refractive index profile were not explicitly stated. Thus we adjusted our parameters A l and B l in (3) to produce a good fit with the reported refractive index. More specifically, we focused on matching slope of Re{n} at the resonance of the cavity ω 0 . This is related to the group velocity, v g /c = 0.0453, which is approximately constant over the spectral bandwidth of the cavity. Table 1 summarizes the transmission and Q-factor of the waveguide for four cases: lossless or lossy media in the computation cell combined with or without EIT medium in the cavity. As reported earlier, the Q-factor of the cavity increases with the introduction of a dispersive medium in both the case of the lossless and lossy media. The transmission and Q-factors calculated agree well with the previously reported values for the lossless cases. The discrepancy in the results in the lossy calculations are due to the fact that the detector locations are likely not to match those of the original paper as well as the method the loss was modelled in the silicon and EIT regions.

Slow light calculations
A distinct advantage of our FDTD implementation is the link with experimental parameters in the calculation of the permittivity. In fact, the symmetry of (3) allows the simulation of a dynamically changing c by changing the coefficients A l and B l accordingly in (4) and (5). In light of this, we conducted numerical experiments similar to the slow light experiments reported in 1999 by Hau et al [2] in which the control field Rabi frequency was decreased while the pulse was contained in the EIT medium, offering a even more dramatic slow down. For this, a very long but narrow 2D FDTD computational domain was used which had periodic BC on the sides and PML BC on the ends of the domain. The background of the domain is vacuum, while the EIT medium  The pulse was allowed to propagate completely into the EIT region before c /ω ab was linearly decreased to 0.086 and held at that value. The decrease in c brought the absorption peaks in Im{ε} closer together and increased the slope of the Re{ε}. At the same time the frequency bandwidth of the pulse was also compressed, since during the reduction of c , the spatial width of the pulse was preserved while the pulse is slowed down. Thus the pulse remained inside the EIT transmission window. The steeper dispersion reduced the group velocity of the pulse from 0.59 v g /c to 0.11 v g /c. The slow down is clearly seen in the plots between time-steps 60 000 and 80 000. Finally, c was increased back to the starting value, and the pulse returned to its original speed and propagated out of the EIT region.

Conclusion
In conclusion, we have successfully modelled a three-level EIT medium with the FDTD method through an exact two-pole representation of the highly dispersive permittivity. The ADE method was employed to incorporate the frequency dependence of ε into the time-stepping algorithm.
Calculations on a dispersive microcavity system using our implementation agreed with previously reported results in which the EIT was modelled as a Lorentz medium. Our code included a time-varying ε which was directly related to the control field Rabi frequency during the FDTD time-stepping. This allowed pulse propagation calculations which mimic actual EIT experiments. We presented a slow light calculation where the control field Rabi frequency was decreased by 30% of the original value, resulting in over a 80% reduction in the v g of the pulse. The control field was then returned to the original value and the pulse returned to the original propagation speed. Through our implementation, one brings the versatility of the FDTD method (e.g. arbitrary structures in the computational domain, a variety of BC, and a dynamic temporal control of the permittivity) to the study of EIT phenomena and opens the door to studies involving heterogeneous structures with highly dispersive elements.