# Cooperative Motions and Topology-Driven Dynamical Arrest in Prime Knots

Hyo Jung Park  Smith College  Anna Lappala<EMAIL_ADDRESS>Massachusetts General Hospital & Department of Molecular Biology, Harvard
University

Knots are entangled structures that cannot be untangled without a cut.
Topological stability of knots is one of the many examples of their important
properties that can be used in information storage and transfer. Knot dynamics
is important for understanding general principles of entanglement as knots
provide an isolated system where tangles are highly controlled and easily
manipulated. To unravel the dynamics of these entangled topological objects,
the first step is to identify the dominant motions that are uniquely guided by
knot structure and its complexity. We identify and classify motions into three
main groups— orthogonal, aligned, and mixed motions, which often act in
unison, orchestrating the complex dynamics of knots. The balance between these
motions is what creates an identifiable signature for every knot. As knot
complexity increases, the carefully orchestrated dynamics is gradually
silenced, eventually reaching a state of topologically driven dynamical
arrest. Depending on their complexity, knots undergo a transition from nearly
stochastic motions to either non-random or even quasiperiodic dynamics before
culminating in dynamical arrest. Here, we show for the first time that
connectivity alone can lead to a topology-driven dynamical arrest in knots of
high complexity. Unexpectedly, we noticed that some knots undergo cooperative
motions as they reach higher complexity, uniquely modulating conformational
patterns of a given knot. Together, these findings demonstrate a link between
topology and dynamics, presenting applications to nanoscale materials.

##  I Introduction

Knots are formally defined as closed curves embedded in $\mathbb{R}^{3}$ that
cannot be untangled without a cut. Prime knots are a subset of closed curves
that cannot be decomposed, i.e. expressed as a connected sum of two non-
trivial knots. Prime knots contain numerous subcategories, and two of those
are the focus of this manuscript—torus and twist knots (Fig. 1). These knots
have intrinsic symmetries, forming braid-like crossings along a circular path.
The only difference between torus and twist knots is that twist knots also
contain self-links—loop-like structures that are not present in torus knots
(compare Fig. 1(a) and (b)). The notation used throughout the manuscript, such
as $3_{1}$ or $10_{124}$, reflects the crossing number, in this case 3 and 10,
respectively. This is the number of times a knot crosses itself, and a higher
number of crossings indicates higher complexity. Subscript values 1 and 124
are based on tabulated knot configurations, and are used to identify a knot
with the same number of crossings, but different configurations.

Knots have played an important role in the understanding and controlling the
dynamics and mixing of fluids such as liquid crystals [1, 2], plasma [3, 4]
and even quantum fluids [5, 6, 7]. Topological stability of knots is one of
the reasons for the growing interest in the field [8, 9, 10], with
applications in information storage from quantum- to biological scales [11,
12]. Knots provide an excellent opportunity to study the intricate
relationship between topology and dynamics in various physical and biological
systems [13, 14, 15, 16]. It is well known that chain entanglement strongly
affects dynamical, thermal and mechanical properties of polymers, including
biopolymers such as DNA [17, 18, 19, 20]. Knots have been used in nanofluidic
applications and nanopore sequencing, serving as invaluable tools for genomic
analysis [21, 22]. While many physical properties of knots can be elegantly
captured in experiments, the details of the relationship between knotting,
entanglement and dynamics are not well understood.

In this work, Molecular Dynamics simulations were used to investigate the link
between topology and dynamics of torus and twist knots. Detailed analysis of
simulation trajectories revealed that knots undergo ‘coherent’ motions, in
which different parts of the knot move cooperatively producing complex motion
patterns. Two types of cooperative motions were identified—“breathing” and
tangential motions. Principal component analysis allowed us to decompose those
dynamics into fundamental ones which we termed orthogonal, tangential, and
mixed motions. Based on principal component analysis, free energy landscapes
(FEL) and correlation matrices were generated. Together with van Hove
correlation functions, they demonstrated how knot complexity can lead to what
we call “topological dynamical arrest”, whereby knots with a high number of
crossings display frozen-in dynamics purely as a result of topology. Lastly,
autocorrelation functions (ACF) of the trajectories were computed, which
showed that dynamical trends in decreasing knot size—as knots evolve from
completely random dynamics to dynamical arrest—are entirely determined by knot
complexity and topology. An attempt to disentangle how various parameters
govern and give rise to unique knot dynamics, this work presents a link
between the motions, dynamical arrest, knot complexity and topology.

Figure 1: A subset of studied prime knots with respective dominant motion
patterns based on principal component analysis of simulated trajectories
depicted as arrows. (a) Torus knots $3_{1}$, $5_{1}$, $7_{1}$ and $9_{1}$. The
trefoil knot $3_{1}$ is also a twist knot. (b) Twist knots $4_{1}$, $5_{2}$,
$7_{2}$ and $9_{2}$. Arrows represent the eigenvectors of the first principal
component (PC1) derived from simulation trajectories, and the chain colors
represent the amplitude of motion (blue-small to red-large).(c) An example
representation of orthogonal, aligned, and mixed motions from simulation
trajectories.

##  II Results and Discussion

### Breathing, tangential, and basis motions

Polymer chains, entangled or not, experience different types of motions as a
result of self-interactions as well as interactions with the solvent. Here,
knotted polymer structures are studied, and to identify motions unique to
knots of a given subgroup and complexity, principal component analysis was
performed. Principal components (PC) represent dominant motions (fluctuations)
extracted from simulation trajectories. The strongest of those—PC1 and PC2—can
be classified into specific types of motion and used to characterize the
overall dynamics of a knot. For this classification, the PC eigenvectors were
used as fluctuation vectors (FV) to construct correlation matrices, which
display correlations between the root-mean-squared (RMS) fluctuations of
individual particles of a knot. The appearance of correlation matrices is
closely related to the behavior of FVs. Examining those matrices thus allows
comparison of the dynamics represented by different PCs. As a result, three
general types of motions were identified—orthogonal, aligned, and mixed (Fig.
1)—which underlie the aforementioned breathing and tangential categories of
cooperative motions.

Orthogonal motions are dynamics in which the FVs of individual particles are
largely perpendicular to the knot strand. These motions are generally
characterized by square patterns in correlation maps as seen in that of
$4_{1}$ (Fig. 2(c,ii)). Such maps show even-sized boxes with alternating blue
and red colors, which correspond to positive and negative correlations,
respectively. The square patterns may also appear to be smoother as in the
correlation map of $3_{1}$ (Fig. 2(a,i)). Since orthogonal motions generate
in- and outward fluctuations of the knot strand, we argue that orthogonal
dynamics account for breathing motion, which is the recursive tightening and
loosening of the crossing regions in a knot.

Aligned motions, on the other hand, describe motions in which the FVs of
adjacent particles point in the same direction, that is, are “aligned” with
each other. These motions in general coexist with smooth patterns in
correlation maps where the same color continues in the oblique direction (see
$5_{1}$ in Fig. 2(a,ii)). In low complexity knots, if the geometry permits,
aligned and orthogonal dynamics cooperate in propelling tangential motion
along the knot, causing the knot strands to smoothly slither or pass by each
other. It is aligned motions that give directionality to tangential dynamics
while orthogonal motions open up tighter parts of the knot to allow the
slithering of the knot strands. Lastly, mixed motions, which are neither
orthogonal nor aligned (or a combination of both), are shown in the PC3 FV
plot of $5_{1}$ (Fig. 1(c)).

Often, these three types of general motions are present in knots with varying
intensities. The torus knot $3_{1}$, for instance, is dominated by orthogonal
(breathing) motions in its PC1 dynamics, and yet in PC3, the varying
orientations of the aligned FVs match the geometry of the knot, describing the
knot’s tangential dynamics (Fig. 1(c)). This contrasts with $5_{1}$, in which
the hierarchy of motions is reversed, with PC1 describing tangential motions
(again, with a proper geometry) and PC3 the breathing. In both $3_{1}$ and
$5_{1}$, mixed motions appear in fifth or higher PCs. Different compositions
of these basis motions create complex patterns of dynamics, giving each knot a
unique dynamical signature.

Figure 2: Analysis of torus and twist knot correlation matrices and their
respective free energy contour profiles. (a) correlation matrices of torus
knots i) $3_{1}$ ii) $5_{1}$ iii) $7_{1}$, iv) $8_{19}$, v)$9_{1}$, vi)
$10_{124}$, showing correlated motion patterns in red and anti-correlated
motions in blue. (b) Free energy profiles of torus knots derived from
principal component analysis of simulation trajectories. i) $3_{1}$ ii)
$5_{1}$ iii) $7_{1}$, iv) $8_{19}$, v)$9_{1}$, vi) $10_{124}$. Blue regions
represent favourable regions of the phase space that are readily explored
within a simulation trajectory and red regions remain unexplored. (c)
correlation matrices of twist knots i) $3_{1}$ ii) $4_{1}$ iii) $5_{2}$, iv)
$6_{1}$, v)$7_{2}$, vi) $8_{1}$, showing correlated motion patterns in red and
anti-correlated motions in blue. (d) Free energy profiles of twist knots
derived from principal component analysis of simulation trajectories. i)
$3_{1}$ ii) $4_{1}$ iii) $5_{2}$, iv) $6_{1}$, v)$7_{2}$, vi) $8_{1}$. Blue
regions represent favourable regions of the phase space that are readily
explored within a simulation trajectory and red regions remain unexplored.

### Comparison of torus and twist knots

Torus and twist knots share common characteristics in their dynamics. First,
the tighter parts of the knot tend to have tangential motions, while the
looser parts tend to exhibit breathing motions. As evident in FV plots (see
Fig. S2), they become more tangential to the strand in the vicinity of
crossings, and more perpendicular near the loops or self-links of the knots.

