Quantum sensing using imbalanced counter-rotating Bose--Einstein condensate modes

A quantum device for measuring two-body interactions, scalar magnetic fields and rotations is proposed using a Bose--Einstein condensate (BEC) in a ring trap. We consider an imbalanced superposition of orbital angular momentum modes with opposite winding numbers for which a rotating minimal atomic density line appears. We derive an analytical model relating the angular frequency of the minimal density line rotation to the strength of the non-linear atom-atom interactions and the difference between the populations of the counter-propagating modes. Additionally, we propose a full experimental protocol based on direct fluorescence imaging of the BEC that allows to measure all the quantities involved in the analytical model and use the system for sensing purposes.


I. INTRODUCTION
Pushing the limits of sensing technologies is one of the main challenges in modern physics, opening the door to highprecision measurements of fundamental constants as well as applications in many different areas of science. Specifically, the development of highly-sensitive compact magnetic field sensors enables from detecting extremely weak biologically relevant signals to localize geological structures or archaeological sites [1]. In this context, superconducting quantum interference devices (SQUIDs) [2,3] and atomic [4][5][6][7][8][9][10][11] and nitrogen-vacancy diamonds [12,13] magnetometers are the three main approaches that allow achieving, in a non-invasive way, unprecedented sensitivity to extremely small magnetic fields.
In particular, the extraordinary degree of control of ultracold atomic systems [14,15] makes them ideal platforms for precision measurements [16]. There are basically two types of ultracold atomic magnetometers depending on whether the magnetic field drives the internal or the external degrees of freedom of the atoms. The former are typically based on the detection of the Larmor spin precession of optically pumped atoms while the latter encode the magnetic field information in the spatial density profile of the matter wave. Atomic magnetometers with Bose-Einstein condensates (BECs) have been investigated, for instance, by using stimulated Raman transitions [17], probing separately the different internal states of a spinor BEC after free fall [18], or measuring the Larmor precession in a spinor BEC [19][20][21][22][23]. In the latter case, sensitivity can be increased by probing spin-squeezed states [24]. In [25], the possibility of taking profit of Feshbach resonances to use a two-component BEC as a magnetometer was also outlined. Ultracold atomic magnetometers based on detecting density fluctuations in a BEC due to the deformation of the trapping potential have also been demonstrated [26][27][28].
In this article, we propose to use a BEC trapped in a two-dimensional (2D) ring potential for measuring with high sensitivity non-linear interactions, scalar magnetic fields and rotations. We consider an imbalanced superposition of counter-rotating Orbital Angular Momentum (OAM) modes, whose spatial density distribution presents a minimal line. A weak two-body interaction between the atoms of the BEC leads to a rotation of the minimal atomic density line whose angular frequency is directly related to the strength of such interactions. This phenomenon is somehow reminiscent of the propagation of gray solitons, which originate in repulsively interacting BECs due to a compensation between the kinetic and mean field interaction energies. In this case, however, the minimal density line appears for attractive, repulsive or even non-interacting BECs, and is a consequence of the interference between the counterpropagating modes that takes place due to the circular geometry of the system. The rest of the paper is organized as follows. In section II we describe the physical system and we derive an analytical expression that accounts for the rotation of the line of minimal density. In section III, we take profit of this expression to propose a full experimental protocol to measure the interaction strength, which is proportional to the s-wave scattering length. Far from the resonant field or with a dilute enough BEC, the relation between the scattering length and the applied magnetic field given by Feshbach resonances could be exploited to use the system as a novel type of scalar magnetometer. We also outline the possibility of using the system as a rotation sensor. Finally, in section IV we summarize the main conclusions. In appendix A we derive the general equations that govern the dynamics of a BEC carrying OAM in a ring potential, and in appendix B we give further details about the experimental implementation of the measurement protocol.

