# Spin-Current Instability at a Magnetic Domain Wall in a Ferromagnetic
Superfluid:
a Generation Mechanism of Eccentric Fractional Skyrmions

Hiromitsu Takeuchi<EMAIL_ADDRESS>http://hiromitsu-
takeuchi.appspot.com  Department of Physics and Nambu Yoichiro Institute of
Theoretical and Experimental Physics (NITEP),
Osaka City University, Osaka 558-8585, Japan

###### Abstract

Spinful superfluids of ultracold atoms are ideal for investigating the
intrinsic properties of spin current and texture because they are realized in
an isolated, nondissipative system free from impurities, dislocations, and
thermal fluctuations. This study theoretically reveals the impact of spin
current on a magnetic domain wall in spinful superfluids. An exact wall
solution is obtained in the ferromagnetic phase of a spin-1 Bose–Einstein
condensate with easy-axis anisotropy at zero temperature. The bosonic-
quasiparticle mechanics analytically show that the spin current along the wall
becomes unstable if the spin-current velocity exceeds the criteria, leading to
complicated situations because of the competition between transverse magnons
and ripplons. Our direct numerical simulation reveals that this system has a
mechanism to generate an eccentric fractional skyrmion, which has a fractional
topological charge, but its texture is not similar to that of a meron. This
mechanism is in contrast to the generation of conventional skyrmions in easy-
axis magnets [S. K. Kim and Y. Tserkovnyak, Phys. Rev. Lett. 119, 047202
(2017)]. The theoretical findings can be examined in the same situation as in
a recent experiment on ultracold atoms and indicate unexplored similar
phenomena in different physical systems, such as chiral superfluids and
superconductors, magnets, spintronics, and particle physics.

###### pacs:

Valid PACS appear here

Membranal topological defects, called domain walls (DWs), are formed when a
discrete symmetry is spontaneously broken in a phase transition. They occur in
every branch of physics, ranging from condensed matter fields, such as
magnetism hubert2008magnetic  and optics gilles2017polarization , to cosmology
and field theory 2006Vachaspati . DWs are also important in engineering
because they appear in electronic and spintronic devices and magnetic
materials, e.g., twisted nematic liquid crystals chandrasekhar_1992 ,
submicrometer ferromagnetic structures atkinson2003magnetic , and
ferromagnetic metals tatara1997resistivity .

DWs in fermionic superfluids and superconductors with internal orbital and
spin degrees of freedom have attracted increasing interest. Chiral symmetry
breaking has been confirmed in superfluid 3He-A PhysRevLett.109.215301 ;
ikegami2013chiral ; ikegami2015observation , and chiral domain structures have
been observed in a slab geometry PhysRevLett.120.205301 . Superfluid 3He-B is
also capable of domain formation in a slab PhysRevLett.122.085301 , and an
exotic DW that terminates on a half-quantum vortex as a nexus
volovikPhysRevResearch.2.023263 ; Zhang2020PhysRevResearch.2.043356  has been
realized in porous media makinen2019half . As fermionic quasiparticles form
anomalous bound states at a DW or interface Aoki2005_PhysRevLett.95.075301 ;
Nurakawa2009_PhysRevLett.103.155301 ; murakawa2010surface ; Okuda2012_JPCM ;
Zheng2016_PhysRevLett.117.195301 ; Zheng2017_PhysRevLett.118.065301  due to a
spatial modulation of the order parameter, the bound quasiparticles in
superfluids and superconductors cause unconventional responses of currents
matsumoto1999quasiparticle ; Ashby2009PhysRevB.79.224509 ;
Serban2010PhysRevLett.104.147001 ; Sauls2011PhysRevB.84.214509 ;
Bouhon2014PhysRevB.90.220511 ; higashitani2014spin  and magnetic fields
TanakaPhysRevB.72.140503 ; Yokoyama2011PhysRevLett.106.246601 ;
Higashitani2013PhysRevLett.110.175301 ; Asano2015PhysRevB.92.224508 .

In contrast, DWs in bosonic superfluids and the properties associated with
quasiparticle-bound states are less apparent. Recently, composite defects of
nematic-spin DWs and half-quantum vortices have been experimentally realized
through the symmetry breaking phase transition in spin-1 Bose–Einstein
condensates (BECs) of 27Na atoms kangPhysRevLett.122.095301 . The mass-current
field around a composite defect is theoretically shown to form an anomalous
elliptic structure in equilibrium TakeuchiPhysRevLett.126.195302 . Magnetic-
domain formation has been realized in a strongly ferromagnetic spinor BEC of
7Li atoms Huh2020PhysRevResearch.2.033471 . Spinor or multicomponent
superfluids are suitable systems for exploring the intrinsic effect of spin
currents because they are realized in isolated systems free from impurities
and dislocations. A spin current (countersuperflow) has been experimentally
shown to be intrinsically unstable without energy dissipation
Hamner2011PhysRevLett.106.065302 ; KimPhysRevLett.119.185302  when the spin-
current velocity exceeds the criteria Law2001PhysRevA.63.063612 ;
Zhu2015PhysRevA.91.023633 . While these spin-current instabilities occur in
the bulk without the involvement of DWs, the following fundamental questions
have not been answered yet: How do spin currents influence a magnetic DW? What
is the role of quasiparticle-bound states?

In this study, we theoretically investigate the properties of magnetic DWs
with spin currents in the ferromagnetic phase of spin-1 BECs. The DWs in this
system are classified into antiferromagnetic (AF)-core and broken axisymmetry
(BA)-core DWs according to the local magnetization at the wall (Fig. 1). The
bosonic-quasiparticle mechanics reveal that a spin current causes spin-
intrinsic instability in an AF-core DW, causing the condensation of the
transverse magnon in a bound state to form a BA-core DW. The instability
dynamics have several channels to different results owing to the competition
between the magnon condensation and Kelvin–Helmholtz instability (KHI).
Finally, our direct numerical simulation reveals that the instability at a BA-
core DW causes eccentric fractional skyrmions, which are topological
structures distinct from conventional skyrmions and merons.

Basic formulation.— Spin-1 BECs of dilute Bose gases at very low temperatures
are described in the mean field theory, where the boson field operator
$\hat{\Psi}_{m}^{(\dagger)}$ is replaced by the macroscopic wave function
$\Psi_{m}^{(*)}$ in the many-body Hamiltonian. Here, $\Psi_{m}({\bm{r}},t)$ is
the complex scalar field of the $\left|m\right>$ Zeeman ($m=0,\pm 1$)
component. The dynamics in a uniform system are described by the Lagrangian
function of $\vec{\Psi}=[\Psi_{+1},\Psi_{0},\Psi_{-1}]^{\rm T}$ in the
Gross–Pitaevskii model kawaguchi2012spinor  as ${\cal L}=\int
d^{3}x\left(\sum_{m}i\hbar\Psi_{m}^{*}\partial_{t}\Psi_{m}-{\cal G}\right)$,
with ${\cal
G}(\vec{\Psi})=\frac{\hbar^{2}}{2M}\sum_{m}({\bm{\nabla}}\Psi_{m}^{*})\cdot({\bm{\nabla}}\Psi_{m})+{\cal
U}$ and

$\displaystyle{\cal
U}=\frac{c_{n}}{2}n^{2}+\frac{c_{s}}{2}{\bm{s}}^{2}-\sum_{m}(\mu+mp-m^{2}q-U_{m})\left|\Psi_{m}\right|^{2}.$
(1)

Here, we used the total density $n=\sum_{m}\left|\Psi_{m}\right|^{2}$ and the
spin density ${\bm{s}}=[s_{x},s_{y},s_{z}]^{\rm
T}=\sum_{mm^{\prime}}\Psi_{m}^{*}\left(\check{\bm{\sigma}}\right)_{mm^{\prime}}\Psi_{m^{\prime}}$
with the spin-1 matrices
$\check{\bm{\sigma}}=[\check{\sigma}_{x},\check{\sigma}_{y},\check{\sigma}_{z}]^{\rm
T}$. The Lagrange multipliers $\mu$ and $p$ are associated with the
conservation of the particle number and longitudinal magnetization,
respectively. We neglect the external potential $U_{m}$, unless otherwise
noted. Spin-1 BECs have four different phases: polar (P), AF, ferromagnetic
(F), and BA phases kawaguchi2012spinor . Magnetic DWs are realized in the F
phase with ferromagnetic interaction $-c_{n}<c_{s}<0$ and a negative quadratic
Zeeman shift $q<0$.

The ground state $\vec{\Psi}=\vec{\Phi}_{\rm F\pm}$ in the F phase for
$p\gtrless 0$ has a magnetization $s_{z}=\pm n$ and is represented as the F±
state, $\Psi_{\pm 1}=\sqrt{\frac{\mu-q\pm p}{c_{n}+c_{s}}}e^{i\theta_{\rm G}}$
and $\Psi_{\mp 1}=\Psi_{0}=0$. Here, $\theta_{\rm G}=const.$ is the global
phase according to the $U(1)$-symmetry breaking. The energy densities ${\cal
G}(\vec{\Phi}_{\rm F\pm})={\cal U}_{\rm F}=-\frac{1}{2}\frac{(\mu-q\pm
p)^{2}}{c_{n}+c_{s}}$ of these states have the same value for $p=0$,
corresponding to the spontaneous breaking of the discrete symmetry with
respect to the spin inversion $s_{z}\leftrightarrow-s_{z}$. Magnetic DWs are
topological defects associated with discrete symmetry breaking, in addition to
quantized vortices with $U(1)$-symmetry breaking. The DW separates the domains
in the two ground states with opposite magnetizations. The bulk ordered state
is suppressed in the core of the DW, in which the macroscopic wave functions
vary continuously between the two states. The core of a topological defect in
spin-1 BECs has been revealed to exhibit complicated states PhysRevA.86.013613
; PhysRevLett.112.075301 ; PhysRevA.93.033633 ; weiss2019controlled ;
liu2020phase ; underwood2020properties ; PhysRevLett.125.170401 ;
PhysRevLett.125.030402 ; PhysRevResearch.3.L012003 ;
TakeuchiPhysRevLett.126.195302 ; 2021Takeuchi_PhysRevA.104.013316 ;
Katsimiga_2021 . Similarly, there are different possibilities for the local
states at the wall in our system.

To approach the problem systematically, we introduce a general rule on the
mass and spin currents along a DW. The rule states that the current velocity
${\bm{V}}_{m}~{}(\bot\hat{\bm{x}})$ of the $m$ component parallel to a flat DW
must satisfy

$\displaystyle{\bm{V}}_{m}={\bm{V}}_{+}+m{\bm{V}}_{-}.$ (2)

with ${\bm{V}}_{\pm}=\frac{1}{2}({\bm{V}}_{+1}\pm{\bm{V}}_{-1})$. The center-
of-mass velocity ${\bm{V}}_{+}$ is not essential for the considered spin-
current instability. The current velocity rule (2) is obtained by substituting
$\Psi_{m}({\bm{r}},t)=\Phi_{m}({\bm{r}})=f_{m}(x)e^{i\frac{\hbar}{M}{\bm{V}}_{m}\cdot{\bm{r}}}$
into the equation of motion SM1 . Here, the real field $f_{m}$ obeys the
coupled equations

$\displaystyle 0=H_{m}f_{m}+c_{s}(f_{+1}+f_{-1})^{2-m^{2}}f_{0}^{1+m^{2}}$ (3)

