Realizing and optimizing an atomtronic SQUID

We demonstrate how a toroidal Bose-Einstein condensate with a movable barrier can be used to realize an atomtronic SQUID. The magnitude of the barrier height, which creates the analogue of an SNS junction, is of crucial importance, as well as its ramp-up and -down protocol. For too low of a barrier, the relaxation of the system is dynamically suppressed, due to the small rate of phase slips at the barrier. For a higher barrier, the phase coherence across the barrier is suppressed due to thermal fluctuations, which are included in our Truncated Wigner approach. Furthermore, we show that the ramp-up protocol of the barrier can be improved by ramping up its height first, and its velocity after that. This protocol can be further improved by optimizing the ramp-up and ramp-down time scales, which is of direct practical relevance for on-going experimental realizations.


I. INTRODUCTION
The advancement of cold atom technology, and the level of control that can be achieved in such systems, has motivated the question if it can be used to emulate electronic circuitry, and possibly move beyond its features, Ref. [1]. While a realization of, say, the equivalent of electrons moving in a semiconducting material is an interesting direction in itself, it is particularly intriguing to capitalize on the specific features of cold atom systems, such as long-range phase coherence in Bose-Einstein condensates. This motivates to realize systems inspired by superconducting circuitry. Experimentally, an interesting starting point, and a remarkable achievement of its own, is the realization of Bose-Einstein condensates (BECs) in toroidal geometries, Refs. [2][3][4][5].
In this paper, we study the equivalent of an electronic SQUID. The condensate wave function is the equivalent of the superconducting wave function, and a potential barrier, at which the condensate density is suppressed, replaces the SNS interface. This barrier is then moved at a constant speed, which imitates a non-zero magnetic flux through the ring. This setup can also be seen as a stirring experiment, testing the superfluid properties of condensates, see Refs. [6,7]. Other theoretical studies of the ring geometry were reported in Refs. [8].
We demonstrate that a regime of a controlled and effective realization of an atomtronic SQUID exists for sufficiently large potential barrier heights, for realistic temperatures. For small barrier heights, the dynamical relaxation of the condensate to the ground state phase winding is suppressed, because phase slips at the barrier occur only at a very small rate. For larger barrier heights, the phase coherence across the barrier is suppressed because of thermal fluctuations in the bulk of the ring. This results in a smoothed-out response of the rotation, approaching a linear dependence on the stirring velocity, rather than a quantized, step-like response, that is characteristic for a SQUID. Furthermore, we consider two seemingly similar ramp-up processes for the SQUID operation: to either first ramp up the barrier to full speed and then ramp up the barrier height, or use the reverse  order. Interestingly, we find that the latter results in noticeably less heating, and that the new ground state is reached more quickly. We then discuss how this protocol can be further improved by different choices for the ramp-up and -down time scales. This paper is organized as follows: In Sect. II we de-scribe our simulation method, in Sect. III we discuss the phase slip dynamics. In Sect. IV we compare the different barrier ramp-up scenarios, and in Sect. V we conclude.

