1. Introduction

The physical description and even the existence of evanescent gain, or wave amplification that occurs within a passive medium upon total internal reflection (TIR) from a gainy medium, has been the subject of controversy for over 40 years. It has long been observed that an index-guided propagating mode in a dielectric optical waveguide with lossless core and gainy cladding will experience growth [1]. Consider instead the related scenario involving a wave propagating towards a dielectric interface with incident angle above the critical angle θ_c. Conventional optical theory predicts the reflection coefficient magnitude Γ<1 for TIR regardless of whether the medium on the other side of the interface is gainy or lossy [2]. Thus, conventional theory predicts that propagating modes in a waveguide with gainy cladding and lossless core will experience loss (barring some other gain mechanism such as reflections from the outer edge of the cladding) in apparent conflict with the experimentally observed amplification of such modes. In this paper we show that, contrary to conventional theory, fields incident upon a gainy medium beyond the critical angle have reflection coefficients with magnitude greater than unity.

As shown in Fig. 1, we consider a region composed of two dielectric half spaces, where the dielectric interface is constant in x. We label the media filling these half spaces as medium 1 and medium 2, where medium 1 is to the left of the interface and relates to the lossless core of the optical waveguide, and medium 2 is to the right and relates to the waveguide’s cladding. To extract information purely concerning reflection from the interface between the two media and avoid contributions from back reflections within medium 2, we eliminate the back interface corresponding to the outside of the cladding so that medium 2 extends to infinity in the x direction. Both media are nonmagnetic, with dielectric constants ε _r1 and ε _r2 and conductivities σ ₁ and σ ₂. Consistent with a positive index-step optical waveguide, and to allow for TIR, we require ε _r2<ε _r1. Medium 1 is lossless, so σ ₁=0. Medium 2 is either lossy if σ ₂>0 or gainy if σ ₂<0. The conductivity σ ₂ is taken to be a constant independent of frequency.

 

Fig. 1. System under consideration: dielectric half spaces. Medium 1 corresponds to the lossless core of a positive index-step optical waveguide, and medium 2 corresponds to either the gainy or lossy cladding.

Download Full Size | PDF

Assuming exp(jωt) temporal dependence, the spatial dependence for the incident and reflected fields are given by exp(-j k _i·r) and exp(-j k _r·r), respectively, where k _i=k _x1 â _x+k _y1 â _y=k ₁(cosθ_i â _x+sinθ _i â _y). The reflected-field wave vector differs from that of the incident field only in the sign of the x-component so that |k _i|=|k _r|=k ₁=ω(ε _r1 ε ₀ µ ₀)^1/2. Spatial dependence in the second medium is given by exp((g–j k _t)·r). The orientation of the gain vector g corresponds to the direction of increasing amplitude. Owing to continuity of the tangential fields across the interface, k _y2=k _y1 and g has only an x component.

The s-polarized Fresnel reflection coefficient can now be expressed as

Alternatively, using the notation of [2], one can write the reflection coefficient as

where

Here n ₁ is the index of refraction of medium 1, and κ ₂ is related to the imaginary part of the complex index of refraction of medium 2, n̂ ₂, as follows:

A change in the sign of σ ₂ changes only the sign of κ ₂ in Eq. (5). Because κ ₂ is squared everywhere it appears in Eqs. (3)–(4), only its magnitude affects the value of R so that the Fresnel reflection coefficient as given by Eqs. (2)–(5) is not affected by the sign of σ ₂. This conventional application of optical theory shows evanescent gain not to occur without some new interpretation of the theory.

