# Splitting of Dirac cones in HgTe quantum wells: Effects of crystallographic
orientation, interface-, bulk-, and structure-inversion asymmetry

M. V. Durnev  G. V. Budkin  S. A. Tarasenko  Ioffe Institute, 194021 St.
Petersburg, Russia

###### Abstract

We develop a microscopic theory of the fine structure of Dirac states in
$(0lh)$-grown HgTe/CdHgTe quantum wells (QWs), where $l$ and $h$ are the
Miller indices. It is shown that bulk, interface, and structure inversion
asymmetry causes the anticrossing of levels even at zero in-plane wave vector
and lifts the Dirac state degeneracy. In the QWs of critical thickness, the
two-fold degenerate Dirac cone gets split into non-degenerate Weyl cones. The
splitting and the Weyl point positions dramatically depend on the QW
crystallographic orientation. We calculate the splitting parameters related to
bulk, interface, and structure inversion asymmetry and derive the effective
Hamiltonian of the Dirac states. Further, we obtain an analytical expression
for the energy spectrum and discuss the spectrum for (001)-, (013)- and
(011)-grown QWs.

##  I Introduction

Heterostructures containing band-inverted compound HgTe are in the focus of
modern research in solid state physics. Depending on heterostructure design,
particularly the thickness of HgTe layer, they host a variety of phases
including the phases of three-dimensional and two-dimensional (2D) topological
insulators, 2D gapless semiconductor, 2D semimetal, etc Qi2011 ; Kvon2020 . Of
special interest is the 2D gapless phase with linearly-dispersion Dirac
fermions that is realized in HgTe quantum wells (QWs) of critical thickness,
i.e., at the point of trivial insulator – topological insulator transition
Buttner2011 ; Tarasenko2015 .

In the model of centrosymmetric heterostructure, the Dirac cones in
HgTe/CdHgTe QWs are 2-fold degenerate yielding the 4-fold degenerate Dirac
point at $\bm{k}=0$ in the QW of critical thickness Buttner2011 . Here,
$\bm{k}$ is the in-plane wave vector. Bulk inversion asymmetry (BIA) related
to the lack of an inversion center in host zinc-blend crystal, interface
inversion asymmetry (IIA) related to anisotropy of chemical bonds at
interfaces, and structure inversion asymmetry (SIA) lift the Dirac state
degeneracy Dai2008 ; Konig2008 ; Winkler2012 ; Weithofer2013 ; Tarasenko2015 ;
Orlita2011 ; Zholudev2012 ; Olbrich2013 ; Minkov2016 ; Durnev2016 . This
splitting is contributed by canonical $\bm{k}$-linear Rashba Rashba1960 ;
Vasko1979 ; Bychkov1984  and Dresselhaus Dresselhaus1955 ; Dyakonov1986 ;
Pikus1988 ; Rashba1988  terms as well as the term lifting the 4-fold
degeneracy at $\bm{k}=0$ Dai2008 ; Konig2008 ; Winkler2012 ; Weithofer2013 ;
Tarasenko2015 . The anticrossing gap at $\bm{k}=0$ was found to be quite large
in (001)-grown QWs and originates mostly from light-hole–heavy-hole mixing at
the QW interfaces Tarasenko2015 ; Minkov2016 .

Many experiments, however, are being carried out on HgTe/CdHgTe structures
grown along low-symmetry crystallographic directions, such as [013] and [012],
see Refs. Zholudev2012, ; Olbrich2013, ; Minkov2016, ; Dantscher2015, ;
Dantscher2017, ; Minkov2017, . The choice of crystallographic orientations is
dictated by technology: MBE growth of HgTe and CdHgTe layers on low-symmetry
(lattice-mismatch) GaAs surface enables one to obtain high-quality structures
Dvoretsky2020 . This motivates theoretical studies of low-symmetry QWs
Minkov2017 ; Raichev2012 ; Budkin2022 .

Here, we develop a microscopic theory of the fine structure of Dirac states in
HgTe/CdHgTe QWs taking account IIA, BIA, and SIA coupling. We show that the
energy spectrum in the QW of the critical thickness dramatically depends on
the QW crystallographic orientation and calculate the splitting parameters. We
explore the class of $(0lh)$-grown QWs, where $l$ and $h$ are the Miller
indices, and study how the fine structure evolves from (001)- to (013)-, and
(011)-grown QWs.

##  II Fine structure of Dirac states

The Dirac states in HgTe/CdHgTe QWs of critical and close-to-critical
thickness are formed from the electron-like $\left|E1,\pm 1/2\right\rangle$
and heavy-hole $\left|H1,\pm 3/2\right\rangle$ states Bernevig2006 ;
Gerchikov1989 . The corresponding basis functions at $\bm{k}=0$ have the form

$\displaystyle\left|E1,+1/2\right\rangle=f_{1}(z)\left|\Gamma_{6},+1/2\right\rangle+f_{4}(z)\left|\Gamma_{8},+1/2\right\rangle\,,$
$\displaystyle\left|H1,+3/2\right\rangle=f_{3}(z)\left|\Gamma_{8},+3/2\right\rangle\,,$
$\displaystyle\left|E1,-1/2\right\rangle=f_{1}(z)\left|\Gamma_{6},-1/2\right\rangle+f_{4}(z)\left|\Gamma_{8},-1/2\right\rangle\,,$
$\displaystyle\left|H1,-3/2\right\rangle=f_{3}(z)\left|\Gamma_{8},-3/2\right\rangle\,,$
(1)

where $\bm{k}=(k_{x},k_{y})$ is the in-plane wave vector, $f_{j}(z)$
($j=1,3,4$) are the envelope functions, $z$ is the growth direction,
$\left|\Gamma_{6},m\right\rangle$ ($m=\pm 1/2$) and
$\left|\Gamma_{8},m\right\rangle$ ($m=\pm 1/2,\pm 3/2$) are the Bloch
amplitudes of the $\Gamma_{6}$ and $\Gamma_{8}$ bands, respectively, in the
Brillouin zone center. We consider $(0lh)$-oriented QWs and use the QW-related
coordinate frame $x\parallel[100]$, $y\parallel[0h\bar{l}]$, and
$z\parallel[0lh]$.

The effective $4\times 4$ Hamiltonian, which describes the coupling of the
basis states (II) and formation of the Dirac-like spectrum, can be derived in
the $\bm{k}$$\cdot$$\bm{p}$ theory, see Sec. III. Taking into account bulk,
structure, and interface inversion asymmetry in $(0lh)$-grown QWs, one can
present the effective Hamiltonian in the form

$H=H_{0}+H_{\rm IIA/BIA}+H_{\rm SIA}\,,$ (2)

where

$H_{0}=\left(\begin{array}[]{cccc}\delta&\mathrm{i}Ak_{+}&0&0\\\
-\mathrm{i}Ak_{-}&-\delta&&0\\\ 0&0&\delta&-\mathrm{i}Ak_{-}\\\
0&0&\mathrm{i}Ak_{+}&-\delta\end{array}\right)$ (3)

is the $\bm{k}$-linear Bernevig-Hughes-Zhang Hamiltonian (2D Dirac
Hamiltonian) Bernevig2006 , $A$ is a parameter determining the velocity of
Dirac fermions, $\delta$ is the energy distance between the $E1$ and $H1$
subbands in the absence of mixing.

Interface inversion asymmetry related to anisotropy of chemical bonds at the
interfaces and bulk inversion asymmetry related to the lack of inversion
center in host crystal lead to a mixing of the basis states. This mixing at
$\bm{k}=0$ is described by the Hamiltonian

$H_{\rm IIA/BIA}=\left(\begin{array}[]{cccc}0&-\eta\sin
2\theta&0&\mathrm{i}\gamma\cos 2\theta\\\ -\eta\sin
2\theta&0&\mathrm{i}\gamma\cos 2\theta&0\\\ 0&-\mathrm{i}\gamma\cos
2\theta&0&\eta\sin 2\theta\\\ -\mathrm{i}\gamma\cos 2\theta&0&\eta\sin
2\theta&0\end{array}\right),$ (4)

where $\eta$ and $\gamma$ are mixing parameters, $\theta={\rm arctan}(l/h)$ is
the angle between the QW growth direction $[0lh]$ and the $[001]$ axis. The
angles $\theta=0$, ${\rm arctan}(1/3)\approx 0.321$, and $\pi/4$ correspond to
$[001]$, $[013]$, and $[011]$ growth directions, respectively.

Structure inversion asymmetry in $(0lh)$-oriented QWs grown from cubic
materials also mixes the basis states at $\bm{k}=0$, which is described by the
Hamiltonian

$H_{\rm SIA}=\left(\begin{array}[]{cccc}0&\mathrm{i}\chi\sin
4\theta&0&\zeta\sin^{2}2\theta\\\ -\mathrm{i}\chi\sin
4\theta&0&\zeta\sin^{2}2\theta&0\\\
0&\zeta\sin^{2}2\theta&0&\mathrm{i}\chi\sin 4\theta\\\
\zeta\sin^{2}2\theta&0&-\mathrm{i}\chi\sin 4\theta&0\end{array}\right)\,.$ (5)

The mixing parameters $\chi$ and $\zeta$ are nonzero if both structure
inversion asymmetry and the cubic shape of lattice unit cells are taken into
account. Note also that $H_{\rm SIA}$ vanishes in (001)-grown QWs.

The parameters $\eta$, $\gamma$, $\chi$, and $\zeta$ are calculated in Sec.
III in the framework of the 6-band $\bm{k}$$\cdot$$\bm{p}$ model. An
estimation for HgTe/Cd0.7Hg0.3Te QWs with the critical thickness
${d_{c}\approx 6.7}$ nm gives $\eta,\gamma\sim 5$ meV and $\zeta,\chi\sim 0.1$
meV in the electric field $E_{z}=15$ kV/cm.

