Decay of the relative phase domain wall into confined vortex pairs: the case of a coherently coupled bosonic mixture

TL;DR

This paper predicts how a relative phase domain wall in a coherently coupled Bose-Einstein condensate decays through vortex pair formation or snake instability, depending on the Rabi coupling strength, supported by numerical simulations.

Contribution

It introduces a detailed analysis of decay mechanisms of phase domain walls in coupled BECs, highlighting the role of Rabi coupling strength and providing numerical predictions for experimental setups.

Findings

Small $\Omega$ causes vortex-antivortex pair formation and wall fragmentation.

Large $\Omega$ leads to snake instability and vortex pair confinement.

Numerical simulations match experimental sodium gas configurations.

Abstract

A domain wall of relative phase in a flattened harmonically-trapped Bose-Einstein condensed mixture of two atomic hyperfine states, subject to a stationary Rabi coupling of intensity $\Omega$, is predicted to decay through two different mechanisms. For small values of $\Omega$ the instability has an energetic nature and is associated with the formation of a vortex-antivortex pair of the same atomic hyperfine states, whose motion inside the trap causes the emergence of magnetization, the bending of the domain wall and its consequent fragmentation. For large values of $\Omega$ the domain wall instead undergoes a dynamic snake instability, caused by the negative value of its effective mass and results in the fast fragmentation of the wall into smaller domain walls confining vortex pairs of different atomic species. Numerical predictions are given by solving the time-dependent Gross-Pitaevskii equation in experimentally available configurations of mixtures of sodium atomic gases.

Equations

$${\mathcal{E}}_{\mathrm{GP}}=$$

$$+g_{12}|{\Psi }_1{|}^2|{\Psi }_2{|}^2-\frac{1}{2}\hslash \Omega \left({\Psi }_1^{\ast }{\Psi }_2+{\Psi }_2^{\ast }{\Psi }_1\right),$$

$${\varphi }_{1,2}(x)=\pm 2\arctan (\exp (x/{\xi }_{\Omega })),$$

$$m^{\ast }=m\frac{\sqrt{3}{\pi }^2}{4}{\xi }_{s}n\frac{\sqrt{\Omega {\Omega }_c}}{{\Omega }_c-\Omega },$$

$$\underset{q\to 0}{\lim }\frac{\omega (q)}{q}=\sqrt{\frac{E_{\mathrm{DW}}^{(v=0)}}{m^{\ast }}},$$

$$\eta =\frac{n_1-n_2}{n}=-\frac{m^{\ast }v}{\hslash n}\text{mbox}{sech}^2\left(\frac{\Omega }{3{\Omega }_c{\xi }_s}(x-vt)\right),$$

$$U\frac{\partial \eta }{\partial \zeta }-2\sin \eta \frac{\partial \eta }{\partial \zeta }\frac{\partial {\varphi }_R}{\partial \zeta }+\cos \eta \frac{{\partial }^2{\varphi }_R}{\partial {\zeta }^2}={\kappa }^2\sin {\varphi }_R,$$

$$U\cos \eta \frac{\partial {\varphi }_R}{\partial \zeta }=$$

$$-{\kappa }^2\sin \eta \cos {\varphi }_R,$$

$$\varphi ={\varphi }_1-{\varphi }_2=4\arctan (\exp (\kappa \zeta )).$$

$$-\frac{{\partial }^2\eta }{\partial {\zeta }^2}-\eta \left[{\left(\frac{\partial \varphi }{\partial \zeta }\right)}^2-1-{\kappa }^2\cos \varphi \right]=-U\frac{\partial \varphi }{\partial \zeta },$$

$${\kappa }^2\cos \varphi ={\kappa }^2+{\left({\partial }_{\zeta }\varphi \right)}^2-\frac{6{\kappa }^2}{{\cosh }^2(\kappa \zeta )}$$

$$\hat{\mathcal{O}}\eta =-U\frac{\partial \varphi }{\partial \zeta },$$

$$\hat{\mathcal{O}}=-\frac{{\partial }^2}{\partial {\zeta }^2}-\frac{6{\kappa }^2}{{\cosh }^2(\kappa \zeta )}+1+{\kappa }^2$$

$${\psi }_{\mathcal{O}}(\zeta )=\frac{\sqrt{3\kappa }}{2{\cosh }^2(\kappa \zeta )},$$

$$\eta (\zeta )=\mathcal{C}{\psi }_{\mathcal{O}}(\zeta )+\delta \eta ,$$

$$\mathcal{C}=-\frac{U}{1-3{\kappa }^2}{\int }_{-\infty }^{\infty }d\zeta {\psi }_{\mathcal{O}}^{\ast }(\zeta )\frac{\partial \varphi }{\partial \zeta }=-\frac{\sqrt{3}\pi }{2}\frac{\sqrt{\kappa }U}{1-3{\kappa }^2}.$$

