#  Skyrmion nucleation on the surface of a topological insulator

Daichi Kurebayashi<EMAIL_ADDRESS>School of Physics, The
University of New South Wales, Sydney 2052, Australia  Center for Emergent
Matter Science, RIKEN, Wako 351-0198, Japan  Oleg A. Tretiakov
<EMAIL_ADDRESS>School of Physics, The University of New South Wales,
Sydney 2052, Australia

###### Abstract

Skyrmion nucleation induced by spin-transfer torques at an interface of a
topological insulator and a ferromagnetic insulator is investigated. Due to
strong spin-orbit coupling on a surface of topological insulators, which
enhances the effect of spin torques, efficient manipulation of skyrmions is
expected, and therefore, topological insulators could provide the ideal
platform to achieve high-performance skyrmionic devices. Using micromagnetic
simulations and energetics, we evaluate properties of the skyrmion nucleation
on a surface of topological insulators, such as nucleation time, critical
electric field, and skyrmion numbers. We show that the nucleation time is
inversely proportional to the applied electric field. We also identify the
Gilbert damping and temperature dependencies of the critical field.
Furthermore, we analytically evaluate the effect of the Dzyaloshinskii-Moriya
interaction and demonstrate that the temperature dependence can be explained
by the reduction of a magnon excitation gap due to the self-energy
corrections.

A magnetic skyrmion is a real-space topological object defined by its non-
coplanar spin texture [1, 2, 3, 4, 5, 6, 7]. Because of their topological spin
structure, skyrmions exhibit unique dynamics [8, 9] and transport properties
associated with a nontrivial quantum Berry phase resulting in emergent
electromagnetic fields [10, 11, 12]. On the one hand, a topological insulator
(TI) is a momentum-space topological object characterized by its nontrivial
band structure [13, 14]. Due to this topologically nontrivial structure,
transport phenomena related to the momentum-space Berry curvature, such as the
quantum spin Hall effect [15, 16, 17] and the quantum anomalous Hall effect
[18, 19, 20], are realized and experimentally observed. Another consequence of
the nontrivial bands is the appearance of metallic gapless surface states. The
surface Dirac electrons mediate strong correlations between spin and current
as their spin and momentum have a one-to-one correspondence, known as the
spin-momentum locking. Although both TIs and skyrmions separately have been
recent emergent topics in condensed matter physics, a combination of them
could be an ideal platform to study the interplay of real- and momentum-space
topology.

Recently, in a heterostructure consisting of a TI and a magnetic insulator,
the ferromagnetic (FM) skyrmion formation has been observed by transport
measurements [21, 22, 23]. In addition to the conventional anomalous Hall
effect, an extra Hall signal has been observed, which was attributed to the
topological Hall effect arising from the emergent electromagnetic field of
skyrmions. Beyond that, a real-space observation by a scanning transmission
X-ray microscopy has been made and confirmed the formation of
antiferromagnetically coupled Néel-type skyrmions at a TI interface with a
ferrimagnet [24].

Skyrmions on a TI surface are exhilarating not only from the viewpoint of
emerging physics arising from the interplay of real- and momentum-space
topologies but also due to their prospects for spintronic nanodevices [6, 7].
It has been theoretically proposed that a skyrmion on a TI surface is
accompanied by a nonzero charge density due to the chiral edge states [25, 26,
27, 28]. Because of this additional charge attached to it, a skyrmion can be
manipulated by external electric fields without Ohmic losses from currents.
Another mechanism has been proposed to manipulate skyrmions by utilizing spin-
transfer torques, which are greatly enhanced due to the spin-momentum locking
on the TI surface [29, 30]. Consequently, the dynamics of skyrmions is
expected to be faster, which is also highly favorable for memory applications.

Although the skyrmion dynamics on a TI surface has been intensively
investigated, the missing ingredient for successful applications of skyrmions
in TIs is their nucleation studies. In conventional FMs, skyrmion nucleation
has been explored [31], for example, by employing geometric structures and
local magnetic fluctuations, such as notches, edges, or impurity sites [32,
33, 34, 35, 36, 37, 38] as well as by utilizing local injection of charge and
spin currents [39, 40, 41, 42]. However, the effect of dissipation on the
nucleation process, including the influence of the Gilbert damping and thermal
fluctuations, has not been understood well and therefore requires further
investigation. Moreover, skyrmion nucleation on a TI surface might be
significantly different from the one in conventional FMs because of the spin-
momentum locking. With increasing interest in skyrmionics with TIs, detailed
studies on the skyrmion nucleation process are highly demanded and therefore
are the subject of this paper. We investigate the general properties of
skyrmion nucleation on a TI surface, such as nucleation time and critical
field. We perform micromagnetic simulations at a finite temperature and
demonstrate the Gilbert damping and temperature dependences of the nucleation
process. To give physical understanding, we describe analytically the effect
of temperature based on a self-energy renormalization. By treating the
Dzyaloshinskii-Moriya interaction (DMI) as a perturbation within the random-
phase approximation, we succeed in reproducing the temperature dependence
obtained in our micromagnetic simulations.

Stochastic magnetization dynamics model.– Magnetization dynamics at finite
temperature is analyzed with the stochastic Landau-Lifshitz-Gilbert (sLLG)
equation [43, 44],

$\displaystyle\frac{\partial\bm{n}}{\partial
t}=-\frac{\gamma}{1+\alpha^{2}}\left[\bm{n}\times\bm{B}^{e}+\alpha\bm{n}\times\left(\bm{n}\times\bm{B}^{e}\right)\right]+\bm{T},$
(1)

where $\gamma$ is the gyromagnetic ratio, $\alpha$ is the Gilbert damping
constant, $\bm{n}(\bm{r})=\bm{M}(\bm{r})/M_{s}$ is the normalized
magnetization, $M_{s}$ is the saturation magnetization,
$\bm{B}^{e}(\bm{r})=-(1/M_{s})\delta F_{M}/\delta\bm{n}(\bm{r})+\bm{B}^{\rm
th}$ is the effective magnetic field, $\bm{B}^{\rm th}$ is the thermal field,
$F_{M}$ is the magnetic free energy, and $\bm{T}$ is the spin-transfer torque.
On a surface of TIs, the magnetic free energy $F_{M}=F_{\rm ex}+F_{\rm
ani}+F_{\rm Z}+F_{\rm DMI}$, where $F_{\rm ex}=(2J_{0}S^{2}/l_{a})\int
dV\,(\bm{\nabla}\bm{n})^{2}$, $F_{\rm ani}=-K\int dV\,(n_{z})^{2}$, $F_{\rm
Z}=-M_{s}\int dV\,B_{z}n_{z}$, and $F_{\rm DMI}$ are the exchange, anisotropy,
Zeeman, and DMI energies, respectively. Here $J_{0}$ is the exchange constant
between local magnetic moments, $l_{a}$ is the magnetic lattice constant, $S$
is the amplitude of the local moments, $K$ is the easy-axis anisotropy
constant along $z$-axis, and $B_{z}$ is the magnetic field perpendicular to
the film 111Note that we have neglected the effect of demagnetizing field as
it can be included as a small addition to an effective perpendicular magnetic
anisotropy of a thin film.. On a TI surface, electrons mediate anisotropic
exchange interaction, namely the DMI, reflecting an inversion symmetry
breaking. Then the DMI takes the form $F_{\rm
DMI}=(D/\xi_{l})\int\\!dV[n_{z}(\partial_{x}n_{x}+\partial_{y}n_{y})-n_{x}\partial_{x}n_{z}-n_{y}\partial_{y}n_{z}]$,
where $\xi_{l}$ is the penetration length of a TI surface state into the
magnetic insulator,
$D=-J^{2}\left[\Theta\left(|J|+E_{F}\right)-\Theta\left(|J|-E_{F}\right)\right]/(8\pi
v_{F})$ is the DMI constant mediated by the surface Dirac electrons [46],
where $J$ is the s-d coupling constant between electron’s spins and local
moments, $v_{F}$ is the Fermi velocity, $E_{F}$ is the Fermi energy, and
$\Theta(x)$ is the Heaviside function. The effective magnetic field is then
given by

$\displaystyle\\!\\!B^{e}_{x,y}$ $\displaystyle=$
$\displaystyle\frac{4J_{0}S^{2}}{M_{s}l_{a}}\Delta
n_{x,y}+\frac{2D}{M_{s}\xi_{l}}\partial_{x,y}n_{z}+B^{\rm th}_{x,y},$ (2)
$\displaystyle\\!\\!B^{e}_{z}$ $\displaystyle=$
$\displaystyle\frac{4J_{0}S^{2}}{M_{s}l_{a}}\Delta
n_{z}+B_{z}-\frac{2D}{M_{s}\xi_{l}}\bm{\nabla}\cdot\bm{n}+\frac{2K}{M_{s}}n_{z}+B^{\rm
th}_{z}.$ (3)

It is known that the spin-transfer torques on a TI surface that couples to
local magnetic moments are substantially modified due to the strong spin-orbit
coupling:

$\bm{T}=-\frac{e\gamma}{M_{s}\xi_{l}}\left(\alpha_{e}\hat{x}+\beta_{e}\hat{y}\right)\left(\bm{\nabla}\cdot\bm{n}\right)E_{x},$
(4)

