## nextnanomat

**Software documentation**

**Operating System**

**nextnano.cloud**

## nextnano++

**Tutorials**

**Software documentation**

**Operating System**

**nextnano.cloud**

**Tutorials**

nnp:optics

Author: Takuma Sato, nextnano GmbH

In this tutorial we illustrate the optics{} module to demonstrate what nextnano++ can simulate for optoelectronic devices. This module performs a detailed calculation to **optical absorption phenomena**, **using 8 (or 6) band** $\mathbf{k}\cdot\mathbf{p}$ **models**. If you are interested in

- the background physics of this module and how to write the input file, go to
**Part 1** - the simulation results for intersubband transitions, go to
**Part 2** - the simulation results for interband transitions, go to
.**Part 3** - optical absorption in 2D devices,
*(under construction)* - optical absorption in broken-gap structures,
*(under construction)*

This algorithm is implemented based on the following diploma thesis:

- Thomas Eißfeller,
*Linear Optical Response of Semiconductor Nanodevices*, Technische Universität München (2008)

For the physics of optical transition in semiconductors and its application, we refer to

- Shun L. Chuang,
*Physics of Optoelectronic Devices*(Wiley, 1995) - S.M. Sze & Kwok K. Ng,
*Physics of Semiconductor Devices*(Wiley, 2007)

In the k.p analysis of one- (or two-) dimensional structures we have a projection of the Bloch wave vector along translation-invariant directions. We denote them as $\mathbf{k}_\parallel=k_y\hat{y}+k_z\hat{z}$ (1D) and $\mathbf{k}_\parallel=k_z\hat{z}$ (2D). Under envelope function approximation the $\mathbf{k}\cdot\mathbf{p}$ model yields the following equation to determine the confined states in structured directions

\begin{equation} \sum_{\mu=1}^8 H_{\nu\mu}^\mathrm{kp8}(\mathbf{k}_\parallel,\mathbf{r}_\perp)f_{m,\mu}(\mathbf{r}_\perp)=E_m(\mathbf{k}_\parallel)f_{m,\nu}(\mathbf{r}_\perp)\ \ \ (\nu=1,\dots,8), \end{equation}

where the Greek indices label the k.p bands and $m$ denotes eigenvalues, $\mathbf{r}_\perp=x\hat{x}$ in 1D and $\mathbf{r}_\perp=x\hat{x}+y\hat{y}$ in 2D. $H^\mathrm{kp8}$ is the 8$\times$8 matrix whose elements are given by the k.p parameters in the database. $f_{m,\mu}(\mathbf{r}_\perp)$ are the envelopes in the structured directions. The full wave function is given at each $\mathbf{k}_\parallel$ as

\begin{equation} \Psi_n(\mathbf{k}_\parallel,\mathbf{r})=\sum_{\mu=1}^8 F_{m,\mu}(\mathbf{k}_\parallel,\mathbf{r})u_\mu(\mathbf{r}) =\sum_{\mu=1}^8 \frac{e^{i\mathbf{k}_\parallel\cdot\mathbf{r}_\parallel}}{\sqrt{A}}f_{m,\mu}(\mathbf{r}_\perp)u_\mu(\mathbf{r}), \end{equation}

where $u_\mu(\mathbf{r})$ is the Bloch function of the band $\mu$ at $\mathbf{k}=0$ and $A=\int d\mathbf{r}_\parallel$. In general, both the conduction band ($\Gamma$) and valence bands contribute to this full wave function. The spinor composition is exported to `Quantum\spinor_composition`

. After solving this “Schrödinger” equation, the wave function is integrated over a limited region in $\mathbf{k}_\parallel$ space to obtain the charge density, which is used in the quantum-current-Poisson iteration. The region is specified under quantum{} as

