# Spirals and skyrmions in antiferromagnetic triangular lattices

Wuzhang Fang  Aldo Raeliarijaona  Po-Hao Chang  Alexey A. Kovalev  Kirill D.
Belashchenko  Department of Physics and Astronomy and Nebraska Center for
Materials and Nanoscience, University of Nebraska-Lincoln, Lincoln, Nebraska
68588-0299, USA

(August 27, 2024)

###### Abstract

We study realizations of spirals and skyrmions in two-dimensional
antiferromagnets with a triangular lattice on an inversion-symmetry-breaking
substrate. As a possible material realization, we investigate the adsorption
of transition-metal atoms (Cr, Mn, Fe, or Co) on a monolayer of
$\text{MoS}_{2}$, $\text{WS}_{2}$, or $\text{WSe}_{2}$ and obtain the
exchange, anisotropy, and Dzyaloshinskii-Moriya interaction parameters using
first-principles calculations. Using energy minimization and parallel-
tempering Monte-Carlo simulations, we determine the magnetic phase diagrams
for a wide range of interaction parameters. We find that skyrmion lattices can
appear even with weak Dzyaloshinskii-Moriya interactions, but their stability
is hindered by magnetic anisotropy. Our results suggest that
$\text{{Cr}/MoS}_{2}$, $\text{{Fe}/MoS}_{2}$, and $\text{{Fe}/WSe}_{2}$
interfaces can host spin spirals formed from the 120$\degree$
antiferromagnetic states. Our results further suggests that for other
interfaces, such as $\text{{Fe}/MoS}_{2}$, the Dzyaloshinskii-Moriya
interaction is strong enough to drive the system into a three-sublattice
skyrmion lattice in the presence of experimentally feasible external magnetic
field.

##  I Introduction

Antiferromagnetic spintronics [1] is an actively developing field with many
potential applications as it explores materials with unique properties, i.e.,
materials associated with lack of stray fields and ultrafast dynamics in THz
range [2]. Concepts of topology have influenced the development of condensed
matter physics, e.g., as can be seen in realizations of magnetic skyrmions in
ferromagnets [3; 4; 5]. Realizations of skyrmions in collinear
antiferromagnets (AFMs) are also being explored [6; 7; 8; 9; 10]. Studies of
three-sublattice systems generalize above ideas to noncollinear
antiferromagnets and demonstrate possibilities for realizations of various
skyrmion phases [11; 12; 13; 14; 15]. It has recently been demonstrated that
the topological charge in AFMs can be also fractionalized [15; 10] which can
be interpreted as a formation of merons in square crystals of vortices and
antivortices [16; 17; 18; 19].

Materials engineering has been used for realizations of room temperature
skyrmions in ferromagnets [20; 21; 22]. In this work, we explore possibilities
of integrating transition metal dichalcogenides (TMDs) with magnetic
transition metals to realize noncollinear AFMs capable of hosting skyrmion
lattices. Studies have demonstrated utility of TMDs for realizations of spin-
valves [23; 24] and TMD-based spintronic devices [25; 23; 24; 26]. Contact
between metals and TMDs can alter the electronic properties of the system and
affect the efficiency of TMD-based devices [27; 25; 28]. It is thus crucial to
study various proximity effects associated with the metal/TMD interfaces [29].
Due to the layered nature of TMDs, capping them with magnetic atoms seems to
be natural in order to bring about novel properties. In a recent study of
Fe/Ir(111) and TMD interface, it has been found that transition-metal atoms in
the TMD monolayer acquire a metallic character and spin imbalance due to its
chemisorption on Fe/Ir [27]. The large spin-orbit coupling (SOC) in monolayer
TMDs [30] combined with structural symmetry breaking in magnetic
heterostructures can potentially lead to sizable Dzyaloshinskii-Moriya
interaction (DMI) [31; 32]. It is thus crucial to investigate magnetic
interactions in TMD-based magnetic systems and determine such phenomenological
parameters as the strength of exchange interaction, the strength of DMI, and
the single-ion magnetic anisotropy. This will also determine the feasibility
of realizing various skyrmion phases in TMD-based triangular AFMs.

A very strong DMI in a triangular AFM can lead to appearance of a three-
sublattice skyrmion crystal [11]; however, it is not clear whether such strong
DMI can be realized in real materials. It has been suggested that a relatively
large DMI, $D>0.2J$ with $J$ being the exchange energy, is required for
stabilization of skyrmions in triangular AFMs [33; 34], which puts a strong
constraint on realizations of three-sublattice skyrmons. Numerically, it has
been demostrated that skyrmion lattices in triangular AFM can be stabilized
even for weak DMI, $D<0.2J$ [35]. In real materials, the presence of SOC will
likely also lead to the apearance of the single-ion magnetic anisotropy.
Despite the existing studies of three-sublattice skyrmions in triangular AFMs,
the analysis of such skyrmions for the experimentaly feasible DMI and in the
presence of single-ion magnetic anisotropy is still lacking.

In this paper, we study realizations of magnetic spirals and skyrmions in
triangular AFMs in the presence of DMI, comparable to or smaller than the
critical DMI strength, $D_{c}=0.2J$ [33], and include the effect of single-ion
magnetic anisotropy. As a possible material realization, we consider one layer
of magnetic transition metal X, such that X={Co, Cr, Fe, Mn}, adsorbed onto a
monolayer of transition metal dichalcogenides ($\text{MoS}_{2}$,
$\text{WS}_{2}$, or $\text{WSe}_{2}$). From ab initio calculations, we
evaluate parameters describing magnetic interactions, namely the exchange
interaction, the uniaxial anisotropy, and the DMI, by computing the difference
in total energy of magnetic configurations. We find that {Co, Cr, Fe,
Mn}/MoS2, $\text{{Mn}/WS}_{2}$, and $\text{{Fe}/WSe}_{2}$ can be described as
AFMs with uniaxial anisotropy and DMI. Uisng Monte Carlo simulations, we study
phase diagrams of triangular AFMs in the presence of DMI and magnetic
anisotropy, and assess the feasibility of realizing magnetic spiral and
skyrmion phases. Our Monte Carlo simulations indicate that skyrmion lattices
will appear even for relatively weak DMI below $D_{c}=0.2J$, but their
stability will be hindered by magnetic anisotropy. We find that in the absence
of magnetic field the magnetic ground state of $\text{{Cr}/MoS}_{2}$,
$\text{{Fe}/MoS}_{2}$, and $\text{{Fe}/WSe}_{2}$ is a spin spiral on three
sublattices of the triangular lattice while the ground state of
$\text{{Mn}/MoS}_{2}$, $\text{{Mn}/WS}_{2}$, and $\text{{Co}/MoS}_{2}$ is the
120$\degree$ Néel type configuration. Our analysis further suggests that
$\text{{Fe}/MoS}_{2}$ can potentially host a three-sublattice skyrmion crystal
in the presence of experimentally feasible magnetic fields.

##  II Methods

In this section, we briefly describe methods used in analyzing triangular
AFMs. The phenomenological parameters for several material candidates are
established by employing first-principles calculations. The zero temperature
phase diagrams are obtained by minimizing the energy for a wide range of
phenomenological parameters. The finite temperature phase diagrams are
obtained by employing the feedback-optimized parallel tempering Monte Carlo
simulations.

### First-principles calculations

First-principles calculations are performed using the Vienna ab initio
simulation package (VASP) [36] and open source package for Material eXplorer
(OpenMX) codes [37]. In VASP calculations, we use the generalized gradient
approximation [38] and projector augmented wave potentials [39; 40]. A
$27\times 27\times 1$ k-point mesh is used to sample the Brillouin zone. The
dense mesh of k-point is necessary to converge to the correct magnetic state
and to get accurate total energy. The plane-wave cutoff energy is 500 eV and
energy convergence criteria are $10^{-5}$ eV for structure relaxations and
$10^{-6}$ eV for total energy calculations. The structures are fully relaxed
until the magnitude of the force on the individual atom is less than 0.01
eV/Å. A thick vacuum layer of 17 Å along the out-of-plane direction is used to
minimize the interaction due to the periodic condition. To describe the
interactions between localized d electrons of transition atoms accurately, we
take into account the on-site Coulomb interactions using the LDA+U method
implemented in VASP [41] and OpenMX [42]. The effective on-site Coulomb
parameters $U$ is calculated for each transition metal atom using the linear
response method, as implemented in VASP, and the effective exchange parameter
$J$ is chosen as 0.9 eV. After determining the adsorption sites and relaxed
structures using VASP, we use OpenMX for the calculations of physical
properties of our systems. OpenMX uses pseudo-atomic localized orbitals [43]
as basis functions. Band structures have been compared to the ones calculated
using VASP to determine the basis sets used in our calculations.

### Model Hamiltonian

The monolayer of adsorbed transition atoms on top of $\text{MoS}_{2}$,
$\text{WS}_{2}$ or $\text{WSe}_{2}$ form a triangular lattice. The neighboring
triangles are not equivalent due to the existence of the transition metal, and
this asymmetry results in the Dzyaloshinskii-Moriya interactions (DMI) [44;
45]. The model Hamiltonian to describe the magnetic interactions between the
adsorbed transition atoms is given by:

$\displaystyle H$ $\displaystyle=$
$\displaystyle\sum_{\braket{ij}}{J\,\mathbf{S}_{i}\cdot\mathbf{S}_{j}+\tilde{\mathbf{D}}_{ij}\cdot(\mathbf{S}_{i}\times\mathbf{S}_{j})}$
(1) $\displaystyle-$
$\displaystyle\tilde{K}\sum_{i}{(S_{i}^{z})^{2}}-\tilde{h}\sum_{i}{S_{i}^{z}},$

where $i$ and $j$ are site indices with $\braket{ij}$ indicating the nearest
neighbors; for each site we introduce a spin variable $\mathbf{S}_{i}$ with
$|\mathbf{S}_{i}|=1$, $J$ describes the exchange interaction, $\tilde{K}$ is
the single-ion anisotropy constant, $\tilde{h}=\mu B$ describes Zeeman
coupling of the magnetic moment $\mu$ with the external magnetic field $B$,
and
$\tilde{\mathbf{D}}_{ij}=\tilde{\mathbf{D}}_{\perp}+\tilde{\mathbf{D}}_{\parallel}$
is the DMI vectors, as shown in Fig. 1, split into out-of-plane,
$\tilde{\mathbf{D}}_{\perp}$, and in-plane, $\tilde{\mathbf{D}}_{\parallel}$,
contributions. We introduce the reduced parameters
$\mathbf{D}_{ij}=\tilde{\mathbf{D}}_{ij}/J$, $h=\tilde{h}/J$, $K=\tilde{K}/J$,
and $T=k_{B}\tilde{T}/J$, which we will use except when we refer to parameters
determined from first-principles.

### Energy minimization

To find the classical ground state in the absence of magnetic field at zero
temperature, we use energy minimization. This can be achieved by using the
spherical approximation within the Luttinger-Tisza method [46]. We replace the
constraint $|\mathbf{S}_{i}|=S$ by the milder constraint
$\sum_{i}|S_{i}|^{2}=NS^{2}$ where $N$ is the number of lattice sites. As an
alternative approach, we also directly substitute a spiral ansatz in Eq. (1)
and perform numerical minimization, which slightly improves on the spherical
approximation. The model Hamiltonian in Eq. (1) is rewritten using the Fourier
transform (see Appendix A) of the spins
${\mathbf{S}_{\mathbf{r},a}=\sum_{\mathbf{q}}\mathbf{S}_{\mathbf{q},a}e^{i\mathbf{q\cdot
r}}}$ as:

$H=\sum_{\mathbf{q}}{\mathbf{\Psi}_{\mathbf{-q}}H_{\mathbf{q}}\mathbf{\Psi}_{\mathbf{q}}},$
(2)