where $\alpha_{e}=\tau J^{3}{\rm sgn}(E_{F})\Theta(E_{F}^{2}-J^{2})/(8\pi
E_{F}^{2})$ and $\beta_{e}=\tau^{2}J^{2}(E_{F}^{2}-J^{2}){\rm
sgn}(E_{F})\Theta(E_{F}^{2}-J^{2})/(8\pi E_{F}^{2})$ are dimensionless
coefficients [29, 30]. These specific spin-transfer torques occur at the
TI/magnetic insulator interface or in magnetic TI thin films. A detailed
derivation of the spin-transfer torques on the TI surface is given in the
Appendix. Note that we neglected the effect of the spin-orbit torques (SOT) on
the magnetization dynamics in Eq. (1). On the TI surface, the SOT plays an
important role in the magnetic dynamics because of the spin-momentum locking,
however, for the skyrmion nucleation, which requires spatial inhomogeneity,
the SOT does not give rise to any qualitative differences as the SOT are
uniform contribution. Instead, the SOT under the DC current is equivalent to a
static uniform in-plane magnetic field, therefore reducing the energy barrier
between the uniform FM and skyrmion states and thus contributing to the
reduction of the critical field.

For micromagnetic simulations, we discretize a space into a square lattice by
using relations
$\left.\partial_{i}\bm{n}(\bm{r})\right|_{\bm{r}_{j}}\approx(\bm{n}_{\bm{j}+\bm{e}_{i}}-\bm{n}_{\bm{j}-\bm{e}_{i}})/(2l_{a})$,
$\left.\partial_{i}^{2}\bm{n}(\bm{r})\right|_{\bm{r}_{j}}\approx(\bm{n}_{\bm{j}+\bm{e}_{i}}-2\bm{n}_{\bm{j}}+\bm{n}_{\bm{j}-\bm{e}_{i}})/l_{a}^{2}$,
and imposing periodic boundary conditions in $x$\- and $y$-directions. To
numerically simulate finite temperature, the thermal field is defined by
$\bm{B}^{\rm th}=\bm{\eta}\sqrt{2\alpha k_{B}T/(M_{s}\gamma V\Delta t)}$,
where $\bm{\eta}$ is the random vector drawn from a standard normal
distribution, $V$ is the average magnetic-ion volume, and $\Delta t$ is the
simulation time-step. As an integration scheme, we have employed the Heun’s
method with the time-step $\Delta t\sim 3$ fs. Since these simulations are all
stochastic, all numerical data are obtained as the statistical average over 50
independent simulations. For typical parameters of a magnetic TI, we have used
values estimated from the first-principle calculations and experiments:
$v_{F}=2.55\ \rm eV\,\AA$, $J=0.15$ eV, $J_{0}=1.38\times 10^{-23}$ J,
$K=7.25\times 10^{-27}$ J/$\AA^{3}$, $M_{s}=1.16\times 10^{4}$ J/(Tm${}^{3})$,
and $l_{a}=\xi_{l}=8.1$ Å [47, 48].

Figure 1: Magnetic profile at (a) $t=0$ s when no current is applied, (b)
$t=0.2$ ns after current pulse is applied, (c) $t=1$ ns after current pulse is
applied, and (d) after the current pulse is switched off. The color code and
arrows show the $z$\- and in-plane components of magnetization, respectively.
The parameters are $\alpha=0.04$, $T=0.2$K, and $E_{x}=6.5\times 10^{4}$ V/m.

Skyrmion nucleation.– We first examine the skyrmion nucleation with uniform
currents. Figure 1 shows the magnetic profiles under a uniform current pulse.
Here, we chose the parameters as $\alpha=0.04$, $T=0.2$ K, and
$E_{x}=6.5\times 10^{4}$ V/m. At $t=0$ ns, when no current is applied, the
magnetization is along the $z$-axis. After application of the current pulse,
it first develops particle-like small fluctuations as shown in Fig. 1(b).
Eventually, these particle-like fluctuations grow into Néel skyrmions once
their radius reaches the critical one determined by $J$, $D$, and $B_{z}$ [49,
50], otherwise they collapse back into a uniform state. Note that, during an
early stage of the nucleation, a skyrmion - antiskyrmion pair is created due
to the topological number conservation [51, 11]. However, because our DMI
stabilizes only skyrmions, the antiskyrmions quickly decay into the uniform
state [52]. After the current pulse is switched off, only the skyrmions
survive, see Fig. 1 (d). These results clearly show that skyrmions can be
nucleated by uniform current pulses at TI/FM interfaces [53].

Effect of Gilbert damping.– There is always a delay before the first skyrmion
is nucleated. To investigate this nucleation time, we first study its Gilbert
damping $\alpha$ dependence. Since $\alpha$ depends on various factors such as
disorder, it is important to understand how it affects the nucleation process.
The skyrmion nucleation time $t_{\rm n}$ for various $\alpha$ is shown in Fig.
2 (a) as a function of applied field $E_{x}$. We define the $t_{\rm n}$ as a
time before the total skyrmion number, $N_{\rm sk}=\left|\iint dxdy\
\bm{n}\cdot(\partial_{x}\bm{n}\times\partial_{y}\bm{n})\right|/(4\pi^{2}),$
exceeds one. The skyrmion nucleation is absent for small $E_{x}$, i.e., there
is a critical field $E_{c}$ for the nucleation process. In this regime, the
energy dissipation caused by the Gilbert damping exceeds the energy influx due
to the spin-transfer torque 222The spin-transfer torque is derived as a
response to an electric field $E$, instead of current $j$. In order to
determine the energy efficiency, it is required to estimate the current
density, which depends on sample quality and detailed material parameters.
Therefore, a theoretical estimation of energy efficiency is not
straightforward. However, in general, surface states of TI are known to be
less affected by impurity scattering and possess higher mobility than those of
conventional metals. Consequently, the energy loss from the Joule heating is
expected to be smaller., such that the total accumulated energy is
insufficient to nucleate a skyrmion.

In terms of the nucleation time, Fig. 2 (a) shows diverging behavior at
$E_{c}$ and monotonically decreases with $E_{x}$. This can be explained based
on the energy considerations. Since the energy influx per unit time to the
system is linearly proportional to $E_{x}$ as the coupling to electrons is
treated within the linear response theory, the total accumulated energy is
$F_{\rm tot}\propto E_{x}$. The nucleation rate $1/t_{\rm n}$ is, then,
proportional to the energy difference between the total accumulated energy
$F_{\rm tot}$ and the nucleation energy of a single skyrmion $F_{\rm sk}$,
i.e., $1/t_{\rm n}\propto F_{\rm tot}-F_{\rm sk}\sim E_{x}-E_{c}$. Thus, the
nucleation time scales as $t_{\rm n}\propto(E_{x}-E_{c})^{-1}$. Indeed, the
numerical data is in excellent agreement with it, see Fig. 2 (a), when fitted
by,

$t_{n}(\alpha)=A[|E_{x}|-E_{c}(\alpha)]^{-1},$ (5)

where $A$ is a coefficient and $E_{c}(\alpha)$ is the critical field at the
given $\alpha$. As shown in the inset of Fig. 2 (a), the critical field
$E_{c}$ is linear in $\alpha$. This is because the energy dissipation is
linear in $\alpha$ and the energy influx is $\propto E_{x}$. Thus, as the
total accumulated energy is determined as the difference of the energy influx
due to spin-transfer torques and the Gilbert dissipation, the required energy
influx to nucleate a skyrmion should linearly increase with $\alpha$.

Figure 2: (a) The nucleation time and (b) the skyrmion number as a function of
electric field $E_{x}$ are plotted for various Gilbert damping $\alpha$ and
$T=1$ K. The dots are numerical data and solid lines are fitting functions.
The inset of panel (a) shows the Gilbert damping dependence of the critical
field $E_{c}$ in the units of $10^{4}$ V/m.

We also examined the total skyrmion number nucleated after 1 ns pulse as a
function of applied field $E_{x}$ for several $\alpha$. As shown in Fig. 2
(b), the nucleated skyrmion number $N_{\rm sk}$ linearly increases with
$E_{x}$ in the vicinity of the critical field. However, as $E_{x}$ increases
further, $N_{\rm sk}$ deviates from a linear slope and saturates. The
saturation occurs because skyrmions start overlapping and merging into large
domains as the skyrmion density increases. By further increasing $E_{x}$,
skyrmion states are almost destroyed by strong magnon excitations. Time
evolution of $N_{\rm sk}$ at each applied field and corresponding magnetic
profiles are shown in Fig. 3. As seen from Fig. 3 (a), $N_{\rm sk}$ stays
constant after reaching a steady state for the fields below $E_{x}\sim
1.5\times 10^{5}$ m/V. In this regime, all nucleated skyrmions are well
separated, as shown in Fig. 3 (b). On the other hand, for larger $E_{x}$
corresponding to the turbulence regime in Fig. 3 (a), $N_{\rm sk}$ decreases
as $E_{x}$ increases and oscillates in time. Figure 3 (c) shows a typical
magnetic profile in the strong field regime, swirling magnetic structures no
longer survive [53], and larger magnetic domains are formed because of strong
magnon excitations. Note that it is a crossover, not a phase transition,
between the steady nucleation and turbulence regimes.

Figure 3: (a) The skyrmion number $N_{\rm sk}$ as a function of time and
electric field $E_{x}$, for $\alpha=0.04$ and $T=0.2$ K. (b), (c) Magnetic
profile at $t=0.5$ ns and $E_{x}=7\times 10^{4}$ V/m, and (c) $t=0.5$ ns and
$E_{x}=22\times 10^{4}$ V/m.

