# Quantum oscillations in 2D electron gases with spin-orbit and Zeeman
interactions

Denis R. Candido  Department of Physics and Astronomy, University of Iowa,
Iowa City, Iowa 52242, USA  Sigurdur I. Erlingsson  Hamed Gramizadeh
Department of Engineering, Reykjavik University, Menntavegi 1, IS-102
Reykjavik, Iceland  João Vitor I. Costa  Instituto de Física de São Carlos,
Universidade de São Paulo, 13560-970 São Carlos, SP, Brazil  Pirmin J. Weigele
Department of Physics, University of Basel, CH-4056, Basel, Switzerland
Dominik M. Zumbühl  Department of Physics, University of Basel, CH-4056,
Basel, Switzerland  J. Carlos Egues  Instituto de Física de São Carlos,
Universidade de São Paulo, 13560-970 São Carlos, SP, Brazil  Department of
Physics, University of Basel, CH-4056, Basel, Switzerland

###### Abstract

Shubnikov-de Haas (SdH) oscillations are the fingerprint of the Landau and
Zeeman splitting level structure on the resistivity in presence of a moderate
magnetic field before full quantization is manifest in the integer quantum
Hall effect. These oscillations have served as a paradigmatic experimental
probe and tool for extracting key semiconductor parameters such as carrier
density, effective mass $m^{*}$, Zeeman splitting with g-factor $g^{*}$,
quantum scattering time and Rashba $\alpha$ and Dresselhaus $\beta$ spin-orbit
(SO) coupling parameters. Analytical descriptions of the SdH oscillations are
available for some special cases, but the generic case with all three terms
simultaneously present has not been solved analytically so far, seriously
hampering the analysis and interpretation of experimental data. Here, we
bridge this gap by providing an analytical formulation for the SdH
oscillations of 2D electron gases (2DEGs) with simultaneous and arbitrary
Rashba, Dresselhaus, and Zeeman interactions. We use a Poisson summation
formula for the density of states of the 2DEG, which affords a complete yet
simple description of the oscillatory behavior of its magnetoresistivity. Our
analytical and numerical calculations allow us to extract the beating
frequencies, quantum lifetimes, and also to understand the role of higher
harmonics in the SdH oscillations. More importantly, we derive a simple
condition for the vanishing of SO induced SdH beatings for all harmonics in
2DEGs: $\alpha/\beta=[(1-\tilde{\Delta})/(1+\tilde{\Delta})]^{1/2}$, where
$\tilde{\Delta}\propto g^{*}m^{*}$ is a material parameter given by the ratio
of the Zeeman and Landau level splitting. This condition is notably different
from that of the persistent spin helix at $\alpha/\beta=1$ for materials with
large $g^{*}$ such as InAs or InSb. We also predict beatings in the higher
harmonics of the SdH oscillations and elucidate the inequivalence of the SdH
response of Rashba-dominated ($\alpha>\beta$) vs Dresselhaus-dominated
($\alpha<\beta$) 2DEGs in semiconductors with substantial $g^{*}$. We find
excellent agreement with recent available experimental data of Dettwiler et
al. Phys. Rev. X 7, 031010 (2017), and Beukman et al., Phys. Rev. B 96, 241401
(2017). The new formalism builds the foundation for a new generation of
quantum transport experiments and spin-orbit materials with unprecedented
physical insight and material parameter extraction.

##  I Introduction

The spin-orbit (SO) interaction couples the orbital and spin degrees of
freedom, not only forms the basis for a range of spin related effects such as
the spin Hall effect D’Yakonov and Perel (1971); Dyakonov and Perel (1971);
Hirsch (1999); Landisman and Connors (2005) and the persistent spin helix
Schliemann  _et al._ (2003a); Bernevig  _et al._ (2006a); Fu  _et al._ (2016),
but also underlies the physical mechanisms of new phases of matter, e.g.,
topological insulators, quantum spin Hall materials Kane and Mele (2005);
Bernevig and Zhang (2006); Bernevig  _et al._ (2006b), and Majorana Kitaev
(2001); Fu and Kane (2009); Candido  _et al._ (2018), Dirac and Weyl fermions
Pal (2011). Accordingly, advancing techniques and methods to measure and
extract SO couplings from experimental data are crucial for the development of
these fields.

Figure 1: Magnetoresistivity for a) pure Rashba $\alpha=7.0$ meV nm and b)
pure Dresselhaus $\beta=7.0$ meV nm with ${m^{*}=0.019m_{o}}$ and
$n_{2D}=3.3\times 10^{11}$ cm-2 from Ref. Gilbertsson  _et al._ , 2008. The
curves in a) and b) are not the same due to the large $g$–factor $g^{*}=-34$.
The insets display the normalized FFT including the 2nd harmonic. The presence
of beating nodes in $\delta\rho_{xx}$ are clearly visible in a) the
fundamental and b) the 2nd harmonic, see Fig. 8. The condition for the absence
of beatings (single peak for each harmonic) is $\alpha=7.0$ meV nm and
$\beta=5.0$ meV nm, shown in c), but not $\alpha=\beta=7.0$ meV nm, the
persistant spin helix case, shown in d), clearly exhibiting a beating (here a
splitting of 2nd harmonic peak).

Shubnikov-de Haas (SdH) oscillationsShubnikov and De Haas (1930); Shubnikov
and de Haas (1930) are among the best techniques to probe simultaneously spin-
and charge-related quantities associated to electrons in semiconductors,
including effective masses, gyromagnetic ratios, quantum scattering times,
densities and SO couplings. Most recently, they have been crucial to the study
and understanding of new materials, as for example, 2D-materials, transition
metal dichalcogenides, van der Waals heterostructures Smoleński  _et al._
(2019); Kormányos  _et al._ (2015); Cui  _et al._ (2015); Xu  _et al._ (2016);
Cui  _et al._ (2017); Rhodes  _et al._ (2019); Masseroni  _et al._ (2023);
Slizovskiy  _et al._ (2023), and also materials hosting new phases of matter
e.g., topological insulators Ferreira  _et al._ (2022), unconventional
superconductivity Cao  _et al._ (2018a) and correlated insulator behavior Cao
_et al._ (2018b). It has also been used to establish the presence of nodal-
lines Hu  _et al._ (2016), Berry’s phase Murakawa  _et al._ (2013); Datta  _et
al._ (2019), and different topology of Fermi surfaces Alexandradinata  _et
al._ (2018). SdH oscillations are magneto-oscillations in the resistivity and
originate from the sequential crossings of the discrete Landau Levels (LLs)
through the Fermi energy. Without SO coupling and in the low-field regime, the
period of the SdH oscillations can be related to the density of the electron
gasIhn (2010). In the presence of SO interaction, on the other hand, the
energy spectrum changes dramatically thus leading to additional frequencies in
the magnetoresistivity and hence beatings, Figs. 1(a). This was first
theoretically described semiclassically by Das et. al., Das  _et al._ (1989).
In the so-called Onsager’s picture, different sub-bands possess different
Sommerfeld quantized orbits (playing the role of the LLs), which cross the
Fermi energy with different frequencies in $B^{-1}$. The spin-split bands give
rise to two distinct oscillating frequencies in the magnetotransport. The
standard experiment relies on Fourier analyzing the measured SdH oscillations.
An experimental method introduced in Refs. Nitta  _et al._ , 1997; Engels  _et
al._ , 1997; Schäpers  _et al._ , 1998 has often been used to estimate the
strength of the Rashba coupling via the splitting of the Fourier frequency
peaks. However, these methods have been criticized for not accounting for the
Zeeman splitting (through the g-factor $g^{*}$) nor for the additional
Dresselhaus SO coupling (Gilbertsson  _et al._ , 2008).

There have been some attempts to analyze the SdH oscillations taking into
account both $\alpha$, $\beta$ and $g^{*}$. However, these mostly involved
qualitative comparison with the energy spectrum of pure Rashba and pure
DresselhausAkabori  _et al._ (2006, 2008). In Ref. Yang and Chang, 2006, fully
numerical calculations of magneto-oscillations were performed but for
relatively high magnetic fields and low electron densities, far away from the
regime of recent experimental works Beukman  _et al._ (2017). Moreover, it was
realized that in the absence of the Zeeman interaction, important features are
absent. More specifically, without accounting for the spin mixing generated by
the magnetic field (via the Zeeman interaction), predictions become imprecise
Winkler  _et al._ (2000), and even fail to describe phenomena such as magnetic
inter-subband scattering Raikh and Shahbazyan (1994) and magnetic breakdown
Averkiev  _et al._ (2005). In general, full quantum mechanical numerics are
generally done in order to check agreement with experiments, which are neither
very practical nor elucidate much of the physics happening in those systems
Beukman  _et al._ (2017); Cimpoiasu  _et al._ (2019). Finally, all the
previous works have neglected the influence of higher harmonics, recently seen
experimentally Dettwiler  _et al._ (2017).

Here, we present a detailed investigation of SdH oscillations in the presence
of SO couplings of both Rashba $\alpha$ and Dresselhaus $\beta$ types and
Zeeman interaction with g-factor $g^{*}$. Our main result is the derivation,
for the first time in the literature, of a simple analytical expression for
the SdH oscillations in the presence of simultaneous arbitrary couplings
$\alpha$ and $\beta$ in addition to $g^{*}$. We note that earlier analytical
descriptions of SdH magnetoresistivity oscillations considered particular
cases, namely, when either only one of the parameters $\alpha$, $\beta$ or
$g^{*}$ was nonzero or any two of these parameters were nonzero, with the
exceptions ($\alpha\neq\beta\neq 0$, $g^{*}=0$) and ($\alpha=\beta$,
$g^{*}\neq 0$).

Interestingly, our analytical formula generalizes previous resultsAverkiev
_et al._ (2005) for $g^{*}=0$ and predicts a new condition for the vanishing
of the SdH magneto-oscillation beatings in all harmonics [e.g., Figs. 1(c)] in
Rashba-Dresselhaus coupled 2DEGs with substantial Zeeman splittings, namely,

$\frac{\alpha}{\beta}=\sqrt{\frac{1-\tilde{\Delta}}{1+\tilde{\Delta}}},$ (1)

where $\tilde{\Delta}$ is a material parameter given by the ratio between the
Zeeman splitting and the Landau-level spacing. As we discuss later on, Eq. (1)
is not associated with a conserved quantity in our system; this contrasts with
the persistent-spin-helix condition $\alpha=\beta$, which predicts spin
conservation along particular axes Schliemann  _et al._ (2003a); Bernevig  _et
al._ (2006a); Fu  _et al._ (2016). We stress that this case with
$\alpha=\beta$ and generic $g^{*}\neq 0$ leads to beating in the frequency
spectrum of our system, Figs. 1(d), as opposed to our new condition in Eq. 1.
As we discuss below, our numerical and analytical approaches show excellent
agreement with available data from Refs. Beukman  _et al._ , 2017 and
Dettwiler  _et al._ , 2017.

Our approach combines a semiclassical formulation for the resistivity of 2DEGs
with a trace formula for the density of states (DOS) in a quantizing magnetic
field. The trace formula expresses the DOS using the usual Poisson summation
formulaBrack and Bhaduri (1997). This formulation brings out the oscillatory
part of the DOS, thus allowing us to clearly identity the higher harmonics of
the SdH oscillations. It enables us to conveniently separate the frequency
scales into “fast” and “slow” oscillations thus allowing for a clearer
interpretation of the underlying physical phenomena, e.g., the slow beating
SdH oscillations due to the SO coupling.

Our main results for the oscillatory part of magnetoresistivity
$\delta\rho_{xx}(1/B)$ and its frequency spectra $I(f)$ [panel insets] are
show in Fig. 1. For pure Rashba [$\alpha\neq 0$, $\beta=0$, Fig. 1a)] and pure
Dresselhaus [$\alpha=0$, $\beta\neq 0$, Fig. 1b)], but non-zero Zeeman term
($g^{*}\neq 0$), the frequency spectra, as usual, show two main peaks, which
correspond to the first two Fourier components of $\delta\rho_{xx}(1/B)$.
These two cases, however, exhibit a marked contrast: while the pure Rashba
shows a peak splitting at the fundamental frequency, the pure Dresselhaus
exhibits a peak splitting in the second harmonic. As we explain in detail in
Sec. V.4, this contrasting behavior arises from the interplay between the
Zeeman and SO interactions, which makes the SdH magneto-responses with nonzero
g-factors $g^{*}$ inequivalent for Rashba-dominated ($\alpha>\beta$) vs.
Dresselhaus-dominated ($\alpha<\beta$) 2DEGs. For $g^{*}=0$, the pure Rashba
and pure Dresselhaus cases give identical results.

Figure 1(c) illustrates our prediction in Eq. (1) thus showing no peak
splitting in the frequency spectra – at any harmonic – when this condition is
satisfied. To emphasize this condition emulates a situation with no SO
coupling (i.e., no beating), we plot in Fig. 1(c) the $\alpha=\beta=0$ (with
$g^{*}\neq 0$) case [dashed curve in 1(c)], which shows complete overlap with
the case satisfying Eq. (1). In contrast and for completeness, Fig. 1d) shows
the $\alpha=\beta\neq 0$ case with $g^{*}\neq 0$, which exhibits peak
splitting in the second harmonic.

We have applied our analytical description to low-density GaAs-based quantum
wells for which there are experimental dataDettwiler  _et al._ (2017) showing
several harmonics in the SdH magneto-oscillations. Figure 2 shows the
excellent agreement obtained, thus illustrating that our semiclassical
formulas can satisfactorily capture the higher harmonics of the SdH
oscillations. Moreover, we have applied our analytical approach to low-density
InSb-based 2DEGsGilbertsson  _et al._ (2008); Akabori  _et al._ (2008) where,
unlike GaAs-based 2DEGs, a strong SO coupling manifests itself as beatings in
the measured SdH oscillations, and find good agreement. We have also
implemented a detailed numerical calculation for high-density InAs-based 2DEGs
for which an analytical description is not adequate. Here again we find very
good agreement with available dataBeukman  _et al._ (2017) and are able to
extract SO coupling parameters.

Next (Sec II), we present a description of the Hamiltonian of our system. In
Sec. III we discuss how to obtain the “$F$–function”, the central quantity in
our formulation, from the Landau-quantized energy spectrum of our system and
its connection with the density of states (DOS). The formalism for obtaining
the Shubnikov-de Haas oscillations in terms of the Poisson summation formula
and the F-function is described in Sec. IV. Finally, in Sec. V we present and
analyze different particular cases of SdH oscillations and, more important,
derive the new condition in Eq. (1) for the complete absence of beatings (all
harmonics) in the SdH oscillations, for 2DEGs with non-zero Rashba,
Dresselhaus, and Zeeman couplings. The appendices present relevant details of
our theoretical formulation.

##  II 2DEG Hamiltonian

Our starting point is the Hamiltonian for a 2DEG confined in a quantum well
($xy$ plane) grown along the [001] crystallographic direction, taken as $z$
axis. In the presence of a perpendicular external magnetic field
$\textbf{B}=(0,0,B)$ and both Rashba(Bychkov and Rashba, 1984) and
DresselhausDresselhaus, 1955 spin orbit interactions, the Hamiltonian reads

$\displaystyle{\cal H}$ $\displaystyle=$
$\displaystyle\frac{1}{2m^{*}}\left(\Pi_{x}^{2}+\Pi_{y}^{2}\right)+\frac{1}{2}{g^{*}\mu_{B}B}{}\sigma_{z}$
(2) $\displaystyle+$
$\displaystyle\frac{\alpha}{\hbar}\left(\Pi_{x}\sigma_{y}-\Pi_{y}\sigma_{x}\right)+\frac{\beta}{\hbar}\left(\Pi_{x}\sigma_{x}-\Pi_{y}\sigma_{y}\right),$

where $g^{*}$ is the g-factor, $m^{*}$ is effective mass,
$\bm{\Pi}=\textbf{p}-q\textbf{A}$ is the canonical momentum, $q$ is the
electric charge, $\mu_{B}$ is the Bohr magneton, $\hbar$ the reduced Planck’s
constant and $\sigma_{x}$, $\sigma_{y}$, $\sigma_{z}$ denote the usual Pauli
matrices. The parameters $\alpha$ and $\beta$ denote the linear-in-$k$ Rashba
and Dresselhaus SO couplings, respectively. The $\beta$ coupling includes a
density dependent correction arising from the cubic Dresselhaus term. As we
discuss later on [Sec. VI.1], our numerical results will account for the full
cubic Dresselhaus term.

Let us introduce the annihilation and creation operators associated to the
Landau level $\left|n\right\rangle$

$\displaystyle a$ $\displaystyle=$
$\displaystyle\frac{\ell_{c}}{\sqrt{2}\hbar}\left(\Pi_{x}-i\zeta\Pi_{y}\right),$
(3) $\displaystyle a^{\dagger}$ $\displaystyle=$
$\displaystyle\frac{\ell_{c}}{\sqrt{2}\hbar}\left(\Pi_{x}+i\zeta\Pi_{y}\right),$
(4)

obeying $\left[a,a^{\dagger}\right]=1$,
${a\left|n\right\rangle=\sqrt{n}\left|n-1\right\rangle}$,
${a^{\dagger}\left|n\right\rangle=\sqrt{n+1}\left|n+1\right\rangle}$,
$\zeta=-{\rm sign}\left(qB\right)$, with the magnetic length and the center of
the Landau orbit denoted by $\ell_{c}=\sqrt{\frac{\hbar}{\left|qB\right|}}$
and $y_{0}=\frac{ek_{x}}{|qB|}$, respectively. In this work, we have $q=-e$,
where $e>0$ is the absolute value of the elementary electronic charge, and we
choose $B>0$, yielding $\zeta=1$. Using Eqs. (3) and (4), our Hamiltonian [Eq.
(2)] becomes

$\displaystyle{\cal H}$ $\displaystyle=$
$\displaystyle\hbar\omega_{c}(a^{\dagger}a+1/2)+\frac{\Delta}{2}\sigma_{z}-\frac{i\alpha}{\sqrt{2}\ell_{c}}(a^{\dagger}\sigma_{-}-a\sigma_{+})$
(5) $\displaystyle+$
$\displaystyle\frac{\beta}{\sqrt{2}\ell_{c}}(a^{\dagger}\sigma_{+}+a\sigma_{-}),$

where the cyclotron frequency is $\omega_{c}=eB/m^{*}$,
${\Delta=g^{*}\mu_{B}B}$, which inherits its sign from $g^{*}$, and
$\sigma_{\pm}=\sigma_{x}\pm i\sigma_{y}$, with $\sigma_{x}$ and $\sigma_{y}$
denoting Pauli matrices. We now perform the canonical transformation
$\tilde{{\cal H}}={\cal U}{\cal H}{\cal U}^{\dagger}$ with ${\cal
U}=e^{-i\frac{\pi}{4}\left(\frac{\sigma_{z}}{2}+a^{\dagger}a\right)}$, which
yields

$\displaystyle{\cal U}\sigma_{\pm}{\cal U}^{\dagger}$ $\displaystyle=$
$\displaystyle\sigma_{\pm}e^{\mp i\frac{\pi}{4}},$ (6) $\displaystyle{\cal
U}\sigma_{z}{\cal U}^{\dagger}$ $\displaystyle=$ $\displaystyle\sigma_{z},$
(7) $\displaystyle{\cal U}a^{\dagger}{\cal U}^{\dagger}$ $\displaystyle=$
$\displaystyle e^{i\frac{\pi}{4}}a^{\dagger},\qquad$ (8)

and finally

$\displaystyle\frac{{\cal\tilde{H}}}{\hbar\omega_{c}}$ $\displaystyle=$
$\displaystyle(a^{\dagger}a+1/2)+\frac{\tilde{\Delta}}{2}\sigma_{z}+\alpha_{B}(a^{\dagger}\sigma_{-}+a\sigma_{+})$
(9) $\displaystyle+\beta_{B}(a^{\dagger}\sigma_{+}+a\sigma_{-}),$