The Hamiltonians $H_{\rm IIA/BIA}$ (4) and $H_{\rm SIA}$ (5) depend on the QW
growth direction defined by the $\theta$ angle. Straightforward
diagonalization of the Hamiltonian (2) yields four dispersion branches

$\displaystyle E_{1,4}$ $\displaystyle=$
$\displaystyle\mp\sqrt{\delta^{2}+\gamma_{\theta}^{2}+\eta_{\theta}^{2}+\zeta_{\theta}^{2}+\chi_{\theta}^{2}+A^{2}k^{2}+2AK}\>,$
$\displaystyle E_{2,3}$ $\displaystyle=$
$\displaystyle\mp\sqrt{\delta^{2}+\gamma_{\theta}^{2}+\eta_{\theta}^{2}+\zeta_{\theta}^{2}+\chi_{\theta}^{2}+A^{2}k^{2}-2AK}\>,$
(6)

where $k^{2}=k_{x}^{2}+k_{y}^{2}$,

$K=\sqrt{(\gamma_{\theta}^{2}+\zeta_{\theta}^{2})k^{2}+\left(\chi_{\theta}k_{x}+\eta_{\theta}k_{y}\right)^{2}}\>,$
(7)

and

$\displaystyle\gamma_{\theta}$ $\displaystyle=$ $\displaystyle\gamma\cos
2\theta\>,~{}~{}~{}\eta_{\theta}=\eta\sin 2\theta\>,$
$\displaystyle\zeta_{\theta}$ $\displaystyle=$
$\displaystyle\zeta\sin^{2}2\theta\>,~{}~{}\chi_{\theta}=\chi\sin 4\theta\>.$
(8)

In the following subsections we analyze the fine structure of Dirac states in
QWs with different crystallographic orientations for different mixing
mechanisms.

###  II.1 Interface and bulk inversion asymmetry

Dirac states in HgTe/CdHgTe QWs with symmetric confinement potential are
described by the Hamiltonian $H=H_{0}+H_{\rm IIA/BIA}$ with $H_{0}$ and
$H_{\rm IIA/BIA}$ given by Eqs. (3) and (4), respectively. In this case, Eq.
(II) yields

$\displaystyle E_{1,4}$ $\displaystyle=$
$\displaystyle\mp\sqrt{\delta^{2}+\gamma_{\theta}^{2}+\eta_{\theta}^{2}+A^{2}k^{2}+2A\sqrt{\gamma_{\theta}^{2}k^{2}+\eta_{\theta}^{2}k_{y}^{2}}}\>,$
$\displaystyle E_{2,3}$ $\displaystyle=$
$\displaystyle\mp\sqrt{\delta^{2}+\gamma_{\theta}^{2}+\eta_{\theta}^{2}+A^{2}k^{2}-2A\sqrt{\gamma_{\theta}^{2}k^{2}+\eta_{\theta}^{2}k_{y}^{2}}}\>.$

Figure 1:  Energy spectra of Dirac states in (001)-, (013)- and (011)-grown
HgTe/CdHgTe QWs of critical thickness with interface- and bulk-inversion
asymmetry included. The spectra are calculated after Eq. (II.1) for $\delta=0$
and $\gamma=\eta$. Color decodes the peudospin projection onto the QW normal
(see text for details): blue and red correspond to $\sigma_{z}=-1$ and
$\sigma_{z}=+1$, respectively, whereas purple corresponds to $\sigma_{z}=0$.
Figure 2:  Energy spectra of Dirac states in (001)-, (013)- and (011)-grown
HgTe/CdHgTe QWs of close-to-critical thickness with interface- and bulk-
inversion asymmetry included. The spectra are calculated after Eq. (II.1) for
$\gamma=\eta=2\delta$. Color decodes the peudospin projection onto the QW
normal: blue and red correspond to $\sigma_{z}=-1$ and $\sigma_{z}=+1$,
respectively, whereas purple corresponds to $\sigma_{z}=0$.

Figure 1 shows the energy spectra of Dirac states in (001), (013), and (011)
QWs of critical thickness ($\delta=0$) with the IIA/BIA term included. The
(001) orientation corresponds to $\theta=0$. In this case, the energy spectrum
consists of two non-degenerate (Weyl) cones shifted vertically (along the
energy axis) with respect to each other Tarasenko2015 . The Weyl points are
located at $\bm{k}=0$ and the energies $E=\pm\gamma$. The (011) orientation
corresponds to $\theta=\pi/4$. In such QWs, the IIA/BIA interaction splits the
Dirac cone into two Weyl cones shifted along $k_{y}$ with respect to each
other. The Weyl points are located at $\bm{k}=(0,\pm\eta/A)$. The spectrum in
(013) and general $(0lh)$ QWs is an intermediate case between the spectra in
(001) and (011) structures. Now, there are four Weyl points in the energy
spectrum. Two points are located at
$\bm{k}=(0,\pm\sqrt{\gamma_{\theta}^{2}+\eta_{\theta}^{2}}/A)$ and zero
energy, while the other two points are at $\bm{k}=0$ and the energies
$E=\pm\sqrt{\gamma_{\theta}^{2}+\eta_{\theta}^{2}}$.

The color in Fig. 1 decodes the projection of pseudospin onto the QW normal
$\sigma_{z}$ defined by
$\sigma_{z}=|c_{1}|^{2}+|c_{2}|^{2}-|c_{3}|^{2}-|c_{4}|^{2}$, where $c_{j}$
are the coefficients of decomposition of a wave function $\psi$ over the basis
functions (II). It illustrates the relative contribution of the “spin-up”
($\left|E1,+1/2\right\rangle$ and $\left|H1,+3/2\right\rangle$) and “spin-
down” ($\left|E1,-1/2\right\rangle$ and $\left|H1,-3/2\right\rangle$) blocks
in the given state $\psi$. In (001) QWs, the Weyl cones are formed by the
“spin-up” and “spin-down” blocks in equal portions and $\sigma_{z}=0$ (purple
color) for all eigen states. In contrast, the split Weyl cones in (011) QWs
are formed by pure “spin-up” and “spin-down” states and characterized by
$\sigma_{z}=+1$ (red color) and $\sigma_{z}=-1$ (red color) pseudospin
projections.

The energy spectra of (001), (013), and (011) QWs of close-to-critical
thickness (with the gap $2|\delta|$) are shown in Fig. 2. In (001) QWs, the
spectrum is given by $E=\pm\sqrt{\delta^{2}+(A|\bm{k}|\pm\gamma)^{2}}$ and the
band extrema are situated at the loop with $|\bm{k}|=|\gamma/A|$. In (011)
QWs, the spectrum has the form
$E=\pm\sqrt{\delta^{2}+A^{2}k_{x}^{2}+(Ak_{y}\pm\eta)^{2}}$ and consists of
the branches with the $\sigma_{z}=\pm 1$ pseudospin projections. In general
case of $(0lh)$ orientation, e.g., (013), the band extrema are situated at the
points $\bm{k}=(0,\pm\sqrt{\gamma_{\theta}^{2}+\eta_{\theta}^{2}}/A)$. The
contours of constant energy are toric sections, in particular, at
$|E|>\sqrt{\delta^{2}+\gamma_{\theta}^{2}+\eta_{\theta}^{2}}$ the isoenergy
contours are ovals elongated along $k_{y}$.

###  II.2 Structure inversion asymmetry

Here, we study the influence of structure inversion asymmetry on the energy
spectrum of Dirac states. For this purpose, we consider the Hamiltonian
$H=H_{0}+H_{\rm SIA}$ with $H_{0}$ and $H_{\rm SIA}$ given by Eqs. (3) and
(5), respectively. In this case, Eq. (II) yields

$\displaystyle E_{1,4}$ $\displaystyle=$
$\displaystyle\mp\sqrt{\delta^{2}+\zeta_{\theta}^{2}+\chi_{\theta}^{2}+A^{2}k^{2}+2A\sqrt{\zeta_{\theta}^{2}k^{2}+\chi_{\theta}^{2}k_{x}^{2}}}\>,$
$\displaystyle E_{2,3}$ $\displaystyle=$
$\displaystyle\mp\sqrt{\delta^{2}+\zeta_{\theta}^{2}+\chi_{\theta}^{2}+A^{2}k^{2}-2A\sqrt{\zeta_{\theta}^{2}k^{2}+\chi_{\theta}^{2}k_{x}^{2}}}\>.$

Figure 3:  Energy spectra of Dirac states in (001)-, (013)- and (011)-oriented
QWs of critical thickness with structure-inversion asymmetry included. The
spectra are calculated after Eq. (II.2) for $\delta=0$ and $\zeta=\chi$. Color
decodes the peudospin projection onto the QW normal: blue and red colors
correspond to $\sigma_{z}=-1$ and $\sigma_{z}=+1$, respectively, whereas
purple color corresponds to $\sigma_{z}=0$.

Figure 3 shows the energy dispersions given by Eq. (II.2) for (001), (013),
and (011) QWs of critical thickness ($\delta=0$). In (001) QWs, structure
inversion asymmetry does not contribute to the mixing of the basis states at
$\bm{k}=0$. As a result, the point $\bm{k}=0$ remains four-fold degenerate
just as it is in the Bernevig-Hughes-Zhang model. In QWs of other
orientations, the SIA interaction lifts the four-fold degeneracy at $\bm{k}=0$
and splits the Dirac cone. Generally, there are four Weyl points in the energy
spectrum located at
$\bm{k}=(\pm\sqrt{\zeta_{\theta}^{2}+\chi_{\theta}^{2}}/A,0)$ and zero energy
and at $\bm{k}=0$ and the energies
$E=\pm\sqrt{\zeta_{\theta}^{2}+\chi_{\theta}^{2}}$, respectively.
Interestingly, the SIA interaction in (011) QWs splits the Dirac cone in a way
similar to the IIA/BIA interaction does in (001) QWs.