with $H_{m}=-\frac{\hbar^{2}}{2M}\frac{{\rm d}^{2}}{{\rm
d}x^{2}}+c_{n}n+mc_{s}s_{z}-\tilde{\mu}+m\tilde{p}+m^{2}\tilde{q}$ and the
hydrostatic variables 2021Takeuchi_PhysRevA.104.013316 ,
$\tilde{\mu}=\mu-\frac{1}{2}M{\bm{V}}_{+}^{2}$,
$\tilde{p}=p-M{\bm{V}}_{+}\cdot{\bm{V}}_{-}$, and
$\tilde{q}=q+\frac{1}{2}M{\bm{V}}_{-}^{2}$. As the hydrostatic pressure in a
domain with $s_{z}=\pm n$ is given by
$P_{\pm}=\frac{1}{2}\frac{\left(\tilde{\mu}-\tilde{q}\pm\tilde{p}\right)^{2}}{c_{n}+c_{s}}$,
the condition $\tilde{p}=0$ is satisfied for a flat DW in the pressure
equilibrium with $P_{+}=P_{-}$:

$\displaystyle n(x\to\pm\infty)\to n_{\rm
F}=\frac{\tilde{\mu}-\tilde{q}}{c_{n}+c_{s}},$ (4)

or equivalently $s_{z}(x\to\pm\infty)\to\pm n_{\rm F}$.

Figure 1:  Schematic profiles of AF-core and BA-core DWs. The curves show the
density $f_{m}^{2}~{}(m=0,\pm 1)$ and amplitude of the spin-current velocity
$v_{\rm spin}$ along the DWs for ${\bm{V}}_{+}=0$. The arrows represent the
spin texture.

Under the boundary condition (4) for $-1<\frac{c_{s}}{c_{n}}<0$ and
$\tilde{q}<0$, we obtain two types of solutions, classified by the local
ordered states at the center $x=0$: the AF-core DW with the local AF state
[$f_{+1}(0)=f_{-1}(0)$ and $f_{0}(0)=0$], and the BA-core DW with the local BA
state [$f_{+1}(0)=f_{-1}(0)$ and $f_{0}(0)\neq 0$] (Fig. 1). The spin-current
velocity ${\bm{v}}_{\rm spin}\equiv\frac{{\bm{j}}_{\rm spin}}{n}$, defined by
the spin-current density ${\bm{j}}_{\rm spin}=\sum_{m}m{\bm{j}}_{m}$, where
${\bm{j}}_{m}=\frac{\hbar}{M}{\rm Im}(\Psi_{m}^{*}{\bm{\nabla}}\Psi_{m})$, is
computed as ${\bm{v}}_{\rm
spin}=\frac{f_{+1}^{2}-f_{-1}^{2}}{n}{\bm{V}}_{+}+\frac{f_{+1}^{2}+f_{-1}^{2}}{n}{\bm{V}}_{-}$,
which reduces to ${\bm{V}}_{-}$ at the center of an AF-core DW with
$|f_{+1}(0)|=|f_{-1}(0)|$. The energy density at the wall is generally higher
than that in the bulk; thus, $n(x=0)$ is smaller than $n_{\rm F}$.

The theoretical correspondence between binary BECs and spin-1 BECs with
$\Psi_{0}=0$ can be efficiently utilized to conduct our analytical
computations. According to the correspondence, the DW solution in binary BECs
PhysRevA.91.033615  gives an exact solution of the AF-core DW for
$c_{s}=-c_{n}/2$ SM1 ,

$\displaystyle f_{\pm 1}=f_{\pm 1}^{\rm ex}=\frac{\sqrt{n_{\rm
F}}}{2}\left[1\pm\tanh\left(\frac{x}{\xi}\right)\right]$ (5)

with $\xi=\frac{\hbar}{\sqrt{M(\tilde{\mu}-\tilde{q})}}$. This solution is
practical because the spin-1 BECs of 7Li atoms have $c_{s}\approx-c_{n}/2$
Huh2020PhysRevResearch.2.033471 . An exact DW solution is also available for
the BA phase Yu2021PhysRevResearch.3.023043 ; Yu2021anomalous .

Bosonic-quasiparticle mechanics— We first discuss the spin-current instability
at an AF-core DW analytically. The linear stability is investigated by solving
the eigenvalue problem, derived by linearizing the equation of motion with
respect to the Bogoliubov modes
$\delta\Psi_{m}({\bm{r}},t)=\Psi_{m}({\bm{r}},t)-\Phi_{m}({\bm{r}})=e^{i\Theta_{m}}\left[u_{m}(x)e^{i{\bm{k}}\cdot{\bm{r}}-i\omega
t}-v_{m}(x)^{*}e^{-i{\bm{k}}\cdot{\bm{r}}+i\omega^{*}t}\right]$ with
${\bm{k}}\bot\hat{\bm{x}}$. As the instability is induced by modes with
nonpositive excitation energy, we investigate the behavior of low-energy
excitations. For $f_{0}=0$, we have two independent eigenvalue equations for
the vector fields $\vec{u}_{0}=(u_{0},v_{0})^{\rm T}$ and
$(u_{+1},u_{-1},v_{+1},v_{-1})^{\rm T}$. The problem with the latter is the
same as the quantum KHI in segregated binary BECs
Takeuchi2010PhysRevB.81.094517 , where the lowest-energy excitations of
ripplons, the quanta of ripple waves on a DW, are described by the low-energy
effective theory SM1 .

The distinction from binary BECs is caused by the excitations of
$\vec{u}_{0}$, called transverse magnons. The magnons obey the eigenvalue
equation

$\displaystyle\hbar\tilde{\omega}\vec{u}_{0}=\begin{pmatrix}h_{0}+h_{\rm
p}&h_{\rm p}\\\ -h_{\rm p}&-h_{0}-h_{\rm p}\end{pmatrix}\vec{u}_{0}$ (6)

with $\tilde{\omega}=\omega-{\bm{k}}\cdot{\bm{V}}_{+}$,
$h_{0}=\frac{\hbar^{2}{\bm{k}}^{2}}{2M}+H_{0}+c_{s}(f_{+1}+f_{-1})^{2}$, and
$h_{\rm p}=-2c_{s}f_{+1}f_{-1}$. We do not explicitly consider the center-of-
mass motion causing the term $-{\bm{k}}\cdot{\bm{V}}_{+}$ in our isolated
system with Galilean invariance.

Transverse magnons form bound states at the wall as follows. In the classical
limit of $-\frac{\hbar^{2}}{2M}\frac{{\rm d}^{2}}{{\rm
d}x^{2}}\to\frac{P_{x}^{2}}{2M}$ in the semiclassical theory of bosonic-
quasiparticle mechanics WKB ; takeuchi2018doubly , we obtain the classical
energy $E_{\rm classic}(P_{x},x)$ as the eigenvalue for ${\bm{k}}=0$ as
$E_{\rm
classic}^{2}=\left(\frac{P_{x}^{2}}{2M}+V_{+}\right)\left(\frac{P_{x}^{2}}{2M}+V_{-}\right)$
with the semiclassical momentum $P_{x}$ in the $x$ direction, and
$V_{\pm}(x)=(c_{n}+c_{s})n\pm 2c_{s}f_{+1}f_{-1}-\tilde{\mu}$. The classical
energy can be lower than the energy gap $|\tilde{q}|$ of the magnon in the
bulk; $E_{\rm classic}(0,x)=V_{\rm eff}(x)\leq|\tilde{q}|$, where we use
$V_{\rm
eff}(x)\equiv\sqrt{\left[(c_{s}+c_{n})n-\tilde{\mu}\right]^{2}-4c_{s}^{2}f_{+1}^{2}f_{-1}^{2}}$
and the inequalities $n<n_{\rm F}$ and $f_{+1}f_{-1}\neq 0$ at the wall. As
the magnons must have an energy higher than $|\tilde{q}|$ to propagate in the
bulk, the bulk is the classically forbidden region for magnons with $E_{\rm
classic}<|\tilde{q}|$, forming bound states at the wall.

The bound magnons can have even nonpositive or imaginary eigenvalues,
indicating instability leading to the spontaneous excitation and condensation
of the magnons. The stability was evaluated using the criterion $V_{\rm
eff}=0$ in the classical limit. The semiclassical theory is applicable for the
spin-1 BECs of 87Rb and 41K with a small negative $c_{s}/c_{n}$, where the
spatial variations of $n$ and $V_{\rm eff}$ are small, with $n\approx n_{\rm
F}$. For $c_{s}\to 0$, the system becomes dynamically unstable when $E_{\rm
classic}(0,x)$ is imaginary, with $\tilde{q}\geq\tilde{q}_{\rm
C}\approx\frac{c_{s}}{c_{n}}\tilde{\mu}\to 0$, where we use $V_{\rm
eff}(0)\approx\tilde{q}^{2}-\left(\frac{c_{s}}{c_{n}+c_{s}}\right)^{2}(\tilde{\mu}-\tilde{q})^{2}$.

To provide a quantitative understanding beyond the semiclassical perspective,
we focus on the excitations at the exact DW solution (5) in 7Li condensates
with $c_{s}=-\frac{c_{n}}{2}$. The perturbation theory of bosonic-
quasiparticle mechanics Skryabin2000PhysRevA.63.013602 ;
Lundh2006PhysRevA.74.063620 ; takeuchi2018doubly ;
2021Takeuchi_PhysRevA.104.013316  is extended to calculate the criterion
$\tilde{q}_{\rm C}(c_{s}=-\frac{c_{n}}{2})$. Regarding $h_{\rm p}$ as the
perturbation in Eq. (6), the unperturbed solution obeys

$\displaystyle E_{\pm}\phi_{\pm}=\pm h_{0}\phi_{\pm}$ (7)

where we use the eigenvalue $\hbar\tilde{\omega}=E_{\pm}$ and eigenvector
$\vec{u}_{0}=\vec{w}_{\pm}$ with $\vec{w}_{+}=(\phi_{+},0)^{\rm T}$ and
$\vec{w}_{-}=(0,\phi_{+})^{\rm T}$ for $h_{\rm p}=0$. This equation is
solvable as a single-particle problem in the Pöschl–Teller potential
LandauQuantum , yielding the bound states with $\pm
E_{\pm}=\frac{\hbar^{2}k^{2}}{2M}-\frac{\tilde{\mu}}{2}-\frac{\tilde{q}}{2}$
and $\phi_{\pm}=\frac{1}{\sqrt{2\xi}}\rm{sech}\left(\frac{z}{\xi}\right)$.

As the low-energy excitation is localized in the core, the perturbed solution
is constructed by a combination of the bound states in the two-mode
approximation as $\vec{u}_{0}=C_{+}\vec{w}_{+}+C_{-}\vec{w}_{-}$ TwoApprox ,
and we obtain

$\displaystyle\hbar\tilde{\omega}=\hbar\tilde{\omega}_{\pm}=\pm\sqrt{\left(\frac{\hbar^{2}k^{2}}{M}-\tilde{\mu}-\tilde{q}\right)\left(\frac{\hbar^{2}k^{2}}{4M}+\frac{\tilde{\mu}-7\tilde{q}}{12}\right)}.$
(8)

The system becomes dynamically unstable when the plus and minus branches of
Eq. (8) collide at $\tilde{q}=\tilde{q}_{\rm
C}(c_{s}=-\frac{c_{n}}{2})=-\tilde{\mu}$, with $k=0$. The dispersion (8) is
consistent with our numerical results near the critical point SM1 .

