# Topological equivalence of stripy states and skyrmion crystals

X. R. Wang  [<EMAIL_ADDRESS>Xu-Chong Hu  [  Zhou-Zhou Sun  [

###### Abstract

Stripy states, consisting of a collection of stripy spin textures, are the
precursors of skyrmion crystals (SkXs). Common belief is that stripy states
and SkXs are topologically unconnected, and transitions between SkXs and
stripy states are phase transitions. Here, we show that both stripy states and
SkXs are skyrmion condensates and they are topologically equivalent. By
gradually tuning the stripe whose width goes from smaller than to larger than
skyrmion-skyrmion separation, the structure of a skyrmion condensate
transforms smoothly and continuously from various stripy phases, including
helical states and mazes, to crystals, showing that stripy states are
topologically connected to SkXs.

###### keywords:

Nano-magnetism, chiral magnets, topological physics, skyrmion crystals,
helical states

The Hong Kong University of Science and Technology] Physics Department, The
Hong Kong University of Science and Technology (HKUST), Clear Water Bay,
Kowloon, Hong Kong [HKUST Shenzhen Research Institute] HKUST Shenzhen Research
Institute, Shenzhen 518057, China The Hong Kong University of Science and
Technology] Physics Department, The Hong Kong University of Science and
Technology (HKUST), Clear Water Bay, Kowloon, Hong Kong [HKUST Shenzhen
Research Institute] HKUST Shenzhen Research Institute, Shenzhen 518057, China
South China Business College] South China Business College, Guangdong
University of Foreign Studies, Guangzhou 510545, China [HKUST Shenzhen
Research Institute] HKUST Shenzhen Research Institute, Shenzhen 518057, China
SkX

Magnetic skyrmions have attracted much attention in recent years for their
academic interest and potential applications 1, 2, 3, 4, 5, 8, 6, 7. Various
aspects of magnetic skyrmions have been extensively and intensively studied
11, 28, 21, 12, 9, 10, 13, 16, 17, 18, 19, 23, 24, 25, 27, 22, 26, 14, 15, 20,
30, 32, 29, 31, 33, 34, 36, 37, 35, 39, 38, 41, 43, 42, 44, 40, 45, 46. Vast
experiments 47 showed that various stripy phases, including the helical states
and mazes, accompany skyrmion crystals (SkXs). General belief is that stripy
states and SkXs are topologically unconnected 47, 48, and transformations
between them are phase transitions. Testing this belief motivates current
study.

For a perpendicularly magnetized chiral film described by magnetization
$\bf{M}$, the Heisenberg exchange stiffness constant $A$, anisotropic
coefficient $K$, and the Dzyaloshinskii-Moriya interaction (DMI) coefficient
$D$, $\kappa\equiv(\pi^{2}D^{2})/(16AK)=1$, an important quantity, separates
isolated skyrmions ($\kappa<1$) from condensed stripe skyrmions ($\kappa>1$)
in the absence of a magnetic field 49, 50, 51, 52. In this paper, we study how
the morphology of a collection of stripe skyrmions changes as the stripe
widths gradually increase. Our results show the topological equivalence of
stripy states and SkXs and raise the question about the notion of the first-
order phase transition in SkX formations.

We first investigate how morphology of 50 skyrmions, in films of
$400\,$nm$\,\times 400\,$nm$\,\times 0.5\,$nm, changes with stripe width.
Previous study 49 showed that one nucleation domain in a chiral magnetic film
of $\kappa=(\pi D)^{2}/(16AK)>1$ develops into a stripe skyrmion. Thus, 50
nucleation domains of $m_{z}=1$ embedded in the $m_{z}=-1$ background of
chiral magnetic films can be used to generate 50 skyrmions. The 50 domains of
6 nm in diameter each are initially distributed in films randomly as shown in
Fig. 1(a1). Figure 1(a2) shows a stable structure of 50 skyrmions for Film 1
with $A/D=1\,$nm (a2). All stripes have the same width satisfying
$L=a(\kappa)A/D$ with $a(\kappa=2)$=6.5 49. $2L=13\,$nm for Film 1 is much
smaller than average skyrmion-skyrmion separation $d_{ss}=\sqrt{400\times
400/50}\simeq 57\,$nm. Since each nucleation domain has a large enough space
to grow independently, the final stable structure is a dense maze formed by 50
ramified/non-ramified stripe skyrmions. Figures 1(a3) is a stable structure
for Film 3 of $A/D=2\,$nm and the dense maze becomes an irregularly arranged
42 stripe skyrmions and 8 circular-like skyrmions. Figures 1(a4)-(a6) are
stable structures for Film 5 with $A/D=3\,$nm (a4); Film 7 with $A/D=4\,$nm
(a5); and Film 9 with $A/D=5\,$nm (a6). Figures 1(a4)-(a5) are the mixtures of
20 circular skyrmions and 30 stripe skyrmions when $L=13.02\,$nm (a4), 40
circular skyrmions and 10 stripe skyrmions when $L=26.04\,$nm (a5),
respectively. SkXs of 50 circular skyrmions form when $2L$ is larger than 57
nm as shown in Fig. 1(a6). These results demonstrate the importance of
$L/d_{ss}$. $L/d_{ss}=0.5$ roughly separates SkXs from mixture of stripe
skymrions and circular skyrmions, as shown in the $L/d_{ss}$ phase diagram.

Figure 1: (a1) Initial configuration of randomly distributed 50 circular
nucleation domains of $m_{z}=1$. The diameter of each domain is $6\,$nm.
(a2)-(a6) The stable spin structures with different $A/D$, starting from
initial configuration (a1). $A/D=1\,$nm (a2), $2\,$nm (a3), $3\,$nm (a4),
$4\,$nm (a5), $5\,$nm (a6), corresponding to Films 1, 3, 5, 7, and 9 in Tab.
1. $\kappa=2$ for all films. $L/d_{ss}$ shows the phase-diagram. (b1)-(b3) The
time dependences of skyrmion number $Q$ (the right y-axis and the blue curves)
and total magnetic energy $E$ (the left y-axis and the red curves).
$A/D=1\,$nm, $2\,$nm (b1); $3\,$nm, $4\,$nm (b2); $5\,$nm (b3). The time is in
the logarithmic scale. All $Q$ reaches the expected 50 within 1 ps, and $E$
approaches its minimum in nanoseconds.

Interestingly, both positive and negative skyrmion charges, whose density is
defined as
$\rho=\frac{1}{4\pi}\mathbf{m}\cdot(\partial_{x}\mathbf{m}\times\partial_{y}\mathbf{m})$,
appear in a stripe skyrmion while circular skyrmions carry only positive
charges. The colours in Figs. 1(a2)-(a6) encode skyrmion charge (positive in
blue and negative in red) distribution. Indeed, each small domain of initial
zero skyrmion number becomes a skyrmion of $Q=\int\rho\rm dx\rm dy=1$ within 1
picoseconds. Each pattern in Figs. 1(a2)-(a6) consists of 50 skyrmions,
evident from the time-dependence of $Q(t)$ in Figs. 1(b1)-(b3). $Q(t)$ reaches
$50$ within picoseconds. The total energy $E$ monotonically decreases with
time and approaches its minimal values in nanoseconds for all films as shown
in Figs. 1(b1)-(b3) (the left y-axis and the red curves).

In the second set of simulations, starting from the stable maze structure of
Fig. 1(a2) with 50 ramified stripe skyrmions for Film 1 of $A/D=1\,$nm as the
initial configuration, we consecutively increase $A$ by $0.5\,\rm pJ/m$ and
decrease $K_{u}$ by proper values to keep $\kappa=2$ every $20\,$ns to
simulate Film 2 to Film 10 listed in Tab. 1. Figures 2(a1)-(a6) are the stable
structures for Film 2, 3, 4, 5, 7, and 9, respectively. Each of them comes
from the stable structure of its predecessor. The time-dependences of $E$ and
$Q$ are plotted in Fig. 2(b), showing that systems reach their equilibrium
states (constant $E$ and $Q$) within a few nanoseconds with new set of
parameters. $Q=50$ holds all the time due to the topological nature of the
quantity.

Figure 2: (a1) Stable structure of Film 2 developed from the stable structure
of Film 1 of Fig. 1(a2). (a2) Stable structure of Film 3 of $A/D=2\,$nm
developed from (a1). (a3) Stable structure of Film 4 of $A/D=2.5\,$nm from
(a2). (a4) Stable structure of Film 5 of $A/D=3\,$nm from (a3). (a5) Stable
structure of Film 7 of $A/D=4\,$nm. (a6) Stable structure of Film 9 of
$A/D=5\,$nm. $\kappa=2$ is for all films. Insets: Stable structures at
temperature of 10$\,K$ for Film 3 (a2); Film 4 (a3); Film 5 (a4); Film 7 (a5);
Film 9 (a6). $L/d_{ss}$ shows the phase-diagram. (b) The blue line ($0\,$K)
and the green dash line ($10\,$K) are the time evolution of skyrmion number
$Q$ (the right y-axis). The time is in the logarithmic scale (up x-axis). The
red curve is the total energy $E$ (the left y-axis) as a function of time in
the normal scale (down x-axis). $E$ for structures in (a1)-(a6) are marked by
cross. $Q$ maintains at 50 for both zero temperature and $10\,K$. $E$
approaches a constant in nanoseconds after tuning the parameters. The inset is
total magnetic energy of the stable structure as a function of $D^{2}/A$. The
symbols are the simulation results and the blue curve is the fit to
$E=c(D^{2}/A)$ with $c=(-7.8\times 10^{3}\,nm^{3})$.

For a given spin structure, the total magnetic energy is, according to Eq. (4)
(see Methods),

$\displaystyle E$
$\displaystyle=d\int\int\\{A|\nabla\mathbf{m}|^{2}+D[m_{z}\nabla\cdot\mathbf{m}$
(1)
$\displaystyle\quad-(\mathbf{m}\cdot\nabla)m_{z}]+K(1-m_{z}^{2})\\}\mathrm{d}x\mathrm{d}y$
$\displaystyle=\frac{\pi^{2}D^{2}d}{16A}\int\int\\{|\frac{4A}{\pi
D}\nabla\mathbf{m}|^{2}+\frac{4}{\pi}[m_{z}\frac{4A}{\pi
D}\nabla\cdot\mathbf{m}$ $\displaystyle\quad-(\mathbf{m}\cdot\frac{4A}{\pi
D}\nabla)m_{z}]+\frac{1}{\kappa}(1-m_{z}^{2})\\}\mathrm{d}x\mathrm{d}y$
$\displaystyle=\kappa
KdL^{2}\int\int\\{|\nabla\mathbf{m}|^{2}+\frac{4}{\pi}[m_{z}\nabla\cdot\mathbf{m}$
$\displaystyle\quad-(\mathbf{m}\cdot\nabla)m_{z}]+\frac{1}{\kappa}(1-m_{z}^{2})\\}\mathrm{d}x\mathrm{d}y.$