In QW structures with a gap (not shown), the band extrema are located at
$\bm{k}=(\pm\sqrt{\zeta_{\theta}^{2}+\chi_{\theta}^{2}}/A,0)$, in particular,
at $\bm{k}=0$ in (001) QWs and at the loop with $|\bm{k}|=|\zeta/A|$ in (011)
QWs.

###  II.3 Interplay of IIA/BIA and SIA

In real QW structures, all types of asymmetry, including bulk, interface, and
structure inversion asymmetry, are present. The dispersion branches in that
case are given by the general Eq. (II). Figure 4 shows the energy spectra of
Dirac states in such QWs with asymmetric confinement potential and grown along
different crystallographic orientations.

In (001) QWs, the splitting of the Dirac cone at $\bm{k}=0$ is determined by
the BIA/IIA term and the energy spectrum coincides with the one shown in Fig.
1. The spectrum of asymmetric (011) QWs is qualitatively similar to the
spectrum of symmetric $(013)$ QWs. There are four Weyl points: two of them
located at $\bm{k}=(0,\pm\sqrt{\zeta_{\theta}^{2}+\eta_{\theta}^{2}}/A)$ and
zero energy and the other two located at $\bm{k}=0$ and the energies
$E=\pm\sqrt{\zeta_{\theta}^{2}+\eta_{\theta}^{2}}$.

The spectrum of a general $(0lh)$-grown QW with asymmetric confinement
potential is shown in the central panel in Fig. 4. The Weyl points at zero
energy are located at the wave vectors

$\bm{k}=\pm\sqrt{\frac{\gamma_{\theta}^{2}+\eta_{\theta}^{2}+\zeta_{\theta}^{2}+\chi_{\theta}^{2}}{A^{2}\left(\eta_{\theta}^{2}+\chi_{\theta}^{2}\right)}}\left(\chi_{\theta},\eta_{\theta}\right)\>.$
(11)

Interestingly, the position of these points in the $\bm{k}$ space is not
pinned to a specific in-plane direction. The angle between the line connecting
the Weyl points and the $k_{x}$ axis, $\arctan(\eta_{\theta}/\chi_{\theta})$,
depends on the SIA parameter $\chi$ and, therefore, can be controlled by an
external electric field applied along the QW normal, e.g., by gate voltage.

Figure 4:  Energy spectra of Dirac states in (001)-, (013)- and (011)-grown
HgTe/CdHgTe QWs of critical thickness with interface-, bulk-, and structure-
inversion asymmetry included. The spectra are calculated after Eq. (II) for
$\delta=0$ and $\gamma=\eta=\zeta=\chi$. Color decodes the peudospin
projection onto the QW normal: blue and red correspond to $\sigma_{z}=-1$ and
$\sigma_{z}=+1$, respectively, whereas purple corresponds to $\sigma_{z}=0$.

##  III 6-band kp theory

In this section, we calculate the energy spectrum and derive the parameters of
the effective Hamiltonian (2) using the extended 6-band
$\bm{k}$$\cdot$$\bm{p}$ theory. The conduction-band and valence-band states in
HgTe/CdHgTe structure are mainly formed from the $\Gamma_{6}$ and $\Gamma_{8}$
bands, which are 2-fold and 4-fold degenerate, respectively, at $\bm{k}=0$ in
the bulk crystal Novik2005 ; Gerchikov1989 ; Bernevig2006 . Taking into
account $\bm{k}$$\cdot$$\bm{p}$ mixing, deformation interaction, and interface
mixing, we present the corresponding 6-band Hamiltonian in the form

$\mathcal{H}=\mathcal{H}_{\rm kp}+\mathcal{H}_{\rm def}+\mathcal{H}_{\rm
int}=\left(\begin{array}[]{c c}\mathcal{H}_{66}&\mathcal{H}_{68}\\\
\mathcal{H}^{\dagger}_{68}&\mathcal{H}_{88}\end{array}\right),$ (12)

where $\mathcal{H}_{66}$ is the 2$\times$2 matrix in the
$\left|\Gamma_{6},+1/2\right\rangle$, $\left|\Gamma_{6},-1/2\right\rangle$
basis, $\mathcal{H}_{88}$ is the 4$\times$4 matrix in the
$\left|\Gamma_{8},+3/2\right\rangle$,
$\left|\Gamma_{8},+1/2\right\rangle,\left|\Gamma_{8},-1/2\right\rangle,\left|\Gamma_{8},-3/2\right\rangle$
basis, $\mathcal{H}_{68}$ is the 2$\times$4 matrix, which couples the
$\Gamma_{6}$ and $\Gamma_{8}$ blocks, and $\mathcal{H}^{\dagger}_{68}$ is the
Hermitian conjugate matrix.

The isotropic version of the Hamiltonian $\mathcal{H}_{\rm kp}$ is the well-
known 6-band Kane Hamiltonian which is often used to model the conduction and
valence bands in narrow gap III-V semiconductors. This model, however, is
simplified and does not take into account the real symmetry ($T_{d}$ point
group) of the zinc-blende lattice. The latter is essential to describe the
fine structure of Dirac states. Therefore, we go beyond the isotropic Kane
Hamiltonian and make use of the so-called extended Kane model Winkler_book ,
which takes into account the cubic shape and inversion asymmetry of the
lattice.

The extended Kane Hamiltonian $\mathcal{H}_{\rm kp}$ is constructed by the
methods of group representation theory BirPikus . In the cubic axes
$x^{\prime}\parallel[100]$, $y^{\prime}\parallel[010]$, and
$z^{\prime}\parallel[001]$, the $\mathcal{H}_{66}^{\rm kp}$ and
$\mathcal{H}_{88}^{\rm kp}$ blocks to the second order in the wave vector
$\bm{k}$ have the form

$\mathcal{H}_{66}^{\rm kp}=U_{6}+\frac{\hbar^{2}k^{2}}{2m_{e}^{\prime}}\>,$
(13) $\mathcal{H}_{88}^{\rm
kp}=U_{8}+\frac{\hbar^{2}}{2m_{0}}\left[-\left(\gamma_{1}^{\prime}+\frac{5}{2}\gamma_{2}^{\prime}\right)k^{2}+2\gamma_{2}^{\prime}(\bm{J}\bm{k})^{2}+\right.\\\
\left.+2(\gamma_{3}^{\prime}-\gamma_{2}^{\prime})\sum_{i\neq
j}\\{J_{i}J_{j}\\}_{s}k_{i}k_{j}\right]+\frac{4\varkappa_{0}}{\sqrt{3}}\bm{V}\bm{k}\>,$
(14)

where $U_{6}$ and $U_{8}$ are the energies of the $\Gamma_{6}$ and
$\Gamma_{8}$ bands at $\bm{k}=0$,
$\bm{k}=(k_{x^{\prime}},k_{y^{\prime}},k_{z^{\prime}})$ is the wave vector,
$\gamma_{1}^{\prime}$, $\gamma_{2}^{\prime}$, $\gamma_{3}^{\prime}$, and
$m_{e}^{\prime}$ are the contributions to the Luttinger parameters and the
effective mass, respectively, from remote bands and free electron dispersion,
$\bm{J}=(J_{x^{\prime}},J_{y^{\prime}},J_{z^{\prime}})$ is the vector composed
of the momentum-3/2 matrices,
$\bm{V}=(V_{x^{\prime}},V_{y^{\prime}},V_{z^{\prime}})$, where
$V_{x^{\prime}}=\\{J_{x^{\prime}},J_{y^{\prime}}^{2}-J_{z^{\prime}}^{2}\\}$
and the other components of $\bm{V}$ are derived by the cyclic permutation of
the subscripts, $\\{A,B\\}_{s}=(AB+BA)/2$ is the symmetrized product of the
operators $A$ and $B$, and $\varkappa_{0}$ is a band parameter.

To construct the $\mathcal{H}_{68}$ block, one notes that the direct product
$\Gamma_{6}\times\Gamma_{8}^{*}$ is decomposed into the irreducible
representations $\Gamma_{3}+\Gamma_{4}+\Gamma_{5}$. The sets
$\\{k_{x^{\prime}},k_{y^{\prime}},k_{z^{\prime}}\\}$ and
$\\{k_{y^{\prime}}k_{z^{\prime}},k_{x^{\prime}}k_{z^{\prime}},k_{x^{\prime}}k_{y^{\prime}}\\}$
transform according to the vector representation $\Gamma_{5}$ whereas the pair
$\\{2k_{z^{\prime}}^{2}-k_{x^{\prime}}^{2}-k_{y^{\prime}}^{2},\sqrt{3}(k_{x^{\prime}}^{2}-k_{y^{\prime}}^{2})\\}$
transforms according to the $\Gamma_{3}$ representation. The combinations that
transform according to the pseudo-vector representation $\Gamma_{4}$ are cubic
in $\bm{k}$ and are not considered. Thus, the $\mathcal{H}_{68}^{\rm kp}$
block is given by Winkler_book

