Nonlinear Waves in an Experimentally Motivated Ring-shaped Bose-Einstein Condensate Setup

TL;DR

This paper constructs and analyzes various stationary soliton states in a ring-shaped Bose-Einstein condensate, exploring their stability, dynamics, and experimental creation methods within a nonlinear Schrödinger framework.

Contribution

It systematically identifies and characterizes a wide range of soliton and vortex states in a ring-shaped BEC using bifurcation analysis and stability studies.

Findings

Multiple stable and unstable soliton configurations identified.

Unstable states can lead to persistent vortical dynamics.

Phase-imprinting can create these states in experiments.

Abstract

We systematically construct stationary soliton states in a one-component, two-dimensional, repulsive, Gross-Pitaevskii equation with a ring-shaped target-like trap similar to the potential used to confine a Bose-Einstein condensate in a recent experiment [Eckel, et al.\ {\em Nature} {\bf 506}, 200 (2014)]. In addition to the ground state configuration, we identify a wide variety of excited states involving phase jumps (and associated dark solitons) inside the ring. These configurations are obtained from a systematic bifurcation analysis starting from the linear, small atom density, limit. We study the stability, and when unstable, the dynamics of the most basic configurations. Often these lead to vortical dynamics inside the ring persisting over long time scales in our numerical experiments. To illustrate the relevance of the identified states, we showcase how such dark-soliton configurations (even the unstable ones) can be created in laboratory condensates by using phase-imprinting techniques.

Equations

$$\displaystyle i\hbar\partial_{t}\Phi=\Big{[}-\frac{\hbar^{2}}{2m}\nabla^{2}+g_{\rm 2D}|\Phi|^{2}+V(\boldsymbol{r})\Big{]}\Phi\,,$$

$$V=\frac{1}{2}m\left({\omega }_x^2{x}^2+{\omega }_y^2{y}^2\right),$$

$$i{\partial }_t\Phi =-\frac{1}{2}{\nabla }^2\Phi +|\Phi {|}^2\Phi +V(\mathbf{r})\Phi ,$$

$$V(r)=\{\begin{array}{ll}1-A{e}^{-\frac{(r-r_{\text{ring}}{)}^2}{s_{\text{ring}}^2}}-e^{-\frac{(r-r_{\text{disk}}{)}^2}{s_{\text{disk}}^2}}& r\ge r_{\text{disk}}\\ -A{e}^{-\frac{(r-r_{\text{ring}}{)}^2}{s_{\text{ring}}^2}}& r<r_{\text{disk}},\end{array}$$

$$\Phi (\mathbf{r},t)=\varphi (\mathbf{r})e^{-i\mu t},$$

$$-\frac{1}{2}{\nabla }^2\varphi +|\varphi {|}^2\varphi +\left[V(x,y)-\mu \right]\varphi =0.$$

$$N={\iint }_{-\infty }^{+\infty }{\left|\varphi (x,y)\right|}^{2}dxdy.$$

$$\varphi (\mathbf{r})={\varphi }_0(\mathbf{r})+\left[a(\mathbf{r})e^{i\omega t}+b^{\star }(\mathbf{r})e^{-i{\omega }^{\star }t}\right],$$

$$-\omega \left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{cc}{A}_{11}& A_{12}\\ -A_{12}^{\star }& -A_{11}\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right),$$

$$A_{11}$$

$$A_{12}$$

Full text

I Introduction

Atomic Bose-Einstein condensates (BECs) Pethick and Smith (2008); Pitaevskii and Stringari (2003); Kevrekidis et al. (2008); Bagnato et al. (2015) offer an ideal testing ground for confronting theoretical models of nonlinear matter waves with experimental data. Since their experimental realization, there have been tremendous advances Bagnato et al. (2015); Cornell and Wieman (2002); Dalfovo et al. (1999); Ketterle (2002); Leanhardt et al. (2002) in trapping, guiding, manipulating and controlling BECs. For instance, recent advances in all-optical trapping Henderson et al. (2009); Gaunt and Hadzibabic (2012); Pasienski and DeMarco (2008) have produced confined atomic clouds with temperatures at the nanokelvin scale. All-optical trapping, in turn, has enabled the strength of the atom-atom interactions in atomic gas BECs to be tuned to any desired value over many orders of magnitude Pollack et al. (2009) by adjusting an external magnetic field through the phenomenon of the Feshbach resonance Inouye et al. (1998). This enables a wide range of experiments to be conducted because the properties of BECs —as well as the nature of their effective nonlinearity— crucially depend on the strength and sign of these interactions.