A. Physical system
We consider a BEC formed by N atoms of mass m confined in the z direction by a harmonic potential of frequency ω z and in the perpendicular plane by an annular trap of radial frequency ω and radius R. We study the system in the limit of strong confinement along the z direction; ω z ω. Under this assumption, in the limit a z a s n 2 1, where a s is the s-wave scattering length, n 2 the two-dimensional density of the BEC and a z = /(mω z ) the harmonic oscillator length along the z direction, the three-dimensional (3D) Gross-Pitaevskii equation (GPE) can be restricted to the x − y plane by considering the profile for the BEC order parameter along the z direction as a Gaussian of width a z , which corresponds to its ground state along this direction [52]. In doing so, the 3D two-body interaction parameter g 3 = (N 4π 2 a s )/m is transformed to its two-dimensional (2D) form g 2 = (N √ 8π 2 a s )/(ma z ) (note that in these expressions we have taken the BEC wave function to be normalized to 1). Thus, the 2D GPE that we will use to describe the system reads where V (r) = 1 2 mω 2 (r − R) 2 is the potential created by the ring. Furthermore, by expressing the distances in units of σ = mω , the energies in units of ω and time in units of 1/ω, we arrive at the following dimensionless form of the 2D GPE, which is the one that we will use throughout the paper where all quantities are now expressed in terms of the above defined units and the dimensionless non-linear interaction parameter is given by The system supports stationary states with a well-defined total OAM l and positive or negative winding number, which we denote as |l, ± . The OAM eigenstates have the wave functions where f l (r) is the corresponding radial part of the wave function.

B. Dynamics in the weakly interacting regime
Let us consider as initial state an imbalanced superposition of the |1, + and |1, − states, with n 1± ≡ p 1+ − p 1− being the population imbalance. Such state could be realized for instance by preparing the BEC in the ground state of the ring, imprinting a 2π round phase and momentarily breaking the cylindrical symmetry of the potential to induce a coupling between the degenerate states of positive and negative circulation [57,58] or by directly transferring OAM with a laser beam [59].
Due to parity reasons, the non-linear term in the GPE can only couple OAM states with odd total OAM l, see appendix A for a more detailed justification. Thus, we can write the total wave function at any time t as Ψ( r, t) = l odd β=± a lβ (t)φ lβ ( r).
Since we focus on the weakly interacting regime, we consider that the only higher energetic states with a relevant role in the dynamics are |3, + and |3, − . In order to simplify the forthcoming analytical expressions, we assume that the radial part of the wave functions are the ones of the ring potential ground state, i.e. we take f l (r) = f 0 (r) in Eq. (4). This is an excellent approximation as long as the width of the the density profile of the BEC along the radial direction is much smaller than the radius of the ring, which is always the case in the weakly interacting regime. The time evolution of the probability amplitudes a l± (t) (l = 1, 3) is obtained by substituting (5) into the GPE (2) (see appendix A for details) where the four-state model (FSM) Hamiltonian reads where ρ i±j± ≡ a i± a * j± with i, j = 1, 3 are the density matrix elements, µ l the chemical potential of the l = 1, 3 OAM states, Hφ l± = µ l φ l± , and U = g 2d d 2 r|f 0 (r)| 4 . From these parameter definitions, the validity condition of the weakly interacting regime reads (µ 3 − µ 1 ) ≡ ∆ U . Within this regime, Fig. 2(a) shows a typical temporal evolution of the populations of all the OAM states involved in the dynamics considering as initial state an imbalanced superposition of the |1, + and |1, − states. The continuous lines have been obtained by solving with a high order Runge-Kutta method the FSM, Eq. (6), and the insets show the comparison with the results obtained by a full numerical integration of the 2D GPE (points). We have performed this integration using a standard Crank-Nicolson algorithm in a space-splitting scheme [60], i.e., we have introduced the Trotter decomposition e iH(x,y)∆t ≈ e iH(x)∆t e iH(x)∆t , where ∆t is the discrete time step, that we have taken to be ∆t = 10 −3 . The grid used for the simulations has a spatial discretization width ∆x = 2.4 × 10 −3 and a total of 1000 points in each dimension. For all the populations, we find an excellent agreement between the results obtained with the two different methods, with relative discrepancies typically on the order of 10 −2 . Despite the fact that the populations of the different OAM states present only very small fluctuations, the initial state is not in general a stationary state of the system because the minimum appearing in the density profile rotates at a constant speed. This fact can be appreciated in Fig. 2(b), where the density profile is shown for different times. At t = 0, the density profile has a minimum density line at x = 0, and as time marches on this line rotates in the x − y plane. The fact that the minimum density line rotates means that there is a timedependent relative phase α(t) between the a 1+ (t) and a 1− (t) coefficients, so that the state of the system evolves in This phase difference is due to the non-linear interaction, and can be understood as a consequence of the presence of off-diagonal terms in the FSM Hamiltonian (7). In order to determine the time dependence of α, in Fig. 2(c) we plot the temporal evolution of the real part of the coherence We observe that it oscillates harmonically, which means that α evolves linearly with time. The oscillation frequency of the coherence corresponds to the rotation frequency of the minimum density line.
From the FSM, we can obtain the oscillation frequency of ρ 1+1− by solving the von Neumann equation iρ = [H FSM , ρ]. After assuming ρ 1+1+ = p 1+ and ρ 1−1− = p 1− to be constant and neglecting all terms O(a 2 3± (t)), we arrive at a linear system of three coupled differential equations The characteristic frequencies k of the system of equations (8) are obtained by solving the eigenvalue equation Since U ∆ in the weakly interacting regime, the term proportional to p 1+ p 1− U 2 can be neglected in front of the others. The three eigenvalues that are obtained after solving Eq. (9) are imaginary. The eigenmode associated to the eigenvalue of lowest modulus k 0 has a predominant component of ρ 1+1− (t), allowing us to write ρ 1+1− (t) ≈ ρ 1+1− (0)e k0t . Thus, the rotation frequency of the nodal line is Ω FSM = − i 2 k 0 , where the subscript indicates that the rotation frequency has been obtained in the context of the FSM. In the limit ∆ Ω F SM , the rotation frequency of the nodal line is given by Note that, although the l = 3 states are nearly not populated during the dynamical evolution, the parameter ∆, which contains the chemical potential µ 3 , plays a significant role in the expression of the rotation frequency (10). Thus, these states must be taken into account for an accurate description of the dynamics of the system.

