# Theory for Cd3As2 thin films in the presence of magnetic fields

M. Smith  Materials Science Division, Argonne National Laboratory, Lemont,
Illinois 60439, USA  V. L. Quito  Ames National Laboratory, Ames, Iowa 50011,
USA  Department of Physics and Astronomy, Iowa State University, Ames, Iowa
50011, USA  São Carlos Institute of Physics, University of São Paulo, PO Box
369, 13560-970, São Carlos, SP, Brazil.  A.A. Burkov  Department of Physics
and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada
Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada
P. P. Orth  Ames National Laboratory, Ames, Iowa 50011, USA  Department of
Physics and Astronomy, Iowa State University, Ames, Iowa 50011, USA
Department of Physics, Saarland University, 66123 Saarbrücken, Germany  I.
Martin  Materials Science Division, Argonne National Laboratory, Lemont,
Illinois 60439, USA

###### Abstract

We present a theory for thin films of the Dirac semimetal Cd3As2 in the
presence of magnetic fields. We show that, above a critical thickness,
specific subbands $n$ of thin film Cd3As2 are in a quantum spin Hall insulator
regime and study their response to in- and out-of-plane magnetic fields. We
find that sufficiently large in-plane Zeeman fields drive the system toward a
2D Dirac semimetal regime, provided the field is directed perpendicular to a
high-symmetry mirror plane. For other directions, we find the Dirac points to
be weakly gapped. We further investigate how the system responds to finite
out-of-plane field components, both starting from the quantum spin Hall regime
at small in-plane fields and from the 2D Dirac semimetal regimes at larger in-
plane fields, addressing recent experimental observations in [A. C. Lygo _et
al._ , Phys. Rev. Lett. 130 046201 (2023)] and [B. Guo _et al._ , Phys. Rev.
Lett. 131, 046601(2023)].

##  I Introduction

The materials search for quantum spin Hall (QSH) insulators has been a major
topic in condensed matter physics ever since their theoretical proposal [1, 2,
3, 4, 5, 6, 7, 8]. Due to their intriguing and robust transport properties,
which arise from topologically protected helical edge states and can be
controlled by electric and magnetic fields, these materials are interesting
both from a fundamental perspective as well as for technological applications,
e.g., for a topological field-effect transistor [9]. The QSH effect has thus
far been realized in various systems such as semiconductor quantum wells [10,
11, 12], heavy-element analogs of graphene such as bismuthene [13] and
germanene [14], and transition metal dichalcogenide monolayers [9, 15, 16].
Another promising route to realizing QSH states is to start from three-
dimensional (3D) bulk topological insulators (TIs) or semimetals (TSMs) and
attempt to reach the QSH regime by tuning quantum confinement in thin films or
heterostructures [17, 18, 19, 20]. In TIs, the very thin film limit is
characterized by a hybridization between the top and the bottom 2D surface
states, creating a gapped low energy spectrum [2, 21, 22, 23, 24]. In
contrast, in a TSM with surface normal directed along the axis that separates
the bulk nodes, there exists no top and bottom surface states to hybridize,
and instead, the bulk spectrum is gapped due to quantum confinement [19, 18].
This situation has recently been experimentally realized in thin films of
Cd3As2 [25, 26], which is a prototypical 3D Dirac TSM in the bulk [27]. It was
driven to the thin film limit in epitaxial $\approx 20$ nm thick films that
were grown along the $[001]$ crystallographic direction, which is parallel to
the separation direction between the bulk Dirac nodes. Various transport
properties were measured in the presence of magnetic fields [25, 26];
explaining these is a primary motivation for our work. Thus, here we
theoretically consider Cd3As2 thin films that terminate along the [001]
crystallographic direction and investigate their band structure and edge
states in the presence of a magnetic field. We characterize their dependence
on film thickness, chemical potential, and magnetic field strength and
direction, a crucial step to understanding Cd3As2 thin film properties.

In thin films of Cd3As2 the transverse momentum quantization leads to a
reorganization of the 3D Dirac dispersion into 2D subbands. The subbands are
generally gapped, although the gaps can be small, making the effects of
external magnetic fields potentially large. Recent magnetotransport
measurements in a perpendicular magnetic field along the $[001]$ direction
[25] revealed a magnetic field dependence of the Landau level spectrum that is
in agreement with that of a 2D TI [10]. Subsequent experiments in tilted
magnetic fields with both in- and out-of-plane components showed an odd-
integer quantum Hall effect [26]. This was attributed to a transition to a 2D
Dirac semimetal phase, with the two Dirac cones appearing along momentum space
directions perpendicular to the in-plane magnetic field.

Motivated by these recent experimental findings, we develop and analyze a
model for the low-energy electronic states in $[001]$ thin films in the
presence of magnetic fields pointing along general directions. Specifically,
we construct a low-energy $k\cdot p$ band structure model that considers the
contributions of the different types of electrons close to the Fermi level. We
demonstrate that the $s$ and relevant $p$ orbitals experience in-plane
magnetic fields in a very distinct way: due to the strong spin-orbit coupling,
the relevant $p$ orbitals show a cubic Zeeman effect, while the $s$ orbitals
experience the standard linear Zeeman coupling. This distinct coupling
manifestly conserves the $C_{4z}$ symmetry of the system and differentiates
our model from recent work in Ref. [28], which considers a linear coupling of
magnetic field to the relevant $p$ orbitals. We show below that the two models
have distinct behaviors as a function of in-plane field orientation.
Nevertheless, both models similarly yield a Dirac semimetal phase for fields
pointed in the $[100]$ direction, which Ref. [28] refers to as a 2D Weyl
semimetal. Ref. [28] does not consider the effect of rotations of the in-plane
field and focuses on non-trivially gapping out the 2D Dirac nodes via the
application of lattice strain. Here, we show that the cubic $p$-orbital Zeeman
effect for in-plane magnetic fields leads to a 2D Dirac semimetal phase only
when the field is applied in a direction perpendicular to a mirror plane and
naturally leads to a gapped spectrum for fields rotated away from these high
symmetry directions.

Specifically, we derive the proper coupling of the relevant $p$ orbitals to
the magnetic field, showing it is consistent with all symmetry operations,
which are explicitly given in Appendix B. We then investigate the band
structure as a function of in-plane magnetic field. For specific field
directions (perpendicular to a mirror plane), we find that there exists a
critical field strength such that the system transitions from a gapped to a
gapless band structure with two Dirac points appearing along the line
perpendicular to the applied field (_i.e._ , along the mirror line). As the
field is rotated away from these high symmetry directions, the Dirac points
become gapped. Interestingly, applying an out-of-plane Zeeman field can close
and reopen these Dirac gaps, leading to a change of the Chern number in the
system. We further use this model to explore the orbital effect of the
perpendicular component of the magnetic field when the system is driven into
the 2D Dirac semimetal regime by an in-plane field. We find there is an
emergent quantum number for the low-lying Landau levels, the Dirac states,
which causes an additional two-fold Landau level degeneracy. These findings
provide a plausible explanation for the recent quantum Hall measurements
reported in Ref. [26]. Lastly, we discuss the evolution of the surface states
in the presence of an in-plane Zeeman field, showing that the conventional
edge states are coupled by the in-plane magnetic field and thus gapped and
eventually disappear for large enough in-plane field strength.

The remainder of this paper is structured as follows: Section II establishes
the effective low-energy theory for Cd3As2 thin films with surface normal
along the [001] direction. The resulting model is generally applicable to
Dirac TSMs with a surface normal that is parallel to the direction separating
the Dirac points. Section III discusses the electronic response to in-plane
Zeeman fields, first only considering the Zeeman coupling of the $s-$orbitals
and then deriving the smaller (nonlinear) Zeeman effect in the relevant
$p-$orbitals. Section IV discusses the surface states of the thin film with
and without in-plane Zeeman coupling. Section V considers the effect of an
out-of-plane orbital magnetic field. We first discuss a purely out-of-plane
field and then include a possibly large in-plane field component. Finally, in
Section VI, we summarize our findings and predictions.

##  II Model for (001) thin films

Here, we derive a low-energy model of Cd3As2 thin films, grown such that the
surface normal points along the $[001]$ direction. We start from a low-energy
four-band model for bulk Cd3As2. The low-energy bands are located around the
$\Gamma$ point in the Brillouin zone and include the $\ket{S_{1/2},m_{J}=\pm
1/2}$ and $\ket{P_{3/2},m_{J}=\pm\ 3/2}$ states [27]. The relevant point group
constraining the model is $D_{4h}$ (or $4/mmm$), and up to cubic order in
momentum, it reads

$\displaystyle H_{0}$
$\displaystyle(\mathbf{k})=\epsilon_{0}(\mathbf{k})\tau_{0}s_{0}+\mathcal{M}(\mathbf{k})\tau_{z}s_{0}$
$\displaystyle+A\left(k_{x}\tau_{x}s_{z}-k_{y}\tau_{y}s_{0}\right)+C_{3a}k_{x}k_{y}\left(k_{x}\tau_{y}s_{0}-k_{y}\tau_{x}s_{z}\right)$
$\displaystyle+C_{3b}\left(k_{x}^{3}\tau_{x}s_{z}-k_{y}^{3}\tau_{y}s_{0}\right)+C_{3c}k_{z}^{2}\left(k_{y}\tau_{y}s_{0}-k_{x}\tau_{x}s_{z}\right)$
$\displaystyle+C_{3d}(k_{x}^{2}-k_{y}^{2})k_{z}\tau_{x}s_{x}-C_{3e}k_{x}k_{y}k_{z}\tau_{x}s_{y}\,.$
(1)

Here, $s_{i}$ and $\tau_{i}$ are Pauli matrices in angular momentum
($m_{J}=\pm$) and orbital ($S,P$) space, respectively, and we have defined

$\displaystyle\epsilon_{0}(\mathbf{k})$
$\displaystyle=C_{0}+C_{1}k_{z}^{2}+C_{2}k_{\parallel}^{2}$ (2a)
$\displaystyle\mathcal{M}(\mathbf{k})$
$\displaystyle=M_{0}-M_{z}k_{z}^{2}-M_{xy}k_{\parallel}^{2}$ (2b)

with $k_{\parallel}^{2}=k_{x}^{2}+k_{y}^{2}$ and $M_{xy},M_{z}>0$, consistent
with [29]. The $k_{z}$ direction points along the fourfold rotation axis,
corresponding to the $[001]$ direction. A detailed derivation of this
effective model is provided in Appendix A. Up to quadratic order, it agrees
with the Hamiltonian derived in Refs. [27, 29, 25]. The four-dimensional basis
is obtained by projecting the two $s$ and six $p$ bands onto the total angular
momentum basis $J$ and keeping the low energy states $\ket{S_{1/2},m_{J}=\pm
1/2}$ and $\ket{P_{3/2},m_{J}=\pm\ 3/2}$. The four other $p$-orbital states
lie approximately $250$ meV higher in energy [27]. The low-energy model for a
film of thickness $L$ is constructed by imposing open boundary conditions
along the $[001]$ direction. The Hamiltonian at zero in-plane momentum,
$\mathbf{k}_{\parallel}=0$, reads

$\displaystyle
H_{0}|_{\mathbf{k}_{\parallel}=0}=(C_{0}-C_{1}\partial_{z}^{2})+(M_{0}+M_{z}\partial_{z}^{2})s_{0}\tau_{z}$
(3)

Let us focus on the $m_{J}=+$ subspace spanned only by the $\ket{S_{1/2},1/2}$
and $\ket{P_{3/2},3/2}$ states as the energy levels of the opposite angular
momentum, $m_{J}=-$, are degenerate. Since $H_{\mathbf{k}_{\parallel}=0}$ is
diagonal, the orbital structure of the eigenstates read $\psi_{+}=(1,0)^{T}$
and $\psi_{-}=(0,1)^{T}$, _i.e._ they localized on the $S$ and $P$ states,
respectively. Using the trial wavefunctions $\Psi_{\beta
n}(z)=A\psi_{\beta}e^{\lambda_{\beta n}z}$, and imposing boundary conditions
such that the wavefunctions vanish at $z=\pm L/2$, we find the wavefunction
for $l=0,1,2,\ldots$ and $\beta=\pm$ is given by

$\displaystyle\Psi_{\beta n}(z)$
$\displaystyle=\sqrt{\frac{2}{L}}\psi_{\beta}\begin{cases}\sin(2l\pi\frac{z}{L}),\quad
n=2l\\\ \cos((2l+1)\pi\frac{z}{L}),\quad n=2l+1\,.\end{cases}$ (4)

The associated energies at $\mathbf{k}_{\parallel}=0$ read

$\displaystyle E_{\beta n}$ $\displaystyle=C_{0}+\beta M_{0}+\left(\frac{\pi
n}{L}\right)^{2}\left(C_{1}-\beta M_{z}\right).$ (5)

Here $n=1,2,3,\ldots$ is the sub-band index, and $\beta=\pm$ denotes the
orbital index. Thus, the subbands for $[001]$ thin films are found by setting
$k_{z}=\pi n/L$. The Hamiltonian for a specific subband takes the form of Eq.
(1) with

$\displaystyle\epsilon_{0}\left(\mathbf{k}\right)$
$\displaystyle\rightarrow\epsilon_{n}\left(\mathbf{k}_{\parallel}\right)=C_{0}+C_{1}\left(\frac{\pi
n}{L}\right)^{2}+C_{2}k_{\parallel}^{2}$ (6a)
$\displaystyle\mathcal{M}\left(\mathbf{k}\right)$
$\displaystyle\rightarrow\mathcal{M}_{n}\left(\mathbf{k}_{\parallel}\right)=M_{0}-M_{z}\left(\frac{\pi
n}{L}\right)^{2}-M_{xy}k_{\parallel}^{2}\,.$ (6b)

The total Hamiltonian thus involves 4 different bands for each subband $n$. In
the following, we use the notation $\mathbf{k}\equiv\mathbf{k}_{\parallel}$
unless specified otherwise. Importantly, the topological character of the
bands is determined by

$\mathcal{M}_{n}(0)\equiv M_{n}=M_{0}-M_{z}\left(\frac{n\pi}{L}\right)^{2}\,.$
(7)

