# Enhanced and stable spin Hall conductivity in a disordered time-reversal and
inversion symmetry broken topological insulator thin film

Siamak Pooyan  Mir Vahid Hosseini<EMAIL_ADDRESS>Department of
Physics, Faculty of Science, University of Zanjan, Zanjan 45371-38791, Iran

###### Abstract

We consider a disordered topological insulator thin film placed on the top of
a ferromagnetic insulator with a perpendicular exchange field $M$ and
subjected to a perpendicular electric field. The presence of ferromagnetic
insulator causes that bottom surface states of the topological insulator thin
film become spin polarized and the electric field provides a potential
difference $V$ between the two surface states, resulting in breaking of time-
reversal and inversion symmetry in the system. Using Kubo formalism and
employing the first Born approximation as well as the self-consistent Born
approximation, we calculate the spin Hall conductivity. We find that for small
values of $V$, a large spin conductivity can be generated through large values
of $M$ away from the charge neutrality point. But for large values of $V$, the
spin conductivity can be promoted even with small values of $M$ around the
charge neutrality point. The effect of vertex corrections and the stability of
the obtained large spin conductivity against disorders are also examined.

##  I Introduction

Topological insulators (TIs) have attracted a lot of attentions from
theoretical viewpoint and potential applications Topo1 ; Topo2 ; Topo3 ; Topo4
; Topo5 ; Topo6 . An interesting feature of these matters is that, in a
topologically nontrivial phase, edge or surface states of the system exhibit
metallic feature, due to intersecting the Fermi level, while bulk states are
an ordinary insulator at the Fermi level Topobound1 ; Topobound2 . Nontrivial
topology of bulk states can be related to the appearance of surface states,
resulting in the bulk-edge correspondence. The topology can be supported by
certain symmetries of the system. This causes surface states become robust
against perturbations respecting essential symmetries of the system
manifesting symmetry-protected topological states.

When the thickness of a 3D TI decreases so that it becomes five to ten
quintuple layers thickness1 ; thickness2 , topological states of the opposite
surfaces can hybridize together providing an interesting opportunity for
applications. Because of the hybridization, a gap opens in the surface
spectrum and carriers behave like massive Dirac fermions in a thin film of TIs
thickness3 ; thickness4 . This feature is expected to lead to the spin Hall
effect in TI thin films SHEinThinFilm .

The spin Hall effect spinHall  can occur due to a large spin-orbit coupling in
time-reversal symmetric TIs. As such, a high spin conductivity would be
expected in TIs owing to the prohibition of back scatterings. TIs have also
been studied in a magnetic field manifesting the magnetoelectric coupling
magnetoelectric  and the magnetoresistance magnetoresistance . Moreover, in a
strong magnetic field perpendicular to the surface of a TI, the charge Hall
effect can coexist with the spin Hall effect charge-spin-Hall . Recently, the
study of magnetic properties of TI thin film has become one of the hot topics
in this issue. It has been shown that TI thin films can exhibit the giant
magneto-optical Kerr effect GiantMagneto , the topological magnetoelectric
effect TopoMagEle , and the giant magnetoresistance giantMagResi  with large
spin Hall angles magnetoTranLarge .

On the other hand, disorders, including impurities and defects, are
ubiquitous, particularly, in solid state materials, e.g., TIs. A competition
between charged impurity scattering and short-range scattering in TIs has been
studied ImpScattTI . In the presence of nonmagnetic and magnetic disorders,
respectively, a large out-of-plane and an in-plane magnetoresistance have been
found on the surface of 3D TIs that is proximity-coupled to a ferromagnetic
insulator (FI) SpinTransportDisorFI-TI . In the thin TI film, surface
electrons can screen a disorder potential that is larger than the
hybridization gap ImpSceenThinTIFilm . However, the spin Hall effect is a
fragile state and its experimental measurement is a challenging task due to
ubiquitous impurities and imperfections in samples spinHall . So, it is
interesting to promote the spin Hall conductivity with a stable character
against disorders paving the way in spintronic applications.

In this respect, while in most of previous cases a magnetic exchange field is
applied to the whole system, including both surfaces of the TI thin film, it
is interesting to know what happens if the magnetization affects on only one
of the surfaces of thin film. In the present work, we calculate the dc spin
conductivity in a one-surface-magnetized TI thin film by the Kubo formalism
with random nonmagnetic potential disorders. To induce the magnetization in a
one of the surfaces of TI thin film, one can attach a ferromagnet to a surface
of the thin film. This also can provide a potential difference between the top
and bottom surfaces. We use the first Born approximation (fBA) and the self-
consistent Born approximation (SCBA) to treat nonmagnetic disorders,
respectively, analytically and numerically. We calculate the self-energies,
density of states (DOS), spin conductivity, and vertex-corrected velocity
function. Interestingly, we find that the spin Hall conductivity can be
enhanced for a large (small) exchange field and small (large) potential
difference in a large (small) chemical potential. Also, it is shown that the
promoted spin Hall conductivity remains survived at high enough impurity
strengths that is a key requirement in spintronics.

The rest of the paper is organized as follows. In Sec. II, we present model
and theory for the surface states of TI thin films with a magnetic exchange
field in a one surface and a potential difference between the two surfaces. In
Sec. III, disorders are modeled by nonmagnetic impurities on the two surfaces
and included in self-energies using the fBA and the SCBA. The effect of
disorders on DOS is also studied. In Sec. IV, the Kubo formula is used to
obtain spin Hall conductivity in the system. Section V is devoted to studying
the effect of vertex corrections. We summarize in Sec. VI.

##  II Model

Figure 1: (Color online) Schematic of a TI thin film placed on the top of a FI
in the presence of an electric field. The surfaces of the film are on the
$xy$-plane and both the polarization of FI and the electric field are along
the $z$ direction.

We consider a TI thin film, having two surface states localized at the
opposite surfaces, placed on the top of a FI in the presence of an electric
field, see Fig. 1. The low-energy effective Hamiltonian of the system can be
written as EffecHamil1 ; EffecHamil2 ; EffecHamil3

$H=rk^{2}+(H_{R}-V)\tau_{z}+\Delta\tau_{x}-M\sigma_{z}\tau_{-},$ (1)

with

$H_{R}=v_{Fk}(k_{x}\sigma_{y}-k_{y}\sigma_{x}),$ (2)

where $v_{Fk}=v_{F}(1+sk^{2})$ with $v_{F}$ being Fermi velocity and $s$
characterizes the next order correction to the Fermi velocity.
$k=\sqrt{k_{x}^{2}+k_{y}^{2}}$ with $k_{x}=k\cos\phi$ and $k_{y}=k\sin\phi$
being the in-plane momentum components. The Pauli matrices $\bm{\sigma}$ and
$\bm{\tau}$ act on the spin and the surface spaces, respectively.
$\tau_{\pm}=(\tau_{0}{\pm}\tau_{z})/2$ with $\tau_{0}$ being a unit matrix.
$r=1/(2m)$ is the inverse mass term originates from asymmetry between the
electron and hole bands EffecHamil1 . Also, $V$ is the potential difference
between the top and bottom surfaces, and $\Delta$ is the hybridization matrix
element between top and bottom surface states Hybrid . $M$ is the exchange
field due to FI being applied to the bottom surface of TI. In this system, the
inversion and time-reversal symmetry are broken by the electric and exchange
fields. These fields are assumed to be perpendicular to the surface of TI
film, i.e., along the normal direction $z$. The effect of an in-plane exchange
field SpinTransportDisorFI-TI  which shifts the position of the Dirac points
oppositely has been studied in the TI thin film leading to the giant
anisotropic magnetoresistance at low dopings AnisoMagneResTITF .

Figure 2: (Color online) Energy spectra as a function of the $k$ for (a) $M=0$
and $V=0$, (b) $M=0.4$ and $V=0$, (c) $M=0$ and $V=0.4$, (d) $M=0.4$ and
$V=0.4$, (e) $M>V$ with $M=0.4$ and $V=0.2$, and (f) $M<V$ with $M=0.2$ and
$V=0.4$. Here $\Delta=0.2$, $s=1$, and $r=1$.

Although Hamiltonian (1) cannot be diagonalized analytically, one can obtain
the energy spectra numerically. In special cases, however, the energy spectra
can be obtained as

$E_{lp}=rk^{2}+l\sqrt{V^{2}+v_{Fk}^{2}k^{2}+\Delta^{2}+2pv_{Fk}kV},$ (3)

for $M=0$ and $V\neq 0$ and

$E_{lp}=rk^{2}+l\sqrt{\frac{M^{2}}{2}+v_{Fk}^{2}k^{2}+\Delta^{2}+p\frac{\sqrt{M^{4}+4M^{2}\Delta^{2}}}{2}},$
(4)

for $M\neq 0$ and $V=0$. Here, $l=\pm$ stand for conduction and valence bands
and $p=\pm$ indicate different subbands. Throughout the paper, we take
$v_{F}/a$ as the unit of energy with $a$ being the lattice constant as the
unit of length.

The energy spectra of Hamiltonian (1) as a function of the $k$ are depicted in
Fig. 2 for different values of the $M$ and $V$. For $M=V=0$, the band
structure is gapped and the doubly degenerate conduction and valence bands,
respectively, have a minimum and a maximum at $k=0$, see Fig. 2(a). As shown
in 2(b), for $M\neq 0$ and $V=0$, the band degeneracy is lifted for the states
near $k=0$. Also, interestingly, the gap between the conduction and valence
band decreases and tends to zero when $M\gg\Delta$. This is in contrast to the
case where the exchange field $M$ is applied to both surfaces of TI thin films
EffecHamil2 . On the contrary, for $M=0$ and $V\neq 0$, the double degeneracy
of bands breaks except at state $k=0$, see Fig. 2(c). Furthermore, the minimum
(maximum) of lower conduction (upper valence) band splits into two minimums
(maximums) located at finite $k$. For finite values of both the $M$ and $V$,
the band structure is shown in Fig. 2(d), Fig. 2(e), and Fig. 2(f),
respectively, when $M=V$, $M>V$, and $M<V$. When $M\leq V$, the degeneracy
point $k=0$ in the conduction bands can be lifted [see Fig. 2(d) and Fig.
2(f)] while the valence bands remain almost intact and look similar to the
case $M=0$ and $V\neq 0$ [see Fig. 2(c)]. Moreover, when $M>V$, the two
minimums (maximums) of lower conduction (upper valence) band located at finite
$k$ begin to merge together providing a single minimum (maximum) at $k=0$ [see
Fig. 2(e)]. As a result, the band structure can be engineered via the combined
effect of both $M$ and $V$ that can be exploited in the following.

##  III Pointlike impurities

In order to investigate the effect of disorder/impurity, we consider the
potential of identical pointlike impurities in the form

$V_{imp}(\mathbf{r})=u_{0}\sigma_{0}\tau_{0}\sum_{i}\delta(\mathbf{r}-\mathbf{R}_{i}),$
(5)

where $u_{0}$ is the strength of scattering potential, $\delta(\mathbf{r})$ is
the Dirac delta function, and $\mathbf{R}_{i}$ are the coordinates of randomly
and equally distributed nonmagnetic impurities on the two surfaces of thin
film. We assume that impurity correlations are Gaussian,

$\displaystyle\langle V_{imp}(\mathbf{r})\rangle$ $\displaystyle=$
$\displaystyle 0,$ (6) $\displaystyle\langle
V_{imp}(\mathbf{r}_{1})V_{imp}(\mathbf{r}_{2})\rangle$ $\displaystyle=$
$\displaystyle n_{i}u^{2}_{0}\delta(\mathbf{r}_{1}-\mathbf{r}_{2}),$ (7)

where $n_{i}$ is the impurity density and $\langle\cdot\cdot\cdot\rangle$ is
the average over space and impurity configurations.

We define the retarded and advanced disorder averaged Green functions as

$G^{\pm}=\left[\mu-H-\Sigma^{\pm}\pm i\eta\right]^{-1},$ (8)

where $\Sigma^{\pm}$ are the self-energies, $\mu$ is the chemical potential,
and $\eta\rightarrow 0^{+}$. Using Dyson’s series Dyson , in the Born
approximation, one can expand the impurity-averaged Green’s function (8) and
gets

$G^{\pm}=\frac{G_{0}^{\pm}}{1+\Sigma^{\pm}G_{0}^{\pm}},$ (9)

where

$G^{\pm}_{0}=\left[\mu-H\pm i\eta\right]^{-1},$ (10)

are the bare retarded and advanced Green’s functions corresponding to the
Hamiltonian (1), that, in special cases, can be expressed as

$\displaystyle G_{\pm}^{0}=$
$\displaystyle\frac{1}{S_{\pm}}[\chi_{\pm}\left(g+V\left(V+2H_{R}\right)\right)-\Delta\left(g+V\left(V+2H_{R}\right)\right)\tau_{x}$
(11)
$\displaystyle+\left(\tilde{V}-\left(g-V^{2}\right)H_{R}\right)\tau_{z}],$

with

$\displaystyle S_{\pm}=$
$\displaystyle\left(\Delta^{2}-\chi_{\pm}^{2}\right)\left(\left(\Delta^{2}-\chi_{\pm}^{2}\right)+2\left(V^{2}+v_{Fk}^{2}k^{2}\right)\right)$
(12) $\displaystyle+\left(V^{2}-v_{Fk}^{2}k^{2}\right)^{2},$

for $V\neq 0$ and $M=0$ and

$\displaystyle G_{0}^{\pm}=$
$\displaystyle\frac{1}{C_{\pm}}[\chi_{\pm}g-{\Delta}\left(g+M\chi_{\pm}\sigma_{z}\right)\tau_{x}+M{\Delta}H_{y}\tau_{y}$
(13) $\displaystyle-
gH_{R}\tau_{z}-\lambda\tau_{+}+M\left(v_{Fk}^{2}k^{2}-\chi_{\pm}^{2}\right)\sigma_{z}\tau_{-}],$

with

$\displaystyle C_{\pm}=$ $\displaystyle
v_{Fk}^{2}k^{2}\left(v_{Fk}^{2}k^{2}+M^{2}-2\left(\chi_{\pm}^{2}-\Delta^{2}\right)\right)$
(14)
$\displaystyle+\left(\chi_{\pm}^{2}-\Delta^{2}\right)^{2}-M\chi_{\pm}^{2},$

for $V=0$ and $M\neq 0$ where

