Bose-Einstein condensate confined in a 1D ring stirred with a rotating delta link

TL;DR

This paper investigates the stationary states, stability, and adiabatic manipulation of a Bose-Einstein condensate in a 1D ring with a rotating delta potential, revealing complex soliton and vortex structures.

Contribution

It provides analytical solutions and stability analysis for the condensate with a rotating delta, introducing new methods to control solitons and vortices in ring traps.

Findings

Energy spectrum exhibits swallowtail structures.

Analytical wave functions derived using Jacobi elliptic functions.

Proposed adiabatic cycles enable creation of vortices with arbitrary angular momentum.

Abstract

We consider a Bose-Einstein condensate with repulsive interactions confined in a 1D ring where a Dirac delta is rotating at constant speed. The spectrum of stationary solutions in the delta comoving frame is analyzed in terms of the nonlinear coupling, delta velocity, and delta strength, which may take positive and negative values. It is organized into a set of energy levels conforming a multiple swallowtail structure in parameter space, consisting in bright solitons, gray and dark solitonic trains, and vortices. Analytical expressions in terms Jacobi elliptic functions are provided for the wave functions and chemical potentials. We compute the critical velocities and perform a Bogoliubov analysis for the ground state and first few excited levels, establishing possible adiabatic transitions between the stationary and stable solutions. A set of adiabatic cycles is proposed in which gray solitons and vortices of arbitrary quantized angular momenta are obtained from the ground state by setting and unsetting a rotating delta.

Tables

Table 1. Table 1: Relation between the adiabatic regions, the Jacobi function with which the solutions in the region are computed, and the constraints of the function at the boundaries. The regions are defined according to the swallowtail (ST) structure of Fig.  2 and through the lines bounding them. In the first swallowtail at α < 0 𝛼 0 \alpha<0 , these lines are described by the union of the various points P i subscript 𝑃 𝑖 P_{i} as in Fig.  9 , and where P ∞ ≡ ( α → − ∞ , Ω = 1 2 ) subscript 𝑃 formulae-sequence → 𝛼 Ω 1 2 P_{\infty}\equiv(\alpha\to-\infty,\Omega=\frac{1}{2}) .

Region		$J$	Bounds
Bottom 1st ST	$P_1-P_2-P_3$	dn	$P_1-P_2$, Eq. (44)
			$P_2-P_3$, $m=1$
	$P_3-P_4-P_5-P_{\infty }$	dc	$P_3-P_4-P_5$, $m=0$
			$P_5-P_{\infty }$, $\gamma =0$
Bottom 1st ST $\alpha >0$		$\widetilde{\mathrm{dn}}$	Eq. (44)
Top 1st ST		dn	$P_1-P_2$, Eq. (44)
			$P_2-P_6$, $m=1$
Top 1st ST $\alpha >0$		$\widetilde{\mathrm{dn}}$	Eq. (44)
Bottom 2nd ST		$\widetilde{\mathrm{dn}}$	Eq. (44)
Bottom 2nd ST $\alpha >0$		dn	Eq. (44)

Table 2. Table 2: Expressions of α 𝛼 \alpha and Ω Ω \Omega for the critical points P i = ( α i , Ω i ) subscript 𝑃 𝑖 subscript 𝛼 𝑖 subscript Ω 𝑖 P_{i}=(\alpha_{i},\Omega_{i}) , i = 1 , ⋯ , 7 𝑖 1 ⋯ 7 i=1,\cdots,7 .

$P_i$	${\alpha }_{P_i}$	${\Omega }_{P_i}$
$P_1$	$0$	$\sqrt{\frac{g}{2\pi }+\frac{1}{4}}$
$P_2$		${\Omega }_{P_2}=\Omega (k_2,1)$
$P_3$	$0$	$\sqrt{\frac{g}{2\pi }}$
$P_4$	$-g$	$0$
$P_5$		$\frac{1}{2}$
$P_6$		$\frac{1}{2}$
$P_7$	$0$	$\sqrt{\frac{g}{2\pi }+1}$

Equations

$$-\frac{1}{2}{\varphi }^{\prime \prime }(\theta )+g|\varphi (\theta ){|}^2\varphi (\theta )=$$

$$\varphi (0)-e^{i2\pi \Omega }\varphi (2\pi )=$$

$${\varphi }^{\prime }(0)-e^{i2\pi \Omega }{\varphi }^{\prime }(2\pi )=$$

$${\int }_0^{2\pi }d\theta |\varphi (\theta ){|}^2=$$

$$r_J^2(\theta )=$$

$${\beta }_J(\theta )=$$

$${\tilde{\Omega }}_n=\sqrt{\frac{g}{2\pi }+\frac{n^2}{4}},$$

$$\Psi =e^{-i\mu t}(\varphi +u{e}^{-i\omega t}-v^{\ast }{e}^{i{\omega }^{\ast }t}),$$

$$-\frac{1}{2}{u}^{\prime \prime }+2g|\varphi {|}^{2}u-\mu u-g{\varphi }^{2}v=wu,$$

$$\frac{1}{2}{v}^{\prime \prime }-2g|\varphi {|}^{2}v+\mu v+g{{\varphi }^{\ast }}^{2}u=wv.$$

