Guided modes in a slab of crystalline silicon

In this section, a method of obtaining propagation properties of light trapped in waveguides is explained and a modesolver algorithm developed.

Electromagnetic introduction and total internal reflection

The propagation constants of electromagnetic radiation confined by total internal reflection to a silicon slab of finite thickness can be found by imposing boundary conditions on the plane-wave solution of the Helmholtz equation. A plane wave's electric field in space and time can be expresed as

\begin{equation} \vec{E}(\vec{x},t)=\vec{E}_0\exp\left[-i\left(k\vec{x}-\omega t\right)\right] \end{equation}

,

where $\vec{x}=\left(x,y,z\right)^T$, $k=\frac{2\pi n}{\lambda_0}$ with free-space wavelength $\lambda_0$ and refractive index $n$, and $\omega=\frac{2\pi c}{\lambda_0}$ with speed of light $c$, and $\vec{E_0}$ describing the amplitude of the field.

The waveguide

If a two-dimensional silicon slab, henceforth also referred to as waveguide, be infinitely extended along the $z$-axis and extends to $\pm \frac{h}{2}$ in the $y$-direction and the remaining area where $|y|>\frac{h}{2}$ be occupied by air, total internal reflection gives rise to electro-magnetic waves propagating along the $z-$direction that are confined to the silicon layer and only have exponentially decaying amplitudes into the surrounding air. A waveguide of this kind is depicted below.

title

Fig. 1: Waveguide and geometric consideration of plane-wave propagation profiles. Here, $n_1$ corresponds to the refractive index of silicon, while $n_2=n_3=1$ for air.

The smallest angle $\theta$ for which total internal reflexion occurs for a wave of given frequency $\omega$ (or equivalently free-space wavelength $\lambda_0$), is given by Snell's law of refraction;

$n_\mathrm{2}\sin(\theta_\mathrm{2})=n_\mathrm{1}\sin(\theta_\mathrm{1})$. [Add here the angle $\theta_2$ to the figure above]

$\rightarrow \theta_1=\arcsin\left(\frac{n_2}{n_1}\right)\equiv\theta_c$.

The propagation constant

Consider a dielectric medium, infinite in both $x$ and $z$, and of height $h$ in the direction of $y$, surrounded by air for all $|y|>h/2$ as depicted in fig. 1 (a). For a wave propagating in the $z$-direction inside the waveguide with its normal field vectors oriented at an angle $\theta$ towards the interface. The portion of the field in the $y$-direction is bound by total internal reflection, the boundary condition for the $y$-component of the field, defined in fig.\,1(b) as $k_y=k_0n_1\cos(\theta)$, is deduced as follows:

Starting from the lower interface, the wave travels a distance of $h$ to the upper interface and experiences a phase change of $\phi_t=k_yh=k_0n_1\cos(\theta)h$ during propagation. The total internal reflection at the interface induces a phase change $\phi_i$, given by the Fresnel coefficients. After reflection at the upper boundary, another distance $h$ is covered towards the lower interface and upon the second reflection and corresponding phase change $\phi_i$ there, one period is completed and the starting condition recovered. If the total phase change corresponding to a period is an integer multiple $m$ of $2\pi$, the boundary condition \begin{equation} \label{eq:boundarycond} 2\pi\cdot m=2\phi_t+2\phi_i \end{equation} is satisfied and a standing wave confined to the waveguide is observed.\ Such waves are termed \textit{modes} with mode number $m$, where the projection of wavevector $k$ inside the waveguide onto the interface $k_z=k_0n_1\sin\theta_1$ is known as the propagation constant\,$\beta$. Also referred to as the mode's wavevector\,$k_m$, $k_m=\beta=k_z=k_0n_1\sin(\theta_1)$ is an invariant for a wave across all layers in a multi-layer structure\,\cite{Ghatak1987}. A \textit{mode} is also said to have an effective refractive index $n_m=n_1\sin(\theta_1)$, such that \begin{equation} k_m=\beta=k_0n_m=k_0n_1\sin(\theta_1). \end{equation} The modes of a waveguide are further classified by their polarization. Since the phase shift $\phi_i$ experienced by a wave upon reflection at an interface depends on the orientation of the electric and magnetic field vectors relative to that interface, solutions to the boundary condition are slightly different if the electric field is oriented parallel or perpendicular to the plane of incidence.