As knot complexity increases, knot’s motions become more confined due to their
complex topology, eventually culminating in dynamical arrest in both types of
knots. The higher the crossing number, the narrower the explored basins of
FELs, and highly complex knots such as $8_{1}$ and $10_{124}$ undergo a
transition into a state where the structure is too tightly intertwined to
experience any motion at all.

In general, torus knots (Fig. 2(b)) have highly symmetric, circular FELs,
while twist knots (Fig. 2(d))—except $3_{1}$, which is both torus and
twist—have asymmetric, broad and level ones. In addition, unlike twist knots,
whose correlation maps along the first two PCs show relatively clear square
patterns, those of torus knots lack such patterns (compare Fig. 2(a) and (c)).

Torus and twist knots also exhibit different trends as knot complexity
increases: For twist knots, the explored areas of the FELs remain relatively
constant, as well as the sizes of “squares” in the correlation maps (see Fig.
2(c) and (d)). This is in contrast to what is observed in torus knots, in
which the areas of FELs and the patterns of correlation maps vary
significantly. It is apparent that the presence of self-links in twist knots
is the source of these different behaviors, demonstrating how even a
relatively minor change in knot architecture can lead to considerable changes
in their dynamics.

Figure 3: Conformational analysis of knots of varying complexity provides
insights into details of their structure and dynamics. (a) The number of
accessible conformational states derived from free energy profiles for
persistence lengths $(l_{p})$ 1 and 50$\sigma$ for knots of varying complexity
shows a complex nonlinear relationship. (b) Average curvature / crossing
number as a function of complexity follows a power law relationship. (c) Loop
size (Number of monomers / number of crossings) as a function of curvature per
loop shows a linear relationship. (d) Number of conformations as a function
loop size follows a linear relationship. (e) Length-scale grows exponentially
with loop size. (f) The number of clusters derived from simulation
trajectories as a function of loop size follows a linear curve.

### The link between knot complexity, motion lengthscales, and dynamical
arrest