$$u(0)-e^{i2\pi \Omega }u(2\pi )=$$

$$u^{\prime }(0)-e^{i2\pi \Omega }{u}^{\prime }(2\pi )=$$

$$v(0)-e^{-i2\pi \Omega }v(2\pi )=$$

$$v^{\prime }(0)-e^{-i2\pi \Omega }{v}^{\prime }(2\pi )=$$

$$u(\theta )=e^{-i\Omega \theta }\tilde{u}(\theta ),$$

$$v(\theta )=e^{i\Omega \theta }\tilde{v}(\theta ),$$

$$i\hslash {\partial }_{t_L}{\psi }_L=$$

$$+\frac{\alpha }{2}\delta ({\theta }_L-\Omega t_L){\psi }_L,$$

$${\psi }_L(\Omega t_L,t_L)=$$

$$\frac{m{R}^2\alpha }{{\hslash }^2}\psi (\Omega t_L,t_L)=$$

$$-({\partial }_{{\theta }_L}\psi ){|}_{{\theta }_L=\Omega t_L+2\pi }.$$

$$\theta ={\theta }_L-\Omega t_L,\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{\partial }_{\theta }={\partial }_{{\theta }_L},$$

$$t=t_L,\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}{\partial }_t={\partial }_{t_L}+\Omega {\partial }_{{\theta }_L},$$

$$\psi (\theta ,t)=e^{\frac{i}{\hslash }\left(\frac{1}{2}m{R}^2{\Omega }^2{t}_L-mR\Omega {\theta }_L\right)}{\psi }_L({\theta }_L,t_L),$$

$$i{\partial }_t\psi =$$

$$\psi (0,t)=$$

$$r_J^2(\theta )=$$

$${\beta }_J^{\prime }(\theta )=$$

$${\theta }_j=\pi -\frac{j}{k}K(m).$$

$$r_{\mathrm{dn}}(\theta )=$$

$$r_{\mathrm{dc}}(\theta )=$$

$${\eta }_{\mathrm{dn}}=$$

$$+E[\mathrm{JA}(k{\theta }_j,m),m],$$

$${\eta }_{\mathrm{dc}}=$$

$${\gamma }_J=$$

$$\times \sqrt{g-2\pi k^2(1-m)+k{\eta }_J},$$

$${\mu }_J=$$

$${\mu }_{\alpha =0}=$$

$$+\frac{6(n-1{)}^2}{\pi }K(m)E(m)].$$

$${\tilde{\Omega }}_n=\sqrt{\frac{g}{2\pi }+\frac{n^2}{4}},$$

Full text

I Introduction

Bose-Einstein condensates (BECs) constrained in annular traps provide a way to study various phenomena related to superfluidity, including persistent currents and their decay, phase slips, and critical velocity ryu07 ; moulder12 ; beattie13 ; wright13pra . Persistent currents can be created experimentally by the application of artificial gauge fields dalibard11 , or by rotating a localized, tunable repulsive barrier around the ring ramanathan11 ; piazza09 ; piazza13 ; wright13prl . Through the latter method, hysteresis between different circulation states was observed in eckel14 . Within the context of atomtronics, a ring condensate is a key atomic circuit element. It has demonstrated its capability as a superconducting quantum interference device ryu13 , entailing the possibility of high precision measurements and applications in quantum information processing hallwood10 ; schenke11 ; amico14 .

In the view of a better control of BECs, phase transitions have been analyzed in different ring settings within the mean field approach. They were first studied in a ring under a rotational drive kanamoto09 , and then through the interplay between rotation and symmetry breaking potentials or rotating lattice rings. One lattice site was studied in fialko12 , a double well in li12 , and a more general unified approach of a ring lattice in munoz19 , all involving the possibility to adiabatically connect different quantized states such as persistent currents or solitons.

By solving the Gross-Pitaevski equation (GPE), various works have studied the energy diagram and metastability of BECs in rings with a rotating defect baharian13 ; munoz15 ; kunimi18 . In the case of the 1D GPE, stationary solutions can be found through the inverse scattering method or by directly integrating and writing them in terms of Jacobi functions. These solutions have been analyzed under box and periodic boundary conditions carr002 , under a rotational drive kanamoto09 , and under some specific constant potentials seaman05 . The flow past an obstacle in the form of a Dirac delta was studied perturbatively in hakim97 ; pavloff02 . In cominotti14 a 1D ring with a rotating Dirac delta was analyzed for some specific rotations, strengths, and nonlinearities.

In this paper, we study a repulsive BEC in a 1D ring where a Dirac delta link is rotating at constant speed. The use of analytical solutions, expressed in terms of Jacobi functions, allows us to compute the stationary wave functions and chemical potentials for the ground state and an arbitrary number of excited energy levels. The obtained energy diagram, depending on the delta velocities and strengths, both attractive and repulsive, is analyzed as a function of the coupling strength. This diagram entails a series of critical velocities which, together with a Bogoliubov analysis, lay out the distribution of stable and metastable states in parameter space, and which adiabatic transitions between them are possible. Within these transitions, we propose a few adiabatic cycles in which excited solitonic and vortex states are produced by setting and unsetting a rotating delta.