As $M_{xy}>0$, for positive $M_{n}>0$ the bands are inverted and topologically
non-trivial, while for negative $M_{n}<0$, the bands un-invert and become
topologically trivial. Considering the film thicknesses reported in Refs. [25,
26], we find that for a thickness $L\sim 19$ nm it is the $n=2$ subband which
is most relevant. It exhibits a topological gap $M_{n=2}\sim 5.5$ meV and the
$n=1,3$ subbands are $\sim 20$ meV away. The relatively small topological gap
of the $n=2$ subband and the large effective $g$-factor of the Zeeman coupling
make thin film Cd3As2 easily tunable by in-plane Zeeman fields. As will be
shown below, the addition of an in-plane Zeeman field causes the topological
gap for the $m_{J}=\pm$ states to be different, leading to a gapless regime
for specific Zeeman field directions.

##  III In-plane Magnetic Field Effects

In this Section, we show the effects of adding in-plane magnetic fields
$\mathbf{B}=B(\cos\phi,\sin\phi,0)$. We start with a derivation of the Zeeman
coupling to a magnetic field in the $xy$ plane.

###  III.1 Microscopic derivation of the Hamiltonian

The magnetic field $\mathbf{B}$ couples to the $s$ and $p$ states according to
the microscopic Zeeman Hamiltonian

$H_{Z}=g_{s}\mu_{B}\mathbf{B}\cdot\left(\mathbf{L}+2\mathbf{S}\right)+g_{p}\mu_{B}\mathbf{B}\cdot\left(\mathbf{L}+2\mathbf{S}\right).$
(8)

Here, $\mathbf{L}$ and $\mathbf{S}$ are the orbital and spin angular momentum
operators, $\mu_{B}$ is the Bohr magneton, and $g_{s},g_{p}$ are the
$g$-factors associated with $s$ and $p$ states. Given the large spin-orbit
coupling in the system, it is convenient to work with the total angular
momentum $J$. The states close to the Fermi surface are the doublets of
$J=3/2$, $\ket{P_{3/2},\pm\frac{3}{2}}$, and of $J=1/2$, $\ket{S_{1/2},\pm
1/2}$ [27]. The $s$ electrons of the Kramers pair $\ket{S_{1/2},\pm 1/2}$
couple to the in-plane magnetic field as
$H_{Z}^{(S)}=g_{s}\mu_{B}\mathbf{B}\cdot\mathbf{s}$, where
$\mathbf{s}=(s_{x},s_{y},s_{z})$ is the vector of the Pauli matrices. The low
energy $p$ states $\ket{P_{3/2},\pm 3/2}$, however, do not couple linearly to
the in-plane component of the field. Instead, the in-plane magnetic field
couples the $\ket{P_{3/2},\pm\frac{3}{2}}$ only in the cubic order. This can
be understood as follows: terms in the Hamiltonian linear in the applied in-
plane field can only change the azimuthal quantum number $m_{J}$ by $1$, and
therefore, the components $m_{J}=\pm 3/2$ are not coupled to linear order in
the field; cubic coupling is, however, possible, involving transition via
$\ket{P_{1/2},\pm\frac{1}{2}}$ and $\ket{P_{3/2},\pm\frac{1}{2}}$ virtual
states. The explicit calculations deriving this coupling are performed in
Appendix B.

To lowest order in an in-plane magnetic field
$\mathbf{B}\equiv\mathbf{B}_{\parallel}=(B_{x},B_{y},0)$, the effective
coupling of the $\ket{P_{3/2},\pm\frac{3}{2}}$ states reads

$\displaystyle
H_{Z}^{(P)}(\mathbf{B}_{\parallel})=\frac{g^{3}_{p}\mu_{B}^{3}}{\Delta^{2}}$
$\displaystyle\Bigl{[}B_{x}\bigl{(}B_{x}^{2}-3B_{y}^{2}\bigr{)}s_{x}$
$\displaystyle\quad-B_{y}\bigl{(}B_{y}^{2}-3B_{x}^{2}\bigr{)}s_{y}\Bigr{]}\,.$
(9)

Here, $\Delta$ is the energy scale of the higher bands at the $\Gamma$ point
(see the derivation in Appendix B). Thus, the in-plane Zeeman field terms for
the effective low-energy model are given by

$\displaystyle H($
$\displaystyle\mathbf{B}_{\parallel})=g_{s}\mu_{B}\left(\frac{\tau_{0}+\tau_{z}}{2}\right)\mathbf{B}_{\parallel}\cdot\mathbf{s}+\frac{g_{p}^{3}\mu_{B}^{3}}{\Delta^{2}}\left(\frac{\tau_{0}-\tau_{z}}{2}\right)$
$\displaystyle\times\left[B_{x}\left(B_{x}^{2}-3B_{y}^{2}\right)s_{x}-B_{y}\left(B_{y}^{2}-3B_{x}^{2}\right)s_{y}\right].$
(10)

It is convenient to parameterize the field in the $xy$ plane as
$B_{x}=B_{\parallel}\cos\phi$ and $B_{y}=B_{\parallel}\sin\phi$. Returning to
Eq. (1), we simplify our calculation by choosing coefficients
$C_{3a}=3C_{3b}=3A_{1}$,$C_{3c}=0$, and $C_{3e}=-C_{3d}/2=-2\eta$, which
simplifies Eq. (1) to its atomic limit form (see Appendix D for further
details). Finally, the full Hamiltonian, including the in-plane field, is
found from Eq. (1) combined with Eq (10). Under the assumption that the
finite-$k$ terms do not lead to significant changes in the wavefunctions along
the $z$ direction, we project the Hamiltonian into states of Eq. (4), leading
to the effective Hamiltonian

$\displaystyle H$
$\displaystyle=\epsilon_{n}(\mathbf{k})\tau_{0}s_{0}+\mathcal{M}_{n}(\mathbf{k})\tau_{z}s_{0}+A(k_{x}\tau_{x}s_{z}-k_{y}\tau_{y}s_{0})$
$\displaystyle+3A_{1}k_{x}k_{y}\left(k_{x}\tau_{y}s_{0}-k_{y}\tau_{x}s_{z}\right)+A_{1}\left(k_{x}^{3}\tau_{x}s_{z}-k_{y}^{3}\tau_{y}s_{0}\right)$
$\displaystyle+\eta\frac{n\pi}{L}\left[(k_{x}^{2}-k_{y}^{2})\tau_{x}s_{x}+2k_{x}k_{y}\tau_{x}s_{y}\right]$
$\displaystyle+\frac{\tau_{0}+\tau_{z}}{2}\bm{h}_{s}\cdot\bm{s}+\frac{\tau_{0}-\tau_{z}}{2}\frac{h_{p}^{3}}{\Delta^{2}}\left[\cos(3\phi)s_{x}+\sin(3\phi)s_{y}\right].$
(11)

Here, we labeled $g_{s}\mu_{B}\mathbf{B}_{\parallel}\equiv\bm{h}_{s}$ and and
$g_{p}\mu_{B}\mathbf{B}_{\parallel}\equiv\bm{h}_{p}$.

###  III.2 Effect of $s$-orbital Zeeman Field

It is one of the main goals of this work to show the importance of including
the magnetic field effect on the $p$ orbitals to properly capture features of
the low energy bands. As $h_{p}^{3}/\Delta^{2}\ll h_{s}$, however, we can
first ignore it to find the eigenvalues (up to quadratic order in $k$)

$\displaystyle E_{\alpha\beta n}$
$\displaystyle=\epsilon_{n}(\mathbf{k})+\alpha\frac{h_{s}}{2}+\beta\sqrt{A^{2}k^{2}+\left(\mathcal{M}_{n}(\mathbf{k})+\alpha\frac{h_{s}}{2}\right)^{2}}.$
(12)

Here $\beta=\pm$ denotes the $S$ and $P$ orbital character of each band at
$\Gamma$, respectively, while $\alpha=\pm$ denotes the two states in the
Kramers doublet of a given orbital. There are two distinct regimes that
display different behavior at the $\Gamma$ point, where the energies simplify
to

$\displaystyle E_{\alpha\beta n}(\mathbf{k}=0)$
$\displaystyle=\epsilon_{n}+\alpha\frac{h_{s}}{2}+\beta\left|M_{n}+\alpha\frac{h_{s}}{2}\right|$
(13)

In the weak field regime, when $\left|M_{n}\right|>\frac{h_{s}}{2}$, the
system is gapped at $\Gamma$. Focusing on the interesting case of a
topological gap, $M_{n}>0$, we find

$\displaystyle E_{\alpha\beta n}(\mathbf{k}=0)$
$\displaystyle=\epsilon_{n}+\alpha(1+\beta)\frac{h_{s}}{2}+\beta M_{n}\,.$
(14)

The $s$ bands ($\beta=+$) are spin-split, whereas the $p$ bands ($\beta=-$)
are degenerate at $\Gamma$.

In the strong field regime, when $h_{s}>2M_{n}$, the $\alpha=-$ bands have
been un-inverted as the gap for the $\alpha=-$ s-orbital band at $\Gamma$ is
given by $M_{n}-\frac{h_{s}}{2}<-M_{n}$. The un-inversion of the $\alpha=-$
bands leads to a quadratic band touching between the $p-$orbital bands. As we
will see below, the inclusion of the $p-$orbital Zeeman field resolves the
quadratic band touching into two Dirac points perpendicular to the applied
field, provided the field points perpendicular to a mirror plane. We estimate
that for sample sizes $L=19$ nm the relevant subband $n=2$ has a topological
gap of size $M_{n=2}\approx 5.5$ meV such that the strong field regime is
reached by Zeeman fields $h_{s}\sim 11$ meV. For a $g-$factor of $g_{s}=12$
[30] this corresponds to fields $B_{\parallel}\sim 16$ T.

Figure 1: Band structure plots for thin film Cd3As2 in the low and high field
regimes when setting $h_{p}=0$. The dispersion is given by Eq. (12) with model
coefficients taken from Ref. [29] (also listed at the end of Appendix E) (a)
Low-field dispersion for $h_{s}=M_{n}/2$ with non-degenerate subbands. The
bands which arise from the $s-$orbitals at the $\Gamma$ point are colored red,
and the ones arising from the $p-$orbitals are colored black. (b) High-field
dispersion for $h_{s}=5M_{n}/2$, where the $\alpha=-$ bands have un-inverted,
resulting in a quadratic band touching of the $p-$bands at $\Gamma$ for
$h_{p}=0$.

###  III.3 Inclusion of the $p-$orbital Zeeman field

The inclusion of the $p-$orbital field is only relevant in the large field
limit where the $\alpha=-$ bands have been un-inverted, and the low energy
theory is controlled by the bands with $p$-orbital character near the $\Gamma$
point. Importantly, as $|\mathbf{k}|$ is increased, there is a significant
contribution to these bands from the $s$-orbitals. It is the main goal of this
section to show that, by including the $p$-orbital field and considering the
terms of cubic order in $k$, the low-energy Dirac fermions are gapped out for
specific field orientations.

By integrating out the bands at higher energies (See Appendix D for further
details), we arrive at an effective Hamiltonian for the lowest energy bands,

$\displaystyle\mathcal{H}_{p}=-M_{n}+\left(\begin{array}[]{cc}\mathcal{B}_{2}&\mathcal{B}_{1}\\\
\mathcal{B}_{1}^{*}&-\mathcal{B}_{2}\end{array}\right)$ (17)

where, up to quintic order,

$\displaystyle\mathcal{B}_{1}$
$\displaystyle=\biggl{(}\frac{A^{2}k_{-}^{2}}{h_{s}}+\frac{2AA_{1}k^{2}k_{+}^{2}}{h_{s}}\biggr{)}e^{-i\phi}$
$\displaystyle\qquad-\frac{\eta^{2}k_{-}^{4}}{h_{s}}\left(\frac{n\pi}{L}\right)^{2}e^{i\phi}+\frac{h_{p}^{3}}{\Delta^{2}}e^{-3i\phi}$
(18a) $\displaystyle\mathcal{B}_{2}$ $\displaystyle=-\frac{2\eta
k^{3}}{h_{s}}\left(\frac{n\pi}{L}\right)\left[A\cos(\theta-\phi)+A_{1}k^{2}\cos(5\theta-\phi)\right]$
(18b)

Here, $\theta$ and $\phi$ are defined via
$\mathbf{k}=k(\cos\theta,\sin\theta,0)$ and
$\mathbf{B}=B(\cos\phi,\sin\phi,0)$, respectively, and we have defined
$k_{-}=k_{x}-ik_{y}$. We first consider the eigenvalues up to
$\mathcal{O}(k^{2})$, which are given by

$\displaystyle E_{p\pm}$
$\displaystyle=-M_{n}\pm\left|\frac{A^{2}k^{2}}{h_{s}}e^{-2i(\theta-\phi)}+\frac{h_{p}^{3}}{\Delta^{2}}\right|.$
(19)

We can see that, to $\mathcal{O}(k^{2})$, two Dirac points are created
perpendicular to the applied field direction at
$\theta_{c\sigma}=\phi+\sigma\pi/2$ and $k_{c0}=\sqrt{h_{s}h_{p}^{3}}/(\Delta
A)$, where $\sigma=\pm$ denotes the two Dirac points.

The inclusion of terms $\mathcal{O}(k^{3})$ lowers the artificial symmetry
$C_{\infty,z}$ down to $C_{4z}$, and the Dirac fermions become gapped for
field directions which are not perpendicular to a high symmetry plane [22].
Further, $k_{c0}$ is shifted slightly to $k_{c}$, which is the solution of
$|\mathcal{B}_{1}|=0$ for $\theta=\theta_{c\sigma}$. This shift, however, is
small due to the smallness of $h_{p}^{3}/\Delta^{2}$ in comparison to $h_{s}$
and therefore $k_{c}$ is well approximated by $k_{c0}$. (See Appendix D for
more details) We may model the Dirac fermions by expanding $\mathcal{B}_{1,2}$
about $\mathbf{k}_{c}=k_{c}(\cos\theta_{c\sigma},\sin\theta_{c\sigma})$,
keeping the lowest order non-zero contributions in $k_{c}$, to get