$\displaystyle g=v_{Fk}^{2}k^{2}+\Delta^{2}-\chi_{\pm}^{2},$ (15)
$\displaystyle\tilde{V}=V\left(V^{2}+\Delta^{2}-\chi_{\pm}^{2}-v_{Fk}^{2}k^{2}\right),$
$\displaystyle\chi_{\pm}=\left(rk^{2}-\mu\pm i\eta\right),$ $\displaystyle
H_{y}=v_{Fk}\left(k_{x}\sigma_{x}+k_{y}\sigma_{y}\right),$
$\displaystyle\lambda=M\left(M\left(\chi_{\pm}+H_{R}\right)+\Delta^{2}\sigma_{z}\right).$

The self-energies $\Sigma^{\pm}$ can also be defend as selfEnergy

$\Sigma^{\pm}=\left\langle
V_{\mathrm{imp}}G^{\pm}V_{\mathrm{imp}}\right\rangle.$ (16)

The above relation, using Eq. (5), can be read as

$\Sigma^{\pm}=n_{i}u_{0}^{2}\int\frac{d^{2}k}{(2\pi)^{2}}G^{\pm}.$ (17)

The self-energies $\Sigma^{\pm}$ can be written in the matrix structure,

$\Sigma^{\pm}=\left(\begin{array}[]{cccc}\Sigma^{\pm}_{00}&\Sigma^{\pm}_{0x}&\Sigma^{\pm}_{0y}&\Sigma^{\pm}_{0z}\\\
\Sigma^{\pm}_{x0}&\Sigma^{\pm}_{xx}&\Sigma^{\pm}_{xy}&\Sigma^{\pm}_{xz}\\\
\Sigma^{\pm}_{y0}&\Sigma^{\pm}_{yx}&\Sigma^{\pm}_{yy}&\Sigma^{\pm}_{yz}\\\
\Sigma^{\pm}_{z0}&\Sigma^{\pm}_{zx}&\Sigma^{\pm}_{zy}&\Sigma^{\pm}_{zz}\\\
\end{array}\right),$ (18)

where $\Sigma^{\pm}_{i,j}$ with $i,j=0,x,y,z$ are the elements of self-
energies. Using (9), (10), and (17), it is easy to show that the nonzero
elements are selfEnergy1 ; selfEnergy2

$\displaystyle\Sigma_{ii}^{\pm}=n_{i}u_{0}^{2}\int\frac{d^{2}k}{(2\pi)^{2}}\frac{S_{ii}^{\pm}}{D_{+}D_{-}},\quad
i=0,x,y,z,$ (19)
$\displaystyle\Sigma_{0x}^{\pm}=n_{i}u_{0}^{2}\int\frac{d^{2}k}{(2\pi)^{2}}\frac{S_{0x}^{\pm}}{D_{+}D_{-}},$
$\displaystyle\Sigma_{yz}^{\pm}=n_{i}u_{0}^{2}\int\frac{d^{2}k}{(2\pi)^{2}}\frac{S_{yz}^{\pm}}{D_{+}D_{-}},$

where $D_{\pm}=Det(G^{\pm})$ and $S^{\pm}_{i,j}$ with $i,j=0,x,y,z$ are given
in Appendix A. Here, the self-energies $\Sigma_{0x}=\Sigma_{x0}$ and
$\Sigma_{yz}=\Sigma_{zy}$. The self-energies (19) can be calculated by the
SCBA method. In the SCBA, with an initial guess of the self-energies
$\Sigma^{\pm}$, one can determine a new value for them. This process can be
done iteratively until the difference between successive values of the self-
energy becomes smaller than a desired value. On the other hand, in the fBA, we
replace $G^{\pm}$ by $G_{0}^{\pm}$ in Eq. (17). The different components of
self-energies in this method are given by

$\displaystyle\Sigma_{00}^{\pm}=n_{i}u_{0}^{2}\int\frac{dk}{2\pi}k\frac{\chi_{\pm}\left(g+V^{2}\right)+\tilde{V}+M^{2}(\chi_{\pm}-V)-M\Delta^{2}}{\Delta^{4}+2\Delta^{2}\left(V^{2}+v_{Fk}^{2}k^{2}-\chi_{\pm}^{2}\right)+\left(v_{Fk}^{2}k^{2}-(V-\chi_{\pm})^{2}\right)\left(M^{2}+v_{Fk}^{2}k^{2}-(V+\chi_{\pm})^{2}\right)},$
(20)
$\displaystyle\Sigma_{xx}^{\pm}=n_{i}u_{0}^{2}\int\frac{dk}{2\pi}k\frac{\chi_{\pm}\left(g+V^{2}\right)-\tilde{V}-\tilde{M}}{\Delta^{4}+2\Delta^{2}\left(V^{2}+v_{Fk}^{2}k^{2}-\chi_{\pm}^{2}\right)+\left(v_{Fk}^{2}k^{2}-(V-\chi_{\pm})^{2}\right)\left(M^{2}+v_{Fk}^{2}k^{2}-(V+\chi_{\pm})^{2}\right)},$
$\displaystyle\Sigma_{yy}^{\pm}=n_{i}u_{0}^{2}\int\frac{dk}{2\pi}k\frac{\chi_{\pm}\left(g+V^{2}\right)+\tilde{V}+M^{2}(\chi_{\pm}-V)+M\Delta^{2}}{\Delta^{4}+2\Delta^{2}\left(V^{2}+v_{Fk}^{2}k^{2}-\chi_{\pm}^{2}\right)+\left(v_{Fk}^{2}k^{2}-(V-\chi_{\pm})^{2}\right)\left(M^{2}+v_{Fk}^{2}k^{2}-(V+\chi_{\pm})^{2}\right)},$
$\displaystyle\Sigma_{zz}^{\pm}=n_{i}u_{0}^{2}\int\frac{dk}{2\pi}k\frac{\chi_{\pm}\left(g+V^{2}\right)-\tilde{V}+\tilde{M}}{\Delta^{4}+2\Delta^{2}\left(V^{2}+v_{Fk}^{2}k^{2}-\chi_{\pm}^{2}\right)+\left(v_{Fk}^{2}k^{2}-(V-\chi_{\pm})^{2}\right)\left(M^{2}+v_{Fk}^{2}k^{2}-(V+\chi_{\pm})^{2}\right)},$
$\displaystyle\Sigma_{0x}^{\pm}=\Sigma_{x0}=-n_{i}u_{0}^{2}\int\frac{dk}{2\pi}k\frac{\Delta\left(g+V^{2}-M(\chi_{\pm}-V)\right)}{\Delta^{4}+2\Delta^{2}\left(V^{2}+v_{Fk}^{2}k^{2}-\chi_{\pm}^{2}\right)+\left(v_{Fk}^{2}k^{2}-(V-\chi_{\pm})^{2}\right)\left(M^{2}+v_{Fk}^{2}k^{2}-(V+\chi_{\pm})^{2}\right)},$
$\displaystyle\Sigma_{yz}^{\pm}=\Sigma_{zy}=-n_{i}u_{0}^{2}\int\frac{dk}{2\pi}k\frac{\Delta\left(g+V^{2}+M(\chi_{\pm}-V)\right)}{\Delta^{4}+2\Delta^{2}\left(V^{2}+v_{Fk}^{2}k^{2}-\chi_{\pm}^{2}\right)+\left(v_{Fk}^{2}k^{2}-(V-\chi_{\pm})^{2}\right)\left(M^{2}+v_{Fk}^{2}k^{2}-(V+\chi_{\pm})^{2}\right)}.$

Here, we have introduced

$\tilde{M}=M\left(\chi_{\pm}^{2}+V^{2}-2\chi_{\pm}V-v_{Fk}^{2}k^{2}\right).$
(21)

For $r=0$, $|\mu|>\Delta$, and in the limit $M\ll 1$ and $V\ll 1$, one can
obtain short analytical approximated expressions to Eqs. (20) as

$\displaystyle\Sigma_{00}^{\pm}\approx\frac{\gamma_{b}}{\pi}\left(2V-\frac{M\Delta^{2}}{\Delta^{2}-\mu^{2}}+(V-\mu)\log\left(\frac{\mu^{2}-\Delta^{2}}{\mu^{2}-\Delta^{2}+(k_{c}v_{F})^{2}}\right)\pm
i\pi sgn(\mu)(\mu-V)\right),$ (22)
$\displaystyle\Sigma_{xx}^{\pm}\approx\frac{\gamma_{b}}{\pi}\left(-2V-\frac{M(\Delta^{2}-2V\mu)}{\Delta^{2}-\mu^{2}}-(M+V+\mu)\log\left(\frac{\mu^{2}-\Delta^{2}}{\mu^{2}-\Delta^{2}+(k_{c}v_{F})^{2}}\right)\pm
i\pi sgn(\mu)(M+V+\mu)\right),$
$\displaystyle\Sigma_{yy}^{\pm}\approx\frac{\gamma_{b}}{\pi}\left(2V+\frac{M\Delta^{2}}{\Delta^{2}-\mu^{2}}+(V-\mu)\log\left(\frac{\mu^{2}-\Delta^{2}}{\mu^{2}-\Delta^{2}+(k_{c}v_{F})^{2}}\right)\pm
i\pi sgn(\mu)(\mu-V)\right),$
$\displaystyle\Sigma_{zz}^{\pm}\approx\frac{\gamma_{b}}{\pi}\left(-2V+\frac{M(\Delta^{2}-2V\mu)}{\Delta^{2}-\mu^{2}}+(M-V-\mu)\log\left(\frac{\mu^{2}-\Delta^{2}}{\mu^{2}-\Delta^{2}+(k_{c}v_{F})^{2}}\right)\pm
i\pi sgn(\mu)(-M+V+\mu)\right),$
$\displaystyle\Sigma_{0x}^{\pm}=\Sigma_{x0}^{\pm}\approx\frac{-\gamma_{b}}{\pi}\left(\frac{\Delta
M(V-\mu)}{\left(\Delta^{2}-\mu^{2}\right)}-\Delta\log\left(\frac{\mu^{2}-\Delta^{2}}{\mu^{2}-\Delta^{2}+(k_{c}v_{F})^{2}}\right)\pm
i\pi sgn(\mu)\Delta\right),$
$\displaystyle\Sigma_{yz}^{\pm}=\Sigma_{zy}^{\pm}\approx\frac{-\gamma_{b}}{\pi}\left(-\frac{\Delta
M(V-\mu)}{\left(\Delta^{2}-\mu^{2}\right)}-\Delta\log\left(\frac{\mu^{2}-\Delta^{2}}{\mu^{2}-\Delta^{2}+(k_{c}v_{F})^{2}}\right)\pm
i\pi sgn(\mu)\Delta\right),$

where $k_{c}=\pi/a$ is the wave vector cutoff, $sgn(x)$ is the Sign function,
and $\gamma_{b}=n_{i}u_{0}^{2}/\left(4v_{F}^{2}\right)$ is the impurity
parameter. In the gap region, $|\mu|<\Delta$, due to the absence of states,
the self-energies take exponentially small values disorderedScatGap1 ;
disorderedScatGap2 . Generally, the self-energies (19) can be written as

$\Sigma^{\pm}=\Sigma^{\prime}\pm i\Gamma,$ (23)

where $\Sigma^{\prime}$ is real part and $\Gamma$ is imaginary part. This can
also be seen explicitly from Eqs. (22).

Figure 3: (Color online) Real part of the self-energies as a function of the
$\mu$ in the fBA and SCBA methods for (a) $M=0$ and $V=0$, (b) $M=0.4$ and
$V=0$, (c) $M=0$ and $V=0.4$, (d) $M=0.4$ and $V=0.4$, (e) $M>V$ with $M=0.4$
and $V=0.2$, and (f) $M<V$ with $M=0.2$ and $V=0.4$. Here $\Delta=0.2$,
$\gamma_{b}=0.2$, $s=1$, and $r=1$. Figure 4: (Color online) Imaginary part of
the self-energies as a function of the $\mu$ in the fBA and SCBA methods for
(a) $M=0$ and $V=0$, (b) $M=0.4$ and $V=0$, (c) $M=0$ and $V=0.4$, (d) $M=0.4$
and $V=0.4$, (e) $M>V$ with $M=0.4$ and $V=0.2$, and (f) $M<V$ with $M=0.2$
and $V=0.4$. Here $\Delta=0.2$, $\gamma_{b}=0.2$, $s=1$, and $r=1$.

The real and imaginary parts of self-energies versus $\mu$ with both the fBA
and the SCBA methods are evaluated in Figs. 3 and 4 for different values of
the $M$ and $V$. From both the figures, one can see that the obtained
solutions of fBA usually correspond to those of SCBA method for relativity
small $\gamma_{b}=0.2$. Both the real and imaginary parts of the self-energies
change their values at low energies while at high energies often take small
constant values. The diagonal components of self-energies are equal as well as
the off-diagonal ones in the absence of $M$ and $V$, i.e.,
$\Sigma_{00}=\Sigma_{xx}=\Sigma_{yy}=\Sigma_{zz}$ and
$\Sigma_{0x}=\Sigma_{yz}$ [see Figs. 3(a) and 4(a)]. This results from the
presence of both time-reversal and inversion symmetries. Also, for $M=0$ and
$V\neq 0$ the components are $\Sigma_{xx}=\Sigma_{zz}$,
$\Sigma_{00}=\Sigma_{yy}$, and $\Sigma_{0x}=\Sigma_{yz}$ [see Figs. 3(c) and
4(c)]. This is because of preserving time-reversal symmetry. On the other
hand, for $M\neq 0$ and $V=0$, the different components take different values
[see Figs. 3(b) and 4(b)]. This indicates that, in our model, the $M$ breaks
not only the time-reversal symmetry but also the inversion one. Furthermore,
the real and imaginary parts of $\Sigma_{xx}$ usually gets large amplitudes
compared with the others.

In the following, we study the disordered DOS disorderedDOS  using the Green’s
function technique. To do so, we start with

$D=-\frac{1}{\pi}\int\frac{kdk}{2\pi}\operatorname{Tr}\operatorname{Im}G^{+}.$
(24)

Plugging (9) and (10) into the above equation, one gets

$\displaystyle D=$
$\displaystyle\frac{-1}{\pi^{2}}\operatorname{Im}\int\frac{kdk}{D_{+}D_{-}}[M^{2}\left(\chi_{+}-V+i\left(\Gamma_{00}+\Gamma_{xx}\right)\right)$
(25) $\displaystyle+$ $\displaystyle
2\left(V^{2}\\!+\\!v_{Fk}^{2}\\!k^{2}\\!-\\!\left(\chi_{+}\\!+\\!i\left(\Gamma_{yy}\\!-\\!\Gamma_{xx}\right)\right)^{2}\\!-\\!\left(\Gamma_{0x}-i\Delta\right)^{2}\right)$
$\displaystyle\times$
$\displaystyle\left(\chi_{+}+i\left(\Gamma_{00}+\Gamma_{zz}\right)\right)].$

