[7]G^ #1,#2_#3,#4(#5 #6| #7)

# Noncollinear twisted RKKY interaction on the optically driven SnTe(001)
surface

Mohsen Yarmohammadi<EMAIL_ADDRESS>Department of Physics,
The University of Texas at Dallas, Richardson, Texas 75080, USA  Marin Bukov
<EMAIL_ADDRESS>Max Planck Institute for the Physics of Complex
Systems, Nöthnitzer Str. 38, 01187 Dresden, Germany  Faculty of Physics, St.
Kliment Ohridski University of Sofia, 5 James Bourchier Blvd, 1164 Sofia,
Bulgaria  Michael H. Kolodrubetz<EMAIL_ADDRESS>Department of Physics,
The University of Texas at Dallas, Richardson, Texas 75080, USA

###### Abstract

The nontrivial spin texture on the (001) surface of topological crystalline
insulator SnTe hosts exotic scientific importance and spintronic applications.
Here, we study the effects of weak Floquet optical driving on the Ruderman-
Kittel-Kasuya-Yosida (RKKY) interaction between two magnetic impurities on a
doped SnTe(001) surface. Due to peculiar spin-orbit hybridization, we find a
noncollinear twisted RKKY interaction comprising XYZ-Heisenberg, symmetric in-
plane, and asymmetric Dzyaloshinskii-Moriya (DM) terms. We see that
contributions from the $z$ ($x$)-component of the XYZ-Heisenberg (DM)
interaction are dominant for most parameters. The interactions, including DM
terms that are responsible for interesting spin textures, require doping in
most cases. We propose to modify the interactions in situ via optical control
of band structure, and thereby doping. A notable aspect of this control
protocol is breaking of electron-hole symmetry, which stems from the DM
interaction.

##  I Introduction

For the past decade, topological insulators (TIs) have been at the forefront
of condensed matter physics [1, 2]. Due to the strong spin-orbit coupling
(SOC) [3, 4], topological phases of matter contain gapped bulk and gapless
surface states, protected by the time-reversal symmetry. Notably, quantum
anomalous Hall effect (QAHE) has been observed in magnetically doped TIs [5,
6, 7, 8], which is a topological effect with potential for low-power
information processing associated with dissipationless chiral edge transport
[9]. However, the QAHE becomes weak in a wide variety of transition metal-
doped TIs [10, 11, 12, 13, 7] if one increases the density of surface
carriers, mainly due to weakening of the magnetic order as a prerequisite for
anomalous Hall physics.

Although the mechanism of this phenomenon is still under discussion [14], it
was suggested that the indirect Ruderman-Kittel-Kasuya-Yosida (RKKY) exchange
coupling between magnetic ions/impurities plays a significant role [15, 16]
among other common types of magnetism because it is tightly connected to the
density of itinerant electrons in the host material [17, *PhysRev.106.893,
*PhysRev.96.99]. Depending on the spin-orbit structures of topological
materials, various types of RKKY mechanisms [20, 21, 22] and eventually
different weak/strong regimes of QAHE can emerge. Furthermore, emergent
spintronic applications need the RKKY mechanism to find the proper material
with improved adjustability of spin alignments for logic magnetic devices [23,
24]. It is the latter purpose that we mainly aim at in this work, and we only
conceptually (i.e., without additional calculations) state its relationship
with the QAHE.

It has been shown that nontrivial RKKY interactions can appear on the surface
of three-dimensional TIs [16, 25, 26, 27] as well as at the edge of two-
dimensional TIs [28, 29, 30]. In this work, we proceed with the topological
crystalline insulators (TCIs) [31, 32, 33, 34, 35] because, in addition to the
time-reversal symmetry, crystal symmetries protect gapless surface states as
well. This, in turn, covers a broad range of materials rather than specific TI
compounds. Moreover, in contrast to an odd number of Dirac cones in strong
TIs, TCIs provide an even number of cones, useful for spintronics and
valleytronics [23, 24, 36]. We choose to focus on SnTe(001) surface as one of
the well-known TCIs with a nontrivial spin texture [37] to see how its
magnetic features can go beyond the ones in TIs to be reflected in spintronics
and QAHE physics.

In order to understand novel features in these materials, we need the ability
to tune the properties of the surface. The best way to do so is through a gap
opening at Dirac cones because it suppresses the large density of the doped
surface around the Fermi energy, first, and second, it makes the magnetic
order strong enough through the suppression of various phase alterations of
ferromagnetic (FM) and antiferromagnetic (AFM) interactions. Among various
ways for gap opening [38], we proceed with optical driving using Floquet-Bloch
states [39, 40, 41, 42] because by tuning the intensity and frequency of the
radiation, one can dynamically tune the photon-dressed bands, and eventually
modulate the properties of an irradiated system in a controlled way.

Motivated by various light-induced phenomena [43, 44, 45, 46, 47, 48], we
consider the interaction of a weak (to avoid the heating issue) circularly
polarized light with the SnTe(001) surface, which opens a gap isotropically at
all four Dirac cones, and thereby modulates the RKKY interaction.
Subsequently, assuming that particles fill the Floquet-Bloch bands up to a
given chemical potential when the weak Floquet drive is turned on, we pose the
question: how does the interplay between doping and optical driving influence
the RKKY interaction for both spintronics and QAHE treatment? To answer this
question, we reformulate the Heisenberg RKKY interaction [49, 50] of hybrid
states of sublattices in SnTe(001) for the bare spin-orbital states of each
sublattice to take into account the effect of a nontrivial spin texture. This
results in noncollinear twisted spin alignments including XYZ-Heisenberg,
symmetric in-plane, and asymmetric DM interactions, implying that magnetic
moments lie on the surface with increased adjustability. The asymmetric DM
interaction is substantial for spintronic applications [51, 52, 53].
Interestingly, the Floquet drive tends to treat the QAHE by suppressing the
density of doped states and RKKY amplitudes; however, we note that topological
features of the SnTe(001) surface are not essential for our study.

This paper is structured as follows. In Sec. II, we present the continuum
Hamiltonian of pristine SnTe(001) surface and construct the main building
block of optically driven states. In Sec. III, we turn to the RKKY theory and
derive the Floquet-Green’s functions for a weak and off-resonant drive. The
numerical results are presented in Sec. IV. We discuss potential extensions in
Sec. V and finally we end the paper with a conclusion in Sec. VI.

##  II Hamiltonian model

###  II.1 Pristine $\rm SnTe$(001)

In this section, we intorduce the continuum Hamiltonian model [33, 35, 34, 32,
54, 38] to describe each coaxial Dirac cone in the low-energy limit of
SnTe(001) surface. Two cones $\Lambda_{x}$ and $\Lambda^{\prime}_{x}$ are
located along the direction $X_{1}-\Gamma-X_{1}$ of the projected surface
Brillouin zone (SBZ), while two cones $\Lambda_{y}$ and $\Lambda^{\prime}_{y}$
are located along the direction $X_{2}-\Gamma-X_{2}$, as shown in the left
side of Fig. 1. The pristine Hamiltonian of the $\Lambda_{x}$ point is given
by ($\hbar=1$) [38, 55, 49, 54, 50, 56]

$\mathcal{H}_{\Lambda_{x}}=\tilde{\eta}_{x}p_{x}\sigma_{y}-\tilde{\eta}_{y}p_{y}\sigma_{x}\,,$
(1)

where $p_{x}=k_{x}-\Lambda_{x}$ is measured from the Dirac cone at
$\Lambda_{x}=(\sqrt{n^{2}+\delta^{2}}/\eta_{x},0)$, $p_{y}=k_{y}$,
$\tilde{\eta}_{x}=(\delta/\sqrt{n^{2}+\delta^{2}})\eta_{x}$, and
$\tilde{\eta}_{y}=\eta_{y}$ [55]. Here, $\eta_{x}=3.53$ eV.Å and
$\eta_{y}=1.91$ eV.Å refer to the bare Fermi velocities along the $x$\- and
$y$-direction, respectively. The Pauli matrices are given by $\sigma_{x}$ and
$\sigma_{y}$. The intervalley scattering parameters $n=55$ meV and $\delta=40$
meV are, respectively, responsible for the high-energy Dirac cones and the
formation of two copies of Dirac cones with opposite chiralities in the four-
band model (not shown here) [32, 38, 54, 49]. Although they correctly describe
intervalley scattering at the lattice scale of SnTe(001), they vary sample to
sample since they arise from non-universal effects such as surface roughness.
Indeed, the parameter $n$ is the threshold energy below which the low-energy
model is valid.

Figure 1: The low-energy spectrum of Dirac fermions on the SBZ of SnTe(001)
surface, left side, and the optically driven spectrum, right side, with
induced gap $\Delta=2a^{2}/\Omega$ from the light with the ac field $a$ and
the frequency $\Omega$. Other Dirac cones follow from $C_{2}$ and $C_{4}$
rotational symmetries.

The rotational symmetry $C_{2}$ protects the coaxial Dirac cones resulting in
$\mathcal{H}_{{\Lambda^{\prime}_{x,y}}}=\mathcal{H}_{{\Lambda_{x,y}}}$, while
the perpendicular cones are formed by the $C_{4}$ symmetry. This, in turn,
implies that
$\mathcal{H}_{\Lambda_{y}}={}\tilde{\eta}_{y}p_{x}\sigma_{y}-\tilde{\eta}_{x}p_{y}\sigma_{x}$.
It should be stressed that the valence and conduction bands of the energy
spectrum refer to the hybridized $p$-orbitals of the cation Sn2+ and anion
Te2- sublattices, namely $|1\rangle=(|\uparrow,\rm Te\rangle+|\downarrow,\rm
Sn\rangle)/\sqrt{2}$ and $|2\rangle=(|\uparrow,\rm Sn\rangle+|\downarrow,\rm
Te\rangle)/\sqrt{2}$ [37, 38, 55, 54].

Before turning to the irradiated SnTe(001) surface, we make one useful
simplification, which will help later to find analytical RKKY interactions.
Because RKKY theory is mainly based on low-momenta response, we neglect the
anisotropy feature of the surface Dirac cones originating from
$\tilde{\eta}_{x}\simeq 2.08$ eV.Å and $\tilde{\eta}_{y}=1.91$ eV.Å. Hence, we
set $v_{\rm F}=(\tilde{\eta}_{x}+\tilde{\eta}_{y})/2\simeq 2$ eV.Å.

###  II.2 Irradiated $\rm SnTe$(001)