$\displaystyle\mathcal{H}_{\sigma}$
$\displaystyle=-M_{n}+\left(\begin{array}[]{cc}\sigma
m(\phi)&\frac{-2iA^{2}k_{c}}{h_{s}}\delta k_{-}e^{-2i\phi}\\\
\frac{2iA^{2}k_{c}}{h_{s}}\delta k_{+}e^{2i\phi}&-\sigma
m(\phi)\end{array}\right)\,,$ (22)

where $m(\phi)=m_{0}\sin(4\phi)$, and $m_{0}=\frac{2\eta
A_{1}k^{5}_{c}}{h_{s}}\left(\frac{n\pi}{L}\right)$ is the mass of the Dirac
fermion. We also define $\delta k_{\pm}=\delta k_{x}\pm i\delta k_{y}$ with
$\delta k_{x,y}$ being the deviation of the momentum away from $\bm{k}_{c}$.
Importantly, the mass $m(\phi)$ is zero when the field is applied
perpendicular to a mirror plane of the relevant point group $D_{4h}$, as
derived in Ref. [22]. Since the system obeys $D_{4h}$ point symmetry, this
corresponds to the field being applied along the $x$ or the $y$ axis as well
as along the lines $x=\pm y$. As the field is rotated through a plane where
$m(\phi)$ vanishes, the mass changes sign, switching the chirality of the two
flavors of Dirac fermions ($\sigma=\pm$), but keeping the total Chern number
of the system unchanged.

Figure 2: 3D Plot of the bands from Eq. (17) when $\phi=0$, showing the two
band touchings which give rise to the emergent Dirac states. We use
$h_{s}=2.5M_{n}$, $h_{p}^{3}/\Delta^{2}=h_{s}/20$, $A_{1}=\eta=100$
meV$\cdot$nm3, and other model parameters are taken from Ref. [29].

###  III.4 Anomalous Hall Effect for tilted field

We now discuss the effect of tilting the Zeeman field out of the $xy$-plane
such that it acquires a finite $B_{z}$ component. As the basis states are
$\alpha=\pm$ states an out-of-plane Zeeman field enters in the form
$h_{z}=2g_{p}\mu_{B}B_{z}s_{z}$ (see also Appendix B). The Dirac Hamiltonian
becomes

$\displaystyle\mathcal{H}_{\sigma}$
$\displaystyle=-M_{n}+\left(\begin{array}[]{cc}\sigma
m(\phi)+h_{z}&\frac{-2iA^{2}k_{c}}{h_{s}}\delta k_{-}e^{-2i\phi}\\\
\frac{2iA^{2}k_{c}}{h_{s}}\delta k_{+}e^{2i\phi}&-\sigma
m(\phi)-h_{z}\end{array}\right).$ (25)

For sufficiently large $h_{z}$ such that $|h_{z}|>|m(\phi)|$, the gap of one
Dirac fermion can change sign, leading to a change of the Chern number of the
system by $\pm 1$ and causing an anomalous Hall effect. The Chern number of
the system as a function of $h_{z}$ is shown in Fig. 3.

Figure 3: The Chern number of the system as a function of out-of-plane Zeeman
field $h_{z}$. When the applied field becomes larger than the mass of each
Dirac point, $|h_{z}|>|m|$, the gap of one of the Dirac points changes sign,
leading to a change in the Chern number of the system.

##  IV Destruction of helical Surface States by in-plane field

We now return to the case of purely in-plane fields, neglecting the orbital
magnetic field, and investigate the evolution of the non-trivial helical edge
states in the presence of an in-plane Zeeman field. Without any applied
magnetic field, the effective model for thin film Cd3As2 is that of a QSH
insulator, and we, therefore, expect non-trivial helical edge states to arise
due to the inverted subbands. After the in-plane field is turned on, we show
that the helical edge states become coupled, resulting in a gapping of the
surface states originating for each inverted subband. Because for reasonable
field strengths $h_{p}^{3}/\Delta^{2}\ll h_{s}$, we will neglect the
$p$-orbital Zeeman field in the following.

###  IV.1 Derivation of surface state solutions

We impose open boundary conditions in the $y$ direction and consider a sample
that is a semi-infinite plane for $y\geq 0$. We thus replace
$k_{y}\rightarrow-i\partial_{y}$ in Eq. (11). For simplicity, we only keep
terms up to $\mathcal{O}(k^{2})$ and neglect the un-inverted band structure
$\epsilon_{n}(\mathbf{k})$, as including these terms would not alter our
conclusions. We then use the ansatz

$\Psi_{\alpha\beta
n}(k_{x},\lambda_{n},y)=\psi_{\alpha\beta}(k_{x},\lambda_{n})e^{\lambda_{n}y}\,,$
(26)

where $\psi_{\alpha\beta}(k_{x},\lambda)$ is a four-spinor, and
$\alpha,\beta,n$ denote spin, orbital, and subband indices. Since the sample
is located at $y\geq 0$, we look for solutions with $\lambda_{n}<0$. The
eigenenergies and eigenvectors as a function of $k_{x}$ and $\lambda$ are
found as

$\displaystyle E_{\alpha\beta n}(k_{x},-i\lambda)$
$\displaystyle=\alpha\frac{h_{s}}{2}+\beta\sqrt{A^{2}k^{2}+\Bigl{(}\mathcal{M}_{n}(\bm{k})+\alpha\frac{h_{s}}{2}\Bigr{)}^{2}}$
(27a) $\displaystyle\psi_{\alpha\beta}(k_{x},-i\lambda)$
$\displaystyle=\mathcal{N}_{\alpha\beta}(\bm{k})\begin{pmatrix}-\alpha
e^{-i\phi}\frac{E_{\alpha\beta
n}(\mathbf{k})+\mathcal{M}_{n}(\mathbf{k})}{Ak_{+}}\\\ -\alpha
e^{-i\phi}\frac{k_{-}}{k_{+}}\\\ -\frac{E_{\alpha\beta
n}(\mathbf{k})+\mathcal{M}_{n}(\mathbf{k})}{Ak_{+}}\\\ 1\end{pmatrix}\,.$
(27b)

Here, $k^{2}=k_{x}^{2}-\lambda^{2}$, $k_{\pm}=k_{x}\pm\lambda$, and
$\mathcal{N}_{\alpha\beta}(\mathbf{k})$ is a normalization factor. We can
express $\lambda$ as a function of energy by solving $E_{\alpha\beta
n}(k_{x},-i\lambda)=\epsilon$ and obtain four solutions

$\displaystyle\lambda_{\alpha\zeta}(\epsilon)$
$\displaystyle=-\sqrt{k_{x}^{2}+\frac{\tilde{B}_{\alpha}}{2}+\zeta\sqrt{\tilde{B}_{\alpha}^{2}/4-\tilde{C}_{\alpha}(\epsilon)}}$
(28)

where $\alpha,\zeta=\pm$ and

$\displaystyle\tilde{B}_{\alpha}$
$\displaystyle=\frac{A^{2}-2M_{xy}\left(M_{n}+\alpha\frac{h_{s}}{2}\right)}{M_{xy}^{2}},$
(29a) $\displaystyle\tilde{C}_{\alpha}(\epsilon)$
$\displaystyle=\frac{\left(M_{n}+\alpha\frac{h_{s}}{2}\right)^{2}-\left(\epsilon-\alpha\frac{h_{s}}{2}\right)^{2}}{M_{xy}^{2}}.$
(29b)

The wavefunction is then

$\displaystyle\Psi(k_{x},y)$
$\displaystyle=\sum_{\alpha\zeta\beta}\tilde{\mathcal{N}}_{\alpha\zeta\beta}\psi_{\alpha\beta}[k_{x},-i\lambda_{\alpha\zeta}(\epsilon)]e^{\lambda_{\alpha\zeta}(\epsilon)y}$
(30)

where $\tilde{\mathcal{N}}_{\alpha\zeta\beta}$ are constants and the spinor
$\psi_{\alpha\beta}(k_{x},-i\lambda_{\alpha\zeta})$ is the eigenvector
associated with $E_{\alpha\beta n}$ and given by Eq. (27b). The boundary
condition is $\Psi(k_{x},\lambda_{n},y=0)=0$, _i.e._ the wavefunction vanishes
on the surface. Without loss of generality we set
$\tilde{\mathcal{N}}_{\alpha\zeta\beta}=\tilde{\mathcal{N}}_{\alpha\zeta}\delta_{\beta+}$,
as for a given $\alpha$ the two $\beta=\pm$ eigenstates give the same
solution. As shown in detail in Appendix E, the boundary conditions can be
expressed as a determinant of some matrix $\mathcal{A}$ via
$\det\mathcal{A}=0$, where the columns of $\mathcal{A}$ contain the spinors
$\psi_{\alpha+}$. We then solve the boundary conditions for $\epsilon$ as a
function of $k_{x}$, which are plotted at various in-plane field values to
yield the dispersions shown in Figs. 4 and Fig. 5.

###  IV.2 Zero field results

In the absence of an in-plane Zeeman field, $h_{s}=0$, we observe linearly
dispersing helical edge states, as expected for a QSH insulator [31]. Because
the $\alpha=\pm$ subspaces are degenerate, we may simplify our boundary
conditions for $\alpha=+$ to

$\displaystyle\frac{E_{+\beta
n}(k_{x},-i\lambda_{++})+\mathcal{M}_{n}(k_{x},-i\lambda_{++})}{E_{+\beta
n}(k_{x},-i\lambda_{+-})+\mathcal{M}_{n}(k_{x},-i\lambda_{+-})}=\frac{k_{x}+\lambda_{++}}{k_{x}+\lambda_{+-}}$
(31)

which we solve numerically for the energy $\epsilon$ as a function of $k_{x}$.
Note that the energy $\epsilon$ enters into Eq. (31) quadratically, and the
correct root is chosen such that the $\beta$ surface state connects to the
$\beta$ bulk band. The $\alpha=-$ surface state can then be obtained via the
application of time reversal. As seen in Fig. 4, the $\alpha=+$ surface state
exhibits $v_{x}>0$ and the $\alpha=-$ surface state exhibits $v_{x}<0$, .

Figure 4: Gapless helical surface states on a single Cd3As2 thin film surface
in the absence of an in-plane Zeeman field. The $\alpha=+$ surface state is
colored red, and the $\alpha=-$ surface state is colored blue, with the bulk
bands colored black. Model parameters are taken from Ref. [29]. We use a film
thickness of $L=19$ nm and show the $n=2$ bulk subbands with energies $\pm
M_{n=2}=\pm 5.5$ meV at $\mathbf{k}=0$.

It should be noted that for the surface states to exist one must have

$\displaystyle\frac{A^{2}}{M_{xy}^{2}}>4\frac{M_{n}}{M_{xy}}>0.$ (32)

in order for $\lambda_{\alpha\zeta}$ to have a real component. The
coefficients used for the numerical evaluation of the surface state dispersion
are the same as used in Ref. [29] (See Appendix E for details).

###  IV.3 Gapping of surface states for nonzero in-plane fields