For $r=s=0$, in the limit $\Delta\ll 1$ the above integral can be approximated
as

$\displaystyle D$
$\displaystyle\approx\frac{-2}{\pi^{2}k_{c}^{4}v_{F}^{6}}[5\mu
V^{2}\left(2\Delta\Gamma_{0x}-3\mu\left(\Gamma_{00}+\Gamma_{xx}\right)\right)$
(26)
$\displaystyle+k_{c}^{2}v_{F}^{2}\left(2\Delta\mu\Gamma_{0x}-3\left(\Gamma_{00}+\Gamma_{xx}\right)\left(V^{2}+\mu^{2}\right)\right)$
$\displaystyle+2k_{c}^{4}v_{F}^{4}\Gamma_{00}\log\left(k_{c}\right)],$

for $M=0$ and $V\neq 0$ and

$\displaystyle D$
$\displaystyle\approx\frac{-1}{\pi^{2}k_{c}^{4}v_{F}^{6}}[M^{2}\left(3\mu^{2}\left(\tilde{\Gamma}_{1}+\tilde{\Gamma}_{2}\right)-4\mu\Delta\Gamma_{0x}\right)$
(27)
$\displaystyle+k_{c}^{2}v_{F}^{2}\left(4\mu\Delta\Gamma_{0x}+(M^{2}-6\mu^{2})\left(\tilde{\Gamma}_{1}+\tilde{\Gamma}_{2}\right)\right)$
$\displaystyle+4k_{c}^{4}v_{F}^{4}\tilde{\Gamma}_{1}\log\left(k_{c}\right)],$

for $M\neq 0$ and $V=0$, where $\tilde{\Gamma}_{1}=\Gamma_{00}+\Gamma_{yy}$
and $\tilde{\Gamma}_{2}=\Gamma_{xx}+\Gamma_{zz}$.

Figure 5: (Color online) Dependence of DOS on the $\mu$ in the fBA and SCBA
methods with $\gamma_{b}=0.2$ for different values of M and (a) $V=0$, (b)
$V=0.2$, (c) $V=0.4$, and (d) $V=0.6$.

Numerical evaluation of (25) as a function of the $\mu$ with different values
of the $M$ and $V$ is plotted in Fig. 5 for the fBA and SCBA. The DOSs of
conduction and valence bands are not equal. The DOS of conduction band is
larger than that of the valence band. As the $M$ increases, the DOS increases
and, at the same time, the gap decreases. The gap will be vanished for large
values of the $M$. For the small values of $V$ ($V=0$ and $0.2)$, in the
valence band, the DOS is almost constant and there is a maximum value of the
DOS in the conduction band. With the increases of $V$ ($V=0.4$ and $0.6$), the
coherent peaks at both the edges of gap become pronounced. Furthermore, for
small $V$, the DOS calculated by the fBA and SCBA are almost the same [see
Figs. 5(a) and 5(b)]. But, for large $V$, the SCBA results deviate from fBA
ones around the gap [see Figs. 5(c) and 5(d)].

##  IV Spin conductivity within Kubo formula

In the linear response theory, a spin current can respond to an applied
voltage through a spin conductivity that can be calculated by the Kubo formula
Dyson . The spin conductivity can be decomposed into three terms as spincond1
; spincond2 ; spincond3

$\sigma_{\alpha\beta}^{\gamma}=\sigma_{\alpha\beta}^{I\gamma}+\sigma_{\alpha\beta}^{II\gamma}+\sigma_{\alpha\beta}^{III\gamma},$
(28)

where the first and second terms include the contribution of states at the
Fermi level and the third one contains the contribution of states below the
Fermi level. Here, the indices $\alpha$ and $\beta$ stand for the coordinates
$x$ and $y$, and $\gamma$ indicates the direction of spin.

The spin conductivity due to Fermi surface states, at zero temperature, can be
written as spincond1

$\sigma_{\alpha\beta}^{I\gamma}=\frac{e}{4\pi}\int\frac{d^{2}k}{(2\pi)^{2}}\operatorname{Tr}\left[j_{\alpha}^{\gamma}G^{+}\tilde{v}_{\beta}G^{-}\right],$
(29)
$\sigma_{\alpha\beta}^{II\gamma}=-\frac{e}{8\pi}\int\frac{d^{2}k}{(2\pi)^{2}}\operatorname{Tr}\left[j_{\alpha}^{\gamma}G^{+}\tilde{v}_{\beta}G^{+}+j_{\alpha}^{\gamma}G^{-}\tilde{v}_{\beta}G^{-}\right],$
(30)

where $j_{\alpha}^{\gamma}=\left\\{\sigma_{\gamma},v_{\alpha}\right\\}/4$ is
the current operator with $v_{\alpha}=\partial H/\partial k_{\alpha}$ being
the velocity operator, $\tilde{v}_{\alpha}$ is the velocity-vertex function.
Note that the terms including only retarded or advanced Green’s functions,
i.e., $\left<j_{\alpha}^{\gamma}G^{+}v_{\beta}G^{+}\right>$ or
$\left<j_{\alpha}^{\gamma}G^{-}v_{\beta}G^{-}\right>$, can be neglected in the
weak scattering limit spincond2 ; zerothterm0 ; zerothterm1 . The spin
conductivity of filled states below the Fermi level is given by spincond1 ;
spincond2 ; spincond3

$\displaystyle\sigma_{\alpha\beta}^{III\gamma}=$
$\displaystyle\frac{e}{8\pi}\int\frac{d^{2}k}{(2\pi)^{2}}\int_{-\infty}^{\mu}f(E)dE$
(31)
$\displaystyle\times\operatorname{Tr}\left[j_{\alpha}^{\gamma}G^{+}\tilde{v}_{\beta}\frac{dG^{+}}{dE}-j_{\alpha}^{\gamma}\frac{dG^{+}}{dE}\tilde{v}_{\beta}G^{+}+\mathrm{c.c.}\right],$

where $f(E)$ is the Fermi-Dirac function, and c.c denotes complex conjugate.
We first neglect the effect of vertex corrections, in the next section, their
effect will be investigated.

In the SCBA, having obtained the self-energies self-consistently and the
impurity-averaged Green’s function, we put $G^{\pm}$ into Eqs. (29)-(31) and
calculate the spin conductivities. But in the fBA, we replace $G^{\pm}$ by
$G_{0}^{\pm}$ in Eq. (17) and substitute the self-energies (17) into Eqs.
(29)-(31) and calculate the spin conductivities. Note that for Hamiltonian
(1), although the spin conductivity tensor $\sigma_{\alpha\beta}^{III\gamma}$
can take finite values in the special case $M=V=0$ spincond3 , but it vanishes
in the cases ($M=0,V\neq 0$), ($M\neq 0,V=0$), or ($M\neq 0,V\neq 0$), due to
breaking the inversion symmetry. Then the term
$\sigma_{\alpha\beta}^{I\gamma}$ is the only considerable one that survives in
the system. To be more specific, in what follows, we focus on the
$z$-component of transverse spin conductivity. Using Eq. (29) and employing
Eqs. (9), (11), and (13), the spin conductivity $\sigma_{xy}^{Iz}$ in the
limit $\Delta\ll 1$, $r\ll 1$, and $s=0$ can be approximated as

$\displaystyle\sigma_{xy}^{Iz}$ $\displaystyle\approx\frac{4\sigma_{0}^{z}\pi
r}{k_{c}^{4}v_{F}^{6}}[V^{2}\mu\left(9\mu\left(\Gamma_{00}+\Gamma_{xx}\right)-4\Delta\Gamma_{0x}\right)$
(32)
$\displaystyle+2k_{c}^{2}v_{F}^{2}\left(\Gamma_{00}+\Gamma_{xx}\right)\left(V^{2}+\mu^{2}\right)-2k_{c}^{4}v_{F}^{4}\log\left(k_{c}\right)],$

for $M=0$ and $V\neq 0$ and

$\displaystyle\sigma_{xy}^{Iz}$ $\displaystyle\approx-\frac{2\sigma_{0}^{z}\pi
r}{k_{c}^{4}v_{F}^{6}}[M^{2}\tilde{\Gamma}_{2}\left(2k_{c}^{2}v_{F}^{2}+3\mu^{2}\right)$
(33)
$\displaystyle+2k_{c}^{2}v_{F}^{2}\left(\tilde{\Gamma}_{1}+\tilde{\Gamma}_{2}\right)\left(k_{c}^{2}v_{F}^{2}\log\left(k_{c}\right)-\mu^{2}\right)],$

for $M\neq 0$ and $V=0$, where $\sigma_{0}^{z}=e/16\pi^{3}$.

Figure 6: (Color online) Density plot of the spin conductivity
$\sigma_{xy}^{Iz}$ in the fBA method as functions of the $M$ and $\mu$ with
the disorder parameter $\gamma_{b}=0.2$ for (a) $V=0$, (b) $V=0.2$, (c)
$V=0.4$, (d) $V=0.6$, (e) $V=0.8$, and (f) $V=1$.

In Fig. 6, the density plot of the spin conductivity $\sigma_{xy}^{Iz}$ as
functions of the $M$ and $\mu$ is depicted in the fBA for various values of
the $V$. All plots show that the spin conductivity via the conduction band
states can have larger values than those of the valence band arising from the
large DOS in the conduction band and more dispersive feature of the conduction
band. For $V=0$, as shown in Fig. 6(a), there is a global gap around $\mu=0$
with the largest bandgap width at $M=0$. Also, a large spin conductivity
occurs for large $|M|$ in the conduction band away the charge neutrality
point, $\mu=0$. Interestingly, as shown in Fig. 6(b) for $V=0.2$ a region with
large values of the spin conductivity splits into two parts so that one part
shifts towards $\mu=0$ and the other one shifts oppositely as $V$ increases,
see Figs. 6(c)-6(f). At the same time, a considerable spin conductivity
emerges in the top of valence band for small $M$. As a result, this changes
the global gap into the partial gap and a large spin conductivity can take
place at low dopings with small amplitudes of the $M$.

Figure 7: (Color online) Density plot of the spin conductivity
$\sigma_{xy}^{Iz}$ in the fBA method as functions of the $V$ and $\mu$ with
the disorder parameter $\gamma_{b}=0.2$ for (a) $M=0$, (b) $M=0.2$, (c)
$M=0.4$, (d) $M=0.6$, (e) $M=0.8$, and (f) $M=1$.

The spin conductivity density plot as functions of $V$ and $\mu$ is depicted
for various values of the $M$ in Fig. 7. For $M=0$, as can be seen from Fig.
7(a), there is a partial gap between conduction and valence bands and the spin
conductivity $\sigma_{xy}^{Iz}$ takes moderate values near the edges of the
gap with smallest values near $V=0$. As $M$ increases, the gap decreases, in
particular, for small $V$ and the spin conductivity $\sigma_{xy}^{Iz}$ near
the lower edge of the gap begin to vanish [see Figs. 7(b)-7(f)]. Moreover, the
spin conductivity $\sigma_{xy}^{Iz}$ increases not only near the upper edge of
gap but also, interestingly, at high dopings for small values of the $V$.

Figure 8: (Color online) The spin conductivity $\sigma_{xy}^{Iz}$ in the fBA
and SCBA methods as a function of the $\mu$ for different values of M with (a)
$V=0$, (b) $V=0.2$, (c) $V=0.4$, and (d) $V=0.6$. Here, the disorder parameter
$\gamma_{b}=0.2$.

In Fig. 8, the spin conductivity $\sigma_{xy}^{Iz}$ is depicted in terms of
the $\mu$ for various values of the $M$ and $V$ in both the fBA and SCBA
methods. For $V=0$, see Fig. 8(a), with $M=0$, for small impurity parameter
$\gamma_{b}=0.2$ the fBA and SCBA have the same diagrams. Moreover, with the
increase of $M$, the solutions of SCBA slightly deviate from those of the fBA.
This is because of lifting the degeneracy at low energies, as already
discussed. As $V$ increases, asymmetric coherent peaks at the edges of band
gap appear, see Figs. 8(b)-8(d). One also realizes that the deviation between
both solutions becomes more pronounced at low energies. This can be attributed
to lifting the degeneracy of more states when the $V$ turns on providing more
scattering state. In particular, the difference of both the solutions in the
conduction band is considerable than that in the valence band. Because, the
Fermi surfaces of conduction band is smaller than those of valence band
increasing more scattering processes. Furthermore, at high energies, the
difference between the solutions of the fBA and SCBA decreases, because the
impurities affect on low energy states.

Figure 9: (Color online) Density plot of the spin conductivity
$\sigma_{xy}^{Iz}$ in the SCBA method as functions of the disorder parameter
$\gamma_{b}$ and $\mu$ for (a) $M=0$ and $V=0$, (b) $M>V$ with $M=1.5$ and
$V=0.2$, and (c) $M<V$ with $M=0.2$ and $V=1$. Here $\Delta=0.2$, s=1, and
r=1.

The dependence of the spin conductivity $\sigma_{xy}^{Iz}$ on the disorder
parameter $\gamma_{b}$ and $\mu$ in the SCBA method is shown in Fig. 9 for
different values of the $M$ and $V$. In the absence of magnetic exchange field
and potential difference, i.e., $M=0$ and $V=0$, the small values of spin
conductivity spoil as $\gamma_{b}$ increases [see Fig. 9(a)]. But,
interestingly, for either $M>V$ with large $\mu$ [see Fig. 9(b)] or $M<V$ with
small $\mu$ [see Fig. 9(c)], the obtained high spin conductivity can sustain
even in the large values of $\gamma_{b}$ providing a stable spin Hall
conductivity. Note that we have calculated the $\sigma_{xy}^{IIz}$ versus the
disorder parameter $\gamma_{b}$ and $\mu$ (not shown). Its values is
negligible compared to the $\sigma_{xy}^{Iz}$. So, the patterns of Fig. 9 do
not change by including the $\sigma_{xy}^{IIz}$.

##  V Vertex corrections

The velocity-vertex function $\tilde{v}_{\alpha}$ satisfies the self-
consistent equation zerothterm0 ; zerothterm1