Metastable spin structures depend only explicitly on $\kappa=(\pi
D)^{2}/(16AK)$, and inexplicitly on the length scale of $L_{1}=4A/(\pi D)$.
Since the integral area is shrunk by $L_{1}^{2}$, the total magnetic energy is
proportional to $\kappa K\propto D^{2}/A$ in the case of constant $\kappa$.
Symbols in the inset of Fig. 2(b) are the total magnetic energies of
metastable states for $A/D=1\,\rm nm$ up to $5\,\rm nm$. The solid curve is
the fit to $E(A/D)=cD^{2}/A$, showing excellent agreement with simulations.
The energy deviates slightly from $cD^{2}/A$ for larger $A/D$ such that $L$
($L_{1}$) is comparable to skyrmion-skyrmion separation and skyrmion
interactions becomes important.

The spikes in $E$-curve come from sudden energy increase when parameters are
modified. The energy of the new stable structure is higher (less negative) as
$A$ increases. Compare Fig. 2(a2) with Fig. 1(a3) for Film 3, Fig. 2(a4) with
Fig. 1(a4) for Film 5, and Fig. 2(a5) with Fig. 1(a5) for Film 7, each pair
has similar, but not identical, stable structures, revealing the fact of a
large number of metastable structures for a collection of stripe skyrmions.
Which metastable structure that a mixture of circular and stripe skyrmions
ends up depends on the history of the system, see more discussions on this
issue in Supporting Information.

Mumax3 simulations 53 at $T=10\,K$, much lower than the Curie temperatures
$T_{c}$ ($84\,$K$\sim 290\,$K for Film 2 to Film 9), are used to demonstrate
the thermal stability, see Supporting Information for the detail of
determining $T_{c}$. The insets of Fig. 2(a1)-(a6) are the metastable
structures of Film 2 at 20ns (a1), Film 3 at 40ns (a2), Film 4 at 60ns (a3),
Film 5 at 80 ns (a4), Film 7 at 120 ns (a4), and Film 9 at 160 ns (a6). $Q=50$
hold throughout 180 ns and morphology of the final states are similar as those
at zero temperature for all films as shown in Fig. 2(b). However, one distinct
difference is that skyrmion surfaces are more irregular and rough at finite
temperatures, agreeing with experiments 54, which may be understood from the
negative skyrmion formation energy. Negative formation energy means negative
surface tension such that such a surface can hardly resist external
perturbations and deformations.

Figure 3:  (a1) Initial configuration of 50 circular nucleation domains of
$m_{z}=1$ arranged orderly in two columns. The diameter of each domain is
$6\,$nm. (a2)-(a6) Stable structures for $A/D=1\,$nm (a2); $2\,$nm (a3);
$3\,$nm (a4); $4\,$nm (a5); and $5\,$nm (a6). All films have $\kappa=2$. (b)
The time evolution of skyrmion number $Q$ (the right y-axis and the blue
curve). Time is in the logarithmic scale (the top x-axis). The red curve is
the total energy $E$ (the left y-axis) as a function of time in the linear
scale (the bottom x-axis). $Q$ reaches 50 within 1 ps and keeps this value
afterword even when the parameters are tuned to different values. $E$ for
structures in (a2)-(a6) are marked by cross. The inset is the $A/D$ dependence
of $G_{6}-G_{2}$. (c1)-(c4) The red and the blue points are the z-component of
magnetization along the green lines labeled by circled $n$ in (a2)-(a5). The
green lines are the fit to our approximate profile. The cyan arrows indicate
the stripe widths of $m_{z}<0$ or $m_{z}>0$.

In contrast to a large number of metastable structures in stripe skyrmion
condensate, variations of structures in SkXs are limited to a global
translational transformation of all skyrmions. This explains why SkXs in Figs.
2(a6) and 1(a6) are almost the same although their kinetic and dynamical
process are not the same. To further verify the independence of SkXs on their
initial configurations, we repeat the simulations for Film 1, starting with
the same 50 nucleation domains as those in Fig. 1 but arranged in two well-
ordered columns as shown in Fig. 3(a1). Figure 3(a2) is the stable helical
state after $20\,$ns evolution with 50 stripe skyrmions ($Q(t)=50$ as shown in
Fig. 3(b)). Stripes tilt about $56^{\circ}$ from the $y$-direction. The stripe
tilt comes from film tiling by 25 stripes of width $L\simeq 6.5\,$nm. The
tilted angle $\theta$ satisfies $2L/(\sin\theta)=400/25$. For $L=6.5\,$nm,
$\theta=\cos^{-1}(13/16)\simeq 54^{\circ}$, very close to the observed value
of $56^{\circ}$.

After consecutively changing model parameters listed in Tab. 1 for Film 3 to
10 every 20 ns, Figs. 3(a3)-(a6) show stable structures for Film 3 with
$A/D=2\,$nm at $60\,$ns (a3), Film 5 with $A/D=3\,$nm at $100\,$ns (a4), Film
7 with $A/D=4\,$nm at $140\,$ns (a5), and Film 9 with $A/D=5\,$nm at $180\,$ns
(a6). Figure 3(b) shows time-dependences of $E$ and $Q$. The spikes in
$E$-curve are from the sudden change of model parameters. $Q$ maintains the
value of 50 after the first a few picoseconds while $E$ reaches its minimal
value in nanoseconds every time after parameter changes. The number of stripe
(circular) skyrmions decreases (increases) as $L$ increases. All skyrmions in
Fig. 3(a6) are circular when $L$ is larger than the average skyrmion-skyrmion
separation.

To demonstrate that structure transformation from stripy states to SkXs is
smooth, we introduce a local orientation order parameter, widely used in
studying triangular lattices 55, 56, 57, for $i$’th skyrmion

$\Psi_{n}(i)=\sum_{NN}e^{in\theta_{ij}}.$ (2)

$\theta_{ij}$ is the angle of the direction from the center of skyrmion $i$ to
the center of skyrmion $j$ with respect to the x-axis, and sum is over all
nearest neighbouring skyrmions $j$. For a lattice with orientation order of a
triangular lattice, $\Psi_{6}(i)$ is close to 1 while it is a small complex
number for a lattice without orientation order. For an imperfect lattice
structure, averaged $\Psi_{n}$ (over all skyrmions), called $G_{n}$,

$G_{n}=\frac{1}{N}\sum_{i}\Psi_{n}(r_{i})\Psi_{n}^{*}(r_{i}),$ (3)

can distinguish one skyrmion structure from another 55, 56, 57. $G_{6}$ for a
SkX in a perfect triangular lattice is 1. Since $G_{6}$ is also close to 1 for
ordered stripe states such as those in Figs. 3(a2)-(a3), we use $G_{6}-G_{2}$,
which is close to 0 for a helical state and 1 for a SkX, as an order parameter
to distinguish the two structures. $G_{6}-G_{2}$ undergoes a sudden transition
for a first-order phase transition and a smooth variation for a non-phase
transition. The inset of Fig. 3(b) shows how $G_{6}-G_{2}$ varies with $A/D$
from 0 for an ordered stripe skyrmion state, via maze and mixture of stripes
and circular skyrmions, to 1 for SkXs. The smooth variation of $G_{6}-G_{2}$
further supports the assertion of no first-order phase transition from helical
states to SkXs. It should be pointed out that our SkXs, similar to those in
the literature, are not in perfect triangular lattice and $G_{6}-G_{2}$ is
smaller than 1.

Previous studies 49, 50, 51, 52 showed that stripe skyrmion structure follows
$m_{z}(x)=\frac{\sinh^{2}(L/2w)-\sinh^{2}(x/w)}{\sinh^{2}(L/2w)+\sinh^{2}(x/w)}$
for $-L/2<x<L/2$ and
$m_{z}(x)=-\frac{\sinh^{2}(L/2w)-\sinh^{2}(x/w-L/w)}{\sinh^{2}(L/2w)+\sinh^{2}(x/w-L/w)}$
for $L/2<x<3L/2$, where $L$ and $w$ are the stripe width and wall thickness.
$x=0$ refers to stripe center where $m_{z}=1$. Figures 3(c1)-(c4) are the
distributions of $m_{z}$ along lines labelled by the green ⓝ in Figs.
3(a2)-(a5). The solid curves are the fits to the theoretical spin profile with
$L=6.51\,$nm and $w=1.73\,$nm for (c1); $L=13.02\,$nm and $w=3.46\,$nm for
(c2); $L=19.53\,$nm and $w=5.19\,$nm for (c3); $L=26.04\,$nm and $w=6.92\,$nm
for (c4). $L$ agrees perfectly with $L=6.51A/D$ for $\kappa=2$. Data from
different stripes falling onto the same curve demonstrates that stripes,
building blocks of stripy pattern, are identical. Although the initial
configurations in Figs. 1, 2, and 3 are different, and the morphologies of the
mixtures of stripe and circular skyrmions are sensitive to the initial
configurations and their dynamical paths, the structures of SkXs do not depend
on the dynamical paths and depend only on the ratio of the stripe width and
skyrmion-skyrmion separation.

To further substantiate the assertion above, we carry out similar simulations
as those for Fig. 2, but with 100 randomly distributed nucleation domains as
the initial configuration. The average skyrmion-skyrmion separation is reduced
to $d_{ss}=\sqrt{400\times 400/100}=40\,$nm. Thus, we should expect the
transformation from maze to the mixture of stripe and circular skyrmions, and
to a SkX at smaller $A/D$ than those in Fig. 2. As shown in Fig. 4, this is
indeed the cases. (a1) is the initial configuration of randomly distributed
100 circular nucleation domains. (a2) is the stable structure of Film 1 of
$A/D=1\,$nm and $\kappa=2$ after $20\,$ns evolution. Clearly, it is a maze of
100 skyrmions, evident from $Q=100$ within a few picoseconds of evolution, and
individual skyrmion is an irregular ramified stripe, or a curved or short
stripe. After consecutive changes of $A$ and $K$ to simulate Films 3 to 10
given in Tab. 1 every $20\,$ns, Figs. 4(a3)-(a9) are their stable structures
at $t=40\,$ns; $60\,$ns; $80\,$ns; $100\,$ns; $120\,$ns; $140\,$ns; and
$160\,$ns respectively for $A/D=1.5\,$nm (a3); $A/D=2\,$nm (a4); $2.5\,$nm
(a5); $3\,$nm (a6), 3.5 nm (a7), 4 nm (a8), and 4.5 nm (a9) with stripe width
$L=9.76\,$nm, $13.02\,$nm, $16.27\,$nm, $19.53\,$nm, $22.78\,$nm, $26.04\,$nm,
$29.29\,$nm, $32.55\,$nm. Indeed, 100 skyrmions form SkXs when $A/D>3\,$nm as
shown in Figs. 4(a7)-(a9), smaller than $4\,$nm in Fig. 2. For $A/D=3\,$nm,
almost all skyrmions are circular except 9 very short stripe skyrmions (a6).
Figure 4(b) is the time dependences of total energy $E$ and $Q$ to show the
steady states. $Q(t)$ demonstrates clearly the topological feature of the
quantity and its resilience to environmental changes.