Effect of temperature.– Next, we examine the temperature effects on the
nucleation phenomenon. The nucleation time for various temperatures is
presented in Fig. 4. One can notice that the temperature $T$ only affects the
critical field, see $E_{c}(\alpha)$ in Eq. (5), whereas the functional form of
$E_{x}$ is hardly modified. The critical field linearly decreases with $T$, as
shown in the inset of Fig. 4. Phenomenologically, the linear dependence on $T$
can be understood as follows. The thermal fluctuations supply the energy $\sim
k_{B}T$ to the system, where $k_{B}$ is the Boltzmann constant. Due to this
additional contribution, the energy required to create a skyrmion reduces
linearly with $T$. We note that the critical field vanishes around $T\sim 8$
K; above this temperature, skyrmions are nucleated even without the spin-
transfer torques due to the thermal fluctuations.

Figure 4: The nucleation time as a function of electric field $E_{x}$ for
various temperatures and $\alpha=0.04$. The dots are numerical data, while the
solid curves are fitting functions given by Eq. (5). The inset shows the
temperature dependence of the critical field $E_{c}$.

Although phenomenological energy considerations explain the temperature effect
on the nucleation process, we move one step further and try to explain this
phenomenon in terms of ferromagnetic magnon excitations. Introducing the
Holstein-Primakoff representation, the free energy of the system can be
transformed to the magnonic Hamiltonian as

$\displaystyle\hat{H}_{m}$ $\displaystyle=$
$\displaystyle\sum_{k}a_{\bm{k}}^{\dagger}\left(J_{k}+\tilde{K}+\tilde{B}-t_{k}\right)a_{\bm{k}}$
(6)
$\displaystyle+\sum_{\bm{k},\bm{q}}\left(D_{\bm{q}}a^{\dagger}_{\bm{q}}a^{\dagger}_{\bm{k}-\bm{q}/2}a_{\bm{k}+\bm{q}/2}+D^{*}_{\bm{q}}a^{\dagger}_{\bm{k}+\bm{q}/2}a_{\bm{k}-\bm{q}/2}a_{\bm{q}}\right),$

where $a_{\bm{k}}$ is the magnon annihilation operator with the wave number
$k$, $J_{k}=8J_{0}S^{2}\sum_{i=x,y}(1-\cos k_{j})$ is the exchange energy,
$\tilde{K}=2Kl_{a}^{3}$ is the easy-axis anisotropy,
$\tilde{B}=M_{s}B_{z}l_{a}^{3}$ is the Zeeman energy,
$D_{\bm{k}}=2iDl_{a}(\sin k_{x}+i\sin k_{y})$ is the DMI, and
$t_{\bm{k}}=4eE_{x}l_{a}(\beta_{e}\sin k_{y}+\alpha_{e}\sin k_{x})$ describes
the effect of the spin-transfer torque. Note that we have retained the linear
in $E_{x}$ term and neglected the higher-order terms. The first term gives a
single-particle magnon dispersion as it already has the bilinear form in
magnon operators. On the other hand, the DMI in the second line of Eq. (6)
contains three-magnon operators and has to be treated perturbatively.

The full Green’s function is given by
$G_{\bm{k}}(D)=\braket{\phi_{k}\bar{\phi}_{k}}=\frac{\int
D(\bar{\phi},\phi)\phi_{k}\bar{\phi}_{k}e^{-S[\phi,\bar{\phi}]}}{\int
D(\bar{\phi},\phi)e^{-S[\phi,\bar{\phi}]}}$, where
$S[\phi,\bar{\phi}]=\sum_{k}\bar{\phi}_{\bm{k}}(-i\omega_{n})\phi_{k}+\mathcal{H}_{m}[\phi,\bar{\phi}]$
is the imaginary-time action, $\phi_{k=(\omega_{n},\bm{k})}$ is the eigenvalue
of the magnon operator $a_{\bm{k}}$, and $\mathcal{H}_{m}$ is the Hamiltonian
in the magnon coherent states basis. Within the random phase approximation,
the Green’s function
$G_{k}\approx(-i\omega_{n}+J_{k}+\tilde{K}+\tilde{B}-t_{k}-\Sigma_{k})^{-1}$,
where $\Sigma_{k}$ is the self-energy induced by the DMI. The real part of the
self-energy modifies the magnon dispersion, while the imaginary part gives a
finite lifetime. Therefore, the effective magnon dispersion including the
effect of DMI takes the form

$\omega_{\rm eff}(\bm{k})=J_{\bm{k}}+\tilde{K}+\tilde{B}-t_{\bm{k}}-{\rm
Re}\left[\Sigma_{k}\right].$ (7)

As depicted in Fig. 5, the self-energy has two contributions
$\Sigma_{k}=\Sigma_{k,d}+\Sigma_{k,p}$ with

$\displaystyle\Sigma_{k,d}$ $\displaystyle=$
$\displaystyle-4|D_{\bm{k}}|^{2}\sum_{\bm{q}}\frac{f_{B}(\omega_{\bm{q}})-f_{B}(\omega_{\bm{k}+\bm{q}})}{\omega_{\bm{q}}-\omega_{\bm{k}+\bm{q}}+i0_{+}},$
(8) $\displaystyle\Sigma_{k,p}$ $\displaystyle=$ $\displaystyle
4\sum_{\bm{q}}|D_{\bm{q}}|^{2}\frac{1+f_{B}(\omega_{\bm{q}})+f_{B}(\omega_{\bm{k}-\bm{q}})}{\omega_{\bm{q}}+\omega_{\bm{k}-\bm{q}}-i0_{+}},$
(9)

where $f_{B}$ is the Bose-Einstein distribution and
$\omega_{\bm{k}}=J_{\bm{k}}+\tilde{K}+\tilde{B}-t_{\bm{k}}$ is the bare magnon
dispersion. Equations (8) and (9) correspond to the magnon density-density
response function, $\Sigma_{d}$, and pair-correlation function, $\Sigma_{p}$.
From Eqs. (8) and (9), one notices that the real part of the self-energy is
always positive, namely, the magnon correlations always reduce the magnon
excitation gap min($\omega_{\rm eff}$). As the Bose-Einstein distribution
function can be expanded as $f_{B}(\omega)\approx k_{B}T/\omega$, the
reduction of the gap min($\omega_{\rm eff}$) due to the self-energy
contributions is linear in temperature. Then, because the magnon instability
is a precursor of the skyrmion nucleation and it occurs when the magnon
excitation gap collapses, we conclude that the critical field linearly
decreases with $T$. Note that this expansion is valid when a single-particle
magnon excitation gap, $\omega_{\rm gap}$, is smaller than $k_{B}T$. In
typical magnetically-doped TIs, the perpendicular magnetic anisotropy is $\sim
10^{-6}$ eV [48], which corresponds to $\sim 0.1$ K. Because the condition
$\omega_{\rm gap}<k_{B}T$ is always met within the temperature range of our
simulations, the linear-$T$ dependence observed in the inset of Fig. 4 can be
explained analytically based on the magnon excitations. Note that the
instability and the nucleation of magnetic texture is a highly nonequilibrium
process, and equilibrium analysis based on the self-energy corrections may
fail to completely explain dynamical properties. We also neglected the magnon
excitation due to the electric field, which may affect the behavior in a low-
temperature regime. However, it is still capable of qualitatively describing
statistical properties such as the average nucleation time and critical field.

Figure 5: The Feynman diagrams contributing to the magnon self-energy. (a) and
(b) show the density-density and the pair correlations, respectively.

Discussion.– Let us review potential candidate materials for magnetic TIs. In
this study, we used material parameters for a typical magnetic TI, such as Cr-
and V-doped $\rm(Bi,Sb)_{2}Te_{3}$ [47, 48]. The Curie temperature of this
group of materials, however, is $T_{c}\approx 15\sim 30\rm K$, and therefore
the skyrmion state is only stable at the cryogenic temperatures, thus
hindering some device applications. In order to overcome this limitation,
achieving high $T_{c}$ is particularly desirable. In a similar group of
materials, such as $\rm MnSb_{2}Te_{4}$ [55], a higher Curie temperature
$T_{c}\sim 50$K is observed. Furthermore, the increase of $T_{c}$ by 60K due
to the exchange bias effect is observed in a superlattice structure of Dy-
doped $\rm Bi_{2}Te_{3}$ and Cr-doped $\rm Sb_{2}Te_{3}$ [56]. As another
approach to realizing high $T_{c}$ magnetic TIs, bilayer systems consisting of
a TI and magnetic insulator with a higher $T_{c}$ are expected to be more
commercially realizable candidates. It is reported that a TI grown on YIG by
molecular-beam epitaxy exhibits the anomalous Hall effect, i.e., a signature
of magnetically-coupled surface states, up to 300K [57, 58]. Thus, with recent
advances in nanotechnology, there is a real prospect that room-temperature
magnetic TIs will be within the reach in the near future.

