# Characterisation of plastic scintillator paddles and
lightweight MWPCs for the MID subsystem of ALICE 3

Ruben Alfaro  Mauricio Alvarado Hernández  Gyula Bencédi  Juan Carlos
Cabanillas Noris  Marco Antonio Díaz Maldonado  Carlos Duarte Galvan  Arturo
Fernández Téllez  Gergely Gábor Barnaföldi  Ádám Gera  Varlen Grabsky  Gergő
Hamar  Gerardo Herrera Corral  Ildefonso León Monzón  Josué Martínez García
Mario Iván Martínez Hernandez  Jesús Eduardo Muñoz Méndez  Richárd Nagy
Rafael Ángel Narcio Laveaga  Antonio Ortiz  Mario Rodríguez-Cahuantzi
Solangel Rojas Torres  Timea Szollosova  Miguel Enrique Patiño Salazar  Jared
Pazarán García  Hector David Regules Medel  Guillermo Tejeda Muñoz  Paola
Vargas Torres  Dezső Varga  Róbert Vértesi  Yael Antonio Vasquez Beltran
Carlos Rafael Vázquez Villamar  Irandheny Yoval Pozos  Instituto de Física,
Universidad Nacional Autónoma de México, Mexico  Instituto de Ciencias
Nucleares, Universidad Nacional Autónoma de México, Mexico  HUN-REN Wigner
Research Centre for Physics, Hungary  Instituto Tecnológico de Culiacán,
Tecnológico Nacional de México  Universidad Autónoma de Sinaloa, Mexico
Facultad de Ciencias Físico Matemáticas, Benemérita Universidad Autónoma de
Puebla, Mexico  Centro de Investigación y Estudios Avanzados del IPN, Mexico
Czech Technical University in Prague, Prague, Czech Republic

###### Abstract

The ALICE collaboration is proposing a completely new detector, ALICE 3, for
operation during the LHC Runs 5 and 6. One the ALICE 3 subsystems is the Muon
IDentifier detector (MID), which has to be optimised to be efficient for the
reconstruction of $J/\psi$ at rest (muons down to $p_{\rm T}\approx 1.5$
GeV/$c$) for $|\eta|<1.3$. Given the modest particle fluence expected in the
MID of a few Hz/cm2, technologies like plastic scintillator bars ($\approx 1$
m length) equipped with wavelength-shifting fibers and silicon photomultiplier
readout, and lightweight Multi-Wire Proportional Chambers (MWPCs) are under
investigation. To this end, different plastic scintillator paddles and MWPCs
were studied at the CERN T10 test beam facility. This paper reports on the
performance of the scintillator prototypes tested at different beam energies
(from 0.5 GeV/$c$ up to 6 GeV/$c$) and positions (horizontal, vertical, and
angular scans). The MWPCs were tested at different energies (from 0.5 GeV/$c$
to 10 GeV/$c$) and beam intensities, their efficiency and position resolutions
were verified beyond the particle rates expected with the MID in ALICE 3.

###### keywords:

Upgrade, Scintillators, SiPM, MWPC, ALICE experiment

††journal: JINST

##  1 Introduction

ALICE (A Large Ion Collider Experiment) is one of the four main experiments at
the CERN Large Hadron Collider (LHC) and it is mainly devoted to heavy-ion
physics [1]. Its main goal is to achieve the detailed characterisation of the
strongly interacting quark–gluon plasma (QGP) that is theoretically predicted
in lattice QCD, and that has been observed in relativistic heavy-ion
collisions by experiments at RHIC and at the LHC [2, 3].

In spite of the ambitious upgrade and physics programs to be completed in the
next eight years, crucial questions will remain open after the ALICE data
taking in the LHC Runs 3 and 4. These include the deeper understanding of the
rich phenomenology of QCD matter by connecting parton-energy loss, collective
flow, hadronisation and electromagnetic radiation in a unified description.
The successful completion of this task requires the construction of a new
detector with increased rate capabilities, improved vertexing and tracking
over wide momentum and rapidity ranges. This project, named ALICE 3, is being
optimized to have excellent capabilities for heavy-flavour and electromagnetic
radiation studies [4]. The ALICE 3 detector is expected to take data during
the LHC Run 5 (2035-2038) and Run 6 (2040-2041). The core of ALICE 3 will be a
full silicon tracking device covering eight units of pseudorapidity, a Time-
Of-Flight (TOF) system consisting of two cylindrical layers with a time
resolution of 20 ps, a Ring Imaging Cherenkov (RICH) detector that will allow
the ${\rm e}/\pi$ separation up to 2 GeV/$c$, and finally an electromagnetic
calorimeter covering the barrel acceptance and one forward direction. All
these detectors will be located inside a superconducting magnet, providing a
magnetic field of $B=2$ T. A Forward Conversion Tracker (FCT) is under
consideration to measure photons down to transverse momenta of $p_{\rm
T}\approx 2$ MeV/$c$. For muon identification, the Muon IDentifier (MID)
detector will consist of a cylindrical shape steel absorber of thickness
$\approx 70$ cm, surrounding the magnet, and of MID chambers around the
absorber allowing the muon tagging by matching hits in the muon detectors with
tracks from the tracker. The absorber thickness is being optimized for the
efficient reconstruction of $J/\psi$ emitted at rest at $\eta=0$, which
requires the muon identification with a transverse momentum as low as $p_{\rm
T}=1.5$ GeV/$c$.

Regarding the MID subsystem, muon chambers with a granularity of
$\Delta\varphi\Delta\eta=0.02\times 0.02$ are proposed [4]. The moderate
charged-particle rate in the MID subsystem of a few Hz/cm2 for pp and Pb–Pb
collisions at an interaction rate of 24 MHz and 93 kHz, respectively, together
with the required granularity of 50-60 mm pad size open the possibility of
using well known technologies like Multi-Wire Proportional Chambers (MWPCs)
and Resistive Plate Chambers (RPCs). However, the use of scintillator paddles,
i.e. scintillator bars equipped with wave-length shifting fibers (WLS) and
silicon photomultiplier (SiPM) readout, are considered as the baseline option
in the ALICE 3 LoI due to their simplicity [4]. This option considers two
layers of crossed scintillator paddles (5 cm wide) with a 10 cm gap between
two layers. Both MWPC and scintillators represent a cost effective solution to
cover a detection area of about $360$ m2.

In this paper, the performance of different options of plastic scintillators
(FNAL-NICADD [5], ELJEN model EJ-208 [6] and U-70 IHEP, PROTVINO) and
photosensors are reported. For MWPCs, two sets of chambers from Wigner RCP
($512\times 512\,{\rm mm}^{2}$, wires: 8mm, pad: 4 mm) and ($768\times
768\,{\rm mm}^{2}$, wires: 12 mm, pad: 12 mm) are tested with high and
moderate position resolution, measuring their performance at beam conditions,
as well as quantifying the performance in high-particle rate environment.

##  2 Experimental setup

The measurements were performed at the T10 beamline at the East Area of the
CERN PS, which underwent a recent overhaul [7] allowing a convenient beam
testing activity. Three scintillator prototypes and two sets of MWPCs were
tested. Figure 1 shows a sketch of the experimental setup and the coordinate
system. The beam particles travel along the $z$ axis, the horizontal axis
mentioned hereinafter corresponds to the $x$ axis, whereas the vertical axis
corresponds to the $y$ axis. The trigger signal was provided by the
coincidence among three scintillator paddles denoted as TRIG. $X1$, TRIG. $Y1$
and TRIG. $X2$. The scintillators are coupled to Hamamatsu R3478
photomultipliers, resulting in high efficiency and high signal-to-noise ratio
with a time resolution of about 1 ns.

The MWPCs operated in a self-trigger mode, which, despite its limited time
resolution, facilitated the detection of all particles from the beam due to
their expansive coverage area. The detected particles include those that
missed the trigger scintillators due to their small detection area ($\approx
3\times 2$ cm2). Thus the real beam profile could be measured.

Figure 1: Sketch of the experimental setup. The trigger detectors are labeled
as TRIG. $X1$, TRIG. $Y1$ and TRIG. $X2$. The beam travels from right to left
along the $z$ coordinate that is perpendicular to the array of detectors
tested at T10. The detector areas are given in terms of the length along the
horizontal axis ($x$) times the length along the vertical axis ($y$).

###  2.1 Plastic scintillator prototypes

The low-cost extruded plastic scintillator from FNAL-NICADD [5] was included
among the prototypes. Other prototypes based on the relatively low-cost
scintillator manufactured by ELJEN (model EJ-208) and plastic scintillator
manufactured at U-70 IHEP, PROTVINO (which was used for the cosmic-ray
detector of ALICE [Fernández:957722]) were also tested. All of them had
approximately one meter length. The FNAL-NICADD and ELJEN EJ-208 prototypes
were constructed and tested with cosmic muons at UNAM, and the PROTVINO
prototype was constructed and tested at BUAP. Table 1 summarizes the main
features of the scintillation bars as well as the model of the SiPMs that were
used in the prototypes. A brief description of the plastic scintillator
prototypes can be found below.

Table 1: Characteristics of the plastic scintillator bars. The width (w), length (l) and thickness (t) are specified, as well as the number of fibers and SiPM model. Company | Bar dimensions | No. of fibers | Fiber diameter | SiPM model
---|---|---|---|---
| (${\rm w}\times{\rm l}\times{\rm t}$) |  |  |
FNAL-NICADD | $4.0\times 100.0\times 1.0\,{\rm cm}^{3}$ | 1 | 1.5 mm | Hamamatsu S13360-6050CS
ELJEN EJ-208 | $5.0\times 100.0\times 1.0\,{\rm cm}^{3}$ | 0 | – | Hamamatsu S13360-6050CS
PROTVINO | $5.0\times 100.0\times 1.0\,{\rm cm}^{3}$ | 1 | 1.5 mm | SensL serie C $6\times 6\,{\rm mm}^{2}$

Given the relatively large area to be covered, an attractive alternative is
the low-cost scintillator produced at the FNAL-NICADD facility [5], which has
been used in several experiments (see e.g. Ref. [8, 9]). This plastic
scintillator has a light emission peak at 420 nm [10]. Since the attenuation
length of the plastic scintillator is short (a few tens of cm), the light
produced by the particle interaction has to be collected, reemitted, and
transported to the photodetectors by wave-length shifting (WLS) fibers. The
prototype used the Kuraray WLS fiber model Y-11(200) that poses a long-
attenuation length ($\approx 400$ cm) with an emission (absorption) peak at
476 nm (375 nm) [11]. The fiber is mounted inside the $2$ mm groove on the
large surface of the scintillator coupled only with air. A Hamamatsu SiPM of
type S13360-6050CS was coupled at one side of the fiber mentioned above. The
SiPM has a photosensitive area of $6.0\times 6.0$ mm2, pixel pitch of 50
$\mu$m, its detection spectral response is found in the wavelength interval
270-900 nm, with a photon detection efficiency of 40% at peak sensitivity
wavelength (450 nm). The SiPM was connected to the module SP5600 manufactured
by CAEN that serves to provide the bias voltage (57 V) and an analog output
connector.

Another alternative that was tested is the relatively low-cost plastic
scintillator produced by ELJEN (EJ-208). This scintillator has a long
attenuation length of about 4 m and therefore no WLS fibers were considered in
this prototype. This plastic scintillator has a light emission peak at 435 nm
with rise (decay) time of 1 ns (3.3 ns). The coupling with the SiPM
(S13360-6050CS) was done with an optical grease (ELJEN) having high
transparency and a refractive index close to the refractive index of the
plastic scintillator ($n=1.58$). The CAEN module SP5600 was also used in this
prototype that serves to provide the bias voltage (57 V) and an analog output
connector.

A plastic scintillator from PROTVINO was also explored in order to see whether
its performance fulfill the requirements of MID. The bar was equipped with two
SiPMs at each side of the bar and one Kuraray WLS fiber model Y-11(200). Two
SiPMs were needed otherwise the efficiency at large distance with respect to
SiPM was observed to drop down to 80%. Data for beam momentum of 6 GeV/c were
used. The SiPMs used in this prototype were produced by SensL (C-series) with
an active area of $6.0\times 6.0$ mm2 operated at 30 V. The fiber was directly
connected to the SiPMs without any optical grease. Its photon detection
efficiency is about 40% at peak sensitivity wavelength ($\approx 420$ nm). A
custom-made front-end electronics card was used.

Data acquisition was done using standard VME and NIM modules. Three
scintillation paddles indicated as TRIG. $X1$, TRIG. $Y1$, and TRIG. $X2$ in
Fig. 1 were used for triggering, their output signals were connected to a NIM
leading edge discriminator (CAEN N840). Two of the signals were processed with
a coincidence module (LeCroy 622A) to define the two-fold coincidence. The
two-fold coincidence was defined by TRIG. $X1$ $\wedge$ TRIG. $Y1$, while the
three-fold one was defined by the coincidence among TRIG. $X1$ $\wedge$ TRIG.
$Y1$ $\wedge$ TRIG. $X2$.