$\mathcal{H}_{68}^{{\rm
kp}{\dagger}}=\left(\begin{array}[]{cc}\frac{\mathrm{i}}{\sqrt{2}}(B_{+}k_{+}k_{z^{\prime}}-Pk_{-})&\frac{\mathrm{i}}{3\sqrt{2}}B_{-}(2k_{z^{\prime}}^{2}-k_{\parallel}^{2})\\\
&\\\
\sqrt{\frac{2}{3}}(\mathrm{i}Pk_{z^{\prime}}+B_{+}k_{x^{\prime}}k_{y^{\prime}})&\frac{\mathrm{i}}{\sqrt{6}}(B_{+}k_{+}k_{z^{\prime}}-Pk_{-})\\\
-\frac{\mathrm{i}}{\sqrt{6}}B_{-}(k_{y^{\prime}}^{2}-k_{x^{\prime}}^{2})&\\\
&\\\
\frac{\mathrm{i}}{\sqrt{6}}(Pk_{+}+B_{+}k_{-}k_{z^{\prime}})&\sqrt{\frac{2}{3}}(\mathrm{i}Pk_{z^{\prime}}+B_{+}k_{x^{\prime}}k_{y^{\prime}})\\\
&+\frac{\mathrm{i}}{\sqrt{6}}B_{-}(k_{y^{\prime}}^{2}-k_{x^{\prime}}^{2})\\\
&\\\
\frac{\mathrm{i}}{3\sqrt{2}}B_{-}(k_{\parallel}^{2}-2k_{z^{\prime}}^{2})&\frac{\mathrm{i}}{\sqrt{2}}(Pk_{+}+B_{+}k_{-}k_{z^{\prime}})\end{array}\right),$
(15)

where $P=\mathrm{i}(\hbar/m_{0})p_{cv}$ is the Kane matrix element, $B_{\pm}$
are band parameters, $k_{\pm}=k_{x^{\prime}}\pm\mathrm{i}k_{y^{\prime}}$, and
$k_{\parallel}^{2}=k_{x^{\prime}}^{2}+k_{y^{\prime}}^{2}$. Note that the
definition of $P$ and $B_{\pm}$ in Eq. (15) differs by the factor of
$\mathrm{i}$ from that in Ref. Winkler_book, .

The extended Kane Hamiltonian given by the blocks (13)-(15) reflects the real
symmetry of the zinc-blende lattice including its cubic shape and the lack of
space inversion center. The isotropic centrosymmetric approximation
corresponds to $\gamma_{2}^{\prime}=\gamma_{3}^{\prime}$, $\varkappa_{0}=0$,
and $B_{\pm}=0$. The nonzero difference
$\gamma_{2}^{\prime}-\gamma_{3}^{\prime}$ takes into account the cubic
anisotropy of the unit cell whereas the nonzero parameters $\varkappa_{0}$ and
$B_{\pm}$ reflect bulk inversion asymmetry. The $B_{-}$ parameter couples the
functions with the opposite spin projections and, hence, is expected to be
smaller than $B_{+}$. We neglect this parameter in the following calculations.

The parameters of the effective $6$-band Hamiltonian $\mathcal{H}_{\rm kp}$
can be expressed via the coupling parameters and the energy gaps in multi-band
$\bm{k}$$\cdot$$\bm{p}$ theory. The results of such calculations in the
14-band $\bm{k}$$\cdot$$\bm{p}$ model Pikus1988 ; Jancu2005 ; Durnev2014 ,
which includes the $\Gamma_{7}$ valence band and the remote
$\Gamma_{8}^{\prime}$ and $\Gamma_{7}^{\prime}$ conduction bands in addition
to the considered $\Gamma_{6}$ and $\Gamma_{8}$ bands, are summarized in Tab.
1.

$m_{0}/m_{e}^{\prime}$ | $\gamma_{1}^{\prime}$ | $\gamma_{2}^{\prime}$ | $\gamma_{3}^{\prime}$ | $\varkappa_{0}$ | $B_{+}$ | $B_{-}$
---|---|---|---|---|---|---
$\dfrac{2m_{0}P^{2}}{3\hbar^{2}(E_{g}+\Delta)}$ | $\dfrac{4m_{0}Q^{2}}{3\hbar^{2}(E_{g}+E_{g}^{\prime})}$ | $-\dfrac{m_{0}Q^{2}}{3\hbar^{2}(E_{g}+E_{g}^{\prime})}$ | $\dfrac{m_{0}Q^{2}}{3\hbar^{2}(E_{g}+E_{g}^{\prime})}$ | 0 | $\dfrac{QP^{\prime}(E_{g}+2E_{g}^{\prime})}{E_{g}^{\prime}(E_{g}+E_{g}^{\prime})}$ | $-\dfrac{QP^{\prime}\Delta^{\prime}(E_{g}^{2}+2E_{g}E_{g}^{\prime}+2E_{g}^{\prime 2})}{2E_{g}^{\prime 2}(E_{g}+E_{g}^{\prime})^{2}}$
Table 1: Parameters of the 6-band Hamiltonian $\mathcal{H}_{\rm kp}$
calculated in the 14-band $\bm{k}$$\cdot$$\bm{p}$ model, which includes the
$\Gamma_{7}$ and $\Gamma_{8}$ valence bands, the $\Gamma_{6}$ conduction band,
and the $\Gamma_{7}^{\prime}$ and $\Gamma_{8}^{\prime}$ remote conduction
bands Pikus1988 ; Jancu2005 ; Durnev2014 . Here, $E_{g}$ and $E_{g}^{\prime}$
are the gaps at $\bm{k}=0$ between $\Gamma_{6}$ and $\Gamma_{8}$ and between
$\Gamma_{7}^{\prime}$ and $\Gamma_{6}$, respectively, $\Delta$ and
$\Delta^{\prime}$ are the spin-orbit splittings of the valence band and the
remote conduction band, $P$, $P^{\prime}$, and $Q$ are the
$\Gamma_{7,8}-\Gamma_{6}$, $\Gamma^{\prime}_{7,8}-\Gamma_{6}$, and
$\Gamma_{7,8}-\Gamma^{\prime}_{7,8}$ coupling parameters, respectively. We
assume that $P^{\prime}\ll P$, $\Delta^{\prime}\ll E_{g}^{\prime}$ and also
neglect $\bm{k}$-independent spin-orbit coupling between the $\Gamma_{7,8}$,
and $\Gamma^{\prime}_{7,8}$ bands.

Layers in epitaxial HgTe/CdHgTe structures are typically strained because of
considerable mismatch (of about 0.3$\,\%$) between HgTe and CdTe lattice
constants. The 6-band strain Hamiltonian $\mathcal{H}_{\rm def}$ can be
constructed in a way similar to the $\bm{k}$$\cdot$$\bm{p}$ Hamiltonian
Pikus1988 ; BirPikus . Such a procedure yields the diagonal blocks

$\mathcal{H}_{66}^{\rm def}=\Xi_{c}{\rm Tr}\epsilon\>,$ (16)
$\mathcal{H}_{88}^{\rm def}=\left(a+\frac{5}{4}b\right){\rm
Tr}\epsilon-b\sum_{i,j}\\{J_{i}J_{j}\\}_{s}\epsilon_{ij}\\\
+\left(b-\frac{d}{\sqrt{3}}\right)\sum_{i\neq
j}\\{J_{i}J_{j}\\}_{s}\epsilon_{ij}\>,$ (17)

where $\Xi_{c}$ is the $\Gamma_{6}$-band deformation potential, $a$, $b$, and
$d$ are the $\Gamma_{8}$-band deformation potentials, $\mathcal{H}_{88}^{\rm
def}$ is the Bir-Pikus Hamiltonian, $\epsilon$ is the strain tensor, and the
off-diagonal blocks

$\mathcal{H}_{68}^{{\rm
def}{\dagger}}=\left(\begin{array}[]{cc}-\Xi_{cv}\dfrac{\epsilon_{y^{\prime}z^{\prime}}-\mathrm{i}\epsilon_{x^{\prime}z^{\prime}}}{\sqrt{2}}&\mathrm{i}\Xi_{cv}^{\prime}\dfrac{2\epsilon_{z^{\prime}z^{\prime}}-\epsilon_{x^{\prime}x^{\prime}}-\epsilon_{y^{\prime}y^{\prime}}}{3\sqrt{2}}\\\
&\\\
\sqrt{\dfrac{2}{3}}\,\Xi_{cv}\epsilon_{x^{\prime}y^{\prime}}&-\Xi_{cv}\dfrac{\epsilon_{y^{\prime}z^{\prime}}-\mathrm{i}\epsilon_{x^{\prime}z^{\prime}}}{\sqrt{6}}\\\
+\mathrm{i}\Xi_{cv}^{\prime}\dfrac{\epsilon_{x^{\prime}x^{\prime}}-\epsilon_{y^{\prime}y^{\prime}}}{\sqrt{6}}&\\\
&\\\
\Xi_{cv}\dfrac{\epsilon_{y^{\prime}z^{\prime}}+\mathrm{i}\epsilon_{x^{\prime}z^{\prime}}}{\sqrt{6}}&\sqrt{\dfrac{2}{3}}\,\Xi_{cv}\epsilon_{x^{\prime}y^{\prime}}\\\
&-\mathrm{i}\Xi_{cv}^{\prime}\dfrac{\epsilon_{x^{\prime}x^{\prime}}-\epsilon_{y^{\prime}y^{\prime}}}{\sqrt{6}}\\\
&\\\
\mathrm{i}\Xi_{cv}^{\prime}\dfrac{\epsilon_{x^{\prime}x^{\prime}}+\epsilon_{y^{\prime}y^{\prime}}-2\epsilon_{z^{\prime}z^{\prime}}}{3\sqrt{2}}&\Xi_{cv}\dfrac{\epsilon_{y^{\prime}z^{\prime}}+\mathrm{i}\epsilon_{x^{\prime}z^{\prime}}}{\sqrt{2}}\end{array}\right),$
(18)

where $\Xi_{cv}$ and $\Xi_{cv}^{\prime}$ are the inter-band deformation
potentials. Note that $\Xi_{cv}$ and $\Xi_{cv}^{\prime}$ vanish in
centrosymmetric crystals. The values of $\Xi_{cv}$ and $\Xi_{cv}^{\prime}$ for
HgTe and CdTe are not known. We expect $\Xi_{cv}^{\prime}$ to be much smaller
than $\Xi_{cv}$ and neglect it below.