III. QUANTUM SENSING PROTOCOL
A. Sensing of two-body interactions Recalling that the parameter U of the FSM Hamiltonian (7) is given by U = g 2d d 2 r|f 0 (r)| 4 ≡ g 2d I and assuming that we are in the regime of validity of the FSM, Eq. (10) allows us to express the interaction parameter g 2d as where Ω is the observed frequency of rotation of the nodal line. The relation (11) constitutes the basis to use the physical system under consideration as a quantum sensing device. By determining the parameters appearing on the right hand side, one can infer the value of g 2d and thus, from Eq. (3), either the s-wave scattering length or the number of atoms forming the BEC. In Fig. 3(a), we plot Ω as a function of g 2d for different values of n 1± , computed using (10) (continuous lines) and the full numerical integration of the 2D GPE (points), showing an excellent agreement between the two methods for low non-linearities and population imbalances. For g 2d < 4, Fig. 3(b) shows the relative error δΩ ΩGPE , where Ω GPE is the rotation frequency of the nodal line obtained from the GPE and δΩ = |Ω FSM − Ω GPE |, as a function of the ab initio values of n 1± and g 2d in the numerical simulation, finding a maximum relative error of 10 −2 . Since all the treatment developed so far is valid for low values of g 2d , this sensing device could be used for dilute BECs. The rotation frequency of the minimum density line, Ω, can be measured by direct imaging in real time of the density distribution of the BEC. If the coherence time of the BEC is τ , in order for this measurement to be possible the condition Ωω > ∼ 1/τ must be fulfilled, since otherwise the rotation would be so slow that it could not be appreciated during the time that the experiment lasts. The upper limit of observable relevant values of Ω is imposed by the regime of validity of the model. If the interaction is too large, the assumptions of the FSM model are no longer valid and it is thus not possible to relate the rotation frequency of the nodal line to the non-linear interaction parameter using (11). The rest of parameters appearing on the right hand side of (11) can be determined experimentally from fluorescence images of the BEC. In Appendix A we design a specific protocol to measure the population imbalance n 1± , the integral of the radial wave function I, and the chemical potential difference ∆. Note that currently there are different approaches to measure the s−wave scattering length of ultracold atoms [61] such as those based on photoassociation spectroscopy, ballistic expansion, and collective excitations. Our proposal constitutes an alternative to these approaches where all the unknowns can be directly inferred from fluorescence images of the BEC. However, the limit g 2d < 4 obtained for the configuration discussed in Fig. 3 implies that for a BEC of, e.g., 10 4 atoms of 23 Na, with a trapping frequency ω z of a few hundreds of Hz, the maximum s−wave scattering length that could be measured with high precision, e.g., with a relative error of 10 −2 , would be few times the Bohr radius.