Critical spin-current velocities— The instability can be interpreted from the
perspective of spin-current instability in the AF phase. The spin current in
the AF state is dynamically unstable when the spin-current velocity is higher
than $V_{\rm AF}=\sqrt{-2\frac{q}{M}}$ Zhu2015PhysRevA.91.023633 , leading to
the nucleation of the $m=0$ component owing to the collisional annihilation of
the $m=\pm 1$ components. A similar interpretation is applicable to the local
AF state at the wall; the AF-core DW becomes unstable when the spin-current
velocity $V_{-}=|{\bm{V}}_{-}|$ at the center exceeds

$\displaystyle V_{\rm mag}=\sqrt{\frac{2\tilde{q}_{\rm C}}{M}+V_{\rm
AF}^{2}}~{}~{}~{}(q\leq\tilde{q}_{\rm C}).$ (9)

This interpretation supports the semiclassical theory because Eq. (9) is
consistently asymptotic to $V_{\rm AF}$ for $\tilde{q}_{\rm C}\to 0$, which is
close to the AF phase with $c_{s}\to 0$.

Formula (9) suggests that the AF-core DWs are always unstable for
$\tilde{q}>\tilde{q}_{\rm C}$. Then, the energy of an AF-core DW is higher
than that of a BA-core DW, at which the local BA state is formed with
$s_{\bot}=\sqrt{s_{x}^{2}+s_{y}^{2}}\neq 0$. This is in contrast to the spin-
current instability in the AF phase KimPhysRevLett.119.185302 , where a
“quench” of spin current effectively causes the transition from the AF phase
($\tilde{q}<0$) to the P phase ($\tilde{q}>0$) with $c_{s}>0$. In our system
with ferromagnetic interaction $c_{s}<0$, the spin-current quench causes a
transition from the local AF state to the local BA state with transverse
magnetization. As the magnetization is energetically preferred for negative
$c_{s}$, the criterion $\tilde{q}_{\rm C}/\tilde{\mu}$ is a monotonically
increasing function of $\frac{c_{s}}{c_{n}}$.

The possibility of BA-core formation creates different situations of spin-
current instability. We describe a complex situation, e.g., if a small
magnetic-field gradient is present with $U_{m}(x)\propto-mx$, where the
critical spin-current velocity for ripplon excitation is $V_{\rm
rip}(\alpha)\propto\alpha^{1/4}$ with the tension $\alpha$ of a DW
Takeuchi2010PhysRevB.81.094517 . In particular, the tension $\alpha_{\rm BA}$
of a BA-core DW is lower than that of an AF-core DW ($\alpha_{\rm AF}$), and
thus, $V_{\rm rip}$ decreases when a BA-core DW is formed after magnon
condensation. The KHI does not occur if $V_{-}<V_{\rm rip}(\alpha_{\rm BA})$
after the BA-core formation, whereas it does for $V_{-}>V_{\rm
rip}(\alpha_{\rm BA})$. On the other hand, when ripplons rather than magnons
grow dominantly in the early stage with $V_{-}>V_{\rm rip}(\alpha_{\rm AF})$,
quantized vortices are nucleated by following the same scenario as the KHI in
binary BECs Takeuchi2010PhysRevB.81.094517 . Then, a BA-core DW will not be
realized if the local spin-current velocity at the DW becomes lower than
$V_{\rm mag}$; otherwise, the magnon condensation occurs following the same
scenario as above.

Skyrmion generation.— Finally, we present an anomalous phenomenon observed in
the numerical experiments of the spin-current instability at a BA-core DW. The
spiraling spin texture along a flat DW [Fig. 2(a)] is identical to the texture
along a magnetic DW in the presence of a spin current in a magnet
Kim2017PhysRevLett.119.047202 , where the spin current is metastable, leading
to the generation of skyrmions with a unit topological charge. However, our
spin-current instability can generate skyrmions with fractional charges.

To demonstrate this anomaly effectively, we consider the simplest case of a
small Weber number with $U_{m}=0$, where the flutter-finger pattern appears in
the early stage of the two-dimensional KHI Kokubo2021_PhysRevA.104.023312
[Fig. 2(b)]. An important distinction becomes apparent after the vortices are
nucleated from the fingertips of DW waves SM1 . The vortex nucleation causes
numerous spin singularities along the wall at which the AF state (${\bm{s}}=0$
and $\Psi_{0}=0$) occurs locally [Fig. 2(c)]. Surprisingly, such singularities
exist for a long time and survive even on the DW loops released into the bulk
[see a DW loop on the far-right side of Fig. 2(d)].

Figure 2:  Dynamics of the spin-current instability at a BA-core DW for
$(c_{s}/c_{n},\tilde{q}/\tilde{\mu},V_{\pm 1}\tau/\xi)=(-0.5,-0.05,\pm 0.368)$
with $\tau=\frac{\hbar}{\tilde{\mu}-\tilde{q}}$. The arrows show the spin
density ${\bm{s}}$. The distributions of $s_{z}$ and $n_{0}=|\Psi_{0}|^{2}$
are represented by the arrow color and background grayscale, respectively. (a)
$t/\tau=0$: The wall is flat in the initial state. (b) $t/\tau=200$: The KHI
makes the flutter-finger pattern. (c) $t/\tau=625$: The vortices are released
as DW loops, and then spin singularities appear along the DWs at which spin
density vanishes locally. (d) $t/\tau=1100$: A spin singularity survives even
on a released loop, forming a fractional skyrmion. The spin texture of
skyrmions is schematically illustrated for (e) $N_{\rm s}=0$, (f) $N_{\rm
s}=1/2$, and (g) $N_{\rm s}=1$, together with the distribution of local
ordered states (BA, AF, F±).

To explain the topological peculiarity of this phenomenon, we consider an
isolated $m=\mp 1$ domain immersed in a sea of $m=\pm 1$ domains and compute
the Mermin–Ho relation Mermin1976_PhysRevLett.36.832.2  between the skyrmion
charge $N_{\rm s}=\frac{1}{4\pi}\int_{\cal
S}dxdy{\bm{s}}\cdot(\partial_{x}{\bm{s}}\times\partial_{y}{\bm{s}})$
Skyrme1961_Proc.R.Soc.A  and the vortex winding number $N_{\rm
v}=\frac{M}{2\pi\hbar}\int_{\cal S}dxdy({\bm{\nabla}}\times{\bm{v}}_{\rm
mass})_{z}=\frac{M}{2\pi\hbar}\oint_{\cal C}d{\bm{r}}\cdot{\bm{v}}_{\rm mass}$
kawaguchi2012spinor .

$\displaystyle N_{\rm s}=\frac{1}{2}N_{\rm v}$ (10)

Here, ${\cal S}$ is the surface enclosed by a closed contour ${\cal C}$ that
surrounds the isolated domain, and ${\bm{v}}_{\rm
mass}=\frac{\hbar}{Mn}\sum_{m}{\rm Im}(\Psi_{m}^{*}{\bm{\nabla}}\Psi_{m})$ is
the mass-current velocity. Although conventional skyrmions have an integer
$N_{\rm s}$ under the boundary condition
${\bm{s}}(|{\bm{r}}|\to\infty)\parallel\pm\hat{\bm{z}}$, we may recognize the
existence of a fractional skyrmion with $N_{\rm s}=1/2$ corresponding to the
unit vortex charge $N_{\rm v}=1$. We typically observed isolated domains with
$N_{\rm v}=0,1,2$ in the numerical simulations, and the corresponding
skyrmions with $N_{\rm s}=0,1/2,1$ are illustrated in Fig. 2(e-g). For
example, the solid angle covered by the spin texture in the $m=+1$ domain at
the lower-left side of Fig. 2(c) is $4\pi$ (i.e., $N_{\rm s}=1$), which can be
understood from the fact that ${\bm{s}}$ winds the equator once along the DW
loop. In contrast, for an isolated domain with $N_{\rm s}=\frac{1}{2}$ on the
far-right side of Fig. 2(d), ${\bm{s}}$ winds half of the equator along the DW
by avoiding a spin singularity.

Discussion.— Our fractional skyrmions with spin singularity are in contrast to
the meron or Mermin–Ho texture Parts1995_PhysRevLett.75.3320 ;
volovik2003universe ; girvin1999quantum ; Moon1995_PhysRevB.51.5138 ;
Kasamatsu2004_PhysRevLett.93.250406 ; kasamatsu2005vortices ;
Ezawa2011_PhysRevB.83.100408 ; Kharkov2017_PhysRevLett.119.207201 ;
gao2019creation . The former exists as an isolated object, whereas fractional
charges in the latter are embedded in a periodic texture or realized only in
pairs under the same boundary condition, e.g., recent observations in
antiferromagnets gao2020fractional ; jani2021antiferromagnetic . Fractional
skyrmions may be assumed unstable because spin singularity increases the
energy cost SV_3He-A . However, the cost is suppressed by the local AF
(nematic-spin) order at the singularity. Furthermore, the splitting of an
integer skyrmion with $N_{\rm s}=1$ into two isolated fractional skyrmions
with $N_{\rm s}=1/2$ is preferred with respect to kinetic energy since two
isolated vortices with $N_{\rm v}=1$ have lower kinetic energy than a vortex
with $N_{\rm v}=2$ in the first approximation pethick2008bose . Energetics are
of fundamental importance because the relationship between the skyrmion charge
and the vortex state with off-centered or eccentric spin singularity in Fig.
2(f) reminds us of that between the vortex winding rule and the non-
axisymmetic vortex TakeuchiPhysRevLett.126.195302 ;
2021Takeuchi_PhysRevA.104.013316 , which is related to the stability of the
eccentric skyrmions. As spin singularities have been numerically observed to
be nucleated via the inhomogeneous growth of the transverse-magnon
condensation (not shown here), our theoretical predictions are expected to be
examined if the ferromagnetic condensates are quenched into the deep F phase
in the experiment Huh2020PhysRevResearch.2.033471 , where domain-coarsening
dynamics are agitated by the KHI, leading to a quantum anomaly in the dynamic
scaling law TakeuchiPhysRevA.97.013617 ; takeuchi2016domain ;
Bourges2017_PhysRevA.95.023616 .

###### Acknowledgements.

This study was supported by JSPS KAKENHI Grant Numbers JP18KK0391, JP20H01842,
and 20H01843 and in part by the OCU “Think globally, act locally ” Research
Grant for Young Scientists through the hometown donation fund of Osaka City.