II. SIMULATION METHOD
Because the SQUID dynamics is dominated by the dynamics at the barrier, and because of the low density at the barrier, and therefore the low mean-field energy, it is imperative to include thermal fluctuations of the system. We include both thermal fluctuations, and the lowest order of quantum fluctuations, within a Truncated Wigner approximation (TWA), see e.g. [9]. The approach that we use is closely related to the one of Ref. [10]. We describe the system with the Hamiltonian where m is the atomic mass and g is the interaction strength, for which we use the approximation g = 4πa s 2 /m, with a s being the s-wave scattering length. The external potential, V (r, t) = V tr (r) + V bar (r, t) consists of the trapping potential, V tr (r) = , with θ being the azimuthal angle. ρ 0 is the radius to the ring. Within the TWA, the operators,Ψ are replaced by classical fields, which are propagated according to the equations of motion of this Hamiltonian. The initial condition of these classical fields are generated from the Wigner distribution of the initial state.
In order to carry out the calculations, we discretize the real-space description by introducing a lattice approximation and work in cylindrical coordinates. We replace the continuous wave function, ψ(ρ, θ, z) by a discrete wavefunction,ψ ijk =ψ(ρ i , θ j , z k ), with the mapping , and ρ 0 = lN θ /(2π) is the radius of the ring at the trap minimum. We emphasize that the discretization length is not constant, but increases linearly with increasing distance from the center axis of the ring. Therefore, the curvature of the ring geometry is fully taken into account. In this representation, |ψ ijk | 2 corresponds to the number of atoms per unit cell, where the volume of the unit cell is l 3 ρ i /ρ 0 . For the lattice size chosen here, the volume of the unit cell varies from 0.6l 3 for i = 0 to 1.4l 3 for i = 16.
In this representation, the equations of motion take the form where J = 2 /(2ml 2 ) is the tunneling energy, and the interaction term is given by On the lattice, the external potential is given by As mentioned above, we initialize the dynamics by sampling from the initial Wigner distribution. For the initial state, we choose a non-interacting Bose gas in the toroidal trap, of zero temperature. After the initialization, the interactions are turned on slowly to generate the desired interacting ensemble in the trap [12]. This process of turning on the interaction results in a nonzero temperature that is comparable or larger than the mean-field energy of the system [13]. This temperature is measured by weakly coupling harmonic oscillators to the current as described in Ref. [10]. For examples discussed in this paper, the temperature of the atomic cloud after initialization is approximately 4.0 J = 43 nK [14].
Throughout this paper, we use a lattice with the dimensions N ρ = 17, N θ = 126, and N z = 5, and N = 50000 atoms. The trapping frequencies are ω ρ /J = 0.5 and ω z /J = 2.5. We propagate and average over 72 initial states and the length scale of the barrier is l b /l = 3. We set U 0 /J = 0.07, which corresponds to a length scale l = 0.99µm, a time scale ∆t = /J = 0.71 ms and an energy scale J = 10.8 nK for sodium atoms. The healing length in the bulk of the system is ξ ≈ 0.9µm. This length is small compared to system size in the radial direction and comparable to the system size in the z-direction, which puts the system in the dimensional cross-over regime between two and three dimensions.
We now consider the following experiment: Starting at t 1 = 100∆t, we ramp the barrier magnitude from 0 to the final barrier height V b over a time period of 200∆t, while keeping the barrier stationary. Then at t 2 = 300∆t, the barrier is accelerated to its final stirring frequency, ω s over 200∆t. At t 3 = 500∆t, the atomic cloud is stirred at constant frequency for a time period of 1200∆t at maximum barrier height. At t 4 = 1700∆t the barrier height is ramped down over the time 200∆t while continuing to stir the atoms at the frequency ω s . We refer to this procedure as protocol 1.
Additionally, we compare this protocol to the process of turning on the barrier height to V b , while stirring  at constant frequency ω s , starting at t 1 and reaching the maximum barrier height at t 2 . The atomic cloud is stirred for 1400∆t before ramping down the barrier, while continuing to stir at constant frequency. This second protocol is reminiscent of the experiments reported in Ref. [2]. As we discuss below, the first protocol is the preferable protocol for large stirring frequencies, because it avoids exciting phase slips before the barrier has reached its maximum height. We elaborate on this protocol further in Sect. IV by considering different ramping time scales.
In Fig. 1, we illustrate the dynamics of the system by depicting the density n(ρ, θ) in the z = 0 layer, and the current j θ (ρ, θ) along the azimuthal direction of the ring, at four different times, for a barrier height of V b /J = 1.8 and stirring frequency ω s /ω 0 = 3.8 for the first protocol, with ω 0 = /(mρ 2 0 ). At the first time, t = 300∆t, the trap potential is fully ramped up, as is visible in the density depletion of the ring shaped condensate, but is still stationary. Here, the phase winding is still zero, and both positive and negative current fluctuations are visible. We note again that these fluctuations are predominantly of thermal origin. At the second time, t = 410∆t, one phase slip has occurred. Now the current has acquired a preferred direction, as is immediately visible in Fig. 1 (d). At the third time, t = 1000∆t, a second phase slip has occurred, and the magnitude of the current has increased, Fig. 1 (f). Finally, at time t = 2000∆t, the barrier has been ramped down, so the density of the condensate does not display a depleted region. However, as is visible from the current in Fig. 1 (h), the stirring of this ring shaped condensate has imparted a finite current circulating in the ring.
For a system with a complex order parameter, such as a condensate of atoms or a superconductor, this current is related to the well-defined phase. In the condensed phase, this results in a quantization of the phase winding around a ring geometry. We determine the phase along the central line along the ring, which tracks the maximal density for a given azimuthal angle, and in the z = 0 plane. We calculate the phase at the angle θ via The phase winding is the sum of the phase differences around the ring, n ∆φ = ( θ δφ(θ)) /2π, where the phase difference between two points, is between −π and π. The phase winding is calculated in each individual realization and then averaged over the realizations to generate the average phase winding.