Complex index of refraction and Beer-Lambert's law

In absorptive media, the refractive index is in fact a complex number $\tilde{n}=n+i\kappa$, of which the real part corresponds to the normal refractive index $n$ and the imaginary part $\kappa$ describes the amplitude attenuation according to Beer-Lambert's law:

\begin{equation} \frac{dE}{dz}=\sqrt{\frac{I}{I_0}}=\exp\left[-k_0\kappa z\right]=\exp\left[-k''z\right], \end{equation}

where the complex wave vector $\tilde{k}=\tilde{n}k_0=(n+i\kappa)k_0=k'+ik''$ was introduced.

Let us reconsider the propagation constant $k_m$ for the complex case.

$\tilde{k}_m=k_z'+ik_z''=\beta +i\frac{\alpha}{2}=k_0n_1\sin\theta_1+ik_z''$.

From the sketch in fig.$\,1b$, it is evident that the distance $d$ travelled by the wave per unit length in $z$ is given by $d(z)=z/\sin(\theta_1) \,\,\,\forall\, \theta_1\,\, \in\, (0,\frac{\pi}{2})$, such that the effective imaginary part of the mode's refractive index is given by

$K=\kappa/\sin\theta_1$,

and consequentially

$k_z''=k_0K=k_0\frac{\kappa}{\sin\theta_1}=k_0\frac{n}{N}\kappa$.

This definition will be referred to as "big K", while $K=\kappa$ is called "small k".

The algorithm for $K_m$ values

From the preceeding it becomes clear that a mode can be uniquely characterized by two parameters, e.g. the free space wavelength $\lambda_0$ and the propagation constant $\beta$. Furthermore, since $\beta=k_0(\lambda_0)\cdot n(\lambda_0)\cdot sin(\theta_1)$, only $\lambda_0$ and $\sin(\theta_1)$ are truly independent variables.

$\sin$($\theta_1$) can be determined from the boundary condition stated above, since the remaining parameters are defined by choosing $\lambda_0$ alone.

$2\pi\cdot m=2\phi_t+2\phi_i$

$\rightarrow \phi_t=\frac{kh}{sin(\theta_1)}=\pi m -\phi_i$

$\rightarrow sin(\theta_1)=\frac{kh}{\pi m -\phi_i}$

The proposed algorithm actually computes the value of $m$ by probing sequential values of $\theta_1$ and determining the deviation of $m$ from a true integer value under consideration of a desired precision. It can be tested after a short introduction to the modesolver library below.

The modesolver library

Let's import the library which only has the Modedata class: (Enter the Code-containing field by marking it and pressing "ENTER" and then compile by pressing "CTRL"+"ENTER")

Use of tabulated data for refractive index values

The library is capable of importing the refractive index from any source, e.g. refractiveindex.info, if the data is provided in two separate text files for real and imaginary part, where the filename consists of the material and the suffix "R" or "I" for real- and imaginary, respectively. The files should each contain two columns with wavelength in microns and the respective value in column 1 and 2, respectively. If the filename cannot be found, generic values for air will be used (real=1, imaginary=0).

Example: Si: Two files: siliconR.dat and siliconI.dat.

Conntent:siliconR.dat

[lambda] [n] ( The first row will be skipped )

0.255 10.4

0.256 10.7

... ...

The importing function is called when the Modedata object is created and it looks for the tabulated data files which are given as names of the layers.

Let's create a Modedata object:

If the tabulated values cannot be found, generic ones are used with n=1, and $\kappa$=0.

What can the Modedata object do?

Let's calculate the Modedata

Let's plot the generated data to see what it looks like:

Small k VS Big K...

The question was raised if the extinction coefficient of a mode needs to take into account the propagation angle, since the absorption is per distance parallel to the waveguide, but the mode itself may zig-zag and thus have an effective pathlength > than the parallel distance. The alternate version of kappa can be activated by setting OBJ.smk=0, as shown below:

Let's compare a few points for their extinction coefficient values!

Great Stuff!! Now to make more precise calculations without waiting, we can import precompiled data:

Import / Export Modedata

Import and compare the Modesolver output of LUMERICAL

Metal Reflector, 1nm steps

The Lumerical Eigenmode Solver does not only output the guided modes, but all radiation that may interact with the waveguide. In order to filter for the waveguided modes, we discard all propagation constants $\beta$ that are smaller that $k_0$. Additionally, the data is brought to the same form as the output of our simple modesolver. This is done with the ImportLumerical function, which needs knowledge about the output file layout, i.e. columns per mode, col of beta and alpha. Please provide Thickness afterwards. Only with Lumerical input, another variable OBJ.all exists which holds all mode data, including non TIR modes. You can use these instead by issuing: OBJ.beta=OBJ.all

Metal Reflector, 10um thick silicon

Compare Lumerical Modes With Modedata generator

Calculating the expected absorption

Let us recover eq.$\,$(2.33) of the Master Thesis, which stems from a paper by Yu et al: Fundamental limit of nanophotonic light trapping in solar cells. The absorption spectrum for a single resonance of a waveguide is given in by:

\begin{equation} A(\omega)=\frac{\gamma_i\gamma_e}{\left(\omega-\omega_0\right)^2 +\left(\gamma_i+N\gamma_e\right)^2/4}, \,\,\,\,(1) \end{equation}

where $\gamma_i$ describes intrinsic absorption, $\gamma_e$ leakage to other channels (other resonance or free space radiation), $\omega_0$ is the frequency of the resonance, $N$ is the number of channels that are coupled with coupling constant $\gamma_e$.

Integration of the spectrum yields the spectral cross-section, which has unit of frequency. The interpretation as given by Yu et al is as follows: "For an incident spectrum with bandwidth $\Delta\omega\gg\sigma$, a resonance contributes an additional $\sigma/\Delta\omega$ to the spectrally-averaged absorption coefficient."

\begin{equation} \sigma=\int_{-\infty}^{\infty} A(\omega)d\omega%\leq\frac{2\pi\gamma_i}{N}. \,\,\,\, (2) \end{equation}

If we add the contributions from all individual resonances $m$ in the frequency range $[\omega,\omega+\Delta\omega]$, the upper limit for absorption is found in the equation below:

\begin{equation} A_\mathrm{T}=\frac{1}{\Delta\omega}\sum_m\frac{\sigma_\mathrm{max}}{N}=\frac{1}{\Delta\omega}\sum_m\frac{2\pi\gamma_{i,m}}{N}. \,\,\,\, (3) \end{equation}

The number of resonances M is given by (Yu et al., S10):

\begin{equation} M=\sum\limits_{m=1}^{m_\mathrm{tot}}\frac{2\pi\omega n_m^2}{c_0^2}\frac{L^2}{(2\pi)^2}\delta\omega. \,\,\,\,(4) \end{equation}

If the complex refractive index ($n+i\kappa$) for a medium is known, its intrinsic loss rate $\gamma_i$ can be calculated as follows(From Yu et al.):

\begin{equation} \gamma_i=\alpha_0\frac{c}{n}=2 k''\frac{c}{n}=\frac{4\pi}{\lambda}\kappa\frac{c}{n}. \,\,\,\, (5) \end{equation}