When the Zeeman field is small, $h_{s}<2M_{n}$, the bulk bands are still
inverted, but the helical edge states are now coupled by $h_{s}$. This opens
up a gap in the edge state spectrum on the order of $h_{s}$, which is
consistent with the explicit breaking of time-reversal symmetry. For
$h_{s}<\frac{2}{3}M_{n}$ there remain two (gapped) surface state solutions
(See Fig. 5(a). As the field becomes larger, $h_{s}\geq\frac{2}{3}M_{n}$, the
$S$-orbital $\alpha=-$ band is pushed through the upper surface state which
merges with the bulk. This leaves only a single gapped surface state solution
in a small region of momentum space, as shown in Fig. 5(b).

Figure 5: Gapped surface states on a single surface of thin film Cd3As2 in the
presence of a small in-plane magnetic field, $h_{s}<2M_{n}$. The surface
states are shown in red and the bulk bands in black. (a) The surface states
and bulk bands for $h_{s}=M_{n}/2$. (b) The surface state at $h_{s}=M_{n}$.
where the upper surface state has already merged with the bulk leaving only
the lower surface state solution present. Both panels use the model parameters
from Ref. [29].

When the field is increased to $h=2M_{n}$, the remaining surface state merges
with the bulk bands, and the $\alpha=-$ bulk bands un-invert, and for fields
$h_{s}>2M_{n}$, no surface state solutions arise.

##  V Orbital effect of Out-Of-Plane Magnetic Field

We now include the orbital effects of an out-of-plane magnetic field field,
which reorganizes bands into Landau levels. Experimentally, a crossing of the
zeroth particle and hole Landau levels as a function of thickness at an out-
of-plane field value $B_{c}$ has been reported [25]. In subsection V.1, we use
our model to calculate $B_{c}$ in terms of k.p parameters, and in subsection
V.2 we consider the case of a large in-plane Zeeman field and an out-of-plane
orbital field and calculate the Landau level spectrum.

###  V.1 Strictly out-of-plane field results

We first consider a purely out-of-plane field $\mathbf{B}=B\hat{\mathbf{z}}$.
Using the Hamiltonian in Eq. (1) to order $\mathcal{O}(k^{2})$, we introduce a
magnetic field via the Peierls substitution,
$\mathbf{k}\rightarrow-i\bm{\nabla}+e\mathbf{A}$. We work in the Landau gauge
$\mathbf{A}=By\hat{\mathbf{x}}$. We then define raising and lowering operators
in the conventional way

$\displaystyle
a^{\dagger}=\frac{l_{B}}{\sqrt{2}}\left(k_{x}-\nabla_{y}+\frac{y}{l_{B}^{2}}\right),$
(33a) $\displaystyle
a=\frac{l_{B}}{\sqrt{2}}\left(k_{x}+\nabla_{y}+\frac{y}{l_{B}^{2}}\right),$
(33b)

where $l_{B}^{2}=1/(|e|B)$ with electron charge $e$. The Hamiltonian in Eq.
(1), up to quadratic order in k, takes the form

$\displaystyle H$
$\displaystyle=\left[M_{n}-\omega\left(a^{\dagger}a+\frac{1}{2}\right)\right]\tau_{z}s_{0}$
$\displaystyle+\frac{A}{\sqrt{2}l_{B}}\left[\tau_{x}s_{z}\left(a+a^{\dagger}\right)+i\tau_{y}s_{0}\left(a-a^{\dagger}\right)\right],$
(34)

where $\omega=2M_{xy}/l_{B}^{2}$ and $M_{n}$ is defined in Eq. (7). Here, we
have neglected the band-diagonal contribution $\epsilon_{n}(\mathbf{k})$,
which, as shown in Appendix F introduces particle-hole asymmetry but does not
change the magnetic field value when the particle and hole Landau levels
cross.

Using the trial wavefunctions $\psi_{0}=\left(0,0,\ket{0},0\right)^{T}$ and
$\psi_{0}=\left(0,\ket{0},0,0\right)^{T}$, for the zeroth Landau levels, we
find the particle and hole zeroth Landau level energies,

$\displaystyle E_{0\pm}$
$\displaystyle=\pm\left|M_{n}-\frac{\omega}{2}\right|.$ (35)

For each subband $n$, there is a crossing of the particle and hole zeroth
Landau level energies at a critical magnetic field $B_{c}(n)$, which is a
function of the thickness of the sample,

$\displaystyle
B_{c}(n)=\frac{1}{|e|M_{xy}}\left[M_{0}-M_{z}\left(\frac{n\pi}{L}\right)^{2}\right]$
(36)

For the $n=2$ subbands, which are most relevant for experimental sample
thicknesses of $18-20$ nm reported in Refs. [25, 26] and using realistic
$k\cdot p$ model parameters from Ref. [29] (see also Appendix E), we find the
critical field $B_{c}(n=2)\approx 14$ T for $L=18$ nm thick samples. For a
slightly larger thickness $L=20$ nm, we find $B_{c}(n=2)\approx 38$ T.
Experiment [25] reports a critical magnetic field of $\approx 10$ T for $20$
nm thick films, in reasonable agreement with the model estimate.

###  V.2 Tilted field results

In the presence of a strong in-plane Zeeman coupling ($h_{s}>2M_{n}$), we have
shown previously that the low-energy physics is controlled by the bands
arising from the $p$-orbitals at $\Gamma$. The $s$-orbital bands at the
$\Gamma$ point are split off in energy by the large Zeeman field
$h_{s}>2M_{n}$. We thus include an out-of-plane orbital field to the effective
Hamiltonian $\mathcal{H}_{p}$ in Eq. (17), in which we keep only
$\mathcal{O}(k^{2})$, to find

$\displaystyle\mathcal{H}_{p}$
$\displaystyle=\left(\begin{array}[]{cc}-M_{n}+\omega(a^{\dagger}a+\frac{1}{2})&\alpha\left(a^{\dagger}\right)^{2}+\frac{h_{p}^{3}}{\Delta^{2}}e^{-3i\phi}\\\
\alpha^{*}a^{2}+\frac{h_{p}^{3}}{\Delta^{2}}e^{3i\phi}&-M_{n}+\omega(a^{\dagger}a+\frac{1}{2})\end{array}\right)$
(39)

where $\alpha=2A^{2}e^{-i\phi}/(l_{B}^{2}h_{s})$. In the following, we assume
that the in-plane field is oriented along the $x$-axis and set $\phi=0$ unless
otherwise stated. Since this direction is perpendicular to a mirror plane,
there exist emergent gapless Dirac states at $\mathbf{k}_{c}=\sigma
k_{c}\hat{\mathbf{y}}$ as described in Sec. III.3.

Figure 6: Schematic plot of the Landau levels arising from $\mathcal{H}_{p}$
in Eq. (39) and their transition from being twofold degenerate near the two
Dirac points (green) to exhibiting no additional degeneracy (red). The Dirac
points occur at $\epsilon=-M_{n}$ and the transition from double to single
degenerate Landau levels occurs at $\epsilon\sim-M_{n}+h_{p}^{3}/\Delta^{2}$
for the particle Landau levels.

As schematically shown in Fig. 6, there are two distinct regimes: for low-
lying Landau levels (shown in green), there is an emergent new quantum number
$\sigma$, which denotes the Dirac node and gives rise to a twofold degeneracy
of the Landau levels. Above some Landau level index $\nu_{c}$, the Dirac
states are no longer relevant, and this twofold Landau level degeneracy is no
longer present. We consider here only the positive $\nu$ – the negative $\nu$
can be obtained using particle-hole symmetry of our model (neglecting
$\epsilon(\mathbf{k})$). The effect of including the particle-hole symmetry
breaking term $\epsilon(\mathbf{k})$ is discussed in detail in Appendix F.

We proceed by first finding analytical solutions for the two cases,
$\nu>\nu_{c}$ and $\nu<\nu_{c}$, and then show numerical results demonstrating
the crossover and estimate the crossover index $\nu_{c}$.

####  V.2.1 High energy Landau Levels, $\nu>\nu_{c}$

We first consider high energy states and solve Eq. (39) for large Landau Level
indices $\nu$ such that $\alpha\nu\gg h_{p}^{3}/\Delta^{2}$ and we can neglect
the $p-$orbital Zeeman field $h_{p}$. In this case, the trial wavefunction

$\displaystyle\psi_{\nu}$
$\displaystyle=\left(\begin{array}[]{c}\tilde{C}_{1,\nu}\ket{\nu}\\\
\tilde{C}_{2,\nu}\ket{\nu-2}\end{array}\right),$ (42)

is an eigenstate of $\mathcal{H}_{p}$ [Eq. (39)] with eigenvalues

$\displaystyle
E_{\pm}(\nu)=-M_{n}+\omega\left(\nu-\frac{1}{2}\right)\pm\sqrt{\omega^{2}+|\alpha|^{2}\nu(\nu-1)}.$
(43)

We see that the eigenvalues are characteristic of Landau levels in a quadratic
band touching [32]; schematically, in Fig. 6 they are shown in red.

####  V.2.2 Low energy Landau Levels, $\nu<\nu_{c}$

At low energies, the physics is dominated by the Dirac states, which occur at
$\mathbf{k}_{c}=\pm\frac{\sqrt{h_{p}^{3}h_{s}}}{A\Delta}\hat{\mathbf{y}}$.
Using the expansion of Eq. (22) about these points, $\mathbf{k}=\delta
k_{x}\hat{\mathbf{x}}+(\sigma k_{c}+\delta k_{y})\hat{\mathbf{y}}$, where
$\sigma=\pm$ denotes the two band touching points, and inserting the Peierls
substitution for $\delta k$, the Hamiltonian is

$\displaystyle\mathcal{H}_{\sigma}$
$\displaystyle=\left(\begin{array}[]{cc}-M_{n}&-2i\frac{A^{2}}{h_{s}}k_{c}\left(\delta
k_{x}-i\delta k_{y}\right)\\\ 2i\frac{A^{2}}{h_{s}}k_{c}\left(\delta
k_{x}+i\delta k_{y}\right)&-M_{n}\end{array}\right).$ (46)

Inserting the orbital field via Peierls substitution for the deviations
$\delta\mathbf{k}$, we get the Hamiltonian of each Dirac state

$\displaystyle\mathcal{H}_{\sigma}$
$\displaystyle=\left(\begin{array}[]{cc}-M_{n}&\Gamma a^{\dagger}\\\
\Gamma^{*}a&-M_{n}\end{array}\right).$ (49)

where

$\Gamma=-i\frac{2^{3/2}A^{2}k_{c}}{h_{s}l_{B}}e^{-2i\phi}.$ (50)

The wavefunction is $\psi_{\nu}=\left(-i\ket{\nu},\ket{\nu-1}\right)^{T},$
while the eigenenergies read

$\displaystyle E_{\sigma\pm}(\nu)$
$\displaystyle=-M_{n}\pm\left|\Gamma\right|\sqrt{|\nu|}\,$ (51)

and we see that the Landau level spectrum about each Dirac point is identical,
yielding a two-fold degeneracy for each Landau level, as indicated by the
green Landau levels in Fig. 6.

Notably, if the in-plane field is rotated away from the high-symmetry
directions, a finite mass term $m(\phi)$ arises. If we also now include the
out-of-plane Zeeman field $h_{z}$, the Hamiltonian for the Landau levels due
to a single Dirac fermion is given by

$\displaystyle\mathcal{H}_{\sigma}$
$\displaystyle=-M_{n}+\left(\begin{array}[]{cc}\sigma m(\phi)+h_{z}&\Gamma
a^{\dagger}\\\ \Gamma^{*}a&-\sigma m(\phi)-h_{z}\end{array}\right)$ (54)

The wavefunction for the zeroth Landau level is $\psi_{0}=(\ket{0},0)^{T}$,
giving for the zeroth Landau level energies

$\displaystyle E_{\sigma}(0)$ $\displaystyle=-M_{n}+\sigma m(\phi)+h_{z}.$
(55)

We see that an out-of-plane Zeeman field $h_{z}$ shifts the zeroth Landau
levels of both Dirac fermions in the same direction, as shown in Fig. 7.
However, the mass $m(\phi)$ term splits the zeroth Landau levels of the two
Dirac fermions if the in-plane component of the field does not point in a high
symmetry direction. We note that in experiment, Ref. [26], the presence of the
$\nu=0$ plateau for large in-plane fields indicates the splitting of the
zeroth Landau levels. If the origin of that splitting is the finite $m(\phi)$,
then we expect that this plateau should disappear or become harder to see when
the in-plane component of magnetic field is pointing in a high-symmetry
direction. However, the splitting may also be due to the appearance of
correlated, “emergent” spontaneously valley-polarized states. Future
experiment should help to resolve which of the possibilities is realized in
practice.

####  V.2.3 Crossover regime

We have established that the Hamiltonian in Eq. (39) exhibits twofold
degenerate Landau levels for small $\nu$, one coming from each band touching
point at $\mathbf{k}=\sigma k_{c}\hat{\mathbf{y}}$, with eigenvalues given by
Eq. (51). For large $\nu$, we found singly degenerate Landau levels as one can
eventually neglect the $p-$orbital Zeeman field, and the Landau levels lose
this twofold degeneracy, with energies given by Eq. (43). The crossover
between these two regimes can be approached numerically by considering a
finite number of Landau levels coupling via Eq. (39). Taking a total of $1000$
levels into account, we show the resulting Landau level spectrum in Fig. 7.
The spectrum is particle-hole symmetric around energy $-M_{n}$, as we are
neglecting the band diagonal $\epsilon_{n}(\mathbf{k})$ as well as $\omega$,
and the index $\nu=500$ corresponds to the first Landau level above $-M_{n}$
(lowest green line in Fig. 6). Clearly, the Landau levels exhibit a crossover
between being doubly degenerate at low $|\nu-500|<\nu_{c}$ to being singly
degenerate at larger $|\nu-500|>\nu_{c}$. In the figure, the crossover occurs
at $|\nu-500|=8$.

Figure 7: $1000$ Landau levels are coupled via Eq. (39) and the resulting
Hamiltonian is diagonalized numerically. The Landau levels plotted are near
$\epsilon=-M_{n}\approx-5.5$ meV, denoted by the red line. The Landau level
energies are calculated at $l_{B}=200$ nm, $h_{s}=5M_{n}/2$,
$h_{p}^{3}/\Delta^{2}=h_{s}/20$, and the values of the other coefficients are
taken from Ref. [29]. We have also neglected $\omega$ as it only introduces
particle-hole asymmetry. (a) Landau level spectrum in the absence of $h_{z}$.
The Landau levels are numbered such that $\nu=501$ occurs at energy
$-M_{n}+|\Gamma|$ (corresponding to $\nu=1$ in Eq. (51)), and $\nu<500$
correspond to Landau levels with energies at or below $-M_{n}$. (b) Landau
level spectrum in the presence of $h_{z}=2g_{p}B_{z}$, where we use
$g_{p}=1/20$ to emphasize the shift of the zeroth Landau levels of the Dirac
fermions in the presence of $h_{z}$.

The crossover between these two regimes occurs when the Landau level energy is
on the order of $E\approx-M_{n}+\frac{h_{p}^{3}}{\Delta^{2}}$ as this is the
energy of the $\alpha=-$ $p$-orbital band at the $\Gamma$ point, and
illustrated in Fig. 6. We can, therefore, estimate the number of Landau levels
with a twofold degeneracy, $\nu_{c}$, using Eq. (51) to obtain

$\displaystyle\nu_{c}$ $\displaystyle\approx
2\frac{h_{p}^{6}}{\Delta^{4}|\Gamma|^{2}}=2\frac{h_{p}^{3}h_{s}l_{B}^{2}}{2^{3}\Delta^{2}A^{2}}\,.$
(56)

The factor of two in front accounts for the double degeneracy. For the
parameters used in Fig. 7, we find $\nu_{c}\approx 7.7$, which is in good
agreement with the numerical results.

Ref. [26] observes quantum Hall plateaus in a magnetic field with a fixed
magnitude $|\bm{B}|=\sqrt{B_{\perp}^{2}+B_{\parallel}^{2}}$ tilted out of the
plane at varying angles $\theta_{tilt}$,
$\bm{B}=(B_{\parallel}\cos\theta_{tilt},0,B_{\perp}\sin\theta_{tilt})$. For
smaller tilt angles, which correspond to larger in-plane magnetic fields, the
quantum Hall plateaus are observed at odd integer filling factors
$\nu=1,3,5,...$ while for fields with a larger tilt angle quantum Hall
plateaus are observed at every integer filling factor. As the magnetic field
magnitude is kept fixed while the field is rotated out of the plane, both
$h_{s}$, $h_{p}$, and $l_{B}$ decrease, decreasing $\nu_{c}$ in agreement with
our model. For in plane magnetic field value $B_{\parallel}=16$ T, which
corresponds to the high field regime, $g_{s}=12$,
$h_{p}^{3}/\Delta^{2}=h_{s}/20$, and $A=111.6$ meV$\cdot$nm from Ref. [29], we
calculate that $\nu_{c}<2$ for $B_{\perp}<6.8$ T, in good agreement with
experiment which sees the transition between these two regimes occur between
$B_{\perp}=3.62$ T and $4.79$ T.

##  VI Conclusions

In this article, we have derived a low-energy $k\cdot p$ model for $[001]$
thin film Cd3As2 in the presence of a magnetic field applied in an arbitrary
direction. At zero field and above a critical thickness, specific subbands $n$
of thin film Cd3As2 are in a quantum spin Hall insulator regime. Importantly,
we find that the $p$-orbital Zeeman field enters the Hamiltonian only to cubic
order in the applied field, while it is linear for the $s$-orbitals. This
distinction is important to understand why the model enters a 2D Dirac
semimetal phase in the presence of a sufficiently large in-plane Zeeman field
($h_{s}>2M_{n}$) that points perpendicular to a high symmetry mirror plane.
The nodal points appear on a mirror plane and become gapped when the applied
in-plane field is rotated away from the high-symmetry directions in plane. The
fact that the system is initially gapped for small in-plane fields but enters
a gapless regime at larger fields (if applied perpendicular to a mirror plane)
can explain the experimental observation of a decrease in longitudinal
resistance for sufficiently large in-plane fields [26]. From the realistic
model parameters, for sample sizes $L=19$ nm with an effective $g-$factor of
$g_{s}=12$ we estimate the large field regime to begin for the $n=2$ subbands
at $B_{\parallel}\sim 16$ T. The energy gap follows a $\sin(4\phi)$ behavior
as the field is rotated through the plane, which we predict to be observable
when fitting the resistance at finite temperatures to an Arrhenius law.
Further, we have investigated the fate of the 1D helical edge states of the
zero field QSH insulator in the presence of an in-plane magnetic field.
Neglecting the small $p$-orbital Zeeman field, we showed how they gap out and
evolve in the presence of an in-plane field that couples the $s$-orbitals.
While they gap out at infinitesimal fields, consistent with the explicit
breaking of time-reversal symmetry, we further showed their complete removal
as the in-plane field $h_{s}$ is increased. We further study the inclusion of
a perpendicular orbital magnetic field and show that, in the absence of an in-
plane Zeeman field, a crossing of the zeroth Landau levels as a function of
film thickness occurs. This behavior is consistent with the experimental
observations reported in Ref. [25]. Lastly, we investigated how the system
responds to a perpendicular orbital magnetic field in the 2D Dirac semimetal
regime and find that the resulting Landau level spectrum exhibits a twofold
degeneracy up to a Landau level index $\nu_{c}\propto
h_{p}^{3}h_{s}l_{B}^{2}/\Delta^{2}A^{2}$, where $A$ sets the Fermi velocity in
the bulk model and $\Delta$ is the energy separation of the low-energy
manifold from higher bands. The double degeneracy arises from the emergence of
a new quantum number associated with the two 2D Dirac nodes that appear at
large in-plane fields. The Landau level spectrum crosses over to be singly
degenerate at larger energies for $\nu>\nu_{c}$. We present analytical results
for the Landau level spectrum of the two regimes and show numerical results
for the crossover. These results can explain recent quantum Hall measurements
presented in Ref. [26], which reported an initial jump of the quantum Hall
plateaus by two conductance quanta $(e^{2}/h$) at small out-of-plane fields
before returning to plateaus changing by one for larger out-of-plane fields.

##  VII Acknowledgments

We thank S. Stemmer for useful discussions that motivated this work. We also
acknowledge discussions with O. Heinonen, J.X. Zhu, and N. Fokkens. This work
was supported by the Center for the Advancement of Topological Semimetals
(CATS), an Energy Frontier Research Center funded by the U.S. Department of
Energy (DOE) Office of Science (SC), Office of Basic Energy Sciences (BES),
through the Ames National Laboratory under contract DE-AC02-07CH11358.
Research at Perimeter Institute is supported in part by the Government of
Canada through the Department of Innovation, Science and Economic Development
and by the Province of Ontario through the Ministry of Economic Development,
Job Creation and Trade.

##  Appendix A Derivation of the k.p model

In this Appendix, we derive the effective $k\cdot p$ model up to cubic order,
complementing the derivations presented before, up to second order [29]. Up to
quadratic order, it matches the model of Refs. [25, 29, 27]. We consider the
space spanned by the states
$\left\\{\left|S_{1/2},1/2\right\rangle,\left|P_{3/2},3/2\right\rangle,\left|S_{1/2},-1/2\right\rangle,\left|P_{3/2},-3/2\right\rangle\right\\}$.
This basis order corresponds to $s_{i}\otimes\tau_{j}$ ordering. In this
basis, the matrix structure of the symmetry transformations generating all
operations in the point group $D_{4h}$ read

$\displaystyle C_{2x}$ $\displaystyle=-is_{x}\tau_{z},$ (57) $\displaystyle
C_{2y}$ $\displaystyle=-is_{y}\tau_{0},$ (58) $\displaystyle P$
$\displaystyle=s_{0}\tau_{z},$ (59) $\displaystyle T$
$\displaystyle=is_{y}\tau_{0}K,$ (60) $\displaystyle C_{4z}$
$\displaystyle=\frac{1}{\sqrt{2}}(s_{0}\tau_{z}-is_{z}\tau_{0})$ (61)

Here, $P$ denotes spatial inversion, $T$ denotes time-reversal, and $K$
complex conjugation. These symmetry operations can be checked as follows:

  * •

$C_{2x}$ and $C_{2y}$: There are two ways of finding the operator. The first
is to rewrite the states in terms of spherical harmonics. For instance,
$\left|P_{3/2},3/2\right\rangle$ comes from $Y_{1,1}$ combined with a spin
half pointing up. Another way is purely in terms of the total angular
momentum. For that, we write down the operator for $J=1/2$ and for $J=3/2$.
For instance, for $J=3/2$, we need
$C_{2x}=\exp\left(-i\frac{2\pi}{2}J_{x}\right)$, where $J_{x}$ is the $x$
component of the spin $3/2$. We project onto the states of $M_{J}=\pm 3/2$. We
repeat the procedure for the $J=1/2$ case. The final operator is the direct
sum of the two.

  * •

$P$: The $p$ orbitals are odd under inversion, while the s orbitals are even;

  * •

$T$: Time reversal acts the usual way as these four states consist of two
Kramers doublets;

  * •

$C_{4z}$: we can use the total angular momentum operators, $J=1/2$ and
$J=3/2$, and compute how they rotate;

Finally, the mirror symmetry operations are:

$\displaystyle M_{xz}$ $\displaystyle=P\circ C_{2y}=-is_{y}\tau_{z},$ (62)
$\displaystyle M_{(x+y)z}$ $\displaystyle=C_{4z}\circ
M_{xz}=\frac{i}{\sqrt{2}}(-s_{x}\tau_{z}+s_{y}\tau_{0})$ (63) $\displaystyle
M_{yz}$ $\displaystyle=C_{4z}^{2}\circ M_{xz}=-is_{x}\tau_{0}$ (64)
$\displaystyle M_{(x-y)z}$ $\displaystyle=C_{4z}^{3}\circ
M_{xz}=-\frac{i}{\sqrt{2}}(s_{x}\tau_{z}+s_{y}\tau_{0}).$ (65)

In general, the Hamiltonian can be expressed as

$\displaystyle H$
$\displaystyle=\sum_{i,j}\sum_{\mathbf{k}}H_{ij}(\mathbf{k})c_{i\mathbf{k}}^{\dagger}c_{j\mathbf{k}}.$
(66)

We call $U_{g}$ the unitary part of the symmetry, acting on the orbital and
spin indices, and $R_{g}$ the representation of the symmetry in momentum
space. The symmetry constraints impose that

$\displaystyle U_{g}H\left(R_{g}^{-1}\mathbf{k}\right)U_{g}^{-1}$
$\displaystyle=H\left(\mathbf{k}\right).$ (67)

The $k\cdot p$ Hamiltonian satisfying the symmetries can be computed using
algebraic subroutines that incorporate all the symmetry operations. We use the
package MagneticKP for Mathematica [33]. Up to order $k^{3}$, we arrive at Eq.
(1) of the main text.

##  Appendix B Microscopic derivation of the Zeeman and spin-orbit terms

In this Appendix, we derive the effects of adding a magnetic field along a
general direction to thin films of Cd3 As2. This Appendix is organized in the
following way. First, in Sec. B.1 we list the different bases: the spherical
basis, the Cartesian one, and the one that diagonalizes the total angular
momentum. In Sec. B.2, we write down the spin-orbit coupling term on the
different bases. In Sec. B.3, we derive the terms in the Hamiltonian
corresponding to a magnetic field.

###  B.1 The different basis and how to relate them

Throughout this Appendix, we are going to change basis several times. First,
we write down all the operators on the spherical basis, ordered according to

$\displaystyle\left\\{\left|\uparrow\right\rangle,\left|\downarrow\right\rangle\right\\}\otimes\left\\{\left|s\right\rangle,\left|1,1\right\rangle,\left|1,0\right\rangle,\left|1,-1\right\rangle\right\\}$
$\displaystyle=\left\\{\left|s,\uparrow\right\rangle,\left|1,1,\uparrow\right\rangle,\ldots\right\\}\,.$
(68)

Here, $\ket{l,m}$ are the indices of the spherical harmonics $Y_{lm}$ denoting
the orbital angular momentum numbers. Another useful basis is the one with the
$s$ and $p_{x,y,z}$ orbitals,

$\left\\{\left|\uparrow\right\rangle,\left|\downarrow\right\rangle\right\\}\otimes\left\\{\left|s\right\rangle,\left|p_{x}\right\rangle,\left|p_{y}\right\rangle,\left|p_{z}\right\rangle\right\\}=\left\\{\left|s,\uparrow\right\rangle,\left|p_{x},\uparrow\right\rangle,\ldots\right\\}.$
(69)

The basic transformations come from writing down the spherical harmonics in
terms of the $p$ orbitals,

$\displaystyle Y_{1,-1}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(p_{x}-ip_{y}\right),$ (70)
$\displaystyle Y_{1,1}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(-p_{x}-ip_{y}\right),$ (71)
$\displaystyle Y_{0,0}$ $\displaystyle=p_{z}.$ (72)

Here, we do not explicitly write the spin index; the transformations are the
same for the up and down components. Conversely,

$\displaystyle p_{x}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(Y_{1,-1}-Y_{1,1}\right),$ (73)
$\displaystyle p_{y}$
$\displaystyle=\frac{i}{\sqrt{2}}\left(Y_{1,-1}+Y_{1,1}\right),$ (74)
$\displaystyle p_{z}$ $\displaystyle=Y_{1,0}.$ (75)

The rotation matrix from the spherical to the $p_{x,y,z}$ basis can be read
off from Eqs. (70), (71), and (72),

$U=\left(\begin{array}[]{cccccccc}1&0&0&0&0&0&0&0\\\
0&-\frac{1}{\sqrt{2}}&\frac{i}{\sqrt{2}}&0&0&0&0&0\\\ 0&0&0&1&0&0&0&0\\\
0&\frac{1}{\sqrt{2}}&\frac{i}{\sqrt{2}}&0&0&0&0&0\\\ 0&0&0&0&1&0&0&0\\\
0&0&0&0&0&-\frac{1}{\sqrt{2}}&\frac{i}{\sqrt{2}}&0\\\ 0&0&0&0&0&0&0&1\\\
0&0&0&0&0&\frac{1}{\sqrt{2}}&\frac{i}{\sqrt{2}}&0\end{array}\right).$ (76)

The first $4\times 4$ block corresponds to spin up, while the second one
corresponds to spin down. Finally, since spin-orbit coupling is large, we will
work with the total angular momentum basis

$\displaystyle\mathcal{H}$
$\displaystyle=\left\\{\left|S_{1/2},1/2\right\rangle,\left|S_{1/2},-1/2\right\rangle,\left|P_{3/2},3/2\right\rangle,\left|P_{3/2},1/2\right\rangle,\right.$
$\displaystyle\left.\left|P_{3/2},-1/2\right\rangle,\left|P_{3/2},-3/2\right\rangle,\left|P_{1/2},1/2\right\rangle,\left|P_{1/2},-1/2\right\rangle\right\\}.$
(77)

Later, we will consider the low-energy modes, in which we truncate the Hilbert
space, $\mathcal{H}_{trunc}$, spanned by the four states (note the
$\tau_{j}\otimes s_{i}$ ordering here):

$\left\\{\left|S_{1/2},1/2\right\rangle,\left|S_{1/2},-1/2\right\rangle,\left|P_{3/2},3/2\right\rangle,\left|P_{3/2},-3/2\right\rangle\right\\}.$
(78)

To change basis to the total angular momentum basis, we have to write down the
$J$ multiplets in terms of the original states. The $S_{1/2}$ are trivially
found from the $s$ orbitals,

$\displaystyle\left|S_{1/2},1/2\right\rangle$
$\displaystyle=\left|s,\uparrow\right\rangle,$ (79)
$\displaystyle\left|S_{1/2},-1/2\right\rangle$
$\displaystyle=\left|s,\downarrow\right\rangle.$ (80)

The four $P_{3/2}$ are constructed from the $L=1$, $S=1/2$ states as

$\displaystyle\left|P_{3/2},3/2\right\rangle$
$\displaystyle=\left|m_{\ell}=1,\uparrow\right\rangle,$
$\displaystyle=\frac{1}{\sqrt{2}}\left(-\left|p_{x},\uparrow\right\rangle-i\left|p_{y},\uparrow\right\rangle\right),$
(81) $\displaystyle\left|P_{3/2},1/2\right\rangle$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\sqrt{2}\left|m_{\ell}=0,\uparrow\right\rangle+\left|m_{\ell}=1,\downarrow\right\rangle\right),$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\sqrt{2}\left|p_{z},\uparrow\right\rangle-\frac{1}{\sqrt{2}}\left|p_{x},\downarrow\right\rangle-\frac{i}{\sqrt{2}}\left|p_{y},\downarrow\right\rangle\right),$
(82) $\displaystyle\left|P_{3/2},-1/2\right\rangle$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\left|m_{\ell}=-1,\uparrow\right\rangle+\sqrt{2}\left|m_{\ell}=0,\downarrow\right\rangle\right),$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\frac{1}{\sqrt{2}}\left|p_{x},\uparrow\right\rangle-\frac{i}{\sqrt{2}}\left|p_{y},\uparrow\right\rangle+\sqrt{2}\left|p_{z},\downarrow\right\rangle\right),$
(83) $\displaystyle\left|P_{3/2},-3/2\right\rangle$
$\displaystyle=\left|m_{\ell}=-1,\downarrow\right\rangle,$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\left|p_{x},\downarrow\right\rangle-i\left|p_{y},\downarrow\right\rangle\right),$
(84)