where we have introduced the real valued, dimensionless quantities
$\alpha_{B}=\frac{\alpha}{\sqrt{2}\hbar\omega_{c}\ell_{c}}$,
$\beta_{B}=\frac{\beta}{\sqrt{2}\hbar\omega_{c}\ell_{c}}$ and
$\tilde{\Delta}=\frac{\Delta}{\hbar\omega_{c}}=\frac{g^{*}m^{*}}{2m_{0}}$.

Analytical solutions for the above Hamiltonian [Eq. (9)] can be found for the
cases with either pure Rashba or pure DresselhausBychkov and Rashba, 1984;
Winkler, 2003. The specific cases of $\alpha=\pm\beta$ turn out to be of great
physical interest, where persistent spin helix (PSH)Schliemann  _et al._
(2003a); Bernevig  _et al._ (2006a); Dettwiler  _et al._ (2017) and persistent
Skyrmion lattice (PSL) Fu  _et al._ (2016) were predicted. Interestingly, the
case with $\alpha=\pm\beta$ maps to the Rabi model in quantum optics and was
recently solved exactly (Braak, 2011). The exact solution relies on obtaining
zeros of a transcendental function. Moreover, previous studies of the Rabi
model have important implications for our system. For instance, we have shown
that the Rabi parity symmetry (Casanova  _et al._ , 2010; Braak, 2011) remains
valid in our problem for arbitrary $\alpha$ and $\beta$ (See Appendix B). This
enables us to separate the Hilbert space in two subspaces with different
parities, which can be individually analyzed and compared. As for general
couplings $\alpha$ and $\beta$, similar systems have been studied before in
the framework of Landau levels, using either variational (Hartree-Fock)
methods (Schliemann  _et al._ , 2003b), second order perturbation (Zarea and
Ulloa, 2005; Wang and Vasilopoulos, 2005) or obtaining the spectrum in terms
of solutions of transcendental equations (Zhang, 2006). A perturbation scheme
based on 4th order Schrieffer-Wolff transformation has also been used to find
approximate analytical solutions (Erlingsson  _et al._ , 2010). However, we
are unaware of any exact analytical solution for general Rashba, Dresselhaus
and Zeeman coupling.

##  III $F$-function and its connection with the energy spectrum and DOS

For our 2DEG in the presence of perpendicular magnetic field, the low magnetic
field regime corresponds to having a very large number of Landau levels below
the Fermi energy $\varepsilon_{F}$ (taken as constant and equal to its zero-
field value), i.e., many occupied states. The system is thus assumed to be far
away from the integer quantum Hall regime where few Landau levels are occupied
and the effects of electron-electron interaction become important Ihn (2010).
Let us denote the eigenenergies of our dimensionless Hamiltonian Eq. (9) by
$\varepsilon_{n,s}$, where $n\in\mathbb{N}_{0}$ represents the LL number and
$s=\pm$ represents a pseudo-spin associated to the presence of two spin-split
bands (due to the Zeemann and SO interactions). With this notation, the
density of states (DOS) reads

$\mathcal{D}(\varepsilon,B)=\frac{{{\color[rgb]{0,0,0}\tilde{D}}}}{A}\sum_{n,s}\delta(\varepsilon-\hbar\omega_{c}\varepsilon_{n,s}),$
(10)

where $\tilde{D}=A/2\pi\ell_{c}^{2}$ is the LL degeneracy and $A$ the 2DEG
area. This LL degeneracy is the same for all 2DEGs studied here in the
presence or absence of Zeemann and SO interactions.

As we show in the next section, the magnetotransport properties of the system
can be determined by the Landau levels sequentially crossing
$\varepsilon_{F}$. The rate at which these crossings happen will determine a
periodic behavior of the magnetotransport properties of the system as the
magnetic field is varied. In order to describe this periodicity, we introduce
the $F$-function Ihn (2010) (see Appendix A for details), which is defined by
the relation

$\varepsilon_{n,s}(B)=\varepsilon\leftrightarrow n=F_{s}(\varepsilon,B).$ (11)

The $F$-function gives the Landau level index $n$ of the state that has energy
$\varepsilon$ and pseudo-spin $s$ at magnetic field $B$.Brack and Bhaduri
(1997); Tarasenko (2002) It is important to notice that the equation for $n$
[Eq. (11)] can also return non-integer values for $n$. In such cases the
$F$-function provides a measure of how close a Landau level $n$ is to the
energy $\varepsilon$, for a given pseudo-spin $s$ and magnetic field $B$.

Since one can relate transport phenomena with the density of states, we
rewrite the DOS of our system in a way that highlights its oscillatory
behavior dependence on both $\alpha$ and $\beta$. First we introduce the
$F_{s}$ function into Eq. (10)

$\displaystyle\mathcal{D}(\varepsilon,B)$
$\displaystyle\approx\frac{m^{*}}{2\pi\hbar^{2}}\sum_{n,s}\delta(n-F_{s}(\varepsilon/\hbar\omega_{c},B)),$
(12)

which neglects terms $\mathcal{O}[(\alpha
m^{*}\ell_{c})^{2}/\hbar]+\mathcal{O}[(\beta m^{*}\ell_{c})^{2}/\hbar)]$. This
holds for typical values of $\alpha$, $\beta$, $m^{*}$ and small magnetic
fields $B\lesssim 1$T. Using the Poisson summation formula
$\sum_{n=0}^{\infty}\delta(n-F_{s})=1+2\sum_{l=1}^{\infty}\cos(2\pi lF_{s})$
and defining the relevant quantities

$\mathcal{F}_{\pm}=\frac{1}{2}(F_{+}\pm F_{-}),$ (13)

we obtain

$\displaystyle\frac{\mathcal{D}(\varepsilon,B)-{{\color[rgb]{0,0,0}2}}\mathcal{D}_{0}}{{{\color[rgb]{0,0,0}2}}\mathcal{D}_{0}}$
$\displaystyle\approx$ $\displaystyle 2\sum_{l=1}^{\infty}\cos(2\pi
l\mathcal{F}_{+})\cos(2\pi l\mathcal{F}_{-}),$ (14)

where $\mathcal{D}_{0}=\frac{m^{*}}{2\pi\hbar^{2}}$ is the (constant) density
of states per spin for the 2DEG at zero magnetic field (see Appendix A for
details). As we are going to see later, $\mathcal{F}_{+}$ contains the fast
oscillations with respect to $1/B$, which is is proportional to the electron
density $n_{2D}$. On the other hand, $\mathcal{F}_{-}$ contains the slow
oscillations that are determined by the spin-orbit coupling terms, $\alpha$
and $\beta$. Moreover, the fast oscillations arising from the terms with $l>1$
correspond to the higher harmonics, and have be seen in experiments Dettwiler
_et al._ (2017).

##  IV SdH oscillations in the magnetoresisitivity

As already mentioned, the oscillations in the magnetoresistivity as a function
of the magnetic field are called SdH oscillations(Ihn, 2010). They appear as a
consequence of the sequential depopulation of the LLs when the magnetic field
is increased. For low magnetic fields where multiple LL are occupied, i.e.,
far from the integer quantum Hall regimeIhn (2010), a semi-classical
description of the magneto-oscillations can be used.

In experiments, the measurement of the SdH oscillations is accessed via the
longitudinal differential resistivity. In general, the resistivity tensor is
defined as the inverse matrix of the conductivity tensor,

$\bm{\rho}=\left(\begin{array}[]{cc}\sigma_{xx}&\sigma_{xy}\\\
\sigma_{xy}&\sigma_{xx}\end{array}\right)^{-1}.$ (15)

The relevant magnetoresistivity component for us is

$\rho_{xx}=\frac{\sigma_{xx}}{\sigma_{xx}^{2}+\sigma_{xy}^{2}},$ (16)

where

$\sigma_{xx(xy)}(B,T)=\int
d\varepsilon\left(-\frac{df_{0}(\varepsilon)}{d\varepsilon}\right)\sigma_{xx(xy)}(B,\varepsilon,T=0),$
(17)

where $f_{0}(\varepsilon)$ is the Fermi-Dirac distribution. Using a semi-
classical approach, we account for the magnetic field dependence of the
conductivity via the electron scattering time $\tau(\varepsilon,B)$, which is
proportional to the DOS ${\cal D}(\varepsilon,B)$ via Fermi’s golden rule.
Accordingly, up to linear order on the deviation of the DOS, we obtain

$\tau(\varepsilon,B)\approx\tau_{0}\left(\varepsilon\right)\left[1-\frac{{\cal
D}(\varepsilon,B)-{\cal D}_{0}(\varepsilon)}{{\cal
D}_{0}(\varepsilon)}\right],$ (18)

with ${\cal D}_{0}(\varepsilon)={\cal D}(\varepsilon,B=0)$ and
$\tau_{0}(\varepsilon)=\tau(\varepsilon,B=0)$. Using the Drude semi-classical
equations for the frequency-independent currentIhn (2010), the normalized
longitudinal resistivity reads

$\displaystyle\delta\rho_{xx}(B)$
$\displaystyle=\frac{\rho_{xx}(B)-\rho_{xx}(B=0)}{\rho_{xx}(B=0)}$ (19)
$\displaystyle=\int
d\varepsilon\left(-\frac{df_{0}(\varepsilon)}{d\varepsilon}\right)\frac{{\cal
D}(\varepsilon,B)-{\cal D}_{0}(\varepsilon)}{{\cal D}_{0}(\varepsilon)}.$ (20)

For the DOS in the presence of Landau level broadening due to scattering
processes, the relation in Eq. (10) is replaced by

$\mathcal{D}(\varepsilon,B)=\frac{{{\color[rgb]{0,0,0}\tilde{D}}}}{A}\sum_{n,s}L_{\Gamma}(\varepsilon-\hbar\omega_{c}\varepsilon_{n,s}),$
(21)

where $L_{\Gamma}(x)$ describes the broadening function, e.g., Lorentzian or
Gaussian, and $\Gamma$ is parameter defining the broadening of the levels (see
Appendix A for details). After applying the Poisson summation formula, we
obtain a result that resembles Eq. (14), apart from the appearance of the the
cosine Fourier transform of $L_{\Gamma}(x)$, denoted with
$\tilde{L}_{\Gamma}(x)$,

$\frac{{\cal D}(\varepsilon,B)-{\cal D}_{0}(\varepsilon)}{{\cal
D}_{0}(\varepsilon)}\approx
2\sum_{l=1}^{\infty}\tilde{L}_{\Gamma}\left(l\frac{\Gamma}{\hbar\omega_{c}}\right)\cos(2\pi
l\mathcal{F}_{-})\cos(2\pi l\mathcal{F}_{+}).$ (22)

The so-called Dingle factor $\tilde{L}_{\Gamma}(x)$ Ihn (2010) sets the limit
of validity of the semi-classical approximation, i.e., that the oscillatory
part of the resistivity should be much smaller than the constant term. It also
gives the regime where it is valid to consider only the lowest harmonic.
Higher harmonics have been observed in magnetoresistivity measurements
Dettwiler  _et al._ (2017) in GaAs-based 2DEGs. The $\mathcal{F}_{-}$ function
can be related to the envelope of the SdH oscillations. The general form of
the temperature-dependent normalized resistivitity reads

$\displaystyle\delta\rho_{xx}(B,T)$ $\displaystyle=2\sum_{l=1}^{\infty}\int
d\varepsilon{\tilde{L}_{\Gamma}\left(l\frac{\Gamma}{\hbar\omega_{c}}\right)}\left(-\frac{df_{0}(\varepsilon,T)}{d\varepsilon}\right)$
(23) $\displaystyle\times\cos(2\pi l\mathcal{F}_{-})\cos(2\pi
l\mathcal{F}_{+}).$

Even though we only consider the zero-temperature limit in the present work,
for completeness, below we present the temperature-dependence of
$\delta\rho_{xx}(B,T)$ valid in the relevant parameter range considered in
this work and for all the systems studied here. As show in Appendix G, we find

$\displaystyle\delta\rho_{xx}(B,T)$ $\displaystyle\approx$ $\displaystyle
2\sum_{l=1}^{\infty}{\tilde{L}_{\Gamma}\left(l\frac{\Gamma}{\hbar\omega_{c}}\right)}{\cal
A}_{l}(T)$ (24) $\displaystyle\times\left.\cos(2\pi l\mathcal{F}_{-})\cos(2\pi
l\mathcal{F}_{+})\right|_{\varepsilon=\varepsilon_{F}},$

where the temperature-dependent coefficient

${\cal
A}_{l}(T)=\frac{2\pi^{2}lk_{B}T/\hbar\omega_{c}}{\sinh\left(2\pi^{2}lk_{B}T/\hbar\omega_{c}\right)}$
(25)

accounts for the temperature dependence of the SdH oscillations. In the limit
of vanishing $\alpha$ and $\beta$ this reduces to the result in Ref. Laikhtman
and Altshuler, 1994, and in the case of both vanishing $\beta$ and broadening
($\Gamma=0$) gives the result in Ref. Winkler, 2003 [Eq. (9.28)]. Here we
assume that $\varepsilon$ is close to the zero magnetic field Fermi energy
$\varepsilon_{F}=\hbar^{2}k_{F}^{2}/2m$.

A widely used method to extract spin-orbit couplings and electronic densities
is to analyze the oscillations by calculating the quantity

$I(f)=\left|\int_{B_{2}^{-1}}^{B_{1}^{-1}}d\left(\frac{1}{B}\right)\frac{\rho_{xx}(B)-\rho_{xx}(B_{1})}{\rho_{xx}(B_{1})}e^{i2\pi
f/B}\right|^{2},$ (26)

which defines the power spectrum of the normalized magneto-resistivity with a
trivial background value $\rho_{xx}(B=B_{1})$ removed. Note that $B_{1}$
should be small enough such that the semiclassical regime of a constant
$\rho_{xx}(B\rightarrow 0)$ is reached.

Figure 2: The power spectrum $I(f)$ for $\delta\rho_{xx}$ measurements on a
GaAs 2DEG in Ref. Dettwiler  _et al._ , 2017 obtained using Eq. (26). The
calculated results used Eq. (30) with one fitting parameter: $\tau_{q}$. The
inset shows the magneto-resistivity data and the corresponding calculated
$\delta\rho_{xx}$.

In Fig. 2 the power spectrum is shown for data from Fig. S11a in Ref.
Dettwiler  _et al._ , 2017, where magnetoresistivity SdH oscillations were
measured in a GaAs 2DEG over a magnetic field interval $[0.20,1.5]$ T. The
power spectrum shows a SdH peak at $f\approx 10.5$ T (the fundamental
frequency), and higher harmonics are clearly visible at $21.0$ T and $31.5$ T,
corresponding to the first and second harmonic, respectively. The experimental
data was fitted with Eq. (30) with one fit parameter: $\tau_{q}$. The
resulting fit matches very well the harmonics of the SdH signal. To acccount
for the small background shift in the experimental data as seen in the inset a
more elaborate modeling of the data would be required. The fitting was done
using six harmonics, and resulted in $\tau_{q}=0.8$ ps, using standard GaAs
parameters $m=0.067{m_{0}}$ and $g^{*}=-0.44$. Note that we have used Eq.(30),
which does not include SO coupling, for our fitting procedure here. This is
justifiable because GaAs-based 2DEGs have relatively small SO couplings, not
accessible via SdH measurements. Weak anti-localization measurements can
access the SO parameter in these systemsBeukman  _et al._ (2017). However,
GaAs-based 2DEGs have relatively high mobilities thus making it possible to
see many harmonics.

##  V Results and discussions

In this section we present the energy spectrum, $F$–function and
magnetoresistivity SdH oscillations for different parameter regimes of our
Hamiltonian, Eq. (9). Additionally, we discuss in detail the interpretation of
the SdH oscillations within the trace formula description (e.g., contribution
of higher harmonics) and show how to extract relevant spin-orbit couplings
from it. The results are presented in order of simplicity, i.e., from the
simplest to the more complex case.

###  V.1 Landau levels with only Zeeman interaction

In the presence of Zeeman and no Rashba and Dresselhaus SO couplings, i.e.,
$\alpha=\beta=0$, the eigenenergies of our Hamiltonian [Eq. (9)] are given by

$\frac{\varepsilon_{n,s}}{\hbar\omega_{c}}=n+\frac{1}{2}+\frac{\tilde{\Delta}}{2}s.$
(27)

with $n\in\mathbb{N}_{0}$ and $s=1$ ($s=-1$) representing the pure spin state
$\left|\uparrow\right\rangle$ ($\left|\downarrow\right\rangle$). In Fig. 3 we
plot the four energy levels corresponding to $n=4,5$ and $s=\pm 1$, along with
$\varepsilon_{F}/\hbar\omega_{c}$, using the following InSb QW parameters from
Refs. Gilbertsson  _et al._ , 2008; Akabori  _et al._ , 2008:
$m^{*}=0.019m_{o}$, $g^{*}=-34$ and electron density $n_{2D}=3.3\times
10^{-3}$ nm-2. For these parameters, the ordering of the energies obeys
$\varepsilon_{n+1,-1}>\varepsilon_{n+1,1}>\varepsilon_{n,-1}>\varepsilon_{n,1}$.
Figure 3 shows how successive levels cross the Fermi energy as a function of
the magnetic field. This, in turn, will reflect on the oscillations of the
resistivity once for $\varepsilon_{F}\approx\varepsilon_{n,s}$, an increase on
the conductivity will happen due to the resonance condition between the
corresponding LL and the Fermi energy.

From the energy expressions above [Eq. (27)], we can obtain the $F$–functions
through Eq. (11), namely,

${F_{s}\left(\varepsilon\right)=\frac{\varepsilon}{\hbar\omega_{c}}-\frac{\tilde{\Delta}}{2}s-\frac{1}{2}},\text{
with }\frac{dF_{s}(\varepsilon)}{d\varepsilon}=\frac{1}{\hbar\omega_{c}},$
(28)

yielding the fast and slow components [Eq. (13)]

$\mathcal{F}_{+}\left(\varepsilon,B\right)=\frac{\varepsilon}{\hbar\omega_{c}}-\frac{1}{2},\qquad\mathcal{F}_{-}\left(\varepsilon,B\right)=-\frac{\tilde{\Delta}}{2}.$
(29)

At $\varepsilon=\varepsilon_{F}$ these can be expressed (to a very good
approximation) as $\mathcal{F}_{+}=\frac{hn_{2D}}{2e}\frac{1}{B}-\frac{1}{2}$
and $\mathcal{F}_{-}=-\frac{g}{4}\frac{m^{*}}{m_{0}}$, where we assume that
$n_{2D}=\frac{k_{F}^{2}}{2\pi}$ is the 2DEG electron density at $B=0$.

Figure 3: Landau levels $n=4$, 5 [Eq. (27)] as a function of $1/B$ for a 2DEG
with only Zeeman interaction and no SO couplings ($\alpha=\beta=0$). The
dotted line shows $\varepsilon_{F}/\hbar\omega_{c}$. Here, we use
$m^{*}=0.019m_{o}$, $g^{*}=-34$ and $n_{2D}=3.3\times 10^{-3}$ nm-2 for InSb-
based wellsGilbertsson  _et al._ (2008); Akabori  _et al._ (2008).

The corresponding resistivity can now be determined through Eq. (22) and reads

$\delta\rho_{xx}(B)=2\sum_{l=1}^{\infty}e^{-l\pi\frac{\hbar/\tau_{q}}{\hbar\omega_{c}}}\frac{2\pi^{2}lk_{B}T/\hbar\omega_{c}}{\sinh\left(2\pi^{2}lk_{B}T/\hbar\omega_{c}\right)}\cos\left[2\pi
l\left(\frac{f^{\rm SdH}}{B}-\frac{1}{2}\right)\right]\cos\left(\pi
l{g^{*}}\frac{m^{*}}{m_{0}}\right),$ (30)