This paper is organized as follows. In the next section, II, we introduce the theoretical model, defining the GPE and boundary conditions in the Dirac delta comoving frame, and provide a method to compute the spectrum. The results are in Sec. III, in which we illustrate the main features of the spectrum (III.1), its stability (III.2), and its dependence on the nonlinearity (III.3). In Sec. IV, we propose a set of adiabatic paths to excite the condensate. We conclude this paper in Sec. V. Mathematical details are found in the Appendices.

II Theoretical model

We consider a BEC at zero temperature in a tightly transverse annular trap in which a Dirac delta link is rotating at constant speed. The point-like potential is chosen instead of a finite one such that analytical solutions can be obtained, with the view that the results may not qualitatively change with respect to a very peaked Gaussian. Considering only stationary solutions, and within the mean field approach, we can determine the condensate wave function in the delta comoving frame, $\phi(\theta)$, by the 1D Gross-Pitaevskii equation. Then, $\phi(\theta)$ is constrained by delta boundary conditions and normalization,

[TABLE]

where $g>0$ is the reduced 1D coupling, $\mu$ the chemical potential, $\theta\in[0,2\pi)$, and $\Omega$ and $\frac{\alpha}{2}$ the velocity and strength of the delta link, see App A. Here and in the rest of the paper we use units $\hbar=R=M=1$, $R$ being the radius of the ring and $M$ the mass of the atoms. Renormalizing a wave function $\phi(\theta)\to\sqrt{N}\phi(\theta)$ amounts to a rescaling of $g\to g\,N$. We choose to fix the normalization and study how the spectrum depends on $g$.

Any solution of Eq. (1), $\phi(\theta)=r(\theta)e^{i\,\beta(\theta)}$, can be written in closed form in terms of a Jacobi elliptic function seaman05 . In particular, the density $\rho(\theta)\equiv r(\theta)^{2}$ depends linearly on the square of one of the twelve Jacobi functions ($J$), and the phase $\beta(\theta)$ is fixed by the density,

[TABLE]

where $k$ is the frequency and $m\in(0,1)$ the elliptic modulus, which generalizes the trigonometric and hyperbolic functions into the Jacobi ones. The constants $A$ and $B$, the shift $\theta_{j}$, $\gamma=\rho(\theta)\beta^{\prime}(\theta)$, a constant representing the current, and $\alpha$, $\Omega$ and $\mu$, are fixed by Eqs. (1)-(4) in terms of $k$ and $m$. The spectrum $\mu(\alpha,\Omega)$ is thus given in parametric form, $(\alpha(k,m),\Omega(k,m),\mu(k,m))$. Running $k$ and $m$ in a systematic way in the three possible solutions allowed by Eqs. (1)-(4), dn, $\tilde{{\rm dn}}$, and dc (see App. B) we obtain $\mu(\alpha,\Omega)$ as a series of surfaces which fold onto each other —energy levels which cross and are degenerate at specific lines $\Omega_{cr}(\alpha)$. Any solution found for a specific $\alpha$ and $\Omega$, also satisfies Eqs. (1)-(4) with $\Omega\to\Omega\pm integer$ and $\beta(\theta)\to\pm\beta(\theta)$. To obtain the complete spectrum, we shift and mirror the obtained spectrum to $\Omega\to\Omega\pm integer$.

III Static properties

Our goal is to analyze the possible stable and adiabatic changes of the condensate as one varies the strength and velocity of the Dirac delta. For this we first study the structure of the spectrum $\mu(\alpha,\Omega)$, i.e. the regions in the $\alpha-\Omega$ plane in which stationary solutions exist for the ground and first excited states, and how the chemical potential depends on $\alpha$ and $\Omega$ (Sec. III.1). Then we analyze whether the solutions at each region are stable or metastable against a perturbation through a Bogoliubov analysis (Sec. III.2). The results in Sec. III.1 and III.2 are analyzed and illustrated for $g=10$. In Sec. III.3 we study how they depend on $g$.

III.1 Spectrum

The spectrum $\mu(\alpha,\Omega)$ consists in a set of surfaces which cross and merge at certain boundaries, as the ones plotted in Fig. 1 and their symmetric versions at $\Omega\to\Omega+integer$. Five sections of $\mu(\alpha,\Omega)$ with constant $\alpha$ are plotted in Fig. 2. Except for the ground state at $\alpha<0$, they present a set of concatenated swallowtail (ST) shapes. For any given $\alpha$, each series of swallowtail diagrams represents a set of stationary solutions continuously connected among them through the parameters $k$ and $m$, or through the velocity $\Omega$ and chemical potential $\mu$. These type of energy diagrams are characteristic of hysteresis and were analyzed in the context of ring condensates in mueller02 . We organize and label them, and the surfaces that they constitute, according to their position, ordered from lower to higher energy. For each diagram we distinguish a bottom and a top part, both merging at the tip of the swallowtail. In the following we present their general features. The structure of the spectrum is analyzed more thoroughly in App. B.

