Magnetic soliton molecules in binary condensates

TL;DR

This paper investigates magnetic soliton molecules in binary Bose-Einstein condensates, analyzing their formation, stability, and collision dynamics through analytical and numerical methods, and suggests experimental realization is feasible.

Contribution

It provides the first detailed analysis of bound states of magnetic solitons, including their interaction potential, dissociation energy, and collision behaviors in binary condensates.

Findings

Bound states of magnetic solitons can form stable molecules.

Collision dynamics include dipole-like interactions between bound states.

Analytical results agree well with numerical simulations.

Abstract

Two-component Bose-Einstein condensates in the miscible phase can support polarization solitary waves, known as magnetic solitons. By calculating the interaction potential between two magnetic solitons, we elucidate the mechanisms and conditions for the formation of bound states -- or molecules -- and support these predictions with dynamical simulations. We analytically determine the dissociation energy of bound states consisting of two oppositely polarized solitons and find good agreement with full numerical simulations. Collisions between bound states -- either with other bound states or with individual solitons -- produce intriguing dynamics. Notably, collisions between a pair of bound states exhibit a dipole-like behavior. We anticipate that such bound states, along with their rich collision dynamics, are within reach of current experimental capabilities.

Equations

$${\Psi }_i(\mathbf{x})={\psi }_i(x)\frac{\exp [-(y^2+z^2)/2l^2]}{l\sqrt{\pi }},$$

$$i\hslash {\partial }_t{\psi }_1$$

$$i\hslash {\partial }_t{\psi }_2$$

$${\partial }_{\tau }{\varphi }_A$$

$${\partial }_{\tau }\theta$$

$$E_{\mathrm{SE}}=\frac{n\hslash c_s}{4}\int d\zeta \left[{\left(\frac{\partial \theta }{\partial \zeta }\right)}^2+{\sin }^2\theta {\left(\frac{\partial {\varphi }_A}{\partial \zeta }\right)}^2+{\cos }^2\theta \right].$$

$$\cos \theta$$

$${\varphi }_A$$

$$\cos \theta (\zeta ,\tau =0)$$

$${\varphi }_A(\zeta ,\tau =0)$$

$$E_{\text{MS}}^{\text{uni}}(U)=n\hslash c_s\sqrt{1-U^2},$$

$$E_{\text{RPDW}}^{\text{uni}}(b)=2n\hslash c_s/b.$$

$$U=\pm \sqrt{1-\frac{1}{b^2}}.$$

$$m_{\text{eff}}=\frac{1}{U{c}_s^2}\frac{d{E}_{\text{MS}}^{\text{uni}}}{dU}=-\frac{n\hslash }{c_s}\frac{1}{\sqrt{1-U^2}}.$$

$$V=\text{sign}(m_{\text{eff}})\left[E_{\text{2MS}}(\Delta x)-E_0\right],$$

$$E^{\text{uni}}=E_{\text{RPDW}}^{\text{uni}}(b)-2E_{\text{MS}}^{\text{uni}}(U=0)=2n\hslash c_s\left(\frac{1}{b}-1\right).$$

$$E=E_{\mathrm{RPDW}}(b)-E_0.$$

$$E_{\text{MS}}=\text{sign}(m_{\text{eff}})\left(E_0^{2L}/2-E_0\right),$$

Full text

I Introduction

Stabilized by a balance between dispersion and nonlinearity, solitons are localized collective excitations in elongated systems, relevant for a broad range of systems ranging from water in narrow channels [1, 2] to optical fibers [3], plasmas [4] and conducting polymers [5, 6]. Moreover, ultracold Bose-Einstein condensates (BECs) provide an especially appealing platform for the investigation of solitons, due to the clean and highly controllable experimental apparatuses available. As classic examples, BECs with defocusing nonlinearities can support dark solitons [7, 8, 9, 10], while bright solitons may form in the presence of focusing interactions [11, 12].