where $H_{\mathbf{q}}$ is the Fourier transformed Hamiltonian, and the vector
$\mathbf{\Psi_{q}}$ is given by:

$\mathbf{\Psi_{q}}=\\{\mathbf{S}_{\mathbf{q},1},\mathbf{S}_{\mathbf{q},2},\mathbf{S}_{\mathbf{q},3}\\}.$
(3)

Within the Luttinger-Tisza formalism, the ground state is determined by
solving the eigenvalue problem
$H_{\mathbf{q}}\Psi_{\mathbf{q}}=\varepsilon_{\mathbf{q}}\Psi_{\mathbf{q}}$
where the constraint leads to the normalization of the eigenvectors [46]. To
improve this procedure, we also numerically minimize Eq. (2) by directly
substituting a spiral magnetic structure.

Figure 1: (a) Three adsorption sites (H, T${}_{\text{S}}$ and
T${}_{\text{Mo}}$) on top of monolayer $\text{MoS}_{2}$. Mo and S atoms are
represented by blue and orange spheres, respectively; (b) the same unit cell
with the Dzyaloshinskii-Moriya vector depicted by the black arrows centered on
the bond of neighboring, interacting atoms. The convention taken for the
direction on bonds is counterclockwise in a triangular plaquette with Mo atom
at its center. The vectors a and b are the primitive vectors of the hexagonal
Bravais lattice.

### Parallel tempering Monte Carlo

To obtain the phase diagram and perform simulated annealing, we use feedback-
optimized parallel tempering Monte Carlo simulations with Metropolis updates
[47; 48]. We used up to 300 spin configurations (replicas), simultaneously
simulated at different temperatures. The temperatures are chosen according to
optimization algorithm improving the exchanges of states of different replicas
[49], see Appendix B for details. The system consists of $N=L^{2}$ sites for
$L=75$ with periodic boundaries (in some cases we considered larger systems
increasing the systems size to $L=150$). The system was equilibrated using
runs of $10^{6}-10^{7}$ Monte-Carlo sweeps (attempts per spin) and
consequently the same number of sweeps was used to accumulate statistics. By
considering different lattice sizes, we observe finite size effects; however,
we observe that our results capture all important features of the phase
diagram given the large size of the system. The magnetic ground state is
characterized using the real space spin texture in combination with the spin
structure factor defined as:

$\displaystyle
S_{\perp}(\mathbf{q})=\frac{1}{N}\bigg{\langle}\sum_{\alpha=x,y}{\left|\sum_{\mathbf{r}}{S^{\alpha}_{\mathbf{r}}\exp^{-i\mathbf{q}\cdot\mathbf{r}}}\right|^{2}}\bigg{\rangle},$
(4) $\displaystyle
S_{\parallel}(\mathbf{q})=\frac{1}{N}\bigg{\langle}\left|\sum_{\mathbf{r}}{S^{z}_{\mathbf{r}}\exp^{-i\mathbf{q}\cdot\mathbf{r}}}\right|^{2}\bigg{\rangle},$
(5)

where $\langle\rangle$ stands for statistical averaging, $N$ is the number of
spins, and the perpendicular or parallel components to the plane normal are
considered. In addition, we introduce the sublattice topological charge [50]
by summing contributions from all elementary triangles of sublattice $\alpha$:

$Q_{\alpha}=\frac{1}{2\pi}\sum_{\Delta}\operatorname{Arg}(A_{\alpha,\Delta}+iB_{\alpha,\Delta}),$
(6)

where

$\displaystyle A_{\alpha,\Delta}$ $\displaystyle=$
$\displaystyle\mathbf{S}_{\alpha,1}\cdot\left(\mathbf{S}_{\alpha,2}\times\mathbf{S}_{\alpha,3}\right),$
(7) $\displaystyle B_{\alpha,\Delta}$ $\displaystyle=$ $\displaystyle
1+\mathbf{S}_{\alpha,1}\cdot\mathbf{S}_{\alpha,2}+\mathbf{S}_{\alpha,2}\cdot\mathbf{S}_{\alpha,3}+\mathbf{S}_{\alpha,1}\cdot\mathbf{S}_{\alpha,3},$
(8)

with $\mathbf{S}_{\alpha,1}$, $\mathbf{S}_{\alpha,2}$, and
$\mathbf{S}_{\alpha,3}$ being spins located on vertices of elementary triangle
of lattice $\alpha$. The total topological charge over all three sublattices
is defined as $Q=\sum_{\alpha=1}^{3}Q_{\alpha}$. To identify the phase
boundaries, we calculate the specific heat:

$\displaystyle
C_{\text{V}}=\frac{1}{T^{2}}\left(\braket{E^{2}}-\braket{E}^{2}\right),$ (9)

and the topological susceptibility:

$\chi_{{\pm}}=\frac{\langle(Q_{\pm})^{2}\rangle-\langle
Q_{\pm}\rangle^{2}}{\langle Q_{\pm}\rangle T},$ (10)

where we define skyrmion-like, $Q_{+}$, and antiskyrmion-like, $Q_{-}$,
topological charge [51; 52], k${}_{\text{B}}$ is the Boltzmann constant, and
$T$ is the reduced temperature. Note that only triangles with positive charge
density contribute to $Q_{+}$ and only triangles with negative charge density
contribute to $Q_{-}$, so that $Q=Q_{+}-Q_{-}$. In the presentation of our
results, we use $\chi_{-}$ as the $\chi_{-}$ peak is higher than the
$\chi_{+}$ peak. According to our results, the positions of both peaks are
almost identical.

Figure 2: Phase diagram as a function of DMI and effective anisotropy of AFM
triangular lattice at zero temperature in the absence of magnetic field. The
horizontal axis indicates the in-plane DMI and the vertical axis indicates the
effective uniaxial anisotropy. The plot shows the ground states of (Cr, Fe,
Co, $\text{{Mn)}/MoS}_{2}$, $\text{{Mn}/WS}_{2}$ and $\text{{Fe}/WSe}_{2}$.
The colors mark the different phases, i.e., the light blue region corresponds
to the Néel phase and the pink region is the spiral phase. The plotting range
was chosen to distinctly display all 6 cases in Table 1. Figure 3: A snapshot
of spin spiral in a triangular AFM obtained by simulated annealing for
$K_{\text{eff}}=0.4$, $D_{\parallel}=0.2$, $h=0$, and $T=0.02$ on a $75\times
75$ triangular lattice with periodic boundary conditions. The letters on the
left side (A, B, C) indicate the plotted sublattice. The colors indicate the
orthogonal to the plane component of the spins.

##  III Results and discussion

In this section, we first discuss results of first-principles calculations
establishing adsoption sites of transition metal atoms on top of the monolayer
TMD. These calculations further establish phenomenological parameters
describing magnetic interactions for several material candidates. Finally, we
discuss zero and finite temperature phase diagrams for a wide range of
phenomenological parameters describing triangular AFMs. Particular emphasis is
given to realizations of SkX in the presence of weak DMI and magentic
anisotropy.

### Adsorption sites, electronic structure, and magnetic ordering

As shown in Fig. 1, three possible adsorption sites of transition-metal atoms
on monolayer TMD are denoted as H (hollow site), T${}_{\text{S}}$ (on top of
S), and T${}_{\text{TM}}$ (on top of TM, where TM is either Mo or W). We
consider the case where all adsorption sites of a given type are occupied. The
calculated $U$ values and preferred adsorption sites are listed in Table 1.

Table 1: Parameters in the model Hamiltonian. Atom | Adsorption site |  $U$ (eV) |  $J$ (meV) | T${}_{\text{N}}$ (K) |  $\tilde{\text{D}}_{\perp}$ (meV) |  $\tilde{\text{D}}_{\parallel}$ (meV) |  $K$ (meV) |  $K_{\text{eff}}$ (meV)
---|---|---|---|---|---|---|---|---
Cr/MoS2 | On top of S | 4.9 | 37.04 | 430 | 0.003 | 0.08 | 0.18 | 0.194
Mn/MoS2 | On top of S | 4.4 | 19.26 | 224 | 0.39 | 0.10 | -0.05 | 1.952
Fe/MoS2 | Hollow | 4.4 | 15.30 | 178 | -0.21 | 1.05 | 1.65 | 0.539
Co/MoS2 | On top of S | 4.9 | 3.30 | 38 | 0.008 | 0.03 | -0.11 | -0.066
Fe/WSe2 | Hollow | 4.4 | 12.45 | 144 | 0.68 | 1.63 | 0.86 | 4.42
Mn/WS2 | On top of S | 4.4 | 23.04 | 267 | 1.03 | 0.24 | -0.88 | 6.23

The band structures and densities of states (DOS) in the ferromagnetic state
for each system are shown in Fig. 4. Red color represents Mo or W character;
blue and green depict majority and minority-spin states of the magnetic
transition-metal atom. In Co/MoS2 and Cr/MoS2, the bands crossing the Fermi
level have a mixed character deriving from the hybridization of the
conduction-band states of monolayer MoS2 with the states of Co or Cr. These
bands have a strongly $k$-dependent exchange splitting reflecting the varying
weight of the Co or Cr states. In systems with Mn the Fermi level lies near
the bottom of additional TMD-derived bands, while in systems with Fe there are
almost pure minority-spin Fe bands at the Fermi level.

Figure 4: Projected band structures and densities of states of (a) Co/MoS2 (b)
Cr/MoS2 (c) Mn/MoS2 (d) Mn/WS2 (e) Fe/MoS2 (f) Fe/WSe2. The red, blue, and
green colors in the band-structure plots represent the partial weights of Mo
or W states, the majority-spin states of the magnetic atom, and the minority-
spin states of the magnetic atom, respectively. Black lines in the DOS plots
show the total DOS; blue and green, respectively, majority-spin and minority-
spin partial DOS for the magnetic atom.

Next, we determine the preferred magnetic ordering of the adsorbed atoms on
the adsorption sites. Ferromagnetic and $120\degree$ in-plane
antiferromagnetic orderings are considered in our calculations. Total energies
calculated with these two orderings suggest that the $120\degree$ in-plane
antiferromagnetic ordering is the preferred magnetic state with the lowest
total energy.

To extract the parameters of the model Hamiltonian (1), we perform constrained
DFT calculations with suitably specified non-collinear spin configurations,
which are described in Appendix C. The resulting parameters are listed in
Table 1.

### Spiral classical ground state at zero temperature

In the absence of DMI and external magnetic field, the ground state is AFM
with a 120$\degree$ three-sublattice magnetic structure [53; 54; 11]. The
spins in each sublattice separately are aligned and the angle between
magnetizations from each sublattice is 120$\degree$. In the presence of DMI,
Eq. (2) has minimum at nonzero wave vector $\mathbf{k}_{\text{min}}$. The
ground state now contains three spirals with wave vector
$\mathbf{k}_{\text{min}}$ set on each sublattice of the triangular lattice. It
is convenient to introduce the effective anisotropy K${}_{\text{eff}}$ that
combines the effect of the uniaxial anisotropy $K$ and D⟂ in Eq. (1), where
K${}_{\text{eff}}>0$ prefers the easy axis alignment, whereas
K${}_{\text{eff}}<0$ prefers the easy plane alignment. By considering the
energy difference between the out-of-plane and in-plane coplanar 120$\degree$
configurations, we can derive the expression of $K_{\text{eff}}$ to be:

$K_{\text{eff}}=K+D_{\perp}(3\sqrt{3}).$ (11)

The phase diagram resulting from energy minimization in Eq. (2) is depicted in
Fig. 2. By comparing energies for different magnetic structures, we
distinguish two types of ground states: the three-sublattice Néel states
(marked as the light blue region in Fig. 2) in which the spins on the same
sublattice are collinear and a spiral state (marked as the pink region in Fig.
2). To confirm the phase diagram in Fig. 2, we also performed simulated
annealing calculations. The resulting ground states have been characterized by
observing the real space spin textures and the spin structure factors in the
reciprocal space in Figs. 3 and 5.