These advances have led to more stable, easier to use experimental settings and high-precision measurements of coherent structures in BECs. In a plethora of experiments, matter-wave dark Frantzeskakis (2010) and bright Strecker et al. (2003); Kh. Abdullaev et al. (2005); Khaykovich et al. (2002); Strecker et al. (2002); Cornish et al. (2006) solitons have been realized in single- and multi-component BECs with repulsive or attractive interatomic interactions, respectively. For example, bright solitons have been formed in ultracold 7Li gas Strecker et al. (2002); Khaykovich et al. (2002) as well as during the collapse of 85Rb condensates Cornish et al. (2006). Dark solitons have been studied in 87Rb condensates Burger et al. (1999); Anderson et al. (2001); Becker et al. (2008); Weller et al. (2008) and in sodium BECs Denschlag et al. (2000); Dutton et al. (2001). Furthermore, coupled dark-bright solitons have been engineered in 87Rb condensates using phase-imprinting methods Becker et al. (2008) or generated during superfluid-superfluid counterflow Hamner et al. (2011); Middelkamp et al. (2011a). Finally, matter wave gap solitons Eiermann et al. (2004); Morsch and Oberthaler (2006) have been produced in BECs trapped in light-induced periodic potentials.

At the theoretical level, and for sufficiently low temperatures, static and dynamical properties of BECs have been quite successfully modeled by an effective mean-field equation known as the Gross-Pitaevskii equation (GPE) Pethick and Smith (2008); Pitaevskii and Stringari (2003); Kevrekidis et al. (2015). The GPE is tantamount to a (cubic) nonlinear Schrödinger (NLS) equation with the addition of the external potential that confines the BEC. The $(2+1)$-dimensional version of the fully 3D equation reads, in terms of physical units, as

[TABLE]

where $\Phi(\boldsymbol{r},t)$ is the macroscopic BEC wavefunction, $\nabla^{2}$ is the Laplacian in $\boldsymbol{r}=(x,y)$, $m$ is the atomic mass and $g_{\rm 2D}$ describes the effective 2D strength of the atom-atom interaction. The effective 2D coupling constant $g_{\rm 2D}$ is given by $g_{\rm 2D}=g/(\sqrt{2\pi}a_{z})=2\sqrt{2\pi}\hbar a_{z}\omega_{z}a$, where $\omega_{z}$ is the harmonic trapping strength in the transverse direction, with $a_{z}$ being its corresponding harmonic oscillator length. The 3D coupling constant is $g=4\pi\hbar^{2}a/m$, where $a$ is the $s$-wave scattering length.

In the following, we set $g_{\rm 2D}>0$, that is the nonlinearity in the GPE is chosen to be defocusing Kevrekidis et al. (2015); Ablowitz et al. (2004); Sulem and Sulem (2007) which models a repulsive interatomic interaction, as is the case, e.g., in 87Rb. Multiple stationary dark-soliton states can emerge when the repulsion between dark solitons is counterbalanced by the inclusion of a trapping potential $V(\boldsymbol{r})$ in Eq. (1). The existence and formation of nonlinear patterns in BECs crucially depend on the chosen form for the applied trapping potential $V(\boldsymbol{r})$. The traditionally used magnetic traps can be adequately modeled by an harmonic external potential of the form Dalfovo et al. (1999); Kevrekidis and Frantzeskakis (2004)

[TABLE]

where, for generality, the trap frequencies $\omega_{x}$ and $\omega_{y}$ along the $x$- and $y$-direction can be chosen to be different. Static and dynamical properties of matter-wave dark solitons have been investigated in great detail in model (1) with the parabolic confining potential (2) and higher-dimensional analogues thereof. For example, dark soliton stripes and multivortex states such as vortex dipoles, tripoles, and quadrupoles have been found Möttönen et al. (2005); Pietilä et al. (2006) and their existence, stability and dynamics have been discussed in detail in the literature Middelkamp et al. (2010); Middelkamp et al. (2011b).