In conclusion, we have comprehensively studied the current-induced skyrmion
nucleation on a surface of a topological insulator. As a system, we have
considered a heterostructure consisting of a TI and an ultrathin ferromagnetic
insulator (applicable to the recent experimental realizations [21, 22, 23,
24]) and employed the specific spin-transfer torque realized on a TI surface
coupling to magnetic local moments. By solving the stochastic LLG equation, we
have determined the critical field $E_{c}$ for skyrmion nucleation and found
that the nucleation time $t_{n}\propto(E_{x}-E_{c})^{-1}$, where $E_{c}$ has
been found to be proportional to the Gilbert damping $\alpha$. These results
suggest that with the advances in nanotechnology, when much cleaner samples
with lower damping will be available, faster skyrmion nucleation with lower
fields may be achieved. Furthermore, we have investigated the temperature
dependence of the nucleation and have observed the temperature effects on the
critical field, which linearly decreases as temperature rises. We have given a
quantum microscopic description for this linear dependence based on magnon
excitations. The self-energy corrections due to the Dzyaloshinskii-Moriya
interaction provide linear in temperature reduction of the magnon excitation
gap and the critical field. Our results give a phenomenological understanding
of the skyrmion nucleation on a surface of a topological insulator and open
doors for developing TI-based skyrmionic memory and logic nanodevices.

## Acknowledgments

The authors are grateful to N. Nagaosa for insightful discussions. O.A.T.
acknowledges the support from the Australian Research Council (Grant No.
DP200101027), the Cooperative Research Project Program at the Research
Institute of Electrical Communication, Tohoku University (Japan), and the
NCMAS grant.

*

##  Appendix A Spin-transfer torque on the surface of a TI

In the appendix, we present a detailed derivation of spin-transfer torques on
the surface of a topological insulator. The Hamiltonian describing the surface
states of a topological insulator is given by

$\displaystyle H$ $\displaystyle=$ $\displaystyle
H_{0}-ev_{F}\sum_{\bm{k}}\psi_{\bm{k}}^{\dagger}\left[\bm{A}(t)\times\bm{\sigma}\right]_{z}\psi_{\bm{k}}-J\sum_{\bm{r}}\bm{n}(\bm{r})\cdot\psi_{\bm{r}}^{\dagger}\bm{\sigma}\psi_{\bm{r}},$

where
$H_{0}=\sum_{\bm{k}}\psi_{\bm{k}}^{\dagger}\left[-v_{F}\left(k_{x}\sigma_{y}-k_{y}\sigma_{x}\right)-J\sigma_{z}\right]\psi_{\bm{k}}$
describes a massive Dirac Fermions with a helical spin texture [47],
$\psi_{\bm{k}}$ is the Fermionic annihilation operator at momentum $\bm{k}$,
$v_{F}$ is the Fermi velocity, $\bm{\sigma}$ is the Pauli matrix describing
real spin degrees of freedom, $\bm{A}=(A_{x},A_{y})$ is the electromagnetic
vector potential, $J$ is a s-d exchange constant, and
$\bm{n}(\bm{r})=\hat{\bm{z}}+\delta\bm{n}(\bm{r})$ is magnetization of local
moments. Here we take the saturated magnetization direction along
$\hat{\bm{z}}$ axis and assume $\delta\bm{n}(\bm{r})$ is a small magnetic
fluctuation around $\hat{\bm{z}}$ axis, namely
$\hat{\bm{z}}\cdot\delta\bm{n}(\bm{r})=0$ and $|\delta\bm{n}(\bm{r})|<<1$.

By treating the electromagnetic field and the magnetic fluctuation,
$\delta\bm{n}(\bm{r})$, as perturbations, a non-equilibrium spin accumulation
induced by an electric field applied along the $x$-direction is evaluated as

$\displaystyle\braket{\bm{\sigma}(\bm{r},t)}$ $\displaystyle=$
$\displaystyle-i\ {\rm
Tr}\left.\left[\bm{\sigma}G^{<}(\bm{r},\bm{r}^{\prime};t,t)\right]\right|_{\bm{r}^{\prime}\rightarrow\bm{r}},$
(11)

where
$G^{<}(\bm{r},\bm{r};t,t)=i\braket{\psi^{\dagger}_{\bm{r}^{\prime}}(t^{\prime})\psi_{\bm{r}}(t)}$
is the lesser Green’s function on the Keldysh contour and the trace is over
the spin indices [59]. Turning into momentum space, the contour ordered
Green’s function obeys the Dyson equation given as

$\displaystyle G_{\bm{k},\bm{k}^{\prime}}(t,t^{\prime})$ $\displaystyle=$
$\displaystyle g_{\bm{k}}(t-t^{\prime})\delta_{\bm{k},\bm{k}^{\prime}}$
$\displaystyle-ev_{F}\sum_{\bm{k}}\int_{C}dt_{1}\
g_{\bm{k}}(t-t_{1})A_{x}(t_{1})\sigma_{y}G_{\bm{k},\bm{k}^{\prime}}(t_{1},t^{\prime})$
$\displaystyle-J\sum_{\bm{k},\bm{q}}\int_{C}dt_{1}\
g_{\bm{k}}(t-t_{1})\delta\bm{n}(\bm{q})\cdot\bm{\sigma}G_{\bm{k}-\bm{q},\bm{k}^{\prime}}(t_{1},t^{\prime}),$

where
$G_{\bm{k},\bm{k}^{\prime}}=\sum_{\bm{k},\bm{k}^{\prime}}\exp(-i\bm{k}\cdot\bm{r}+i\bm{k}^{\prime}\cdot\bm{r}^{\prime})G(\bm{r},\bm{r}^{\prime})$
and $g_{\bm{k}}=[i\partial_{t}-H_{0}]^{-1}$ is an unperturbed Green’s
function. By expanding the contour-ordered Green’s function up to the first
order in the external field and the magnetic fluctuation,
$\delta\bm{n}(\bm{r})$, one can obtain that

$\displaystyle G_{\bm{k},\bm{k}^{\prime}}(t,t^{\prime})$
$\displaystyle\approx$ $\displaystyle-
ev_{F}J\delta_{\bm{k}-\bm{q},\bm{k}^{\prime}}\sum_{\bm{k},\bm{q}}\int_{C}dt_{1}\int_{C}dt_{2}\delta\bm{n}(\bm{q})\cdot\left[g_{\bm{k}}(t-t_{1})A_{x}(t_{1})\sigma_{y}g_{\bm{k}}(t_{1}-t_{2})\bm{\sigma}g_{\bm{k}-\bm{q}}(t_{2}-t^{\prime})\right.$
(13) $\displaystyle\hskip
147.95424pt\left.+g_{\bm{k}+\bm{q}}(t-t_{1})\bm{\sigma}g_{\bm{k}}(t_{1}-t_{2})A_{x}(t_{2})\sigma_{y}g_{\bm{k}}(t_{2}-t^{\prime})\right].$

Note that, in Eq. (13), we have neglected the terms irrelevant to the non-
equilibrium spin accumulation. By substituting Eq. (13) into Eq. (11), the
spin accumulation is given by

$\displaystyle\braket{\sigma_{\alpha}(\bm{r},t)}$ $\displaystyle=$
$\displaystyle
iev_{F}J\sum_{\bm{k},\bm{q}}e^{i\bm{q}\cdot\bm{r}}A_{x}(t_{1})\delta
n_{\beta}(\bm{q}){\rm
Tr}\left[\int_{C}dt_{1}\int_{C}dt_{2}\sigma_{\alpha}g_{\bm{k}}(t-t_{1})\sigma_{y}g_{\bm{k}}(t_{1}-t_{2})\sigma_{\beta}g_{\bm{k}-\bm{q}}(t_{2}-t)\right.$
(14) $\displaystyle\hskip
108.12054pt\left.+\int_{C}dt_{1}\int_{C}dt_{2}\sigma_{\alpha}g_{\bm{k}+\bm{q}}(t-t_{2})\sigma_{\beta}g_{\bm{k}}(t_{2}-t_{1})\sigma_{y}g_{\bm{k}}(t_{1}-t)\right]^{<}.$

With use of the Langreth theorem,

$\displaystyle\left[\int_{C}dt_{1}\int_{C}dt_{2}A(t-t_{1})B(t_{1}-t_{2})C(t_{2}-t^{\prime})\right]^{<}$
$\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}\left[A^{r}(t-t_{1})B^{r}(t_{1}-t_{2})C^{<}(t_{2}-t^{\prime})\right.$
(15)
$\displaystyle\left.+A^{r}(t-t_{1})B^{<}(t_{1}-t_{2})C^{a}(t_{2}-t^{\prime})+A^{<}(t-t_{1})B^{a}(t_{1}-t_{2})C^{a}(t_{2}-t^{\prime})\right],$

where superscripts $(a,r)$ denote the advanced and retarded functions, the
non-equilibrium spin accumulation is obtained as

$\displaystyle\braket{\sigma_{\alpha}(\bm{r},t)}=\sum_{\bm{q}}\int\frac{d\Omega}{2\pi}e^{i\bm{q}\cdot\bm{r}-i\Omega
t}\ \chi_{\alpha\beta}(\bm{q},\Omega)A_{x}(\Omega)\delta n_{\beta}(\bm{q}).$
(16)

Here $\chi_{\alpha\beta}(\bm{q},\Omega)$ is a spin susceptibility defined by