The two-fold coincidence signal served as the primary trigger to initiate data
acquisition, while the three-fold coincidence was digitised for offline
analysis. The signals from the SiPMs were fully digitised using the CAEN V1742
digitizer with 32 channels, enabling the simultaneous recording of signals
from all prototypes. The digitizer was configured with a sampling rate of 1
GS/$s$. A TTI PL330DP power supply provided the voltage (30 V) to SensL SiPMs.
The output signals from these SiPMs were amplified tenfold before being
connected to the digitizer using a variable fast amplifier (CAEN V974). The
data of the 32 channels were transferred to a computer with the VME bridge
(CAEN V1718) for further processing and analysis using standard CAEN software.

###  2.2 MWPC prototypes

The MID subsystem requires high efficiency ($\approx 100\%$) and lightweight
detectors. The so-called “CCC” (closed-cathode chamber [12]) MWPC used for
similar purposes in the VHMPID upgrade project of ALICE [13] has been
investigated in the past. The purpose of the present beam test was to verify
the compatibility of similar MWPC structures with the MID requirements. For
this reason, a more standard but still lightweight and low-cost MWPC using
contemporary materials and design concept [14] was built. Similar MWPC
detectors are broadly applied for cosmic muon imaging, where high reliability
and robust structures are required [15].

Eight MWPC detectors spaced over more than 3 m length, defining very precise
individual beam particles, were tested. Two detector sizes were used, labeled
as “50 cm” and “80 cm” prototypes. The detectors have a gas gap of 18–22 mm,
an anode wire spacing of 8 mm for “50 cm” and 12 mm for “80 cm” prototypes,
with two-dimensional projective readout achieved by field wires parallel to
the anode wires, as well as perpendicular “pick-up” wires [15]. Having 8 MWPCs
provided sufficiently precise and redundant track information. A “track” is
defined if at least 7 out of the 8 available detectors form a straight line
within the position resolution of the detectors. Individual detector
efficiencies are typically above 95%, which ensures high efficiency for track
identification.

The detector system was controlled and data were read out using a DAQ system
designed for cosmic muons [14], however, it proved to be sufficiently fast to
operate up to the highest relevant beam intensity. The dead time of the DAQ
system was about 0.1 ms, leading to 10 kHz maximum event rate for recording.

The beam profile at any plane perpendicular to the beam direction was
determined extrapolating the individual tracks to the beam source. The beam
profile was found to be strongly energy and position dependent, and slightly
intensity dependent. Figure 2 shows an example for the $x$ position of the
first MWPC. The left panel shows the beam profile (number of particles along
the $x$ axis) with triggering only on the MWPC detectors, which is the “real”
beam profile reaching the detectors. The right-hand side panel shows that the
two-fold coincidence signal used for efficiency measurements gives a much
cleaner beam definition. Apparently the beam widens at lower energies as
expected [7], but it is conveniently usable for the measurements over the full
energy range. The triggered beam profile is relevant for efficiency
measurements, however the complete beam profile must be considered when the
beam intensity is determined at any specific detector location.

Figure 2: Horizontal beam profiles (number of particles along the $x$ axis)
measured with MWPCs at various beam momenta. The beam profile is shown
considering the self-triggered MWPC mode (left panel) and triggered mode using
the scintillation paddles (right panel).

##  3 Results

###  3.1 Plastic scintillator bars

To illustrate the quality of data, the charge ($Q$) distributions obtained
with FNAL-NICADD prototype are shown in Fig. 3 for different beam momenta as
well as for different distances between the hit and the SiPM. The charge was
obtained integrating the pulse within a time gate of 300 ns. The main peak of
the signal exhibits no or little variation with respect to the distance, due
to the WLS fiber. A model based on a convolution of a Landau and a Gaussian
distribution is fitted to the charge distributions. The mean of the Gaussian
and the most-probable value (MPV) of the the Landau function are constrained
to be the same. The Gaussian distribution serves to describe detector effects
which smear the charge distribution. The measured charge (MPV) is around 18 pC
that correspond to $\approx 40$ photoelectrons. Figure 3 shows a clear and
consistent signal-to-pedestal separation as a function of hit position and
beam momentum. A similar good pedestal/signal separation is seen for ELJEN
EJ-208 and PROTVINO prototypes. Although, ELJEN EJ-208 scintillator prototype
gives a light output of up to 160 photoelectrons, the main disadvantage of
this plastic scintillator is the cost that is 4 times that of FNAL-NICADD.
Last but not least, the variation of charge (MPV) with time gate is found to
be below 10%, and it is stable with the beam momentum.

Figure 3: Charge distributions measured for the FNAL-NICADD prototype. Results
are displayed for different beam momenta, as well as for three different
distances between the hit and the SiPM of $\approx 5$ cm (left), $\approx 50$
cm (middle), and $\approx 95$ cm (right).

Light collection is expected to depend on the mean path length ($\langle
L\rangle$) traveled by the particle in the scintillator. The data presented in
this paper consider hadrons crossing the scintillator bars at normal incidence
(incidence angle of $\sim 90^{\circ}$), therefore they traverse the 1 cm
thickness of the bar. In experiments at colliders like LHC, only particles
coming from the interaction point at pseudorapidity equal to zero ($\eta=0$)
would satisfy this condition. However, the acceptance of the MID subsystem is
expected to be around $|\eta|<1.24$, therefore particles at $\eta=1.24$
(incidence angle of $\sim 32^{\circ}$) might reach the detector. Particles
reaching MID at smaller incidence angles ($0-30^{\circ}$) originated from
beam-induced background could also be observed. It is therefore useful to
study the detector response as a function of the incidence angle. This study
was conducted by rotating the prototypes covering the interval between
$10^{\circ}$ and $60^{\circ}$. Figure 4 shows the charge (MPV) as a function
of the incidence angle and path length. As expected, the charge value shows a
modest variation ($\sim 1.7$) for incidence angles within
$30^{\circ}-60^{\circ}$, which is roughly the interval covered by particles
from the interaction point. For smaller angles, the charge exhibits a steep
increase with decreasing the incidence angle. For angles of about $\sim
10^{\circ}$, the charge is four times that achieved at normal incidence.

Figure 4: Charge (MPV) as a function of the incidence angle and path length
($\langle L\rangle$). Results for FNAL-NICADD and ELJEN EJ-208 are displayed.
The statistical errors are smaller than the markers, and therefore, are not
visible.

The efficiency is determined considering events that are triggered by the
triple coincidence TRIG. $X1$ $\wedge$ TRIG. $Y1$ $\wedge$ TRIG. $X2$. It is
obtained as the ratio of the number of events registered by the prototype
(measured charge above the threshold) divided by the total number of events.
In the case of PROTVINO prototype, the detector signal comes from any of the
two SiPMs that were placed at each side of the bar. The efficiencies as a
function of the hit position along the bar length are shown in Fig. 5. Results
for hadron momentum of 6 GeV/$c$ are presented. The efficiencies are above 98%
along the length of the bar and the detector efficiency is found to be quite
uniform. For completeness, a similar measurement was performed along the bar
width (vertical scan) at mid position (hit position at 50 cm from the SiPM).
Given the size of the active area of the trigger, $2\times 2$ cm2, the
efficiency measurement was restricted to $\pm 1.00$ cm ($\pm 1.5$ cm) around
the center of the bar for FNAL-NICADD (ELJEN EJ-208) prototype. The right
panel of Fig. 5 shows that within those intervals the efficiency is above 98%
for both prototypes.

Figure 5: Left: efficiency as a function of the distance between the SiPM and
the hit position (horizontal scan). Right: efficiency along the width of the
bar (vertical scan). The statistical errors are smaller than the markers, and
therefore, are not visible.

The time resolution of scintillator prototypes is mostly determined by the
length of the scintillator bar ($\approx 1$ m) and is about 1.7 ns and 1.5 ns
for FNAL-NICADD and ELJEN EJ-208 prototypes, respectively. This time
resolution is already enough to operate the detector in the expected LHC
conditions during Runs 5 and 6. In the case of ALICE 3 the expected
interaction rate for pp collisions is of around 24 MHz. Assuming an
integration time of the ALICE 3 tracker of about 100 ns [4], a pileup of
around 2-4 pp collisions is expected. The measured time resolution in our
prototypes would be more than enough for the MID subsystem.

###  3.2 MWPC

Similar “CCC-type” detectors were thoroughly characterised in test beams
earlier [12] and field measurements for muon imaging [16]. Detailed studies on
MWPC with cosmic muons have been reported [14], therefore this paper primarily
addresses the intensity dependence of the MWPC performance.

The beam test conditions allowed the measurements to span over a broad range
of beam intensities, by adjusting the collimator settings, and by changing the
position along the diverging beamline. Particularly, beam rates two orders of
magnitude above the expected rate in the ALICE 3 MID could be safely reached
during the tests.

Intensity is defined as particle rate divided by the beam spot area. The beam
area could be determined assuming a Gaussian beam profile,
$A=\pi\sigma_{x}\sigma_{y}$ where $\sigma$ is the width of the beam spot. The
beam profile is shown in Fig. 2 at the respective position of the detector
under study. The particle rate was derived directly from the count rate
measured using default beam-instruments of T10 and our DAQ. The estimation of
the intensity has a considerable systematic error, around 30% as averaging on
a mostly-Gaussian distribution, to which the main contribution is from the
non-Gaussian shape of the beam.

Figure 6 shows the measured efficiency as a function of the estimated beam
intensity. The efficiency is defined as the probability to find a valid hit in
a given chamber, under the condition that using the other 7 chambers, a
straight line track can be formed using at least 6 chambers.

Figure 6: MWPC detection efficiency measured as a function of beam intensity
for beam momentum of 5 GeV/$c$.

The measured efficiency is around 98% at beam intensity of 10 Hz/cm2, and
decreases only slightly until a steep drop around 1 kHz/cm2. Since the top
beam intensities were saturating the DAQ with total beam rate well above 10
kHz/cm2, it is possible that the efficiency drop is connected to the readout
electronics, rather than the detector amplification structure.

The key conclusion is that the efficiency is above 95% and constant for
intensities up to 100 Hz/cm2, well above the intensities expected at the ALICE
3 MID.

The position information from the MWPC chambers is acquired from wires running
parallel and perpendicular to the anode wires, called “pads” and “field”
wires, respectively [15]. In order to simplify the readout system, only a
discriminated signal was registered from each wire, which results in a
position resolution in the order of the segmentation (8 mm and 12 mm for the
50 cm and 80 cm versions, respectively). The segmentation is considerably
better than that of the scintillator bars (5 cm), which seems sufficient for
the purpose of the MID detector.

The measurement of position resolution in a specific chamber uses the
information provided by the eight layers considered in the setup. When a track
is identified (with at least seven valid hits along a straight line), the
chamber under study is removed from the tracking, the track is further fitted
considering the remaining (six or seven) hit points, and then the difference
between the measured point in the chamber under study, and the track
intersection is determined.

The position resolution is expected to degrade at high intensities and at
lower energies. Figure 7 shows the measured resolution as a function of the
beam intensity.

Figure 7: MWPC position resolution as a function of beam intensity.
Measurements were performed for beam momentum of 5 GeV/$c$.

There is a clear broadening of the resolution toward high intensities. The
resolution is about 10-30% worse for intensities around 1000 Hz/cm2, where the
efficiency drop was observed.

Regarding the energy dependence, one can expect the degradation of the
apparent resolution due to the scattering of the particles along the beam
line. This means that the observed broadening is not due to the change of the
intrinsic resolution, but rather indicative of the quality of the beamline.
Figure 8 shows the measured position resolution as a function of beam
momentum. Particularly at the lowest momenta, 0.5 and 1 GeV/$c$, the
broadening amounts to about 40-60%. Out of the total length of the tracking
system (more than 2.5 m), 5 mm corresponds to merely 2 mrad in angular spread.

Figure 8: MWPC position resolution as a function of beam momentum.

The systematic measurements presented above indicate that up to beam
intensities of 100 Hz/cm2, similarly to what was observed for the efficiency,
there is a very small degradation of position resolution, and that the angular
spread of the beam particles was sufficiently low for a meaningful measurement
down to $p=2$ GeV/$c$.

##  4 Summary and outlook

The ALICE 3 experiment entails a large-area muon detector, with high and
uniform efficiency, and excellent time resolution, and running with moderate
occupancy and particle fluence. Scintillator, MWPC, and RPC based options are
all taken into account, and two of these technologies were tested at the CERN
T10 test beam facility.