## References

  * (1) Alex Hubert and Rudolf Schäfer.  Magnetic domains: the analysis of magnetic microstructures.  Springer Science & Business Media, 2008.
  * (2) Marin Gilles, P-Y Bony, Jocelin Garnier, Antonio Picozzi, Massimiliano Guasoni, and Julien Fatome.  Polarization domain walls in optical fibres as topological bits for data transmission.  Nature photonics, 11(2):102–107, 2017.
  * (3) T. Vachaspati.  Kinks and Domain Walls: An Introduction to Classical and Quantum Solitons.  Cambridge University Press,Cambridge, 2006.
  * (4) S. Chandrasekhar.  Liquid Crystals.  Cambridge University Press, 2 edition, 1992.
  * (5) Del Atkinson, Dan A Allwood, Gang Xiong, Michael D Cooke, Colm C Faulkner, and Russell P Cowburn.  Magnetic domain-wall dynamics in a submicrometre ferromagnetic structure.  Nature materials, 2(2):85–87, 2003.
  * (6) Gen Tatara and Hidetoshi Fukuyama.  Resistivity due to a domain wall in ferromagnetic metal.  Physical Review Letters, 78(19):3773, 1997.
  * (7) P. M. Walmsley and A. I. Golov.  Chirality of Superfluid ${}^{3}\mathrm{He}$-$A$.  Phys. Rev. Lett., 109:215301, Nov 2012.
  * (8) H Ikegami, Y Tsutsumi, and K Kono.  Chiral symmetry breaking in superfluid 3He-A.  Science, 341(6141):59–62, 2013.
  * (9) Hiroki Ikegami, Yasumasa Tsutsumi, and Kimitoshi Kono.  Observation of intrinsic Magnus force and direct detection of chirality in superfluid 3He-A.  Journal of the Physical Society of Japan, 84(4):044602, 2015.
  * (10) J. Kasai, Y. Okamoto, K. Nishioka, T. Takagi, and Y. Sasaki.  Chiral Domain Structure in Superfluid ${}^{3}\mathrm{He}\text{$-$}A$ Studied by Magnetic Resonance Imaging.  Phys. Rev. Lett., 120:205301, May 2018.
  * (11) Lev V. Levitin, Ben Yager, Laura Sumner, Brian Cowan, Andrew J. Casey, John Saunders, Nikolay Zhelev, Robert G. Bennett, and Jeevak M. Parpia.  Evidence for a Spatially Modulated Superfluid Phase of ${}^{3}\mathrm{He}$ under Confinement.  Phys. Rev. Lett., 122:085301, Feb 2019.
  * (12) G. E. Volovik and K. Zhang.  String monopoles, string walls, vortex skyrmions, and nexus objects in the polar distorted B phase of 3He.  Phys. Rev. Research, 2:023263, Jun 2020.
  * (13) K. Zhang.  One-dimensional nexus objects, network of Kibble-Lazarides-Shafi string walls, and their spin dynamic response in polar-distorted $B$-phase of ${}^{3}\mathrm{He}$.  Phys. Rev. Research, 2:043356, Dec 2020.
  * (14) JT Mäkinen, VV Dmitriev, Jaakko Nissinen, Juho Rysti, GE Volovik, AN Yudin, Kuang Zhang, and VB Eltsov.  Half-quantum vortices and walls bounded by strings in the polar-distorted phases of topological superfluid 3He.  Nature communications, 10(1):1–8, 2019.
  * (15) Y. Aoki, Y. Wada, M. Saitoh, R. Nomura, Y. Okuda, Y. Nagato, M. Yamamoto, S. Higashitani, and K. Nagai.  Observation of Surface Andreev Bound States of Superfluid ${}^{3}\mathrm{He}$ by Transverse Acoustic Impedance Measurements.  Phys. Rev. Lett., 95:075301, Aug 2005.
  * (16) S. Murakawa, Y. Tamura, Y. Wada, M. Wasai, M. Saitoh, Y. Aoki, R. Nomura, Y. Okuda, Y. Nagato, M. Yamamoto, S. Higashitani, and K. Nagai.  New Anomaly in the Transverse Acoustic Impedance of Superfluid ${}^{3}\mathrm{He}\mathrm{\text{$-$}}\mathbf{B}$ with a Wall Coated by Several Layers of ${}^{4}\mathrm{He}$.  Phys. Rev. Lett., 103:155301, Oct 2009.
  * (17) Satoshi Murakawa, Yuichiro Wada, Yuta Tamura, Masahiro Wasai, Masamichi Saitoh, Yuki Aoki, Ryuji Nomura, Yuichi Okuda, Yasushi Nagato, Mikio Yamamoto, et al.  Surface Majorana cone of the superfluid 3He B phase.  Journal of the Physical Society of Japan, 80(1):013602, 2010.
  * (18) Okuda Y and Nomura R.  Surface Andreev bound states of superfluid3He and Majorana fermions.  Journal of Physics: Condensed Matter, 24(34):343201, aug 2012.
  * (19) P. Zheng, W. G. Jiang, C. S. Barquist, Y. Lee, and H. B. Chan.  Anomalous Damping of a Microelectromechanical Oscillator in Superfluid ${}^{3}\mathrm{He}$-B.  Phys. Rev. Lett., 117:195301, Nov 2016.
  * (20) P. Zheng, W. G. Jiang, C. S. Barquist, Y. Lee, and H. B. Chan.  Critical Velocity in the Presence of Surface Bound States in Superfluid ${}^{3}\mathrm{He}\text{$-$}\mathrm{B}$.  Phys. Rev. Lett., 118:065301, Feb 2017.
  * (21) Masashige Matsumoto and Manfred Sigrist.  Quasiparticle states near the surface and the domain wall in apx$\pm$ipy-wave superconductor.  Journal of the Physical Society of Japan, 68(3):994–1007, 1999\.
  * (22) Phillip E. C. Ashby and Catherine Kallin.  Suppression of spontaneous supercurrents in a chiral $p$-wave superconductor.  Phys. Rev. B, 79:224509, Jun 2009.
  * (23) I. Serban, B. Béri, A. R. Akhmerov, and C. W. J. Beenakker.  Domain wall in a chiral $p$-wave superconductor: A pathway for electrical current.  Phys. Rev. Lett., 104:147001, Apr 2010.
  * (24) J. A. Sauls.  Surface states, edge currents, and the angular momentum of chiral $p$-wave superfluids.  Phys. Rev. B, 84:214509, Dec 2011.
  * (25) Adrien Bouhon and Manfred Sigrist.  Current inversion at the edges of a chiral $p$-wave superconductor.  Phys. Rev. B, 90:220511(R), Dec 2014.
  * (26) Seiji Higashitani.  Spin Current as a Manifestation of Surface Odd-Frequency Pairing in Superfluid 3He.  Journal of the Physical Society of Japan, 83(7):075002, 2014.
  * (27) Y. Tanaka, Y. Asano, A. A. Golubov, and S. Kashiwaya.  Anomalous features of the proximity effect in triplet superconductors.  Phys. Rev. B, 72:140503(R), Oct 2005.
  * (28) Takehito Yokoyama, Yukio Tanaka, and Naoto Nagaosa.  Anomalous Meissner Effect in a Normal-Metal–Superconductor Junction with a Spin-Active Interface.  Phys. Rev. Lett., 106:246601, Jun 2011.
  * (29) S. Higashitani, H. Takeuchi, S. Matsuo, Y. Nagato, and K. Nagai.  Magnetic Response of Odd-Frequency $s$-Wave Cooper Pairs in a Superfluid Proximity System.  Phys. Rev. Lett., 110:175301, Apr 2013.
  * (30) Yasuhiro Asano and Akihiro Sasaki.  Odd-frequency cooper pairs in two-band superconductors and their magnetic response.  Phys. Rev. B, 92:224508, Dec 2015.
  * (31) Seji Kang, Sang Won Seo, Hiromitsu Takeuchi, and Y. Shin.  Observation of Wall-Vortex Composite Defects in a Spinor Bose-Einstein Condensate.  Phys. Rev. Lett., 122:095301, Mar 2019.
  * (32) Hiromitsu Takeuchi.  Quantum Elliptic Vortex in a Nematic-Spin Bose-Einstein Condensate.  Phys. Rev. Lett., 126:195302, May 2021.
  * (33) SeungJung Huh, Kyungtae Kim, Kiryang Kwon, and Jae-yoon Choi.  Observation of a strongly ferromagnetic spinor Bose-Einstein condensate.  Phys. Rev. Research, 2:033471, Sep 2020.
  * (34) C. Hamner, J. J. Chang, P. Engels, and M. A. Hoefer.  Generation of Dark-Bright Soliton Trains in Superfluid-Superfluid Counterflow.  Phys. Rev. Lett., 106:065302, Feb 2011.
  * (35) Joon Hyun Kim, Sang Won Seo, and Y. Shin.  Critical Spin Superflow in a Spinor Bose-Einstein Condensate.  Phys. Rev. Lett., 119:185302, Oct 2017.
  * (36) C. K. Law, C. M. Chan, P. T. Leung, and M.-C. Chu.  Critical velocity in a binary mixture of moving Bose condensates.  Phys. Rev. A, 63:063612, May 2001.
  * (37) Qizhong Zhu, Qing-feng Sun, and Biao Wu.  Superfluidity of a pure spin current in ultracold Bose gases.  Phys. Rev. A, 91:023633, Feb 2015.
  * (38) Yuki Kawaguchi and Masahito Ueda.  Spinor Bose–Einstein Condensates.  Physics Reports, 520(5):253–381, 2012.
  * (39) Justin Lovegrove, Magnus O. Borgh, and Janne Ruostekoski.  Energetically stable singular vortex cores in an atomic spin-1 bose-einstein condensate.  Phys. Rev. A, 86:013613, Jul 2012.
  * (40) Justin Lovegrove, Magnus O. Borgh, and Janne Ruostekoski.  Energetic Stability of Coreless Vortices in Spin-1 Bose-Einstein Condensates with Conserved Magnetization.  Phys. Rev. Lett., 112:075301, Feb 2014.
  * (41) Justin Lovegrove, Magnus O. Borgh, and Janne Ruostekoski.  Stability and internal structure of vortices in spin-1 Bose-Einstein condensates with conserved magnetization.  Phys. Rev. A, 93:033633, Mar 2016.
  * (42) Lauren S Weiss, Magnus O Borgh, Alina Blinova, Tuomas Ollikainen, Mikko Möttönen, Janne Ruostekoski, and David S Hall.  Controlled creation of a singular spinor vortex by circumventing the Dirac belt trick.  Nature communications, 10(1):1–8, 2019.
  * (43) I-Kang Liu, Shih-Chuan Gou, and Hiromitsu Takeuchi.  Phase diagram of solitons in the polar phase of a spin-1 Bose-Einstein condensate.  Phys. Rev. Research, 2:033506, Sep 2020.
  * (44) Andrew P. C. Underwood, D. Baillie, P. Blair Blakie, and H. Takeuchi.  Properties of a nematic spin vortex in an antiferromagnetic spin-1 Bose-Einstein condensate.  Phys. Rev. A, 102:023326, Aug 2020.
  * (45) Stefan Lannig, Christian-Marcel Schmied, Maximilian Prüfer, Philipp Kunkel, Robin Strohmaier, Helmut Strobel, Thomas Gasenzer, Panayotis G. Kevrekidis, and Markus K. Oberthaler.  Collisions of Three-Component Vector Solitons in Bose-Einstein Condensates.  Phys. Rev. Lett., 125:170401, Oct 2020.
  * (46) X. Chai, D. Lao, Kazuya Fujimoto, Ryusuke Hamazaki, Masahito Ueda, and C. Raman.  Magnetic Solitons in a Spin-1 Bose-Einstein Condensate.  Phys. Rev. Lett., 125:030402, Jul 2020.
  * (47) Xiao Chai, Di Lao, Kazuya Fujimoto, and Chandra Raman.  Magnetic soliton: From two to three components with SO(3) symmetry.  Phys. Rev. Research, 3:L012003, Jan 2021.
  * (48) Hiromitsu Takeuchi.  Phase diagram of vortices in the polar phase of spin-1 Bose-Einstein condensates.  Phys. Rev. A, 104:013316, Jul 2021.
  * (49) G C Katsimiga, S I Mistakidis, P Schmelcher, and P G Kevrekidis.  Phase diagram, stability and magnetic properties of nonlinear excitations in spinor Bose–Einstein condensates.  New Journal of Physics, 23(1):013015, jan 2021.
  * (50) See Supplemental Material for details of the current velocity rule, correspondence of the domain-wall solutions, low-energy effective theory, excitation spectrum of the bound transeverse magnons, instability dynamics of a BA-core DW, and numerical methods. .
  * (51) Joseph O. Indekeu, Chang-You Lin, Nguyen Van Thu, Bert Van Schaeybroeck, and Tran Huu Phat.  Static interfacial properties of Bose-Einstein-condensate mixtures.  Phys. Rev. A, 91:033615, Mar 2015.
  * (52) Xiaoquan Yu and P. B. Blakie.  Dark-soliton-like magnetic domain walls in a two-dimensional ferromagnetic superfluid.  Phys. Rev. Research, 3:023043, Apr 2021.
  * (53) Xiaoquan Yu and PB Blakie.  Anomalous dynamics of magnetic field-driven propagating magnetic domain walls.  arXiv preprint arXiv:2104.12967, 2021.
  * (54) Hiromitsu Takeuchi, Naoya Suzuki, Kenichi Kasamatsu, Hiroki Saito, and Makoto Tsubota.  Quantum Kelvin-Helmholtz instability in phase-separated two-component Bose-Einstein condensates.  Phys. Rev. B, 81:094517, Mar 2010.
  * (55) This approximation is a bosonic-quasiparticle version of the Wentzel-Kramers-Brillouin (WKB) approximation in quantum mechanics, see also takeuchi2018doubly . .
  * (56) Hiromitsu Takeuchi, Michikazu Kobayashi, and Kenichi Kasamatsu.  Is a doubly quantized vortex dynamically unstable in uniform superfluids?  Journal of the Physical Society of Japan, 87(2):023601, 2018.
  * (57) Dmitry V. Skryabin.  Instabilities of vortices in a binary mixture of trapped Bose-Einstein condensates: Role of collective excitations with positive and negative energies.  Phys. Rev. A, 63:013602, Dec 2000.
  * (58) Emil Lundh and Halvor M. Nilsen.  Dynamic stability of a doubly quantized vortex in a three-dimensional condensate.  Phys. Rev. A, 74:063620, Dec 2006.
  * (59) Problem 4 in Sec. 23 of L. D. Landau and E. Lifshitz.  Quantum mechanics: Non-relativistic theory, 3rd ed., course of theoretical physics vol. 3.  (Pargamon, New York,1999).
  * (60) The two mode approximation gives an exact solution with $C_{+}=-C_{-}$ for $\tilde{\omega}=k=0$ at the critical point $\tilde{q}=-\tilde{\mu}$. .
  * (61) Se Kwon Kim and Yaroslav Tserkovnyak.  Magnetic Domain Walls as Hosts of Spin Superfluids and Generators of Skyrmions.  Phys. Rev. Lett., 119:047202, Jul 2017.
  * (62) Haruya Kokubo, Kenichi Kasamatsu, and Hiromitsu Takeuchi.  Pattern formation of quantum kelvin-helmholtz instability in binary superfluids.  Phys. Rev. A, 104:023312, Aug 2021.
  * (63) N. D. Mermin and Tin-Lin Ho.  Circulation and Angular Momentum in the $A$ Phase of Superfluid Helium-3.  Phys. Rev. Lett., 36:832–832, Apr 1976.
  * (64) T. H. R. Skyrme.  A non-linear field theory.  Proc. R. Soc. A, 260:127–138, 1961.
  * (65) Ü. Parts, J. M. Karimäki, J. H. Koivuniemi, M. Krusius, V. M. H. Ruutu, E. V. Thuneberg, and G. E. Volovik.  Phase Diagram of Vortices in Superfluid 3 $He-\mathit{A}$.  Phys. Rev. Lett., 75:3320–3323, Oct 1995.
  * (66) Grigory E Volovik.  The universe in a helium droplet, volume 117.  Oxford University Press on Demand, New York, 2003.
  * (67) Steven M Girvin.  The quantum Hall effect: novel excitations and broken symmetries.  In Aspects topologiques de la physique en basse dimension. Topological aspects of low dimensional systems, pages 53–175. Springer, 1999\.
  * (68) K. Moon, H. Mori, Kun Yang, S. M. Girvin, A. H. MacDonald, L. Zheng, D. Yoshioka, and Shou-Cheng Zhang.  Spontaneous interlayer coherence in double-layer quantum Hall systems: Charged vortices and Kosterlitz-Thouless phase transitions.  Phys. Rev. B, 51:5138–5170, Feb 1995.
  * (69) Kenichi Kasamatsu, Makoto Tsubota, and Masahito Ueda.  Vortex Molecules in Coherently Coupled Two-Component Bose-Einstein Condensates.  Phys. Rev. Lett., 93:250406, Dec 2004.
  * (70) Kenichi Kasamatsu, Makoto Tsubota, and Masahito Ueda.  Vortices in multicomponent Bose–Einstein condensates.  International Journal of Modern Physics B, 19(11):1835–1904, 2005\.
  * (71) Motohiko Ezawa.  Compact merons and skyrmions in thin chiral magnetic films.  Phys. Rev. B, 83:100408(R), Mar 2011.
  * (72) Y. A. Kharkov, O. P. Sushkov, and M. Mostovoy.  Bound States of Skyrmions and Merons near the Lifshitz Point.  Phys. Rev. Lett., 119:207201, Nov 2017.
  * (73) Ningbo Gao, S-G Je, M-Y Im, Jun Woo Choi, Masheng Yang, Qin-ci Li, TY Wang, S Lee, H-S Han, K-S Lee, et al.  Creation and annihilation of topological meron pairs in in-plane magnetized films.  Nature communications, 10(1):1–9, 2019.
  * (74) Shang Gao, H Diego Rosales, Flavia A Gómez Albarracín, Vladimir Tsurkan, Guratinder Kaur, Tom Fennell, Paul Steffens, Martin Boehm, Petr Čermák, Astrid Schneidewind, et al.  Fractional antiferromagnetic skyrmion lattice induced by anisotropic couplings.  Nature, 586(7827):37–41, 2020.
  * (75) Hariom Jani, Jheng-Cyuan Lin, Jiahao Chen, Jack Harrison, Francesco Maccherozzi, Jonathon Schad, Saurav Prakash, Chang-Beom Eom, Ariando Ariando, Thirumalai Venkatesan, et al.  Antiferromagnetic half-skyrmions and bimerons at room temperature.  Nature, 590(7844):74–79, 2021.
  * (76) Similar objects have been discussed in superfluid 3He-A; ‘SV’ (singular vortex) in the phase diagram of vortices Parts1995_PhysRevLett.75.3320 ; lounasmaa1999vortices  (see also seppala1983evidence ; Seppala1984_PhysRevLett.52.1802 ; Volovic1984_PhysRevB.29.6090 ; Simola1987_PhysRevLett.58.904 ; fetter1987vortex ) .
  * (77) Christopher J Pethick and Henrik Smith.  Bose–Einstein condensation in dilute gases.  Cambridge university press, 2008.
  * (78) Hiromitsu Takeuchi.  Domain-area distribution anomaly in segregating multicomponent superfluids.  Phys. Rev. A, 97:013617, Jan 2018.
  * (79) Hiromitsu Takeuchi.  Domain size distribution in segregating binary superfluids.  Journal of Low Temperature Physics, 183(3):169–174, 2016.
  * (80) Andréane Bourges and P. B. Blakie.  Different growth rates for spin and superfluid order in a quenched spinor condensate.  Phys. Rev. A, 95:023616, Feb 2017.
  * (81) Olli V Lounasmaa and Erkki Thuneberg.  Vortices in rotating superfluid 3He.  Proceedings of the National Academy of Sciences, 96(14):7760–7767, 1999.
  * (82) Harri Kalervo Seppälä and GE Volovik.  Evidence for nonsingular vorticity in the Helsinki experiments on rotating 3 He-A.  Journal of low temperature physics, 51(3):279–290, 1983.
  * (83) H. K. Seppälä, P. J. Hakonen, M. Krusius, T. Ohmi, M. M. Salomaa, J. T. Simola, and G. E. Volovik.  Continuous Vortices with Broken Symmetry in Rotating Superfluid ${}^{3}\mathrm{He}$-$A$.  Phys. Rev. Lett., 52:1802–1805, May 1984.
  * (84) V. Z. Vulovic, D. L. Stein, and A. L. Fetter.  NMR of textures in rotating ${}^{3}\mathrm{He}-A$.  Phys. Rev. B, 29:6090–6095, Jun 1984.
  * (85) J. T. Simola, L. Skrbek, K. K. Nummila, and J. S. Korhonen.  Two different vortex states in rotating ${}^{3}\mathit{A}$ observed by use of negative ions.  Phys. Rev. Lett., 58:904–907, Mar 1987.
  * (86) Alexander L Fetter.  Vortex structures in rotating 3 He-A.  Journal of low temperature physics, 67(3):145–153, 1987.
  * (87) Tomoya Isoshima, Kazushige Machida, and Tetsuo Ohmi.  Quantum vortex in a spinor Bose-Einstein condensate.  Journal of the Physical Society of Japan, 70(6):1604–1610, 2001\.
  * (88) R. A. Barankov.  Boundary of two mixed Bose-Einstein condensates.  Phys. Rev. A, 66:013612, Jul 2002.
  * (89) Bert Van Schaeybroeck.  Interface tension of Bose-Einstein condensates.  Phys. Rev. A, 78:023624, Aug 2008.
  * (90) Joseph O. Indekeu, Chang-You Lin, Nguyen Van Thu, Bert Van Schaeybroeck, and Tran Huu Phat.  Static interfacial properties of Bose-Einstein-condensate mixtures.  Phys. Rev. A, 91:033615, Mar 2015.
  * (91) Hiromitsu Takeuchi and Kenichi Kasamatsu.  Nambu-Goldstone modes in segregated Bose-Einstein condensates.  Phys. Rev. A, 88:043612, Oct 2013.
  * (92) Shohei Watabe, Yusuke Kato, and Yoji Ohashi.  Excitation transport through a domain wall in a Bose-Einstein condensate.  Phys. Rev. A, 86:023622, Aug 2012.
  * (93) P. K. Kundu and I. M. Cohen.  Fluid Mechanics, 4th ed. Academic Press, New York, 2008.
  * (94) Makoto Tsubota, Michikazu Kobayashi, and Hiromitsu Takeuchi.  Quantum hydrodynamics.  Physics Reports, 522(3):191–238, 2013.
  * (95) Anthony James Leggett et al.  Quantum liquids: Bose condensation and Cooper pairing in condensed-matter systems.  Oxford university press, Oxford, 2006.