III. PHASE SLIP DYNAMICS
In Fig. 2 we show the time evolution of the average phase winding along the center of the toroidal trap for three different barrier heights and three different stirring frequencies. For the smallest barrier height of V b /J = 1, and for all three rotation frequencies ω s , the phase winding does not relax to the new ground state. Instead, on the time scales of a typical experiment, it remains in a long-lived, metastable state. On the other hand, for large barrier heights such as V b /J = 2.4, the coupling across the barrier is small compared to the temperature, so that the dynamics are oscillatory and noisy for the length of the experiment. Here, the dynamics of the experiment results in a distribution of phase windings, rather than -essentially -a single phase winding number. However, we observe that for intermediate barrier heights such as V b /J = 1.8, the oscillations damp out in a realistic time.
We note that the magnitude the phase slip rate, and its dependence on the system properties war discussed in Ref. [10]. In particular, it was discussed that the phase slip rate is consistent with the scaling τ −1 ph ∼ exp(−E b /k B T ), i.e. an Arrhenius law. The energy scale E b is controlled by the barrier height and width. This strong, exponential dependence is reflected in the behavior we observe here, where the relaxation dynamics of the system is nearly suppressed at V b /J = 1.0, while at around V b /J ≈ 2.4 the barrier is so high that any coherence across the barrier is suppressed.
To understand the origin of the oscillatory behavior we show the time evolution of the density along the ring. We define radially projected density n 1D (θ) = dρ dzn(ρ, θ, z). We emphasize that while this quantity has the dimensions of a one-dimensional density we do not imply that the dynamics can be reduced to that of a one-dimensional system. In fact, as we had discussed in Ref. [10], the phase slips that occur in the system are due to a non-trivial process of a vortex traversing the barrier region. We merely use this integrated density as a convenient way to depict the system evolution.
In Figure 3 (a) and (c ) the time evolution of the radially projected density is shown for V b /J = 1.8 and V b /J = 2.4, respectively, and for ω s /ω 0 = 3.4. In (b) and (d), the same evolution is shown in the reference frame of the barrier. As is visible in this figure, the oscillations in the average phase winding are due phonon pulses generated during the initialization and acceleration of the barrier. The density waves observed in the density travel at the speed of sound, which is approximately v s = 1.52l/∆t or 2.1mm/s. With this velocity, the period of the oscillations of the average phase winding is of the order of 2T = 166∆t, corresponding to traveling back and forth along the circumference of the ring. This matches the period of the oscillations observed in Fig. 2. We observe that the oscillations for V b /J = 1.8 damp out during the time of the simulation, whereas for V b /J = 2.4 they do not. This is due to the reflection of the phonon pulse at the barrier. For the high barrier this reflection is essentially complete, and the pulse travels back and forth with only little damping. For the intermediate barrier height, the reflection is partial, which results in dephasing and damping.
Motivated by this observation, we suggest ways to minimize the detrimental effect of these phonons pulses and the resulting oscillations in the phase winding in Sect. IV, as a step towards improving the SQUID operation of the condensate ring.
The phase winding that emerges from this dynamical evolution is shown in Fig. 4. We depict the final, average phase winding after the barrier has been ramped down, as a function of the stirring frequency for several barrier heights, in the main panel. As indicated, the idealized, fully quantized, and fully relaxed phase winding shows steps at ω s /ω 0 = n + 1/2, with n being an integer. This is shown as a black line. For too small of a barrier height, such as V b /J = 1, the system remains dynamically trapped at a smaller phase winding than for the fully relaxed system. As the barrier height is increased, the behavior of the phase winding approaches a smoothed-out step-like behavior, for V b /J = 1.8 and V b /J = 2.4. We note that this smoothed out step-like response in the phase winding is comparable to the results that were obtained experimentally in Refs. [2]. In the inset of Fig. 4, we plot the average phase winding, timeaveraged over four phonon periods, [t 4 − 4T, t 4 ], before the barrier is ramped down. This is the phase winding that the system relaxes to at longer times, after the oscillatory behavior has damped down, with the stirring on. We note that for V b /J = 1.0 and V b /J = 1.8 the timeaveraged phase winding is very close to the phase winding shown in the main panel, because the oscillatory behavior has damped out on the time scale of the experiment. However, for V b /J = 2.4, the phase winding that is shown in the inset is smoothed out to an almost linear behavior. This indicates that the classical limit of this response is nearly reached, as expected for a fully disconnected ring. We therefore conclude that this limit should be observed in experiment for high barrier potentials, at long stirring times and instantaneous ramp-down of the barrier. Furthermore, we conclude that the re-emergence of the step-like behavior after the ramp-down of the barrier, is due to the non-zero time of this ramp-down, during which the system develops a well-defined phase winding. This motivates our proposal to increase this ramp-down time in Sect. IV. In Fig. 5 (a) and (b) we elaborate on the behavior shown in Fig. 4. We show the distribution of the phase winding, as a function of time, which goes beyond the expectation value of this distribution that was shown in Fig. 4. In Fig. 5 (a)  We note that the width of this distribution is controlled by the long range phase fluctuations along the quasi-1D geometry of the ring-shaped condensate. As discussed in Ref. [11], the single particle correlation func-tion along the ring falls off exponentially, with a length scale l φ = 2 N 0 /(πmρ 0 k B T ). N 0 is the number of condensed atoms. This length scale has to be compared to the circumference of the ring, L c = 2πρ 0 ≈ 124.4µm. For N 0 ≈ N , and N being the total atom number, and for T = 50 nK, we have l φ = 336µm. This results in a ratio L c /l φ that is smaller than 1. However, this does suggest, that the phase coherence across the barrier can be further stabilized by increasing this ratio. This ratio can also be written as where λ T is the thermal de Broglie wavelength, λ T = (2π 2 /(mk B T )) 1/2 , and n 2D = N/(πρ 2 0 ) is a the density of a hypothetical system of area πρ 2 0 with N atoms. Another way of stating the meaning of this ratio is that it describes the magnitude of the phase difference across the barrier, (∆φ) 2 ∼ L c /l φ . Therefore, the figure of merit that determines if the elongated 3D condensate is phase coherent along its extended 1D axis, is of the form of an inverse 2D phase space density. The optimal regime is that of low temperatures, and ring condensate with a small radius and high density. If, on the other hand, one wants to explore the effective 1D regime of a phase-fluctuating condensate, this figure of merit has to be increased above 1.