$$E$$

$$\displaystyle+\sin^{2}\eta+2\kappa^{2}(1-\cos\eta\cos\phi_{R})\bigg{]}\;,$$

$$E^{(2)}=\frac{n\hslash c_s}{4}\int d\zeta \left[{\left(\frac{\partial \eta }{\partial \zeta }\right)}^2-{\eta }^2{\left(\frac{\partial \varphi }{\partial \zeta }\right)}^2+{\eta }^2+{\kappa }^2{\eta }^2\cos \varphi \right].$$

$$E^{(2)}=$$

$$\frac{3{\pi }^2}{16}\frac{n\hslash \kappa v^2}{c_s}\frac{1}{1-3{\kappa }^2},$$

$$m^{\ast }=\underset{v\to 0}{\lim }\frac{1}{v}\frac{dE(v)}{dv},$$

$$m^{\ast }=m\frac{3{\pi }^2}{4}\frac{{\xi }_{s}n\kappa }{1-3{\kappa }^2}.$$

$$\omega (q)=i\gamma (q)$$

$$\mathcal{L}=$$

$$-\frac{1}{2}\left[(\nabla {\varphi }_R{)}^2+(\nabla {\varphi }_T{)}^2+(\nabla \eta {)}^2-2\sin \eta \nabla {\varphi }_R\nabla {\varphi }_T\right]$$

$$+{\kappa }^2\cos \eta \cos {\varphi }_R.$$

$$\frac{\partial \mathcal{L}}{\partial \mathcal{F}}-\nabla \frac{\partial \mathcal{L}}{\partial (\nabla \mathcal{F})}-\frac{d}{d\tau }\frac{\partial \mathcal{L}}{\partial ({\partial }_{\tau }\mathcal{F})}=0,$$

$$\frac{\partial \eta }{\partial \tau }+\frac{{\partial }^2\chi }{\partial {\zeta }^2}-q^2{\xi }_s^2\chi -\nabla (\eta \nabla {\varphi }_T)={\kappa }^2\cos {\varphi }_R$$

$$\frac{\partial \chi }{\partial \tau }+\frac{{\partial }^2\eta }{\partial {\zeta }^2}-q^2{\xi }_s^2\eta +\nabla \chi \nabla {\varphi }_T+\nabla {\varphi }_R\nabla \beta -\eta ={\kappa }^2\eta \cos {\varphi }_R$$

$$\frac{{\partial }^2\beta }{\partial {\zeta }^2}-q^2{\xi }_s^2\beta -\nabla (\eta \nabla {\varphi }_R)=0,$$

$$\eta (\vec{r},\tau )=$$

$${\varphi }_R(\vec{r},\tau )=$$

$${\varphi }_T(\vec{r},\tau )=$$

Full text

The study of topological defects in ordered phases encompasses different fields of physics from quantum fluids and superconductors to cosmology. Of particular interest is the case of composite defects due to the presence of multiple symmetry breaking Bunkov and Godfrin (2000). In this respect ultra-cold gases mixtures are emerging as one of the most suitable platforms for studying the dynamical behaviour of topological defects and in particular their formation and decay. The easiest system showing interesting composite defects is most probably a binary Bose-Bose mixture with an interconversion term, realized via external coherent Rabi coupling, between two hyperfine components Matthews et al. (1999a); Zibold et al. (2010); Nicklas et al. (2011, 2015). The intriguing features of such a gas are related to the controlled – via the tunability of the Rabi coupling – breaking of the conservation of the relative atom number. Only the total atom number is conserved, due to the $U(1)$ symmetry related global phase (see, e.g., Abad and Recati (2013) and reference therein).

An interesting solitonic solution is represented by the sine-Gordon domain wall of the relative phase, explored by Son and Stephanov Son and Stephanov (2002) for Bose-Einstein condensed gases and by Tanaka Tanaka (2001) for two-band superconductors. The relative phase domain wall (hereafter simply called domain wall) is characterized by an asymptotic $2\pi$ jump of the relative phase between the two condensates, whose space gradient is localized in a narrow region fixed by the strength of the Rabi coupling, which makes it energetically costly. Such a solitonic configuration is particularly relevant because it provides non trivial constraints on the topology of quantized vortices. Indeed, in the presence of coherent coupling between the two condensates, a quantized vortex in one of the two components of the mixture, also known as half-quantized vortex, has always a domain wall attached to it. As a consequence isolated half-quantized vortices cannot exist, but are rather paired in order to screen the relative phase change. This situation is reminiscent of the problem of quarks in particle physics, which cannot exist as individual particles, but are confined by strings. The analogy with the physics of quarks and strings is further reinforced by the linear increase of the confining potential with the distance between the vortices in a pair and the possibility for the domain wall to connect either a vortex and an antivortex in the same atomic species or vortices with the same angular momentum in different atomic species (such pairs are named mesons and baryons, respectively, in the recent work by Eto and Nitta Eto and Nitta (2018)).