However, in recent years there has been increasing research activity in exploring different choices (specifically non-parabolic ones) for the external trapping potential in Eq. (1). Examples of trapping configurations recently used in BEC experiments include: double Milburn et al. (1997); Capuzzi and Hernández (1999); Holthaus (2001); Shin et al. (2005); Wang et al. (2005); Albiez et al. (2005); Zibold et al. (2010), and more-well (such as four-well Wang et al. (2009)) potentials, box potentials Pitaevskii and Stringari (2003), optical lattice potentials Pitaevskii and Stringari (2003); Adhikari and Muruganandam (2003); Choi and Niu (1999), or magnetic quadrupole trap combined with an optical dipole trap Lin et al. (2009), among many others.

In this article, we wish to explore the existence and stability of localized states in the two-dimensional (2D) GPE (1) with a ring-shaped trapping potential and repulsive interatomic interactions. A key feature of our work is the identification of a wide variety of nonlinear states in this system including ones bearing different numbers of phase jumps and associated dark solitons. The bifurcation analysis of such stationary solutions is complemented by the corresponding stability analysis, and the dynamical evolution of potentially unstable configurations. Equally importantly, phase imprinting protocols are utilized in suitably crafted numerical experiments in order to illustrate the potential of such states towards being realized in recently considered experimental setups.

More specifically, our considerations are tailored the recent experimental setup of atomtronic systems Pepino et al. (2009); Seaman et al. (2007), that are confined, neutral, ultracold atomic gases which exhibit behavior analogous to semiconductor electronic devices and circuits. In atomtronics, ring BECs are used Ramanathan et al. (2011); Wright et al. (2013); Murray et al. (2013) to realize atomic-gas analogs of superconducting quantum interference devices (SQUIDs). In Ref. Ramanathan et al. (2011), a closed-loop atom circuit was implemented for the first time in a ring-shaped confining potential. Rf SQUIDs Clarke and Braginski (2004) have been created Wright et al. (2013) in ring BECs by rotating a weak link (a localized region of reduced superfluid density) around the ring-shaped condensate. A rotating weak link was used to drive phase slips which changed the circulation around the ring and simulations, based on the GPE, showed how the circulation of the ring BEC can be probed by measuring the distribution of hole areas in time-of-flight images Murray et al. (2013). We also note in passing that ring-shaped BECs have been recently argued Eckel et al. (2017) as an interesting laboratory testbed for cosmological physics.

The article is structured as follows. In Sec. II we briefly review some of the properties of the GPE in $(2+1)$ dimensions and introduce the chosen ring-shaped trapping potential. For a detailed discussion of the existence and stability analysis of steady-state solutions in the 2D GPE with repulsive interactions we refer the interested reader to the reviews and textbooks Kevrekidis et al. (2015); Ablowitz et al. (2004); Kevrekidis et al. (2008). Our numerical results are reported in Sec. III. Finally, in Sec. IV, we summarize our conclusions and discuss possible directions for further work.

II Model and Methodology

To simplify our numerical calculations, we rewrite Eq. (1) in its well-known dimensionless form Kevrekidis et al. (2015, 2008)

[TABLE]

where $\Phi=\Phi(x,y)$ is the 2D wavefunction and $\nabla^{2}$ is the Laplacian in $\boldsymbol{r}=(x,y)$. Equation (3) is obtained from Eq. (1) by averaging (integrating) along the $z$-direction and rescaling space coordinates by the the transverse oscillator length $a_{z}$ and time by $\omega_{z}^{-1}$. Then, the density $|\Phi|^{2}$, length, time and energy are respectively measured in units of $(2\sqrt{2\pi}aa_{z})^{-1}$, the harmonic oscillator length $a_{z}=\sqrt{\hbar/(m\omega_{z})}$, the inverse trap frequency $\omega_{z}^{-1}$ and energy $\hbar\omega_{z}$.

We choose an external trapping potential as experimentally obtained from a fit provided by NIST experimentalists corresponding to a ring-shaped channel of mean radius $r_{\text{ring}}$ together with a central well of radius $r_{\text{disk}}$. Stationary ground-state condensates filling this potential (see Fig. 1) consist of a central disk surrounded by a ring thus motivating the names $r_{\text{disk}}$ and $r_{\text{ring}}$ Mathew et al. (2015); Eckel et al. (2014). This potential has the flexibility to be either a ring-plus-disk or just a ring. In the case where a ring is present, the disk can be used as a phase reference to detect phase variations in the ring caused by, e.g., stirring. Specifically, the fitted potential from the experiments takes the radial form:

[TABLE]