Multi-component BECs allow one to investigate the rich physics offered by vector solitons, non-dispersive waves spanning multiple coupled components 111It is worth noting that, although vector solitons may not always be true solitons in the strict mathematical sense, they possess key soliton characteristics, and we will refer to them as solitons for consistency with the literature.. For the case of two distinguishable components [14, 15, 16, 17, 18, 17]—a binary condensate—experiments have observed dark-bright [19] and dark-antidark solitons [20, 21, 22, 23, 24], as well as the closely related magnetic solitons [25, 26] (see related theory [27, 28]), which occur in balanced BEC mixtures close to the miscible-immiscible transition.

Soliton bound states—composite structures of multiple solitons—were predicted by Öhberg and Santos using pairs of magnetic solitons of opposite polarization [29]222In Ref. [29], magnetic solitons are referred to as kink-antikink solitons. (also see more recent works [31, 32]). Closely related bound structures were later created in experiments by accelerating two overlapping condensates in opposite directions to create counterflows [20, 21]. Despite being predicted relatively early, much remains unknown about magnetic soliton bound states and their underlying physical mechanisms.

In this paper, we theoretically uncover the properties and underlying mechanisms of bound states formed by magnetic solitons in two-component BECs. To build intuition, Figs. 1(a) and (b) show the density and phase, respectively, of a pair of at-rest solitons with opposite polarization. Since each soliton has zero velocity, its core exhibits a $\pi$ phase jump in one component, while the other remains phase-flat. Notice also that magnetic solitons manifest predominantly in the spin density $s(x)=n_{1}(x)-n_{2}(x)$, whereas the total density $n(x)=n_{1}(x)+n_{2}(x)$ remains nearly constant. When this pair is evolved in time, Fig. 1(c) shows that the constituent solitons attract and pass through one another, forming an oscillating bound state, whereas solitons with like polarization repel each other [Fig. 1(d)]. These contrasting behaviors motivate our systematic study of soliton bound states. We numerically calculate the interaction potential in both cases, elucidating the mechanism for bound state formation. We further show that the potential flattens anharmonically with increasing separation, consistent with the diverging oscillation period at larger excitation amplitudes. We also derive the bound state dissociation energy analytically and confirm quantitative agreement with dynamical simulations. Finally, we investigate collisions involving soliton bound states, revealing dynamics that range from nearly elastic bounces to soliton swapping, analogous to substitution reactions.

II Formalism

II.1 System and methods

We consider two coupled BECs trapped in an elongated geometry, where both components experience the same harmonic potential along the $y$ and $z$ directions, with isotropic trapping frequencies $\omega_{y}=\omega_{z}\equiv\omega$, and are unconfined along $x$. The wavefunction of component $i$ is described by the separable ansatz,

[TABLE]

with characteristic length $l=\sqrt{\hbar/m\omega}$, where the mass $m$ is assumed to be the same for both components. Ansatz (1) is substituted into the two-component Gross-Pitaevskii equation (GPE) [33]. The radial part is then integrated out, and subtracting the radial zero-point energy yields the quasi-one-dimensional (quasi-1D) GPEs:

[TABLE]

with linear densities $n_{j}=|\psi_{j}(x)|^{2}$ and interaction strengths $g_{ij}=4\pi a_{ij}\hbar^{2}/m$, where $a_{ij}$ is the $s$-wave scattering length between components $i$ and $j$, and $g_{12}=g_{21}$.

Equation (2) is solved with a 4th-order Runge-Kutta method, and stationary solutions are found using imaginary time evolution. To create zero-velocity solitons, we enforce the phase of the order parameter to exhibit a $\pi$ step function for the appropriate component at the desired position for every imaginary time step, similar to the approach used in Refs. [34, 35] to create solitons in one-component systems. We validate this method by evolving the converged solutions in real time. For both isolated at-rest solitons and bound-state pairs produced at stationary points of their inter-soliton potential, we observe no change during the evolution, confirming their stationary state character. This method is also useful for creating zero-velocity soliton pairs that are not necessarily stationary states. For example, the initial state ($t=0$) in Fig. 1(a) and 1(a) was created in this way.

II.2 Parameters and considerations for experiments