B. Sensing of magnetic fields
Assuming that the total number of atoms of the BEC N and the trapping frequency in the z direction ω z are precisely known quantities, Eqs. (11) and (3) together with the protocols to measure n 1± , I and ∆ allow to determine the scattering length a S at zero magnetic field. Alternatively, if the scattering length is a known quantity, the measurements of Ω, I and ∆ can be used to determine n 1± through the aforementioned relations.
If the scattering length depends somehow on the modulus of the external magnetic field B, turning on the field will be translated into a variation of Ω. Thus, the system could be used as a scalar magnetometer by relating changes on the frequency of rotation of the minimal line to variations of the modulus of the magnetic field. Taking into account that I and ∆ are almost independent of g 2d and thus of B in the regime of interaction strengths for which the model is valid, combining Eqs. (3) and (10) we can evaluate the sensitivity that this magnetometer would have as Since we must have U ∆ in order for the model to be valid, we can define a threshold limit for the sensitivity by taking U/∆ = 1 in (12). Defining the aspect ratio Λ ≡ ω z /ω and changing the differentials in (12) by finite increments, we find the following upper threshold for the sensitivity in magnetic field variations ∆B th as a function of the change in the rotation frequency of the nodal line From Eq. (13), we observe that the sensitivity is improved by having a large number of condensed particles and a strong dependence of the scattering length on the magnetic field modulus. However, since the parameter g 2d ∝ N a s needs to be small in order for the model to be valid, it is also required that the scattering length takes small values.
In the presence of a Feshbach resonance, the scattering length depends of the magnetic field modulus as whereã S is the background scattering length, B 0 is the value of B at resonance and δ is the width of the resonance. Thus, by placing the magnetic field close to the resonant value B 0 , one could in principle meet both the requirement that the scattering length is small and that it depends strongly on the magnetic field modulus. However, in most cases this procedure would have the inconvenience that close to a Feshbach resonance the three-body losses are greatly enhanced, limiting the lifetime of the BEC and hindering the measurement procedure. Nevertheless, some atomic species such as 85 Rb [62], 133 Cs [63], 39 K [64] or 7 Li [65] have been reported to form BECs that are stable across Feshbach resonances, so they could be potential candidates for using the system as a magnetometer. Additionally, the BECs formed by these species have lifetimes on the order of a few seconds. Taking into account that the trapping frequency ω, in units of which Ω is expressed, is typically of the order of a few hundreds of Hz for ring-shaped traps, and considering typical values of Ω shown in figure 3 (a), in International System units Ω ∼ 1Hz. This means that in the typical time that an experiment would last, τ ∼ 1s, the minimum density line would perform some complete laps. Under the reasonable assumption that the fluorescence imaging system could resolve angular differences on the order of ∼ 0.1 rad, incrementals in the rotation frequency on the order of 10 −2 Hz could be measured. Thus, in the dimensionless units of Eq. (13), sensitivites on the order of ∆Ω ∼ 10 −4 could be achieved. These atomic species have, however, the drawback that they typically form BECs with a low number of particles, which limits the sensitivity to magnetic fields. Although it is outside of the scope of this paper to give accurate values of the sensitivities that could be achieved with this apparatus, making use of Eqs. (13) and (14), and considering the experimental parameters reported in [54], we have estimated that, in principle, this magnetometer would allow to measure changes in the magnetic field on the order of a few pT at a bandwidth of 1 Hz.
As a last remark, we point out that after measuring the scattering length, far from the resonant field B 0 , if the line of minimal density rotates at a constant speed the relation (14) can be inverted to infer the absolute value of the magnetic field.