IV. OPTIMIZING THE BARRIER PROTOCOL
As a last point, we compare the two ramp-up protocols, mentioned above, as well as different ramping times for the barrier. In Fig. 5 (b) we show the case in which the barrier height is ramped up first, while remaining stationary, and then its velocity. This is the case which has been discussed thoughout this paper. In Fig. 5 (c) we show the case in which the velocity is always at its final magnitude, and the height is ramped up from zero. We see that in the latter case the response is visibly more oscillatory and noisy. Here, the phase slips begin to occurs around t = 200∆t, when the barrier is still well below its maximum height. Each phase slip generates phononic excitations which are released into the bulk of the condensate, as visible in Fig. 5 (c). Similar processes were observed in Ref. [10]. Additionally, we observe that it takes longer for the system to relax to a single phase winding. This suggests that the density at the barrier should be minimized when the phase slips occur to reduce undesirable excitations. This is illustrated in Fig. 6. Here, we show the total energy per particle during the evolution of these two protocols. A larger magnitude of this quantity will result in a higher temperature of the system after it has thermalized, and therefore can be used as a measure for how well the SQUID is implemented. As visible, in this protocol additional undesired excitations are created, along with the desired phase slips, which leads to longer relaxation  , the barrier is stationary while it is ramped up, from t1 = 100∆t to t2 = 300∆t and then accelerated at maximum height, from t2 = 300∆t to t3 = 500∆t (protocol 1). In (c), the barrier stirs with a constant frequency while the magnitude is ramped up, from t1 = 100∆t to t2 = 300∆t (protocol 2).
times and additional heating of the system. These effects are more pronounced at higher stirring frequencies than at lower stirring frequencies and are expected to play a bigger role as the stirring frequency is further increased. Again, we observe that the preferable operation of the SQUID consists of first ramping up the barrier height, and then the barrier velocity. Next we address the influence of the ramp times on the SQUID operation. As a key example, we focus on the first step of the phase winding, around ω s /ω 0 ≈ 1/2. In Fig. 7, we consider three different time sequences. The sequence that we discussed up to here, is shown in panel (c ). Additionally, we consider two other sequences, shown in (d) and (e). The sequence in (d) features a slow ramp down, and the sequence in (e) features both a slow ramp-up and a slow ramp down. For stirring frequencies near the first step, these two sequences both lead to an improvement: In panel (f) we show the resulting phase winding. The phase winding increases more steeply for the time   Figure 6: We depict the dynamics of the energy for sweep frequency ωs/ω0 = 3.8 and barrier height, V b /J = 1.8. In protocol 1 the barrier is stationary while it is ramped up, from t1 = 100∆t to t2 = 300∆t and then accelerated at maximum height, from t2 = 300∆t to t3 = 500∆t. In protocol 2, the barrier stirs at constant frequency while the magnitude is ramped up, from t1 = 100∆t to t2 = 300∆t. E0 is the total energy at t=0.
sequences with slower ramp-up and ramp-down. Furthermore, in panel (g) we show the resulting increase of the energy of the system, which is again improved for the slower ramp times. For these comparatively small stirring frequencies, the improvement is primarily due to the slower ramp-down. For higher stirring frequencies, the slower ramp-up time leads to a further improvement of the operation, because the phonon pulse that is created by the barrier acceleration is reduced. We give an indication for this behavior in panel (b). Among the three phase winding evolutions that are shown, the protocol that is shown in panel (e) is the least oscillatory. Furthermore, we point out that the slow ramp-down time also allows for the sloshing motion of the system to damp down. The detrimental feature of this motion is not due to the magnitude of the phase fluctuations, but because it can the lead the system to arrive at a phase winding other than the ground state one, if the barrier is ramped down fast. As mentioned above, this phonon motion damps out for intermediate values of the barrier, but not for high barriers. If the ramp-down is done sufficiently slow, this provides enough time for the oscillations to damp out, while the barrier is at intermediate values.