Here, the focus is on the connection between structure in general as well as
specific structural elements such as loops, and dynamics, and how both respond
to topological changes, such as knot complexity. It is evident that knots with
a high number of crossings (higher complexity) tend to also be denser– but are
local densities the only factor that contributes to the inevitable dynamical
arrest in highly complex knots? The simple answer is it is not– and while high
density does strongly contribute to dynamical arrest, the extent of it also
depends on the size of loops, chain flexibility and local geometry. In
general, the number of accessible conformations reduces non-linearly as knot
complexity increases. (Fig. 3(a), and as expected, the number of accessible
states is decreased with an increase in the persistence length of the chain. A
careful look into trajectory-averaged curvature per loop as a function of
complexity follows a power law, suggesting that highly curved loop segments
are mostly present in low complexity knots. In fact, loop size is linearly
dependent on curvature per loop, confirming that on average,larger knots have
highly curved loops. (Figs. 3(b, c). Perhaps unsurprisingly, this leads to a
linear relationship between the number of energetically accessible
conformations and loop size. (Fig. 3(d) Fourier transforms of correlation maps
were examined in order to study how the correlated motions evolve as a
function of knot complexity. For each PC1 correlation map, the largest
wavelength (length-scale) of blue-to-red inversion was calculated from the
map’s Fourier transform (Fig. 3(e)). As expected in both types of knots, a
general decrease in lengthscale was observed as the crossing number increased.
This trend stems from the fixed size of knots (all knots are $N=31$ monomers),
which results in decreased spacing between crossing regions, or “length of
loops,” making the correlation patterns more granular towards higher knot
complexity. Twist knots have constant lengthscales from $4_{1}$ to $7_{2}$, in
agreement with the previous observation that the correlation maps of those
four knots have nearly identical sizes of squares. A notable feature is the
sudden drop in lengthscale after $7_{2}$, which contrasts the relatively
smooth decrease in lengthscale in torus knots. This reflects the sharp
transition into dynamical arrest in twist knots that is observed in their FELs
(Fig. 2(d)). Evidently, all these unique characteristics of twist knots—the
sharp drop and the flat region in lengthscale as well as the nearly identical
correlation maps of the four twist knots—stem from the presence of self-links,
which may sustain the variety or strength of dynamics against the decreasing
length of loops or give rise to similar dynamics in the four twist knots of
different complexities. Interestingly, lengthscale scales exponentially with
loop size (Fig. 3(e)).

Cluster analysis of simulation trajectories also illustrate these interesting
trends: the number of dynamical clusters decreases as a function of knot
complexity, with torus knots following a linear decline and twist knots mostly
only having two identifiable clusters. In this context, clusters signify types
of observed motions derived by grouping similar basins in the FELs (see
Section IV for details on cluster analysis). Like the number of conformational
states, number of clusters scales linearly with loop size.

All the data presented so far reveal the key characteristics of dynamical
arrest in high-complexity knots: highly restrictive free energy profiles and
small correlation lengthscales. In order to verify if those are indeed marks
of dynamical arrest, van Hove self-correlation functions were computed for all
knots, and it is clear that the RMS profiles for dynamically arrested knots
such as $8_{1}$, $9_{2}$, and $10_{124}$ are stable over long times, as shown
in Fig. 4(f), 5 and Fig. S2. These findings suggest that topology alone can
lead to dynamical arrest in knotted structures.

Figure 4: Van Hove self-correlation function for torus knots at different time
frames of the simulation trajectory. (a) $3_{1}$ (b) $5_{1}$ (c) $7_{1}$, (d)
$8_{19}$, (e)$9_{1}$, (f) $10_{124}$. Figure 5: A phase diagram representing a
relationship between the number of energetically feasible conformational
states, crossing number (complexity) and persistence length. Three major
phases are observed: a liquid-like ergodic phase (red), a glass-like
dynamically arrested state (blue) and an intermediate state that involves
cooperative motions.

### Controlling randomness: topology-dependent dynamical trends in varying
knot size

Figure 6: Decoupling density-induced dynamical arrest from geometry-induced
one. (a) Spherically confined knots $3_{1}$ and $7_{6}$ showed higher RMS
displacements as opposed to $8_{1}$ and $10_{124}$. Knot $7_{6}$ showed a
unique profile that corresponds to different conformational states. (b) A
histogram based on RMSD values from (a) for $3_{1}$, showing at least 2 peaks.
(c) As in (b), $7_{6}$ showing at least 4 peaks. (d) As in (b), $8_{1}$,
showing at least 2 peaks. (e) As in (b), $10_{124}$, showing only one peak.(f)
Autocorrelation functions of some knots and the dynamical trends observed in
varying knot size. Collective autocorrelation functions of $3_{1}$ with size
$N=11$, of the $6_{1}$ knot with size $N=27$, and the $7_{6}$ knot with size
$N=31$. Each is an example illustrating complete randomness, non-randomness,
and quasiperiodicity, respectively. (g) A schematic of the transition from
complete randomness to dynamical arrest as knot size decreases. $3_{1}$ does
not lose randomness until it reaches a state of dynamical arrest (Route Red:
a). For knots with crossing number of $>4$, randomness gradually decreases as
knot size shrinks and non-random patterns appear. Only a subset of high-
complexity knots exhibit quasiperiodicity in their collective ACFs.

It may be tempting to assume that dynamical arrest in higher complexity knots
is simply triggered by higher contact densities. In order to decouple the
effects of local density from geometry, the structures were spherically
confined to reach the same effective volume, demonstrating that knots of
varying complexities display different dynamics as a result of geometry alone.
To confirm this, a number of structures ranging from a linear (unknotted)
chain to a knot with 10 crossings were evaluated under spherical confinement
to control for density, with the idea that if density was the only factor
influencing the dynamics of knots, root mean squared displacements would be
the same for all the different knot topologies. Yet our results (Fig. 6)
demonstrate that this is not the case, and knots of higher complexity display
a dramatic decrease in mobility, confirming the hypothesis that topology plays
a critical role in the initiation of this new kind of purely topology-driven
dynamical arrest. Autocorrelation analysis of knot trajectories showed
interesting dynamical trends as knot size was varied. Autocorrelation
functions (ACF) capture deviations from standard random dynamics, and for each
knot, both the autocorrelation of the collective structure (“collective ACF,”
which takes into account all the $3N$ coordinates in the knot) and the
autocorrelation of individual coordinates (“individual ACFs”) were computed
(see Section. IV for more detail). The results revealed that randomness in
knot dynamics is highly controllable by varying knot size. When the size of a
knot is large enough that the entire structure is loose, both the knot’s
individual and collective ACFs are flat with values nearly equal to zero,
implying complete randomness in the knot’s dynamics (see Fig. S4—ACFs of
$3_{1}$ $N=31$ and some large high-complexity knot). In the lowest-complexity
knot ($3_{1}$), randomness is sustained regardless of knot size (see $3_{1}$
in Fig. 6). However, knots with crossing number of 4 or higher begin to
diverge from complete randomness as their size decreases and show significant
non-random patterns in their trajectories. In some of these higher-complexity
knots, periodic trends emerge at some particular knot sizes, which we termed
“quasiperiodicity” Examples of this general non-randomness and
quasiperiodicity are show in Fig. 6 (a) ($6_{1}$ and $7_{6}$) (see S5 for more
examples). Quasiperiodicity (Fig. 6 (g, Route c) is a subset of general non-
randomness (Route B) and characterized by a vanishing semi-sinusoidal curve in
the knots’ collective ACFs (see Fig. 6 (f) $7_{6}$). Interestingly, it only
emerges in some torus knots $7_{1}$, $8_{19}$, and $10_{124}$ and in certain
high-complexity knots such as $7_{2}$, $7_{4}$, and $7_{6}$, which are neither
torus nor twist. Upon closer examination, it was noticed that those knots
generally have ring-shaped FELs (see Fig. 2(b) and Fig. S6 (FELs of “other”
knots)). A striking feature of these results is that the “fate” of dynamical
trends in varying knot size is solely determined by knot complexity and
topology: knots must have an average crossing number of at least 4 or 5 in
order to exhibit non-randomness. The type of non-randomness depends on
topological structures, with only a subset of high-complexity knots allowed to
have quasiperiodic behaviors. This is another demonstration that topology
alone can affect the dynamics of entangled objects.

##  III Conclusion

In this work, the dynamics of prime knots was studied, focusing on torus- and
twist knots due to their intrinsically symmetrical structures. Results of this
work revealed a strong connection between knot dynamics, topology and knot
complexity. Dominant motions in torus and twist knots were studied in detail,
and three types of motion patterns were identified: “breathing”, aligned, and
a combination of the two. These motions work in unison to generate complex
patterns in the dynamics of individual knots. It was also found that
randomness in knot dynamics can be modulated by varying knot size. As knot
size shrinks, knots exhibit either complete randomness, non-randomness, or
quasiperiodicity depending on knot complexity and topology, before reaching a
state of dynamical arrest. It was observed that some knots undergo cooperative
motions as they reach higher complexity, uniquely modulating conformational
patterns of a given knot. Intriguingly, knots of high complexity were found to
undergo a so called topology-driven dynamical arrest that manifests in a
dramatic reduction of explored configuration space accompanied by a shift in
typical patterns of motion in knots. Understanding periodic and non-random
dynamics as well as topology-guided dynamical arrest will enable future
studies of entangled systems not limited to knots, leading to the control and
regulation of local dynamics of entangled structures, including biological
systems such as proteins, RNA and chromosomes.

##  IV Methods

Initial configurations of prime knots (Fig. 1(a)) were generated based on
models from KnotPlot software [23]. Coordinates were interpolated using three-
dimensional parametric curves in order to build knots of desired size. Each
curve was scaled to fit 31 particles (a number that produces stable
simulations for all studied knots), and along the curves knots were
constructed particle by particle, ensuring that all bond lengths are within
$\pm 0.08\sigma$ and the inter-atomic distances between non-bonded pairs at
least $0.72\sigma$, where $\sigma$ is the monomer radius in Lennard-Jones
units. Molecular Dynamics simulations described in this work were performed
using LAMMPS software [24]. From the simulation trajectories, for each knot
the data matrix $X$ of shape [$3N,m$] was obtained, where $N$=31, and $m$ is
the number of frames in the trajectory. The transitions and rotations along
the trajectory were removed before $X$ was constructed, as the focus was on
the internal dynamics of the knots [25]. Principal component analysis of the
data matrix was performed to extract the most dominant fluctuations from
Molecular Dynamics (MD) trajectories. The correlation matrix was first
calculated

$\displaystyle Q=\frac{1}{m}\overline{X}\overline{X}^{T},$ (1)

where $\overline{X}$ is the mean-centered data matrix obtained by subtracting
means from the rows of $X$. Eigendecomposition of $Q$

$\displaystyle Q=V\Lambda V^{-1}$ (2)

yields the orthonormal eigenvectors in the columns of $V$ and the eigenvalues
in the diagonal of $\Lambda$. The principal components (PC) are these
eigenvectors $\bm{v_{\mathrm{i}}}$ and each describe the direction of a
particular deviation from the average conformation. Each of their
corresponding eigenvalues $\lambda_{i}$ indicates the portion of the total
variance of $X$ that is explained by that PC. A large eigenvalue means that
the direction of that specific deviation is frequently explored over the
entire trajectory. Thus, by sorting the PCs in a descending order of their
eigenvalues and analyzing the first few dominant PCs, one can identify the
strongest motions that are present in the knots. In this work the PCs are
referred to as PC1 (the first PC), PC2 (the second), and so forth. Projections
along principal components, free energy landscapes, correlation maps, and
fluctuation vector plots were examined to characterize the motions of the
knots. Cluster analysis, van Hove self-correlation functions, and
autocorrelation functions were used in studying the relationship between those
motions, knot complexity, and dynamical arrest. Fig. S2 shows the projections
of conformations on the plane formed by two PCs. Each projection contains $m$
scores (shown as dots), and each score provides information of how much the
two PCs account for that particular conformation. Values of $z_{ij}$ are
obtained by taking the dot product of the $i$th frame and $j$th eigenvector:

$\displaystyle z_{\mathrm{ij}}=\bm{x_{\mathrm{i}}}\cdot\bm{v_{\mathrm{j}}},$
(3)

where $\bm{x_{\mathrm{i}}}$ is the $i$th column of $\overline{X}$, and
$\bm{v_{\mathrm{j}}}$ the $j$th column of $V$. The calculation of scores can
be easily achieved via matrix algebra:

$\displaystyle Z=\overline{X}^{T}V.$ (4)

In this case, the scores are the elements of the score matrix $Z$. Free energy
landscapes (FEL) were computed from the distributions of projections along
principal components based on the formula

$\displaystyle\Delta G(r_{1},r_{2})=-k_{\mathrm{B}}T\ln{P(r_{1},r_{2})},$ (5)

where $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ the temperature of the
simulation, and $P(r_{1},r_{2})$ the estimated probability density function
along the axes $r_{1}$ and $r_{2}$ of some two PCs. $P(r_{1},r_{2})$ was
estimated (Fig. 1(c)) by fitting a score plot using a two-dimensional Gaussian
kernel function with the improved Sheather Jones bandwidth [26]. After
obtaining $P(r_{1},r_{2})$, one can plot the three-dimensional free energy
landscapes and also their free energy contour maps. Correlation maps display
the correlations between the root-mean-squared (RMS) fluctuations of
individual particles in a knot. The correlation between particles $i$ and $j$
is defined as

$\displaystyle
c_{\mathrm{ij}}=\frac{<\Delta\bm{R_{\mathrm{i}}}\cdot\Delta\bm{R_{\mathrm{j}}}>}{\sqrt{<\Delta\bm{R_{\mathrm{i}}}^{2}>\cdot<\Delta\bm{R_{\mathrm{j}}}^{2}>}},$
(6)

where $\Delta\bm{R_{\mathrm{i}}}$ and $\Delta\bm{R_{\mathrm{j}}}$ are the
fluctuation vectors. These correlations $c_{ij}$ comprise the $N\times N$
correlation map $C$. Since the directions of dominant conformational changes
are described by the PCs, one can use the PC eigenvectors as the fluctuation
vectors $\Delta\bm{R}$ of the particles. Thus, the correlation map $C$ of each
PC can be found, and for this calculation the algorithm developed by Eyal et
al.[27] was used. Visualizing these correlations as heat maps provides a means
to examine and compare the specific correlated motions described by different
PCs. Fluctuation vector (FV) plots show the PC eigenvectors superimposed onto
initial conformations. By examining how the fluctuation vectors are aligned
with the knot strands, one can determine the types of dynamics the PCs
represent, such as orthogonal, aligned, and mixed motions. One can also
deconstruct each PC to study its contribution to the overall cooperative
dynamics—tangential and breathing. This can be achieved by projecting the PC
eigenvectors onto the Frenet frame $(\mathbf{T},\mathbf{N},\mathbf{B})$ along
the knot, where the tangential vector $\mathbf{T}$ describes tangential
dynamics, and the normal $\mathbf{N}$ and binormal $\mathbf{B}$ vectors
breathing dynamics. To reproduce the overall motions along a particular Frenet
vector, one can use equations of the form

$\displaystyle\bm{v}_{\mathrm{fluc}}=\sum_{i}\bm{v_{\mathrm{i}}}\cdot\bm{\mathrm{F}}\sqrt{\lambda_{\mathrm{i}}}.$
(7)

Here, $\sqrt{\lambda_{\mathrm{i}}}$ —the square root of the $i$th
eigenvalue—describes the standard deviation along the corresponding PC axis,
whose direction is given by the eigenvector $\bm{v_{\mathrm{i}}}$. The Frenet
vector $\mathbf{F}$ should be replaced by either $\mathbf{T}$, $\mathbf{N}$,
or $\mathbf{B}$. This Frenet method is especially useful if cooperative
motions depend on more than one basis motion types (e.g. tangential dynamics
are due to both orthogonal and aligned motions). The correlation maps for
breathing dynamics obtained from this method do generally exhibit square
patterns, in agreement with the results obtained from individual PCs (Fig.
2(c)). Such comparison, however, cannot be made between the correlation maps
of aligned and tangential dynamics, as the former dynamics do not always match
the latter. Cluster analysis was used to group motions of a knot into clusters
based on their similarities. We clustered frames in the simulation trajectory,
i.e. the columns of the mean-centered data matrix $\overline{X}$, by using the
average-linkage algorithm. We specifically chose this algorithm among other
options such as complete linkage or K-means, since average linkage produces
more consistent results for MD trajectories. For the clustering criterion,
correlations between trajectory frames were used, and for determining the
number of clusters, the so-called “elbow method” was employed. Clusters can be
plotted on the PC1-PC2 plane to visualize their locations in phase space,
which is achieved by coloring scores in the score plot based on cluster
indices. The number of identified clusters serves as a measure of the
diversity of motions in a knot, allowing the study of how a certain parameter
such as knot size and bond stiffness affects the knot’s dynamics. In this
work, this number was used to demonstrate the link between knot complexity and
dynamical arrest: knots become tighter as the crossing number increases, which
confines knot dynamics to fewer conformations. This is exactly reflected in
their lower numbers of identified clusters, signifying less variety in their
motions. The van Hove correlation function gives the average probability
density of finding a particle at position $\mathbf{r^{\prime}}+\mathbf{r}$ at
time $t^{\prime}+t$, given that there was a particle at position
$\mathbf{r^{\prime}}$ at time $t^{\prime}$. In other words, it describes the
correlation between the “presence” of a particle in some position and the
presence of another (or the same) particle in some position at some reference
time [28]. Normally the reference time $t^{\prime}$ is taken to be zero. For a
system of $N$ particles, the van Hove correlation function can be written as

$\displaystyle
G(\mathbf{r},t)=\frac{1}{N}\langle\sum_{i=1}^{N}\sum_{j=1}^{N}\delta^{(3)}(\mathbf{r}-\mathbf{r_{i}}(t)+\mathbf{r_{j}}(0))\rangle,$
(8)

where $\mathbf{r_{i}}$ and $\mathbf{r_{j}}$ represent particle positions,
$\delta^{(3)}(\cdot)$ a three-dimensional Dirac delta function, and
$\langle\cdot\rangle$ an ensemble average (an average over identical systems).
At time $t=0$, the function reduces to

$\displaystyle G(\mathbf{r},0)$
$\displaystyle=\frac{1}{N}\langle\sum_{i=1}^{N}\sum_{j=1}^{N}\delta^{(3)}(\mathbf{r}-\mathbf{r_{i}}(0)+\mathbf{r_{j}}(0))\rangle$
(9) $\displaystyle=\delta^{(3)}(\mathbf{r})+\frac{1}{N}\sum_{i\neq
j}\langle\delta^{(3)}(\mathbf{r}-\mathbf{r_{i}}(0)+\mathbf{r_{j}}(0))\rangle$
(10) $\displaystyle=\delta^{(3)}(\mathbf{r})+g(\mathbf{r}),$ (11)

where $g(\mathbf{r})$ is the usual radial distribution function [28]. In the
case of systems composed of distinguishable particles, the function $G$ can be
separated into two parts:

$\displaystyle G(\mathbf{r},t)$
$\displaystyle=\frac{1}{N}\langle\sum_{i=1}^{N}\delta^{(3)}(\mathbf{r}-\mathbf{r_{i}}(t)+\mathbf{r_{i}}(0))\rangle+\frac{1}{N}\langle\sum_{i\neq
j}^{N}\delta^{(3)}(\mathbf{r}-\mathbf{r_{i}}(t)+\mathbf{r_{j}}(0))\rangle$
(12) $\displaystyle:=G_{s}(\mathbf{r},t)+G_{d}(\mathbf{r},t).$ (13)

The “self” part $G_{s}(\mathbf{r},t)$ gives the distribution of displacements
of each particle from their initial position over time, and the “distinct”
part $G_{d}(\mathbf{r},t)$ describes the distribution of relative
displacements between two different particles at different times. Both $G_{s}$
and $G_{d}$ can be calculated from the histograms of corresponding
displacements obtained from MD trajectories. In this work, Gaussian kernel
density functions with improved Sheather Jones bandwidth were used to fit the
histograms and approximate smooth $G_{s}$ and $G_{d}$ functions [26]. As
discussed in Section II, the self part $G_{s}$ (Fig. 4) demonstrated dynamical
arrest in high-complexity knots.

The autocorrelation function (ACF) is a function of delay or lag $\tau$ that
gives the correlation between a time series $X$ and a delayed copy of itself:

$\displaystyle
R(\tau)=\frac{\sum_{t=1}^{m-\tau}(X(t)-\overline{X})(X(t+\tau)-\overline{X})}{\sum_{t=1}^{m}(X(t)-\overline{X})^{2}},$
(14)

where $m$ is the total time over which $X$ was recorded (which in our case is
the number of frames in a MD trajectory), and $\overline{X}$ is the time-
average of $X$. In essence, the ACF is a measure of the self-similarity of a
data set over different delay times $\tau$ [29]. The magnitude and shape of
ACF allow one to infer the dynamical trend of the data, i.e. whether the data
is random, periodic, autoregressive, etc. An essentially random time series
like white noise has an ACF that is equal to zero for all lags $\tau>0$. For
this reason, the standard error of the white noise is used to determine
whether a given time series $X$ has statistically significant autocorrelation
values [29]. $X$ is said to have a significant autocorrelation value at lag
$\tau$ if $R(\tau)$ lies outside the confidence interval given by

$\displaystyle\pm\frac{\sqrt{2}z_{p}}{\sqrt{m}},$ (15)

where $z_{p}$ is the cumulative distribution function of the standard normal
distribution (also known as the error function), and $p=0.95$ ($95\%$) the
confidence level, i.e. the probability with which a new normal random variable
will lie within the interval. Autocorrelation values that fall outside this
interval suggest non-randomness in the data. In this work, two kinds of ACFs
were computed for each knot, which we named a “collective” ACF and
“individual” ACFs. A collective ACF is autocorrelation that takes into account
all $3N$ coordinates in a knot (each coordinate serves as a time series). The
collective ACF is useful in that it encodes the dynamical behavior of the
entire knot; however, as it is the ACF of a multivariate data set, it is
tricky to determine the confidence interval for the data. Thus, individual
ACFs (the ACFs of individual coordinates) were computed in order to verify if
an interesting trend shown in the collective ACF is indeed statistically
significant. Also, for completeness, ACFs of the other prime knots were
computed (not only torus or twist knots), such as $5_{2}$, $6_{3}$, $7_{2}$,
etc.

## References

  * Smalyukh [2020] I. I. Smalyukh, Reports on Progress in Physics 83, 106601 (2020).
  * Alexander [2019] G. P. Alexander, Science 365, 1377 (2019).
  * Smiet  _et al._ [2017] C. B. Smiet, A. Thompson, P. Bouwmeester, and D. Bouwmeester, New Journal of Physics 19, 023046 (2017).
  * Smiet  _et al._ [2019] C. B. Smiet, H. J. de Blank, T. A. de Jong, D. N. L. Kok, and D. Bouwmeester, Journal of Plasma Physics 85, 10.1017/s0022377818001290 (2019).
  * Hall  _et al._ [2016] D. S. Hall, M. W. Ray, K. Tiurev, E. Ruokokoski, A. H. Gheorghe, and M. Möttönen, Nature Physics 12, 478 (2016).
  * Kleckner  _et al._ [2016] D. Kleckner, L. H. Kauffman, and W. T. M. Irvine, Nature Physics 12, 650 (2016).
  * Barenghi [2007] C. F. Barenghi, Milan Journal of Mathematics 75, 177 (2007).
  * Oberti and Ricca [2016] C. Oberti and R. L. Ricca, Journal of Knot Theory and Its Ramifications 25, 1650036 (2016).
  * Schaufelberger [2020] F. Schaufelberger, Communications Chemistry 3, 10.1038/s42004-020-00433-7 (2020).
  * Fielden  _et al._ [2017] S. D. P. Fielden, D. A. Leigh, and S. L. Woltering, Angewandte Chemie International Edition 56, 11166 (2017).
  * Bonesteel  _et al._ [2005] N. E. Bonesteel, L. Hormozi, G. Zikos, and S. H. Simon, Physical Review Letters 95, 10.1103/physrevlett.95.140503 (2005).
  * Lim and Jackson [2015] N. C. H. Lim and S. E. Jackson, Journal of Physics: Condensed Matter 27, 354101 (2015).
  * Ghrist [1998] R. W. Ghrist, Chaos, Solitons & Fractals 9, 583 (1998).
  * Boi [2005] L. Boi, in  _Geometries of Nature, Living Systems and Human Cognition_ (World Scientific, 2005) pp. 203–278.
  * Michieletto  _et al._ [2022] D. Michieletto, Y. A. G. Fosado, E. Melas, M. Baiesi, L. Tubiana, and E. Orlandini, Nucleic Acids Research 50, 4659 (2022).
  * Smrek  _et al._ [2021] J. Smrek, J. Garamella, R. Robertson-Anderson, and D. Michieletto, Science Advances 7, 10.1126/sciadv.abf9260 (2021).
  * Michieletto  _et al._ [2014] D. Michieletto, D. Marenduzzo, E. Orlandini, G. P. Alexander, and M. S. Turner, Soft Matter 10, 5936 (2014).
  * Michieletto  _et al._ [2017] D. Michieletto, N. Nahali, and A. Rosa, Physical Review Letters 119, 10.1103/physrevlett.119.197801 (2017).
  * Signorini  _et al._ [2021] L. F. Signorini, C. Perego, and R. Potestio, The Journal of Chemical Physics 155, 115101 (2021).
  * Rosa  _et al._ [2012] A. Rosa, M. D. Ventra, and C. Micheletti, Physical Review Letters 109, 10.1103/physrevlett.109.118301 (2012).
  * Klotz  _et al._ [2017] A. R. Klotz, V. Narsimhan, B. W. Soh, and P. S. Doyle, Macromolecules 50, 4074 (2017).
  * Amin  _et al._ [2018] S. Amin, A. Khorshid, L. Zeng, P. Zimny, and W. Reisner, Nature Communications 9, 10.1038/s41467-018-03901-w (2018).
  * Scharein [1998] R. G. Scharein,  _Interactive Topological Drawing_ , Ph.D. thesis, Department of Computer Science, The University of British Columbia (1998).
  * Plimpton  _et al._ [2020] S. Plimpton, A. Kohlmeyer, A. Thompson, S. Moore, and R. Berger, LAMMPS stable release 29 October 2020 (2020).
  * David and Jacobs [2014] C. C. David and D. J. Jacobs, Principal component analysis: A method for determining the essential dynamics of proteins, in  _Protein Dynamics_, Vol. 1084, edited by D. R. Livesay (Humana Press, 2014) p. 193–226.
  * Odland [2018] T. Odland, tommyod/kdepy: Kernel density estimation in python (2018).
  * Eyal  _et al._ [2015] E. Eyal, G. Lum, and I. Bahar, Bioinformatics 31, 1487 (2015).
  * Van Hove [1954] L. Van Hove, Phys. Rev. 95, 249 (1954).
  * Vilela  _et al._ [2013] M. Vilela, N. Halidi, S. Besson, H. Elliott, K. Hahn, J. Tytell, and G. Danuser, in  _Fluorescence Fluctuation Spectroscopy (FFS), Part B_, Methods in Enzymology, Vol. 519, edited by S. Y. Tetin (Academic Press, 2013) pp. 253–276.


# Unraveling the dynamics of magnetization in topological insulator-
ferromagnet heterostructures via spin-orbit torque

Taekoo Oh  [<EMAIL_ADDRESS>Naoto Nagaosa  [<EMAIL_ADDRESS>
###### Abstract

Spin-orbit coupling stands as a pivotal determinant in the realm of condensed
matter physics. In recent, its profound influence on spin dynamics opens up a
captivating arena with promising applications. Notably, the topological
insulator-ferromagnet heterostructure has been recognized for inducing spin
dynamics through applied current, driven by spin-orbit torque. Building upon
recent observations revealing spin flip signals within this heterostructure,
our study elucidates the conditions governing spin flips by studying the
magnetization dynamics. We establish that the interplay between spin-
anisotropy and spin-orbit torque plays a crucial role in shaping the physics
of magnetization dynamics within the heterostructure. Furthermore, we
categorize various modes of magnetization dynamics, constructing a
comprehensive phase diagram across distinct energy scales, damping constants,
and applied frequencies. This research not only offers insights into
controlling spin direction but also charts a new pathway to the practical
application of spin-orbit coupled systems.

###### keywords:

American Chemical Society, LaTeX

RIKEN]RIKEN Center for Emergent Matter Science (CEMS), Wako, Saitama 351-0198,
Japan RIKEN]RIKEN Center for Emergent Matter Science (CEMS), Wako, Saitama
351-0198, Japan IR,NMR,UV

##  1 Keywords

Spintronics, Magnetization dynamics, Topological insulator, Spin-orbit torque.

##  2 Introduction

Spin-orbit coupling (SOC), acknowledged as one of the fundamental interactions
in materials 1, reveals a fertile landscape within condensed matter physics. A
standout illustration of its impact is evident in the burgeoning field of
spintronics, where SOC-driven spin dynamics unveils a promising horizon for
practical applications. Accordingly, recent attention has been drawn to
topological insulator-ferromagnet (TI-FM) heterostructures, as evinced by
notable reports 2, 3, 4, 5, 6, 7. The surge in interest stems from their
remarkable efficiency in catalyzing spin dynamics, primarily attributed to the
spin-orbit torque (SOT) induced by an applied current 8, 9.

These heterostructures exhibit intriguing phenomena that closely follow its
intricate dynamics of spins. For instance, an emergent inductance from the
motive force in spiral magnets 10, 11, 12 was extended to the TI-FM
heterostructures 13, 14, 15. Nonreciprocal transport phenomena associated with
the SOC was explored 16, 17.The significant increase of Curie temperature was
reported as well 18, 19, 20, 21. Notably, it has been acknowledged that the
SOT at the TI-FM interface is robust enough to flip magnetization and induce
the sign change of Hall Effect 22, 23, 24, 25, 26. Such properties hold
promise for applications in spintronics devices, particularly in utilizing
topological SOTs for magnetic memories 27, 28, 29, 30, 31, 32..

Motivated by such intriguing phenomena due to the SOC and SOT, this study
explores magnetization dynamics in TI-FM heterostructures with an applied
current, aiming to unveil the conditions for magnetization flip. Considering
both spin-anisotropy and SOT, we establish the equilibrium state of
magnetization and develop a model describing its dynamics under direct current
(DC) or alternating current (AC). In the absence of damping, we identify
oscillating, faltering, and flipping modes for DC, with the latest inducing
magnetization flip. The choice between modes is determined by the relative
strengths of SOT to spin-anisotropy. With the introduction of damping, the
occurrence of spin flip hinges on the duration elapsed until the crossover
from flipping to oscillating modes takes place.

For AC, on the other hand, we discern three final states—adiabatic,
resonating, and chaotic. We classify five distinct modes in the adiabatic
state at low frequencies, which is derived from DC modes. By providing the
phase diagrams of adiabatic modes, we illustrate that the modes are determined
by SOT, spin-anisotropy, driving frequency, and initial driving phase. Lastly,
we explore their transition to resonating and chaotic states at higher
frequencies in the viewpoint of the Fourier transform, revealing that the
periodic array of peaks gives rise to the chaotic state. By overhauling the
dynamics, we provide insights into the complex dynamics of magnetization in
TI-FM heterostructures.

Figure 1: Magnetization dynamics of TI-FM heterostructure. (a) The schematics
of the TI-FM heterostructure. (b) The equilibrium state of magnetization in
the different relative strengths of SOT $x=\sigma/K$.

##  3 Results

###  3.1 Model and its equilibrium

The physical configuration is illustrated in Fig. 1(a), where a current flows
along the $x$-axis, the itinerant electron spin aligns with the $y$-axis, and
the ferromagnet’s magnetization lie along the $+z$-axis at $t=0$. The polar
and azimuthal angles with respect to the $y$-axis are denoted as $\theta$ and
$\phi$. We represent the magnetization as
$\vec{m}=m(\sin\theta\sin\phi,\cos\theta,\sin\theta\cos\phi)$ and the
itinerant electron spin as $\vec{\sigma}=\sigma(t)(0,1,0)$. For simplicity, we
exclusively focus on the dynamics of the single $\vec{m}$, and ignore that of
$\vec{\sigma}$. This is justified by the stronger spin-momentum locking of the
TI surface state compared to anisotropy 15, 33, 34, 35, 3, 36. For instance,
Ni has the easy-axis anisotropy energy about $2.7~{}\mu$eV$/$atom in the bulk
37 and $4.0~{}$meV$/$atom in the monolayer 38, while the spin-orbit coupling
of Bi2Se3 is estimated as $\sim 1.2~{}$eV$/$atom. Accordingly, we set our unit
of energy to be $1.0~{}\mu$eV $\sim 2.4~{}$GHz and the unit of time to be
$(2.4~{}$GHz$)^{-1}\approx 0.42~{}$ns. The spin flip is then denoted as the
switching of spin direction between $+z$ and $-z$ half space.

Two crucial potential energies emerge: $V_{A}=-\frac{K}{2}m_{z}^{2}$ $(K>0)$,
representing magnetic anisotropy in the ferromagnet, and
$V_{S}=-\gamma\vec{m}\cdot\vec{\sigma}$ signifying the Rashba-field-like 39,
40 (or domain-damping-like) 41 SOT. The equilibrium state of magnetization is
obtained, showcased in Fig. 1(b) for fixed $\vec{\sigma}$. Normalizing
$m=\gamma=1$, key energy scales become the anisotropy energy $K$ and SOT
$\sigma$. The dimensionless parameter $x=\frac{\sigma}{K}$ is defined. As $x$
increases from 0, the equilibrium states initially align along $\pm z$-axis
and gradually rotate towards the $y$-axis. In terms of angles,
$\sin\theta=x,\phi=\pm\pi/2$. Before reaching the $y$-axis, the equilibrium
states are twofold degenerate. For $x\geq 1$, the degenerate equilibrium
states converge at the $y$-axis.

To describe magnetization dynamics out of equilibrium, the Lagrangian density
is considered:

$\displaystyle\mathcal{L}=\mathcal{L}_{B}-V_{A}-V_{S},~{}~{}\mathcal{L}_{B}=m(1-\cos\theta)\dot{\phi}.$
(1)

Note that $\mathcal{L}_{B}$ is the Berry phase term, making polar and
azimuthal angles conjugate 42, 43. Thus, as magnetization rotates from the
$z$-axis to its equilibrium state, deviations from the $yz$-plane induce
complex motion. Introducing Gilbert damping through the Rayleigh dissipation
function 44, $\mathcal{R}=\frac{\eta}{2}\frac{\partial\vec{m}}{\partial
t}\cdot\frac{\partial\vec{m}}{\partial t}$, with the dimensionless damping
constant $\eta$, yields the equations of motion:

$\displaystyle(1+\eta^{2})\dot{\theta}=$
$\displaystyle-\frac{K}{2}\sin\theta\sin 2\phi+\frac{\eta K}{2}\sin
2\theta\cos^{2}\phi-\eta\sigma(t)\sin\theta,\text{ and}$ (2)
$\displaystyle(1+\eta^{2})\dot{\phi}=$
$\displaystyle-K\cos\theta\cos^{2}\phi+\sigma(t)-\frac{\eta K}{2}\sin 2\phi.$
(3)

The initial conditions are set at $\theta=\pi/2$ and $\phi=0$. This system of
equations governs physics for both DC and AC scenarios. The numerical
computation is performed by Mathematica.

It should be noted that the field-like SOT $\sigma(t)$ precesses the
magnetization around $y$-axis while the damping-like SOT $\eta\sigma(t)$ cants
the magnetization toward $y$-axis. This is consistent with the Landau-
Lifshitz-Gilbert equation. It is noteworthy that the relation between the
damping constant $\eta$ and the Gilbert damping constant $\alpha$ is
$\alpha=\eta\gamma$, where $\gamma$ is the gyromagnetic ratio 44. We already
set $\gamma=1$, resulting in $\alpha=\eta$. As a consequence, the field-like
SOT is responsible for the spin flip, whereas the damping-like SOT deters it.

###  3.2 Magnetization dynamics under DC

In the exploration of magnetization dynamics under DC, we initiate with a
straightforward scenario where $\alpha=0$ and $\sigma(t)=0$ for $t<0$ and
$\sigma(t)=\sigma$ for $t\geq 0$. Lacking damping, we anticipate permanent
magnetization precession, described by the simplified equations:

$\displaystyle\dot{\theta}=-K\sin\theta\sin\phi\cos\phi,~{}~{}\dot{\phi}=-K\cos\theta\cos^{2}\phi+\sigma.$
(4)

These equations are intricate to solve analytically, but insight can be gained
by examining two limiting cases, $x\ll 1$ and $x\gg 1$ with $\alpha=0$. In the
first limit, we find oscillatory behavior:

$\displaystyle\theta=\frac{\pi}{2}+x(\cos(Kt)-1),~{}~{}\phi=x\sin(Kt).$ (5)

Both $\theta$ and $\phi$ oscillate with amplitude $x$ and frequency $K$,
termed the oscillating mode. In the second limit:

$\displaystyle\theta=\frac{\pi}{2}+\frac{1}{4x}(\cos(2\sigma
t)-1),~{}~{}\phi=\sigma t.$ (6)

Here, $\theta$ oscillates with amplitude $1/(4x)$ and frequency $2\sigma$,
while $\phi$ monotonically increases. This implies the spin precession about
the $y$-axis and a spin flip occurring at period of $\pi/\sigma$, termed the
flipping mode. Comparing the frequency of $\theta$ dynamics in each mode, we
speculate that the transition from the oscillating to the flipping modes might
occur at $x=1/2$ where $K=2\sigma$.

The guess can be verified by numerical computations of frequency in units of
$K$ for $\sigma(t)=\sigma$ for $t\geq 0$ and $\alpha=0$ in Fig. 2(a). The
frequency variation exhibits $\omega_{DC}\approx K=1$ for $x\ll 1$ and
$\omega_{DC}\approx 2\sigma$ for $x\gg 1$. The singularity at $x=1/2$
indicates the expected transition. Only at this transition point, the spins
are faltering between $z$ axis and $xy$-plane in period, so we call this
faltering mode. For DC, we focus on the oscillating and flipping modes, as the
faltering mode requires a fine tuning.

Figure 2: The magnetization dynamics under DC. (a) Computed frequencies of
oscillating $(x<0.5)$, faltering $(x=0.5)$, and flipping modes $(x>0.5)$ under
DC without damping in varying $x=\sigma/K$. (b) Dynamics of magnetization at
$x=0.6,\alpha=0.05$ (lower) and $x=0.6,\alpha=0.1$ (upper) in time.
Chromatically, the red indicates $+z$, the green indicates in-plane, and the
blue indicates $-z$ directions of magnetization. The direction inward the
paper is the $y$ direction. (c) A 2D phase diagram of the final states in
varying $x$ and $\alpha$. The dark blue indicates the even-flipped final
state, and the yellow indicates the odd-flipped final state.

When damping $\alpha$ is introduced under DC, the magnetization eventually
attains an equilibrium state. For $x\geq 1$, the equilbrium state aligns with
$y$-axis, so the spin flip is absent. However, since two equilibrium states
exists for $x<1$, the selection between them occurs, which can flip the spin.
The selection is determined by both $x$ and $\alpha$. Figure 2(b) displays two
examples of dynamics leading to an even-flipped ($x=0.6,\alpha=0.05$) and odd-
flipped ($x=0.6,\alpha=0.1$) final state, respectively. The magnetization in
the even-flipped state stays in $+z$ half space, while that in the odd-flipped
state stays in $-z$ half space. In both examples, owing to the damping effect,
the crossover from the flipping to the oscillating mode occurs after a
specific time $\tau_{c}$. In Fig. 2(b), for the former, after $\tau_{c}\approx
10.8$ ns, while for the latter, $\tau_{c}\approx 4.2$ ns. Empirically, we find
that $\tau_{c}\propto e^{4.53x}\alpha^{-1}$ for small $\alpha$. [See
Supporting Information (SI).] These show that the variation in $\tau_{c}$
caused by the interplay of $x$ and $\alpha$ serves as a determining factor for
the final state.

We further delve into the relation of the final state with $x$ and $\alpha$ by
a 2D phase diagram in Fig. 2(c). When $x<1/2$, exclusively the even-flipped
final state manifests, whereas for $x\geq 1/2$, both odd-flipped and even-
flipped states emerge, forming a fan-like configuration. It should be noted
that the fan-like configuration is not clear in $\alpha<0.01$ due to the
resolution. The exclusive appearance of even-flipped state for $x<1/2$ is due
to the absence of spin flip in the oscillating mode. However, in the case of
$x\geq 1/2$, the spin flip can occur multiple times before the crossover from
the flipping to the oscillating modes takes place. If the crossover time
exceeds half of the flipping mode period, the spin flips to $-z$ half space.
With a longer crossover time surpassing a full flipping mode period, the spin
flips twice, returning to $+z$ half space. By extending the crossover time
further, the spin flips repeatedly. Even (odd) numbers of spin flips lead to
an even-flipped (odd-flipped) final state, giving rise to the distinctive fan-
like feature in $x\geq 1/2$. This pattern persists until $x=1$, where the
equilibrium state converges to $y$-axis.

###  3.3 Magnetization dynamics under AC

Based on above results, we here explore the magnetization dynamics under AC,
which is given by $\sigma(t)=\sigma\sin(\omega_{\text{AC}}t+\delta)$ for
$t\geq 0$ and $\sigma(t)=0$ for $t<0$. $\omega_{\text{AC}}$ is the driving
frequency, and $\delta$ is the initial phase of AC. The initial state of the
magnetization is again aligned with $+z$-axis. Owing to the damping $\alpha$,
the relative strength of SOT $x$, and driving frequency $\omega_{\text{AC}}$,
the initial state overcomes the irregular dynamics and transits to the final
states after some time. We identify three distinct final states in Fig. 3(a):
adiabatic, resonating, and chaotic states. Both adiabatic and resonating
states are the steady states coming after the decay of initial state. At each
time $t$, $\sigma(t)$ determines an equilbrium state as shown in Fig. 1(b).
The adiabatic state denotes that the magnetization mostly adheres to the
equilibrium state at each time. The resonating state denotes that the
magnetization dynamics is periodic but detach from the equilibrium state at
each time. On the other hand, the chaotic state is evolved from the initial
state, retaining its irregularity, in which the magnetization never has a
periodic motion.

We primarily observe the adiabatic states for low-frequency or strong damping
regime. We observe five distinct modes in the adiabatic states: I) an even-
flipped oscillating mode, II) an odd-flipped oscillating mode, III) an even-
flipped faltering mode, IV) an odd-flipped faltering mode, and V) a
periodically flipping mode. These modes are depicted in Fig. 3(b). The even-
flipped (odd-flipped) oscillating mode or Mode I (II) is the oscillation of
magnetization within $+z$ ($-z$) half space. The even-flipped (odd-flipped)
faltering mode or Mode III (IV) is the repeated faltering of magnetization
between $xy$-plane and $+z$ ($-z$) half space. The periodically flipping mode
or Mode V is the repeated magnetization flip. Notably, Modes I and II (or III
and IV) are almost identical, but differ in the number of spin flips before
decaying to the adiabatic state. Only Mode V shows the continuous spin flip in
time in its adiabatic state. As one can expect from their names, the adiabatic
modes originate from DC modes. Modes I and II come from oscillating mode,
Modes III and IV come from the faltering mode, and Mode V comes from the
flipping mode.

