format: text

text:
# Nonequilibrium Green’s function theory for transport and gain properties of quantum cascade structures

## I Introduction

Quantum cascade (QC) structures are semiconductor heterostructures grown with a complicated sequence of alternating layers of different semiconductor materials and with varying thicknesses. This sequence of layers is repeated many times, up to tens or even over a hundred periods. Figure 1 shows an example of the conduction-band line-up in a QC structure. These structures form the basis of a new type of semiconductor laser,Fai94b  in which the laser light emission occurs through intersubband or interminiband transitions (in most cases within the conduction band) rather than interband transitions. These lasers have a great variety of designs, and a recent review is given in Ref. Gma01, . Until recently,Wan01  all quantum cascade laser (QCL) structures were designed so that each period in the structure contains an active region in which the lasing transition occurs, and a separate injector region. The injector acts as a reservoir of electrons for injection into the active region of the next stage. It also acts as a collector of electrons from the preceding active region. The direction of the electron flow is indicated in Fig. 1, and the electron flow is seen to resemble a cascade as the electrons move from one stage to the next when a bias is applied. Hence, the electron transport through a QCL structure is a complicated interplay between relaxation through light emission and scattering in the active region, and transmission through tunneling and scattering in the injector region.Iot01

Initially, most theoretical investigationsIot01 ; Tor00 ; Har99a  of transport in QC structures have focused on the role played by scattering processes in determining transport properties and the dynamics of the electron distributions in these structures. Very recently, however, the first theoretical investigations of quantum transport in these structures which have considered or incorporated quantum effects beyond a Boltzmann equation approach have been reported.Wac01 ; Iot01a A theoretical study of quantum transport may be treated using the density matrix formalism,Kuh98  or with a nonequilibrium Green’s function (NGF) approach.Kad62 ; Hau96 ; Lak97 ; Wac02

In the work reported here, we extend the NGF theory described in Ref. Wac02,  to the study of quantum transport in QC structures. Very early results from this investigation have been reported in Refs. Wac01,  and Lee02, . In this paper, we present further and more detailed results from this study. In addition, we have extended the work further to the problem of evaluating the gain or absorption spectra of these structures,Wac02a  and we also report and discuss results of this work here.

In the following section, we describe the theoretical formulation that we use to derive the transport properties, and then the linear optical response of QC structures. In Sec. III, we apply this theory to example QC structures, and describe the results obtained. The last section contains a summary and conclusion.

## II Theoretical Formulation

### II.1 Basis states and Hamiltonian

We model the QC structure as a periodic superlattice structure, in which each period contains $`N_s`$ semiconductor layers with varying thicknesses. The Hamiltonian $`\widehat{H}`$ which we use to model this system may be separated into two parts: $`\widehat{H}=\widehat{H}_o+\widehat{H}_{\mathrm{scatt}}`$. $`\widehat{H}_o`$ contains the superlattice potential and a static electric field $``$ applied in the growth direction, i.e., $`\widehat{H}_o=\widehat{H}_{\mathrm{SL}}+\widehat{H}_{}`$. The Hamiltonian $`\widehat{H}_{\mathrm{scatt}}`$ describes the scattering processes included in the theory.

The Hamiltonian is expressed in a set of basis states, $`\mathrm{\Psi }_{𝐤,\alpha }(𝐫,z)=(e^{i𝐤.𝐫}/\sqrt{𝒜})\psi _\alpha (z)`$, which we assume separable, although this is an approximation when the effective mass is position dependent. In the plane of the semiconductor layers, the basis functions behave as plane waves. The normalization constant $`𝒜`$ is the sample area in this plane. The position vector $`𝐫`$ and in-plane wave vector k are two-dimensional (2D) vectors. In the growth direction, here labeled $`z`$, there are several possible choices (see Ref. Wac02,  for a discussion) for the functions $`\psi _\alpha (z)`$. Although the physical results obtained from the theory should be independent of the specific choice of these functions, there are different advantages or drawbacks attached to a given choice. For instance, as we discuss next, one choice may be more suited to expediting the numerical computation, while another choice may more easily allow the extraction of physical information in a form that can be compared with experimental measurements. Possible choices are (i) Bloch functions which are eigenstates of the bare superlattice potential $`\widehat{H}_{\mathrm{SL}}`$, and are spatially extended across the whole structure. (ii) Wannier-Stark states which are eigenstates of $`\widehat{H}_o`$, i.e., of the superlattice potential and the applied bias. These eigenstates are metastableGlu02  (non-Hermitian) with complex energies. They are often treated approximately as stationary states (their metastable nature is neglected), and this leads to an ambiguity in the definition of these states depending on how this approximation is made. (iii) Wannier functions, which should not be confused with Wannier-Stark functions, may be constructed so that they are spatially well localized. In particular, they may be constructed as eigenstates of the position operator. In symmetric superlattice structures, these eigenstates are maximally localized.Koh59  The spatial localization of the Wannier states enables one to setup a picture of transport in which scattering occurs from one spatial region of the structure to another. The Wannier functions do not depend on bias, and this is computationally advantageous since coupling matrix elements between Wannier functions may be calculated only once at zero bias, and the same matrix elements may then be used at any applied bias. For these reasons, the theory presented here is formulated in the Wannier function basis (see Appendix A), which we derive as eigenstates of the position operator. A disadvantage in using this basis is that the Wannier functions are not energy eigenstates of the structure, and it is therefore difficult to obtain a physical interpretation in the energy picture or domain. For instance, there is no one-to-one correspondence between an optical transition in the structure and a transition between a pair of Wannier functions. This difficulty may be partially circumvented, however, at a later stage in the calculation, by transforming the results obtained into the Wannier-Stark basis. We emphasize again that the basis choice is not in itself an important issue, i.e., the physical content of the theory should not depend on this choice, but a suitable choice of basis can facilitate the numerical computation (e.g., Wannier states), or more easily allow physical interpretation (e.g., Wannier-Stark states).

Expressing the Hamiltonian in the Wannier basis we obtain

$`\widehat{H}_{\mathrm{SL}}`$ $`=`$ $`{\displaystyle \underset{n,\nu }{}}{\displaystyle \underset{𝐤,𝐬}{}}[E_{\nu 𝐤}\widehat{a}_{n,𝐤,s}^\nu \widehat{a}_{n,𝐤,s}^\nu `$ (1)
$`+T_1^\nu (\widehat{a}_{n+1,𝐤,s}^\nu \widehat{a}_{n,𝐤,s}^\nu +\widehat{a}_{n1,𝐤,s}^\nu \widehat{a}_{n,𝐤,s}^\nu )],`$
where the index $`n`$ labels a period in the superlattice, and the index $`\nu `$ labels a Wannier function $`\psi _\nu (z)`$ within a period. $`\widehat{a}_{n,𝐤,s}^\nu `$ and $`\widehat{a}_{n,𝐤,s}^\nu `$ are creation and annihilation operators for an electron with in-plane wave vector k, and spin index $`s`$, in the $`\nu `$th Wannier level, in period $`n`$. As stated earlier, the Wannier functions are not eigenstates of $`\widehat{H}_{\mathrm{SL}}`$, and hence $`T_l^\nu `$ represents the off-diagonal couplings between Wannier levels in different periods, and $`E_{\nu 𝐤}`$ represents the diagonal elements of $`\widehat{H}_{\mathrm{SL}}`$ in this basis. We keep only terms in $`T_1^\nu `$, i.e., we consider only couplings between adjacent periods. The next-nearest-neighbor couplings $`T_2^\nu `$ are two or more orders of magnitude smaller. The Hamiltonian $`H_{}`$, due to the electric field $``$, is written as

$`\widehat{H}_{}`$ $`=`$ $`{\displaystyle \underset{n,\nu ,\mu }{}}{\displaystyle \underset{𝐤,s}{}}\{eR_0^{\mu \nu }\widehat{a}_{n,𝐤,s}^\mu \widehat{a}_{n,𝐤,s}^\nu ned\delta _{\mu \nu }\widehat{a}_{n,𝐤,s}^\mu \widehat{a}_{n,𝐤,s}^\mu `$ (2)
$`eR_1^{\mu \nu }[\widehat{a}_{n+1,𝐤,s}^\mu \widehat{a}_{n,𝐤,s}^\nu +\widehat{a}_{n,𝐤,s}^\nu \widehat{a}_{n+1,𝐤,s}^\mu ]\},`$
where $`R_l^{\mu \nu }=𝑑z\psi _\mu ^{}(zld)z\psi _\nu (z)`$. $`d`$ is the length of one period and $`e<0`$ is the electron charge.

In the scattering Hamiltonian $`\widehat{H}_{\mathrm{scatt}}`$, we include interface roughness, impurity, and electron-phonon scattering processes. Both acoustic phonons and longitudinal optical (LO) phonons are considered. The electron–acoustic-phonon scattering rates are far smaller than those of the other scattering processes (Table 1), but we include these to provide a channel for small energy transfers between the electrons and lattice. This opens up a route for the carrier-lattice system to move towards thermal equilibrium, particularly at low bias when the electrons may have insufficient energy to emit LO-phonons.

In the NGF theory, the scattering processes expressed in $`\widehat{H}_{\mathrm{scatt}}`$ are treated in the form of self-energies which are described below. We have also included electron-electron scattering in a mean-field approximation so that the interaction between electrons appears as an additional single-particle potential in $`\widehat{H}_{\mathrm{scatt}}`$, i.e., we replace $`\widehat{H}_{\mathrm{scatt}}`$ with $`\widehat{H}_{\mathrm{scatt}}^{}=\widehat{H}_{\mathrm{scatt}}+\widehat{H}_{\mathrm{MF}}`$. $`\widehat{H}_{\mathrm{MF}}`$ is written as

$$\widehat{H}_{\mathrm{MF}}=\underset{m\mu ,n\nu }{}\underset{𝐤,s}{}[V_{\mathrm{MF}}]_{m\mu ,n\nu }\widehat{a}_{m,𝐤,s}^\mu \widehat{a}_{n,𝐤,s}^\nu ,$$
(3)
where $`[V_{\mathrm{MF}}]_{m\mu ,n\nu }=𝑑z\psi _{m\mu }^{}(z)V_{\mathrm{MF}}(z)\psi _{n\nu }(z)`$. The evaluation of the mean-field potential $`V_{\mathrm{MF}}(z)`$ is described in the following subsection.

### II.2 Quantum transport equations and self-energies

To determine the transport properties of a QC structure under applied bias, we need first to obtain a description of the nonequilibrium stationary state of the system. This information is contained in the nonequilibrium Green’s functions: the retarded Green’s function $`𝐆^{\mathrm{ret}}(E)`$ and the correlation function $`𝐆^<(E)`$ (bold type G, and $`𝚺`$ below, indicates a matrix in, e.g., the Wannier basis). The quantum transport equations, which these functions obey, can be derivedHau96 ; Wac02  from the Hamiltonian $`\widehat{H}`$ and have the following form. The Dyson equation for the matrix element $`G_{\alpha _1\alpha _2,𝐤}^{\mathrm{ret}}(E)`$ is written as

$$EG_{\alpha _1\alpha _2,𝐤}^{\mathrm{ret}}(E)\underset{\beta }{}\left[(\widehat{H}_o+\widehat{H}_{\mathrm{MF}})_{\alpha _1\beta ,𝐤}+\mathrm{\Sigma }_{\alpha _1\beta ,𝐤}^{\mathrm{ret}}(E)\right]G_{\beta \alpha _2,𝐤}^{\mathrm{ret}}(E)=\delta _{\alpha _1\alpha _2},$$
(4)
where $`\alpha _1`$, $`\alpha _2`$, and $`\beta `$ are general indices that include both the period and Wannier level indices, e.g. $`\alpha _1(m,\mu )`$. The correlation function $`G_{\alpha _1\alpha _2,𝐤}^<(E)`$ is obtained from the Keldysh relation:

$$G_{\alpha _1\alpha _2,𝐤}^<(E)=\underset{\beta \beta ^{}}{}G_{\alpha _1\beta ,𝐤}^{\mathrm{ret}}(E)\mathrm{\Sigma }_{\beta \beta ^{},𝐤}^<(E)G_{\beta ^{}\alpha _2,𝐤}^{\mathrm{adv}}(E),$$
(5)
where $`𝐆^{\mathrm{adv}}(E)=[𝐆^{\mathrm{ret}}(E)]^{}`$.

The self-energies $`𝚺^{\mathrm{ret}}(E)`$ and $`𝚺^<(E)`$ originate from the scattering processes contained in $`\widehat{H}_{\mathrm{scatt}}`$ (i.e., excluding $`\widehat{H}_{\mathrm{MF}}`$), and are evaluated within the self-consistent Born approximation. Assuming that the diagonal parts of the Green’s functions and the self-energies dominate (because the basis states are spatially localized), the self-energy for interface roughness or impurity scattering has the form

$$\mathrm{\Sigma }_{\alpha _1\alpha _1,𝐤}^{<,\mathrm{rough}/\mathrm{imp}}(E)=\underset{\beta ,𝐤^{}}{}|V_{\alpha _1\beta }^{\mathrm{rough}/\mathrm{imp}}(𝐤𝐤^{})|^2G_{\beta \beta ,𝐤^{}}^<(E).$$
(6)
$`V_{\alpha _1\beta }^{\mathrm{rough}}(𝐤^{}𝐤)`$ and $`V_{\alpha _1\beta }^{\mathrm{imp}}(𝐤^{}𝐤)`$ are matrix elements for interface roughness and impurity scattering. In the former case, the angle brackets denote an average over all possible distributions of variations in the interface thickness, and in the latter case these brackets denote an average over all possible impurity distributions. The equation for $`𝚺^{\mathrm{ret},\mathrm{rough}/\mathrm{imp}}(E)`$ is obtained by making the replacements: $`𝚺^{<,\mathrm{rough}/\mathrm{imp}}(E)𝚺^{\mathrm{ret},\mathrm{rough}/\mathrm{imp}}(E)`$, and $`𝐆^<(E)𝐆^{\mathrm{ret}}(E)`$.

For optical or acoustic phonon scattering, the self-energies are