Because of the high energy cost, it is preferable for long domain walls to decay into smaller ones connecting vortex pairs Son and Stephanov (2002); Gallemí et al. (2016); Tylutki et al. (2016); Eto and Nitta (2018). This is expected to be the relevant decay mechanism for small values of the Rabi coupling. If instead the Rabi coupling is large, domain walls become dynamically unstable Son and Stephanov (2002) and exhibit snake instability. Understanding the mechanism of fragmentation of the domain wall is of crucial importance in view of future experiments aimed to investigate the formation and the dynamics of confined vortex pairs in ultracold atomic gases.

In this Letter we first provide a thorough discussion of the equilibrium and dynamical properties of the domain wall in uniform matter. Then, by carrying out a fully-numerical analysis based on the solution of the time-dependent Gross-Pitaevskii equations (GPEs), we explore the decay mechanisms of a domain wall in a mixture trapped by a harmonic potential, suggesting a natural protocol for possible experiments. As we show in this Letter, the decay of a domain wall into shorter domain walls can be monitored not only by measuring the relative phase of the gas, but also through the total density and the magnetization profiles.

We consider a zero temperature mixture of two hyperfine states of a bosonic atomic isotope, which are coherently-coupled through a Rabi term $\Omega$. Defining $(\Psi_{1},\Psi_{2})$ the order parameter of the mixture, the grand canonical energy functional takes the form:

[TABLE]

where $m$ is the atomic mass, $\mu$ is the chemical potential and $V_{\rm ext}$ an external trapping potential. The interaction strengths $g_{ij}=4\pi\hbar^{2}a_{ij}/m$ and $g_{12}=4\pi\hbar^{2}a_{12}/m$ are given in terms of the $s$-wave scattering lengths, $a_{ij}$. In the following we consider the specific case of the two hyperfine levels $|F=1,m_{F}=\pm 1\rangle$ of 23Na for which $a_{11}=a_{22}=54.54a_{B}$ and $a_{12}=50.78a_{B}$ (where $a_{B}$ is the Bohr radius) and we define $g=g_{11}=g_{22}$. The last equality means that the system is $\mathbb{Z}_{2}$ symmetric and the ground state is characterised by $n_{1}=n_{2}=n/2$ 111The $\mathbb{Z}_{2}$ symmetry is spontaneously broken if $\delta g$ is enough negative (see, e.g., Abad and Recati (2013)). A peculiarity of sodium atoms is that the two $m_{F}=\pm 1$ states are miscible ($\delta g=g-g_{12}>0$) and that the value $\delta g/g=0.07$ is sufficiently large to allow for the observability of the effects discussed in this Letter. These important features satisfied by 23Na have already proven to be important to observe magnetic effects Gallemí et al. (2018); Kang et al. (2019) and spin superfluidity Kim et al. (2017); Fava et al. (2018).

As already mentioned, the two coupled GPEs (see Supplementary Material, Eq.(D1)) derived by minimizing the energy functional (1) admit, in uniform matter, a stationary solitonic solution characterized by a peculiar behavior of the relative phase between the two condensates. Indeed under the assumption that the densities are homogeneous and equal, the GPEs can be reduced to the simpler sine-Gordon equation $\hbar\nabla^{2}\phi=m\Omega\sin\phi$ Son and Stephanov (2002) for the relative phase $\phi=\phi_{1}-\phi_{2}$ between the two condensates. This equation admits the following 1D solitonic solution Tabor (1989)

[TABLE]

which connects the values [math] and $2\pi$ of the relative phase $\phi$ as $x$ goes from $-\infty$ to $\infty$. Such a solution generates a counterflow current, accumulated in a small region of the relative phase gradient of size $\xi_{\Omega}=\sqrt{\hbar/m\Omega}$ (Rabi healing length), fixed by $\Omega$. The solution (2) of the sine-Gordon equation is actually an approximation to the solution of the coupled GPEs, which still exhibits an equivalent phase pattern, but with a density dip near the wall, as a consequence of the compressibility of the gas. The dip increases with the increase of the Rabi coupling Usui and Takeuchi (2015), as shown in Fig. 1 (a), and in the case of 23Na, the central density vanishes for $\Omega=20.57\,\Omega_{c}$, reaching the dark soliton state Tsuzuki (1971).

