Solitons in inhomogeneous gauge potentials: integrable and nonintegrable dynamics

TL;DR

This paper introduces an exactly integrable nonlinear model for spinor solitons in space-dependent gauge potentials, explores the effects of Zeeman splitting on soliton dynamics, and demonstrates a controllable transition between integrable and nonintegrable regimes.

Contribution

The work presents a new integrable model for spinor solitons in inhomogeneous gauge potentials and analyzes how Zeeman splitting controls the transition between integrable and nonintegrable dynamics.

Findings

At zero Zeeman splitting, solitons move without scattering in a disordered SOC landscape.

Nonzero Zeeman splitting causes strong scattering and nonintegrable behavior.

Large Zeeman splitting restores integrability, enabling control over soliton dynamics.

Abstract

We introduce an exactly integrable nonlinear model describing the dynamics of spinor solitons in space-dependent matrix gauge potentials of rather general types. The model is shown to be gauge equivalent to the integrable system of vector nonlinear Schr\"odinger equations known as the Manakov model. As an example we consider a self-attractive Bose-Einstein condensate with random spin-orbit coupling (SOC). If Zeeman splitting is also included, the system becomes nonintegrable. We illustrate this by considering the random walk of a soliton in a disordered SOC landscape. While at zero Zeeman splitting the soliton moves without scattering along linear trajectories in the random SOC landscape, at nonzero splitting it exhibits strong scattering by the SOC inhomogeneities. For a large Zeeman splitting the integrability is recovered. In this sense the Zeeman splitting serves as a parameter controlling the crossover between two different integrable limits.

Equations

$$i\frac{\partial \mathbf{\psi }}{\partial t}=\frac{1}{2}{\mathrm{P}}^2(x)\mathbf{\psi }+\frac{\Omega }{2}{\sigma }_3\mathbf{\psi }-\left({\mathbf{\psi }}^{†}\mathbf{\psi }\right)\mathbf{\psi }.$$

$$i{\bm{u}}_{t}+\frac{1}{2}{\bm{u}}_{xx}+({\bm{u}}^{\dagger}{\bm{u}}){\bm{u}}=\Omega\kappa(x){\bm{u}},\quad\kappa=\left(\begin{array}[]{cc}s_{3}&\kappa_{3}\\ \kappa_{3}^{*}&-s_{3}\end{array}\right),$$

$$\displaystyle\mbox{\boldmath$\varphi$\unboldmath}_{x}=U\mbox{\boldmath$\varphi$\unboldmath}-i\lambda\mbox{\boldmath$\varphi$\unboldmath}{\cal E},\quad\mbox{with}\,\,U=i(\lambda{\cal E}+{\cal A}+{\cal U})$$

$${\cal A}=\left(\begin{array}[]{ccc}a_{3}&a_{1}+ia_{2}&0\\ a_{1}-ia_{2}&-a_{3}&0\\ 0&0&0\end{array}\right)\!,\,\,{\cal U}=\left(\begin{array}[]{ccc}0&0&\psi_{1}^{*}\\ 0&0&\psi_{2}^{*}\\ \psi_{1}&\psi_{2}&0\end{array}\right).$$

$$\frac{dX}{dt}=\Pi ,\frac{d\Pi }{dt}=i\frac{\Omega }{2\Vert \mathbf{\psi }{\Vert }^2}\int {\mathbf{\psi }}^{†}[{\sigma }_3,\mathrm{A}]\mathbf{\psi }dx(\equiv F(t))$$

$$i\frac{\partial \mathbf{\psi }}{\partial t}=\frac{1}{2}{\mathbf{P}}^2(\mathbf{r})\mathbf{\psi }-\left({\mathbf{\psi }}^{†}\mathbf{\psi }\right)\mathbf{\psi }.$$

$$\mathbf{P}=-i\nabla +\mathbf{A}(\mathbf{r})$$

$$\mathbf{P}{\xi }_{\mathbf{k}}(\mathbf{r})=\mathbf{k}{\xi }_{\mathbf{k}}(\mathbf{r})$$

$$\displaystyle-i\nabla\mbox{\boldmath$\phi$\unboldmath}_{1}+{\bm{A}}({\bm{r}})\mbox{\boldmath$\phi$\unboldmath}_{1}=0.$$

$$i\frac{\partial A_x}{\partial y}-i\frac{\partial A_y}{\partial x}+[A_x,A_y]=0,$$

$$\displaystyle\nabla(\mbox{\boldmath$\phi$\unboldmath}_{1}^{\dagger}\mbox{\boldmath$\phi$\unboldmath}_{1})=\nabla(\mbox{\boldmath$\phi$\unboldmath}_{2}^{\dagger}\mbox{\boldmath$\phi$\unboldmath}_{2})=0$$

$$\displaystyle\mbox{\boldmath$\phi$\unboldmath}_{1}^{\dagger}\mbox{\boldmath$\phi$\unboldmath}_{2}=0.$$