As a result, early observations of gain in the propagating modes of an optical waveguide with lossless core and gainy cladding incited numerous explanations; a small sampling of these are described here. Romanov and Shakhidzhanov (R&S) posited that gain should be possible upon TIR from a gainy medium and forbidden for incidence below the critical angle [3]. Shortly after these results were published, Lebedev, Volkov and Kogan (LV&K) produced experimental values for the reflection coefficient that conflicted with the predictions of the R&S theory [4]. Several years later, Callary and Carniglia presented an alternate theory which predicted gain upon reflection for incidence below the critical angle, as well as above [5]. Soon after, Cybulski and Carniglia proposed a theory to explain the LV&K result, but did not find agreement [6].

The question of the mechanism of evanescent gain has been the subject of a more recent surge of interest. Citing the lack of consensus in the scientific community and the broad range of predicted values for gain upon reflection given in the literature, Lee, Kim, Noh and Jhe asserted that bulk material properties cannot be used to describe evanescent gain [7]. The Lee et al. work insisted that a full quantum mechanical approach must be taken to accurately describe this phenomenon, and it did not make the connection to classical theory. More recently Fan, Dogariu and Wang reexamined the theory as given by R&S, and used numerical methods to predict the Goos-Hänchen shift that occurs upon reflection from gainy and lossy media [8]. Overall, past research has generally agreed that gain can occur upon reflection from a gainy medium and that the amount of gain depends on the incident angle; however, a vast disparity exists in predictions for the circumstances under which gain will occur, and the amount of gain that will occur upon reflection.

Because u ₂ and ν ₂ are squared in Eqs. (3) and (4), their definition contains an inherent sign ambiguity. The presence of this sign choice produces two solutions to Maxwell’s equations which are mathematically valid. Which solution is physically valid? The sign choice affects the resulting Fresnel reflection coefficient magnitude Γ. Conventionally the choice of signs for u ₂ and ν ₂ correspond to taking k _x2 to be positive and g_x to be negative for σ ₂≠0. This yields phase propagation away from the interface and exponential amplitude decay away from the interface in medium 2.

These sign choices stem from the conventional branch cut in the complex plane of system parameters, which assumes that the reflection coefficient varies continuously with incident angle. Is this assumed branch cut correctly defined? Historically, neither theoretical analysis nor experiment has provided a definitive answer, as proven by the enduring lack of consensus within the scientific community on this point. Recent advances in computational electromagnetics techniques and an increase in computing power permit an alternate analysis using the finite-difference time-domain (FDTD) method [9]. A full-wave time-domain Maxwell’s equations solver, FDTD allows for precise control of system parameters impossible in experiment and makes none of the underlying assumptions that are required in theoretical analyses. FDTD is ideally suited to examine this system and answer the question of sign choice.

In this paper we present a new examination of the existence and properties of evanescent gain, using first-principles FDTD analysis. First we model the interaction of a localized wave (a Gaussian beam) with the gainy half space. Using the FDTD-computed reflection coefficients as a benchmark, we define a new branch cut in the complex plane of system parameters. The choice of this branch cut is made specifically so that the analytical description of the reflection coefficient for the localized wave, which relies on a branch cut definition, achieves agreement with the FDTD results. We believe that for a medium with constant conductivity the outcome of this first-principles investigation is a conclusive description of the circumstances under which evanescent gain will occur and the elimination of ambiguities in the Fresnel equations.

2. FDTD analysis of Gaussian beam reflection

To examine the interaction of a localized wave with this dielectric interface we use a 2D FDTD grid (defined as a TM_z grid with non-zero E_z, H_x, and H_y field components). The computational domain is shown in Fig. 2. For these simulations ε _r1=4 and ε _r2=2. Both media are nonmagnetic, as before. Again, medium 1 is taken to be lossless, and medium 2 has conductivity σ ₂=±500 S/m, where we take the upper sign for a lossy medium and the lower sign for a gainy medium. The incident field is a Gaussian beam that is pulsed in time with a Gaussian envelope and carrier frequency f ₀=3.33×10¹⁴ Hz. In medium 1 there are 20 samples per wavelength at the center frequency, and the Courant number is 0.95/√2 in free space. We use the conventional definition of the critical angle, so that here θ_c=45°.