A major issue concerns the stability of the domain wall. Dark solitons Muryshev et al. (1999); Pitaevskii and Stringari (2016), as well as dark-bright solitons Busch and Anglin (2001), are dynamical unstable, except for one dimensional geometries. As first pointed out by Son and Stephanov Son and Stephanov (2002), the relative phase domain wall is instead dynamically stable for $\Omega$ smaller then a critical value $\Omega_{c}$, even in 3D configurations. Under the assumption that $\hbar\Omega,\delta gn\ll gn$, the critical value takes the simple expression $\Omega_{c}=n\delta g/3$ Son and Stephanov (2002). The (transverse) dynamical stability of the domain wall is related to the sign of its effective mass, which fixes the energy change of the configuration when the soliton moves with small velocity $v$ according to $E_{DW}(v)=E_{DW}^{(v=0)}+m^{*}v^{2}/2$, where $E_{DW}^{(v=0)}=A^{-1}\int d\vec{r}\,^{3}(\mathcal{E}_{\rm GP}[\Psi_{\rm DW}]-\mathcal{E}_{\rm GP}[\Psi_{\rm GS}])$ is the energy cost per unit area ($A$) to create a domain wall at rest. Here $\Psi_{\rm GS}$ is the order parameter of a uniform mixture and $\Psi_{\rm DW}$ the one in presence of the domain wall. Analytical expression for $m^{*}$ can be obtained in two limiting cases. For $\Omega\ll\Omega_{c}$, one gets Qu et al. (2017) $m^{*}=8m\xi_{s}n\sqrt{\Omega/\Omega_{c}}$, where $\xi_{s}=\hbar/\sqrt{2m\delta gn}$ is the spin healing length. In the case $\Omega\simeq\Omega_{c}$ we find that the effective mass per unit area (see Supplementary material, Section II) in the incompressible regime ($\hbar\Omega\ll\mu$) can be approximated as

[TABLE]

showing that the effective mass diverges at $\Omega_{c}$ and becomes negative for $\Omega>\Omega_{c}$.

It was shown in Kamchatnov and Pitaevskii (2008) for dark solitons that the low-momentum dispersion for the transverse excitations can be written as

[TABLE]

where $q$ is the in-plane wave vector of the excitation. Equation (4) holds under very general assumptions and can be applied also to the domain wall. In particular one sees that for $\Omega<\Omega_{c}$, the positiveness of the effective mass results into a real excitation frequency. For larger values of $\Omega$, the effective mass is instead negative, leading to a purely imaginary excitation frequency, which indicates the occurrence of a dynamic instability. In Fig. 1 (c) we report the behavior of the imaginary component of the transverse excitations for $\Omega>\Omega_{c}$ as a function of the dimensionless parameter $q\xi$, with $\xi=\hbar/\sqrt{mgn}$, obtained by solving the Bogoliubov equations (see Supplementary material, section IV) in the case of 23Na. The spectrum is reminiscent of the dark soliton one Muryshev et al. (1999), and converges to it when the central density becomes zero. In Fig. 1(d) we also compare the value of the effective mass calculated numerically using Eq. (4) with the prediction of the analytic formula (3). The agreement is very good in a wide range of values of the Rabi coupling.

As we will discuss in the second part of the Letter domain walls are most likely generated at a finite velocity. In Ref. Qu et al. (2017) it was shown that a moving domain wall is charaterized by finite magnetization localized near the wall. Both the density dip, discussed above, and the magnetization are important features for the experimental identification of the soliton. The direct measurement of the relative phase is in fact of difficult experimental access Matthews et al. (1999b), although the relative flow, i.e., the gradient of $\phi$, could be extracted via spin selective spatially resolved Bragg scattering. This technique has been very recently applied in Ref. Kang et al. (2019) to identify topological composite objects in a spin-1 Bose gas.

Again in the limit $\hbar\Omega,\delta gn\ll gn$ and for small velocity, an analytical expression for the local magnetization $\eta$ can be derived (see Supplementary material, Section I):

[TABLE]

which is shown in Fig. 1(b) in units of $\eta/U$ where $U\equiv v/c_{s}$ is the velocity in units of the spin speed of sound $c_{s}=\hbar/2m\xi_{s}$ calculated for $\Omega=0$.