Regarding plastic scintillators, the option that offers good performance on
light-yield output (around 40 photoelectrons), good time resolution ($<2$ ns),
and that represents a moderate-cost solution is the extruded FNAL-NICADD
scintillator bar equipped with WLS fiber and only one SiPM at one side of the
bar. Other possibilities were tested, in particular ELJEN EJ-208 scintillator
bar also exhibits a good performance but the cost per bar is 4 times that of
the FNAL-NICADD option. For the PROTVINO scintillator bar to be competitive,
it had to be equipped with WLS fiber and two SiPMs. However, having two
sensors per bar in the MID detector would increase the cost of the readout
electronics. The other issue is that the plastic is more than 40 years old and
its technical specifications are missing.

Concerning the tested MWPCs, the detectors were functioning reliably during
the test beam period, and the results offer new detailed information with
control on beam intensity, beam energy, and geometry which is matched with the
MID structure. The tested MWPC type is competitive, having high efficiency and
excellent position resolution even beyond the required particle fluence. Its
moderate intrinsic time resolution shall be checked against the expected rate
and occupancy in ALICE 3 MID simulations.

## Acknowledgments

Support for this work has been received from CONAHCyT under the Grants CF No.
2042, CB A1-S-22917 and A1-S-13525; PAPIIT-UNAM under the project IG100524;
VIEP-BUAP 2023/213; from Ministry of Education, Youth, and Sports of the Czech
Republic, grant number LM2023040; as well as from Hungarian National Research,
Development and Innovation Office under the Grants: OTKA-FK-135349, ELKH-KT-
SA-88/2021, TKP2021-NTKA-10, and 2021-4.1.2-NEMZ-KI-2022-00018/11223-23. This
project has received funding from the European Unions Horizon Europe research
and innovation programme under grant agreement No 101057511. Authors
acknowledge the ALICE Collaboration and T10 CERN team for facilities, software
and organisation of the test beam.

## References

  * [1] K. Aamodt et al., The ALICE experiment at the CERN LHC, JINST 3 (2008) S08002.  doi:10.1088/1748-0221/3/08/S08002.
  * [2] W. Busza, K. Rajagopal, W. van der Schee, Heavy Ion Collisions: The Big Picture, and the Big Questions, Ann. Rev. Nucl. Part. Sci. 68 (2018) 339–376.  arXiv:1802.04801, doi:10.1146/annurev-nucl-101917-020852.
  * [3] R. Bala, I. Bautista, J. Bielcikova, A. Ortiz, Heavy-ion physics at the LHC: Review of Run I results, Int. J. Mod. Phys. E 25 (07) (2016) 1642006.  arXiv:1605.03939, doi:10.1142/S0218301316420064.
  * [4] Letter of intent for ALICE 3: A next-generation heavy-ion experiment at the LHC (11 2022).  arXiv:2211.02491.
  * [5] A. Pla-Dalmau, A. D. Bross, K. L. Mellott, Low-cost extruded plastic scintillator, Nucl. Instrum. Meth. A 466 (2001) 482–491.  doi:10.1016/S0168-9002(01)00177-2.
  * [6] Eljen technology, https://eljentechnology.com/products/plastic-scintillators, accessed: 2023-11-8.
  * [7] J. Bernhard, The New CERN Low-energy Facilities for Neutrino Detector Tests, J. Phys.: Conf. Ser. 1486 (2020) 012127.
  * [8] K. S. McFarland, MINERvA: A Dedicated neutrino scattering experiment at NuMI, Nucl. Phys. B Proc. Suppl. 159 (2006) 107–112.  arXiv:physics/0605088, doi:10.1016/j.nuclphysbps.2006.08.073.
  * [9] D. G. Michael et al., The magnetized steel and scintillator calorimeters of the minos experiment, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 596 (2) (2008) 190–228.  doi:https://doi.org/10.1016/j.nima.2008.08.003.
  * [10] O. A. Grachov, T. M. Cormier, A. Pla-Dalmau, A. Bross, V. Rykalin, Study of new FNAL-NICADD extruded scintillator as active media of large EMCal of ALICE at LHC, in: 11th International Conference on Calorimetry in High-Energy Physics (Calor 2004), 2004, pp. 49–54.  arXiv:physics/0405028, doi:10.1142/9789812701978_0007.
  * [11] Kuraray, http://kuraraypsf.jp/psf/ws.html, accessed: 2023-11-8.
  * [12] D. Varga et al., Close cathode chamber: Low material budget MWPC, Nucl. Instrum. Methods Phys. Res. A 698 (2016) 11–18.  doi:https://doi.org/10.1016/j.nima.2012.09.025.
  * [13] T. V. Acconcia et al., A Very High Momentum Particle Identification Detector, Eur. Phys. J. Plus 129 (2014) 91.  arXiv:1309.5880, doi:10.1140/epjp/i2014-14091-5.
  * [14] D. Varga et al., Detector developments for high performance muography applications, Nucl. Instrum. Methods Phys. Res. A 958 (2020) 162236.  doi:10.1016/j.nima.2019.05.077.
  * [15] D. Varga et al., High Efficiency Gaseous Tracking Detector for Cosmic Muon Madiography, Adv. High Energy Phys. 2016 (2016) 1962317.  doi:https://doi.org/10.1155/2016/1962317.
  * [16] L. Oláh, G. G. Barnaföldi, G. Hamar, H. G. Melegh, G. Surányi, D. Varga, Ccc-based muon telescope for examination of natural caves, Geoscientific Instrumentation, Methods and Data Systems 1 (2) (2012) 229–234.  doi:10.5194/gi-1-229-2012.


# Embedded Topological Semimetals

Saavanth Velury  Taylor L. Hughes  Department of Physics and Institute of
Condensed Matter Theory, University of Illinois at Urbana-Champaign, IL 61801,
USA

###### Abstract

Topological semimetals, such as Dirac, Weyl, or line-node semimetals, are
gapless states of matter characterized by their nodal band structures and
surface states. In this work, we consider layered (topologically trivial)
insulating systems in $D$ dimensions that are composed of coupled multi-layers
of $d$-dimensional topological semimetals. Despite being nominal bulk
insulators, we show that crystal defects having co-dimension $(D-d)$ can
harbor robust lower dimensional topological semimetals embedded in a trivial
insulating background. As an example we show that defect-bound topological
semimetals can be localized on stacking faults and partial dislocations.
Finally, we propose how an embedded topological Dirac semimetal can be
identified in experiment by introducing a magnetic field and resolving the
relativistic massless Dirac Landau level spectrum at low energies in an
otherwise gapped system.

###### pacs:

Valid PACS appear here

††preprint: APS/123-QED

##  I Introduction

Over the past decade, the classification of gapped, non-interacting
topological phases of matter has been extensively developed. Starting with
strong topological insulators protected by only time-reversal, particle-hole,
and/or chiral symmetriesKane and Mele (2005a, b); Bernevig  _et al._ (2006a);
Schnyder  _et al._ (2008); Kitaev (2009); Qi  _et al._ (2009); Ryu  _et al._
(2010); Hasan and Kane (2010); Qi and Zhang (2011), the classification has
grown to include weak topological insulators (those protected by translation
symmetries)Fu  _et al._ (2007); Moore and Balents (2007); Roy (2009) as well
as those protected by crystalline symmetriesFu (2011); Fang  _et al._ (2012);
Chiu  _et al._ (2013); Fang  _et al._ (2013); Jadaun  _et al._ (2013); Slager
_et al._ (2013); Benalcazar  _et al._ (2014); Shiozaki and Sato (2014);
Kruthoff  _et al._ (2017); Song  _et al._ (2017); Po  _et al._ (2017); Bradlyn
_et al._ (2017); Cano  _et al._ (2018); Po  _et al._ (2018); Song  _et al._
(2018); Ono  _et al._ (2020); Po (2020); Song  _et al._ (2020); Elcoro  _et
al._ (2020); Ono  _et al._ (2021).Mirroring this development in topological
insulators, the classification of topological semimetals has seen recent rapid
developmentYang and Nagaosa (2014); Bradlyn  _et al._ (2016); Watanabe  _et
al._ (2016); Wu  _et al._ (2021) building upon earlier work that explored
their exotic transport propertiesHosur  _et al._ (2012); Hosur and Qi (2013);
Wang  _et al._ (2017), robust anomalous surface statesWan  _et al._ (2011);
Matsuura  _et al._ (2013); Xu  _et al._ (2015); Bian  _et al._ (2016a);
Kargarian  _et al._ (2016); Wieder  _et al._ (2020); Obakpolor and Hosur
(2021), and electromagnetic responsesShi  _et al._ (2007); Zyuzin and Burkov
(2012); Vazifeh and Franz (2013); Haldane (2014); Ramamurthy and Hughes (2015,
2017); Armitage  _et al._ (2018); Robredo  _et al._ (2021); Gioia  _et al._
(2021); Pal  _et al._ (2021).

Unlike their insulating counterparts, topological semimetals are gapless in
nature, and their topology is characterized by their symmetry/topology
protected band degeneracies (nodes) near the Fermi level. In general,
topological semimetals can be classified as Dirac, Weyl, or line-node based on
the dispersion of their nodes, and what symmetries are required to protect the
degeneracies. This classification has become more refined through
consideration of other features such as the number of degeneracies, the
resulting topology of the Fermi surface, their distribution throughout the
Brillouin zone (BZ), etc. Additionally, the progress made in the theoretical
examination of topological semimetals has been aided by their experimental
discovery in real materials. Some notable examples include the three-
dimensional Dirac semimetals Na3BiLiu  _et al._ (2014a) and Cd3As2Liu  _et
al._ (2014b); Neupane  _et al._ (2014), the Weyl semimetal TaAsXu  _et al._
(2015), and line-node semimetals PbTaSe2Bian  _et al._ (2016b) and TlTaSe2Bian
_et al._ (2016a). There have also been many material predictions made for
topological semimetals based on first-principles calculationsXu  _et al._
(2011); Burkov and Balents (2011); Wang  _et al._ (2012); Young  _et al._
(2012); Steinberg  _et al._ (2014); Bulmash  _et al._ (2014); Chen  _et al._
(2017) that will likely lead to new experimental discoveries.

The connection between symmetry-protected topology and defect sensitivity is
fundamental to the study of topological insulatorsTeo and Hughes (2017).
Previous work has unveiled a rich set of phenomena associated to the
generation of stable fermionic bound states localized on various types of
topological defects, e.g., vortices, dislocations, or disclinationsRan  _et
al._ (2009); Teo and Kane (2010); Teo and Hughes (2013); Benalcazar  _et al._
(2014); Queiroz  _et al._ (2019); Li  _et al._ (2020); Roy and Juričić (2021);
Nag and Roy (2021). Interestingly, a recent workTuegel  _et al._ (2019)
provided evidence of a nominally topologically trivial insulator capable of
hosting topological features on crystal defects if it is formed from lower-
dimensional layers of topologically non-trivial insulators. While the total
topology of the entire layered system is trivial, the topology of isolated
layers can be revealed with stacking faults or partial dislocations, which
hence give rise to protected surface states at the boundaries of the defect.
This concept was dubbed an embedded topological insulator (ETI) - topological
insulators embedded as defects in (topologically) trivial insulating
environments of greater dimensionality.

It was subsequently shown in later worksKhalaf (2018); Matsugatani and
Watanabe (2018); Trifunovic and Brouwer (2019) that higher-order topological
insulators (HOTIs) can be smoothly deformed to ETIs (and vice versa) via a
procedure that preserves the bulk gap and global crystalline symmetry at the
expense of breaking translational and local crystalline symmetries. Given a
$d$-dimensional HOTI hosting an $n^{\text{th}}$-order gapless boundary mode
($d>n$), an ETI can be generated by trivializing the regions of the HOTI
adjacent to this boundary mode, resulting in a gapped topological phase with
codimension $n-1$ embedded in a trivial bulk, which is an ETI. A specific
example of this is considered inMatsugatani and Watanabe (2018) for a 3D
inversion-symmetric HOTI composed of a coupled layers of Chern insulators with
alternating Chern number $C=\pm 1$, with a single Chern insulator of $C=+1$
containing the inversion center. The Chern layer containing the inversion
center hosts the chiral edge mode that corresponds to the hinge mode of this
HOTI. One can smoothly adjust the inter-layer coupling to dimerize the layers
above and below the inversion center, resulting in gapless edge modes with
opposite chirality becoming gapped out and effectively trivializing the bulk
halves neighboring the inversion center. The residual Chern insulator
containing the inversion center is precisely a 2D embedded Chern insulator.

In light of these results, it is natural to ask if topological semimetals can
emerge in a trivial, bulk insulating environment through the introduction of
defects, i.e., can we have robust, embedded topological semimetals? A
precursor to this question was answered in the context of topological band
insulators (TBIs) in Ref. Slager  _et al._ , 2016, which demonstrated that
grain boundaries in TBIs can host two-dimensional time-reversal symmetry
protected semimetals. From this, it is natural to ask if other types of
defects, such as stacking faults or dislocations, can give rise to novel
embedded topological (semi)metallic phases that would be unexpected in a bulk
insulator. Indeed, since transport is a global measurement, such defect-
generated conductivity would show up as seemingly spurious transport signals
in otherwise insulating system.