Each energy level —top or bottom part of a swallowtail diagram— is symmetric with respect to $\Omega=\frac{l}{2}$ and bounded by a pair of critical velocities $\frac{l}{2}\pm\Omega_{cr}(\alpha)$, with $l$ an integer. Solutions corresponding to a level centered at a link velocity $\Omega=\frac{l}{2}$, are a boost from those in the analogous level at $\Omega=0$ or $\Omega=\frac{1}{2}$. The $n$th set of swallowtail diagrams, $n=1$ being the bottom one, entails densities with $n$ depressions, considering both the valleys characteristic of the Jacobi functions, and the downward kinks in the case of $\alpha>0$. We distinguish three types of solutions, depending on the depth of the depressions: vortex states, dark solitons, and gray solitons.

Vortex states. The red parabolas of plot (c) ($\alpha=0$) in Fig. 2 centered at $l$ represent the chemical potential of vortex states of angular momentum $l$ as observed from the frame moving at $\Omega$. The minima of $\mu(\Omega)$ are at $\Omega=l$, where the observer is comoving with the vortex.

Dark solitons. At precisely $\Omega=\frac{l}{2}$, solutions consist in dark solitonic trains (except for the ground state), see the black dots in Fig. 2. In the delta comoving frame, the dark solitonic trains are stationary, with zero current and constant phase, except for a phase jump of $\pi$ at each zero in the density. They correspond to the minima of the energy spectrum $\mu(\Omega)$ for any particular fixed $\alpha$. In the lab frame, they comove with the condensate and the delta link at $\Omega=\frac{l}{2}$.

Gray solitons. Solutions with velocities that depart from $\Omega=\frac{l}{2}$, consist in gray solitonic trains, with shallower waves and faster currents the larger $|\Omega-\frac{l}{2}|$. At $\alpha=0$, gray solitonic trains with $n$ depressions become completely flat and merge with vortex states at $\Omega=\frac{l}{2}\pm|\tilde{\Omega}_{n}-\frac{n}{2}|$, with

[TABLE]

see App. B or carr002 . At $\alpha\neq 0$, the rotational symmetry is broken, and a pair of critical velocities, corresponding to the tips of the swallowtails, limit the range of $\Omega$ for which stationary solutions exist. In particular, the width of the bottom part of the $n$th swallowtail centered at $\frac{l}{2}$ monotonously decreases from $\frac{l}{2}\pm|\tilde{\Omega}_{n}-\frac{n-1}{2}|$, at $\alpha=0$, to $\frac{l}{2}\pm\frac{1}{2}$ at $|\alpha|\to\infty$. The condensate in the ring is therefore able to sustain stationary solutions, consisting in gray solitonic trains comoving with the kink, up to a certain stirring velocity —relative to $\frac{l}{2}$—, which decreases with the strength of the delta link.

A sample density and phase for each colored region of Fig. 1 (and the corresponding ones at $\alpha<0$), are plotted in Fig. 3. The densities corresponding to the bottom of first swallowtail surface at $\alpha<0$ entail an upward kink. This kink becomes higher and more peaked as the delta potential becomes more attractive, and can be understood as a bright soliton. In contrast, for the rest of levels, as $\alpha\to\pm\infty$, the density at the delta position becomes zero. The densities corresponding to the first swallowtail diagram have one depression, which for $\alpha<0$ consists of one valley and for $\alpha>0$ a downward kink. Similarly, the four plots corresponding to the second swallowtail levels have two depressions. For $\alpha>0$, one of these depressions is also understood as the downward kink imposed by the delta. From these plots it can also be inferred the relation between the the phase and the density, $\beta^{\prime}(\theta)=\frac{\gamma}{\rho(\theta)}$, which implies higher phase gradients (velocities) for lower densities.

III.2 Metastability

A Dirac delta with fixed strength and rotating at a constant speed allows an infinite set of solutions organized in chemical potential levels. The stability of these solutions can be studied by adding a small perturbation to the stationary wave function

[TABLE]

and analyzing how it evolves. Replacing this function in the time-dependent Gross-Pitaevskii equation (Eq. (23)) and linearizing in $u$ and $v$, we obtain the corresponding Bogoliubov system of equations bogoliubov47 ,

[TABLE]

Both perturbations, $u$ and $v$, must satisfy the boundary conditions separately, which read,

[TABLE]

We solve this system of equations, by changing variables to

[TABLE]

and then expanding $\tilde{u}$ and $\tilde{v}$ in an orthonormal basis, thus converting it into a matrix eigenvalue problem, and by the Direct and Arnoldi methods integrated in the differential solvers in Mathematica, see App. E for more details.

Both methods yield the same eigenvalues, and analyzing whether they are real or complex, we split each region in stable and unstable parts.

The bottom of the first swallowtail and both upper parts are found completely stable and unstable, respectively. In contrast, the lower part of the second swallowtail is stable except for regions of metastability in form of stripes in the plane $\alpha-\Omega$ at $\alpha<0$, as shown in Fig. 4.

III.3 Dependence on nonlinearity

The results presented in the previous subsections have been illustrated with the nonlinearity fixed to $g=10$. In Fig. 5 we show the same adiabatic regions for $g=20,10,5$, and $1$. The first and third columns correspond to the bottom part of first and second swallowtail diagrams (as in Fig. 2), while the middle column represents the top part of the first one. Larger nonlinearity implies a greater span and overlap of all the levels. As $g$ decreases, and for any fixed $\alpha$, the tail part of the swallowtail diagrams becomes smaller, vanishing at $g=0$: both the region in the middle column and the overlaps of shifted regions in the others decrease in size. A condensate with larger nonlinearity is therefore able to sustain stationary solutions for faster stirring velocities, while in the linear limit, $g=0$, the critical velocities are independent of the delta strength and fixed to $\Omega=\frac{l}{2}\pm\frac{1}{2}$.

