# Modes in Dielectric Waveguides

PHYSICAL PICTURE

To develop a qualitative sense of how light propagates within a waveguide, we have two main pictures: plane waves and the angular spectrum picture.

(1) Plane waves: If we think of our plane wave with wave fronts perpendicular to the direction of propagation, then we have two pictures: $k$ pointing up and $k$ pointing down. Thus, the sum of these components yields a plane wave that propagates straight down the waveguide, i.e. $k_{tot} = \beta$. The spatial period is equal to wavelength of the mode, which is $2\pi/\beta$.

(2) Angular spectrum: Recalling that the angular spectrum means that we can describe any field as a continuum of plane waves, we consider the picture where the electric field profile is a Gaussian transverse to the direction of propagation i.e. its amplitude is in the x- and y-plane with propagation in the z-direction. We can redraw this picture as a “bunch” of plane waves propagating directly down the fiber that add in the transverse direction to be a Gaussian.

WAVEVECTOR MODEL OF SLAB WAVEGUIDE

We can consider the case first of a symmetric slab waveguide, meaning that the upper and lower guiding layers have the same refractive index i.e. $n_1 = n_3$. This may occur in silica glass waveguides or in semiconductor laser heterostructures. The goal here in this section is to establish that only certain $k$-vectors that satisfy the phase conditions of the waveguide will propagate. The selection of the waves is dependent on the thickness of the layer where the wave is propagating, and on the refractive indices of $n_2$ and $n_1 = n_3$. We can derive the exact conditions by looking at the following figure:

As we might imagine, if the guiding layer thickness is small, then only a few modes will be allowed to propagate, and vice versa.

MODE STRUCTURE

Although the wavevector picture is not entirely accurate, as we would need to use EM theory in order to get the actual field profiles, we can nevertheless press forward with this idea in order to get a qualitative idea about what the field profiles will be. This perhaps will then make the EM theory picture less shocking.

ELECTROMAGNETIC ANALYSIS OF SLAB WAVEGUIDE

Alright, here we are at the precipice of understanding waveguide modes via EM theory! To begin, we have from before that the wave equation is:

\begin{equation} \nabla \times (\nabla \times E) = \nabla(\nabla\cdot{E})-\nabla^2E = -\mu\epsilon\frac{\partial^2E}{\partial{t}^2} = -\frac{\epsilon_r}{c^2}\frac{\partial^2E}{\partial{t}^2} \end{equation}

Before, in the assumption that the medium was homogeneous ($\epsilon$ and $\mu$ constant over all space) with no free charges ($\rho = 0$), Gauss’s Law yielded the result that $\nabla\cdot{E}=0$. This in turn allowed us to derive the usual form of the wave equation. However, in this case, the dielectric constant $\epsilon$ is not constant over all space due to the fact that we have different materials as part of our waveguide. The dielectric constant thus becomes spatially dependent on the x-direction, $\epsilon = \epsilon(x)$. If we then apply Gauss’s Law, we have:

\begin{align} \nabla \cdot D = 0 & = \nabla\cdot(\epsilon{E}) \\ & = \epsilon(\nabla\cdot{E})+(\nabla\epsilon)\cdot{E} \\ & = \epsilon\nabla\cdot{E}+\epsilon_0(\frac{d\epsilon_r}{dx})E_x\end{align}

Rearranging to find $\nabla\cdot{E}$:

\begin{equation} \nabla\cdot{E} = -\frac{1}{\epsilon}\frac{d\epsilon}{dx}E_x = -E_x\frac{d}{dx}ln\epsilon \end{equation}

Substituting this into the wave equation then yields:

\begin{equation} \nabla(-E_x\frac{d}{dx}ln\epsilon)-\nabla^2E = -\frac{n^2}{c^2}\frac{\partial^2{E}}{\partial^2{t}} \end{equation}