Figure 4: (a1) The initial configuration with 100 randomly distributed
nucleation domains of $m_{z}=1$ in the films. The diameter of each domain is
$6\,$nm. (a2)-(a9) The stable structures for $A/D=1\,$nm at $t=20\,$ns (a2);
$1.5\,$nm at $t=40\,$ns (a3); $2\,$nm at $t=60\,$ns (a4); $2.5\,$nm at
$t=80\,$ns (a5); $3\,$nm at $t=100\,$ns (a6); $3.5\,$nm at $t=120\,$ns (a7);
$4\,$nm at $t=140\,$ns (a8); $4.5\,$nm at $t=160\,$ns (a9). $\kappa=2$ is for
all films. (b) The time-dependences of skyrmion number $Q$ (the right y-axis
and the blue curve) and the total energy $E$ (the left y-axis and the red
curves). Time $t$ is in the linear scale for $E$ (the lower x-axis) and in in
the logarithmic scale (the top x-axis) for $Q$. $E$ for structures in
(a2)-(a9) are marked by cross. $Q$ reaches 100 within 1 ps and remains this
value even when model parameters are tuned. $E$ approaches its minimal value
in nanoseconds after model parameters are tuned.

Our results can either be experimentally tested or be used to explain
experiments because magnetic material parameters depend on the temperature.
Experiments 58, 59 showed that $A$ and $D$ depends on $M$ very differently
although mean field theories predict the same dependence of $A$ and $D$ on
$M$. Of course, both $L$ and $\kappa$ change with the temperature in realities
58, 59, in contrast to the assumption of constant $\kappa$ here. Since
skyrmion structures are determined by the ratio of $L$ and skyrmion-skyrmion
separation, not on $\kappa$ as long as it supports condensed skyrmion phase
($\kappa>1$) such that whether $\kappa$ is a constant or not is not essential.
Beside the temperature, other external knobs such as electrical fields,
magnetic fields, as well as strains can also change $L$ 60, 61, 62, 63, 64.

Physics reported here occurs in the bulk, and should not depend on the
boundary conditions. This is also what observed in experiments 54. However, it
is known that helical spin structures ended at a sample edge can be very
complicated. Skyrmion numbers of such spin textures are usually not integer.
In the Supporting Information, we discuss possible features of textures at
sample edges.

Looking forward, our assertion of no first-order phase transition may not be
applicable to SkX formation from the closely packed isolated skyrmions. The
similarities and differences between the two types of SkXs should be an
interesting issue that deserves a careful study. SkXs from a collection of
isolated skyrmions due to the skyrmion-skyrmion interaction has many
similarities with the formation of atomic crystals through atom-atom
interactions. It is well known that the melting or crystallization of an
atomic crystal to or from a liquid or an amorphous state is the first-order
phase transition involving latent heat. Of course, unlike real crystals in
which atoms cannot be destroyed, skyrmion number does not need to be conserved
in SkX formation.

In conclusion, a collection of skyrmions transforms into a SkX as the width of
stripe skyrmions becomes larger than the skyrmion-skyrmion separation, no
matter whether the initial configuration of skyrmions are ordered or randomly
distributed in a film and whether the width increases gradually or suddenly.
Since the SkX comes from the skyrmion-skyrmion repulsion, it explains why
individual skyrmion in a SkX is not a truly circular spin texture. Our finding
does not support the notion of the first-order phase transition for SkX
formation, but current results cannot rule out the second-order or higher-
order phase transitions.

##  1 Methods

Our model is a ferromagnetic film of thickness $d$ in $xy$ plane with the
total magnetic energy $E$ consisting of the exchange energy $E_{\mathrm{ex}}$,
the DMI energy $E_{\mathrm{DM}}$, the anisotropy energy $E_{\mathrm{an}}$, and
the magnetic dipolar energy $E_{\mathrm{d}}$

$E=E_{\mathrm{ex}}+E_{\mathrm{DM}}+E_{\mathrm{an}}+E_{\mathrm{d}},$ (4)

where $E_{\mathrm{ex}}=Ad\iint|\nabla\mathbf{m}|^{2}\mathrm{d}S$,
$E_{\mathrm{DM}}=Dd\iint[m_{z}\nabla\cdot\mathbf{m}-(\mathbf{m}\cdot\nabla)m_{z}]\mathrm{d}S$,
$E_{\mathrm{an}}=K_{u}d\iint(1-m_{z}^{2})\mathrm{d}S$, and
$E_{\mathrm{d}}=\mu_{0}M_{\rm s}d\iint\mathbf{H}_{\rm
d}\cdot\mathbf{m}\mathrm{d}S$. $M_{\rm s}$ and $\mathbf{m}$ are the saturation
magnetization and the unit direction of $\mathbf{M}$, respectively. $K_{u}$ is
the perpendicular magneto-crystalline anisotropy and $\mu_{0}$ is the vacuum
permeability. Magnetization dynamics is governed by the Landau-Lifshitz-
Gilbert (LLG) equation,

$\frac{\partial\mathbf{m}}{\partial t}=-\gamma\mathbf{m}\times\mathbf{H}_{\rm
eff}+\alpha\mathbf{m}\times\frac{\partial\mathbf{m}}{\partial t},$ (5)

where $\gamma$ and $\alpha$ are gyromagnetic ratio and the Gilbert damping
constant, respectively. $\mathbf{H}_{\rm
eff}=\frac{2A}{\mu_{0}M_{s}}\nabla^{2}\mathbf{m}+\frac{2K_{u}}{\mu_{0}M_{s}}m_{z}\hat{z}+\mathbf{H}_{\rm
d}+\mathbf{H}_{\rm DM}+\mathbf{h}$ is the effective field including the
exchange field, crystalline anisotropy field, demagnetization field, DMI
field, and a temperature-induced random magnetic field of magnitude
$h=\sqrt{2\alpha k_{B}T/(\mu_{0}M_{s}\gamma\Delta V\Delta t)}$ where $\Delta
V$, $\Delta t$, and $T$ are the cell volume, time step, and the temperature,
respectively. Demagnetization field can be included in the effective
anisotropy $K=K_{u}-\mu_{0}M_{s}^{2}/2$ when the film thickness $d$ is much
smaller than the exchange length 13. The thickness of our samples are
$0.5\,$nm which is indeed much smaller than the exchange length of
$\sqrt{2A/(\mu_{0}M_{s}^{2})}>60\,$nm, where $\mu_{0}=4\pi\times 10^{-7}\,{\rm
N/A^{2}}$ is vacuum permeability.

In the absence of energy sources such as an electric current and the heat
bath, the LLG equation describes a dissipative system whose energy can only
decrease 65. Thus, the steady state solutions of LLG equation with proper
initial magnetization distributions are the stable/metastable spin textures of
Eqs. (4). This is an efficient way to find the structure of a skyrmion
condensate. Mumax3 package 53 is used to search various metastable spin
textures in films of $400\,\rm nm\times\,400\,nm\times\,0.5\,nm$ with periodic
boundary conditions in both $x$ and $y$-directions. The mesh size is $1\,\rm
nm\times\,1\,nm\times\,0.5\,nm$ so that no change in simulation results is
detected when a smaller mesh size is used. In this study, we consider 10
different films with the same $M_{s}=0.15{\,}\rm MA/m$ and $D=1\,\rm
mJ/m^{2}$. $A$ and $K_{u}$ can vary and are used to tune the stripe width
while $\kappa=2$ are the same for all 10 films. The values of $A$ and $K_{u}$,
as well as the effective $K$ are listed in Tab. 1. $A/D$ of the 10 films
varies from 1 nm to 5.5 nm. Also a large $\alpha=0.5$ is used to speed up our
simulations.

Table 1: The value of material parameters, $A$, $K_{u}$, and $K$ for 10 different films.  Films | $A(\,\rm pJ/m\,)$ | $K_{u}(\,\rm MJ/m^{3}\,)$ | $K(\,\rm MJ/m^{3}\,)$
---|---|---|---
1 | 1.0 | 0.322 | 0.308
2 | 1.5 | 0.220 | 0.205
3 | 2.0 | 0.168 | 0.154
4 | 2.5 | 0.137 | 0.122
5 | 3.0 | 0.116 | 0.102
6 | 3.5 | 0.102 | 0.088
7 | 4.0 | 0.091 | 0.077
8 | 4.5 | 0.082 | 0.068
9 | 5.0 | 0.075 | 0.061
10 | 5.5 | 0.070 | 0.056

##  2 ASSOCIATED CONTENTS

Supporting Information

The Supporting Information is available free of charge on the ACS Publication
website at DOI:

  * •

Supporting Information: Discussion of spin structures at sample boundary,
damping dependence of metastable state, and determination of the Curie
temperatures. (PDF)

##  3 AUTHOR INFORMATION

Corresponding Author
*E-mail<EMAIL_ADDRESS>Phone: +852-23587488. Fax: +852-23581652.

Notes
The authors declare no competing financial interest.

##  4 ACKNOWLEDGMENTS

This work is supported by the National Key Research and Development Program of
China (grant No. 2020YFA0309600), the NSFC (grant No. 11974296), and Hong Kong
RGC Grants (No. 16300522, 16301619 and 16302321).