Having characterized the basic properties of the domain wall in uniform matter we are now ready to discuss the possible realization and the study of the different decay mechanisms when the mixture is trapped in a harmonic potential. In order to provide a quantitative description of the domain wall decay we rely on the solution of time-dependent coupled GPEs. We consider $N=8\times 10^{5}$ atoms in a quasi-2D geometry confined by a cylindrically symmetric harmonic trap $V_{\rm ext}=m(\omega^{2}r^{2}+\omega_{z}^{2}z^{2})/2$, with $\omega=2\pi\times 8\,\mbox{Hz}$ and $\omega_{z}=2\pi\times 320\,\mbox{Hz}$ 222In a quasi-2D geometry, the scattering lengths are renormalized to $a_{i}/\sqrt{2\pi l_{z}^{2}}$ with $l_{z}=\sqrt{\hbar/m\omega_{z}}$. By evolving the GPE in imaginary time, we first prepare the condensate in its ground state in which each atom is in a 50/50 superposition of the two hyperfine internal levels. The cloud has the typical inverted parabola density profile and we define for later convenience the central critical Rabi coupling as $\Omega_{c}^{0}=\delta gn(0)/3$ with $n(0)$ the density in the center of the trap.

A relative phase profile in the form of Eq. (2), which crosses the whole cloud at the center of the trap, is then imprinted onto the condensate wave functions and let to evolve in real time by using the GPEs. Notice that this initial condition should be easily realised in current experiments. Since the phase imprinting does not change the total density, the initial configuration is not the proper domain wall solution and a local dip is rapidly produced (see first and fourth column of Fig. 2). The noise created in this first stage of the evolution will play a crucial role in driving the subsequent dynamics of the wall, which depends very much on whether the value of $\Omega$ is smaller or larger than $\Omega_{c}^{0}$. The two cases give rise to very different scenarios for the evolution of the domain wall, both concerning the nature of the instability and the corresponding time scales.

When $\Omega<\Omega_{c}^{0}$, a local density approach suggests that in the region where $\Omega<\Omega_{c}(n(r))$ the domain wall is dynamically stable. For sufficiently small values of $\Omega$ the region of dynamical instability is confined to a narrow region near the surface (which is shown as a thin line) and we find that a domain wall confining a vortex-antivortex pair in the same component is formed after a while. The vortices are located symmetrically at a distance $R$ from the center of the trap, near the surface. This pair does not carry angular momentum, and starts moving with linear velocity, causing the bending of the wall and the emergence of a visible local magnetization (second column of Fig. 2). Such a bending is in agreement with the bending already reported in Ref. Eto and Nitta (2018) in uniform matter. When the bending becomes pronounced the domain wall breaks down into a certain number of vortex-vortex and antivortex-antivortex pairs and, depending on the actual value of $\Omega$, an additional vortex-antivortex pair, that rapidly moves towards the cloud boundary, where gets absorbed or breaks into other pairs. A typical frame of such rich dynamics is shown in the third column of Fig. 2.

When $\Omega>\Omega_{c}^{0}$ a vortex-vortex pair appears, which starts rotating, dragging the edges of the wall (see Fig. 2, fifth column). However, the domain wall is dynamically unstable similarly to the homogeneous case and the noise caused by the phase imprinting procedure triggers the decay, leading to the snaking of the domain wall. The simulation shows that on short times this snaking nucleates multiple domain walls of smaller size which confine pairs of vortices of equal vorticity, belonging to different atomic species (sixth column of Fig. 2). These smaller objects are more stable since the minimum length required to exhibit dynamic instability is given by $\pi/q_{0}$, where $q_{0}$ is the threshold value of the wave vector below which the Bogoliubov excitation frequency becomes imaginary (see Fig. 1(c) and Supplementary Material, Sections III and IV).

As reported in Fig. 2, it is evident that the domain wall breaking for $\Omega<\Omega_{c}^{0}$ is almost an order of magnitude slower than for $\Omega>\Omega_{c}^{0}$, reflecting the energetic nature of the instability. With the chosen values of the trap fequencies the decay time varies from 1 to 10 seconds in the two cases reported in Fig. 2. Even assuming that such time scale is reduced by one or two orders of magnitude in a real experiment, such a phenomenon should be still measurable.

For completeness in the first two rows of Fig. 2, we also show the behavior of the total density and of the magnetization for the same temporal sequence. These quantities are measurable experimentally (typical values of the dip reported in Fig. 2 correspond to a contrast of about 5% in the central region) and can therefore directly probe the nature of the decay and of the corresponding dynamical evolution of the vortex pairs, given the strong correlation between the density dip, the magnetized wall arising at $\Omega<\Omega_{c}^{0}$ and the kink in the relative phase.