where $f^{\rm SdH}=\frac{hn_{2D}}{2e}$ and we have assumed a Lorentzian form
for the $L_{\Gamma}$ broadening. For small magnetic fields, both effective
mass and $g$-factor nominal values do not depend on the magnetic field Lei
_et al._ (2020). As a result, the $1/B$-dependence of the resistivity in a
2DEG with only Zeeman coupling, displays oscillations with frequencies
multiple of $f^{\rm{SdH}}$, and absence of beating. This can be seen from Fig.
4, where we plot $\delta\rho_{xx}(B)$ vs $1/B$ for the harmonics $l=1,2,3$ and
clearly see oscillations with the respective frequencies
$f^{\rm{SdH}}$,$2f^{\rm{SdH}}$, and $3f^{\rm{SdH}}$. The solid (dotted) curves
correspond to $g^{*}=-34$ and $m^{*}=0.019m_{o}$ ($g^{*}=0$ and
$m^{*}=0.019m_{o}$) Gilbertsson  _et al._ (2008); Akabori  _et al._ (2008).
Note that the higher harmonics have smaller resistivity amplitudes. This
occurs due to the Dingle factor $\propto e^{-l/B}$, which suppresses the
higher harmonic components.

Figure 4: Magnetoresistivity deviation $\delta\rho_{xx}(B)$ as a function of
$1/B$ for a 2DEG with only Zeeman coupling and no SO couplings. The lowest
curve corresponds to $\delta\rho_{xx}(B)$ and the curves labelled by $l$ are
the individual frequency components in Eq. (30). The solid (dashed) line
corresponds to $g^{*}=-34$ ($g^{*}=0$), ${m^{*}=0.019m_{o}}$ and
$n_{2D}=3.3\times 10^{-3}$ nm-2. These parameters are for InSb-based wells
Gilbertsson  _et al._ (2008); Akabori  _et al._ (2008).

We should stress that the effects of the Zeeman coupling within the plot of
$\delta\rho_{xx}(B)$ are not immediately obvious. For instance, it can be seen
that for $g^{*}=0$ and $g^{*}\neq 0$, the corresponding
$\delta\rho_{xx}^{l=1}(B)$ (blue curves depicting the first harmonic) only
differ from themselves by the amplitude of the oscillation. For
$\tilde{\Delta}=-0.323$, $\cos(2\pi\tilde{\Delta}/2)$ is smaller than one,
thus yielding a reduction of the total amplitude for $g^{*}\neq 0$ as compared
to $g=0$. As a consequence, the presence of Zeeman is not readily evident from
the oscillations of $\delta\rho_{xx}^{l=1}(B)$. Conversely, the Zeeman is only
manifested within the resistivity when one considers many harmonics, as we
discuss below.

The definition of DOS in Eq. (21) gives broadened Landau levels separated by
$\hbar\omega_{c}$ , which are in turn spin split by the Zeeman term
$\tilde{\Delta}$ [See Eq. (27) and Fig. 3]. This spin splitting can only be
seen in the resistivity [Eq. (30)] when the contributions from the first and
second harmonics, $\cos(2\pi f^{\rm{SdH}}/B-\pi)\cos(2\pi\tilde{\Delta}/2)$
and $\cos(4\pi f^{\rm{SdH}}/B-2\pi)\cos(4\pi\tilde{\Delta}/2)$, respectively,
have opposite signs. For the parameters of Fig. 3 $\tilde{\Delta}=-0.323$ the
Zeeman term significantly affects the maximum of the resistivity. This can be
seen in Fig. 4, where the resistivities associated to harmonics $l=1$ and
$l=2$ (blue and cyan solid curves, respectively), interfere in a destructive
way, producing the double-peak feature in the total resistivity (purple solid
lines), characteristic of the incipient spin splitting in such data. We
emphasize, however, that this feature can be absent depending on the
broadening of the energy levels (due to the overlap of the spin-split levels).
This is the reason why the double-peak feature is not seen on the other
maximum peaks.

Although the $g^{*}$-factor term does not depend explicitly on magnetic field,
it can manifest itself in the magneto-oscillations. More specifically, Zeeman-
only effects can have a pronounced effect on the magneto-oscillations,
controlling the amplitude and sign of how subsequent harmonics are added,
either constructively or destructively, before being damped by the quantum
life time. Furthermore, it is important to say that the Zeeman can give rise
to interesting features and affect drastically the understanding of the
magneto-oscillations. For instance, if one could engineer a material 111InAs
has the value of $\tilde{\Delta}\approx 0.24$, but one would have a symmetric
structure to minimize the spin-orbit contribution such that
$\tilde{\Delta}=\frac{g^{*}}{2}\frac{m^{*}}{m_{o}}=0.5+m$ with
$m\in\mathbb{Z}$, then the main weight of the resistivity would be due to the
second harmonic with SdH frequently $2f^{\rm SdH}$ as
$\cos(l\pi\tilde{\Delta})=0$ for $l=1$.

###  V.2 Landau Levels with Zeeman and Rashba interactions

We now analyze the case where we have the presence of both Zeeman and Rashba
terms, i.e., $\tilde{\Delta}\neq 0$, $\alpha\neq 0$ and no Dresselhaus
coupling $\beta=0$ in Eq. (9). In the spin basis
$\left\\{\left|\uparrow\right\rangle,\thinspace\left|\downarrow\right\rangle\right\\}$,
the corresponding Hamiltonian assumes the following matrix form

$\frac{{\cal\tilde{H}}}{\hbar\omega_{c}}=\left(\begin{array}[]{cc}a^{\dagger}a+\frac{1}{2}+\frac{\tilde{\Delta}}{2}&2\alpha_{B}a\\\
2\alpha_{B}a^{\dagger}&a^{\dagger}a+\frac{1}{2}-\frac{\tilde{\Delta}}{2}\end{array}\right).$
(31)

Interestingly, the operator ${\cal N}_{+}=a^{\dagger}a+\sigma_{z}/2$ commutes
with the Hamiltonian above, i.e., ${[{\cal\tilde{H}},{\cal N}_{+}]=0}$, and
hence ${\cal\tilde{H}}$ and ${\cal N}_{+}$ share the same eigenstates. Hence
we have ${{\cal
N}_{+}\left|n,\uparrow\right\rangle=\left(n+1/2\right)\left|n,\uparrow\right\rangle}$
and ${{\cal
N}_{+}\left|n+1,\downarrow\right\rangle=\left(n+1/2\right)\left|n+1,\downarrow\right\rangle}$
, i.e., for $n\in\mathbb{N}$, $\left|n,\uparrow\right\rangle$ and
$\left|n+1,\downarrow\right\rangle$ are degenerate with respect to the
operator ${\cal N}_{+}$, except for $\left|0,\downarrow\right\rangle$ with
corresponding energy
$\frac{\varepsilon_{0,\downarrow}}{\hbar\omega_{c}}=\frac{1}{2}(1-\tilde{\Delta})$.
As a consequence, a linear combination of $\left|n,\uparrow\right\rangle$ and
$\left|n+1,\downarrow\right\rangle$ is also an eigenstate of our Hamiltonian
Eq. (31). This motivates us to rewrite the total Hamiltonian as a direct sum
of $2\times 2$ block Hamiltonians in the basis
$\left\\{\left|n,\uparrow\right\rangle,\thinspace\left|n+1,\downarrow\right\rangle\right\\}$(${\cal\tilde{H}}_{\left|n,\uparrow\right\rangle;\left|n+1,\downarrow\right\rangle}$),
in addition to the non-degenerate decoupled Hamiltonian
(${\cal\tilde{H}}_{\left|0,\downarrow\right\rangle}$), namely

${\cal\tilde{H}}={\cal\tilde{H}}_{\left|0,\downarrow\right\rangle}\oplus\bigoplus_{n=0}^{\infty}{\cal\tilde{H}}_{\left|n,\uparrow\right\rangle;\left|n+1,\downarrow\right\rangle},$
(32)

with
${\cal\tilde{H}}_{\left|0,\downarrow\right\rangle}=\varepsilon_{0,\downarrow}$
and

${\cal\tilde{H}}_{\left|n,\uparrow\right\rangle;\left|n+1,\downarrow\right\rangle}=\hbar\omega_{c}\left(\begin{array}[]{cc}n+\frac{1}{2}+\frac{\tilde{\Delta}}{2}&2\alpha_{B}\sqrt{n+1}\\\
2\alpha_{B}\sqrt{n+1}&n+1+\frac{1}{2}-\frac{\tilde{\Delta}}{2}\end{array}\right).$
(33)

The diagonalization of the Hamiltonian Eq. (33) yields energies

$\displaystyle\frac{\varepsilon_{n,s}}{\hbar\omega_{c}}=$
$\displaystyle\left(n+\frac{1}{2}+\frac{s}{2}\right)$ (34)
$\displaystyle-\frac{s}{2}\frac{1-\tilde{\Delta}}{|1-\tilde{\Delta}|}\sqrt{\left(1-\tilde{\Delta}\right)^{2}+16\alpha_{B}^{2}\left(n+\frac{1}{2}+\frac{s}{2}\right)},$

with $s=\pm$ and $n\in\mathbb{N}_{0}$, which already incorporates the energy
of the decoupled state $\left|0,\downarrow\right\rangle$,
$\varepsilon_{0,-}\equiv\varepsilon_{0,\downarrow}$
$(\varepsilon_{0,+}\equiv\varepsilon_{0,\downarrow})$ if
${1-\tilde{\Delta}>0}$ ${(1-\tilde{\Delta}<0)}$.

Figure 5: Landau levels $n=4$, 5 [Eq. (34)] as a function of $1/B$ for a 2DEG
with non-zero Zeeman and Rashba interactions but no Dresselhaus coupling
($\beta=0$). The dotted line denotes $\varepsilon_{F}/\hbar\omega_{c}$. The
parameters here are $\alpha=10$ meV nm, $m^{*}=0.019m_{o}$, $g^{*}=-34$ and
$n_{2D}=3.3\times 10^{-3}$ nm-2 for InSb-based wellsGilbertsson  _et al._
(2008); Akabori  _et al._ (2008). The dashed lines show the corresponding
levels for $\alpha=0$.

These LLs are plotted in Fig. 5 as a function of $1/B$ for parameters
$\alpha=10$ meV nm, $m^{*}=0.019m_{o}$ and $g^{*}=-34$ Gilbertsson  _et al._
(2008); Akabori  _et al._ (2008). Due to the spin-orbit coupling, the energy
levels $\frac{\varepsilon_{l,s}}{\hbar\omega_{c}}$ are no longer equidistant,
and their separation changes as function of $1/B$. On this scale, the energy
dispersion appears linear in $1/B$. In fact, for $\tilde{\Delta}<0$
($\tilde{\Delta}>0$) the spin-splitting is enhanced (suppresses) relative to
the case with $\alpha=0$ (See Fig. 3). This can be seen through the expansion
of the term $(1-\tilde{\Delta})^{2}$ within the square root of Eq. (34),
yielding $-2\tilde{\Delta}$, which enhances the Zeeman splitting in the
presence of Rashba SO coupling.Akabori  _et al._ (2008)

Table 1: Definitions of the Zeeman and SO-related quantities used is this work.  Zeeman $\left(g^{*}\right)$ | $\omega_{c}=\frac{eB}{m^{*}}$ | $\ell_{c}=\sqrt{\frac{\hbar}{eB}}$ | $\Delta={g^{*}\mu_{B}B}$ | $\tilde{\Delta}=\frac{\Delta}{\hbar\omega_{c}}=\frac{g^{*}m^{*}}{2m_{0}}$
---|---|---|---|---
Rashba $\left(\alpha\right)$ | $\varepsilon_{R}=\frac{\alpha^{2}m^{*}}{2\hbar^{2}}=\frac{\hbar^{2}k_{R}^{2}}{2m^{*}}$ | $\alpha_{B}=\frac{\alpha}{\sqrt{2}\hbar\omega_{c}\ell_{c}}$ | $\frac{\varepsilon_{R}}{\hbar\omega_{c}}=\alpha_{B}^{2}$ |
Dresselhaus $\left(\beta\right)$ | $\varepsilon_{D}=\frac{\beta^{2}m^{*}}{2\hbar^{2}}=\frac{\hbar^{2}k_{D}^{2}}{2m^{*}}$ | $\beta_{B}=\frac{\beta}{\sqrt{2}\hbar\omega_{c}\ell_{c}}$ | $\frac{\varepsilon_{D}}{\hbar\omega_{c}}=\beta_{B}^{2}$ |
SO parameters | $\gamma=\alpha_{B}+\beta_{B}$ | $\delta=\alpha_{B}-\beta_{B}$ | $\Omega=\frac{2\varepsilon_{R}/\hbar\omega_{c}}{1-\tilde{\Delta}}+\frac{2\varepsilon_{D}/\hbar\omega_{c}}{1+\tilde{\Delta}}$ | $\Lambda=\frac{2\varepsilon_{R}/\hbar\omega_{c}}{1-\tilde{\Delta}}-\frac{2\varepsilon_{D}/\hbar\omega_{c}}{1+\tilde{\Delta}}$

Accordingly, for this case we obtain

$\displaystyle{\cal F}_{+}(\varepsilon,B)$
$\displaystyle=\frac{\varepsilon}{\hbar\omega_{c}}-\frac{1}{2}+2\alpha_{B}^{2},$
(35) $\displaystyle{\cal F}_{-}(\varepsilon,B)$
$\displaystyle=-\frac{1}{2}+\frac{1}{2}\frac{1-\tilde{\Delta}}{|1-\tilde{\Delta}|}\sqrt{(1-\tilde{\Delta})^{2}+16\alpha_{B}^{2}\left(\alpha_{B}^{2}+\frac{\varepsilon}{\hbar\omega_{c}}\right)}.$
(36)

Differently from the results in the previous section, here both
$\mathcal{F}_{\pm}$ functions depend on the magnetic field. As a consequence,
we will have more complex oscillations in $\rho_{xx}(B)$ as compared to the
case without Rashba coupling [Fig. 4].

In Fig. 6, we plot the total differential magneto-resistivity
$\delta\rho_{xx}(B)$, and the independent contributions from harmonics $l=1,2$
and $l=3$. Here we use $\alpha=10$ meV nm, $m^{*}=0.019m_{o}$, $g^{*}=-34$ and
$n_{2D}=3.3\times 10^{-3}$ nm-2 Gilbertsson  _et al._ (2008); Akabori  _et
al._ (2008). Similarly to the case with $\alpha=0$ (dashed line in Fig. 6),
here we also see oscillations for the $l=1,2,3$ harmonics with respective
frequencies $f^{\rm{SdH}},2f^{\rm{SdH}}$ and $3f^{\rm{SdH}}$. However, for
$l=1$ we observe beating, which can be expected as both ${\cal
F}_{-}(\varepsilon,B)$ and ${\cal F}_{+}(\varepsilon,B)$ frequencies now
depend on $1/B$. More specifically, this beating appears here because in the
magnetic range considered we have $2\pi l\mathcal{F}_{-}(B)=\frac{\pi}{2}$,
which leads to a node in $\delta\rho_{xx}$ as
$\delta\rho_{xx}\propto\cos\left[2\pi l\mathcal{F}_{-}(B)\right]$. Note that
this only occurs for $l=1$, since for higher harmonics this condition is not
satisfied. Due to the larger amplitude of the harmonic $l=1$, this beating is
also seen in the total magneto-resistivity.

Figure 6: Magnetoresistivity deviation $\delta\rho_{xx}(B)$ as a function of
$1/B$ for a 2DEG with Zeeman and Rashba interactions and no Dresselhaus
coupling ($\beta=0$). The lowest curve corresponds to $\delta\rho_{xx}(B)$ and
the curves labelled by $l$ are the individual frequency components (i.e.,
harmonics) in Eq. (30). The solid lines are obtained with $g^{*}=-34$,
$\alpha=10$ meV nm, ${m^{*}=0.019m_{o}}$ and $n_{2D}=3.3\times 10^{-3}$ nm-2
Gilbertsson  _et al._ (2008); Akabori  _et al._ (2008); the dotted lines show
the corresponding $\alpha=0$ case.

###  V.3 Landau Levels with Zeeman and Dresselhaus interaction

In the case of Zeeman with pure Dresselhaus, i.e., $\tilde{\Delta}\neq 0$,
$\alpha=0$ and $\beta\neq 0$, the Hamiltonian Eq. (9) in the spin basis is
given by

$\frac{{\cal\tilde{H}}}{\hbar\omega_{c}}=\left(\begin{array}[]{cc}a^{\dagger}a+\frac{1}{2}-\frac{\tilde{\Delta}}{2}&2\beta_{B}a^{\dagger}\\\
2\beta_{B}a&a^{\dagger}a+\frac{1}{2}+\frac{\tilde{\Delta}}{2}\end{array}\right).$
(37)

Differently from the case of pure Rashba, here the operator ${\cal
N}_{-}=a^{\dagger}a-\sigma_{z}/2$ commutes with the Hamiltonian above. For
this case we have ${{\cal
N}_{-}\left|n,\downarrow\right\rangle=\left(n+1/2\right)\left|n,\downarrow\right\rangle}$
and ${{\cal
N}_{-}\left|n+1,\uparrow\right\rangle=\left(n+1/2\right)\left|n+1,\uparrow\right\rangle}$,
i.e., for $n\in\mathbb{N}$, $\left|n,\downarrow\right\rangle$ and
$\left|n+1,\uparrow\right\rangle$ are degenerate with respect to the operator
${\cal N}_{-}$, except for the state $\left|0,\uparrow\right\rangle$ with
corresponding energy
$\frac{\varepsilon_{0,\uparrow}}{\hbar\omega_{c}}=\frac{1}{2}(1+\tilde{\Delta})$.
As a consequence, a linear combination of $\left|n,\downarrow\right\rangle$
and $\left|n+1,\uparrow\right\rangle$ is also an eigenstate of our
Hamiltonian. Therefore, differently from the previous case here the
Hamiltonian reads,

${\cal\tilde{H}}={\cal\tilde{H}}_{\left|0,\uparrow\right\rangle}\oplus\bigoplus_{n=0}^{\infty}{\cal\tilde{H}}_{\left|n,\downarrow\right\rangle;\left|n+1,\uparrow\right\rangle},$
(38)

with
${\cal\tilde{H}}_{\left|0,\uparrow\right\rangle}=\varepsilon_{0,\uparrow}$ and

Figure 7: Landau levels $n=4$, 5 [Eq. (40)] as a function of $1/B$ for a 2DEG
with Zeeman and Dresselhaus interations but no Rashba coupling ($\alpha=0$).
The dotted line denotes $\varepsilon_{F}/\hbar\omega_{c}$. The parameters here
are $\beta=10$ meV nm, $m^{*}=0.019m_{o}$, $g^{*}=-34$ and $n_{2D}=3.3\times
10^{-3}$ nm-2 Gilbertsson  _et al._ (2008); Akabori  _et al._ (2008) dashed
lines show the corresponding levels for $\beta=0$.

${\cal\tilde{H}}_{\left|n,\downarrow\right\rangle;\left|n+1,\uparrow\right\rangle}=\hbar\omega_{c}\left(\begin{array}[]{cc}n+\frac{1}{2}-\frac{\tilde{\Delta}}{2}&2\beta_{B}\sqrt{n+1}\\\
2\beta_{B}\sqrt{n+1}&n+1+\frac{1}{2}+\frac{\tilde{\Delta}}{2}\end{array}\right).$
(39)

The diagonalization of the Hamiltonian Eq. (39) yields energies

$\displaystyle\frac{\varepsilon_{n,s}}{\hbar\omega_{c}}=$
$\displaystyle\left(n+\frac{1}{2}-\frac{s}{2}\right)$ (40)
$\displaystyle+\frac{s}{2}\frac{1+\tilde{\Delta}}{|1+\tilde{\Delta}|}\sqrt{\left(1+\tilde{\Delta}\right)^{2}+16\beta_{B}^{2}\left(n+\frac{1}{2}-\frac{s}{2}\right)},$