Interfaces in heterostructures introduce additional mechanics of Bloch state
coupling. In zinc-blende structures, interface inversion asymmetry related to
the anisotropy of chemical bonds leads to light-hole–heavy-hole mixing
Aleiner1992 ; Ivchenko1996 . For an arbitrary crystallographic orientation of
the interface, the mixing can be described by the Hamiltonian

$\mathcal{H}_{88}^{\rm
int}=\frac{\hbar^{2}t_{l-h}}{\sqrt{3}a_{0}m_{0}}\delta(\bm{r}\cdot\bm{n}+r_{\rm
int})\sum_{i}\\{J_{i}J_{i+1}\\}_{s}n_{i+2}\>,$ (19)

where $t_{l-h}$ is a dimensionless mixing parameter, $a_{0}$ is the lattice
constant, $\bm{n}=(n_{x^{\prime}},n_{y^{\prime}},n_{z^{\prime}})$ is the unit
vector directed along the interface normal, say from CdHgTe to HgTe,
$\bm{r}\cdot\bm{n}+r_{\rm int}=0$ is the equation of the interface plane, and
$r_{\rm int}$ is the distance between the interface and the coordinate origin.

To describe the electron and hole states in $(0lh)$-oriented QWs we rotate the
Hamiltonian (12) in the reference frame relevant to the QW. The transition
from the reference frame $(x^{\prime},y^{\prime},z^{\prime})$ to $(x,y,z)$
corresponds to the rotation around the $x$ axis by the angle $\theta$. Under
this rotation, the basis functions $\left|\Gamma_{6},m\right\rangle$ and
$\left|\Gamma_{8},m\right\rangle$ transform as the basis functions of the
$1/2$ and $3/2$ angular momentum, respectively. Therefore, the 6-band
Hamiltonian (12) in the QW reference frame assumes the form

$\mathcal{H}_{\rm xyz}=R^{-1}\mathcal{H}R\,,$ (20)

where

$R=\begin{pmatrix}R_{6}&0\\\ 0&R_{8}\end{pmatrix}\,,$ (21)

$R_{6}=\exp(is_{x}\theta)$ and $R_{8}=\exp(iJ_{x}\theta)$ are the $2\times 2$
and $4\times 4$ rotation matrices, respectively, and $s_{x}=\sigma_{x}/2$. The
wave vector and strain tensor components are transformed as

$\displaystyle k_{x^{\prime}}$ $\displaystyle=k_{x}\>,$ $\displaystyle
k_{y^{\prime}}$ $\displaystyle=k_{y}\cos\theta+k_{z}\sin\theta\>,$
$\displaystyle k_{z^{\prime}}$ $\displaystyle=k_{z}\cos\theta-
k_{y}\sin\theta\>,$ (22)

and

$\displaystyle\epsilon_{x^{\prime}x^{\prime}}$
$\displaystyle=\epsilon_{xx}\>,$
$\displaystyle\epsilon_{x^{\prime}y^{\prime}}$
$\displaystyle=\epsilon_{xy}\cos\theta+\epsilon_{xz}\sin\theta\>,$
$\displaystyle\epsilon_{x^{\prime}z^{\prime}}$
$\displaystyle=\epsilon_{xz}\cos\theta-\epsilon_{xy}\sin\theta\>,$
$\displaystyle\epsilon_{y^{\prime}y^{\prime}}$
$\displaystyle=\epsilon_{yy}\cos^{2}\theta+\epsilon_{zz}\sin^{2}\theta+\epsilon_{yz}\sin
2\theta\>,$ $\displaystyle\epsilon_{y^{\prime}z^{\prime}}$
$\displaystyle=\epsilon_{yz}\cos
2\theta+(1/2)(\epsilon_{zz}-\epsilon_{yy})\sin 2\theta\>,$
$\displaystyle\epsilon_{z^{\prime}z^{\prime}}$
$\displaystyle=\epsilon_{zz}\cos^{2}\theta+\epsilon_{yy}\sin^{2}\theta-\epsilon_{yz}\sin
2\theta\>.$ (23)

We assume that the strain in the QW is caused by mismatch between the lattice
constant $a$ of the buffer layer and the lattice constant $a_{0}$ of the
unstrained QW material. Then, the in-plane strain in the QW is determined from
the lattice matching condition and given by
$\epsilon_{xx}=\epsilon_{yy}=a/a_{0}-1$. The other strain components are found
from the elastic energy minimization and given by Dantscher2015

$\displaystyle\epsilon_{zz}$
$\displaystyle=\dfrac{c_{11}^{2}+2c_{11}(c_{12}-c_{44})+c_{12}(-3c_{12}+10c_{44})-(c_{11}+3c_{12})(c_{11}-c_{12}-2c_{44})\cos
4\theta}{-c_{11}^{2}-6c_{11}c_{44}+c_{12}(c_{12}+2c_{44})+(c_{11}+c_{12})(c_{11}-c_{12}-2c_{44})\cos
4\theta}\epsilon_{xx}\>,$ $\displaystyle\epsilon_{yz}$
$\displaystyle=\dfrac{(c_{11}+2c_{12})(c_{11}-c_{12}-2c_{44})\sin
4\theta}{-c_{11}^{2}-6c_{11}c_{44}+c_{12}(c_{12}+2c_{44})+(c_{11}+c_{12})(c_{11}-c_{12}-2c_{44})\cos
4\theta}\epsilon_{xx}\>,\;\;\;\epsilon_{xz}=0\>,\;\;\;\epsilon_{xy}=0\>,$ (24)

where $c_{11}$, $c_{12}$, and $c_{44}$ are the elastic constants. Note that
$c_{11}-c_{12}-2c_{44}=0$ in the model of isotropic elastic medium and, hence,
$\epsilon_{yz}\neq 0$ due to cubic shape of the unit cell.

Figure 5:  Dispersions in (013)-grown HgTe/Cd0.7Hg0.3Te QW structure with the
QW thickness $w=6.7$ nm for the in-plane directions (a) $\bm{k}\parallel
x\parallel[100]$ and (b) $\bm{k}\parallel y\parallel[03\bar{1}]$. The
dispersions are calculated numerically in the 6-band $\bm{k}$$\cdot$$\bm{p}$
theory using HgTe and CdTe parameters from Refs. Novik2005, ; Dantscher2015, ;
Tarasenko2015,  and the electric field $E_{z}=15$ kV/cm.

Figure 5 shows cross sections of the energy spectrum calculated numerically
for asymmetric (013)-grown HgTe/Cd0.7Hg0.3Te QW of critical thickness.
Splitting of Dirac states at $\bm{k}=0$ and the in-plane anisotropy of the
energy spectrum are readily seen.

To obtain the parameters of the effective Hamiltonian (2) we solve the
Schrödinger equation $\mathcal{H}_{\rm xyz}^{(\rm iso)}\Psi=E\Psi$ for zero
in-plane wave vector, where $\mathcal{H}_{\rm xyz}^{(\rm iso)}$ is the
isotropic part of the Hamiltonian (20), and find the functions $\left|E1,\pm
1/2\right\rangle$ and $\left|H1,\pm 3/2\right\rangle$. Then, we project the
Hamiltonian $(\mathcal{H}_{\rm xyz}-\mathcal{H}_{\rm xyz}^{(\rm iso)})$ onto
the basis states $\left|E1,\pm 1/2\right\rangle$ and $\left|H1,\pm
3/2\right\rangle$. This procedure yields the effective Hamiltonian (2) with
$A=P/\sqrt{2}\int f_{1}(z)f_{3}(z)dz$,

$\displaystyle\eta$ $\displaystyle=$
$\displaystyle\frac{\hbar^{2}t_{l-h}}{2m_{0}a_{0}}\left[f_{3}(w/2)f_{4}(w/2)-f_{3}(-w/2)f_{4}(-w/2)\right]$
$\displaystyle+$ $\displaystyle\frac{1}{\sqrt{2}}\int
dzf_{1}(z)\Xi_{cv}\left[\epsilon_{yz}\cot
2\theta+\frac{\epsilon_{zz}-\epsilon_{yy}}{2}\right]f_{3}(z)$ $\displaystyle-$
$\displaystyle\frac{1}{2\sqrt{2}}\int
dzf_{1}(z)\partial_{z}B_{+}\partial_{z}f_{3}(z)\>,$ $\displaystyle\gamma$
$\displaystyle=$
$\displaystyle\frac{\hbar^{2}t_{l-h}}{2m_{0}a_{0}}\left[f_{3}(w/2)f_{4}(w/2)-f_{3}(-w/2)f_{4}(-w/2)\right]\>,$
$\displaystyle\chi$ $\displaystyle=$
$\displaystyle\frac{\sqrt{3}\hbar^{2}}{4m_{0}}\int
dzf_{4}(z)\partial_{z}(\gamma_{2}^{\prime}-\gamma_{3}^{\prime})\partial_{z}f_{3}(z)$
$\displaystyle-$ $\displaystyle\frac{1}{2}\int
dzf_{3}(z)\left\\{\frac{\left[(d+\sqrt{3}b)+(d-\sqrt{3}b)\cos
4\theta\right]\epsilon_{yz}}{\sin 4\theta}\right.$
$\displaystyle\left.+(d-\sqrt{3}b)\frac{\epsilon_{zz}-\epsilon_{yy}}{2}\right\\}f_{4}(z)\>,$
$\displaystyle\zeta$ $\displaystyle=$
$\displaystyle\frac{\sqrt{3}\hbar^{2}}{4m_{0}}\int
dzf_{4}(z)\partial_{z}(\gamma_{2}^{\prime}-\gamma_{3}^{\prime})\partial_{z}f_{3}(z)$
$\displaystyle+$ $\displaystyle\frac{1}{2}\int
dzf_{3}(z)(\sqrt{3}b-d)\left[\epsilon_{yz}\cot
2\theta+\frac{\epsilon_{zz}-\epsilon_{yy}}{2}\right]f_{4}(z)\>,$