In this work, we answer these questions by extending the idea presented in
Ref. Tuegel  _et al._ , 2019. Namely, we introduce the notion of an embedded
topological semimetal (ETSM). The main focus of our work is on layered
(topologically trivial) insulating systems in $D$-dimensions that are composed
of coupled multilayers of $d$-dimensional topological semimetals. We will
consider the effects of two types crystalline defects: stacking faults and
partial dislocations. We demonstrate that when these defects are introduced in
these nominal bulk insulators, defect-bound lower dimensional topological
semimetals can emerge, embedded in a trivial insulating background.
Furthermore, we show that when these defects host 2D topological Dirac
semimetals their presence can be detected by observing the characteristic
Landau level spectrum for massless Dirac fermions.

This article is organized as follows. In Sec. II, we initially build up the
concept of an ETSM by first presenting a model for embedded metals occurring
in a $D=2$-dimensional (2D) insulating system composed of layered
$d=1$-dimensional (1D) spinless electron chains coupled to each other. By
introducing two stacking faults at a set distance apart from each other, the
system transitions from being a bulk insulator having metallic behavior near
the defects. Then, in Sec. III, we introduce our definitive example of an
ETSM: an embedded topological Dirac semimetal. The construction involves a
$D=3$-dimensional (3D) topologically trivial insulating system composed of
layered 2D topological Dirac semimetals. We demonstrate that introducing a
stacking fault or partial dislocation in this system results in a topological
Dirac semimetal in an otherwise trivially gapped system. In Sec. IV, we
discuss how an ETSM can be experimentally determined by introducing a magnetic
field and resolving the relativistic massless Dirac Landau level spectrum at
low energies in an otherwise gapped system. Finally, in Sec. V, we offer
concluding remarks.

##  II Embedded Metals in (1+1)-D

We begin by introducing the model of embedded metals in a 2D insulating
system. This section serves as a warmup to the concept of an ETSM. Although
the embedded metals that we explain in detail in this section are
topologically trivial and therefore, not ETSMs, it serves as an example of the
importance of translational-symmetry breaking required to realize the embedded
nature of these systems that we will discuss throughout this work.
Specifically, we consider stacking faults introduced at certain locations on
the real-space lattice and directly correspond them to the specific bands in
the bulk band structure that intersect the Fermi level. We will show that
these metallic bands arise directly from the exposed lower dimensional
metallic layers where the stacking faults are introduced, which we identify as
the embedded metals.
First, we will describe the pristine model when no stacking faults are
present. This model consists of 1D spinless electron chains extended in the
$x$-direction, and stacked and coupled to each other along the $y$-direction.
Each one-dimensional chain is described by the Bloch Hamiltonian

$\begin{split}h_{\text{1D}}(k_{x})=\pm 2t\cos(k_{x}),\end{split}$ (1)

where $t$ is the strength of the nearest-neighbor hopping, and we have set the
lattice constant $a=1$. Layering these chains on top of each other, we place
two chains with opposite hoppings $+t$ and $-t$ within each unit cell along
the $y$-direction, with a total of $N_{y}$ unit cells in the $y$-direction.
This is illustrated in Fig. 1. We denote the chain with hopping $-t$ with the
label $A$ and the other with hopping $+t$ with the label $B$.

Figure 1: Illustration of the stacking fault procedure. The pristine bulk
insulator described by the Bloch Hamiltonian given by Eq. (2) is shown above
for $N_{y}=8$ unit cells (indicated by the dashed outlines). Each unit cell
contains two spinless 1D electron chains with hopping $-t$ (denoted by the
horizontal red colored bonds, labeled $A$ within each unit cell) and $+t$
(denoted by the horizontal blue colored bonds, labeled $B$ within each unit
cell). The intra-cell coupling strength is denoted by $\gamma$ (vertical green
colored bonds in each unit cell) and the inter-cell coupling strength is
denoted by $\delta$ (vertical orange colored bonds between unit cells). Two
chains are extracted, one of type $A$ and the other of type $B$ as
highlighted. Figure 2: Band structure of the pristine bulk insulator given by
Eq. (2) shown in (a), for $\gamma=0.75$, $\delta=0.5$, and $t=1$ for
$N_{y}=50$ unit cells ($100$ spinless electron chains), with periodic
boundaries in the $y$ direction. In (b), the band structure of an embedded
metal after the stacking faults are introduced at $n_{y}=6^{\text{th}}$ (a $B$
chain is extracted, the effect on the band structure is highlighted by the
green bands) and at $N_{y}+1-n_{y}=45^{\text{th}}$ unit cells (an $A$ chain is
extracted, the effect on the band structure is highlighted by the red bands).
The exposed chains are coupled with strength $\lambda=\delta=0.5$.

The Bloch Hamiltonian describing this model is given by
$H(\mathbf{k})=\vec{\mathbf{d}}(\mathbf{k})\cdot\vec{\bm{\sigma}}$ with

$\begin{split}\vec{\mathbf{d}}=(\gamma+2\delta\cos(k_{y}),2\delta\sin(k_{y}),-2t\cos(k_{x})),\end{split}$
(2)

where $\gamma$ is the intra-cell hopping between the chains, $\delta$ is the
inter-cell hopping, and $\vec{\bm{\sigma}}=(\sigma_{x},\sigma_{y},\sigma_{z})$
are the Pauli matrices expressed in the wire basis within each unit cell
(i.e., $A$ and $B$). The pristine Hamiltonian has a spinless time-reversal
symmetry given by $T=K$ where $K$ is the complex conjugation operator
($TH(\mathbf{k})T^{-1}=H^{*}(\mathbf{k})=H(-\mathbf{k})$). The energy spectrum
of (2) is given by
$E(\mathbf{k})=\pm|\vec{\mathbf{d}}|=\pm\sqrt{(\gamma+\delta\cos(k_{y}))^{2}+(\delta\sin(k_{y}))^{2}+(-2t\cos(k_{x}))^{2}}$.
Thus, for nonzero $\gamma$ and $\delta$, the spectrum of the Hamiltonian is
gapped for all $|\gamma|\neq|\delta|$, and therefore a bulk insulator. The
band structure of the pristine bulk insulator is illustrated in Fig. 2 (a).

Now we will consider the effects of introducing stacking faults, as
illustrated in Fig. 1. We generate stacking faults by removing a spinless
electron chain of type $A$ and another of type $B$ at two separate locations,
which will introduce two inequivalent stacking faults. A chain of type $B$ is
removed from the $n_{y}^{\text{th}}$ unit cell where $1\leq n_{y}<N_{y}/2$,
and a chain of type $A$ is removed from the $(N_{y}+1-n_{y})^{\text{th}}$ unit
cell. The two chains that neighbor each of the removed lines are then coupled
with strength $\lambda$. The resulting energy spectrum after this procedure is
illustrated in Fig. 2 (b). We find two metallic dispersing bands that lie in
the insulating gap, one arising from, and localized at, each stacking fault.
There are gapless Fermi points ($E_{F}=0$) at $k_{x}=\pm\frac{\pi}{2}$.
Intuitively, this occurs because the exposed $A$ and $B$ chains in the
$n_{y}^{\text{th}}$ and $(N_{y}+1-n_{y})^{\text{th}}$ unit cells respectively
each have energy dispersion given by $E_{\text{1D}}^{(\pm)}(k_{x})=\pm
2t\cos(k_{x})$, with the exposed $A$ chain carrying the positive dispersion
and the exposed $B$ chain carrying the negative one. The metallic dispersions
persist for any value of $\lambda\in[0,\max(\gamma,\delta)]$, not just
$\lambda=\delta$. However, the Fermi velocity of the band dispersion around
$k_{x}=\pm\frac{\pi}{2}$, decreases as $\lambda$ is increased. Additional
plots demonstrating this, and a $k\cdot p$ expansion of the dispersion around
$k_{x}=\pm\frac{\pi}{2},$ are provided in Appendix A. Therefore, the exposed
$A$ and $B$ chains lying on the stacking faults are denoted as the embedded
metals of the system. Thus, this example is a clear illustration that a
metallic system can be realized from a bulk insulator through stacking fault
defects that expose the underlying embedded metals within the system.

##  III Embedded Topological Dirac Semimetal in (2+1)-D

Having established how a bulk insulator can be turned into a metallic system
by introducing stacking faults, we now turn our attention $D=3$ dimensions.
Our paradigmatic example is the embedded topological Dirac semimetal (ETDSM),
which we will establish in this section. We will consider a 3D nominal bulk
insulator constructed from 2D topological Dirac semimetal layers. First, we
will establish the technical details of the model, including the layer
Hamiltonian and its phase diagram, and then the full model Hamiltonian built
from coupling the individual layers. Then, we establish how the stacking fault
is introduced into our system, and illustrate the change in the bulk band
structure, which transforms from being gapped to becoming gapless with robust
Dirac nodes. We identify these Dirac nodes as those arising from the exposed
2D topological Dirac semimetal layer lying on the stacking fault, which is the
ETDSM. Finally, we will characterize the topology of the resulting ETSM by
analyzing the Berry phase characteristics across the BZ. In this discussion,
we note that the symmetries required for topological protection of the Dirac
nodes of the ETDSM are the same symmetries required to protect the Dirac nodes
of the individual topological Dirac semimetal layers.

Figure 3: Phase diagram of the 2D layer Hamiltonian given by Eq. (4) for model
parameters $m$ and $t_{y}$. The phases of this model include the ordinary
insulator (OI) phase, weak topological insulator (WTI) phase, Dirac semimetal
1 (DSM 1) phase, and Dirac semimetal 2 phase (DSM 2) phase. The DSM 1 phase
contains one pair of Dirac nodes as opposed to two pairs of Dirac nodes in the
DSM 2 phase.

###  III.1 Layer Hamiltonian and Phase Diagram

Each of the 2D layers comprising our full 3D model are adapted from a square
lattice Dirac semimetal model from Ref. Ramamurthy and Hughes, 2015, the
details of which we briefly review in this subsection. The 2D Dirac semimetal
layer itself is constructed from coupling 1D topological insulator chains
having Bloch Hamiltonian:

$\begin{split}H_{\text{1DTI}}(k_{x})=\sin(k_{x})\sigma_{y}+(1-m-\cos(k_{x}))\sigma_{z},\end{split}$
(3)

where $m$ is a parameter controlling the insulating phases of the model, and
the Pauli matrices act on spinless orbitals labeled by $1$ and $2$
respectively. This model has an inversion symmetry given by $I=\sigma_{z}$
($IH_{\text{1DTI}}(k_{x})I^{-1}=H_{\text{1DTI}}(-k_{x})$), which can protect
an inversion symmetric topological insulator phaseHughes  _et al._ (2011). The
important topological invariant that characterizes this model is
$\mathbb{Z}_{2}$ in nature, and is captured by the bulk charge polarization
$P$. For $m<0$ or $m>2$, the bulk charge polarization is given by $P=0\hskip
2.84544pt(\text{mod}\hskip 1.42271pte),$ and for $0<m<2$, $P=\frac{e}{2}\hskip
2.84544pt(\text{mod}\hskip 1.42271pte)$. Hence, when $0<m<2$, Eq. (3) is said
to be in its obstructed atomic limit phaseBradlyn  _et al._ (2017), whereas
when $m<0$ or $m>2$, it is said to be in its trivial phase.

To construct a single layer of our full three-dimensional insulator, we stack
several copies of (3) along the $y$-direction, and introduce a coupling
strength $t_{y}$ along the stacking direction. This leads to the following
Bloch Hamiltonian for each layer:

$\begin{split}H_{\text{2DWTI}}(k_{x},k_{y})=H_{\text{1DTI}}(k_{x})-t_{y}\cos(k_{y})\sigma_{z}.\end{split}$
(4)

For $0<m<2$, as long no solution exists for (i) $\cos(k_{y})=-\frac{m}{t_{y}}$
and (ii) $\cos(k_{y})=\frac{2-m}{t_{y}}$, the system is a weak topological
insulator. If for any value of $(m,k_{x},k_{y})$, a solution exists to either
(i) or (ii), the system becomes a Dirac semimetal with a single pair of Dirac
nodes in the BZ (DSM 1 in Fig. 3). Similarly, if solutions exist for both (i)
and (ii), then the system is a Dirac semimetal with two pairs of Dirac nodes
in the BZ (DSM 2 in Fig. 3). These results can be summarized in a phase
diagram for (4) shown in Fig. 3. In both Dirac semimetal phases, the nodes are
stabilized by a composite $TI$ symmetry (where $T=K$ is spinless time-reversal
symmetry). For simplicity, we will consider the regime where $|m|<|t_{y}|$ or
$|2-m|<|t_{y}|$ so that each layer is in the DSM 1 phase.

###  III.2 Model Hamiltonian and Stacking Fault Procedure

