Flemish Strings of Magnetic Solitons and a Non-Thermal Fixed Point in a One-Dimensional Antiferromagnetic Spin-1 Bose Gas

TL;DR

This paper demonstrates that thermalization in a one-dimensional antiferromagnetic spin-1 Bose gas occurs via a non-thermal fixed point, involving the annihilation of Flemish-string magnetic soliton bound states, with potential experimental implications.

Contribution

It introduces the concept of Flemish-string bound states of magnetic solitons and their role in non-thermal fixed points during thermalization in a spin-1 Bose gas.

Findings

Thermalization proceeds via a non-thermal fixed point.

Flemish-string bound states are key to soliton annihilation.

Potential experimental setups are discussed.

Abstract

Thermalization in a quenched one-dimensional antiferromagnetic spin-1 Bose gas is shown to proceed via a non-thermal fixed point through annihilation of Flemish-string bound states of magnetic solitons. A possible experimental situation is discussed.

Equations

$$C(\mathbf{x},t)=t^{\gamma }f(\mathbf{x}{t}^{-\beta }),$$

$${\hat{H}}_{\mathrm{BH}}=$$

$$\displaystyle\mathcal{N}^{\rm AF}_{\mu\nu}=\frac{N}{2M}\left(\begin{array}[]{ccc}1-{\rm cos(2\alpha)}&-{\rm sin(2\alpha)}&0\\ -{\rm sin(2\alpha)}&1-{\rm cos(2\alpha)}&0\\ 0&0&2\end{array}\right),$$

$$C_{\mathrm{N}}(j,t):=\frac{1}{M}\sum _{k=0}^{M-1}⟨{\hat{N}}_{xy,j+k}{\hat{N}}_{xy,k}⟩(t).$$

$$i\hslash \frac{\partial }{\partial t}{b}_{m,j}=$$

$$i\frac{\partial }{\partial t}{\psi }_m(x,t)=-J^{\prime }\frac{{\partial }^2}{\partial x^2}{\psi }_m(x,t)-q^{\prime }{m}^2{\psi }_m(x,t)+U_0^{\prime }\rho (x,t){\psi }_m(x,t)+U_2^{\prime }\sum _{n=-1}^1\mathbf{S}(x,t)\cdot (\mathbf{s}{)}_{mn}{\psi }_n(x,t),$$

$$i\frac{\partial }{\partial t}{\psi }_1(x,t)=$$

$$i\frac{\partial }{\partial t}{\psi }_{-1}(x,t)=$$

$$\left(\begin{array}{c}{\psi }_1(x,t)\\ {\psi }_{-1}(x,t)\end{array}\right)=\sqrt{{\rho }_{\mathrm{tot}}(x,t)}{e}^{i\Phi (x,t)/2}\left(\begin{array}{c}\cos (\theta (x,t)/2)e^{-i\varphi (x,t)/2}\\ \sin (\theta (x,t)/2)e^{i\varphi (x,t)/2}\end{array}\right)$$

$$\frac{\partial }{\partial t}{\theta }^{\prime }+A\frac{\partial }{\partial x}{\theta }^{\prime }+B({\theta }^{\prime }{)}^2\frac{\partial }{\partial x}{\theta }^{\prime }+C\frac{{\partial }^3}{\partial x^3}{\theta }^{\prime }=0,$$

$$A=\sqrt{2J^{\prime }\bar{\rho }{U}_2^{\prime }},$$

$$B=-\frac{3A(8U_0^{\prime }+10U_2^{\prime })}{8(U_0^{\prime }-U_2^{\prime })},$$

$$C=-\frac{J^{\prime 3}}{2A},$$

$$i\frac{\partial }{\partial t}{\psi }_1(x,t)=$$

$$i\frac{\partial }{\partial t}{\psi }_{-1}(x,t)=$$

$$S_{z,\mathrm{tot}}=\frac{1}{M}\sum _{j=0}^{M-1}{S}_{z,j}.$$

$$i\hslash \frac{\partial }{\partial t}{\psi }_m(\mathbf{r},t)$$

$$\rho (\mathbf{r},t)=\sum _{m=-1}^1|{\psi }_m(\mathbf{r},t){|}^2,$$

$$\mathbf{S}(\mathbf{r},t)=\sum _{m,n=-1}^1{\psi }_m(\mathbf{r},t{)}^{\ast }(\mathbf{s}{)}_{mn}{\psi }_n(\mathbf{r},t),$$

$$V_{\mathrm{trap}}(\mathbf{r})=\frac{1}{2}M({\omega }_z^2{z}^2+{\omega }_r^2{r}^2),$$

$$c_0=\frac{4\pi {\hslash }^2(a_0+2a_2)}{3M},$$

$$c_2=\frac{4\pi {\hslash }^2(a_2-a_0)}{3M}.$$

$${\psi }_m(\mathbf{r},t)={\varphi }_m(z,t)F(r),$$

$$i\hslash \frac{\partial }{\partial t}{\varphi }_m(z,t)$$

$$\displaystyle c_{j}^{\prime}=c_{j}\Biggl{(}\int dxdy|F(r)|^{4}\Biggl{)}^{-1}~{}~{}(j=0,2).$$