where $w$ is the QW thickness and the strain tensor $\epsilon$ is given by Eq.
(III). The strain-induced contributions to $\eta$, $\chi$, and $\zeta$ depend
on the $\theta$ angle. Nevertheless, the major effect of QW crystallographic
orientation on the fine structure of Dirac states is captured by Eqs. (4) and
(5). Note that, in QWs with symmetric confinement, $\chi=\zeta=0$ since
$f_{3}(z)$ and $f_{4}(z)$ have opposite parities.

Finally, we estimate the coupling parameters $\gamma$, $\eta$, $\zeta$, and
$\chi$ for HgTe/Cd0.7Hg0.3Te QWs with the thickness $w=6.7$ nm using the band
parameters, elastic constants, and interface mixing from Refs. Novik2005, ;
Dantscher2015, ; Tarasenko2015, . The estimation shows that the dominant
contribution to $\gamma$ and $\eta$ of approximately $5$ meV comes from IIA
with $t_{l-h}\sim 1$. The strain contribution to $\eta$ is of the order of 1
meV for the interband deformation potential $\Xi_{cv}\approx-1$ eV and the
strain tensor components $\epsilon_{yy}=3\cdot 10^{-3}$,
$\epsilon_{zz}=-4\cdot 10^{-3}$ and $\epsilon_{yz}=10^{-3}$ calculated after
Eq. (III). The BIA contribution to $\eta$ is $\sim 0.2$ meV for the parameter
$B_{+}\approx 0.4\hbar^{2}/m_{0}$ estimated following Tab. 1 for
$E_{g}^{\prime}=4.5$ eV Lu1989 , $Q=P$, and $P^{\prime}=0.1~{}P$. The $\chi$
and $\zeta$ parameters are related to SIA and estimated for the electric field
$E_{z}=15$ kV/cm2. The contribution to $\chi$ and $\zeta$ from the cubic
warping of the energy spectrum determined by
$\gamma_{2}^{\prime}-\gamma_{3}^{\prime}$ is of the order of $0.1$ meV. The
strain contributions to $\chi$ and $\zeta$ are about $-0.1$ meV and $-0.2$
meV, respectively. Note that SIA-related coupling parameters scale with the
electric field $E_{z}$ and can be larger.

##  IV Summary

To summarize, we have studied theoretically the fine structure of Dirac states
in HgTe/CdHgTe QWs of critical and close-to-critical thicknesses. Taking into
account bulk, interface, and structure inversion asymmetry in these zinc-
blende-type QWs, we have derived the effective Hamiltonian that describes the
splitting of Dirac states in the class of $(0lh)$-oriented QWs from a unified
standpoint. The Hamiltonian contains four parameters of splitting at zero in-
plane wave vector. We have calculated these parameters in the 6-band extended
Kane model, which takes into account the lack of inversion center and the
cubic shape of the host crystal lattice, the elastic strain in the QW, and the
heavy-hole–light-hole mixing at the QW interfaces. Further, we have derived an
analytical expression for the energy spectrum of the Dirac states as a
function of the growth direction and studied how the spectrum evolves from
(001)- to (013)- and (011)-grown QWs. In general case, the spectrum is
anisotropic and in the QWs of critical thickness contains four Weyl points.
The positions of the Weyl points depend on the QW crystallographic orientation
and structure inversion asymmetry and can be controlled by an external
electric field applied along the QW normal, e.g., by gate voltage.

###### Acknowledgements.

This work was supported by the Russian Science Foundation (project
22-12-00211). G.V.B. acknowledges the support from the “BASIS” foundation.

## References

  * (1) X.-L. Qi and S. C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * (2) Z. D. Kvon, D. A. Kozlov, E. B. Olshanetsky, G. M. Gusev, N. N. Mikhailov, and S A Dvoretsky, Topological insulators based on HgTe, Phys.-Usp. 63, 629 (2020).
  * (3) B. Büttner, C. X. Liu, G. Tkachov, E. G. Novik, C. Brüne, H. Buhmann, E. M. Hankiewicz, P. Recher, B. Trauzettel, S. C. Zhang, and L. W. Molenkamp, Single valley Dirac fermions in zero-gap HgTe quantum wells, Nat. Phys. 7, 418 (2011).
  * (4) S. A. Tarasenko, M. V. Durnev, M. O. Nestoklon, E. L. Ivchenko, J.-W. Luo, and A. Zunger, Split Dirac cones in HgTe/CdTe quantum wells due to symmetry-enforced level anticrossing at interfaces, Phys. Rev. B 91, 081302(R) (2015).
  * (5) X. Dai, T. L. Hughes, X. L. Qi, Z. Fang, and S. C. Zhang, Helical edge and surface states in HgTe quantum wells and bulk insulators, Phys. Rev. B 77, 125319 (2008).
  * (6) M. König, H. Buhmann, L.W. Molenkamp, T. L. Hughes, C.-X. Liu, X. L. Qi, and S. C. Zhang, The quantum spin Hall effect: theory and experiment, J. Phys. Soc. Jpn. 77, 031007 (2008).
  * (7) R. Winkler, L. Y. Wang, Y. H. Lin, and C. S. Chu, Robust level coincidences in the subband structure of quasi-2D systems, Solid State Commun. 152, 2096 (2012).
  * (8) L. Weithofer and P. Recher, Chiral Majorana edge states in HgTe quantum wells, New J. Phys. 15, 085008 (2013).
  * (9) M. Orlita, K. Masztalerz, C. Faugeras, M. Potemski, E. G. Novik, C. Brüne, H. Buhmann, and L. W. Molenkamp, Fine structure of zero-mode Landau levels in HgTe/HgxCd1-xTe quantum wells, Phys. Rev. B 83, 115307 (2011).
  * (10) M. Zholudev, F. Teppe, M. Orlita, C. Consejo, J. Torres, N. Dyakonova, M. Czapkiewicz, J. Wróbel, G. Grabecki, N. Mikhailov, S. Dvoretskii, A. Ikonnikov, K. Spirin, V. Aleshkin, V. Gavrilenko, and W. Knap, Magnetospectroscopy of two-dimensional HgTe-based topological insulators around the critical thickness, Phys. Rev. B 86, 205420 (2012).
  * (11) P. Olbrich, C. Zoth, P. Vierling, K.-M. Dantscher, G. V. Budkin, S. A. Tarasenko, V. V. Belkov, D. A. Kozlov, Z. D. Kvon, N. N. Mikhailov, S. A. Dvoretsky, and S. D. Ganichev, Giant photocurrents in a Dirac fermion system at cyclotron resonance, Phys. Rev. B 87, 235439 (2013).
  * (12) G. M. Minkov, A.V. Germanenko, O. E. Rut, A. A. Sherstobitov, M. O. Nestoklon, S. A. Dvoretski, and N. N. Mikhailov, Spin-orbit splitting of valence and conduction bands in HgTe quantum wells near the Dirac point, Phys. Rev. B 93, 155304 (2016).
  * (13) M. V. Durnev and S. A. Tarasenko, Magnetic field effects on edge and bulk states in topological insulators based on HgTe/CdHgTe quantum wells with strong natural interface inversion asymmetry, Phys. Rev. B 93, 075434 (2016).
  * (14) E. I. Rashba, Properties of semiconductors with an extremum loop. 1. Cyclotron and combinational resonance in a magnetic field perpendicular to the plane of the loop, Sov. Phys. Solid. State 2, 1109 (1960).
  * (15) F. T. Vasko, Spin splitting in the spectrum of two-dimensional electrons due to the surface potential, JETP Lett. 30, 541 (1979).
  * (16) Y. A. Bychkov and E. I. Rashba, Properties of a 2D electron gas with lifted spectral degeneracy, JETP Lett. 39, 78 (1984).
  * (17) G. Dresselhaus, Spin-orbit coupling effects in zinc blende structures, Phys. Rev. 100, 580 (1955).
  * (18) M. I. D’yakonov and V. Y. Kachorovskii, Spin relaxation of conduction electrons in noncentrosymmetric semiconductors, Sov. Phys. Semicond. 20, 110 (1986).
  * (19) G. E. Pikus, V. A. Maruschak, and A. N. Titkov, Spin splitting of energy bands and spin relaxation of carriers in cubic III-V crystals, Sov. Phys. Semicond. 22, 115 (1988).
  * (20) E. I. Rashba and E. Y. Sherman, Spin-orbital band splitting in symmetric quantum wells, Phys. Lett. A 129, 175 (1988).
  * (21) K.-M. Dantscher, D. A. Kozlov, P. Olbrich, C. Zoth, P. Faltermeier, M. Lindner, G. V. Budkin, S. A. Tarasenko, V. V. Bel’kov, Z.D. Kvon, N. N. Mikhailov, S. A. Dvoretsky, D. Weiss, B. Jenichen, and S. D. Ganichev, Cyclotron-resonance-assisted photocurrents in surface states of a three-dimensional topological insulator based on a strained high-mobility HgTe film, Phys. Rev. B 92, 165314 (2015).
  * (22) K.-M. Dantscher, D. A. Kozlov, M. T. Scherr, S. Gebert, J. Bärenfänger, M. V. Durnev, S. A. Tarasenko, V. V. Bel’kov, N. N. Mikhailov, S. A. Dvoretsky, Z. D. Kvon, J. Ziegler, D. Weiss, and S. D. Ganichev, Photogalvanic probing of helical edge channels in two-dimensional HgTe topological insulators, Phys. Rev. B 95, 201103(R) (2017).
  * (23) G. M. Minkov, V. Ya. Aleshkin, O. E. Rut, A. A. Sherstobitov, A. V. Germanenko, S. A. Dvoretski, and N. N. Mikhailov, Valence band energy spectrum of HgTe quantum wells with an inverted band structure, Phys. Rev. B 96, 035310 (2017).
  * (24) S. Dvoretsky, N. Mikhailov, D. Ikusov, V. Kartashev, A. Kolesnikov, I. Sabinina, Y. G. Sidorov, and V. Shvets, The growth of CdTe layer on GaAs substrate by MBE, in Methods for Film Synthesis and Coating Procedures (IntechOpen, 2020).
  * (25) O. E. Raichev, Effective Hamiltonian, energy spectrum, and phase transition induced by in-plane magnetic field in symmetric HgTe quantum wells, Phys. Rev. B 85, 045310 (2012).
  * (26) G. V. Budkin and S. A. Tarasenko, Spin splitting in low-symmetry quantum wells beyond Rashba and Dresselhaus terms, Phys. Rev. B 105, L161301 (2022).
  * (27) B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Quantum spin Hall effect and topological phase transition in HgTe quantum wells, Science 314, 1757 (2006).
  * (28) L. G. Gerchikov and A. V. Subashiev, Non-monotonic behavior of the energy gap in the film made of a gapless semiconductor, Sov. Phys. Semicond. 23, 1368 (1989).
  * (29) E. G. Novik, A. Pfeuffer-Jeschke, T. Jungwirth, V. Latussek, C. R. Becker, G. Landwehr, H. Buhmann, and L. W. Molenkamp, Band structure of semimagnetic Hg1-yMnyTe quantum wells, Phys. Rev. B 72, 035321 (2005).
  * (30) R. Winkler, Spin-Orbit Coupling Effects in Two-Dimensional Electron and Hole Systems, (2003).
  * (31) G. L. Bir and G. E. Pikus, Symmetry and Strain-induced Effects in Semiconductors, (Wiley, New York, 1974).
  * (32) J.-M. Jancu, R. Scholz, E. A. de Andrada e Silva, and G. C. L. Rocca, Atomistic spin-orbit coupling and $\bm{kp}$ parameters in $\mathrm{III-V}$ semiconductors, Phys. Rev. B 72, 193201 (2005).
  * (33) M. V. Durnev, M. M. Glazov, and E. L. Ivchenko, Spin-orbit splitting of valence subbands in semiconductor nanostructures, Phys. Rev. B 89, 075430 (2014).
  * (34) I. L. Aleiner and E. L. Ivchenko, Anistropic exchange splitting in type II GaAs/AlAs superlattices, JETP Lett. 55, 692 (1992).
  * (35) E. L. Ivchenko, A. Y. Kaminski, and U. Rössler, Heavy-light hole mixing at zinc-blende (001) interfaces under normal incidence, Phys. Rev. B 54, 5852 (1996).
  * (36) Z. W. Lu, D. Singh, and H. Krakauer, Total-energy study of the equation of state of HgTe and HgSe, Phys. Rev. B 39, 10154 (1989).