##  S1 Current velocity rule

The current velocity rule [Eq. (2)] is a natural extension of the vortex
winding rule for a rotational flow around an axisymmetric vortex
isoshima2001quantum ; TakeuchiPhysRevLett.126.195302 ;
2021Takeuchi_PhysRevA.104.013316 . The rule is applicable when all Zeeman
components have finite population. In the presence of currents along a flat DW
normal to the $x$ axis, the wave function in the stationary state is expressed
as $\Psi_{m}({\bm{r}},t)=f_{m}(x)e^{i\Theta_{m}}$ with the real function
$f_{m}$ and the phase

$\displaystyle\Theta_{m}=\frac{\hbar}{M}{\bm{V}}_{m}\cdot{\bm{r}}+\vartheta_{m},$
(S1)

where we used the current velocity ${\bm{V}}_{m}\bot\hat{\bm{x}}$ of the $m$
component defined by the current density $\mathbf{j}_{m}=\frac{\hbar}{M}{\rm
Im}(\Psi_{m}^{*}{\bm{\nabla}}\Psi_{m})=f_{m}^{2}{\bm{V}}_{m}$. By substituting
this formula into the equation of motion obtained from the Lagrangian, we have
the coupled equations

$\displaystyle\frac{\hbar^{2}}{2M}f_{0}^{\prime\prime}=g_{0}f_{0}+c_{s}f_{0}\left(f^{2}_{+1}+f^{2}_{-1}+2f_{+1}f_{-1}e^{i\delta\Theta}\right)$
(S2) $\displaystyle\frac{\hbar^{2}}{2M}f_{\pm 1}^{\prime\prime}=g_{\pm
1}f_{\pm 1}+c_{s}f^{2}_{0}\left(f_{\pm 1}+f_{\mp 1}e^{-i\delta\Theta}\right)$
(S3)