C. Sensing of rotations
Let us consider the case when the BEC is placed in a reference frame rotating at an angular frequency Ω ext , which is positive (negative) if the rotation is clockwise (counter-clockwise). Now the dynamics is governed by the modified GPE where L z = −i ∂ ∂ϕ is the z component of the angular momentum operator. The ideal instance for using the system under study as a sensor of rotations is the non-interacting limit g 2d = 0. In that case, it can be easily shown that the effect of the external rotation is to make the line of minimal density rotate at an angular speed Ω ext , which can be directly measured in experiments.
In the weakly interacting regime, the system under study can still be used as a sensor of external rotations. In that case, we find that the only difference in the dynamics with respect to the case when there is no external rotation is that the rotation frequency of the nodal line is shifted precisely by a quantity Ω ext . Thus, if g 2d is known and I, n 1± and ∆ are measured using the protocol provied in the appendix A, the system under consideration can be used as a sensing device for external rotations by computing the external rotation as Ω ext = Ω − Ω FSM , where Ω is the rotation frequency of the nodal line observed in the experiment and Ω FSM is given by (10).
The proposed setup constitutes an alternative to the two main lines of development of rotation sensors using ultracold atoms: the atomic-gas analogues of superconducting quantum interference devices (SQUIDs) [41][42][43][44][45][46] and the Sagnac interferometers, for a review see [48]. Gyroscopes based on the Sagnac effect measure a rotation rate relative to an inertial reference frame, based on a rotationally induced phase shift between two paths of an interferometer and the low available atomic fluxes and low effective areas are the main limiting factors of their sensitivity.