where $r_{\text{ring}}$, $A$ and $s_{\text{ring}}$ represent, respectively, the radius, the amplitude and the width of this ring-shaped potential. The experimentally fitted potential parameters correspond to: $r_{\text{ring}}=22.27$ $\mu$m, $r_{\text{disk}}=2.597$ $\mu$m, $s_{\text{ring}}=3.913$ $\mu$m, $s_{\text{disk}}=4.717$ $\mu$m, and $A=0.8206$. Expressed in terms of the dimensionless units of Eq. (3), based on a transverse trap frequency $\omega_{z}/2\pi=500$ Hz, these quantities correspond to: $r_{\text{ring}}=25.304738$, $r_{\text{disk}}=2.95089$, $s_{\text{ring}}=4.446226$, $s_{\text{disk}}=5.3597867$, and $A=0.8206$. A plot of the resulting ring-shaped potential is displayed in Fig. 1. Note that for ease of interpretation, we opt to display in this figure, and all subsequent ones, the spatial dimensions in the original variables, namely in microns.

Let us now construct stationary solutions of Eq. (3) by separating space and time according to

[TABLE]

where $\mu$ is the (dimensionless) chemical potential. Substituting ansatz (5) into the 2D GPE (3) yields the steady-state equation

[TABLE]

Steady-state solutions for Eq. (6) correspond to mono-parametric branches parametrized by the chemical potential $\mu$ which, in turn, fixes the number of BEC atoms in the condensate. This relationship is obtained through the conserved quantity of the GPE corresponding to the (squared) $L^{2}$ norm of the solution:

[TABLE]

Thus, after bringing back the dimensions into Eq. (7), $N$ can be identified with the mass or total number of atoms in the BEC. In what follows we find suitable starting points on a given solution branch and then vary $\mu$ using continuation methods to follow the entire branch possibly leading to bifurcations (when two solution branches collide or when new branches emanate from existing ones) as the chemical potential $\mu$ is varied Charalampidis et al. (2018); Middelkamp et al. (2011b). For given chemical potential $\mu$, we find stationary nonlinear solutions to Eq. (6) by using two different implementations of Newton algorithms. Details on these numerical methods are found in Sec. III.

After having numerically computed solutions, for each chosen value of $\mu$, we proceed to study their instability modes by performing the well-known Bogoliubov-de Gennes (BdG) stability analysis Pitaevskii and Stringari (2003); Pethick and Smith (2008); Kevrekidis et al. (2008). We perturb around a stationary solution $\phi_{0}$ using the perturbation ansatz

[TABLE]

where $(\cdot)^{\star}$ denotes complex conjugation and $\omega$ is a complex eigenfrequency. Linearization of the GPE (3) around the stationary solution $\phi_{0}$ via the ansatz (8) yields the following BdG eigenvalue problem

[TABLE]

where the matrix elements are explicitly given by

[TABLE]

We compute the eigenfunctions $\{a(x,y),b(x,y)\}$ and eigenfrequencies $\omega$ of the BdG eigenvalue problem (9) for a steady-state solution $\phi_{0}$ and for a given value $\mu$ using the $\tt{eigs}$ MATLAB routine Stewart (2002); Lehoucq et al. (1998) and our results are further checked with the Scalable Library for Eigenvalue Problem Computations (SLEPc) Hernandez et al. (2005, 2003); Roman et al. (2017). The BdG stability results are then depicted in terms of the corresponding spectra by plotting the real and imaginary parts of the eigenfrequencies as a function of $\mu$. Recall that for a linearly (neutrally) stable soliton configuration, all eigenfrequencies must be real, that is $\text{Im}(\omega)=0$.

III Numerical Simulations

Our numerical results are based on discretizing the ensuing nonlinear equations —for the dynamics Eq. (3), for the steady states Eq. (6), and for the BdG spectra eigenvalue problem Eq. (9)— on the rectangular, uniform, 2D grid $(x,y)\in[-50,50]$ and with grid spacing $\Delta x=0.2$. The steady-state equation (6) is solved using a Newton-Krylov algorithm Kelley (2003) and then, the obtained states are checked using Newton iterations implemented in the SNES libraries of PETSc Balay et al. (2017a, b, 1997).

In order to pick a suitable initial guess for convergence towards the steady state we use the first few solutions close to the linear limit. The linear limit, corresponding to weak nonlinearities in Eq. (9), may be formally identified with $N\rightarrow 0$. Then, stationary states for larger values of $\mu$ are obtained via numerical continuation by taking as initial guess the configuration calculated at nearby chemical potential values. The numerical results presented below were carried out with the chemical potential $\mu$ varying over the interval $[0,1]$ with steps of $\Delta\mu=0.002$. If not otherwise stated, all configurations depicted here correspond to the chemical potential $\mu=0.9$.