with $\delta\Theta=\Theta_{+1}+\Theta_{-1}-2\Theta_{0}$,
$f_{m}^{\prime\prime}=\frac{{\rm d}^{2}}{{\rm d}x^{2}}f_{m}$, and

$g_{m}=c_{n}n+mc_{s}s_{z}-\mu-mp+m^{2}q+\frac{1}{2}M{\bm{V}}_{m}^{2}.$

When $f_{m}(x)$ is nonzero for all components by satisfying
$f_{+1}f_{0}f_{-1}\neq 0$ at a certain place, $\delta\Theta$ must be an
integer multiple of $\pi$ for satisfying the coupled equations above. We may
set as $\delta\Theta=0$ when the factor $e^{i\pi}=-1$ is included into the
real function $f_{m}$ by changing its sign. In this way, we have
${\bm{V}}_{+1}+{\bm{V}}_{-1}=2{\bm{V}}_{0}$ with
$\vartheta_{+1}+\vartheta_{-1}-2\vartheta_{0}=0$, equivalent to the current
velocity rule. The velocities ${\bm{V}}_{+}$ and ${\bm{V}}_{-}$ are related to
the mass and spin currents at the soliton center ($x=0$).

By substituting Eq. (2) into Eqs. (S2) and (S3) with $\delta\Theta=0$, one
obtains Eq. (3). Equation (3) are equivalent to Eqs. (S2) and (S3) with
$\delta\Theta=0$ and ${\bm{V}}_{m}=0$ when $(\tilde{\mu},\tilde{p},\tilde{q})$
are replaced by $(\mu,p,q)$; the solution $f_{m}(x)$ of the former is equal to
that of the latter when $(\tilde{\mu},\tilde{p},\tilde{q})=(\mu,p,q)$.
Accordingly, we readjust the condition of the hydrostatic variables as
$\tilde{\mu}>0$ and $\tilde{q}<0$ for realizing the F states in the bulk.

##  S2 Domain-wall solution

Figure S1:  (a) Numerical plots of density profiles $f_{+1}^{2}$ (red mark)
and $f_{-1}^{2}$ (blue mark) of an AF-core DW for
$c_{s}/c_{n}=-0.8,-0.5,-0.2,-0.03$. The dimensionless plots are independent of
$\tilde{q}/\tilde{\mu}$. Solid curves shows the exact solution $(f_{\pm
1}^{\rm ex})^{2}$ for $c_{s}/c_{n}=-0.5$ given by Eq. (5). (b) Numerical plots
of $f_{m}^{2}~{}(m=0,\pm 1)$ of the lowest-energy DW for
$\tilde{q}/\tilde{\mu}=-1,-0.5,-0.25,-0.1$ with $c_{s}/c_{n}=-0.5$.

Here, we demonstrate the details of analytical and numerical results of DW
solutions and the tension $\alpha$. The results are summarized in Figs. S1 and
S2.

We first introduce the useful correspondence between binary BECs and spin-1
BECs in the absence of the $m=0$ component, which enable us to obtain some
analytical expressions. After the replacement of
$(\Psi_{+1},\Psi_{-1})\to(\Psi_{1},\Psi_{2})$ and $(\mu+p-q,\mu-
p-q)\to(\mu_{1},\mu_{2})$, the Lagrangian of spin-1 BECs with $\Psi_{0}=0$ is
rewritten as ${\cal L}=\sum_{j=1,2}\int
d^{3}x\left[i\hbar\Psi_{j}^{*}\partial_{t}\Psi_{j}-\frac{\hbar^{2}}{2M}({\bm{\nabla}}\Psi_{j}^{*})\cdot({\bm{\nabla}}\Psi_{j})-{\cal
U}_{j}\right]$ with

$\displaystyle{\cal
U}_{j}=-\mu_{j}|\Psi_{j}|^{2}+\sum_{k=1,2}\frac{g_{jk}}{2}|\Psi_{j}|^{2}|\Psi_{k}|^{2}.$
(S4)

This is just the Lagrangian of binary BECs with the intra- and inter-component
coupling constants, $g=g_{11}=g_{22}=c_{n}+c_{s}$ and
$g_{12}=g_{21}=c_{n}-c_{s}$, respectively. The ferromagnetic interaction
$(0>c_{s}>-c_{n})$ satisfies the immiscible condition of binary BECs,
$0<g<g_{12}$, while the anti-ferromagnetic one $(c_{s}>0)$ does the miscible
condition, $0<|g_{12}|<g$ pethick2008bose .

Joseph et al. pointed out that there is an exact solution of a DW in a
segregated binary BEC with $g_{12}/g=3$. This case corresponds to
$c_{n}/c_{s}=-1/2$ in spin-1 BECs according to the above correspondence. Then
the exact solution of AF-core DW is given by Eq. (5). Figure S1 shows typical
profiles of AF-core and BA-core DWs obtained numerically by solving Eq. (4).
See the section of “Numerical methods” for details of the numerical method. A
stationary solution of a flat DW is obtained under the Neumann boundary
condition at the $x=\pm L_{x}/2$ with the system size $L_{x}=256\xi$. The bulk
density is determined by the pressure balance between the two domains with
opposite magnetization. Since there exists only a single spin component far
from the DW, the hydrostatic pressure $P_{\pm}$ in a domain with $s_{z}=\pm
n_{\pm}$ is computed in a similar manner as a scalar BEC pethick2008bose ;
$P_{\pm}=-\frac{\partial G}{\partial V}=\frac{1}{2}(c_{n}+c_{s})n_{\pm}^{2}$
with the system volume $V$ and the bulk density
$n_{\pm}=\frac{\tilde{\mu}-\tilde{q}\pm\tilde{p}}{c_{n}+c_{s}}$. The pressure
balance $P_{+}=P_{-}=P_{\rm F}$ is realized under the condition $\tilde{p}=0$
or

$\displaystyle p=M{\bm{V}}_{+}\cdot{\bm{V}}_{\rm-}.$ (S5)

Here, the bulk pressure is represented as

$\displaystyle P_{\rm F}=\frac{1}{2}(c_{n}+c_{s})n_{\rm F}^{2}$ (S6)

with the bulk density $n_{\rm F}$ of Eq. (4). It is noted that a flat DW can
be realized even for nonzero $p$ in the presence of currents by satisfying Eq.
(S5). The Galilean invariance allow us to set as ${\bm{V}}_{+}=0$ without loss
of generality in an isolated system and thus we put $p=0$ in this work.

Figure S2:  The tension $\alpha$ (a) and the maximum transverse spin density
$\max{s_{\bot}}$ (b) of the lowest-energy DW for
$c_{s}/c_{n}=-0.03,-0.2,-0.5,-0.8$. The analytical results, $\alpha_{\rm ex}$
[Eq. (S9)] and $\alpha_{\rm weak}$ [Eq. (S10)], for AF-core DWs are compared
with some numerical results. The tension is asymptotic to the critical
behavior $\frac{\alpha}{(\tilde{\mu}-\tilde{q})n_{\rm
F}\xi}\propto\sqrt{|\tilde{q}|/\tilde{\mu}}$ for $\tilde{q}/\tilde{\mu}\to 0$.
The lowest-energy DW is a BA-core DW with $\max{s_{\bot}}>0$ when $\tilde{q}$
is larger than a critical value $\tilde{q}_{\rm C}(c_{s})<0$.

Next, we derive the formula of DW tension, which is important to determine its
static and dynamic properties. The tension $\alpha$ is defined by the excess
energy in the presence of a DW of area $S$ as

$\displaystyle\alpha$ $\displaystyle=$ $\displaystyle S^{-1}\int
d^{3}x\left[{\cal G}-\tilde{\cal U}_{\rm b}\right]$ (S7) $\displaystyle=$
$\displaystyle\int_{-\infty}^{\infty}dx\left[\frac{\hbar^{2}}{2M}\sum_{m}f_{m}^{\prime
2}+U_{\rm pot}\right].$

where $\tilde{\cal U}_{\rm b}=-P_{\rm F}$ is the energy density in the bulk.
The second line is obtained by substituting $\Psi_{m}=f_{m}e^{i\Theta_{m}}$
into the first line and composed of the kinetic energy associated with the
spatial gradient $f_{m}^{\prime}=\frac{{\rm d}}{{\rm d}x}f_{m}$ and the
residual term $U_{\rm pot}$. Equation (3) is represented as
$f_{m}^{\prime}\frac{\hbar^{2}}{M}\frac{{\rm d}^{2}}{{\rm
d}x^{2}}f_{m}^{\prime}=f_{m}^{\prime}\frac{\partial U_{\rm pot}}{\partial
f_{m}}$. The integration and summation of this equation with respect to $x$
and $m$, one obtains $\frac{\hbar^{2}}{2M}\sum_{m}f_{m}^{\prime 2}=U_{\rm
pot}$ with the boundary conditions $f_{m}^{\prime}(\pm\infty)=0$ and $U_{\rm
pot}(f_{m}(\pm\infty))=0$. By using this result, the tension is represented by
the formula

$\displaystyle\alpha=\frac{\hbar^{2}}{M}\sum_{m}\int_{-\infty}^{\infty}dxf_{m}^{\prime
2}=2\int_{-\infty}^{\infty}dxU_{\rm pot}.$ (S8)

The tension $\alpha_{\rm AF}$ of an AF-core DW is evaluated analytically for
some cases. The tension $\alpha_{\rm ex}$ for the exact solution of Eq. (5) is
computed as

$\displaystyle\alpha_{\rm AF}=\alpha_{\rm
ex}=\frac{2}{3}(\tilde{\mu}-\tilde{q})n_{\rm F}\xi~{}~{}~{}(c_{s}=-c_{n}/2).$
(S9)

According to the correspondence with binary BECs, the tension is asymptotic to
$\alpha_{\rm weak}$ ($\alpha_{\rm stong}$) in the weakly (strongly)
segregating limit of $\frac{g_{12}}{g}\to 1$ ($\frac{g_{12}}{g}\to\infty$),
given by Barankov2002_PhysRevA.66.013612 ; Schaeybroeck2008_PhysRevA.78.023624
; Indekeu2015_PhysRevA.91.033615

$\displaystyle\alpha_{\rm AF}\to\alpha_{\rm
weak}=\sqrt{\frac{-c_{s}}{c_{n}+c_{s}}}(\tilde{\mu}-\tilde{q})n_{\rm
F}\xi~{}~{}~{}(c_{s}\to 0),$ (S10) $\displaystyle\alpha_{\rm AF}\to\alpha_{\rm
strong}=\frac{4}{3}(\tilde{\mu}-\tilde{q})n_{\rm F}\xi~{}~{}~{}(c_{s}\to-
c_{n}).$ (S11)