For the purpose of realizing magnetic solitons, we consider the case of a balanced mixture, i.e., $g_{11}=g_{22}\equiv g$, $N_{1}=N_{2}\equiv N/2$ and average 1D density per component $n_{0}=N/2L$. We choose a system size $L=30\,\upmu$m and spin-healing length $\xi_{s}=0.4703\,\upmu$m, where $\xi_{s}=\hbar/\sqrt{2mn_{0}\delta g/\pi l^{2}}$, with $\delta g=g-g_{12}$.

While our results are given in dimensionless units, they can easily be translated into SI units. For example, using the mass of 87Rb, our computational parameters can correspond to $n_{0}=50\,\upmu\rm{m}^{-1}$, $\omega=2\pi\times 4969\,$Hz, and $(a_{11},\,a_{12},\,a_{22})=(100,95,100)\,a_{0}$, where $a_{0}$ is the Bohr radius. Alternatively, with the choice of 23Na, our computational parameters can correspond to $n_{0}=200\,\upmu\rm{m}^{-1}$, $\omega=2\pi\times 8611\,$Hz, and $(a_{11},\,a_{12},\,a_{22})=(54.54,51.81,54.54)\,a_{0}$. These interaction combinations are close to the regimes already realized in magnetic soliton experiments, e.g., see [20, 25], and we have checked that our results are qualitatively representative of a broad range of density and trapping strength combinations.

III Formulation of the bound state dissociation criterion

The order parameter of component $i$ can be parametrized as $\psi_{i}(x)=\sqrt{n_{i}(x)}e^{i\phi_{i}(x)}$, where $\phi_{i}$ is the phase. For spin excitations and dynamics in two-component BECs with $\delta g\ll g$, one can safely assume that the total density $n$ is uniform (uniform-$n$ approximation) and the dynamics are solely characterized by the spin density $s(x)$ and relative phase $\phi_{A}(x)=\phi_{1}(x)-\phi_{2}(x)$ [27]. Introducing the polarization distribution, $\cos\theta(x)=s(x)/n$, the GPEs can be reformulated into the following coupled differential equations [28]

[TABLE]

where $\zeta=x/\xi_{s}$, $\tau=c_{s}t/\xi_{s}$ are the dimensionless coordinate and time, respectively, with $c_{s}=\sqrt{n_{0}\delta g/2m\pi l^{2}}$ the speed of spin sound. This formalism allows the energy of an excitation in the spin density to be evaluated with [28]

[TABLE]

Prior work [27] revealed a traveling-wave solution, the so-called magnetic soliton, of the form $\theta(\zeta,\tau)\equiv\theta(\zeta-U\tau)$ and $\phi_{A}(\zeta,\tau)\equiv\phi_{A}(\zeta-U\tau)$ where $U=v/c_{s}$ is the speed of the soliton. The static magnetic soliton $(U=0)$ takes the following simple form

[TABLE]

where $\Theta(\zeta)$ is the Heaviside step function. We define the $(+)$ solution as spin-up magnetic soliton, since $s(x)>0$ around the core [see soliton on the right in Fig. 1(a)]. Conversely, the $(-)$ solution is referred to as a spin-down soliton. Note that, for $U=0$, the antidark component ($n_{i}(x)>n_{0}$ around the core) has no phase variation and the dark component ($n_{i}(x)<n_{0}$) exhibits a $\pi$ phase jump.

An intriguing fact is that a pair of oppositely polarized magnetic solitons, initially cospatial, can be generated with opposite velocities by imprinting a relative-phase domain wall of the form [36]:

[TABLE]

where $b$ is its characteristic width. Although the spin density is spatially uniform, $\phi_{A}$ can be symmetrically distributed between the two components, revealing the soliton pairs shown in Fig. 2, where three different widths are considered. The $\tau>0$ dynamics of the pair falls into one of two categories, as can be predicted by the following energy consideration. Using Eq. (4), we can obtain the energy of a single magnetic soliton

[TABLE]

and the energy of the relative-phase domain wall

[TABLE]

On the one hand, for large widths $b>1$ (small energies), the energy of the initially imprinted relative-phase domain wall can be equal to the sum of the energies of two widely separated solitons moving in opposite directions with $U$ given by

[TABLE]