Finally, the two $P_{1/2}$ states are also found from the $L=1$, $S=1/2$ as

$\displaystyle\left|P_{1/2},1/2\right\rangle$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\sqrt{2}\left|m_{\ell}=1,\downarrow\right\rangle-\left|m_{\ell}=0,\uparrow\right\rangle\right),$
$\displaystyle=\frac{1}{\sqrt{3}}\left(-\left|p_{x},\downarrow\right\rangle-i\left|p_{y},\downarrow\right\rangle-\left|p_{z},\uparrow\right\rangle\right),$
(85) $\displaystyle\left|P_{1/2},-1/2\right\rangle$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\left|m_{\ell}=0,\downarrow\right\rangle-\sqrt{2}\left|m_{\ell}=-1,\uparrow\right\rangle\right),$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\left|p_{z},\downarrow\right\rangle-\left|p_{x},\uparrow\right\rangle+i\left|p_{y},\uparrow\right\rangle\right).$
(86)

The matrix that changes the basis from $s,p_{x,y,z}$ to the total angular
momentum is immediately found by collecting the results from Eqs. (79) to
(86),

$U_{J}=\left(\begin{array}[]{cccccccc}1&0&0&0&0&0&0&0\\\
0&0&-\frac{1}{\sqrt{2}}&0&\frac{1}{\sqrt{6}}&0&0&-\frac{1}{\sqrt{3}}\\\
0&0&-\frac{i}{\sqrt{2}}&0&-\frac{i}{\sqrt{6}}&0&0&\frac{i}{\sqrt{3}}\\\
0&0&0&\sqrt{\frac{2}{3}}&0&0&-\frac{1}{\sqrt{3}}&0\\\ 0&1&0&0&0&0&0&0\\\
0&0&0&-\frac{1}{\sqrt{6}}&0&\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{3}}&0\\\
0&0&0&-\frac{i}{\sqrt{6}}&0&-\frac{i}{\sqrt{2}}&-\frac{i}{\sqrt{3}}&0\\\
0&0&0&0&\sqrt{\frac{2}{3}}&0&0&\frac{1}{\sqrt{3}}\end{array}\right).$ (87)