with $s=\pm$ and $n\in\mathbb{N}_{0}$, which already incorporates the energy
of the decoupled state $\left|0,\uparrow\right\rangle$,
$\varepsilon_{0,+}\equiv\varepsilon_{0,\uparrow}$
($\varepsilon_{0,-}\equiv\varepsilon_{0,\uparrow}$) if $1+\tilde{\Delta}>0$
($1+\tilde{\Delta}<0$). Here, it is important to notice the opposite sign of
$s$ with respect to Eq. (33). This happens because the pure Dresselhaus
Hamiltonian Eq. (39) has opposite basis ordering of the spin states as
compared to the pure Rashba Hamiltonian Eq. (33). Accordingly, the ${\cal
F}_{\pm}$ functions change slightly and read

$\displaystyle{\cal F}_{+}(\varepsilon,B)$
$\displaystyle=\frac{\varepsilon}{\hbar\omega_{c}}-\frac{1}{2}+2\beta_{B}^{2},$
(41) $\displaystyle{\cal F}_{-}(\varepsilon,B)$
$\displaystyle=\frac{1}{2}-\frac{1}{2}\frac{1+\tilde{\Delta}}{|1+\tilde{\Delta}|}\sqrt{\left(1+\tilde{\Delta}\right)^{2}+16\beta_{B}^{2}\left(\beta_{B}^{2}+\frac{\varepsilon}{\hbar\omega_{c}}\right)}.$
(42)

Due to the: i) similarity of the Dresselhaus expression Eqs. (40), (41) and
(42) to the ones arising from the pure Rashba case, Eqs. (34), (35) and (36);
ii) cosine dependence of the ${\cal F}_{\pm}$ functions within the resistivity
Eq. (22); all the results and equations in the last section also holds here by
making $\alpha_{B}\rightarrow\beta_{B}$,
$\tilde{\Delta}\rightarrow-\tilde{\Delta}$ and $s\rightarrow-s$. This can also
be seen on the level of the Hamiltonian in Eq. (2) where applying the unitary
transformation $W=e^{i\frac{\pi}{2}\sigma_{x}}e^{i\frac{\pi}{4}\sigma_{z}}$
results in

$\displaystyle W\frac{{\cal\tilde{H}}}{\hbar\omega_{c}}W^{\dagger}$
$\displaystyle=$
$\displaystyle(a^{\dagger}a+1/2)+\frac{(-\tilde{\Delta})}{2}\sigma_{z}+\beta_{B}(a^{\dagger}\sigma_{-}+a\sigma_{+}),$
(43) $\displaystyle+\alpha_{B}(a^{\dagger}\sigma_{+}+a\sigma_{-}),$

which is the expected result. This mapping from $(\alpha,\tilde{\Delta})$ to
$(\beta,-\tilde{\Delta})$ has visible consequences on the energy levels. In
Fig. 7 we plot the corresponding LLs [Eq. (40)] as a function of $1/B$ for
parameters $\beta=10$ meV nm, $m^{*}=0.019m_{o}$ and $g^{*}=-34$ Gilbertsson
_et al._ (2008); Akabori  _et al._ (2008). Due to the spin-orbit coupling, the
energy levels $\frac{\varepsilon_{l,s}}{\hbar\omega_{c}}$ are no longer
equidistant, and their separation changes as function of $1/B$. However,
differently from the pure Rashba case, now the Dresselhaus competes with the
Zeeman coupling, even leading to LL-dependent crossings. This can be seen
through the expansion of $(1+\tilde{\Delta})^{2}$ within the square root [Eq.
(40)], which will give rise to $2\tilde{\Delta}<0$, thus suppressing the spin
splitting in the presence of Dresselhaus SO coupling.

Figure 8: Magnetoresistivity deviation $\delta\rho_{xx}(B)$ as a function of
$1/B$ for a 2DEG with Zeeman and Dresselhaus interactions and no Rashba
coupling ($\alpha=0$). The lowest curve corresponds to $\delta\rho_{xx}(B)$
and the curves labelled by $l$ are the individual frequency components in Eq.
(30). The solid lines are calculated for $g^{*}=-34$, $\beta=10$ meV nm,
${m^{*}=0.019m_{o}}$ and $n_{2D}=3.3\times 10^{-3}$ nm-2 Gilbertsson  _et al._
(2008); Akabori  _et al._ (2008); the dotted line shows the corresponding
$\beta=0$ case.

In Fig. 8 we plot the total differential magneto-resistivity
$\delta\rho_{xx}(B)$, and the individual contributions from the harmonics
${l=1,2}$ and $l=3$. We use $\beta=10$ meV nm, $m^{*}=0.019m_{o}$, $g^{*}=-34$
and $n_{2D}=3.3\times 10^{-3}$ nm-2 Gilbertsson  _et al._ (2008); Akabori  _et
al._ (2008). First, similarly to the previous cases, here we can also clearly
see oscillations with frequencies $f^{\rm{SdH}},2f^{\rm{SdH}}$,
$3f^{\rm{SdH}}$. Differently from the previous case with $\alpha=10$ meV nm
and $\beta=0$, now we see no beating for the $l=1$ harmonic but find beating
for $l=2$. This happens as $2\pi l\mathcal{F}_{-}(B)=\frac{\pi}{2}$ – the
condition to observe beating – is only satisfied for $l=2$. Even though the
beating appears within the second harmonic, it is not manifested in the total
differential magneto-resistivity $\delta\rho_{xx}(B)$ for our choice of
parameters. This is due to the smaller oscillation amplitude of $l=2$ with
respect to $l=1$.

###  V.4 Beatings in the SdH oscillations with nonzero Zeeman and in the
presence of either Rashba or Dresselhaus: a unified description

In this section we will discuss more thoroughly the conditions for the
appearance of beatings. The two functions $\mathcal{F}_{+}$ and
$\mathcal{F}_{-}$, Eq. (13), determine the fast and slow component,
respectively, of the SdH oscillations. To highlight this point and its
connection to the power spectrum in Eq. (26), we start by rewriting Eqs. (35)-
(36), and Eqs. (42)- (42) in a unified way

$\displaystyle{\cal F}_{+}(\varepsilon,B)$
$\displaystyle=f_{{R(D)}}^{\mathrm{SdH}}\frac{1}{B}-\frac{1}{2},$ (44)
$\displaystyle{\cal F}_{-}(\varepsilon,B)$
$\displaystyle={\mp}\frac{1}{2}{\pm}{\frac{1}{2}\frac{1\mp\tilde{\Delta}}{|1\mp\tilde{\Delta}|}}\sqrt{(1{\mp}\tilde{\Delta})^{2}+4\left(f_{R(D)}\frac{1}{B}\right)^{2}},$
(45)

where we have introduced the magneto-oscillation frequencies

$\displaystyle f_{{R(D)}}^{\mathrm{SdH}}$
$\displaystyle=\frac{h}{2e}\left(n_{2D}+\frac{k_{R(D)}^{2}}{\pi}\right),$ (46)
$\displaystyle f_{R(D)}$
$\displaystyle=\frac{h}{2e}\sqrt{\frac{2k_{R(D)}^{2}}{\pi}}\sqrt{n_{2D}+\frac{k_{R(D)}^{2}}{2\pi}},$
(47)

where the $R$ $(D)$ index refers to either pure Rasbha (Dresselhaus) case,
with $k_{R}=\frac{m\alpha}{\hbar^{2}}$ ($k_{D}=\frac{m\beta}{\hbar^{2}}$).
Here, the upper (lower) sign refers to the Rashba (Dresselhaus) case.

Figure 9: Frequency response $\sqrt{I(f)}$ for $\alpha=10.0$ meV nm and
$\beta=0.0$ [red curve], and $\alpha=0.0$ and $\beta=10.0$ meV nm [blue
curve]. Other parameters are $m^{*}=0.019m_{o}$, $g^{*}=-34$ and
$n_{2D}=3.3\times 10^{-3}$ nm-2 Gilbertsson  _et al._ (2008); Akabori  _et
al._ (2008). The solid black shows corresponds to no spin-orbit coupling
($\alpha=\beta=0$) and black dashed corresponds to $\alpha=\beta=g=0$.

In the case where $n_{2D}\gg{k_{R(D)}^{2}}/{2\pi}$, and ${f_{R(D)}}/{B}\gg 1$,
the beating frequency takes the standard form
$f_{R(D)}=\frac{h}{2e}\sqrt{{2k_{R(D)}^{2}n_{2D}}/{\pi}}$, in which case
$\tilde{\Delta}$ becomes irrelevant for the magnitude of the beating
frequencyEngels  _et al._ (1997).

The frequency $f_{R(D)}^{\rm{SdH}}$ [Eq. (46)] is the main SdH frequency of
the magneto-resistance oscillations, usually extracted from experiments to
infer the 2D electronic density $n_{2D}$. On the other hand, the frequency
$f_{R(D)}$ [Eq. (47)] is the one allowing for possible beatings in the
magneto-oscillation. As previously discussed in the last two sections, the
presence of beating happens when $2\pi l\mathcal{F}_{-}(B)=\frac{\pi}{2}$ is
satisfied, which depends on the value of both $f_{R(D)}$ and $\tilde{\Delta}$.

The presence or absence of beatings can also be visualized through the power
spectrum defined by Eq. (26). From interference of waves, we know that the
presence of beatings correspond to sum of cosines waves with slightly
different frequencies. Accordingly, the power spectrum for this case would
show two peaks located at slightly different frequencies. In Fig. 9 we plot
$\sqrt{I(f)}$ for $m^{*}=0.019m_{o}$ and $n_{2D}=3.3\times 10^{-3}$ nm-2,
using different spin-orbit parameters and $g$-factor values. For all different
sets of parameters, we always have the presence of two main peaks located at
both $1/B\approx 6.8$ T-1 and $1/B\approx 13.6$ T-1. These correspond to the
main SdH frequencies for the first and second harmonics, $f_{R(D)}^{\rm{SdH}}$
and $2f_{R(D)}^{\rm{SdH}}$, respectively. In the absence of both Rashba,
Dresselhaus and $g$-factor (dashed yellow curve), we observe no beating in the
$\delta\rho_{xx}$ (Fig. 4).

On the other hand, for the case of pure Rashba $\alpha=10$ meV nm with $g=-34$
(solid red curve), the presence of the beating in Fig. 6 is made clear by the
splitting of the peak of the power spectrum around $f=f_{R}^{\rm{SdH}}$ in
Fig. 9. Interestingly, for $\alpha=10$ meV nm with $g=0$ (dashed red curve),
the splitting of the peak is not seen anymore, thus highlighting the important
role of the Zeeman on the visualization of beatings. For the pure Dresselhaus
case with $\beta=10$ meV nm and $g=-34$ (solid blue line), we do not see a
peak splitting at the $f=f_{D}^{\rm{SdH}}$ but rather at
$f=2f_{D}^{\rm{SdH}}$, which is consistent with the presence of the beating
seen on the second harmonic in Fig. 8. Similarly to the pure Rashba case, for
$\beta=10$ meV nm with $g=0$ (dashed blue line), the splitting of the peak is
not seen anymore, corroborating again the role of the Zeeman term on the
presence of beatings.

The apparent “asymmetry” in having peak-splitting for Rashba spin-orbit
coupling but not for Dresselhaus (even when they have same SO strength) can be
understood from the behavior of the $l\mathcal{F}_{-}$-function vs $1/B$,
shown in Fig. 10. As already discussed previously in Secs. V.2 and V.3, the
condition for beating happens when $\cos(2l\pi\mathcal{F}_{-})=0$ or
equivalently, $l\mathcal{F}_{-}=\pm 1/4$ ($\pm 1/4$ plotted as gray lines). In
the case of Rashba (green lines) one has $(1-\tilde{\Delta})>1$, and the
condition for a beating node, $\cos(2l\pi\mathcal{F}_{-})=0$, is reached in
the interval of $1/B$ for $l=1$ (solid purple) (gray circles). In the
Dresselhaus case, $(1+\tilde{\Delta})<1$, such that $l\mathcal{F}_{-}$ for
$l=1$ only crosses $-1/4$ for large values of $1/B$, where the amplitude of
the SdH has already been suppressed. Conversely, $l\mathcal{F}_{-}$ crosses
1/4 for $l=2$ at smaller values of $1/B$, thus guaranteeing the presence of a
beating within the magnetic field range, as shown in Fig. 8.

Figure 10: Plot of $l\mathcal{F}_{-}(B)$ vs $1/B$ for $l=1$ and $l=2$ using
$\alpha=10$ meVnm with $\beta=0$ (purple lines), $\alpha=\beta=0$ (cyan
lines), and $\beta=10$ meVnm with $\alpha=0$ (green lines). The solid gray
lines indicate $\pm 1/4$ and the gray circles indicate where beating nodes
occur. For all curves, we use $m^{*}=0.019$$m_{o}$, $g^{*}=-34$ and
$n_{2D}=3.3\times 10^{-3}$ nm-2, parameters for InSb-based 2DEGsGilbertsson
_et al._ (2008); Akabori  _et al._ (2008).

###  V.5 Landau Levels with simultaneous Zeeman, Rashba and Dresselhaus
interactions: Analytical results

As mentioned earlier, to the best of our knowledge, there are no general
_exact_ analytical results for the energies and SdH oscillations corresponding
to the case with simultaneous and arbitrary Zeeman, Rashba and Dresselhaus
couplings. Therefore, in this section we will outline how to derive an
effective _approximate_ solution that can be used to shed light on
magnetotransport results for materials, e.g. GaAs or InAs, in which all the
three couplings are present. For convenience, we define the sum and difference
of the spin-orbit couplings

$\displaystyle\gamma$ $\displaystyle=$ $\displaystyle\alpha_{B}+\beta_{B},$
(48) $\displaystyle\delta$ $\displaystyle=$
$\displaystyle\alpha_{B}-\beta_{B},$ (49)

[see definitions of $\alpha_{B}$ and $\beta_{B}$ following Eq. (9)] which
allows us to rewrite Eq. (9) as

$\displaystyle\frac{{\cal\tilde{H}}}{\hbar\omega_{c}}$ $\displaystyle=$
$\displaystyle
a^{\dagger}a+\frac{1}{2}+\frac{\tilde{\Delta}}{2}\sigma_{z}+\frac{\gamma+\delta}{2}\left(a^{\dagger}\sigma_{-}+a\sigma_{+}\right)$
(50)
$\displaystyle+\frac{\gamma-\delta}{2}\left(a^{\dagger}\sigma_{+}+a\sigma_{-}\right).$

Note that both the pure Rashba and pure Dresselhaus cases are recovered from
the equation above for $\gamma=\delta$ and $\gamma=-\delta$, respectively.
Next, we define the Hamiltonian for $\gamma=\delta$ and $\gamma=-\delta$

$\displaystyle\frac{{\mathcal{\tilde{H}}}_{\mathrm{\pm}}}{\hbar\omega_{c}}$
$\displaystyle=$ $\displaystyle
a^{\dagger}a+\frac{1}{2}+\frac{\tilde{\Delta}}{2}\sigma_{z}\pm\delta(a^{\dagger}\sigma_{\mp}+a\sigma_{\pm}),$
(51)

which describes the pure Rashba ($+$) and pure Dresselhaus ($-$) cases in the
presence of the Zeeman coupling. As we already discussed in the previous
sections, by defining the operator
$\mathcal{N}_{\pm}=a^{\dagger}a\pm\frac{1}{2}\sigma_{z}$, we obtain
$[{{\cal\tilde{H}}}_{\pm},\mathcal{N}_{\pm}]=0$, so the eigenstates of
${{\cal\tilde{H}}}_{\pm}$ are also eigenstates of $\mathcal{N}_{\pm}$. The
eigenstates of $\mathcal{N}_{+}$ ${({\mathcal{N}}_{-})}$ are then constructed
from the pair $\\{|n,\uparrow\rangle,|n+1,\downarrow\rangle\\}$
($\\{|n,\downarrow\rangle,|n+1,{\uparrow}\rangle\\}$). The above statement is
true except for the decoupled eigenstates $|0,\uparrow\rangle$
($|0,\downarrow\rangle$) with corresponding eigenenergy
$\hbar\omega_{c}(1-\tilde{\Delta})/2$ $[\hbar\omega_{c}(1+\tilde{\Delta})/2]$.
The diagonalization of each two state subspace results in

$\displaystyle\frac{\varepsilon_{n,s}}{\hbar\omega_{c}}$ $\displaystyle=$
$\displaystyle\left(n+\frac{1}{2}+\frac{\delta}{|\delta|}\frac{s}{2}\right)-\frac{\delta}{|\delta|}\frac{s}{2}\left(1-\frac{\delta}{|\delta|}\tilde{\Delta}\right)$
(52) $\displaystyle\times$
$\displaystyle\sqrt{1+\frac{16\delta^{2}}{(1-\frac{\delta}{|\delta|}\tilde{\Delta})^{2}}{\left(n+\frac{1}{2}+\frac{\delta}{|\delta|}\frac{s}{2}\right)}},$

with $s=+$ $(-)$ and $n\in\mathbb{N}_{0}$. Note that this form is valid for
both pure Rashba ($\delta=\gamma>0$) and Dresselhaus ($\delta=-\gamma<0$),
Eqs. (34) and (40), respectively, thus also including the corresponding
decoupled state with the lowest eigenvalues of ${{\cal N}}_{\pm}$. Note that
to recover the pure Zeeman case with no Rashba and Dresselhaus, we should take
$\delta\rightarrow 0$ with $\delta/|\delta|\rightarrow 1$.

When both Rashba and Dresselhaus are present, we can use second order
perturbation theory with respect to $\delta,\gamma\ll 1$ (See Appendix E), to
obtain the approximate eigenvalues of the Hamiltonian in Eq. (50), namely

$\displaystyle\frac{\varepsilon_{n,s}}{\hbar\omega_{c}}$ $\displaystyle=$
$\displaystyle
n+\nicefrac{{1}}{{2}}+s\frac{\tilde{\Delta}}{2}-2s\Lambda(l+\nicefrac{{1}}{{2}})-\Omega$
(53)

where the quantities $\Lambda$ and $\Omega$ are defined as

$\displaystyle\Lambda$ $\displaystyle=$
$\displaystyle\frac{(\gamma^{2}+\delta^{2})\tilde{\Delta}+2\gamma\delta}{(1-\tilde{\Delta}^{2})}=\frac{{2}\frac{\varepsilon_{R}}{\hbar\omega_{c}}}{(1-\tilde{\Delta})}-\frac{{2}\frac{\varepsilon_{D}}{\hbar\omega_{c}}}{(1+\tilde{\Delta})},$
(54) $\displaystyle\Omega$ $\displaystyle=$
$\displaystyle\frac{(\gamma^{2}+\delta^{2})+2\gamma\delta\tilde{\Delta}}{(1-\tilde{\Delta}^{2})}=\frac{{2}\frac{\varepsilon_{R}}{\hbar\omega_{c}}}{(1-\tilde{\Delta})}+{{\color[rgb]{0,0,0}\frac{{2}\frac{\varepsilon_{D}}{\hbar\omega_{c}}}{(1+\tilde{\Delta})},}}$
(55)

where we have introduced ${\varepsilon_{R}}/{\hbar\omega_{c}}=\alpha_{B}^{2}$
and $\varepsilon_{D}/\hbar\omega_{c}=\beta_{B}^{2}$.

Our goal now is to rewrite Eq. (53) in a form that recovers the already
obtained exact results for pure Rashba and pure Dresselhaus cases. First, we
write $\Lambda=\frac{\Lambda}{|\Lambda|}|\Lambda|$ since $\Lambda$ changes
sign depending on the relative strengths of $\alpha$ and $\beta$, similarly to
the sign of $\delta$ that enters into Eq. (52). By adding and subtracting a
term $\frac{s}{2}\,\frac{\Lambda}{|\Lambda|}$ in Eq. (53) and after some
straightforward algebra we obtain

$\displaystyle\frac{\varepsilon_{n,s}}{\hbar\omega_{c}}$ $\displaystyle=$
$\displaystyle\left(n+\frac{1}{2}+\frac{\Lambda}{|\Lambda|}\frac{s}{2}\right)-\frac{\Lambda}{|\Lambda|}\frac{s}{2}{\left(1-\frac{\Lambda}{|\Lambda|}\tilde{\Delta}\right)}$
$\displaystyle\times$
$\displaystyle\left\\{1+\frac{{4}}{1-\frac{\Lambda}{|\Lambda|}\tilde{\Delta}}{\left[|\Lambda|\left(n+\frac{1}{2}\right)+\,\Omega\frac{\Lambda}{|\Lambda|}\frac{s}{2}\right]}{}\right\\}.$