$$F(r)=\{\begin{array}{rl}\sqrt{\frac{2}{\pi R_{\mathrm{eff}}^2}}\sqrt{1-\frac{r^2}{R_{\mathrm{eff}}^2}}& (r<R_{\mathrm{eff}});\\ 0(\mathrm{otherwise})& \end{array}$$

$$c_j^{\prime }=\frac{4c_j}{3\pi R_{\mathrm{eff}}^2}(j=0,2).$$

$$\displaystyle\bm{\Phi}(z)=\left(\begin{array}[]{c}0\\ \Phi_{0}(z)\\ 0\end{array}\right),$$

$$-\frac{{\hslash }^2}{2M}\frac{{\partial }^2}{\partial z^2}{\Phi }_0(z)+\frac{1}{2}M{\omega }_z^2{z}^2{\Phi }_0(z)+c_0^{\prime }|{\Phi }_0(z){|}^2{\Phi }_0(z)=\mu {\Phi }_0(z),$$

$$\displaystyle\phi_{1}(z)=\sum_{j=1}^{N_{{\rm cut},1}}\Bigl{(}R_{1,j}\alpha_{j}(z)+R_{-1,j}^{*}\beta_{j}(z)^{*}\Bigl{)},$$

Full text

Introduction.- Ultracold atomic gases offer an ideal playground for studying universal nonequilibrium dynamics due to their high controllability and isolation from the environment Pol11 ; Eisert15 ; Altman . Indeed, fundamental aspects of universal thermalization dynamics near critical points at equilibrium phase transitions have been studied in terms of the Kibble-Zurek mechanism KZ1 ; KZ2 ; KZ3 ; KZ4 ; KZ5 ; KZ6 and dynamical critical phenomena DPT1 . Even far from critical points, isolated systems are found to exhibit universal thermalization in the decay of turbulence NTF1 ; NTF2 and coarsening dynamics CDinBEC0 ; CDinBEC1 ; CDinBEC2 ; CDinBEC3 ; CDinBEC4 ; CDinBEC5 ; CDinBEC6 ; CDinBEC6_1 ; CDinBEC7 ; CDinBEC8 ; CDinBEC9 ; CDinBEC9_1 ; CDinBEC10 ; CDinBEC11 ; CDinBEC12 ; CDinBEC13 ; CDinBEC14 ; CDinBEC15 ; CDinBEC16 . However, unlike the dynamics near critical points DCP1 ; DCP2 , a unified framework for understanding the universal behavior is yet to be established.

Recently, a non-thermal fixed point (NTFP) NTF1 ; NTF2 ; NTF3 ; NTF4 has attracted growing interest as a universal thermalization scenario in isolated quantum systems. It is expected to provide a unified framework of different nonequilibrium phenomena in diverse systems ranging from cosmology to cold atoms NTF3 ; NTF4 . As illustrated in Fig. 1(a), a system with a strongly quenched initial state undergoes relaxation via a NTFP at which the scale-invariant thermalization dynamics emerges transiently. The NTFP is characterized by a universal scaling of an equal-time correlation function $C({\bm{x}},t)$ as $C({\bm{x}},t)=s^{-\gamma}C({\bm{x}}s^{\beta},ts)$ with an arbitrary parameter $s$ and scaling exponents $\beta$ and $\gamma$. Setting $s=t^{-1}$, we obtain the dynamical scaling:

[TABLE]

where $f({\bm{y}})=C({\bm{y}},1)$ is a universal function. After passing through the NTFP, the system relaxes toward a thermalized or stationary state. This scenario has been studied theoretically in two- and three-dimensional Bose gases NTF_GP1 ; NTF_GP2 ; NTF_GP3 ; NTF_GP4 ; NTF_GP5 ; NTF_GP6 ; NTF_GP7 and a relativistic $O(n)$ model NTF_O1 ; NTF_O2 ; NTF_O3 in the context of wave turbulence S (12, 13); WT3 and coarsening dynamics coar1 ; coar2 . Here weakly interacting waves or topological objects such as vortices and domain walls play essential roles in thermalization processes. Very recently, the universal scaling in Eq. (1) in momentum space was observed in one-dimensional (1D) ultracold atomic gases NTF_exp1 ; NTF_exp2 . However, it remains to be fully understood under what conditions a system shows universal behavior especially in an experimentally controllable manner. For example, a weak quench in Fig. 1(a) generates a small number of elementary excitations and may become thermalized without passing through a NTFP.

In this Letter, we theoretically study quench dynamics in a 1D spin-1 Bose gas with an antiferromagnetic (AF) interaction, unveiling universal thermalization dynamics through a NTFP caused by exotic bound states of twisted magnetic solitons, which we refer to as Flemish strings (see Figs. 1(b) and (c)). Since the spin-exchange interaction is AF, the relaxation dynamics after the quench in Fig. 2(a) is characterized by nematicity Nematic1 ; Nematic2 ; Nematic3 ; Nematic4 . We find that a nematic correlation function transiently exhibits the dynamical scaling (1) with $\beta\simeq 0.32$ and $\gamma\simeq 0.11$. The underlying physics of this universal thermalization through a NTFP is highly nontrivial cooperative soliton dynamics, where two magnetic solitons, which are stable individually magsoli1 ; magsoli2 ; magsoli3 ; magsoli4 , undergo pair-annhilation by forming Flemish strings. While there are a number of studies for quench dynamics of an AF spinor Bose gas CDinBEC9_1 ; quench0 ; quench1 ; quench2 ; quench3 ; quench4 ; quench5 ; quench6 ; quench7 ; quench8 ; quench9 ; quench10 , universal thermalization through Flemish strings has never been reported.

Our results reported here offer the first clear evidence of the soliton-induced NTFP. In Ref. NTF_pre , a NTFP in a 1D Bose gas is theoretically discussed, but the dynamical scaling (1), which is regarded as a hallmark of a NTFP NTF3 , is not confirmed. In Ref. NTF_exp2 , a relation between the observed universal thermalization dynamics and a random grey soliton model in a 1D repulsive Bose gas NTF_pre is discussed. However, the observed exponent $\beta$ deviates by 50% from the exponent obtained by assuming the theoretical model, and the other key exponent $\alpha=\beta+\gamma$ is not explained from the perspective of solitons. Thus, the relation between the NTFP and solitons has remained unclear.

Furthermore, we address an unsolved issue raised in Refs. AF_quench1 ; AF_quench2 , where thermalization dynamics after the same quench as ours in Fig. 2(a) was experimentally studied in trapped 1D AF Bose gases of ${}^{23}{\rm Na}$. The experiments investigate the observed dynamics in terms of coarsening coar1 ; coar2 , but universal aspects of thermalization have remained an unsolved issue. We argue that the universal thermalization dynamics characterized by the NTFP did not emerge in the experiments of Refs. AF_quench1 ; AF_quench2 because the quench was so weak that only a few solitons were created, and show that the NTFP due to Flemish strings should emerge for stronger quench.

Universal thermalization dynamics in a uniform system.- We consider a translation-invariant spin-1 Bose-Hubbard model SBH subject to the periodic boundary condition with the lattice constant $a$ and the number of lattices $M$. We denote the annihilation (creation) operator of spin-1 bosons with magnetic quantum number $m$ at site $j$ as $\hat{b}_{m,j}$ ($\hat{b}^{\dagger}_{m,j}$). Then, the Hamiltonian is given by

[TABLE]

with the hopping parameter $J$, the quadratic Zeeman coefficient $q$, the density interaction coefficient $U_{0}$, and the spin interaction coefficient $U_{2}$. Here, $\hat{\rho}_{j}\coloneqq\sum_{m}\hat{b}_{m,j}^{\dagger}\hat{b}_{m,j}$ is the total particle-number operator and $\hat{S}_{\alpha,j}\coloneqq\sum_{m,n}\hat{b}_{m,j}^{\dagger}({s_{\alpha}})_{mn}\hat{b}_{n,j}~{}(\alpha=x,y,z)$ with the spin-1 spin matrices $(s_{\alpha})_{mn}$ is the spin-vector operator.

We focus on the AF interaction ($U_{2}>0$), so that the mean-field ground-state phase is either polar ($q>0$) or AF ($q<0$) (see Fig. 2(a)). Since both phases are non-magnetic, the spinor order parameter is a second-rank tensor of the spin matrices Nematic1 ; Nematic2 ; Nematic3 described by the nematic tensor $\mathcal{N}_{\mu\nu,j}\coloneqq\langle\hat{N}_{\mu\nu,j}\rangle$ with \hat{N}_{\mu\nu,j}\coloneqq\frac{1}{2}\sum_{l,m,n}\hat{b}_{m,j}^{\dagger}\Bigl{[}(s_{\mu})_{mn}(s_{\nu})_{nl}+(\mu\leftrightarrow\nu)\Bigl{]}\hat{b}_{l,j}. Here, $\langle\cdots\rangle$ means the average over the state vector $|\psi(t)\rangle$. In the AF phase, the nematic tensor becomes

[TABLE]

where $\alpha$ is the azimuthal angle in spin space and $N$ is the number of condensed bosons. When we quench $q$ from the polar phase to the AF phase as shown in Fig. 2(a), the system gets highly excited and starts relaxation.

We employ the truncated Wigner approximation (TWA) TWA0 ; TWA1 ; TWA2 to numerically study the thermalization dynamics of a quenched AF Bose gas with $J/U_{0}=40$, $M=512$, and $N=2\times 10^{4}$, with which the system is in a deep superfluid regime where the TWA should work well TWA1 ; TWA2 . The initial state is chosen to be the Bogoliubov vacuum for the polar phase CDinBEC14 , and we suddenly quench $q$ from $q_{\rm i}=NU_{2}/M$ to $q_{\rm f}=(1-\Gamma)NU_{2}/M$ where $\Gamma$ controls the strength of the quench.

Figure 2(b) shows the spatio-temporal distributions of nematicity $\hat{N}_{xy,j}$ (left) and magnetization $\hat{S}_{z,j}$ (right) obtained from a single trajectory. The typical domain size of $\hat{N}_{xy,j}$ grows in time, while the number of locally magnetized domains, which are magnetic solitons magsoli1 ; magsoli2 ; magsoli3 ; magsoli4 , decreases in time (see Fig. 4(c)). As shown in Fig. 2(c), magnetic solitons can form a Flemish string in which two magnetic solitons form a virtual bound state and evolve in time in a twisted manner. The formation of a Flemish string is essential for thermalization as discussed later.

To investigate the universal scaling in Eq. (1), we calculate the nematic correlation function defined by

[TABLE]

Here, we consider the $xy$-component of $\hat{N}_{\mu\nu,j}$ because Eq. (6) shows that the order parameter of the AF phase can be characterized by $\alpha$ alone. The time evolution of $C_{{\rm N}}(j,t)$ shows that the correlation length $L_{\rm c}(t)$ monotonically increases as shown in Fig. 3 (a). Here, $L_{\rm c}(t)$ is determined from $C_{\rm N}(j=L_{c}(t)/a,t)=0.5C_{\rm N}(j=0,t)$. When the abscissa and the ordinate are normalized by $L_{\rm c}(t)/a$ and $C_{{\rm N}}(0,t)$, respectively, the correlation functions at different times collapse into a single universal curve as shown in Fig. 3(b). This dynamical scaling indicates the emergence of a NTFP.

The universality class of the NTFP can be classified by the exponents $\beta$ and $\gamma$ in Eq. (1), which are derived from power laws $L_{\rm c}(t)\propto t^{\beta}$ and $C_{\rm N}(0,t)\propto t^{\gamma}$ scaling . Figures 3(c) and (d) show $L_{\rm c}(t)$ and $C_{\rm N}(0,t)$ for four sets of parameters $U_{\rm 2}$ and $\Gamma$, from which we obtain $\beta\simeq 0.32$ and $\gamma\simeq 0.11$ by the least-square fit. These exponents are independent of the parameters and thus universal.

These exponents have never been reported in literature. To compare them with the wavenumber representation of previous studies NTF_exp1 ; NTF_exp2 ; NTF_GP1 ; NTF_GP2 ; NTF_GP3 ; NTF_GP4 ; NTF_GP5 ; NTF_GP6 ; NTF_GP7 , we perform the Fourier transformation of Eq. (1) in 1D systems, obtaining $\bar{C}(k,t)=t^{\alpha}g(kt^{\beta})$ with some function $g(x)$. Then, we find an exponent $\alpha=\gamma+\beta$ with $\alpha\simeq 0.43$, which does not agree with the results in Refs. NTF_exp1 ; NTF_exp2 ; NTF_GP1 ; NTF_GP2 ; NTF_GP3 ; NTF_GP4 ; NTF_GP5 ; NTF_GP6 ; NTF_GP7 .

The universality of the thermalization dynamics emerges from two different types of soliton solutions, i.e., magnetic solitons and Flemish strings. The magnetic soliton is a solution of the integrable Gardner equation Gardner1 , which can be derived from application of singular perturbation theory to the continuous classical model of Eq. (2) magsoli1 (see Supplemental Material supl ). As shown in Fig. 4(a), it has a locally magnetized stable object and separates the two different nematic orders, thereby decreasing the nematic correlation length. The thermalization dynamics would be promoted if a single magnetic soliton could be annihilated. However, it is, in fact, stable and cannot disappear spontaneously.

The annihilation of solitons proceeds through that of a Flemish string. This is a soliton solution of the integrable Manakov equation Manakov1 ; Manakov2 ; Manakov3 ; Manakov4 , which is a continuous classical model of Eq. (2) with $U_{2}=0$ (see Supplemental Material supl ). In contrast to magnetic solitons, it does not separate the nematic order as shown in Fig. 4(b). This string is unstable for $U_{2}\neq 0$, so that it eventually vanishes and nematic domains can grow. We have confirmed that Flemish strings are formed through collisions of magnetic solitons, as indicated by the solid curves in Fig. 2(b). The formation of Flemish strings and their subsequent decay constitute the main mechanism for soliton annihilation. The details including the lifetime and the number of Flemish strings are discussed in Supplemental Material supl .

Another numerical evidence for the nematic domain growth due to solitons is the number of magnetic solitons. We calculate $n_{\rm sol}\coloneqq\sum_{j=0}^{M-1}\langle\theta(\hat{S}_{z,j}^{2}-S_{\rm th}^{2})\rangle$ with $S_{\rm th}=0.7\rho_{\rm b}$ for the bulk particle number $\rho_{\rm b}$ soliton_number , where $\theta(x)$ is the unit-step function. This should be proportional to $t^{-\beta}$ if $L_{\rm c}(t)$ obeys the power law $t^{\beta}$. We find $n_{\rm sol}\propto t^{-0.31}$ as shown in Fig. 4(c), which is consistent with $\beta\simeq 0.32$ obtained from $L_{\rm c}(t)$. This confirms the fact that annihilation of magnetic solitons is essential in the universal thermalization dynamics through the NTFP. This universal behavior appears when many solitons are generated. However, for longer times, the universal relaxation stops and $L_{\rm c}(t)$ and $C_{\rm N}(0,t)$ saturate because the number of solitons becomes small.

Analysis of the experimental thermalization dynamics.- Let us examine the experiments AF_quench1 ; AF_quench2 where the thermalization dynamics in trapped 1D AF spinor Bose gases was investigated using density correlation functions. One is naturally led to ask whether or not the observed thermalization dynamics is universal. By numerically solving the spin-1 Gross-Pitaevskii (GP) equation with initial noises in the finite-temperature TWA method TWA1 ; FT_TWA1 ; FT_TWA2 ; FT_TWA3 , we obtain results in good agreement with the experiments supl . We calculate the correlation function $C_{\rm N}(z,t)$ with the parameters used in the experiments as shown in Fig. 5(a) definition_of_N . Although $C_{\rm N}(0,t)$ shows the power-law behavior in a short time ($200\sim 350~{}{\rm ms}$), we find neither the dynamical scaling of $C_{\rm N}(z,t)$ nor the power law of $L_{\rm c}(t)$. We therefore conclude that the universal thermalization dynamics does not emerge in the experiments.

The absence of the universal scaling in Eq. (1) is attributed to a weak quench in the experiments, where $q(t)$ is quenched from $2.8~{}{\rm Hz}$ to $-4.2~{}{\rm Hz}$ exp_comment . The energy scale of this change is given in terms of the spin interaction coefficient $c_{\rm 2}^{\prime}$ of the GP model and the bulk density $\rho_{\rm b,exp}$ by $0.06c_{\rm 2}^{\prime}\rho_{\rm b,exp}$, which is too small to excite a sufficient number of magnetic solitons.

To investigate the universal thermalization dynamics, we need a strong quench whose energy scale is of the order of the spin interaction energy $c_{\rm 2}^{\prime}\rho_{\rm b,exp}$. We numerically quench $q$ from $2.8~{}{\rm Hz}$ to $-117~{}{\rm Hz}$. Figure 5(b) shows the dynamical scaling in Eq. (1) and the power laws with the same exponents as those of the uniform system. Hence, we can confirm the universal thermalization dynamics through the NTFP. These results are consistent with the NTFP scenario in Fig. 1(a), where the strong quench leads to the NTFP but the weak one does not. While the spatial distribution of $\hat{N}_{xy,j}$ has not been observed, $\hat{S}_{z,j}$ can be measured experimentally. Thus, it is possible to estimate the exponent $\beta$ from $n_{\rm sol}$. The inset of Fig. 5(b-3) confirms $n_{\rm sol}\propto t^{-\beta}$, which is an experimentally observable signature of the NTFP.

$Conclusion.-$ By numerically studying the 1D thermalization in an antiferromagnetic spin-1 Bose gas, we find the universal nematic thermalization dynamics through a NTFP with scaling exponents $\beta\simeq 0.32$ and $\gamma\simeq 0.11$ in Eq. (1). We identify the thermalization mechanism to be annihilation of Flemish strings of magnetic soliton pairs. The universal thermalization dynamics is discussed in the experimental settings in Refs. AF_quench1 ; AF_quench2 . In 1D ultracold atomic gases, several different types of solitons other than magnetic solitons have experimentally been created soliton_exp1 ; soliton_exp2 ; soliton_exp3 ; soliton_exp4 ; soliton_exp5 ; soliton_exp6 ; soliton_exp7 ; soliton_exp8 ; soliton_exp9 and their dynamics has extensively been investigated atomic_gas_soliton ; soliton1 ; soliton2 ; soliton3 ; soliton4 ; soliton5 ; soliton6 ; soliton7 ; soliton8 ; soliton9 ; soliton10 ; soliton11 ; soliton12 ; soliton13 ; soliton14 ; soliton15 ; soliton16 . It is interesting to investigate other NTFP universality classes with these solitons.

Supplemental Material for “Flemish Strings of Magnetic Solitons and a Non-Thermal Fixed Point in a One-Dimensional Antiferromagnetic Spin-1 Bose Gas”

We describe integrable models for magnetic solitons and Flemish strings, discuss the main mechanism for magnetic-soliton annihilation, and provide details of our numerical analysis of the experiments of Refs. S (1, 2).

I Integrable models for magnetic solitons and Flemish strings

I.1 Continuum approximation

A classical trajectory in the truncated Wigner approximation (TWA) obeys the following classical equation of motion derived from Eq. (2) in the main text:

[TABLE]

where $b_{m,j}$, $\rho_{j}\coloneqq\sum_{m}b_{m,j}^{*}b_{m,j}$ and $S_{\alpha,j}\coloneqq\sum_{m,n}b_{m,j}^{*}({s_{\alpha}})_{mn}b_{n,j}~{}(\alpha=x,y,z)$ with the spin-1 spin matrices $(s_{\alpha})_{mn}$ are c-numbers corresponding to the annihilation operators $\hat{b}_{m,j}$, the total particle number, and the spin vector, respectively.

When the width of a soliton $l_{\rm s}$ is much larger than the lattice constant $a$, we can apply the continuum approximation to Eq. (S-1). In our numerical calculation for the thermalization dynamics in the main text, we have $l_{\rm s}\sim 9a$. Denoting a macroscopic wave function by $\psi_{m}(x=ja)=b_{m,j}$, we derive the following spinor Gross-Pitaevskii (GP) equation:

[TABLE]

where $J^{\prime}=Ja^{2}/\hbar$, $q^{\prime}=-q/\hbar$, $U_{0}^{\prime}=U_{0}/\hbar$, $U_{2}^{\prime}=U_{2}/\hbar$, $\rho(x,t)=\sum_{m=-1}^{1}|\psi_{m}(x,t)|^{2}$, and $\bm{S}(x,t)=\sum_{n,m=-1}^{1}\psi_{n}(x,t)^{*}({\bm{s}})_{nm}\psi_{m}(x,t)$. Here, we eliminate the chemical potential term proportional to $b_{m,j}$ by applying the global gauge transformation: $\psi(x,t)\rightarrow\psi(x,t){\rm exp}(-iWt)$ with $W=2J+2U_{0}-U_{2}$.

In the quench dynamics considered in the main text, the parameter $q$ is quenched from the polar phase to the antiferromagnetic phase, so that the particle number of the $m=0$ component is very small. Therefore, we can neglect $\psi_{0}(x,t)$ in Eq. (S-2), obtaining

[TABLE]

[TABLE]

These coupled equations have some soliton solutions as explained in the following.

I.2 The Gardner equation

Let us assume that the system has small but finite fluctuations for the amplitude and phase of $\psi_{m}(x,t)$. Then, we can derive an integrable equation called the Gardner equation using a singular perturbation method. We follow the derivation of Ref. S (3). Firstly, we write the wavefunction as

[TABLE]

with the local total particle number $\rho_{\rm tot}=|\psi_{1}|^{2}+|\psi_{-1}|^{2}$, the angle $\theta$ determining the relative particle number, the global phase $\Phi$, and the relative phase $\phi$. Secondly, by employing the singular perturbation method under the assumption that fluctuations of $\rho_{\rm tot}$, $\theta$, $\Phi$, and $\phi$ are small but nonzero, we obtain an effective equation of motion of $\theta^{\prime}=\theta-\pi/2$ given by

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and $\bar{\rho}$ is the average particle number. This equation is a special case of the integrable Gardner equation. It has a magnetic-soliton solution, which is called a polarization Gardner soliton in Ref. S (3). We note that the magnetic-soliton solution can be analytically derived without resort to the singular perturbation method if the local total particle number $\rho_{\rm tot}$ is independent of space S (4).

I.3 The Manakov equation

When we choose $U_{2}^{\prime}=0$, Eqs. (S-3) and (S-4) become

[TABLE]

[TABLE]

These coupled equations are integrable and known as the Manakov equation S (5). This equation has soliton solutions, one of which is a Flemish string. The mathematical expression of the solution is given in Ref. S (6)

II Mechanism for magnetic-soliton annihilation

We discuss the main mechanism for magnetic-soliton annihilation in our universal relaxation. First, we point out the crucial role played by spin conservation and integrability breaking in the spin-1 Bose-Hubbard model, and then explain why formation of Flemish strings through magnetic-soliton collisions and their subsequent decay constitute the main mechanism for magnetic-soliton annihilation. Second, we numerically confirm this annihilation process by counting the number of Flemish strings and nematic domains.

Soliton annihilation requires collision processes due to the stability of a single magnetic soliton. Our numerical calculations show that a single magnetic soliton is quite stable, and that any spontaneous single-soliton annihilation occurs. If we strongly excite the system with energy injection much larger than the spin exchange interaction, magnetic solitons might be broken into large amplitude spin waves. However, we do not inject such high energy into the system in the relaxation dynamics considered here, so that the spontaneous soliton annihilation should not occur. Thus, collision processes are needed for soliton annihilation.

First, let us make a physical argument about the annihilation process in light of spin conservation. Our system respects spin rotational symmetry about the z-axis, so that the spatially averaged magnetization $S_{z,{\rm tot}}$, defined in the following, is conserved:

[TABLE]

Thus, collision processes of two magnetic solitons with the same sign of $S_{z}$ (see Figs. S-1 (a) and (b)) are generally forbidden. As shown below, these processes are indeed not found in our numerical simulations. On the other hand, a collision process with opposite signs of $S_{z}$ as in Fig. S-1(c) can lead to annihilation of two magnetic solitons without breaking spin conservation. Therefore, spin conservation imposes a strong constraint on the soliton-annihilation mechanism.

It is notable that Flemish strings can be formed through the collision process in Fig. S-1(c) before soliton-pair annihilation. This is because our model is close to the integrable Manakov model as described in the main text and this model has an exact Flemish-string solution. However, our model is not exactly integrable due to the weak spin-exchange interaction, and the Flemish string can disappear with a finite lifetime. Thus, this integrability breaking eventually leads to the soliton-pair annihilation. From this discussion, the time scale for the instability of Flemish strings is roughly estimated to be $\hbar/U_{2}$ with the spin-exchange interaction energy $U_{2}$. The quantitative lifetime might be analytically obtained by implementing soliton perturbation calculations S (7), from which we may derive an analytical expression of the lifetime. However, the investigation of this problem is far beyond the scope of the present work and we would like to leave it as a future problem.

Summarizing the above discussion, we conclude that the main mechanism for the soliton annihilation is the collision process between the positive and negative magnetic-solitons that form unstable Flemish strings and the subsequent decay. In this mechanism, the spin conservation and the integrability breaking for the integrable Manakov model play essential roles.

We perform numerical simulations to investigate the annihilation processes of magnetic solitons, and find that two magnetic solitons with opposite signs disappear only through formation of a Flemish string. In the numerical calculation, preparing an initial state with 8 magnetic solitons (8 nematic domains), we solve the classical equation corresponding to the spin-1 Bose-Hubbard model with $U_{0}/J=0.025$ and $U_{2}/J=0.0025$ under a periodic boundary condition. The time is normalized by the characteristic time $\tau=4\hbar/J$. Figure S-2 shows spatio-temporal distributions of nematicity and spin, where there are initially 8 nematic domains divided by magnetic solitons. As time goes by, Flemish strings, indicated by arrows, are formed. Counting the number of nematic domains under the periodic boundary condition, one can confirm that the number decreases via formation of Flemish strings. For example, the domain number is 6 around $t/\tau=3200$ after the first formation of a Flemish string ($t/\tau=2900$, solid arrows) and decreases to 4 around $t/\tau=4400$ after the second formation ($t/\tau=3500$, dashed arrows). The spin amplitudes of Flemish strings become smaller in time and gradually disappear. Here, the time scale in which the amplitude starts to change is of the order of $\hbar/U_{2}\sim 100\tau$, which is consistent with our discussion of the Flemish-string instability based on the integrability breaking caused by the weak spin-exchange interaction. Our numerical calculations show no evidence of other annihilation mechanisms for magnetic solitons.

Finally, we discuss the number of Flemish strings. Although it is difficult to automatically identify the number by a computer, we can count the number by eye and confirm that the nematic order indeed becomes larger as the Flemish strings are formed and subsequently decay as shown in Fig. S-2. For example, one can find that two Flemish strings are formed around $t/\tau=2900,3500$, and that their subsequent decay leads to a decrease in the number of solitons. This is a clear evidence for the growth of nematic order through formation of Flemish strings. The reason for the difficulty of the automatic identification by a computer is that we need to scrutinize the time evolution of magnetic solitons to judge whether they are bound or not but we have yet to understand how to automatically judge twisting of Flemish strings.

III Numerical analysis for the experiments

III.1 Model

We consider a Bose-Einstein condensate of ${}^{23}{\rm Na}$ with mass $M$ in a cigar-shaped harmonic trapping potential. To investigate the quench dynamics in this system, we use the GP equation given by

[TABLE]

with the macroscopic wavefunction $\psi_{m}(\bm{r},t)~{}(m=-1,0,1)$ and the quadratic Zeeman coefficient $q(t)$. The total particle density $\rho({\bm{r}},t)$, the spin density vector $\bm{S}({\bm{r}},t)$, and the trapping potential $V_{\rm trap}(\bm{r})$ are given by

[TABLE]

[TABLE]

[TABLE]

where $r=\sqrt{x^{2}+y^{2}}$ is the radial coordinate. The parameters $c_{0}$ and $c_{2}$ are the strengths of spin-independent and spin-dependent interactions, which are expressed in terms of the s-wave scattering lengths $a_{0}$ and $a_{2}$ for the total-spin 0 and 2 channels by

[TABLE]

[TABLE]

We adopt the following parameters used in the experiments of Refs S (1, 2): $\omega_{z}=2\pi\times 7~{}{\rm Hz}$, $\omega_{r}=2\pi\times 470~{}{\rm Hz}$, and $M=23{\rm u}$, where $\rm u$ is the unified atomic mass unit. We choose $N_{\rm c}=9.4\times 10^{6}$, $a_{0}=48.8a_{\rm B}$, and $a_{2}=51.2a_{\rm B}$ with the Bohr radius $a_{\rm B}$ such that the Thomas-Fermi radii $R_{r}$ and $R_{z}$ in the radial and axial directions, the bulk chemical potential $\mu_{\rm b}$, and the bulk spin interaction energy $\mu_{\rm s}$ agree with the experimentally observed values. Using these parameters, we obtain $R_{r}\sim 5.4~{}{\rm\mu m}$, $R_{z}\sim 360~{}{\rm\mu m}$, $\mu_{\rm b}\sim 354~{}{\rm nK}$, and $\mu_{\rm s}\sim h\times 120~{}{\rm Hz}$, where $h$ is the Planck constant. These values agree well with those observed experimentally S (1, 2).

III.2 Reduction of spatial dimensions to one

Our system has a cigar-shaped configuration, so that we can reduce the three-dimensional system to a one-dimensional (1D) one. Because of tight confinement in the radial direction, we can express the wavefunction as

[TABLE]

where $\phi_{m}(z,t)$ is the one-dimensional wave function with magnetic quantum number $m$ and $F(r)$ is the radial wave function. As described in Ref. S (1), when the energy scale of the quench is smaller than that of the first excited state in the radial direction, we can show that $F(r)$ is independent of time.

Substituting Eq. (S-19) into Eq. (S-13), we obtain the 1D GP equation given by

[TABLE]

where $\rho_{\rm 1d}(z,t)=\sum_{m=-1}^{1}|\phi_{m}(z,t)|^{2}$ and $\bm{S}_{\rm 1d}(z,t)=\sum_{m,n=-1}^{1}\phi_{m}(z,t)^{*}(\bm{s})_{mn}\phi_{n}(z,t)$. The 1D interaction parameters are given by

[TABLE]

Finally, we consider the form of $F(r)$. The radial size of the condensate is $R_{r}\sim 5.4~{}{\rm\mu m}$, which is much larger than the harmonic-oscillator length $a_{r}=\sqrt{\hbar/2M\omega_{r}}\sim 0.7~{}{\rm\mu m}$. Thus, we assume the Thomas-Fermi approximation for $F(r)$:

[TABLE]

with the effective radius $R_{\rm eff}$. As a result, we obtain

[TABLE]

Here, we use $R_{\rm eff}=4.05~{}{\rm\mu m}\sim R_{r}$ which makes the physical quantities, such as the bulk chemical potential, agree well with those of the experiments in Refs. S (1, 2).

III.3 How to generate initial states

To solve Eq. (S-20), we take the polar state as an initial state, which is generated by using the finite-temperature TWA method S (8, 9, 10, 11, 12).

First, we numerically obtain a mean-field ground state by the imaginary-time evolution of Eq. (S-20):

[TABLE]

which satisfies

[TABLE]

where $\mu$ is the chemical potential. Then, the initial state for Eq. (S-20) can be expressed as

[TABLE]

[TABLE]

[TABLE]

where $N_{{\rm cut},m}$ is a cutoff to avoid the ultraviolet divergence, which is unavoidable when the TWA method is applied to continuous models S (8, 9). In our numerical calculation, the cutoff $N_{{\rm cut},m}$ is determined so that $E_{m,j}<\mu$ is satisfied. The functions $\alpha_{j}(z)$, $\beta_{j}(z)$, $u_{j}(z)$, and $v_{j}(z)$ obey the following Bogoliubov-de Gennes equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The quantities $R_{1,j}$, $R_{0,j}$, and $R_{-1,j}$ are random numbers which are generated by the Wigner function defined by

[TABLE]

with $E_{-1,j}=E_{1,j}$.

III.4 Comparison between experimental and numerical results

We perform numerical calculations of Eq. (S-20) with Eqs. (S-27)-(S-36), and compare the results with the experimental ones. Here, we compare the sum of the normalized particle numbers of the $m=1$ and $m=-1$ components and density fluctuations as a function of time. The former is defined by

[TABLE]

where $N_{\rm tot}$ is the total particle number and $N_{m}(t)$ is the particle number of the $m$-component. The latter are defined by

[TABLE]

where ${\rho}_{m}(z)=|\phi_{m}(z)|^{2}$ is the particle density for the $m$-component, and the bracket denotes the ensemble average.

Figure S-3 shows the time evolution of $\mathcal{C}_{\rm sum}(t)$ and its dependence on temperature $T$. The temperature reported in the experiments of Refs. S (1, 2) is $T_{\rm exp}=400~{}{\rm nK}$. We find that our numerical results at $200~{}{\rm nK}$ and $400~{}{\rm nK}$ show better agreement with the experimental data than that at the zero-temperature result in the early stage ($t<30~{}{\rm ms}$), but in the later stage ($t>60~{}$ms) the results at $T=0~{}{\rm nK}$, $200~{}{\rm nK}$, and $400~{}{\rm nK}$ agree qualitatively with the data.

The numerical results of $\mathcal{D}_{m,n}(t)$ are shown in Fig. S-4, where we plot the two results at $T=0~{}{\rm nK}$ and $400~{}{\rm nK}$. In contrast to Fig. S-3, the density fluctuations $\mathcal{D}_{m,n}(t)$ at = $400~{}{\rm nK}$ show the better agreement with the experimental results than those at zero temperature. Compared with the previous numerical results (Fig. 3 in Ref. S (1)), the agreement improves considerably. We suspect that the remaining deviation between our numerical results and the experimental ones may arise from the time-of-flight Stern-Gerlach observation in the experiments because it takes about $30~{}{\rm ms}$, which is comparable with the time scale needed before $\mathcal{C}_{\rm sum}(t)$ and $\mathcal{D}_{m,n}(t)$ in Figs. (S-3) and (S-4) start to grow.

In summary, our numerical calculation (especially for $T=400~{}{\rm nK}$) can describe the experimental result well if we use the finite-temperature TWA method. Since the parameter region is in the deep superfluid regime where the method works well, we expect that our numerics can predict the thermalization dynamics observed in the experiments qualitatively beyond the experimental time scale $\sim 100~{}{\rm ms}$ as well.

III.5 Spatio-temporal distribution of $S_{z}$ and $N_{xy}$ in the quench dynamics

We perform numerical calculations with two quench protocols: weak quench corresponding to the experiments of Refs. S (1, 2) and strong quench. We suddenly change $q$ from $2.8~{}{\rm Hz}$ to $-4.2~{}{\rm Hz}$ in the former case, and from $2.8~{}{\rm Hz}$ to $-117~{}{\rm Hz}$ in the latter case, and examine the spatio-temporal distributions of $S_{{\rm 1d},z}(z,t)$ and $N_{{\rm 1d},xy}(z,t)$ for both cases, which are calculated from a single trajectory of Eq. (S-20). Here, $N_{{\rm 1d},xy}(z,t)$ is given by

[TABLE]

[TABLE]

Figure S-5 illustrates the spatio-temporal distributions for the weak quench. In this case, the energy scale of the quench $\sim 0.06c_{2}^{\prime}\rho_{\rm b,exp}$ is much smaller than the spin-exchange interaction energy $c_{2}^{\prime}\rho_{\rm b,exp}$ at the bulk density $\rho_{\rm b,exp}$, so that only a few magnetic solitons are generated. On the other hand, for strong quench, many magnetic solitons are created as shown in Fig. S-6 because the energy scale for this quench is about $c_{2}^{\prime}\rho_{\rm b,exp}$.

III.6 Discussion on the density correlation function observed in Ref. S (2)

We discuss how the density correlation function observed in Ref. S (2) is related to our universal relaxation dynamics. The experiment is a sequel to Ref. S (1), and studied the relaxation dynamics in an antiferromagnetic phase after quenching the quadratic Zeeman coupling $q$ from positive to negative values. By observing a spatial correlation function for a density distribution, the experiment confirmed that the correlation function exhibits relaxation to a stationary state.

We conclude that our universal relaxation and the power exponents are not directly related to the experimentally observed results because of the following reason. The experiment observed relaxation after the weak quench, and generated only a few magnetic solitons as shown in our numerical calculations in Fig. S-5. Thus, the observed relaxation is essentially different from our universal relaxation with many magnetic solitons, and there is no direct relation between our power exponents and the observed density correlation function.

Finally, let us briefly discuss our system after a long time, although this is beyond the scope of our research. In this case, our system has only a few magnetic solitons, and becomes similar to the experimental situation. Thus, on this time scale, a relation might be established between the experiment and our system. Reference S (2) has reported that the density correlation length decreases in time (the exponent $\beta$ is negative), which implies growth of fine density structures. By taking into account that only a few magnetic solitons are generated, this behavior reminds us of a cascade in wave turbulence S (12, 13), in which nonlinear interactions between waves with different wavenumbers generate finer structures. In the field of ultracold gases, several theoretical and experimental studies have found such turbulent cascades S (14, 15, 16, 17, 18, 19, 20, 21, 22). We thus expect that a signature of a wave turbulent cascade may appear for few soliton cases in experiments or in the long-time regime of our setup. Applying weak wave turbulence theory S (12, 13) to the spin-1 BH model, we could analytically derive several power exponents for density, spin, and nematic correlation functions, which might provide direct relation between the experimental and our theoretical results.