Having all the transformation matrices, we are ready to address how the Zeeman
and the the spin-orbit terms are written in each basis.

###  B.2 The spin-orbit coupling

Before writing down the magnetic field terms, we present the different ways of
writing down the spin-orbit coupling matrix. The spin-orbit coupling can be
easily written on the spherical basis by writing the operators $\mathbf{S}$
for the spin $1/2$ and $\mathbf{L}$ for the spin $1$. The SOC written in
matrix form, with the basis ordered according to Eq. (68), reads

$H_{SO}=\Delta\mathbf{L}\cdot\mathbf{S}=\Delta\left(\begin{array}[]{cccccc}\frac{1}{2}&0&0&0&0&0\\\
0&0&0&\frac{1}{\sqrt{2}}&0&0\\\ 0&0&-\frac{1}{2}&0&\frac{1}{\sqrt{2}}&0\\\
0&\frac{1}{\sqrt{2}}&0&-\frac{1}{2}&0&0\\\ 0&0&\frac{1}{\sqrt{2}}&0&0&0\\\
0&0&0&0&0&\frac{1}{2}\end{array}\right).$ (88)

Rotating to the $p_{x,y,z}$ basis, of Eq. (69), we find, using the rotation
matrix of Eq. (69),

$\tilde{H}_{SO}=U^{\dagger}H_{SO}U=\frac{\Delta}{2}\left(\begin{array}[]{cccccccc}0&0&0&0&0&0&0&0\\\
0&0&-i&0&0&0&0&1\\\ 0&i&0&0&0&0&0&-i\\\ 0&0&0&0&0&-1&i&0\\\ 0&0&0&0&0&0&0&0\\\
0&0&0&-1&0&0&i&0\\\ 0&0&0&-i&0&-i&0&0\\\ 0&1&i&0&0&0&0&0\end{array}\right).$
(89)

It matches the form found in Ref. [29]. Finally, we rotate from the
$p_{x,y,z}$ to the total angular momentum basis,

$\displaystyle H_{SOC}^{\left(J\right)}=U_{J}^{\dagger}\tilde{H}_{SOC}U_{J}$
$\displaystyle=\Delta\left(\begin{array}[]{cccccccc}0&0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0&0\\\ 0&0&\frac{1}{2}&0&0&0&0&0\\\ 0&0&0&\frac{1}{2}&0&0&0&0\\\
0&0&0&0&\frac{1}{2}&0&0&0\\\ 0&0&0&0&0&\frac{1}{2}&0&0\\\ 0&0&0&0&0&0&-1&0\\\
0&0&0&0&0&0&0&-1\end{array}\right).$ (98)

This matrix is diagonal in the total angular momentum basis. It can also be
found simply from
$\mathbf{L}\cdot\mathbf{S}=\frac{1}{2}\left(\mathbf{J}^{2}-\mathbf{L}^{2}-\mathbf{S}^{2}\right)$.
For the $J=3/2$ multiplet (entries 3 to 6), coming from the $p$ electrons
$\left(l=1\right)$ and $s=1/2$, we get $\mathbf{L}\cdot\mathbf{S}=1/2$. For
the cases of $\left(J,L,S\right)=\left(1/2,1,1/2\right)$, we find
$\mathbf{L}\cdot\mathbf{S}=-1$, while for the states coming from $s$ orbitals,
$\left(J,L,S\right)=\left(1/2,0,1/2\right)$, yielding
$\mathbf{L}\cdot\mathbf{S}=0$.

###  B.3 The effects of SOC and magnetic fields

The Zeeman term can be easily written in the spherical basis by simply writing
$g_{i}\mu_{B}\mathbf{B}\cdot\left(\mathbf{L}+2\mathbf{S}\right)$ for all
degrees of freedom, following the basis defined in Eq. (68). Here, $g_{i}$
denotes the $g$-factor for a given orbital. In this subsection, for
compactness, we will set $\mu_{B}=1$, keeping in mind that all terms that
follow are proportional to $\mu_{B}$. The Zeeman term reads

$\displaystyle
H_{Z}=g_{s}\mathbf{B}\cdot\left(\mathbf{L}+2\mathbf{S}\right)+g_{p}\mathbf{B}\cdot\left(\mathbf{L}+2\mathbf{S}\right)$
$\displaystyle=\\!\begin{pmatrix}B_{z}g_{s}&0&0&0&g_{s}(B_{x}-iB_{y})&0&0&0\\\
0&2B_{z}g_{p}&\frac{g_{p}(B_{x}-iB_{y})}{\sqrt{2}}&0&0&g_{p}(B_{x}-iB_{y})&0&0\\\
0&\frac{g_{p}(B_{x}+iB_{y})}{\sqrt{2}}&B_{z}g_{p}&\frac{g_{p}(B_{x}-iB_{y})}{\sqrt{2}}&0&0&g_{p}(B_{x}-iB_{y})&0\\\
0&0&\frac{g_{p}(B_{x}+iB_{y})}{\sqrt{2}}&0&0&0&0&g_{p}(B_{x}-iB_{y})\\\
g_{s}(B_{x}+iB_{y})&0&0&0&-B_{z}g_{s}&0&0&0\\\
0&g_{p}(B_{x}+iB_{y})&0&0&0&0&\frac{g_{p}(B_{x}-iB_{y})}{\sqrt{2}}&0\\\
0&0&g_{p}(B_{x}+iB_{y})&0&0&\frac{g_{p}(B_{x}+iB_{y})}{\sqrt{2}}&-B_{z}g_{p}&\frac{g_{p}(B_{x}-iB_{y})}{\sqrt{2}}\\\
0&0&0&g_{p}(B_{x}+iB_{y})&0&0&\frac{g_{p}(B_{x}+iB_{y})}{\sqrt{2}}&-2B_{z}g_{p}\end{pmatrix}\\!.$
(99)

By rotating to the $p_{x,y,z}$ basis, we find

$\displaystyle\tilde{H}_{Z}=U^{\dagger}H_{Z}U,$
$\displaystyle=\begin{pmatrix}g_{s}B_{z}&0&0&0&g_{s}(B_{x}-iB_{y})&0&0&0\\\
0&g_{p}B_{z}&-ig_{p}B_{z}&ig_{p}B_{y}&0&g_{p}(B_{x}-iB_{y})&0&0\\\
0&ig_{p}B_{z}&g_{p}B_{z}&-ig_{p}B_{x}&0&0&g_{p}(B_{x}-iB_{y})&0\\\
0&-ig_{p}B_{y}g_{p}&ig_{p}B_{x}&g_{p}B_{z}&0&0&0&g_{p}(B_{x}-iB_{y})\\\
g_{s}(B_{x}+iB_{y})&0&0&0&-g_{s}B_{z}&0&0&0\\\
0&g_{p}(B_{x}+iB_{y})&0&0&0&-g_{p}B_{z}&-ig_{p}B_{z}&ig_{p}B_{y}\\\
0&0&g_{p}(B_{x}+iB_{y})&0&0&ig_{p}B_{z}&-g_{p}B_{z}&-ig_{p}B_{x}\\\
0&0&0&g_{p}(B_{x}+iB_{y})&0&-ig_{p}B_{y}&ig_{p}B_{x}&-g_{p}B_{z}\end{pmatrix}.$
(100)

Finally, in the total angular momentum basis, the Zeeman term reads

$\displaystyle H_{B}^{\left(J\right)}=U_{J}^{\dagger}H_{Z}U_{J},$
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}g_{s}\left(\begin{array}[]{cc}B_{z}&B_{x}-iB_{y}\\\
B_{x}+iB_{y}&-B_{z}\end{array}\right)&0\\\
0&g_{p}\left(\begin{array}[]{cccccc}2B_{z}&\frac{2\left(B_{x}-iB_{y}\right)}{\sqrt{3}}&0&0&\frac{\left(B_{x}-iB_{y}\right)}{\sqrt{6}}&0\\\
\frac{2(B_{x}+iB_{y})}{\sqrt{3}}&\frac{2}{3}B_{z}&\frac{4}{3}\left(B_{x}-iB_{y}\right)&0&-\frac{1}{3}\left(\sqrt{2}B_{z}\right)&\frac{\left(B_{x}-iB_{y}\right)}{3\sqrt{2}}\\\
0&\frac{4}{3}\left(B_{x}+iB_{y}\right)&-\frac{2}{3}B_{z}&\frac{2\left(B_{x}-iB_{y}\right)}{\sqrt{3}}&-\frac{\left(B_{x}+iB_{y}\right)}{3\sqrt{2}}&-\frac{1}{3}\left(\sqrt{2}B_{z}\right)\\\
0&0&\frac{2\left(B_{x}+iB_{y}\right)}{\sqrt{3}}&-2B_{z}&0&-\frac{\left(B_{x}+iB_{y}\right)}{\sqrt{6}}\\\
\frac{\left(B_{x}+iB_{y}\right)}{\sqrt{6}}&-\frac{1}{3}\left(\sqrt{2}B_{z}\right)&-\frac{\left(B_{x}-iB_{y}\right)}{3\sqrt{2}}&0&\frac{B_{z}}{3}&\frac{1}{3}\left(B_{x}-iB_{y}\right)\\\
0&\frac{\left(B_{x}+iB_{y}\right)}{3\sqrt{2}}&-\frac{1}{3}\left(\sqrt{2}B_{z}\right)&-\frac{\left(B_{x}-iB_{y}\right)}{\sqrt{6}}&\frac{1}{3}\left(B_{x}+iB_{y}\right)&-\frac{B_{z}}{3}\end{array}\right)\end{array}\right).$
(111)

From the form of the matrix, we can already find many important properties.
Within each subspace of total angular momentum, the field couples as
$g_{J}\mathbf{B}\cdot\mathbf{J}$, with $J$ the total angular momentum and
$g_{J}$ the corresponding Landé g-factor [34]

$g_{J}=\frac{3}{2}+\frac{S\left(S+1\right)-L\left(L+1\right)}{2J\left(J+1\right)}.$
(112)

Within a given angular momentum multiplet, the Zeeman term is proportional to
the total angular momentum $J$ [34]. $J_{x}$ and $J_{y}$ do not couple
$M_{J}=3/2$ with $M_{J}=-3/2$. For this multiplet, therefore, only $J_{z}$
will matter (and, therefore, only $B_{z}$). In fact, for the $J=3/2$ states,
the Landé g-factor is $g_{J=3/2}\left(J=3/2,L=1,S=1/2\right)=4/3$ while for
the $J=1/2$ multiplet of $p$ electrons,
$g_{J=1/2}\left(J=1/2,L=1,S=1/2\right)=2/3$. The diagonal entries of the
matrix are thus $g_{3/2}M_{J}=\frac{4}{3}\times(\pm\frac{3}{2})=\pm 2$,
$g_{3/2}M_{J}=\frac{4}{3}\times(\pm\frac{1}{2})=\pm\frac{2}{3}$, and
$g_{1/2}M_{J}=\frac{2}{3}\times(\pm\frac{1}{2})=\pm\frac{1}{3}$. The $s$
electrons behave purely as spin 1/2 objects, with a Landé factor
$g_{J=1/2}\left(J=1/2,L=0,S=1/2\right)=2$. The diagonal Zeeman coupling is
proportional to $g_{J}M_{J}$ and thus $g_{1/2}=2\times(\pm\frac{1}{2})=\pm 1$.
As for the $s$ electrons, only the spin is coupled, and $g_{s}s=2\times
1/2=1$.

If we confine ourselves to the low-energy states [Eq. (78)] by simply
projecting $H_{B}^{\left(J\right)}$ to those states, we find
$H_{P}=PH_{B}^{\left(J\right)}P$, with a matrix representation