In the case of pure Rashba we have
$\Lambda=\Omega=\frac{\delta^{2}}{1-\tilde{\Delta}}>0$ while for pure
Dresselhaus $\Lambda=-\Omega=-\frac{\delta^{2}}{1+\tilde{\Delta}}<0$; these
neatly reduce to the exact results when using second order Taylor expansion of
Eq. (52). Note that Eq. (LABEL:eq:pertEnSgn) also reproduces the exact result
for when $\alpha=\beta$ and $g^{*}=0$ Tarasenko and Averkiev (2002),
represented here by $\Lambda\rightarrow 0$ with $\Lambda/|\Lambda|\rightarrow
1$, $\tilde{\Delta}=0$, and $\Omega=2\varepsilon_{D/R}/\hbar\omega_{c}$ . The
mathematical expression of Eqs. (34) and (40) motivate us to rewrite Eq.
(LABEL:eq:pertEnSgn) as

$\displaystyle\frac{\varepsilon_{n,s}}{\hbar\omega_{c}}=\left(n+\frac{1}{2}+\frac{\Lambda}{|\Lambda|}\frac{s}{2}\right)-\frac{\Lambda}{|\Lambda|}\frac{s}{2}\frac{1-\frac{\Lambda}{|\Lambda|}\tilde{\Delta}}{|1-\frac{\Lambda}{|\Lambda|}\tilde{\Delta}|}\times\sqrt{\left(1-\frac{\Lambda}{|\Lambda|}\tilde{\Delta}\right)^{2}+8\left(1-\frac{\Lambda}{|\Lambda|}\tilde{\Delta}\right)\left[|\Lambda|\left(n+\frac{1}{2}\right)+\,\Omega\frac{\Lambda}{|\Lambda|}\frac{s}{2}\right]},$
(57)

where we have used $1+\frac{x}{2}\approx\sqrt{1+x}$ 222Here we emphasize that
depending on the values of $\Lambda$ and $\Omega$, the argument of square root
of Eq. (57) becomes negative, thus yielding imaginary energies. This happens
when $x\gg 1$, which violates the assumption of writing Eq. (57). It is
important to note that although $|\Lambda|\ll 1$, $\Lambda$ enters the square
root multiplied by $n$, the Landau level index. This means that for high
enough $n$, the product $|\Lambda|n$ is not necessarily a small quantity.
Accordingly, although the equation above becomes exact for either pure Rashba
or Dresselhaus case, for $\alpha,\beta\neq 0$ Eq. (57) is only valid when
$|\Lambda|n\lesssim 1$, besides $\alpha_{B},\beta_{B},\delta,\gamma\ll 1$
already assumed in Appendix E to obtain Eq. (53).

We reiterate that Eq. (57) satisfies the exact results for (i) the Zeeman-only
case [Eq. (27)], (ii) the pure Rashba plus nonzero $g^{*}$ [Eq. (34)] and
(iii) the pure Dresselhaus plus nonzero $g^{*}$ [Eq. (40)]. The case
$\alpha=\beta$ with $g^{*}=0$, for which there is also an exact
solutionTarasenko and Averkiev (2002), is satisfied to leading order using
$\sqrt{1+x}\approx 1+x/2$ for with $x=8\Omega(s/2)/(1-\tilde{\Delta})\ll 1$.
That is, as mentioned in the previous paragraph, the approximate solution
given by Eq. (LABEL:eq:pertEnSgn) reproduces the exact solution for
$\alpha=\beta$ with $g^{*}=0$Tarasenko and Averkiev (2002).

As in the case of pure Zeeman, Rashba or Dresselhaus, we can now calculate the
$F$-function from Eq. (57). The corresponding results are presented in
Appendix F, and by neglecting SO contributions higher or equal than second
order in the spin-orbit parameters $\Lambda$ and $\Omega$ (or fourth order in
$\gamma$ and $\delta$), we obtain

$\displaystyle{\cal F}_{+}$
$\displaystyle=\frac{\varepsilon}{\hbar\omega_{c}}-\frac{1}{2}+\Omega-\Lambda\tilde{\Delta},$
(58) $\displaystyle{\cal F}_{-}$
$\displaystyle=-\frac{1}{2}\frac{\Lambda}{|\Lambda|}+\frac{1}{2}\frac{\Lambda}{|\Lambda|}\frac{1-\frac{\Lambda}{\left|\Lambda\right|}\tilde{\Delta}}{\left|1-\frac{\Lambda}{\left|\Lambda\right|}\tilde{\Delta}\right|}\sqrt{\left(1-\frac{\Lambda}{\left|\Lambda\right|}\tilde{\Delta}\right)^{2}+8\left|\Lambda\right|\left(1-\frac{\Lambda}{\left|\Lambda\right|}\tilde{\Delta}\right)\left[\frac{\varepsilon}{\hbar\omega_{c}}+\frac{1}{2}\left|\Lambda\right|\left(1-\frac{\Lambda}{\left|\Lambda\right|}\tilde{\Delta}\right)\right]}.$
(59)

It is easy to see that these equations recover all the previous results: pure
Zeeman [Eq. (29)], Zeeman with pure Rashba [Eqs. (35) and (36)], and Zeeman
with pure Dresselhaus [Eqs. (41) and (42)]. Additionally, in the case of
$\Lambda\approx 0$, ${\cal
F}_{-}{{\color[rgb]{0,0,0}\approx-\tilde{\Delta}/2}}$ , which reduces to the
pure Zeeman case. Accordingly, here ${\cal F}_{-}$ becomes independent of $B$
(for $B\lesssim 1$ T), and therefore, we expect the absence of beatings in the
magneto-resistivity, previously seen for both pure Rashba and pure Dresselhaus
cases.

###  V.6 Generalized SdH magneto-resistivity for arbitrary $\alpha$, $\beta$
and $g^{*}$ : new prediction for the absence of beatings.

Using the Eqs. (58) and (59) in Eq. (24), we can derive the magnetoresistivity
$\delta\rho_{xx}(B)$ for the case with arbitrary Rashba and Dresselhaus
couplings and simulta- neous nonzero Zeeman field,

$\displaystyle\delta\rho_{xx}(B)=$ $\displaystyle
2\sum_{l=1}^{\infty}e^{-l\pi\frac{\hbar/\tau_{q}}{\hbar\omega_{c}}}\frac{2\pi^{2}lk_{B}T/\hbar\omega_{c}}{\sinh(2\pi^{2}lk_{B}T/\hbar\omega_{c})}$
$\displaystyle\times\cos\left[2\pi
l\left(\frac{\varepsilon_{F}}{\hbar\omega_{c}}+\frac{2\varepsilon_{R}}{\hbar\omega_{c}}+\frac{2\varepsilon_{D}}{\hbar\omega_{c}}\right)\right]\cos\left\\{\pi
l\sqrt{\left(1-\frac{\Lambda}{\left|\Lambda\right|}\tilde{\Delta}\right)^{2}+16\lambda_{B}^{2}\left(\lambda_{B}^{2}+\frac{\varepsilon_{F}}{\hbar\omega_{c}}\right)}\right\\},$
(60)

with
$\lambda_{B}^{2}=\frac{\left|\Lambda\right|}{2}\left(1-\frac{\Lambda}{\left|\Lambda\right|}\tilde{\Delta}\right)$.
From Eq. (60), we can derive the condition for the absence of beatings for any
$l$ by finding the condition for the second cosine being independent of $1/B$
; this implies $|\Lambda|=0$, which leads to

$\frac{\alpha}{\beta}=\sqrt{\frac{1-\tilde{\Delta}}{1+\tilde{\Delta}}},$ (61)

thus yielding Eq. (1) presented in the introduction. For $\tilde{\Delta}\ll
1$, the above condition is reduced to $\alpha\approx\beta$, corresponding to
the situation where the total SO $k$-dependent effective field becomes
unidirectional Schliemann  _et al._ (2003a); Bernevig  _et al._ (2006a); Fu
_et al._ (2016).

Note that the above condition does not correspond to any fundamental symmetry,
since there is no new conserved quantity in our Hamiltonian with both non-zero
Zeeman ($g^{*}\neq 0$) and Rashba-Dresselhaus couplings. We reiterate that Eq.
(61) is entirely distinct from the persistent-spin-helix condition
$\alpha=\beta$. As shown in Fig. 1(d), the case $\alpha=\beta$ and $g^{*}\neq
0$ does not show peak splitting in the first harmonic but ehxibits beating (or
peak splitting) in the second harmonic. Only when $g^{*}=0$ (no Zeeman) and
$\alpha=\beta$ there are peak splittings absent altogether Tarasenko and
Averkiev (2002); Averkiev  _et al._ (2005).

###  V.7 Beatings for both $\alpha$ and $\beta$ non-zero

In the previous sections, we studied the effect of the Zeeman interaction on
the frequency splitting of the power spectrum peaks, which represents the
beatings in the SdH oscillations. Here we study the interplay of both the
Dresselhaus and Rashba interactions on the beatings of the SdH oscillations.

Similarly to what we did leading up to Eq. (47), we can obtain the effective
beating frequency from the $\mathcal{F}_{-}$–function in Eq. (59) which
results in

$\displaystyle f_{R+D}$
$\displaystyle=\frac{h}{2e}\sqrt{\left|\frac{2k_{R+D}^{2}}{\pi}\left(n_{\mathrm{2D}}+\frac{k_{R+D}^{2}}{2\pi}\right)\right|},$
(62)

where the effective SO momentum is

$k_{R+D}=\frac{m^{*}}{\hbar^{2}}\sqrt{\left(1-\frac{\Lambda}{|\Lambda|}\tilde{\Delta}\right)\left(\frac{\alpha^{2}}{1-\tilde{\Delta}}-\frac{\beta^{2}}{1+\tilde{\Delta}}\right)}.$
(63)

Figure 11: Normalized power spectrum $\sqrt{I(f)}$ for a fixed $\alpha=7.0$
meV nm, and for $\beta=0.0$ to $10$ meV nm, from bottom to top. [red curves],
along with full numerical results [black dashed]. The curve corresponding to
$k_{RD}=0$ with $\beta=5.0$ meV nm, is shown [blue curve]. The gray circles
indicate the frequency splitting in Eq. (62). Other parameters are
$m^{*}=0.019m_{o}$, $g^{*}=-34$, $n_{2D}=3.3\times 10^{-3}$ nm-2, and
$\hbar\tau_{q}^{-1}=1.75$ meV, for InSb-based 2DEGsGilbertsson  _et al._
(2008); Akabori  _et al._ (2008) Figure 12: Normalized power spectrum
$\sqrt{I(f)}$ for a fixed $\beta=5.0$ meV nm, and for $\alpha=0.0$ to $9.0$
meV nm, from bottom to top. [red curves] along with full numerical results
[black dashed]. The curve corresponding to $k_{RD}=0$ with $\alpha=7.0$ meV nm
[blue curve]. The gray circles indicate the frequency splitting in Eq. (62).
Other parameters are $m^{*}=0.019m_{o}$, $g^{*}=-34$, $n_{2D}=3.3\times
10^{-3}$ nm-2, and $\hbar\tau_{q}^{-1}=1.75$ meV, for InSb-based
2DEGsGilbertsson  _et al._ (2008); Akabori  _et al._ (2008)

We start with the pure Rashba case plus Zeeman, $\alpha=7.0$ meV nm and
$g^{*}=-34$. The corresponding power spectrum yields the bottom curve in Fig.
11, similar to the one plotted in Fig. 9. This curve shows two main peaks
representing the first two harmonics, and the presence of a split main peak.
We assume a Lorentzian broadening $\hbar\tau_{q}^{-1}=1.75$ meV. When the
Dresselhaus coupling $\beta$ increases, we see the splitting of the main peak
reduces until it vanishes for $\beta=5.0$ meV nm (the frequency splitting from
Eq. (62) is indicated by the gray circles). The absence of beating is indeed
expected as predicted by the condition
$\beta=\alpha\sqrt{\frac{1+\tilde{\Delta}}{1-\tilde{\Delta}}}=5.0$ meV in Eq.
(61). For larger $\beta$, we see that the splitting of the main peak remains
neglible. However, in the second harmonic a clear splitting opens up. The
condition for having no peak splitting at any harmonics is indeed the
condition in Eq. (61), where the effects of the SO couplings basically
disappear [there are still small SO terms $\varepsilon_{R}$ ,
$\varepsilon_{D}$ in Eq. (60)]. The power spectrum using full numerical
calculations are also shown [black dashed], and for this parameter regime the
analytical and numerical results agree well.

A similar analysis can be done for the case of pure Dresselhaus with Zeeman,
$\beta=5.0$ meV nm and $g^{*}=-34$, see Fig. 12. Here the splitting is not
observed in the main peak, but rather in the second harmonic. As $\alpha$ is
increased from 0.0 to 9.0 meV nm, the splitting in the second harmonic
decreases, and vanishes at $\alpha=7.0$ meV nm, which again corresponds to the
condition in Eq. (61). Despite the good accuracy of the approximate analytical
solution for $\alpha\lesssim 8.0$ meV nm, it starts to deviate from the exact
one (full numerics) for higher values of $\beta$. This happens because for
these values, the combined effect Rashba and Dresselhaus is more pronounced,
producing an anti-crossing between different energy levels (see blue curves
Fig. 14, discussed further below). While the approximate energies obtained
here are always monotonic with respect to $1/B$, around the anti-crossing the
numerical ones are not. Accordingly, our $F$-function calculation will not be
able to fully describe the SdH oscillations and frequencies around the anti-
crossing regions, specifically the approximate solution misses a central peak
that starts developing, which will be discussed in the next section. In terms
of the $F$-function, the occurance of level anticrossings corresponds to
$|{\cal F}_{-}|\approx 1/2$. Since the power spectrum is obtained by
integrating $\delta\rho_{xx}$ over a range of $1/B$, there is no simple
condition determining the validity of the approximate solution. However,
looking at the $\lambda_{B}$ term in Eq. (59) the condition

$8\pi
n_{\mathrm{2D}}\left(\frac{k_{R}^{2}}{1-\tilde{\Delta}}-\frac{k_{D}^{2}}{1+\tilde{\Delta}}\right)l_{c}^{4}\lesssim
1,$ (64)

yields a useful estimate for the $1/B$ values where the Dingle factor has not
suppressed $\delta\rho_{xx}$. Equation (64) generalizes a similar condition
derived in Ref. Erlingsson  _et al._ , 2010. It is also interesting to note
that the analytical result is more accurate for higher harmonics, as the
Dingle-factor helps diminishing the amplitude of the anti-crossing at higher-
fields (see Fig. 14).

##  VI Landau Levels with Zeeman, Rashba and Dresselhaus interactions:
Numerical results

In the previous section, we have derived an approximate analytical result for
the magnetoresistance oscillations in the presence of both Rashba, Dresselhaus
and Zeeman interactions. The assumptions and approximations underlying the
derivation involved the relatively small SO coupling and the low number of
occupied Landau levels. These are satisfied in the low electron density InSb-
based 2DEGs of Refs. Gilbertsson  _et al._ , 2009; Akabori  _et al._ , 2008.
For higher electron density systems (but still with just a singly-occupied
subband at $B=0$), such as the InAs/GaSb wells in Ref. Beukman  _et al._ ,
2017, a numerical approach is needed. Below we outline the numerical
procedure. The numerical approach also allows us to account for the full form
of cubic Dresselhaus term, see Sec. VI.1.

For the case of either pure Rashba or Dresselhaus with Zeeman, the absence of
anti-crossing in the LL spectrum allow us to obtain _exact_ analytical results
for the problem. As we explain below, this does not hold in the presence of
both Rashba and Dresselhaus with the Hamiltonian (in the spin basis) Eq. (9)

$\frac{{\cal\tilde{H}}}{\hbar\omega_{c}}=\left(\begin{array}[]{cc}a^{\dagger}a+\frac{1}{2}+\frac{\tilde{\Delta}}{2}&2\alpha_{B}a+2\beta_{B}a^{\dagger}\\\
2\alpha_{B}a^{\dagger}+2\beta_{B}a&a^{\dagger}a+\frac{1}{2}-\frac{\tilde{\Delta}}{2}\end{array}\right).$
(65)

Therefore, here we calculate the magnetotransport numerically via the
diagonalization of the Hamiltonian above. The $F$-function method used for the
analytical cases can be extended to allow for numerical methods for
calculating the energy spectrum, see App. D.

As opposed to both the pure Rashba and pure Dresselhaus cases, ${\cal
N}_{\pm}$ do not commute with the Hamiltonian above, and therefore, the
diagonal basis cannot be described by any linear combination of the previous
degenerate eigenstates of ${\cal N}_{\pm}$. However, there is still a unitary
operator, ${\cal P}=\exp\left\\{i\pi\left({\cal
N}_{\pm}-\frac{1}{2}\right)\right\\}$ that commutes with this Hamiltonian,
called the _parity_ operator (Casanova  _et al._ , 2010; Braak, 2011), which
is discussed in detail in App. B. The corresponding unitary transformation
gives $\mathcal{P}a\mathcal{P}^{\dagger}=-a$,
$\mathcal{P}a^{\dagger}\mathcal{P}^{\dagger}=-a^{\dagger}$ and
$\mathcal{P}\sigma_{\pm}\mathcal{P}^{\dagger}=-\sigma_{\pm}$, which clearly
makes the Hamiltonian Eq. (9) invariant due to presence of only
$a^{\dagger}a$, $a^{\dagger}\sigma_{\pm}$ and $a\sigma_{\pm}$ terms. The
eigenvalues of $\mathcal{P}$, $\pm 1$, help analyze the energy spectrum
behavior.

To understand the influence on the spectrum of both Rashba and Dresselhaus
contributions, we first recall that in the absence of the latter, the Rashba
term is responsible for coupling $\left|n,\uparrow\right\rangle$ to
$\left|n+1,\downarrow\right\rangle$, for $n\in\mathbb{N}_{0}$ , thus yielding
decoupled $2\times 2$ block diagonal Rashba Hamiltonians (shown by the red
boxes in the Hamiltonian below). When we account for the Dresselhaus
contribution, we obtain a coupling between states
$\left|n,\downarrow\right\rangle$ and $\left|n+1,\uparrow\right\rangle$ for
$n\in\mathbb{N}_{0}$, which belongs to different Rashba blocks. More
specifically, the Dresselhaus term produces a coupling between blocks
$\left\\{\left|n,\uparrow\right\rangle,\left|n+1,\downarrow\right\rangle\right\\}$
and $\left\\{\left|n+\Delta n,\uparrow\right\rangle,\left|n+1+\Delta
n,\downarrow\right\rangle\right\\}$ with $\Delta n=2$, which is indicated by
the blue box in the Hamiltonian below (See App. B). As a consequence, we have
two decoupled orthogonal basis set given by
$\left\\{\left|0^{+}\right\rangle\right\\}=\left\\{\left|n,\uparrow\right\rangle,\left|n+1,\downarrow\right\rangle,\thinspace\ldots\right\\}$
and
$\left\\{\left|0^{-}\right\rangle\right\\}=\left\\{\left|n,\downarrow\right\rangle,\left|n+1,\uparrow\right\rangle,\thinspace\ldots\right\\}$
with $n\in\mathbb{N}_{0}$. Interestingly, these decoupled basis have different
eigenvalues with respect to the parity operator, i.e., ${\cal
P}\left|0^{\pm}\right\rangle=\pm 1\left|0^{\pm}\right\rangle$ and therefore,
represent different parity subspace.

Figure 13: Graphical illustration of the parity subspaces in the matrix
representation of the Hamiltonian Eq. (65). Here we see the Rashba interaction
couples $\left|n-1,\uparrow\right\rangle$ to $\left|n,\downarrow\right\rangle$
(red boxes), while the Dresselhaus interaction couples
$\left|n,\downarrow\right\rangle$ to $\left|n+1,\uparrow\right\rangle$ (blue
boxes).