To model the interaction of the incident wave with a true half space and thereby create reflections from only the single interface, we terminate the grid boundaries with uniaxial perfectly matched layer (PML) absorbing boundary conditions for general (including conductive) media [9]. The PML is designed to absorb outward-propagating waves that reach the grid boundary. The performance of PML in terminating lossy media has been documented extensively in the literature. We validated the performance of the PML in providing essentially a reflection-free termination of the gainy medium by comparing the reflections from the PML with those from a perfect electric conductor (PEC). We found that the PML reflection at normal incidence was 140 dB below that from the PEC. In the case of evanescent decay, we require that the interface be defined far enough from the PML in medium 2 that the full decaying field be represented within the FDTD grid; thus we require at least a few wavelengths of dielectric in medium 2 between the interface and right-hand PML. As the incident radiation approaches and passes through the critical angle, any outgoing radiation in medium 2 is absorbed by the PML and does not reflect back into the main part of the grid to interact with the interface again; this is essential to properly represent the semi-infinite half space.

 

Fig. 2. Snapshots of a pulsed-in-time Gaussian beam launched into the FDTD grid, where red indicates positive values of E_z and blue indicates negative values of E_z. (a) shows the incident pulse, and (b) shows the pulse after reflection from the dielectric interface.

Download Full Size | PDF

2.1. Incident wave source condition

To launch the localized wave in the FDTD grid, we specify E_z fields along a sourcing line (the bottom-most row of grid cells, as shown in Fig. 2) at each FDTD time step. The sourcing line forms an angle θ_i with the desired beam axis as shown in Fig. 3. The source field values are pre-calculated using an analytical description of a Gaussian beam modulated with a Gaussian pulse [10]. By controlling the amplitude and phase of the fields along the sourcing line, we are able to launch a pulse with a chosen direction of propagation into the FDTD grid. Fig. 2(a) shows a snapshot in time of the pulse after it has been launched from the sourcing line. It propagates up and to the right to intersect the interface with incident angle θ_i with respect to the interface normal. As long as we respect the paraxial approximation inherent to the analytical description of the Gaussian beam, and use a reasonable grid resolution in the FDTD simulations, we have great latitude in manipulating the Gaussian beam parameters. Here the beam waist W ₀ is modified through its association with the Rayleigh diffraction length, and the incident angle of the pulse is adjusted to sweep through the critical angle for TIR.

2.2. Reflection coefficient calculation

Figure 2 shows two observation lines in medium 1 along which the E_z and H_x components of the incident and reflected waves are recorded. A discrete Fourier transform is used to obtain the field phasors at the carrier frequency f ₀ from the time-domain fields. The observation lines run through the centers of the Yee grid cells, so that the E_z components are located on the observation lines and the H_x components are positioned a half grid cell above and below the line.

After completion of the time-stepping loop, one set of field phasors (either E_z or H_x) is multiplied by a phase factor so that all field components use a common temporal reference point. This has the effect of removing the stagger in time between the electric and magnetic field that is inherent in the FDTD method. Additionally, the geometric mean is taken of the magnetic fields to either side of the electric fields on the observation line. This serves to spatially collocate the electric and magnetic fields [12]. These unstaggered and collocated fields are used to calculate the y-directed power density at each point along the observation lines. The power density values are then summed to compute the incident and reflected y-directed power [W/m in 2D] that passes across each line. The square root of the ratio of the reflected power to the incident power yields the magnitude of the reflection coefficient.

 

Fig. 3. Gaussian beam time-modulated with a Gaussian pulse. The diagonal line where E_z is pre-calculated corresponds to the bottom-most row of E_z field components in the FDTD grid.

Download Full Size | PDF

2.3. FDTD results for the Gaussian beam