Performing the same metastability analysis we find that the bottom part of the first swallowtail and both upper parts remain completely stable and unstable for all the $g$ tested. The third region of Fig. 5, corresponding to the bottom of the second swallowtail diagram, is also found completely stable at $\alpha>0$, while the metastable stripes at $\alpha<0$ become thiner (thicker) as $g$ decreases (increases).

IV Adiabatic generation of vortex states and excited solitons

All complex solutions found are continuously connected at $\alpha=0$, where gray solitonic trains merge into vortex states as their velocity $\Omega$ departs from $\frac{l}{2}$, with $l$ an integer (see middle panel of Fig. 2). Once a finite delta strength $\frac{\alpha}{2}\neq 0$ is set, a gap between the various swallowtail diagrams appears, and depending on the initial rotational velocity and initial state, different energy levels can be accessed. In particular, setting a delta link in the ground state while rotating at $\Omega_{i}\in(\tilde{\Omega}_{n},\tilde{\Omega}_{n+1})$, $n\geq 1$, a gray solitonic train with $n+1$ depressions is obtained—reaching thus the bottom of the $n+1$th swallowtail diagram. As the velocity is decreased, the valleys become deeper, and at $\Omega=\frac{n}{2}$ the solution turns into a wave train with $n$ dark solitons. If the velocity is further decreased, the density profile becomes a gray solitonic train again, with shallower waves with lower velocity. At any velocity, one can unset the delta link, making the waves shallower as the delta strength decreases back to $\alpha=0$. For $\Omega>\tilde{\Omega}_{n}$, unsetting the link completely flattens the density, and the ground state is recovered. At velocities $\Omega\in(-\tilde{\Omega}_{n}+n,\tilde{\Omega}_{n})$, the final state consists in a gray solitonic train, with deeper waves the closer $\Omega$ is to $\frac{n}{2}$. If one unsets the delta link at precisely $\Omega=\frac{n}{2}$, a dark solitonic train is obtained. In the case in which $\alpha$ is brought back to zero at $\Omega<-\tilde{\Omega}_{n}+n$, the waves also become infinitely shallow, and merge with the vortex state of $n$ quanta.

Cycles in which the final state is a vortex can be concatenated. Once a vortex is produced, any observer can rotate at the same velocity as the condensate current, such that in the comoving frame the vortex is observed as the ground state. The cycle to produce a vortex from the ground state can then be repeated. Let us consider the vortex state of five quanta of angular momenta. It can be produced by setting a delta link rotating at $\Omega_{i}\in(\tilde{\Omega}_{5},\tilde{\Omega}_{6})$, decreasing its velocity to $\Omega_{f}\in(-\tilde{\Omega}_{6}+6,-\tilde{\Omega}_{5}+5)$, and then unsetting the delta. It involves a middle step in which a wave train with five dark solitons is produced. Another possibility is to concatenate five analogous cycles in which $\Omega_{i}\in(\tilde{\Omega}_{1},\tilde{\Omega}_{2})$, and $\Omega_{f}\in(-\tilde{\Omega}_{2}+2,-\tilde{\Omega}_{1}+1)$, each involving the production of only one dark soliton. After each individual cycle, the observer is boosted $\Omega\to\Omega+1$.

Any of these proposed paths should avoid the metastable and non-stationary regions analyzed in the previous section. In Fig. 6, we schematically demonstrate four of such adiabatic paths. They excite the condensate to one dark soliton and a vortex state, each one through the setting and unsetting of a delta link, either repulsive or attractive. In both cases, the delta strength is limited by the line $\Omega_{cr}(\alpha)$ bounding the adiabatic region corresponding to the bottom of the second swallowtail. For the attractive delta, one also needs to avoid the metastable region, which further limits the value of $\alpha$. These adiabatic cycles can be reproduced for all $g$ tested. However, the range of velocities at which the Dirac delta strength must be set, decreases with $g$, $\tilde{\Omega}_{2}-\tilde{\Omega}_{1}=\frac{3}{4}\sqrt{\frac{\pi}{2g}}+\mathcal{O}(g^{-\frac{3}{2}})$. Moreover, in the case of the cycle with an attractive rotating delta, the constrain on the magnitude of the delta potential will depend on the metastability stripes shown in Fig. 5.

For any adiabatic path involving higher swallowtail levels, the corresponding metastability analysis and determination of critical velocity should be performed. Alternatively, one can simulate these cycles by solving the time-dependent Gross-Pitaevskii equation with the Dirac delta link replaced by a Gaussian one. We solve this differential equation using the method of lines, and reproduce the cycles proposed in Fig. 6. Moreover, we also obtain dark solitonic trains and vortex states with up to five solitons or quanta. On the one hand, these simulations validate the results found for the delta link. On the other, they indicate that the structure of the spectrum laid out in the delta case is also able to depict the main features of the spectrum of a BEC stirred with a peaked Gaussian link.