In terms of the spectrum, in the presence of only Rashba SO coupling, we
observe multiple crossing between the Rashba eigenstates
$\\{|n,-\rangle,|n,+\rangle\\}$ for different $n\in\mathbb{N}_{0}$, with
energy given by Eq. (34), obtained through the diagonalization of the Rashba
blocks (red boxes within the Hamiltonian matrix in Fig. 13). This is shown by
the red solid lines in Fig. 14(a) for $\alpha=7.5$ meV nm. In the presence of
Dresselhaus SO coupling, the states $|n,-\rangle$ and $|n+\Delta n,+\rangle$
with $\Delta n\in\mathbb{N}_{\mathrm{odd}}$ belong to the same parity subspace
and adding a Dresselhaus contribution will yield anti-crossing, which open up
gaps in the spectrum (blue curves). Conversely, the decoupling between the
different parity sets, i.e., $|n,-\rangle$ and $|n+\Delta n,+\rangle$ with
$\Delta n\in\mathbb{N}_{\mathrm{even}}$, implies multiple crossing between
their corresponding energy states. These features are shown by the blue curve
in Figs. 14(a) and (c), where we have used $\beta=3.0$ meV nm. Other
parameters are $m^{*}=0.04$, $g^{*}=-12$ and $n_{2D}=17.6\times 10^{-3}$ nm-2.
These parameters are for InAs/GaSb-based (double) quantum wellsBeukman  _et
al._ (2017) in the electron regime. This regime, as emphasized in Ref. Beukman
_et al._ , 2017, corresponds to the configuration in which the GaSb well is
depleted and the system is effectively a single InAs-based asymmetric quantum
well with electrons only. Furthermore, we also observe that the effect of the
Dresselhaus term is to simply shift the crossing point to a different magnetic
field and energy (the crossing-point energy remains constant to lowest order
in $\beta$ but does in general shift for higher values of $\beta$).

The contrasting behavior of crossings vs. anti-crossings has direct
consequences on the $F$-function, which will be analyzed in the next
paragraphs. First we consider the crossing between states $|n,-\rangle$ and
$|n+\Delta n,+\rangle$, with even $\Delta n$ (corresponding to states
belonging to different parity subspaces). The $F$-function are

$\displaystyle\varepsilon_{n,-}(B)=\frac{\varepsilon}{\hbar\omega_{c}}\leftrightarrow
n=F_{-}\left(\frac{\varepsilon}{\hbar\omega_{c}},B;\alpha,\beta\right),$ (66)

and

$\displaystyle\varepsilon_{n+\Delta
n,+}(B)=\frac{\varepsilon}{\hbar\omega_{c}}\leftrightarrow n+\Delta
n=F_{+}\left(\frac{\varepsilon}{\hbar\omega_{c}},B;\alpha,\beta\right),$ (67)

where we have explicitly added their dependence on $\alpha$ and $\beta$. This
results in an $F$-function difference [see Eq. 13] at the crossing
$\varepsilon=\varepsilon_{c}$ and $B=B_{c}$

$\mathcal{F}_{-}\left(\frac{\varepsilon_{c}}{\hbar\omega_{c}},B_{c};\alpha,\beta\right)=\frac{\Delta
n}{2}\in\mathbb{Z}.$ (68)

Note that since the SdH oscillation is dependent on $\mathcal{F}_{\pm}$ in the
form of $\cos(2\pi\mathcal{F}_{\pm})$, we can re-define $\mathcal{F}_{-}$ to
lie within an unit interval, e.g., $\mathcal{F}_{-}\in[-1/2,1/2]$.
Accordingly, integer values of $\mathcal{F}_{-}$ are equivalent to
$\mathcal{F}_{-}=0$ and therefore, the vanishing of ${\cal F}_{-}$ provides
the field values where the crossing happens. The curves for ${\cal F}_{-}$ are
plotted in Fig. 14(b) for the same parameters as in Fig. 14(a). It presents a
sawtooth pattern because values of $|{\cal F}_{-}|>1/2$ are shifted back to
the $[-1/2,1/2]$ interval. The role of the Dresselhaus coupling for these
crossings is evident in Fig. 14(b), where the zeros of $\mathcal{F}_{-}$
remain zeros for any value of $\beta$, but are simply shifted to new values of
magnetic field, open circle moves to open rectangle Fig. 14(b).

Next, we look at the crossing between states belonging to the same parity
subspace, i.e., $|n,-\rangle$ and $|n+\Delta n,+\rangle$ for odd $\Delta n$.
We recall that this crossing only exists for the pure Rashba case, shown in
both Figs. 14(a) and (c). Here the relations in Eqs. (66) and (67) still hold,
the only difference being the value of $\Delta n$, which results in

$\mathcal{F}_{-}\left(\frac{\varepsilon_{c}}{\hbar\omega_{c}},B_{c};\alpha,\beta=0\right)=\frac{\Delta
n}{2}\in\mathbb{Z}+\frac{1}{2}.$ (69)

Adding a non-zero Dresselhaus contribution will couple these states and lead
to an anti-crossing, shown in Figs. 14(a) and (c). The anti-crossing result in
non half-integer values of $\mathcal{F}_{\pm}$ in Eqs. (66) and (67) and will
lead to a rounding of the sawtooth pattern as seen in Fig. (14)(b) (blue
curves).

The conditions in Eqs. (68) and (69) lead to values of
$\cos(2\pi\mathcal{F}_{-})=1$ [filled circle and rectangle in Fig. 14b)] and
$\cos(2\pi\mathcal{F}_{-})=-1$ [open cirlce circle in 14b)], respectively, in
the case of either pure Rashba or pure Dresselhaus. However, when both Rashba
and Dresselhaus are present only the former condition
$\cos(2\pi\mathcal{F}_{-})=1$ holds (crossing of states with opposite parity)
but the latter condition changes such that $\cos(2\pi\mathcal{F}_{-})>-1$ due
to anticrossings of states with same parity eigenvalue [open rectangle in Fig.
14b)]. This, in turn, affects the shape of the magneto-oscillations leading to
an asymmetry in the maximum and minimum values of $\cos(2\pi\mathcal{F}_{-})$.

Figure 14: The energy spectrum for two sets of $(\alpha,\beta)=(7.5,0.0)$ meV
nm [red] and $(\alpha,\beta)=(7.5,3.0)\,$meV nm [blue], along with
$\varepsilon_{F}/\hbar\omega_{c}$ (black dashed), a) around $n=125$ and c)
$n=255$. b) $\mathcal{F}_{c}$ for the same pair of parameters. Note sawtooth
form for pure Rashba [red], and for $(\alpha,\beta)=(7.5,3.0)\,$meV nm [blue]
a rounding, and translation, of the cusps due to level anticrossing [solid
circles]. Other parameters are $m^{*}=0.04$, $g^{*}=-12$ and
$n_{2D}=17.6\times 10^{-3}$ nm-2, for InAs-based quantum wellsBeukman  _et
al._ (2017).  Figure 15: Magnetooscillations for four different parameter
values, including pure Rashba, then different combinations of
$(\alpha,\beta)$. The anti-crossings in the spectrum complicates the beating
behavior, which eventually vanishes for around $(\alpha,\beta)=(4.5,3.0)\,$meV
nm. Other parameters are $m^{*}=0.04$, $g^{*}=-12$ and $n_{2D}=17.6\times
10^{-3}$ nm-2.

In Fig. 15 this asymmetry is visible in the magneto-osillations. Here we
assume Gaussian broadening with $B_{q}=0.50$ T which forms an envelope (black
dashed curve). The lowest curve is the pure Rashba
$(\alpha,\beta)=(7.5,0.0)\,$meV nm and there all maximas intersect the
envelope [black circles]. The curve for $(\alpha,\beta)=(7.5,3.0)\,$meV nm
shows that only some maxima intersect the envelope, the other maximas
correspond to $\cos(2\pi\mathcal{F}_{-})>-1$ do not (black circle). This is a
direct consequence of the anti-crossing in the spectrum in Fig. 14. The curves
for $(\alpha,\beta)=(5.5,3.0)\,$meV nm and $(4.5,3.0)\,$meV nm show how the
anti-crossing becomes larger, eventually leading to an absence of beatings.
This can also be seen in the frequency spectrum shown in Fig. 16, for the
$f\approx f^{\mathrm{SdH}}$ peak. The lowest curve (blue) corresponds to
$(\alpha,\beta)=(7.5,3.0)\,$meV nm where the spectrum shows well separated
peaks. However, as the strength of the Rashba coupling is decreased all the
way down to $\alpha=0.5$ meV nm for a fixed value of $\beta=3.0$ meV nm a
central peak develops and for $\alpha$ between 4.5 and 1.5 meV nm, the two
split peaks are barely visible.

Figure 16: Power spectrum for af fixed $\beta=3.0$ meV nm for $\alpha=7.5$
down to 0.5 meV nm. Other parameters are $m^{*}=0.04$, $g^{*}=-12$ and
$n_{2D}=17.6\times 10^{-3}$ nm-2, from Ref. Beukman  _et al._ , 2017.

###  VI.1 Extracting $\alpha$ and $\beta$ from SdH data

The magneto-oscillations can be thought of as a fingerprint of the sample
parameters, including Fermi energy $\varepsilon_{F}$, effective mass $m^{*}$,
$g^{*}$, and $\alpha$ and $\beta$. To better capture the influence of the
spin-orbit couplings for higher electron density, the full form of the
Dresselhaus interaction will be used. For non-zero magnetic fields, this
corresponds to having Dresselhaus SO term in Eq. (9) replaced with

$\displaystyle\left[\frac{1}{\sqrt{2}\hbar\omega_{c}\ell_{c}}\left(\beta_{1}-\gamma\frac{a^{\dagger}a}{2\ell_{c}^{2}}\right)a^{\dagger}\sigma_{+}+\frac{\gamma}{2\ell_{c}^{2}}a^{3}\sigma_{-}\right]+\mathrm{h.c.},$
(70)

where $\beta_{1}=\gamma\langle k_{z}^{2}\rangle$, $\gamma$ is material-
dependent parameter describing the SO interaction due to bulk inversion
asymmetry, and $\langle k_{z}^{2}\rangle$ is the expectation value of the
$z$-component of the square of momentum operator (divided by $\hbar$), see
App. C for details of full Dresselhaus coupling. Note that $\beta$ in Eq. (2)
is assumed to include the first harmonic of the cubic Dresselhaus Fu  _et al._
(2016); Dresselhaus (1955), which makes it linearly dependent on the electron
density. For instance, if the potential confining the 2DEG is assumed to be an
infinite well of width $d_{\mathrm{QW}}$ then $\langle
k_{z}^{2}\rangle=\pi^{2}/d_{\mathrm{QW}}^{2}$. To model the magnetoresistance
data we start from Eq. (22), which features $(i)$ a sum over higher harmonics
, $(ii)$ rapid oscillations coming from $\mathcal{F}_{+}$, and $(iii)$ damping
due to Landau level broadening $\tilde{L}_{\Gamma}$. The analysis introduced
in the previous section was based on the study of the properties of
$\cos(2\pi\mathcal{F}_{-})$, which forms an envelope on top of the rapid
oscillations. Note that in the case having both Rashba and Dresselhaus
coupling the rapid oscillations are still dominated by the normal SdH
oscillations, i.e.

$\displaystyle\mathcal{F_{+}}(B)$ $\displaystyle=$
$\displaystyle-\frac{1}{2}+\frac{\varepsilon_{F}}{\hbar\omega_{c}}\left(1+\mathcal{O}\left(\frac{\varepsilon_{R}}{\varepsilon_{F}},\frac{\varepsilon_{D}}{\varepsilon_{F}}\right)\right)$
(71) $\displaystyle\approx$
$\displaystyle-\frac{1}{2}+\frac{f^{\mathrm{SdH}}}{B},$

so the SO coupling does not affect the rapid oscillations. The resulting
lowest harmonic form of the magneto-resistivity is

$\delta\rho_{xx}(B)=-2\tilde{L}_{\Gamma}(B)\cos(2\pi
l\mathcal{F}_{-}(B))\cos\left(2\pi\frac{f_{\mathrm{SdH}}}{B}\right),$ (72)

which can be fitted to available data.

Figures 17-19 show the experimental data from Ref. Beukman  _et al._ , 2017
for InAs/GaSb quantum wells in the electron regime) along with our theoretical
fits [Eq. (72)]. We focus on the experimental curves 1, 5 and 10, of Fig. S4
of Ref. Beukman  _et al._ , 2017 that we label as C1, C5 and C10 in Fig.
17-19. The data was fitted to $\delta\rho_{xx}(B)$ in Eq. (72), where
$\mathcal{F}_{-}$ was calculated numerically. For the fitting we consider both
the Dresselhaus coupling in Eq. (9) [black dashed lines], and also with the
full Dresselhaus term in Eq. (70) [solid red lines]. The black dots are
reference points extracted from the data, which are used in the fitting of
$\tilde{L}_{\Gamma}(B)\cos(2\pi\mathcal{F}_{-})$. The best fittings were
produced by assuming Gaussian broadening, namely.

$\tilde{L}_{\Gamma}(B)=\exp\left(-2\pi^{2}\frac{\Gamma^{2}}{(\hbar\omega_{c})^{2}}\right)=\exp\left(-\frac{B_{q}^{2}}{B^{2}}\right),$
(73)

where $B_{q}=\sqrt{2}\pi\frac{m^{*}\Gamma}{\hbar e}$ and $\Gamma$ is a
constant Landau level broadening.

Figure 17: The black dots are reference points for curve C1, solid black line.
The black dashed curve is the linear Dresselhaus result and solid red curve
full Dresselhaus result. Parameter values from fitting are shown in the inset.
Other parametersBeukman  _et al._ (2017) are $m^{*}=0.019$, $g=-12$ and
$n_{\mathrm{2D}}=0.0176$ nm-2. Figure 18: Similar to Fig. (17), but for curve
C5, solid black line. The black dashed curve represents the linear Dresselhaus
result, while the solid red curve, the full cubic Dresselhaus term. Extracted
fitting parameters are shown in the inset. Other parametersBeukman  _et al._
(2017) are $m^{*}=0.019$, $g=-12$ and $n_{\mathrm{2D}}=0.0176$ nm-2.  Figure
19: Similar to Fig. (17), but for curve C10, solid black line. The black
dashed curve represents the linear Dresselhaus result, which fails to fit the
data. However, the full cubic Dresselhaus term (solid red curve) results in a
good fit. Extracted fitting parameters are shown in the inset. Other
parametersBeukman  _et al._ (2017) are $m^{*}=0.019$, $g=-12$ and
$n_{\mathrm{2D}}=0.0176$ nm-2.

For curve C1 in Fig. 17 the fitting with linear Dresselhaus yields values
$\alpha=7.2$ meV nm and $\beta=3.0$ meV nm. On the other hand, for fitting to
the full model we obtain $\alpha=7.6$ meV nm, and $\gamma=85$ meV nm3. We see
that both fits produce equally good curves fitting the experimental data
points, with comparable values for the extracted Rashba SO coupling. This
indicates that when the Rashba coupling dominates the cubic Dresselhaus term
($a^{3}$-term in Eq. (70)), fitting the data with the addition of the cubic
term does not strongly affect the result. The results for curve C5 in Fig. 18
behave similarly, i.e. we find fitted values of the Rashba coefficient,
$\alpha=6.7$ meV nm for the linear Dresselhaus with $\beta=2.8$ meV nm, and
$\alpha=6.3$ meV nm for the full cubic Dresselhaus, with $\gamma=82$ meV nm3.

However, the story is different for the curve C10 shown in Fig. 19. Here the
value of Rashba and Dresselhaus coupling are closer, and then the details of
the linear vs. cubic Dresselhaus become relevant. Indeed, the linear
Dresselhaus model fitting yields $\alpha=5.5$ meV nm and $\beta=2.5$ meV nm
while the cubic fit gives $\alpha=4.9$ meV nm. More importantly the error in
the linear fit is quite high, and the fit [black dashed curve] fails to
describe the data points. However, the cubic model gives a good fit , with
$\gamma=80$ meV nm3. This clearly shows the importance of the cubic
contributions in samples with high density, where the Rashba and Dresselhaus
contributions are of comparable magnitudes.

The fit results in Fig. 17-19 were done for $\langle
k_{z}^{2}\rangle=\pi^{2}/d_{\mathrm{QW}}^{2}$ where $d_{\mathrm{QW}}=12.5$ nm
Beukman  _et al._ (2017). To fully model the sample a self-consistent Poisson-
Schrödinger calculation is requiredDettwiler  _et al._ (2017); Fu  _et al._
(2020, 2016), which is beyond the scope of this work. We can however use
different values of $\langle k_{z}^{2}\rangle$, which indirectly emulate self-
consistent potential details, i.e. increasing the value of $\langle
k_{z}^{2}\rangle$ suggests a stronger confinement in the InAs quantum well,
and decreased value of $\langle k_{z}^{2}\rangle$ would correspond to
wavefunctions being less localized in the InAs quantum well.

In Fig. 20 the values of $\alpha$, $\beta_{1}$, and $\beta$ are shown as a
function of $\langle k_{z}^{2}\rangle$ from
0.75$\frac{\pi^{2}}{d_{\mathrm{QW}}^{2}}$ to
1.25$\frac{\pi^{2}}{d_{\mathrm{QW}}^{2}}$. The data from the three curves are
indicated by different forms: C1: circle, C5: triangle, and C10: square. For
each value of $\langle k_{z}^{2}\rangle$, specific values of $\alpha$
$\beta_{1}$, and $\beta$ are obtained from the fit. The fit results for
$\alpha$ and $\beta$ for each curve remain relatively insensitive to $\langle
k_{z}^{2}\rangle$-variations. Note that as $\langle k_{z}^{2}\rangle$ varies
$\beta_{1}$ changes quite rapidly via the fitted value of $\gamma$. This is to
be expected since lower values $\langle k_{z}^{2}\rangle$, correspond to the
electron leaking out the InAs quantum well $\gamma$ into the GaSb, which has a
higher bulk value of $\gamma$. For higher values of $\langle k_{z}^{2}\rangle$
the system becomes more strongly confined in the InAs quantum well and the
value of $\gamma$ should tend to the value corresponding to bulk InAs.

The fact that the values of $\alpha$ and $\beta$ change only slightly as
function of $\langle k_{z}^{2}\rangle$, as can be seen in Fig. 20, has
important consequences on the fitting proceedure. For this reason a fitting
with $\gamma$ and $\langle k_{z}^{2}\rangle$ both being independent fitting
parameters can not be performed, since if $\beta_{1}=\gamma\langle
k_{z}^{2}\rangle$ is the dominant contribution to the Dresselhaus couplings
then there are multiple (infinite) solutions to the equation $\gamma\langle
k_{z}^{2}\rangle=\mathrm{const.\ }$ and fitting the data with $\gamma$ and
$\langle k_{z}^{2}\rangle$ independent will not converge Beukman  _et al._
(2017).

Figure 20: The spin-orbit parameters that result from the fitting as a
function of $\langle k_{z}^{2}\rangle$. Other parameters are $m^{*}=0.019$,
$g=-12$ and $n_{\mathrm{2D}}=0.0176$ nm-2. The three different symbols
represent different curves: curve 1: circle, 5: triangle, and 10: square.

##  VII Summary

We investigated the SdH magneto-oscillations in the resistivity $\rho_{xx}$ of
2DEGs in the presence of spin-orbit (Rashba-Dresselhaus) and Zeeman couplings.
We used a semiclassical approach for the resistivity combined with a Poisson
summation formula for the Landau-quantized DOS. Our approach allows for an
intuitive separation of the slow and fast quantum oscillations in terms of
“F-functions”, central quantities in our description, essentially being the
inverse functions of the spin-resolved Landau-level structure of the system.
We study a variety of exact cases such as the pure Zeemann, pure Dresselhaus,
and pure Rashba cases – all of which provide analytical expressions for the
magnetoresistivity.