The initial conditions discussed so far, based on the phase imprinting of the domain wall is not the only possibility that one can explore theoretically. In order to check the general validity of the instability mechanisms, we have considered two additional routes for generating the decay. In the first case we initially create an infinitely long domain wall corresponding to the exact stationary solution of the GPEs. Being exact, such a metastable configuration would not evolve in time unless a small amount of noise is added to the system. We have analyzed this possibility by evolving in time the wave function $\beta(\vec{r}\,)\Psi_{\rm DW}(\vec{r}\,)$, where $\Psi_{\rm DW}(\vec{r}\,)$ is the exact domain wall solution and $\beta(\vec{r})$ is a random number in the range $(1-\varepsilon,1+\varepsilon)$, $\varepsilon$ being the noise strength ($1\gg\varepsilon>0$). The presence of additional noise triggers the dynamical instability for $\Omega>\Omega_{c}^{0}$ and it is indispensable for $\Omega<\Omega_{c}^{0}$ in order to favour the transition to lower energy configurations characterized by a domain wall of smaller size connecting a vortex pair, which is protected by an energy barrier due to the presence of the trap (see Supplementary Material, Section V). Alternatively instead of adding noise to the initial configuration, we produce a fast ramping of $\Omega$, which sets the system out of equilibrium. Such a protocol would allow to improve the visibility of the magnetized strings and the density dip. In both cases, very similar dynamics take place.

In conclusion, we have investigated the decay mechanisms exhibited by a relative phase domain wall in a coherently-coupled condensate, with special emphasis to the role played by the external harmonic trapping. When the Rabi frequency is larger than the critical value above which domain walls are dynamically unstable, the snake instability mechanism triggers the decay, leading to the breaking of the domain wall into confined vortex pairs. In contrast, when the Rabi frequency is below the critical one the instability of the soliton has an energetic nature and additional excitations are required to allow the entrance of the vortices from the border of the trapped atomic cloud. The novel scenarios discussed in this Letter are expected to open challenging perspectives to be explored in future experiments on coherently coupled atomic gases to obtain composite topological objects Kang et al. (2019) and mimic pair creation and string breaking.

Note added:

During the final writing of this Letter we have become aware of a new theoretical work which addresses the problem of the dynamic instability of the relative phase domain walls in uniform matter Ihara and Kasamatsu (2019), in general agreement with some of our findings.

Acknowledgements.

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No. 641122 “QUIC”, the FIS$\hbar$ project of the Istituto Nazionale di Fisica Nucleare and the Provincia Autonoma di Trento. We would like to thank useful discussions with Gabriele Ferrari and Giacomo Lamporesi.

In the following Supplementary Material we derive some useful results for the magnetization, the effective mass and the frequency of the elementary excitations associated with the domain wall near the critical point $\Omega_{c}$, which have been discussed in the main text. In Section E we also derive explicit results for the barrier encountered by a vortex pair before entering a harmonically trapped atomic cloud.

Appendix A Velocity-induced magnetization

In this section we calculate the magnetization of a slowly moving domain wall. Let us define $\kappa^{2}=\Omega/3\Omega_{c}=\xi_{s}/\xi_{\Omega}$, where $\xi_{s}=\hbar/\sqrt{2m\delta gn}$ and $\xi_{\Omega}=\sqrt{\hbar/m\Omega}$ are the spin and the Rabi healing lengths respectively (see main text) and introduce $U=v/c_{s}$ the velocity of the domain wall in units of the spin sound speed $c_{s}=\sqrt{\delta gn/2m}$, and $n\sin\eta=n_{1}-n_{2}$ the magnetization of the atomic cloud. In the following analysis we always ignore the density dip exhibited by the domain wall and make us of the Eqs. (15) and (16) of Ref. Qu et al. (2017):

[TABLE]

and

[TABLE]

where $\zeta=(x-vt)/\xi_{s}$ and $\phi_{R}$ is the relative phase of the domain wall. We will solve the above equations considering small deviations from the Son-Stephanov stationary solution

[TABLE]

The perturbative calculation requires $\phi_{R}-\phi=\chi\ll 1$ and $\eta\ll 1$. It is easy to show that to lowest order, $\eta\sim U$ and $\chi\sim U^{2}$. Thus, one can simplify Eq. (7) to the form

[TABLE]

where $\partial_{\zeta}\phi=2\kappa/\cosh(\kappa\zeta)$. Using the identity

[TABLE]

we can rewrite Eq. (9) as:

[TABLE]

where the Pöschl-Teller-like operator

[TABLE]

has the lowest eigenvalue $1-3\kappa^{2}$ with its normalized eigenvector

[TABLE]

corresponding to the bound solution of even parity. The exact solution of Eq. (11) can be presented in the form

[TABLE]

where $\delta\eta$ is the contribution of the delocalized states, orthogonal to $\psi_{\mathcal{O}}$ . By multiplying both members of Eq. (11) by $\psi_{\mathcal{O}}^{*}=$ $\psi_{\mathcal{O}}$, and integrating over $\zeta$ from $-\infty$ to $\infty$, the proportionality constant $\mathcal{C}$ is found to be