IV. CONCLUSIONS
We have studied the dynamics of an imbalanced superposition of the two degenerate counter-rotating l = 1 OAM modes of a weakly interacting BEC trapped in a 2D ring potential. We have found that the non-linear interaction induces a time-dependent phase difference between these two modes which leads to a rotation of the line of minimal atomic density of the BEC. The derived few state model provides a simple analytical dependence between the rotation frequency and the non-linear parameter which, for low non-linearities, perfectly matches with the ab initio numerical simulations. The measurement of the rotation frequency allows to use the system as a quantum sensor of two-body interactions, scalar magnetic fields and rotations. The theoretical treatment exposed in this work can also be extended to a regime of higher interactions, where higher OAM modes are excited and a myriad of new physical scenarios opens up.
We start by considering the expansion of a general state of the BEC in terms of the OAM modes where m ∈ Z is an index that corresponds to the multiplication of the indices l and β in the expression of the OAM modes of the main text (4). Like in the main text, we have assumed that the radial parts of the OAM modes correspond to the lowest vibrational state of the ring, i.e., f m (r) = f 0 (r) ∀m. Since the OAM wave functions are normalized to unity, |φ m (r, ϕ)| 2 rdrdϕ = 1, the amplitudes in the expansion (A1) fulfill the constraint m |a m (t)| 2 = 1. Substitution of the wave function (A1) into the 2D GPE (2) yields (we drop the explicit dependences on t and r) From the expression (A2), an equation of motion for each of the amplitudes can be found by multiplying both sides by φ * l and integrating over the whole 2D space where we have defined the quantity U ≡ g 2d rdrdϕ |f 0 (r)| 4 and we have taken profit of the fact that the OAM modes φ l (r, ϕ) are eigenstates of the time-independent 2D GPE with eigenvalue equal to their chemical potential µ l , i.e., ∇ 2 2 + 1 2 (r − R) 2 + g 2d |φ l (r, ϕ)| 2 φ l (r, ϕ) = µ l φ l (r, ϕ). From Eq. (A3), one can see that the term U m =m a m a * m a (l+m −m) , which appears due to the presence of the non-linear interaction term in the GPE, introduces coupling between different OAM modes. In this paper, we have considered initial states of the form Ψ(0) = ( √ p 1+ φ 1 (r, ϕ) + √ p 1− φ −1 (r, ϕ)); p 1+ + p 1− = 1. Thus, since a l (0) = 0 ∀l = 1, −1, the only higher order OAM states that will initially be directly coupled will be those in which the non-linear part of (A3) has a term such as a 1 a * −1 a 1 or any other combination of l = ±1 amplitudes. The only modes that have terms of this type are those of OAM l = ±3. Higher odd OAM modes, l = ±5, ±7, ..., are subsequently populated through coupling terms that contain combinations of amplitudes of lower odd OAM modes. However, for this particular form of the initial state, the modes with an even value of the OAM, l = 0, ±2, ..., cannot be excited because in their dynamical equations the terms in m =m a m a * m a (l+m −m) always contain at least one even OAM amplitude and, since the lowest even modes l = 0, ±2 are not directly coupled to the l = ±1 modes, none of the even modes will be populated during the time evolution. This justifies the expansion of the wave function in terms of only odd OAM modes of Eq. (5).
The FSM Eqs. (6) and (7) have been obtained directly from Eqs. (A3) by truncating the basis at the l = ±3 OAM modes and rewriting the resulting four non-linear coupled equations in a more compact matrix form.

Population imbalance
The population imbalance between the two l = 1 states can be determined from the density profile per particle at any time t, which can be obtained by fluorescence imaging. Since the wave function is given by Ψ( r, t) = f 0 (r)( √ p 1+ e iϕ + √ p 1− e −i(ϕ+Ωt) ), its density profile reads Thus, the atom density has a minimum at ϕ = π/2 − Ωt and a maximum at ϕ = −Ωt. Let us now consider the two integration regions A 1 and A 2 shown in Fig. 4(a), which are arcs of radius ρ and angle 2θ centred around the maximum and minimum of intensity, respectively. The integrals of |Ψ| 2 over A 1 and A 2 can be performed numerically and, for sufficiently small θ, they yield approximately Thus, combining (B2a) and (B2b) one can determine the product of populations as which, together with the constraint p 1+ + p 1− = 1, allows to determine the population imbalance from a fluorescence image.
The chemical potential of the angular momentum states can be decomposed into its kinetic, potential and interaction contributions. Since one can assume that the wave functions take the form φ l± ( r) = f 0 (r)e ±ilϕ , the potential and interaction contributions will be the same regardless of l, while the kinetic contribution is given by Thus, the chemical potential difference is only due to the difference in the centrifugal terms of the kinetic energy From Eq. (B1), one can see that the integral (B7) can be numerically performed after determining f 2 0 (r) from a fluorescence image as f 2 0 (r) = |Ψ(r,ϕ=−Ωt,t)| 2

1+2
√ p1+p1− . In order to check the accuracy of the proposed experimental protocol, we have computed g 2d using Eq. (11) and determining all the parameters on the right hand side following the above described numerical procedures, and later on comparing with the ab initio used value of g 2d in the simulation. In Fig. 4(b) we plot the relative error δg 2d g 2d committed as a function of the ab initio values of g 2d and p 1+ . In the region g 2d ≈ 1 and n 1± ≈ 0.6, the relative error is minimal and it reaches very low values, on the order of 10 −5 . The maximum value of the relative error is about 10%, and is found for low values of n 1± . In our simulations, we have used a grid of dimensions 24 × 24 and 1000 points in each spatial direction. With higher grid precision, the relative error committed with the proposed protocol could prove to be even lower.