## References

  * 1 Bogdanov, A. N.; Rößler, U. K. Chiral Symmetry Breaking in Magnetic Thin Films and Multilayers. Phys. Rev. Lett. 2001, 87 (3), 037203. DOI:10.1103/PhysRevLett.87.037203
  * 2 Rößler, U. K.; Bogdanov, A. N.; Pfleiderer, C. Spontaneous skyrmion ground states in magnetic metals. Nature 2006, 442 (7104), 797-801. DOI:10.1038/nature05056
  * 3 Mühlbauer, S.; Binz, B.; Jonietz, F.; Pfeiderer, C.; Rosch, A.; Neubauer, A.; Georgii, G.; Böni, P. Skyrmion lattice in a chiral magnet. Science 2009, 323 (5916), 915-919. DOI:10.1126/science.1166767
  * 4 Yu, X. Z.; Onose, Y.; Kanazawa, N.; Park, J. H.; Han, J. H.; Matsui, Y.; Nagaosa, N.; Tokura, Y. Real-space observation of a two-dimensional skyrmion crystal. Nature 2010, 465 (7300), 901-904. DOI:10.1038/nature09124
  * 5 Yu, X. Z.; Kanazawa, N.; Onose, Y.; Kimoto, K.; Zhang, W. Z.; Ishiwata, S.; Matsui, Y.; Tokura, Y. Near room-temperature formation of a skyrmion crystal in thin-films of the helimagnet FeGe. Nat. Mater. 2011, 10 (2), 106-109. DOI:10.1038/nmat2916
  * 6 Fert, A.; Cros, V.; Sampaio, J. Skyrmions on the track. Nat. Nanotech. 2013, 8 (3), 152-156. DOI:10.1038/nnano.2013.29
  * 7 Sampaio, J.; Cros, V.; Rohart, S.; Thiaville, A.; Fert, A. Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures. Nat. Nanotech 2013, 8 (11), 839-844. DOI:10.1038/nnano.2013.210
  * 8 Nagaosa, N.; Tokura, Y. Topological properties and dynamics of magnetic skyrmions. Nat. Nanotech. 2013, 8 (12), 899-911. DOI:10.1038/nnano.2013.243
  * 9 Onose, Y.; Okamura, Y.; Seki, S.; Ishiwata, S.; Tokura, Y. Observation of magnetic excitations of skyrmion crystal in a helimagnetic insulator $\mathrm{Cu}_{2}\mathrm{OSeO}_{3}$. Phys. Rev. Lett. 2012, 109 (3), 037603. DOI:10.1103/PhysRevLett.109.037603
  * 10 Park, H. S.; Yu, X.; Aizawa, S.; Tanigaki, T.; Akashi, T.; Takahashi, Y.; Matsuda, T.; Kanazawa, N.; Onose, Y.; Shindo, D.; Tonomura, A.; Tokura, Y. Observation of the magnetic flux and three-dimensional structure of skyrmion lattices by electron holography. Nat. Nanotech. 2014, 9 (5), 337-342. DOI:10.1038/nnano.2014.52
  * 11 Krause, S.; Wiesendanger, R. Spintronics: Skyrmionics gets hot. Nat. Mater. 2016, 15 (5), 493-494. DOI:10.1038/nmat4615
  * 12 Rohart, S.; Thiaville, A. Skyrmion confinement in ultrathin film nanostructures in the presence of Dzyaloshinskii-Moriya interaction. Phys. Rev. B 2013, 88 (18), 184422. DOI:10.1103/PhysRevB.88.184422
  * 13 Wang, X. S.; Yuan, H. Y.; Wang, X. R. A theory on skyrmion size. Commun. Phys. 2018, 1 (1), 31. DOI:10.1038/s42005-018-0029-0
  * 14 Gong, X.; Yuan, H. Y.; Wang, X. R. Current-driven skyrmion motion in granular films. Phys. Rev. B 2020, 101 (6), 064421. DOI:10.1103/PhysRevB.101.064421
  * 15 Gong, X.; Jing, K. Y.; Lu, J.; Wang, X. R. Skyrmion pinning by disk-shaped defects. Phys. Rev. B 2022, 105 (9), 094437. DOI:10.1103/PhysRevB.105.094437
  * 16 Lenov, A. O.; Monchesky, T. L.; Romming, N.; Kubetzka, A.; Bogdanov, A. N.; Wiesendanger, R. The properties of isolated chiral skyrmions in thin magnetic films. New J. Phys. 2016, 18 (6), 065003. DOI:10.1088/1367-2630/18/6/065003
  * 17 Braun, H.-B. Fluctuations and instabilities of ferromagnetic domain-wall pairs in an external magnetic field. Phys. Rev. B 1994, 50 (22), 16485. DOI:10.1103/PhysRevB.50.16485
  * 18 Siemens, A.; Zhang, Y.; Hagemeister, J.; Vedmedenko, E. Y.; Wiesendanger, R. Minimal radius of magnetic skyrmions: statics and dynamics. New. J. Phys. 2016, 18 (4), 045021. DOI:10.1088/1367-2630/18/4/045021
  * 19 Simon, E.; Palotás, K.; Rózsa, L.; Udvardi, L.; Szunyogh, L. Formation of magnetic skyrmions with tunable properties in PdFe bilayer deposited on Ir(111). Phys. Rev. B, 2014, 90 (9), 094410. DOI:10.1103/PhysRevB.90.094410
  * 20 Karhu, E. A.; Rößler, U. K.; Bogdanov, A. N.; Kahwaji, S.; Kirby, B. J.; Fritzsche, H.; Robertson, M. D.; Majkrzak, C. F.; Monchesky, T. L. Chiral modulations and reorientation effects in MnSi thin films. Phys. Rev. B 2012, 85 (9), 094429. DOI:10.1103/PhysRevB.85.094429
  * 21 Iwasaki, J.; Mochizuki. M.; Nagaosa, N. Universal current-velocity relation of skyrmion motion in chiral magnets. Nat. Commun. 2013, 4 (1), 1463. DOI:10.1038/ncomms2442
  * 22 Lin, S.-Z.; Reichhardt, C.; Batista, C. D.; Saxena, A. Particle model for skyrmions in metallic chiral magnets: Dynamics, pinning, and creep. Phys. Rev. B 2013, 87 (21), 214419. DOI:10.1103/PhysRevB.87.214419
  * 23 Koshibae, W.; Nagaosa, N. Theory of current-driven skyrmions in disordered magnets. Sci. Rep. 2018, 8 (1), 6238. DOI:10.1038/s41598-018-24693-5
  * 24 Juge, R.; Je, S.-G.; Chaves, D. de S.; Pizzini, S.; Buda-Prejbeanu, L. D.; Aballe, L.; Foerster, M.; Locatelli, A.; Mentes, T. O.; Sala, A.; Maccherozzi, F.; Dhesi, S. S.; Auffret, S.; Gautier, E.; Gaudin, G.; Vogel, J.; Boulle, O. Magnetic skyrmions in confined geometries: Effect of the magnetic field and the disorder. J. Magn. Magn. Mater. 2018 (1), 455, 3. DOI:10.1038/s41598-018-24693-5
  * 25 Kim, J.-V.; Yoo, M.-W. Current-driven skyrmion dynamics in disordered films. Appl. Phys. Lett. 2017, 110 (13), 132404. DOI:10.1063/1.4979316
  * 26 Hoshino, S.; Nagaosa, N. Theory of the magnetic skyrmion glass. Phys. Rev. B 2018, 97 (2), 024413. DOI:10.1103/PhysRevB.97.024413
  * 27 Woo, S.; Song, K. M.; Zhang, X.; Zhou, Y,; Ezawa, M.; Liu, X.; Finizio, S.; Raabe, J.; Lee, N. J.; Kim, S.-I.; Park, S.-Y.; Kim, Y.; Kim, J.-Y.; Lee, D.; Lee, O.; Choi, J. W.; Min, B.-C.; Koo, H. C.; Chang, J. Current-driven dynamics and inhibition of the skyrmion Hall effect of ferrimagnetic skyrmions in GdFeCo films. Nat. Commun. 2018, 9 (1), 959. DOI:10.1038/s41467-018-03378-7
  * 28 Woo, S.; Litzius, K.; Krüger, B.; Im, M.-Y.; Caretta, L.; Richter, K.; Mann, M.; Krone, A.; Reeve, R. M.; Weigand, M.; Agrawal, P.; Lemesh, I.; Mawass, M.-A.; Fischer, P.; Kläui, M.; Beach, G. S. D. Observation of room-temperature magnetic skyrmions and their current-driven dynamics in ultrathin metallic ferromagnets. Nat. Mater. 2016, 15 (5), 501. DOI:10.1038/nmat4593
  * 29 Yuan, H. Y.; Wang, X. R. Skyrmion Creation and Manipulation by Nano-Second Current Pulses. Sci. Rep. 2016, 6 (1), 22638. DOI: 10.1038/srep22638
  * 30 Zhou, Y.; Ezawa, M. A. Reversible conversion between a skyrmion and a domain-wall pair in a junction geometry. Nat. Commun. 2014, 5 (1), 4652. DOI:10.1038/ncomms5652
  * 31 Jiang, W.; Upadhyaya, P.; Zhang, W.; Yu, G.; Jungfleisch, M. B.; Fradin, F. Y.; Pearson, J. E.; Tserkovnyak, Y.; Wang, K. L.; Heinonen, O.; te Velthuis, S. G. E.; Hoffmann, A. Blowing magnetic skyrmion bubbles. Science 2015, 349 (6245), 283-286. DOI:10.1126/science.aaa1442
  * 32 Li, J.; Tan, A.; Moon, K. W.; Doran, A.; Marcus, M. A.; Young, A. T.; Arenholz, E.; Ma, S.; Yang, R. F.; Hwang, C.; Qiu, Z. Q. Tailoring the topology of an artificial magnetic skyrmion. Nat. Commun. 2014, 5 (1), 4704. DOI:10.1038/ncomms5704
  * 33 Heinze, S.; Bergmann, K. von; Menzel, M.; Brede, J.; Kubetzka, A.; Wiesendanger, R.; Bihlmayer, G.; Blügel, S. Spontaneous atomic-scale Magnetic skyrmion lattice in two dimensions. Nat. Phys. 2011, 7 (9), 713-718. DOI:10.1038/nphys2045
  * 34 Dürrenfeld, P.; Xu, Y.; Åkerman, J.; Zhou, Y. Controlled skyrmion nucleation in extended magnetic layers using a nanocontact geometry. Phys. Rev. B 2017, 96 (5), 054430. DOI:10.1103/PhysRevB.96.054430
  * 35 Vidal-Silva, N.; Riveros, A.; Escrig, J. Stability of Néel skyrmions in ultra-thin nanodots considering Dzyaloshinskii-Moriya and dipolar interactions. J. Magn. Magn. Mater. 2017, 443, 116-123. DOI:10.1016/j.jmmm.2017.07.049
  * 36 Wilson, M. N.; Butenko, A. B.; Bogdanov, A. N.; Monchesky, T. L. Chiral skyrmions in cubic helimagnet films: The role of uniaxial anisotropy. Phys. Rev. B 2014, 89 (9), 094411. DOI:10.1103/PhysRevB.89.094411
  * 37 Romming, N.; Kubetzka, A.; Hanneken, C.; von Bergmann, K.; Wiesendanger, R. Field-Dependent Size and Shape of Single Magnetic Skyrmions. Phys. Rev. Lett. 2015, 114 (17), 177203. DOI:10.1103/PhysRevLett.114.177203
  * 38 Du, H.; Che, R.; Kong, L.; Zhao, X.; Jin, C.; Wang, C.; Yang, J.; Ning, W.; Li, R.; Jin, C.; Chen, X.; Zang, J.; Zhang, Y.; Tian, M. Edge-mediated skyrmion chain and its collective dynamics in a confined geometry. Nat. Commun. 2015, 6 (1), 8504. DOI:10.1038/ncomms9504
  * 39 Romming, N.; Hanneken, C.; Menzel, M.; Bickel, J. E.; Wolter, B.; von Bergmann, K.; Kubetzka, A.; Wiesendanger, R. Writing and deleting single magnetic skyrmions. Science 2013, 341 (6146), 636-639. DOI:10.1126/science.1240573
  * 40 Reichhardt, C.; Ray, D.; Olson Reichhardt, C. J. Collective Transport Properties of Driven Skyrmions with Random Disorder. Phys. Rev. Lett. 2015, 114 (21), 217202. DOI:10.1103/PhysRevLett.114.217202
  * 41 Yuan, H. Y.; Wang, X. S.; Yung, Man-Hong; Wang, X. R. Wiggling skyrmion propagation under parametric pumping. Phys. Rev. B 2019, 99 (1), 014428. DOI:10.1103/PhysRevB.99.014428
  * 42 W. Jiang,; Zhang, X.; Yu, G.; Zhang, W.; Wang, X.; Jungfleisch, M. B.; Pearson, J. E.; Cheng, X.; Heinonen, O.; Wang, K. L.; Zhou, Y.; Hoffmann, A.; te Velthuis, S. G. E. Direct observation of the skyrmion Hall effect. Nat. Phys. 2017, 13 (2), 162–169. DOI:10.1038/nphys3883
  * 43 Jaiswal, S.; Litzius, K.; Lemesh, I.; Büttner, F.; Finizio, S.; Raabe, J.; Weigand, M.; Lee, K.; Langer, J.; Ocker, B.; Jakob, G.; Beach. G. S. D.; Kläui, M. Investigation of the Dzyaloshinskii-Moriya interaction and room temperature skyrmions in W/CoFeB/MgO thin films and microwires. Appl. Phys. Lett. 2017, 111 (2), 022409. DOI:10.1063/1.4991360
  * 44 Litzius, K.; Lemesh, I.; Krüger, B.; Bassirian, P.; Caretta, L.; Richter, K.; Büttner, F.; Sato, K.; Tretiakov, O. A.; Förster, J.; Reeve, R. M.; Weigand, M.; Bykova, I.; Stoll, H.; Schtz, G.; Beach, G. S. D.; Kläui, M. Skyrmion Hall effect revealed by direct time-resolved X-ray microscopy. Nat. Phys. 2017, 13 (2), 170-175. DOI:10.1038/nphys4000
  * 45 Reichhardt, C.; Olson Reichhardt, C. J. Noise fluctuations and drive dependence of the skyrmion Hall effect in disordered systems. New J. Phys. 2016, 18 (9), 095005. DOI:10.1088/1367-2630/18/9/095005
  * 46 Sampaio, J.; Cros, V.; Rohart, S.; Thiaville A.; Fert, A. Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures. Nat. Nanotech. 2013, 8 (11), 839-844. DOI:10.1038/nnano.2013.210
  * 47 Back, C.; Cros, V.; Ebert, H.; Everschor-Sitte, K.; Fert, A.; Garst, M.; Ma, T. P.; Mankovsky, S.; Monchesky, T. L. ; Mostovoy, M.; Nagaos, N.; Parkin, S. S. P.; Pfleiderer, C.; Reyren, N.; Rosch, A.; Taguchi, Y.; Tokura, Y.; von Bergmann, K.; Zang, J. D. The 2020 skyrmionics roadmap. J. Phys. D: Appl. Phys. 2020, 53 (36), 363001. DOI:10.1088/1361-6463/ab8418
  * 48 Karube, K.; White, J. S.; Reynolds, N.; Gavilano, J. L.; Oike, H.; Kikkawa, A.; Kagawa, F.; Tokunaga, Y.; R$\rm\phi$nnow, H. M.; Tokura, Y.; Taguchi, Y. Robust metastable skyrmions and their triangular–square lattice structural transition in a high-temperature chiral magnet. Nat. Mat. 2016 (12), 15, 1237–1242. DOI:10.1038/nmat4752
  * 49 Wang, X. R.; Hu, X. C.; Wu, H. T. Stripe skyrmions and skyrmion crystals. Commun. Phys. 2021, 4 (1), 142. DOI:10.1038/s42005-021-00646-9
  * 50 Wu, H. T.; Hu, X. C.; Jing, K. Y.; Wang, X. R. Size and profile of skyrmions in skyrmion crystals. Commun. Phys. 2021, 4 (1), 210. DOI:10.1038/s42005-021-00716-y
  * 51 Wu, H. T.; Hu, X. C.; Wang, X. R. Nematic and smectic stripe phases and stripe-SkX transformations. Science China 2022, 65 (4), 247512. DOI:10.1007/s11433-021-1852-8
  * 52 Hu, X. C.; Wu, H. T.; Wang, X. R. A theory of skyrmion crystal formation. Nanoscale 2022, 14 (20), 7516-7529. DOI:10.1039/D2NR01300B
  * 53 Vansteenkiste, A.; Leliaert, J.; Dvornik, M.; Helsen, M.; Garcia- Sanchez F.; Waeyenberge, F. B. V. The design and verification of MuMax3. AIP Adv. 2014, 4 (10), 107133. DOI:10.1063/1.4899186
  * 54 Lin, T.; Wang, C.; Qiu, Z.; Chen C.; Xing, T.; Sun, L.; Liang J.; Wu, Y.; Shi, Z.; Lei, N. Magnetic triangular bubble lattices in bismuth-doped yttrium iron garnet. Chinese Physics B 2023, 32, 027505. doi:10.1088/1674-1056/aca604
  * 55 Kosterlitz, J. M.; Thouless, D. J. Ordering, metastability and phase transitions in two-dimensional systems. J. Phys. C 1973, 6, 1181. DOI:10.1088/0022-3719/6/7/010
  * 56 Halperin, B. I.; Nelson, D. R. Theory of two-dimensional melting. Phys. Rev. Lett. 1978, 41, (2) 121. DOI:10.1103/PhysRevLett.41.121
  * 57 Nelson, D. R.; Halperin, B. I. Dislocation-mediated melting in two dimensions. Phys. Rev. B 1979, 19 (5), 2457. DOI:10.1103/PhysRevB.19.2457
  * 58 Zhang, Y.; Kong, X.; Xu, G.; Jin, Y.; Jiang, C. J.; Chai, G. Z. Direct observation of the temperature dependence of the Dzyaloshinskii–Moriya interaction. J. Phys. D: Appl. Phys. 2022, 55 (19), 195304. DOI:10.1088/1361-6463/ac50cd
  * 59 Ham, W. S.; Ishibashi, M.; Kim, K.; Kim, S. ; Ono, T. Observation of temperature-dependent Dzyaloshinskii–Moriya interaction within the 50–300 K range. J. Phys. D: Appl. Phys. 2022, 61 (2), 02090. DOI:10.35848/1347-4065/ac46b0
  * 60 Schott, M.; Ranno, L.; Béa, H.; Baraduc, C.; Auffret, S.; Bernand-Mantel, A. Electric field control of interfacial Dzyaloshinskii-Moriya interaction in Pt/Co/Al$O_{x}$ thin films. J. Magn. Magn. Mater. 2021, 520, 167122. DOI:10.1016/j.jmmm.2020.167122
  * 61 Srivastava, T.; Schott, M.; Juge, R.; Křižáková, V.; Belmeguenai, M.; Roussigné, Y.; Bernand-Mantel, A.; Ranno, L.; Pizzini, S. ; Chérif, S.; Stashkevich, A.; Auffret, S.; Boulle, O.; Gaudin, G.; Chshiev, M.; Baraduc, C.; Béa, H. Large-voltage tuning of Dzyaloshinskii–Moriya interactions: A route toward dynamic control of skyrmion chirality. Nano Lett. 2018, 18 (8), 4871-4877. DOI:10.1021/acs.nanolett.8b01502
  * 62 Deger, C. Strain-enhanced dzyaloshinskii–moriya interaction at Co/Pt interfaces. Sci. Rep. 2020, 10 (1), 12314. DOI:10.1038/s41598-020-69360-w
  * 63 Gareeva, Z. V.; Zvezdin, A. K.; Kalyakin, L. A.; Gareev, T. T. Dzyaloshinskii–Moriya interaction, epitaxial strains, phase transitions in multiferroics with cycloidal structure. J. Magn. Magn. Mater. 2020, 515, 167255. DOI:10.1016/j.jmmm.2020.167255
  * 64 Udalov, O. G.; Beloborodov, I. S. Strain-dependent Dzyaloshinskii-Moriya interaction in a ferromagnet/heavy-metal bilayer. Phys. Rev. B 2020, 102 (13), 134422. DOI:10.1103/PhysRevB.102.134422
  * 65 Wang, X. R.; Yan, P.; Lu, J.; He, C. Magnetic field driven domain-wall propagation in magnetic nanowires. Ann. Phys. (NY) 2009, 324 (8), 1815-1820. DOI:10.1016/j.aop.2009.05.004