In this section, we present the expression for the effective two-band
Hamiltonian model discussed above for the case when the SnTe(001) Dirac cones
are weakly driven by a circularly polarized light. To do so, the time-periodic
vector potential $\vec{A}(t)=A_{0}[\sin(\Omega t),\cos(\Omega t)]$, where
$A_{0}=E_{0}/\sqrt{2}\Omega$ ($E_{0}$ and $\Omega$ are the amplitude and
frequency of light, respectively) is chosen to drive the Hamiltonian of the
$\Lambda_{x}$ point (and other Dirac cones) through the minimal coupling
scheme as

$\mathcal{H}_{\Lambda_{x}}(t)={}\begin{pmatrix}0&-{\tt i}v_{\rm
F}p_{-}-ae^{+{\tt i}\Omega t}\\\ \\\ +{\tt i}v_{\rm F}p_{+}-ae^{-{\tt i}\Omega
t}&0\end{pmatrix},$ (2)

where $p_{\pm}=p_{x}\pm{\tt i}p_{y}$ and $a=eA_{0}v_{\rm F}$. The drive
induces transitions between the eigenstates. Since the Hamiltonian is periodic
in time, $\mathcal{H}_{\Lambda_{x}}(t)=\mathcal{H}_{\Lambda_{x}}(t+T)$ with
$T=2\pi/\Omega$, the valid theory to address the time-periodic Hamiltonians is
the Floquet-Bloch’s theorem [57, 42]: a complete set of solutions of the time-
dependent Schrödinger equation
$\mathcal{H}_{\Lambda_{x}}(t)|\psi(t)\rangle={\tt
i}\partial_{t}|\psi(t)\rangle$ reads as $|\psi_{f,p}(t)\rangle={}e^{-{\tt
i}\varepsilon_{f,p}}|\mathcal{F}_{f,p}(t)\rangle$, where the integer Floquet
index $f$ classifies different sidebands with energies $\varepsilon_{f,p}$.
Accordingly, the Fourier series
$|\mathcal{F}_{f,p}(t)\rangle=\sum^{\infty}_{j=-\infty}e^{-{\tt i}j\Omega
t}|\mathcal{F}^{j}_{f,p}\rangle$ describes the time-periodic Floquet states
$|\mathcal{F}_{f,p}(t)\rangle=|\mathcal{F}_{f,p}(t+2\pi/\Omega)\rangle$. It
should be mentioned that $f=0$ maps the quasienergies to the first Floquet
zone $[-\Omega/2,\Omega/2]$. Then, it is straightforward to deduce the Floquet
Hamiltonian $\mathcal{H}_{{\rm
F},\Lambda_{x}}=\mathcal{H}_{\Lambda_{x}}(t)-{\tt i}\partial_{t}$ in the
Fourier domain [58] $\sum^{\infty}_{j=-\infty}\mathcal{H}_{{\rm
F},mj}|\mathcal{F}^{j}_{f,p}\rangle=\varepsilon_{f,p}|\mathcal{F}^{m}_{f,p}\rangle$.
Having the driving period $T=2\pi/\Omega$, the matrix elements of the Floquet
Hamiltonian can be found via $\mathcal{H}_{{\rm
F},mj}={}\mathcal{H}_{mj}-j\Omega\delta_{m,j}$, where $j$ is the photon number
and [59, 60, 61]

$\mathcal{H}_{mj}={}\frac{1}{T}\int^{T}_{0}\mathcal{H}_{\Lambda_{x}}(t)e^{{\tt
i}(m-j)\Omega t}dt\,.$ (3)

The Floquet Hamiltonian can be numerically diagonalized to study the exact
dynamics governed by the application of light. However, for sufficiently weak
electric field intensities, i.e., when the field intensity $a$ is small
compared to the frequency $\Omega$, we may approximate the Floquet Hamiltonian
using a high-frequency expansion explained below. In the high-frequency limit,
the Floquet sidebands do not overlap; thus, there are no gaps due to band
crossings so the hybridization between different Floquet bands becomes weak,
meaning that interband transitions of electrons are suppressed.

Thus, using a perturbative approach known as the van Vleck inverse-frequency
expansion, we obtain the effective Floquet Hamiltonian in powers of
$\Omega^{-1}$, $\mathcal{H}_{{\rm
F},\Lambda_{x}}\simeq{}\mathcal{H}_{\Lambda_{x}}+\frac{\left[\mathcal{H}_{-1},\mathcal{H}_{+1}\right]}{\Omega}$
[41, 62, 41, 63], which reads as:

$\displaystyle\mathcal{H}_{{\rm F},\Lambda_{x}}={}$
$\displaystyle\begin{pmatrix}+\Delta/2&-{\tt i}v_{\rm F}\,p\,e^{-{\tt
i}\varphi_{p}}\\\ \\\ +{\tt i}v_{\rm F}\,p\,e^{+{\tt
i}\varphi_{p}}&-\Delta/2\end{pmatrix}\,,$ (4)

where a gap $\Delta={}2a^{2}/\Omega$ opens up at the Dirac cone through a two-
photon process similar to that in Refs. [39, 64, 65, 59, 61]. Using the
$C_{2}$ and $C_{4}$ symmetries, one can obtain the Floquet Hamiltonian at
other Dirac cones. Diagonalizing the above Hamiltonian leads to a gapped Dirac
dispersion $\mathcal{E}_{\Lambda_{x}}(p)={}\pm\,\sqrt{v^{2}_{\rm
F}p^{2}+\Delta^{2}/4}$, as shown in the right panel of Fig. 1.

In addition to optical driving, there are various ways to open the gap at the
Dirac cones [38], namely the exchange magnetization, the Zeeman term without
an external magnetic field, and proximity coupling to a ferromagnet. However,
as discussed in the introduction, the Floquet drive allows us to tune the
properties of the dispersion for various purposes by adjusting the intensity
and frequency of light. Moreover, further photon processes in
$\mathcal{H}_{{\rm F},\Lambda_{x}}$ can be taken into account by including
higher-order terms [59, 59, 61].

For our purposes, the Floquet drive will be used to tune the RKKY interaction
for both spintronics and eventually QAHE physics (indirectly through the
density of surface carriers). We now proceed to the RKKY interaction in
irradiated SnTe(001).

##  III RKKY interaction in irradiated $\rm SnTe$(001)

As explained in the introduction, the indirect exchange coupling between two
magnetic impurities or two localized spins $\vec{S}_{1}$ and $\vec{S}_{2}$
mediated by the host itinerant electron spins $\vec{s}$ is described by RKKY
theory [17, *PhysRev.106.893, *PhysRev.96.99]. In our case, the host itinerant
electrons are the optically driven ones on the SnTe(001) surface given by Eq.
(4). According to this theory, the magnetic moments are treated as immobile
defects on the lattice sites $\vec{R}_{1}$ and $\vec{R}_{2}$, and the
interaction Hamiltonian is given by $\mathcal{H}_{\rm
int}=J\sum_{i}\vec{S}_{i}\cdot\vec{s}_{i}$, where $J$ denotes the bare
exchange energy. Thus, we consider the total Hamiltonian for SnTe(001) surface
including a contact interaction with the magnetic centers, namely,
$\mathcal{H}=\mathcal{H}_{{\rm F},\Lambda_{x}}+\mathcal{H}_{\rm int}$. Using
second-order perturbation theory, one finds [17, *PhysRev.106.893,
*PhysRev.96.99]

$\mathcal{H}^{\alpha\beta}_{\rm
RKKY}={}J^{2}\sum_{l,j}\,S^{l\alpha}_{1}\,\chi^{\alpha\beta}_{lj}(\,\vec{R}\,)\,S^{j\beta}_{2}\,,$
(5)

with lattice sites $\alpha$ and $\beta$, and the impurity separation
$\vec{R}=\vec{R}_{2}-\vec{R}_{1}$. In the above equation, the spin
susceptibility $\chi^{\alpha\beta}_{lj}(\,\vec{R}\,)$ can be obtained from the
retarded Green’s functions in the spin-space for different flavors
$\\{l,j\\}\in\\{x,y,z\\}$ [66, 67, 68]:

$\chi^{\alpha\beta}_{lj}(\,\vec{R}\,)=-\frac{2}{\pi}\,\Im\int^{\mathcal{E}_{\mathrm{F}}}_{-\infty}d\mathcal{E}\,\mathbb{F}^{\alpha\beta}_{lj}(\mathcal{E},\vec{R})\,,$
(6)

where $\mathcal{E}_{\rm F}$ is the Fermi energy and

$\mathbb{F}^{\alpha\beta}_{lj}(\mathcal{E},\vec{R})=\mathrm{Tr}\left[\sigma_{l}G^{\alpha\beta}(\mathcal{E},\vec{R})\sigma_{j}G^{\beta\alpha}(\mathcal{E},-\vec{R})\right]\,.$
(7)

We set the temperature to zero in the present work, so $\mathcal{E}_{\rm
F}=\mu$, where $\mu$ is the chemical potential. The role of doping is then
characterized by $\mu>0$ and $\mu<0$ referring to electron and hole doping,
respectively. We assume that the electrons fill all energies up to $\mu$
within the Floquet-modified dispersion because the optical drive considered is
sufficiently weak and off-resonant that it may be turned on continuously,
resulting in a new (quasi zero temperature) steady state. Note that this is
fundamentally different from many other Floquet systems which are
experimentally constrained to pump-probe non-equilibrium states to avoid
heating, for which the electron filling is quite different [69]. In Eq. (7),
we have

$G^{\alpha\beta}(\mathcal{E},\vec{R})=\begin{pmatrix}G^{\uparrow\uparrow}_{\alpha\beta}\,(\mathcal{E},\vec{R})&&G^{\uparrow\downarrow}_{\alpha\beta}\,(\mathcal{E},\vec{R})\\\
\\\
G^{\downarrow\uparrow}_{\alpha\beta}\,(\mathcal{E},\vec{R})&&G^{\downarrow\downarrow}_{\alpha\beta}\,(\mathcal{E},\vec{R})\end{pmatrix}\,,$
(8)

which denotes the real-space non-interacting Green’s function of sublattices
in spin-space, given by the photo-dressed Dirac bands. However, we must sum
over the directions $X_{1}-\Gamma-X_{1}$ and $X_{2}-\Gamma-X_{2}$ on the
SnTe(001) surface, as well as the pair of Dirac cones at $X_{1}$ and $X_{2}$
points, to cover the entire SBZ. Thus, we have

$\displaystyle G^{\sigma\sigma^{\prime}}_{\alpha\beta}(\mathcal{E},\vec{R})=$
$\displaystyle\frac{1}{S_{\mathrm{SBZ}}}\int
d^{2}p\,e^{\texttt{i}\vec{p}\cdot\vec{R}}\Big{[}e^{\texttt{i}\vec{X}_{1}\cdot\vec{R}}G^{\sigma\sigma^{\prime}}_{\alpha\beta}(\vec{p}+\vec{X}_{1},\mathcal{E})$
(9)
$\displaystyle+e^{\texttt{i}\vec{X}_{2}\cdot\vec{R}}G^{\sigma\sigma^{\prime}}_{\alpha\beta}(\vec{p}+\vec{X}_{2},\mathcal{E})\Big{]}\,,$