On the other hand, when the initial width is narrow, $b<1$, the energy of the generated soliton pair is too large to form well-separated solitons. Instead, they quickly evolve into a pair of static magnetic solitons ($U=0$) at finite separation, before again moving towards one another, thus forming an oscillatory bound state (analogous to an excited molecule). This counterintuitive behavior may be understood by noting that the energy of a soliton decreases with increasing velocity [see Eq. (7)], which corresponds to a negative effective mass 333For comparison, it is interesting to note that while dark solitons also have negative effective masses [39], certain dark-bright solitons have been predicted to dynamically oscillate between negative and positive effective mass [40].

[TABLE]

The boundary between the two regimes is set by the critical width $b_{\text{c}}=1$, which corresponds to the critical bound-state dissociation energy $E_{\text{c}}=2n\hbar c_{s}$.

IV Numerical results

IV.1 Soliton-soliton interaction potential

To shed light on the dynamics of magnetic soliton pairs, we define the intersoliton potential as

[TABLE]

where $E_{\text{2MS}}$ is the total system energy for two solitons at rest separated by $\Delta x=x_{2}-x_{1}$, with $x_{i}$ denoting the position of soliton $i$. The reference energy $E_{0}$ corresponds to the total system energy for maximally separated solitons at rest, so that $V\to 0$ when the solitons are far apart. The sign prefactor in Eq. (11) is included to make the relationship between the force and acceleration more intuitive. Since magnetic solitons have a negative effective mass, $\text{sign}(m_{\text{eff}})=-1$ [see Eq. (10)], the acceleration $a=F/m_{\text{eff}}$ is actually opposite in direction to the force $F$. However, this sign difference is absorbed into our definition of $V$ and bound states thus occur within potential minima.

Figure 3(a) shows $V$ for two solitons with opposite polarizations. The blue points are calculated using the GPEs following the approach outlined in Sec. II.1. Importantly, there is a minimum at $\Delta x=0$, indicating that the solitons experience an attraction toward zero separation. This behavior is qualitatively different from that of real molecular potentials, where the minima occur at finite separation. We have confirmed that zero-separation bound states remain stationary under real time evolution of the GPEs (not shown). In contrast, Fig. 3(b) shows $V$ for a pair of solitons with like polarization. In the limit $\Delta x\to 0$, $V$ approaches the energy of a single magnetic soliton in the system (horizontal dashed line; see Appendix A for details). To understand this, recall that the wavefunction phase has two $\pi$ jumps in the same component. Taking the limit $\Delta x\to 0$ but excluding $\Delta x=0$, the region where the density is forced to zero approaches a single point, just as for a single soliton. For increasing $|\Delta x|$, $V$ decreases monotonically, indicating that like-sign solitons repel each other and thus cannot form bound states [see Fig. 1(d)].

A semianalytic prediction for the intersoliton potential can be made based on the uniform-$n$ approximation outlined in Sec. III. The soliton interaction energy is obtained from the difference between the energy of a relative-phase domain wall of width $b<1$ and that of two isolated magnetic solitons at zero velocity, i.e.,

[TABLE]

The energy of the relative phase domain wall is equivalent to that of two solitons at rest separated by $\Delta x$, provided that $b<1$, because these two configurations correspond to different phases of the same oscillation. To determine the relation between $b$ and $\Delta x$, we evolve relative phase domain walls under the uniform-$n$ equations (3) and extract $\Delta x$ as the maximum soliton separation at rest after a quarter of an oscillation period. The resulting semianalytic prediction of the potential $V^{\text{uni}}=\text{sign}(m_{\text{eff}})E^{\text{uni}}$ as a function of $\Delta x$ is plotted in Fig. 3(a) as solid lines, demonstrating good agreement with the GPE results for moderate soliton separations. Note, however, that $V^{\text{uni}}$ diverges for $\Delta x\ll\xi_{s}$, because the uniform-$n$ constraint significantly overestimates the kinetic energy when $b\ll 1$.