The tension $\alpha_{\rm BA}$ of a BA-core DW is computed by numerically
minimizing $\alpha$, or equivalently, solving Eq. (3). Figure S2 shows the
$\tilde{q}$-dependence of $\alpha_{\rm BA}$ for some typical values of $c_{s}$
together with the plot of the peak value $\max s_{\bot}$ of the transverse
spin density at the DW. The tension $\alpha_{\rm BA}$ approaches $\alpha_{\rm
AF}$ for $\tilde{q}\to\tilde{q}_{\rm C}$, at which a BA-core DW becomes an AF-
core DW with $\max s_{x}=0$.

##  S3 Low-energy effective theory

We introduce the low-energy effective theory to describe the hydrodynamic
aspect of a DW, whose theoretical framework is not affected explicitly by the
internal structure of DWs. A similar theory has been applied to a DW in binary
BECs for describing quantum Kelvin-Helmholtz instability (KHI)
Takeuchi2010PhysRevB.81.094517  and Nambu-Goldstone (NG) modes
TakeuchiPhysRevA.88.043612 . Here, we extend the theory to spin-1 BECs.

The low-energy effective theory is constructed by considering the degrees of
freedom associated with the NG modes and neglecting higher-energy Bogoliubov
excitations. We first consider Bogoliubov excitations that propagate in the
bulk far from the DW. An excitation in the bulk is described as a collective
fluctuation $\delta\vec{\Psi}=\vec{\Psi}-\vec{\Phi}_{\rm
F}=\vec{u}-\vec{v}^{*}$ with the Bogoliubov coefficients $\vec{u}$ and
$\vec{v}$. The excitations can be classified by the coefficients, and we have
three types of excitations; e.g. in a magnetic domain with $s_{z}=+n$, (i)
$\vec{u}\propto[1,0,0]^{\rm T}$ and $\vec{v}\propto[1,0,0]^{\rm T}$, (ii)
$\vec{u}\propto[0,1,0]^{\rm T}$ and $\vec{v}=0$, and (iii)
$\vec{u}\propto[0,0,1]^{\rm T}$ and $\vec{v}=0$. The mode (i) is a gapless
mode whose excitation spectrum $\varepsilon(k)$ does not have an energy gap
for zero wave number ($k=0$): $\varepsilon(0)=0$. This mode is called phonon,
corresponding to the NG mode associated with the $U(1)$ symmetry breaking.
Other two modes have energy gaps in the excitation spectrum: (ii)
$\varepsilon(0)=|p|-q>0$ and (iii) $\varepsilon(0)=2|p|-2c_{s}n>0$. The
Bogoliubov coefficients of (ii) and (iii) are composed only by the $m=0$ or
$m=-1$ component while those of (i) by the $m=1$ component. Therefore we
neglect the fluctuations of $m=0$ and $m=-1$ components in the low-energy
effective theory. Similarly, in a domain with $s_{z}=-n$, the degrees of
freedom associated with the $m=0$ and $m=+1$ components are neglected
approximately.

There exists another NG mode, called the ripplon, corresponding to the
vibration of a DW. According to the above consideration, we describe the state
in a domain with $s_{z}=\pm n$ by a single wave function $\Psi_{\pm 1}$ in the
effective theory. Representing the position of the DW between two magnetic
domains by the single-valued function $x=\eta(y,z,t)$, we introduce the
effective Lagrangian,

$\displaystyle{\cal L}_{\rm eff}=\int{\rm d}y{\rm
d}z\left(\int^{+\infty}_{\eta}{\cal P}_{+1}dx+\int^{\eta}_{-\infty}{\cal
P}_{-1}dx\right)-\alpha S$

with the area $S$ and the tension $\alpha$ of the DW. The Lagrangian density
${\cal P}_{\pm 1}$ in the domain with $s_{z}=\pm n$ is obtained by neglecting
the contributions from the $m=0$ and $m=\mp 1$ components in the original
Lagrangian ${\cal L}$; ${\cal P}_{\pm 1}=i\hbar\Psi_{\pm
1}^{*}\partial_{t}\Psi_{\pm
1}-\frac{\hbar^{2}}{2M}\sum_{m}({\bm{\nabla}}\Psi_{\pm
1}^{*})\cdot({\bm{\nabla}}\Psi_{\pm 1})-{\cal U}_{\pm 1}$ with ${\cal U}_{\pm
1}=\frac{1}{2}(c_{n}+c_{s})|\Psi_{\pm 1}|^{4}-(\mu\pm p-q)|\Psi_{\pm 1}|^{2}.$

By writing as $S=\int{\rm d}x{\rm
d}y\sqrt{1+(\partial_{y}\eta)^{2}+(\partial_{z}\eta)^{2}}\approx\int{\rm
d}x{\rm
d}y\left[1+\frac{1}{2}(\partial_{y}\eta)^{2}+\frac{1}{2}(\partial_{z}\eta)^{2}\right]$
for a small fluctuation from a stationary state of $\eta=0$, one obtains the
equation of motion with respect to $\eta$ from the effective Lagrangian,

$\displaystyle{\cal P}_{+1}(\eta)-{\cal
P}_{-1}(\eta)=\alpha\left(\partial_{y}^{2}+\partial_{z}^{2}\right)\eta.$ (S12)

This equation is an analog of the Young-Laplace equation in fluid mechanics,
where the Laplace pressures ${\cal P}_{+1}-{\cal P}_{-1}$ balances the effect
of the shape of the wall represented by the right hand side of Eq. (S12). In
the stationary state of $\eta=0$ the pressures in the two domains are equal
with ${\cal P}_{+1}={\cal P}_{-1}=P_{\rm F}$. The phonon is decoupled with the
ripplon in the first approximation. This is because the low-energy phonon can
propagate through the wall without disturbance, associated with the anomalous
tunneling effect of a DW PhysRevA.86.023622 . Hence, the translational motion
of a wall is treated independently from the bulk fluctuation and we mainly
focus the former as a vibration mode of a DW.

The vibration of a DW is determined by Eq. (S12) under a proper boundary
condition on the wall. The boundary condition can be different between our
system and binary BECs while the theoretical treatment of the bulk wave
function is identical to each other. A possible mechanism to make the
difference is the spin interaction. The number of particles is conserved for
each component in binary BECs, which demands the kinematic boundary condition
2008Kundu , whereas the population transfer occurs between different spin
components via the spin interaction in spinor BECs. However, we neglect the
transfer in the effective theory for the following reason. The mechanism of
the transfer originally comes from the operator $\propto
c_{s}\hat{\Psi}_{0}^{\dagger}\hat{\Psi}_{0}^{\dagger}\hat{\Psi}_{+1}\hat{\Psi}_{-1}$
and its Hermitian conjugate in the spin interaction in the many-body
Hamiltonian. The particle transfer can happen together with a pair
annihilation or creation of the $m=0$ component, resulting in a fluctuation of
the wave function $\Psi_{0}$. Such a process is forbidden at the DW in the
low-energy effective theory because it can cost an extra energy by changing
the internal structure of the wall, except for the case of the ripplon
excitation to induce the transition to a BA-core DW. Moreover, the ripplon is
the NG mode associated with the spontaneous breaking of translational
symmetry, and the fluctuation due to a ripplon excitation makes a
translational shift, expressed by $\delta\vec{\Psi}\propto\eta\frac{{\rm
d}}{{\rm d}x}\vec{\Psi}$ by conserving the population of each component. We
hence apply the same boundary condition as the one used in segregated binary
BECs.

Once the correspondence between binary BECs and our system is established on
the level of the effective theory, we can use the result of the dispersion
relation for ripplon excitations in binary BECs. Here, the roles of the first
and second components in Ref. Takeuchi2010PhysRevB.81.094517 ;
TakeuchiPhysRevA.88.043612 ; tsubota2013quantum  are played by the $m=-1$ and
$m=+1$ components in our system, respectively. The dispersion $\omega_{\rm
rip}({\bm{k}})$ of a ripplon mode
$\eta\propto\cos({\bm{k}}\cdot{\bm{r}}-\omega_{\rm rip}t)$ is given by

$\displaystyle\omega_{\rm
rip}({\bm{k}},{\bm{V}}_{+},{\bm{V}}_{-})={\bm{k}}\cdot{\bm{V}}_{+}+\omega_{\rm
KH}({\bm{k}},{\bm{V}}_{-})$ (S13)

with the KH dispersion

$\displaystyle\omega_{\rm
KH}({\bm{k}},{\bm{V}}_{-})=\sqrt{\frac{\alpha|{\bm{k}}|^{3}}{2Mn_{\rm
F}}-\left({\bm{k}}\cdot{\bm{V}}_{-}\right)^{2}}.$ (S14)

The first term on the right hand side of Eq. (S13) comes from the so-called
Doppler effect induced by the center-of-mass motion related to the Landau
instability of ripplon Takeuchi2010PhysRevB.81.094517 .

Equation (S14) is regarded as the dispersion relation of capillary waves in
fluid dynamics. We can realize also an analogue of the gravity-capillary waves
by introducing external potentials, which plays the role of gravity to
stabilize the interface at the equilibrium position $\eta=0$
Takeuchi2010PhysRevB.81.094517 . Such potentials are realized by imposing a
magnetic field gradient in the $x$ direction, which is represented by the
potential felt by the $m$ components as $U_{m}(x)=-m\frac{F}{2n_{\rm F}}x$.
Considering a two-dimensional system for ${\bm{V}}_{+}=V_{+}\hat{\bm{y}}$,
${\bm{V}}_{-}=V_{-}\hat{\bm{y}}$ and ${\bm{k}}=k\hat{\bm{y}}$ under the
potential, the dispersion relation of Eq. (S13) is replaced by

$\displaystyle\omega_{\rm rip}=kV_{+}+|k|\sqrt{\frac{F+\alpha
k^{2}}{2|k|Mn_{\rm F}}-V_{-}^{2}}.$ (S15)

The dynamic instability happens when the frequency becomes complex for
$\frac{F+\alpha k^{2}}{2|k|Mn_{\rm F}}-V_{\rm-}^{2}<0$. The instability is
realized when the spin-current velocity $V_{-}$ at the DW exceeds the critical
value

$\displaystyle V_{\rm rip}=\sqrt{\frac{\sqrt{F\alpha}}{Mn_{\rm F}}}.$ (S16)

This result implies that ripplon excitations cause the dynamic instability for
an arbitrary spin-current velocity without magnetic-field gradient ($F=0$).

##  S4 Excitation spectrum of bound transverse magnons

Figure S3:  Bogoliubov excitation spectrum of the 1st and 2nd lowest-energy
excitations at an AF-core DW for $c_{s}/c_{n}=-1/2$ and $V_{+}=0$. Circles
show the real part (black circle) and the imaginary part (red circle) obtained
numerically by solving the BdG equation. Solid curves represents the
analytical plots of ${\rm Re}~{}\tilde{\omega}_{+}$ (green) and ${\rm
Im}~{}\tilde{\omega}_{+}$ (orange) in Eq. (8) and ${\rm Re}~{}\omega_{\rm KH}$
(blue) in Eq. (S14) with $\alpha=\alpha_{\rm ex}$.

Here, we derive Eq. (8) in the two-mode approximation based on the
perturbation theory of bosonic-quasiparticle mechanics, which is a natural
extension of the perturbation theory in quantum mechanics for single-particle
problem. In this analogy, the eigenvalue equation of Eq. (6) corresponds to
the ‘Schrödinger equation’ in quantum mechanics.

We consider the perturbation expansion in the following

$\displaystyle\hbar\tilde{\omega}\vec{u}_{0}=(\hat{h}_{0}+\delta\hat{h})\vec{u}_{0}$
(S17)