with $\\{\sigma,\sigma^{\prime}\\}=\\{\uparrow,\downarrow\\}$, where
$S_{\mathrm{SBZ}}$ is the area of the entire SBZ. The integral is dominated by
low momentum, $p\ll|\Lambda_{x,y}|$, justifying our choice to treat each Dirac
cone independently in the low-energy limit of two-band model.

It is necessary to point out that in the Hamiltonian model, a superposition of
sublattices given by the eigenstates $|1\rangle$ and $|2\rangle$ appears,
rather than the bare sublattices $|{\rm Te/Se},\uparrow/\downarrow\rangle$.
Therefore, one needs to express the real-space Green’s functions of the hybrid
states in terms of the bare sublattices considering the corresponding spins.
Hence, we rewrite the above Green’s function as [see Appendix A]

$G^{\sigma\sigma^{\prime}}_{\alpha\beta}(\mathcal{E},\vec{R})=\sum_{\gamma,\eta}u^{*}_{\alpha,\gamma,\sigma}u_{\beta,\gamma,\sigma^{\prime}}\mathcal{G}_{\gamma\eta}(\mathcal{E},\vec{R})\,,$
(10)

in which the summation runs over hybrid states $\\{|1\rangle,|2\rangle\\}$ and
$u_{\rm Te,1,\uparrow}=u_{\rm Te,2,\downarrow}=u_{\rm Sn,1,\downarrow}=u_{\rm
Sn,2,\uparrow}=1/\sqrt{2}$. To calculate
$\mathcal{G}_{\gamma\eta}(\mathcal{E},\vec{R})$, we first calculate the
reciprocal-space Green’s functions around all Dirac cones within the weak-
driven approximation. Therefore, using the relation
$\mathcal{G}^{\Lambda_{x}}(p,\mathcal{E})={}\left[\mathcal{E}+{\tt
i}o^{+}-\mathcal{H}_{{\rm F},\Lambda_{x}}\right]^{-1}$ with infinitesimal
parameter $o^{+}\to 0$ for the $\Lambda_{x}$ point as well as using the
expressions $d^{2}p=p\,dp\,d\varphi_{p}$ and
$\exp\big{[}\texttt{i}\vec{p}\cdot\vec{R}\big{]}=\exp\left[\texttt{i}\,p\,R\,\cos\left(\varphi_{p}-\varphi_{R}\right)\right]$
with $\varphi_{p}=\tan^{-1}\left(p_{y}/p_{x}\right)$, we obtain

$\displaystyle\mathcal{G}_{\gamma\eta}(\mathcal{E},\vec{R})={}\frac{1}{S_{\mathrm{SBZ}}}\int^{p_{\rm
c}}_{0}\int^{2\pi}_{0}d^{2}p{}e^{\texttt{i}\,\vec{p}\cdot\vec{R}}\,\mathcal{G}_{\gamma\eta}(p,\mathcal{E}).$
(11)

We set the cutoff $p_{\rm c}\to\infty$, as the long-range interactions are
dominated by low momentum modes within this integral. Hence, for the
$\Lambda_{x}$ point, the real-space Green’s functions of hybrid states
$|1\rangle$ and $|2\rangle$ read as

$\displaystyle\mathcal{G}_{11}(\mathcal{E},\vec{R})={}$
$\displaystyle\frac{-2\pi\left[\mathcal{E}+{\tt
i}o^{+}+\Delta/2\right]}{S_{\mathrm{SBZ}}\,v^{2}_{\rm F}}\mathbb{K}_{0}(-{\tt
i}\widetilde{\mathcal{E}}\,R/v_{\rm F})\,,$ (12a)
$\displaystyle\mathcal{G}_{12}(\mathcal{E},\vec{R})={}$
$\displaystyle\frac{+2\pi\,{\tt i}\,e^{-{\tt
i}\varphi_{R}}\widetilde{\mathcal{E}}}{S_{\mathrm{SBZ}}\,v^{2}_{\rm
F}}\mathbb{K}_{1}(-{\tt i}\widetilde{\mathcal{E}}\,R/v_{\rm F})\,,$ (12b)
$\displaystyle\mathcal{G}_{21}(\mathcal{E},\vec{R})={}$
$\displaystyle\frac{-2\pi\,{\tt i}\,e^{+{\tt
i}\varphi_{R}}\widetilde{\mathcal{E}}}{S_{\mathrm{SBZ}}\,v^{2}_{\rm
F}}\mathbb{K}_{1}(-{\tt i}\widetilde{\mathcal{E}}\,R/v_{\rm F})\,,$ (12c)
$\displaystyle\mathcal{G}_{22}(\mathcal{E},\vec{R})={}$
$\displaystyle\frac{-2\pi\left[\mathcal{E}+{\tt
i}o^{+}-\Delta/2\right]}{S_{\mathrm{SBZ}}\,v^{2}_{\rm F}}\mathbb{K}_{0}(-{\tt
i}\widetilde{\mathcal{E}}\,R/v_{\rm F})\,,$ (12d)

where $\widetilde{\mathcal{E}}=\sqrt{(\mathcal{E}+{\tt
i}o^{+})^{2}-\Delta^{2}/4}$ and $\mathbb{K}_{0,1}$ denote modified Bessel
functions. For the $\Lambda^{\prime}_{x}$, $\Lambda_{y}$, and
$\Lambda^{\prime}_{y}$ points, we have the exact same expressions since the
gap is opened isotropically at all Dirac cones by the optical driving. The
above elements fulfil the relations
$\mathcal{G}_{11/22}(\mathcal{E},-\vec{R})=\mathcal{G}_{11/22}(\mathcal{E},\vec{R})$
and
$\mathcal{G}_{12/21}(\mathcal{E},-\vec{R})=-\mathcal{G}_{12/21}(\mathcal{E},\vec{R})$.

Turning back to Eq. (5), one would rewrite the RKKY Hamiltonian as

$\mathcal{H}_{\rm
RKKY}={}-\frac{2J^{2}}{\pi}\text{Im}\int_{-\infty}^{\mu}d\mathcal{E}\sum_{lj}S^{l}_{1}S^{j}_{2}\,\mathbb{F}_{lj}\,.$
(13)

The components $\mathbb{F}^{\rm\alpha\,\alpha}_{lj}$ and
$\mathbb{F}^{\rm\alpha\,\beta}_{lj}$ for the spins on the same and different
sublattices, respectively, are presented in Appendix B. Accordingly, the
Hamiltonians read as

$\displaystyle\mathcal{H}^{\alpha\,\alpha,\alpha_{\rm s}}_{\rm RKKY}=$
$\displaystyle\sum_{i}\mathcal{J}^{\alpha\,\alpha}_{i}S^{i}_{1}S^{i}_{2}+\vec{\mathcal{J}}_{\rm
DM}^{\alpha_{\rm s}}\cdot(\vec{S}_{1}\times\vec{S}_{2})$
$\displaystyle+\alpha_{\rm
s}\mathcal{J}^{\alpha\,\alpha}_{xy}[S^{x}_{1}S^{y}_{2}+S^{y}_{1}S^{x}_{2}]\,,$
(14a) $\displaystyle\mathcal{H}^{\alpha\,\beta,\alpha_{\rm d}}_{\rm RKKY}=$
$\displaystyle\sum_{i}\mathcal{J}^{\alpha\,\beta}_{i}S^{i}_{1}S^{i}_{2}+\vec{\mathcal{J}}_{\rm
DM}^{\alpha_{\rm d}}\cdot(\vec{S}_{1}\times\vec{S}_{2}),$ (14b)

where we have $\vec{\mathcal{J}}^{\alpha_{\rm s}}_{\rm
DM}=(\mathcal{J}^{\alpha\,\alpha}_{{\rm DM},x},\alpha_{\rm
s}\mathcal{J}^{\alpha\,\alpha}_{{\rm DM},y},0)$, and
$\vec{\mathcal{J}}^{\alpha_{\rm d}}_{\rm
DM}=(\mathcal{J}^{\alpha\,\beta}_{{\rm DM},x},0,\alpha_{\rm
d}\mathcal{J}^{\alpha\,\beta}_{{\rm DM},z})$. When the impurities are placed
on the TeTe (SnSn) sublattices, we use $\alpha_{\rm s}=+1~{}(-1)$, while
$\alpha_{\rm d}=+1~{}(-1)$ for the TeSn (SnTe) sublattices.

The first terms in both Hamiltonians of Eq. (14) couple the same spin
directions with different exchange strengths {$\mathcal{J}_{x}$,
$\mathcal{J}_{y}$, $\mathcal{J}_{z}$}, which gives the XYZ-Heisenberg
interaction. Furthermore, due to the intrinsic SOC of the SnTe(001) surface,
the symmetry of spin space is expected to be broken in response to the
impurities, leading to off-diagonal components of the RKKY Hamiltonian (the
last terms) [70, 71]. While the asymmetric second terms resemble a DM
interaction, the third term in Eq. (14a) gives rise to an in-plane symmetric
interaction. The large DM interaction (which we will show in the following)
compared to the one in TI thin films [72] can lead to interesting phenomena
such as spin Nernst effect, the appearance of the nontrivial topology, and
topological spin excitations [73, 74, 75]. The DM interaction has been also
explored extensively in recent years for spintronic applications [52, 53, 51,
76]. These off-diagonal terms enable the full RKKY Hamiltonian to provide
noncollinear twisted alignment between the magnetic impurities.

Let us define the form factor between Dirac points and $\\{X_{1}$ and
$X_{2}\\}$ points from Eq. (9)

$f(\vec{R})=2\,[e^{{\tt i}X_{1}R_{x}}\cos(\Lambda_{x}R_{x})+e^{{\tt
i}X_{2}R_{y}}\cos(\Lambda_{y}R_{y})]\,.$ (15)

For clarity, we work with the normalized exchange interactions,
$\mathcal{C}\,J^{2}\,|f(\vec{R})|^{2}$, where $\mathcal{C}=4\pi/S^{2}_{\rm
SBZ}v^{4}_{\rm F}$ thus,
$\widetilde{\mathcal{J}}(R)=\mathcal{J}(R)/\mathcal{C}\,J^{2}\,|f(\vec{R})|^{2}$
suppresses the beating type of oscillations from the multiple surface Dirac
cones [49], as shown in Fig. C.1 of Appendix C. With tthis, we show the
results for a single Dirac cone in what follows. The modulations captured in
$|f(\vec{R})|^{2}$ are different than those of graphene and other systems [77,
78, 79, 80] and can be useful to modulate the couplings on medium-length
scales $1/\Lambda_{x,y}$. Crucially, such a length scale only exists in TCIs
and is potentially tunable via, e.g., strain.