$`\mathrm{\Sigma }_{\alpha _1\alpha _1,𝐤}^{\mathrm{ret},\mathrm{phon}}(E)`$ $`=`$ $`{\displaystyle \underset{\beta ,𝐤^{}}{}}|V_{\alpha _1\beta }^{\mathrm{phon}}(𝐤,𝐤^{})|^2`$ (7)
$`\times \{[f_B(E_{\mathrm{phon}})+1]G_{\beta \beta ,𝐤^{}}^{\mathrm{ret}}(EE_{\mathrm{phon}})`$
$`+f_B(E_{\mathrm{phon}})G_{\beta \beta ,𝐤^{}}^{\mathrm{ret}}(E+E_{\mathrm{phon}})`$
$`+{\displaystyle \frac{1}{2}}\left[G_{\beta \beta ,𝐤^{}}^<(EE_{\mathrm{phon}})G_{\beta \beta ,𝐤^{}}^<(E+E_{\mathrm{phon}})\right]`$
$`+i{\displaystyle }{\displaystyle \frac{dE_1}{2\pi }}G_{\beta \beta ,𝐤^{}}^<(EE_1)(𝒫\left\{{\displaystyle \frac{1}{E_1E_{\mathrm{phon}}}}\right\}𝒫\left\{{\displaystyle \frac{1}{E_1+E_{\mathrm{phon}}}}\right\})\},`$
and

$`\mathrm{\Sigma }_{\alpha _1\alpha _1,𝐤}^{<,\mathrm{phon}}(E)={\displaystyle \underset{\beta ,𝐤^{}}{}}|V_{\alpha _1\beta }^{\mathrm{phon}}(𝐤,𝐤^{})|^2`$
$`\times \{f_B(E_{\mathrm{phon}})G_{\beta \beta ,𝐤^{}}^<(EE_{\mathrm{phon}})`$
$`+[f_B(E_{\mathrm{phon}})+1]G_{\beta \beta ,𝐤^{}}^<(E+E_{\mathrm{phon}})\},`$ (8)
where $`V^{\mathrm{phon}}`$ represents the interaction with optical or acoustic phonons, and $`E_{\mathrm{phon}}`$ represents the energy of the optical or acoustic phonon. $`f_B(E)=1/[\mathrm{exp}(E/k_BT)1]`$ is the equilibrium phonon distribution at a temperature $`T`$. This is the only place where the lattice temperature appears in this theory. Further detail about the evaluation of the scattering matrix elements and the self-energies is given in Appendix B. For simplicity in the numerical evaluation, we neglect the last line containing the principle value terms in Eq. (7). To further expedite the numerical computation, we use momentum-independent scattering matrix elements $`V_{\alpha \beta }^{\mathrm{phon}/\mathrm{rough}/\mathrm{imp}}(𝐤_{\mathrm{typ}},𝐤_{}^{}{}_{\mathrm{typ}}{}^{})`$, employing a representative momentum $`𝐤_{\mathrm{typ}}`$ (see Appendix B).

As mentioned in the previous section, we include also electron-electron scattering in the form of a mean-field potential which we determine by solving Poisson’s equation

$$\frac{d^2V_{\mathrm{MF}}(z)}{dz^2}=\frac{e}{ϵ_s}[\rho _{\mathrm{dope}}(z)+\rho _e(z)],$$
(9)
with periodic boundary conditions $`V_{\mathrm{MF}}(z+d)=V_{\mathrm{MF}}(z)`$. $`ϵ_s`$ is the absolute static background permittivity. $`\rho _{\mathrm{dope}}(z)`$ is the dopant density profile in the structure, and the electron density profile

$$\rho _e(z)=\mathrm{\hspace{0.33em}2}e\underset{\alpha \beta }{}\underset{𝐤}{}\frac{dE}{2\pi }iG_{\alpha \beta ,𝐤}^<(E)\psi _\beta ^{}(z)\psi _\alpha (z),$$
(10)
with the factor of 2 from summation over the spin index.

To determine the Green’s functions, $`𝐆^{\mathrm{ret}}(E)`$ and $`𝐆^<(E)`$, we loop iteratively between the quantum transport equations, Eqs. (4) and (5), and the equations for the self-energies and mean-field potential Eqs. (6) – (10), until we reach a self-consistent solution for these equations. To solve these equations, we assume that the system is an infinite periodic structure with period $`d`$. The boundary conditions $`G_{\alpha \beta }(E)=G_{\alpha ^{}\beta ^{}}(E+ed)`$ are applied between each period to model the bias drop per period in the structure, where $`\alpha ^{}`$ and $`\beta ^{}`$ are shifted one period to the left of $`\alpha `$ and $`\beta `$. As starting conditions for the iterative loop, we assume the electrons are in a Fermi distribution with temperature $`T`$, and that the density of states in each Wannier subband is a simple stepfunction. However, once a converged solution was obtained for a given applied voltage, we used this solution as a starting point for the calculation at the next bias to facilitate the convergence process. To test for convergence, we compared the differences in the $`𝐤`$-integrated Green’s functions $`[𝐆_{\alpha \alpha }^{\mathrm{ret},<}(E)]_{i+1}[𝐆_{\alpha \alpha }^{\mathrm{ret},<}(E)]_i`$, and carrier density $`[n_e]_{i+1}[n_e]_i`$, evaluated in two successive iteration steps, with a given tolerance value. A further test for convergence could be carried out by starting the calculation at different bias points, e.g., starting at zero bias and increasing the voltage to a high bias point, then restarting the calculation at a high bias point and decreasing the voltage to a low value.

### II.3 Transport properties

Having determined $`𝐆^{\mathrm{ret}}(E)`$ and $`𝐆^<(E)`$, information about the stationary state of the system, such as excitation energies, energy renormalizations, broadenings or lifetimes, density of states, populations, and distribution functions can be extracted. Examples of such information and their application are given in Sec. III.1. We can also evaluate the current density using the rate of change of the position operator $`\widehat{z}`$

$$J(t)=\frac{e}{𝒱}\frac{d\widehat{z}}{dt}=\frac{e}{𝒱}\frac{i}{\mathrm{}}[\widehat{H},\widehat{z}].$$
(11)
$`J(t)`$ is the average current density flowing in the $`z`$ direction, in system volume $`𝒱`$. Recalling the separation of the Hamiltonian $`\widehat{H}`$ into two parts, $`\widehat{H}_o`$ and $`\widehat{H}_{\mathrm{scatt}}^{}`$, we can determine two separate contributions to the current $`J(t)`$. The first contribution comes from $`\widehat{H}_o`$,

$$J_o(t)=\frac{e}{𝒱}\frac{i}{\mathrm{}}[\widehat{H}_o,\widehat{z}]=\frac{e}{\mathrm{}𝒱}\mathrm{Tr}\{[\widehat{H}_o,\widehat{z}]𝐆^<(t,t)\},$$
(12)
which we evaluate in the energy representation

$$J_o=\frac{2e}{\mathrm{}𝒱}\underset{\alpha \beta }{}\underset{𝐤}{}\frac{dE}{2\pi }[\widehat{H}_o,\widehat{z}]_{\alpha \beta }G_{\beta \alpha ,𝐤}^<(E).$$
(13)
The second contribution is from $`\widehat{H}_{\mathrm{scatt}}^{}`$,

$`J_{\mathrm{scatt}}^{}(t)`$ $`=`$ $`J_{\mathrm{scatt}}(t)+J_{\mathrm{MF}}(t)`$ (14)
$`=`$ $`{\displaystyle \frac{e}{𝒱}}{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{scatt}}+\widehat{H}_{\mathrm{MF}},\widehat{z}],`$
where we evaluate $`J_{\mathrm{MF}}`$ similarly to $`J_o`$, i.e., Ref. foot.jmf,

$$J_{\mathrm{MF}}=\frac{2e}{\mathrm{}𝒱}\underset{\alpha \beta }{}\underset{𝐤}{}\frac{dE}{2\pi }[\widehat{H}_{\mathrm{MF}},\widehat{z}]_{\alpha \beta }G_{\beta \alpha ,𝐤}^<(E).$$
(15)
The current contribution driven by $`H_{\mathrm{scatt}}`$ is written in the energy representation as (see Appendix C)

$`J_{\mathrm{scatt}}`$ $`=`$ $`{\displaystyle \frac{2e}{\mathrm{}𝒱}}{\displaystyle \underset{\alpha \beta \gamma }{}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{dE}{2\pi }[G_{\beta \gamma ,𝐤}^<(E)\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{\mathrm{adv}(\alpha )}(E)+G_{\beta \gamma ,𝐤}^{\mathrm{ret}}(E)\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{<(\alpha )}(E)]z_{\gamma \beta }}`$ (16)
$`z_{\alpha \gamma }[\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{<(\beta )}(E)G_{\gamma \alpha ,𝐤}^{\mathrm{adv}}(E)+\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{\mathrm{ret}(\beta )}(E)G_{\gamma \alpha ,𝐤}^<(E)].`$
The superscript notation $`(\alpha )`$ in, for example, the self-energy $`\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{<(\alpha )}(E)`$ indicates that we take only the part of the self-energy which depends on $`G_{\alpha \alpha }(E)`$. The factor of 2 in Eqs. (13), (15), and (16) is from the spin index summation.

We note here, that although we have divided the current into separate contributions driven by $`\widehat{H}_o`$, $`\widehat{H}_{\mathrm{MF}}`$, and $`\widehat{H}_{\mathrm{scatt}}`$, this should not be taken to mean that the scattering processes play no role in driving, for example, the current $`J_o`$. Although, $`\widehat{H}_{\mathrm{scatt}}`$ does not appear explicitly in Eq. (12), it is instrumental in determining $`𝐆^<`$, and, hence, implicitly influences the current $`J_o`$. The scattering processes provide channels for energy dissipation from the electronic system, so that as the electrons descend downwards through the electric potential their kinetic energy does not increase but remains constant, and a steady current flow is maintained through the structure.

### II.4 Gain and absorption spectra

We can also use the information contained in $`𝐆^{\mathrm{ret}}(E)`$ and $`𝐆^<(E)`$ as a starting point to evaluate the gain or absorption spectra of the structure. To do this, we consider the linear response of the stationary state described by $`𝐆^{\mathrm{ret}}(E)`$ and $`𝐆^<(E)`$ to a small perturbation from an optical field. From this response we determine the susceptibility $`\chi (\omega )`$, and from Maxwell’s equationsHau93  we obtain the gain coefficient

$$g(\omega )\frac{\omega }{c}\frac{\mathrm{Im}[\chi (\omega )]}{n_B},$$
(17)
where $`n_B\sqrt{ϵ_s}`$ is the background refractive index. The applied optical field

$$𝑬(𝐫,t)=𝐞_z\frac{d\omega }{2\pi }(\omega )e^{ik(\omega )yi\omega t},$$
(18)
which propagates in the $`y`$ direction is included as a small time-dependent perturbation to the Hamiltonian: $`\widehat{H}\widehat{H}+\delta \widehat{V}(t)`$, where the perturbation

$$\delta \widehat{V}(𝐫,t)=\frac{e}{2m_e(z)}\left[\widehat{𝐩}𝐀(𝐫,t)+𝐀(𝐫,t)\widehat{𝐩}\right]+e\varphi (𝐫,t),$$
(19)
with the momentum operator $`\widehat{𝐩}`$, and $`m_e(z)`$ is the spatially dependent effective mass. The vector potential $`𝐀(𝐫,t)`$ and scalar potential $`\varphi (𝐫,t)`$ are related to the optical field through

$`𝐀(𝐫,t)`$ $`=`$ $`{\displaystyle \frac{d\omega }{2\pi }\frac{(\omega )}{i\omega }e^{ik(\omega )yi\omega t}𝐞_z},`$
$`\varphi (𝐫,t)`$ $`=`$ $`0,`$ (20)
in the Coulomb gauge. A discussion of the use of different gauges in this problem is given in Ref. Wac02b, . For the results presented here, we apply the long-wavelength approximation, in which we neglect the spatial variation of the optical field across the structure. Hence, in this gauge, the perturbation is written in the energy representation as

$$\delta V_{\alpha \beta }(\omega )=i\frac{e(\omega )}{\omega }\left[\frac{\widehat{p}_z}{m_e(z)}\right]_{\alpha \beta }=\frac{e(\omega )}{\mathrm{}\omega }[\widehat{H}_o,\widehat{z}]_{\alpha \beta }.$$
(21)
The linear changes in the Green’s functions and self-energies caused by the additional term $`\delta \widehat{V}(t)`$ in the Hamiltonian represent the linear response of the nonequilibrium stationary state to the applied optical field. These changes may be written asWac02b

$$\delta 𝐆^{\mathrm{ret}}(\omega ,E)=𝐆^{\mathrm{ret}}(E+\mathrm{}\omega )[\delta 𝐕(\omega )+\delta 𝚺^{\mathrm{ret}}(\omega ,E)]𝐆^{\mathrm{ret}}(E)$$
(22)
$$\delta 𝐆^{\mathrm{adv}}(\omega ,E)=𝐆^{\mathrm{adv}}(E+\mathrm{}\omega )[\delta 𝐕(\omega )+\delta 𝚺^{\mathrm{adv}}(\omega ,E)]𝐆^{\mathrm{adv}}(E),$$
(23)
and