$$\displaystyle\bm{\psi}({\bm{r}},t)=u_{1}({\bm{r}},t)\mbox{\boldmath$\phi$\unboldmath}_{1}({\bm{r}})+u_{2}({\bm{r}},t)\mbox{\boldmath$\phi$\unboldmath}_{2}({\bm{r}})$$

$$i{\mathbf{u}}_t+\frac{1}{2}{\nabla }^2\mathbf{u}+({\mathbf{u}}^{†}\mathbf{u})\mathbf{u}=0$$

$$A_x=\mathbf{a}(x)\cdot \mathbf{\sigma },A_y=\mathbf{b}(x)\cdot \mathbf{\sigma }$$

$$\frac{d\mathbf{b}}{dx}=2\mathbf{a}\times \mathbf{b}.$$

$$\displaystyle\mbox{\boldmath$\phi$\unboldmath}_{1}({\bm{r}})={\tilde{}\mbox{\boldmath$\phi$\unboldmath}}_{1}(x)e^{-iqy}$$

$$\displaystyle\frac{1}{i}\frac{d{\tilde{}\mbox{\boldmath$\phi$\unboldmath}}_{1}}{dx}+A_{x}(x){\tilde{}\mbox{\boldmath$\phi$\unboldmath}}_{1}=0,\quad A_{y}(x){\tilde{}\mbox{\boldmath$\phi$\unboldmath}}_{1}=q{\tilde{}\mbox{\boldmath$\phi$\unboldmath}}_{1},$$

$$H_0=-\hslash ({\Omega }_1|0⟩⟨1|+{\Omega }_2|0⟩⟨2|+{\Omega }_3|0⟩⟨3|)+h.c.$$

$$\displaystyle\begin{array}[]{c}\displaystyle{\Omega_{1}=\Omega e^{i\Theta({\bm{r}})}e^{ik_{x}x+ik_{y}y}}\\[8.53581pt] \displaystyle{\Omega_{2}=\Omega e^{i\Theta({\bm{r}})}e^{-ik_{x}x+ik_{y}y}}\\[8.53581pt] \displaystyle{\Omega_{3}=2\Omega e^{ikz}}\end{array}$$

$$e^{i\tilde{k}\xi }{e}^{i\Theta (\xi ,\eta ,\zeta )}={\int }_{k_0-\delta k}^{k_0+\delta k}d{k}_{\perp }{e}^{-i\sqrt{k^2-k_{\perp }^2}\xi }{\int }_{-\pi }^{\pi }d\nu \mathcal{A}(\nu ,k_{\perp })e^{i{k}_{\perp }(\eta \cos \nu +\zeta \sin \nu )},$$

$$|D_1⟩=\frac{1}{\sqrt{2}}\left\{e^{-i{k}_{x}x-i{k}_{y}y}|1⟩-e^{i{k}_{x}x-i{k}_{y}y}|2⟩\right\},$$

$$|D_2⟩=\frac{1}{\sqrt{3}}\left\{e^{-i\Theta (\mathbf{r})}{e}^{-i{k}_{x}x-i{k}_{y}y}|1⟩+e^{-i\Theta (\mathbf{r})}{e}^{i{k}_{x}x-i{k}_{y}y}|2⟩-e^{-ikz}|3⟩\right\}.$$

$$|\Psi ⟩={\Psi }_1(\mathbf{r})|D_1⟩+{\Psi }_2(\mathbf{r})|D_2⟩.$$

$$H_{\mathrm{lin}}=\frac{1}{2M}{\left(\frac{\hslash }{i}\nabla -A\right)}^2+V_{\mathrm{tot}}(\mathbf{r}),$$

$$A_{mn}=i\hslash ⟨D_m(\mathbf{r})|\nabla D_n(\mathbf{r})⟩.$$

$$U_{ij}=\frac{{\hslash }^2}{2M}\left(⟨\nabla D_i|\nabla D_j⟩+\sum _{l=1}^2⟨D_i|\nabla D_l⟩⟨D_l|\nabla D_j⟩\right).$$

$$A=\alpha {\sigma }_1{e}^{i{\sigma }_3\theta (x)},\alpha =\sqrt{\frac{2}{3}}{a}_{\perp }{k}_x,$$

$$U=\frac{{\hslash }^2}{6M}\text{mbox}diag\left(\frac{k_x^2}{2},\frac{k_y^2}{3}-\frac{k^2}{6}\right).$$

Full text