For numerical purposes, we first define
$\mathcal{E}-\mu=\mathcal{E}^{\prime}$, owing to the electron-hole symmetry
[79, 67, 68, 66], to set the upper limit of the integral in Eq. (13) to zero.
Then, we use the imaginary energy representation $\mathcal{E}^{\prime}+{\tt
i}o^{+}={\tt i}\omega$ to get rid of the numerical singularities and tight
tolerances in the real energy integration. After simple algebra calculations
and defining $\widetilde{\omega}^{2}=(\omega-{\tt i}\mu)^{2}+(\Delta^{2}/4)$,
the components of RKKY Hamiltonians for the first configuration (impurities on
the same sublattices) read as:

$\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{x}(R)=$
$\displaystyle+\text{Im}\int_{o^{+}}^{\infty}{\tt
i}d\omega\,\widetilde{\omega}^{2}\big{[}\mathbb{K}^{2}_{0}(\widetilde{\omega}\,R/v_{\rm
F})$
$\displaystyle\qquad+\cos(2\varphi_{R})\mathbb{K}^{2}_{1}(\widetilde{\omega}\,R/v_{\rm
F})\big{]},$ (16a)
$\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{y}(R)=$
$\displaystyle+\text{Im}\int_{o^{+}}^{\infty}{\tt
i}d\omega\,\widetilde{\omega}^{2}\big{[}\mathbb{K}^{2}_{0}(\widetilde{\omega}\,R/v_{\rm
F})$
$\displaystyle\qquad-\cos(2\varphi_{R})\mathbb{K}^{2}_{1}(\widetilde{\omega}\,R/v_{\rm
F})\big{]},$ (16b)
$\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{z}(R)=$
$\displaystyle+\text{Im}\int_{o^{+}}^{\infty}{\tt
i}d\omega\,\widetilde{\omega}^{2}\big{[}\mathbb{K}^{2}_{0}(\widetilde{\omega}\,R/v_{\rm
F})$ $\displaystyle\qquad+\mathbb{K}^{2}_{1}(\widetilde{\omega}\,R/v_{\rm
F})\big{]}-\frac{\Delta^{2}}{2}\text{Im}\int_{o^{+}}^{\infty}{\tt
i}d\omega\mathbb{K}^{2}_{0}(\widetilde{\omega}\,R/v_{\rm F}),$ (16c)
$\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{xy}(R)=$
$\displaystyle+\sin(2\varphi_{R})\text{Im}\int_{o^{+}}^{\infty}{\tt
i}d\omega\,\widetilde{\omega}^{2}\,\mathbb{K}^{2}_{1}(\widetilde{\omega}\,R/v_{\rm
F}),$ (16d) $\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{{\rm
DM},x}(R)=$
$\displaystyle-2\sin(\varphi_{R})\,\text{Im}\int_{o^{+}}^{\infty}d\omega\,\widetilde{\omega}\sqrt{\widetilde{\omega}^{2}-\Delta^{2}/4}$
$\displaystyle\hskip
28.45274pt\,\times\mathbb{K}_{0}(\widetilde{\omega}\,R/v_{\rm
F})\,\mathbb{K}_{1}(\widetilde{\omega}\,R/v_{\rm F})\,,$ (16e)
$\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{{\rm DM},y}(R)=$
$\displaystyle-2\cos(\varphi_{R})\,\text{Im}\int_{o^{+}}^{\infty}d\omega\,\widetilde{\omega}\sqrt{\widetilde{\omega}^{2}-\Delta^{2}/4}$
$\displaystyle\hskip
28.45274pt\,\times\mathbb{K}_{0}(\widetilde{\omega}\,R/v_{\rm
F})\,\mathbb{K}_{1}(\widetilde{\omega}\,R/v_{\rm F})\,.$ (16f)

Due to the symmetrical form of the states $|1\rangle$ and $|2\rangle$, the $x$
and $y$ components of the XYZ-Heisenberg interaction differ only by a sign in
the off-diagonal components of the real-space Green’s functions. Moreover, the
DM components are shifted by the $\pi/2$ phase.

For the second configuration (impurities on different sublattices), we find

$\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\beta}_{x}(R)={}$
$\displaystyle+\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{x}(R)\,,$ (17a)
$\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\beta}_{y}(R)={}$
$\displaystyle-\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{y}(R)\,,$ (17b)
$\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\beta}_{z}(R)={}$
$\displaystyle-\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{z}(R)\,,$ (17c)
$\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\beta}_{{\rm DM},x}(R)={}$
$\displaystyle-\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{{\rm DM},x}(R)\,,$
(17d) $\displaystyle\widetilde{\mathcal{J}}^{\alpha\,\beta}_{{\rm
DM},z}(R)={}$
$\displaystyle-\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{xy}(R)\,,$ (17e)

stemming from the spin-orbit structure of the model, and not from the spatial
symmetry of the SnTe(001) surface.

The magnetic impurities can reside on the same, Fig. 2(a), or different, Fig.
2(b), sublattices. Other configurations, such as impurities on the bonds or in
the center of the unit cell, can be simply obtained by adding together the
results for the RKKY interaction on the nearby sites. For example, if one
magnetic impurity is located halfway between a Sn site at location
$a_{0}\hat{y}/2$ and Te site at $-a_{0}\hat{y}/2$ ($a_{0}$ is the lattice
constant), then the interaction with an impurity located at a distance
$\vec{R}$ away on a Sn site is given by $\mathcal{H}^{\rm bond}_{\rm
RKKY}=\mathcal{H}^{\rm SnSn}_{\rm
RKKY}(\vec{R}-a_{0}\hat{y}/2)+\mathcal{H}^{\rm TeSn}_{\rm
RKKY}(\vec{R}+a_{0}\hat{y}/2)$. Note that the microscopic coupling $J$ for
this case will be different from the site-located impurity, and will depend on
the sublattice. Generalizations to other impurity positions are
straightforward.

Figure 2: A simple schematic for the position of magnetic impurities on (a)
the same and (b) different sublattices of the SnTe(001) surface. Cation Sn2+,
anion Te2-, and the local magnetic impurities are, respectively, shown by the
black, red, and blue spheres.

In general, for the sublattice setup, we continue only with $\alpha_{\rm
s}=+1$ in Eq. (14) henceforth for the magnetic impurities on the same TeTe
sublattices. A sign change in $\widetilde{\mathcal{J}}^{\alpha\alpha}_{{\rm
DM},y}$ and $\widetilde{\mathcal{J}}^{\alpha\alpha}_{xy}$ gives rise to the
results for magnetic impurities on the same SnSn sublattices; for different
TeSn sublattices, one can simply use Eq. (17). Finally, the sign of
$\widetilde{\mathcal{J}}^{\alpha\beta}_{{\rm DM},z}$ should be swapped for the
results on different SnTe sublattices.

##  IV Results and discussion

The role of the chemical potential $\mu$ and the optical gap $\Delta$ in
various terms of the RKKY interaction form the main messages of the present
work for both spintronic applications and QAHE physics. The Fermi sea can be
drained out if one keeps increasing the driving strength for a fixed impurity
separation $R$ and direction $\varphi_{R}$. This, in turn, manipulates the
surface states involved in the RKKY coupling. It is very important to control
the chemical potential to bring the surface states into the so-called low-
energy regime. This can be tuned with, e.g., a back gate from the thin film
substrates [81, 82].

It is necessary to mention that the Floquet parameter $a/\Omega$ mainly tunes
the light-matter interaction effect on the SnTe(001) surface, which should be
small in the weak drive regime (to avoid heating of the sample). It is
worthwhile commenting that the alignment of magnetic moments can be influenced
by both $a$ and $\Omega$. In Fig. D.1 of Appendix D, we systematically address
this matter. As soon as we turn on the light, independent of the symmetric and
asymmetric contributions to the RKKY interaction, magnetic moments talk to
each other only beneath the line $a\propto\sqrt{\Omega}$. The bandwidth of our
model is given by $2n=110$ meV. To have a broader range of light intensity for
tuning the RKKY interaction (less zero response), we continue with an off-
resonant high-frequency $\Omega=1$ eV $\propto 270$ THz much larger than the
bandwidth. Thus, the gap of the system can reach $\Delta\simeq 23$ K [83],
i.e., $a\simeq 47$ meV. However, we vary intensity $a$ up to 200 meV (optical
gap $\Delta$ up to 80 meV) in what follows.

From the Eqs. (16a)-(16f) and Eqs. (17a)-(17e), one can observe the
periodicity of interactions in $\varphi_{R}$ – except
$\varphi_{R}$-independent $\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{z}(R)$ –
which provides $\mathcal{H}_{\rm RKKY}(\varphi_{R})=\mathcal{H}_{\rm
RKKY}(\varphi_{R}+\pi)$. Due to the square lattice of the host SnTe(001)
surface, three angles $\varphi_{R}=0$, $\pi/4,$ and $\pi/2$ between the
magnetic impurities lead to special RKKY interactions. A quick analysis
results in $\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{xy}(R)=0$ for both $0$
and $\pi/2$ and $\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{{\rm
DM},x/y}(R)=0$, respectively, for $0$ and $\pi/2$. While for $\pi/4$, we find
$\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{x}(R)=\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{y}(R)$
and $\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{{\rm
DM},x}(R)=\widetilde{\mathcal{J}}^{\alpha\,\alpha}_{{\rm DM},y}(R)$. But, for
noncollinear twisted alignments of magnetic moments on the surface with more
adjustability, we continue with an in-plane angle, e.g., $\varphi_{R}=\pi/3$.
We set the chemical potential and the optical gap at $\mu=27.5$ meV and
$\Delta=45$ meV, respectively, where necessary. In the driven-doped surface,
this set helps to have effective interactions because of $\mu>\Delta/2$ (see
interpretation below).