The phase diagram in Fig. 2 shows that for 3 of the 6 compounds in Table 1,
$\text{{Cr}/MoS}_{2}$, $\text{{Fe}/MoS}_{2}$, and $\text{{Fe}/WSe}_{2}$, the
magnetic ground state is a spin spiral. As an example, in Figs. 3 and 5(a) we
show the real space spin textures and the spin structure factors for the case
$K_{\text{eff}}=0.4$, $D_{\parallel}=0.2$. Note that additional (substantially
smaller in height) peaks in Fig. 5(a) arise due to the fact that the spin
structure factor is calculated by performing the statistical averaging at low
but finite temperature, $T=0.02$. The dominant peak clearly identifies the
spiral phase in Fig. 5(a).

Figure 5: Structure factor $S_{\parallel}(\mathbf{q})$ of a triangular AFM
from Monte Carlo simulations. The vertical and horizontal axes represent the
reciprocal space coordinates $k_{\text{x}}$ and $k_{\text{y}}$ in units of
$1/L$. The plot is for sublattice A; those for B and C are almost identical.
(a) The plot corresponds to spiral phase realizable for $K=0.4$,
$D_{\parallel}=0.2$, $h=0$, and $T=0.02$. (b) The plot corresponds to SkX
phase realizable for $K=0$, $D_{\parallel}=0.2$, $h=3.1$, and $T=0.02$.

Figure 6: (a) and (b) Specific heat and topological susceptibility as a
function of temperature for different values of the external magnetic field,
$D_{\parallel}=0.2$, $K=0$. (c) and (d) Specific heat and topological
susceptibility as a function of temperature for different values of the
external magnetic field, $D_{\parallel}=0.2$, $K=0.1$.

### Antiferromagnetic skyrmion crystals

In this section, we show that it is possible to stabilize antiferromagnetic
skyrmion crystals (SkXs) at a finite temperature and in the presence of a
magnetic field. To this end, we consider a $75\times 75$ triangular lattice
with periodic boundary conditions and perform parallel tempering Monte Carlo
simulations. To identify the phase transitions, we calculate the specific
heat, the topological susceptibility, and the total topological charge
according to Eqs. (9), (10), and (6), respectively. Some of the curves used
for identifying phase transitions are shown in Fig. 6 where we plot the
specific heat and topological susceptibility as a function of temperature. In
Figs. 6(a) and (c), we plot the specific heat and topological susceptibility
for $D_{\parallel}=0.2$, $K=0$ and for 4 values of the external field,
$h=0,\,1.6,\,5,\,6.6$. In Figs. 6(b) and (d), we plot the specific heat and
topological susceptibility for $D_{\parallel}=0.2$, $K=0.1$ and for 4 values
of the external field, $h=0,\,1.8,\,4,\,5$. The locations of phase transitions
are identified by tracking the peaks in the specific heat and topological
susceptibility. It is necessary to account for the finite size effects and
scaling relations to pinpoint the exact locations of phase transitions. Even
though we do not perform the scaling analysis, by comparing different system
sizes we conclude that the chosen system size is large enough and allows us to
correctly describe the schematic of the phase diagrams for our systems.

Figure 7: A snapshot of SkX in a triangular AFM obtained by simulated
annealing for $K=0$, $D_{\parallel}=0.2$, $h=3.1$, and $T=0.02$ on a $75\times
75$ triangular lattice with periodic boundary conditions. The letters on the
left side (A, B, C) indicate the plotted sublattice. The colors indicate the
orthogonal to the plane component of the spins.

Figure 8: Magnetic field–temperature (h–T) phase diagram of a triangular
lattice AFM obtained by parallel tempering Monte Carlo simulations on a
lattice of $75\times 75$ atoms. The red circles correspond to peaks of the
specific heat and the blue squares correspond to peaks of the topological
susceptibility. (a) Phases SP, SkX, and VL are realized when
$D_{\parallel}=0.2$ and $K=-0.05$. (b) Phases SP, SkX, and VL are realized
when $D_{\parallel}=0.2$ and $K=0$. (c) Phases SP, UUD, SkX, and VL are
realized when $D_{\parallel}=0.2$ and $K=0.05$. (d) Phases SP, Y, UUD, SkX,
and VL are realized when $D_{\parallel}=0.2$ and $K=0.1$. The color indicates
the total topological charge in a region of $75\times 75$ spins.

The results of our calculations are shown in Fig. 8. In the absence of DMI,
our results agree with the known results [55; 53; 54] (not shown). In a small
magnetic field, a 120$\degree$ three-sublattice state is canted, forming a Y
state with one spin aligned by the magnetic field and two others forming an
angle with the opposite direction. At fields $h>3$, we observe a coplanar V
canted state where two spins are aligned with each other and form V with the
third spin. Note that at low temperatures, the degeneracy between the
collinear and coplanar states is lifted by the order-from-disorder effects
[56]. At $T=0$, $h=3$ and in the intermediate region, the system orders into
UUD state with two spins aligned along the magnetic field and the third spin
aligned in the opposite direction (see Fig. 1 in Ref. [53] for representations
of Y, V, and UUD states).

In the presence of DMI, we identify regions in the phase diagram corresponding
to the spiral, antiferromagnetic SkX, and vortical-like (VL) states (see Figs.
8(a), 8(b), 8(c), and 8(d)). Our results mostly agree with Ref. [11] except
for the presence of VL phase. According to our results, antiferromagnetic SkX
phase appears even for relatively small strength of DMI, $D_{\parallel}=0.2$,
as can be seen in Fig. 8. At low temperatures and low magnetic fields the
system prefers the spiral state as a natural extension of the behavior at
$T=0$ in Fig. 2. Increasing the magnetic field or temperature drives system
into three-sublattice SkX phase. In Fig. 5(b), we show the spin structure
factor for SkX phase realizable for $K=0$, $D_{\parallel}=0.2$, $h=3.1$, and
$T=0.02$. A snapshot of spin texture corresponding to SkX phase is shown in
Fig. 7. The SP-SkX transition is identified by peaks in the specific heat and
the topological susceptibility. The SkX-VL transition is identified by peaks
in the topological susceptibility (see Fig. 8). The topological charge can
also be used to identify the transitions at lower temperatures as can be seen
from color representation of the topological charge in Fig. 8. The exact
location of the transition to paramagnetic (P) phase is harder to identify as
peaks for the specific heat and the topological susceptibility do not
coincide, and the average topological charge is nonzero in P phase. This
behavior has been reported recently for realizations of SkX in ultrathin films
Pd/Fe/Ir(111) [52]. Due to the presence of a well defined temperature-driven
phase transition from a spiral phase, accompanied by a sharp rise in
topological charge, we predict that SkX phase will persist even up to very low
magnetic fields, where $h=0.1$ is the lowest field considered in Fig. 8. As
expected, the topological charge vanishes at $h=0$.

We further explore the effect of the easy axis ($K>0$) and easy plane ($K<0$)
magnetic anisotropies. In Fig. 8(a), we plot the phase diagram for $K=-0.05$
which at $h=0$, $T=0$ corresponds to the spiral phase in Fig. 2. We identify
regions corresponding to the SP, SkX, and VL phases, with SkX occupying a
larger portion of the phase diagram compared to Fig. 8(b) corresponding to
$K=0$. In Fig. 8(c), we plot the phase diagram for $K=0.05$. We identify
regions corresponding to the SP, SkX, and VL phases; however, the region of
SkX phase is shrunk due to appearance of additional UUD region. The UUD region
becomes larger as we increase the easy axis anisotropy to $K=0.1$, as can be
seen in Fig. 8(d). In addition, for $K=0.1$ we identify a pocket corresponding
to Y phase. Overall, the presence of easy axis magnetic anisotropy is not
beneficial for realization of the SkX phase. On the other hand, the easy plane
magnetic anisotropy can be beneficial for driving a spiral phase into skyrmion
phase by applying magnetic field.

Of the systems listed in Table 1, $\text{{Fe}/MoS}_{2}$ and
$\text{{Fe}/WSe}_{2}$ seem to be the most promising candidates for realizing
the SkX phase thanks to the presence of strong DMI. However, large magnetic
anisotropy in $\text{{Fe}/WSe}_{2}$ will tend to suppress this phase,
according to the phase diagrams in Fig. 8. By exploiting the high-temperature
SkX regions in Fig. 8, one can apply magnetic fields that are much smaller
compared to the exchange scale given by $J/\mu_{B}$. For $\text{{Fe}/MoS}_{2}$
this field could be on the order of 20 T. In order to obtain skyrmions in a
three-sublattice antiferromagnet, it is desirable to find materials with
weaker exchange interaction and smaller magnetic anisotropy. However, a weak
easy plane magnetic anisotropy can be beneficial for stabilizing the SkX
phase.

##  IV Conclusion

We have studied the phase diagram of triangular AFM in the presence of DMI and
magnetic anisotropy. Our results based on parallel tempering Monte Carlo
simulations suggest feasibility of realizing spiral and skyrmion phases in
triangular AFMs. The spiral phase appears at low temperatures and weak
magnetic fields due to the presence of DMI. By adding magnetic field and by
raising the temperature, one can drive the system into the skyrmion lattice
phase. Our Monte Carlo simulations indicate that both spirals and skyrmion
lattices can appear even for relatively weak DMI, but their stability will be
hindered by magnetic anisotropy. However, a weak easy plane magnetic
anisotropy can be beneficial for stabilizing the skyrmion phase. We observe
that skyrmion lattices can be stabilized below the threshold
$D_{\parallel}=0.2$ suggested in Refs [33; 34].

As a realization of triangular AFMs, we have considered TMD-based
antiferromagnetic triangular lattices. We have used first-principles
computations to determine the parameters of magnetic interactions in monolayer
TMDs ($\text{MoS}_{2}$,$\text{WS}_{2}$, and $\text{WSe}_{2}$) onto which 4
different transition metals (Co, Mn, Cr, and Fe) are adsorbed. We have shown
that there are preferred adsorption sites, according to the strength of the
Hubbard-U potential, and that all the compounds are antiferromagnets with
nonzero uniaxial anisotropy and DMI. According to our phase diagrams, at low
temperatures the ground state of $\text{{Cr}/MoS}_{2}$, $\text{{Fe}/MoS}_{2}$,
and $\text{{Fe}/WSe}_{2}$ is a three-sublattice spiral. By further lowering
anisotropy using materials engineering, TMD-based antiferromagnetic triangular
lattices, such as $\text{{Fe}/WSe}_{2}$, can potentially exhibit skyrmion
lattices. On the other hand, we predict that $\text{{Fe}/MoS}_{2}$ can
potentially host a three-sublattice skyrmion crystal in the presence of
experimentally feasible magnetic fields.

###### Acknowledgements.

AAK and AR were supported by the U.S. Department of Energy, Office of Science,
Basic Energy Sciences, under Award No. DE-SC0021019. WF, PHC, and KDB were
supported by the National Science Foundation (NSF) through the Nebraska
Materials Research Science and Engineering Center (MRSEC) (grant no.
DMR-1420645) and grant No. DMR-1916275. This research used resources of the
National Energy Research Scientific Computing Center (NERSC), a U.S.
Department of Energy Office of Science User Facility operated under Contract
No. DE-AC02-05CH11231. Part of this work was also completed utilizing the
Holland Computing Center of the University of Nebraska, which receives support
from the Nebraska Research Initiative.

## Appendix A: Fourier transformed Hamiltonian for energy minimization

The elements of the Fourier transform of the Hamiltonian in Eq. (2) are given
by the following matrix:

$H_{q}=\begin{pmatrix}A_{1}&\text{m}_{1,2}(\mathbf{k})&\text{m}_{1,3}(\mathbf{k})\\\
\text{m}^{*}_{1,2}(\mathbf{k})&A_{2}&\text{m}_{2,3}(\mathbf{k})\\\
\text{m}^{*}_{1,3}(\mathbf{k})&\text{m}^{*}_{2,3}(\mathbf{k})&A_{3}\\\
\end{pmatrix}.$ (12)

The $A_{\alpha}=\text{diag}(0,0,-\text{K})$ is the single-ion anisotropy for
sublattice atom $\alpha$ and the submatrix $\text{m}_{\alpha,\beta}$ denotes
the interaction between neighboring spins with sublattice indices $\alpha$,
$\beta=\\{1,2,3\\}$. Using $k$ and $l$ as the Cartesian coordinates indices,
these interaction submatrices read:

$m^{kl}_{1,2}(\mathbf{k})=(1+e^{i\mathbf{k}\cdot(\mathbf{e}_{1}+\mathbf{e}_{2})}+e^{i\mathbf{k}\cdot(2\mathbf{e}_{1}-\mathbf{e}_{2})})\delta^{kl}+D_{\parallel}\sum_{n}(\mathbf{e}_{1}+\mathbf{e}_{2}e^{i\mathbf{k}\cdot(\mathbf{e}_{1}+\mathbf{e}_{2})}+\mathbf{e}_{3}e^{i\mathbf{k}\cdot(2\mathbf{e}_{1}-\mathbf{e}_{2})})^{n}\varepsilon^{kln};$
(13)
$m^{kl}_{1,3}(\mathbf{k})=(1+e^{i\mathbf{k}\cdot(\mathbf{e}_{1}+\mathbf{e}_{2})}+e^{i\mathbf{k}\cdot(-\mathbf{e}_{1}+2\mathbf{e}_{2})})\delta^{kl}+D_{\parallel}\sum_{n}(\mathbf{e}_{2}+\mathbf{e}_{1}e^{i\mathbf{k}\cdot(\mathbf{e}_{1}+\mathbf{e}_{2})}+\mathbf{e}_{3}e^{i\mathbf{k}\cdot(-\mathbf{e}_{1}+2\mathbf{e}_{2})})^{n}\varepsilon^{kln};$
(14)
$m^{kl}_{2,3}(\mathbf{k})=(1+e^{i\mathbf{k}\cdot(\mathbf{e}_{1}-2\mathbf{e}_{2})}+e^{i\mathbf{k}\cdot(2\mathbf{e}_{1}-\mathbf{e}_{2})})\delta^{kl}+D_{\parallel}\sum_{n}(\mathbf{e}_{3}+\mathbf{e}_{2}e^{i\mathbf{k}\cdot(\mathbf{e}_{1}-2\mathbf{e}_{2})}+\mathbf{e}_{1}e^{i\mathbf{k}\cdot(2\mathbf{e}_{1}-\mathbf{e}_{2})})^{n}\varepsilon^{kln},$
(15)

with $D_{\parallel}$ being the strength of the DMI in reduced units,
$\delta^{kl}$ being the Kronecker delta symbol, $\mathbf{e}_{1}=\\{1,0,0\\}$,
$\mathbf{e}_{2}=\\{\frac{1}{2},\frac{\sqrt{3}}{2},0\\}$ and
$\mathbf{e}_{3}=\mathbf{e}_{1}-\mathbf{e}_{2}$, and $\varepsilon^{kln}$ being
the Levi-Civita tensor.

## Appendix B: Feedback-optimized parallel tempering Monte Carlo

The parallel tempering Monte Carlo algorithm simulates a set of noninteracting
replicas at temperatures ${T_{1},\dots},T_{M}$ where $M$ is the number of
replicas [48; 47; 57; 58]. After performing a fixed number of sweeps
(typically a few), the algorithm suggests a swap of replicas at neighbouring
temperatures, $T_{i}$ and $T_{i+1}$, accepting the swap with probability,

$p(E_{i},T_{i}\rightarrow E_{i+1},T_{i+1})=\min\\{1,\exp(\Delta\beta\Delta
E)\\},$ (16)

where we introduce the difference between inverse temperatures,
$\Delta\beta=1/(T_{i+1})-1/(T_{i})$, and the difference in energy of the two
replicas, $\Delta E=E_{i+1}-E_{i}$. The algorithm leads to diffusion of
replicas in temperature space. Visits of the high temperature regions, where
relaxation happens faster, facilitate relaxation at lower temperatures. Thus,
the system can efficiently relax in the presence of complicated energy
landscapes. In this process, it is important to maximize the number of round
trips between the lowest and the highest temperature, $T_{1}$ and $T_{M}$. The
number of round trips will strongly depend on the choice of temperature set,
${T_{1},\dots},T_{M}$. A temperature set can be chosen according to a
geometric progression. This, however, will result in suppressed exchanges
between temperatures at a phase transition where specific heat diverges.

In this work, we use approach suggested in Ref. [49] where the temperature set
is optimized by measuring the diffusion of replicas. In particular, a replica
assumes index “up” after visiting $T_{1}$ and a replica assumes index “down”
after visiting $T_{M}$. The index is rewritten from “up” to “down” every time
the “up” replica visits $T_{M}$, and the index is rewritten from “down” to
“up” every time the “down” replica visits $T_{1}$. After each sweep, we update
quantities $n_{\text{up}}(T_{i})$ and $n_{\text{down}}(T_{i})$ for each
$T_{i}$ by adding $1$ to $n_{\text{up}}(T_{i})$ if the replica at $T_{i}$ has
index “up” and by adding $1$ to $n_{\text{down}}(T_{i})$ if the replica at
$T_{i}$ has index “down”. The quantity of interest, $f(T)$, evaluated for each
temperature characterizes the fraction of replicas which have visited either
$T_{1}$ or $T_{M}$:

$f(T_{i})=\frac{n_{\text{up}}(T_{i})}{n_{\text{up}}(T_{i})+n_{\text{down}}(T_{i})},$
(17)

where $n_{\text{up/down}}(T_{i})$ counts the total number of “up”/“down”
walkers throughout the whole Monte Carlo simulation for a temperature $T_{i}$.
Using function $f(T_{i})$ recalculated for each temperature set, we
iteratively maximize the current of replicas in temperature space [49]. For an
optimal temperature set, the function $f(T_{i})$ will decrease at a constant
rate as a function of the temperature index, as seen in Fig. 9(a).

Figure 9: Fraction of temperatures $f(T)$ moving from the lowest to the
highest temperature as a function of the temperature index in (a), and as a
function of temperature in (b). The calculation corresponds to parameters:
$D_{\parallel}=0.2$, $K=0.1$, and $h=0.8$. Figure 10: Spin configurations for
extracting parameters in the model Hamiltonian: (a) and (b) in-plane coplanar
120° spin configurations with an opposite chirality. (c) an out-of-plane
coplanar 120° spin configuration. (d) and (e) coplanar 120° spin
configurations in a $3\times 1$ unit cell.

In Fig. 9, we give an example of $f(T_{i})$ calculated before (a geometric
progression) and after the optimization procedure. To increase the temperature
diffusion, temperatures in a temperature set are shifted toward the region of
a phase transition, as can be seen in Fig. 9(b).

## Appendix C: Extracting parameters from model Hamiltonian

To extract parameters of Hamiltonian (1), we perform several constrained DFT
calculations with specified non-collinear spin configurations. To account for
the $120\degree$ in-plane antiferromagnetic ordering, we use a
$\sqrt{3}\times\sqrt{3}$ unit cell (three transition atoms per unit cell) to
perform the calculations. The nearest-neighbor exchange parameter $J$ is
obtained without spin-orbit interaction, which also results in vanishing
anisotropy and DMI. We tilt the spins out of plane by several small angles.
The energy difference with respect to the in-plane spin configuration as a
function of the tilting angle $\theta$ is $\Delta
E\approx\frac{27}{2}J\theta^{2}$ which allows us to determine the exchange
parameter by fitting with the energy as a function of the tilting angle
$\theta$. The single-ion anisotropy $K$, and in-plane and out-of-plane DMI
components can be determined by taking the energy difference between suitably
specified spin configurations. For out-of-plane DMI component, we use two in-
plane coplanar 120° spin configurations with opposite chirality as shown in
Figs. 10(a) and (b). The energy difference between these two spin
configurations is $\Delta E=9\sqrt{3}\tilde{\text{D}}_{\perp}$. For single-ion
anisotropy $K$, we use one in-plane and one out-of-plane coplanar 120° spin
configurations as shown in Figs. 10(a) and (c). The energy difference between
these two spin configurations is $\Delta
E=\frac{9}{2}\sqrt{3}\tilde{\text{D}}_{\perp}+\frac{3}{2}K$. For in-plane DMI
component, we notice that it is not possible to obtain it in a
$\sqrt{3}\times\sqrt{3}$ unit cell as the contribution to the total energy
from in-plane DMI component is canceled out by symmetry regardless of how we
specify the spin configurations. To break this symmetry, a $3\times 1$ unit
cell (three transition atoms per unit cell) is used for the calculation of the
in-plane DMI component. We use two spin configurations with opposite chirality
as shown in Figs. 10(d) and (e). The energy difference between these two spin
configurations is $\Delta E=\frac{9}{2}\sqrt{3}\tilde{\text{D}}_{\parallel}$.