To construct our 3D bulk insulator, we first layer copies of (4) along the
$z$-direction such that each unit cell in the $z$-direction consists of two
layers of (4). We denote these two layers as $A$ and $B$ which have model
parameters $(m_{A},t_{y,A})$ and $(m_{B},t_{y,B})$ respectively. In this limit
when the layers are decoupled the Bloch Hamiltonian for the pristine 3D system
is

$\begin{split}&H(\mathbf{k})=\sin(k_{x})(\mathbf{1}\otimes\sigma_{y})\\\
&+(1-m-\cos(k_{x})-t_{y}\cos(k_{y}))(\mathbf{1}\otimes\sigma_{z}),\end{split}$
(5)

where the Pauli matrices $\tau_{i}$ ($i=x,y,z$) represent the layers within
each unit cell. The Bloch Hamiltonian in (LABEL:eq:3Dpristine1) maintains the
spinless time-reversal symmetry present in (4), as well as the inversion
symmetry within the $xy$ plane, i.e., $C_{2z}$ symmetry represented by the
operator $C_{2z}=\mathbf{1}\otimes\sigma_{z}$. When the model parameters
$(m,t_{y})$ are tuned to the DSM 1 phase as stated in Subsec. III.1, each
layer is a Dirac semimetal with a pair of Dirac nodes in the bulk BZ. In order
to gap out these Dirac nodes and turn (LABEL:eq:3Dpristine1) into an
insulator, we introduce a coupling term $H_{\text{coupling}}(\mathbf{k})$ that
couples the layers of (LABEL:eq:3Dpristine1):

$\begin{split}&H_{\text{coupling}}(\mathbf{k})=(\gamma_{z}+\delta_{z}\cos(k_{z}))(\tau_{y}\otimes\sigma_{x})\\\
&+(\delta_{z}\sin(k_{z}))(\tau_{x}\otimes\sigma_{x}),\end{split}$ (6)

where $\gamma_{z}$ and $\delta_{z}$ are the intra-cell and inter-cell
couplings respectively. We note that in the presence of the coupling term the
model has broken $T$ and $C_{2z},$ but the composite $C_{2z}T$ symmetry is
maintained. Hence, when each individual layer of our 3D model is tuned to the
DSM 1 phase in (LABEL:eq:3Dpristine1), the coupling Eq. (LABEL:eq:3Dpristine2)
gaps the entire 3D system to form a trivial insulator. The open boundary
spectrum of the 3D model with the coupling term included is illustrated in
Fig. 4 (a). The bulk, along with the typical Dirac semimetal edge states, are
completely gapped out, indicating that the system is a trivial insulator in
the absence of a stacking fault.

Figure 4: Band structure of the pristine 3D bulk insulator given by Eqns.
(LABEL:eq:3Dpristine1) and (LABEL:eq:3Dpristine2) shown in (a), for $m=2$,
$t_{y}=1$ with $\gamma_{z}=0.05$ and $\delta_{z}=0.5$ with open boundaries in
the $x$ direction and periodic boundaries in the $z$ direction, for $N_{z}=14$
unit cells ($28$ layers in total) and $N_{x}=28$ sites. The bulk band
structure after the stacking fault procedure is performed is shown in (b)
(periodic boundaries in both the $x$ and $z$ directions). The layer in the
$7^{\text{th}}$ unit cell ($14^{\text{th}}$ layer) was extracted and the
exposed layers are rejoined with coupling strength $\delta_{z}=0.5$. Dirac
nodes appear in the bulk at $k_{y}=\pm\frac{\pi}{2}$. Figure 5: Illustration
of stacking fault procedure (Volterra construction) in the layered
topologically trivial 3D insulator. In (a), the pristine 3D insulator is
shown. The $B$ layer to be extracted is indicated by the yellow outline. Once
this stacking fault is introduced, the resulting structure is shown in (b).

We now introduce a stacking fault into this system by extracting a single
layer from the $\left(\frac{N_{z}}{2}\right)^{\text{th}}$ unit cell (i.e., the
$N_{z}^{\text{th}}$ layer) (Fig. 5). The two exposed layers are then coupled
with a coupling strength $\delta_{z}$, and the resulting band structure is
shown in Fig. 4(b), which features Dirac nodes at $k_{y}=\pm\frac{\pi}{2}$.
Therefore, by introducing a stacking fault, we have explicitly broken the unit
cell structure of the pristine trivial insulator and generated an embedded
Dirac semimetal. As we will see below it is also possible for other types of
metallic features to emerge when partial dislocation defects are introduced
into this trivial insulator (see Subsec. III.4), but first we discuss the
stability of the Dirac nodes in the ETSM.

###  III.3 Characterizing the Topology of the ETSM

In this subsection, we will demonstrate the topological nature of the embedded
Dirac semimetal by calculating the Zak-Berry phase for a system with and
without a stacking fault. Since we will consider the presence of a defect, we
treat the $z$-direction in real-space and have an effective 2D BZ spanned by
$(k_{x},k_{y}).$ For each fixed value of $k_{y}$ in this 2D BZ we can consider
a quasi-1D system extended along the $k_{x}$ direction. Hence for each fixed
$k_{y}$ we can consider the Zak-Berry phase along $k_{x}$ by integrating the
adiabatic connection along $k_{x}$Zak (1989):

$\begin{split}\theta_{B}(k_{y})=\int\limits_{-\pi}^{\pi}dk_{x}\hskip
2.84544pt\text{tr}\hskip 1.42271pt\mathcal{A}_{x}(k_{x}),\end{split}$ (7)

where

$\begin{split}\mathcal{A}_{x}(k_{x})=\langle
u_{m}(k_{x},k_{y})|i\nabla_{k_{x}}|u_{n}(k_{x},k_{y})\rangle,\end{split}$ (8)

is the non-Abelian adiabatic connection derived from the occupied states
(i.e., states with energy less than the Fermi-level $E_{F}=0$). For each value
of $k_{y}$ this quasi-1D Berry phase is quantized by the $C_{2z}T$ symmetry
and can take values of only $0$ or $\pi$ if $k_{y}$ does not pass through a
gapless Dirac point (which we avoid in our plots). We compute $\theta_{B}$
along a discretized grid of $k_{x}$ points at each $k_{y}$ value (excluding
the values where the Dirac nodes occur at $k_{y}=\pm\frac{\pi}{2}$). To do
this, we use the discretized gauge invariant Wilson loop expression

$\begin{split}&\theta_{B}(k_{y})=\text{Im}\hskip
1.42271pt\text{log}\prod\limits_{i=1}^{N_{\text{grid},x}}\text{det}\hskip
2.84544pt\mathcal{U}_{i}(k_{y})\\\ &[\mathcal{U}_{i}(k_{y})]_{mn}=\langle
u_{m}(k_{x,i},k_{y})|u_{n}(k_{x,i+1},k_{y})\rangle.\end{split}$ (9)

Figure 6: Plot of the discretized (Zak)-Berry phase given by Eq.
(LABEL:eq:discretizedZakBerryphase) as a function of $k_{y}$, for a
discretized 2D BZ of $N_{\text{grid},x}=100$ points along $k_{x}$ by
$N_{\text{grid},y}=500$ points along $k_{y}$. Identical model parameters were
used as in Fig. 4. The plot in (a) is when the system is a trivial insulator
with no stacking fault and (b) corresponds to the system when the stacking
fault is introduced. The dashed vertical lines indicate discontinuities where
the Berry phase cannot be computed, which is precisely at the location of the
Dirac nodes, $k_{y}=\pm\frac{\pi}{2}$.

A plot of the calculated Berry phase is shown in Fig. 6 as $k_{y}$ is tuned
from $-\pi$ to $\pi.$ For Fig. 6 (a) we consider a system without a stacking
fault and the Berry phase is constant across the $k_{y}$ BZ. For Fig. 6 (b) we
show the results for a single stacking fault and find that the Berry phase
remains quantized, but it changes its value by a quantized amount precisely at
the locations of the Dirac nodes, $k_{y}=\pm\frac{\pi}{2}$. This result can be
understood as follows. In computing the Berry phase across the BZ, we
effectively treated the BZ as a family of 1D insulators parametrized by
$k_{y}$. That is, each value of $k_{y}$ (except at $k_{y}=\pm\frac{\pi}{2}$)
represents a 1D insulator oriented along the $x$-direction having an effective
$y\to-y$ inversion symmetry generated by $C_{2z}T$. So as $k_{y}$ is tuned
within $\left(-\frac{\pi}{2},\frac{\pi}{2}\right)$, the system is effectively
in a trivial insulator phase since each 1D insulator has Berry phase
$\theta_{B}=0$. Likewise, the system is in a topological insulator (obstructed
atomic limit) phase when $k_{y}$ is varied within
$\left(-\pi,-\frac{\pi}{2}\right)$ and $\left(\frac{\pi}{2},\pi\right]$ since
each 1D insulator has Berry phase $\theta_{B}=\pi$. Therefore, each Dirac node
appearing in momentum space represents a transition between a trivial phase
and a topological phase based on the Berry phase characteristics, making this
an ETDSM. The key quantized Berry phase jump that illustrates the topological
nature of the 2D Dirac cones is protected by the $C_{2z}T$ symmetry, and this
symmetry also provides the local momentum-space stability of the defect-
localized Dirac node degeneracy. We also characterize the topological nature
of the embedded Dirac semimetal using the single-particle entanglement
spectrum and by computing the entanglement Berry phase in analogy to Ref.
Tuegel  _et al._ , 2019, which is discussed in Appendix B.

###  III.4 Topological Metallic States Induced by a Partial Dislocation

In addition to stacking faults, partial dislocations (which can be
heuristically thought of as terminated stacking faults) provide another
example of a crystalline defect that can give rise to metallic states when
introduced to a trivial insulating system. In particular, we show that
metallic states that are topologically protected can emerge in a trivial
background when a partial dislocation is applied to our 3D model. To implement
the partial dislocation, we remove only half of the layer located in the
$\left(\frac{N_{z}}{2}\right)^{\text{th}}$ unit cell (i.e., the
$N_{z}^{\text{th}}$ layer) and the layers above and below are recoupled.
Specifically, we perform the cut parallel to the $y$-direction so that
translational invariance is preserved in this direction and $k_{y}$ remains a
good quantum number. This procedure is illustrated in Fig. 7 (a) and (b), and
the resulting energy spectrum is shown in Fig. 8 (b) which can be compared
with the defectless spectrum shown in Fig. 8 (a). In the system with a partial
dislocation we find a set of gapless modes that connect between the embedded
Dirac nodes at $k_{y}=\pm\frac{\pi}{2}$ and are located in position space on
edge of the cleaved layer, i.e., on the partial dislocation line. These modes
emerge in the system since the layer itself is a Dirac semimetal and these are
the flatband edge states that connect the Dirac nodes. Hence, we have found
that both stacking faults and partial dislocations give rise to embedded
metallic states that emerge in the bulk trivial insulator.

Figure 7: Illustration of the partial dislocation procedure. (a) and (b) show
how a partial dislocation is introduced in the pristine 3D topologically
trivial insulator. A layer is cut in half parallel to the $y$-direction so
that $k_{y}$ remains a good quantum number, and the layers above and below are
recoupled. The band structure of the trivial 3D insulator before the partial
dislocation is illustrated in (a) for periodic boundaries along the $z$
direction and open boundaries along the $x$ direction. As expected, the bulk
is gapped and the edge states are gapped out. The bulk band structure after
the partial dislocation is performed is shown in (b). Figure 8: Illustration
of how topologically protected gapless modes can emerge from an ordinary
insulator. As shown in (a) before the partial dislocation is created, the bulk
is gapped and the edge states are gapped out. The bulk band structure after
the partial dislocation is performed is shown in (b) and was plotted for
periodic boundaries in both the $z$ and $x$ directions. Gapless modes
intervene between $k_{y}=-\frac{\pi}{2}$ and $k_{y}=\frac{\pi}{2}$ that
traverse the boundary of the BZ. Identical model parameters were used to
construct the plots as in Fig. 4.

##  IV Landau Level Spectrum of Embedded Topological Dirac Semimetal

