Open Access
Add to BibSonomyAbstract
Wave guidance is an important aspect of light trapping in thin film photovoltaics making it important to properly model the effects of loss on the field profiles. This paper derives the full-field solution for electromagnetic wave propagation in a symmetric dielectric slab with finite absorption. The functional form of the eigenvalue equation is identical to the lossless case except the propagation constants take on complex values. Additional loss-guidance and anti-guidance modes appear in the lossy model which do not normally exist in the analogous lossless case. An approximate solution for the longitudinal attenuation coefficient αz is derived from geometric optics and shows excellent agreement with the exact value. Lossy mode propagation is then explored in the context of photovoltaics by modeling a thin film solar cell made of amorphous silicon.
© 2011 Optical Society of America
1. Introduction
As photovoltaic device thicknesses continue to shrink and nanostructured designs emerge, the effects of discretization of their optical mode spectrum are becoming critical to understand. It has thus become necessary to explore the properties of guided mode propagation in lossy materials. Unlike the field of communications, where loss is something to be minimized and is often therefore neglected, the goal of photovoltaic structures is to absorb as much light as possible within specific layers. However, when feature sizes reach the order of 1.0 μm or less, geometric optics no longer provides an accurate description of field propagation at optical wavelengths. It is therefore an important goal to better understand the problem of wave guidance in lossy dielectric films. The solution to this problem serves as a useful model for thin-film photovoltaic devices and reveals many interesting insights into the nature of lossy guidance. It also serves as a useful benchmark from which to evaluate the reliability of low-loss approximations against their exact, analytical solutions.
The theory of guided wave propagation has a long and rich history with many comprehensive references [1–5]. Even so, the problem of guided propagation in lossy media is not nearly as well-understood as the lossless case. Until recently, such problems were normally solved through the use of perturbation theory applied to an ideal lossless model. Such methods are useful for describing guided waves in low-loss structures, but fail to account for the effects of a high-loss material on the exact field solutions. Some efforts have been made to derive more exact solutions to the propagation constants in simple models [6, 7], though the scope of such analysis is currently very limited. Tremendous progress was eventually made in the solution to multilayer guiding structures [8], but the generality of such a solution tends to obscure much of the underlying physics within the model. The problem of wave guidance in amplifying media has been extensively analyzed [9], and has proved very useful in the field of semiconductor lasers. This problem is mathematically analogous to that of wave guidance in a lossy substrate, which is frequently encountered in the field of thin-film photovoltaics. The goal of this paper is to examine the problem of lossy waveguide propagation, and do so by utilizing the familiar formalism of field matching at planar boundaries found in most graduate-level references. Such analysis reveals useful insight into the behavior of such physical models, while also serving as a template for solving variations on the same basic structure. We shall also show that simple approximations to the exact solutions can provide very high accuracy without the need for heavy complication of the classic formulations.
The model structure is depicted in Figure 1 and consists of a dielectric slab with thickness 2h placed between two half-space dielectrics. The dielectric slab, also called the film region, is lossy and may therefore be characterized by a complex index of refraction nf = n + jκ. For simplicity, we shall assume that the cladding layers are both lossless and symmetric, such that they may both be defined by the purely real index of refraction nc. Further complications to the structure may include losses or asymmetries in the cladding layers but will not be considered in this work.
Fig. 1 A dielectric slab with thickness 2h is surrounded by two cladding layers. The slab is a lossy material with complex index of refraction nf = n + jκ. The cladding layers are both lossless and symmetric with real-valued index nc.
Throughout this paper, we shall adhere to the convention of a time-dependent phasor notation with the form of Re{ejkxe−jωt}. This is consistent with the majority of the optics literature but contrary to much of the RF and microwave literature. It also means that a positive value for the extinction coefficient κ implies a lossy material rather than a gain material. When calculating the time derivative of a phasor, we may further make the substitution d/dt = −jω. The ultimate field solutions can likewise assume two possible polarization states, which are the transverse-electric (TE) and transverse-magnetic (TM) polarizations. For the TE case, the electric field intensity E is defined to be polarized along the ŷ-direction. For the TM case, the magnetic flux density H is defined to be ŷ-polarized.
2. Eigenvalue equations
Like all propagation problems in a uniform medium, the functional form of the field solution must satisfy the vector Helmholtz equation. We begin by assuming TE polarization on the electric field, though the derivation follows a nearly parallel procedure for the TM case. The two-dimensional electric field must therefore satisfy [1]
For the case of a lossy film region, the propagation constant k = kf is a complex number that may be written as kf = k0nf. The free-space wavenumber is then given by k0 = 2π/λ0, where λ0 is the free-space wavelength of excitation for the model. The functional expression for E(x,z) can take on many forms, and the only condition for a “correct” choice is the constraint that the solution satisfy the boundary conditions imposed by the geometry of the model.When inside the film region (|x| ≤ h), we shall assume that the electric field can be expressed in the following form:
This expression represents two uniform plane waves propagating through a lossy medium. The constants βx and βz are the transverse and longitudinal phase constants, while αx and αz are the transverse and longitudinal attenuation coefficients, respectively. Both waves possess a forward component along the +ẑ-direction and opposing components along x̂. The ± sign is indicative of a choice between either even (+) or odd (−) symmetry along x. E0 is an arbitrary complex constant that determines the overall intensity and phase of the electric field.The functional form of the field solution can be simplified by defining the complex-valued wavenumbers kx = βx + jαx and kz = βz + jαz. Plugging back into Equation (2) therefore leads to
This expression is identical to the field profiles typically assumed for the lossless case, with the only change being the assumption of complex values for the wavenumbers. The procedure is therefore directly analogous to the lossless case found in most references, which we review in Appendix A. The resultant TE eigenvalue equations are therefore given as Similarly, for the TM case, the eigenvalue equations are Together, Equations (4) through (7) determine the allowed values for the transverse complex wavenumber of the even and odd modes under both polarizations.3. Eigenvalue solutions
Because the eigenvalue equations are transcendental, they can only be solved through iterative methods. We therefore begin by defining the residual function f(kx) as the error in the eigenvalue equation with respect to kx. For example, with the even TE case, this is simply
while for the odd case we have Next, we define the misfit function ϕ as the squared norm of the residual: Note that the asterisk {*} denotes the complex conjugate operation. Our goal is to solve for all complex values of kx such that ϕ(kx) = 0. We shall also enforce the condition βx > 0. This task is accomplished using standard nonlinear optimization techniques [15], and we opted for the steepest descent method with linear line search outlined in Appendix B.As a demonstration case, let us assume TE polarization with nf = 2.0 + j0.5, nc = 1.5, and a normalized film thickness of h/λ0 = 0.5. The logarithmic power of the misfit function, 10 log10 ϕ, is plotted in Figure 2(a) for the even modes. Mode solutions manifest as zeros in ϕ (dark blue regions) with the M = 0 and M = 2 modes indicated. The initial trial solution is chosen by the M = 2 mode for an equivalent lossless system (βxλ0 = 7.27) and indicated by the “X” mark. After applying the steepest decent method, the exact lossy mode solution converges to a value of kxλ0 = 7.89 + j0.77, indicated by the “O” mark on the M = 2 mode. Solving Equations (A.4) and (A.12) (see Appendix A) then produces the full set of propagation constants. Expressed in terms of electrical length, these are
Figure 2(b) shows the electric fields along both x and z generated from Equation (A.7) and normalized to unit amplitude. The horizontal lines indicate the waveguide boundaries. We can also observe the exponential decay in field strength along z, which we naturally expect from propagation through a lossy film.Fig. 2 (a) Logarithmic power of the misfit function, 10 log10 ϕ, for the even modes. Mode solutions manifest as zeros in ϕ (dark blue regions) with the M = 0 and M = 2 modes indicated. Model parameters are nf = 2 + j0.5, nc = 1.5 and h/λ0 = 0.5. The “X” mark indicates the initial trial solution (kx,0λ0 = 7.27). The “O” mark indicates the lossy solution for the M = 2 mode (kxλ0 = 7.89 + j0.77). (b) Normalized electric field profile along both x and z (V/m). Horizontal bars indicate the waveguide boundaries.
Figure 3(a) shows the first four mode profiles for the same model. These mode profiles are generally analogous to the same solutions for a lossless system and propagate in the same manner. An exception to this is the M = 3 mode which actually does not exist as a viable solution for a lossless slab. This is a key difference between lossy and lossless waveguides, where lossy systems possess an equal or greater number of eigenmodes. These extra modes are referred to as “loss guided” modes [9], since they only occur in lossy films. The loss-guided mode also has an effective index of kz/k0 = 1.14 + j0.6, which is less than the real index of either the film or the cladding layers.
Fig. 3 (a) Electric field profiles (normalized to unit amplitude) for the first four modes of the lossy waveguide example. (b) M = 2 profile under increasing values of κ. Vertical bars indicate the waveguide boundaries.
Figure 3(b), demonstrates the change in the M = 2 mode as the extinction coefficient is incremented from κ = 0.01 to κ = 2. It is useful to note that as the film region grows more lossy, its field profile within the slab experiences relatively little perturbation with respect to the lossless case. The only significant change is that higher loss tends to result in greater field confinement within the slab.
Another apparent aspect of the field solution is that γ possesses a real-valued component. Under lossless conditions, γ would be purely imaginary, thereby producing strict exponential decay to the fields in the cladding region. When the film is lossy, the real component to γ causes an additional phase oscillation with respect to x on top of the decay fields. Using the diagram in Figure 4, we can understand this behavior by depicting the evanescent fields as rays that leave the waveguide at some point z and then re-enter at some point z + Δz (a common description for the Goos-Hanchen effect [10]). However, because the film is lossy, the ray R1 carries greater field intensity than the ray R2. This means that, for |x| > h, the time-averaged Poynting’s vector will always contain some nonzero x-component that points toward the film region. Visually speaking, this manifests as a “skewing” of the fields in the cladding region as the constant phase fronts are pushed off-normal with respect to the waveguide. The effect is somewhat subtle for the case in Figure 2(b) but very dramatic for the case in Figure 6(b).
Fig. 4 Ray diagram of the Goos-Hanchen effect for a lossy film. The ray R1 is more intense than the ray R2, leading to a net flow of power into the film at any given point along z. Consequently, the time-averaged Poynting’s vector S in the cladding region has an x-component that points toward the film.
4. Branch cut solutions
A special phenomenon worth noting is the existence of a branch cut in the misfit ϕ due to the square-root function within the residual. The cut occurs when the imaginary component of the radical is set to zero, since this is the point where the square-root of a complex value switches phase during numerical computation. Carrying out this calculation therefore leads us to
The branch point itself is found by setting the argument of the radical to zero, which occurs atAn example branch cut is highlighted by the dashed curves in Figure 5 using the same model parameters as those in Figure 2. Although ϕ maps to two unique values along every point in the complex plane, only one at a time can be rendered on a single 2D graph. The default mapping is shown in Figure 5(a) using the positive value of the radical. An alternative mapping is found by switching the sign of the square-root function in Equations (4)–(7) and then re-deriving ϕ accordingly. This function is shown in Figure 5(b) and reveals a potential set of new zeros in ϕ. However, because of the negative sign on the radical, the imaginary component for γ generally switches sign as well. Such solutions imply a field profile that diverges as |x| → ∞ and are therefore not physically realizable as guided modes.
Fig. 5 (a) Positive and (b) negative solutions to the misfit function with respect to the complex radical in Equation (4). The negative misfit reveals a new set of solutions to the eigenvalue equation, though such solutions are not physically admissible.
Fig. 6 (a) Log-power of the even misfit using nf = 2.0 + j0.1, nc = 2.25, and h/λ0 = 0.5, demonstrating an anti-guidance mode at kxλ0 = 2.8 + j0.77. (b) Normalized field profile of anti-guidance mode. The “skewing” effect on the evanescent fields is much more dramatic in these modes.
The the origins of loss-guided modes can also be understood from Figure 5. As κ grows in value, the branch cut moves further away from the origin to expose zeros in ϕ that otherwise would have been disallowed. Because of this trend, lossy waveguides will always possess an equal or greater number of viable mode solutions than their lossless counterparts. It also allows us to visualize the origins of anti-guidance, which violates the typical restrictions of total internal reflection (TIR) [9]. As an example, Figure 6 shows a waveguide solution using nf = 2.0 + j0.1, nc = 2.25, and h/λ0 = 0.5. Although the cladding index is greater than the real part of the film index, there still exists a mode solution at kxλ0 = 2.8 + j0.77. The anti-guided mode also has an effective index of kz/k0 = 1.95 + j0.07, which is lower than either the film or cladding. The transition for this mode occurs at κ ≈ 0.03, below which the branch cut overlaps the zero in ϕ and removes it as a physically viable solution. Low-loss and loss-less waveguides therefore do not possess guided mode solutions for nc > nf, which serves to reinforce the conventional TIR condition.
5. Longitudinal attenuation and low-loss approximation
Another aspect of lossy mode propagation is that mode attenuation along the z-direction increases with mode number. This behavior is demonstrated in Figure 7(a) using nf = 2.5+ j0.01, nc = 1.5, h = 1.5 μm and λ0 = 1.0 μm. This prediction can be corroborated through numerical simulation with the finite-difference time-domain method [11]. Because the mode profiles change very little under low-loss conditions, we may excite the lossy film with its equivalent lossless mode and then measure the simulated power loss over some fixed distance. This result is plotted in Figure 7(a) by the red markers and agrees extremely well with the analytical computations.
Fig. 7 (a) Longitudinal attenuation coefficient (αz) versus mode number using nf = 2.5 + j0.01, nc = 1.5, h = 1.5 μm and λ0 = 1.0 μm. Exact computations (black) are compared against numerical simulation (red) and the low-loss approximation (blue). (b) Ray propagation in the film according to geometric optics. High-order modes (black) experience greater path-length than low-order modes (red) for a given displacement along z. Some fraction of the overall path-length is also spent within the cladding region due to the evanescent field propagation. Since this region is lossless, it does not contribute to the ray attenuation.
Such behavior in mode attenuation may be visually understood from a geometric optics perspective by considering the ray paths illustrated in Figure 7(b). Because the high-order modes (black) propagate at steeper angles than low-order modes (red), each ray must propagate through more material for the same displacement along z. For the case of low-loss dielectrics (κ << 1), simple trigonometry suggests that
where θ = tan−1(βx/βz) is the elevation angle of the ray and α0 = k0κ is the intrinsic attenuation coefficient of the film. The coefficient Γ is a correction factor that accounts for the fraction of total power in the film region relative to the total power in the mode itself. Because the cladding layers are not lossy, the evanescent fields do not experience any Ohmic power loss and therefore need to be accounted for. A rigorous discussion of this problem is presented in [12] for the case of amplifying media, but naturally carries over to the case of lossy media. The parameter Γ can therefore be set to the modal confinement factor given by where ℰm(x) is the normalized TE field profile along x for the mth mode. These mode profiles can be filled by the solution for an equivalent lossless slab, which is valid as long as the lossless field profiles generally match the lossy case. Using Figure 3 as a guide, we surmise that this will generally be true for κ << 1.Computations from Equation (13) are also plotted in Figure 7(a) and agree very well with the exact values. The utility of applying the low-loss approximation is that it avoids the added tedium of solving Equation (4) for complex wavenumbers. It also presents an intuitively satisfying picture based on purely ray-tracing arguments. For the case of TM modes, the normalized field profile ℰm(x) can be taken from the y-component of the magnetic field for similar results [12]. In general, the net effect is that high-order modes will attenuate more rapidly than low-order modes. The only exception to this trend occurs for modes that propagate very close to the critical angle of the slab. In such cases, the evanescent fields will penetrate far into the cladding region, thereby reducing mode confinement. The result is a small value for Γ that offsets the large secθ dependence.
6. Applications to thin-film photovoltaics
One simple way to model a photovoltaic cell is to consider a lossy dielectric slab backed by a perfect electrical conductor (PEC) as shown in Figure 8(a). Neglecting any surface oxide layers or anti-reflective coatings, such a model is a reasonable representation of a typical thin-film solar cell. We can then mathematically solve for the eigenmodes by extracting only the odd-symmetry solutions from a symmetric dielectric slab with twice the width.
Fig. 8 Configuration for a lossy dielectric waveguide backed by a PEC ground plane. Solutions are equivalent to the odd modes of a symmetric dielectric slab with twice the width.
As an example, Figure 9(a) shows the TE field distribution for the M = 4 mode of an h = 500 nm film of amorphous silicon at λ0 = 600 nm (nf = 4.6 + j0.3, nc = 1) [13]. Figure 9(b) computes αz for all TE modes of the device. The solutions using the low-loss approximation are also indicated by the blue markers at each data point. Clearly, the low-loss approximation produces a very accurate measure of the mode attenuation. As expected, we can also see that high-order modes are generally more lossy than low-order modes, with the M = 6 mode possessing an absorption coefficient that is roughly twice that of the fundamental M = 0 mode. However, the lossy waveguide also possesses an extra five eigenmodes in addition to the seven that exist for the equivalent lossless model. These loss-guided modes cannot be derived from any low-loss approximation and are strictly unique to the lossy waveguide model. It is also apparent that such modes possess dramatically higher values for αz than do the classical modes.
Fig. 9 (a) Electric field profile of the M = 4 mode for thin (h = 500 nm) film of amorphous silicon at λ0 = 600 nm. Indices are given by nf = 4.6 + j0.3 and nc = 1. (b) Longitudinal attenuation coefficient versus mode number. The low-loss approximation is applied to the first seven modes, but does not exist for the extra five in the lossy model.
An example profile of the M = 7 mode is illustrated in Figure 10(a). When viewed along the transverse axis, such modes generally behave in the familiar fashion one should expect from a guided mode. In particular, there are eight characteristic peaks and valleys, as well as a strong evanescent decay in the cladding region. Figure 10(b) shows the same mode when viewed along the longitudinal axis. From this perspective, the mode is almost completely attenuated within a single wavelength cycle. This behavior is characteristic for the additional modes in a lossy dielectric waveguide, which tend to have very high values for αz. Such results may have useful implications for light-trapping applications where the ultimate goal is to maximize light absorption in a finite film [14].
Fig. 10 (a) Electric field profile of the M = 7 mode from the amorphous silicon model. (b) Full 2D profile, showing dramatic longitudinal absorption of the loss-guided mode.
7. Conclusion
This paper derives the full-field solution to guided wave propagation in a symmetric dielectric slab with finite loss. The resultant eigenvalue equations were shown to be identical to the lossless case, but with complex wavenumbers instead of purely real-valued ones. Solutions were computed along the complex plane using an iterative nonlinear inversion scheme. Lossy waveguides also possess a greater number of eigenmode solutions than do their lossless counterparts. These “loss-guided” modes behave much like their conventional counterparts, but with much stronger attenuation in z. It is also possible to demonstrate eigenmodes that violate the traditional TIR requirement nf < nc, producing a form of “anti-guidance” in the waveguide. Using an argument based purely on geometric optics, we have verified a simple low-loss approximation for the longitudinal attenuation coefficient αz. This approximation produces very accurate results when κ << 1 without any need for the complex-valued eigenmode solver. Finally, we then applied this theory to a simplified model of a thin-film photovoltaic cell made from amorphous silicon. This was accomplished by modeling a symmetric dielectric slab with twice the width and solving only for the odd-symmetry modes. It is our intention that this theory will be used to refine the understanding of light-trapping in thin-film photovoltaics, as well as further understanding of any related applications that may be modeled using lossy dielectric waveguides.
Appendix A: Full-field solutions
The exact, full-field solution follows the same general procedure typically outlined in graduate-level textbooks [1–5]. Beginning with the TE case, the assumed functional form of the electric field profile satisfies
The cases for even and odd symmetry are usually managed separately by noting that each case can be specifically expressed using Note that we have taken the liberty of lumping any constant coefficients into the E0 term during the substitution. If we plug either of these expressions back into the Helmholtz equation, we arrive at the dispersion relationThe next step is to assume a functional form for the field profile within the cladding region (|x| ≥ h). In this case, we shall assume that E may be written as
Note that because of the phase-matching condition between planar boundaries, the z-component to the wavenumber is identical in both the slab and cladding regions. The constant coefficient C will be determined later after enforcing continuity at the boundaries. The plus or minus sign is again determined by imposing even or odd symmetry. The complex propagation constant γ = γr + jγi then satisfies the dispersion relation where kc = k0nc is the intrinsic wavenumber of the cladding region.The total electric field through all space may now be expressed as a piecewise function. For the even case, this takes on the form of
with a similar expression for the odd case. We can also solve for H(x,z) by applying Faraday’s law, ∇ × E = jωμ0H. For the case of even symmetry, this gives us with again a similar expression for the case of odd symmetry.The final step is to apply boundary conditions by enforcing continuity along the tangential field components at x = ±h. Starting with the even case, continuity on the E-field implies
Applying continuity on the z-component to the H-field then gives us Now divide Equation (A.10) by Equation (A.9) to find To handle the γ term, we combine Equations (A.4) and Equation (A.6) to arrive at The final result after substitution is Equation (4). Applying the same process to odd symmetry or the TM polarization produces Equations (5)–(7).Appendix B: Steepest descent method with linear line search
Many possible inversion schemes are readily capable of solving the eigenvalue equations given in this paper and a comprehensive reference on various inversion methods can be found in [15]. For this problem we have opted for the steepest descent method with linear line search and is summarized below.
The first step is to compute the complex-valued derivative of the residual with respect to kx. For example, the even TE case produces
with similar expressions for the odd case and TM polarization. Starting with a trial solution kx,0, we then execute the following algorithm: This algorithm is repeated until ϕ(kx,n) falls below some given threshold of error tolerance (say, 10−9). For the case of relatively low-loss dielectrics, a good choice for kx,0 is the mode solution to an equivalent system without any loss. It is also important to bear in mind that the line-search algorithm applies a linear approximation to the derivative of the residual when determining the optimal step size un. As a result, un can often be over-estimated, thus leading to a divergence in the search algorithm. This can be alleviated by ensuring that each step of the algorithm enforces the descent condition, which is written as Whenever a new iteration fails to satisfy this condition, a very simple remedy is to reduce un by a factor of 1/2. The iteration attempt is then repeated until Equation (B.8) finally holds, and the algorithm can then continue normally.References and links
1. C. A. Balanis, Advanced Engineering Electromagnetics (Wiley, New York, NY, 1989).
2. J. A. Kong, Electromagnetic Wave Theory (EMW Publishing, Cambridge, MA, 2000).
3. G. H. Owyang, Foundations of Optical Waveguides (Elsevier, New York, NY, 1981).
4. C. R. Pollock and M. Lipson, Integrated Photonics (Kluwer Academic Publishers, Boston, MA, 2003).
5. A. W. Snyder and J. D. Love, Optical Waveguide Theory (Chapman and Hall, New York, NY, 1983).
6. S. Miyanaga and H. Fujiwara, “Effects of absorption on the propagation constants of guided modes in an asymmetric slab optical waveguide,” Opt. Commun. 64, 31–35 (1987). [CrossRef]
7. K. Fujii, “Dispersion relations of TE-mode in a 3-layer slab waveguide with imaginary-part of refractive index,” Opt. Commun. 171, 245–251 (1999). [CrossRef]
8. T. D. Visser, H. Blok, and D. Lenstra, “Modal analysis of a planar waveguide with gain and losses,” IEEE J. Quantum Electron. 31, 1803–1810 (1995). [CrossRef]
9. A. E. Siegman, “Propagating modes in gain-guided optical fibers,” J. Opt. Soc. Am. A 20, 1617–1628 (2003). [CrossRef]
10. A. W. Snyder and J. D. Love, “Goos-Hanchen shift,” Appl. Opt. 15, 236–238 (1973). [CrossRef]
11. Lumerical Solutions Inc., http://www.lumerical.com/.
12. T. D. Visser, H. Blok, and D. Lenstra, “Confinement factors and gain in optical amplifiers,” IEEE J. Quantum Electron. 33, 1763–1766 (1997). [CrossRef]
13. G. K. M. Thutupalli and S. G. Tomlin, “The optical properties of amorphous and crystalline silicon,” J. Phys. C Solid State 10, 467–477 (1977). [CrossRef]
14. Z. Yu, A. Raman, and S. Fan, “Fundamental limit of nanophotonic light trapping in solar cells,” Proc. Natl. Acad. Sci. U.S.A. 107, 17491–17496 (2010). [CrossRef] [PubMed]
15. J. A. Snyman, Practical Mathematical Optimization: An Introduction to Basic Optimization Theory and Classical and New Gradient-Based Algorithms (Applied Optimization) (Springer, 2005). [PubMed]