Further insights into the dynamical properties and stability of the found steady states can be obtained by perturbing these solutions with the eigenvectors, computed in the BdG linearization analysis (9), and studying their temporal evolution. To simulate the time evolution based on Eq. (3), we employ a fourth order Runge-Kutta integrator in time with a second-order finite differences used for the discretization of the spatial derivatives.

III.1 States bifurcating from the linear limit

The most basic steady state is given by the ground state. For our system with the potential given in Eq. (4), the ground state emerging from the linear limit simply corresponds to a localized “hump” of atoms that populate the central well of the potential (see panels (a) in Fig. 2). The corresponding particle number (or mass) for the ground state branch as a function of the chemical potential is depicted in Fig. 3 (see line denoted by GS). It is interesting to note that the ground state does not populate the ring of the external potential. In fact the ring does not get populated until $\mu$ reaches $\mu\simeq\mu_{\text{crit}}^{(0)}=0.313$

For $\mu\geq\mu_{\text{crit}}^{(0)}$ a new state emerges from the linear, $N\simeq 0$, limit that starts filling the ring with atoms (see panels (b) in Fig. 2). This ring-shaped solution would correspond to the ground state if the central well was absent. The mass for this ring state is depicted in Fig. 3 (see line denoted by 0S). Since this ring state could be considered as a quasi-1D periodic line of density, it is possible to think about the configurations stemming from its excited states.

For instance, in an infinite 1D line density, in the absence of external potential, the repulsive GPE admits a dark soliton solution Kevrekidis et al. (2008, 2015) corresponding to the first excited state. In the case of the ring line density, the wavefunction necessarily has to be periodic along the ring. This topological constraint restricts the number of dark solitons that can be excited along the ring to be an even number. With an even number of dark solitons along the ring, periodic boundaries are automatically satisfied.

We show these $n$-soliton steady-state solutions in Fig. 2 for $n=0$ (the ring state without any solitons), $n=2$ (a pair of dark solitons), $n=4$ (two dark soliton pairs), and $n=6$ (three dark soliton pairs). Note that, due to symmetry, in the steady state all the dark solitons must be equidistant from each other along the periodic ring. The particle numbers corresponding to these $n$-soliton solutions are depicted in Fig. 3. Note that the $n$-dark soliton solutions, populating the ring, bifurcate from the linear limit ($N\simeq 0$) and are independent of the ground state that populates the central well.

III.2 States bifurcating from the ground state

We also explored states bifurcating from the ground state. In particular, at $\mu\approx 0.560$ a double-ring solution bifurcates away from the ground state. This double-ring (see panels (a) in Fig. 4) contains the ground state populating the central well coupled to two out-of-phase rings, that populate the ring portion of the external potential, as can be seen in the top panel of Fig. 4(a) —depicting the real part of the solution— where the phase difference between the inner and outer rings is evident. Namely, this state effectively contains a ring dark soliton Theocharis et al. (2003); Frantzeskakis (2010); Kevrekidis et al. (2015) inside the outside ring channel. Figure 6 depicts the BdG spectra for the ground state and the double-ring state as a function of $\mu$. As expected, the ground state is always (neutrally) stable. However, as it is clear from the figure, the double ring is unstable since its inception. It is relevant to note that this has been recently demonstrated to be generically the case due to their azimuthal undulations in the presence of an external radial potential with the quadrupolar undulations representing the first among such spatial modes that becomes unstable Kevrekidis et al. (2017).

Interestingly, there exist further states bifurcating in turn from the double-ring solution. These states, depicted in panels (b)–(d) in Fig. 4, correspond to the double-ring with out-of-phase “petals” along the azimuthal direction. The bifurcation progression of the double-ring from the ground state and, subsequently, the states bifurcating from the double-ring is more evidently portrayed in Fig. 5 that depicts the particle numbers for these solutions as a function of $\mu$. It is relevant to note that, apparently, configurations with higher number of petals bifurcate first from the double-ring. This bifurcation cascade continues beyond what is shown in Fig. 4 (where only the first couple of bifurcating branches are depicted).