More importantly, from our unified and general formulation we also derive, for
the first time, an analytical solution for the case with arbitrary Rashba and
Dresselhaus couplings and simultaneous non-zero Zeeman coupling ($g^{*}\neq
0$).Interestingly, this allows us to derive a unique new condition for the
vanishing of the SO-induced beatings in the SdH signals:
$\alpha/\beta=[(1-\tilde{\Delta})/(1+\tilde{\Delta})]^{1/2}$, where
$\tilde{\Delta}=g^{*}m^{*}/2m_{0}$ (i.e., ratio (Zeeman
energy)/$\hbar\omega_{c}$). This new condition does not correspond to any
conserved quantity in our Hamiltonian, unlike the persistent-spin-helix
condition $\alpha=\beta$ which is associated with the conservation of spin
along some particular axes. We emphasize that our new condition precludes
beatings in all harmonics of the quantum oscillations.

We have applied our analytical formulation to describe low-density data for
SdH oscillations showing many harmonics in GaAs-based 2DEGs (see SM in Ref.
Dettwiler  _et al._ , 2017) and found an excellent agreement, Fig. 2. We have
also applied our theory to low-density InSb-based 2DEGsGilbertsson  _et al._
(2008); Akabori  _et al._ (2008). In addition, we have also developed a
detailed numerical calculation for high-density InAs-based 2DEGs, in which an
analytical description is not satisfactory. We find excellent agreement with
available data for high-density InAs-based 2DEGs Beukman  _et al._ (2017);
Dettwiler  _et al._ (2017). We have also pointed out an inequivalence between
the Rashba-dominated + Zeeman vs the Dresselhaus-dominate + Zeeman cases, with
only the former showing beatings. This follows from a distinct interplay
between the SO and Zeeman terms in these two regimes. We hope our detailed
study and unified general formulation will stimulate furhter experimental
investigations aiming at verifying our theoretical predictions.

##  VIII Acknowledgments

The authors would like to thank Arjan Beukman and Leo Kouvenhowen for sharing
experimental data from Ref. Beukman  _et al._ , 2017 that was used for
fitting. We also thank Thomas Schaepers and Makoto Kohda for useful
discussions. The authors acknowledge funding from the Reykjavik University
Research Fund, the São Paulo Research Foundation (FAPESP) Grants No.
2016/08468-0 and No. 2020/00841-9, Conselho Nacional de Pesquisas (CNPq),
Grants No. 306122/2018-9 and 301595/2022-4, the Swiss Nanoscience Institute
(SNI), the NCCR SPIN and grant no. 179024 of the Swiss NSF, the Georg H.
Endress Foundation, the EU H2020 European Microkelvin Platform EMP (grant no.
824109) and FET TOPSQUAD (grant no. 862046).

##  Appendix A Density of states and F-Functions

Here we follow closely the discussion (and notation) in Sec. 3.2.2 of the book
Semiclassical Physics by Brack and Bhaduri Brack and Bhaduri (1997).

For simplicity, we first consider the case with a discrete spectrum
$\varepsilon_{n}=f(n)$, $n=0,1,2,...$ in which each level has a degeneracy
$d_{n}=\tilde{D}(n)$, with $\tilde{D}(n)$ being an arbitrary function of $n$.
Later on we will account for a (pseudo) spin index. As an example, we note
that for the usual 2DEG Landau levels (LLs) (in the absence of both Zeeman or
SO interaction), $\varepsilon_{n}=\hbar\omega_{c}(n+1/2)$ and
$d_{n}=AeB/h=\tilde{D}(n)$ (A: area of the 2DEG, $e>0$); in this case,
$d_{n}=\tilde{D}(n)$ denotes the LL degeneracy and is independent of $n$. This
same Landau degeneracy holds in the presence of Zeeman and SO interactions.
For later convenience, we define $D(n)=\tilde{D}(n)/A$ to be the level
degeneracy per unit area [e.g., for LLs $D(n)=n_{LL}(B)=eB/h$]. As in Ref.
Brack and Bhaduri, 1997, let $f(n)$ be an arbitrary monotonic function with a
differentiable inverse $f^{-1}(x)=F(x)$, the relevant “F-function” in our
discussion. In this case, because $f^{-1}(f(x))=x=f(f^{-1}(x))$ it follows
that $n=F(\varepsilon_{n})$. Next we define the DOS of our system and relate
it the to the F-function, which ultimately allows us to calculate the
oscillatory part of the DOS relevant for our semiclassical transport
calculation.

###  A.1 Density of states without LL broadening

Quite generally we can define the DOS of our system as,

$g(\varepsilon)=\frac{1}{A}\sum_{n=0}^{\infty}\tilde{D}(n)\delta(\varepsilon-\varepsilon_{n}).$
(74)

Note that the above DOS is defined per area and energy. In Ref. Brack and
Bhaduri (1997) the DOS is defined just per energy. Motivated by the property
$\delta[{{\color[rgb]{0,0,0}y}}(x)]=\frac{1}{|{{\color[rgb]{0,0,0}y}}^{\prime}(x_{0})|}\delta(x-x_{0})$
where $x_{0}$ denotes a root of ${{\color[rgb]{0,0,0}y}}(x)$, i.e.,
${{\color[rgb]{0,0,0}y}}(x_{0})=0$ and
${{\color[rgb]{0,0,0}y}}^{\prime}(x)=d{{\color[rgb]{0,0,0}y}}(x)/dx$, we
define $h(\varepsilon)=n-F(\varepsilon)$, which obeys $h(\varepsilon_{n})=0$
as $n=F(\varepsilon_{n})$ by construction. We can then write

$\delta[h(\varepsilon)]=\delta(n-F(\varepsilon))=\frac{1}{|F^{\prime}(\varepsilon_{n})|}\delta(\varepsilon-\varepsilon_{n})=\frac{1}{|F^{\prime}(\varepsilon)|}\delta(\varepsilon-\varepsilon_{n}),$
(75)

or

$\delta(\varepsilon-\varepsilon_{n})=|F^{\prime}(\varepsilon)|\delta(n-F(\varepsilon)).$
(76)

Substituting (76) into (74), we find

$\displaystyle g({{\color[rgb]{0,0,0}\varepsilon}})$ $\displaystyle=$
$\displaystyle
D(\varepsilon)|F^{\prime}(\varepsilon)|\sum_{n=0}^{\infty}\delta(n-F(\varepsilon)),$
(77)

where $D(\varepsilon)\equiv D(F(\varepsilon))$. Noting that

$\sum_{n=0}^{\infty}\delta(\varepsilon-n)=\sum_{l=-\infty}^{\infty}e^{2\pi
il\varepsilon},\qquad(\varepsilon>0)$ (78)

we can straightforwardly cast (77) in the form

$g(\varepsilon)=D(\varepsilon)|F^{\prime}(\varepsilon)|\sum_{l=-\infty}^{\infty}e^{2\pi
ilF(\varepsilon)}.$ (79)

Now we introduce the (pseudo) spin index $s=\pm 1$ by adding a subscript $s$
to all quantities [except $D(\varepsilon)$ for it is not (pseudo) spin
dependent]. This $s$ index accounts for the spin-dependent Zeeman and SO
interactions in our 2DEG. With this new index, the DOS in Eq. (78), viewed as
per spin now, becomes

$g_{s}(\varepsilon)=D(\varepsilon)|F_{s}^{\prime}(\varepsilon)|\sum_{l=-\infty}^{\infty}e^{2\pi
ilF_{s}(\varepsilon)},$ (80)

or

$g_{s}(\varepsilon)=D(\varepsilon)|F_{s}^{\prime}(\varepsilon)|\left\\{1+2\sum_{l=1}^{\infty}\cos[2\pi
lF_{s}(\varepsilon)]\right\\}.\qquad(\varepsilon>0)$ (81)

By summing over $s$, we obtain the total DOS,

$g(\varepsilon)=D(\varepsilon)\sum_{s}|F_{s}^{\prime}(\varepsilon)|\left\\{1+2\sum_{l=1}^{\infty}\cos[2\pi
lF_{s}(\varepsilon)]\right\\}.$ (82)

For the systems investigated in our work, $F_{s}^{\prime}(\varepsilon)\simeq
1/\hbar\omega_{c}$. This is actually exact for the Zeeman-only case, see Eq.
(28), main text, but only approximate in the presence of SO interaction [see
Eq. (111)]. In this case and
$D(\varepsilon)|F_{s}^{\prime}(\varepsilon)|=\frac{m^{*}}{2\pi\hbar^{2}}$, we
find

$g(\varepsilon)\simeq\frac{m^{*}}{\pi\hbar^{2}}\left\\{1+\sum_{l=1}^{\infty}\left(\cos[2\pi
lF_{+}(\varepsilon)]+\cos[2\pi lF_{-}(\varepsilon)]\right)\right\\}.$ (83)

Using the identity,

$\cos a+\cos b=2\cos[(a+b)/2]\cos[(a-b)/2],$ (84)

we can rewrite Eq. (83) as

$g(\varepsilon)\simeq\frac{m^{*}}{\pi\hbar^{2}}\left\\{1+\sum_{l=1}^{\infty}{{\color[rgb]{0,0,0}2}}\cos[2\pi
l\mathcal{F}_{+}(\varepsilon)]\cos[2\pi
l\mathcal{F}_{-}(\varepsilon)]\right\\},$ (85)

where

$\mathcal{F}_{\pm}(\varepsilon)=\frac{1}{2}[F_{+}(\varepsilon)\pm
F_{-}(\varepsilon)].$ (86)

To regain the DOS notation in the main text, we now make
$g(\varepsilon)\rightarrow\mathcal{D}(\varepsilon,B)$ and use
$\mathcal{D}_{0}=\frac{m^{*}}{2\pi\hbar^{2}}$. Hence, Eq. (85) becomes

$\mathcal{D}(\varepsilon,B)\simeq
2\mathcal{D}_{0}\left\\{1+{{\color[rgb]{0,0,0}2}}\sum_{l=1}^{\infty}\cos[2\pi
l\mathcal{F}_{+}(\varepsilon)]\cos[2\pi
l\mathcal{F}_{-}(\varepsilon)]\right\\},$ (87)

or

$\frac{\mathcal{D}(\varepsilon,B)-{{\color[rgb]{0,0,0}2}}\mathcal{D}_{0}}{{{\color[rgb]{0,0,0}2}}\mathcal{D}_{0}}\simeq
2\sum_{l=1}^{\infty}\cos[2\pi l\mathcal{F}_{+}(\varepsilon)]\cos[2\pi
l\mathcal{F}_{-}(\varepsilon)],$ (88)

which is Eq. (14) in the main text.

###  A.2 Density of states including Landau level broadening

We can account for LL broadening in the DOS calculation by considering
Lorentzian or Gaussian functions as particular representations of the ideal
$\delta$ functions describing the discrete levels. We consider a simple
phenomenological description which assumes that all LLs have the same spin-
independent broadening $\Gamma$.

####  A.2.1 Lorentzian DOS case

Here we take the delta function representing the DOS of a single LL as,

$\delta(\varepsilon-\varepsilon_{n})=\lim_{\Gamma\rightarrow
0}\frac{1}{\pi}\frac{\Gamma/2}{(\varepsilon-\varepsilon_{n})^{2}{{\color[rgb]{0,0,0}+}}(\Gamma/2)^{2}}=\lim_{\Gamma\rightarrow
0}L_{\Gamma}(\varepsilon-\varepsilon_{n}),$ (89)

where

$L_{\Gamma}(\varepsilon)=\frac{1}{\pi}\frac{\Gamma/2}{\varepsilon^{2}{{\color[rgb]{0,0,0}+}}(\Gamma/2)^{2}},$
(90)

with

$\int_{-\infty}^{\infty}L_{\Gamma}(\varepsilon)d\varepsilon=1.$ (91)

Note that

$\int_{-\infty}^{\infty}L_{\Gamma}(\varepsilon)e^{-2\pi
il\varepsilon}d\varepsilon=e^{-\Gamma\pi|l|}=\tilde{L}_{\Gamma}(k),$ (92)

where $\tilde{L}_{\Gamma}(k)$ is the Fourier transform (FT) of
$L_{\Gamma}(\varepsilon)$ and $l\in\mathbb{Z}$. Using the shifting property of
FTs, it follows that the FT of $L_{\Gamma}(\varepsilon-x)$ is $e^{-2\pi
ikx}\tilde{L}_{\Gamma}(k)$. Generalizing Eq. (74) for Lorentzian-broadened
levels we have (we will add a subindex $s$ later on)

$g(\varepsilon)=\lim_{\Gamma\rightarrow
0}\sum_{n=0}^{\infty}D(n)L_{\Gamma}(\varepsilon-\varepsilon_{n}),$ (93)

which we can rewrite as,

$g(\varepsilon)=\lim_{\Gamma\rightarrow
0}\sum_{n=0}^{\infty}\int_{-\infty}^{\infty}{D(n)}L_{\Gamma}(\varepsilon-x)\delta(x-\varepsilon_{n})dx.$
(94)

Considering that $D(n)$ is independent of $n$ and using Eq. (78) with the
replacement $\varepsilon\rightarrow F(\varepsilon)$, we obtain

$g(\varepsilon)=\lim_{\Gamma\rightarrow
0}\int_{-\infty}^{\infty}D(F(x))|F^{\prime}(x)|\sum_{l=-\infty}^{\infty}e^{2\pi
ilF(x)}L_{\Gamma}(\varepsilon-x)dx.$ (95)

Since $L_{\Gamma}(\varepsilon-x)=L_{\Gamma}(x-\varepsilon)$ is peaked at
$x=\varepsilon$, it is convenient to expand $F(x)$ around this point. Then
$g(\varepsilon)$ becomes

$\displaystyle g(\varepsilon)=\lim_{\Gamma\rightarrow
0}\int_{-\infty}^{\infty}D[F(\varepsilon)+|F^{\prime}(\varepsilon)|(x-\varepsilon)][|F^{\prime}(\varepsilon)|+$
$\displaystyle|F^{\prime\prime}(\varepsilon)|(x-\varepsilon)]\sum_{l=-\infty}^{\infty}e^{2\pi
il\left[F(\varepsilon)+|F^{\prime}(\varepsilon)|(x-\varepsilon)\right]}L_{\Gamma}(\varepsilon-x)dx.$
(96)

Neglecting the contribution $|F^{\prime\prime}(\varepsilon)|(x-\varepsilon)$
[as a matter of fact, this contribution vanishes identically in the limit
$L_{\Gamma}(x-\varepsilon)\rightarrow\delta(x-\varepsilon)$, because
$\int_{-\infty}^{\infty}f(x)\delta(x-x_{0})dx=f(x_{0})$], we have

$\displaystyle g(\varepsilon)=D(\varepsilon)\lim_{\Gamma\rightarrow
0}|F^{\prime}(\varepsilon)|\sum_{l=-\infty}^{\infty}e^{2\pi
ilF(\varepsilon)}\times$ $\displaystyle\int_{-\infty}^{\infty}e^{2\pi
il|F^{\prime}(\varepsilon)|(x-\varepsilon)}L_{\Gamma}(x-\varepsilon)d(x-\varepsilon).$
(97)

Using Eq. (92), we can write

$g(\varepsilon)=D(\varepsilon)\lim_{\Gamma\rightarrow
0}|F^{\prime}(\varepsilon)|\sum_{l=-\infty}^{\infty}e^{2\pi
ilF(\varepsilon)}\tilde{L}_{\Gamma}(l|F^{\prime}(\varepsilon)|).$ (98)

or

$g(\varepsilon)=D(\varepsilon)\lim_{\Gamma\rightarrow
0}|F^{\prime}(\varepsilon)|\sum_{l=-\infty}^{\infty}e^{2\pi
ilF(\varepsilon)}e^{-\Gamma\pi|lF^{\prime}(\varepsilon)|},$ (99)

where have used,

$\tilde{L}_{\Gamma}(l|F^{\prime}(\varepsilon)|)=e^{-\Gamma\pi|lF^{\prime}(\varepsilon)|}$
(100)

As before [Eq. (80)], we can rewrite Eq. (99) by adding a subindex $s$ to
obtain the LL-broadened DOS per spin

$g_{s}(\varepsilon)=D(\varepsilon)|F_{s}^{\prime}(\varepsilon)|\left\\{1+2\sum_{l=1}^{\infty}\cos[2\pi
lF_{s}(\varepsilon)]e^{-\Gamma\pi l|F_{s}^{\prime}(\varepsilon)|}\right\\}.$
(101)

In the above we have dropped the $\lim_{\Gamma\rightarrow 0}$, since a real
system has a finite $\Gamma$. Interestingly, the broadened DOS in Eq. (101)
can be obtained directly from the case without broadening [Eq. (81)] by simply
multiplying the oscillating components (harmonics) in the latter by the
exponential (“Dingle”) factor $e^{-\Gamma\pi l|F_{s}^{\prime}(\varepsilon)|}$.

Here again, for the systems of interest here
$F_{s}^{\prime}(\varepsilon)\simeq 1/\hbar\omega_{c}$ and the exponential
factor in Eq. (101) becomes

$e^{-\Gamma\pi l|F_{s}^{\prime}(\varepsilon)|}=e^{-\pi
l\Gamma/\hbar\omega_{c}},$ (102)

where $\Gamma\equiv\hbar/\tau_{q}$, $\tau_{q}$ is the quantum lifetime of the
LL. Summing over the (pseudo) spin index $s$ Eq. (101) becomes

$g(\varepsilon)=\frac{m^{*}}{\pi\hbar^{2}}\left\\{1+{{\color[rgb]{0,0,0}2}}\sum_{l=1}^{\infty}\cos[2\pi
l\mathcal{F}_{+}(\varepsilon)]\cos[2\pi
l\mathcal{F}_{-}(\varepsilon)]e^{-{{\color[rgb]{0,0,0}\frac{\pi
l\Gamma}{\hbar\omega_{c}}}}}\right\\}.$ (103)

In the notation of the main text we have

$\frac{\mathcal{D}(\varepsilon,B)-{{\color[rgb]{0,0,0}2}}\mathcal{D}_{0}}{{{\color[rgb]{0,0,0}2}}\mathcal{D}_{0}}\simeq
2\sum_{l=1}^{\infty}\cos[2\pi l\mathcal{F}_{+}(\varepsilon)]\cos[2\pi
l\mathcal{F}_{-}(\varepsilon)]e^{-{{\color[rgb]{0,0,0}\frac{\pi
l\Gamma}{\hbar\omega_{c}}}}},$ (104)

which is the Eq. (22) of the main text, but written for the Lorentzian
broadening case.

####  A.2.2 Gaussian DOS case

The Gaussian-broadened case can be treated similarly by considering the delta
function representation

$\delta(\varepsilon-\varepsilon_{n})=\lim_{\Gamma\rightarrow
0}\frac{1}{\sqrt{2\pi}\Gamma}e^{-\frac{(\varepsilon-\varepsilon_{n})^{2}}{2\Gamma^{2}}}.$
(105)

From this we can evaluate the integral in Eq. (92) which results in the
Gaussian version of Eq. (100):

$\tilde{L}_{\Gamma}\left(l|F^{\prime}(\epsilon)|\right)=e^{-2\pi^{2}(l|F^{\prime}(\epsilon)|)^{2}\Gamma^{2}}.$
(106)

This reduces to Eq. (73) for $l=1$ (fundamental frequency) and
$|F^{\prime}(\epsilon)|=1/\hbar\omega_{c}$.

####  A.2.3 Calculating the F-function and its derivative
$F^{\prime}(\varepsilon)$

Here we illustrate the calculation of $F_{s}(\varepsilon)$ and its derivative
with respect to $\varepsilon$, $F^{\prime}(\varepsilon)$, in the presence of
SO interaction. For simplicity, we consider the pure Rashba case (no Zeeman).
To determine the F-functions we need to invert
${{\color[rgb]{0,0,0}\varepsilon_{n,s}}}=\varepsilon$, where