## References

  * Baltz  _et al._ [2018] V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, Rev. Mod. Phys. 90, 015005 (2018).
  * Olejník  _et al._ [2018] K. Olejník, T. Seifert, Z. Kašpar, V. Novák, P. Wadley, R. P. Campion, M. Baumgartner, P. Gambardella, P. Němec, J. Wunderlich, J. Sinova, P. Kužel, M. Müller, T. Kampfrath, and T. Jungwirth, Sci. Adv. 4, eaar3566 (2018).
  * Rößler  _et al._ [2006] U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Nature 442, 797 (2006).
  * Muhlbauer  _et al._ [2009] S. Muhlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Boni, Science 323, 915 (2009).
  * Yu  _et al._ [2010] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Nature 465, 901 (2010).
  * Barker and Tretiakov [2016] J. Barker and O. A. Tretiakov, Phys. Rev. Lett. 116, 147203 (2016).
  * Zhang  _et al._ [2016] X. Zhang, Y. Zhou, and M. Ezawa, Sci. Rep. 6, 24795 (2016).
  * Göbel  _et al._ [2017] B. Göbel, A. Mook, J. Henk, and I. Mertig, Phys. Rev. B 96, 060406 (2017).
  * Zarzuela  _et al._ [2019] R. Zarzuela, S. K. Kim, and Y. Tserkovnyak, Phys. Rev. B 100, 100408 (2019).
  * Li and Kovalev [2020] B. Li and A. A. Kovalev, Phys. Rev. Lett. 125, 257201 (2020).
  * Rosales  _et al._ [2015] H. Rosales, D. Cabra, and P. Pujol, Physical Review B 92, 214439 (2015).
  * Leonov and Mostovoy [2015] A. O. Leonov and M. Mostovoy, Nature Communications 6, 8275 (2015).
  * Díaz  _et al._ [2019] S. A. Díaz, J. Klinovaja, and D. Loss, Phys. Rev. Lett. 122, 187203 (2019).
  * Kamiya and Batista [2014] Y. Kamiya and C. D. Batista, Phys. Rev. X 4, 011023 (2014).
  * Gao  _et al._ [2020] S. Gao, H. D. Rosales, F. A. G. Albarracín, V. Tsurkan, G. Kaur, T. Fennell, P. Steffens, M. Boehm, P. Čermák, A. Schneidewind, E. Ressouche, D. C. Cabra, C. Rüegg, and O. Zaharko, Nature 586, 37 (2020).
  * Lin  _et al._ [2015] S.-Z. Lin, A. Saxena, and C. D. Batista, Phys. Rev. B 91, 224407 (2015).
  * Güngördü  _et al._ [2016] U. Güngördü, R. Nepal, O. A. Tretiakov, K. Belashchenko, and A. A. Kovalev, Phys. Rev. B. 93, 064428 (2016).
  * Kovalev and Sandhoefner [2018] A. A. Kovalev and S. Sandhoefner, Frontiers in Physics 6, 00098 (2018).
  * Yu  _et al._ [2018] X. Z. Yu, W. Koshibae, Y. Tokunaga, K. Shibata, Y. Taguchi, N. Nagaosa, and Y. Tokura, Nature 564, 95 (2018).
  * Moreau-Luchaire  _et al._ [2016] C. Moreau-Luchaire, C. Moutafis, N. Reyren, J. Sampaio, C. A. F. Vaz, N. V. Horne, K. Bouzehouane, K. Garcia, C. Deranlot, P. Warnicke, P. Wohlhüter, J.-M. George, M. Weigand, J. Raabe, V. Cros, and A. Fert, Nature Nanotechnology 11, 444 (2016).
  * Woo  _et al._ [2016] S. Woo, K. Litzius, B. Krüger, M.-Y. Im, L. Caretta, K. Richter, M. Mann, A. Krone, R. M. Reeve, M. Weigand, P. Agrawal, I. Lemesh, M.-A. Mawass, P. Fischer, M. Kläui, and G. S. D. Beach, Nature Materials 15, 501 (2016).
  * Boulle  _et al._ [2016] O. Boulle, J. Vogel, H. Yang, S. Pizzini, D. de Souza Chaves, A. Locatelli, T. O. Menteş, A. Sala, L. D. Buda-Prejbeanu, O. Klein, M. Belmeguenai, Y. Roussigné, A. Stashkevich, S. M. Chérif, L. Aballe, M. Foerster, M. Chshiev, S. Auffret, I. M. Miron, and G. Gaudin, Nature Nanotechnology 11, 449 (2016).
  * Wang  _et al._ [2015] W. Wang, A. Narayan, L. Tang, K. Dolui, Y. Liu, X. Yuan, Y. Jin, Y. Wu, I. Rungger, S. Sanvito, and F. Xiu, Nano letters 15, 5261 (2015).
  * Wu  _et al._ [2015] H.-C. Wu, C. Ó. Coileáin, M. Abid, O. Mauit, A. Syrlybekov, A. Khalid, H. Xu, R. Gatensby, J. J. Wang, H. Liu, L. Yang, G. S. Duesberg, H.-Z. Zhang, M. Abid, and I. V. Shvets, Scientific Reports 5, 15984 (2015).
  * Allain  _et al._ [2015] A. Allain, J. Kang, K. Banerjee, and A. Kis, Nature materials 14, 1195 (2015).
  * Safeer  _et al._ [2019] C. K. Safeer, J. Ingla-Aynés, F. Herling, J. H. Garcia, M. Vila, N. Ontoso, M. R. Calvo, S. Roche, L. E. Hueso, and F. Casanova, Nano Letters 19, 1074 (2019).
  * Caciuc  _et al._ [2018] V. Caciuc, N. Atodiresei, and S. Blügel, Physical Review Materials 2, 084001 (2018).
  * Zhang  _et al._ [2019] W. Zhang, P. K. J. Wong, X. Zhou, A. Rath, Z. Huang, H. Wang, S. A. Morton, J. Yuan, L. Zhang, R. Chua, S. Zeng, E. Liu, F. Xu, Ariando, D. H. C. Chua, Y. P. Feng, G. van der Laan, S. J. Pennycook, Y. Zhai, and A. T. S. Wee, ACS Nano 13, 2253 (2019).
  * Žutić  _et al._ [2019] I. Žutić, A. Matos-Abiague, B. Scharf, H. Dery, and K. Belashchenko, Materials Today 22, 85 (2019).
  * Latzke  _et al._ [2015] D. W. Latzke, W. Zhang, A. Suslu, T.-R. Chang, H. Lin, H.-T. Jeng, S. Tongay, J. Wu, A. Bansil, and A. Lanzara, Physical Review B 91, 235202 (2015).
  * Dzyaloshinsky [1958] I. Dzyaloshinsky, J. Phys. Chem. Solids 4, 241 (1958).
  * Moriya [1960a] T. Moriya, Phys. Rev. 120, 91 (1960a).
  * Osorio  _et al._ [2017] S. A. Osorio, H. D. Rosales, M. B. Sturla, and D. C. Cabra, Phys. Rev. B 96, 024404 (2017).
  * Mohylna and Žukovič [2020] M. Mohylna and M. Žukovič, Acta Physica Polonica A 137, 616 (2020).
  * Mohylna  _et al._ [2021] M. Mohylna, J. Buša. Jr., and M. Žukovič, Journal of Magnetism and Magnetic Materials 527, 167755 (2021).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Ozaki [2003] T. Ozaki, Phys. Rev. B 67, 155108 (2003).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Blöchl [1994] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Liechtenstein  _et al._ [1995] A. I. Liechtenstein, V. I. Anisimov, and J. Zaanen, Phys. Rev. B 52, R5467 (1995).
  * Han  _et al._ [2006] M. J. Han, T. Ozaki, and J. Yu, Phys. Rev. B 73, 045110 (2006).
  * Ozaki and Kino [2004] T. Ozaki and H. Kino, Phys. Rev. B 69, 195113 (2004).
  * Dzyaloshinskii [1958] I. E. Dzyaloshinskii, J. Phys. Chem. Solids 4, 241 (1958).
  * Moriya [1960b] T. Moriya, Phys. Rev. 120, 91 (1960b).
  * Luttinger and Tisza [1946] J. Luttinger and L. Tisza, Physical Review 70, 954 (1946).
  * Swendsen and Wang [1986] R. H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 57, 2607 (1986).
  * Hukushima and Nemoto [1996] K. Hukushima and K. Nemoto, Journal of the Physical Society of Japan 65, 1604 (1996).
  * Katzgraber  _et al._ [2006] H. G. Katzgraber, S. Trebst, D. A. Huse, and M. Troyer, Journal of Statistical Mechanics: Theory and Experiment 2006, P03018 (2006).
  * Berg and Lüscher [1981] B. Berg and M. Lüscher, Nuclear Physics B 190, 412 (1981).
  * Lau and Dasgupta [1989] M.-h. Lau and C. Dasgupta, Phys. Rev. B 39, 7212 (1989).
  * Böttcher  _et al._ [2018] M. Böttcher, S. Heinze, S. Egorov, J. Sinova, and B. Dupé, New Journal of Physics 20, 103014 (2018).
  * Gvozdikova  _et al._ [2011] M. Gvozdikova, P. Melchy, and M. Zhitomirsky, Journal of Physics: Condensed Matter 23, 164209 (2011).
  * Seabra  _et al._ [2011] L. Seabra, T. Momoi, P. Sindzingre, and N. Shannon, Physical Review B 84, 214418 (2011).
  * Kawamura and Miyashita [1985] H. Kawamura and S. Miyashita, Journal of the Physical Society of Japan 54, 4530 (1985).
  * Villain  _et al._ [1980] J. Villain, R. Bidaux, J.-P. Carton, and R. Conte, Journal de Physique 41, 1263 (1980).
  * Marinari and Parisi [1992] E. Marinari and G. Parisi, Europhysics Letters (EPL) 19, 451 (1992).
  * Lyubartsev  _et al._ [1992] A. P. Lyubartsev, A. A. Martsinovski, S. V. Shevkunov, and P. N. Vorontsov-Velyaminov, The Journal of Chemical Physics 96, 1776 (1992).


when at least three HLC hits are recorded in the DeepCore fiducial volume, the
$\rho$ and $z$ positions of the first HLC hit are likely to be near the
interaction point.

The second method is VertexGuess, which is the position of the first hit DOM
in a cleaned hit series. For a neutrino signal event, the interaction happens
in the DeepCore fiducial volume, whereas an atmospheric muon is expected to
leave early hits in the veto region. Therefore, the variable VertexGuessZ
provides a quick guess at an early stage in the event selection process for
the $z$ position of the interaction point.

###  A.2 Charge Variables

The event selection of Analyses $\mathcal{A}$ and $\mathcal{B}$ use the
following charge information from a cleaned hit series to identify downward-
going atmospheric muon events.

First, NAbove200 is the integrated charge from hit DOMs above $z=-200$ m and
$2~{}\mu$s before the DeepCore trigger. Compared to an upward-going neutrino,
a downward-going atmospheric muon is more likely to hit DOMs in the upper part
of the detector and deposit a larger total charge in the veto region before
trigger time.

The second charge variable is the charge ratio (QR) and has several distinct
implementations. The variable QR6 is the fraction of accumulated charges from
cleaned hits during the first 600 ns after the DeepCore trigger with respect
to the total accumulated charge from all cleaned hits; that is,

$QR6=\frac{1}{Q_{\text{total}}}\sum_{i}Q_{i}$ (9)

for $0\,\rm{ns}<t(Q_{i})<600\,\rm{ns}$. A contained neutrino event deposits
more charge in a shorter time scale than a through-going muon, so a QR6 value
closer to zero indicates a potential atmospheric muon event. Similarly, the
variable QR3 is calculated with a tighter time window of 300 ns instead of 600
ns. Further, to reduce the impact from noise hits, the charge ratio variables
C2QR6 and C2QR3 are calculated as described above for QR6 and QR3,
respectively, but with the two initial hits in the cleaned hit series removed.

###  A.3 Veto Regions

Veto variables are used in both analyses to identify and remove atmospheric
muon backgrounds. These variables define veto regions, event-by-event, based
on the positions and/or photon arrival times of the hit DOMs.

First used in Ha (2011) and optimized for DeepCore atmospheric oscillation
searches in Euler (2014), the Veto Identified Causal Hits (VICH) algorithm
uses the cleaned hit closest to the trigger time as a reference hit and
determines a veto region in which the hits may be causally connected with the
reference hit. Fig. 18 shows the causal veto region (in red), which is defined
by

  * •

$\Delta r/c<2.5\ \rm{\mu s}$,

  * •

$\Delta r/c<-\frac{2}{3}\Delta t+\frac{1}{3}\ \rm{\mu s}$,

  * •

$\Delta t-0.15\ \rm{\mu s}<\Delta r/c<\Delta t+1.85\ \rm{\mu s}$,

where $\Delta r$ and $\Delta t$ are the distance and photon arrival time
between a given hit and the reference hit, respectively, and $c$ is the speed
of light in vacuum. The width of the veto region accounts for a reasonable
amount of scattering in the ice combined with the typical time scale of a GeV-
scale neutrino event in IceCube (1 $\mu$s). Restrictions of the red region
close to the time and position of the reference hit are made to allow
variations between the interaction time, the first trigger, and subsequent
hits that could occur in a neutrino event. With the defined veto region, the
VICH variable is the integrated charge from cleaned hits that lie in the veto
region. The more charge from cleaned hits deposited inside the veto region,
the more likely that the event is caused by an atmospheric muon.

Figure 18: Definition of veto region based on the VICH algorithm. The four
lines surrounding the red region define the causal veto volume, event-by-
event, based on the hit closest to the trigger time. The more total charge
from cleaned hits deposited inside the veto region, the more likely that the
event is caused by an atmospheric muon.