# Two-Step Knowledge Distillation for Tiny Speech Enhancement

###### Abstract

Tiny, causal models are crucial for embedded audio machine learning
applications. Model compression can be achieved via distilling knowledge from
a large teacher into a smaller student model. In this work, we propose a novel
two-step approach for tiny speech enhancement model distillation. In contrast
to the standard approach of a weighted mixture of distillation and supervised
losses, we firstly pre-train the student using only the knowledge distillation
(KD) objective, after which we switch to a fully supervised training regime.
We also propose a novel fine-grained similarity-preserving KD loss, which aims
to match the student’s intra-activation Gram matrices to that of the teacher.
Our method demonstrates broad improvements, but particularly shines in adverse
conditions including high compression and low signal to noise ratios (SNR),
yielding signal to distortion ratio gains of 0.9 dB and 1.1 dB, respectively,
at -5 dB input SNR and 63$\times$ compression compared to baseline.

Index Terms—  speech enhancement, knowledge distillation, tinyML, model
compression

##  1 Introduction

In recent years, deep neural network (DNN) models have become a common
approach to the speech enhancement (SE) problem, due to their performance [1,
2, 3]. However, large, powerful models are often unsuitable for resource-
constrained platforms, like hearing aids or wearables, because of their memory
footprint, computational latency, and power consumption [2, 4, 5, 6]. To meet
these constraints, audio TinyML research tends to focus on designing model
architectures with small numbers of parameters, using model compression
techniques to reduce the size of large models, or both [4, 5, 6, 7].

Pruning is a popular method for reducing the size of DNN models for SE [4, 5,
6, 8]. The goal of pruning is to remove weights least contributing to model
performance. In its simplest form, this can be performed post-training by
removing weights with the lowest magnitudes. Online pruning, where the model
is trained and pruned concurrently, builds on post-training pruning by
exposing the model to pruning errors while training, allowing it to adapt to
this form of compression noise [4]. Unstructured pruning of individual weights
can yield impressive model size reduction with little performance sacrifice,
but corresponding savings in computational throughput are not possible without
hardware support for sparse inference, which is unusual in embedded hardware.
Structured pruning of blocks of weights and/or neurons is often designed with
broader hardware compatibility in mind, but the performance drop tends to be
larger than for the unstructured case [6].