$\displaystyle\chi_{\alpha\beta}(\bm{q},\Omega)$ $\displaystyle=$
$\displaystyle iev_{F}J\ \int\frac{d\omega}{2\pi}\sum_{\bm{k}}{\rm
Tr}\left[\sigma_{\alpha}g_{\bm{k}}(\omega+\Omega)\sigma_{y}g_{\bm{k}}(\omega)\sigma_{\beta}g_{\bm{k}-\bm{q}}(\omega)+\sigma_{\alpha}g_{\bm{k}+\bm{q}}(\omega)\sigma_{\beta}g_{\bm{k}}(\omega)\sigma_{y}g_{\bm{k}}(\omega-\Omega)\right]^{<}$
(17) $\displaystyle\approx$
$\displaystyle\frac{i\Omega}{2\pi}ev_{F}J\sum_{\bm{k}}\left[\sigma_{\alpha}g^{r}_{\bm{k}}\sigma_{y}g_{\bm{k}}^{a}\sigma_{\beta}g_{\bm{k}-\bm{q}}^{a}+\sigma_{\alpha}g^{r}_{\bm{k}+\bm{q}}\sigma_{\beta}g_{\bm{k}}^{r}\sigma_{y}g_{\bm{k}}^{a}\right],$

where $g_{\bm{k}}^{r,a}\equiv g_{\bm{k}}^{r,a}(\omega=0)=\left[-H_{0}\pm
i\eta\right]^{-1}$, $\eta$ is the damping rate induced by impurity scattering,
and we have expanded the function with respect to the frequency up to the
first order. Here we have considered the small scattering regime
$(\eta<<E_{F},\ J)$ and only retained the dominant contributions which include
both the retarded and advanced Green’s functions. By expanding the Green’s
function with respect to the external momentum, $\bm{q}$, as
$g_{\bm{k}+\bm{q}}\approx
g_{\bm{k}}-v_{F}\epsilon_{\gamma\rho}q_{\gamma}g_{\bm{k}}\sigma_{\rho}g_{\bm{k}}+\mathcal{O}(q^{2})$,
one can obtain that

$\displaystyle\chi_{\alpha\beta}(\bm{q},\Omega)$ $\displaystyle\approx$
$\displaystyle\frac{i\Omega}{2\pi}ev_{F}^{2}J\sum_{\bm{k}}\epsilon_{\gamma\rho}\left[\sigma_{\alpha}g^{r}_{\bm{k}}\sigma_{y}g_{\bm{k}}^{a}\sigma_{\beta}g_{\bm{k}}^{a}\sigma_{\rho}g_{\bm{k}}^{a}\right.$
(18)
$\displaystyle\left.-\sigma_{\alpha}g^{r}_{\bm{k}}\sigma_{\rho}g^{r}_{\bm{k}}\sigma_{\beta}g_{\bm{k}}^{r}\sigma_{y}g_{\bm{k}}^{a}\right]q_{\gamma}.$

Note that we have dropped the $\mathcal{O}(q^{0})$-terms because they are the
uniform contributions, namely corresponding to the spin-orbit torque, which is
well known and there is no need to reproduce in this Appendix. After momentum
integration in Eq. (18), the non-equilibrium spin accumulation corresponding
to the spin-transfer torques is obtained by

$\displaystyle\braket{\bm{\sigma}_{\perp}(\bm{r},t)}$ $\displaystyle=$
$\displaystyle\bm{C}(E_{F})\bm{\nabla}\cdot\delta\bm{n}(\bm{r})E_{x},$ (19)

where $\bm{C}(E_{F})=\frac{eJ}{8\pi
E_{F}^{2}}\left[\tau^{2}(E_{F}^{2}-J^{2})\hat{x}-\tau JM_{z}\hat{y}\right]{\rm
sgn(E_{F})}\Theta(E_{F}^{2}-J^{2})$, $\tau=\frac{1}{2\eta}$ is the electron
scattering lifetime, ${\rm sgn}(x)$ is the sign function, and $\Theta(x)$ is
the Heaviside function. Finally, the spin-transfer torques are given by

$\displaystyle\bm{T}$ $\displaystyle=$
$\displaystyle-\frac{\gamma_{0}J}{M_{S}\xi_{l}}\hat{\bm{z}}\times\braket{\bm{\sigma}_{\perp}(\bm{r},t)},$
(20) $\displaystyle T_{x}$ $\displaystyle=$
$\displaystyle-\frac{e\gamma_{0}\alpha_{e}}{M_{S}\xi_{l}}\left(\bm{\nabla}\cdot\delta\bm{n}\right)E_{x},$
(21) $\displaystyle T_{y}$ $\displaystyle=$
$\displaystyle-\frac{e\gamma_{0}\beta_{e}}{M_{S}\xi_{l}}\left(\bm{\nabla}\cdot\delta\bm{n}\right)E_{x},$
(22)

where $\alpha_{e}=\frac{\tau J^{3}}{8\pi E_{F}^{2}}{\rm
sgn}(E_{F})\Theta(E_{F}^{2}-J^{2})$ and
$\beta_{e}=\frac{\tau^{2}J^{2}(E_{F}^{2}-J^{2})}{8\pi E_{F}^{2}}{\rm
sgn}(E_{F})\Theta(E_{F}^{2}-J^{2})$ are dimensionless coefficients
characterizing the spin-transfer torque mediated by the Dirac electrons [29,
30]. Note that Eq. (21) and (22) correspond to the adiabatic and non-adiabatic
torques, respectively. In Fig. 6, the Fermi energy dependence of the spin-
transfer torque coefficients, $\alpha_{e}$ and $\beta_{e}$, are shown. The
results indicate that the non-adiabatic contribution dominates over the
adiabatic one. It is known that the non-adiabatic contribution becomes
significant when the spin-relaxation is large in the system [60]. Due to the
strong spin-orbit coupling on a surface of topological insulators, it is
reasonable to have substantial non-adiabatic contributions.

Figure 6:  The coefficients characterizing the spin-transfer torques on the
surface of the topological insulator. The material parameters are taken as
$\frac{1}{2\tau}=\eta_{0}|E_{F}|$, $\eta_{0}=0.1$, $v_{F}=2.55\,{\rm
eV\,\AA}$, and $J=0.03$ eV.