The oscillation period $T$ of excited bound states, shown in the inset to Fig. 3(a), can serve as a useful experimental probe. We calculate $T$ by evolving the at-rest solitons from the main plot of Fig. 3(a) as initial states in time-dependent GPE simulations. The period exhibits anharmonicity, diverging as the oscillation amplitude $\Delta x_{\text{max}}$ increases and the bound state dissociation energy is approached. The large-amplitude case highlighted by the magenta circles corresponds to the dynamics in Fig. 1(c). For reference, taking the mass of 87Rb (see Sec. II.2) yields a period of $T\approx 60.6\,$ms for this case.

IV.2 Critical binding energy

A relative-phase domain wall is equivalent to a pair of counterpropagating, oppositely polarized solitons. The soliton speed increases as $1/b$ decreases, while the energy decreases. In Fig. 4, we exploit this by performing dynamical GPE simulations to determine the critical energy below which a soliton pair becomes unbound.

Relative-phase domain walls of width $b_{0}$ [Eq. (6)] are first imprinted on otherwise uniform BECs. These are then subjected to a short imaginary-time evolution ($\Delta t=312i/\omega$) to relax the density near the soliton cores, which need not be uniform. The resulting states, which serve as our initial conditions for real-time evolution, differ slightly in width from $b_{0}$; we extract the actual $b$ by fitting their phase profiles [Eq. (6)]. The fit errors are small, with error bars in Fig. 4 generally smaller than the symbols (see Appendix B). The system energy $E_{\rm RPDW}$ is calculated for each state, with the reference energy of two maximally separated solitons subtracted to obtain

[TABLE]

The states are then evolved in real time to determine whether they remain bound or become unbound.

Figure 4 shows the energy of relative phase domain walls [Eq. (13)] as a function of inverse width. Squares (circles) denote soliton pairs that evolve to be unbound (bound) in real-time simulations. Two example simulations are plotted as insets, one on each side of the binding threshold, i.e., $b^{-1}=0.9$ and $b^{-1}=1.1$. The analytic uniform-$n$ prediction [Eq. (12)] agrees quantitatively with the GPE simulations [Eq. (13)] near the threshold. At larger $1/b$ (slow-moving solitons), deviations arise from the constraint of uniform total density (cf. Sec. IV.1).

IV.3 Collisions involving bound states

A well-known property of solitons is their ability to survive collisions unscathed. Because this work does not focus on the Manakov limit, which requires $\Delta g=0$, our system is not perfectly integrable, and we therefore do not expect perfectly elastic collisions. Nevertheless, it is of interest to test whether magnetic soliton bound states can remain intact after collisions. In realistic experimental settings, isolated bound states are likely to occur in excited configurations. Accordingly, in Fig. 5 we investigate collisions involving excited bound states under various conditions.

Since solitons of like polarization repel and those of opposite polarization attract, an excited bound state behaves like a dipole oscillating with period $T$. The time-dependent polarity of these dipoles is illustrated in Fig. 5 by black arrows pointing from spin-down (red) to spin-up (blue) solitons within each bound state. Figure 5(a) shows a collision between two tightly bound states, each with an initial $T=31.04\,$ms. In this case, the polarities of the two bound states remain opposite throughout the collision, leading to inter bound-state repulsion and a clear bounce, with both bound states remaining intact. Evidence for bound state energy loss is seen in the slightly increasing oscillation periods over time. The energy lost is converted to phonons, which are not visible with this color scheme.

Figure 5(b) shows two tightly bound states (initially $T=34.2\,$ms) with the same polarity (black arrows pointing in the same direction). In this case, the bound-state dipoles tend to attract each other, yet both bound states nevertheless survive the collision. This attractive interaction is only weakly visible due to the tight binding and the relatively fast center-of-mass motions of the bound states. By contrast, the attractive interaction between same-polarity bound states has striking consequences for two loosely bound states (initially $T=128.9\,$ms), as shown in Fig. 5(c). Here, the inner solitons—having opposite polarization—not only cross each other but are actually torn from their respective bound states and form a new bound state, while the two outer solitons become unbound. Finally, in the case of a single soliton colliding with a bound state (initially $T=36.59\,$ms), Fig. 5(d) illustrates a type of substitution reaction, in which the newcomer replaces one of the bound state’s solitons and ejects it in the process.

V Discussion