Figure 3: The modes in the adiabatic state under AC. (a) The classification
and comparison of final states. The final states transit with the driving
frequency of AC, from adiabatic to resonating, and to chaotic state in
sequence. While adiabatic and resonating states have periodicity, chaotic
state does not. While adiabatic state adheres to the equilibrium state at each
$t$, others do not. (b) Schematics of distinct AC adiabatic modes after decay.
(c-d) 2D phase diagrams of adiabatic modes in $x$ and $\omega_{\text{AC}}$
with $\alpha=0.03$, at (c) $\delta=0$ and (d) $\delta=\pi/2$. (e-f) 2D phase
diagrams of adiabatic modes in $\delta$ and $\omega_{\text{AC}}$ at (e)
$x=0.6$ and (f) $x=0.8$. The violet denotes Mode I, the red denotes Mode II,
the blue denotes Mode III, the yellow denotes Mode IV, and the green denotes
Mode V.

The modes are chosen by the interplay of $x$, $\omega_{\text{AC}}$, $\delta$,
and $\alpha$. Both $x$ and $\omega_{\text{AC}}$ underscore their importance
after decaying to adiabatic states while $\delta$ and $\alpha$ play a pivotal
role during the decay process. Primary investigation is performed by phase
diagrams for the modes in $x$ and $\omega_{\text{AC}}$ at $\alpha=0.03$ with
different $\delta$ presented in Figs. 3(c,d). For $x<1$, Modes I and II
appear, while for $x\geq 1$, Modes III, IV, and V manifest. This transition
occurs around $x\sim 1$ since $x=1$ is where two equilibrium states meet at
the $y$-axis. Specifically, for $x<1$, as the equilibrium position does not
reach $y$-axis at any time, the magnetization oscillates within either $+z$ or
$-z$ half space. For $x\geq 1$, however, as the equilibrium reaches $y$-axis,
the magnetization either falters on the $xy$-plane or flips its position
between $+z$ and $-z$ half spaces periodically.