Gauge invariance having its origins in the theory of electromagnetism, is known to be a general principle playing crucial role in almost any field history . One of its applications, intensively discussed nowadays are the synthetic gauge fields and potentials. Such potentials of practically arbitrary form can be created in atomic systems, illuminated by proper combination of the laser beams Ruseckas . In this way it is possible to emulate in systems of neutral atoms analogs of electric electric and magnetic magnetic fields, as well as the spin-orbit coupling (SOC). The latter technique has recently made possible to engineer spin-orbit coupled Bose-Einstein condensates (SO-BECs) Nature . Two important properties, the tunability of the SOC in atomic systems tun-SO1 ; tun-SO2 ; tun-SO3 ; tun-SO4 , as well as the intrinsic nonlinearity of SO-BECs stemming from inter-atomic interactions, have stimulated extensive studies of soliton dynamics in BECs with inhomogeneous SOC. In particular, the interactions of one-dimensional (1D) solitons in SO-BEC with a localized coupling defect have been studied in 1D KarKonZez14 and in 2D Malomed_2Ddefect settings, and the propagation of soliton in a BEC with inhomogeneous helicoidal SOC was addressed KarKon2017 .

The gauge invariance is also known to be a powerful tool of generating and studying nonlinear integrable systems FadTak . In particular, two integrable nonlinear equations describing different physical phenomena, may be found to be gauge equivalent, i.e. reducible to each other by a gauge transformation. Also, gauge transformation can be used to generate continuous cont-gauge and discrete KonChubVaz-1 ; KonChubVaz-2 integrable models with inhomogeneous coefficients departing from the homogeneous ones.

It was found previously KarKon2017 that if a quasi-1D BEC with equal inter– and intra–component interactions has a helicoidal structure, and no other potentials or Zeeman splitting is present, the dynamics of soliton is reduced to the soliton of the exactly integrable Manakov model Manakov which is a system of nonlinearly coupled SU(2) invariant nonlinear Schrödinger (NLS) equations. A similar result was pointed out in KarKonZez14 for the case of a particular inhomogeneous gauge potential. The possibility of generalizing these results for arbitrary potentials still remains an open question.

In the present Letter we prove that coupled NLS equations with $x-$dependent matrix Hermitian gauge potential of a general type, is an integrable model, which is gauge equivalent to the Manakov system. The inclusion of the Zeeman splitting, quantified below by the field $\Omega$, makes the system non-integrable at finite values of $\Omega$, while its integrability is restored in the limit $\Omega\to\infty$. The strength of the Zeeman field is a parameter describing a crossover between two different integrable limits.

To study the above mentioned crossover, we address the evolution of a matter soliton in a BEC with a random SOC (for recent study of random SOC in linear systems see randSOC ). This is the second goal of this Letter. In particular, we show that the gauge transformation in the integrable case effectively separates random evolution of the pseudo-spin and deterministic evolution of the soliton envelope with random initial conditions, similarly to transformation between random and regular SOC explored in the linear theory gauge-linear . The Zeeman field couples these dynamical processes leading to anomalous diffusion of a soliton.

Let us consider a 1D Gross-Pitaevskii equation (GPE) describing a spinor $\bm{\psi}(x,t)=(\psi_{1},\psi_{2})^{T}$ ($T$ stands for transpose) in the presence of a random gauge potential ${\rm A}(x)$ and of a Zeeman coupling $\Omega\sigma_{3}/2$:

[TABLE]

Here ${\rm P}=-i\partial_{x}+{\rm A}(x)$ is a generalized momentum, and $\sigma_{1,2,3}$ are the Pauli matrices. Inter- and intra-species interactions are assumed to be attractive and equal. The units are chosen to make the atomic mass $M=\hbar=1$, and $\bm{\psi}$ is normalized to have nonlinear coefficient equal to 1.

We impose three constraints on the $x$-dependent gauge field, requiring it to be Hermitian ${\rm A}={\rm A}^{\dagger}$, as needed for Hermiticity of the generalized momentum, to anti-commute with the time reversal operator for spin $1/2$ particles ${\cal T}=i\sigma_{2}{\cal K}$, where ${\cal K}$ is the complex conjugation: ${\rm A}{\cal T}+{\cal T}{\rm A}=0$, and to have determinant equal to a constant, characterizing the SOC strength: $\det{\rm A}=-\alpha^{2}$. The first two requirements define the general form of the gauge field: {\rm A}(x)=\mbox{\boldmath\sigma\unboldmath}{\bm{a}}, where ${\bm{a}}=(a_{1},a_{2},a_{3})$ is a real vector and \mbox{\boldmath\sigma\unboldmath}=(\sigma_{1},\sigma_{2},\sigma_{3}) is the Pauli matrix vector. Next we consider an eigenstate ${\bm{\xi}}_{k}(x)$ of the ${\rm P}$ operator: ${\rm P}{\bm{\xi}}_{k}=k{\bm{\xi}}_{k}$, $k$ being the eigenvalue. We also define \mbox{\boldmath\phi\unboldmath}_{1}(x)=e^{-ikx}{\bm{\xi}}_{k} and \mbox{\boldmath\phi\unboldmath}_{2}={\cal T}\mbox{\boldmath\phi\unboldmath}_{1}. It is straightforward to verify that {\rm P}\mbox{\boldmath\phi\unboldmath}_{j}=0 and \mbox{\boldmath\phi\unboldmath}_{j}^{\dagger}\mbox{\boldmath\phi\unboldmath}_{j} is $x$-independent. Thus the spinors \mbox{\boldmath\phi\unboldmath}_{1,2}(x) make up an orthonormal basis in the spinor subspace: \mbox{\boldmath\phi\unboldmath}_{i}^{\dagger}\mbox{\boldmath\phi\unboldmath}_{j}=\delta_{ij}. The vectors \mbox{\boldmath\phi\unboldmath}_{1} and \mbox{\boldmath\phi\unboldmath}_{2} describe opposite pseudo-spin distributions \tilde{\mbox{\boldmath\sigma\unboldmath}}(x) and -\tilde{\mbox{\boldmath\sigma\unboldmath}}(x), where \tilde{\mbox{\boldmath\sigma\unboldmath}}=(\tilde{\sigma}_{1},\tilde{\sigma}_{2},\tilde{\sigma}_{3}) with \tilde{\sigma}_{j}(x)=\frac{1}{2}\mbox{\boldmath\phi\unboldmath}_{1}^{\dagger}\sigma_{j}\mbox{\boldmath\phi\unboldmath}_{1} $(j=1,2,3)$.