In Sec. 2.1 we defined the incident angle of the Gaussian beam as the angle of the beam axis with respect to the interface normal. In reality the spectral content of any single beam results in incidence at a range of angles, and the chosen single “θ_i” in our calculations is an approximation for the weighted average of the full spectrum of incident angles. (For true incidence with a single θ_i we would use plane wave analysis.) Each angular component of the incident field has a reflection coefficient that depends on the component angle with respect to the interface normal; the sum of all such field components produces the beam reflection coefficient.

The angular spread of a Gaussian beam is inversely related to its spatial width with a plane wave serving as the limit of an infinitely wide beam with a single incident angle. The various curves with data markers in Fig. 4 show the FDTD-computed reflection coefficient magnitude at the carrier frequency as a function of incident angle for several beam widths. Clearly the incident field angular spread strongly affects the resulting beam reflection coefficient near the critical angle. We now have some explanation for the disparity in reported experimental results: gain upon reflection is strongly dependent on the angular spread of the incident beam, especially near the critical angle.

As will be discussed in more detail in the next section, the solid black line in Fig. 4 is the analytical Fresnel reflection coefficient where we have taken k _x2 to be negative for incident angles greater than θ_c. Figure 4 shows that as the width of the beam waist increases, the FDTD beam results approach the analytical plane wave results obtained using this non-conventional sign choice. The mathematical justification for this sign choice is discussed below.

3. Theoretical analysis: a new branch cut in the complex plane

As stated previously, the spatial dependence in medium 2 is given by

 

Fig. 4. Comparison of three FDTD-computed beam reflection coefficient magnitudes for a gainy half space with a critical angle of 45°. The three different beam waists W ₀ are given in terms of λ, the wavelength at the carrier frequency in medium 1. The plane wave reference curve is computed using Fresnel theory, taking k _x2 to be negative for θ_i>θ_c.

Download Full Size | PDF

Importantly, if g_x and k _x2 are known unambiguously, then from Eq. (1) the reflection coefficient is also known unambiguously.

Applying the wave equation to any field component in medium 2 yields the dispersion equation

Solving for (−g_x+jk _x2) yields

where $k_0=\omega \sqrt{{\mu }_0{\epsilon }_0}$ . Let the complex quantity ξ represent the term in parentheses so that

In Eq. (9) the real and imaginary parts of ξ are grouped in parentheses. The real part ξ_r is only a function of the incident angle θ_i for given material parameters, while the imaginary part ξ_i is only a function of the conductivity σ ₂ for a given frequency. We can now write

Since k ₀ is real, to uniquely determine g_x and k _x2 our task is reduced to determining the proper branch cut for ξ ^1/2. (Thus, rather than determining which of the four possible combinations of signs for g_x and k _x2 is correct, we must determine one branch cut.)

Regardless of the conductivity, we define the critical angle θ_c to be

At the critical angle the real part of ξ is zero. We will say that the field is evanescent in medium 2 if θ_i>θ_c(ξ_r>0). Conversely, the field is said to be propagating in medium 2 if θ_i<θ_c(ξ_r<0). Again we note that we adopt the terms evanescent or propagating independent of the conductivity (since only the incident angle affects the real part of ξ).

 

Fig. 5. The (θ_i,σ ₂)-plane. Each point in the plane maps to a unique value of ξ as given by Eq. (9). The point A corresponds to (θ_i=θ_c,σ ₂=0), i.e., ξ=0 where the field is neither propagating nor decaying in the x direction. The incident angle θ_i is between zero and π/2. The conductivity σ ₂ can take on any real value. The four quadrants are distinguished by whether they representing a propagating or evanescent field in either a lossy or a gainy medium. These quadrants also correspond to the different signs for the real and imaginary parts of ξ.

Download Full Size | PDF

At this point it is instructive to visualize a two-dimensional plane in which the independent and orthogonal coordinates are the incident angle and the conductivity in medium 2, i.e., the “(θ_i,σ ₂)-plane,” which is shown in Fig. 5. We consider incident angles θ_i between 0 and π/2 while the conductivity can take on any real value.