$\omega_{\text{AC}}$ also plays a role in determining the modes. Primarily,
$\omega_{\text{AC}}$ chooses the faltering and flipping modes when $x\geq 1$,
since it determines the time duration $\tau_{y}$ that the equilibrium state at
each time stays at the $y$-axis. One can obtain the duration by finding the
maximum $\tau_{y}$ satisfying $\sigma(t)\geq 1$ in
$t\in[t_{0},t_{0}+\tau_{y}]$. In the case of sinusoidal $\sigma(t)$, this can
be expressed as $\tau_{y}=\frac{1}{\omega_{\text{AC}}}(\pi-2\zeta)$, where
$\zeta=\arcsin(1/x)\in(0,\pi/2]$. For $\tau_{y}$, as the equilibrium state is
at the $y$-axis, the magnetization modulates near the $y$-axis. After
$\tau_{y}$, the equilibrium position is divided again and deviates away from
the $y$-axis. Then, depending on its modulated position, the magnetization
chooses one of the equilibrium state, which leads to either faltering or
flipping modes. Additionally, the window of $x$ for every mode at higher
$\omega_{\text{AC}}$ is opened up wider than that at lower
$\omega_{\text{AC}}$, showing the fan-like feature in Figs. 3(c,d). Unlike the
faltering mode under DC, the window of Modes III and IV, originating from the
faltering mode, expands to the finite range of $x$. Lastly, it is noteworthy
that the phase of $\omega_{AC}=0$ in Figure 3(d) is the same as the phase in
$\alpha=0.03$ line of Fig. 2(c). This substantiates that adiabatic modes can
be derived from DC modes.