## References

  * Pokrovsky [1979] V. Pokrovsky, Properties of ordered, continuously degenerate systems, Adv. Phys. 28, 595 (1979).
  * Bogdanov and Yablonskii [1989] A. N. Bogdanov and D. A. Yablonskii, Thermodynamically stable “vortices” in magnetically ordered crystals. The mixed state of magnets, Sov. Phys. JETP 68, 101 (1989).
  * Mühlbauer  _et al._ [2009] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer, A. Rosch, A. Neubauer, R. Georgii, and P. Böni, Skyrmion lattice in a chiral magnet, Science 323, 915 (2009).
  * Yu  _et al._ [2010a] X. Z. Yu, Y. Onose, N. Kanazawa, J. H. Park, J. H. Han, Y. Matsui, N. Nagaosa, and Y. Tokura, Real-space observation of a two-dimensional skyrmion crystal, Nature 465, 901 (2010a).
  * Finocchio  _et al._ [2016] G. Finocchio, F. Büttner, R. Tomasello, M. Carpentieri, and M. Kläui, Magnetic skyrmions: from fundamental to applications, J. Phys. D Appl. Phys. 49, 423001 (2016).
  * Fert  _et al._ [2017] A. Fert, N. Reyren, and V. Cros, Magnetic skyrmions: advances in physics and potential applications, Nat. Rev. Mater. 2, 17031 (2017).
  * Göbel  _et al._ [2021] B. Göbel, I. Mertig, and O. A. Tretiakov, Beyond skyrmions: Review and perspectives of alternative magnetic quasiparticles, Phys. Rep. 895, 1 (2021).
  * Jiang  _et al._ [2017] W. Jiang, X. Zhang, G. Yu, W. Zhang, X. Wang, M. B. Jungfleisch, J. E. Pearson, X. Cheng, O. Heinonen, K. L. Wang, Y. Zhou, A. Hoffmann, and S. G. E. te Velthuis, Direct observation of the skyrmion Hall effect, Nat. Phys. 13, 123 (2017).
  * Litzius  _et al._ [2017] K. Litzius, I. Lemesh, B. Krüger, P. Bassirian, L. Caretta, K. Richter, F. Büttner, K. Sato, O. A. Tretiakov, J. Förster, R. M. Reeve, M. Weigand, I. Bykova, H. Stoll, G. Schütz, G. S. D. Beach, and M. Kläui, Skyrmion Hall effect revealed by direct time-resolved X-ray microscopy, Nat. Phys. 13, 170 (2017).
  * Schulz  _et al._ [2012] T. Schulz, R. Ritz, A. Bauer, M. Halder, M. Wagner, C. Franz, C. Pfleiderer, K. Everschor, M. Garst, and A. Rosch, Emergent electrodynamics of skyrmions in a chiral magnet, Nat. Phys. 8, 301 (2012).
  * Nagaosa and Tokura [2013] N. Nagaosa and Y. Tokura, Topological properties and dynamics of magnetic skyrmions, Nat. Nanotechnol. 8, 899 (2013).
  * Akosa  _et al._ [2018] C. A. Akosa, O. A. Tretiakov, G. Tatara, and A. Manchon, Theory of the topological spin Hall effect in antiferromagnetic skyrmions: Impact on current-induced motion, Phys. Rev. Lett. 121, 097204 (2018).
  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * Kane and Mele [2005] C. L. Kane and E. J. Mele, ${Z}_{2}$ topological order and the quantum spin Hall effect, Phys. Rev. Lett. 95, 146802 (2005).
  * Bernevig and Zhang [2006] B. A. Bernevig and S.-C. Zhang, Quantum spin Hall effect, Phys. Rev. Lett. 96, 106802 (2006).
  * König  _et al._ [2007] M. König, S. Wiedmann, C. Brüne, A. Roth, H. Buhmann, L. W. Molenkamp, X.-L. Qi, and S.-C. Zhang, Quantum spin Hall insulator state in HgTe quantum wells, Science 318, 766 (2007).
  * Yu  _et al._ [2010b] R. Yu, W. Zhang, H.-J. Zhang, S.-C. Zhang, X. Dai, and Z. Fang, Quantized anomalous Hall effect in magnetic topological insulators, Science 329, 61 (2010b).
  * Chang  _et al._ [2013] C.-Z. Chang, J. Zhang, X. Feng, J. Shen, Z. Zhang, M. Guo, K. Li, Y. Ou, P. Wei, L.-L. Wang, Z.-Q. Ji, Y. Feng, S. Ji, X. Chen, J. Jia, X. Dai, Z. Fang, S.-C. Zhang, K. He, Y. Wang, L. Lu, X.-C. Ma, and Q.-K. Xue, Experimental observation of the quantum anomalous Hall effect in a magnetic topological insulator, Science 340, 167 (2013).
  * Checkelsky  _et al._ [2014] J. G. Checkelsky, R. Yoshimi, A. Tsukazaki, K. S. Takahashi, Y. Kozuka, J. Falson, M. Kawasaki, and Y. Tokura, Trajectory of the anomalous Hall effect towards the quantized state in a ferromagnetic topological insulator, Nat. Phys. 10, 731 (2014).
  * Yasuda  _et al._ [2016] K. Yasuda, R. Wakatsuki, T. Morimoto, R. Yoshimi, A. Tsukazaki, K. S. Takahashi, M. Ezawa, M. Kawasaki, N. Nagaosa, and Y. Tokura, Geometric Hall effects in topological insulator heterostructures, Nat. Phys. 12, 555 (2016).
  * Chen  _et al._ [2019] J. Chen, L. Wang, M. Zhang, L. Zhou, R. Zhang, L. Jin, X. Wang, H. Qin, Y. Qiu, J. Mei, F. Ye, B. Xi, H. He, B. Li, and G. Wang, Evidence for magnetic skyrmions at the interface of ferromagnet/topological-insulator heterostructures, Nano Lett. 19, 6144 (2019).
  * Li  _et al._ [2021] P. Li, J. Ding, S. S. L. Zhang, J. Kally, T. Pillsbury, O. G. Heinonen, G. Rimal, C. Bi, A. DeMann, S. B. Field, W. Wang, J. Tang, J. S. Jiang, A. Hoffmann, N. Samarth, and M. Wu, Topological Hall effect in a topological insulator interfaced with a magnetic insulator, Nano Lett. 21, 84 (2021).
  * Wu  _et al._ [2020] H. Wu, F. Groß, B. Dai, D. Lujan, S. A. Razavi, P. Zhang, Y. Liu, K. Sobotkiewich, J. Förster, M. Weigand, G. Schütz, X. Li, J. Gräfe, and K. L. Wang, Ferrimagnetic skyrmions in topological insulator/ferrimagnet heterostructures, Adv. Mater. 32, 2003380 (2020).
  * Nomura and Nagaosa [2010] K. Nomura and N. Nagaosa, Electric charging of magnetic textures on the surface of a topological insulator, Phys. Rev. B 82, 161401 (2010).
  * Hurst  _et al._ [2015] H. M. Hurst, D. K. Efimkin, J. Zang, and V. Galitski, Charged skyrmions on the surface of a topological insulator, Phys. Rev. B 91, 060401 (2015).
  * Andrikopoulos  _et al._ [2016] D. Andrikopoulos, B. Sorée, and J. De Boeck, Skyrmion-induced bound states on the surface of three-dimensional topological insulators, J. Appl. Phys. 119, 193903 (2016).
  * Araki  _et al._ [2016] Y. Araki, A. Yoshida, and K. Nomura, Universal charge and current on magnetic domain walls in Weyl semimetals, Phys. Rev. B 94, 115312 (2016).
  * Sakai and Kohno [2014] A. Sakai and H. Kohno, Spin torques and charge transport on the surface of topological insulator, Phys. Rev. B 89, 165307 (2014).
  * Kurebayashi and Nagaosa [2019] D. Kurebayashi and N. Nagaosa, Theory of current-driven dynamics of spin textures on the surface of a topological insulator, Phys. Rev. B 100, 134407 (2019).
  * Heinrich [2021] B. Heinrich, Skyrmion birth at the notch, Nat. Nanotechnol. 16, 1051 (2021).
  * Iwasaki  _et al._ [2013] J. Iwasaki, M. Mochizuki, and N. Nagaosa, Current-induced skyrmion dynamics in constricted geometries, Nat. Nanotechnol. 8, 742 (2013).
  * Zhou and Ezawa [2014] Y. Zhou and M. Ezawa, A reversible conversion between a skyrmion and a domain-wall pair in a junction geometry, Nat. Commun. 5, 4652 (2014).
  * Everschor-Sitte  _et al._ [2017] K. Everschor-Sitte, M. Sitte, T. Valet, A. Abanov, and J. Sinova, Skyrmion production on demand by homogeneous dc currents, New J. Phys. 19, 092001 (2017).
  * Stier  _et al._ [2017] M. Stier, W. Häusler, T. Posske, G. Gurski, and M. Thorwart, Skyrmion–anti-skyrmion pair creation by in-plane currents, Phys. Rev. Lett. 118, 267203 (2017).
  * Büttner  _et al._ [2017] F. Büttner, I. Lemesh, M. Schneider, B. Pfau, C. M. Günther, P. Hessing, J. Geilhufe, L. Caretta, D. Engel, B. Krüger, J. Viefhaus, S. Eisebitt, and G. S. D. Beach, Field-free deterministic ultrafast creation of magnetic skyrmions by spin–orbit torques, Nat. Nanotechnol. 12, 1040 (2017).
  * Leonov and Inoue [2018] A. O. Leonov and K. Inoue, Homogeneous and heterogeneous nucleation of skyrmions in thin layers of cubic helimagnets, Phys. Rev. B 98, 054404 (2018).
  * Büttner  _et al._ [2021] F. Büttner, B. Pfau, M. Böttcher, M. Schneider, G. Mercurio, C. M. Günther, P. Hessing, C. Klose, A. Wittmann, K. Gerlinger, L.-M. Kern, C. Strüber, C. von Korff Schmising, J. Fuchs, D. Engel, A. Churikova, S. Huang, D. Suzuki, I. Lemesh, M. Huang, L. Caretta, D. Weder, J. H. Gaida, M. Möller, T. R. Harvey, S. Zayko, K. Bagschik, R. Carley, L. Mercadier, J. Schlappa, A. Yaroslavtsev, L. Le Guyarder, N. Gerasimova, A. Scherz, C. Deiter, R. Gort, D. Hickin, J. Zhu, M. Turcato, D. Lomidze, F. Erdinger, A. Castoldi, S. Maffessanti, M. Porro, A. Samartsev, J. Sinova, C. Ropers, J. H. Mentink, B. Dupé, G. S. D. Beach, and S. Eisebitt, Observation of fluctuation-mediated picosecond nucleation of a topological phase, Nat. Mater. 20, 30 (2021).
  * Tchoe and Han [2012] Y. Tchoe and J. H. Han, Skyrmion generation by current, Phys. Rev. B 85, 174416 (2012).
  * Sampaio  _et al._ [2013] J. Sampaio, V. Cros, S. Rohart, A. Thiaville, and A. Fert, Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures, Nat. Nanotechnol. 8, 839 (2013).
  * Hrabec  _et al._ [2017] A. Hrabec, J. Sampaio, M. Belmeguenai, I. Gross, R. Weil, S. M. Chérif, A. Stashkevich, V. Jacques, A. Thiaville, and S. Rohart, Current-induced skyrmion generation and dynamics in symmetric bilayers, Nat. Commun. 8, 15765 (2017).
  * Finizio  _et al._ [2019] S. Finizio, K. Zeissler, S. Wintz, S. Mayr, T. Weßels, A. J. Huxtable, G. Burnell, C. H. Marrows, and J. Raabe, Deterministic field-free skyrmion nucleation at a nanoengineered injector device, Nano Lett. 19, 7246 (2019).
  * Garanin [1997] D. A. Garanin, Fokker-Planck and Landau-Lifshitz-Bloch equations for classical ferromagnets, Phys. Rev. B 55, 3050 (1997).
  * Evans  _et al._ [2012] R. F. L. Evans, D. Hinzke, U. Atxitia, U. Nowak, R. W. Chantrell, and O. Chubykalo-Fesenko, Stochastic form of the Landau-Lifshitz-Bloch equation, Phys. Rev. B 85, 014433 (2012).
  * Note [1] Note that we have neglected the effect of demagnetizing field as it can be included as a small addition to an effective perpendicular magnetic anisotropy of a thin film.
  * Wakatsuki  _et al._ [2015] R. Wakatsuki, M. Ezawa, and N. Nagaosa, Domain wall of a ferromagnet on a three-dimensional topological insulator, Sci. Rep. 5, 13638 (2015).
  * Liu  _et al._ [2010] C.-X. Liu, X.-L. Qi, H. Zhang, X. Dai, Z. Fang, and S.-C. Zhang, Model Hamiltonian for topological insulators, Phys. Rev. B 82, 045122 (2010).
  * Wang  _et al._ [2015] W. Wang, F. Yang, C. Gao, J. Jia, G. D. Gu, and W. Wu, Visualizing ferromagnetic domains in magnetic topological insulators, APL Mater. 3, 083301 (2015).
  * Bogdanov and Hubert [1994] A. Bogdanov and A. Hubert, Thermodynamically stable magnetic vortex states in magnetic crystals, J. Magn. Magn. Mater. 138, 255 (1994).
  * Bessarab  _et al._ [2019] P. F. Bessarab, D. Yudin, D. R. Gulevich, P. Wadley, M. Titov, and O. A. Tretiakov, Stability and lifetime of antiferromagnetic skyrmions, Phys. Rev. B 99, 140411 (2019).
  * Tretiakov and Tchernyshyov [2007] O. A. Tretiakov and O. Tchernyshyov, Vortices in thin ferromagnetic films and the skyrmion number, Phys. Rev. B 75, 012408 (2007).
  * Potkina  _et al._ [2020] M. N. Potkina, I. S. Lobanov, O. A. Tretiakov, H. Jónsson, and V. M. Uzdin, Stability of long-lived antiskyrmions in the Mn-Pt-Sn tetragonal Heusler material, Phys. Rev. B 102, 134430 (2020).
  * [53] See Supplementary Material at [url will be inserted by publisher] for 1) the movie of the skyrmion nucleation by uniform current pulse and 2) the movie of the turbulent regime.
  * Note [2] The spin-transfer torque is derived as a response to an electric field $E$, instead of current $j$. In order to determine the energy efficiency, it is required to estimate the current density, which depends on sample quality and detailed material parameters. Therefore, a theoretical estimation of energy efficiency is not straightforward. However, in general, surface states of TI are known to be less affected by impurity scattering and possess higher mobility than those of conventional metals. Consequently, the energy loss from the Joule heating is expected to be smaller.
  * Wimmer  _et al._ [2021] S. Wimmer, J. Sánchez-Barriga, P. Küppers, A. Ney, E. Schierle, F. Freyse, O. Caha, J. Michalička, M. Liebmann, D. Primetzhofer, M. Hoffman, A. Ernst, M. M. Otrokov, G. Bihlmayer, E. Weschke, B. Lake, E. V. Chulkov, M. Morgenstern, G. Bauer, G. Springholz, and O. Rader, Mn-rich MnSb2Te4: A topological insulator with magnetic gap closing at high Curie temperatures of 45–50 K, Adv. Mater. 33, 2102935 (2021).
  * Liu  _et al._ [2020] J. Liu, A. Singh, Y. Y. F. Liu, A. Ionescu, B. Kuerbanjiang, C. H. W. Barnes, and T. Hesjedal, Exchange bias in magnetic topological insulator superlattices, Nano Lett. 20, 5315 (2020).
  * Fanchiang  _et al._ [2018] Y. T. Fanchiang, K. H. M. Chen, C. C. Tseng, C. C. Chen, C. K. Cheng, S. R. Yang, C. N. Wu, S. F. Lee, M. Hong, and J. Kwo, Strongly exchange-coupled and surface-state-modulated magnetization dynamics in Bi2Se3/yttrium iron garnet heterostructures, Nat. Commun. 9, 223 (2018).
  * Pereira  _et al._ [2020] V. M. Pereira, S. G. Altendorf, C. E. Liu, S. C. Liao, A. C. Komarek, M. Guo, H.-J. Lin, C. T. Chen, M. Hong, J. Kwo, L. H. Tjeng, and C. N. Wu, Topological insulator interfaced with ferromagnetic insulators: $\mathrm{B}{\mathrm{i}}_{2}\mathrm{T}{\mathrm{e}}_{3}$ thin films on magnetite and iron garnets, Phys. Rev. Mater. 4, 064202 (2020).
  * Tatara  _et al._ [2008] G. Tatara, H. Kohno, and J. Shibata, Microscopic approach to current-driven domain wall dynamics, Phys. Rep. 468, 213 (2008).
  * Kohno  _et al._ [2006] H. Kohno, G. Tatara, and J. Shibata, Microscopic calculation of spin torques in disordered ferromagnets, J. Phys. Soc. Jpn. 75, 113706 (2006).