Figure 3 shows the RKKY terms as a function of impurity separation $R/a_{0}$.
For the pristine (undriven-undoped) surface, $\mu=0$ and $\Delta=0$, which
leads to real $\widetilde{\omega}=\omega$, the analytical solution of Eqs.
(16a)-(16d) is already known [49, 50, 80, 67, 68, 66] and can
straightforwardly be found with the help of Eq. (23) in Appendix E such that
the non-zero interactions decay as $R^{-3}$, see Fig. 3(a). Under the same
conditions, the DM terms in Eqs. (16e) and (16f) vanish
$\widetilde{\mathcal{J}}^{\alpha\alpha}_{{\rm DM},x/y}(R)=0$ due to the
absence of imaginary energy. The physical reason for a vanishing DM
interaction stems from the presence of inversion symmetry between sublattices
(absence of $\mu$) for which asymmetric alignment of spins from the interplay
of diagonal and off-diagonal spin-orbital states does not take place. The
presence of DM interaction due to doping has been experimentally confirmed in
TI thin films [52]. Interestingly, only
$\widetilde{\mathcal{J}}^{\alpha\alpha}_{x}(R)$ forms an AFM order since the
response from different spin-orbital states, $\\{|1\rangle$ and $|2\rangle\\}$
or $\\{|2\rangle$ and $|1\rangle\\}$, are three times larger than the same
spin-orbital states, $\\{|1\rangle$ and $|1\rangle\\}$ or $\\{|2\rangle$ and
$|2\rangle\\}$, resulting in
$\widetilde{\mathcal{J}}^{\alpha\alpha}_{x}(R)<0$.

Figure 3: RKKY interactions on the SnTe(001) surface at $\varphi_{R}=\pi/3$ as
a function of impurity separation $R$ when the surface is (a) pristine
(undriven-undoped), (b) undriven-doped, (c) driven-undoped, and (d) driven-
doped. We set $\mu=27.5$ meV and $\Delta=45$ meV, respectively, for the
chemical potential and optical gap. It can be observed that the asymmetric DM
term exists only on the doped surface. Furthermore, the Floquet drive
suppresses all terms due to the reduced density of mediating states.

For the undriven-doped surface, $\mu\neq 0$ and $\Delta=0$, which leads to
imaginary $\widetilde{\omega}=\omega-{\tt i}\mu$, one needs Eq. (24) in
Appendix E. It is obvious that the short-range RKKY interactions are given by
the power-law $R^{-3}$, while for the long-range couplings, oscillatory Meijer
functions [84, 85] emerge, as presented in Fig. E.1. Our numerical data in
Fig. 3(b) confirm this dependence. In contrast to the pristine surface, DM
interactions wake up when the surface is doped for the same reason of
inversion symmetry breaking between sublattices.

For the driven-undoped surface, $\mu=0$ and $\Delta\neq 0$, which leads to
real $\widetilde{\omega}=\sqrt{\omega^{2}+\Delta^{2}/4}$, with the help of Eq.
(25) in Appendix E, we find short-range RKKY interactions $\propto R^{-3}$,
while for the long-range interactions, one needs to use $\lim_{\Delta R/v_{\rm
F}\to\infty}\mathbb{K}_{n}(\Delta R/v_{\rm F})\simeq\sqrt{\pi v_{\rm
F}/2\Delta R}\exp(-\Delta R/v_{\rm F})$ and the Laplace method to obtain RKKY
interactions $\propto(R/\Delta)^{-3/2}\exp(-2\Delta R/v_{\rm F})$ [66], which
are again in agreement with the numerical results in Fig. 3(c). For the reason
mentioned before about the doping effect on the DM interaction, one still
expects a vanishing DM response for the driven-undoped surface.

For the driven-doped surface in Fig. 3(d), we consider the $\mu>\Delta/2$
regime where the chemical potential is outside the optical gap so the states
added by doping play role in RKKY interactions (for $\mu<\Delta/2$ the DM
terms vanish). Although it is not easy to find analytical expressions for the
driven-doped surface, we expect a modified version of Meijer functions for the
long-range coupling as the interactions behave like the undriven-doped
surface. However, we believe that the power-law $R^{-3}$ is still valid for
short-range couplings. As both driving and doping are present, the drastic
change of oscillations for all RKKY components is evident. This is a direct
consequence of the optical gap for which the corresponding density of
mediating states in the RKKY interaction decreases, and thus, most of the
doped states are washed out. For this reason, the oscillations are accompanied
by smaller amplitudes compared to the undriven-doped surface. This is mainly
where our novel contribution pays off in QAHE physics since magnetic
characters focus more on purely positive (FM) and negative (AFM) signs.

It is also worth mentioning that various FM and AFM characters emerge for
symmetric XYZ-Heisenberg $\widetilde{\mathcal{J}}_{i}$ and in-plane
$\widetilde{\mathcal{J}}_{xy}$ components depending on the interplay between
driving and doping, as well as between the impurity separation and direction.
But, for the change of signs in DM components, one needs to take into account
the chirality. When DM is positive, the coupling between magnetic impurities
favors clockwise rotation, meaning that it tends to align right handedly more
naturally than left. Similarly, a negative DM interaction favors counter-
clockwise rotation for magnetic moments.

Figure 4: RKKY interactions for (a) $\mu=13.75$ meV and (b) $\mu=27.5$ meV at
$R/a_{0}=50$ and $\varphi_{R}=\pi/3$ as a function of the optical gap $\Delta$
on the SnTe(001) surface. Independent of the chemical potential, RKKY
responses vanish for $\Delta\geq 2|\mu|$ once the chemical potential lies in
the gap, due to their exponential decay with $R$.

To gain further insights into the gap dependency of interactions, we proceed
to plot the RKKY couplings against the optical gap $\Delta$ for $\mu=13.75$
meV, Fig. 4(a), and $\mu=27.5$ meV, Fig. 4(b), at $R/a_{0}=50$ and
$\varphi_{R}=\pi/3$. As soon as we turn on the Floquet drive, all interactions
start to emerge due to the presence of $\mu$. It is noteworthy that once the
gap edge lies at the chemical potential level, i.e., when $\Delta=2|\mu|$,
they start to vanish because the doped states start to lie in the gap without
any significant density of mediating states. If we dope the surface stronger,
Fig. 4(b), RKKY responses become stronger compared to Fig. 4(a). This can be
easily understood by taking into account that the density of carriers on the
surface with $\mu=27.5$ meV is larger than at $\mu=13.75$ meV; thus, more
itinerant carriers act as mediators for the indirect exchange coupling between
impurities.

Next, we systematically explore the RKKY interactions of driven-doped
SnTe(001) surface at fixed direction $\varphi_{R}=\pi/3$ and impurity
separation $R/a_{0}=50$. In Fig. 5, at first glance, one would observe that
the amplitude of oscillations increases as more surface states (larger $\mu$)
are involved in the exchange coupling. However, the gap opening reduces the
rate of involved states and for a certain chemical potential $\mu$, RKKY
interactions are non-zero only below $\Delta=2|\mu|$. For example, at
$\mu=\pm\,27.5$ meV, one obtains zero RKKY responses at $\Delta\geq 50$ meV,
as confirmed in Figs. 5(a)-(f), implying that there is no effective mediating
state anymore on the surface (at least, within our approximations) so that the
RKKY responses approach zero. The red and blue colors highlight the character
of FM and AFM (clockwise and counter-clockwise spin rotations) for symmetric
(asymmetric) RKKY components, while the yellow color displays nearly zero
responses.

Figure 5: RKKY interactions at $\varphi_{R}=\pi/3$ and $R/a_{0}=50$ as a
function of both optical gap $\Delta$ and chemical potential $\mu$ on the
driven-doped SnTe(001) surface. Red (blue) color refers to FM (AFM) order in
(a)-(d), while they refer to clockwise (contour-clockwise) spin rotation in
the DM terms of (e) and (f). The yellow color depicts nearly zero RKKY
response. The dominant contribution to the RKKY interaction belongs to the $z$
($x$)-component of symmetric (asymmetric) interaction, i.e.,
$\widetilde{\mathcal{J}}^{\alpha\alpha}_{z}$
($\widetilde{\mathcal{J}}^{\alpha\alpha}_{{\rm DM},x}$). As expected, no
response appears for $\Delta\geq 2|\mu|$ due to the vanishing density of
mediating states.

Additionally, the dominant oscillation amplitude belongs to the $z$-component
of the XYZ-Heisenberg, though the $x$-component of the DM term is the next
prevalent among others. As mentioned before, DM interaction enables the system
to be applicable for emergent spintronic applications [52, 53, 51, 76]. This
does not occur on the surface of conventional TI thin films [72]. Therefore,
our large DM compatible with $\widetilde{\mathcal{J}}^{\alpha\alpha}_{z}(R)$
is a very useful feature of TCIs.

It is also worth exploring the intrinsic electron-hole symmetry on the driven-
doped SnTe(001) surface. The electron-hole symmetry is valid for all symmetric
interactions, XYZ-Heisenberg and in-plane $xy$, independent of doping and
optical gap (also independent of the impurity separation and position of
impurities – not shown here). For the DM interaction at a given direction, one
immediately extracts the following relation

$\widetilde{\mathcal{J}}^{\alpha\alpha}_{{\rm
DM},x/y}(-\mu,\Delta)=-\widetilde{\mathcal{J}}^{\alpha\alpha}_{{\rm
DM},x/y}(+\mu,\Delta)\,,$ (18)

which has already been confirmed experimentally on the surface of TI thin
films [52], as also shown in Figs. 5(e) and 5(f). It has an opposite sign when
the chemical potential cuts through the upper or lower band. We emphasize that
the presence of inversion symmetry between sublattices and the corresponding
spin-orbital states is the main origin of the above relationship.

For the experimental perspectives, one would extract the following critical
Floquet driving amplitude $E^{\rm c}_{0}$ (considering an arbitrary driving
frequency $\Omega$) in the low-energy regime of SnTe(001) surface above which
a nearly zero RKKY coupling is achieved:

$E^{\rm c}_{0}={}\frac{\sqrt{n}}{ev_{\rm F}}\Omega^{3/2}\,.$ (19)

For instance, in our case, for the selected light frequency of $\Omega\simeq
1$ eV, the Fermi velocity of $v_{\rm F}\simeq 2$ eV. Å, and $n=55$ meV, one
finds $E^{\rm c}_{0}\simeq 10^{9}$ V/m, which is nowadays easily realizable in
experiment [86, 87, 88]. Having this information, we would suggest using a
light intensity below $E^{\rm c}_{0}$ in experiments so tuning the RKKY
coupling with light can play role in determining the prerequisites for
spintronic applications and QAHE. Experimentally, single-atomic magnetometry
and magnetotransport through scanning tunneling microscopy and angle-resolved
photoemission spectroscopy can measure the RKKY interaction [89, 90, 91].

##  V Outlook

Before closing the discussion, let us briefly discuss possible avenues for
upcoming research. The present study for the “weak” light-induced RKKY
interaction, which still reports applicable findings, can also be addressed
for the “strong” driving regime. First, in the strong driving regime, the full
Floquet Hamiltonian is required for further photon processes [59]; second, the
nonlinear corrections should be solved numerically [59, 61]. The extended
continuous strong driving must inevitably lead to heating, which without
remediation would destroy the sample. We can propose a heat sink, e.g., a
high-quality Al [92], which maintains a low system temperature despite the
Floquet drive. Beyond the weak driving regime, our calculations could also be
done for a “linearly” polarized light if one considers the $n$th-order Bessel
functions for the corrections up to infinite orders [93].