Now follows the characterization of the channel number $N$, which is assumed the same for all resonances. In the supporting information, the number of channels is approximated as an area fraction in k-space,

$N=\frac{A}{A_m}$,

where the total area $A$ is calculated from the interval of k-vectors $[k_1,k_2]$ that can be excited by the HUD-pattern,

$A=\pi\left(k_1^2-k_2^2\right)$,

and $A_m$ is the area occupied by a single resonance at

$A_m=(\Delta k_m)^2$.

Thus, \begin{equation} N=\frac{\pi\left(k_2^2-k_1^2\right)}{(\Delta k_m)^2}. \,\,\,\, (6) \end{equation} Substitution of $M,N$, and $\gamma_i$ into (2) yields:

\begin{equation} \tilde{A}=\frac{1}{\Delta\omega}\sum_\mathrm{m}\frac{4\pi\alpha_m k_m}{\left(k_2^2-k_1^2\right)}. \,\,\,\,(7) \end{equation}

Next, the analytical solutions for a lambertian scatterer are used to calculate the absorption in the slab. The formula is taken from M.A. Green's "Lambertian Light Trapping in Textured Solar Cells and Light-Emitting Diodes: Analytical Solutions".

There are several equations derived by Green that are more or less approximate depending on the magnitude of the absorption coefficient $\alpha$ times thickness $W$, since also the radial distribution of light intensity from a lambertian scatterer inside the absorber depends on the product $\alpha W$. The average pathlength of a light ray approaches $W$ for strong absorption $\alpha W>5$, but can increase up to $4W$ in the weak absorptive limit.

The product $\alpha W$ from Green can be identified with $\tilde{A}$ from (7) in our case, since it is unitless and the influence of the thickness is to determine the number of modes that take part in the sum over m.

We'll look at the magnitude of $A_T$ across the spectral range below and decide for the appropriate optical pathlength enhancement factor $f_p$ there. The total absorption $A_{TT}$ is then defined as:

\begin{equation} A_{T}(\omega)=\frac{1-\exp[-f_p \tilde{A}(\omega)]}{1-(1-n_m^{-1})\exp[-f_p \tilde{A}(\omega)]}. \,\,\,\, (8) \end{equation}

In order to arrive at the spectrally averaged integrated absorption, summation over $A_{T}(\omega)$ is performed, where each contribution is weighed by the flux of the incident radiation at that wavelength $F_\mathrm{AM1.5}(\omega)$, such that the integrated absorption:

\begin{equation} IA = \frac{1}{\int_\omega F(\omega) d\omega} \int_\omega A_{T}(\omega) F(\omega) d\omega . \,\,\,\, (9) \end{equation}

The algorithm for calculating the expected absorption

There are three code fields following: Starting with the

  1. Initialization, next come the

  2. definitions of subroutines / equations and functions, and finally some examplary

  3. plotfunctions for visualization.

The code blocks are best evaluated after one-another with CTRL+ENTER.

Initialization

Now follow the definitions of the subroutines and equations

Quick explanation of the functions above

extractKL(k1,k2,lamb)

#Extract all modes with certain wavelength: lamb and propagation constant: k1 <= beta <= K2 and return in separate Array

eq7a(lamb,k1,k2) #Equation 7, version as above for specific lambda

eq7b(lamb,k1,k2) #Equation 7, version with N_m from Femius' equation (cp. Thesis sec. 2.2) for specific lambda

eq8a(alpha) #Equation 8 with f_p=1, alpha the argument is supposed to be the return value of eq7

eq8b(alpha) #Equation 8 with f_p=4, alpha the argument is supposed to be the return value of eq7

IA(k1,k2,eq7,eq8)

#Spectrally averaged integrated absorption found through convolution of (8) with air mass 1.5 spectrum

Plotfunctions for the estimated absorption

Plot for some k1, k2 values absorption, IA, or IA for an array of k1, k2 values

Eq. (7) vs. $\alpha_0\cdot W$