# Simulation of angular resolution of a new electromagnetic sampling
calorimeter

Junlee Kim<EMAIL_ADDRESS>Eun-Joo Kim<EMAIL_ADDRESS>Young Jun Kim
Jung Keun Ahn  Gei Youb Lim  Division of Science Education, Jeonbuk National
University, Jeonju 54896, Korea  Department of Physics, Korea University,
Seoul 02841, Korea  Institute of Particle and Nuclear Studies (IPNS), High
Energy Accelerator Research Organization (KEK), Tsukuba 305-0801, Japan

###### Abstract

We report the simulation results for angular resolution of an electromagnetic
(EM) sampling calorimeter with photons in the range of 100 MeV to 2 GeV. The
simulation model of the EM calorimeter consists of alternating layers of an
1-mm-thick Pb plate and a 5-mm-thick plastic scintillator plate. The
scintillator plates are alternately segmented into horizontal and vertical
strips. In this study, we obtained energy deposits in individual strips using
Geant4 simulations and reconstructed the incident photon angles using XGBoost
with gradient-boosted decision trees. The performance of the angle
reconstruction depends on both the detector configuration and accuracy of
machine learning. The angular resolution so obtained can be expressed as
$0.24\oplus 1.24/\sqrt{E_{\gamma}}$, where $E_{\gamma}$ is the incident photon
energy in GeV. This energy dependence is consistent for different incident
angles in the range of 15∘ to 40∘.

###### keywords:

Electromagnetic Calorimeter , Geant4 , XGBoost

††journal: Nuclear Instruments and Methods in Physics Research A

##  1 Motivation

Electromagnetic (EM) calorimeter plays an important role in the experimental
studies of nuclear and particle physics [1, 2, 3]. Recently, sampling
calorimeters consisting of alternating layers of an absorber that generates EM
showers and an active medium that provides signals have become popular,
especially in large-scale high-energy experiments. In sampling calorimeters,
the energy deposited in the active medium fluctuates because the active layers
are interleaved with the absorber layers. This so-called sampling fluctuation
dominates the energy resolution. This sampling fluctuation can be improved by
optimizing the thickness ratio of the absorber to the active medium. A
3-m-long cylindrical calorimeter constructed of alternating layers of Pb and
plastic scintillator plates is an example of a sampling calorimeter [4].

In addition to the alternating layer structure, the segmented structure of the
active medium facilitates the measurement of the lateral distributions of the
EM shower along the incident photon direction. Energy deposits in the active
medium segmented in both the longitudinal and lateral directions allow for the
deduction of the incident direction of the incoming photon. A limitation in
the angular resolution, especially for low-energy photons, results from the
stochastic behavior of EM shower development [5]. We performed a simulation
study to deduce the incident angle of the photons, which can be a starting
point for designing a detector configuration to provide precise angle
measurements for photons.

This study describes a machine-learning approach based on XGBoost (XGB), which
is an optimized distributed gradient-boosting library [6]. With machine-
learning algorithms under the gradient-boosting framework, we can deduce the
direction of the incident photons entering a sampling calorimeter with
alternating layers of Pb plate and plastic scintillator strips in the
horizontal and vertical directions.

##  2 Electromagnetic shower simulation

Figure 1:  Schematic of the sampling calorimeter model consisting of 105
alternating layers of a Pb and segmented scintillator plate in 35 strips. Each
scintillator plate is oriented alternately in horizontal and vertical
directions. See text for details.

The sampling calorimeter was designed as a block consisting of alternating
layers of an 1-mm-thick Pb absorber and a 5-mm-thick polyvinyltoluene-based
plastic scintillator. The plastic scintillator is segmented into 15-mm-wide
strips, which are alternately oriented in the vertical and horizontal
directions, as shown in Figure 1. It has a cross-section of size 525 $\times$
525 mm2 and accommodates 105 alternating layers of 630 mm length, which
corresponds to 20 radiation lengths (20$X_{0}$) that are sufficiently long to
absorb full photon energy in the range of 0.1 to 2 GeV.

Figure 2:  Event display of simulated energy deposit patterns for an 1-GeV
photon entering the calorimeter ($\theta=$ 0) in (a) $xz$\- and (b)
$yz$-planes.

The detector response to the incident photons was simulated using Geant4 (ver.
4.10.06) with standard EM sub-packages [8]. The direction normal to the
detector surface defines the $z$-axis. The photon direction was defined by the
polar angle ($\theta$) with respect to the $z$-axis. Figure 2 illustrates the
simulated energy deposit patterns in each strip for an 1-GeV photon at normal
incidence in the $xz$\- and $yz$-planes. Each segmented region shown in Figure
2 represents a channel.

##  3 Reconstruction of incident angles

The incident angles of the photons are reconstructed using XGB [6], which maps
a feature dataset of the energy deposits in each scintillator strip onto a
target variable, the incident angle of a photon. Training data were carefully
prepared using the Geant4 simulation such that the input datasets were
representative of the detector response of the real sampling calorimeters. To
minimize the bias in the training process related to the dataset, the incident
angles were uniformly generated at the detector surface in the angular range
of 0 $<\theta<$ 50∘ and 0 $<\varphi<$ 360∘, where $\varphi$ denotes the
azimuthal angle. The number of training samples was $10^{5}$ considering the
limited computing resources. To test the reconstruction of the incident
angles, we generated photons at a fixed incident angle $\theta$ with a known
incident energy.