The velocities $\tilde{\Omega}_{n}=\sqrt{\frac{1}{2\pi}\frac{g}{g_{nat}}+\frac{n^{2}}{4}}\,\Omega_{nat}$, constraining the adiabatic cycles, have been presented in natural units, where $\Omega_{nat}=\frac{\hbar}{MR^{2}}=1$ and $\alpha_{nat}=g_{nat}=\frac{\hbar^{2}}{MR}=1$. For a condensate of 87Rb atoms, with mass $M=86.909$ u and ring traps of radiuses $R_{1}=20{\,\rm\mu m}$ to $R_{2}=100{\,\rm\mu m}$, the natural velocities range from $\Omega_{nat}^{(1)}=1.83$ rad/s to $\Omega_{nat}^{(2)}=0.073$ rad/s. Taking $g=10\,\tilde{g}$, the thresholds speeds to access the first excited states are $\tilde{\Omega}_{n}=2.48$, 2.95, 3.59 rad/s for the smaller radius, and $\tilde{\Omega}_{n}=0.099$, 0.118, 0.143 rad/s for the larger one. In the linear limit, $g=0$, $\tilde{\Omega}_{n}=\frac{n}{2}\Omega_{nat}$.

V Conclusions

The spectrum of a 1D ring condensate with a Dirac delta rotating at constant speed has been analyzed in terms of the nonlinearity $g$, the delta velocity $\Omega$, and the delta strength $\alpha/2$. Analytical expressions are provided for the wave function, the current, and the chemical potential. For a fixed $g$ and $\alpha$, the dependence of the chemical potential on the delta velocity, $\mu(\Omega)$, consists of a series of swallowtail diagrams. These diagrams can be organized from smaller to larger energies, and each one can be split into a bottom and a top part. The lowest diagram at $\alpha<0$ is an exception, and consists only of the bottom part, with a more complex structure depending on a set of critical points $P_{i}$. As the magnitude of the delta strength increases, the sizes of the tails in each diagram decrease, while the energy gap among each diagram becomes larger. At $\alpha=0$, the top parts of the diagrams merge with the bottom parts of the immediate upper ones. The spectrum $\mu(\alpha,\Omega)$ thus consists in a multiple swallowtail 3D structure, each region providing a range of delta strengths and velocities which can be varied adiabatically to access different solitonic solutions. These solutions consist in gray or dark solitonic trains, where the number of depressions in each train has been related to the position of the swallowtail. In particular an odd (even) number of dark solitons comove with the condensate at $\Omega=\frac{l}{2}$, where $l$ is an odd (even) integer.

In order to support the possible adiabatic processes allowed by the Gross-Pitaevskii spectrum, we have analyzed the metastability of each solution for the first two swallowtail diagrams. The top parts are found unstable and the bottom ones stable, except for the bottom of the second diagram at $\alpha<0$, which presents a series of metastable stripes in parameter space.

We have proposed a method to produce dark and gray solitons, and vortex states of arbitrary quantized angular momentum, by controlling the stirring velocity and strength of the potential. The method consists in setting and unsetting a Dirac delta potential while rotating it around the condensate at certain velocities. In particular, as a rotating observer sets a delta link at $\Omega_{i}\in(\tilde{\Omega}_{n},\tilde{\Omega}_{n+1})$, the bottom part of the $n+1$th swallow tail is reached, and a gray solitonic train with $n+1$ depressions is produced. The cycles in the parameter space defined by $\alpha$ and $\Omega$ corresponding to the various production processes are constrained by the width of the swallowtail diagrams, and also by the metastable regions. These adiabatic paths are qualitatively reproduced by solving the time-dependent Gross-Pitaevskii equation in which a finite width Gaussian link is rotating at constant speed.

Appendix A Gross-Pitaevskii equation and boundary conditions

The evolution of the condensate wave function in the Lab frame, $\psi_{L}(\theta_{L},t_{L})$, is governed by the 1D Gross-Pitaevskii equation,

[TABLE]

where $M$ is the atomic mass, $R$ the radius of the ring, $\theta_{L}\in(0,2\pi)$ and $t_{L}$ the angular and time coordinates in the lab frame, $\alpha/2$ and $\Omega$ the magnitude and angular velocity of the Dirac delta, and $g>0$ the reduced 1D coupling strength. The circular topology imposes continuity conditions in the wave function,

[TABLE]

and the Dirac delta constrains its derivatives through boundary conditions. These are obtained by integrating Eq. (17) in a small contour around the delta, $\theta_{L}\in(\Omega t_{L}-\epsilon,\Omega t_{L}+\epsilon)$, and taking the limit $\epsilon\to 0$,

[TABLE]

We change variables to the delta rotating frame leggett73 ,

[TABLE]

and use units $\hbar=M=R=1$. Then eqs. (17), (18), and (19) become,

[TABLE]

For a stationary solution, $\psi(\theta,t)=e^{-i\,\mu t}\phi(\theta)$, where $\mu$ is the chemical potential, these equations result in Eqs. (1)-(3).

Appendix B Solutions

To obtain the solutions of Eq. (1) we write the wave function as $\phi(\theta)=r(\theta)e^{i\,\beta(\theta)}$. Separating into real and imaginary parts, and integrating, the density and phase take the general form,