We now shift our focus to discussing a viable method for identifying an ETDSM
in an experimental setting. In particular, we propose that an ETDSM can be
readily identified by resolving its relativistic massless Dirac Landau level
(LL) spectrum and the associated integer QHE sequence. In order to resolve the
LL spectrum, the system with the stacking fault must be placed in a magnetic
field with a non-vanishing component normal to the plane of the stacking
fault. For simplicity we will consider a magnetic field along the
$z$-direction, and use a Landau gauge $\mathbf{A}=Bx\hat{y}.$ We take the
magnetic field strength to be $B=\phi\frac{\Phi_{0}}{a^{2}},$ where
$\Phi_{0}=h/e$ is the flux quantum, and the lattice constant is denoted by $a$
along the $x$ and $y$ directions (we set both $\Phi_{0}$ and $a$ to unity). We
let $\phi=\frac{p}{q}$ be a rational numberHofstadter (1976) where $p$ gives
the number of integer flux quanta per magnetic unit cell of size $q$. The
technical details of the stacking fault Hamiltonian in a magnetic field are
included in Appendix C.
To correctly obtain the Dirac LL spectrum from a lattice model we must choose
a large enough $q$ to recover the continuum resultBernevig  _et al._ (2006b,
2007). Diagonalizing this Hamiltonian for $\phi=1/31$ produces the Landau
levels shown in Fig. 9. Fig. 9 (a) shows the LL spectrum for the insulator
given by Eqns. (LABEL:eq:3Dpristine1) and (LABEL:eq:3Dpristine2) with no
stacking fault, while Fig. 9 (b) shows the LL spectrum for the ETDSM when the
stacking fault is introduced. We see a clear zeroth LL (formed by four
degenerate bands at $E=0$) and non-uniform LL spacing in Fig. 9 (b) which is
expected in a system with 2D Dirac points. These zero-energy bands are a key
feature that distinguishes the trivial insulator from the ETDSM and signifies
the presence of bulk Dirac nodes. To summarize, we find that when the ETDSM is
subject to a magnetic field, the Dirac LL spectrum can be successfully
resolved at low energies. Consequently, given the consistency of our numerical
results presented in this section with those found in Refs. Bernevig  _et al._
, 2006b, 2007, for a graphene lattice model, we expect the ETDSM to exhibit an
associated relativistic integer QHE sequence, which could potentially be an
experimentally accessible signature. We note that a similar proposal of
measuring the unconventional QHE sequence was made in Ref. Slager  _et al._ ,
2016 to identify time-reversal symmetry protected semimetals along grain
boundaries in TBIs.

##  V Conclusion

In conclusion, we have investigated the viability of realizing new,
potentially topological, (semi)-metallic phases in systems via defect
insertions in ordinary insulators. We established the notion of ETSMs as
topological semimetals localized on defects embedded in trivial insulating
environments, by considering topologically trivial insulating systems in $D$
dimensions composed of coupled multi-layers $d$-dimensional topological
semimetals. Starting with the simplest example in $D=2$ dimensions, we
considered an ordinary insulator composed of layered $d=1$ spinless electron
chains, and introduced two stacking faults, which unveiled an embedded metal
with metallic Fermi points in its band structure. Extending this to $D=3$
dimensions, we considered an insulator consisting of coupled layers of 2D
Dirac semimetals. Introducing a stacking fault in this system resulted in an
ETDSM, with visible Dirac nodes in the bulk band structure. We also showed the
partial dislocations can bind the expected edge states of a 2D Dirac
semimetal. Finally, we discussed how an ETSM can be identified in experiment
by resolving the relativistic massless Dirac LL spectrum at low energies.
Although embedded Weyl and line-node semimetals were not explicitly considered
in this work, we note that the construction laid out here for embedded Dirac
semimetals applies to these systems straightforwardly for $D=4$ dimensions.
For example, an embedded Weyl semimetal can be realized in an topologically
trivial insulator in $D=4$ dimensions composed of 3D layers, which themselves
are coupled multi-layers of 2D ($d=2$) Chern insulators, by applying a
stacking fault. However, because these systems require higher dimensionality
beyond $D=3$ dimensions, their physical relevance will be limited to
engineered materials having additional synthetic dimensions.

Figure 9: Landau Level spectrum of (a) the trivial insulator before the
stacking fault is introduced and (b) the ETDSM when the stacking fault is
introduced. For each plot, $N_{x}=62$ sites and $q=31$ as required by magnetic
translational invariance, with $m=2$, $t_{y}=1$, $\gamma_{z}=0.5$, and
$\delta_{z}=0.15$ with $N_{z}=28$ layers. The Dirac LL spectrum is realized at
this value of $q$ (e.g., low values of $B$) as shown in (b). The LL spacing
changes from the spacing found in relativistic LL to the uniform spacing found
in non-relativistic LL. Furthermore, a four-fold degenerate set of zero energy
bands appears in the ETDSM spectrum of (b), and (c) is a zoomed-in plot of the
relativistic LL in (b) illustrating the zero energy bands. (d) is a plot of
the zero mode wavefunction probability distribution along the stacking
direction. The zero modes are localized along the layer with the ETDSM after
the stacking fault is introduced, as shown in (d).

Finally, it warrants further examination whether embedded topological phases
can be related to other types of higher-order topological phases. As mentioned
in the introduction, it has been shown previously inKhalaf (2018); Matsugatani
and Watanabe (2018); Trifunovic and Brouwer (2019) that HOTIs can be smoothly
deformed to ETIs. However, it may be possible that higher-order topological
semimetals (HOTSMs)Lin and Hughes (2018); Wang  _et al._ (2020); Ghorashi  _et
al._ (2020); Simon  _et al._ (2021) can also harbor ETIs. For example, the 3D
higher-order Dirac semimetal discussed inLin and Hughes (2018) and the higher-
order Weyl semimetal discussed in Ghorashi  _et al._ (2020) relies on coupled-
layer constructions of 2D quadrupole models. If one considers a dimerized
coupling structure along the stacking direction of the quadrupole models that
breaks the translational symmetry but preserves crystalline symmetries such as
$C_{4}$ and/or mirror symmetries, this can result in an embedded quadrupole
insulator. For the ETSMs considered in this work, it may be viable to employ a
similar type of argument to relate ETSMs to HOTIs or HOTSMs, which we leave to
future study.

##  VI Acknowledgments

We thank Penghao Zhu and Ze-Min Huang for insightful discussions, and Nick
Abboud for some useful initial discussions on this work. S.V. is supported by
the NSF Graduate Research Fellowship Program under Grant No. DGE - 1746047.
T.L.H. thanks ARO under the MURI Grant W911NF2020166.