Using the basis \mbox{\boldmath\phi\unboldmath}_{1,2}(x) one can write the solution of (1) as \bm{\psi}(x,t)=u_{1}(x,t)\mbox{\boldmath\phi\unboldmath}_{1}(x)+u_{2}(x,t)\mbox{\boldmath\phi\unboldmath}_{2}(x), and verify that the “envelope” spinor ${\bm{u}}=(u_{1},u_{2})^{T}$ solves the equation

[TABLE]

where \kappa_{3}(x)=\frac{1}{2}\mbox{\boldmath\phi\unboldmath}_{1}^{\dagger}\sigma_{3}\mbox{\boldmath\phi\unboldmath}_{2} describes the coupling of the envelope components. This allows us to interpret the vectors \mbox{\boldmath\phi\unboldmath}_{1,2} and ${\bm{u}}$ respectively as the pseudo-spin and the soliton “degrees” of freedom, which are coupled by the Zeeman field when $\Omega\neq 0$. In the absence of Zeeman field, $\Omega=0$, Eq. (2) describes deterministic evolution of the envelope ${\bm{u}}$. We note, that even in the case of deterministic initial condition for the field $\bm{\psi}_{0}(x)=\bm{\psi}(x,0)$, the initial conditions for the fields $u_{1,2}(x,t)$ are random: they are defined by the projections of $\bm{\psi}_{0}(x)$ on \mbox{\boldmath\phi\unboldmath}_{1,2}(x).

Remarkably, the described separation of spinor and nonlinear degrees of freedom can be performed also in 2D and 3D cases for spatially dependent non-Abelian gauge potentials, whose components are related by zero-curvature conditions suppl .

For $\Omega=0$ Eq. (1) is gauge equivalent to the Manakov model  Manakov , hence it is exactly integrable. Indeed, for $\Omega=0$, Eq. (1) is obtained from the compatibility condition $U_{t}-V_{x}+[U,V]=0$ of the eigenvalue problems:

[TABLE]

and \mbox{\boldmath\varphi\unboldmath}_{z}=V\mbox{\boldmath\varphi\unboldmath}-i\lambda^{2}\mbox{\boldmath\varphi\unboldmath}{\cal E}, with $V=i\lambda^{2}{\cal E}+i\lambda{\cal U}+\frac{1}{2}{\cal E}{\cal U}_{x}-\frac{i}{2}{\cal E}{\cal U}^{2}-\frac{i}{2}{\cal E}[{\cal A},{\cal U}]$, where $\varphi$ is a $3$ matrix, $\lambda$ is the spectral parameter, ${\cal E}=$diag$(1,1,-1)$,

[TABLE]

When all $a_{j}=0$ we recover the $UV-$representation of the Manakov model Doctor .

Turning now to the opposite limit of large Zeeman splitting and performing the rotation ${\bm{\Psi}}=e^{-i\Omega\sigma_{3}t/2}{\bm{\psi}}$ KarKonZez14 , the GPE (1) for the spinor ${\bm{\Psi}}$ preserves its original form, but now without Zeeman field and with time-dependent gauge potential ${\rm A}_{\Omega}(t)=a_{3}\sigma_{3}+e^{i\Omega\sigma_{3}t}(a_{1}\sigma_{1}+a_{2}\sigma_{2})$. At $\Omega\to\infty$, the last components of ${\rm A}_{\Omega}(t)$ become rapidly oscillating and their average effect on the dynamics vanishes. This corresponds to the rotating wave approximation with SOC being a perturbation with respect to the Zeeman field. In this case ${\rm A}_{\Omega}(t)\to a_{3}\sigma_{3}$ and thus the model again becomes exactly integrable (although different from the limit of zero Zeeman splitting). This limit can be also viewed as the nonlinear analog of Paschen-Back effect Paschen-Back , which for an atom with a random coupling (although not having a gauge structure) was discussed in MoSherKon17 .