As concerns the stability of the bifurcating states, it is important to stress that the double-ring solution is unstable since its emergence from the ground state around $\mu\simeq 0.560$ and, therefore, all the subsequent bifurcating states from the double-ring inherit the instability from their double ring “ancestor” and are thus always unstable as well. Furthermore, it is interesting to note that the the first few instabilities seen in the BdG spectrum of the double ring (bottom-right panel in Fig. 6) coincide with the critical mass values corresponding to the emergence of the different bifurcating states from the double-ring configuration. Another way to state this in the language of dynamical systems is that these multi-petal states are emerging via supercritical pitchfork (symmetry breaking) bifurcations, leading to the further destabilization of the radially symmetric state via the emergence of a wide variety of azimuthally modulated ones.

Finally, in order to monitor the evolution of instabilities for the double-ring, we depict in Fig. 7 the dynamical destabilization of the double-ring. In this case, we perturb the double-ring profile with an eigenvector picked from the third instability in the BdG spectra (see bottom-right panel in Fig. 6). The wave form involving the relevant wavenumber is clearly dynamically amplified and eventually destroys the ring like structure in favor of one that bears the periodicity of the imposed perturbation.

III.3 States bifurcating from the $n$-dark soliton configurations

In a similar manner as we identified bifurcating states from the ground state and subsequently from the double-ring in the previous section, we now follow the bifurcating states from the $n$-dark soliton solutions and their associated phenomenology. For instance, Fig. 8 depicts, alongside the 2-soliton solution, its first three bifurcating states. In this case the bifurcating states pertain to excitations of the central well of the external potential. These central well excitations correspond to azimuthal, out-of-phase, “multi-petal” configurations. Figure 9 depicts the particle numbers for these configurations. In particular, the bottom panel displays the particle number difference $\Delta N$ between the central excited states and the original 2-soliton solution. When this diagnostic departs from zero, it signals the emergence of a bifurcation of a new branch from a previously existing one (with $\Delta N=0$). As shown in Fig. 10, similarly to the bifurcating states from the 2-soliton configuration, we were able to identify bifurcating states from the 4-, 6-, and 8-soliton configurations.

Now that we have identified the 2-soliton solution and its bifurcating states, let us briefly discuss the ensuing stability as a function of $\mu$. In Fig. 11 we depict the BdG spectra of the 2-soliton state (top row) together with the spectra of its first bifurcating states (bottom row) —profiles for these configurations for $\mu=0.7$ are depicted in the first two columns of Fig. 8. The BdG spectrum for the 2-soliton configuration indicates that this profile is (neutrally) stable for $\mu<0.702$. For larger values of the chemical potential (not shown here) other instabilities arise, however we do not consider them here given their much weaker growth rates.

For instance, Fig. 12 shows the long time evolution of the 2-soliton ground state heavily perturbed with an eigenvector picked from the second instability in the BdG spectrum. We observe that, when perturbed, the 2-soliton configuration develops two pairs of vortices which travel inside the ring. The vortex nature of these traveling localized solutions becomes apparent in the phase plots (see bottom row of panels) and the corresponding 2$\pi$ winding at the vortex locations. To guide the reader we have included (red) arrows that indicate the direction of motion for the vortices. As time progresses one of the vortices in each vortex pair gets “absorbed” by the edge of the ring ($t\approx 50$) leaving only two vortices of opposite charge to run along the ring. The two vortices travel towards each other, then reverse direction, move again towards each other, bounce off again etc. The vortices are found to move back and forth for a prolonged time before they annihilate for longer times (not shown here) and the configuration settles down to a slightly perturbed ring without any apparent vortices in the bulk of the ring.

Here, we omit the time evolution of instabilities corresponding to the higher excited states of the 2-soliton configuration since they do not provide any new insights into the dynamical properties. In all cases, vortices are found to travel back and forth inside the ring. For the excited states of the 2-soliton configuration, we also observe that vortices are created in the central portion of the cloud. However, those might be less relevant for experiments as the density is low there and the vortices are more tightly packed.

For completeness, we depict in Fig. 13 the BdG stability spectra for the 4-, 6- and 8-soliton solutions. As it was the case for the 2-soliton configuration, the $n$-soliton configurations are also stable for $\mu<0.702$ and the spectra are quite similar. This is straightforward to understand as the corresponding dark solitons are placed relatively far away from each other along the ring and, therefore, their mutual interaction is (exponentially) weak and thus not very noticeable when dealing with a handful of solitons. Nonetheless, higher-order excited states including a large number of dark solitons will correspond to relatively shorter mutual separations leading to stronger interactions and modifications of the stability spectra. We defer the study of such cases to future publications.