On the other hand, $\delta$ acts as a switch to turn on the spin flip during
the decay process to adiabatic states. Specifically, comparing Figs. 3(b) to
(c), Modes II and IV barely appear when $\delta=0$, which means that the spin
flip is turned off near $\delta=0$. This happens because the average of
$\sigma(t)$ during the decay time is small, so the spin cannot flip before
decaying the adiabatic state. This can be supported by Figs. 3(e-f). In Fig.
3(e), we present phase diagrams for the modes in $\delta$ and
$\omega_{\text{AC}}$ at $x=0.6$. The transition from Modes I to II occurs once
during the increase of $\delta$. This happens because the spin flip is turned
on by finite $\delta$. Increasing $x$, the spin flips more times just as in DC
case, so the repeated transition between Modes I and II can also be observed
in Fig. 3(f). We should note that the increase of $\alpha$ reduces the decay
time. Above arguments holds also for the triangular wave instead of sinusoidal
wave. [See SI.]

Figure 4: The transition from adiabatic to resonating and chaotic states under
AC. (a) 2D phase diagrams of the final states in $\omega_{\text{AC}}$ and
$\alpha$. The adiabatic state is in dark blue, the resonating state is in
green, and the chaotic state is in dark red. The top panel is at $x=0.3$ and
the bottom panel is at $x=0.7$. (b) The schematics of change in Fourier
transform in for each state. (c) The corresponding examples to each state of
(b) at $x=0.3$ and $\alpha=0$.