[TABLE]

This expression diverges for $\kappa=1/\sqrt{3}$ (the critical point), which suggests that the contribution to the magnetization (14) that comes from the delocalized states can be indeed neglected. Notice that in order to fulfill the condition $\eta\ll 1$ the velocity of the soliton should satisfy the condition ${U\ll 1-3\kappa^{2}}$.

Appendix B Effective mass of a domain wall

In order to calculate the effective mass we start from the general expression of the energy of the domain wall per unit area, given in Ref. Qu et al. (2017):

[TABLE]

where we have again neglected the presence of the density dip. Keeping only terms of order $U^{2}$ (that is $\eta^{2}$ or $\chi$, the terms linear in $\chi$ being vanishing), we can obtain the second order term in the expansion of the energy as a function of the velocity, hereafter called $E^{(2)}$:

[TABLE]

By inserting in Eq. (17) the Son-Stephanov expression (8) for the relative phase $\phi$ of the moving domain wall and performing partial integration, one can rewrite (17) as

[TABLE]

where we have used the ansatz (14) with $\delta\eta=0$ and the approximate result (15) for the constant $\mathcal{C}$.

Thus, the effective mass per unit area, defined by:

[TABLE]

takes the simple insightful form:

[TABLE]

This is the line shown as a red dashed line in Fig. 1(d) in the main text. The limit of validity of this equation is exactly the same as in the calculation of the magnetization (see previous section), being based on the approximated solution (13), yielding result (15) for the relevant constant $\mathcal{C}$.

Appendix C Imaginary frequencies of the excitation spectrum

In order to get a good approximation for the full $q$-dependence of the imaginary frequency:

[TABLE]

that arises at $\Omega>\Omega_{c}$, one has to consider not only the fluctuations of the magnetization and the relative phase, but also the fluctuations of the total phase $\phi_{T}=\phi_{1}+\phi_{2}$, which are particularly relevant at large wave vectors.

To this aim, we will start from the Lagrangian previously introduced in Ref. Qu et al. (2017), including not only the magnetization and the relative phase but also the total phase:

[TABLE]

One can then obtain the corresponding equations of motion for each of the variables $\mathcal{F}=\{\eta,\phi_{R},\phi_{T}\}$, using the Lagrange equation:

[TABLE]

from which one can obtain the following set of the equations:

[TABLE]

where we have expressed the fluctuations as:

[TABLE]

By working to order $U$ (or equivalently $\eta$), the previous equations take the simplified form:

[TABLE]

Let us now solve the first equation, by writing $\chi(\zeta)=\mathcal{D}\psi_{\mathcal{P}}(\zeta)+\delta\chi$, where $\mathcal{D}$ is a proportionality coefficient (the final result will not depend on its value) and

[TABLE]

is the only bound eigenvector of the Pöschl-Teller operator

[TABLE]

whose eigenvalue is [math]. Notice that the zero-mode solution (29) has a simple physical meaning. This is a change of the Son-Stephanov static solution (8) after an infinitesimal small displacement along $\zeta$ direction.

The function $\delta\chi$ is a combination of the other eigenstates that are, by construction, orthogonal to $\psi_{\mathcal{P}}$. Analogously to what we did in Section A for the magnetization, we can multiply both members of Eq. (26) by $\psi_{\mathcal{P}}^{*}$ and integrate over all space to get:

[TABLE]

where

[TABLE]

The other two equations are coupled. In order to decouple them, we first transform Eq. (27) with help of Eq. (10). One gets:

[TABLE]

Now let us multiply this equation on $\psi_{\mathcal{O}}^{*}$ and integrate over the space. The result is

[TABLE]

where we have ignored smaller $\delta\chi$ and $\delta\eta$ terms. This equation still depends on $\beta$, but it can be further simplified if we take into account that because of the fast decreasing of the functions $\psi_{\mathcal{O}}$ and $\partial_{\zeta}\phi$ as $|\zeta|\to\infty$, the integrals of the right hand side of Eq. (34), are relevant only when $\zeta\lesssim\kappa^{-1}$. On the other hand, it is easy to show that the term proportional to $q^{2}$ in Eq. (28) is important only on distances larger than $q^{-1}\xi_{s}^{-1}$. However, near the critical point the instability takes place only for $q<q_{0}\propto\sqrt{3\kappa^{2}-1}\ll 1$ (see below). This means that in the integration (34) one can calculate $\beta$ neglecting the term proportional to $q^{2}$ in Eq. (28), which leads to:

[TABLE]

Integration of Eq. (35) yields