Finally, we depict in Fig. 14 the corresponding dynamical evolution for the $n$-soliton profiles for $n=4$, 6, and 8, when perturbed with eigenvectors picked from the first instability in their BdG spectra. Note that in all cases the dynamics tends to lead to the disintegration of the dark solitons (through collisions and/or splitting into vortex pairs that in turn get “absorbed” by the periphery of the ring). Eventually, and potentially after long transient stages, the evolution settles into a perturbed ring structure without dark solitons or vortices in its bulk.

III.4 Phase imprinting of $n$-dark soliton states

We now explore the especially important —in terms of a practical implementation— possibility of seeding in the experiment some of the excited state configurations that we described above. In particular, we are interested in the experimental possibility of initializing configurations that bear $n$-dark solitons and let them evolve to study their interactions and collisional dynamics. For that purpose, we start with the ring steady state depicted in Fig. 2(b). As mentioned above, this solution exists for $\mu\geq 0.313$ and it is stable for $\mu<0.702$ and therefore it is a good candidate to be attainable in a physical experiment. Then, by using a phase imprinting technique, e.g., by shining laser light on one half of the condensate for a short period of time Denschlag et al. (2000); Burger et al. (1999, 2002); Schulte et al. (2002), whereby half of the ring’s phase is shifted by $\pi$ with respect to the other half, it is possible to generate an initial condition that has the correct phase profile of a 2-dark soliton state. Such scenarios with multiple phase jumps have been previously used in quasi-1d settings in order to examine the effectively 1d interaction of dark solitary waves Stellmer et al. (2008).

We have tested that this technique is successful at seeding $n$-dark solution solutions for chemical potentials below the instability threshold around $\mu\approx 0.7$ (results not shown here). However, as we are interested not only in seeding steady states in the experiments, but also in observing the potentially unstable dynamics of these $n$-dark soliton solutions. In that light, we focus our attention here on phase imprinting $n$-dark soliton solutions past their stability threshold (i.e., $\mu$ a bit larger than 0.7). This is precisely what is depicted in Fig. 15 where the initial condition (first column of panels) corresponds the ring steady state with a phase imprinting such that the phase of the left half is $+\pi$ while the phase of the right half is [math]. As can be observed from the figure, after an initial period of adjustment ($t<300$), where the imprinted phase forces the dark soliton nucleation, a pair of dark solitons on opposite sides of the ring is formed. This configuration corresponds to a slightly perturbed 2-dark soliton state. This state, being unstable for the chosen value of $\mu$ as per the discussion in the previous sections, evolves in a manner akin to the one depicted in Fig. 12. Namely, the dark solitons start moving and colliding along the ring.

This phase-imprinting technique can be straightforwardly generalized to higher number of dark solitons by imprinting the appropriate phase. For instance, by imprinting a phase difference across the horizontal axis and then doing the same across the vertical axis, one is left with the appropriate phase to nucleate the $4$-dark soliton state. This case is depicted in Fig. 16 whose dynamical evolution in now similar to the one depicted in the first row of panels in Fig. 14. It is relevant to mention that the dynamics of the unstable $n$-dark soliton eventually leads to a perturbed ground state as the dark solitons destabilize towards the formation of vortex pairs, which in turn scatter and ultimately get absorbed by the periphery of the ring. It is natural to expect that as the ring gets thinner and more quasi-one-dimensional the relevant states will be progressively stabilized against such transverse undulations and the associated breakup towards vortex dipoles Brand and Reinhardt (2002).

III.5 A Zoo of More Exotic States

In addition to the states we constructed from the linear limit, there also exist states which bifurcate from the ground and dipole states and their excitations. Appropriate initial guesses for these states have been constructed by using the well-known ground and dipole ansätze for solutions of Eq. (6) in the presence of a harmonic external potential. For instance, as depicted in Fig. 17, there is a plethora of states bifurcating from the ground state. All of the states presented in this figure are real and pertain the combination of an $n$-dark soliton solution (populating the ring) coupled to a phase-less hump of mass localized in the central well (namely, the remnant of the ground state of the system). We have checked that all of these states are actually unstable (results not shown here). Similarly, as depicted in Fig. 18, it is possible to find more families of purely real solutions corresponding to the combination of, again, an $n$-dark soliton solution (populating the ring) but now coupled to the first excited state of the ground state (namely, the dipole consisting of a plus-minus hump at the center of the cloud). We have also checked that all of these states are actually unstable (results not shown here). This process can be extended for higher excited states of the ground state coupled to the $n$-dark soliton configuration on the ring.