We move on to high-$\omega_{\text{AC}}$ and low-$\alpha$ regime to investigate
ferromagnetic magnetoresonance (FMR). One could expect that the resonance
frequency is closely related to the frequencies $\omega_{\text{DC}}$ in Fig.
2(a). In fact, when $\omega_{\text{AC}}$ approaches $\omega_{\text{DC}}$, the
oscillating amplitude becomes larger, and eventually the adiabatic state
changes to the resonating state and to the chaotic state sequentially. By
checking the final state after $t\sim 1~{}$ms, we present 2D phase diagrams in
$\omega_{\text{AC}}$ and $\alpha$ at $x=0.3$ and $0.7$ in Fig. 4(a). Both
diagrams show similar behavior. It begins with the adiabatic state at low
frequencies, confronts the transition to the resonating state, and reaches the
chaotic state at a certain range of $\omega_{\text{AC}}$. For $x=0.7$, the
window of chaotic states is opened up widely as SOT becomes large.

We further illuminate the transition of states by Fourier transform. The
transition is illustrated in Fig. 4(b), where the change in $|\phi(\omega)|$
by $\omega_{\text{AC}}$ is schematically presented. Here, $\phi(\omega)$
corresponds to the Fourier transform of azimuthal angle $\phi(t)$ from $t=0$
to $t=0.21~{}$ms ($5\times 10^{5}$ unit times). Figure 4(c) showcases
representative examples of Fourier transforms at $x=0.3$ and $\alpha=0$. For
adiabatic and resonating states, four distinct peaks are observed in the top
panels: i) a main peak at $\omega_{1}$ derived from $\omega_{\text{DC}}$, ii)
another main peak driven by AC at $\omega_{2}=\omega_{\text{AC}}$, iii)
induced peaks from i) and ii) at
$\omega_{3}=\omega_{\text{1}}+n(\omega_{\text{1}}-\omega_{\text{AC}})$, and
iv) subpeaks at $\omega_{4}=|\omega_{1,2,3}\pm 2n\omega_{\text{AC}}|$
($n\in\mathbb{N}$). Notably, as the time interval is elongated from 0.21 ms,
Peaks i and iii decay while Peaks ii and iv gain intensity. [See SI.] This
means that Peaks i and iii are related to the decay process while Peaks ii and
iv are related to the stable state after decaying. Increasing
$\omega_{\text{AC}}$ to the higher frequency, the change only occurs to the
distance between peaks, as the frequency of Peak i decreases and that of Peak
ii increases. Thus, adiabatic and resonating states are indistinguishable
solely by Fourier transform. This is substantiated by the phase diagrams in
Fig. 4(a), where the transition from adiabatic and resonating states is not
distinctly delineated. This behavior is consistent for all modes in Fig. 3(b).
Near the chaotic state, the peaks form array in a period of
$\omega_{\text{AC}}-\omega_{1}$ as shown in the middle panels. As
$\omega_{\text{AC}}$ increases more, chaos sets in, destroying all peaks and
being $\phi(\omega)\propto 1/\omega$ as shown in the bottom panels. This
behavior does not change although the range of time is elongated.

##  4 Discussion

We address here about the typical values of parameters in the realistic
systems. The typical value of ferromagnetic anisotropy is $1-10~{}\mu$eV/atom
26, that of the current density in the experiments is $\sim
10^{7}~{}$A$\cdot$cm-2, and that of the itinerant spin polarization by Rashba-
Edelstein Effect is estimated about $10^{-4}\hbar$ per unit cell 45. The
typical value of Gilbert damping constant is $\sim 10^{-3}-10^{-2}$ 46, and
that of the unit of AC frequency is estimated about $1-10$ GHz.

So far, we discuss the magnetization dynamics at the interface of TI-FM
heterostructure. Under DC, we observe the oscillating and flipping modes
without damping, which is determined by the relative strength of anisotropy
and SOT. With damping, the number of spin flips is determined by the duration
of the crossover from flipping to the oscillating mode. Under AC, we observe
five distinct modes in the low frequency regime and their evolution to the
resonating and chaotic states in the high frequency regime. Although we mainly
discuss the TI-FM heterostructure due to its high efficiency, our work can be
applied to the general systems with strong Rashba spin-orbit coupling. This is
because the assumption underlying in our work is only that the current carries
finite spin due to the Rashba-Edelstein Effect. Our work offers insights into
the spin control by spin-orbit coupling, underscoring the practical aspects of
the world of spintronics.

We thank Wataru Koshibae for the fruitful discussions. This work was supported
by JST, CREST Grant Number JPMJCR1874, Japan.

This material is available free of charge via the Internet at
https://pubs.acs.org/

$\bullet$ The dependence of crossover time from flipping to oscillating mode
on $x$ and $\alpha$, the dependence of decay time to adiabatic modes on
$\alpha$, the comparison between triangular and sinusoidal waves, and the time
evolution of peaks in Fourier transform.

See pages - of Supp.pdf