Another effective way to identify atmospheric muon background is to count the
number of corridor DOMs. As shown in Fig. 1, the IceCube and DeepCore strings
are arranged in a roughly triangular lattice in the horizontal plane. An
atmospheric muon interacting in DeepCore can potentially come from a corridor
(shown as the purple arrow in the figure) leaving no detected light in the
veto region. To identify these muons, known corridor regions are studied.
Given an interaction vertex, the corridor algorithm first finds direct hits
from a cleaned hit series. Direct hits are hits due to minimally-scattered
photons in the ice between their emission point and their detection in the DOM
(the procedure to identify direct hits can be found in Aartsen et al. (2015)).
The algorithm then looks through the closest IceCube strings along the known
corridors and counts the number of direct hits on those strings, which is
defined as the number of corridor DOMs.

###  A.4 Center of Gravity

The Center of Gravity (CoG) is a parameter that measures space-time
correlations between assumed signal hits in DeepCore and likely veto hits in
the surrounding IceCube DOMs. For a total of $N$ hits, the average position
$\vec{x}_{\text{CoG}}$ is given by

$\vec{x}_{\text{CoG}}=\frac{1}{N}\sum_{i}^{N}\vec{x}_{i},$ (10)

where $\vec{x}_{i}$ is the position of the $i^{\text{th}}$ hit DOM relative to
the center of IceCube. Then, an average CoG time $t_{\text{CoG}}$ is
calculated by assuming a simple cascade hypothesis with light propagating from
the CoG without scattering;

$t_{\text{CoG}}=\frac{1}{N}\sum_{i}^{N}\bigg{(}t_{i}-\frac{|\vec{x}_{i}-\vec{x}_{\text{CoG}}|}{c_{\text{ice}}}\bigg{)}.$
(11)

Here, $t_{i}$ is the photon arrival time at the $i^{\text{th}}$ hit DOM, and
$c_{\text{ice}}$ is the speed of light in ice. The same calculations can be
applied to obtain the average position and time from a specific group of
cleaned hits.

CoG is often used to check for causality. The CoG position and time are
calculated from the cleaned hits inside the DeepCore fiducial volume. Then,
for each cleaned hit in the veto region, a causally connected veto hit is
identified if its $\vec{x}$ and $t$ satisfy

$0.25\text{
m/ns}\leq\frac{|\vec{x}_{\text{CoG}}-\vec{x}|}{t_{\text{CoG}}-t}\leq 0.4\text{
m/ns}.$ (12)

This requirement is used to reject muon tracks entering the fiducial volume
after leaving hits in the veto region.

CoG is also used to evaluate event topology by dividing the time-sorted and
cleaned hit series equally into four quartiles. Quartile 1 (Q1) consists of
the earliest hits; quartile 4 (Q4) the latest hits. The CoG position and time
are calculated for the hits in each quartile. The following three variables,
based on these CoG quartiles, are used to identify background events.

First, to identify atmospheric muons, a z-travel variable is defined as the
vertical distance between the CoG $z$ position of Q1 and that from all cleaned
hits in the fiducial volume; that is, z-travel $\equiv
z_{\text{all}}-z_{\text{Q1}}$. For a downward-going atmospheric muon event,
its earlier hits tend to have a $z_{\text{Q1}}$ position above the average
$z_{\text{all}}$, resulting in a negative z-travel. Similarly, an upward-going
neutrino event usually has a positive z-travel.

Second, the spatial separation between Q1 and Q4 is also used for rejecting
muon background events. The Q1-Q4 separation is defined as
$|\vec{x}_{\text{Q4}}-\vec{x}_{\text{Q1}}|$. Because atmospheric muons usually
travel long distances across the detector, they often have a larger spatial
separation between their earlier and later hits compared to neutrino events.

Third, to discriminate noise triggers from physics triggers, the space-time
interval $\Delta s^{2}$ between Q1 and Q4 is also used. By definition, $\Delta
s^{2}\equiv|\vec{x}_{\text{Q4}}-\vec{x}_{\text{Q1}}|^{2}-(ct_{\text{Q4}}-ct_{\text{Q1}})^{2}$.
For an event caused by random detector noise, its $\Delta s^{2}$ will appear
either too time-like or too space-like compared to an event due to a neutrino
interaction.

In addition to hit causality and event topology, a charge-weighted CoG can
also estimate the size of an event by determining the standard deviations of
$z$ position ($\sigma_{z}$) and photon arrival time ($\sigma_{t}$) from all
cleaned hits. An atmospheric muon tends to produce hits across the detector
for a longer period of time, so its $\sigma_{z}$ and $\sigma_{t}$ are
typically larger than for a neutrino event.

###  A.5 Quick Track Reconstructions

At lower selection levels with high event rates, computationally-inexpensive
reconstruction algorithms are often used to provide a rough estimate of the
event parameters related to the interacting particle. These parameters include
the particle’s speed, direction, and point of interaction. The following two
quick algorithms assume a track event hypothesis and are used in both Analyses
$\mathcal{A}$ and $\mathcal{B}$.

The improved LineFit, or iLineFit, is based on the LineFit reconstruction.
Assuming an infinitely-long muon track, the simple LineFit algorithm
analytically minimizes a least-squares fit of cleaned hits with respect to the
event parameters. The iLineFit then takes into account effects such as random
detector noise and the scattering and absorption properties of the ice. The
fitting procedure is described in Aartsen et al. (2014c).

The second reconstruction is a likelihood-based single photoelectron fit with
eleven seeds (SPEFit11). Given a track-like event with a set of event
parameters, the likelihood between the expected and observed photon arrival
times is determined for each DOM. The total likelihood from all DOMs is
minimized with respect to the event parameters. This fit runs iteratively from
eleven different starting orientations to avoid falling into local minima.
More information is found in Ahrens et al. (2004).

##  Appendix B Common Event Selection Criteria

This section discusses the basic event filtering at the early stages of
selection processes. These early selection cuts are mostly identical between
Analyses $\mathcal{A}$ and $\mathcal{B}$. The first requirements (Levels 1 and
2) rely on trigger conditions, whereas the next selection criteria (Level 3)
depends largely on veto algorithms.

###  B.1 Common Level 1 and 2

Levels 1 and 2 include the standard online triggering and filtering, both of
which rely on charges recorded by the PMTs in all DOMs. In IceCube, the charge
that a DOM records is measured in effective photoelectron units (p.e.), and
the calibration and characteristics of a PMT are described in Abbasi et al.
(2010b). When a DOM’s PMT exceeds a 0.25 p.e. threshold, an incident hit is
detected, and the DOM is known as a “hit DOM.”

When several nearby DOMs are hit, a local coincidence (LC) occurs, which
indicates a potential neutrino signal event. In particular, two LC types are
of interest: LC1 coincidences between two nearest-neighbor DOMs on a string
and LC2 coincidences between two next-nearest-neighbor DOMs on a string. If
either the LC1 or LC2 condition is met, the LC is called hard local
coincidence (HLC), and the initial hit is a HLC hit. For each recorded HLC, a
full digitization readout is performed Aartsen et al. (2017). If a DOM is hit
with no coincidence from its neighbors, the hit only results in a charge and
time stamp readout for the DOM instead of a full waveform readout.

Both the online trigger (Level 1) and online filter (Level 2) are performed at
the South Pole Aartsen et al. (2017). Level 1 is a simple multiplicity trigger
(SMT) that requires at least three HLC hits within 2.5 $\mu$s among the
DeepCore DOMs; this trigger condition is known as SMT3. The Level 2 online
filter looks for causally-connected hits in the veto and DeepCore regions
using the CoG algorithm discussed in Appendix A.4. For a given event, if the
CoG algorithm identifies one or more causally connected veto hits, the event
is likely caused by an atmospheric muon and is thereby rejected.

###  B.2 Common Level 3

In general, the goal of the Level 3 selection is to remove events that are
triggered by random detector noise and atmospheric muons. Most algorithms at
Level 3 rely on hit information inside and outside the extended DeepCore
volume. This region is defined to include DOMs that are 2100 m below the
surface on all strings except the strings in the outermost three layers of
IceCube. Compared to the standard DeepCore fiducial volume defined in Fig. 1,
the extended DeepCore volume contains five more IceCube strings.

Two algorithms are used to identify events triggered by random detector noise.
First, the NoiseEngine algorithm looks for directionality among hits. It
starts by removing isolated hits and determines a time window that maximizes
the number of cleaned hits. For each cleaned hit within the time window, it is
connected to all other hits if the hit pair satisfies a space-time correlation
window. A map of all possible hit pairs is produced, and they are projected
onto a binned HEALPix sphere. If more than three pairs land in a single
HEALPix bin, then hits are directional, and the event is unlikely to be noise-
triggered. The second algorithm uses charge and hit information from cleaned
hits within a dynamic time window. Noise-triggered events tend to have no more
than two cleaned hits with a total charge less than 2 p.e.

Two selection criteria are applied to quickly identify candidate neutrino
signal events in DeepCore. First, the total amount of charge recorded from a
set of cleaned hits on DOMs inside the extended DeepCore fiducial volume must
be greater than zero. The second selection variable is the VertexGuessZ
described in Appendix A.1. An event passes if its first cleaned hit has a $z$
position below $-120$ m, which is the top of the extended DeepCore fiducial
volume.

To identify obvious atmospheric muon events, two simple charge variables are
defined to study light deposited in the veto region. As discussed in Appendix
A.2, an atmospheric muon tends to have a higher value of NAbove200 than a
neutrino signal event. Thus, an event is rejected if its NAbove200 is $\geq
12$ p.e. A second method uses the CoG algorithm discussed in Appendix A.4. The
algorithm is used to identify causally-connected hits in the veto region, and
their charges are summed. If the total charge is $\geq 7$ p.e., the event is
likely caused by an atmospheric muon and rejected.

The remaining two common Level 3 selection criteria are also charge-related.
The first variable is the ratio of total charge outside the extended DeepCore
fiducial region to that inside. The calculation is performed on a cleaned hit
series, and events with a ratio smaller than 1.5 are kept. Second, the charge
ratio C2QR6 discussed in Appendix A.2 is used. A typical background muon event
has a lower value of C2QR6, so only events with C2QR6 greater than 0.4 are
kept.

With these quick and simple Level 3 criteria applied, the event rate is
diminished by a factor of 20 to roughly 1 Hz. Both Analyses $\mathcal{A}$ and
$\mathcal{B}$ share the same Level 3 selection criteria above, and an
additional cut at Level 3 is applied in Analysis $\mathcal{A}$ as discussed in
Appendix C.1.

Type | Filtering | Analysis $\mathcal{A}$ | Analysis $\mathcal{B}$
---|---|---|---
| Total | DeepCore | L3 | L4 | L5 | L6 | L7 | L7 LE | L3 | L4 | L5 | L6 | L6 LE
atm. $\mu$ | 991000 | 9180 | 970 | 50.5 | 4.10 | 0.443 | 0.100 | 0.092 | 1310 | 44.7 | 0.163 | 0.0297 | 0.0259
Noise | 35900 | 8120 | 284 | 12.0 | 1.80 | 0.102 |  $<$0.001 |  $<$0.001 | 292 | 0.0006 |  $<$0.0006 | 0.0003 |  $<$0.0003
$\nu_{e}$ | 1.84 | 1.72 | 1.26 | 0.783 | 0.544 | 0.362 | 0.325 | 0.194 | 1.29 | 0.278 | 0.180 | 0.149 | 0.126
$\nu_{\mu}$ | 11.3 | 6.36 | 4.76 | 2.50 | 1.63 | 1.01 | 0.676 | 0.552 | 4.93 | 1.02 | 0.558 | 0.396 | 0.342
$\nu_{\tau}$ | 0.293 | 0.270 | 0.206 | 0.134 | 0.103 | 0.074 | 0.051 | 0.045 | 0.210 | 0.052 | 0.038 | 0.031 | 0.028
MC Total | 1030000 | 17300 | 1260 | 65.9 | 8.18 | 1.99 | 1.15 | 0.884 | 1608 | 46.1 | 0.94 | 0.61 | 0.52
Data | 1150000 | 19100 | 1090 | 68.6 | 7.42 | 1.84 | 0.87 | 0.715 | 1981 | 34.9 | 0.844 | 0.504 | 0.432
Table 3:  The event rate in mHz for the common filtering and the subsequent
event selection levels for Analyses $\mathcal{A}$ and $\mathcal{B}$,
respectively. After the final selection level, the analyses only include
events with energies in the region from 5.6 GeV to 56 GeV which is denoted as
“LE”.