$\tilde{v}_{\alpha}=v_{\alpha}+n_{i}u_{0}^{2}\int\frac{d^{2}k}{(2\pi)^{2}}G^{R}\tilde{v}_{\alpha}G^{A}.$
(34)

To ensure about the true form of vertex function matrix structure, we
calculate the first-order correction of velocity-vertex function
$\tilde{v}_{\alpha}$. Using iteration, the first-order correction to velocity-
vertex function reads

$\tilde{v}_{\alpha}^{(1)}=n_{i}u_{0}^{2}\int\frac{d^{2}k}{(2\pi)^{2}}G^{R}v_{\alpha}G^{A},$
(35)

in the case $M=0$ and $V\ll 1$ for $r=s=0$, it can be obtained as

$\displaystyle\tilde{v}_{x}^{(1)}=\frac{\gamma_{b}v_{F}}{\pi}[$
$\displaystyle\left(\frac{8\Gamma_{xx}\Gamma_{0x}\Delta\mu}{3\left(\Delta^{2}-\mu^{2}\right)^{2}}-1\right)\sigma_{y}\tau_{z}$
(36)
$\displaystyle+\frac{2\left(\Gamma_{0x}\mu+\Gamma_{00}\Delta\right)}{3\left(\Delta^{2}-\mu^{2}\right)}\sigma_{y}\tau_{y}$
$\displaystyle+\frac{4V\Gamma_{00}\Gamma_{0x}\left(\Delta\sigma_{y}\tau_{0}-\mu\sigma_{y}\tau_{x}\right)}{3\left(\Delta^{2}-\mu^{2}\right)^{2}}],$

$\displaystyle\tilde{v}_{y}^{(1)}=-\frac{\gamma_{b}v_{F}}{\pi}[$
$\displaystyle\left(\frac{8\Gamma_{xx}\Gamma_{0x}\Delta\mu}{3\left(\Delta^{2}-\mu^{2}\right)^{2}}-1\right)\sigma_{x}\tau_{z}$
(37)
$\displaystyle+\frac{2\left(\Gamma_{0x}\mu+\Gamma_{00}\Delta\right)}{3\left(\Delta^{2}-\mu^{2}\right)}\sigma_{x}\tau_{y}$
$\displaystyle+\frac{4V\Gamma_{00}\Gamma_{0x}\left(\Delta\sigma_{x}\tau_{0}-\mu\sigma_{x}\tau_{x}\right)}{3\left(\Delta^{2}-\mu^{2}\right)^{2}}].$

In the equations above, for weak scattering limit the first terms will be
dominated, so one can expand the velocity vertex as

$\displaystyle\tilde{v}_{x}=v_{x}+\delta\sigma_{y}\tau_{z},$ (38)
$\displaystyle\tilde{v}_{y}=v_{y}-\delta\sigma_{x}\tau_{z}.$

Plugging Eqs. (38) into Eq. (34), one can obtain the correction $\delta$ to
the velocities numerically. Note that for the case $V=0$ and $M\ll 1$, we
examined that the matrix structure of Eqs. (38) remains the same.

Figure 10: (Color online) The vertex correction $\delta$ in the fBA and SCBA
methods as a function of the $\mu$ for different values of $M$ with (a) $V=0$,
(b) $V=0.2$, (c) $V=0.4$, and (d) $V=0.6$. Here, the disorder parameter is
$\gamma_{b}=0.2$.

In Fig. 10, the velocity correction $\delta$ is evaluated as a function of the
$\mu$ in both the fBA and SCBA methods with different values of the $M$ and
$V$. At high energies, the velocity correction $\delta$ has larger values than
those for low energies almost independent of the $M$. Moreover, at low
energies, there is a dip such that the wide of dip becomes narrowed as $M$
increases. For $V=0$, the dips are centered at $\mu=0$ as shown in Fig. 10(a).
Also, as $V$ increases the center of dips deviates from zero depending on the
$M$ and, at the same time, the $\delta$ exhibits a strong sudden change in
$\mu>0$ as shown in Figs. 10(b)-10(d).

##  VI Summary

We considered a disordered 2D topological insulator thin film, having two
surface states, subjected to a magnetic exchange field and a potential
difference. The magnetic exchange field is applied perpendicular to a one of
the surface states and the potential difference is between the two surfaces.
Using the fBA and SCBA, we calculated both self-energies and DOS analytically
and numerically. It is found that, in the model, a large spin Hall
conductivity can be reached for either large magnetic exchange fields and
small potential differences at high dopings or small magnetic exchange fields
and large potential differences at low dopings. The promoted spin Hall effect
is also examined with respect to the impurity parameter exhibiting a good
stability against scattering of charge impurities. Also, the effect of vertex
corrections is investigated changing the velocity at low energies.

##  Appendix A Some notations

We have introduced $S^{\pm}_{i,j}$ in Eq. (19) of the main text as

$\displaystyle S_{00}^{\pm}=$ $\displaystyle
V\left(V^{2}-\left(v_{Fk}^{2}k^{2}+\left(\chi_{\pm}+\Sigma_{00}\right)^{2}-\left(\Delta+\Sigma_{0x}\right)^{2}\right)\right)-M\left(\left(V-\chi_{\pm}-\Sigma_{yy}\right)+\left(\Delta+\Sigma_{0x}\right)^{2}\right)$
$\displaystyle+\left(\chi_{\pm}+\Sigma_{zz}\right)\left(V^{2}-\left(-v_{Fk}^{2}k^{2}+\left(\chi_{\pm}+\Sigma_{xx}\right)^{2}-\left(\Delta+\Sigma_{yz}\right)^{2}\right)\right),$
$\displaystyle S_{xx}^{\pm}=$
$\displaystyle-V\left(V^{2}-\left(v_{Fk}^{2}k^{2}+\left(\chi_{\pm}+\Sigma_{00}\right)^{2}-\left(\Delta+\Sigma_{0x}\right)^{2}\right)\right)+M\left(v_{Fk}^{2}k^{2}+\left(V-\chi_{\pm}-\Sigma_{yy}\right)\right)$
$\displaystyle+\left(\chi_{\pm}+\Sigma_{zz}\right)\left(V^{2}-\left(-v_{Fk}^{2}k^{2}+\left(\chi_{\pm}+\Sigma_{yy}\right)^{2}-\left(\Delta+\Sigma_{yz}\right)^{2}\right)\right),$
$\displaystyle S_{yy}^{\pm}=$ $\displaystyle
V\left(V^{2}-\left(v_{Fk}^{2}k^{2}+\left(\chi_{\pm}+\Sigma_{00}\right)^{2}-\left(\Delta+\Sigma_{0x}\right)^{2}\right)\right)-M\left(\left(V-\chi_{\pm}-\Sigma_{yy}\right)-\left(\Delta+\Sigma_{0x}\right)^{2}\right)$
$\displaystyle+\left(\chi_{\pm}+\Sigma_{zz}\right)\left(V^{2}-\left(-v_{Fk}^{2}k^{2}+\left(\chi_{\pm}+\Sigma_{xx}\right)^{2}-\left(\Delta+\Sigma_{yz}\right)^{2}\right)\right),$
$\displaystyle S_{zz}^{\pm}=$
$\displaystyle-V\left(V^{2}-\left(v_{Fk}^{2}k^{2}+\left(\chi_{\pm}+\Sigma_{00}\right)^{2}-\left(\Delta+\Sigma_{0x}\right)^{2}\right)\right)-M\left(v_{Fk}^{2}k^{2}-\left(V-\chi_{\pm}-\Sigma_{yy}\right)\right)$
$\displaystyle+\left(\chi_{\pm}+\Sigma_{zz}\right)\left(V^{2}-\left(-v_{Fk}^{2}k^{2}+\left(\chi_{\pm}+\Sigma_{yy}\right)^{2}-\left(\Delta+\Sigma_{yz}\right)^{2}\right)\right),$
$\displaystyle S_{0x}^{\pm}=$
$\displaystyle-\left(\Delta+\Sigma_{0x}\right)\left(v_{Fk}^{2}k^{2}+V^{2}-M\left(-V+\chi_{\pm}+\Sigma_{00}\right)-\left(\chi_{\pm}+\Sigma_{xx}\right)^{2}+\left(\Delta+\Sigma_{0x}\right)^{2}\right),$
$\displaystyle S_{yz}^{\pm}=$
$\displaystyle-\left(\Delta+\Sigma_{yz}\right)\left(v_{Fk}^{2}k^{2}+V^{2}+M\left(-V+\chi_{\pm}+\Sigma_{yy}\right)-\left(\chi_{\pm}+\Sigma_{zz}\right)^{2}+\left(\Delta+\Sigma_{yz}\right)^{2}\right).$

## References

  * (1) M.Z. Hasan and C.L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * (2) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * (3) Y. Ando, J. Phys. Soc. Jpn. 82, 102001 (2013).
  * (4) J.C.Y. Teo and C.L. Kane, Phys. Rev. B 82, 115120 (2010).
  * (5) C.-K. Chiu, J.C.Y. Teo, A.P. Schnyder, and S. Ryu, Rev. Mod. Phys. 88, 035005 (2016).
  * (6) A.P. Schnyder, S. Ryu, A. Furusaki, and A.W.W. Ludwig, Phys. Rev. B 78, 195125 (2008).
  * (7) J.E. Moore, Nature (London) 464, 194 (2010).
  * (8) X.-L. Qi and S.-C. Zhang, Phys. Today 63, 33 (2010).
  * (9) G. Zhang, H. Qin, J. Teng, J. Guo, Q. Guo, X. Dai, Z. Fang, and K. Wu, Appl. Phys. Lett. 95, 053114 (2009).
  * (10) H. Peng, K. Lai, D. Kong, S. Meister, Y. Chen, X. Qi, S. Zhang, Z. Shen, and Y. Cui, Nat. Mater. 9, 225 (2010).
  * (11) Y. Zhang, K. He, C.-Z. Chang, C.-L. Song, L.-L. Wang, X. Chen, J.-F. Jia, Z. Fang, X. Dai, W.-Y. Shan et al., Nat. Phys. 6, 584 (2010).
  * (12) D. Kim, P. Syers, N.P. Butch, J. Paglione, and M.S. Fuhrer, Nat. Commun. 4, 2040 (2013).
  * (13) S.K. Chong, L. Liu, T.D. Sparks, F. Liu, and V.V. Deshpande, arXiv:2004.04870.
  * (14) J. Sinova, S. O. Valenzuela, J. Wunderlich, C. Back, and T. Jungwirth, Rev. Mod. Phys. 87, 1213 (2015).
  * (15) A. LaForge, A. Frenzel, B. Pursley, T. Lin, X. Liu, J. Shi, and D. Basov, Phys. Rev. B 81, 125120 (2010).
  * (16) A. Taskin, K. Segawa, and Y. Ando, Phys. Rev. B 82, 121302R (2010).
  * (17) D.-H. Lee, Phys. Rev. Lett. 103, 196804 (2009).
  * (18) W.-K. Tse and A.H. MacDonald, Phys. Rev. Lett. 105, 057401 (2010).
  * (19) T. Morimoto, A. Furusaki, and N. Nagaosa, Phys. Rev. B 92, 085113 (2015).
  * (20) R.-L. Chu, J. Li, J.K. Jain, and S.-Q. Shen, Phys. Rev. B 80, 081102 (2009).
  * (21) Y. Fan, P. Upadhyaya, X. Kou, M. Lang, S. Takei, Z. Wang, J. Tang, L. He, L.-T. Chang, M. Montazeri, G. Yu,W. Jiang, T. Nie, R.N. Schwartz, Y. Tserkovnyak, and K.L. Wang, Nat. Mater. 13, 699 (2014).
  * (22) D. Culcer, E.H. Hwang, T.D. Stanescu, and S. Das Sarma, Phys. Rev. B 82, 155457 (2010).
  * (23) T. Chiba, S. Takahashi, and G.E.W. Bauer, Phys. Rev. B 95, 094428 (2017).
  * (24) Y. Huang and B.I. Shklovskii, Phys. Rev. B 103, 165409 (2021).
  * (25) C.-X. Liu, X.-L. Qi, H.J. Zhang, X. Dai, Z. Fang, and S.-C. Zhang, Phys. Rev. B 82, 045122 (2010).
  * (26) W. Luo and X.-L. Qi, Phys. Rev. B 87, 085431 (2013).
  * (27) A.T. Lee, M.J. Han, K. Park, Phys. Rev. B 90, 155103 (2014).
  * (28) A.A. Zyuzin, M.D. Hook, and A.A. Burkov, Phys. Rev. B 83, 245428 (2011).
  * (29) A. G. Moghaddam, A. Qaiumzadeh, A. Dyrdał, and J. Berakdar, Phys. Rev. Lett. 125, 196801 (2020).
  * (30) G. D. Mahan, Many-Particle Physics, (Kluwer Academic, New York, 2000).
  * (31) A. Altland and B.D. Simons, Condensed Matter Field theory (Cambridge, 2010).
  * (32) G. Rickaysen, Green’s Functions and Condensed Matter (Academic, New York, 1980).
  * (33) S. Doniach and E. H. Sondheimer, Green’s Functions for Solid State Physicists (Imperial College, London, 1998).
  * (34) P.M. Ostrovsky, I.V. Gornyi, and A.D. Mirlin, Phys. Rev. B 74, 235443 (2006).
  * (35) M.V. Hosseini, J. Phys. Condens. Matter 29, 315502 (2017).
  * (36) B.Y.-K. Hu, E.H. Hwang, and S. Das Sarma, Phys. Rev. B 78, 165411 (2008).
  * (37) M.-F. Yang and M.-C. Chang, Phys. Rev. B 73, 073304 (2006).
  * (38) D. Kodderitzsch, K. Chadova, and H. Ebert, Phys. Rev. B 92, 184415 (2015).
  * (39) R.S. Akzyanov, Phys. Rev. B 100, 045403 (2019).
  * (40) N.A. Sinitsyn, J.E. Hill, H. Min, J. Sinova, and A.H. MacDonald, Phys. Rev. Lett. 97, 106804 (2006).
  * (41) N.A. Sinitsyn, A.H. MacDonald, T. Jungwirth, V.K. Dugaev, and J. Sinova, Phys. Rev. B 75, 045315 (2007).


# Maximum likelihood (ML) estimators for scaled mutation parameters with a
strand symmetric mutation model in equilibrium

Claus Vogl<EMAIL_ADDRESS>Lynette Caitlin Mikula  lcm29@st-
andrews.ac.uk Department of Biomedical Sciences, Vetmeduni Vienna,
Veterinärplatz 1, A-1210 Wien, Austria  Centre for Biological Diversity,
School of Biology, University of St. Andrews, St Andrews KY16 9TH, UK