[TABLE]

and substitution into Eq. (34) yields to the cancellation of the two terms in the right hand side. Finally, we obtain the simple equation:

[TABLE]

In conclusion, Eq. (37) can be combined with Eq. (31) to give the final result:

[TABLE]

for the imaginary part of the frequency holding for $\Omega>\Omega_{c}$.

Figure 3 (a) shows the imaginary part of the eigenfrequencies for $\Omega$ larger but close to $\Omega_{c}$, as a function of the wave vector $q$, for different values of the Rabi frequency. The solid lines with symbols correspond to the results obtained from the numerical solution of the Bogoliubov equations (see section D), and the dashed lines with the same colors correspond to the curves given by Eq. (38). Panel (b) shows the value of $q_{0}$, which is the threshold point at which the imaginary frequencies vanish at finite $q$, as a function of $\Omega/\Omega_{c}$. Crosses correspond to numerics and the red dashed line corresponds to the analytical prediction given by Eq. (38).

Appendix D Bogoliubov equations

From the energy functional (Eq. (1) in the main text), one can obtain the corresponding coupled Gross-Pitaevskii equations for the two wave functions $\psi_{i}$ as: $i\hbar\partial_{t}\psi_{i}=\delta\mathcal{E}/\delta\psi_{i}^{*}$. Static solutions (like the domain wall), obey $\psi_{i}=\psi_{i}\exp(i\mu t/\hbar)$, and thus can be computed by solving:

[TABLE]

Solutions of the previous equations can be directly seen in Figs. 1(c) in the main text and 3(a) in the Supplementary Material in the absence of the external trapping ($V_{\mathrm{ext}}=0$).

From the Gross-Pitaevskii equations, one derives the Bogoliubov equations by analyzing the fluctuations on the following static solution $\psi_{i}\to\psi_{i}+u_{i}\exp(i\omega t/\hbar)+v^{*}_{i}\exp(-i\omega t/\hbar)$, where $u_{i}(\vec{r}\,)$ and $v_{i}(\vec{r}\,)$ correspond to the amplitudes of the fluctuations. For the domain wall solution, both amplitudes can be factorized by symmetry reasons, as $u_{i}(\vec{r}\,)=\exp(iq\rho)\,u_{i}(x)$ and $v_{i}(\vec{r}\,)=\exp(iq\rho)\,v_{i}(x)$, respectively. By introducing the equilibrium state $\psi(x)$, numerically obtained by solving the Gross-Pitaevskii equation, the Bogoliubov equations

[TABLE]

can be solved by diagonalizing the matrix

[TABLE]

with the operators $B_{i}$ defined by :

[TABLE]

and the chemical potential given by

[TABLE]

Diagonalization of this matrix allows us to obtain the eigenfrequencies $\omega$, which turn out to be real for $\Omega<\Omega_{c}$ for any value of $q$ and imaginary for $\Omega>\Omega_{c}$ and $q$ below $q_{0}$, which depends on $\Omega$. The dependence of the $\mbox{Im}(\omega)$ as a function of $q$ is very similar to that appearing for the dark soliton, with a slope at $q\to 0$ increasing with $\Omega$. As a consequence, domain walls are dynamically unstable for $\Omega>\Omega_{c}$.

Appendix E Energy of a domain wall as a function of its length in a trap

In this section we study the dependence of the energy of a domain wall as a function of its length localized along one of the diametral axis of an isotropic harmonic trap. The authors of Ref. Eto and Nitta (2018) have found that in uniform matter, the critical length above which the domain wall is unstable, is much longer for the vortex-vortex than for the vortex-antivortex configuration. Under the presence of the harmonic trap, a new effect appears, playing a crucial role in the decay of the domain wall.

The space-dependence of the density in fact gives a strong modification of the energy as a function of the length, which is specially important when the gradient of the density is sharp, i.e. near the boundary. In practice the new effect results in an energy barrier that prevents the vortices that live at the edge of the condensate to enter inside the central region of higher density, unless a thermal or quantum tunnelling mechanism allows for such an entrance.

In Fig. 4 we show the dependence of the energy as a function of the length of a domain wall connecting two vortices in a vortex-vortex configuration in a gas of $8\times 10^{5}$ 23Na atoms confined in a harmonic trap with a trapping frequency equal to $\omega=2\pi\times 8\,\mbox{Hz}$, and $\Omega=1.2\,\Omega_{c}$. The figure clearly shows the existence of an energy barrier when the length of the domain wall approaches the size of the condensate. It can be understood in terms of the energy cost of introducing the core of a vortex into the higher-density region. A similar picture has been also calculated for the vortex-antivortex configuration, for which we find essentially identical results.