The point A in Fig. 5 corresponds to (θ_i=θ_c,σ ₂=0)—a wave incident on a lossless medium at the critical angle. At this point ξ=0 so that the field in medium 2 would be neither propagating nor decaying in the x direction. As will be discussed in a moment, the branch cut (which is really pertinent to the ξ-plane and not the (θ_i,σ ₂)-plane) will be defined as the line segment (θ_i=θ_c,σ ₂<0) which corresponds to the heavy line drawn below the point A.

The ξ-plane is shown on the left side of Fig. 6. The corresponding quadrants from the (θ_i,σ ₂)-plane carry over directly as shown by the quadrant labels. The branch cut corresponds to the negative imaginary axis. Given the branch cut shown in Fig. 6, an arbitrary point in the ξ-plane can be written as

Returning to g_x and k _x2, we have

Thus one can express g_x and k _x2 as

The mapping of the ξ-plane to the (g_x,k _x2)-plane is shown in Fig. 6.

The quadrant of a point is dictated by the value of θ_ξ while the branch cut dictates the valid range of θ_ξ. Given this, the gain g_x and phase constant k _x2 take on the following values:

 

Fig. 6. The mapping of the ξ-plane to the (g_x,k _x2)-plane. Points in the (g_x,k _x2)-plane are obtained from the square root of points in the ξ -plane scaled by $\omega \sqrt{{\mu }_0{\epsilon }_0}$ . In the (g_x,k _x2)- plane points below and to the right of the heavy diagonal line are not realizable. The sector labels shown in the (g_x,k _x2)-plane are taken from the corresponding quadrant labels in the ξ-plane.

Download Full Size | PDF

Referring back to Fig. 5, the dividing line between quadrants I and IV is given by (θ_c<θ_i≤π/2,σ ₂=0). For these values of θ_i and σ ₂ the field in medium 2 is purely evanescent and there is no loss or gain. Along this line k _x2 is zero while g_x is non-zero and it varies continuously to either side. This line essentially corresponds to the branch cut that has been advocated in the past (which yields 0≤θ_ξ<2π). Were one to use this line as the branch cut there would be a discontinuity in g_x as a function of conductivity for a fixed angle.

No matter where the branch cut is placed, as one varies either the conductivity or the incident angle a discontinuity will and must exist in either g_x or k _x2 or in both. For the branch cut we propose, there is a discontinuity in both g_x and k _x2 (a change in sign of both quantities) as one passes through the critical angle for negative conductivity. Such a discontinuity is not unprecedented as the Goos-Hänchen shift is absent for angles less than the critical angle but present for any angle greater than the critical angle [11]. The shift is either there or it is not. It does not vary continuously and the discontinuity occurs at the critical angle. In the same spirit, and consistent with the results obtained from the FDTD simulations, we have kept the critical angle as the branch cut even when gain is present.

Note that region IV is the only region where phase propagation is directed back toward the interface. The gain vector is also directed toward the interface. This produces a reflection coefficient with a magnitude greater than unity (this is only true for region IV). Thus, a slight change in incident angle, i.e., a change that transitions from the propagating to the evanescent side of the branch cut, will result in the magnitude of the reflection coefficient changing from a value less than unity to one greater than unity. This “switch flipping” will happen discontinuously as shown for the plane wave curve in Fig. 4. This discontinuity is consistent with all the FDTD simulations we have performed.

3.1. Theoretical analysis of Gaussian beam reflection

The pulsed incident wave can be decomposed into its frequency-domain components and, for any given frequency, the corresponding Gaussian beam can be represented as a continuous spectrum of plane waves as depicted in Fig. 7. The χζ-coordinate system is such that the central beam axis is aligned with the χ axis thus forming the angle θ_i between the x and χ axes.

 