At $\Omega>0$, the system is not integrable, but a wavepacket obeys the Ehrenfest theorem Ehrenfest

[TABLE]

which is written in terms of the soliton center of mass position $X=\|\bm{\psi}\|^{-2}\int_{-\infty}^{\infty}\bm{\psi}^{\dagger}x\bm{\psi}dx$, where the norm $\|\bm{\psi}\|^{2}=\int_{-\infty}^{\infty}\bm{\psi}^{\dagger}\bm{\psi}\,dx$ is a conserved quantity, and of the integral momentum of the soliton $\Pi=\|\bm{\psi}\|^{-2}\int_{-\infty}^{\infty}\bm{\psi}^{\dagger}{\rm P}(x)\bm{\psi}\,dx$ which is a conserved quantity in both the integrable limits discussed above.

To explore the crossover between the integrable limits we consider a soliton in a BEC with SOC of the form ${\rm A}(x)=\alpha\sigma_{1}e^{i\sigma_{3}\theta(x)}$, where $\theta(x)$ is a random function. The experimental feasibility of the model stems from different scales of the wavelength of the laser beams producing the SOC (typically below one micron), and of the random potential variations, which is about 10 $\mu$m for a 1D condensate of a transverse width of a few microns. The random field can be produced by spatially modulated beams, as shown in suppl for an example of a tripod scheme Ruseckas . Use of monochromatic quasi-nondiffracting beams nondiff1 allows for designing practically arbitrary spatial modulations nondiff5 on basis of algorithms developed in nondiff3 . For alternative possibilities of producing prescribed gauge potentials see e.g. GaugeReview .

For the numerical simulations we choose the random function $\theta(x)=2\pi f(x)/f_{\textrm{max}}\in[-2\pi,2\pi]$, where $f(x)=\sum_{j=-n}^{n}r_{j}e^{-(x-j-r_{j})^{2}/2}$, with $r_{j}\in[-0.5,0.5]$ being a uniformly random distribution with zero average value, $\langle\theta(x)\rangle=0$ (angular brackets stand for statistical averaging). The initial condition in all simulations was chosen in the form of a wavepacket with only lower state populated $\bm{\psi}_{0}(x)=e^{ivx}\mbox{sech}(x)(0,1)^{T}$, which corresponds to $\|\bm{\psi}\|^{2}=2$, $X(0)=0$, $\Pi(0)=v$, and $F(0)=0$. The evolution of such state was obtained by solving Eq. (1) for long times (up to $t_{f}=10^{3}$) for each realization of $\theta(x)$, and the subsequent averaging was performed over $10^{3}$ realizations of $\theta(x)$.

The evolution of the averaged atomic density of the dominant $\psi_{2}$ component is illustrated in Fig. 1. For each realization of the random function $\theta(x)$ the excited soliton moves as a localized object that does not spread, i.e. the $\psi_{1}$ component always accompanies the dominant $\psi_{2}$ component, and moves along the same trajectory in the $(x,t)$ plane. Such individual trajectories are resolvable in the averaged density distributions. The dynamics in panels (a) to (d) shows the crossover between the two integrable limits of $\Omega=0$ and $\Omega\to\infty$ (the transition to the latter limit is obvious already at $\Omega\sim 1$). A peculiarity of this system is that even in the integrable limit $\Omega=0$ we observe a small divergence of the linear trajectories from the central one, indicated by the dashed line [Fig. 1(a)]. This reflects the fact that the eigenvalue problem (3) is random even for deterministic initial conditions, i.e. solitons generated by the same initial condition ${\bm{\psi}}_{0}$ in different realizations of the gauge potential acquire randomly distributed parameters, including random velocities (concentrated in a narrow interval around $v$). Already for small $\Omega\sim 0.1$, when the integrability is lost [panels (b) and (c)], one observes a considerable scattering of solitons by inhomogeneities of the SOC landscape that in many cases may lead to the inversion of the soliton velocity. This scattering occurs due to random perturbation in the right-hand side of (2). In terms of the “Newtonian” dynamics of the soliton (5), this is the effect of the time dependent force $F(t)$ stemming from the noncommutativity of the gauge and Zeeman fields. Scattering becomes much weaker at $\Omega\sim 1$ [Fig. 1(d)].

To characterize the statistical properties of the evolution dynamics, we studied the average soliton displacement and the mean squared displacement (MSD). In all realizations of SOC landscape the integral soliton center of mass $X(t)$ practically coincides with the position of the soliton maximum $x_{\rm m}(t)$ defined through the relation $S_{0}(t)\equiv\bm{\psi}^{\dagger}(x_{\rm m},t)\bm{\psi}(x_{\rm m},t)=\max_{x}[\bm{\psi}^{\dagger}(x,t)\bm{\psi}(x,t)]$. However, definition for the MSD based on the position of maximum $d=\langle x_{\rm m}^{2}\rangle-\langle x_{\rm m}\rangle^{2}$ is much more accurate than the integral one, because it disregards radiation emitted by soliton interacting with random potential. For the above reasons, below we use averaged quantities based on the position of soliton maximum $x_{\rm m}(t)$.