To address the low-frequency limit beyond our high-frequency expansion, band
crossings should also be considered in the theory; this requires another
expansion of the effective Floquet Hamiltonian and extra terms will be
generated. Additionally, the finite temperature effects can also be explored
for both short- and long-range couplings within the finite-temperature self-
consistent field approximation [94, 95].

Note that the drive used to tune the RKKY interaction also opens up a
topological gap in the surface Dirac cones. With four Dirac cones of identical
chirality, it will correspond to Chern number 2, suggesting that TCIs may also
be a useful playground for light-induced anomalous quantum Hall physics.
Unlike existing two-dimensional systems, such as graphene [88], we propose to
realize the gapped TCI surface states by continuous, moderate amplitude
optical driving, which allowed the realization of an effective Fermi surface
earlier in the paper, Sec IV. This also implies the quasi-static topology of
the surface electronic states, which will be an interesting topic for future
exploration. However, unlike graphene, this topological state is complicated
by the presence of gapless edge states on the other [non-(001)] surfaces as
well as the bottom (001) surface of the material. We leave an analysis of
topological transport in this geometry for future work.

##  VI Conclusions

Exploring the physics of magnetically doped topological systems in tuning QAHE
as well as for the spintronics community has triggered interest in condensed
matter physics. Recently, experimental photonic platforms have also become
increasingly urgent for tuning the properties of materials. To contribute to
these, we have employed the isotropic optically driven continuum model for
Dirac cones with nontrivial spin textures on the doped SnTe(001) surface as a
well-known TCI to investigate the (quasi)out-of-equilibrium physics of the
RKKY interaction between two magnetic impurities. This, in turn, aims at
providing spintronic applications as well as making the TCI-based QAHE strong.
To tune the features, in particular, we have focused on the weak driving
effects and off-resonant regime between the light and bandwidth of the Dirac
cones using the van Vleck inverse-frequency expansion. We make use of the bare
spin-orbital states of each sublattice, resulting in noncollinear twisted RKKY
interaction comprising of XYZ-Heisenberg, symmetric in-plane, and asymmetric
DM terms.

Preliminary analyses highlighted the $z$ ($x$)-component of the XYZ-Heisenberg
(DM) interaction as the first (second) dominant contribution to the total RKKY
interaction. Moreover, depending on the position of magnetic impurities,
chemical potential, and the optical gap, the driven-doped SnTe(001) surface
reaches various modulations of FM and AFM characters for the symmetric RKKY
terms as well as of the clockwise and counter-clockwise spin rotations for the
asymmetric DM ones. A systematic analysis of the RKKY coupling on the
interplay between doping and driving demonstrated that the efficient range of
chemical potential in controlling the amplitudes belongs to the strengths
outside the optical gap. The optical gap (Floquet drive) in the presence of
doping leads to a reduction of the RKKY amplitudes because of the decreased
total density of mediating states. This is where we propose a Floquet drive to
make the TCI-based QAHE strong when the surface is doped. Alongside, we find
non-zero DM interactions only for the doped surface because of inversion
symmetry breaking between sublattices. Moreover, to evaluate the intrinsic
electron-hole symmetry of the system, we found that the symmetry is broken
because of the DM term.

Providing reasonable light intensity and frequency compatible with the low-
energy Dirac bands of the SnTe(001) surface, insights from the present work
are discussed for feasible experimental observations. Finally, we would
mention that the large DM term on the surface of TCIs is highly desirable for
spintronic applications, highlighting their additional usefulness compared to
TIs due to the presence of multiple Dirac cones.

## Acknowledgments

M.Y. greatly thanks Wang-Kong Tse for helpful discussions. This work was
supported by the National Science Foundation through award numbers
DMR-1945529, PHY-1607611, NSF PHY1748958, and the Welch Foundation through
award number AT-2036-20200401. This project was funded by The University of
Texas at Dallas Office of Research and Innovation through the SPIRe program.
M.B. was supported by the Marie Skłodowska-Curie grant agreement No 890711
(until 01.09.2022).

##  Appendix A Components of the real-space Green’s function

The real-space Green’s functions of hybrid states in terms of bare sublattices
considering corresponding spins requires the following definition
$G^{\sigma\sigma^{\prime}}_{\alpha\beta}(\mathcal{E},\vec{R})=\sum_{\gamma,\eta}u^{*}_{\alpha,\gamma,\sigma}u_{\beta,\gamma,\sigma^{\prime}}\mathcal{G}_{\gamma\eta}(\mathcal{E},\vec{R})$,
where components are given by [in the following we omit $(\dots)$, so
$\mathcal{G}_{\gamma\eta}$ means
$\mathcal{G}_{\gamma\eta}(\mathcal{E},\vec{R})$]

$\displaystyle G^{\uparrow\uparrow}_{\rm Te\,Te}={}$
$\displaystyle\,\mathcal{G}_{11}/2\,,\quad G^{\uparrow\downarrow}_{\rm
Te\,Te}={}\,\mathcal{G}_{12}/2\,,\quad$ (20a) $\displaystyle
G^{\downarrow\uparrow}_{\rm Te\,Te}={}$
$\displaystyle\,\mathcal{G}_{21}/2\,,\quad G^{\downarrow\downarrow}_{\rm
Te\,Te}={}\,\mathcal{G}_{22}/2\,,$ (20b) $\displaystyle
G^{\uparrow\uparrow}_{\rm Te\,Sn}={}$
$\displaystyle\,\mathcal{G}_{12}/2\,,\quad G^{\uparrow\downarrow}_{\rm
Te\,Sn}={}\,\mathcal{G}_{11}/2\,,\quad$ (20c) $\displaystyle
G^{\downarrow\uparrow}_{\rm Te\,Sn}={}$
$\displaystyle\,\mathcal{G}_{22}/2\,,\quad G^{\downarrow\downarrow}_{\rm
Te\,Sn}={}\,\mathcal{G}_{21}/2\,,$ (20d) $\displaystyle
G^{\uparrow\uparrow}_{\rm Sn\,Te}={}$
$\displaystyle\,\mathcal{G}_{21}/2\,,\quad G^{\uparrow\downarrow}_{\rm
Sn\,Te}={}\,\mathcal{G}_{22}/2\,,\quad$ (20e) $\displaystyle
G^{\downarrow\uparrow}_{\rm Sn\,Te}={}$
$\displaystyle\,\mathcal{G}_{11}/2\,,\quad G^{\downarrow\downarrow}_{\rm
Sn\,Te}={}\,\mathcal{G}_{12}/2\,,$ (20f) $\displaystyle
G^{\uparrow\uparrow}_{\rm Sn\,Sn}={}$
$\displaystyle\,\mathcal{G}_{22}/2\,,\quad G^{\uparrow\downarrow}_{\rm
Sn\,Sn}={}\,\mathcal{G}_{21}/2\,,\quad$ (20g) $\displaystyle
G^{\downarrow\uparrow}_{\rm Sn\,Sn}={}$
$\displaystyle\,\mathcal{G}_{12}/2\,,\quad G^{\downarrow\downarrow}_{\rm
Sn\,Sn}={}\,\mathcal{G}_{11}/2\,.$ (20h)

##  Appendix B Components of the spin susceptibility

The spin susceptibility is described by
$\chi^{\alpha\beta}_{lj}(\,\vec{R}\,)=-\frac{2}{\pi}\,\Im\int^{\mu}_{-\infty}d\mathcal{E}\,\mathbb{F}^{\alpha\beta}_{lj}(\mathcal{E},\vec{R})$
with the following expressions for the first configuration (impurities on the
same sublattices)

$\displaystyle\mathbb{F}^{\rm\alpha\,\alpha}_{xx}={}$
$\displaystyle\frac{|f(\vec{R})|^{2}}{4}\left[2\,\mathcal{G}_{11}\mathcal{G}_{22}-\mathcal{G}^{2}_{12}-\mathcal{G}^{2}_{21}\right]\,,$
(21a) $\displaystyle\mathbb{F}^{\rm\alpha\,\alpha}_{yy}={}$
$\displaystyle\frac{|f(\vec{R})|^{2}}{4}\left[2\,\mathcal{G}_{11}\mathcal{G}_{22}+\mathcal{G}^{2}_{12}+\mathcal{G}^{2}_{21}\right]\,,$
(21b) $\displaystyle\mathbb{F}^{\rm\alpha\,\alpha}_{zz}={}$
$\displaystyle\frac{|f(\vec{R})|^{2}}{4}\left[2\,\mathcal{G}_{12}\mathcal{G}_{21}+\mathcal{G}^{2}_{11}+\mathcal{G}^{2}_{22}\right]\,,$
(21c) $\displaystyle\mathbb{F}^{\rm\alpha\,\alpha,\alpha_{\rm s}}_{xy/yx}={}$
$\displaystyle{\tt i}\alpha_{\rm
s}\,\frac{|f(\vec{R})|^{2}}{4}\left[\mathcal{G}^{2}_{21}-\mathcal{G}^{2}_{12}\right]\,,$
(21d) $\displaystyle\mathbb{F}^{\rm\alpha\,\alpha,\alpha_{\rm s}}_{xz/zx}={}$
$\displaystyle\pm\alpha_{\rm
s}\,\frac{|f(\vec{R})|^{2}}{4}\left[\mathcal{G}_{11}+\mathcal{G}_{22}\right]\left[\mathcal{G}_{21}-\mathcal{G}_{12}\right]\,,$
(21e) $\displaystyle\mathbb{F}^{\rm\alpha\,\alpha}_{yz/zy}={}$
$\displaystyle\mp{\tt
i}\frac{|f(\vec{R})|^{2}}{4}\left[\mathcal{G}_{11}+\mathcal{G}_{22}\right]\left[\mathcal{G}_{12}+\mathcal{G}_{21}\right]\,.$
(21f)

For the second configuration (impurities on different sublattices), we acheive