$`\delta 𝐆^<(\omega ,E)`$ $`=`$ $`𝐆^{\mathrm{ret}}(E+\mathrm{}\omega )\delta 𝐕(\omega )𝐆^<(E)+𝐆^<(E+\mathrm{}\omega )\delta 𝐕(\omega )𝐆^{\mathrm{adv}}(E)`$ (24)
$`+𝐆^{\mathrm{ret}}(E+\mathrm{}\omega )\delta 𝚺^{\mathrm{ret}}(\omega ,E)𝐆^<(E)+𝐆^{\mathrm{ret}}(E+\mathrm{}\omega )\delta 𝚺^<(\omega ,E)𝐆^{\mathrm{adv}}(E)`$
$`+𝐆^<(E+\mathrm{}\omega )\delta 𝚺^{\mathrm{adv}}(\omega ,E)𝐆^{\mathrm{adv}}(E).`$
Note that $`\delta 𝐆^{\mathrm{adv}}(\omega ,E)[\delta 𝐆^{\mathrm{ret}}(\omega ,E)]^{}`$. The expressions for $`\delta 𝚺(E)`$ take the same form as in Eqs. (6) – (8) where the functions $`𝐆(E)`$ are replaced by $`\delta 𝐆(E)`$. For $`\delta 𝚺^{\mathrm{adv}}(E)`$, the expression for $`\delta 𝚺^{\mathrm{ret}}(E)`$ is used with the replacement $`\delta 𝐆^<(E)\delta 𝐆^<(E)`$ and $`\delta 𝐆^{\mathrm{ret}}(E)\delta 𝐆^{\mathrm{adv}}(E)`$. $`\delta 𝚺^{</\mathrm{ret}}(E)`$ must be evaluated self-consistently with $`\delta 𝐆^{</\mathrm{ret}}(E)`$, and therefore another iterative loop must be carried out, similar to the earlier procedure to determine $`𝐆^{</\mathrm{ret}}(E)`$ and $`𝚺^{</\mathrm{ret}}(E)`$. From the change in the Green’s functions and self-energies, and using Eqs. (11) – (16), we obtain the change in current density

$$\delta J(\omega )=\delta J_o(\omega )+\delta J_{\mathrm{MF}}(\omega )+\delta J_{\mathrm{scatt}}(\omega ),$$
(25)
wherefoot.jmf

$`\delta J_o(\omega )`$ $`=`$ $`{\displaystyle \frac{2e}{\mathrm{}𝒱}}{\displaystyle \underset{\alpha \beta }{}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{dE}{2\pi }[\widehat{H}_o,\widehat{z}]_{\alpha \beta }\delta G_{\beta \alpha ,𝐤}^<(\omega ,E)}+[\delta V,\widehat{z}]_{\alpha \beta }G_{\beta \alpha ,𝐤}^<(\omega ,E),`$ (26)
$`\delta J_{\mathrm{MF}}(\omega )`$ $`=`$ $`{\displaystyle \frac{2e}{\mathrm{}𝒱}}{\displaystyle \underset{\alpha \beta }{}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{dE}{2\pi }[\widehat{H}_{\mathrm{MF}},\widehat{z}]_{\alpha \beta }\delta G_{\beta \alpha ,𝐤}^<(\omega ,E)},`$ (27)
and

$`\delta J_{\mathrm{scatt}}(\omega )`$ $`=`$ $`{\displaystyle \frac{2e}{\mathrm{}𝒱}}{\displaystyle }{\displaystyle \frac{dE}{2\pi }}{\displaystyle \underset{\alpha \beta \gamma }{}}{\displaystyle \underset{𝐤}{}}[\delta G_{\beta \gamma ,𝐤}^<(\omega ,E)\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{\mathrm{adv}(\alpha )}(E)+\delta G_{\beta \gamma ,𝐤}^{\mathrm{ret}}(\omega ,E)\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{<(\alpha )}(E)`$ (28)
$`+G_{\beta \gamma ,𝐤}^<(E+\mathrm{}\omega )\delta \mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{\mathrm{adv}(\alpha )}(\omega ,E)+G_{\beta \gamma ,𝐤}^{\mathrm{ret}}(E+\mathrm{}\omega )\delta \mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{<(\alpha )}(\omega ,E)]z_{\gamma \beta }`$
$`{\displaystyle \underset{\alpha \beta \gamma }{}}{\displaystyle \underset{𝐤}{}}z_{\alpha \gamma }[\delta \mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{<(\beta )}(\omega ,E)G_{\gamma \alpha ,𝐤}^{\mathrm{adv}}(E)+\delta \mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{\mathrm{ret}(\beta )}(\omega ,E)G_{\gamma \alpha ,𝐤}^<(E)`$
$`+\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{<(\beta )}(E+\mathrm{}\omega )\delta G_{\gamma \alpha ,𝐤}^{\mathrm{adv}}(\omega ,E)+\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{\mathrm{ret}(\beta )}(E+\mathrm{}\omega )\delta G_{\gamma \alpha ,𝐤}^<(\omega ,E)],`$
where the factor 2 is again from the spin index summation. From the change in current density $`\delta J(\omega )`$, we obtain the polarization $`\delta P(\omega )=i\delta J(\omega )/\omega `$ induced in the material, and hence the complex susceptibility

$$\chi (\omega )=\frac{\delta P(\omega )}{ϵ_0(\omega )}=\frac{i}{ϵ_0}\frac{\delta J(\omega )}{\omega (\omega )},$$
(29)
and from Eq. (17) we obtain the gain coefficient $`g(\omega )`$.

## III Application to quantum cascade structures

We have applied the theory presented above to some example QCL structures reported in the literature. These structures were grown in the GaAs/Al<sub>x</sub>Ga<sub>1-x</sub>As material system, but this theory is also applicable to other material systems. The examples we consider are \[A\] the first GaAs/Al<sub>x</sub>Ga<sub>1-x</sub>As QCL structureSir98  reported by Sirtori et al., with $`x=0.33`$, and \[B\] a structurePag01  reported by Page et al. with a similar layer design to structure A but with a higher Al content in the barrier, $`x=0.45`$, resulting in higher barriers. This QCL structure operates in pulsed mode at room temperature. The lasing transition in structures A and B occurs between quantum well subbands. Other parameters for these structures are summarized in Table 2. In this section, we present and discuss results obtained for these structures.

### III.1 Nonequilibrium stationary state

In this subsection, we describe some of the information about the nonequilibrium stationary state of a QC structure under an applied bias that can be extracted from the Green’s functions. In Fig. 2, we show some examples of the diagonal elements of $`\mathrm{Im}[𝐆^{\mathrm{ret}}(E)]`$ for $`𝐤=0`$ plotted as a function of the energy parameter $`E`$. Figure 2(a) shows two examples of this function evaluated in the Wannier basis. These examples were evaluated for structure A at a bias of 0.2 V/period. Each curve represents a Wannier state in one period of the QC structure. The Wannier states are not energy eigenstates of the system, and if a Wannier state is expressed as a linear superposition of the energy eigenstates, the position of the peaks in, for instance, the curve labeled $`\nu =3`$ represent the energies of the eigenstates which comprise this superposition.

It is possible also to express the Green’s functions in the Wannier-Stark basis. From the diagonalization of the Hamiltonian $`\widehat{H}_o`$ expressed in the Wannier basis we obtain the transformation matrix between the Wannier and Wannier-Stark basis.foot.ws  Applying this transformation matrix to $`𝐆^{\mathrm{ret}}(E)`$ we obtain the example curves shown in Fig. 2(b). Each of these curves (diagonal elements of $`\mathrm{Im}[𝐆^{\mathrm{ret}}(E)]`$ for $`𝐤=0`$ in the Wannier-Stark basis) represents a Wannier-Stark state. In contrast to the curves shown in Fig. 2(a) with their complicated structure of peaks, each curve in the Wannier-Stark basis consists of a single, broadened peak. This indicates that at this bias, the Wannier-Stark states approximate closely the true energy eigenstates of the system. The positions of these peaks $`E_i`$ are shifted downwards by around 10 to 15 meV in comparison to the eigenenergies obtained from the diagonalization of $`\widehat{H}_o`$. This shift represents the renormalization of the Wannier-Stark energy levels due to the scattering processes described in $`\widehat{H}_{\mathrm{scatt}}^{}`$. The broadening of the peaks also originates from these scattering processes, and the width of the peaks gives a measure of the decay rate $`\mathrm{\Gamma }_i`$ of these levels due to scattering. In particular, we take $`\mathrm{\Gamma }_i=\mathrm{\Delta }E_{\mathrm{fwhm}}`$, where $`\mathrm{\Delta }E_{\mathrm{fwhm}}`$ is the full width at half maximum of the single peak in the function $`\mathrm{Im}[G_{ii,𝐤=0}^{\mathrm{ret}}(E)]`$, in the Wannier-Stark basis. In Sec. III.3, we will use this information to analyze the gain spectra obtained from the NGF theory.

The curves in Fig. 2 were evaluated for $`𝐤=0`$. If instead we integrate over the in-plane wave vector $`𝐤`$, we obtain the curves shown in Fig. 3(a). These staircaselike curves are the density of states (DOS) in two example Wannier subbands (the two Wannier states considered in Fig. 2(a) are the $`𝐤=0`$ states in these subbands). The onset of each step in the staircase corresponds to the position of each peak in the corresponding curve in Fig. 2(a). If we carried out the $`𝐤`$ integral in the Wannier-Stark basis, we would obtain a single step instead of a staircaselike structure, corresponding to the single peak (for a given curve) seen in Fig. 2(b). If we sum the DOS curves for all Wannier subbands within a period, we obtain Fig. 3(b), which is the total DOS per period of the structure. The result in Fig. 3(b) is basis independent if we sum over contributions from all basis states within the energy range shown.

We consider now the lesser correlation function $`𝐆^<(E)`$. At equal times, $`G_{\nu \nu ,𝐤}^<(t)=i\widehat{a}_{\nu 𝐤}^{}(t)\widehat{a}_{\nu 𝐤}(t)=i\widehat{n}_{\nu 𝐤}(t)`$, where $`\widehat{n}_{\nu 𝐤}(t)`$ is the occupation of a state k in subband $`\nu `$. If we move into the energy representation, and sum over the contributions from all states in all subbands in a period, we obtain the population $`n(E)=2_{\nu 𝐤}G_{\nu \nu ,𝐤}^<(E)`$ as a function of energy within the period. Figure 4(a) shows the normalised population per period per unit energy for several applied voltages. Dividing these curves by the total DOS curves for the corresponding bias gives the distribution function as a function of energy, i.e.,

$$f(E)=\frac{_{\nu 𝐤}G_{\nu \nu ,𝐤}^<(E)}{2_{\nu 𝐤}\mathrm{Im}[G_{\nu \nu ,𝐤}^{\mathrm{ret}}(E)]}.$$
(30)
The distribution function at 0.26 V/period in structure A is shown in Fig. 4(b). In this curve, the jagged points between 0 and 0.2 eV are LO-phonon replicas, and the onset of a small population inversion at around 0.2 eV is also seen.

### III.2 Transport

We consider now the transport properties of the QC structures described above. Figure 5 shows a plot of current density against V/period, calculated from Eqs. (11) – (16). The current density for structure B is reduced below structure A because of the higher barrier heights in B. The separate current contributions $`J_o`$ and $`J_{\mathrm{scatt}}+J_{\mathrm{MF}}`$ from Eqs. (13) – (16) are also shown. In both the structures, the main contribution is from $`J_{\mathrm{scatt}}`$ ($`J_{\mathrm{MF}}`$ is negligible in comparison with $`J_{\mathrm{scatt}}`$). This finding is in contrast to the behavior in simple superlattice structures where the dominant contribution to the current is from $`J_o`$. This difference could perhaps be explained by the far greater number of subbands within one period of a QC structure (in comparison to one period of a simple symmetric superlattice structure), in which the envelope functions $`\psi _\alpha (z)`$ are spatially displaced across the period. This opens up the possibility for many more scattering transitions which facilitates the transport of carriers across the structure through scattering processes. Comparing A and B, we see that $`J_o`$ is almost identical for both structures. The current $`J_{\mathrm{scatt}}`$ in B is around half that in A so that the reduction in total current in B compared to A arises from a decrease in the transport driven by the scattering Hamiltonian $`\widehat{H}_{\mathrm{scatt}}`$.

To compare the theoretical results with experimental data reported in the literature, we show in Fig. 6 the results for A and B plotted against similar scales as the experimental plots shown in Fig. 3 in Ref. Sir98, , and Fig. 4 in Ref. Pag01, . Comparing the results for structure A with the data given in Ref. Sir98, , we see the voltage trend in both curves agree well. In the theory, the voltage tends to be lower for a given current than in the experiment, but this may be accounted for by an additional series resistance in the cladding layers. For structure B, we see again that the voltage trends agree well with the data shown in Ref. Pag01, . In particular, when comparing the trends at different temperatures, we find that in both the theory and experiment the current is higher for a given voltage at the higher temperature. In both theory and experiment, the curves at both temperatures converge at a certain current density. Unlike the experiment, the two curves do not cross at the point of convergence in the theory, but diverge beyond this point. The theoretical curves reproduce the sharp break in the voltage trend at around 23 kA/cm<sup>2</sup> as seen in the experiment, and the theory shows that this break corresponds to the onset of a region of negative differential resistivity (NDR). NDR is not seen explicitly in the experimental data in Ref. Pag01,  because of an artefact of the measurement system.foot.pagendr

The break in the voltage trend is attributedPag01  by Page et al. to the breaking of the level alignment in the injector region and the active region causing an interruption in the current flow. To verify this, we show the wave functions and alignment of the injector level i and upper laser level u (in the Wannier-Stark basisfoot.ws ) at 77 K, for different applied voltages, in Fig. 7. Figures 7(a) – (c) correspond to the voltages marked a, b, and c in Fig. 6(b). Figure 6(b) shows that as the voltage increases from 10 V, the current-voltage characteristic passes through a peak at 13 V, and then beyond this point the current decreases to a minimum at around 16.5 V. Comparing this behavior with Fig. 7, we see that the change in voltage, and the behavior of the I-V characteristic, is accompanied by a shift in the position of level i from below to above level u. Level i passes through the resonance with level u which causes the peak in the I-V characteristic seen at 13 V. The population percentage in these levels is also shown in Fig. 7, and we see that away from the resonance position, $`50`$% of the population is in the injector level, with only $`10`$% in the upper laser level. Close to resonance, however, we see that the population is more evenly distributed between the two levels, with $`30`$% in the injector level, and 20% in the upper laser level. The wave functions in the near-resonance case (b) are less well-localized than in the off-resonance cases (a) and (c), and tend to spread across both the injector and active region. This increases the overlap of the two wave functions and facilitates the transfer of carriers between the injector and active region.