$H_{P}=\begin{pmatrix}g_{s}B_{z}&g_{s}(B_{x}-iB_{y})&0&0\\\
g_{s}(B_{x}+iB_{y})&-g_{s}B_{z}&0&0\\\ 0&0&2g_{p}B_{z}&0\\\
0&0&0&-2g_{p}B_{z}\end{pmatrix}.$ (113)

Here, $P$ is the projector onto the low-energy subspace. We find that only the
$z$ component of the magnetic field couples to the $J=3/2$, $M_{J}=\pm 3/2$
states. This is expected as, within the $J=3/2$ subspace, the Zeeman term is
$g_{J}\mathbf{B}\cdot\mathbf{J}$ and the $x$ and $y$ components cannot couple
two states with $M_{z}$ differing by $3$.

##  Appendix C Low-energy effective coupling of p-electrons to in-plane
fields: perturbation theory

In this Appendix, we consider how the components of the magnetic field
parallel to the $xy$ plane couple to the $p$ electrons by performing a
perturbative calculation. We can already infer that we need a process that
couples the states of $M_{J}=\pm 3/2$. Such processes need to change $M_{z}$
by 3 units, which requires a cubic coupling with the magnetic field. We now
demonstrate this by an explicit calculation.

We will be focusing on states at the $\Gamma$ point. In what follows, we
neglect the effects of crystal fields that preserve the azimuthal quantum
numbers and, therefore, do not alter our conclusions. We call $E_{s}$ the
energy of the s orbitals and $\Delta_{1}$, $\Delta_{2}$ and $\Delta_{3}$ the
three different energies of the different $p$ doublets. Using the basis as
defined in Eq. (77), the Hamiltonian $H_{0}$ is diagonal, with energy $E_{s}$
for the s-orbitals and energies $\Delta_{1},\Delta_{2}$ for the $p$-orbitals
of $J=3/2$ and azimuthal numbers $\pm 3/2$ and $\pm 1/2$, respectively. The
energy of the p-orbitals of $J=1/2$ is $\Delta_{3}$.

These masses are entering to produce differences in the energy of the bands,
to zeroth order, and will lead to the $s$ states $\left|S_{1/2},\pm
1/2\right\rangle$ and the states $\left|P_{3/2},\pm 3/2\right\rangle$ close to
the Fermi level. We decompose the Hamiltonian as

$H=H_{0}+H_{B}^{\left(J\right)},$ (114)

where the perturbation $H_{B}^{\left(J\right)}$ is the term with magnetic
field, Eq. (111). We are going to follow the approach of Ref. [35]. We define
the resolvent operator

$R=\frac{Q}{E_{0}-H_{0}},$ (115)

with $Q=1-P$, where $P$ projects onto the four-dimensional set of states
$\\{\ket{S_{1/2},\pm 1/2},\ket{P_{3/2},\pm 3/2}\\}$, and $E_{0}$ the
eigenenergies of $H_{0}$. With the projector $P$ and the resolvent $R$, the
different corrections of the Hamiltonian can be found in a straightforward way
as [35]

$\displaystyle\Delta H^{\left(1\right)}$
$\displaystyle=PH_{B}^{\left(J\right)}P,$ (116) $\displaystyle\Delta
H^{\left(2\right)}$
$\displaystyle=PH_{B}^{\left(J\right)}RH_{B}^{\left(J\right)}P,$ (117)
$\displaystyle\Delta H^{\left(3\right)}$
$\displaystyle=PH_{B}^{\left(J\right)}RH_{B}^{\left(J\right)}RH_{B}^{\left(J\right)}P-PH_{B}^{\left(J\right)}R^{2}H_{B}^{\left(J\right)}PH_{B}^{\left(J\right)}P.$
(118)

All the terms $\Delta H^{\left(i\right)}$ are effective $4\times 4$ models,
written in the subspace spanned by the states listed in Eq. (78).

The first correction, $\Delta H^{\left(1\right)}$, is equivalent to neglecting
the effects of the high-energy bands and leads to the projected Hamiltonian of
Eq. (113). In what follows, we will consider only the effects on the $P_{3/2}$
states as the higher-order contributions for the $s$ states vanish.

The term $\Delta H^{\left(2\right)}$ is an isotropic contribution that adds to
the diagonal components and reads

$\Delta
H^{\left(2\right)}=u^{\left(2\right)}\left|\mathbf{B}_{\parallel}\right|^{2}\left(\begin{array}[]{cc}1&0\\\
0&1\end{array}\right).$ (119)

Here, $\mathbf{B}_{\parallel}=(B_{x},B_{y},0)$ and

$u^{\left(2\right)}=\frac{9\Delta_{1}-\Delta_{2}-8\Delta_{3}}{6\left(\Delta_{1}-\Delta_{2}\right)\left(\Delta_{1}-\Delta_{3}\right)}.$
(120)

Finally, the third-order correction reads $\Delta H^{\left(3\right)}=\Delta
H_{\perp}^{\left(3\right)}+\Delta H_{\parallel}^{\left(3\right)}$, with

$\Delta
H_{\parallel}^{\left(3\right)}=u_{\parallel}^{\left(3\right)}\left(\begin{array}[]{cc}0&\left(B_{x}-iB_{y}\right)^{3}\\\
\left(B_{x}+iB_{y}\right)^{3}&0\end{array}\right).$ (121)

and $\Delta H_{\perp}^{\left(3\right)}=u_{\perp}^{\left(3\right)}s_{z}$
Explicitly,

$\displaystyle u_{\parallel}^{\left(3\right)}$
$\displaystyle=\frac{\left(9\Delta_{1}-\Delta_{2}-8\Delta_{3}\right)\left(3\Delta_{1}+\Delta_{2}-4\Delta_{3}\right)}{18\left(\Delta_{1}-\Delta_{2}\right)^{2}\left(\Delta_{1}-\Delta_{3}\right)^{2}},$
(122) $\displaystyle u_{\perp}^{\left(3\right)}$
$\displaystyle=\frac{\left(3\Delta_{1}+\Delta_{2}-4\Delta_{3}\right)^{2}}{18\left(\Delta_{1}-\Delta_{2}\right)^{2}\left(\Delta_{1}-\Delta_{3}\right)^{2}}.$
(123)

Decomposing in terms of the $\mathbf{s}$ and $\bm{\tau}$ matrices, we find

$\displaystyle\Delta H^{\left(3\right)}$
$\displaystyle=\mu_{B}^{3}g_{p}\left(\frac{\tau_{0}-\tau_{z}}{2}\right)\left\\{B_{x}\left(B_{x}^{2}-3B_{y}^{2}\right)s_{x}\right.$
$\displaystyle\left.-B_{y}\left(B_{y}^{2}-3B_{x}^{2}\right)s_{y}+r\left|\mathbf{B}_{\perp}\right|^{2}B_{z}s_{z}\right\\},$
(124)

$r=u_{\perp}^{\left(3\right)}/u_{\parallel}^{\left(3\right)}$ a dimensionless
number of order one.

Adding all contributions, we find an effective coupling with the $\mathbf{B}$
field shown in Eq. (10) of the main text.

##  Appendix D Derivation of $\mathcal{H}_{p}$

In this Appendix, we explicitly show that the specific relations between the
third-order coefficients in the $k\cdot p$ model that we use in the main text
to simplify the resulting expressions do not change any of our conclusions. In
other words, it is not necessary to make assumptions about the coefficients
$C_{3a,b,c,d,e}$ as is done in the main text, and it will be shown that for
arbitrary coefficients, the form of the effective Hamiltonian
$\mathcal{H}_{p}$ in the high field limit takes the same form as in the main
text. To begin, our full Hamiltonian is given by

$\displaystyle\mathcal{H}$
$\displaystyle=\left(\begin{array}[]{cccc}M_{n}&A_{+}&h_{s}e^{-i\phi}&F\\\
A_{-}&-M_{n}&F&\tilde{h}_{p}e^{-3i\phi}\\\
h_{s}e^{i\phi}&F^{*}&M_{n}&-A_{-}\\\
F^{*}&\tilde{h}_{p}e^{3i\phi}&-A_{+}&-M_{n}\end{array}\right)$ (129)

in the basis
$(\ket{S,1/2},\ket{P_{3/2},3/2},\ket{S,-1/2},\ket{P_{3/2},-3/2})$, and we use
the shorthand $\tilde{h}_{p}=h_{p}^{3}/\Delta^{2}$ and

$\displaystyle A_{\pm}$ $\displaystyle=\tilde{A}k_{\pm}\mp
iC_{3a}k_{x}k_{y}k_{\mp}+C_{3b}(k_{x}^{3}\pm ik_{y}^{3}),$ (130a)
$\displaystyle F$
$\displaystyle=iC_{3e}k_{x}k_{y}k_{z}+C_{3d}k_{z}(k_{x}^{2}-k_{y}^{2}),$
(130b)

where we have conveniently written $\tilde{A}=A-C_{3c}k_{z}^{2}$ and use the
shorthand $k_{z}=n\pi/L$. As the bands with $p$-orbital character at $\Gamma$
have energies $\epsilon\approx-M_{n}$, we may take the $s$-orbital subspace of
the Schrödinger equation and solve for the $s$-orbital wavefunctions in terms
of the $p$-orbital wavefunctions in the limit $h_{s}\gg M_{n}$:

$\displaystyle s_{1/2}$
$\displaystyle\approx\frac{-F^{*}p_{3/2}+A_{-}p_{-3/2}}{h_{s}e^{i\phi}}$
(131a) $\displaystyle s_{-1/2}$
$\displaystyle\approx\frac{-A_{+}p_{3/2}-Fp_{-3/2}}{h_{s}e^{-i\phi}}$ (131b)

Substituting these solutions back into the equations for the $p$-orbitals we
obtain a low energy theory for the relevant bands near $\epsilon=-M_{n}$,
$\mathcal{H}_{p}$, with

$\displaystyle\mathcal{B}_{1}$
$\displaystyle=\frac{A_{-}^{2}}{h_{s}}e^{-i\phi}-\frac{F^{2}}{h_{s}}e^{i\phi}+\tilde{h}_{p}e^{-3i\phi},$
(132a) $\displaystyle\mathcal{B}_{2}$
$\displaystyle=-\left[\frac{F^{*}A_{-}}{h_{s}}e^{-i\phi}+\frac{FA_{+}}{h_{s}}e^{i\phi}\right].$
(132b)

$\mathcal{B}_{1}$ to $\mathcal{O}(k^{2})$ is given by

$\displaystyle\mathcal{B}_{1}^{(0)}$
$\displaystyle=\frac{\tilde{A}^{2}k_{-}^{2}}{h_{s}}e^{-i\phi}+\tilde{h}_{p}e^{-3i\phi}.$
(133)

$\mathcal{B}_{1}^{(0)}=0$ occurs at
$k_{c0}=\sqrt{\tilde{h}_{p}h_{s}/\tilde{A}^{2}}$ and
$\theta_{c\sigma}=\phi\pm\pi/2$. Inclusion of higher order in $k$ terms
slightly shifts $k_{c}$ away from $k_{c0}$ and gives rise to a mass term which
arises from $\mathcal{B}_{2}$, as we will show below. For coefficients in the
atomic limit $\mathcal{B}_{1}$ is

$\displaystyle\mathcal{B}_{1}e^{3i\phi}$
$\displaystyle=\frac{1}{h_{s}}\left[A^{2}k^{2}e^{-2i(\theta-\phi)}+2AA_{1}k^{4}e^{2i(\theta+\phi)}\right.$
$\displaystyle+\left.A_{1}^{2}k^{6}e^{2i(3\theta+\phi)}-\eta^{2}k_{z}^{2}k^{4}e^{-4i(\theta-\phi)}\right]+\tilde{h}_{p},$
(134)

the magnitude of which determines $k_{c}$ when
$\theta=\theta_{c\sigma}=\phi+\sigma\pi/2$. Thus, $\mathcal{B}_{1}$ evaluated
at $\theta=\theta_{c\sigma}$ is

$\displaystyle\mathcal{B}_{1}e^{3i\phi}$
$\displaystyle=\frac{-1}{h_{s}}\left[A^{2}k^{2}+2AA_{1}k^{4}\cos
4\phi+A_{1}^{2}k^{6}\cos 8\phi-\eta^{2}k_{z}^{2}k^{4}\right]$
$\displaystyle+\tilde{h}_{p}-i\frac{2A_{1}k^{4}}{h_{s}}\left(A+A_{1}k^{2}\cos
4\phi\right)\sin 4\phi$ (135)

In general, the equation for $k_{c}$ using Eq. (D) is a 6th order polynomial
for $k_{c}^{2}$. However, as we limit our Hamiltonian Eq. (129) to
$\mathcal{O}(k^{3})$, we should only keep up to $\mathcal{O}(k^{4})$ terms in
$\mathcal{B}_{1}$, and the resulting equation for $k_{c}$ is determined by

$\displaystyle\left[A^{2}k_{c}^{2}+2AA_{1}k_{c}^{4}\cos
4\phi-\eta^{2}k_{z}^{2}k_{c}^{4}-h_{s}\tilde{h}_{p}\right]^{2}$
$\displaystyle+4A_{1}^{2}A^{2}k_{c}^{8}\sin^{2}4\phi=0.$ (136)

which is a quartic equation for $k_{c}^{2}$ and thus has a general solution.
Due to the smallness of $\tilde{h}_{p}$ compared to $h_{s}$, however, $k_{c}$
is well approximated by $k_{c0}$. We can expand $\mathcal{B}_{1}$ up to lowest
order in $k_{c}$, giving us

$\displaystyle\mathcal{B}_{1}$
$\displaystyle\approx-\frac{2iA^{2}k_{c}}{h_{s}}\delta k_{-}e^{-2i\phi}.$
(137)

The general form of $\mathcal{B}_{2}$ after expressing $\mathbf{k}$ in polar
coordinates is