$\displaystyle\mathbb{F}^{\alpha\,\beta/\beta\,\alpha}_{xx}={}$
$\displaystyle\frac{|f(\vec{R})|^{2}}{4}\left[2\,\mathcal{G}_{11}\mathcal{G}_{22}-\mathcal{G}^{2}_{12}-\mathcal{G}^{2}_{21}\right]\,,$
(22a) $\displaystyle\mathbb{F}^{\alpha\,\beta/\beta\,\alpha}_{yy}={}$
$\displaystyle\frac{|f(\vec{R})|^{2}}{4}\left[-2\,\mathcal{G}_{11}\mathcal{G}_{22}-\mathcal{G}^{2}_{12}-\mathcal{G}^{2}_{21}\right]\,,$
(22b) $\displaystyle\mathbb{F}^{\alpha\,\beta/\beta\,\alpha}_{zz}={}$
$\displaystyle\frac{|f(\vec{R})|^{2}}{4}\left[-2\,\mathcal{G}_{12}\mathcal{G}_{21}-\mathcal{G}^{2}_{11}-\mathcal{G}^{2}_{22}\right]\,,$
(22c) $\displaystyle\mathbb{F}^{\rm\alpha\,\beta,\alpha_{\rm d}}_{xy/yx}={}$
$\displaystyle\pm{\tt i}\alpha_{\rm
d}\,\frac{|f(\vec{R})|^{2}}{4}\left[\mathcal{G}^{2}_{12}-\mathcal{G}^{2}_{21}\right]\,,$
(22d) $\displaystyle\mathbb{F}^{\rm\alpha\,\beta,\alpha_{\rm d}}_{xz/zx}={}$
$\displaystyle\alpha_{\rm
d}\,\frac{|f(\vec{R})|^{2}}{4}\left[\mathcal{G}_{11}+\mathcal{G}_{22}\right]\left[\mathcal{G}_{12}-\mathcal{G}_{21}\right]\,,$
(22e) $\displaystyle\mathbb{F}^{\rm\alpha\,\beta}_{yz/zy}={}$
$\displaystyle\pm{\tt
i}\frac{|f(\vec{R})|^{2}}{4}\left[\mathcal{G}_{11}+\mathcal{G}_{22}\right]\left[\mathcal{G}_{12}+\mathcal{G}_{21}\right]\,.$
(22f)

##  Appendix C Beating oscillation of RKKY interaction

Figure C.1: Beating type of RKKY oscillations due to the multiple Dirac cones
on the SnTe(001) surface, Eq. (15), for impurity separations $R/a_{0}=10$ and
30 as a function of the polar angle $\varphi_{R}$.

Here, we aim at showing the beating type of oscillations from multiple Dirac
cones on the SnTe(001) surface, see Eq. 15 for $|f(\vec{R})|^{2}$, as shown in
Fig. C.1. The momenta contributed to this beating function are
$X_{1/2}\pm\Lambda_{x/y}$ and $2\Lambda_{x/y}$ on the SBZ [49].

##  Appendix D Floquet engineering of magnetic moments alignment on the $\rm
SnTe$(001) surface

Figure D.1: Magnetic phase diagram of RKKY interactions on the driven-doped
SnTe(001) surface at $R/a_{0}=50$ and $\varphi_{R}=\pi/3$. Depending on the
interplay between light intensity $a$ and frequency $\Omega$, various FM/AFM
and clockwise/contour-clockwise characters can be seen for symmetric and
asymmetric components, respectively. Also,
$\widetilde{\mathcal{J}}^{\alpha\alpha}_{z}(R)>\widetilde{\mathcal{J}}^{\alpha\alpha}_{{\rm
DM},x}(R)$ holds true.

An important motivation for Fig. D.1 presented here originates from the
spintronics community where the quest for the control of magnetic phases is
still active. Among the fundamental open questions is whether the magnetic
moments alignment on the SnTe(001) surface with a nontrivial spin texture is
adjustable through the interplay between the intensity and frequency of light.
This is somehow important for the Floqeut engineering of surface. As can be
seen, phase diagrams for symmetric and asymmetric contributions to the RKKY
interaction highlight the characters of FM/AFM (for symmetric contributions in
Figs. D.1(a)-(d)) as well as clockwise/contour-clockwise (for asymmetric
contributions in Figs. D.1(e)-(f)) in different controlled ways, characterized
by red and blue colors. In a certain regime above the line
$a\propto\sqrt{\Omega}$, the RKKY interaction vanishes. Getting away from the
alignments, we have the relation
$\widetilde{\mathcal{J}}^{\alpha\alpha}_{z}(R)>\widetilde{\mathcal{J}}^{\alpha\alpha}_{{\rm
DM},x}(R)$ between the strongest symmetric and asymmetric RKKY components.

##  Appendix E Solution of integrals in Eq. (16) for undriven-undoped,
undriven-doped, and driven-undoped $\rm SnTe$(001) surface

Figure E.1: Oscillatory Meijer functions, Eq. (24), versus $x=\mu R/v_{\rm F}$
for undriven-doped SnTe(001) surface, resulting in oscillatory RKKY
interactions in Fig. 3(b).

In this part, we present the analytical expressions for the integrals in Eq.
(16). For the pristine surface, we use the following relation for
$n=\\{0,1\\}$:

${\rm Im}\int^{\infty}_{o^{+}}{\tt
i}d\omega\,\omega^{2}\mathbb{K}^{2}_{n}(\omega R/v_{\rm
F})=+(2n+1)\frac{\pi^{2}v^{3}_{\rm F}}{32R^{3}}\,.$ (23)

For the undriven-doped surface, the following relations are required to
understand the short- and long-range responses:

$\displaystyle\mathbb{K}_{n}({\tt i}u)={}(-1)^{n+1}\frac{\pi}{2}e^{{\tt
i}n\pi/2}\big{[}\mathbb{Y}_{n}(u)+{\tt i}\mathbb{J}_{n}(u)\big{]}\,,$ (24a)
$\displaystyle{\rm Im}\int^{\infty}_{o^{+}}{\tt i}d\omega[\omega-{\tt
i}\mu]^{2}\mathbb{K}^{2}_{n}([\omega-{\tt i}\mu]R/v_{\rm F})={}$
$\displaystyle+(2n+1)\frac{\pi^{2}v^{3}_{\rm
F}}{32R^{3}}-(-1)^{n}\frac{\pi^{2}\mu^{3}}{4\sqrt{\pi}}\underbrace{\mathbb{G}^{\,*,2}_{1,2}\lparen\begin{smallmatrix}4\\\
-\tfrac{1}{2},\tfrac{1}{2}\end{smallmatrix}|\,0,n,-\tfrac{3}{2},-n\rparen{\frac{\mu^{2}R^{2}}{v^{2}_{\rm
F}}}}_{\mathcal{M}_{n}(\mu R/v_{\rm F})},$ (24b)

where $\mathcal{M}_{n}(x)$ is the special Meijer function [84, 85], as shown
in Fig. E.1, which leads to oscillatory RKKY interactions in Fig. 3(b).

For the driven-undoped surface, we use the following mathematical identity
[96]

$\displaystyle\int_{-\infty}^{0}x^{j-1}dx\,\mathbb{K}_{n}(ax)\mathbb{K}_{m}(ax)={}\frac{2^{j-3}}{a^{j}\Gamma(j)}\Gamma([j+n+m]/2)$
(25)
$\displaystyle\quad\Gamma([j+n-m]/2)\Gamma([j-n+m]/2)\Gamma([j-n-m]/2)\,,$

where $\Gamma(\dots)$ is the special Gamma function.