In the experimental data at 233 K, the break in the voltage trend is less pronounced, and this was attributedPag01  to the presence of electrons with higher energy leaking into continuum states, and maintaining the current flow. In the theoretical results, shown in Fig. 6(b), there is a small increase in current at 233 K, but it is difficult to make a quantitative comparison with the experimental data because of the uncertainty in the experimental data in the NDR region.foot.pagendr  Figure 8 shows the electron distribution functions $`f(E)`$ at 77 K and 233 K, in one period of the structure at 14.4 V, i.e., just beyond the peak in the I-V characteristic in Fig. 6(b). The thin horizontal line marks the conduction-band offset or barrier height. The peak in the distributions at around 0.2 eV corresponds to the population inversion in the upper laser level. The population inversion at 233 K is about a third less than that at 77 K. In the high-energy tails, the distribution is slightly larger at 233 K than at 77 K. The difference in the distribution functions is not large, however, and although this difference may result in more electron leakage into the continuum at 233 K, it is not clear at this point if this is a sufficiently large effect to elucidate the experimental results.

### III.3 Gain and absorption spectra

In this section, we apply the theory described in Sec. II.4 to structure A. We distinguish here between the gain coefficient $`g(\omega )`$ evaluated using $`\delta J(\omega )`$ from Eq. (25), and the gain coefficient $`g_o(\omega )`$ which is evaluated neglecting $`\delta J_{\mathrm{scatt}}(\omega )`$ in Eq. (25), and neglecting terms containing $`\delta 𝚺^{</\mathrm{ret}/\mathrm{adv}}`$ in $`\delta 𝐆^<(\omega ,E)`$ \[Eq. (24)\]. Hence, $`g_o(\omega )`$ does not depend on changes in the self-energies, and it is simpler to evaluate since its evaluation does not require a further self-consistent calculation.

In Fig. 9, we show the gain coefficient $`g(\omega )`$ calculated with the NGF theory for different applied voltages, between 0 to 0.2 V/period. With zero applied bias, there is a broad absorption ranging from around 120 to 140 meV. As the applied voltage increases, the absorption gradually decreases in this range. There is a corresponding slow increase in absorption between around 80 to 100 meV. At around 0.18 V/period, a positive gain begins to appear in the range 120 to 140 meV, and this gain increases further as the voltage is further increased. These results were obtained in the Wannier basis and as stated in Sec. II.1, there is a problem in interpreting the origin of the different features in the spectra because these cannot be related to transitions between Wannier states. Later in this section, we discuss a way round this problem by considering a transformation to Wannier-Stark states.

Before this we compare the results obtained for $`g(\omega )`$ from the more rigorous NGF theory, with the gain $`g_o(\omega )`$ which is calculated neglecting terms from $`\delta J_{\mathrm{scatt}}`$. This comparison is shown in Fig. 10. We see that the simpler theory gives rise to gain curves which are offset to higher values (i.e., absorption is reduced and gain is increased) in comparison to the more rigorous theory.

To determine the origin of the different features in the spectra, we consider a different approach for evaluating the gain coefficient. In this approach, we transform the Green’s functions $`𝐆^{</\mathrm{ret}}`$ obtained in the Wannier basis to the Wannier-Stark basis ($`𝐆_{\mathrm{WS}}^{</\mathrm{ret}}`$) as described earlier in Sec. III.1. From the diagonal elements of $`G_{\mathrm{WS},ii}^<`$, we obtain the populations $`n_i`$ of the Wannier-Stark levels, and from the positions and widths of the single peaks in $`\mathrm{Im}[G_{\mathrm{WS},ii}^{\mathrm{ret}}]`$, we obtain the Wannier-Stark level energies, $`E_i`$, for $`𝐤=0`$, and the corresponding decay rates, $`\mathrm{\Gamma }_i`$. We treat each pair of Wannier-Stark subbands in a simple two-band model, in which the optical response due to subbands $`i`$ and $`j`$ is given by,Hau93

$$\mathrm{Im}[\chi _{ij}(\omega )]=\frac{2\pi }{ϵ_0𝒱}\underset{𝐤}{}|d_{ij}|^2(f_{i,𝐤}f_{j,𝐤})\delta (\mathrm{}\omega +E_{i,𝐤}E_{j,𝐤}),$$
(31)
where $`d_{ij}=e𝑑z\psi _i^{}(z)z\psi _j(z)`$ is the dipole matrix element between the Wannier-Stark functions $`\psi _i(z)`$ and $`\psi _j(z)`$, and $`f_{i,𝐤}`$ and $`f_{j,𝐤}`$ are the nonequilibrium distribution functions in subbands $`i`$ and $`j`$. The factor 2 comes from summing over the spin index. Assuming parabolic bands and replacing the $`\delta `$ function with a Lorentzian to model the broadening of the levels we obtain

$$\mathrm{Im}[\chi (\omega )]=\frac{\pi }{ϵ_0d}\underset{\genfrac{}{}{0pt}{}{ij}{(E_j>E_i)}}{}|d_{ij}|^2(n_in_j)_{ij}(\omega ),$$
(32)
with the Lorentzian $`_{ij}(\omega )=(\mathrm{\Gamma }_{ij}/2\pi )/[(\mathrm{}\omega \mathrm{\Delta }E_{ji})^2+(\mathrm{\Gamma }_{ij}/2)^2]`$, where $`\mathrm{\Gamma }_{ij}=\mathrm{\Gamma }_i+\mathrm{\Gamma }_j`$, and $`\mathrm{\Delta }E_{ji}=E_jE_i`$ (these parameters are defined in Sec. III.1). $`n_i`$ ($`n_j`$) is the 2D carrier density in subband $`i`$ ($`j`$). In Eq. (32), the contributions from all transitions with $`E_j>E_i`$ are included. Using Eq. (32) with Eq. (17) gives us an estimate of the gain coefficient $`g_{WS}(\omega )`$ within the Wannier-Stark picture.

Disadvantages of this simpler approach, compared to the approach described in Sec. II.4, are, firstly, that the Wannier-Stark levels become increasingly delocalized at low bias and an increasingly large number of basis states must be used to obtain a good approximation to the energy eigenstates as the applied voltage decreases. Thus, this approach becomes impracticable to apply at low bias. In addition, the function $`\mathrm{Im}[G^{\mathrm{ret}}(E)]`$ is not always a simple Lorentzian as seen in Fig. 2(b) for each Wannier-Stark level, but may, for some levels, have additional peaks or shoulders indicating that these levels do not approximate well a single eigenstate of the system but are superpositions of the eigenstates. This behavior occurs mainly at low bias, but can also occur when two Wannier-Stark levels are very close in energy which gives a double-peaked structure to $`\mathrm{Im}[G^{\mathrm{ret}}(E)]`$, e.g., for the two levels in panel (b) in Fig. 7. In these cases, it is not possible to use this simple model to estimate the gain. Another drawback of this model is that by only making use of the diagonal elements of the Green’s functions, $`G_{\mathrm{WS},ii}^<`$ and $`G_{\mathrm{WS},ii}^{\mathrm{ret}}`$, we have discarded information on quantum effects that are contained in the offdiagonal elements.

In Fig. 11, we compare the gain coefficient $`g(\omega )`$ from the NGF theory with the gain coefficient $`g_{WS}(\omega )`$ from the simple model in the Wannier-Stark basis. Spectra are shown for a bias of 0.12 – 0.2 V/period. We find that the simpler model reproduces well the main features of the NGF theory. As in the previous figure, we see that for the simpler theory in Fig. 11, the gain curves are again offset to higher values compared to the more rigorous theory.

We consider now the origin of the different features in the spectra. This can be most simply understood if we look at the gain spectra obtained in the Wannier-Stark basis. For clarity, only these spectra are shown again in Fig. 12, and for a wider bias range, 0.1 to 0.22 V/period, than seen in Fig. 11. These spectra contain transitions between all possible pairs of levels with energy separation $`\mathrm{\Delta }E_{ji}`$ in the energy range shown on the $`x`$ axis. At the lowest bias shown, i.e., 0.1 V/period, we see that the structure is almost transparent with only a very small absorption for the energy range shown. As the applied voltage increases, absorption increases in the range 0.08 to 0.11 eV. Simultaneously, gain appears and increases in the range 0.11 to 0.16 eV. The gain peak shifts to higher energies as the voltage increases.

To explain these features, we consider the curve at 0.2 V/period in Fig. 12. In Fig. 13, we show the main Wannier-Stark energy levels involved in the transitions which give rise to the absorption and gain features seen at this bias. From an inspection of the transitions contributing to the sum in Eq. (32), we find as expected that the gain feature originates mainly from the transition (marked B in the figure) between the upper laser level (labeled 2) and the lower laser level (1). There is also some contribution to the gain from transitions between an injector level (2’) and the lower laser level (1), and between other injector levels and the upper or lower laser level (not shown). The absorption feature below 0.1 eV is due mainly to transitions (marked A in the figure) from the upper laser level (2) to levels in the continuum (3 and 3’), with additional contributions also from transitions between injector levels. The relative population in each level is also shown in the figure, and levels 2 and 2’ have a greater fraction of the population compared to levels 1, 3, and 3’. This supports the attribution of gain to transition B, and absorption to transition A. The largest fraction of the population $``$ 70% is in the lowest injector level (not shown in the figure). Level 0 contains around 2% of the population.

When we consider the relative populations in these levels for different applied voltages (in the range 0.1 – 0.22 V/period as shown in Fig. 12), we find that the population in the lower laser level remains at around 1% at all voltages. The population in the upper laser level, however, changes from around 0.2% at 0.1 V/period to over 10% at 0.22 V/period. Hence, there is a transition to a population inversion between the two levels as the bias increases. The population in level 0 ranges from around 20% at 0.1 V/period to around 2 to 3% above 0.16 V/period. The largest proportion of the population remains in the injector region, particularly in the lowest injector level which forms a reservoir of around 60 to 70% of the carriers at all applied voltages.

We also note here that the broad gain or absorption features seen in Fig. 12 are not due to transitions between a single pair of levels but contain contributions from many transitions. For example, in Fig. 14, we show the different Lorentzian contributions (thin lines) to the gain curve (thick line) at 0.22 V/period. Although the main contribution to the absorption feature (thin line labeled A) originates from the upper laser level – continuum transition, and the main contribution to the gain feature (thin line labeled B) arises from the upper – lower laser level transition, there is also a substantial contribution from other (energetically) neighboring transitions. The figure also shows that the inhomogeneous broadening due to the contributions from different transitions, and the broadening due to scattering, observed in the linewidth of each individual Lorentzian, are comparable in size. At higher applied bias the number of contributing transitions decreases, e.g, at 0.3 V/period, the contributions to the gain feature come mainly from only two transitions.

Finally, in the last set of results, we consider structure B. Experimental data showing the light output vs. current density is shown in Fig. 4 of Ref. Pag01, . At 77 K, the light output is seen to increase with increasing current density until it reaches a maximum at around 22 kAcm<sup>-2</sup>, corresponding to the peak in the I-V characteristic seen in Fig. 6(b). The light output goes to zero for higher current densities. A similar behavior is seen at 233 K but with a much reduced light output. Figure 15 shows the gain coefficient $`g_{WS}(\omega )`$ calculated at 77 K for structure B at different applied voltages around the peak of the I-V characteristic shown in Fig. 6(b). The gain feature increases as the applied voltage or current density increases, reaching a maximum around the peak of the I-V characteristic, and then decreases beyond this point. A similar behavior is also seen at 233 K except that the gain is reduced at the higher temperature. For comparison, a gain curve calculated at 233 K, at the peak current density 22 kAcm<sup>-2</sup>, is also shown in Fig. 15. The gain feature in this curve is much smaller than the gain feature in the corresponding curve (at 22 kAcm<sup>-2</sup>) at 77 K. Hence, the behavior of the calculated gain curves correlate well with the light output vs. current density curves measured in the experiment.

## IV Conclusion

We have applied an NGF theory to obtain a description of the nonequilibrium stationary state of QC structures under an applied bias. Using this information, we evaluate the current-voltage characteristic of example QC structures reported in the literature. The theoretical results are quantitatively close to experimental I-V data, and reproduce well the trends seen in the data. In addition, we determine two contributions to the current density. The first contribution $`J_o`$ is driven by $`\widehat{H}_o`$ which is the Hamiltonian for the superlattice potential with applied bias. The other contribution $`J_{\mathrm{scatt}}^{}`$ is driven by $`\widehat{H}_{\mathrm{scatt}}+\widehat{H}_{\mathrm{MF}}`$ which describes the scattering processes in the structure. We find that, in contrast to simple superlattice structures, $`J_{\mathrm{scatt}}^{}`$ is the main contribution to the current in the QC structures we consider.

In addition, we have extended the theory to determine the linear response of the nonequilibrium stationary state of these structures to a small applied optical perturbation. This enables us to evaluate the linear susceptibility and hence the gain or absorption spectra of the structure. We compare the spectra obtained using a more rigorous NGF theory in which the changes $`\delta 𝐆`$, $`\delta 𝚺`$, and $`\delta J_{\mathrm{scatt}}`$ due to the optical perturbation are considered, to simpler models in which (i) only $`\delta 𝐆`$ is considered, or (ii) by summing over transitions in a simple two-band model (summing over different pairs of bands) with Lorentzian broadened levels. We find that the simpler models result in spectra which are offset to higher values than in the more rigorous theory. We have also made a detailed analysis of the origin of the different gain and absorption features in the spectra, and of the redistribution of population within the Wannier-Stark levels as the applied voltage changes. The gain and absorption features correlate well with the distribution of population within these levels, which is determined from the NGF theory.

###### Acknowledgements.

The authors thank the Deutsche Forschungsgemeinschaft for financial support through Grant No. FOR394. Discussions with A. Amann, H. Page, J. Schlesner, and M. Woerner are gratefully acknowledged.

## Appendix A

The construction of the Wannier function basis used in the calculations is described in this appendix. As the first step, we solve the one-dimensional Schrödinger equation

$$\left[\frac{\mathrm{}^2}{2m_e(z)}\frac{d^2}{dz^2}+V(z)\right]\psi (z)=E\psi (z),$$
(33)
for the envelope functions $`\psi (z)`$. The spatially dependent superlattice potential $`V(z)`$ and the effective mass $`m_e(z)`$ are assumed constant in each semiconductor layer, i.e., $`V(z)=V_o`$ and $`m_e(z)=m_b`$ in the barriers, and $`V(z)=0`$ and $`m_e(z)=m_w`$ in the wells. Equation (33) can be solved with a transfer matrix method (for a textbook discussion see, e.g., Ref. Yu99, ). In this approach, the envelope function in a semiconductor layer $`j`$ is written as