In contrast to pruning, knowledge distillation (KD) adopts a different
framework. The goal of KD is to utilize a strong pre-trained teacher model to
guide the training of the smaller student [9, 10, 11]. The underlying
assumption is that the pre-trained teacher network offers additional useful
context compared to the ground truth data by itself. Unlike pruning, this
process does not involve modifying the student network from its original dense
form, which reduces the complexity of the deployment process. In this work, we
focus on KD due to its above-outlined benefits over pruning.

KD methods have been applied to various classification tasks in the audio
domain [12, 13]. However, KD has not been extensively explored for causal low-
latency SE, which often requires tiny networks (sub-100k parameters) optimized
for low-resource wearable devices, such as hearing aids [5, 6]. So-called
response-based KD approaches use the pre-trained teacher model’s outputs to
train a student network [14, 15]. However, distillation can be further
facilitated by taking advantage of intermediate representations of the two
models, not just their outputs [10]. One common obstacle in such feature-based
KD is the dimensionality mismatch between teacher and student activations due
to the model size difference. To alleviate this issue, [16] proposed aligning
intermediate features, while [17] used attention maps to do so. The latter was
applied in the context of SE in [18] using considerably large, non-causal
student models intended for offline applications. In [19], the authors
addressed the dimensionality mismatch problem for the causal SE models by
using frame-level Similarity Preserving KD [20] (SPKD). SPKD captures the
similarity patterns between network activations for different training
examples and aims to match those patterns between the student and the frozen
pre-trained teacher models. The authors of [19] also introduced fusion blocks,
analogous to [21], to distill relationships between consecutive layers.

Here, we show that the efficacy of conventional KD methods is limited for
tiny, causal SE models. To improve distillation efficacy, we first extend the
method from [19] by computing SPKD for each bin of the latent representations,
corresponding to the time frame (as in [19]) but also the frequency bin of the
input, thus providing more resolution for KD loss optimization. The proposed
extension outperforms other similarity-based KD methods which we also explore.
Second, we hypothesize that matching a large teacher model might be
challenging for small student models and thus may lead to suboptimal
performance. Inspired by [22], we propose a novel two-step framework for
distilling tiny SE models. In the first step, the student is pre-trained using
only the KD criterion to match the activation patterns of the teacher, with no
additional ground truth supervision. The goal of this unsupervised KD pre-
training is to make the student similar to the teacher prior to the main
training. Then, the pre-trained student model is further optimized in a
supervised fashion and/or using KD routines. We find that pre-training using
the proposed SPKD method at the level of the individual bin of the latent
representation, followed by fully supervised training yields superior
performance compared to other distillation approaches utilizing weighted
mixtures of KD and supervised losses. We report the performance of our method
across various student model sizes, input mixture signal-to-noise ratios
(SNRs), and finally, assess the similarity between the activation patterns of
the teacher and distilled student.

Fig. 1: (a) Distillation process overview (b) Self-Similarity Gram matrices
computation. (c) Flow matrices computation ($\bigotimes$ denotes matrix
multiplication). Note that, for clarity, transpositions and matrix
multiplications are applied only to the last two dimensions of each tensor.

##  2 Methods

###  2.1 Model architecture

Our backbone architecture for the exploration of tiny SE KD is the
Convolutional Recurrent U-Net for SE (CRUSE) topology [7]. However, note that
the methodology developed here can, in principle, be applied to any other
architecture. The CRUSE model operates in the time-frequency domain and takes
power-law compressed log-mel spectrograms (LMS) as input. The LMS is obtained
by taking the magnitude of the complex short-time Fourier transform (STFT,
512/256 samples frame/hop size), processing it through a Mel-space filterbank
(80 bins, covering 50-8k Hz range) and finally compressing the result by
raising it to the power of $0.3$. The model output is a real-valued time-
frequency mask bounded within the range $(0,1)$ through the sigmoid activation
of the final block. The mask is applied to the noisy model input and
reconstituted into the time domain using the inverse STFT and the noisy input
phase.

The model comprises four encoder/decoder blocks and a bottleneck with grouped
GRU units (4 groups), reducing the computational complexity compared to a
conventional GRU layer with the same number of units [23]. The encoder/decoder
blocks are composed of 2D convolution/transpose convolution layers with (2, 3)
kernels (time, frequency) and (1, 2) strides, followed by cumulative layer
normalization [24] and leaky ReLU activation ($\alpha$ = 0.2). To further
reduce the model complexity, skip connections between the encoder and decoder
used in classic U-Net are replaced with 1x1 convolutions, whose outputs are
summed into the decoder inputs. We enforce the model’s frame-level causality
by using causal convolutions and causal cumulative layer norms. The total
algorithmic latency of the model is 32 ms (single STFT frame) [2].

In our experiments, both teacher and student are CRUSE models and their sizes
are adjusted by changing the number of units in each block. In particular, the
teacher uses {32, 64, 128, 192} encoder/decoder channels and 960 bottleneck
GRU units, resulting in 1.9M parameters, and 13.34 MOps/frame (i.e. the number
of operations required to process a single STFT frame). Our default student
uses {8, 16, 32, 32} encoder/decoder channels and 160 bottleneck GRU units
resulting in 62k parameters (3.3% of the teacher), and 0.84 MOps/frame (6.3%
of the teacher).

###  2.2 Self-similarity local knowledge distillation

Inspired by previous work [19, 20], we address the issue of dimensionality
mismatch between teacher and student models by computing similarity-based
distillation losses. The method captures and compares the relationship between
batch items at each layer output, between teacher and student (Fig. 1a, ${\cal
L}_{KD}^{local}$). We refer to this relationship as the self-similarity Gram
matrix $\textbf{G}_{x}$.

Self-similarity matrices (Fig. 1b) can be computed for an example network
latent activation X of shape [$b$, $c$, $t$, $f$], where $b$ \- batch size,
$c$ \- channel, $t$ \- activation width (corresponding to the input time
dimension), $f$ \- activation height (corresponding to the input frequency
dimension), as shown in Fig. 1b. The original implementation from [20]
involves reshaping X to [$b$, $ctf$] and matrix multiplying it by its
transpose $\textbf{X}^{T}$ to obtain the [$b$, $b$] symmetric self-similarity
matrix G. Analogously, this operation can be performed for each $t$ or $f$
dimension independently with resulting $\textbf{G}_{t/f}$ matrices of size
[$t/f$, $b$, $b$]. Such an increase in granularity improved the KD performance
in [19]. Here, we obtain even more detailed intra-activation Gram matrices by
considering each $(t,f)$ bin separately, resulting in the $\textbf{G}_{tf}$
self-similarity matrix with shape [$t$, $f$, $b$, $b$].

Finally, the local KD loss is computed using self-similarity matrices
$\textbf{G}_{x}$ of any kind $x$ obtained from teacher $T$ and student $S$ as:

${\cal
L}_{KD}^{local}=\frac{1}{b^{2}}\sum_{i}\left\lVert\textbf{G}_{x}^{T_{i}}-\textbf{G}_{x}^{S_{i}}\right\rVert^{2}_{F},$
(1)

where $i$ is the layer index and $\left\lVert\right\rVert^{2}_{F}$ is the
Frobenius $l_{2}$ norm.

###  2.3 Information flow knowledge distillation

The above-outlined local similarity losses can be extended to capture
relationships between activations of subsequent layers of the teacher and
student models (Fig. 1a, ${\cal L}_{KD}^{flow}$). The method is inspired by
the Flow of Solution Procedure (FSP) matrices introduced in [22] and aims to
not only match local similarity between the teacher and student in the
corresponding layers but also global inter-layer relations.

We propose two versions of flow matrices between layers $i$ and $j$ in our
model (Fig. 1c). The first one, $\textbf{G}_{t}^{i\rightarrow j}$, leverages
$\textbf{G}_{t}$ self-similarity matrices. Thereby each self-similarity block
shares the $t$-dimension and thus the interaction between the layers’ self-
similarity can be captured by performing matrix multiplication of
$\textbf{G}_{t}^{i}$ and transposed $\textbf{G}_{t}^{j}$ (both sized [$t$,
$b$, $b$]) for each time frame $t$.

The second version leverages $\textbf{G}_{tf}$ self-similarity matrices.
However, the $f$ dimension in our model changes for each block due to the
strided convolutions. To quantify the relationship between layers $i$ and $j$
of different dimensions we reshaped $\textbf{G}_{tf}^{i/j}$ to the size of
[$t$, $b$, $f_{i/j}$, $b$]. Then for each time-batch-item pair ($t$,$b$), we
obtain a [$f_{i/j}$, $b$] sub-matrix, which can be matrix multiplied with its
transpose to obtain the flow matrix $\textbf{G}_{tf}^{i\rightarrow j}$ of size
[$t$, $b$, $f_{i}$, $f_{j}$].

We define the loss similarly to Eq. 1 by comparing the teacher
$\textbf{G}^{T_{i\rightarrow j}}_{x}$ flow matrix with the student
$\textbf{G}^{S_{i\rightarrow j}}_{x}$ flow matrix, of the same kind $x$, for
every 2-layer-combination ($i$, $j$):

${\cal
L}_{KD}^{flow}=\frac{1}{b^{2}}\sum_{i}\sum_{j>i}\left\lVert\textbf{G}^{T_{i\rightarrow
j}}_{x}-\textbf{G}^{S_{i\rightarrow j}}_{x}\right\rVert^{2}_{F}$ (2)

###  2.4 Training objective and two-step KD

We use phase-sensitive spectrum approximation (PSA) [25], with clean speech as
a target, as the supervised portion ${\cal L}_{PSA}$ of the total loss. The KD
portion ${\cal L}_{KD}$ of the total loss does not use the ground truth
objective but instead, features obtained from the pre-trained, frozen teacher
model. In particular, ${\cal L}_{KD}$ can match model outputs (i.e. response
distillation, analogous to [15]), $\textbf{G}_{x}$ (${\cal L}_{KD}^{local}$),
or $\textbf{G}^{i\rightarrow j}_{x}$ (${\cal L}_{KD}^{flow}$) matrices
introduced in Sections 2.2 and 2.3, respectively. ${\cal L}_{PSA}$ and ${\cal
L}_{KD}$ are mixed using $\gamma$ coefficient to form the total loss:

${\cal L}=\gamma{\cal L}_{KD}+(1-\gamma){\cal L}_{PSA}$ (3)