###### Abstract

With the multiallelic parent-independent mutation-drift model, the equilibrium
proportions of alleles are known to be Dirichlet distributed. A special case
is the biallelic model, in which the proportions are beta distributed. A
sample taken from these models is then Dirichlet-multinomially or beta-
binomially distributed, respectively. Maximum likelihood (ML) estimators for
the mutation parameters of the biallelic parent-independent mutation model are
available via an expectation maximization algorithm. Assuming small scaled
mutation rates, the distribution of a sample of size $M$ can be expanded in a
Taylor series of first order. Then the ML estimators for the two parameters in
the biallelic model can be expressed using the site frequency spectrum. In
this article, we go beyond parent-independent mutation and analyse a strand-
symmetric mutation model with six scaled mutation parameters that deviates
from parent independent mutation and, generally, from detailed balance. We
derive ML estimators for these six parameters assuming mutation-drift
equilibrium and small scaled mutation rates. This is the first time that ML
estimators are provided for a mutation model more complex than parent-
independent mutation.

###### keywords:

strand-symmetric mutation, mutation-drift model , scaled mutation parameters ,
maximum likelihood inference , expectation-maximization algorithm.

††journal: Theoretical Population Biology

##  1 Introduction

With the parent-independent mutation-drift model, theoretical results for the
estimation of mutation parameters are available: The proportions of alleles of
a biallelic model in equilibrium are beta distributed [18]. Data for the
inference of parameters are usually in the form of site frequency spectrum
(also called allele frequency spectrum) data taken from one population. Such
data are a vector of allele frequencies for samples of $M$ haploid individuals
at a total of $L$ sites (or loci) whereby there are $L_{0},L_{1},\dots,L_{M}$
sites with $y=0,1,\dots,M$ alleles of the focal type. Given a sample of size
$M$ and a binomial sampling distribution conditional on the allele proportion,
one obtains a beta-binomial distribution for the allele frequencies. Maximum
likelihood (ML) estimators can be constructed from the beta-binomial
distribution via an expectation maximization algorithm, where a root of a
polynomial of the order of the sample size needs to be evaluated at each
iteration [16]. With a multiallelic parent-independent mutation model, a
Dirichlet-multinomial distribution for the allele frequencies follows
analogously.

In the limit of small scaled mutation rates, it is convenient to reparametrize
the two parameters of the beta distribution for equilibrium allele proportions
with the mutation bias $\alpha$ and the overall mutation rate $\theta$.
Expanding the beta-binomial distribution in a Taylor series up to first order
of the scaled mutation rate $\theta$ then results in a relatively simple
equilibrium sample distribution [16]:

${\Pr}(y{\,|\,}M,\alpha,\theta\to
0)=\begin{cases}(1-\alpha)\big{(}1-\alpha\theta
H_{M-1}\big{)}+O(\theta^{2})&\text{for $y=0$,}\\\
\alpha(1-\alpha)\theta\,\frac{M}{y(M-y)}+O(\theta^{2})&\text{for $1\leq
y\leq(M-1)$,}\\\ \alpha\big{(}1-(1-\alpha)\theta
H_{M-1}\big{)}+O(\theta^{2})&\text{for $y=M$}\,.\end{cases}$ (1)

whereby $H_{m}$ denotes the $m$th harmonic number. Note that since $H_{M-1}$
approaches infinity logarithmically, the sample size must stay below a limit
of approximately $M<e^{\frac{1}{\max(1-\alpha,\alpha)\theta}}$. Furthermore,
simulations have shown that the mutation rates for the first order Taylor
series expansion must be $2\alpha(1-\alpha)\theta\leq 0.025$ [17]; this upper
boundary for the scaled mutation rate was further substantiated by Schrempf
and Hobolth [14]. It is higher than the scaled mutation rate of most
eukaryotes [12]. With this approximation, simple ML estimators given site-
frequency spectrum data are obtained:

$\hat{\alpha}=\frac{L_{M}+L_{p}}{L}$ (2)

and

$\hat{\theta}=\frac{1}{2\alpha(1-\alpha)}\frac{L_{p}}{H_{M-1}L}\,,$ (3)

respectively [16].

Vogl and Clemente [17] proposed a Moran model, in which mutations are
restricted to monomorphic states. In this so-called boundary-mutation Moran
model, the same equilibrium sample distribution is obtained as with the first
order Taylor series expansion in formula (1). Schrempf and Hobolth [14] were
able to explicitly derive the approximate neutral multiallelic stationary
distribution of allele frequencies as the equilibrium state of a discrete
boundary mutation Moran model with a completely general mutation model.
Starting from a Wright-Fisher model and passing to the diffusion limit, Burden
and Griffiths [1] showed that a Taylor series expansion in the scaled mutation
rate can be used to obtain the stationary distribution for mutation models
more complex than the parent-independent mutation model. This is remarkable
because the equilibrium distribution of such models is not known. Below, we
provide an alternative proof.

Our main results, however, are ML estimators for the scaled mutation
parameters of a strand-symmetric mutation model. The two strands of the DNA
double helix are held together in anti-parallel orientation by hydrogen bonds.
Assuming symmetry between the two strands, the twelve mutation rate parameters
between the four bases adenine, thymine, cytosine, and guanine ($A$, $T$, $C$,
and $G$) are reduced to six because the transition rates of complementary
bases are identical, e.g., the transition rate from $A$ to $C$ is identical to
that from $T$ to $G$. For these six parameters, we provide estimators using
expectation-maximization (EM) algorithms that are guaranteed to converge to
the global maximum.

##  2 The strand-symmetric mutation-drift model

At each site, we assume $K=4$ alleles corresponding to the four bases. The
evolutionary process is described by a continuous time Markov chain with
instantaneous transition rate matrix $\mathbf{Q}$ of $dim(\mathbf{Q})=4$x$4$
and corresponding forward transition probability matrix $\mathbf{P}$. Under
the most general mutation model, the off-diagonal entries of $\mathbf{Q}$ are
the strictly positive transition rates among the four bases $A$, $T$, $G$, and
$C$:

$\mathbf{Q}=\tiny\begin{pmatrix}-a_{12}-a_{13}-a_{14}&a_{12}&a_{13}&a_{14}\\\
a_{21}&-a_{21}-a_{23}-a_{24}&a_{23}&a_{24}\\\
a_{31}&a_{32}&-a_{31}-a_{32}-a_{34}&a_{34}\\\
a_{41}&a_{42}&a_{43}&-a_{41}-a_{42}-a_{43}\\\ \end{pmatrix}\,.$ (4)

Thus the general transition rate matrix has twelve parameters. Finding the
unique stationary state for a general transition rate matrix $\mathbf{Q}$
requires scaling the individual entries to ensure they are $\leq 1$. Applying
the Perron-Frobenius theorem to the stochastic matrix given by the sum of this
scaled transition rate matrix and the identity matrix yields a dominant
eigenvalue of zero. All further eigenvalues have a negative real component.
Provided the process is irreducible, the eigenvector associated with the zero-
valued eigenvalue corresponds to the unique stationary distribution of the
process.

With the strand-symmetric mutation model, similar conclusions can be drawn
without invoking the Perron-Frobenius theorem. The number of parameters of
$\mathbf{Q}$ is reduced to six because the transition rates between bases are
identical on complementary DNA strands. The system is thereby simplified and
$\mathbf{Q}$ can be written as in Lobry [11] (note that our $\mathbf{Q}$ is
the transpose of Lobry’s):

$\mathbf{Q}=\begin{pmatrix}-a-c-e&a&c&e\\\ a&-a-c-e&e&c\\\ b&d&-b-d-f&f\\\
d&b&f&-d-b-f\end{pmatrix}$ (5)

. Again $a,..,f$ are strictly positive and the rows of $\mathbf{Q}$ sum to
zero. The leading eigenvalue of $\mathbf{Q}$ is $\lambda_{1}=0$ with the
associated normed eigenvector

$v_{1}=\begin{pmatrix}\frac{d+b}{2(b+c+d+e)}&\frac{d+b}{2(b+c+d+e)}&\frac{c+e}{2(b+c+d+e)}&\frac{c+e}{2(b+c+d+e)}\end{pmatrix}\,.$
(6)

This vector corresponds to the stationary distribution [11]. Clearly, the
probabilities of the four bases at stationarity follow Chargaff’s second
parity rule: For each of the two DNA strands, the expectations of the
proportions of $A$ and $T$ are identical, as are those of $C$ and $G$.
Chargaff’s second parity rule has been shown to hold for all types of double-
stranded DNA except organellar DNA [13].

The second eigenvalue [11] is $\lambda_{2}=-(b+d+c+e)$ with associated
eigenvector

$v_{2}=\begin{pmatrix}1&1&-1&-1\end{pmatrix}\,.$

This eigenvector contrasts the sum of $A$ and $T$ with the sum of $C$ and $G$.

The next two eigenvalues ($\lambda_{3},\lambda_{4}$) may have complex
components that introduce a probability flow into the system [11]. However, it
must be noted that if all transition rates are similar in magnitude, this
potentially imaginary term fluctuates around zero. The stationary distribution
given by the leading eigenvector is stable despite the complex eigenvalues
because the real components of $\lambda_{2,r;3,r;4,r}$ are strictly negative
[11]. Provided the individual transition rates are indeed similar in
magnitude, the real components of $\lambda_{3,r}$ and $\lambda_{4,r}$ are
equal to $-a-f+\frac{\lambda_{2}}{2}$ and are of the same order of magnitude
as $\lambda_{2}$ (The sign in front of the fraction is incorrect in Lobry’s
listing of eigenvalues; his further arguments indicate that this is merely a
typo). Thus the evolutionary process converges exponentially towards a
stationary state at a rate that depends on $\lambda_{2,3,4}$.

On the basis of Lobry’s work, we can be confident of a unique stationary
distribution given the transition rates among alleles in the strand-symmetric
mutation model except in degenerate and thus biologically meaningless cases.

##  3 The stationary distribution of the strand-symmetric mutation-drift model
with small scaled mutation rates

In this section, we derive the stationary distribution for a sample of size
$M$ taken from the general multiallelic boundary-mutation Moran model in the
limit of small scaled overall mutation rates. It is assumed that the mutation
matrix gives rise to a unique stationary distribution. This is a similar
approach to that of Burden and Griffiths [1] and recovers the full stationary
distribution for the proportions of alleles (their formulas 7-9 in Theorem 1)
[see also 14].

We reparametrize the general mutation matrix as follows:

$\mathbf{Q}=\theta\cdot\begin{pmatrix}-\alpha_{1}&\alpha_{12}&\alpha_{13}&\alpha_{14}\\\
\alpha_{21}&-\alpha_{2}&\alpha_{23}&\alpha_{24}\\\
\alpha_{31}&\alpha_{32}&-\alpha_{3}&\alpha_{34}\\\
\alpha_{41}&\alpha_{42}&\alpha_{43}&-\alpha_{4}\\\ \end{pmatrix}\,,$ (7)

where $\alpha_{k}=\sum_{l=1,l\neq k}^{K}\alpha_{kl}$.

Assuming a unique stationary distribution for $\mathbf{Q}$ and recalling the
duality between the Moran model and Kingman’s coalescent [6, chapt. 2], we can
use a sampling algorithm proposed by Stephens and Donnelly [15] to build a
genealogical realization of a Moran model of sample size $M$ forwards in time.
This sampling algorithm is essentially an urn sampling process. Hoppe [9]
first drew the connection between Polya-like urn processes and Kingman’s
coalescent; Donnelly and Kurtz [5] proposed an algorithm by which one can
obtain the realization of a sample from a particle Fleming-Viot process at
stationarity forwards in time (The Moran model falls within the attraction of
the Fleming-Viot processes). The Stephens-Donnelly algorithm is the continuous
version of such particle sampling algorithms:

  1. 1.

Start with a sample of size $m=1$. Randomly select an allele of type $k$ from
$1\leq k\leq K$ with probability $\Pr(k{\,|\,}\vec{\alpha})$ taken from the
stationary distribution characterized by $\mathbf{Q}$. Split immediately into
two lineages of this type.

  2. 2.

At this point there are $m$ lineages in the ancestry. Wait for an
exponentially distributed time with rate $m(m-1+\theta)$ and then select an
ancestral lineage at random. Split it with probability $m(m-1+\theta)$,
otherwise introduce a mutation.

  3. 3.

Hit a predetermined stopping criterion at sample size $M+1$, then go back to
the time when there were $M$ ancestral lineages.

##  4 Sampling Paths

###  4.1 Taylor Expansion

Given a sample of alleles drawn with the Stephens-Donnelly algorithm, one can
explicitly calculate the probabilities of the possible sampling paths that
lead to the observed sample.

A monomorphic sample of size $M$ is either the result of a sampling path that
consists purely of splitting ancestries or of one that has an even number of
reversible mutations (which we also take to include the unlikely case of four
mutations through all variants back to the original allele). In the former
case, the sampling path can be written as follows:

$\begin{split}f_{0}(\theta)&=\Pr(y_{k}=M\mid\vec{\alpha},M,\theta,\text{no
mutation})\\\
&=\Pr(k\mid\vec{\alpha})\cdot\frac{1}{1+\alpha_{k}\theta}\cdot\frac{2}{2+\alpha_{k}\theta}\cdots\frac{M-1}{(M-1)+\alpha_{k}\theta}\,.\end{split}$
(8)

With $\theta\ll 1$ and $\alpha_{i}$, $\alpha_{ij}$ of order one or smaller, a
Taylor expansion of the full sampling path at $\theta=0$ is to first order:

$\Pr(y_{k}=M\mid\vec{\alpha},M,\text{no
mutation})=\Pr(k\mid\vec{\alpha})-\Pr(k\mid\vec{\alpha})\cdot\sum_{m=2}^{M}\frac{\alpha_{k}\theta}{m-1}+\mathcal{O}(\theta^{2})\,.$
(9)

Here the full derivation: We have

$f_{0}(\theta=0)=\Pr(k\mid\vec{\alpha})$

and