$$\psi _j(z)=A_je^{ik_j(E)(zz_j)}+B_je^{ik_j(E)(zz_j)},$$
(34)
where $`z_j`$ labels the position of interface $`j`$, and $`k_j(E)=\sqrt{2m_j(EV_j)}/\mathrm{}`$ with $`m_j`$ and $`V_j`$ the mass and potential in that layer. Applying continuity conditions

$`\psi _j`$ $`=`$ $`\psi _{j+1},`$ (35)
$`{\displaystyle \frac{1}{m_j}}{\displaystyle \frac{d\psi _j}{dz}}`$ $`=`$ $`{\displaystyle \frac{1}{m_{j+1}}}{\displaystyle \frac{d\psi _{j+1}}{dz}},`$ (36)
at the interface $`j+1`$, gives

$$\left(\genfrac{}{}{0pt}{}{A_{j+1}}{B_{j+1}}\right)=_j(E)\left(\genfrac{}{}{0pt}{}{A_j}{B_j}\right),$$
(37)
with

$$_j(E)=\frac{1}{2}\left(\begin{array}{cc}\left(1+\frac{m_{j+1}k_j}{m_jk_{j+1}}\right)e^{ik_j(z_{j+1}z_j)}& \left(1\frac{m_{j+1}k_j}{m_jk_{j+1}}\right)e^{ik_j(z_{j+1}z_j)}\\ \left(1\frac{m_{j+1}k_j}{m_jk_{j+1}}\right)e^{ik_j(z_{j+1}z_j)}& \left(1+\frac{m_{j+1}k_j}{m_jk_{j+1}}\right)e^{ik_j(z_{j+1}z_j)}\end{array}\right).$$
(38)
If a single period $`d`$ of the structure contains $`M`$ semiconductor layers, the Bloch condition $`\psi _q(z+d)=e^{iqd}\psi _q(z)`$ implies

$$\left(\genfrac{}{}{0pt}{}{A_{M+1}}{B_{M+1}}\right)=\underset{j=1}{\overset{M}{}}_j(E)\left(\genfrac{}{}{0pt}{}{A_1}{B_1}\right)=e^{iqd}\left(\genfrac{}{}{0pt}{}{A_1}{B_1}\right).$$
(39)
For a given value of $`q`$, only certain values of $`E`$ allow the solution of Eq. (39), and this defines the miniband structure $`E^\nu (q)`$. For each $`q`$ value, we determine $`E^\nu (q)`$ numerically by looking for the zeroes of the determinant of the matrix product $`_{j=1}^M_j(E)`$. These zeroes were found by stepping through the energy $`E`$, and comparing the signs of the determinant for two consecutive values of $`E`$ (separated by $`\mathrm{\Delta }E`$). When these signs are opposite this sets the first coarse bracketing of the zero position. This position is then further refined by halving this interval and successive intervals up to a hundred times, while comparing the determinant signs of the interval endpoints at each iteration. Once $`E^\nu (q)`$ is determined, we can also obtain $`A^\nu (q)=A_1`$ and $`B^\nu (q)=B_1`$, and the Bloch functions $`\psi _q^\nu (z)`$, from Eqs. (34), (37), and (39). In the calculations reported here, $`E^\nu (q)`$ was evaluated for 500 $`q`$ points for each miniband, and $`\mathrm{\Delta }E`$ was set to 5 meV. The Bloch functions were evaluated on a position grid with 1500 points per period. In structure A (B), there are eight (nine) minibands below the conduction-band offset. For the calculations here, we included one miniband above the continuum, so the band index $`\nu `$ runs from 1 to 9 (10) for structure A (B).

In the next step, we construct the Wannier functions (associated with miniband $`\nu `$)

$$\psi _\nu ^W(znd)=\sqrt{\frac{d}{2\pi }}_{\pi /d}^{\pi /d}𝑑qe^{inqd}\psi _q^\nu (z)$$
(40)
from a superposition of the Bloch functions in miniband $`\nu `$. The Wannier functions constructed by this superposition are not unique, and can be very different depending on the phase of the Bloch functions which can be chosen arbitrarily for each value of $`q`$. Ideally, we would like the Wannier functions to be as spatially localized as possible to reduce the interaction matrix elements between Wannier functions in different periods. As a first step to construct Wannier functions with this property, we first fix the phase of the Bloch functions such that these functions are real at some arbitrary position $`x_s`$. Different values of $`x_s`$ can be tested, and the resulting Wannier functions checked to see if they are fairly well-localized, e.g., within one period of the structure. We then express the position operator $`\widehat{z}`$ in the basis made up of this initial set of Wannier functions. Finally, we diagonalize the resulting matrix representation of $`\widehat{z}`$, and the resulting eigenfunctions give us the required set of Wannier functions associated with a given miniband. This process is repeated for each miniband.

The Wannier states associated with a given miniband are degenerate, and their energy expectation values lie at the center of the miniband. Their wave functions are spatially displaced from each other, with a separation given by the period of the structure, i.e, there is one Wannier state in each period of the structure. To set up the matrix representation of $`\widehat{z}`$ for a given miniband, we use Wannier states in 11 neighboring periods. The large number of periods used in this construction is necessary to improve the numerical accuracy of the result.

In the transport and gain calculations, however, we consider only couplings and matrix elements between Wannier functions in the same period, and nearest-neighbor periods. Hence, the matrices representing the Green’s functions and self-energies are constructed with 27 basis states (for structure A) from three periods.

## Appendix B

This appendix describes in more detail the evaluation of the self-energies and scattering matrix elements in Eqs. (6) – (8). The self-energies contain summations over $`𝐤^{}`$ of the form $`{\displaystyle \underset{𝐤^{}}{}}|V_{\alpha \beta }(𝐤,𝐤^{})|^2\{G.F.\}`$ where $`V_{\alpha \beta }`$ is a generic matrix element representing interface roughness, impurity, or phonon scattering, and $`\{G.F.\}`$ represents a Green’s function $`G_{\beta \beta ,𝐤^{}}(E^{})`$, with $`E^{}=E`$ or $`E^{}=E\pm E_{\mathrm{phon}}`$, and including, for the case of phonon scattering, a phonon distribution factor. Taking the summation to the continuous limit leads to

$$\underset{𝐤^{}}{}|V_{\alpha \beta }(𝐤,𝐤^{})|^2\{G.F.\}=\frac{𝒜}{(2\pi )^2}_0^{\mathrm{}}dk^{}k^{}_0^{2\pi }d\theta |V_{\alpha \beta }(k,k^{},\theta )|^2\{G.F.\}.$$
(41)
We carry out the angle integration assuming the Green’s function term does not depend on the angle, and defining the angle-integrated quantity $`[f(k^{})]_\theta ={\displaystyle _0^{2\pi }}𝑑\theta f(k^{},\theta )`$ we obtain

$`{\displaystyle \underset{𝐤^{}}{}}|V_{\alpha \beta }(𝐤,𝐤^{})|^2\{G.F.\}`$ $`=`$ $`{\displaystyle \frac{𝒜}{(2\pi )^2}}{\displaystyle }dk^{}k^{}[|V_{\alpha \beta }(k,k^{})|^2]_\theta \{G.F.\}`$ (42)
$`=`$ $`{\displaystyle \frac{𝒜}{(2\pi )^2}}{\displaystyle \frac{m_e}{\mathrm{}^2}}{\displaystyle }dE_k^{}[|V_{\alpha \beta }(k,k^{})|^2]_\theta \{G.F.\}`$
$`=`$ $`{\displaystyle \frac{1}{2\pi }}{\displaystyle \frac{𝒜\rho _o}{2}}[|V_{\alpha \beta }(k_{\mathrm{typ}},k_{\mathrm{typ}}^{})|^2]_\theta {\displaystyle }dE_k^{}\{G.F.\}.`$
The integral over momentum $`k^{}`$ has been transformed to an integral over energy $`E_k^{}`$, and the density of states $`\rho _o=m_e/(\pi \mathrm{}^2)`$. In the last line of Eq. (42) we assume the matrix element $`V_{\alpha \beta }(k,k^{})`$ is slowly varying compared to the Green’s function term and can be taken out of the integral. It is evaluated at fixed momenta $`k_{\mathrm{typ}}`$ and $`k_{\mathrm{typ}}^{}`$, and the choice of these momenta is described below.

From a comparison of the factors outside the integral, in the last line of Eq. (42), with a scattering or transition probability rate derived from Fermi’s golden rule we can definefoot.fgrrate  a fictitious scattering rate (in energy units)

$$\gamma _{\alpha \beta }=\frac{𝒜\rho _o}{2}\left[|V_{\alpha \beta }(k_{\mathrm{typ}},k_{\mathrm{typ}}^{})|^2\right]_\theta .$$
(43)
With this definition we can rewrite Eqs. (41) and (42) as

$$\underset{𝐤^{}}{}|V_{\alpha \beta }(𝐤,𝐤^{})|^2\{G.F.\}=\frac{\gamma _{\alpha \beta }}{2\pi }dE_k^{}\{G.F.\}.$$
(44)
Estimates of $`\gamma `$ for different scattering processes are given in Table 1.

As stated after Eq. (8), and earlier in this appendix, we evaluate the scattering matrix elements $`V_{\alpha \beta }^{\mathrm{phon}/\mathrm{rough}/\mathrm{imp}}(𝐤,𝐤^{})`$ using fixed momenta $`𝐤_{\mathrm{typ}}`$ and $`𝐤_{}^{}{}_{\mathrm{typ}}{}^{}`$ (with the corresponding energies $`E_{\mathrm{typ}}`$ and $`E_{\mathrm{typ}}^{}`$) to accelerate the numerical computation. To fix $`E_{\mathrm{typ}}`$, we consider the energy dependence of the scattering matrix elements. For LO-phonon scattering, there is an energy threshold for phonon emission because of energy conservation and the fixed phonon energy $`E_{\mathrm{LO}}`$. The scattering matrix element is maximum at the emission threshold, and decreases monotonically with increasing energy above this threshold. To obtain an estimate of the scattering matrix element, which lies between the higher values near threshold, and the lower values far above threshold, we set $`E_{\mathrm{typ}}`$ one LO-phonon energy above the LO-phonon threshold. We then fix $`E_{\mathrm{typ}}^{}=E_{\mathrm{typ}}+\mathrm{\Delta }E_{\alpha \beta }E_{\mathrm{LO}}`$, where $`\mathrm{\Delta }E_{\alpha \beta }=[\widehat{H}_o]_{\alpha \alpha }[\widehat{H}_o]_{\beta \beta }`$ (see Fig. 16). To test the sensitivity of the results to the value chosen for $`E_{\mathrm{typ}}`$, we have carried out runs with other values of $`E_{\mathrm{typ}}`$, e.g., at threshold, $`\frac{1}{2}E_{\mathrm{LO}}`$ above threshold, and $`2E_{\mathrm{LO}}`$ above threshold. The calculated current density changes by at most $`9`$% (at some bias points with $`E_{\mathrm{typ}}`$ at threshold) but in most cases the change is around 5% or much less (1 – 2%). This tends to support the assumption that, for the LO-phonon process, the results are not very sensitive to the specific value of $`E_{\mathrm{typ}}`$. For impurity and interface roughness scattering, we set $`E_{\mathrm{typ}}=15`$ meV, and $`E_{\mathrm{typ}}^{}=E_{\mathrm{typ}}+\mathrm{\Delta }E_{\alpha \beta }`$. Test runs were also carried out for $`E_{\mathrm{typ}}=1`$, 7, and 30 meV. The results are more sensitive to the value of $`E_{\mathrm{typ}}`$. For $`E_{\mathrm{typ}}=7`$ and 30 meV, the calculated current density changed by at most 15%. For $`E_{\mathrm{typ}}=1`$ meV, near the bottom of the subband, the difference was much larger, ranging from 10 to 50%. The value $`E_{\mathrm{typ}}=15`$ meV was chosen as a value that lies near the centre of the distributions in each subband, to give an estimate of the average scattering matrix element.

With the momenta fixed, the scattering matrix elements can be taken outside the integrals in Eqs. (6) – (8), and the self-energies depend only on integrals of the form \[as shown in Eq. (42)\]

$$G_{\alpha \alpha }^{\mathrm{ret}}(E)=_0^{\mathrm{}}𝑑E_k^{}G_{\alpha \alpha ,𝐤^{}}(E),$$
(45)
and

$$\mathrm{Im}[G_{\alpha \alpha }^<(E)]=_0^{\mathrm{}}𝑑E_k^{}\mathrm{Im}[G_{\alpha \alpha ,𝐤^{}}^<(E)].$$
(46)
Note that the diagonal elements of $`𝐆^<`$ are pure imaginary. A problem arises in evaluating $`G_{\alpha \alpha }^{\mathrm{ret}}(E)`$ when a momentum-independent scattering matrix element is used because this leads to a divergence in the integral as $`E_k^{}\mathrm{}`$. To deal with this problem, a cutoff energy for the upper limit of the integral is used in the numerical integration. The following subsections give more detail concerning the evaluation of the matrix elements for the different scattering processes (see, also, Refs. Lak97, , Che93, , Rob93, ).

### B.1 Interface roughness and impurity scattering

To derive matrix elements for interface roughness scattering we consider an interface $`j`$ located at $`z=z_j`$ with thickness fluctuations $`\xi _j(𝐫)`$ of the order of one monolayer. We assume correlation functions for the fluctuations given by

$`\xi _j(𝐫)`$ $`=`$ $`0`$ (47)
$`\xi _j(𝐫)\xi _j^{}(𝐫^{})`$ $`=`$ $`\delta _{j,j^{}}\stackrel{~}{f}(|𝐫𝐫^{}|)\text{with}\stackrel{~}{f}(r)=\eta ^2\mathrm{exp}\left({\displaystyle \frac{r}{\lambda }}\right)`$ (48)
$`\eta `$ denotes the root-mean-square of the roughness height and $`\lambda `$ is a typical island size. Correlations between neighboring interfaces are neglected. The angle brackets \[as stated after Eq. (6)\] denote an average over different distributions of thickness fluctuations.