[TABLE]

where $J$ is one of the 12 Jacobi functions and $A$, $B$, $k$, $\theta_{j}$, $m$, and $\gamma_{J}$ are constants. The squares of the six convergent or divergent Jacobi functions are related among themselves linearly and through shifts in $\theta$, and therefore one may consider only a convergent one and a divergent one with general $A$, $B$, and $\theta_{j}$. Eq. (27) represents the stationarity condition, with $\gamma_{J}=r_{J}^{2}\beta_{J}^{\prime}$ the current. The shift $\theta_{j}$ is fixed by the continuity condition, $r(0)=r(2\pi)$: the angular length of the condensate, $2\pi$, has to be equal to an integer number of periods ($j\,T$) plus twice the shift, $2\pi=j\,T+2\theta_{j}$ (see Fig. 7).

The period of $r(\theta)$ is given in terms of the elliptic integral of first kind ($K(m)$), $T=\frac{2K(m)}{k}$, and therefore

[TABLE]

Since we take $k$ and $m$ as a parameters, $j$ can be fixed to $j=0,1$. Eqs. (1) and (4) fix $A$, $B$, $\gamma_{J}$ and $\mu_{J}$ in terms of $\theta_{j}$, $k$, and $m$. Using the Jacobi functions dn and dc as the convergent and divergent independent solutions, respectively, the amplitudes read,

[TABLE]

where

[TABLE]

with $E$ the elliptic integral of the second kind, ${\rm JA}$ the Jacobi amplitude, sc the Jacobi function, and where dc allows only for $j=0$ and $k<K(m)/\pi$ in order to be convergent in $\theta\in[0,2\pi)$. The phases $\beta_{{\rm dn}}$ and $\beta_{{\rm dc}}$ are then integrated from the Eq. (27).

For the types of Jacobi functions chosen, $J={\rm dn},{\rm dc}$, the corresponding currents and chemical potentials read,

[TABLE]

This leaves the frequency $k$ and elliptic modulus $m$ as the only free parameters. They constrain $\alpha$ and $\Omega$ through the boundary conditions in Eqs. (2) and (3). $k$ and $m$ are either real, $k>0$, $m\in[0,1]$, or, in the case of real solutions, may also take complex values with $|m|=1$ and $k\propto 1/\sqrt{1+m}$. For the real solutions (with general real boundary conditions) we refer to perezobiol19 . The elliptic modulus is further constrained by the condition that $\gamma_{J}\in\mathcal{R}$, which is satisfied when and odd number of radicants in Eq. 33 are positive. Note that the transition to an even number of radicants being negative, where $\gamma_{J}$ is not real, happens at $\gamma_{J}=0$.

In the case of $\alpha=0$, the wave functions are plane waves, where $\mu_{n}=\frac{g}{2\pi}+\frac{1}{2}(\Omega+n)^{2}$, or solitonic trains, which have, for each level $n$,

[TABLE]

Both expressions correspond to periodic boundary conditions (solved in Ref. carr002 ) and coincide at $m=0$ and $k=\frac{n-1}{\pi}K(0)=\frac{n-1}{2}$. These values determine the critical velocity,

[TABLE]

that bounds both regions $n$ and $n+1$ at $\alpha=0$.

For a similar treatment of the Jacobi functions and a complete derivation of the solutions see e.g. carr002 . The main difference between carr002 and our work is that in carr002 $k$ and $m$ are not taken as parameters to account for a phase jump and a kink, but adjusted such that periodic boundary conditions are obtained.

Appendix C Computation of the spectrum

All solutions $\phi=r\,e^{i\,\beta}$ satisfying Eqs. (1)-(4) can be obtained by running $k$ and $m$ in their allowed ranges in any of the three Jacobi functions, two convergent and one divergent, and which we label as

[TABLE]

Then, for each $k$ and $m$, the delta strength, velocity, and the chemical potential are obtained from Eqs. (2), (3), and (1),

[TABLE]

A sample of the spectrum $\mu(\alpha,\Omega)$ produced this way is shown in plot (a) of Fig. 8. Its projection on the $\alpha-\Omega$ plane is potted on panel (b) of the same figure. In this plane, lines parametrized by $k$ and fixed $m$, ($\alpha_{m}(k)$,$\Omega_{m}(k)$), and lines parametrized by $m$ and fixed $k$, ($\alpha_{k}(m)$,$\Omega_{k}(m)$), are tangent at the degeneracy line,

[TABLE]

see bottom-right panel of Fig. 8. Therefore, any degeneracy line may be computed by solving

[TABLE]

Appendix D Spectral structure