## References

  * Hasan and Kane [2010] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 (2010).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Rev. Mod. Phys. 83, 1057 (2011).
  * Zhang  _et al._ [2009] H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C. Zhang, Nature Physics 5, 438 (2009).
  * Moore and Balents [2007] J. E. Moore and L. Balents, Phys. Rev. B 75, 121306 (2007).
  * Yu  _et al._ [2010] R. Yu, W. Zhang, H.-J. Zhang, S.-C. Zhang, X. Dai, and Z. Fang, Science 329, 61 (2010).
  * Chang  _et al._ [2013] C.-Z. Chang _et al._ , Science 340, 167 (2013).
  * Kou  _et al._ [2014] X. Kou _et al._ , Phys. Rev. Lett. 113, 137201 (2014).
  * Checkelsky  _et al._ [2014] J. G. Checkelsky, R. Yoshimi, A. Tsukazaki, K. S. Takahashi, Y. Kozuka, J. Falson, M. Kawasaki, and Y. Tokura, Nature Physics 10, 731 (2014).
  * Chi and Moodera [2022] H. Chi and J. S. Moodera, APL Materials 10, 090903 (2022).
  * Zhou  _et al._ [2006] Z. Zhou, Y.-J. Chien, and C. Uher, Phys. Rev. B 74, 224418 (2006).
  * Zhang  _et al._ [2012] D. Zhang _et al._ , Phys. Rev. B 86, 205127 (2012).
  * Hor  _et al._ [2010] Y. S. Hor _et al._ , Phys. Rev. B 81, 195203 (2010).
  * Chang  _et al._ [2015] C.-Z. Chang, W. Zhao, D. Y. Kim, H. Zhang, B. A. Assaf, D. Heiman, S.-C. Zhang, C. Liu, M. H. W. Chan, and J. S. Moodera, Nature Materials 14, 473 (2015).
  * Peixoto  _et al._ [2016] T. R. F. Peixoto _et al._ , Phys. Rev. B 94, 195140 (2016).
  * Wang  _et al._ [2015] J. Wang, B. Lian, and S.-C. Zhang, Phys. Rev. Lett. 115, 036805 (2015).
  * Liu  _et al._ [2009] Q. Liu, C.-X. Liu, C. Xu, X.-L. Qi, and S.-C. Zhang, Phys. Rev. Lett. 102, 156603 (2009).
  * Kasuya [1956] T. Kasuya, Progress of Theoretical Physics 16, 45 (1956).
  * Yosida [1957] K. Yosida, Phys. Rev. 106, 893 (1957).
  * Ruderman and Kittel [1954] M. A. Ruderman and C. Kittel, Phys. Rev. 96, 99 (1954).
  * Parhizgar  _et al._ [2013] F. Parhizgar, H. Rostami, and R. Asgari, Phys. Rev. B 87, 125401 (2013).
  * Klinovaja and Loss [2013] J. Klinovaja and D. Loss, Phys. Rev. B 87, 045422 (2013).
  * Dmitrienko  _et al._ [2010] V. E. Dmitrienko, E. N. Ovchinnikova, J. Kokubun, and K. Ishida, JETP Letters 92, 383 (2010).
  * Awschalom and Flatté [2007] D. D. Awschalom and M. E. Flatté, Nature Physics 3, 153 (2007).
  * Hanson  _et al._ [2007] R. Hanson, L. P. Kouwenhoven, J. R. Petta, S. Tarucha, and L. M. K. Vandersypen, Rev. Mod. Phys. 79, 1217 (2007).
  * Garate and Franz [2010] I. Garate and M. Franz, Phys. Rev. B 81, 172408 (2010).
  * Abanin and Pesin [2011] D. A. Abanin and D. A. Pesin, Phys. Rev. Lett. 106, 136802 (2011).
  * Efimkin and Galitski [2014] D. K. Efimkin and V. Galitski, Phys. Rev. B 89, 115431 (2014).
  * Gao  _et al._ [2009] J. Gao, W. Chen, X. C. Xie, and F.-c. Zhang, Phys. Rev. B 80, 241302 (2009).
  * Kurilovich  _et al._ [2016] P. D. Kurilovich, V. D. Kurilovich, and I. S. Burmistrov, Phys. Rev. B 94, 155408 (2016).
  * Yevtushenko and Yudson [2018] O. M. Yevtushenko and V. I. Yudson, Phys. Rev. Lett. 120, 147201 (2018).
  * Fu [2011] L. Fu, Phys. Rev. Lett. 106, 106802 (2011).
  * Hsieh  _et al._ [2012] T. H. Hsieh, H. Lin, J. Liu, W. Duan, A. Bansil, and L. Fu, Nature Communications 3, 982 (2012).
  * Tanaka  _et al._ [2012] Y. Tanaka, Z. Ren, T. Sato, K. Nakayama, S. Souma, T. Takahashi, K. Segawa, and Y. Ando, Nature Physics 8, 800 (2012).
  * Dziawa  _et al._ [2012] P. Dziawa _et al._ , Nature Materials 11, 1023 (2012).
  * Xu  _et al._ [2012] S.-Y. Xu _et al._ , Nature Communications 3, 1192 (2012).
  * Kundu  _et al._ [2016] A. Kundu, H. A. Fertig, and B. Seradjeh, Phys. Rev. Lett. 116, 016802 (2016).
  * Wang  _et al._ [2013a] Y. J. Wang, W.-F. Tsai, H. Lin, S.-Y. Xu, M. Neupane, M. Z. Hasan, and A. Bansil, Phys. Rev. B 87, 235317 (2013a).
  * Ezawa [2014] M. Ezawa, Phys. Rev. B 89, 195413 (2014).
  * Oka and Aoki [2009] T. Oka and H. Aoki, Phys. Rev. B 79, 081406 (2009).
  * Oka and Kitamura [2019] T. Oka and S. Kitamura, Annual Review of Condensed Matter Physics 10, 387 (2019).
  * Kitagawa  _et al._ [2011] T. Kitagawa, T. Oka, A. Brataas, L. Fu, and E. Demler, Phys. Rev. B 84, 235108 (2011).
  * Platero and Aguado [2004] G. Platero and R. Aguado, Physics Reports 395, 1 (2004).
  * Nasu [2004] K. Nasu,  _Photoinduced Phase Transitions_ (World Scientific, 2004).
  * Yan and Wang [2016] Z. Yan and Z. Wang, Phys. Rev. Lett. 117, 087402 (2016).
  * Narayan [2016] A. Narayan, Phys. Rev. B 94, 041409 (2016).
  * Ezawa [2017] M. Ezawa, Phys. Rev. B 96, 041205 (2017).
  * Kundu and Seradjeh [2013] A. Kundu and B. Seradjeh, Phys. Rev. Lett. 111, 136402 (2013).
  * Yang  _et al._ [2021] Z. Yang, Q. Yang, J. Hu, and D. E. Liu, Phys. Rev. Lett. 126, 086801 (2021).
  * Yarmohammadi and Cheraghchi [2020] M. Yarmohammadi and H. Cheraghchi, Phys. Rev. B 102, 075411 (2020).
  * Cheraghchi and Yarmohammadi [2021] H. Cheraghchi and M. Yarmohammadi, Scientific Reports 11, 5273 (2021).
  * Fert  _et al._ [2013] A. Fert, V. Cros, and J. Sampaio, Nature Nanotechnology 8, 152 (2013).
  * Liu  _et al._ [2017] C. Liu, Y. Zang, W. Ruan, Y. Gong, K. He, X. Ma, Q.-K. Xue, and Y. Wang, Phys. Rev. Lett. 119, 176809 (2017).
  * Wakatsuki  _et al._ [2015] R. Wakatsuki, M. Ezawa, and N. Nagaosa, Scientific Reports 5, 13638 (2015).
  * Liu  _et al._ [2013] J. Liu, W. Duan, and L. Fu, Phys. Rev. B 88, 241303 (2013).
  * Serbyn and Fu [2014] M. Serbyn and L. Fu, Phys. Rev. B 90, 035402 (2014).
  * Okada  _et al._ [2013] Y. Okada _et al._ , Science 341, 1496 (2013).
  * Grifoni and Hänggi [1998] M. Grifoni and P. Hänggi, Physics Reports 304, 229 (1998).
  * Shirley [1965] J. H. Shirley, Phys. Rev. 138, B979 (1965).
  * Busl  _et al._ [2012] M. Busl, G. Platero, and A.-P. Jauho, Phys. Rev. B 85, 155449 (2012).
  * Tsuji  _et al._ [2008] N. Tsuji, T. Oka, and H. Aoki, Phys. Rev. B 78, 235124 (2008).
  * Ke  _et al._ [2020] M. Ke, M. M. Asmar, and W.-K. Tse, Phys. Rev. Research 2, 033228 (2020).
  * Kitagawa  _et al._ [2010] T. Kitagawa, E. Berg, M. Rudner, and E. Demler, Phys. Rev. B 82, 235114 (2010).
  * Ezawa [2013] M. Ezawa, Phys. Rev. Lett. 110, 026603 (2013).
  * Calvo  _et al._ [2011] H. L. Calvo, H. M. Pastawski, S. Roche, and L. E. F. F. Torres, Applied Physics Letters 98, 232103 (2011).
  * Zhou and Wu [2011] Y. Zhou and M. W. Wu, Phys. Rev. B 83, 245436 (2011).
  * Kogan [2013] E. Kogan, Graphene 2, 8 (2013).
  * Sherafati and Satpathy [2011a] M. Sherafati and S. Satpathy, Phys. Rev. B 83, 165425 (2011a).
  * Sherafati and Satpathy [2011b] M. Sherafati and S. Satpathy, Phys. Rev. B 84, 125416 (2011b).
  * Rudner and Lindner [2020] M. S. Rudner and N. H. Lindner, Nature Reviews Physics 2, 229 (2020).
  * Imamura  _et al._ [2004] H. Imamura, P. Bruno, and Y. Utsumi, Phys. Rev. B 69, 121303 (2004).
  * Zhu  _et al._ [2011] J.-J. Zhu, D.-X. Yao, S.-C. Zhang, and K. Chang, Phys. Rev. Lett. 106, 097201 (2011).
  * Shiranzaei  _et al._ [2017] M. Shiranzaei, H. Cheraghchi, and F. Parhizgar, Phys. Rev. B 96, 024413 (2017).
  * Cheng  _et al._ [2016] R. Cheng, S. Okamoto, and D. Xiao, Phys. Rev. Lett. 117, 217202 (2016).
  * Kim  _et al._ [2016] S. K. Kim, H. Ochoa, R. Zarzuela, and Y. Tserkovnyak, Phys. Rev. Lett. 117, 227201 (2016).
  * Chen  _et al._ [2018] L. Chen, J.-H. Chung, B. Gao, T. Chen, M. B. Stone, A. I. Kolesnikov, Q. Huang, and P. Dai, Phys. Rev. X 8, 041028 (2018).
  * Sampaio  _et al._ [2013] J. Sampaio, V. Cros, S. Rohart, A. Thiaville, and A. Fert, Nature Nanotechnology 8, 839 (2013).
  * Asmar and Tse [2019] M. M. Asmar and W.-K. Tse, Phys. Rev. B 100, 014410 (2019).
  * Dugaev  _et al._ [2006] V. K. Dugaev, V. I. Litvinov, and J. Barnas, Phys. Rev. B 74, 224438 (2006).
  * Saremi [2007] S. Saremi, Phys. Rev. B 76, 184430 (2007).
  * Black-Schaffer [2010] A. M. Black-Schaffer, Phys. Rev. B 81, 205416 (2010).
  * Chen  _et al._ [2010] J. Chen _et al._ , Phys. Rev. Lett. 105, 176602 (2010).
  * Yeats  _et al._ [2017] A. L. Yeats, P. J. Mintun, Y. Pan, A. Richardella, B. B. Buckley, N. Samarth, and D. D. Awschalom, Proceedings of the National Academy of Sciences 114, 10379 (2017).
  * Schoenlein  _et al._ [2000] R. W. Schoenlein, S. Chattopadhyay, H. H. W. Chong, T. E. Glover, P. A. Heimann, C. V. Shank, A. A. Zholents, and M. S. Zolotorev, Science 287, 2237 (2000).
  * Luke [1969] Y. Luke,  _The Special Functions and Their Approximations_, ISSN (Elsevier Science, 1969).
  * Fields [1972] J. L. Fields, Mathematics of Computation 26, 757 (1972).
  * Wang  _et al._ [2013b] Y. H. Wang, H. Steinberg, P. Jarillo-Herrero, and N. Gedik, Science 342, 453 (2013b).
  * Nuske  _et al._ [2020] M. Nuske, L. Broers, B. Schulte, G. Jotzu, S. A. Sato, A. Cavalleri, A. Rubio, J. W. McIver, and L. Mathey, Phys. Rev. Research 2, 043408 (2020).
  * McIver  _et al._ [2020] J. W. McIver, B. Schulte, F.-U. Stein, T. Matsuyama, G. Jotzu, G. Meier, and A. Cavalleri, Nature Physics 16, 38 (2020).
  * Hindmarch and Hickey [2003] A. T. Hindmarch and B. J. Hickey, Phys. Rev. Lett. 91, 116601 (2003).
  * Khajetoorians  _et al._ [2012] A. A. Khajetoorians, J. Wiebe, B. Chilian, S. Lounis, S. Blügel, and R. Wiesendanger, Nature Physics 8, 497 (2012).
  * Zhou  _et al._ [2010] L. Zhou, J. Wiebe, S. Lounis, E. Vedmedenko, F. Meier, S. Blügel, P. H. Dederichs, and R. Wiesendanger, Nature Physics 6, 187 (2010).
  * Yarmohammadi  _et al._ [2021] M. Yarmohammadi, C. Meyer, B. Fauseweh, B. Normand, and G. S. Uhrig, Phys. Rev. B 103, 045132 (2021).
  * Hausinger and Grifoni [2010] J. Hausinger and M. Grifoni, Phys. Rev. A 81, 022117 (2010).
  * Park  _et al._ [2018] S. Park, H. Min, E. H. Hwang, and S. Das Sarma, Phys. Rev. B 98, 064425 (2018).
  * Min  _et al._ [2017] H. Min, E. H. Hwang, and S. Das Sarma, Phys. Rev. B 95, 155414 (2017).
  * Prudnikov  _et al._ [1986] A. Prudnikov, I. Brychkov, I. Brychkov, and O. Marichev,  _Integrals and Series: Special functions_, Integrals and Series (Gordon and Breach Science Publishers, 1986).