The averaged displacement and MSD are shown in Fig. 2. Panel (a) shows variation with time of the averaged displacement in the crossover between the two integrable limits (the displacement first rapidly decreases at $\Omega\sim 0.1$, curve 2, but then gradually increases with grows of $\Omega$, see curves 3 and 4). The effect of the Zeeman field $\Omega$ is illustrated in Fig. 2(b), where a deep minimum appears in the displacement computed at $t_{f}=10^{3}$ obtained for two different Zeeman fields. This minimum, observed when the Zeeman field and the strength of the gauge field are of the same order, $\Omega\sim\alpha$, corresponds to a parameter range where the impact of the effective force $F(t)$ on the soliton propagation is strongest. In Fig. 2(c) we show the anomalous diffusion of the soliton (recall that the parameter $d$ characterizes the deviation of trajectories of the soliton motion from mean trajectory, i.e. in a sense this is a measure of the width of the averaged patterns from Fig. 1). In the integrable limit $\Omega=0$, the curve 1 represents an exact parabola, because now both $\langle x_{\rm m}^{2}\rangle$ and $\langle x_{\rm m}\rangle^{2}$ scale as $t^{2}$, with the coefficients of the proportionality being determined by the distribution of the discrete spectrum of the eigenvalue problem (3). Much stronger diffusion is observed in the nonintegrable limit at weak Zeeman field (curves 2 and 3). Interestingly, the anomalous diffusion becomes weaker with the increase of $\Omega$, and it may be even lower than diffusion at $\Omega=0$. This is also obvious from Fig. 1(d), where the width of the pattern becomes relatively narrow. This is the effect of the fast rotations, leading to zero effective gauge potential ${\rm A}_{\Omega}(t)\to 0$ (see above) for the chosen model of SOC. Thus in our system MSD is also nonmonotonic function of the Zeeman field: it is very small in two integrable limits and acquires maximal values in the crossover regime [Fig. 2(d)].

It follows from (5) that for sufficiently small $\Omega$ one can estimate $F(t)\sim\alpha\Omega$, i.e. by fixing a nonzero Zeeman field and increasing the SOC strength one results in stronger effect of the random gauge potential on the soliton. The decay of the force at large $\Omega$ is due to fast oscillations in the integrand in (5), corresponding to the limit of rotating wave approximation. Quantitatively this is illustrated in Fig. 3, where the average soliton displacement rapidly decreases to zero (due to increasing dispersion of the soliton trajectories) and by a sharp maximum of the MSD in the region where $\alpha\sim\Omega$ [cf. Fig. 2 (d)]

Turning to the (pseudo-)spinor characteristics we define {\bm{s}}(t)={S_{0}^{-1}(t)}\bm{\psi}^{\dagger}(x_{\rm m},t)\mbox{\boldmath\sigma\unboldmath}\bm{\psi}(x_{\rm m},t) with ${\bm{s}}(t)$ being always on the Bloch sphere: $|{\bm{s}}|=1$. The choice of the initial $\bm{\psi}_{0}(x)$ for numerical simulations corresponds to the “pure” state soliton bearing the spin: ${\bm{s}}(0)=(0,0,-1)$. Due to random time-dependence of the direction of ${\bm{s}}$, determined by the realization of the random gauge field, the ensemble-averaged $\langle s_{3}\rangle$ undergoes a relatively fast relaxation [see the example in Fig. 4 (a)], characterized by time Glazov of $\tau_{\rm s}\sim 1/\alpha^{2}\sqrt{\langle\Pi^{2}\rangle}\zeta,$ [cf. Eq. (5)] with $\zeta$ being the correlation length of the ${\rm A}(x)-$field. For the chosen model parameters $\zeta\sim 1$ and $v=\alpha=0.3,$ we obtain $\tau_{\rm s}\sim 30,$ in a good agreement with Fig. 4(a). The maximal relaxation of the initial spin is achieved in the integrable limit at zero Zeeman splitting [Fig. 4(b)]. A specific feature of the nonintegrable regime, shown in Fig. 1, is the decrease in $\sqrt{\langle\Pi^{2}\rangle}$ with time, slowing the relaxation down, as can be seen in Fig. 4(a). Such a behavior is a consequence of the “independent” deterministic dynamics of the soliton center of mass, described by ${\bm{u}}$ in (2) at $\Omega=0$ and stochastic dynamics of the soliton pseudo-spin ${\bm{s}}.$ Increasing the Zeeman field results in restoring the pure character of the soliton spin state, which is observed in Fig. 4 already at $\Omega\gtrsim 0.3$. After a short interval of growth of $\langle s_{3}\rangle$, in the interval $0.04\lesssim\Omega\lesssim 0.3$, increasing of the Zeeman splitting results in gradual decrease of $\langle s_{3}\rangle$.