Fig. 9 shows the projections of $\mu(\alpha,\Omega)$ in the $\alpha-\Omega$ plane for the ground state and first excited levels. The ground and first excited state for $\alpha<0$ present a more complex structure, which we analyze in this Appendix. Its bounds are determined by the limits of the elliptic modulus $m$ —0, 1, or the ones fixed by $\gamma_{J}=0$—, or by Eq. (44), and conform the set of lines uniting the points $P_{1}-(\cdots)-P_{5}$ and $P_{1}-P_{2}-P_{6}$. Each segment of these bounds is listed in Table 1 together with the function and limit of $m$ they represent. The line $P_{2}-(\cdots)-P_{5}$ and its continuation at $\Omega=\frac{1}{2}$ is the only curve not defining the degeneracy of two energy levels. The points $P_{i}$, $i=1,7$, themselves can be further determined. Points $P_{1}$, $P_{3}$, and $P_{7}$ correspond to $\alpha=0$ and velocities $\tilde{\Omega}_{1}$, $\tilde{\Omega}_{0}$, and $\tilde{\Omega}_{2}$, respectively. The values of $P_{4}$, $P_{5}$ are constrained by their current and elliptic modulus being zero, $\gamma_{J}=m=0$, while $P_{6}$ has $\gamma_{J}=0$ and $m\to 1$. Point $P_{2}$ is obtained by minimizing $\alpha$ with $m=1$. All these constrains fix $P_{i}$ to the values shown in Table 2.

The dependence of the spectrum on $g$ can be analyzed quantitatively through the expressions in Table 2 for the points $P_{i}$. $P_{1}$, $P_{3}$, $P_{4}$, and $P_{7}$ are given in analytical form, and $P_{2}$, $P_{5}$, and $P_{6}$ are plotted in Fig. 10. $\alpha_{P_{2}}$ approaches zero as $g$ increases, and the structures bounded by the lines $P_{1}-P_{2}-P_{3}$ and $P_{1}-P_{2}-P_{6}$ (middle panel in Fig. 9) vanish in the limit $g\to\infty$. In contrast, $|\alpha_{P_{4}}|$, $|\alpha_{P_{5}}|$ and $\Omega_{P_{3}}$ increase with $g$, and the region bounded by these points grows at large interactions. At $g=-\frac{2}{\pi}+2\pi k_{cr}^{2}\simeq 0.280$, where $k_{cr}\simeq 0.382$ solves $\pi k_{cr}\tanh(k_{cr}\pi)=1$, both Eqs. for $k_{2}$ and $k_{6}$ in Table 2 are satisfied and points $P_{2}$ and $P_{6}$ coincide. In the limit $g\to 0$, $P_{6}$ approaches $P_{5}$ at $\alpha=-\frac{2}{\pi}$, $k_{5}$ and $k_{6}$ tend to zero, and $\Omega_{P_{1}}=\frac{1}{2}$: the parts of the region merge into a flat band. In general, at $g=0$, all levels turn into regions spanning $\Omega\in[l,l+1]$, where all the solutions are stable (see App. F for the linear solutions).

Appendix E Bogoliubov analysis

The differential equations and boundary conditions constraining $\tilde{u}$ and $\tilde{v}$ are

[TABLE]

[TABLE]

In order to turn this system of equations into a linear eigenvalue problem, we expand $\tilde{u}$ and $\tilde{v}$ in an orthonormal basis. This basis does not consist in periodic plane waves, since the derivatives must be discontinuous according to the delta conditions, but in the solutions of Eqs. (1)-(4) with $g=0$ and $\Omega=0$. Imposing these constraints on exponential and trigonometric functions, we obtain the basis,

[TABLE]

with $n$ a positive integer, and where the element $s_{0}(\theta)$ is only used for $\alpha<0$. $\tilde{u}$ and $\tilde{v}$ expanded in this set of functions solve Eqs. (47) and (49), and Eqs. (48) and (50) are satisfied as long as $k_{0}$ and $k_{n}$ are the solutions of, respectively,

[TABLE]

Appendix F Linear limit

Solutions of Eqs. (1)-(4) with $g=0$ can be found analytically proceeding analogously to App. B and replacing Jacobi functions by trigonometric and hyperbolic ones. They read,

[TABLE]

where,

[TABLE]

and the frequency $k$ is real and the current $\gamma$ positive. Note that $r_{c}$ and $r_{ch}$ solutions are related by $k\to i\,k$. For a given $k$, $\gamma$ is limited by the square roots in $A_{c}$, $A_{ch}$ being real. The phases, $\alpha$ and $\Omega$, are computed according to Eqs. (27), (40), and (41), respectively, where now $k$ and $\gamma$ are taken as parameters, and the chemical potentials read,

[TABLE]

The spectrum consists in a series of layered levels, each one spanning all $\alpha$ and $\Omega\in[\frac{n}{2},\frac{n}{2}+\frac{1}{2}]$, given that $|k|\in[\frac{n}{2},\frac{n}{2}+\frac{1}{2}]$, where $n=0,1,2,$ etc., see Fig. 11. Adding a perturbation to these solutions in the form of Eq. (8) must satisfy the same linear equations with $\mu\to\mu\pm w$, as stated by Eqs. (9) and (10) with $g=0$. The solutions only satisfy the boundary conditions for real eigenvalues, and therefore the frequencies $w$ are not imaginary and all the solutions stable.

Acknowledgments

This work was supported by the Japan Ministry of Education, Culture, Sports, Science and Technology under the Grant number 15K05216. We thank Muntsa Guilleumas, Bruno Juliá-Díaz and Ivan Morera for fruitful discussions and a careful reading of the paper. The article has also greatly benefited from ideas and comments from the Atomtronics 2019 community.