## References

  * Witczak-Krempa et al. 2014 Witczak-Krempa, W.; Chen, G.; Kim, Y. B.; Balents, L. Correlated quantum phenomena in the strong spin-orbit regime. _Annu. Rev. Condens. Matter Phys._ 2014, _5_ , 57–82
  * Hasan and Kane 2010 Hasan, M. Z.; Kane, C. L. Colloquium: topological insulators. _Reviews of modern physics_ 2010, _82_ , 3045
  * Qi and Zhang 2011 Qi, X.-L.; Zhang, S.-C. Topological insulators and superconductors. _Reviews of Modern Physics_ 2011, _83_ , 1057
  * Pi et al. 2010 Pi, U. H.; Won Kim, K.; Bae, J. Y.; Lee, S. C.; Cho, Y. J.; Kim, K. S.; Seo, S. Tilting of the spin orientation induced by Rashba effect in ferromagnetic metal layer. _Applied Physics Letters_ 2010, _97_
  * Liu et al. 2011 Liu, L.; Moriyama, T.; Ralph, D.; Buhrman, R. Spin-torque ferromagnetic resonance induced by the spin Hall effect. _Physical review letters_ 2011, _106_ , 036601
  * Kim et al. 2013 Kim, J.; Sinha, J.; Hayashi, M.; Yamanouchi, M.; Fukami, S.; Suzuki, T.; Mitani, S.; Ohno, H. Layer thickness dependence of the current-induced effective field vector in Ta— CoFeB— MgO. _Nature materials_ 2013, _12_ , 240–245
  * Che et al. 2020 Che, X.; Pan, Q.; Vareskic, B.; Zou, J.; Pan, L.; Zhang, P.; Yin, G.; Wu, H.; Shao, Q.; Deng, P., et al. Strongly surface state carrier-dependent spin–orbit torque in magnetic topological insulators. _Advanced Materials_ 2020, _32_ , 1907661
  * Brataas et al. 2012 Brataas, A.; Kent, A. D.; Ohno, H. Current-induced torques in magnetic materials. _Nature materials_ 2012, _11_ , 372–381
  * Shao et al. 2021 Shao, Q.; Li, P.; Liu, L.; Yang, H.; Fukami, S.; Razavi, A.; Wu, H.; Wang, K.; Freimuth, F.; Mokrousov, Y., et al. Roadmap of spin–orbit torques. _IEEE transactions on magnetics_ 2021, _57_ , 1–39
  * Nagaosa 2019 Nagaosa, N. Emergent inductor by spiral magnets. _Japanese Journal of Applied Physics_ 2019, _58_ , 120909
  * Yokouchi et al. 2020 Yokouchi, T.; Kagawa, F.; Hirschberger, M.; Otani, Y.; Nagaosa, N.; Tokura, Y. Emergent electromagnetic induction in a helical-spin magnet. _Nature_ 2020, _586_ , 232–236
  * Kitaori et al. 2021 Kitaori, A.; Kanazawa, N.; Yokouchi, T.; Kagawa, F.; Nagaosa, N.; Tokura, Y. Emergent electromagnetic induction beyond room temperature. _Proceedings of the National Academy of Sciences_ 2021, _118_ , e2105422118
  * Ieda and Yamane 2021 Ieda, J.; Yamane, Y. Intrinsic and extrinsic tunability of Rashba spin-orbit coupled emergent inductors. _Physical Review B_ 2021, _103_ , L100402
  * Yamane et al. 2022 Yamane, Y.; Fukami, S.; Ieda, J. Theory of Emergent Inductance with Spin-Orbit Coupling Effects. _Physical review letters_ 2022, _128_ , 147201
  * Araki and Ieda 2023 Araki, Y.; Ieda, J. Emergence of inductance and capacitance from topological electromagnetism. _Journal of the Physical Society of Japan_ 2023, _92_ , 074705
  * Yasuda et al. 2016 Yasuda, K.; Tsukazaki, A.; Yoshimi, R.; Takahashi, K.; Kawasaki, M.; Tokura, Y. Large unidirectional magnetoresistance in a magnetic topological insulator. _Physical review letters_ 2016, _117_ , 127202
  * Yasuda et al. 2017 Yasuda, K.; Tsukazaki, A.; Yoshimi, R.; Kondou, K.; Takahashi, K.; Otani, Y.; Kawasaki, M.; Tokura, Y. Current-nonlinear Hall effect and spin-orbit torque magnetization switching in a magnetic topological insulator. _Physical review letters_ 2017, _119_ , 137204
  * Alegria et al. 2014 Alegria, L.; Ji, H.; Yao, N.; Clarke, J.; Cava, R. J.; Petta, J. R. Large anomalous Hall effect in ferromagnetic insulator-topological insulator heterostructures. _Applied Physics Letters_ 2014, _105_
  * Katmis et al. 2016 Katmis, F.; Lauter, V.; Nogueira, F. S.; Assaf, B. A.; Jamer, M. E.; Wei, P.; Satpati, B.; Freeland, J. W.; Eremin, I.; Heiman, D., et al. A high-temperature ferromagnetic topological insulating phase by proximity coupling. _Nature_ 2016, _533_ , 513–516
  * Wang et al. 2020 Wang, H.; Liu, Y.; Wu, P.; Hou, W.; Jiang, Y.; Li, X.; Pandey, C.; Chen, D.; Yang, Q.; Wang, H., et al. Above room-temperature ferromagnetism in wafer-scale two-dimensional van der Waals Fe3GeTe2 tailored by a topological insulator. _ACS nano_ 2020, _14_ , 10045–10053
  * Ou et al. 2023 Ou, Y.; Mirzhalilov, M.; Nemes, N. M.; Martinez, J. L.; Rocci, M.; Akey, A.; Ge, W.; Suri, D.; Wang, Y.; Ambaye, H., et al. Enhanced Ferromagnetism in Monolayer Cr2Te3 via Topological Insulator Coupling. _arXiv preprint arXiv:2312.15028_ 2023,
  * Fan et al. 2014 Fan, Y.; Upadhyaya, P.; Kou, X.; Lang, M.; Takei, S.; Wang, Z.; Tang, J.; He, L.; Chang, L.-T.; Montazeri, M., et al. Magnetization switching through giant spin–orbit torque in a magnetically doped topological insulator heterostructure. _Nature materials_ 2014, _13_ , 699–704
  * Mellnik et al. 2014 Mellnik, A.; Lee, J.; Richardella, A.; Grab, J.; Mintun, P.; Fischer, M. H.; Vaezi, A.; Manchon, A.; Kim, E.-A.; Samarth, N., et al. Spin-transfer torque generated by a topological insulator. _Nature_ 2014, _511_ , 449–451
  * Han et al. 2017 Han, J.; Richardella, A.; Siddiqui, S. A.; Finley, J.; Samarth, N.; Liu, L. Room-temperature spin-orbit torque switching induced by a topological insulator. _Physical review letters_ 2017, _119_ , 077702
  * Wang et al. 2017 Wang, Y.; Zhu, D.; Wu, Y.; Yang, Y.; Yu, J.; Ramaswamy, R.; Mishra, R.; Shi, S.; Elyasi, M.; Teo, K.-L., et al. Room temperature magnetization switching in topological insulator-ferromagnet heterostructures by spin-orbit torques. _Nature communications_ 2017, _8_ , 1364
  * Mogi et al. 2021 Mogi, M.; Yasuda, K.; Fujimura, R.; Yoshimi, R.; Ogawa, N.; Tsukazaki, A.; Kawamura, M.; Takahashi, K. S.; Kawasaki, M.; Tokura, Y. Current-induced switching of proximity-induced ferromagnetic surface states in a topological insulator. _Nature communications_ 2021, _12_ , 1404
  * Yokoyama et al. 2010 Yokoyama, T.; Zang, J.; Nagaosa, N. Theoretical study of the dynamics of magnetization on the topological surface. _Physical Review B_ 2010, _81_ , 241410
  * Garate and Franz 2010 Garate, I.; Franz, M. Inverse spin-galvanic effect in the interface between a topological insulator and a ferromagnet. _Physical review letters_ 2010, _104_ , 146802
  * Pesin and MacDonald 2012 Pesin, D.; MacDonald, A. H. Spintronics and pseudospintronics in graphene and topological insulators. _Nature materials_ 2012, _11_ , 409–416
  * Araki and Ieda 2021 Araki, Y.; Ieda, J. Intrinsic torques emerging from anomalous velocity in magnetic textures. _Physical Review Letters_ 2021, _127_ , 277205
  * Wu et al. 2021 Wu, H.; Chen, A.; Zhang, P.; He, H.; Nance, J.; Guo, C.; Sasaki, J.; Shirokura, T.; Hai, P. N.; Fang, B., et al. Magnetic memory driven by topological insulators. _Nature communications_ 2021, _12_ , 6251
  * Yamanouchi et al. 2022 Yamanouchi, M.; Araki, Y.; Sakai, T.; Uemura, T.; Ohta, H.; Ieda, J. Observation of topological Hall torque exerted on a domain wall in the ferromagnetic oxide SrRuO3. _Science advances_ 2022, _8_ , eabl6192
  * Daalderop et al. 1990 Daalderop, G.; Kelly, P.; Schuurmans, M. First-principles calculation of the magnetocrystalline anisotropy energy of iron, cobalt, and nickel. _Physical Review B_ 1990, _41_ , 11919
  * König et al. 2008 König, M.; Buhmann, H.; W. Molenkamp, L.; Hughes, T.; Liu, C.-X.; Qi, X.-L.; Zhang, S.-C. The quantum spin Hall effect: theory and experiment. _Journal of the Physical Society of Japan_ 2008, _77_ , 031007
  * Cullity and Graham 2011 Cullity, B. D.; Graham, C. D. _Introduction to magnetic materials_ ; John Wiley & Sons, 2011
  * Xu et al. 2020 Xu, Y.-K.; Li, H.; He, B.-G.; Cheng, Z.-P.; Zhang, W.-B. Electronic structure and magnetic anisotropy of single-layer rare-earth oxybromide. _ACS omega_ 2020, _5_ , 14194–14201
  * Fritsche et al. 1987 Fritsche, L.; Noffke, J.; Eckardt, H. A relativistic treatment of interacting spin-aligned electron systems: application to ferromagnetic iron, nickel and palladium metal. _Journal of Physics F: Metal Physics_ 1987, _17_ , 943
  * Gay and Richter 1986 Gay, J.; Richter, R. Spin anisotropy of ferromagnetic films. _Physical review letters_ 1986, _56_ , 2728
  * Liu et al. 2012 Liu, L.; Lee, O.; Gudmundsen, T.; Ralph, D.; Buhrman, R. Current-induced switching of perpendicularly magnetized magnetic layers using spin torque from the spin Hall effect. _Physical review letters_ 2012, _109_ , 096602
  * Manchon et al. 2019 Manchon, A.; Železnỳ, J.; Miron, I. M.; Jungwirth, T.; Sinova, J.; Thiaville, A.; Garello, K.; Gambardella, P. Current-induced spin-orbit torques in ferromagnetic and antiferromagnetic systems. _Reviews of Modern Physics_ 2019, _91_ , 035004
  * Kim 2019 Kim, B. S. Skyrmions and Hall transport. _Journal of Physics: Condensed Matter_ 2019, _31_ , 383001
  * Auerbach 2012 Auerbach, A. _Interacting electrons and quantum magnetism_ ; Springer Science & Business Media, 2012
  * Tatara and Fukuyama 2014 Tatara, G.; Fukuyama, H. Phasons and excitations in skyrmion lattice. _Journal of the Physical Society of Japan_ 2014, _83_ , 104711
  * Gilbert 2004 Gilbert, T. L. A phenomenological theory of damping in ferromagnetic materials. _IEEE transactions on magnetics_ 2004, _40_ , 3443–3449
  * Chen 2019 Chen, W. Edelstein and inverse Edelstein effects caused by the pristine surface states of topological insulators. _Journal of Physics: Condensed Matter_ 2019, _32_ , 035809
  * Barati et al. 2013 Barati, E.; Cinal, M.; Edwards, D.; Umerski, A. Calculation of Gilbert damping in ferromagnetic films. EPJ Web of Conferences. 2013; p 18003