# Markov cohort state-transition model: A multinomial distribution
representation

Rowan Iskandar  Center of Excellence in Decision-Analytic Modeling and Health
Economics Research,sitem-insel, Switzerland  Brown University, Providence, RI,
USA  Institute of Social and Preventive Medicine, University of Bern,
Switzerland  Cassandra Berns  Center of Excellence in Decision-Analytic
Modeling and Health Economics Research,sitem-insel, Switzerland  Brown
University, Providence, RI, USA

###### Abstract

Markov cohort state-transition models have been the standard approach for
simulating the prognosis of patients or, more generally, the life trajectories
of individuals over a time period. Current approaches for estimating the
variance of a Markov model using a Monte Carlo sampling or a master equation
representation are computationally expensive and analytically difficult to
express and solve. We introduce an alternative representation of a Markov
model in the form of a multinomial distribution. We derive this representation
from principles and then verify its veracity in a simulation exercise. This
representation provides an exact and fast approach to compute the variance and
a way to estimate transition probabilities in a Bayesian setting.

##  1 Introduction

Markov cohort state-transition models have been used in cost-effectiveness
analyses and decision-analytic modeling studies to simulate the life
trajectories of a patient or a group of patients following the receipt of a
health intervention [Beck and Pauker(1983)]. Briefly, a Markov cohort state-
transition model is a recursive matrix formula that calculates the average
number of individuals in each state using a transition probability matrix and
a initial distribution of individuals across the states. This formulation
introduced by Beck and Pauker [Beck and Pauker(1983)], and its subsequent
tutorial [Sonnenberg and Beck(1993)], however, does not fully capture the
stochastic nature of the Markov model since it only represents the average
behavior of the stochastic process [Iskandar(2018)]. A more useful quantity is
the probability distribution of all individuals at all times. To obtain this
quantity, Iskandar [Iskandar(2018)] explicated the stochastic process
underlying a Markov model in the form of a well-known time-evolution equation
of a probability function, i.e., a master equation, whose solution is the
probability distribution of interest [Iskandar(2018)]. Although the concept of
a master equation and its solution provide a complete characterization of a
Markov model, including higher-order moments, such a formulation requires some
understanding of stochastic processes. This brief report introduces an
alternative representation of a Markov model in the form of a multinomial
distribution that is intuitive and hence more accessible to practitioners. We
start with a description of a Markov cohort state-transition model. Then, we
formally derive a multinomial distribution representation. Lastly, we conduct
a numerical exercise to verify its equivalence to a canonical approach for
simulating cohort models.

##  2 Markov cohort state-transition model

A cohort of $n_{0}$ individuals where each individual (indexed by
$j=1,\ldots,n_{0}$) follows a Markov chain on a finite set of mutually
exclusive and completely exhaustive $s$ health states:
$\mathbb{S}=\\{S_{1},S_{2},\ldots,S_{s}\\}$, is defined as a stochastic
process, $\\{\mathbf{N}(t)\\}_{t\geq t_{0}}$, where $t$ and $t_{0}$ are time
and the initial time, respectively. $\mathbf{N}(t)$ represents a random vector
of individuals across $s$ states (state-configuration) :

$\displaystyle\mathbf{N}(t)=[N_{1}(t)\;N_{2}(t)\;\ldots\;N_{s}(t)],$ (1)

in which $\mathbf{N}(t)=\mathbf{n}\in\mathbb{N}^{s}$ ($\mathbb{N}$ is the set
of non-negative integers). The probability of observing a particular state-
configuration $\mathbf{N}(t)=\mathbf{n}\in\mathbb{N}^{s}$ at time $t$ or
$Pr[\mathbf{N}(t)=\mathbf{n}]\coloneqq G(\mathbf{n},t)$ (the notation
$Pr[\cdot]$ denotes the probability of observing the random event
$\\{\cdot\\}$ at time $t$) provides the complete information about the cohort
model. Our goal is to provide a multinomial distribution representation of
this $Pr(\mathbf{N}(t)=\mathbf{n})$.

##  3 Multinomial distribution representation

To derive the multinomial representation, we start by deriving the stochastic
process for each individual and aggregate the individual processes by
appealing to the multinomial distribution. Individual process. For each $j$-th
individual, we define a $1\times s$ random unit vector of health state
occupancy, $\mathbf{Y}_{j}(t)$, where the value of the $i$-th element of
$\mathbf{Y}_{j}(t)$, i.e., $Y_{j,i}(t)$, represents whether the individual is
in state $S_{i}$: $1$ (occupied) or $0$ (not occupied), i.e.,
$\mathbf{Y}_{j}(t)\in\\{0,1\\}^{s}$ and $\sum\limits_{i=1}^{s}Y_{j,i}(t)=1$
(an individual must be in one and only one state at any time). The dynamics of
$\mathbf{Y}_{j}(t)$ is governed by the allowed transitions at discrete time
steps ($t_{z}$; $z=1,\ldots$) between all pairs of states in $\mathbb{S}$ and
the intensities of these transitions. We define an $s\times s$ stochastic
matrix of transition probabilities governing the intensities of the
transitions among states in $\mathbb{S}$, i.e.,
$\mathbf{P}_{\tau}(t_{z})=\left[p_{kl}(t_{z})\right]_{1\leq k,l\leq s}$, where
$\tau$ denotes the time step (Markov cycle). Each $p_{kl}(t)$ has the usual
interpretation of the probability of an individual transitioning from states
$S_{k}$ to $S_{l}$ in one time step $\tau$ (Markov cycle) at time $t_{z}$. We
are interested in deriving the probability that the individual $j$ occupies a
particular state at any time $t_{z}$ given an initial condition. Such a
probability is also a transition probability and can be derived by iteratively
propagating the initial probability distribution of state occupancy (denoted
by the $1\times s$ vector $\mathbf{p}(t_{z})$) at the initial time, i.e.,
$\mathbf{p}_{j}(t_{0})$, forward in time using the transition probability
matrix:

$\displaystyle\mathbf{p}_{j}(t_{z})=\mathbf{p}_{j}(t_{0})\prod_{u=1}^{z-1}\mathbf{P}_{\tau}^{j}(t_{u})$
(2)

To derive the probability distribution of $\mathbf{Y}_{j}$, we note that there
are $s$ possible realizations of $\mathbf{Y}_{j}$. The probability of a
realization, i.e., the probability of an individual in some state $k$, is
equal to $[\mathbf{p}_{j}(t_{z})]_{k}$, where $[\cdot]_{k}$ denotes the $k$-th
element of a vector (e.g., $[\mathbf{y}_{j}(t_{Z})]_{k}=y_{j,k}(t_{z})$). We
write probability distribution of $\mathbf{Y}_{j}(t_{z})$ as follows:

$\displaystyle
Pr\left(\mathbf{Y}_{j}(t_{z})=\mathbf{y}_{j}(t_{z})\right)=\prod\limits_{k=1}^{s}\left(\left[\mathbf{p}_{j}(t_{z})\right]_{k}\right)^{[\mathbf{y}_{j}(t_{Z})]_{k}}$
(3)