In the chosen dimensions of our slab to have propagation purely in the x-z direction, the derivatives of the y-component of the electric field are zero. With this in mind, we can then propose another separable solution to the new wave equation of the form:

\begin{equation} E(x,z,t) = E(x)e^{i(\omega{t}-\beta{z})} \end{equation}

Now then, we have to think back to our definitions of s- and p-polarization, as the form of the final wave equation and solution will be different in each case. We might have expected this, in fact, as we discussed that s- and p-polarization distinctions come into play whenever we have light incident on a surface. Now, for s-polarization, the electric field is perpendicular to the field (“senkrecht”), so in this case the perpendicular polarization direction would be in the $\hat{y}$ direction. Thus, the electric field is $E = E_y\hat{y}$ and the wave equation is given as:

\begin{equation} \frac{\partial^2E_y(x)}{\partial{x}^2} – \beta^2E_y(x) = -\frac{n_j^2\omega^2}{c^2}E_y(x) \end{equation}

Then, for P-polarization, the magnetic field is in the $\hat{y}$ direction, so the electric field is in the $\hat{x}$ and $\hat{z}$ directions. Thus, for the $\hat{x}$ component, we have:

\begin{equation} \frac{\partial^2E_x(x)}{\partial{x}^2}-\beta^2E_x(x)+\frac{d}{dx}(E_x\frac{d}{dx}ln\epsilon) = -\frac{n_j^2\omega^2}{c^2}E_x(x) \end{equation}

For this problem, we can make the weak guiding approximation, which means that the discontinuity in $\epsilon$ can be neglected and thus that the derivative term $\frac{d}{dx}(E_x\frac{d}{dx}ln\epsilon) \rightarrow 0$, since the derivative would be zero over a constant value of $\epsilon$. We can then define the vacuum wavevector $k_0$ as $k_0^2=\omega^2/c^2$ and rewrite the wave equation as:

\begin{equation} \frac{\partial^2{E(x)}}{\partial{x^2}} + (n_j^2k_0^2-\beta^2)E*x( = 0 \end{equation}

where j = 1, 2, 3 for each of the different regions of the waveguide. The form of the solution to this equation depends on whether $\beta^2$ is greater or less than $n_jk_0^2$. We can define a parameter “$h$” as:

\begin{equation} h_j = n_j^2k_0^2 – \beta^2 \end{equation}

and consider the two cases for “$h$”:

Case I: $h^2 > 0$, which corresponds to $h^2 = (n_j^2k_0^2-\beta^2)(n_j^2k_0^2-\beta^2) = n_j^4k_0^4 – 2\beta^2n_j^2k_0^2 + \beta^2 > 0$. This then means that

\begin{equation} n_j^2k_0^2 > \beta^2 \end{equation}

which we can see makes sense with our above answer that $h^2 > 0$ by comparing the sizes of the terms. We then have:

\begin{equation} \frac{\partial^2E}{\partial{x}^2} + h^2E = 0 \end{equation}

which implies an oscillatory form of the solution such as the general form $E = Acos(hx) + Bsin(hx)$.

Case II: $h^2 < 0$, which corresponds to complex $h$ and thus $\beta^2 > n^2k_0^2$ because then $h = \sqrt{n_j^4k_0^4 – 2\beta^2n_j^2k_0^2 + \beta^2}$ would only be imaginary with these relative values of $\beta$ and $n^2k_0^2$. Hence we have:

\begin{equation} \frac{\partial^2E}{\partial{x^2}} – \gamma^2E = 0 \end{equation}

which has solutions of the form E = Ae^{\gamma{x}} + Be^{-\gamma{x}}.

We can now find some qualitative solutions to the wave equation by comparing the sizes of $\beta$ and $n_j^2k_0^2$:

Case (a): $\beta^2 > n_2^2k_0^2 > n_3^2k_0^2 > n_1^2k_0^2$. In this case, from our above qualitative form of the solution for $\beta^2 > n_j^2k_0^2$ of $Ae^{\gamma{x}} + Be^{-\gamma{x}}$, the solution would be exponential in all three layers since $\beta$ is greater than the other $k$-vectors. The specific form of the solution looks like plot (a) in the figure above, which would be obtained from applying the appropriate boundary conditions. Here we can see an unfortunate physical interpretation of this solution: as $|x| \rightarrow \infty$, the field E(x) also goes to infinity, which is of course physically not possible as that would imply infinite energy. This case is thus not possible, even though it solves the wave equation.

Case (b,c): $k_0n_1 < k_0n_3 < \beta < k_0n_2$. In this case, the $\beta$-vector is still greater than the vectors in the two cladding regions, so these have exponential solutions. However, in the core, the solution is sinusoidal. In order to provide physical solutions, the exponential solutions must be decaying. The plotted solution for this case is shown in cases (b) and (c) above, and it can be seen that these solutions form confined, or guided modes.

Case (d): $k_0n_1 < \beta < k_0n_3 < k_0n_2$. This case corresponds to an evanescent wave in the region 1, but oscillating wave solutions in regions 2 and 3. This is often called a substrate mode, since it penetrates into the substrate as it propagates.

Case (e): $\beta < k_0n_1 < k_0n_3 < k_0n_2$. This solution yields oscillatory solutions in all regions, which are called radiation modes. These are unconfined modes.

Now, we can see what happens when $\beta = k_0n_1$.

\begin{equation} sin\theta_i = \frac{\beta}{k_0n_2} = \frac{k_0n_1}{k_0n_2} = \frac{n_1}{n_2} \end{equation}

Thus we can see that the transition between the radiation modes and the guided modes is the same as the transition between partially and totally internally reflected waves. This is the same for each of the solutions with guided modes, as in cases (b) and (c): the wave is being totally internally reflected at each interface.

We are now ready to find the explicit solutions by actually applying the boundary conditions. We will begin with the case of TE waves (s-polarization) in a symmetric waveguide to make our lives easier.

SYMMETRIC SLAB WAVEGUIDE

We first consider TE modes with $n_1 = n_3 < n_2$, such that the propagating wave is guided since the core region has a higher refractive index than the cladding regions. The solutions to the symmetric slab waveguide are either symmetric or anti-symmetric, which are each pictured in the image below:

Now the solution for the top cladding region is given as:

\begin{equation} E = E_1e^{-\gamma{x}},~~\gamma = \sqrt{\beta^2-n_1^2k_0^2} > 0 \end{equation}

which we found because we would have an exponentially decaying field in this region for a guided mode, which corresponds to the conditions given for $\beta$ and $n_1^2k_0^2$. Then, for a guided more, the core region should have an oscillatory solution given as:

\begin{equation} E = E_scos(hx) + E_asin(hx), h = \sqrt{n_2^2k_0^2-\beta^2}>0 \end{equation}

where we can see that the cosine term would produce the symmetric mode and the sine term the antisymmetric term, so the corresponding mode amplitudes are labeled as such. Finally, for the bottom cladding layer, we again have a decaying exponential given as:

\begin{equation} E = E_3e^{\gamma{x}}, \gamma = \sqrt{\beta^2-n_1^2k_0^2} > 0 \end{equation}

where the exponential term is positive because “x” will be negative in the bottom cladding layer, so the decaying part has been taken care of. We can now apply the boundary conditions and match E and $\partial{E}/\partial{x}$ at each interface.

Now, first considering the symmetric modes, we can match the boundary conditions at $x = a$ using the $E_scoshx$ term from the core region solution as follows:

\begin{align} E_s cos(ha) & = E_1e^{-\gamma{a}} \label{E}\\ -hE_ssinha & = -\gamma{E_1}e^{-\gamma{a}} \label{del_E}\end{align}

where we can see that the first equation comes from matching “E” and the second term from matching $\partial{E}/\partial{x}$. The modal profiles are determined by $h$ and $\gamma$ and not the field amplitudes, so we can divide Equation \eqref{del_E} by Equation \eqref{E} to obtain:

\begin{equation} h tan ha = \gamma \label{sym}\end{equation}

We can do the same procedure for the antisymmetric case:

\begin{align} E_asinhx & = E_1e^{-\gamma{a}} \\ hE_asinha & = -\gammaE_1e^{-\gamma{a}} \end{align}

and again dividing:

\begin{equation} -htan(ha) = \gamma \label{anti-sym}\end{equation}

Using Equations \eqref{sym} and \eqref{anti-sym}, we can find the conditions for the symmetric and the antisymmetric modes. However, first we will also scale the mode parameters $h$ and $\gamma$ to the waveguide dimensions by multiplying by the radius of the slab, “$a$”. We then have:

\begin{align} tan(ha) & = \frac{\gamma{a}}{ha},~symmetric~modes \\ & = tan(ha) & = -\frac{ha}{\gamma{a}},~antisymmetric~modes \end{align}

We now recall that $\gamma$ and $h$ are connected via the $\hat{z}$ propagation constant:

\begin{equation} \gamma^2 = \beta^2 – n_1^2k_0^2 \end{equation}

\begin{equation} h^2 = n_2^2k_0^2 – \beta^2 \end{equation}

Then adding these can get rid of the propagation constant $\beta$:

\begin{equation} h^2 + \gamma^2 = n_2^2k_0^2 – n_1^2k_0^2 \end{equation}

Adding back in “a” to get this to a dimensionless form (since “a” has units of length that would cancel with the $\lambda$ in $k_0 = \frac{2\pi}{\lambda}$):

\begin{equation} a^2h^2 + a^2\gamma^2 = (n_2^2-n_1^2)a^2k_0^2 = V^2 \end{equation}

This last relation gives us a very important property: the V-number! Note that the only free parameters, or the only part that matters, is the product $a^2k_0^2 = \frac{4\pi^2{a^2}}{\lambda^2}$, since this scales the wave guide radius to the wavelength of light. Actually, for solid-state cavity design, scaling the radius of the waveguide to the wavelength of light of the junction is a common starting point. The solutions for which modes propagate thus need the symmetric and antisymmetric equations as well as the V-number equation.

The same procedure would follow for TM modes. In this case, we would have $H = H_y\hat{y}$ (note that we do not use the electric field in the $\hat{x}$ and $\hat{z}$ directions in this case, as this would be rather difficult to solve at the boundaries and I am not a monster). The continuity conditions applied here are continuity in the transverse $H$ ie $H_y$ and in the transverse $E$ i.e. $E_z$, although we could also do the normal $D$ component through its relation to the electric field $D = \epsilon{E}$. The result is:

\begin{equation} tan(ha) = (\frac{n_2}{n_1})^2\frac{\gamma{a}}{ha},~symmetric~modes \end{equation}

\begin{equation} tan(ha) = -(\frac{n_1}{n_2})^2\frac{ha}{\gamma{a}},~antisymmetric~modes \end{equation}

TRANSFER MATRIX METHOD

Another approach, and in some cases with much more layers, a better approach is to use the transfer matrices that we derived in the previous section to determine which modes will propagate. This may seem somewhat hard to imagine as being an adept technique, as we will be comparing the layers longitudinally while the light propagates through transversely, as shown in the figure. However, we can see that this will be useful by recalling that for a multilayer system, we can relate the left and right propagating fields as:

\begin{equation} \begin{bmatrix} E_R^+ \\ E_R^- \end{bmatrix} = M\begin{bmatrix} E^+_L \\ E_L^- \end{bmatrix} \end{equation}

We can also recall that we can find the transmission coefficient as:

\begin{equation} t = S_{1,2} = \frac{1}{M_{22}} = \begin{bmatrix} E_j^+ \\ E_i^+ \end{bmatrix}\end{equation}

Now we will consider that if our modes are guided, i.e. there is no transmission to the left or right of the dielectric boundary, then $1/M_{22} = 0$ which means that the transmission coefficient $t \rightarrow \infty$. This divergence of the transmission coefficient even without an input field can be thought of as the presence of the evanescent field outside of the boundary of the waveguide. We already know from the first subsection of this section that when we have TIR, there is an evanescent field along the boundary surface of the interface plane. Thus, since waveguiding works vai TIR, then this case must corresponds to guided modes! For radiation modes, there is transmission on either side of the dielectric stack as the fields penetrate into the surrounding layers, so the transmission coefficient would converge i.e. $M_{22} \neq 0$ and the transmitted amplitude would be $t = S_{1,2}$. Hence, for the transfer matrix approach, we can just find when $M_{22} = 0$.

Since we found before the transfer matrix across a boundary, we can just write the results for what the transfer matrix across the waveguide should be here. For a symmetric slab waveguide with TE modes, we have:

\begin{equation} M = \frac{1}{4}\begin{bmatrix} 1+ \frac{i\gamma}{h} & 1 – \frac{i\gamma}{h} \\ 1- \frac{i\gamma}{h} & 1 + \frac{i\gamma}{h} \end{bmatrix}\begin{bmatrix} e^{-i2ah} & 0 \\ 0 & e^{i2ah} \end{bmatrix}\begin{bmatrix} 1+ \frac{i\gamma}{h} & 1 – \frac{i\gamma}{h} \\ 1- \frac{i\gamma}{h} & 1 + \frac{i\gamma}{h} \end{bmatrix} \end{equation}

which we can see is for TE modes by the form of the outer two matrices and that the waveguide is symmetric from the symmetry of these two matrices as well. Carrying out the matrix multiplication, we have:

\begin{equation} M_{22} = cos(2ah) + \frac{1}{2}(\frac{h}{\gamma} – \frac{\gamma}{h})sin(2ah) \end{equation}

With $M_{22} = 0$ for guided modes, we then have:

\begin{equation} cot(2ah) = -\frac{1}{2}(\frac{\gamma{a}}{ha}-\frac{ha}{\gamma{a}}) \label{origin}\end{equation}

From trigonometry, we can rewrite the cotangent function as:

\begin{equation} cot(2ah) = \frac{1}{2}[\frac{1}{tan(ha)} – tan(ha)] \label{new}\end{equation}

Now, then, comparing Equations \eqref{origin} and \eqref{new}, we can find two solutions:

\begin{equation} tan(ha) = \frac{\gamma{a}}{ha},~symmetric~modes \end{equation}

\begin{equation} tan(ha) = -\frac{ha}{\gamma{a}},~antisymmetric~modes \end{equation}

which are of course the same solutions we had before! Woohoo!

Now, we have two equations that relate $h$ and $\gamma$, so we can use these to the relate the two as:

\begin{equation} hatan(ha) = \sqrt{V^2-a^2h^2} \end{equation}

which doesn’t have an analytical solution, and so the equations must be solved by plotting both the V-number equation and the symmetric/antisymmetric modes with $\gamma{a}$ versus $ha$. This is shown in the figure below: James F. Lotspeich, “Explicit General Eigenvalue Solutions for Dielectric Slab Waveguides,” Appl. Opt. 14, 327-335 (1975)

We can also recall the conditions that:

\begin{align} \beta & > k_0n_2 = \frac{n_2\omega}{c},~no~modes~are~possible \\ \beta & < k_0n_1 = \frac{n_1\omega}{c},~have~radiation~modes \\ \frac{n_1\omega}{c} & < \beta < \frac{n_2\omega}{c},~have~discrete~confined~modes \end{align}

We can then see these dispersion relations plotted below!