## References

  * Kane and Mele (2005a) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 (2005a).
  * Kane and Mele (2005b) C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 (2005b).
  * Bernevig  _et al._ (2006a) B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 314, 1757 (2006a).
  * Schnyder  _et al._ (2008) A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Phys. Rev. B 78, 195125 (2008).
  * Kitaev (2009) A. Kitaev, AIP Conference Proceedings 1134 (2009), 10.1063/1.3149495.
  * Qi  _et al._ (2009) X.-L. Qi, T. L. Hughes, S. Raghu, and S.-C. Zhang, Phys. Rev. Lett. 102, 187001 (2009).
  * Ryu  _et al._ (2010) S. Ryu, A. P. Schnyder, A. Furusaki, and A. W. W. Ludwig, New J. Phys. 12, 065010 (2010).
  * Hasan and Kane (2010) M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang (2011) X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Fu  _et al._ (2007) L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 98, 106803 (2007).
  * Moore and Balents (2007) J. E. Moore and L. Balents, Phys. Rev. B 75, 121306(R) (2007).
  * Roy (2009) R. Roy, Phys. Rev. B 79, 195322 (2009).
  * Fu (2011) L. Fu, Phys. Rev. Lett. 106, 106802 (2011).
  * Fang  _et al._ (2012) C. Fang, M. J. Gilbert, and B. A. Bernevig, Phys. Rev. B 86, 115112 (2012).
  * Chiu  _et al._ (2013) C.-K. Chiu, H. Yao, and S. Ryu, Phys. Rev. B 88, 075142 (2013).
  * Fang  _et al._ (2013) C. Fang, M. J. Gilbert, and B. A. Bernevig, Phys. Rev. B 87, 035119 (2013).
  * Jadaun  _et al._ (2013) P. Jadaun, D. Xiao, Q. Niu, and S. K. Banerjee, Phys. Rev. B 88, 085110 (2013).
  * Slager  _et al._ (2013) R.-J. Slager, A. Mesaros, V. Juričić, and J. Zaanen, Nat. Phys. 9, 98 (2013).
  * Benalcazar  _et al._ (2014) W. A. Benalcazar, J. C. Y. Teo, and T. L. Hughes, Phys. Rev. B 89, 224503 (2014).
  * Shiozaki and Sato (2014) K. Shiozaki and M. Sato, Phys. Rev. B 90, 165114 (2014).
  * Kruthoff  _et al._ (2017) J. Kruthoff, J. de Boer, J. van Wezel, C. L. Kane, and R.-J. Slager, Phys. Rev. X 7, 041069 (2017).
  * Song  _et al._ (2017) H. Song, S.-J. Huang, L. Fu, and M. Hermele, Phys. Rev. X 7, 011020 (2017).
  * Po  _et al._ (2017) H. C. Po, A. Vishwanath, and H. Watanabe, Nature Communications 8, 50 (2017).
  * Bradlyn  _et al._ (2017) B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Nature 547, 298 (2017).
  * Cano  _et al._ (2018) J. Cano, B. Bradlyn, Z. Wang, L. Elcoro, M. G. Vergniory, C. Felser, M. I. Aroyo, and B. A. Bernevig, Phys. Rev. B 97, 035139 (2018).
  * Po  _et al._ (2018) H. C. Po, H. Watanabe, and A. Vishwanath, Phys. Rev. Lett. 121, 126402 (2018).
  * Song  _et al._ (2018) Z. Song, T. Zhang, Z. Fang, and C. Fang, Nature Communications 9, 3530 (2018).
  * Ono  _et al._ (2020) S. Ono, H. C. Po, and H. Watanabe, Science Advances 6, eaaz8367 (2020).
  * Po (2020) H. C. Po, Journal of Physics: Condensed Matter 32, 263001 (2020).
  * Song  _et al._ (2020) Z.-D. Song, L. Elcoro, and B. A. Bernevig, Science 367, 794 (2020).
  * Elcoro  _et al._ (2020) L. Elcoro, B. J. Weider, Z. Song, Y. Xu, B. Bradlyn, and B. A. Bernevig, arXiv preprint arXiv:2010.00598 (2020).
  * Ono  _et al._ (2021) S. Ono, H. C. Po, and K. Shiozaki, Phys. Rev. Research 3, 023086 (2021).
  * Yang and Nagaosa (2014) B.-J. Yang and N. Nagaosa, Nature Communications 5, 4898 (2014).
  * Bradlyn  _et al._ (2016) B. Bradlyn, J. Cano, Z. Wang, M. G. Vergniory, C. Felser, R. J. Cava, and B. A. Bernevig, Science 353, aaf5037 (2016).
  * Watanabe  _et al._ (2016) H. Watanabe, H. C. Po, M. P. Zaletel, and A. Vishwanath, Phys. Rev. Lett. 117, 096404 (2016).
  * Wu  _et al._ (2021) L. Wu, F. Tang, and X. Wan, Phys. Rev. B 104, 045107 (2021).
  * Hosur  _et al._ (2012) P. Hosur, S. A. Parameswaran, and A. Vishwanath, Phys. Rev. Lett. 108, 046602 (2012).
  * Hosur and Qi (2013) P. Hosur and X.-L. Qi, Comptes Rendus Physique 14, 857 (2013).
  * Wang  _et al._ (2017) S. Wang, B.-C. Lin, A.-Q. Wang, D.-P. Yu, and Z.-M. Liao, Advances in Physics: X 2, 518 (2017).
  * Wan  _et al._ (2011) X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Phys. Rev. B 83, 205101 (2011).
  * Matsuura  _et al._ (2013) S. Matsuura, P.-Y. Chang, A. P. Schnyder, and S. Ryu, New J. Phys. 15, 065001 (2013).
  * Xu  _et al._ (2015) S.-Y. Xu, C. Liu, S. K. Kushwaha, R. Sankar, J. W. Krizan, I. Belopolski, M. Neupane, G. Bian, N. Alidoust, T.-R. Chang, H.-T. Jeng, C.-Y. Huang, W.-F. Tsai, H. Lin, P. P. Shibayev, F.-C. Chou, R. J. Cava, and M. Z. Hasan, Science 347, 294 (2015).
  * Bian  _et al._ (2016a) G. Bian, T.-R. Chang, H. Zheng, S. Velury, S.-Y. Xu, T. Neupert, C.-K. Chiu, S.-M. Huang, D. S. Sanchez, I. Belopolski, N. Alidoust, P.-J. Chen, G. Chang, A. Bansil, H.-T. Jeng, H. Lin, and M. Z. Hasan, Phys. Rev. B 93, 121113(R) (2016a).
  * Kargarian  _et al._ (2016) M. Kargarian, M. Randeria, and Y.-M. Lu, PNAS 113, 8648 (2016).
  * Wieder  _et al._ (2020) B. J. Wieder, Z. Wang, J. Cano, X. Dai, L. M. Schoop, B. Bradlyn, and B. A. Bernevig, Nature Communications 11, 627 (2020).
  * Obakpolor and Hosur (2021) O. E. Obakpolor and P. Hosur, arXiv preprint arXiv:2108.05380 (2021).
  * Shi  _et al._ (2007) J. Shi, G. Vignale, D. Xiao, and Q. Niu, Phys. Rev. Lett. 99, 197202 (2007).
  * Zyuzin and Burkov (2012) A. A. Zyuzin and A. A. Burkov, Phys. Rev. B 86, 115133 (2012).
  * Vazifeh and Franz (2013) M. M. Vazifeh and M. Franz, Phys. Rev. Lett. 111, 027201 (2013).
  * Haldane (2014) F. D. M. Haldane, arXiv preprint arXiv:1401.0529 (2014).
  * Ramamurthy and Hughes (2015) S. T. Ramamurthy and T. L. Hughes, Phys. Rev. B 92, 085105 (2015).
  * Ramamurthy and Hughes (2017) S. T. Ramamurthy and T. L. Hughes, Phys. Rev. B 95, 075138 (2017).
  * Armitage  _et al._ (2018) N. P. Armitage, E. J. Mele, and A. Vishwanath, Rev. Mod. Phys. 90, 015001 (2018).
  * Robredo  _et al._ (2021) I. Robredo, P. Rao, F. de Juan, A. Bergara, J. L. Manes, A. Cortijo, M. G. Vergniory, and B. Bradlyn, Phys. Rev. Research 3, L032068 (2021).
  * Gioia  _et al._ (2021) L. Gioia, C. Wang, and A. A. Burkov, Phys. Rev. Research 3, 043067 (2021).
  * Pal  _et al._ (2021) H. K. Pal, O. E. Obakpolor, and P. Hosur, arXiv preprint arXiv:2112.14170 (2021).
  * Liu  _et al._ (2014a) Z. K. Liu, B. Zhou, Y. Zhang, Z. J. Wang, H. M. Weng, D. Prabhakaran, S.-K. Mo, Z. X. Shen, Z. Fang, X. Dai, Z. Hussain, and Y. L. Chen, Science 343, 864 (2014a).
  * Liu  _et al._ (2014b) Z. K. Liu, J. Jiang, B. Zhou, Z. J. Wang, Y. Zhang, H. M. Weng, D. Prabhakaran, S.-K. Mo, H. Peng, P. Dudin, T. Kim, M. Hoesch, Z. Fang, X. Dai, Z. X. Shen, D. L. Feng, Z. Hussain, and Y. L. Chen, Nature Materials 13, 677 (2014b).
  * Neupane  _et al._ (2014) M. Neupane, S.-Y. Xu, R. Sankar, N. Alidoust, G. Bian, C. Liu, I. Belopolski, T.-R. Chang, H.-T. Jeng, H. Lin, A. Bansil, F. Chou, and M. Z. Hasan, Nature Communications 5, 3786 (2014).
  * Bian  _et al._ (2016b) G. Bian, T.-R. Chang, R. Sankar, S.-Y. Xu, H. Zheng, T. Neupert, C.-K. Chiu, S.-M. Huang, G. Chang, I. Belopolski, D. S. Sanchez, M. Neupane, N. Alidoust, C. Liu, B. Wang, C.-C. Lee, H.-T. Jeng, C. Zhang, Z. Yuan, S. Jia, A. Bansil, F. Chou, H. Lin, and M. Z. Hasan, Nature Communications 7, 10556 (2016b).
  * Xu  _et al._ (2011) G. Xu, H. Weng, Z. Wang, X. Dai, and Z. Fang, Phys. Rev. Lett. 107, 186806 (2011).
  * Burkov and Balents (2011) A. A. Burkov and L. Balents, Phys. Rev. Lett. 107, 127205 (2011).
  * Wang  _et al._ (2012) Z. Wang, Y. Sun, X.-Q. Chen, C. Franchini, G. Xu, H. Weng, X. Dai, and Z. Fang, Phys. Rev. B 85, 195320 (2012).
  * Young  _et al._ (2012) S. M. Young, S. Zaheer, J. C. Y. Teo, C. L. Kane, E. J. Mele, and A. M. Rappe, Phys. Rev. Lett. 108, 140405 (2012).
  * Steinberg  _et al._ (2014) J. A. Steinberg, S. M. Young, S. Zaheer, C. L. Kane, E. J. Mele, and A. M. Rappe, Phys. Rev. Lett. 112, 036403 (2014).
  * Bulmash  _et al._ (2014) D. Bulmash, C.-X. Liu, and X.-L. Qi, Phys. Rev. B 89, 081106(R) (2014).
  * Chen  _et al._ (2017) R. Chen, H. C. Po, J. B. Neaton, and A. Vishwanath, Nat. Phys. 14, 55 (2017).
  * Teo and Hughes (2017) J. C. Y. Teo and T. L. Hughes, Annual Review of Condensed Matter Physics 8, 211 (2017).
  * Ran  _et al._ (2009) Y. Ran, Y. Zhang, and A. Vishwanath, Nat. Phys. 5, 298 (2009).
  * Teo and Kane (2010) J. C. Y. Teo and C. L. Kane, Phys. Rev. B 82, 115120 (2010).
  * Teo and Hughes (2013) J. C. Y. Teo and T. L. Hughes, Phys. Rev. Lett. 111, 047006 (2013).
  * Queiroz  _et al._ (2019) R. Queiroz, I. C. Fulga, N. Avraham, H. Beidenkopf, and J. Cano, Phys. Rev. Lett. 123, 266802 (2019).
  * Li  _et al._ (2020) T. Li, P. Zhu, W. A. Benalcazar, and T. L. Hughes, Phys. Rev. B 101, 115115 (2020).
  * Roy and Juričić (2021) B. Roy and V. Juričić, Phys. Rev. Research 3, 033107 (2021).
  * Nag and Roy (2021) T. Nag and B. Roy, Nature Communications Physics 4, 157 (2021).
  * Tuegel  _et al._ (2019) T. I. Tuegel, V. Chua, and T. L. Hughes, Phys. Rev. B 100, 115126 (2019).
  * Khalaf (2018) E. Khalaf, Phys. Rev. B 97, 205136 (2018).
  * Matsugatani and Watanabe (2018) A. Matsugatani and H. Watanabe, Phys. Rev. B 98, 205129 (2018).
  * Trifunovic and Brouwer (2019) L. Trifunovic and P. W. Brouwer, Phys. Rev. X 9, 011012 (2019).
  * Slager  _et al._ (2016) R.-J. Slager, V. Juričić, V. Lahtinen, and J. Zaanen, Phys. Rev. B 93, 245406 (2016).
  * Hughes  _et al._ (2011) T. L. Hughes, E. Prodan, and B. A. Bernevig, Phys. Rev. B 83, 245132 (2011).
  * Zak (1989) J. Zak, Phys. Rev. Lett. 62, 2747 (1989).
  * Hofstadter (1976) D. R. Hofstadter, Phys. Rev. B 14, 2239 (1976).
  * Bernevig  _et al._ (2006b) B. A. Bernevig, T. L. Hughes, S.-C. Zhang, H. D. Chen, and C. Wu, International Journal of Modern Physics B 20, 3257 (2006b).
  * Bernevig  _et al._ (2007) B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Solid State Communications 143, 20 (2007).
  * Lin and Hughes (2018) M. Lin and T. L. Hughes, Phys. Rev. B 98, 241103(R) (2018).
  * Wang  _et al._ (2020) H.-X. Wang, Z.-K. Lin, B. Jiang, G.-Y. Guo, and J.-H. Jiang, Phys. Rev. Lett. 125, 146401 (2020).
  * Ghorashi  _et al._ (2020) S. A. A. Ghorashi, T. Li, and T. L. Hughes, Phys. Rev. Lett. 125, 266804 (2020).
  * Simon  _et al._ (2021) S. Simon, M. Geier, and P. W. Brouwer, arXiv preprint arXiv:2109.02664 (2021).
  * Wang  _et al._ (2014) H. Wang, R. Wang, Y. Pan, L. Sheng, B. G. Wang, and D. Y. Xing, EPL 107, 40007 (2014).

##  Appendix A Band Dispersion near Fermi Level vs. Stacking Fault Coupling
Strength for Embedded Metal in (1+1)-D

In this section, we present additional figures of the band dispersion near the
Fermi level as a function of the stacking fault coupling strength. When the
stacking fault procedure is performed, the exposed chains are taken to have
generic coupling strength $\lambda$ where
$0\leq\lambda\leq\max(\gamma,\delta)$. Figs. 10 (a)-(d) provide a full
illustration of the evolution of the band dispersion as $\lambda$ is tuned. As
$\lambda$ is increased, the edges of the conduction and valence bands giving
rise to the nodes at $k_{x}=\pm\frac{\pi}{2}$ approach each other. This is
also accompanied by the fact that the dispersion of the nodes changes as
$\lambda$ approaches $\max(\gamma,\delta)$. Namely, the slope of the
dispersion near the Fermi level (i.e., the Fermi velocity) decreases as
$\lambda$ increases.

Figure 10: Plots of the band dispersion near the Fermi level ($E_{F}=0$)
versus the stacking fault coupling strength $\lambda$. We take $\gamma>\delta$
and use the same parameter values $\gamma=0.75$ and $\delta=0.5$ as in Fig. 2.
(a) Plot of the band dispersion near the Fermi level as a function of both
$\lambda$ and $k_{x}$. The bolded red lines highlight how the bulk band edges
bend towards each other as $\lambda$ is increased. (b)-(d) Plots of the band
structure for various values of $\lambda$. When $\lambda$ is small, the
dispersion of the nodes at $k_{x}=\pm\frac{\pi}{2}$ is linear. In particular,
when $\lambda=0$, numerical calculation of the slope shows that $\frac{\Delta
E}{t\Delta k_{x}}\approx 2.0000636$ (with lattice constant $a=1$) and for
$\lambda=0.03$ (shown in (b)), $\frac{\Delta E}{t\Delta k_{x}}\approx
1.99366$. The dispersion of the nodes remains linear in (c), but the slope of
the dispersion begins to shift noticeably. As $\lambda$ approaches the value
of $\gamma>\delta$, the dispersion of the nodes is flattened. The intermediate
value of $\lambda=\delta$ is illustrated in Fig. 2 b) in the main text.

We can gain some insight into the relationship between the band dispersion
near the Fermi level at $k_{x}=\pm\frac{\pi}{2}$ and the stacking fault
coupling strength $\lambda$ by considering the Hamiltonian given by Eq. (2) in
one of two “dimerized” limits: (i) $\gamma\neq 0$, $\delta=0$ and (ii)
$\gamma=0$, $\delta\neq 0$. Below, we consider the dimerized limit given by
(i), but the results are exactly the same for (ii) by simply replacing
$\gamma$ with $\delta$. Performing a $k\cdot p$ expansion, the dispersion near
$k_{x}=\pm\frac{\pi}{2}$ has the following form,

$\begin{split}E_{\pm}(\delta
k_{x})\approx\pm\left(\alpha(\gamma,\lambda)\delta
k_{x}+\beta(\gamma,\lambda)\delta k_{x}^{3}\right)\end{split}$ (10)

where $\delta k_{x}=k_{x}\pm\frac{\pi}{2}$ is the deviation in momentum from
the location of the nodes and the coefficients $\alpha(\gamma,\lambda)$ and
$\beta(\gamma,\lambda)$ are given by,

$\begin{split}&\alpha(\gamma,\lambda)=\frac{2t(\gamma^{2}-\lambda^{2})}{\gamma^{2}+\lambda^{2}}\\\
\\\
&\beta(\gamma,\lambda)=-\frac{t(\gamma^{8}+2\gamma^{6}\lambda^{2}-192t^{2}\gamma^{2}\lambda^{4}-2\gamma^{2}\lambda^{6}-\lambda^{8})}{3(\gamma^{2}+\lambda^{2})^{4}}\end{split}$
(11)

When the stacking fault coupling is very weak (i.e., $\lambda\ll\gamma$), the
coefficient $\alpha(\gamma,\lambda)$ is much greater than the coefficient
$\beta(\gamma,\lambda)$ and the dispersion near $k_{x}$ is dominated by the
first term in (LABEL:eq:kdotpdispersioncoefficients) and is linear. But as
$\lambda$ is increased and specifically, when $\lambda$ becomes very close in
value to $\gamma$, $\alpha(\gamma,\lambda)$ is significantly smaller than
$\beta(\gamma,\lambda)$. Therefore, at strong stacking fault coupling
$\lambda$ that is comparable to $\gamma$, the second term in (10) dominates
and leads to a cubic dispersion in energy near $k_{x}=\pm\frac{\pi}{2}$. In
fact, at $\lambda=\gamma$, the band dispersion is purely cubic. This leads to
the following relationship between the stacking fault coupling $\lambda$ and
the magnitude of the group velocity near the Fermi level $v_{F}$,

$\begin{split}v_{F}=\left|\frac{\delta E}{\delta
k_{x}}\right|\approx\alpha(\gamma,\lambda)+\beta(\gamma,\lambda)\delta
k_{x}^{2}\end{split}$ (12)

which shows that at weak $\lambda$, the group velocity is constant and at
strong $\lambda$, the group velocity has a quadratic dispersion in momentum.

##  Appendix B Topological Characterization of ETSM from Entanglement Spectrum