quantum{ region{ kp_8band{ k_integration{ relative_size = $r_quantum # size of k||-space in quantum{} (relative to the Brillouin zone) num_points = $N_quantum # number of k|| points where Schödinger eq. is solved num_subpoints = $Nsub_quantum # number of points between k|| points where wave functions and eigenvalues are interpolated } } } }

When 1) Schrödinger equation is solved with k.p method, 2) optics{} flag is present and 3) the specifier `calculate_optics{}`

is present under run{} flag, nextnano++ calculates the absorption spectrum.

optics{ region{ ... # see below for details } } run{ solve_quantum{} calculate_optics{} }

The optical absorption accompanied by excitation of charge carriers (state $n\to m$) in a condensed matter is calculated from Fermi's golden rule [Chuang]. The absorption coefficient has the dimension of (length)$^{-1}$.
\begin{equation}
\alpha(\vec{\epsilon},\omega)=\frac{\pi e^2}{n_sc\varepsilon_0 m_0^2\omega}\frac{1}{V}\sum_{n>m}\sum_{\mathbf{k}_\parallel}
|\vec{\epsilon}\cdot\vec{\pi}_{nm}(\mathbf{k}_\parallel)|^2 (f_m-f_n)\delta(E_n-E_m-\hbar\omega),
\end{equation}
where the first sum runs over bands that fulfill $E_n>E_m$, and $f_m(\mathbf{k}_\parallel)=[1+e^{[E_m(\mathbf{k}_\parallel)-E_F]/k_BT}]^{-1}$ is the occupation of eigenstate $m$. When `optics{ occupation_ignore=yes }`

(default is `no`

), the program assumes
\begin{equation}
\begin{cases}
f_m(\mathbf{k}_\parallel)=0 \ \ \text{if}\ \ m\in\text{conduction band} \\
f_m(\mathbf{k}_\parallel)=1 \ \ \text{if}\ \ m\in\text{valence band}
\end{cases}
\end{equation}
The light polarization $\vec{\epsilon}$ and refractive index $n_s$ are specified in the input file. The refractive index is in general frequency-dependent, but we assume it to be constant and equal to the substrate value.

optics{ region{ polarization{ name="TM" re = [1,0,0] } # in 1D simulation, x is the growth direction polarization{ name="TE" re = [0,1,0] } # complex (circular) polarization is also allowed refractive_index = # (optional) use alternative value for the refractive index (default: substrate value) } }

The core of the optical transition is the **optical matrix elements** $\vec{\epsilon}\cdot\vec{\pi}_{nm}(\mathbf{k}_\parallel)$ because the kinematic momentum operator $\vec{\pi}=(\pi_x,\pi_y,\pi_z)$ couples linearly to the vector potential that accounts for the electromagnetic field. Here $\vec{\pi}$ is the sum of the conventional momentum operator $\mathbf{p}$ and the contribution of spin-orbit interaction. The optical matrix elements are calculated as
\begin{equation}
\vec{\pi}_{nm}(\mathbf{k}_\parallel)=\langle n|\vec{\pi}|m\rangle
=\int d\mathbf{r}
\begin{pmatrix}
F_{n1}^* & \cdots & F_{n8}^*
\end{pmatrix}
\begin{pmatrix}
&&\\
&\vec{\pi}_{\nu\mu}^\mathrm{kp8}&\\
&&
\end{pmatrix}
\begin{pmatrix}
F_{m1} \\ \vdots \\ F_{m8}
\end{pmatrix},
\end{equation}
where the 8$\times$8 matrix representation of the momentum operator, $\vec{\pi}_{\nu\mu}^\mathrm{kp8}$, has been derived using the Hellmann-Feynman theorem extended to the 8 band k.p model up to first order in $\mathbf{k}$ [Eißfeller].
For the analysis of the absorption spectrum, nextnano++ also prints out some fractions of the absorption coefficient formula in the output folder, namely

- occupation
`\Optics\occupation`

$f_m(\mathbf{k}_\parallel)$ - eigenvalue
`\Optics\eigenvalue`

$E_m(\mathbf{k}_\parallel)$ - transition intensity (if output_transitions=yes)
`\Optics\transition`

$T_{nm}(\vec{\epsilon},\mathbf{k}_\parallel)=\frac{2}{m_0}|\vec{\epsilon}\cdot\vec{\pi}_{nm}(\mathbf{k}_\parallel)|^2$ - imaginary part of the dielectric function
`\Optics\imepsilon`

$\mathrm{Im}\varepsilon_{nm}(\vec{\epsilon},\omega)=\frac{m_0}{2\omega^2}\frac{\pi e^2}{m_0^2\varepsilon_0}\frac{1}{V}\sum_{\mathbf{k}_\parallel}T_{nm}(\vec{\epsilon},\mathbf{k}_\parallel)(f_m-f_n)\delta(E_n-E_m-\hbar\omega)$ - absorption coefficient spectrum for each transition
`\Optics\absorption`

$\alpha_{nm}(\vec{\epsilon},\omega)=\frac{\omega}{n_sc}\mathrm{Im}\varepsilon_{nm}(\vec{\epsilon},\omega)$ - total absorption coefficient spectrum
`\Optics\absorption`

$\alpha(\vec{\epsilon},\omega)=\sum_{n>m}\alpha_{nm}(\vec{\epsilon},\omega)$

The following part of the input specifies how much transitions to be taken into account. The setting for `k_integration{}`

is explained in the next section.

optics{ region{ interband = $INTERBAND # yes or no intraband = $INTRABAND # yes or no energy_min = $ENERGY_MIN # minimum energy of the absorption spectrum energy_max = $ENERGY_MAX # maximum energy of the absorption spectrum energy_resolution = $ENERGY_RESOLUTION # energy grid spacing k_integration{ relative_size = $r_optics # size of k||-space in optics{} (relative to the Brillouin zone) num_points = $N_optics # number of k|| points where transition intensities are computed num_subpoints = $Nsub_optics # number of points between k|| points where transition intensity is interpolated } } }

Parameters in `k_integration{}`

in optics{} flag (hereafter $r_\mathrm{opt}, N_\mathrm{opt}, N'_\mathrm{opt}$) specify **the size and resolution of the $\mathbf{k}_\parallel$ space integration in absorption calculation**, $\sum_{\mathbf{k}_\parallel}$. This should not be confused with the specifier `k_integration{}`

in quantum{} flag used for **quantum mechanical charge density integration** (hereafter $r_q, N_q, N'_q$, see Figure 1).

First we discuss the parameters $r_\mathrm{opt}$ and $N_\mathrm{opt}$. The size of k|| space may affect the validity of simulation results. It also determines the simulation load. Here are some hints to determine the appropriate parameter sets:

- In undoped systems, integrating up to $|\mathbf{k}_\parallel|$ that gives in-plane kinetic energy $\hbar^2 k_\parallel^2/2m$ corresponding to $2k_BT$ or $3k_BT$ should be sufficient. Usually $r_\mathrm{opt}=0.3$ is sufficiently large to include all occupied states. In doped systems, it depends on the Fermi energy.
- To see the range of occupied states in k|| space, run a simulation and look at the output
`\Optics\occupation`

. We recommend to check the box “Show grid” on the left panel in Output tab of nextnanomat (see also nextnanomat documentation). This shows the occupation $f_m(\mathbf{k}_\parallel)$ as a function of $\mathbf{k}_\parallel$. Let us consider 1D simulation and suppose you got the following:

where $(r_\mathrm{opt},N_\mathrm{opt})=(0.3,8)$. The horizontal- and vertical axes are $k_y$ and $k_z$, respectively. The area $|k_{y,z}|\leq r_\mathrm{opt}\frac{\pi}{a}$ is shown with the k||-space gridding (thin white lines). The number of k|| points in one direction is $2N_\mathrm{opt}+3$. The occupation profile is not smooth and you might want a higher resolution by increasing the parameter $(r_\mathrm{opt},N_\mathrm{opt})\to(0.3,60)$:

The occupation becomes smooth, but at the same time this significantly increases the number of k points (in 1D simulation, (the number of k points)$\propto (r_\mathrm{opt}N_\mathrm{opt})^2$). Noting that the black region, where occupation is zero, does not contribute to the absorption, you can “zoom in” to the colored region by decreasing $r_\mathrm{opt}$ and $N_\mathrm{opt}$ in such a way that the ratio $r_\mathrm{opt}/N_\mathrm{opt}$ remains constant. This will cut down the irrelevant region without changing the resolution. For example, if you set $(r_\mathrm{opt},N_\mathrm{opt})=(0.05,10)$, you obtain

and this should be sufficient for the k||-space integration.

After tuning the parameters $r_\mathrm{opt}, N_\mathrm{opt}$, we can further optimize the setting regarding to the interpolation. The number of subpoints $N'_\mathrm{opt}$ determines at how many k|| points the transition intensity should be interpolated. Increasing $N'_\mathrm{opt}$ gives $E_m(\mathbf{k}_\parallel)$ of higher resolution and makes the absorption spectrum smooth. Figure 2 shows that this parameter improves the absorption coefficient spectrum.

*To do: investigate spin_degeneracy=yes/no and dipole_approximation = yes/no*

In the following we apply the formalism to several devices. As a first example, we simulate the absorption spectrum of an AlGaAs/GaAs quantum well infrared photodetector (QWIP). The QWIP is based on photoconductivity due to intersubband excitation.

- QWIP_singleQW_GaAs_AlGaAs_nnp.in
- QWIP_singleQW_InAs_AlSb_nnp.in
- QWIP(Gunapala_Levine_Chad-J.Appl.Phys70-1991)_nnp.in

The first example uses the same parameters used in

- FIG. 20 in B.F. Levine, J. Appl. Phys. 74 (8), 15 (1993),

while the third example is based on

- S.D. Gunapala
*et al*., J. Appl. Phys., 70 (1), 1 (1991)

We first illustrate the first example “QWIP_singleQW_GaAs_AlGaAs_nnp.in”. In this example, we simulate optical absorption in single quantum well structure. The following input is required for coupled quantum-current-Poisson simulation:

quantum{ region{ name = "optical_active" no_density = no kp_8band{ num_electrons = $OptNumE num_holes = $OptNumH } } } poisson{} current{} run{ solve_strain{} # strain calculation solve_current_poisson{} solve_quantum{} outer_iteration{} calculate_optics{} # absorption calculation }

The specifier `no_density=no`

lets the program calculate quantum mechanical charge density (default). Current-Poisson equation takes over this value. The bandstructure and wave functions are shown in Figure 3 and 4, respectively.

The output folder `\Optics`

contains simulation results for optical absorption. Let us first check the occupation $f_m(\mathbf{k}_\parallel)$ used in the absorption calculation. When comparing the results `\Optics\occupation`

, please mind the auto scale mode of nextnanomat:

The absorption coefficient for TE ($\vec{\epsilon}=\hat{y}$) and TM ($\vec{\epsilon}=\hat{x}$) light polarization is shown in Figure 6. The energy grid spacing here is `$ENERGY_RESOLUTION`

=0.5meV. For single-band models the peak becomes very sharp unless one introduces phenomenological broadening function such as Lorentzian. In k.p calculation, in contrast, peaks gets broadened because the transition energies, $E_n(\mathbf{k}_\parallel)-E_m(\mathbf{k}_\parallel)$, depends on k||. One can confirm this by comparing the output `\Optics\eigenvalue`

for states m=1 and 2 (not shown). In intersubband transitions the transition energies can be concave downward in k|| space, i.e., $E_n(\mathbf{k}_\parallel)-E_m(\mathbf{k}_\parallel)\propto -k^2$, depending on the masses. In the present case the absorption spectrum has a tail in the region $\hbar\omega<\Delta E$.

The optical transitions between conduction band states (intersubband transitions) in response to TE-polarized light is only allowed when eigenstates have finite spinor components in valence bands. In the present case its large bandgap and small confinement leads to small band-mixing, rendering TE absorption orders of magnitude smaller than TM polarization (Figure 6). As seen in the output `\Quantum\spinor_composition`

, eigenstates contain approximately 98% contribution from conduction band and 2% from valence band.

In the second example “QWIP_singleQW_InAs_AlSb_nnp.in”, single quantum well is narrower and the bandgap is smaller than the first example. The small bandgap and large confinement of the wave function (Figure 7) leads to large band mixing. In fact, the output `\Quantum\spinor_composition`

shows that the ground states in Figure 7 consists of 80.7% of conduction band and 19.3% of valence band contribution.

In the third example “QWIP(Gunapala_Levine_Chad-J.Appl.Phys70-1991)_nnp.in”, we set the bias to zero and impose the periodic boundary condition. The GaAs/Al$_x$Ga$_{1-x}$As supperlattice structure induces miniband states below the barriers, enabling bound-to-continuum absorptions of sub-eV photons. This $\mu\mathrm{m}$-wavelength photodetector works without electron tunneling through the barriers, thereby improving the detectivity [Gunapala]. The band structure `bandedges.dat`

and wave functions `\Quantum\probabilities_shift.dat`

are shown in Figure 9. We have continuum states above the barriers as well as bound states in the superlattice (miniband).

The absorption coefficient is exported to `\Optics\absorption`

. The indices in the filename “*_kp8_TE_m_n.dat” refer to the transition from state m to state n. The files without indices contain the total absorption (sum over all transitions). The total absorption coefficient for TE and TM polarization looks like this:

The peak positions do not depend on polarization, while the peak height is much larger for TM polarization compared to the one for TE. Looking at the absorption spectrum for each transition, we identify which transition contributes to which peak (Figure 11).

Let us look at the eigenvalue and occupation of each state to confirm this result. The eigenvalues of the bound- and continuum states are written in the output `\Quantum\probabilities_shift.dat`

or `\Quantum\energy_spectrum`

. **NOTE: quantum{} uses spin-resolved index for the eigenstates, so there are 80 states in total. In Optics{}, however, two spin-degenerate states are summed up and there are only 40 states. This number (1 to 40) is used in the \Optics output filenames. (To do: examine the specifier spin_degeneracy). For the consistency, we use the latter notation throughout.**
Based on the indices in Figure 11, we identify the first four peaks to the following four different transitions (Figure 12). We have confirmed that the peak energies in Figure 11 are consistent to the energy separation of the corresponding states.

Lastly we check the occupation (Fermi-Dirac distribution) $f_m(\mathbf{k}_\parallel)$. In the output `\Optics\eigenvaluespectrum`

(Figure 13), occupation at k||=0 of $m$-th state, $f_m(\mathbf{k}_\parallel=0)$, is plotted at corresponding eigenvalues $E_m$. The function takes the maximum value at the origin $\mathbf{k}_\parallel=0$. In the present system, $f_1(0)=0.087, f_2(0)=0.077, \dots, f_{10}(0)=0.0148$ for the bound states, whereas $f_m(0)<10^{-4}$ for continuum states ($m\geq 11$). Therefore the initial states in Figure 12 are well occupied and the final states are mostly empty. This enables optical absorption via bound-to-continuum excitation of electrons, thereby realizing a quantum well photodetector with high detectivity.

- AlGaAs_QW_Frankenberger_Simple_nnp.in
- AlGaAs_QW_Frankenberger_Simple_nnp_fast.in
- AlGaAs_QW_Frankenberger_Doping_schottky07_nnp.in
- AlGaAs_QW_Frankenberger_Doping_schottky07_nnp_fast.in

These files are located in the sample files folder. The *fast* examples reduce the computation load by limiting exact solution only to $k=0$ point and computing all other $k$ points in the basis of the $k=0$ wavefunctions (`force_k0_subspace`

; see quantum and optics documentations).

In the input file `AlGaAs_QW_Frankenberger_Simple_nnp.in`

, we consider a single quantum well structure:

The program solves the 8-band k.p model coupled to the Poisson equation to find the eigenstates and compute the absorption coefficient. Figure 15 shows the absorption coefficient spectrum for circularly polarized light ($\vec{\epsilon}=\hat{y}-i\hat{z}$). In contrast to QWIP examples above, peaks have long tails toward higher energy. This is because the transition energies $E_n(\mathbf{k}_\parallel)-E_m(\mathbf{k}_\parallel)$ in interband transitions are concave upward $\sim +k^2$ (here we do not consider Type 2 semiconductors).

The steps of this absorption spectrum are associated with the following interband transitions:

* Note*: In the end of the log file, you find the message “Integration reliable up to —eV”. This tells you up to which energy the absorption spectrum is reliable. Since we only consider the vicinity of the origin $\mathbf{k}_\parallel=0$, the reliable energy interval is bound from above by the energy difference of the initial and final states at the edge of the k||-space considered. The upper limit d [eV] is given by
\begin{equation}
d=\min_{\mathbf{k}_\parallel\in \ \Omega^*\mathrm{\ edge}} |E_n(\mathbf{k}_\parallel)-E_m(\mathbf{k}_\parallel)|,
\end{equation}
where $\Omega^*$ is the region in k||-space specified in

`optics{ region{ k_integration{} } }`

with parameters $r_\mathrm{opt}$ and $N_\mathrm{opt}$. In the present case d=3.2eV, while the calculation is safely performed for the interval [1.4, 1.7] (eV). This message appears only when interband transitions are computed, i.e. when `interband=yes`

and `intraband=no`

in optics{} flag.
In the second input file `AlGaAs_QW_Frankenberger_Doping_schottky07_nnp.in`

, we consider the following structure:

Figure 19 compares the results for different settings for occupation $f_m(\mathbf{k}_\parallel)$. When `optics{ occupation_ignore=yes }`

, valence bands and conduction bands are considered to be fully occupied and fully empty, respectively. When the actual occupation of eigenstates are taken into account, in contrast, optical transitions to conduction band states just above the Fermi energy are prohibited because of the thermal distribution of electrons.

nnp/optics.txt · Last modified: 2021/05/18 17:35 by takuma.sato