In this paper, we investigated the mechanisms underlying the formation of magnetic soliton bound states and their collision dynamics. By calculating the inter-soliton potential, we showed that magnetic solitons with opposite polarization can form bound molecular states, whereas those with like polarization cannot. The analytically determined dissociation energy of such bound states shows excellent agreement with numerical simulations. The dynamics of bound states exhibit a diverging oscillation period as the unbinding threshold is approached. Collisions between excited bound states exhibit dipole-like interactions, leading to attraction or repulsion depending on the relative dipole polarity. We also demonstrated substitution-type reactions, in which one soliton of a bound state is replaced by another that collides with it.

Future studies could investigate bound states in binary dipolar BECs, where long-range interactions and the presence of a roton-maxon spectrum may qualitatively affect the inter-soliton potentials and dynamics. Another interesting direction is the generalization to Rabi-coupled systems, which exhibit rich magnetic-soliton dynamics. Furthermore, one could study the crossover from the quasi-1D to the three-dimensional cigar-shaped regime. While magnetic solitons are expected to be more robust against three-dimensional instabilities than their single-component counterparts—since their size is governed by the spin-healing length, which grows large near the immiscibility threshold—it would be valuable to investigate the two-component analog of snake instabilities and the formation of vortex rings. More generally, dark-antidark solitons also occur in imbalanced binary BECs, where the intraspecies interactions of the components differ [38]. We have confirmed that bound states also arise in such imbalanced systems (up to at least $g_{11}/g_{22}=1.2$; not shown); a systematic study of these bound states represents another interesting research direction. From an experimental perspective, oscillations of bound states, as well as collisions involving magnetic solitons and their molecular counterparts, could provide direct information about the interaction potentials, offering a route to test the predictions of this work.

Acknowledgements

We thank T. Billam, T. Bland, F. Ferlaino, L. Bellinato Giacomelli, A. Madhusudan, N. Masalaeva, N. Parker, E. Poli, and P. Senarath Yapa for fruitful discussions. This research was funded in part by the Austrian Science Fund (FWF) [Grant DOI: 10.55776/P 36850-N]. C. Q. is supported by ACC New Jersey under Contract No. W15QKN-18-D-0040.

Appendix A Numerical calculation of single-soliton energy

Our numerical GPE simulations do not assume a constant total density $n$ and are performed in a finite-sized system. For these reasons, rather than using Eq. (7), we numerically calculate the energy of a single soliton in our system. Due to the periodic boundary conditions, our imprinting protocol does not permit the creation of a single isolated soliton. Instead, we consider a system of length $2L$ containing two maximally separated solitons (with energy $E_{0}^{2L}$), which has the same number of solitons per unit length as a system of size $L$ with only one soliton. The numerically calculated energy associated with a single soliton is then defined as

[TABLE]

where $E_{0}$ is the same reference energy used in Eq. (11). Equation (14) is plotted as the horizontal dashed line in Fig. 3(b). Note that the difference between Eq. (14) and Eq. (7) is approximately 7$\%$ for $U=0$.

Alternatively, recall that for like-sign solitons in the limit $\Delta x\to 0$, the inter-soliton potential $V$ [Eq. (11)] approaches the energy of a single magnetic soliton. In Fig.3(b), a single data point represents $V(\Delta x\to 0)$. This value cannot be obtained using realistic numerical grids. Instead, we employ a numerical workaround: during each step of imaginary-time evolution in a system of length $L$, the density is artificially constrained to zero for only a single grid point. This procedure yields excellent agreement with the single-soliton energy $E_{\text{MS}}$.

Appendix B Fitting of relative phase domain walls to find their widths

In Sec. IV.2, we extract the relative phase domain wall width $b$ by fitting the relative phase $\phi_{A}$ of the initial states to Eq. (6). The fit is performed using three different fitting widths: $x_{b}-x_{a}$, $4(x_{b}-x_{a})$, and $12(x_{b}-x_{a})$, with $x_{a}$ and $x_{b}$ being the positions where $\phi_{A}(x_{a})=0.6\pi$ and $\phi_{A}(x_{b})=1.4\pi$, respectively (see Fig. 6). The final value of $b$ is taken as the mean of the three fits, and the error is given by the standard deviation. An example is shown in Fig. 6.