To conclude, we described the evolution of solitons in inhomogeneous gauge potentials. In the absence of the Zeeman field the model is exactly integrable for arbitrary spatial distributions of the matrix gauge potential. Solitons, and more sophisticated solutions can also be constructed using the Inverse Scattering Technique. We described statistics of solitons affected by random SOC. In the integrable case solitons move with constant velocities which are different for different realization of the SOC. When the Zeeman splitting is large, the system again approaches an integrable limit, although different from the one at zero Zeeman splitting. The crossover between these two integrable limits is characterized by strong interaction of a soliton with the random gauge potential, manifesting itself in a slowing down average motion and strongly anomalous diffusion of solitons. Each soliton carries a pseudo-spin. The dynamics of ensemble-averaged pseudo-spinors is characterized by two temporal scales: the fast relaxation at initial stages, well-described in quasi-linear approximation, and the long-time slow evolution.

I On separation of nonlinear time evolution and linear field distribution in a two-dimensional case with a non-Abelian gauge potential.

The separation of ”linear” and ”nonlinear” dynamics, reported in the main text for one-dimensional (1D) Gross-Pitaevskii equations (GPE) can also be performed in the 2D and 3D cases for specific types of non-Abelian potentials. In order to illustrate this, here we consider 2D GPE

[TABLE]

where

[TABLE]

${\bm{A}}=(A_{x},A_{y})$, with $A_{x,y}$ being $2\times 2$ matrices, is an inhomogeneous non-Abelian gauge potential (it is arbitrary, so far), and the Zeeman splitting is set to zero.

Next we consider the eigenvalue problem

[TABLE]

which by the ansatz {\bm{\xi}}_{\bm{k}}({\bm{r}})=e^{i{\bm{k}}\cdot{\bm{r}}}\mbox{\boldmath\phi\unboldmath}_{1}({\bm{r}}), is reduced to {\bm{P}}\mbox{\boldmath\phi\unboldmath}_{1}=0, or explicitly

[TABLE]

The obtained equation, is not solvable for arbitrary potential ${\bm{A}}({\bm{r}})$. Indeed considering the components of (9) separately one obtains that the condition

[TABLE]

which can be viewed as the zero-curvature condition, must be satisfied. Notice, that for a constant potential, i.e., for $A_{x}=$const and $A_{y}=$const, the solvability condition requires the potential to be Abelian, i.e., to have $[A_{x},A_{y}]=0$.

Thus, we require (10) to be satisfied. Furthermore, like in the main text, we require the gauge potential to be Hermitian, i.e., $A_{x,y}^{\dagger}=A_{x,y}$ and to anti-commute with the time reversal: ${\cal T}{\bm{A}}+{\bm{A}}{\cal T}=0$. Thus the vector function defined by \mbox{\boldmath\phi\unboldmath}_{2}={\cal T}\mbox{\boldmath\phi\unboldmath}_{1}, solves {\bm{P}}\mbox{\boldmath\phi\unboldmath}_{2}=0. Furthermore, it is straightforward to verify that

[TABLE]

and

[TABLE]

Now, one can search for a solution of (1) in a form of the ansatz

[TABLE]

where $u_{1,2}({\bm{r}},t)$ are two unknown functions, which solve the equation

[TABLE]

where ${\bm{u}}=(u_{1},u_{2})^{T}$ and we used property (12).

The above consideration can be straightforwardly generalized to the 3D case, where the solvability condition for (9) requires vanishing of all nondiagonal elements of the curvature tensor.

For the sake of illustration of a (nontrivial) non-Abelian gauge potential for which the separation of nonlinear time evolution and linear spatial spinor field distribution is possible, we consider

[TABLE]

where ${\bm{a}}$ and ${\bm{b}}$ depend only on $x$. Then (10) is reduced to the system of ODEs

[TABLE]

As it is mentioned above, it follows from (10) that two arbitrary stationary (coordinate independent) potentials must be Abelian, for the separation ansatz to be applicable. In particular, the conventional two-dimensional Rashba and Dresselhaus couplings are given by coordinate-independent ${\bm{a}}=(0,1,0),\,{\bm{b}}=(-1,0,0)$ and ${\bm{a}}=(1,0,0),\,{\bm{b}}=(0,-1,0),$ correspondingly, with $|{\bm{a}}\times{\bm{b}}|=1$. Therefore, these interactions cannot be gauged out Aleiner . However, if either ${\bm{a}}=\mathbf{0}$, or ${\bm{b}}=\mathbf{0}$, or the SOC contains the Rashba and Dresselhaus contributions of equal strengths, then ${\bm{a}}\times{\bm{b}}=\mathbf{0},$ and the required gauge transformation is possible Tokatly .

If Eq. (16) is satisfied by given (either coordinate-dependent or independent) vectors ${\bm{a}}$ and ${\bm{b}}$, then for the vector potential given by (15) we can look for a solution of (9) in the form