Equation 3 has the form of a multinomial distribution, where the number of
individuals is equal to one ($n_{0}=1$). Each term in the product corresponds
to the probability of individual $j$ being in one of the $s$ states at time
$t_{z}$.

The expected value of $\mathbf{Y}_{j}(t_{z})$ is equal to
$Pr[\mathbf{Y}_{j}(t_{z})]$ since the expected value of each element of
$\mathbf{Y}_{j}(t_{z})$ is equal to the expected value of an indicator
function of whether the $j$-th individual occupies the corresponding state,
i.e.,
$E\left[Y_{j,k}(t_{z})\right]=\left[\mathbf{p}_{j}(t_{z})\right]_{k}=p_{j,k}(t_{z})$.
The covariance matrix of $\mathbf{Y}_{j}(t_{z})$, denoted by
$\Sigma_{\mathbf{Y_{j}}}(t_{z})$, is calculated using the definition of a
covariance matrix,
$\left[\Sigma_{\mathbf{Y_{j}}}(t_{z})\right]_{kl}=E[Y_{j,k}Y_{j,l}]-E[Y_{j,k}]E[Y_{j,l}]$,
and is equal to:

$\displaystyle\left[\Sigma_{\mathbf{Y_{j}}}(t_{z})\right]_{kl}=\begin{cases}p_{j,k}(t_{z})(1-p_{j,l}(t_{z}))&\mbox{
if }k=l\\\ -p_{j,k}(t_{z})p_{j,l}(t_{z})&\mbox{ if }k\neq l\end{cases}$ (4)

From individuals to a cohort. We relate $\mathbf{N}(t)$ (Equation 1) with
$\mathbf{Y}_{j}(t)$ by using the following relationship:

$\mathbf{N}(t)=\sum\limits_{j=1}^{n_{0}}\mathbf{Y}_{j}(t)$ (5)

where we write a cohort as the sum of its individual members. To derive the
probability distribution of $\mathbf{N}(t)$, we need to associate each element
in $\mathbf{N}(t)$ with the corresponding element in the sum. Each element in
the sum consists of a vector of $1$s and $0$s, i.e., $\mathbf{Y}_{j}(t)$. The
sum is constrained to be $n_{0}$ (a closed cohort). The number of possible
ways that $\mathbf{N}(t)$ can be realized, given that ordering of individuals
within group does not matter, is equal to the number of possible realizations
of the vector
$[\mathbf{Y}_{1}(t)\mathbf{Y}_{2}(t)\ldots\mathbf{Y}_{n_{0}}(t)]$. Since all
individuals are independent and follow the same transition rules and
intensities ($\mathbf{P}_{\tau}^{j}(t)=\mathbf{P}_{\tau}(t)$,
$p_{j,k}(t_{z})=p_{k}(t_{z})=$ for all $k$s, and
$\mathbf{p}_{j}(t_{0})=\mathbf{p}(t_{0})$), each realization of
$[\mathbf{Y}_{1}(t_{z})\mathbf{Y}_{2}(t_{z})\ldots\mathbf{Y}_{n_{0}}(t_{z})]$
occurs with a probability of

$\displaystyle
Pr([\mathbf{Y}_{1}(t_{z})\mathbf{Y}_{2}(t_{z})\ldots\mathbf{Y}_{n_{0}}(t_{z})])$
$\displaystyle=\prod\limits_{j=1}^{n_{0}}\prod\limits_{k=1}^{s}p_{k}(t_{z})^{y_{j,k}(t)}$
$\displaystyle=\prod\limits_{k=1}^{s}p_{k}(t_{z})^{\sum\limits_{j=1}^{n_{0}}y_{j,k}(t)}$
$\displaystyle=\prod\limits_{k=1}^{s}p_{k}(t_{z})^{n_{k}(t)}$ (6)

Multinomial distribution. The number of ways that $\mathbf{N}(t)$ can be
realized is identical to the problem of assigning $n_{0}$ of people into $s$
groups, which is given by the multinomial coefficient:
$\binom{n_{0}}{n_{1}n_{2}\ldots n_{s}}$. Putting all components together, the
probability of $\mathbf{N}(t_{z})$ at time $t_{z}$ is given by the following
multinomial distribution:

$\displaystyle Pr(\mathbf{N}(t_{z}))=\binom{n_{0}}{n_{1}n_{2}\ldots
n_{s}}\prod\limits_{k=1}^{s}p_{k}(t_{z})^{n_{k}(t_{z})}$ (7)

To derive the first two moments, we use the first two moments of the
individual process ($\mathbf{Y}_{j}(t)$) and the relationship between each
individual and a cohort of individuals (Equation 5). By linearity of
expectation, the expected value of $\mathbf{N}(t_{z})$ is equal to:

$\displaystyle E[\mathbf{N}(t_{z})]=n_{0}\mathbf{p}(t_{z})$ (8)

By independence of individuals, the variance of $\mathbf{N}(t_{z})$ is the sum
of $n_{0}$ identical variance of $\mathbf{Y_{j}}(t_{z})$ and is given by:

$\displaystyle\Sigma_{\mathbf{N}}(t_{z})=n_{0}\Sigma_{\mathbf{Y_{j}}}(t_{z})$
(9)

The covariance matrix of $\mathbf{N}(t_{z})$, denoted by
$\Sigma_{\mathbf{N}}(t_{z})$, is calculated using the definition of a
covariance matrix,
$\left[\Sigma_{\mathbf{N}}(t_{z})\right]_{uv}=E[N_{u}N_{v}]-E[N_{u}]E[N_{u}]$,
and is equal to:

$\displaystyle\left[\Sigma_{\mathbf{N}}(t_{z})\right]_{uv}=\begin{cases}n_{0}p_{u}(t_{z})(1-p_{u}(t_{z}))&\mbox{
if }v=u\\\ -n_{0}p_{u}(t_{z})p_{v}(t_{z})&\mbox{ if }v\neq u\end{cases}$ (10)

##  4 Numerical verification

We conduct a simulation study to verify whether the multinomial distribution
(Equation 7) represents the solution to a Markov cohort state-transition
model. We consider a $4$-state model with the allowed transitions and their
probabilities as follows: $Pr[S_{1}\rightarrow
S_{2}]=0.1$,$Pr[S_{1}\rightarrow S_{3}]=0.05$,$Pr[S_{1}\rightarrow
S_{4}]=0.14$, $Pr[S_{2}\rightarrow S_{3}]=0.07$,$Pr[S_{2}\rightarrow
S_{4}]=0.17$, and $Pr[S_{3}\rightarrow S_{4}]=0.11$. The simulation focuses on
comparing the mean (Equation 8) and variance (Equation 10) of the multinomial
distribution with those of a microsimulation. We simulate a cohort of 10000
individuals ($n_{0}=10000$) with $\tau=1$-year. Each individual’s life
trajectory is a realization of a Markov chain based on the given transition
probabilities. We replicate the simulation $1000$ times and calculate the
means and variances of the number of individuals across states at all times.
The results of the population trajectories across the four states are given in
Figure 1.

Figure 1: The time trajectories of individual counts in each state using two
approaches, i.e., microsimulation (blue solid lines) and multinomial
distribution (red dotted lines). For each approach, there are three lines,
i.e., mean, mean+standard deviation, mean-standard deviation. The lines are on
the top of each other indicating no difference between the two approaches.

We observe no difference between the empirical estimates of the mean and
variance from the microsimulation and the multinomial distribution. The code
for the numerical exercise is available under a GNU GPL license and can be
found at https://github.com/rowaniskandar/CM_multinomial.

##  5 Concluding remarks

This study explicates an alternative representation of a Markov cohort state-
transition model in the form of a multinomial distribution. We derive the
equivalent representation by using elementary arguments. First, the derivation
starts with specifying a Markov state-transition model for simulating an
individual following Beck and Pauker [Beck and Pauker(1983)]. We then extend
the model to a cohort of individuals by imposing independence among
individuals and uncover the multinomial distribution. We verify the first two
moments of the derived multinomial distribution of the number of individuals
across health states against using a microsimulation. The formula for the
first moment (Equation 8) is indeed the formula for simulating Markov models
as introduced by Beck and Pauker [Beck and Pauker(1983)]. To estimate the
variance, practitioners often use microsimulations, which may be
computationally expensive since we need to replicate the cohort simulation a
number of times in addition to the individual Monte Carlo runs within each
simulated cohort. In contrast, the multinomial representation provides a more
direct and non-computationally demanding approach, particularly for estimating
the variance. This approach relies only on the formula for the mean, which in
turn depends only on the initial distribution and the time-dependent or time-
invariant translation probability matrix. In addition to the computational
advantage, the multinomial representation provides a convenient way to conduct
Bayesian inference on the transition probabilities [Walley(1996)]. In a
Bayesian setting, the Markov model would be treated as the likelihood and
naturally takes the form of a multinomial distribution using the result of
this study. We can then derive the posterior distribution of the transition
probabilities by utilizing the conjugacy between the Dirichlet and multinomial
distributions. In sum, this study introduces another formulation of a well-
established methodology and reinforces the utility of Markov models further.

## References

  * [Beck and Pauker(1983)] J Robert Beck and Stephen G Pauker.  The markov process in medical prognosis.  _Med Decis Making_ , 3(4):419–458, 1983.
  * [Sonnenberg and Beck(1993)] Frank A Sonnenberg and J Robert Beck.  Markov models in medical decision making: a practical guide.  _Medical decision making_ , 13(4):322–338, 1993\.
  * [Iskandar(2018)] Rowan Iskandar.  A theoretical foundation for state-transition cohort models in health decision analysis.  _PloS one_ , 13(12):e0205543, 2018.
  * [Walley(1996)] Peter Walley.  Inferences from multinomial data: learning about a bag of marbles.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 58(1):3–34, 1996.