$\begin{split}f_{0}^{\prime}(\theta)&=\Pr(k\mid\vec{\alpha})\cdot\bigg{(}-\frac{\alpha_{k}}{(1+\alpha_{k}\theta)^{2}}\bigg{[}\frac{2}{2+\alpha_{k}\theta}\cdot\frac{3}{3+\alpha_{k}\theta}\cdots\frac{(M-1)}{(M-1)+\alpha_{k}\theta}\bigg{]}-\\\
&\qquad-\frac{2\alpha_{k}}{(2+\alpha_{k}\theta)^{2}}\left[\frac{1}{1+\alpha_{k}\theta}\cdot\frac{3}{3+\alpha_{k}\theta}\cdots\frac{(M-1)}{(M-1)+\alpha_{k}\theta}\right]-\cdots\\\
&\qquad\cdots-\frac{(M-1)\alpha_{k}}{((M-1)+\alpha_{k}\theta)^{2}}\bigg{[}\frac{1}{1+\alpha_{k}\theta}\cdot\frac{2}{2+\alpha_{k}\theta}\cdots\frac{(M-2)}{(M-2)+\alpha_{k}\theta}\bigg{]}\bigg{)}\,,\end{split}$

such that

$f_{0}^{\prime}(\theta=0)=-\Pr(k\mid\vec{\alpha})\cdot\sum_{m=2}^{M}\frac{\alpha_{k}}{m-1}\,.$

Note that $m$ can be interpreted as the sample size at which a mutation could
potentially occur.

A sampling path with an even number of mutations that restores the monomorphic
condition is unlikely with low mutation rates. Indeed it is easy to show that,
in this case, probabilities of such sampling paths are at least of second
order in $\theta$ and thus do not contribute to the first order approximation.

Consider now the sampling paths that create polymorphic samples of size $M$:
The simplest of these includes one mutation, e.g., when the sample size $m$ is
equal to $y_{k}+1$:

$\begin{split}f_{1}(\theta)&=\Pr(y_{k}=M-y_{l}{\,|\,}\vec{\alpha},M,\theta,\text{one
mutation at $m=y_{k}+1$})\\\
&=\Pr(k{\,|\,}\vec{\alpha})\cdot\frac{1}{1+\alpha_{k}\theta}\cdot\frac{2}{2+\alpha_{k}\theta}\cdots\frac{y_{k}-1}{y_{k}-1+\alpha_{k}\theta}\cdot\frac{\alpha_{kl}\theta}{m-1+\alpha_{k}\theta}\\\
&\qquad\times\frac{1}{y_{k}+\frac{y_{k}}{y_{k}+1}\alpha_{k}\theta+\frac{1}{y_{k}+1}\alpha_{l}\theta}\cdots\frac{y_{l}-1}{M-1+\frac{y_{k}}{M}\alpha_{k}\theta+\frac{y_{l}}{M-y_{k}}\alpha_{l}\theta}\end{split}$
(10)

A Taylor expansion of the probability of this sampling path at $\theta=0$
yields:

$\Pr(y_{k}=M-y_{l}{\,|\,}\vec{\alpha},M,\text{one
mutation})=\Pr(k{\,|\,}\vec{\alpha})\cdot\frac{\alpha_{kl}}{m-1}\frac{(y_{k}-1)!(y_{m}-1)!}{(M-1)!}+\mathcal{O}(\theta^{2})$
(11)

The derivation is similar to those above: Set

$\begin{split}g_{y}(\theta)&=\Pr(k{\,|\,}\vec{\alpha})\cdot\frac{1}{1+\alpha_{k}\theta}\frac{2}{2+\alpha_{k}\theta}\cdots\frac{y_{k}-1}{y_{k}-1+\alpha_{k}\theta}\\\
&\qquad\times\frac{1}{y_{k}+\frac{y_{k}}{y_{k}+1}\alpha_{k}\theta+\frac{1}{y_{k}+1}\alpha_{l}\theta}\cdots\frac{y_{l}-1}{M-1+\frac{y_{k}}{M}\alpha_{k}\theta+\frac{y_{l}}{M-y_{k}}\alpha_{l}\theta}\,,\end{split}$

such that

$\begin{split}g_{y}(\theta=0)&=\Pr(k{\,|\,}\vec{\alpha})\frac{(y_{k}-1)!(y_{l}-1)!}{(M-1)!}\,.\end{split}$
(12)

Then we have

$\begin{split}f_{1}(\theta)&=\frac{\alpha_{kl}\theta}{(m-1)+\alpha_{k}\theta}g_{y}(\theta)\\\
f_{1}^{\prime}(\theta)&=\bigg{(}\frac{\alpha_{kl}\theta}{(m-1)+\alpha_{k}\theta}\bigg{)}^{\prime}g_{m}(\theta)+\frac{\alpha_{kl}\theta}{(m-1)+\alpha_{k}\theta}g_{m}^{\prime}(\theta)\,.\end{split}$
(13)

Setting $\theta=0$, we get

$\begin{split}f_{1}^{\prime}(\theta=0)&=\frac{\alpha_{kl}}{m-1}g_{m}(\theta=0)+0\\\
&=\Pr(k{\,|\,}\vec{\alpha})\cdot\frac{\alpha_{kl}}{m-1}\frac{(y_{k}-1)!(y_{m}-1)!}{(M-1)!}\,.\end{split}$
(14)

The occurrence of two or three mutations that lead to a polymorphic sample
with two, three, or four segregating alleles is a theoretical possibility.
However, the expansion of the sampling path probabilities would again
approximate to zero for the same reason as in the case of (multiple)
reversible mutations in a monomorphic sampling path.

Note that the mutation may in principle occur at any sample size between
$2\leq m_{k}\leq(y_{k}+1)$. Hence, while the series expansion of the
monomorphic sampling paths accounts for the only possible sampling scenario,
the expansion of a polymorphic sampling path represents one of several
feasible branching structures along the sampling algorithm that result in the
given configuration of alleles. We thus need to sum over these possibilities.

###  4.2 Sum of Ordered Probabilities

In the following subsection, we show that the sampling distribution of a
polymorphic sampling path, conditional on the mutation occurring at sample
size $m$, is beta-binomial. The argument is similar to that in Burden and
Griffiths [1]. These authors use a direct result equivalent to the deFinetti
density of a Polya urn model. We start from the properties of the sampling
path and build up the deFinetti representation.

According to deFinetti’s theorem [3], there exists for every infinite sequence
of Bernoulli random variables a probability distribution $F$ on
$\left[0,1\right]$ so that

$\begin{split}&\Pr(X_{1}=1,\cdots,X_{k}=1,X_{k+1}=0,\cdots,X_{n}=0)=\\\
&\qquad\int_{0}^{1}\nu^{k}(1-\nu^{n-k})dF(\nu)=E\left[U^{k}(1-U^{n-k})\right]\end{split}$
(15)

for a random variable $U$, whereby

$\Pr(U\leq\nu)=F(\nu)$

and the $(X_{i})_{i}\leq 1$ are i.i.d. when conditioned on a $U$ that fulfills
$\Pr(X_{i}=k)=U$. The following also holds:

$\Pr(S_{n}=\sum_{i}^{n}X_{i}=k)=\int_{0}^{1}\binom{n}{k}\nu^{k}(1-\nu^{n-k})dF(\nu)\,.$

The joint distribution of the polymorphic sampling path, which is a finite
number of draws of alleles, does not change if the order of alleles within the
path is changed. We start from an initial allele of type $k$, then there is a
mutational event to allele $l$ at a random time point. At this point, $m-1$
alleles of type $k$ are already in the sample. Each subsequent draw can
increase the number of sampled alleles of either type by one. This is exactly
the specification of the following Polya urn process:

$(1,m-1)\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}\,,$

which in our case runs until we have $y_{l}-1$ alleles of type $l$ and
$y_{k}-m=M-y_{l}-m$ alleles of type $k$. By noting that the probabilities of
specific sampling paths are sequences of beta distributed moments and
recalling Hausdorff’s moment problem [7, 8], it follows that

$U\sim B(1,m-1)\,.$

Therefore, the probability of a process that yields $y$ alleles of type $l$ in
a sample size of $M-m$ (denoted by $S_{M-m}$) is given by the following:

$\begin{split}&\Pr(S_{M-m}{\,|\,}M,m,1)\\\
&=\int_{0}^{1}\binom{M-m}{y_{l}-1}\nu^{y_{l}-1}(1-\nu)^{M-y_{l}-m}(1-\nu)^{m-1}\nu^{0}\frac{\Gamma(m)}{\Gamma(m-1)\Gamma(1)}d\nu\\\
&=\binom{M-m}{y_{l}-1}\frac{\Gamma(M-y_{l}-m+m-1)\Gamma(y_{l})}{\Gamma(M-y_{l}-m+m-1+y_{l})}\frac{\Gamma(m)}{\Gamma(m-1)\Gamma(1)}\\\
&=\binom{M-m}{y_{l}-1}\frac{\Gamma(M-y_{l})\Gamma(y_{l})}{\Gamma(M)}\frac{\Gamma(m)}{\Gamma(m-1)\Gamma(1)}\\\
&=\binom{M-m}{y_{l}-1}\frac{\Gamma(M-y_{l})\Gamma(y_{l})}{\Gamma(M)}\frac{(m-1)!}{(m-2)!(1-1)!}\\\
&=\binom{M-m}{y_{l}-1}(m-1)\frac{\Gamma(M-y_{l})\Gamma(y_{l})}{\Gamma(M)}\,.\end{split}$
(16)

For a given proportion $y_{l}$ of mutant alleles, we need to sum over all
possible mutation points, i.e., $2\leq m\leq(M-y_{l}+1)$. Given that the first
allele sampled at sample size $m=1$ is of type $k$, we have:

$\begin{split}\Pr(y_{l}{\,|\,}M,\vec{\alpha},k,\theta\to
0)&=\sum_{m=2}^{M-y_{l}+1}\binom{M-m}{y_{l}-1}(m-1)\frac{\Gamma(M-y_{l})\Gamma(y_{l})}{\Gamma(M)}\frac{\alpha_{kl}\theta}{m-1}+\mathcal{O}(\theta^{2})\\\
&=\alpha_{kl}\theta\sum_{m=2}^{M-y_{l}+1}\binom{M-m}{y_{l}-1}\frac{\Gamma(M-y_{l})\Gamma(y_{l})}{\Gamma(M)}+\mathcal{O}(\theta^{2})\\\
&=\frac{\alpha_{kl}\theta}{y_{l}}+\mathcal{O}(\theta^{2})\,.\end{split}$ (17)

Line three follows from line two by repeatedly applying the identity:

$\begin{split}\binom{N}{x}&=\binom{N-1}{x-1}+\binom{N-1}{x}\Leftrightarrow\binom{N-1}{x-1}=\binom{N}{x}-\binom{N-1}{x}\,,\end{split}$
(18)

from which it follows that

$\begin{split}&\sum_{m=2}^{M-y_{l}+1}\binom{M-m}{y_{l}-1}\\\
&=\binom{M-2}{y_{l}-1}+\binom{M-3}{y_{l}-1}+\dots+\binom{y_{1}}{y_{l}-1}+\binom{y_{l}-1}{y_{l}-1}\\\
&=\binom{M-1}{y_{l}}-\binom{M-2}{y_{l}}+\binom{M-2}{y_{l}}-\binom{M-3}{y_{l}}+\\\
&\qquad\dots+\binom{y_{l}+1}{y_{l}}-\binom{y_{l}}{y_{l}}+\binom{y_{l}}{y_{l}}-\binom{y_{l}-1}{y_{l}}\\\
&=\binom{M-1}{y_{l}}\,,\end{split}$ (19)

since $\binom{y_{l}-1}{y_{l}}=0$. Note that equation (19) also appears in
Burden and Griffiths [1], but is solved differently.

This provides an alternative proof of Theorem 1 in Burden and Griffiths [1]
and concludes the characterization of the polymorphic sampling paths.

##  5 Stationary Distribution

On slightly reparametrizing the results of the previous two subsections, the
probability of obtaining the overall configuration $y$ of alleles can be
concisely formulated: Given a fixed sample size
$M<e^{\frac{1}{\alpha_{k}\theta}}$, mutation biases of maximal order $1$, and
a small scaled overall mutation rate $\theta=N\mu<0.1$, the stationary allelic
configuration $Y$ to the first order of $\theta$ is the distribution
$\Pr(y\mid M,\vec{\alpha},\theta)$:

$\begin{cases}\Pr(k\mid\vec{\alpha})-\Pr(k\mid\vec{\alpha})\cdot\sum_{k\neq
l}^{K}\alpha_{k}\theta H_{M-1}+\mathcal{O}(\theta^{2})&y_{k}=M,y_{l\neq
k}=0\\\
\Pr(k\mid\vec{\alpha})\frac{\alpha_{kl}\theta}{y_{l}}+\Pr(l\mid\vec{\alpha})\frac{\alpha_{lk}\theta}{y_{k}}+\mathcal{O}(\theta^{2})&1\leq
y_{k},y_{l\neq k}\leq M-1\\\ &\text{and }y_{l}+y_{k}=M\,.\end{cases}$ (20)

All other possibilities have probabilities of at most
$\mathcal{O}(\theta^{2})$ and thus do not contribute to the distribution. This
stationary distribution is equivalent to that derived by Burden and Griffiths
[1], as expected.

##  6 Maximum likelihood (ML) estimators of the scaled mutation rate
parameters

Starting from the general stationary distribution and recalling that we are
certain of a unique stationary distribution of $\mathbf{Q}$ in the case of
strand symmetry (i.e., we can build a sample of size $M$ using the Stephens-
Donnelly algorithm if we assume strand symmetry), we now aim to infer the six
parameters of the strand-symmetric boundary-mutation Moran model given site
frequency spectrum data.

Figure 1: Transition Rates for the Strand Symmetric Model

Recalling the parametrization for $\mathbf{Q}$ in formula (5) as visualized in
Fig 1, we set $\theta=-\lambda_{2}=b+d+c+e$. Then
$\alpha_{1}=\frac{a+c+e}{b+d+c+e}$, $\alpha_{12}=\frac{a}{b+d+c+e}$,
$\alpha_{13}=\frac{c}{b+d+c+e}$, $\alpha_{14}=\frac{e}{b+d+c+e}$, etc.

The full stationary distribution for the strand-symmetric model is then:

$\begin{split}&\Pr(y\mid M,\vec{\alpha},\theta)=\\\
&\begin{cases}\frac{b+d}{2(b+d+c+e)}-\frac{b+d}{2(b+d+c+e)}(a+c+e)H_{M-1}+\mathcal{O}(\theta^{2})&y_{A}=M,y_{T,G,C}=0\\\
\frac{b+d}{2(b+d+c+e)}-\frac{b+d}{2(b+d+c+e)}(a+c+e)H_{M-1}+\mathcal{O}(\theta^{2})&y_{T}=M,y_{A,G,C}=0\\\
\frac{c+e}{2(b+d+c+e)}-\frac{c+e}{2(b+d+c+e)}(b+d+f)H_{M-1}+\mathcal{O}(\theta^{2})&y_{G}=M,y_{A,T,C}=0\\\
\frac{c+e}{2(b+d+c+e)}-\frac{c+e}{2(b+d+c+e)}(b+d+f)H_{M-1}+\mathcal{O}(\theta^{2})&y_{C}=M,y_{A,T,G}=0\\\
\frac{b+d}{2(b+d+c+e)}\frac{a}{y_{A}}+\frac{b+d}{2(b+d+c+e)}\frac{a}{y_{T}}+\mathcal{O}(\theta^{2})&1\leq
y_{A},y_{T}\leq M-1,y_{A}+y_{T}=M\\\
\frac{b+d}{2(b+d+c+e)}\frac{e}{y_{A}}+\frac{c+e}{2(b+d+c+e)}\frac{d}{y_{C}}+\mathcal{O}(\theta^{2})&1\leq
y_{A},y_{C}\leq M-1,y_{A}+y_{C}=M\\\
\frac{b+d}{2(b+d+c+e)}\frac{c}{y_{A}}+\frac{c+e}{2(b+d+c+e)}\frac{b}{y_{G}}+\mathcal{O}(\theta^{2})&1\leq
y_{A},y_{G}\leq M-1,y_{A}+y_{G}=M\\\
\frac{b+d}{2(b+d+c+e)}\frac{c}{y_{T}}+\frac{c+e}{2(b+d+c+e)}\frac{b}{y_{C}}+\mathcal{O}(\theta^{2})&1\leq
y_{T},y_{C}\leq M-1,y_{T}+y_{C}=M\\\
\frac{b+d}{2(b+d+c+e)}\frac{e}{y_{T}}+\frac{c+e}{2(b+d+c+e)}\frac{d}{y_{G}}+\mathcal{O}(\theta^{2})&1\leq
y_{T},y_{G}\leq M-1,y_{T}+y_{G}=M\\\
\frac{c+e}{2(b+d+c+e)}\frac{f}{y_{G}}+\frac{c+e}{2(b+d+c+e)}\frac{f}{y_{C}}+\mathcal{O}(\theta^{2})&1\leq
y_{G},y_{C}\leq M-1,y_{G}+y_{C}=M\,.\end{cases}\end{split}$ (21)

In the next subsections, we will derive ML estimators for the parameter vector
$(a,b,c,d,e,f)$ in the distribution (21) by linear transformation.

###  6.1 Biallelic Model

Grouping the bases $A$ and $T$ as well as $G$ and $C$ together results in a
biallelic model [16]: If one interprets the allelic state $k$ as representing
the grouped state $A-T$ including those sites polymorphic for $A$/$T$ and the
allelic state $l$ as representing the $G-C$ alleles and the polymorphic $G/C$
sites, $\alpha_{kl}$ and $\alpha_{lk}$ become equal to $\alpha_{k}$ and
$\alpha_{l}$ respectively (see Figure 2).

Figure 2: Transition Rates for the Biallelic Strand Symmetric Model

In particular, this biallelic model enables us to infer
$\Pr(k=AT\mid\vec{\alpha})=\frac{b+d}{b+d+c+e}$ and
$\Pr(l=GC\mid\vec{\alpha})=\frac{e+c}{b+d+c+e}$. Using the frequencies of $A$
and $T$, one can then infer $a$ using a symmetric version of the biallelic
model and similarly infer $f$ from the frequencies of $G$ and $C$. The
remaining transition rates can be estimated via the expectation-maximization
algorithm.

Following Vogl [16], we will first find the ML estimators for
$\Pr(k=AT\mid\vec{\alpha})=\frac{b+d}{b+d+c+e}=\beta$ and
$\Pr(l=GC\mid\vec{\alpha})=\frac{e+c}{b+d+c+e}=1-\beta$ using
$\alpha_{13}+\alpha_{14}=1-\beta$ and starting from the following biallelic
version of 20:

$\Pr(y\mid M,\vec{\alpha},\theta)=\begin{cases}\beta-\beta(1-\beta)\theta
H_{M-1}+\mathcal{O}(\theta^{2})&y_{k}=M,y_{l\neq k}=0\\\
\beta(1-\beta)\theta\bigg{(}\frac{1}{y_{l}}+\frac{1}{y_{k}}\bigg{)}+\mathcal{O}(\theta^{2})&1\leq
y_{k},y_{l\neq k}\leq M-1,\\\ &y_{l}+y_{k}=M\\\ (1-\beta)-\beta(1-\beta)\theta
H_{M-1}+\mathcal{O}(\theta^{2})&y_{k}=0,y_{l\neq k}=M\,.\end{cases}$ (22)

The estimators are to be expressed via site frequency data: We assume
knowledge of the allele frequency counts at $L<\infty$ loci indexed as $1\leq
l\leq L$ and further assume that these counts are based on a sample of size
$M$ at each locus (the constant sample size merely simplifies the notation and
calculations and is not a theoretical necessity). $L_{0}$ is the number of
loci at which there are no occurrences of the state $A-T$ (i.e., the number of
sites monomorphic for $G$ or $C$, or $G/C$ polymorphic), and $L_{M}$ is the
number of loci either monomorphic for $A$ or $T$, or $A/T$ polymorphic. The
remaining $L_{y}$ with $L_{p}=\sum_{y=1}^{M-1}L_{y}$ are the observed
polymorphic counts between the $A-T$ and $G-C$ states. We define

$\gamma=\beta(1-\beta)\theta=\frac{(b+d)(e+c)}{b+d+c+e}$

to efficiently handle products of random variables.

The likelihood of the site frequency spectrum then becomes:

$\begin{split}&\Pr(L_{0}\dots L_{M}\mid M,\vec{\alpha},\gamma)=\\\
&\qquad\frac{L!}{\prod_{y=0}^{M}L_{y}!}(\beta-\gamma
H_{M-1})^{L_{M}}(2\beta\gamma H_{M-1})^{L_{p}}((1-\beta)-\gamma
H_{M-1})^{L_{0}}\,,\end{split}$ (23)

whereby $H_{m}$ again denotes the $m$th harmonic number.

Then the log-likelihood can be written as proportional to:

$\begin{split}\log(\Pr(L_{0}\dots L_{M}\mid
M,\vec{\alpha},\gamma))&=const+{L_{M}}\log(\beta-\gamma H_{M-1})\\\
&\qquad+{L_{p}}\log(2\beta\gamma H_{M-1})+{L_{0}}\log((1-\beta)-\gamma
H_{M-1})\,.\end{split}$ (24)

With asymmetric mutation rates between the two classes, we can take the
derivative of the log-likelihood by $\beta$:

$\begin{split}\frac{-L_{0}}{(1-\beta)-\gamma
H_{M-1}}+\frac{L_{M}}{\beta-\gamma H_{M-1}}&\overset{!}{=}0\\\
-L_{0}\beta+L_{0}\gamma H_{M-1}+L_{M}-L_{M}\beta-L_{M}\gamma
H_{M-1}&=0\,.\end{split}$ (25)

It follows that

$\beta=\frac{L_{0}\gamma H_{M-1}+L_{M}(1-\gamma H_{M-1})}{L_{0}+L_{M}}\,,$
(26)

and therefore

$\begin{split}\beta-\gamma H_{M-1}&=\frac{-L_{0}\gamma H_{M-1}-L_{M}\gamma
H_{M-1}+L_{0}\gamma H_{M-1}+L_{M}(1-\gamma H_{M-1})}{L_{0}+L_{M}}\\\
&=\frac{L_{M}(1-2\gamma H_{M-1})}{L_{0}+L_{M}}\,,\end{split}$ (27)

and similarly

$(1-\beta)-\gamma H_{M-1}=\frac{L_{0}(1-2\gamma H_{M-1})}{L_{0}+L_{M}}\,.$
(28)

Taking the derivative of the log-likelihood by $2\gamma H_{M-1}$ and
substituting we get:

$\begin{split}&\log(\Pr(L_{0}\dots L_{M}\mid
M,\vec{\alpha},\gamma)=\frac{-L_{M}(L_{0}+L_{M})\frac{L_{M}}{L_{0}+L_{M}}}{L_{M}(1-2\gamma
H_{M-1})}+\frac{L_{p}}{2\gamma H_{M-1}}\\\
&\qquad-\frac{L_{0}(L_{0}+L_{M})\frac{L_{0}}{L_{0}+L_{M}}}{L_{0}(1-2\gamma
H_{M-1})}\overset{!}{=}0\\\ &-(L_{0}+L_{M})2\gamma H_{M-1}+L_{p}(1-2\gamma
H_{M-1})=0\,.\end{split}$ (29)

It follows that

$\widehat{\gamma}=\frac{L_{p}}{2(L_{0}+L_{M}+L_{p})H_{M-1}}\,.$ (30)

This corresponds to the estimator in formula (36) in [16]. Analogously, we
get:

$\beta-\widehat{\gamma}H_{M-1}=\frac{L_{M}}{(L_{0}+L_{M}+L_{p})}\,.$ (31)

Similarly, we have

$(1-\beta)-\widehat{\gamma}H_{M-1}=\frac{L_{0}}{(L_{0}+L_{M}+L_{p})}$ (32)

respectively yielding

$\widehat{\Pr(k=AT\mid\vec{\alpha})}=\widehat{\beta}=\frac{L_{M}+\frac{L_{p}}{2}}{(L_{0}+L_{M}+L_{p})}\,,$
(33)

which corresponds to the estimator in formula (37) in [16], and

$\widehat{\Pr(k=GC\mid\vec{\alpha})}=1-\widehat{\beta}=\frac{L_{0}+\frac{L_{p}}{2}}{(L_{0}+L_{M}+L_{p})}\,.$
(34)

Thus we provide ML estimators of the parameters of the combined states $A-T$
and $G-C$ in the biallelic model resulting from combining $A-T$ and $C-G$
under stationarity.

###  6.2 Symmetric Biallelic Model

Figure 3: Transition Rates for the Symmetric Biallelic Strand Symmetric Model

Next we recover estimators for the transition rates $a$ and $f$ from symmetric
biallelic systems [see 2]. Let us focus on $a$ (see Figure 3): Recalling
Chargaff’s second parity rule and its correspondence with the stationary
distribution of base frequencies, we can write

$\Pr(k=A\mid\vec{\alpha})=\Pr(k=T\mid\vec{\alpha})=\frac{\Pr(k=AT\mid\vec{\alpha})}{2}=\frac{\beta}{2}=\frac{b+d}{2(b+d+c+e)}.$

Furthermore, we know that

$\gamma=\beta(1-\beta)\theta=\frac{b+d}{b+d+c+e}(c+e)\,,$

from which it follows that

$\frac{\beta}{2}(a+c+e)=\frac{\beta}{2}a+\frac{\gamma}{2}$

Subsequently the distribution (20) can be rewritten for a symmetric biallelic
A-T model using the previously estimated parameters:

$\Pr(y\mid
M,a,\beta,\gamma)=\begin{cases}\frac{\beta}{2}-\frac{\beta}{2}aH_{M-1}-\frac{\gamma}{2}H_{M-1}+\mathcal{O}(\theta^{2})&y_{A}=M,y_{T}=0\\\
&\text{ or }y_{A}=0,y_{T}=M\\\ \beta
a\bigg{(}\frac{1}{y_{A}}+\frac{1}{y_{T}}\bigg{)}+\mathcal{O}(\theta^{2})&1\leq
y_{A},y_{T}\leq M-1,\\\ &y_{A}+y_{T}=M\end{cases}$ (35)

This yields the following log-likelihood of the site frequency spectrum:

$\begin{split}&\log(\Pr(L_{0},\dots,L_{M}\mid
M,a,\beta,\gamma))=const+{L_{M}}\log\bigg{(}\frac{\beta}{2}-\frac{\beta}{2}aH_{M-1}-\frac{\gamma}{2}H_{M-1}\bigg{)}\\\
&\qquad+{L_{p}}\log(\beta
aH_{M-1})+{L_{0}}\log\bigg{(}\frac{\beta}{2}-\frac{\beta}{2}aH_{M-1}-\frac{\gamma}{2}H_{M-1}\bigg{)}\\\
\end{split}$ (36)

Taking the derivative by $a$, one can calculate the ML estimator for $a$ as
follows:

$\begin{split}\frac{-L_{M}\frac{\beta}{2}H_{M-1}}{\frac{\beta}{2}-\frac{\beta}{2}aH_{M-1}-\frac{\gamma}{2}H_{M-1}}+\frac{L_{p}\beta
H_{M-1}}{\beta
aH_{M-1}}-\frac{L_{0}\frac{\beta}{2}H_{M-1}}{\frac{\beta}{2}-\frac{\beta}{2}aH_{M-1}-\frac{\gamma}{2}H_{M-1}}&\overset{!}{=}0\\\
\frac{-L_{M}}{\frac{1}{H_{M-1}}-a-(1-\beta)\theta}+\frac{L_{p}}{a}-\frac{L_{0}}{\frac{1}{H_{M-1}}-a-(1-\beta)\theta}&=0\\\
\frac{L_{p}\bigg{(}\frac{1}{H_{M-1}}-(1-\beta)\theta\bigg{)}}{(L_{M}+L_{p}+L_{0})}&=a\\\
\end{split}$ (37)

From earlier we have

$(1-\beta)\theta=\frac{\gamma}{\beta}\,,$

such that

$\widehat{(1-\beta)\theta}=\frac{\hat{\gamma}}{\hat{\beta}}=\frac{L_{p}}{2(L_{M}+\frac{L_{p}}{2})H_{M-1}}\,,$

and furthermore

$\frac{1}{H_{M-1}}-\widehat{(1-\beta)\theta}=\frac{2L_{M}}{2H_{M-1}(L_{M}+\frac{L_{p}}{2})}\,.$

As a result the ML estimator for $a$ becomes:

$\begin{split}\hat{a}&=\frac{L_{p}L_{M}}{(L_{M}+L_{p}+L_{0})(L_{M}+\frac{L_{p}}{2})H_{M-1}}\\\
\end{split}$ (38)

Analogously, we get

$\begin{split}\hat{f}&=\frac{L_{p}L_{0}}{(L_{M}+L_{p}+L_{0})(L_{0}+\frac{L_{p}}{2})H_{M-1}}\,.\end{split}$
(39)

###  6.3 Expectation-Maximization Algorithm

The remaining transition rates $(b,d,e,c)$ cannot be disentangled and
expressed individually through reparametrization in the same way as $a$ and
$f$. Instead, we will use the expectation-maximization (EM) algorithm [4] to
obtain estimators of these transition rates as sequential updates that cycle
through the parameter pairs $d,e$ and $b,c$, respectively.

Figure 4: Transition Rates for the Diagonally Symmetric Biallelic Strand
Symmetric Model

Let us focus on a biallelic A-C system (see Figure 4) and rewrite 20 to obtain
the appropriate stationary distribution:

$\Pr(y\mid
M,e,d,\beta,\theta)=\begin{cases}\frac{\beta}{2}-\frac{\beta}{2}(a+c+e)H_{M-1}+\mathcal{O}(\theta^{2})&y_{A}=M,y_{C}=0\\\
\frac{\beta}{2}\frac{e}{y_{A}}+\bigg{(}\frac{1-\beta}{2}\bigg{)}\frac{d}{y_{C}}+\mathcal{O}(\theta^{2})&1\leq
y_{A},y_{C}\leq M-1,\\\ &y_{A}+y_{C}=M\\\
\bigg{(}\frac{1-\beta}{2}\bigg{)}-\bigg{(}\frac{1-\beta}{2}\bigg{)}(d+b+f)H_{M-1}+\mathcal{O}(\theta^{2})&y_{C}=M,y_{A}=0\,.\end{cases}$
(40)

Let $L_{i}$ denote the number of polymorphic samples in the frequency
spectrum, whereby $y_{A}=i$ and $y_{C}=M-i$. In order to distinguish between
mutations of different directions, we introduce the auxiliary variable $z_{i}$
($0\leq z_{i}\leq L_{i}$) that counts the mutations from $A$ to $C$.
Conversely, $(L_{i}-z_{i})$ counts the mutations from $C$ to $A$. This yields
the complete data log-likelihood:

$\begin{split}&\Pr(y,z_{i}{\,|\,}L_{0},\dots,L_{i},\dots,L_{M},\dots)=L_{M}\log\bigg{(}\frac{\beta}{2}-\frac{\beta}{2}(a+c+e)H_{M-1}\bigg{)}\\\
&\qquad+\sum_{i=1}^{M-1}z_{i}\log\bigg{(}\frac{\beta}{2}\frac{e}{i}\bigg{)}+\sum_{i=1}^{M-1}(L_{i}-z_{i})\log\bigg{(}\frac{1-\beta}{2}\frac{d}{M-i}\bigg{)}\\\
&\qquad+L_{0}\log\bigg{(}\frac{1-\beta}{2}-\frac{1-\beta}{2}(d+b+f)H_{M-1}\bigg{)}\,.\end{split}$
(41)

The expectation of $z_{i}$ corresponds to the mean of a binomial distribution
with sample size $L_{i}$ and
$p_{i}=\frac{\frac{e}{i}}{\frac{e}{i}+\frac{d}{M-i}}$:

$\operatorname{E}(z_{i}{\,|\,}L_{i},e^{t},d^{t})=L_{i}\frac{\frac{e^{t}}{i}}{\frac{e^{t}}{i}+\frac{d^{t}}{M-i}}\,.$
(42)

The expectation step of the EM-algorithm constitutes taking the expectation of
the full log-likelihood and noting that only the part
$\operatorname{Q}(d,e{\,|\,}d^{t},e^{t})=\operatorname{E}(\Pr(y,z_{i}{\,|\,}L_{0},\dots,L_{i},\dots,L_{M},\dots))$
of the function needs to be maximized:

$\begin{split}&\operatorname{Q}(d,e{\,|\,}d^{t},e^{t})=const+L_{M}\log\bigg{(}1-(a+c+e)H_{M-1}\bigg{)}\\\
&\qquad+\sum_{i=1}^{M-1}L_{i}\frac{\frac{e^{t}}{i}}{\frac{e^{t}}{i}\frac{d^{t}}{M-i}}\log(e)+\sum_{i=1}^{M-1}L_{M-i}\frac{\frac{d^{t}}{M-i}}{\frac{e^{t}}{i}+\frac{d^{t}}{M-i}}\log(d)\\\
&\qquad+L_{0}\log\bigg{(}1-(d+b+f)H_{M-1}\bigg{)}\,.\end{split}$ (43)

To calculate the parameter updates according to the maximization step, we must
take the appropriate derivatives of $\operatorname{Q}(d,e{\,|\,}d^{t},e^{t})$.
For the calculation of $e^{t+1}$, we take the derivative by $e$ and set it to
zero:

$\begin{split}\frac{d}{de}\operatorname{Q}(d,e{\,|\,}d^{t},e^{t})&\overset{!}{=}0\\\
0&=\frac{-L_{M}H_{M-1}}{1-(a+c+e)H_{M-1}}+\sum_{i=1}^{M-1}\frac{\frac{L_{i}\frac{e^{t}}{i}}{\frac{e^{t}}{i}+\frac{d^{t}}{M-i}}}{e}\\\
0&=\frac{-L_{M}H_{M-1}}{1-(a+(1-\beta)\theta)H_{M-1}}+\frac{\sum_{i=1}^{M-1}L_{i}\frac{\frac{e^{t}}{i}}{\frac{e^{t}}{i}+\frac{d^{t}}{M-i}}}{e}\\\
e&=\frac{1-(a+(1-\beta)\theta)H_{M-1}}{L_{M}H_{M-1}}\sum_{i=1}^{M-1}L_{i}\frac{\frac{e^{t}}{i}}{\frac{e^{t}}{i}+\frac{d^{t}}{M-i}}\,,\end{split}$
(44)

where we substituted $(1-\beta)\theta$ for $c+e$ to obtain the third line.

Substituting the ML estimators for $a$ and $(1-\beta)\theta$ in the numerator
of the first factor yields:

$\begin{split}(\hat{a}+\widehat{(1-\beta)\theta})H_{M-1}&=\frac{L_{p}L_{M}}{(L_{M}+L_{p}+L_{0})(L_{M}+\frac{L_{p}}{2})}+\frac{L_{p}}{2(L_{M}+\frac{L_{p}}{2})}\\\
&=\frac{3L_{p}L_{M}+L_{p}L_{p}+L_{p}L_{0}}{2(L_{M}+L_{p}+L_{0})(L_{M}+\frac{L_{p}}{2})}\,.\end{split}$
(45)

Furthermore, we have:

$\begin{split}1-(\hat{a}+\widehat{(1-\beta)\theta})H_{M-1}&=\frac{2L_{0}L_{M}+2L_{M}L_{M}}{2(L_{M}+L_{p}+L_{0})(L_{M}+\frac{L_{p}}{2})}\\\
&=\frac{L_{M}(L_{0}+L_{M})}{(L_{M}+L_{p}+L_{0})(L_{M}+\frac{L_{p}}{2})}\,.\end{split}$
(46)

Therefore, we get:

$\begin{split}\hat{e}=\frac{(L_{0}+L_{M})}{(L_{M}+L_{p}+L_{0})(L_{M}+\frac{L_{p}}{2})H_{M-1}}\sum_{i=1}^{M-1}L_{i}\frac{\frac{e^{t}}{i}}{\frac{e^{t}}{i}+\frac{d^{t}}{M-i}}\,.\end{split}$
(47)

Similarly, one obtains $d^{t+1}$ by taking the derivative of
$Q(d,e{\,|\,}d^{t},e^{t})$ by $d$ and setting it to zero.

The overall iteration scheme is then given by:

$\begin{split}e^{t+1}&=\frac{(L_{0}+L_{M})}{(L_{M}+L_{p}+L_{0})(L_{M}+\frac{L_{p}}{2})H_{M-1}}\sum_{i=1}^{M-1}L_{i}\frac{\frac{e^{t}}{i}}{\frac{e^{t}}{i}+\frac{d^{t}}{M-i}}\\\
d^{t+1}&=\frac{(L_{0}+L_{M})}{(L_{M}+L_{p}+L_{0})(L_{0}+\frac{L_{p}}{2})H_{M-1}}\sum_{i=1}^{M-1}L_{i}\frac{\frac{d^{t}}{i}}{\frac{e^{t}}{i}+\frac{d^{t}}{M-i}}\,.\end{split}$
(48)

Considering a biallelic A-G system and using the EM-algorithm analogously, the
parameter updates for $c$ and $b$ can also be determined:

$\begin{split}c^{t+1}&=\frac{(L_{0}+L_{M})}{(L_{M}+L_{p}+L_{0})(L_{M}+\frac{L_{p}}{2})H_{M-1}}\sum_{i=1}^{M-1}L_{i}\frac{\frac{c^{t}}{i}}{\frac{c^{t}}{i}+\frac{b^{t}}{M-i}}\\\
b^{t+1}&=\frac{(L_{0}+L_{M})}{(L_{M}+L_{p}+L_{0})(L_{0}+\frac{L_{p}}{2})H_{M-1}}\sum_{i=1}^{M-1}L_{i}\frac{\frac{b^{t}}{i}}{\frac{c^{t}}{i}+\frac{b^{t}}{M-i}}\,.\end{split}$
(49)

For both pairs of parameters, cyclical calculation of estimators guarantees
convergence towards a local maximum of the marginal likelihood by properties
of the EM-algorithm. Furthermore:
$\Pr(y{\,|\,}L_{0},\dots,L_{i},\dots,L_{M},\dots)$ is the distribution of
counts of each type of segregating allele. The configurations of the types of
alleles themselves are constructed via Polya-like urn processes. As such, the
marginal likelihood $\Pr(y{\,|\,}L_{0},\dots,L_{i},\dots,L_{M},\dots)$ takes
the form of a Dirichlet-multinomial distribution, which is known to be
unimodal [10]. Therefore the estimators determined via the expectation-
maximization algorithm converge towards the global optimum of the marginal
distribution.

##  7 Summary

An alternative derivation of the distribution of a sample of size $M$ is
derived given a general mutation model under the assumptions of small scaled
mutation rates and mutational equilibrium. Further assuming a standard four
allele DNA model with strand-symmetric mutation on complementary DNA strands
and available site frequency spectrum data, ML estimators for the six scaled
mutation parameters are determined. This is the first time such estimators are
provided for a mutation model more complex than the parent-independent
mutation model.

## Acknowledgments

This article is a condensed version of LCM’s master’s thesis supervised by CV.
CV wants to thank Juraj Bergman and Sandra Peer for discussions. CV’s research
is supported by the Austrian Science Fund (FWF): DK W1225-B20; LCM’s by the
School of Biology at the University of St.Andrews.

## References

## References

  * Burden and Griffiths [2018] Burden, C. and Griffiths, R. (2018).  The stationary distribution of a sample from the Wright-Fisher diffusion model with general small mutation rates.  Journal of Mathematical Biology, 78, 1211–1224.
  * Burden and Tang [2016] Burden, C. and Tang, Y. (2016).  An approximate stationary solution for multi-allele neutral diffusion with low mutation rates.  Theoretical Population Biology, 112, 22–32.
  * de Finetti [1931] de Finetti, B. (1931).  Funzione caratteristica di un fenomeno aleatorio.  Atti della R. Academia Nazionale dei Lincei, Serie 6. Memorie, Classe di Scienze Fisiche, Mathematice e Naturale, 4, 251–299.
  * Dempster et al. [1977] Dempster, A. P., Laird, N. M., and Rubin, D. B. (1977).  Maximum Likelihood from Incomplete Data via the EM Algorithm.  Journal of the Royal Statistical Society. Series B, 39(1), 1–38.
  * Donnelly and Kurtz [1996] Donnelly, P. and Kurtz, T. (1996).  A Countable Representation of the Fleming-Viot Measure-Valued Diffusion.  The Annals of Probability, 24(2), 698–742.
  * Etheridge [2012] Etheridge, A. (2012).  Some Mathematical Models from Population Genetics: Lecture Notes in Mathematics.  Springer Verlag. Berlin, Heidelberg.
  * Hausdorff [1921a] Hausdorff, F. (1921a).  Summationsmethoden und Momentenfolgen 1.  Mathematische Zeitschrift, 9, 74–109.
  * Hausdorff [1921b] Hausdorff, F. (1921b).  Summationsmethoden und Momentenfolgen 2.  Mathematische Zeitschrift, 9, 280–299.
  * Hoppe [1987] Hoppe, F. M. (1987).  The Sampling Theory of Neutral Alleles and an Urn Model in Population Genetics.  Journal of Mathematical Biology, 25, 123–159.
  * Levin and Reeds [1977] Levin, B. and Reeds, J. (1977).  Compound Multinomial Likelihood Functions are Unimodal: Proof of a Conjecture of I.J. Good.  The Annals of Statistics, 5(1), 79–87.
  * Lobry [1995] Lobry, J. (1995).  Properties of a General Model of DNA Evolution Under No-Strand-Bias Conditions.  Journal of Molecular Evolution, 40, 326–330.
  * Lynch et al. [2016] Lynch, M., Ackerman, M., Gout, J., Long, H., Sung, W., Thomas, W., and Foster, P. (2016).  Genetic drift, selection and the evolution of the mutation rate.  Nature, 17, 704–714.
  * Mitchell and Bridge [2006] Mitchell, D. and Bridge, R. (2006).  A test of Chargaff’s second rule.  Biochemical and Biophysical Research Communications, 340(1), 90–94.
  * Schrempf and Hobolth [2017] Schrempf, D. and Hobolth, A. (2017).  An alternative derivation of the stationary distribution of the multivariate neutral Wright-Fisher model for low mutation rates with a view to mutation rate estimation from site frequency data.  Theoretical Population Biology, 114, 88–94.
  * Stephens and Donnelly [2000] Stephens, M. and Donnelly, P. (2000).  Inference in Molecular Population Genetics.  Journal of the Royal Statistical Society, Series B, page Discussion Paper.
  * Vogl [2014] Vogl, C. (2014).  Estimating the scaled mutation rate and mutation bias with site frequency data.  Theoretical Population Biology, 98, 19–27.
  * Vogl and Clemente [2012] Vogl, C. and Clemente, F. (2012).  The allele-frequency spectrum in a decoupled Moran model with mutation, drift, and directional selection, assuming small mutation rates.  Theoretical Population Biology, 81(3), 197–209.
  * Wright [1931] Wright, S. (1931).  Evolution in Mendelian populations.  Genetics, 16, 97–159.