[TABLE]

where

[TABLE]

(these equations are consistent) and $q$ is a constant. Finally, looking for solutions of (2) independent on $y$ we arrive at the 1D Manakov model considered in the main text.

II A scheme for random gauge potentials

In order to describe a possibility of how a random gauge potential can be generated, let us consider a BEC with four-level atoms in a tripod configuration described by the Hamiltonian Ruseckas

[TABLE]

where $|j\rangle$ ($j=1,2,3)$ are the low-energy states and $|0\rangle$ is the excited state which is coupled to the states $|j\rangle$ by the Rabi frequencies $\Omega_{j}$. Consider now a BEC loaded in a cigar-shaped trap, which is long enough along the $x$-direction (say, approximately 200$\,\mu$m length) and has transverse radial width in the $(y,z)-$plane of the order of $a_{\bot}\approx 10\,\mu$m. The coupling of the low-energy states with the excited state is assured by the two counter-propagating laser beams:

[TABLE]

where ${\bm{r}}=(x,y,z)$, $\vartheta$ is a real constant and $\Omega$ is the field amplitude. The beams $\Omega_{1,2}$ propagating along the directions $(\pm\cos\varphi,\sin\varphi,0)$, where $\cos\varphi=k_{x}/\tilde{k}$ and $\sin\varphi=k_{y}/\tilde{k}$ (i.e. $\tilde{k}^{2}=k_{x}^{2}+k_{y}^{2}$), can be created as superposition of nondiffracting beams nondiff1 . For example, one can represent

[TABLE]

where $(\xi,\eta,\zeta)$ are the Euclidian coordinates in the rotated frame, $\xi=x\cos\varphi+y\sin\varphi$, $\eta$ and $\zeta$ are the coordinates in the plane orthogonal to $\xi-$axis, and $k=\omega/c$. In Eq. (24) the angular spectrum ${\cal A}(\nu,k_{\bot})$ is defined in the Fourier domain, on a narrow annular ring of the width $2\delta k$ having central radius $k_{0}=\sqrt{k^{2}-\tilde{k}^{2}}$ ($\nu$ is the angular coordinate).

Engineering of the spectrum ${\cal A}(\nu,k_{\bot})$ using iterative Fourier methods, reminiscent of the methods employed in phase retrieval and image processing algorithms nondiff3 , allows researchers to produce quasi-nondiffracting monochromatic light patterns with any desired phase or intensity distribution in the $(\eta,\zeta)$ plane and characteristic features with scales $\sim 2\pi/k_{0}$ ranging from several to hundreds of microns, as demonstrated in nondiff5 .

Now, the two dark states of $H_{0}$ can be found in the form:

[TABLE]

The spinor wave-function is sought in the form

[TABLE]

Now, in the absence of interactions the evolution of the spinor $\bm{\psi}=\left(\Psi_{1},\Psi_{2}\right)^{T}$ is governed by the Hamiltonian Ruseckas :

[TABLE]

where $M$ is the atomic mass and the vector matrix $A$ (known also as Berry connection) has elements

[TABLE]

The total potential $V_{\rm tot}({\bm{r}})$ consists of two parts: one is the external trap potential $V_{\rm ext}$ which is a matrix if the components are coupled, while another part $U$ is a matrix potential induced by the laser beams (23). It has components Juzel

[TABLE]

Due to quasi-one-dimensionality of the condensate, we are interested only in the distribution of $\theta(x)=\Theta(\hat{{\bf i}}x)$, i.e. in the distribution of $\Theta({\bm{r}})$ along the $x$-axis (at $y=z=0$). The only requirement for the function $\Theta({\bm{r}})$, used so far in (24), is that it must be slowly varying on the scale of the wavelength of the beams $\Omega_{1,2}$, i.e. on the scale $\lambda=2\pi/k$.

Substitution of the dark states (25) and (26) in this formula yields the $x-$component of the dimensionless gauge potential

[TABLE]

i.e. the formula used in the text. Here we neglected the derivative of slowly varying $\theta(x)$.

Considering the matrix $U$ in the same approximation of slowly varying $\theta(x)$, one obtains that this potential is diagonal:

[TABLE]

Thus it can be compensated by the respective constant external potentials for the spinor components.

Including inter-atomic interaction, averaging over the cross-section of the trap in the ($y,z$) plane, and rescaling variables such that the longitudinal coordinate is measured in the units of $a_{\bot}$, while the energy is measured in the units of $\hbar\omega_{\bot}$ (where $\omega_{\bot}$ is the linear harmonic oscillator frequency of the parabolic trap in the transverse direction), one ends up with equation (1) from the main text.

As the final step we take into account that the experimental length scale values of the coupling field are typically hundreds of nanometers. On the other hand, typical transverse scale of the trap is of several microns, while its length can be of a few hundreds of microns. Thus, the suggested beam configuration can create almost arbitrary, in particular random, potentials $\theta(x)$ using monochromatic beams, as describes above.