The Hamiltonian for interface roughness scattering is written as

$$\widehat{H}_{\mathrm{rough}}=\underset{\genfrac{}{}{0pt}{}{𝐤,𝐩}{m\mu ,n\nu }}{}[V_{m\mu ,n\nu }^{\mathrm{rough}}(𝐩)\widehat{a}_{m,\mu }^{}(𝐤+𝐩)\widehat{a}_{n,\nu }(𝐤)+\mathrm{H}.\mathrm{c}.],$$
(49)
with the matrix element

$$V_{m\mu ,n\nu }^{\mathrm{rough}}(𝐩)=\underset{j}{}\frac{1}{A}d^2re^{i𝐩𝐫}\xi _j(𝐫)\mathrm{\Delta }E_c\psi _\mu ^{}(z_jmd)\psi _\nu (z_jnd),$$
(50)
where $`\mathrm{\Delta }E_c`$ is the band offset.

Within the self-consistent Born approximation, the self-energy contribution from interface roughness scattering is written as

$$\mathrm{\Sigma }_{\alpha _1\alpha _2,𝐤}^{<,\mathrm{rough}}(E)=\underset{\beta \beta ^{},𝐤^{}}{}V_{\alpha _1\beta }^{\mathrm{rough}}(𝐤𝐤^{})V_{\beta ^{}\alpha _2}^{\mathrm{rough}}(𝐤^{}𝐤)G_{\beta \beta ^{},𝐤^{}}^{<,\mathrm{rough}}(E).$$
(51)
This equation is more general than Eq. (6) since it includes the offdiagonal contributions. Now we assume that the diagonal parts of $`G_{\beta \beta ^{},𝐤^{}}^{<,\mathrm{rough}}`$ dominate, and we keep only the terms $`\beta =\beta ^{}`$ in the summation. Then we obtain for the matrix element term

$`V_{\alpha _1\beta }^{\mathrm{rough}}(𝐩)V_{\beta \alpha _2}^{\mathrm{rough}}(𝐩)`$ $`=`$ $`V_{m_1\mu _1,n\nu }^{\mathrm{rough}}(𝐩)V_{n\nu ,m_2\mu _2}^{\mathrm{rough}}(𝐩)`$ (52)
$`=`$ $`{\displaystyle \frac{(\mathrm{\Delta }E_c)^2}{𝒜^2}}{\displaystyle \underset{j^{}}{}}{\displaystyle }d^2r_1e^{i𝐩𝐫_1}\xi _j^{}(𝐫_1)\psi _{\mu _1}^{}(z_j^{}m_1d)\psi _\nu (z_j^{}nd)`$
$`\times {\displaystyle \underset{j}{}}{\displaystyle }d^2r_2e^{i𝐩𝐫_2}\xi _j(𝐫_2)\psi _\nu ^{}(z_jnd)\psi _{\mu _2}(z_jm_2d)`$
$`=`$ $`{\displaystyle \frac{(\mathrm{\Delta }E_c)^2}{𝒜^2}}{\displaystyle d^2r_1d^2r_2e^{i𝐩(𝐫_1𝐫_2)}\stackrel{~}{f}(|𝐫_1𝐫_2|)}`$
$`\times {\displaystyle \underset{j}{}}\psi _{\mu _1}^{}(z_jm_1d)|\psi _\nu (z_jnd)|^2\psi _{\mu _2}(z_jm_2d)`$
$`=`$ $`{\displaystyle \frac{(\mathrm{\Delta }E_c)^2}{𝒜^2}}{\displaystyle d^2re^{i𝐩𝐫}\stackrel{~}{f}(|𝐫|)}`$
$`\times {\displaystyle \underset{j}{}}\psi _{\mu _1}^{}(z_jm_1d)|\psi _\nu (z_jnd)|^2\psi _{\mu _2}(z_jm_2d),`$
where we have expanded the general indices $`\alpha _1,\alpha _2,\beta `$ in terms of the period and Wannier level indices \[see after Eq. (4)\], and defined $`𝐩=𝐤^{}𝐤`$. Because of the orthogonality of the wave functions, and because the wave functions extend over many interfaces, the sum over $`j`$ tends to vanish for $`m_1m_2`$ and $`\mu _1\mu _2`$. This suggests replacing this term in the last line with $`F_{m_1n}^{\mu _1\nu }\delta _{m_1,m_2}\delta _{\mu _1,\mu _2}`$ where

$$F_h^{\mu _1\nu }=\underset{j}{}|\psi _{\mu _1}^{}(z_jhd)|^2|\psi _\nu (z_j)|^2.$$
(53)
Applying Eq. (48) we obtain

$$V_{m\mu ,n\nu }^{\mathrm{rough}}(𝐩)V_{n\nu ,m\mu }^{\mathrm{rough}}(𝐩)=\frac{(\mathrm{\Delta }E_c)^2\eta ^2}{\rho _o𝒜E_\lambda }\frac{1}{(1+E_p/E_\lambda )^{3/2}}F_{mn}^{\mu \nu }=(E_p),$$
(54)
with $`E_\lambda =\mathrm{}^2/2m\lambda ^2`$ and $`E_p=\mathrm{}^2p^2/2m=E_k+E_k^{}2\sqrt{E_kE_k^{}}\mathrm{cos}\theta `$, where $`\theta `$ is the angle between $`𝐤`$ and $`𝐤^{}`$. The subscript 1 from the indices $`m_1`$ and $`\mu _1`$ is dropped for simplicity. Substituting this latter result in Eq. (51) gives Eq. (6) which contains only the diagonal terms in the self-energy and the Green’s function. We now follow the procedure outlined in Eqs. (41) – (42) and take the summation over $`𝐤^{}`$ in Eq. (6) to the continuous limit. We observe that both the self-energy and the Green’s function depend only on $`E_k`$ and $`E_k^{}`$, but not on the angle $`\theta `$. Hence, as shown in Eq. (42), the angle integration over the matrix element can be carried out analytically, and we define the angle-integrated quantity (Ref. Gra80, , 2.575):

$`\left[|V_{m\mu ,n\nu }^{\mathrm{rough}}(𝐤𝐤^{})|^2\right]_\theta `$ $`=`$ $`{\displaystyle _0^{2\pi }}𝑑\theta (E_k+E_k^{}2\sqrt{E_kE_k^{}}\mathrm{cos}\theta )`$ (55)
$`=`$ $`{\displaystyle \frac{(\mathrm{\Delta }E_c)^2\eta ^2}{\rho _o𝒜E_\lambda }}F_{mn}^{\mu \nu }{\displaystyle \frac{4}{(ab)\sqrt{a+b}}}E\left(\sqrt{{\displaystyle \frac{2b}{a+b}}}\right),`$
with

$$a=1+\frac{E_k+E_k^{}}{E_\lambda }b=2\frac{\sqrt{E_kE_k^{}}}{E_\lambda }.$$
(56)
$`E(x)`$ is the complete elliptic integral of the second kind which is of order $`\pi /2=E(0)>E(x)>E(1)=1`$ (Ref. Abr70, , 17.3). Therefore we set $`E(x)\pi /2`$, and fixing $`E_k=E_{\mathrm{typ}}`$ and $`E_k^{}=E_{\mathrm{typ}}+\mathrm{\Delta }E_{\alpha \beta }`$ as described earlier, we define the fictitious interface roughness scattering rate \[using Eq. (43)\]

$$\gamma _{(mn),\mu \nu }^{\mathrm{rough}}=\frac{\pi (\mathrm{\Delta }E_c)^2\eta ^2\sqrt{E_\lambda }F_{mn}^{\mu \nu }}{\left[E_\lambda +\left(\sqrt{E_{\mathrm{typ}}}\sqrt{E_{\mathrm{typ}}+\mathrm{\Delta }E_{\alpha \beta }}\right)^2\right]\sqrt{E_\lambda +\left(\sqrt{E_{\mathrm{typ}}}+\sqrt{E_{\mathrm{typ}}+\mathrm{\Delta }E_{\alpha \beta }}\right)^2}}.$$
(57)
Impurity scattering is mediated through the (3D) Coulomb interaction $`V1/q^21/|𝐤𝐤^{}|^2`$. Neglecting the angle-dependent term gives $`|V|^21/(E_k+E_k^{})^2`$. Fixing $`E_k=E_{\mathrm{typ}}`$ and $`E_k^{}=E_{\mathrm{typ}}+\mathrm{\Delta }E_{\alpha \beta }`$ as before, and neglecting factors of 2, we obtain an order-of-magnitude estimate of the impurity scattering rate

$$\gamma _{\alpha \beta }^{\mathrm{imp}}=\gamma _{\mathrm{par}}^{\mathrm{imp}}\frac{\left[𝑑z|\psi _\alpha ^{}(z)\psi _\beta (z)|\right]^2}{(1+\mathrm{\Delta }E_{\alpha \beta }/E_{\mathrm{typ}})^2}.$$
(58)
The parameter $`\gamma _{\mathrm{par}}^{\mathrm{imp}}`$ (see Table 1) is estimated from calculations of impurity scattering in simple superlattice structures.Wac97d

### B.2 LO-phonon

To derive the electron–LO-phonon scattering matrix element we start from the interaction HamiltonianFer89

$$\widehat{H}_{\mathrm{lo}}=\underset{𝐐}{}i\alpha (𝐐)[e^{i𝐐𝐑}\widehat{b}_𝐐^{}e^{i𝐐𝐑}\widehat{b}_𝐐],$$
(59)
with

$$\alpha (𝐐)=\left[\frac{e^2}{2𝒱}\frac{\mathrm{}\omega _{\mathrm{lo}}}{Q^2}\left(\frac{1}{ϵ_{\mathrm{}}}\frac{1}{ϵ_s}\right)\right]^{\frac{1}{2}}.$$
(60)
$`𝐐=(𝐪,q_z)`$ is the 3D phonon wave vector with the 2D in-plane component $`𝐪`$ and the component $`q_z`$ in the growth direction. Similarly, the 3D position vector $`𝐑=(𝐫,z)`$. $`\widehat{b}_𝐐^{}`$ $`(\widehat{b}_𝐐)`$ is the phonon creation (annihilation) operator, $`\mathrm{}\omega _{\mathrm{lo}}=E_{\mathrm{lo}}`$ is the LO-phonon energy, and $`𝒱`$ is the crystal volume. $`ϵ_{\mathrm{}}`$ and $`ϵ_s`$ are the high-frequency and static absolute permittivities.

The matrix element

$$V_{\alpha \beta }^{\mathrm{lo}}=\alpha |\widehat{H}_{\mathrm{lo}}|\beta =\underset{𝐐}{}i\alpha (𝐐)\alpha |e^{i𝐐𝐫}\widehat{b}_𝐐^{}e^{i𝐐𝐫}\widehat{b}_𝐐|\beta ,$$
(61)
is written with initial and final states, $`|\beta `$ and $`|\alpha `$, given by

$`|\beta `$ $`=`$ $`|\mathrm{\Psi }_{𝐤,\beta }(𝐫,z)|n(\omega _{\mathrm{lo}}),`$
$`|\alpha `$ $`=`$ $`|\mathrm{\Psi }_{𝐤^{},\alpha }(𝐫,z)|n(\omega _{\mathrm{lo}})\pm 1.`$
The ket containing $`n(\omega _{\mathrm{lo}})`$ is a phonon number state. The upper (lower) sign in the ket $`|n(\omega _{\mathrm{lo}})\pm 1`$ corresponds to phonon emission (absorption). The ket containing $`\mathrm{\Psi }`$ gives the electron state, and the electron wave function $`\mathrm{\Psi }_{𝐤,\beta }(𝐫,z)=𝒜^{\frac{1}{2}}e^{i𝐤𝐫}\psi _\beta (z).`$ Evaluating Eq. (61) leads to

$$|V_{\alpha \beta }^{\mathrm{lo}}(𝐤,𝐤^{})|^2=|\alpha |\widehat{H}_{\mathrm{lo}}|\beta |^2=\underset{𝐐}{}\alpha ^2(𝐐)[n(\omega _{\mathrm{lo}})+\frac{1}{2}\pm \frac{1}{2}]|M_{\alpha \beta }(q_z)|^2\delta _{𝐤^{},𝐤𝐪},$$
(62)
with

$$M_{\alpha \beta }(q_z)=_0^{L_w}𝑑ze^{iq_zz}\psi _\alpha ^{}(z)\psi _\beta (z).$$

$`L_w`$ is the distance in the $`z`$-direction over which the wave functions $`\psi (k_z,z)`$ extend. Converting the sum over Q in Eq. (62) to an integral, and evaluating the in-plane component $`𝐪`$ with the help of the Kronecker delta gives

$$|V_{\alpha \beta }^{\mathrm{lo}}(𝐤,𝐤^{})|^2=[n(\omega _{\mathrm{lo}})+\frac{1}{2}\pm \frac{1}{2}]\frac{E_{\mathrm{lo}}e^2}{4\pi ϵ_p𝒜}_{\mathrm{}}^{\mathrm{}}𝑑q_z\frac{|M_{\alpha \beta }(q_z)|^2}{q_z^2+p^2},$$
(63)
where $`p=|𝐤𝐤^{}|^2=k^2+k_{}^{}{}_{}{}^{2}2kk^{}\mathrm{cos}\theta `$.

As described earlier in this appendix, the matrix element is evaluated within a self-energy integral, for instance of the form: $`\mathrm{\Sigma }(𝐤,E)=_𝐤^{}|V_{\alpha \beta }(𝐤,𝐤^{})|^2\{G.F\}`$. Following again the procedure shown in Eqs. (41) – (42), we take the summation over $`𝐤^{}`$ in the self-energy to the continuous limit, and we define the angle-integrated quantity