Inspired by [22], we propose a two-step KD approach by separating the student
distillation process into two distinct parts. Step 1: In the first step,
$\gamma=1$ for a fixed set of epochs to solely minimize ${\cal L}_{KD}$. While
excluding the supervised ${\cal L}_{PSA}$ does not contribute to the optimal
objective performance, this step provides strong initial weights for further
student model training. Step 2: After this pre-training step, the student is
further optimized to maximize its objective performance using a fully
supervised loss by setting $\gamma=0$, (${\cal L}_{PSA}$ only) or a weighted
${\cal L}_{KD}$/${\cal L}_{PSA}$ loss obtained by setting $\gamma=0.5$. For
one-step KD using a weighted ${\cal L}_{KD}$/${\cal L}_{PSA}$ loss, we set
$\gamma=0.5$.

###  2.5 Training setup

During training, each epoch consists of 5,000 training steps, with each step
being a 32-item batch of 2-second-long audio clips. We used the Adam optimizer
with $6\cdot 10^{-5}$ learning rate. The teacher model is trained till
convergence to ensure the best performance for subsequent distillations. We
train each student model for a total of 400 epochs (2M steps, 35.5k+ hours of
audio). For student model pre-training (Step 1, Section 2.4), we use 100
epochs with $\gamma=1$, thus excluding supervised term ${\cal L}_{PSA}$ (Eq.
3).

##  3 Results

We use the dataset from the Interspeech 2020 Deep Noise Suppression (MS-DNS)
Challenge [1] for experimentation, which consists of 500+ hours of clean
speech and 100+ hours of noise, all mono clips sampled at 16 kHz. For model
training, we mix the speech and noise at various SNR levels, sampled from a
uniform distribution between -5 and 15 dB. We employ a LUFS-based SNR
calculation for more perceptually relevant mixtures and to de-emphasize the
effects of impulsive noises [26]. To evaluate the trained models we use the
non-reverberant evaluation set consisting of 150 clips of noisy speech samples
and their respective clean references.

We quantify the performance of each model via Signal-to-Distortion ratio (SDR)
[27], wide-band Perceptual Evaluation of Speech Quality (PESQ) [28], and
extended Short-Time Objective Intelligibility (eSTOI) [29]. We also report
scores obtained from DNS-MOS P.835 [30], a DNN mean opinion score (MOS)
estimator showing a good correlation with subjective ratings. All of our
results are reported as improvements over unprocessed noisy inputs ($\Delta$).

Table 1: One-step KD for tiny SE. Output: ${\cal L}_{KD}$ comparing teacher
and student outputs (similar to [15]). $\textbf{G}_{x}$: Feature-based ${\cal
L}_{KD}$ using self-similarity matrix of type $x$ (Fig. 1b). All models are
initialized with the same random weights and use $\gamma=0.5$ (Eq. 3).

|  |  |  | $\Delta\textbf{DNS-MOS}$
---|---|---|---|---
Model |  $\Delta\textbf{SDR}$ (dB) |  $\Delta\textbf{PESQ}$ (MOS) |  $\Delta\textbf{eSTOI}$ (%) | BAK | OVRL | SIG
Teacher | 8.65 | 1.25 | 10.07 | 1.44 | 0.69 | 0.06
Student | 6.34 | 0.75 | 5.82 | 1.27 | 0.55 | -0.02
Distillation |  |  |  |  |  |
Output [15] | 6.35 | 0.75 | 5.59 | 1.33 | 0.56 | -0.03
G [20] | 6.32 | 0.75 | 5.70 | 1.29 | 0.56 | -0.02
$\textbf{G}_{t}$ [19] | 6.50 | 0.77 | 5.95 | 1.33 | 0.55 | -0.04
$\textbf{G}_{f}$ | 6.47 | 0.74 | 6.03 | 1.29 | 0.56 | -0.02
$\textbf{G}_{tf}$ (ours) | 6.68 | 0.77 | 5.99 | 1.36 | 0.57 | -0.04

###  3.1 Self-similarity local KD

Table 1 shows the efficacy of local similarity-based one-step KD approaches
when training student models from scratch. Using teacher output as ${\cal
L}_{KD}$ in Eq. 3 or G similarity [20] does not improve, or even deteriorates
the student performance. $\textbf{G}_{t}$ similarity proposed in [19] provides
0.16 dB SDR improvement over the student alone, alongside the best PESQ score.
Our proposed time-frequency similarity calculation method $\textbf{G}_{tf}$
outperforms $\textbf{G}_{t}$ by doubling its SDR improvement (+0.34 dB, w.r.t.
the student alone) and increasing all other scores. This suggests that
increasing the granularity of the similarity matrix in the ${\cal
L}_{KD}^{local}$ calculation facilitates the KD process and overall improves
the performance of the distilled student model.

Table 2: Two-step KD. Step 1 \- Student pre-training using only ${\cal
L}_{KD}$ ($\gamma=1$) or no pre-training (None). Step 2  \- ${\cal L}_{PSA}$:
student training with only PSA loss ($\gamma=0$; supervised),
$\textbf{G}_{tf}$: Loss from Eq. 3 using $\textbf{G}_{tf}$-based ${\cal
L}_{KD}$ and $\gamma=0.5$ (best from Table 1).

|  |  |  | $\Delta\textbf{DNS-MOS}$
---|---|---|---|---
Model |  $\Delta\textbf{SDR}$ (dB) |  $\Delta\textbf{PESQ}$ (MOS) |  $\Delta\textbf{eSTOI}$ (%) | BAK | OVRL | SIG
Teacher | 8.65 | 1.25 | 10.07 | 1.44 | 0.69 | 0.06
Student | 6.34 | 0.75 | 5.82 | 1.27 | 0.55 | -0.02
Step 1 | Step 2 |  |  |  |  |  |
None | $\textbf{G}_{tf}$ | 6.68 | 0.77 | 5.99 | 1.36 | 0.57 | -0.04
| ${\cal L}_{PSA}$ | 6.46 | 0.78 | 6.07 | 1.29 | 0.56 | -0.02
$\textbf{G}^{i\rightarrow j}_{t}$ | $\textbf{G}_{tf}$ | 6.54 | 0.78 | 5.88 | 1.33 | 0.56 | -0.04
| ${\cal L}_{PSA}$ | 6.54 | 0.79 | 5.87 | 1.33 | 0.57 | -0.02
$\textbf{G}^{i\rightarrow j}_{tf}$ | $\textbf{G}_{tf}$ | 6.76 | 0.80 | 6.06 | 1.33 | 0.57 | -0.03
| ${\cal L}_{PSA}$ | 6.77 | 0.81 | 6.38 | 1.34 | 0.59 | -0.01
$\textbf{G}_{tf}$ | $\textbf{G}_{tf}$ | 6.75 | 0.80 | 6.34 | 1.32 | 0.57 | -0.02

###  3.2 Two-step KD

Table 2 presents the results of the proposed two-step distillation process
described in Section 2.4. We find that using time-preserving flow matrices
$\textbf{G}^{i\rightarrow j}_{t}$as the ${\cal L}_{KD}$ pre-training objective
(Step 1) yields comparable or worse performance than using local similarity
$\textbf{G}_{tf}$ with no pre-training. However, changing the pre-training
objective to $\textbf{G}^{i\rightarrow j}_{tf}$, which captures interactions
between latent representations in greater detail, yields improvement across
nearly all the metrics when paired with $\textbf{G}_{tf}$-based KD as Step 2.
Most interestingly, pre-training the student with $\textbf{G}_{tf}$ criterion
and continuing with only the supervised loss ${\cal L}_{PSA}$ provides
substantial improvements across all the metrics, especially SDR (+0.44 dB,
w.r.t. the student alone) and eSTOI (+0.56%), suggesting improved
intelligibility.

Fig. 2: Block-wise CKA similarity between students and teacher networks,
averaged over the MS-DNS test set. Mean(diag) and Mean(all) denote the average
similarity for the corresponding blocks (diagonal) or all the block
combinations, respectively.

We further investigate our best two-stage KD approach by performing Central
Kernel Alignment (CKA) [31] analysis. In principle, CKA allows us to compare
the similarity between activation patterns across different models in response
to a set of inputs. We use the entire evaluation dataset to probe the models
and compute CKA similarities for each pair of layers (averaged over all the
audio clips). Fig. 2-left presents CKA similarity between the teacher and
student trained independently. Fig. 2-middle compares the teacher to the
student pre-trained with $\textbf{G}_{tf}$ criterion (only Step 1). As
expected, the first step alone increases the similarity between the
corresponding teacher and student layers (diagonal). Finally, Fig. 2-right
shows the best student from Table 2, namely Step 1: $\textbf{G}_{tf}$ ${\cal
L}_{KD}$-only pre-training, and Step 2: fully-supervised ${\cal L}_{PSA}$. The
overall similarity to the teacher decreases but remains much higher than for
the student trained independently. This suggests that a brief pre-training
distillation ($\gamma=1$, no ${\cal L}_{PSA}$) allows the student to develop
its unique solution starting from strong prior knowledge inherited from the
teacher.

Table 3: Evaluation of the two-step KD approach on the MS-DNS dataset remixed
at fixed SNRs. Proposed: two-step KD using $\textbf{G}_{tf}$ pre-training
followed by fully-supervised training (best in Table 2).

|  |  |  |  | $\Delta\textbf{DNS-MOS}$
---|---|---|---|---|---
SNR (dB) | Model |  $\Delta\textbf{SDR}$ (dB) |  $\Delta\textbf{PESQ}$ (MOS) |  $\Delta\textbf{eSTOI}$ (%) | BAK | OVRL | SIG
| Teacher | 14.05 | 0.62 | 19.12 | 2.16 | 1.02 | 0.64
| Student | 10.82 | 0.30 | 10.07 | 1.86 | 0.79 | 0.51
\- 5 | Proposed | 11.73 | 0.35 | 11.61 | 1.98 | 0.81 | 0.47
| Teacher | 12.30 | 0.92 | 17.83 | 1.99 | 0.98 | 0.40
| Student | 9.65 | 0.49 | 10.56 | 1.75 | 0.75 | 0.26
0 | Proposed | 10.23 | 0.56 | 11.51 | 1.84 | 0.79 | 0.25
| Teacher | 10.27 | 1.21 | 13.98 | 1.65 | 0.78 | 0.02
| Student | 7.97 | 0.69 | 8.58 | 1.44 | 0.59 | -0.10
5 | Proposed | 8.43 | 0.76 | 9.32 | 1.51 | 0.62 | -0.09

###  3.3 Impact of the student model size and mixture SNR