V. CONCLUSIONS
In conclusion, we have demonstrated the atomtronic implementation of a SQUID in a toroidal Bose-Einstein condensate, with particular emphasis on the effect of thermal fluctuations. These are included in our Truncated Wigner approximation of a realistic system, as realized in the experiments of Refs. [2,3]. We show that the regime of SQUID operation is viable for sufficiently large barrier heights, which imitates an SNS junction of a solid state SQUID. For too low of a barrier height, the rate of phase slips is too low for the system to reach the equilibrium phase winding number for a given moving barrier speed. For larger barrier heights, the thermal phase fluctuations suppress the coherence across the barrier. The characteristic, step-like behavior of the phase winding is achieved for either intermediate values of the barrier height, or during the ramp-down of the barrier, if this occurs on a sufficiently long time scale. Furthermore, we investigate two ramp-up protocols of the barrier height and the barrier velocity, and show that ramping up the barrier height first, before setting it in motion, results in less heating. An additional improvement can be achieved by increasing both the ramp-up time of the barrier and the ramp-down time. The slower ramp-up results in a reduction of the phonon pulse that is emitted when the barrier is set in motion. The slower ramp-down of the barrier improves the reemergence of an integer phase winding after the stirring. We emphasize that these results and considerations will similarly apply to all atomtronic circuits, that are based on condensate dynamics in non-trivial trap geometries, and are therefore of broad interest to the emerging field of imitating superconducting circuitry with Bose-Einstein condensates.