Furthermore, it is also possible to find rich families of genuinely complex solutions. For instance, as seen in panels (a)–(d) of Fig. 19, it is possible to couple the $n$-dark soliton state with the ground state with a non-trivial phase difference between these two states. In the same vein, as is shown in panels (e)–(h) in Fig. 19, it is possible to couple with a non-trivial relative phase the $n$-dark soliton state with the dipole state at the center of the cloud. We have also checked that all of these states are always unstable (results not shown here).

Finally, it is relevant to mention that non-trivial phase configurations can be constructed by replacing the $n$-dark soliton solutions on the ring by a necklace of $n$-vortex solutions. These more exotic profiles are depicted in panels (i)–(k) in Fig. 19 for the case of 2, 4, and 8 vortices, respectively.

IV Conclusions and Future Challenges

We have studied the stationary and dynamical properties of BEC profiles supported by a ring-shaped potential with a target-like profile that has been used in a number of recent experiments conducted at NIST Eckel et al. (2014); Mathew et al. (2015). By following steady states and their bifurcations from the linear (low atom number) limit, we have obtained a wide range of solution branches (not all of which were shown here) and studied the corresponding stability properties as the chemical potential $\mu$ (cf. atom number) is varied. Importantly, numerous among these states were found to be potentially stable, including states carrying multiple (2-, 4-, 6-, 8-) solitons between the starting point of the respective branches and up to a suitable critical value of the chemical potential. Past this critical $\mu$ value, we studied the ensuing dynamics of the dark solitons around the ring. We typically observed that the dark solitons bounce back-and-forth in the ring until they disappear in a process involving each dark soliton splitting into a vortex pair and then the vortices getting eventually absorbed by the periphery of the ring. This process eventually led to a weakly perturbed (i.e, almost homogeneous) ring void of any dark solitons or vortices that persisted for long times.

In the case of $n$-dark soliton solutions, taking advantage of their spectral stability, we illustrated their potential for experimental realization by using phase-imprinting techniques to seed them in the condensate. We were not only able to seed stable $n$-dark soliton solutions but, equally interestingly, to seed unstable solutions whereby the ensuing dark soliton instability dynamics can be studied.

Additionally, a plethora of states was identified involving a combination of (ground or excited) states supported by the central well of the target-like potential coupled with states supported by the ring channel. The states supported by the central well corresponded to the trivial-phase ground state and its excitations in the form of dipole, quadrupole, etc. states. On the other hand, the ring channel accepts $n$-dark (equidistant) soliton solutions where $n$ is even as the periodicity of the ring enforces an even number of dark solitons. We also followed states that, instead of bifurcating from the linear limit, bifurcate from the ground state of the system (a phase-less hump populating the central well). These states correspond to double-ring, out-of-phase, solutions and “petal”-like patterns around the ring.

It would be interesting to implement the phase-imprinting methodology in the actual experiment as it would naturally allow for the study of dark soliton dynamics and interactions especially so in such an annular setup. The potential control of the spatial width of the annulus and the associated control of the snaking stability of the solitonic structures could play a significant role in the explored dynamics. From the modeling perspective it would be interesting to study the stability and dynamics of steady states bearing a large number of dark solitons. For instance, it is known that a chain of dark solitons can be approximated by a Toda lattice on the solitons’ positions and thus one can create (Toda) solitons riding on a backbone of dark solitons (see Ref. Ma et al. (2016) and references therein). Furthermore, a systematic extension of the present studies considering the vortex patterns in the present setting would naturally complement the present solitonic considerations. Lastly, considering extensions of this type of set up also in higher dimensions and suitable (e.g. toroidal-poloidal) geometries may be particularly interesting and relevant in its own right, as well as an appreciation of which (potentially vortical) patterns may be dynamically stable.

Acknowledgements

Some of the work of M.H. was undertaken as a visiting research scholar at the Department of Mathematics and Statistics, University of Massachusetts, employed by the University of Oldenburg and financially supported by FP7, Marie Curie Actions, People, International Research Staff Exchange Scheme (IRSES-605096). P.G.K. and R.C.G. and M.A.E. gratefully acknowledge the support from the National Science Foundation, under grants PHY-1602994, PHY-1603058, and PHY-1707776.