##  Appendix C Higher Level Selection for Sample $\mathcal{A}$

This section discusses the progressive stages of the event selection for
Analysis $\mathcal{A}$. Optimized for atmospheric $\nu_{\tau}$ analyses, the
event selection focuses on both cascade- and track-like neutrino signatures at
$\mathcal{O}$(10) GeV, while removing as many background events as possible.
Given three years of detector exposure, Analysis $\mathcal{A}$ expects more
than 55,000 neutrino events.

###  C.1 Sample $\mathcal{A}$, Level 3

In addition to the common selection criteria described in Sec. B.2, an extra
criterion is applied at Level 3 based on charges inside the extended DeepCore
fiducial volume and that outside the volume. The total fiducial charge is
simply the sum of all charges inside the volume. For charges outside the
extended DeepCore fiducial region, an algorithm is performed to search for the
largest clusters of hits, and the total veto charge is the sum of all hits in
the cluster. Based on the charges deposited in the two regions, a cut is
applied to remove potential background muon events.

###  C.2 Sample $\mathcal{A}$, Level 4

After removing events that are likely caused by random detector noise and
atmospheric muons, the Analysis $\mathcal{A}$ event selection uses a boosted
decision tree (BDT) Freund and Schapire (1997) to further reduce atmospheric
muon background at Level 4. Six variables are used in training the BDT,
including QR6, C2QR6, and NAbove200 discussed in Appendix A.2, as well as
VertexGuessZ in Appendix A.1.

Figure 19: (top) Distribution of ToIEVal at Level 4 (before cut applied). Each
shaded color represents the stacked histogram from each event type. Black dots
represent the data distribution. MC events are weighted by world-averaged best
fit oscillation parameters. (bottom) Ratio of distribution from data to that
from MC. Black error bars are the statistical fluctuations from data, whereas
shaded red areas are the uncertainties from limited MC statistics.

The remaining two variables are the estimated speed from a quick iLinefit
reconstruction (see Appendix A.5) and the rough event topology. The concept of
moment of inertia in classical mechanics is adapted to describe the overall
shape of an event with $N$ hits, each of which has a charge of $q_{i}$. The
diagonal elements of the moment of inertia are then given by

$\begin{split}I_{x}=\sum_{i}^{N}q_{i}\big{(}y_{i}^{2}+z_{i}^{2}\big{)},\\\
I_{y}=\sum_{i}^{N}q_{i}\big{(}x_{i}^{2}+z_{i}^{2}\big{)},\\\
I_{z}=\sum_{i}^{N}q_{i}\big{(}x_{i}^{2}+y_{i}^{2}\big{)},\end{split}$ (13)

where $x_{i},y_{i},z_{i}$ are the distance of the $i^{\text{th}}$ hit from
String 36 at $z$ = 0 m. For a spherically-shaped cascade-like event, the
numerical values of $I_{x,y,z}$ are similar, while an elongated track-like
muon background is more likely to have a smaller lateral distribution of hits
than longitudinal. Therefore, a variable based on the tensor of inertia
eigenvalue ratio, (ToIEVal), is defined as:

$\text{ToIEVal}=\frac{I_{\text{smallest}}}{I_{x}+I_{y}+I_{z}}.$ (14)

ToIEVal is included in training the Level 4 BDT, and its distribution is shown
in Fig. 19.

The BDT score distribution is shown in Fig. 2, and a score cut is applied to
keep events with a score above 0.04. Compared to Level 3, the muon background
from MC estimates after the BDT score cut is reduced from 970 mHz to 51 mHz,
while $\approx 55$% of all neutrino events are kept.

###  C.3 Sample $\mathcal{A}$, Level 5

A second BDT at Level 5 is trained to further reduce background muon
contamination. After Level 4 cut is applied, over fifteen times more
atmospheric muons and three times more pure-noise events remain compared to
the integrated number of events from all neutrino flavors. Six variables are
used for training the BDT, four of which are the radial position $\rho$ of the
FirstHLC (Appendix A.1), the Q1-Q4 separation and z-travel from the CoG
algorithm (Appendix A.4), and the VICH veto variable (Appendix A.3). The
distribution of VICH is shown in Fig. 20 as an example.

Figure 20: (top) Distribution of VICH at Level 5 (before cut applied). Each
shaded color represents the stacked histogram from each event type. Black dots
represent the data distribution. MC events are weighted by world averaged best
fit oscillation parameters. (bottom) Ratio of distribution from data to that
from MC. Black error bars are the statistical fluctuation from data, whereas
shaded red areas are the uncertainties from limited MC statistics.

The remaining two variables are the accumulated time and a simple zenith angle
estimation. Because atmospheric muons are likely to travel across the
detector, more time is needed, compared to neutrino events, for the hit DOMs
to detect photons from the light source. Thus, the time to accumulate 75% of
the charge from a cleaned hit series is included in the Level 5 BDT training.
Further, because a majority of muon background events are down going, a fast
reconstruction SPEFit11 (see Appendix A.5) is performed on a cleaned hit
series to provide a rough estimate on the zenith angle of the interacting
particle.

The Level 5 BDT score distribution is shown in Fig. 3, and a BDT score cut is
applied to keep events with a score above 0.04. Compared to Level 4, the
background muon rate after the Level 5 BDT score cut is reduced to $\approx 4$
mHz, which is slightly more than a factor of 10. Moreover, $\approx 85$% of
noise-triggered events are rejected.

###  C.4 Sample $\mathcal{A}$, Level 6

Previous selection criteria have reduced the background rate to $\approx 6$
mHz, which is comparable to the total neutrino rate of $>$ 2 mHz. At this
stage, most obvious background muons are rejected, and the remaining
atmospheric muons are more difficult to be identified. Therefore, Level 6
includes two straight cuts to reject events caused by random detector noise
and two more cuts to identify sneaky atmospheric muons.

Figure 21: (top) FillRatio distribution at Level 6 (before cut applied). Each
shaded color represents the stacked histogram from each event type. Black dots
represent the data distribution. The vertical line is the cut value of 0.05,
events below which are rejected. MC events are weighted by world averaged best
fit oscillation parameters. (bottom) Ratio of distribution from data to that
from MC. Black error bars are the statistical fluctuation from data, whereas
shaded red areas are the uncertainties from limited MC statistics.

To identify events caused by random detector noise, an algorithm called Fill-
Ratio is performed to look for the topology of hits in an event around an
estimated vertex position given by FirstHLC (see Appendix A.1). A sphere is
defined around the vertex with a radius 60% larger than the average distance
from all cleaned hits. The pattern of hits in an event can be estimated with a
Fill-Ratio variable defined by the ratio of number of hit DOMs inside the
sphere to the total number of DOMs inside the sphere. Since hits in a pure-
noise event are randomly scattered across the detector, such event tends to
have a longer mean distance away from its vertex and less hit DOMs inside the
sphere. Therefore, pure-noise events often have smaller values of Fill-Ratio
compared to physics events, as shown in Fig. 21. A straight cut is placed at
0.05, events above which are kept.

Figure 22: (top) NChannel distribution at Level 6 (before cut applied). Each
shaded color represents the stacked histogram from each event type. Black dots
represents the data distribution. The vertical line is the cut value of 8,
events below which are rejected. MC events are weighted by world averaged best
fit oscillation parameters. (bottom) Ratio of distribution from data to that
from MC. Black error bars are the statistical fluctuation from data, whereas
shaded red areas are the uncertainties from limited MC statistics.

A second cut placed on NChannel also help further remove noise-triggered
events. NChannel is the number of hit DOMs in a cleaned hit series, and its
distribution is shown in Fig. 22. Most noise-triggered events have less than
eight hit DOMs. Further, the computationally-intensive, likelihood-based
reconstruction algorithm described in Sec. III.3 fits eight parameters. In
order to have enough degrees of freedom for a reasonable fit, a direct cut is
applied on NChannel to remove events with less than eight hit DOMs, reducing
the contribution from pure-noise events by a factor of 20.

Figure 23: Fractional distribution at Level 6 (before cut applied) in vertex
radial $\rho$ and depth $z$ positions by FiniteReco for atmospheric muon
events. Color axis represents the fraction of atmospheric muon with respect to
the total expected MC events. Red lines represent the cut, within which events
are kept.

For atmospheric muons, two straight cuts are applied. The first cut is applied
on the number of corridor DOMs as explained in Appendix A.3. Events with two
or more direct hits found along those corridors are rejected. The second cut
involves a likelihood-based algorithm known as FiniteReco. Given an infinite
track event hypothesis, FiniteReco reconstructs an event based on the
probabilities of the individual DOMs to see a hit or not Euler (2014). It
provides a relatively quick estimate on the interaction vertex of an event.
Figure 23 shows the fractional 2D distribution of radial $\rho$ and depth $z$
positions from FiniteReco for atmospheric muons. Most atmospheric muon events
have vertices located 125 m away from the center of the detector and above
$-200$ m in depth. Therefore, three cuts are applied such that events are
required to have a vertex position with $z\ <-225$ m, $\rho\ <$ 125 m, and $z\
<\ -3\cdot\rho$. The latter cut removes events at the upper edge of the
DeepCore fiducial volume, where atmospheric muons can enter DeepCore through
the dust layer.

After Level 6 cuts are applied, more than 92% and 96% of atmospheric muon and
noise-triggered events are rejected compared to Level 5. Neutrinos contribute
more than 66% of the sample, which opens the opportunity to use a
computational-expensive but more comprehensive reconstruction algorithm
discussed in Sec. III.3.

###  C.5 Sample $\mathcal{A}$, Level 7

The comprehensive reconstruction discussed in Sec. III.3 is performed between
Level 6 and 7. To reduce the remaining background events and to improve
agreement between data and MC, final selection cuts are applied at Level 7.

Figure 24: Fractional distribution at Level 7 (before cut applied) in vertex
radial $\rho$ and depth $z$ positions by the final reconstruction discussed in
Sec. III.3 for atmospheric muon events. Color axis represents the fraction of
atmospheric muon with respect to the total expected MC events per bin. Red
lines represent the cut, within which events are kept.

The first selection criteria is similar to the containment requirement from
the FiniteReco vertex positions at Level 6. With a more sophisticated
reconstruction method described in Sec. III.3, the fitted interaction vertex
of an event is better than the estimates from FiniteReco. Figure 24 shows the
fractional 2D distributions of radial $\rho$ and depth $z$ positions obtained
from the final reconstruction. An extra containment condition is added to
exclude events below $z$ position of $-500$ m, and events outside the red
lines are rejected.

Two final cuts are applied during the development of the $\mathcal{A}$
selection. The first cut is based on a 2D distribution of reconstructed energy
per number of hit DOMs and the root mean square (RMS) of photon arrival times
from a cleaned hit series. The second cut is related to ‘flaring’ DOMs, which
emit light sporadically. The extra light is not simulated in MC, and a
disagreement between data and MC is shown from the distribution of normalized
RMS of total charges (see Fig. 25). Therefore, a cut is placed to remove
events in which the normalized RMS of total charges are above 85%. The above
two cuts only remove 5% of the total events and do not alter any physics
results. Nonetheless, they are applied.