$`[|V_{\alpha \beta }^{\mathrm{lo}}(𝐤,𝐤^{})|^2]_\theta `$ $`=`$ $`C{\displaystyle _0^{2\pi }}𝑑\theta {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑q_z{\displaystyle \frac{|M_{\alpha \beta }(q_z)|^2}{q_z^2+k^2+k_{}^{}{}_{}{}^{2}2kk^{}\mathrm{cos}\theta }}`$ (64)
$`=`$ $`2\pi C{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑q_z{\displaystyle \frac{|M_{\alpha \beta }(q_z)|^2}{\sqrt{(q_z^2+k^2+k_{}^{}{}_{}{}^{2})^24k^2k_{}^{}{}_{}{}^{2}}}}`$
$`=`$ $`{\displaystyle \frac{C}{\rho _o}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑q_z{\displaystyle \frac{|M_{\alpha \beta }(q_z)|^2}{\sqrt{(\frac{\mathrm{}^2q_z^2}{2m}+E_k+E_k^{})^24E_kE_k^{}}}},`$
where $`C=E_{\mathrm{lo}}e^2[n(\omega _{\mathrm{lo}})+\frac{1}{2}\pm \frac{1}{2}]/(4\pi ϵ_p𝒜)`$.

Using this result in Eq. (43), and fixing the energies $`E_k=E_{\mathrm{typ}}`$ and $`E_k^{}=E_{\mathrm{typ}}^{}`$ as described earlier, we can define the scattering rate

$$\gamma _{\alpha \beta }^{\mathrm{lo}}=[n(\omega _{\mathrm{lo}})+\frac{1}{2}\pm \frac{1}{2}]\frac{E_{\mathrm{lo}}e^2}{8\pi ϵ_p}_{\mathrm{}}^{\mathrm{}}𝑑q_z\frac{|M_{\alpha \beta }(q_z)|^2}{\sqrt{(\frac{\mathrm{}^2q_z^2}{2m}+E_{\mathrm{typ}}+E_{\mathrm{typ}}^{})^24E_{\mathrm{typ}}E_{\mathrm{typ}}^{}}}.$$
(65)
### B.3 Acoustic phonon

The phonons are implemented as an artificial optical phonon with a phonon energy $`E_{\mathrm{ac}}`$ which should be smaller than $`k_BT`$ and which should not be commensurable with the optical phonon energy. The matrix element, or equivalently, the fictitious scattering rate $`\gamma _{\alpha \beta }^{\mathrm{ac}}`$, is set to

$$\gamma _{\alpha \beta }^{\mathrm{ac}}=\gamma _{\mathrm{par}}^{\mathrm{ac}}\left[𝑑z|\psi _\alpha ^{}(z)\psi _\beta (z)|\right]^2.$$
(66)
## Appendix C

A derivation of the current contribution $`J_{\mathrm{scatt}}`$ \[Eq. (16)\] is given here. We assume that the scattering Hamiltonian $`\widehat{H}_{\mathrm{scatt}}`$ has the form

$$\underset{\genfrac{}{}{0pt}{}{\alpha \beta }{𝐤,𝐤^{},s}}{}\widehat{O}_{\alpha 𝐤,\beta 𝐤^{}}(t)\widehat{a}_{\alpha 𝐤s}^{}(t)\widehat{a}_{\beta 𝐤^{}s}(t),$$
(67)
where $`\widehat{O}_{\alpha \beta }(t)`$ may be just a scalar time-independent matrix element, i.e. $`\widehat{O}_{\alpha 𝐤,\beta 𝐤^{}}(t)=V_{\alpha 𝐤,\beta 𝐤^{}}`$, as in interface roughness or impurity scattering. Alternatively, as in phonon scattering, it may be a time-dependent operator containing the phonon operators $`\widehat{b}(t)`$ and $`\widehat{b}^{}(t)`$. Then we find

$`J_{\mathrm{scatt}}(t)`$ $`=`$ $`{\displaystyle \frac{e}{𝒱}}{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{scatt}},\widehat{z}]`$ (68)
$`=`$ $`{\displaystyle \frac{e}{𝒱}}{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{\alpha \beta \gamma }{𝐤,𝐤^{},s}}{}}i\widehat{a}_{\alpha 𝐤s}^{}(t)[\widehat{O}_{\alpha 𝐤,\gamma 𝐤^{}}(t)z_{\gamma \beta }z_{\alpha \gamma }\widehat{O}_{\gamma 𝐤,\beta 𝐤^{}}(t)]\widehat{a}_{\beta 𝐤^{}s}(t).`$
In the following, this expression will be expressed in terms of Green’s functions and self-energies. For this purpose, we consider the contour-ordered Green’s function:

$$F_{\genfrac{}{}{0pt}{}{\alpha \beta \gamma }{𝐤,𝐤^{},s}}^c(\tau _1,\tau _2)=i\widehat{T}_c\{\widehat{a}_{\beta 𝐤^{}s}(\tau _1)\widehat{a}_{\alpha 𝐤s}^{}(\tau _2)\widehat{O}_{\alpha 𝐤,\gamma 𝐤^{}}(\tau _2)z_{\gamma \beta }z_{\alpha \gamma }\widehat{O}_{\gamma 𝐤,\beta 𝐤^{}}(\tau _1)\widehat{a}_{\beta 𝐤^{}s}(\tau _1)\widehat{a}_{\alpha 𝐤s}^{}(\tau _2)\},$$
(69)
where $`\widehat{T}_c`$ orders the times arguments $`\tau `$ on the complex contour from $`\tau =\mathrm{}+i0^+`$ to $`\tau =\mathrm{}i0^+`$ (see Ref. Hau96, ). By taking the lesser component $`F_{\genfrac{}{}{0pt}{}{\alpha \beta \gamma }{𝐤,𝐤^{},s}}^<(t,t)`$, one obtains the summand in Eq. (68).

Using the interaction picture with superscript $`D`$, $`F^c`$ is now evaluated according to the standard perturbation expansion of

$`F_{\genfrac{}{}{0pt}{}{\alpha \beta \gamma }{𝐤,𝐤^{},s}}^c(\tau _1,\tau _2)`$ $`=`$ $`i\widehat{T}_c\{\mathrm{exp}({\displaystyle \frac{i}{\mathrm{}}}{\displaystyle }d\tau \widehat{H}^D(\tau ))[\widehat{a}_{\beta 𝐤^{}s}^D(\tau _1)\widehat{a}_{\alpha 𝐤s}^D(\tau _2)\widehat{O}_{\alpha 𝐤,\gamma 𝐤^{}}(\tau _2)z_{\gamma \beta }`$ (70)
$`z_{\alpha \gamma }\widehat{O}_{\gamma 𝐤,\beta 𝐤^{}}(\tau _1)\widehat{a}_{\beta 𝐤^{}s}^D(\tau _1)\widehat{a}_{\alpha 𝐤s}^D(\tau _2)]\}.`$
The exponential term can be expanded as

$$\mathrm{exp}\left(\frac{i}{\mathrm{}}𝑑\tau \widehat{H}^D(\tau )\right)1\frac{i}{\mathrm{}}𝑑\tau \underset{\genfrac{}{}{0pt}{}{\delta ϵ}{𝐩,𝐪,s^{}}}{}\widehat{O}_{\delta 𝐩,ϵ𝐪}(\tau )\widehat{a}_{\delta 𝐩s^{}}^D(\tau )\widehat{a}_{ϵ𝐪s^{}}^D(\tau ).$$
(71)
Substituting this expansion into Eq. (70), and noting that terms containing only a single $`\widehat{O}`$ are zero after averaging, the lowest-order nonvanishing terms give:

$`F_{\genfrac{}{}{0pt}{}{\alpha \beta \gamma }{𝐤,𝐤^{},s}}^c(\tau _1,\tau _2){\displaystyle \frac{1}{\mathrm{}}}{\displaystyle \underset{\delta ϵ}{}}{\displaystyle }d\tau [G_{\beta 𝐤^{}s,\delta 𝐤^{}s}^{c0}(\tau _1,\tau )O_{\delta 𝐤^{},ϵ𝐤}(\tau )G_{ϵ𝐤s,\alpha 𝐤s}^{c0}(\tau ,\tau _2)O_{\alpha 𝐤,\gamma 𝐤^{}}(\tau _2)z_{\gamma \beta }`$
$`z_{\alpha \gamma }O_{\gamma 𝐤,\beta 𝐤^{}}(\tau _1)G_{\beta 𝐤^{}s,\delta 𝐤^{}s}^{c0}(\tau _1,\tau )O_{\delta 𝐤^{},ϵ𝐤}(\tau )G_{ϵ𝐤s,\alpha 𝐤s}^{c0}(\tau ,\tau _2)]`$ (72)
with the bare Green’s functions $`G_{\alpha 𝐤s,\beta 𝐤s}^{c0}(\tau _1,\tau )=i\widehat{T}_c\{\widehat{a}_{\alpha 𝐤s}^D(\tau _1)\widehat{a}_{\beta 𝐤s}^D(\tau )\}`$. Only Green’s functions diagonal in the momentum and spin indices are kept in Eq. (72). In order to be consistent with the perturbation expansion in the Green’s functions, further terms are taken into account, which replace the bare Green’s functions by the full Green’s functions. Then we find