Figure 3:  The relative incident angle for 1-GeV photons generated at
$\theta=$ 10∘. The distribution is fitted with the Gaussian and generalized
Gaussian functions. The generalized Gaussian function gives a better result.

Figure 3 represents the distribution of the relative incident angle
($\Delta\theta$) for 1-GeV photons generated at $\theta=$ 10∘. $\Delta\theta$
represents the radial displacement of the reconstructed incident angle
relative to the true incident direction. The central 98% of the distribution
is fitted with the Gaussian and generalized Gaussian (GG) functions [7]. We
tested two functions, Gaussian and GG, to describe the reconstructed angular
distribution. The GG function, also known as the generalized error
distribution, is expressed as:

$\displaystyle
f(x;\mu,\alpha,\beta)=\frac{\beta}{2\alpha\Gamma(1/\beta)}e^{-(|x-\mu|/\alpha)^{\beta}},$
(1)

where $\mu$ is the mean value. The parameters $\alpha$ and $\beta$ determine
the scale and shape of the distribution, respectively. The variance in the GG
function is given by
$\sigma^{2}\equiv\alpha^{2}\Gamma(3/\beta)/\Gamma(1/\beta)$. The angular
resolution of the incident angle reconstruction is defined as $\sigma$.

Figure 4: Angular resolutions are displayed in terms of combination of
N_estimators and Max. depth. The best angular resolution is obtained with
N_estimators = 300 and Max. depth = 100.

Because the performance of the angle reconstruction of the XGB depends on the
hyperparameters, we scanned the evaluated angular resolution by changing the
hyperparameters. As there were correlations among the hyperparameters, the
scanning process was performed in a five-dimensional space for all possible
combinations. As an example, Figure 4 shows the test results for N_estimators
and Max. depth. N_estimators defines the maximum allowed number of decision
trees to be developed, and Max. depth defines the complexity of the decision-
tree structure. We explored a hyperparameter combination that generated the
best angular resolution. Consequently, the N_estimators and Max. depth were
set to 300 and 100, respectively. Similar scans were also performed for
different hyperparameters, subsample, learning rates, and gamma. Subsample
controls the fraction of total event samples for each boosting procedure, the
learning rate weights a decision tree to be added onto the current model, and
the gamma regulates the evaluation of each decision tree. Table 1 lists the
optimized set of hyperparameters that are used in further studies. During the
optimization process, the fraction of the tail that is not well described by
the GG function is approximately 2%.

Table 1: Hyperparameters of the XGB model Parameter | Function | Default value | Used value
---|---|---|---
N_estimators | The number of decision trees | N.A. | 300
Max. depth | Possible maximum depth of tree structure | 6 | 100
Subsample | Fraction of total data used for a single decision | 1 | 0.8
Learning rate | Step length for calculation | 0.3 | 0.02
Gamma | Requirement on minimum loss function | 0 | 0
Figure 5:  Angular resolutions are deduced in terms of the strip width for
1-GeV photons at $\theta=$ 10∘

Angular resolutions were deduced by varying the strip width from 5 to 35 mm,
as shown in Figure 5. The narrower the strip width, the larger is the size of
the features. Consequently, the angle-reconstruction performance is better
with 10-mm-wide strips than 5-mm-wide strips. However, the wider strips hardly
provide detailed information on the EM shower. The 15-mm-wide strips yield a
good angular resolution of 1.23 $\pm$ 0.01∘.

c0.cmb-0.4cm(a)  c0.cmb-0.4cm(b)

Figure 6:  (a) Reconstructed polar angles ($\theta$) for 1-GeV photons at
$\theta=$ 10∘. The red histogram represents the polar angle distribution for
$E_{\rm{tot}}>$ 10 MeV. (b) Angular resolution as a function of the radiation
length ($X_{0}$).

Focusing only on the angle reconstruction, we designed the detector with the
smallest number of layers. Figure 6 (a) shows the angles reconstructed using
the front 32 layers of the detector, which correspond to 6.2$X_{0}$, for 1-GeV
photons. A fraction of the photons fail to be reconstructed because of
insufficient information gathering caused by deeply penetrating photons
without shower generation. The failed events are represented as a delta
function near 0, and such events are removed by requiring the total energy
deposit to be higher than 10 MeV, which is 1% of the incident energy. The
angular resolution of the front layer is estimated to be 1.30 $\pm$ 0.01∘ when
1.9% of the photons fail to be reconstructed.

c0.cmb-0.4cm(a)  c0.cmb-0.4cm(b)

Figure 7:  (a) Angular resolution as a function of the incident angle for
different incident energies. (b) Angular resolution as a function of the
incident energy for $\theta=$ 25∘; it is fitted to the function $p_{0}\oplus
p_{1}/\sqrt{E_{\gamma}(\mathrm{GeV})}$

Figure 7 (a) shows the angular resolution as a function of incident angle for
different incident energies ($E_{\gamma}$). The angular resolution does not
depend on the incident angle for high incident energies. However, the angular
resolution changes significantly for low incident energies at small incident
angles. Figure 7 (b) shows the angular resolution as a function of the
incident photon energy at $\theta=$ 25∘. The angular resolution fits to the
function $p_{0}\oplus p_{1}/\sqrt{E_{\gamma}(\mathrm{GeV})}$, where $p_{0}$
and $p_{1}$ represent the energy-independent and energy-dependent
contributions, respectively.

Figure 8:  $p_{0}$ and $p_{1}$ as a function of the incident angle

Figure 8 shows the estimated $p_{0}$ and $p_{1}$ as a function of the incident
angle. The average values of $p_{0}$ and $p_{1}$ are estimated to be 0.24∘ and
1.24∘, respectively, for $\theta>$ 10∘. For a small $\theta$, $p_{1}$ is
smaller than that for a larger $\theta$ because the angular resolution is
largely dependent on the incident angle for low incident energies, whereas the
opposite trend can be observed for $p_{0}$.

Figure 9:  The angular resolution for 1-GeV photons at $\theta=$ 10∘ as a
function of the number of training samples with different strip widths
obtained using the first 32 layers or the full layers.

Figure 9 shows the angular resolution for different numbers of training
samples and different detector configurations. All configurations show
improvement in the angular resolution with increasing training samples, which
indicates that the maturity of machine learning is an important factor in
angle reconstruction. With an increase in the number of training events, the
Angular resolution improves most rapidly for the configuration with 5-mm-wide
scintillator strips. In other words, the required training events must be
sufficiently large to correlate a larger number of features. Considering the
limited computing resources, the configuration of 32 layers and 15-mm-wide
strip with $10^{5}$ training events was chosen to train the XGB model. The
training result for such a configuration differs by only 0.2∘ from that of the
setup with full layers, 5-mm-wide strips, and $5\times 10^{5}$ training
events.

##  4 Summary

We report the simulation results of the incident angle reconstruction of EM
sampling calorimeters with photons in the range of 0.1 to 2 GeV. To improve
the angular resolution, we tuned the hyperparameters of XGBoost while
considering the correlating effects among them.

In terms of detector configuration, strips with narrower widths provide better
angular resolution using a finely measured EM shower with the assumption that
the training is fully conducted. However, the configuration with a smaller
strip width causes a larger feature size, which negatively affects the
training depending on the quantity of training data. The angular resolution
may be optimized by considering the two competing effects of the strip width
and quantity of training data. We dedicated the front layers of the detector
to reconstructing the incident angle. We obtain an angular resolution of 1.30
$\pm$ 0.01∘ for the configuration of 32 layers (6.2$X_{0}$) of 15-mm-wide
strips, which is a reduction of 0.1∘ compared to that of the configuration
with full layers (20$X_{0}$) in the case of $10^{5}$ training samples.

The angular resolution depends on the incident energy of the photon, and the
dependence is expressed as $p_{0}\oplus p_{1}/\sqrt{E_{\gamma}}$. Except for
small incident angles ($\theta<$ ∘ deg), the two parameters do not
significantly change up to $40^{\circ}$, and the dependence is estimated to be
0.24 $\oplus$ 1.24${}^{\circ}/\sqrt{E_{\gamma}}$.

## Acknowledgment

This work was supported by the National Research Foundation of Korea (NRF)
(Grants No. 2019R1A2C1084552, 2022R1A5A1030700, and 2020R1A3B2079993) and JSPS
KAKENHI Grant No. JP21H01118.

## References

  * [1] Y. Tajima et al., Nucl. Instrum. Meth. A 592, 261 (2008).
  * [2] Giovanni Mazza et al. (CMS Collaboration), Nucl. Instrum. Meth. A 958 (2020)
  * [3] K. Miyabayashi et al., JINST 9 (2014) 09, P09011
  * [4] R. Murayama et al., Nucl. Instrum. Meth. A 953 (2020)
  * [5] O. Adriani et al., Nucl. Instrum. Meth. A 854 (2017)
  * [6] C. Tianqi et al., Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, (2016) 785.
  * [7] D. Gonzalez et al., IEEE Trans. Image Processing. 4, 485 (2007)
  * [8] S. Agostinelli et al., Nucl. Instr. and Meth. A 506 (2003) 250.