Figure 25: (top) distribution of normalized RMS of total charge at Level 7
(before cut applied). Each shaded color represents the stacked histogram from
each event type. Black dots represent the data distribution. The vertical line
is the cut value of 0.85, events above which are rejected. MC events are
weighted by world averaged best fit oscillation parameters. (bottom) Ratio of
distribution from data to that from MC. Black error bars are the statistical
fluctuation from data, whereas shaded red areas are the uncertainties from
limited MC statistics.

Finally, only events within the analysis histogram ranges stated in Sec. IV
are used for the oscillation analyses. This restricts to events from all sky
with a reconstructed energy between 5.6 GeV to 56 GeV and a reconstructed
track length less than 1,000 m.

In summary, the event rates as a function of event type and selection level is
shown in Table 3. The neutrino rates are the combination of the NC+CC channels
and use the atmospheric neutrino flux predictions from Honda et al. (2015)
with values of $\theta_{23}$ and $\Delta m^{2}$ from Fogli et al. (2012). At
the final level, 92% of the $\mathcal{A}$ sample is neutrino events, while the
contamination from atmospheric muon and noise-triggered events are 8% and 0.1%
respectively.

##  Appendix D Higher Level Selection for Sample $\mathcal{B}$

This section focuses on the event selection method for analysis $\mathcal{B}$.
It features a set of straight cuts and a boosted decision tree to improve the
purity of neutrino events at the final selection level. Based on simulations,
about 40,000 neutrinos are expected given three years of detector exposure.
The resultant sample is also used for the most recent published measurements
of atmospheric neutrino oscillation parameters from the IceCube Collaboration
Aartsen et al. (2018a).

###  D.1 Sample $\mathcal{B}$, Level 4

After the common Level 3 filtering discussed in Appendix B.2, a set of
straight cuts is applied at Level 4 to further remove events due to random
detector noise and atmospheric muons. These selection requirements rely on hit
information from a cleaned hit series, charge information in veto regions, and
an estimated interaction vertex position of an event.

To ensure that enough information are detected in an event, the first two
selection variables are NChannel and RTFiducialQ. NChannel is the number of
hit DOMs in a cleaned pulseseries, and only events with at least eight cleaned
hits are kept. Further, RTFiducialQ is a charge-related variable from an
algorithm, which searches for clusters of cleaned hits in the DeepCore
fiducial region that satisfy space-time correlations. Events with a minimum of
one cluster with at least 7 p.e. are kept.

Three additional selection criteria are placed based on variables from the
Center of Gravity (CoG) algorithm discussed in Appendix A.4. First, the space-
time interval $\Delta s^{2}$ has some power to distinguish events caused by
random detector noise from physics events. Therefore, a cut is applied to keep
events where $\Delta s^{2}$ is between $-(400$ m)2 and 0 m2. The remaining two
cuts depend on the size of an event, which is estimated by the charge-weighted
spread of vertex $z$ position ($\sigma_{z}$) and photon arrival time
($\sigma_{t}$) from a cleaned hit series. To be specific, only events with
$\sigma_{t}\ \leq$ 1,000 ns and 7 m $\leq\ \sigma_{z}\ \leq$ 100 m are kept.

Because an atmospheric muon event can be identified using the charge
information in the veto regions, two veto requirements are placed. The first
variable counts veto charges using the CoG algorithm (see Appendix A.4). This
cut is very similar to the total veto charge requirement at Level 3 where the
veto region is outside the extended DeepCore volume. At Level 4, the same
algorithm is applied to the veto region outside the standard DeepCore fiducial
volume defined in Fig. 1. With a slightly larger veto volume, a tighter cut is
applied at Level 4 to remove events with a total veto charge greater than 5
p.e. The second veto charge requirement is applied on the VICH variable
discussed in Appendix A.3. With a veto volume defined by the estimated point
of interaction, VICH looks for potentially causally-related hits in the event-
by-event veto region. A cut is placed at 7 p.e. to reject potential
atmospheric muon events.

To increase the purity of $\nu_{\mu}$ CC events, a cut is applied on the
number of direct hits. Direct hits are hits due to photons that experience
minimal scattering in the ice between its emission point and detection in the
DOM. When a muon track passes next to a string, the intersection of its
Cherenkov cone with a string results forms a hyperbolic pattern as a function
of the photon direct arrival depth and times. The orientation of this pattern
is determined by the angle between the string and the passing muon track. The
procedure to identify direct hits is explained in Aartsen et al. (2015), and
only events with at least three direct hits are kept.

Finally, four containment conditions are applied to ensure that the point of
interaction in an event is located inside the DeepCore fiducial volume. The
first two cuts are based on the radial $\rho_{\text{HLC}}$ and depth
$z_{\text{HLC}}$ positions from FirstHLC discussed in Appendix A.1. Only
events with $\rho_{\text{HLC}}$ less than 150 m and $z_{\text{HLC}}$ position
between $-475$ m and $-200$ m are kept. In addition, a similar criteria is
placed on the $\rho_{\text{Q1}}$ and $z_{\text{Q1}}$ vertex positions from the
first quartile (Q1) in the CoG algorithm (see Appendix A.4). Only events with
$\rho_{\text{Q1}}$ less than 150 m and $z_{\text{Q1}}$ position between $-475$
m and $-150$ m are kept.

After the above selection criteria, the number of atmospheric muon events is
reduced by more than 95% compared to Level 3. Further, the contamination due
to random detector noise is also significantly dropped.

###  D.2 Sample $\mathcal{B}$, Level 5

At Level 5, a boosted decision tree (BDT) Hocker et al. (2007) is trained to
further reduce the atmospheric muon background. Eleven variables are included
for training the BDT.

Figure 26: (top) Distribution of Q1-Q4 separation at Level 5 (before BDT score
cut applied). Each shaded color represents the stacked histogram from each
event type. Black dots represents the data distribution. MC events are
weighted by world averaged best fit oscillation parameters. (bottom) Ratio of
distribution from data to that from MC. Black error bars are the statistical
fluctuation from data, whereas shaded red areas are the uncertainties from
limited MC statistics.

Three variables for BDT training are obtained from the two quick
reconstruction algorithms described in Appendix A.5. They are the
reconstructed speed from iLineFit and the zenith angles from SPEFit11 and
iLineFit.

The next four BDT variables are related to the hit and charge information from
a cleaned pulseseries, including the charge ratios QR3 and C2QR3 discussed in
Appendix A.2. The number of hit DOMs in the cleaned hit series is also
included. Further, the total charge from a cleaned hit series of an event is
also used for BDT training. Since an atmospheric muon tends to deposit more
charges compared to a neutrino, the background-dominated, high-charge region
in the total charge distribution can help identify atmospheric muons from
neutrino events.

Figure 27: (top) BDT score distribution at Level 5 (before cut applied). Each
shaded color represents the stacked histogram from each event type. Black dots
represents the data distribution. The blue vertical line is the cut value of
0.2, events below which are rejected. MC events are weighted by world averaged
best fit oscillation parameters. (bottom) Ratio of distribution from data to
that from MC. Black error bars are the statistical fluctuation from data,
whereas shaded red areas are the uncertainties from limited MC statistics.

The last four BDT variables are based on the CoG algorithm discussed in
Appendix A.4. One of them is the separation between the first and the last
quartiles of a cleaned hit series, or Q1-Q4 separation. As shown in Fig. 26,
atmospheric muon background tends to have a longer spatial distance between
the two quartiles, compared to neutrino events. The remaining three variables
are reused from the previous level. They are the charge-weighted spread of
vertex $z$ position ($\sigma_{z}$) from all cleaned hits and the estimated
radial $\rho_{\text{Q1}}$ and depth $z_{\text{Q1}}$ positions of the
interaction vertex from the first quartile (Q1) of CoG.

Figure 27 shows the BDT score distribution, and a cut is applied to accept
events with a score above 0.2. Compared to Level 4, 99.9% of the atmospheric
muon background is removed, whereas 58% of all neutrinos is kept after the BDT
score cut is placed.

###  D.3 Sample $\mathcal{B}$, Level 6

In between Level 5 and 6, the comprehensive reconstruction discussed in Sec.
III.3 is performed. At Level 6, two final selection requirements are placed to
further improve the quality of the final sample.

Figure 28: Fraction of $\mu$ events to total at Level 6 (before cut applied)
in the starting radial $\rho$ and depth $z$ positions of reconstructed track
by the final reconstruction discussed in Sec. III.3. Red lines represent the
cut, within which events are kept.

First, final containment criteria is required based on the starting and
stopping positions of a reconstructed track. The track must starts within the
region defined by the red lines in Fig. 28, which further rejects atmospheric
muons coming through the dust layer. In addition, energy resolution can be
improved by ensuring the entire track is contained within the more densely
instrumented DeepCore region. Thus, the stopping radial $\rho_{\text{stop}}$
and depth $z_{\text{stop}}$ positions (see Fig. 29) of the track must satisfy
$\rho_{\text{stop}}\leq$ 150 m and $-500$ m $\leq\ z_{\text{stop}}\ \leq$
$-200$ m respectively.

Figure 29: (top) Distribution of stopping $z_{\text{stop}}$ positions of the
reconstructed track at Level 6 (before cut applied). Each shaded color
represents the stacked histogram from each event type. Black dots represents
the data distribution. The vertical lines are the cut values; events between
the two vertical lines are kept. MC events are weighted by world averaged best
fit oscillation parameters. (bottom) Ratio of distribution from data to that
from MC. Black error bars are the statistical fluctuation from data, whereas
shaded red areas are the uncertainties from limited MC statistics.

Last, to further reduce the contamination from sneaky penetrating atmospheric
muons, the number of corridor DOMs, defined in Appendix A.3, is required to be
less than 2.

Table 3 shows the event rates for each event type at each selection level. The
neutrino rates are the combination of the NC+CC channels and use the
atmospheric neutrino flux predictions from Honda et al. (2015) with world
average best fit values for the oscillation parameters. At the final level,
the rate of atmospheric muons is reduced by a factor of $\approx 10^{8}$, and
neutrino events contribute to $\approx 95$% of the $\mathcal{B}$ sample.

##  Appendix E Nuisance and Physics Parameter Impacts in Analysis
$\mathcal{A}$

This section shows the change in expected event rates for all systematic
uncertainties included in Analysis $\mathcal{A}$. The figures show the
percentage change when the corresponding nuisance parameter is shifted off
nominal by the amount specified in each subcaption.

(a) $\nu_{e}/\nu_{\mu}$ Flux Ratio +5% (b) $\nu/\bar{\nu}$ Flux Ratio
+1$\sigma$ (c) Upward/Horizontal Flux Ratio +1$\sigma$

(d) Neutrino Spectral Index +0.1 (e) $M_{A}^{CCQE}$ +1$\sigma$ (+0.248 GeV)

(f) $M_{A}^{RES}$ +1$\sigma$ (+0.22 GeV)

Figure 30: Changes in the event rates from a subset of the systematic
uncertainties included in the analyses.

(a) $\nu^{NC}$ Normalization +20% (b) $\theta_{23}$ shifted from 40∘ to 45∘

(c) $\Delta m^{2}_{32}$ shifted from 2.526 to 2.778 ($10^{-3}\rm{eV}^{2}$)
(d) Overall Optical Efficiency +10%

(e) Lateral Optical Efficiency +1$\sigma$ (f) Head-on Optical Efficiency +1
(a.u.)

Figure 31: Changes in the event rates from a subset of the systematic
uncertainties included in the analyses.

(a) Bulk Ice Scattering +10% (b) Bulk Ice Absorption +10%

(c) Atmospheric Muon Fraction shifted from 8.1% to 8.9% (d) Muon Spectral
Index +1$\sigma$

(e) Coincident $\mu$+$\nu$ Fraction +10%

Figure 32: Changes in the event rates from a subset of the systematic
uncertainties included in the analyses.