In this section, we provide an alternative means for identifying the
topological characteristics of the embedded Dirac semimetal utilizing the
single-particle entanglement spectrum. Using the same procedure as in Subsec.
III.3, we analyze entanglement Berry phase characteristics across the 2D BZ.

###  B.1 Brief Review of Embedded Topological Insulators (ETIs) and
Entanglement Spectrum (ES)

To motivate this section and understand why the entanglement spectrum is a
useful tool for diagnosing ETSMs, we briefly review the concept of an Embedded
Topological Insulator (ETI) from Ref. Tuegel  _et al._ , 2019. Consider an
ordinary (topologically trivial) environment $\mathcal{M}$ and a topological
insulator (TI) $\mathcal{N}$ embedded in $\mathcal{M}$, and consider the case
where $\text{dim}(\mathcal{M})>\text{dim}(\mathcal{N})$. Nominally, if a
physical cut is made between the TI $\mathcal{N}$ and its trivial environment
$\mathcal{M}$, topologically protected modes will emerge, such as surface
states. However, if $\mathcal{N}$ is coupled to its environment $\mathcal{M}$,
it becomes non-trivial to diagnose the topology of $\mathcal{N}$. At weak
coupling it is possible to maintain these surface modes as long as
$\mathcal{N}$ remains gapped and symmetries are preserved. However, at strong
coupling $\mathcal{N}$ could potentially no longer be gapped, destabilizing
the TI $\mathcal{N}$. In this latter case, the topology of $\mathcal{N}$ can
be determined by disentangling $\mathcal{N}$ from its environment
$\mathcal{M}$ through the use of the entanglement spectrum.
For non-interacting fermionic systems harboring ETIs, the entanglement
spectrum is derived as follows. First, a region $\mathcal{R}$ is identified as
harboring the presumed ETI, and an entanglement cut is performed between
$\mathcal{R}$ and its complementary region $\mathcal{R}^{c}$. This leads to a
1-body reduced density matrix $\rho_{\mathcal{R}}$ is determined by
restricting the projector to the region $\mathcal{R}$ (i.e.,
$\rho_{\mathcal{R}}=P|_{\mathcal{R}}$). The entanglement spectrum (ES) is the
eigenvalue spectrum $\\{\lambda_{i}\\}$ of the reduced density matrix
$\rho_{\mathcal{R}}$ and is distributed over the unit interval
$\\{\lambda_{i}\\}\subset[0,1]$. If the ES is gapped at $\lambda=\frac{1}{2}$,
the set of eigenstates of $\rho_{\mathcal{R}}$ with $\lambda_{i}>\frac{1}{2}$
are identified as having the dominant contribution to correlations within
$\mathcal{R}$ (i.e., “occupied” entanglement states). The spectral projector
is then constructed as

$\begin{split}\Pi_{\mathcal{R}}=\sum\limits_{\lambda_{i}>\frac{1}{2}}|\lambda_{i}\rangle\langle\lambda_{i}|.\end{split}$
(13)

If the topological index computed from $\Pi_{\mathcal{R}}$ is non-zero, then
$\mathcal{R}$ is said to be topologically non-trivial. Furthermore, if any
(symmetry-preserving) entanglement cuts made through $\mathcal{R}$ leads to a
gapless entanglement spectrum, then $\mathcal{R}$ is said to harbor an ETI.

###  B.2 Entanglement Berry Phase Analysis of the Embedded Dirac Semimetal

We now compute the entanglement spectrum of the embedded topological Dirac
semimetal discussed in Sec. III. First, we take the spectral projector to be
$P(k_{x},k_{y})=\sum\limits_{E_{n}(k_{x},k_{y})<E_{F}}|u_{n}(k_{x},k_{y})\rangle\langle
u_{n}(k_{x},k_{y})|$ where the Fermi level is taken to be $E_{F}=0$. Before
proceeding, we note that the band structure for the embedded Dirac semimetal
is gapless precisely at $(k_{x},k_{y})=\left(\pm\pi,\pm\frac{\pi}{2}\right)$,
which means our construction of the spectral projector at these points is not
well-defined. Therefore, we modify our construction of the spectral projector
by specifically avoiding the points where the band structure is gapless, i.e.

$\begin{split}P(k_{x},k_{y})=\sum\limits_{E_{n}(k_{x},k_{y})<E_{F}}|u_{n}(k_{x},k_{y})\rangle\langle
u_{n}(k_{x},k_{y})|.\end{split}$ (14)

The sum is performed over all
$(k_{x},k_{y})\neq\left(\pi,\pm\frac{\pi}{2}\right)$. With this definition of
the spectral projector, we are ensuring that we consider only the eigenstates
with energies strictly below the Fermi level and the calculation takes place
in only the gapped regions of the band structure. This is crucial in
evaluating the entanglement spectrum. From here, we construct the reduced
density matrix $\rho_{\mathcal{R}}$ by restricting $P$ to the region
$\mathcal{R}$ illustrated in Fig. 11 (a). The entanglement cuts are made
directly above and below the layer where the stacking fault is located, which
are above and below the $(N_{z}-1)^{\text{th}}$ layer specifically following
the stacking fault procedure detailed in Subsec. III.2. After computing the
eigenvalue spectrum $\\{\lambda_{i}\\}$ of $P|_{\mathcal{R}}$, we construct
the spectral projector over the entanglement eigenstates according to (13).
The resulting entanglement ‘band structure’ constructed from the eigenvalues
$\\{\lambda_{i}\\}$ is shown in Fig. 11 (b) and is clearly gapless, indicating
that the region $\mathcal{R}$ we chose is harboring an ETSM.
Just as was outlined in Subsec. III.3, we divide the 2D BZ into 1D layers
extended along the $k_{x}$ direction and treat $k_{y}$ as a layering parameter
(we also note that a similar method for computing the entanglement spectrum in
Weyl semimetals was also developed in Ref. Wang  _et al._ , 2014). We evaluate
the entanglement Berry phase and entanglement Berry connection, which we
define as follows,

$\begin{split}\theta_{B,\text{ent}}(k_{y})=\int\limits_{-\pi}^{\pi}dk_{x}\hskip
2.84544pt\text{tr}\hskip 1.42271pt\mathcal{A}_{x}(k_{x})\end{split}$ (15)
$\begin{split}\mathcal{A}_{x}(k_{x})=\langle\lambda_{m}(k_{x},k_{y})|i\nabla_{k_{x}}|\lambda_{n}(k_{x},k_{y})\rangle\end{split}$
(16)

When implementing (15) and (16) in practice, we use the discretized gauge
invariant Wilson loop expression over the entanglement eigenstates, which is
generalized straightforwardly from Eq. (LABEL:eq:discretizedZakBerryphase) in
the main text. The result of this is illustrated in Fig. 11 (c), and the
results are identical to Fig. 6 in the main text. As $k_{y}$ is varied from
$(-\pi,\pi]$, the entanglement Berry phase is quantized throughout and only
undergoes a change in value at $k_{y}=\pm\frac{\pi}{2}$.

Figure 11: Characterizing the topology of the embedded Dirac semimetal from
the entanglement perspective. (a) illustrates the entanglement cuts and the
region $\mathcal{R}$ chosen after the stacking fault procedure is performed
(the layers above and below are faded out to emphasize the layer lying on the
stacking fault after the extraction). In (b), the gapless entanglement ‘band
structure’ is shown. In (c), the entanglement Berry phase as a function of
$k_{y}$ is shown, with the dashed lines indicating discontinuities in which
the phase sharply changes value. These occur precisely at
$k_{y}=\pm\frac{\pi}{2}$, precisely where the gapless points in the
entanglement ‘band structure’ occur. Identical model parameters were used to
construct the plots in (b) and (c) as in the main text.

##  Appendix C Details of Model for ETDSM in a Magnetic Field

In this section, we provide the technical details for the ETDSM model when
placed in a magnetic field. We introduce a magnetic field along the
$z$-direction by working in the Landau gauge $\mathbf{A}=Bx\hat{y}$. The
vector potential is incorporated into the Hamiltonian via the Peierls
substitution,

$\begin{split}&c_{n_{x}^{\prime},n_{y}^{\prime},n_{z}^{\prime}}^{\dagger}c_{n_{x},n_{y},n_{z}}\\\
&\to
c_{n_{x}^{\prime},n_{y}^{\prime},n_{z}^{\prime}}^{\dagger}c_{n_{x},n_{y},n_{z}}\exp\left\\{i\int\limits_{(n_{x},n_{y},n_{z})}^{(n_{x}^{\prime},n_{y}^{\prime},n_{z}^{\prime})}\mathbf{A}\cdot\hskip
1.42271ptd\mathbf{r}\right\\}\end{split}$ (17)

The magnetic field strength $B$ is taken to be
$B=\phi\left(\frac{\Phi_{0}}{a^{2}}\right)$ where $p$ and $q$ are co-prime
integers, $\Phi_{0}$ is the flux quantum, and $a$ is the lattice constant
along the $x$ and $y$ directions (we take $\Phi_{0}=a=1$). The flux per
plaquette in the $x-y$ plane is given by $\phi$, where $p$ prescribes the
integer number of flux quanta per magnetic unit cell of size $q$ along the
$x$-direction. Given the choice of Landau gauge, the standard translational
invariance along the $x$-direction is no longer present. Instead, a magnetic
translational invariance is present in which the Hamiltonian is invariant
under translations of the form $n_{x}\to n_{x}+q$. Standard translational
invariance is maintained along the $y$-direction under this choice of Landau
gauge. Performing the Peierls substitution and taking a Fourier transform
along the $y$ direction, the stacking fault Hamiltonian is now

$\begin{split}H=H_{1}+H_{2}+H_{3}\end{split}$ (18)

where,

$\begin{split}&H_{1}=\sum\limits_{n_{x}}\sum\limits_{k_{y}}\sum\limits_{n_{z}=1}^{N_{z}-1}\left(c_{n_{x},k_{y},n_{z}}^{\dagger}[H_{\text{on-
site}}(n_{x},k_{y},\phi)]c_{n_{x},k_{y},n_{z}}\right.\\\
&\left.\vphantom{\sum\limits_{n_{x}}\sum\limits_{k_{y}}}+c_{n_{x}+1,k_{y},n_{z}}^{\dagger}[H_{+\hat{x}}]c_{n_{x},k_{y},n_{z}}+\text{h.c.}\right)\end{split}$
(19)
$\begin{split}&H_{2}=\sum\limits_{n_{x}}\sum\limits_{k_{y}}\sum\limits_{n_{z}=1}^{\frac{N_{z}}{2}-1}\left(i\gamma_{z}c_{n_{x},k_{y},2n_{z}}^{\dagger}c_{n_{x},k_{y},2n_{z}-1}\right.\\\
&\left.\vphantom{\sum\limits_{n_{x}}\sum\limits_{k_{y}}}+i\delta_{z}c_{n_{x},k_{y},2n_{z}+1}^{\dagger}c_{n_{x},k_{y},2n_{z}}+\text{h.c.}\right)\end{split}$
(20)
$\begin{split}&H_{3}=\sum\limits_{n_{x}}\sum\limits_{k_{y}}\sum\limits_{n_{z}=\frac{N_{z}}{2}}^{N_{z}-1}\left(i\gamma_{z}c_{n_{x},k_{y},2n_{z}+1}^{\dagger}c_{n_{x},k_{y},2n_{z}}\right.\\\
&\left.\vphantom{\sum\limits_{n_{x}}\sum\limits_{k_{y}}}+i\delta_{z}c_{n_{x},k_{y},2n_{z}}^{\dagger}c_{n_{x},k_{y},2n_{z}-1}+\text{h.c.}\right)\end{split}$
(21)

where we have introduced periodic boundaries in the $x$-direction (i.e.,
$n_{x}+N_{x}\equiv n_{x}$) and open boundaries in the $z$ direction. In
(LABEL:eq:stackingfaulthamiltonian21), we have introduced the following
quantities,

$\begin{split}[H_{\text{on-site}}(n_{x},k_{y},\phi)]=(1-m-t_{y}\cos(2\pi\phi
n_{x}-k_{y}))\sigma_{z}\end{split}$ (22)
$\begin{split}[H_{+\hat{x}}]=\frac{i\sigma_{y}-\sigma_{z}}{2}\end{split}$ (23)

Note that in
(LABEL:eq:stackingfaulthamiltonian21)-(LABEL:eq:stackingfaulthamiltonian23),
in order for magnetic translational invariance to hold, we must have $N_{x}$
be divisible by the magnetic unit cell size $q$ (i.e., $N_{x}=lq$ for some
$l\in\mathbb{Z}$). Furthermore, for numerical calculation purposes, we take
the values of $k_{y}$ to be over a discrete grid with spacing
$\frac{2\pi}{N_{y}}$ such that
$k_{y}\in\left\\{-\pi\left(1-\frac{2}{N_{y}}\right),\ldots,\pi\right\\}$.