The MS-DNS evaluation dataset consists of relatively high SNRs between 0 and
19 dB (mean 9.07 dB). To assess the SNR-dependent benefit of the proposed two-
step KD approach, we remix the entire evaluation set to obtain mixtures of the
same speech and noise clips but at fixed SNRs of -5, 0 and 5 dB. In Table. 3,
we observe the inverse relationship between the benefit of our approach and
SNR of the noisy mixtures. In particular, for -5 dB SNR mixtures, our KD
approach improves student performance by approximately 1 dB SDR, 1.5% eSTOI,
and 0.1 DNS-MOS BAK. This is an important observation as tiny SE models (here,
3.3% teacher size) tend to exhibit the most significant performance drop in
the low-SNR cases, compared to their larger counterparts [5].

In Table 4 we showcase the efficacy of the proposed KD framework for various
student sizes using the same teacher. We observe that the smaller the
downstream model, the larger benefit our KD method provides over the student
trained alone. In particular, for the 30k-parameter student ($\sim$1.5%
teacher size), the improvements are the largest with over 1 dB SDR, 0.1 PESQ,
and nearly 2% eSTOI. For model sizes above 200k ($\sim$15% teacher size), the
improvements start to plateau. These findings indicate that our method
provides the largest performance boost for the most resource-constrained
cases, usually deemed as the most challenging [4, 5, 6].

##  4 Conclusions

This work proposes a novel two-step KD protocol for distilling tiny, causal SE
models. No previous KD work has investigated this class of embedded-scale SE.
Our framework consists of two distinct steps: 1. Distilling the student model
using only KD objective using only our proposed fine-grained self-similarity
matrix $\textbf{G}_{tf}$ for computing distillation loss; 2. Training the
model obtained in Step 1 via supervised loss. Our results show that tiny SE
models distilled in this fashion perform better than KD methods utilizing
weighted loss between supervised and distillation objectives. Our experimental
evaluation shows that the proposed two-step KD provides the largest benefits
for low-SNR mixtures and smaller student models. Future work should explore
integrating the proposed two-step KD with pruning and/or quantization to
achieve SE models of even lower complexity and apply the method to other
audio-to-audio problems such as source separation, bandwidth extension, or
signal improvement.

Table 4: Impact of the student model size on the two-step KD performance. OPS:
number of operations per frame at inference time. Proposed: two-step KD using
$\textbf{G}_{tf}$ pre-training followed by fully-supervised training (best in
Table 2).

|  |  |  |  |  | $\Delta\textbf{DNS-MOS}$ |
---|---|---|---|---|---|---|---
Model |  Params / OPS (M) |  $\Delta\textbf{SDR}$ (dB) |  $\Delta\textbf{PESQ}$ (MOS) |  $\Delta\textbf{eSTOI}$ (%) | BAK | OVRL | SIG
Teacher | 1.9 / 13.34 | 8.65 | 1.25 | 10.07 | 1.44 | 0.69 | 0.06
Student |  | 4.42 | 0.50 | 2.59 | 1.21 | 0.47 | -0.07
Proposed | 0.03 / 0.42 | 5.52 | 0.61 | 4.55 | 1.18 | 0.47 | -0.05
Student |  | 6.34 | 0.75 | 5.82 | 1.27 | 0.55 | -0.02
Proposed | 0.06 / 0.84 | 6.77 | 0.81 | 6.38 | 1.34 | 0.59 | -0.01
Student |  | 7.24 | 0.93 | 7.53 | 1.38 | 0.62 | 0.00
Proposed | 0.24 / 2.48 | 7.60 | 0.97 | 7.71 | 1.41 | 0.64 | 0.01
Student |  | 7.51 | 0.99 | 7.97 | 1.39 | 0.63 | 0.01
Proposed | 0.35 / 3.08 | 7.54 | 1.01 | 8.22 | 1.38 | 0.64 | 0.02

## References

  * [1] Chandan K. A. Reddy et al.,  “The interspeech 2020 deep noise suppression challenge: Datasets, subjective testing framework, and challenge results,”  in Proc. Interspeech. ISCA, 2020, pp. 2492–2496.
  * [2] Zhong-Qiu Wang, Gordon Wichern, Shinji Watanabe, and Jonathan Le Roux,  “Stft-domain neural speech enhancement with very low algorithmic latency,”  IEEE TASLP, vol. 31, pp. 397–410, 2022.
  * [3] Bryce Irvin, Marko Stamenovic, Mikolaj Kegler, and Li-Chia Yang,  “Self-supervised learning for speech enhancement through synthesis,”  in ICASSP. IEEE, 2023, pp. 1–5.
  * [4] Jiasi Chen and Xukan Ran,  “Deep learning with edge computing: A review,”  Proceedings of the IEEE, vol. 107, no. 8, pp. 1655–1674, 2019.
  * [5] Igor Fedorov et al.,  “TinyLSTMs: Efficient neural speech enhancement for hearing aids,”  in Proc. Interspeech. ISCA, 2020, pp. 4054–4058.
  * [6] Marko Stamenovic, Nils Westhausen, Li-Chia Yang, Carl Jensen, and Alex Pawlicki,  “Weight, block or unit? exploring sparsity tradeoffs for speech enhancement on tiny neural accelerators,”  in NeurIPS ENLSP, 2021.
  * [7] Sebastian Braun, Hannes Gamper, Chandan KA Reddy, and Ivan Tashev,  “Towards efficient models for real-time deep noise suppression,”  in ICASSP. IEEE, 2021, pp. 656–660.
  * [8] Zeyuan Wei, Li Hao, and Xueliang Zhang,  “Model Compression by Iterative Pruning with Knowledge Distillation and Its Application to Speech Enhancement,”  in Proc. Interspeech. ISCA, 2022, pp. 941–945.
  * [9] Geoffrey Hinton, Oriol Vinyals, and Jeffrey Dean,  “Distilling the knowledge in a neural network,”  in NIPS Deep Learning and Representation Learning Workshop, 2015\.
  * [10] Jianping Gou, Baosheng Yu, Stephen J. Maybank, and Dacheng Tao,  “Knowledge distillation: A survey,”  International Journal of Computer Vision, vol. 129, no. 6, pp. 1789–1819, 2021.
  * [11] Chuanguang Yang, Xinqiang Yu, Zhulin An, and Yongjun Xu,  “Categories of response-based, feature-based, and relation-based knowledge distillation,”  in Advancements in Knowledge Distillation: Towards New Horizons of Intelligent Systems, pp. 1–32. Springer, 2023.
  * [12] Florian Schmid, Shahed Masoudian, Khaled Koutini, and Gerhard Widmer,  “Knowledge distillation from transformers for low-complexity acoustic scene classification,”  in DCASE workshop, 2022.
  * [13] Ji Won Yoon, Beom Jun Woo, Sunghwan Ahn, Hyeonseung Lee, and Nam Soo Kim,  “Inter-kd: Intermediate knowledge distillation for ctc-based automatic speech recognition,”  in SLT workshop. IEEE, 2023, pp. 280–286.
  * [14] Xiang Hao, Shixue Wen, Xiangdong Su, Yun Liu, Guanglai Gao, and Xiaofei Li,  “Sub-band knowledge distillation framework for speech enhancement,”  in Proc. Interspeech. ISCA, 2020, pp. 2687–2691.
  * [15] Sotaro Nakaoka, Li Li, Shota Inoue, and Shoji Makino,  “Teacher-student learning for low-latency online speech enhancement using wave-u-net,”  in ICASSP. IEEE, 2021, pp. 661–665.
  * [16] Adriana Romero, Nicolas Ballas, Samira Ebrahimi Kahou, Antoine Chassang, Carlo Gatta, and Yoshua Bengio,  “Fitnets: Hints for thin deep nets,”  arXiv preprint arXiv:1412.6550, 2014.
  * [17] Sergey Zagoruyko and Nikos Komodakis,  “Paying more attention to attention: Improving the performance of convolutional neural networks via attention transfer,”  in ICLR, 2017.
  * [18] Wooseok Shin, Hyun Joon Park, Jin Sob Kim, Byung Hoon Lee, and Sung Won Han,  “Multi-view attention transfer for efficient speech enhancement,”  in Proc. Interspeech. ISCA, 2022.
  * [19] Jiaming Cheng, Ruiyu Liang, Yue Xie, Li Zhao, Björn Schuller, Jie Jia, and Yiyuan Peng,  “Cross-Layer Similarity Knowledge Distillation for Speech Enhancement,”  in Proc. Interspeech. ISCA, 2022, pp. 926–930.
  * [20] Frederick Tung and Greg Mori,  “Similarity-preserving knowledge distillation,”  in CVPR. IEEE, 2019, pp. 1365–1374.
  * [21] Pengguang Chen, Shu Liu, Hengshuang Zhao, and Jiaya Jia,  “Distilling knowledge via knowledge review,”  in CVPR. IEEE, 2021, pp. 5008–5017.
  * [22] Junho Yim, Donggyu Joo, Jihoon Bae, and Junmo Kim,  “A gift from knowledge distillation: Fast optimization, network minimization and transfer learning,”  in CVPR. IEEE, 2017, pp. 7130–7138.
  * [23] Ke Tan and DeLiang Wang,  “Learning complex spectral mapping with gated convolutional recurrent networks for monaural speech enhancement,”  IEEE TASLP, vol. 28, pp. 380–390, 2020.
  * [24] Yi Luo and Nima Mesgarani,  “Conv-TasNet: Surpassing ideal time–frequency magnitude masking for speech separation,”  IEEE TASLP, vol. 27, no. 8, pp. 1256–1266, 2019.
  * [25] Hakan Erdogan, John R Hershey, Shinji Watanabe, and Jonathan Le Roux,  “Phase-sensitive and recognition-boosted speech separation using deep recurrent neural networks,”  in ICASSP. IEEE, 2015, pp. 708–712.
  * [26] ITU-R,  “ITU-R BS.1770-4: Algorithms to measure audio programme loudness and true-peak audio level,” 2015.
  * [27] Jonathan Le Roux, Scott Wisdom, Hakan Erdogan, and John R Hershey,  “Sdr–half-baked or well done?,”  in ICASSP. IEEE, 2019, pp. 626–630.
  * [28] ITU-R,  “P.862 : Perceptual evaluation of speech quality (PESQ): An objective method for end-to-end speech quality assessment of narrow-band telephone networks and speech codecs,” 2015.
  * [29] Jesper Jensen and Cees H Taal,  “An algorithm for predicting the intelligibility of speech masked by modulated noise maskers,”  IEEE TASLP, vol. 24, no. 11, pp. 2009–2022, 2016.
  * [30] Chandan K. A. Reddy, Vishak Gopal, and Ross Cutler,  “Dnsmos p.835: A non-intrusive perceptual objective speech quality metric to evaluate noise suppressors,”  in ICASSP. IEEE, 2022, pp. 886–890.
  * [31] Simon Kornblith, Mohammad Norouzi, Honglak Lee, and Geoffrey Hinton,  “Similarity of neural network representations revisited,”  in ICML, 2019, pp. 3519–3529.