Fig. 7. Decomposition of an incident Gaussian beam to a spectrum of plane waves. The reflected beam is similar to the incident beam except it is weighted by the reflection coefficient and propagates to the left.

Download Full Size | PDF

The incident field is given by

where the amplitude of the angular spectrum is

and k_ζ and k_χ represent the wave numbers in the ζ and χ directions, respectively. A component of the beam not propagating along the beam axis forms an incident angle θ_b with the x axis. This angle can be greater than or less than the nominal incident angle θ_i.

The reflected beam is obtained from the incident-beam spectrum but a sign change accounts for the fact that the reflected field propagates in the -χ direction and the reflection coefficient weights the angular spectrum. The reflected field is thus given by

Because of the geometry of the FDTD simulation, only wave numbers corresponding to propagating waves with “beam angle” θ_b between 0° and 90° are considered. Although the reflection coefficient changes depending on the branch cut that is used, the form of Eq. (18) is independent of the branch cut. (We note that in practice, unlike the depiction in Fig. 7, our calculations assume the origin of the χζ-coordinate system is collocated with the origin of the xy-coordinate system. Since we are ultimately interested in power and not the fields, the choice of the origin is immaterial.)

From the incident and reflected fields, one can obtain the incident and reflected power. Taking the square root of the ratio of these two yields the beam reflection coefficient.

3.2. Comparison of analytical and FDTD results for the Gaussian beam

The results of the FDTD beam reflection coefficient calculations for both the gainy and lossy cases are shown in Fig. 8. Also shown are the analytically calculated beam reflection coefficients when using the new branch cut. The Fresnel plane wave reflection coefficients calculated assuming either the new branch cut or the conventional branch cut are also plotted. The clear agreement between our FDTD result, which depends purely on Maxwell’s equations, and the analytical result that employs the new branch cut, gives strong support for the use of this branch cut.

 

Fig. 8. Comparison of reflection coefficient magnitudes computed for the incident beam using the analytical and FDTD formulations. (a) shows results for reflection from the lossy half space, and (b) shows results for reflection from the gainy half space.

Download Full Size | PDF

The small difference between the FDTD and analytical data can be attributed to the numerical approximations inherent to the FDTD formulation, as well as the fact that the analytical solution has no curvature in the phase fronts of the Gaussian beam. The phase front curvature affects the incident-angle spectrum which, in turn, affects the value of Γ.

To resolve fine angle-dependent detail in the reflection coefficient calculation we employ the analytic field propagation total-field/scattered-field (AFP-TFSF) technique, a recent advance in FDTD methods which allows for high-accuracy modeling of plane wave interaction with the dielectric interface [13]. Plane wave illumination allows for incoming radiation with a single incident angle, permitting further verification of the precise location of the branch cut in Fig. 5. Further, the analytical field descriptions employed in the AFP-TFSF technique explicitly demonstrate the agreement of Maxwell’s equations with the backwards phase propagation in medium 2 described by the new branch cut. A detailed discussion of the AFP-TFSF technique and plane wave simulations is beyond the scope of this paper.

4. Conclusion

We have conducted rigorous FDTD simulations of the interaction of a pulsed Gaussian beam with a gainy half space both above and below the critical angle. We compared the FDTD-computed reflection coefficient for the beam with the analytically derived reflection coefficient, where the analytical case assumed a new, non-conventional branch cut in the complex plane defined by the critical angle. The excellent agreement that was observed between these two sets of data supports the validity of this branch cut for a medium with constant conductivity.

This new branch cut creates a clear connection between system parameters and sign choice in the Fresnel equations. Once a branch cut is fixed there is no sign ambiguity associated with the square root. Through pulsed beam analysis, FDTD has confirmed the validity of the new branch cut. This research has shown that evanescent gain occurs for incidence above the critical angle and not for incidence below the critical angle. These results agree with the original theory presented by Romanov and Shakhidzhanov as well as the more recent reworking of theory by Fan, Dogariu and Wang.