with $\displaystyle\hat{h}_{0}=\begin{pmatrix}h_{0}&0\\\
0&-h_{0}\end{pmatrix}$ and $\displaystyle\delta\hat{h}=\begin{pmatrix}h_{\rm
p}\ &h_{\rm p}\\\ -h_{\rm p}&-h_{\rm p}\end{pmatrix}$.

The unperturbed solution is given by the equations,

$\displaystyle\begin{cases}\varepsilon u_{0}=h_{0}u_{0}\\\ \varepsilon
v_{0}=-h_{0}v_{0}\end{cases}$ (S18)

with the unperturbed eigenvalue $\varepsilon$. These equations are solved
exactly for the case of $c_{s}=-\frac{c_{n}}{2}$. Then, substituting the exact
DW solution of Eq. (5), Eqs. (S18) is reduced to

$\displaystyle\begin{cases}\epsilon_{+}u_{0}=\left[-\frac{1}{2}\frac{d^{2}}{d(x/\xi)^{2}}-{\rm
sech}^{2}\left(\frac{x}{\xi}\right)\right]u_{0}\\\
\epsilon_{-}v_{0}=\left[-\frac{1}{2}\frac{d^{2}}{d(x/\xi)^{2}}-{\rm
sech}^{2}\left(\frac{x}{\xi}\right)\right]v_{0}\end{cases}$ (S19)

with $\epsilon_{\pm}=\frac{\tilde{\mu}-c_{n}n_{\rm
F}/2-e_{k}\pm\varepsilon}{c_{n}n_{\rm F}/2}$ and
$e_{k}=\frac{\hbar^{2}k^{2}}{2M}$.

The eigenvalue solution is given by solving the single-particle problem in the
Pöschl-Teller potential LandauQuantum ; $\epsilon_{\pm}=-\frac{1}{2}$ and
$u_{0},v_{0}\propto{\rm sech}\left(\frac{x}{\xi}\right)$. Accordingly, we
obtain the original eigenvalues

$\displaystyle\varepsilon=\hbar\omega_{\pm}=\pm(e_{k}-\tilde{\mu}+c_{n}n_{\rm
F}/4).$ (S20)

It is expected that the instability is induced by these bound states and the
perturbed solution could be constructed by a combination of them;

$\displaystyle\vec{u}_{0}=C_{+}\vec{u}_{+}+C_{-}\vec{u}_{-},$ (S21)

where we used the eigen solutions
$(\varepsilon,\vec{u}_{0})=(\hbar\omega_{\pm},\vec{u}_{\pm})$;

$\displaystyle\begin{cases}\vec{u}_{+}=\left[\frac{1}{\sqrt{2\xi}}{\rm
sech}\left(\frac{x}{\xi}\right),0\right]^{\rm T}\\\
\vec{u}_{-}=\left[0,\frac{1}{\sqrt{2\xi}}{\rm
sech}\left(\frac{x}{\xi}\right)\right]^{\rm T}.\end{cases}$ (S22)

According to the formulation of the two-mode approximation, one obtains the
result, equivalent to Eq. (8),

$\displaystyle\hbar\tilde{\omega}$ $\displaystyle=$
$\displaystyle\frac{\tilde{\varepsilon}_{+}+\tilde{\varepsilon}_{-}}{2}\pm\sqrt{\left(\frac{\tilde{\varepsilon}_{+}-\tilde{\varepsilon}_{-}}{2}\right)^{2}+M_{-+}M_{+-}}$
(S23) $\displaystyle=$
$\displaystyle\pm\sqrt{\frac{1}{12}\left(2e_{k}-\tilde{\mu}-\tilde{q}\right)\left(6e_{k}+\tilde{\mu}-7\tilde{q}\right)}$

with $\tilde{\varepsilon}_{\pm}=\hbar\omega_{\pm}+M_{\pm\pm}$. Here,
$M_{\alpha\beta}$ is defined as

$\displaystyle M_{\alpha\beta}=N_{\alpha\alpha}\int
d^{3}x\vec{u}_{\alpha}^{\dagger}\hat{\sigma}_{z}\delta\hat{h}\vec{u}_{\beta}$
(S24)

with $N_{\alpha\beta}=\int
d^{3}x\vec{u}_{\alpha}^{\dagger}\hat{\sigma}_{z}\vec{u}_{\beta}=\pm\delta_{\alpha\beta}$
and $\hat{\sigma}_{z}={\rm diag.}(1,-1)$. The matrix elements
$M_{\alpha\beta}$ were computed by using the normalization condition
$N_{\pm\pm}=\pm 1$ as

$\displaystyle M_{++}=-M_{--}=M_{+-}=-M_{-+}=\frac{c_{n}n_{\rm F}}{6}.$ (S25)

The dispersion relation of Eq. (8) is compared with the results of the
numerical diagonalization of the Bogoliubov-de Gennes (BdG) equation
leggett2006quantum , which is obtained after linearizing the equation of
motion with respect to $u_{m}(x)$ and $v_{m}(x)$. See the section of
“Numerical methods” for details on the method of the numerical
diagonalization.

##  S5 Instability dynamics of a BA-core DW

Here, we explain the detailed information on the numerical result of the spin-
current instability demonstrated in Fig. 2. The images in Fig. 2(a-d) are
magnified snapshots of a numerical simulation. Figure S4 shows the time
evolution of the phases $\arg\Psi_{\pm 1}$ and the transverse spin amplitude
$s_{\bot}$ in the same simulation.

The time evolution of the simulation is obtained by numerically solving the
time-dependent Gross-Pitaevskii (GP) equations of spin-1 BECs
kawaguchi2012spinor , the equation of motion obtained from the Lagrangian. The
system size of the simulation is $256\xi\times 256\xi$ and we imposed the
Neumann and periodic boundary conditions in the $x$ and $y$ directions,
respectively. The instability is triggered by adding a small random seed to
the stationary solution of a BA-core DW as the initial state ($t=0$) of the
time evolution. See the section of “Numerical methods” for details of the
numerical method.

The attached movie files demonstrate the animations of different quantities
made with 51 snapshots from $t/\tau=0$ to $t/\tau=1250$ in the same
simulation. The box size is the same as the system size of the simulation.
Figure 2(a-d) are magnified images of some snapshots in the animation of Movie
S1. Top ($\arg\Psi_{-1}$), middle ($\arg\Psi_{+1}$) and bottom ($s_{\bot}$)
panels in Fig S4 are some snapshots in the animations of Movie S2, Movie S3,
and Movie S4, respectively. The plot area of Fig. 2(a-d) is implied by dashed
squares in the bottom panels in Fig. S4.

Some isolated domains with different topological charges are implied by arrows
in the bottom right in Fig. S4. According to Eq. (10), we can convert the
number $N_{\rm v}$ of branch cuts (jump from $\arg\Psi_{\pm 1}=-\pi$ to
$\pi$), which are terminated by an isolated domain, to the skymion charge
$N_{\rm s}(=2N_{\rm v})$ of the domain. We see that a branch cut passes
through an isolated domain corresponding to the case of $N_{\rm s}=0$ from the
comparison between the distributions of $\arg\Psi_{+1}$ and $s_{\bot}$ at
$t/\tau=1100$. Two branch cuts end at a domain by forming a conventional
skymion with $N_{\rm s}=1$. An isolated domain that terminates a branch cut
corresponds to an eccentric fractional skymion with $N_{\rm s}=1/2$. The spin
singularities appear in the form of cut points of a BA-core domain or white
spots in the distribution of $s_{\bot}$ in the bottom panels of Fig. S4.

Figure S4:  Snapshots of the phases $\arg\Psi_{\pm 1}$ and the transverse spin
density $s_{\bot}$ in the time evolution of Fig. 2. The display range of Fig.
2 is shown by red broken lines in the bottom panels. Black areas in top and
middle panels represent the region of $s_{z}>0$ and $s_{z}<0$, respectively.
The transverse spin density $s_{\bot}$ has a similar distribution to the
background images of $n_{0}$ in Fig. 2. In the later stage of the dynamics, a
lot of ”notches” appear as white spots along DWs in the distribution of
$s_{\bot}$, which correspond to the spin singularities in the spin density.
The number of branch cuts (jump from $\arg\Psi_{\pm 1}=-\pi$ to $\pi$)
terminated by an isolated domain wall is equal to $N_{\rm v}(=N_{\rm s}/2)$. A
spin singularity exists along a closed DW surrounding an isolated domain, at
which a branch cut is terminated. The closed DW with a spin singularity looks
like a “C” mark in the distribution plot of $s_{\bot}$ or $n_{0}$. No
singularity happen if the domain contains even number of terminated branch
cuts. Some examples of domains with $N_{\rm s}=0,1/2,1$ are denoted by red
arrows in the lower right.

##  S6 Numerical methods

Here, we explain the method of numerical simulations used in this work. All
simulations are done by rescaling the length, time, and wave function by
$\xi$, $\tilde{\mu}-\tilde{q}$, and $\sqrt{n_{\rm F}}$, respectively.

First, we describe how to obtain the DW solutions in Fig. S1. The stationary
solution of a magnetic DW is obtained by solving Eq. (3) with the steepest
descent method under the Neumann boundary condition $\left.\frac{{\rm
d}f_{m}}{{\rm d}x}\right|_{x=\pm L_{x}/2}=0$ at the system boundary
$x=\pm\frac{L_{x}}{2}$. The space coordinate is discretized as $x\to
x_{i}=-\frac{L_{x}}{2}+\Delta(i-1/2)~{}(i=0,1,2,...,N_{x}+1)$ with
$N_{x}=1024$ and $\Delta=0.25\xi$. The spatial derivative of $f_{m}$ is
computed with finite difference approximation; $\frac{{\rm d}^{2}f_{m}}{{\rm
d}x^{2}}$ is computed by the central difference of the second order. The
solutions of AF-core DWs is obtained by imposing the condition $f_{0}=0$.

The excitation spectrum in Fig. S3 is obtained by numerically diagonalizing
the full BdG equations with respect to the eigenvector
$[u_{+1}(x),u_{0}(x),u_{-1}(x),v_{+1}(x),v_{0}(x),v_{-1}(x)]^{\rm T}$ around
the DW solution obtained above. The spatial discretization, the spatial
derivative, and the boundary condition for the eigenvector field are done in a
similar way as described above. The numerical diagonalization is performed
with the double precise by using Intel® Fortran Compiler with the Linear
Algebra PACKage (LAPACK).

The result of Fig. 2 (Fig. S4 and the movie files) is obtained by solving the
GP equation in two dimensions. The space coordinate is discretized as
$(x,y)\to(x_{i},y_{j})=(-\frac{L}{2}+\Delta(i-1/2),-\frac{L}{2}+\Delta(j-1/2))~{}(i,j=0,1,2,...,N+1)$
with $N=512$ and $\Delta=0.5\xi$. Here, we impose the periodic boundary
condition in the $y$ direction as
$\Psi_{m}(x_{i},y_{N+1})=\Psi_{m}(x_{i},y_{1})$ and
$\Psi_{m}(x_{i},y_{0})=\Psi_{m}(x_{i},y_{N})$. The spatial derivative is done
in a similar way as described above. The initial state of the time evolution
is prepared by the steepest descent method. The time $t$ is discretized as
$t=\Delta tn_{t}$ with $\Delta t=0.0025\tau$. A small random fluctuation added
to the initial state to seed the instability. The time evolution is computed
by utilizing the Crank-Nicolson method.