$\displaystyle\mathcal{B}_{2}$
$\displaystyle=-\frac{\tilde{A}k^{3}k_{z}}{h_{s}}\left(C_{-}\cos(\theta-\phi)+C_{+}\cos(3\theta+\phi)\right)-\frac{k^{5}k_{z}}{4h_{s}}\left[\left((C_{3a}+C_{3b})C_{3e}+4C_{3b}C_{3d}\right)\cos(\theta-\phi)\right.$
$\displaystyle\left.-(C_{3a}-3C_{3b})C_{+}\cos(3\theta+\phi)+(C_{3a}+C_{3b})C_{-}\cos(5\theta-\phi)\right]$
(138)

with $C_{\pm}=C_{3d}\pm C_{3e}/4$. Substituting in the relation between
coefficients $C_{i}$ used in the main text, Eq. (D) reduces to Eq. (18b). At
$\theta=\phi\pm\pi/2$, we have

$\displaystyle\mathcal{B}_{2}=\pm
m(\phi)=\pm\left[-\frac{\tilde{A}k_{c}^{3}k_{z}}{h_{s}}C_{+}\pm\frac{k^{5}_{c}k_{z}}{4h_{s}}\left((C_{3a}-3C_{3b})C_{+}+(C_{3a}+C_{3b})C_{-}\right)\right]\sin(4\phi)$
(139)

The mass thus vanishes if the field points perpendicular to a mirror plane,
since this mirror then remains a symmetry of the system. Combining time-
reversal with the mirror that contains the field direction is also a symmetry.

##  Appendix E Surface States Dispersion Calculation

In this Appendix, we provide details on the calculation of the helical edge
modes with and without an in-plane magnetic field. The boundary condition for
the surface state is expressed in terms of the eigenvectors in Eq. (27b)

$\displaystyle\psi_{1\alpha}(\bm{k}_{x},\bm{k}_{y})$
$\displaystyle=-\alpha\mathcal{N}_{\alpha+}(k_{x},k_{y})\frac{E_{\alpha+}(k_{x},k_{y})+M(k_{x},k_{y})}{Ak_{+}}$
(140a) $\displaystyle\psi_{2\alpha}(\bm{k}_{x},\bm{k}_{y})$
$\displaystyle=-\alpha\mathcal{N}_{\alpha+}e^{-i\phi}\frac{k_{-}}{k_{+}}$
(140b) $\displaystyle\psi_{3\alpha}(\bm{k}_{x},\bm{k}_{y})$
$\displaystyle=-\mathcal{N}_{\alpha+}(k_{x},k_{y})\frac{E_{\alpha+}(k_{x},k_{y})+M(k_{x},k_{y})}{Ak_{+}}$
(140c) $\displaystyle\psi_{4\alpha}(\bm{k}_{x},\bm{k}_{y})$
$\displaystyle=\mathcal{N}_{\alpha+}(k_{x},k_{y})$ (140d)

the matrix $\mathcal{A}$ whose determinant being zero gives the surface state
dispersion $\epsilon(k_{x})$ is then given by

$\displaystyle\mathcal{A}$
$\displaystyle=\left(\begin{array}[]{cccc}\psi_{1+}(k_{x},-i\lambda_{++})&\psi_{1+}(k_{x},-i\lambda_{+-})&\psi_{1-}(k_{x},-i\lambda_{-+})&\psi_{1-}(k_{x},-i\lambda_{--})\\\
\psi_{2+}(k_{x},-i\lambda_{++})&\psi_{2+}(k_{x},-i\lambda_{+-})&\psi_{2-}(k_{x},-i\lambda_{-+})&\psi_{2-}(k_{x},-i\lambda_{--})\\\
\psi_{3+}(k_{x},-i\lambda_{++})&\psi_{3+}(k_{x},-i\lambda_{+-})&\psi_{3-}(k_{x},-i\lambda_{-+})&\psi_{3-}(k_{x},-i\lambda_{--})\\\
\psi_{4+}(k_{x},-i\lambda_{++})&\psi_{4+}(k_{x},-i\lambda_{+-})&\psi_{4-}(k_{x},-i\lambda_{-+})&\psi_{4-}(k_{x},-i\lambda_{--})\end{array}\right)$
(145)

We then solve $\det\mathcal{A}=0$ numerically to obtain the surface state
dispersion $\epsilon(k_{x})$. We use the parameters from Ref. [29],

$\displaystyle A$ $\displaystyle=111.6\,\text{meV}\,\text{nm}$ (146a)
$\displaystyle M_{0}$ $\displaystyle=28.2\,\text{meV}$ (146b) $\displaystyle
M_{z}$ $\displaystyle=207.2\,\text{meV}\,\text{nm}^{2}$ (146c) $\displaystyle
M_{xy}$ $\displaystyle=133.2\,\text{meV}\,\text{nm}^{2}$ (146d)

and we use $n=2$, $L=19$ nm to model the Cd3As2 thin film setup reported in
Refs. [25, 26].

##  Appendix F Particle-Hole Asymmetry in Landau Level Spectrum

Figure 8:  The $\nu=0-4$ Landau levels plotted as a function of the inverse
magnetic length squared $eB$. The $\nu=0$ levels are shown with thicker
strokes, with $\nu\neq 0$ particle and hole levels colored orange and blue,
respectively. The inclusion of the band-diagonal $\epsilon_{n}(\mathbf{k})$
breaks particle-hole symmetry by introducing $\omega_{2}$, which provides a
positive dispersion in magnetic field for every Landau level. The values of
model coefficients used are from Ref. [29] [see also Eq. (146)].

The inclusion of the band-diagonal term $\epsilon_{n}(\mathbf{k})$ [see Eq.
(6a)] causes an additional contribution to Eq. (34), and the full Hamiltonian
is given by

$\displaystyle H$
$\displaystyle=C_{n}+\omega_{2}\left(a^{\dagger}a+\frac{1}{2}\right)+\left[M_{n}-\omega\left(a^{\dagger}a+\frac{1}{2}\right)\right]\tau_{z}s_{0}$
$\displaystyle+\frac{A}{\sqrt{2}l_{B}}\left[\tau_{x}s_{z}\left(a+a^{\dagger}\right)+i\tau_{y}s_{0}\left(a-a^{\dagger}\right)\right]$
(147)

where $C_{n}=C_{0}+C_{1}\left(\frac{n\pi}{L}\right)^{2}$ and
$\omega_{2}=\frac{2C_{2}}{l_{B}^{2}}$. The $\nu\neq 0$ Landau levels are found
with the trial wavefunctions

$\displaystyle\psi_{\nu+}$
$\displaystyle=\left(0,0,B_{1+}\ket{\nu},B_{2+}\ket{\nu-1}\right)^{T}$ (148)
$\displaystyle\psi_{\nu-}$
$\displaystyle=\left(B_{1-}\ket{\nu-1},B_{2-}\ket{\nu},0,0\right)^{T}$ (149)

giving the $\nu\neq 0$ Landau level energies

$\displaystyle\epsilon_{\alpha\beta}(\nu)$
$\displaystyle=C_{n}-\alpha\frac{\omega}{2}+\nu\omega_{2}$
$\displaystyle\quad+\beta\sqrt{\frac{2A^{2}}{l_{B}^{2}}\nu+\left(\alpha\left(M_{n}-\nu\omega\right)+\frac{\omega_{2}}{2}\right)^{2}}.$
(150)

The inclusion of the band-diagonal term causes two contributions to the Landau
level spectrum: first, a trivial shift of the overall energies by $C_{n}$, but
more importantly, $\omega_{2}$ contributes a positive dispersion to every
Landau level as a function of magnetic field, introducing a particle-hole
asymmetry as shown in Fig. 8

## References

  * Kane and Mele [2005] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005).
  * Bernevig  _et al._ [2006] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006).
  * Bernevig and Zhang [2006] B. A. Bernevig and S.-C. Zhang, Phys. Rev. Lett. 96, 106802 (2006).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Maciejko  _et al._ [2011] J. Maciejko, T. L. Hughes, and S.-C. Zhang, Annual Review of Condensed Matter Physics 2, 31 (2011).
  * Qian  _et al._ [2014] X. Qian, J. Liu, L. Fu, and J. Li, Science 346, 1344 (2014).
  * Cao and Chen [2019] C. Cao and J.-H. Chen, Advanced Quantum Technologies 2, 1900026 (2019).
  * Wang  _et al._ [2019] D. Wang, F. Tang, J. Ji, W. Zhang, A. Vishwanath, H. C. Po, and X. Wan, Phys. Rev. B 100, 195108 (2019).
  * Fei  _et al._ [2017] Z. Fei, T. Palomaki, S. Wu, W. Zhao, X. Cai, B. Sun, P. Nguyen, J. Finney, X. Xu, and D. H. Cobden, Nature Phys 13, 677 (2017).
  * Konig  _et al._ [2007] M. Konig, S. Wiedmann, C. Brune, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Science 318, 766 (2007).
  * Roth  _et al._ [2009] A. Roth, C. Brune, H. Buhmann, L. W. Molenkamp, J. Maciejko, X.-L. Qi, and S.-C. Zhang, Science 325, 294 (2009).
  * Du  _et al._ [2015] L. Du, I. Knez, G. Sullivan, and R.-R. Du, Phys. Rev. Lett. 114, 096802 (2015).
  * Reis  _et al._ [2017] F. Reis, G. Li, L. Dudy, M. Bauernfeind, S. Glass, W. Hanke, R. Thomale, J. Schäfer, and R. Claessen, Science 357, 287 (2017).
  * Bampoulis  _et al._ [2023] P. Bampoulis, C. Castenmiller, D. J. Klaassen, J. van Mil, Y. Liu, C.-C. Liu, Y. Yao, M. Ezawa, A. N. Rudenko, and H. J. W. Zandvliet, Phys. Rev. Lett. 130, 196401 (2023).
  * Tang  _et al._ [2017] S. Tang, C. Zhang, D. Wong, Z. Pedramrazi, H.-Z. Tsai, C. Jia, B. Moritz, M. Claassen, H. Ryu, S. Kahn, J. Jiang, H. Yan, M. Hashimoto, D. Lu, R. G. Moore, C.-C. Hwang, C. Hwang, Z. Hussain, Y. Chen, M. M. Ugeda, Z. Liu, X. Xie, T. P. Devereaux, M. F. Crommie, S.-K. Mo, and Z.-X. Shen, Nature Physics 13, 683 (2017).
  * Wu  _et al._ [2018] S. Wu, V. Fatemi, Q. D. Gibson, K. Watanabe, T. Taniguchi, R. J. Cava, and P. Jarillo-Herrero, Science 359, 76 (2018).
  * Linder  _et al._ [2009] J. Linder, T. Yokoyama, and A. Sudbø, Phys. Rev. B 80, 205401 (2009).
  * Lu  _et al._ [2010] H.-Z. Lu, W.-Y. Shan, W. Yao, Q. Niu, and S.-Q. Shen, Physical Review B 81, 115407 (2010), publisher: American Physical Society.
  * Liu  _et al._ [2010] C.-X. Liu, H. Zhang, B. Yan, X.-L. Qi, T. Frauenheim, X. Dai, Z. Fang, and S.-C. Zhang, Physical Review B 81, 041307 (2010), publisher: American Physical Society.
  * Zyuzin and Burkov [2011] A. A. Zyuzin and A. A. Burkov, Phys. Rev. B 83, 195413 (2011).
  * Zhang  _et al._ [2010] Y. Zhang, K. He, C.-Z. Chang, C.-L. Song, L.-L. Wang, X. Chen, J.-F. Jia, Z. Fang, X. Dai, W.-Y. Shan, S.-Q. Shen, Q. Niu, X.-L. Qi, S.-C. Zhang, X.-C. Ma, and Q.-K. Xue, Nature Phys 6, 584 (2010).
  * Burkov [2018] A. A. Burkov, Phys. Rev. B 97, 165104 (2018).
  * Akzyanov [2019] R. S. Akzyanov, Phys. Rev. B 100, 045403 (2019).
  * Chorsi  _et al._ [2020] H. T. Chorsi, S. Yue, P. P. Iyer, M. Goyal, T. Schumann, S. Stemmer, B. Liao, and J. A. Schuller, Advanced Optical Materials 8, 1901192 (2020).
  * Lygo  _et al._ [2023] A. C. Lygo, B. Guo, A. Rashidi, V. Huang, P. Cuadros-Romero, and S. Stemmer, Phys. Rev. Lett. 130, 046201 (2023).
  * Guo  _et al._ [2023] B. Guo, W. Miao, V. Huang, A. C. Lygo, X. Dai, and S. Stemmer, Phys. Rev. Lett. 131, 046601 (2023).
  * Wang  _et al._ [2013] Z. Wang, H. Weng, Q. Wu, X. Dai, and Z. Fang, Physical Review B 88, 125427 (2013).
  * Miao  _et al._ [2023] W. Miao, B. Guo, S. Stemmer, and X. Dai, arXiv preprint arXiv:2309.15457 (2023).
  * Villar Arribi  _et al._ [2020] P. Villar Arribi, J.-X. Zhu, T. Schumann, S. Stemmer, A. A. Burkov, and O. Heinonen, Phys. Rev. B 102, 155141 (2020).
  * Guo  _et al._ [2022] B. Guo, A. C. Lygo, X. Dai, and S. Stemmer, APL Materials 10, 091116 (2022).
  * Zhou  _et al._ [2008] B. Zhou, H.-Z. Lu, R.-L. Chu, S.-Q. Shen, and Q. Niu, Phys. Rev. Lett. 101, 246807 (2008).
  * McCann and Fal’ko [2006] E. McCann and V. I. Fal’ko, Phys. Rev. Lett. 96, 086805 (2006).
  * Zhang  _et al._ [2023] Z. Zhang, Z.-M. Yu, G.-B. Liu, Z. Li, S. A. Yang, and Y. Yao, Computer Physics Communications 290, 108784 (2023).
  * Cohen-Tannoudji  _et al._ [2019] C. Cohen-Tannoudji, B. Diu, and F. Laloë,  _Quantum Mechanics, Volume 2: Angular Momentum, Spin, and Approximation Methods_ (Wiley, 2019).
  * Lindgren [1974] I. Lindgren, Journal of Physics B: Atomic and Molecular Physics 7, 2441 (1974).