$\displaystyle\frac{\varepsilon_{n,s}}{\hbar\omega_{c}}=$
$\displaystyle\left(n+\frac{1}{2}+\frac{s}{2}\right)$ (107)
$\displaystyle-\frac{s}{2}\frac{1-\tilde{\Delta}}{|1-\tilde{\Delta}|}\sqrt{\left(1-\tilde{\Delta}\right)^{2}+16\alpha_{B}^{2}\left(n+\frac{1}{2}+\frac{s}{2}\right)},$

is the pure Rashba energy, Eq. (33) in the main text. Squaring
$\varepsilon-\tilde{n}\hbar\omega_{c}$, with $\tilde{n}=n+(1+s)/2$, we find

$\displaystyle[\varepsilon-\tilde{n}\hbar\omega_{c}]^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{4}(\hbar\omega_{c}-\Delta)^{2}+4\varepsilon_{R}\hbar\omega_{c}\tilde{n}$
$\displaystyle\varepsilon^{2}-2\varepsilon\hbar\omega_{c}\tilde{n}$
$\displaystyle+$
$\displaystyle\tilde{n}^{2}\hbar^{2}\omega_{c}^{2}=\frac{1}{4}(\hbar\omega_{c}-\Delta)^{2}+4\varepsilon_{R}\hbar\omega_{c}\tilde{n}$
$\displaystyle\tilde{n}^{2}\hbar^{2}\omega_{c}^{2}-(2\varepsilon\hbar\omega_{c}$
$\displaystyle+$ $\displaystyle
4\varepsilon_{R}\hbar\omega_{c})\tilde{n}-\frac{1}{4}(\hbar\omega_{c}-\Delta)^{2}+\varepsilon^{2}=0$
$\displaystyle\tilde{n}^{2}-\left(\frac{2\varepsilon}{\hbar\omega_{c}}+\frac{4\varepsilon_{R}}{\hbar\omega_{c}}\right)\tilde{n}$
$\displaystyle-$
$\displaystyle\left(\frac{1}{2}-\frac{\Delta}{2\hbar\omega_{c}}\right)^{2}+\left(\frac{\varepsilon}{\hbar\omega_{c}}\right)^{2}=0$
(108)

We can easily solve (108) for $\tilde{n}_{s}(\varepsilon)\Rightarrow
n_{s}(\varepsilon)=-(1+s)/2+\tilde{n}_{s}(\varepsilon)=f_{s}^{-1}=F_{s}(\varepsilon)$

$\displaystyle F_{s}(\varepsilon)$ $\displaystyle=$
$\displaystyle-\frac{1+s}{2}+\frac{\varepsilon}{\hbar\omega_{c}}+\frac{2\varepsilon_{R}}{\hbar\omega_{c}}$
$\displaystyle+$ $\displaystyle
s\sqrt{\left(\frac{\varepsilon}{\hbar\omega_{c}}+\frac{2\varepsilon_{R}}{\hbar\omega_{c}}\right)^{2}+\left(\frac{1}{2}-\frac{\Delta}{2\hbar\omega_{c}}\right)^{2}-\left(\frac{\varepsilon}{\hbar\omega_{c}}\right)^{2}}.$

We obtain $F_{\pm}^{\prime}(\varepsilon)$ by differentiating (A.2.3),

$\displaystyle F_{s}^{\prime}(\varepsilon)$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar\omega_{c}}+s\frac{1}{2}\frac{2\left(\frac{\varepsilon}{\hbar\omega_{c}}+\frac{2\varepsilon_{R}}{\hbar\omega_{c}}\right)\frac{1}{\hbar\omega_{c}}-\frac{2\varepsilon}{\hbar^{2}\omega_{c}^{2}}}{\sqrt{\left(\frac{\varepsilon}{\hbar\omega_{c}}+\frac{2\varepsilon_{R}}{\hbar\omega_{c}}\right)^{2}+\left(\frac{1}{2}-\frac{\Delta}{2\hbar\omega_{c}}\right)^{2}-\left(\frac{\varepsilon}{\hbar\omega_{c}}\right)^{2}}},$
or (110) $\displaystyle F_{s}^{\prime}(\varepsilon)$ $\displaystyle=$
$\displaystyle\frac{1}{\hbar\omega_{c}}+s\frac{\frac{2\varepsilon_{R}}{\hbar\omega_{c}}}{\sqrt{\frac{4\varepsilon\varepsilon_{R}}{\hbar^{2}\omega_{c}^{2}}+\left(\frac{2\varepsilon_{R}}{\hbar\omega_{c}}\right)^{2}+\left(\frac{1}{2}-\frac{\Delta}{2\hbar\omega_{c}}\right)^{2}}}.$
(111)

As mentioned earlier, the leading term in $F_{s}^{\prime}(\varepsilon)$ is
$1/\hbar\omega_{c}$. By expanding the above expression, we can easily find
$\mathcal{O}\left(\varepsilon_{R}/\varepsilon_{F}\right)=\mathcal{O}[(\alpha
m^{*}\ell_{c})^{2}/\hbar]$ corrections. The above calculation also holds for
the Dresselhaus case. The general case with simultaneous and arbitrary Rashba
and Dresselhaus couplings lead to the corrections $\mathcal{O}[(\alpha
m^{*}\ell_{c})^{2}/\hbar]+\mathcal{O}[(\beta m^{*}\ell_{c})^{2}/\hbar)]$
mentioned following Eq. (12).

##  Appendix B Orthogonal subspaces $\mathcal{P}$

When both Rashba and Dresselhaus are present neither ${\cal N}_{+}$ nor ${\cal
N}_{-}$ are conserved, i.e. $[{\cal N}_{\pm},{\cal\tilde{H}}]\neq 0$. This
will result in mixing of states, e.g. the pure Rashba states will get couple
to each other when a finite $\beta$ is introduced, and vice versa. However,
there is a conserved quantity that can be constructed from ${\cal N}_{\pm}$ by
defining Casanova  _et al._ (2010); Braak (2011)

$\displaystyle\mathcal{P}_{\pm}=\exp(i\pi(\mathcal{N}_{\pm}+\nicefrac{{1}}{{2}})).$
(112)

Using the definition of ${\cal N}_{+}=a^{\dagger}a+\frac{1}{2}\sigma_{z}$ we
can show that

$\displaystyle\mathcal{P}_{+}$ $\displaystyle=$
$\displaystyle\exp\left(i\pi(a^{\dagger}a+\frac{1}{2}\sigma_{z})\right)=\exp\left(i\pi(a^{\dagger}a-\frac{1}{2}\sigma_{z}+\sigma_{z})\right)$
(113) $\displaystyle=$
$\displaystyle\mathcal{P}_{-}\exp(i\pi\sigma_{z})=-\mathcal{P}_{-},$

where we used $\exp(i\pi\sigma_{z})=-1$. Since $\mathcal{P}_{\pm}$ have
eigenvalue $\pm 1$, we only need to consider
$\mathcal{P}=\mathcal{P}_{+}=-\mathcal{P}_{-}$. First, we look at how the
operator $\mathcal{P}$ affects the operators $a$, and $\sigma_{+}$:

$\displaystyle\mathcal{P}\sigma_{+}\mathcal{P}^{\dagger}$ $\displaystyle=$
$\displaystyle
e^{i\frac{\pi}{2}\sigma_{z}}\sigma_{+}e^{-i\frac{\pi}{2}\sigma_{z}}=e^{i\pi}\sigma_{+}=-\sigma_{+}$
(114) $\displaystyle\mathcal{P}a\mathcal{P}^{\dagger}$ $\displaystyle=$
$\displaystyle e^{i\pi a^{\dagger}a}ae^{-i\pi a^{\dagger}a}=e^{i\pi}a=-a$
(115)

The Hamiltonians in both Eqs. (9) and Eq. (70) contain diagonal terms
($a^{\dagger}a$ and $\sigma_{z}$) that commute with $\mathcal{P}$, and non-
diagonal terms that involve odd power $a,a^{\dagger}$ multiplying
$\sigma_{+},\sigma_{-}$, so then its straightforward to show that
$[\mathcal{H},{\cal P}]=0$. Note that $\mathcal{P}$ is unitary so the
condition $[\mathcal{H},{\cal P}]=0$, can be rewritten as
$\mathcal{P}\mathcal{H}\mathcal{P}^{\dagger}=\mathcal{H}$. Focusing on the
spin-orbit part of Eq. (9) one obtains

$\displaystyle\mathcal{P}\left(\alpha_{B}a^{\dagger}\sigma_{-}+\beta_{B}a^{\dagger}\sigma_{+}\right)\mathcal{P}^{\dagger}+\mathrm{h.c.}$
(116) $\displaystyle=$ $\displaystyle\left(\alpha_{B}{\cal
P}a^{\dagger}\mathcal{P}^{\dagger}\mathcal{P}\sigma_{-}\mathcal{P}^{\dagger}+\beta_{B}\mathcal{P}a^{\dagger}\mathcal{P}^{\dagger}\mathcal{P}\sigma_{+}\mathcal{P}^{\dagger}\right)+\mathrm{h.c.}$
$\displaystyle=$
$\displaystyle\left(\alpha_{B}(-a^{\dagger})(-\sigma_{-})+\beta_{B}(-a^{\dagger})(-\sigma_{+})\right)+\mathrm{h.c.}$
$\displaystyle=$
$\displaystyle\left(\alpha_{B}a^{\dagger}\sigma_{-}+\beta_{B}a^{\dagger}\sigma_{+}\right)+\mathrm{h.c.},$

which shows that $\mathcal{P}\mathcal{H}\mathcal{P}^{\dagger}=\mathcal{H}$,
since the diagonal terms in $\mathcal{H}$ trivially commute with
$\mathcal{P}$.

The practical results of having a diagonal operator $\mathcal{P}$ that
commutes with $\mathcal{H}$ is that the Hamiltonian can be diagonalized using
two separate sets of basis states:

$\displaystyle\mathcal{P}=+1$ $\displaystyle:$
$\displaystyle\quad\\{|0,\uparrow\rangle,|1,\downarrow\rangle,|2,\uparrow\rangle,|3,\downarrow\rangle,|4,\uparrow\rangle,\dots\\}$
$\displaystyle\mathcal{P}=-1$ $\displaystyle:$
$\displaystyle\quad\\{|0,\downarrow\rangle,|1,\uparrow\rangle,|2,\downarrow\rangle,|3,\uparrow\rangle,|4,\downarrow\rangle,\dots\\}$

Diagonalizing $\mathcal{H}$ in either of the $\mathcal{P}=+1$, or $-1$,
subspaces will result in a set of states that all anticross. We can connect
these sets of states to $\mathcal{N}_{+}$ eigenstates

$\displaystyle\mathcal{P}=+1$ $\displaystyle:$
$\displaystyle\quad\\{\overbrace{|0,\uparrow\rangle,|1,\downarrow\rangle}^{|0,+\rangle,|1,-\rangle},\overbrace{|2,\uparrow\rangle,|3,\downarrow\rangle}^{|2,+\rangle,|3,-\rangle},|4,\uparrow\rangle,\dots\\}$
$\displaystyle\mathcal{P}=-1$ $\displaystyle:$
$\displaystyle\quad\\{|0,\downarrow\rangle,\underbrace{|1,\uparrow\rangle,|2,\downarrow\rangle}_{|1,+\rangle,|2,-\rangle},\underbrace{|3,\uparrow\rangle,|4,\downarrow\rangle}_{|3,+\rangle,|4,-\rangle},\dots\\},$

and similarly for the $\mathcal{N}_{-}$ eigenstates

$\displaystyle\mathcal{P}=+1$ $\displaystyle:$
$\displaystyle\quad\\{|0,\uparrow\rangle,\overbrace{|1,\downarrow\rangle,|2,\uparrow\rangle}^{|1,-\rangle,|2,+\rangle},\overbrace{|3,\downarrow\rangle,|4,\uparrow\rangle}^{|3,-\rangle,|4,+\rangle},\dots\\}$
$\displaystyle\mathcal{P}=-1$ $\displaystyle:$
$\displaystyle\quad\\{\underbrace{|0,\downarrow\rangle,|1,\uparrow\rangle}_{|0,-\rangle,|1,+\rangle},\underbrace{|2,\downarrow\rangle,|3,\uparrow\rangle}_{|2,-\rangle,|3,+\rangle},|4,\downarrow\rangle,\dots\\}$

Note that $\mathcal{P}$ also commutes with the cubic Dresselhaus terms as is
discussed in App. C.

##  Appendix C Cubic Dresselhaus

The Hamiltonian in Eq. (9) describes a 2DEG with both Rashba and linear
Dresselhaus. For the numerical part we also include the full cubic Dresselhaus
contribution. Starting from Eq. (6.1) in Ref. Winkler, 2003, and projecting
down to the lowest transverse level results in

$\displaystyle\mathcal{H}_{D3}$ $\displaystyle=$
$\displaystyle\frac{(-\gamma\langle
k_{z}^{2}\rangle)}{\hbar}\left(\left[\frac{1}{2}\Pi_{+}\sigma_{+}\right.\right.$
(117) $\displaystyle\left.\left.-\frac{1}{8\hbar^{2}\langle
k_{z}^{2}\rangle}\\{\Pi_{+}^{2}-\Pi_{-}^{2},\Pi_{-}\\}\right]+\mathrm{h.c.}\right),$

where $\Pi_{\pm}=\Pi_{x}\pm i\Pi_{y}$, and
$\langle\Pi_{z}^{2}\rangle=\hbar^{2}\langle k_{z}^{2}\rangle$. Note that now
the Dresselhaus spin-orbit coupling is parametrized by two parameters $\gamma$
and $\langle k_{z}^{2}\rangle$, while for the linear approximation, only the
single parameter $\beta=(-\gamma)\langle k_{z}^{2}\rangle$ is required. Using
the definition in Eqs. (3) and (4) the full Dresselhaus Hamiltonian becomes

$\displaystyle\mathcal{H}_{D3}$ $\displaystyle=$
$\displaystyle\frac{(-\gamma\langle
k_{z}^{2}\rangle)}{\hbar}\left\\{\left[\left(1-\frac{1}{2\langle
k_{z}^{2}\rangle\ell_{c}^{2}}a^{\dagger}a\right)a^{\dagger}\sigma_{+}\right.\right.$
(118) $\displaystyle\left.\left.+\frac{1}{2\langle
k_{z}^{2}\rangle\ell_{c}^{2}}a^{3}\sigma_{+}\right]+\mathrm{h.c.}\right\\}.$

In the absence of spin-orbit interaction $a^{\dagger}a$ can be replaced by its
eigenvalue $n$, which in turn is related to the ratio of the Fermi energy and
$\hbar\omega_{c}$ (valid for $\varepsilon_{F}\ll\hbar\omega_{c}$)

$\displaystyle\frac{1}{\ell_{c}^{2}}a^{\dagger}a\rightarrow\frac{1}{\ell_{c}^{2}}n$
$\displaystyle\approx$
$\displaystyle\frac{1}{\ell_{c}^{2}}\frac{\varepsilon_{F}}{\hbar\omega_{c}}=\frac{k_{F}^{2}}{2}=\pi
n_{2D}.$ (119)

In the presence of spin-orbit we can still formally rewrite Eq. (118) as

$\displaystyle\mathcal{H}_{D3}$ $\displaystyle=$
$\displaystyle\frac{(-\gamma)\left(\langle
k_{z}^{2}\rangle-\frac{\pi}{2}n_{2D}\right)}{\hbar}\left\\{\left[\frac{\langle
k_{z}^{2}\rangle-\frac{1}{2\ell_{c}^{2}}a^{\dagger}a}{\langle
k_{z}^{2}\rangle-\frac{\pi}{2}n_{2D}}a^{\dagger}\sigma_{+}\right.\right.$
(120) $\displaystyle\left.\left.+\frac{1}{2\ell_{c}^{2}}\frac{1}{\langle
k_{z}^{2}\rangle-\frac{\pi}{2}n_{2D}}a^{3}\sigma_{+}\right]+\mathrm{h.c.}\right\\}.$

The prefactor $-\gamma\left(\langle
k_{z}^{2}\rangle-\frac{\pi}{2}n_{2D}\right)$ is defined as

$\displaystyle\beta$ $\displaystyle=$ $\displaystyle\beta_{1}-\beta_{3}$ (121)
$\displaystyle=$ $\displaystyle\left[(-\gamma)\langle
k_{z}^{2}\rangle\right]-\left[(-\gamma)\frac{\pi}{2}n_{2D}\right],$

which reduces to the traditional definition of $\beta$ for low density samples
as considered in Sec. V.

The parity operator $\mathcal{P}$ introduced in App. B also commutes with the
Hamiltonian in Eq. (118), since the spin-orbit terms involve odd powers of
$a,a^{\dagger}$ multiplied by either $\sigma_{+}$ or $\sigma_{-}$, and the
sign introduced the unitary transformation gets cancelled.

##  Appendix D The numerical procedure for finding the $F$-function

For fixed parameter values, the eigenenergies of the Hamiltonian Eq. (9) take
discrete values. They are obtained numerically by diagonalizing the
Hamiltonian matrix using a large enough set of basis states. Finding the
$F$-function as described in Eq. (11) is equivalent to a root finding problem
for the function

$g_{s}(n)=\varepsilon_{n,s}(B)-\varepsilon_{F}=0.$ (122)

This requires the quantum number $n$ to be a continuous variable. which leads
to a minor modification of the Hamiltonian diagonlization procedure. The
standard diagonalization proceedure is to construct a $2N_{L}$ matrix from
$N_{L}$ harmonic oscillator eigenstates, in addition to the spin degree of
freedom. The Pauli matrices form $2\times 2$ blocks that are coupled by the
ladder operators $a$ and $a^{\dagger}$, leading to block tri-diagonal matrix
with $2\times 2$ block matrices

$\displaystyle h_{l,l}$ $\displaystyle=$
$\displaystyle(l-1)\left[\begin{array}[]{cc}1&0\\\
0&1\end{array}\right]+\left[\begin{array}[]{cc}\frac{1-\tilde{\Delta}}{2}&0\\\
0&\frac{1+\tilde{\Delta}}{2}\end{array}\right]$ (127) $\displaystyle
h_{l,l+1}$ $\displaystyle=$
$\displaystyle\sqrt{l+1}\left[\begin{array}[]{cc}0&2\alpha_{\beta}\\\
2\beta_{B}&0\end{array}\right],$ (130)

where $l$ runs from 1 to $N_{L}$ (number of Landau levels used in the
calculations). To obtain a continuous version of Eqs. (127) and (130) a
variable $x$ is added to the index $l$, resulting in

$\displaystyle h_{l,l}(x)$ $\displaystyle=$
$\displaystyle(l+x-1)\left[\begin{array}[]{cc}1&0\\\
0&1\end{array}\right]+\left[\begin{array}[]{cc}\frac{1-\tilde{\Delta}}{2}&0\\\
0&\frac{1+\tilde{\Delta}}{2}\end{array}\right]$ (135) $\displaystyle
h_{l,l+1}(x)$ $\displaystyle=$
$\displaystyle\sqrt{l+x+1}\left[\begin{array}[]{cc}0&2\alpha_{\beta}\\\
2\beta_{B}&0\end{array}\right].$ (138)

The full block-tridiagonal matrix based on the submatrices in Eqs. (135) and
(138) will then yield a spectrum $\varepsilon_{n+x,s}$, for $x\in[-1,1]$. To
further simply the calculations the basis states can be split into
$\mathcal{P}=\pm 1$ subspaces. Each $\mathcal{P}$-subspace contains ordered
states $\\{\epsilon_{0+x},\epsilon_{1+x},\dots\\}$. For each subspace, one
chooses the two adjecent eigenenergies determined by the condition
$\epsilon_{n+x}<\frac{\varepsilon_{F}}{\hbar\omega_{c}}<\epsilon_{n+x+1}$.
Subsequently the value of $x$ is found by solving $g_{s}(n+x)=0$.

##  Appendix E Perturbation theory and “Bogoliubov-de Gennes Hamiltonian”