$`{\displaystyle \underset{𝐤,𝐤^{}}{}}F_{\genfrac{}{}{0pt}{}{\alpha \beta \gamma }{𝐤,𝐤^{},s}}^c(\tau _1,\tau _2)`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{}}}{\displaystyle \underset{\delta ,𝐤}{}}{\displaystyle }d\tau [G_{\beta 𝐤s,\delta 𝐤s}^c(\tau _1,\tau )\mathrm{\Sigma }_{\delta \gamma 𝐤}^{c(\alpha ,r)}(\tau ,\tau _2)z_{\gamma \beta }`$ (73)
$`z_{\alpha \gamma }\mathrm{\Sigma }_{\gamma \delta 𝐤}^{c(\beta ,l)}(\tau _1,\tau )G_{\delta 𝐤s,\alpha 𝐤s}^c(\tau ,\tau _2)]`$
where

$$\mathrm{\Sigma }_{\delta \gamma ,𝐤}^{c(\alpha ,r)}(\tau ,\tau _2)=\underset{ϵ𝐤^{}}{}O_{\delta 𝐤,ϵ𝐤^{}}(\tau )G_{ϵ𝐤^{}s,\alpha 𝐤^{}s}^c(\tau ,\tau _2)O_{\alpha 𝐤^{},\gamma 𝐤}(\tau _2)$$
(74)
denotes the part of the self-energy which exhibits $`O_{\alpha \gamma }`$ on the right-hand side, and

$$\mathrm{\Sigma }_{\gamma \delta ,𝐤}^{c(\beta ,l)}(\tau _1,\tau )=\underset{ϵ𝐤^{}}{}O_{\gamma 𝐤,\beta 𝐤^{}}(\tau _1)G_{\beta 𝐤^{}s,ϵ𝐤^{}s}^c(\tau _1,\tau )O_{ϵ𝐤^{},\delta 𝐤}(\tau )$$
(75)
with $`O_{\gamma \beta }`$ on the left-hand side. In the second term in Eq. (73) we have exchanged the dummy indices $`\delta `$ and $`ϵ`$, and in the first term we have exchanged $`𝐤`$ and $`𝐤^{}`$.

For diagonal self-energies, which depend on diagonal Green’s functions (see Sec. II.2), one has $`\mathrm{\Sigma }_{\delta \gamma ,𝐤}^{c(\alpha ,r/l)}=\delta _{\delta \gamma }\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{(\alpha )}`$ where $`\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{(\alpha )}`$ is defined after Eq. (16). Using Langreth rulesLan76a  and changing to the energy representation leads to Eq. (16) which is given again here for reference:

$`J_{\mathrm{scatt}}`$ $`=`$ $`{\displaystyle \frac{2e}{\mathrm{}𝒱}}{\displaystyle \underset{\alpha \beta \gamma }{}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{dE}{2\pi }\underset{\mathrm{I}}{\underset{}{\left(G_{\beta \gamma ,𝐤}^<(E)\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{\mathrm{adv}(\alpha )}(E)+G_{\beta \gamma ,𝐤}^{\mathrm{ret}}(E)\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{<(\alpha )}(E)\right)z_{\gamma \beta }}}}`$ (76)
$`\underset{\mathrm{II}}{\underset{}{z_{\alpha \gamma }\left(\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{<(\beta )}(E)G_{\gamma \alpha ,𝐤}^{\mathrm{adv}}(E)+\mathrm{\Sigma }_{\gamma \gamma ,𝐤}^{\mathrm{ret}(\beta )}(E)G_{\gamma \alpha ,𝐤}^<(E)\right)}}.`$
We can obtain some insight into this expression for $`J_{\mathrm{scatt}}`$ if we consider only the diagonal terms of the above equation, i.e., we set $`\gamma =\beta `$ in part I, and $`\gamma =\alpha `$ in part II, to obtain

$`{\displaystyle 𝑑E\mathrm{I}}`$ $``$ $`{\displaystyle 𝑑E\left(\underset{\mathrm{Ia}}{\underset{}{G_{\beta \beta }^<\mathrm{\Sigma }_{\beta \beta }^{\mathrm{adv}(\alpha )}}}+\underset{\mathrm{Ib}}{\underset{}{G_{\beta \beta }^{\mathrm{ret}}\mathrm{\Sigma }_{\beta \beta }^{<(\alpha )}}}\right)z_{\beta \beta }},`$ (77)
$`{\displaystyle 𝑑E\mathrm{II}}`$ $``$ $`{\displaystyle 𝑑Ez_{\alpha \alpha }\left(\underset{\mathrm{IIa}}{\underset{}{\mathrm{\Sigma }_{\alpha \alpha }^{<(\beta )}G_{\alpha \alpha }^{\mathrm{adv}}}}+\underset{\mathrm{IIb}}{\underset{}{\mathrm{\Sigma }_{\alpha \alpha }^{\mathrm{ret}(\beta )}G_{\alpha \alpha }^<}}\right)},`$ (78)
where for brevity we neglect the index $`𝐤`$, and $`\mathrm{\Sigma }`$ and $`G`$ are understood to be functions of $`E`$. We observe that each term in the integrals above is a product of a Green’s function and a self-energy, i.e., of the form $`G\mathrm{\Sigma }`$ or $`\mathrm{\Sigma }G`$. These terms can be interpreted as scattering rates. In particular, Eq. (77) contains information about scattering rates into and out of the state $`\beta `$, and Eq. (78) describes scattering into and out of the state $`\alpha `$. To be more specific, we interpret the terms containing $`\mathrm{\Sigma }^{\mathrm{ret}}`$ or $`\mathrm{\Sigma }^{\mathrm{adv}}`$ as a scattering-out rate, e.g, the term $`\mathrm{Ia}`$, $`𝑑EG_{\beta \beta }^<\mathrm{\Sigma }_{\beta \beta }^{\mathrm{adv}(\alpha )}`$, can be interpreted as a rate $`\mathrm{\Gamma }_{\beta \alpha }^{\mathrm{out}}`$ for scattering out of state $`\beta `$ into state $`\alpha `$. Similarly, term $`\mathrm{IIb}`$, $`𝑑E\mathrm{\Sigma }_{\alpha \alpha }^{\mathrm{ret}(\beta )}G_{\alpha \alpha }^<`$, is interpreted as a rate $`\mathrm{\Gamma }_{\alpha \beta }^{\mathrm{out}}`$ for scattering out of state $`\alpha `$ into state $`\beta `$. On the other hand, the terms containing $`\mathrm{\Sigma }^<`$ are interpreted as scattering-in rates. Thus, term $`\mathrm{Ib}`$ describes the rate $`\mathrm{\Gamma }_{\alpha \beta }^{\mathrm{in}}`$ for scattering into state $`\beta `$ from $`\alpha `$, and term $`\mathrm{IIa}`$ describes the rate $`\mathrm{\Gamma }_{\beta \alpha }^{\mathrm{in}}`$ for scattering into state $`\alpha `$ from $`\beta `$. We note here that $`\mathrm{\Gamma }_{\beta \alpha }^{\mathrm{out}}=\mathrm{\Gamma }_{\beta \alpha }^{\mathrm{in}}<0`$ and $`\mathrm{\Gamma }_{\alpha \beta }^{\mathrm{in}}=\mathrm{\Gamma }_{\alpha \beta }^{\mathrm{out}}>0`$. Combining Eqs. (77) and (78) gives

$`{\displaystyle 𝑑EI}+II`$ $`=`$ $`\left(\mathrm{\Gamma }_{\beta \alpha }^{\mathrm{out}}+\mathrm{\Gamma }_{\alpha \beta }^{\mathrm{in}}\right)z_{\beta \beta }z_{\alpha \alpha }\left(\mathrm{\Gamma }_{\beta \alpha }^{\mathrm{in}}+\mathrm{\Gamma }_{\alpha \beta }^{\mathrm{out}}\right)`$ (79)
$`=`$ $`\left(\mathrm{\Gamma }_{\beta \alpha }^{\mathrm{out}/\mathrm{in}}+\mathrm{\Gamma }_{\alpha \beta }^{\mathrm{in}/\mathrm{out}}\right)(z_{\beta \beta }z_{\alpha \alpha }).`$
Thus, this expression is the product of the distance $`|z_{\beta \beta }z_{\alpha \alpha }|`$ with the net transfer rate ($`\mathrm{\Gamma }_{\beta \alpha }^{\mathrm{out}/\mathrm{in}}+\mathrm{\Gamma }_{\alpha \beta }^{\mathrm{in}/\mathrm{out}}`$) between state $`\beta `$ and $`\alpha `$, and we can interpret this as a velocity or charge transfer rate from e.g., $`z_{\beta \beta }`$ to $`z_{\alpha \alpha }`$, i.e., a current flow from $`z_{\beta \beta }`$ to $`z_{\alpha \alpha }`$. (The direction of current flow depends on the net transfer rate.)


text:
# CPT-2002/P.4416x CERN-TH/2002-350x BUHEP-02-41 𝐵_𝐾 from quenched overlap QCDWork supported in part by EU HPP contract HPRN-CT-2000-00145 and US DOE grant DE-FG02-91ER40676. We thank Boston University for use of its supercomputer facilities and the Gauge Connection for its gauge configurations.

## 1 Introduction

$`K^0\overline{K}^0`$ mixing induces indirect CP violation in $`K\pi \pi `$ decays and is governed, in the standard model, by the matrix element

$$\overline{K}^0|𝒪_{\mathrm{\Delta }S=2}(\mu )|K^0=\frac{16}{3}M_K^2F_K^2\times B_K(\mu ),$$
(1)
where $`𝒪_{\mathrm{\Delta }S=2}`$ is the four-quark operator given below (Eq. (2)) and where $`B_K`$ parametrizes deviation from the vacuum saturation approximation.

The benefit of using overlap fermions  to compute this matrix element is their exact $`SU(N_f)_L\times SU(N_f)_R`$ flavor-chiral symmetry at finite lattice spacing . This implies that the $`\mathrm{\Delta }S=2`$ operator is simple to construct on the lattice (unlike with staggered fermions) and does not mix with chirally dominant operators  (as it does with Wilson fermions). Overlap fermions further guarantee full $`𝒪\left(a\right)`$-improvement. On a different note, this exploratory computation tests the feasibility of calculating matrix elements of four-quark operators with overlap fermions.

## 2 Simulation details

We use Neuberger’s overlap action

$$=a^4\underset{x}{}\overline{q}\left[\left(1\frac{1}{2\rho }am_q\right)D+m_q\right]q,$$

with

$$D=\frac{\rho }{a}\left(1+X/\sqrt{X^{}X}\right)\mathrm{and}X=D_W\rho /a,$$

where $`D_W`$ is the massless Wilson-Dirac operator.

The renormalized $`\mathrm{\Delta }S=2`$ operator of Eq. (1) is related to the bare lattice operator through

$$𝒪_{\mathrm{\Delta }S=2}(\mu )=Z_{\mathrm{\Delta }S=2}(a\mu ,g_0)𝒪_{\mathrm{\Delta }S=2}^{bare}(g_0),$$
(2)
where the lattice operator with correct chiral properties is

$$𝒪_{\mathrm{\Delta }S=2}^{bare}=[\overline{s}\gamma _\mu (1\gamma _5)\widehat{d}][\overline{s}\gamma _\mu (1\gamma _5)\widehat{d}],$$

with $`\widehat{d}=\left(1aD/2\rho \right)d`$.

The computation is performed with the Wilson gauge action in the quenched approximation at $`\beta =6.0`$ ($`a^1(r_0)2.1\text{GeV}`$) on a $`16^3\times 32`$ lattice with a statistics of 80 configurations and with $`\rho =1.4`$. The pseudoscalar mesons simulated are composed of degenerate quarks with masses $`m_q=0.04,0.055,0.07,0.085,0.1`$ corresponding to a range from $`m_s/21.3m_s`$. Details regarding the implementation and the inversion of the overlap operator can be found in .

The three-point function used to extract the $`\mathrm{\Delta }S=2`$ matrix element and its asymptotic time behavior are

$$\underset{\mathrm{𝐱𝐲}}{}J_0^L(x)𝒪_{\mathrm{\Delta }S=2}^{bare}(0)J_0^L(y)\stackrel{Ty_0>x_01}{}$$
(3)
$$\frac{|K|J_0^L|0|^2}{(2M_K)^2}\overline{K}|𝒪_{\mathrm{\Delta }S=2}^{bare}|K\{\mathrm{e}^{M_K(y_0Tx_0)}+R_K$$

$$\times (\mathrm{e}^{M_K(x_0+y_0)\delta x_0}+\mathrm{e}^{M_K(x_0+y_02T)+\delta (y_0T)})\},$$

with $`J_\mu ^L=\overline{d}\gamma _\mu (1\gamma _5)\widehat{s}`$. We allow for a time-reversed contribution because our lattice is rather short in the time direction and this contribution can be as large as 3-4% for our lightest quark mass at $`x_0=Ty_0=8a`$, a point which is included in our fits. In Eq. (3) $`\delta `$ measures the finite-volume shift in the energy of the two kaon state. This shift varies from about 2% to 5% in going from our heaviest to our lightest mass and is consistent with zero within roughly three to one standard deviations. $`R_K`$ is clearly intimately related to $`0|𝒪_{\mathrm{\Delta }S=2}^{bare}|KK`$, but finite volume effects can be significant and we choose not to pursue the study of this time-reversed matrix element. The fit that we retain for determining the forward matrix element is to times $`x_0/a`$ and $`y_0/a`$ in the intervals $`68`$ and $`2426`$ with $`y_0>x_0`$. To obtain $`|K|J_0^L|0|`$ and $`M_K`$ we compute two-point functions of left-handed currents which we symmetrize in time and fit in the time range $`515`$.

The motivation for using left-handed currents is that they eliminate spurious zero-mode contributions which appear in finite volume and in the quenched approximation because the overlap operator satisfies an index theorem , as a good Dirac operator should. In our $`J_\mu ^L`$ correlators, topological zero modes, $`\varphi _0(x)`$, contribute as

$$\gamma _\mu (1\gamma _5)\frac{\varphi _0(x)\varphi _0^{}(y)}{m_q}\gamma _\mu (1\gamma _5),$$

which vanishes because of the chirality of these modes. The use of left-handed currents is proposed in  where it is suggested that the weak interactions of Goldstone bosons be studied in the $`ϵ`$-regime of Gasser and Leutwyler.

## 3 Chiral behavior of $`B_K`$

We obtain $`B_K^{bare}`$ from

$$B_K^{bare}=\frac{3}{8}\frac{\overline{K}|𝒪_{\mathrm{\Delta }S=2}^{bare}|K}{\overline{K}|\overline{J}_0^L|00|\overline{J}_0^L|K},$$

where $`\overline{J}_0^L`$ is the conjugate of $`J_0^L`$. Quenched chiral perturbation theory (Q$`\chi `$PT) predicts

$$B_K=B\left[16\left(\frac{M}{4\pi F}\right)^2\mathrm{ln}\frac{M}{\mathrm{\Lambda }_B}+b\left(\frac{M}{4\pi F}\right)^4\right],$$
(4)
where $`B`$, $`F`$ and $`M`$ are the leading order values of the $`B`$-parameter, leptonic decay constant and pseudoscalar mass, respectively. $`\mathrm{\Lambda }_B`$ is a scale which corresponds to $`𝒪\left(p^4\right)`$ contact terms in the effective theory and the term proportional to $`b`$ is added to parametrize likely higher order contributions. We fit this functional form to our results for $`B_K^{bare}`$, with $`M=M_K`$, $`F=F_\pi ^{phys}`$ or $`F_K^{phys}`$ and $`b0`$ or not. The fit parameters are $`B^{bare}`$, $`\mathrm{\Lambda }_B`$ and $`b`$ when it is not fixed to zero. 111Note that the definition of $`b`$ changes with $`F`$. The results of these fits are shown in Fig. 1. We find that the mass-dependence of $`B_K`$ is consistent with expectations from $`\chi PT`$, though results at lower masses and with better statistics would be necessary to confirm the coefficient of the logarithm. Until such results become available, we choose not to quote a number in the chiral limit. At the physical point $`M_K=M_K^{phys}`$, however, where the value of the $`B`$-parameter obtained is insensitive to the choice of functional form, these considerations are not important. 222The result for $`B_K`$ given below is obtained from the chiral fit with $`F=F_\pi ^{phys}`$ and $`b0`$.

## 4 Non-perturbative renormalization

We perform all renormalizations non-perturbatively in the RI/MOM scheme la . Thus we fix gluon configurations to Landau gauge and compute numerically approriate, amputated forward quark Green’s functions with legs of momentum $`p\sqrt{p^2}`$. Then we use the fact that with overlap fermions the renormalization constants $`Z_{J^L}=Z_V`$ and $`Z_{\mathrm{\Delta }S=2}=Z_{VV+AA}`$ to define

$$Z_{B_K}^{\mathrm{RI}}(ap,g_0)=\frac{\mathrm{\Gamma }_V(pa)^2}{\mathrm{\Gamma }_{VV+AA}(pa)},$$

where $`\mathrm{\Gamma }_𝒪`$ is the value of the non-perturbative, amputated Green’s function of operator $`𝒪`$ projected onto the spin-color structure of $`𝒪`$.

We find that in a range from approximatively 1.5 to 2 GeV the $`p`$-dependence of $`Z_{B_K}^{\mathrm{RI}}`$ matches 2-loop running obtained by combining the results of  and , with $`N_f=0`$ and $`\alpha _s`$ taken from . This is shown in Fig. 2 where we plot, vs $`p`$, $`Z_{B_K}^{\mathrm{RGI}}`$ which we define as the ratio of $`Z_{B_K}^{\mathrm{RI}}`$ to the 2-loop running expression. Also shown in Fig. 2 are the non-perturbative $`Z_{B_K}^{\mathrm{RI}}(p)`$ and its 1-loop counterpart $`Z_{B_K}^{\mathrm{RI},\mathrm{PT}}(p)`$, obtained by expanding in $`\alpha _s`$ the ratio of $`Z_{\mathrm{\Delta }S=2}^{\mathrm{RI},\mathrm{PT}}(p)`$ to $`(Z_A^{\mathrm{RI},\mathrm{PT}})^2`$. The latter are taken from  and are properly matched. The value of $`Z_{B_K}^{\mathrm{RGI}}`$ that we use is the result of fitting our data to a constant, in the momentum range shown in the plot. The renormalization constant in the $`\overline{\mathrm{MS}}\mathrm{NDR}`$ scheme is then computed through multiplying this value by the 2-loop running expression of  with $`N_f`$ and $`\alpha _s`$ chosen as above.

## 5 Conclusion

Though costly, it is entirely feasible to compute weak matrix elements with overlap fermions, at least in the quenched approximation. The advantages of this approach are the chiral symmetry and the $`𝒪\left(a\right)`$-improvement.

Our preliminary results are $`Z_{B_K}^{\mathrm{RGI}}=1.29(1)`$ for the renormalization constant and, for the $`B`$-parameter:

$$B_K^{\mathrm{RGI}}=0.85(10),B_K^{\mathrm{NDR}}(2\text{GeV})=0.61(7),$$

where the errors do not yet include an estimate of systematic uncertainties. For the $`B`$-parameter, agreement with continuum-limit world averages (see e.g. ) based on the staggered result of  is excellent. Agreement with the results of , obtained with a different formulation of overlap fermions, perturbative matching and larger quark masses is good. Finally, our result is respectively a half and one standard deviations above the results of  and , both obtained using domain-wall fermions with a finite fifth dimension.