Implementation of an atomtronic SQUID in a strongly confined toroidal condensate

We investigate the dynamics of an atomtronic SQUID created by two mobile barriers, moving at two different, constant velocities in a quasi-1D toroidal condensate. We implement a multi-band truncated Wigner approximation numerically, to demonstrate the functionality of a SQUID reflected in the oscillatory voltage-flux dependence. The relative velocity of the two barriers results in a chemical potential imbalance analogous to a voltage in an electronic system. The average velocity of the two barriers corresponds to a rotation of the condensate, analogous to a magnetic flux. We demonstrate that the voltage equivalent shows characteristic flux-dependent oscillations. We point out the parameter regime of barrier heights and relaxation times for the phase slip dynamics, resulting in a realistic protocol for atomtronic SQUID operation.

An important and compelling phenomenon is the quantum interference of currents that flow through two Josephson junctions. When they are brought to interference the resulting current gives access to the accumulated relative phase. The relative phase is due to the magnetic flux enclosed by the currents. These devices couple to the magnetic field via the electric charge of the Cooper pairs. This was first observed in superconducting quantum interference devices (SQUIDs) [33]. SQUIDs have applications in quantum sensing [34] and information processing [35].
The field of atomtronics aims to implement analogues of electronic devices such as SQUIDs. Superflow in toroidal atom condensates allows us to realize atom analogs of dc-and rf-SQUIDS using a weak link [21,36], and dc-SQUIDs using two weak links [37,38]. Recently, an atomic dc-SQUID was realized in a voltage state at a constant bias current [39], demonstrating the quantum interference of currents.
In this paper we propose an experimental setup to measure the voltage-flux dependence of a dc-SQUID, which is created by two mobile barriers in a quasi-1D ring condensate. We study the system dynamics using a multiband simulation approach suitable for toroidal condensates in quasi-1D setting. Within a truncated Wigner approximation, we take quantum fluctuations into ac-count and extend the theoretical description beyond the Gross-Pitaevskii equation [39]. We simulate the presence of a magnetic field by letting the barriers rotate around the ring. This motion excites phase slips in the condensate, similar to how a penetrating flux induces a phase winding in a superconducting loop. We show that the average phase winding depends on the barrier height and results in a step-like flux dependence for large barrier heights. We then superimpose another barrier motion which lets the barriers approach each other with a nonzero relative velocity. This simulates a bias current flowing through the SQUID. Operating the SQUID with a bias current results in the formation of a density imbalance, which is analogous to a voltage in an electronic circuit and is a distinct feature of dc-SQUIDs in the resistive regime. We analyze the dependence of the voltage on the applied magnetic flux and find characteristic oscillations at the periodicity of the flux quantum. This demonstrates quantum interference effects and is analogous to the voltage-flux relation based on electronic SQUID models. Our work builds on a previous study on an atomtronic SQUID in a 3D toroidal setting [40] and experimental work utilizing a mobile barrier based SQUID without magnetic flux [38].
This paper is organized as follows. In Sec. II we describe the multi-band simulation approach to simulate the dynamics of a quasi-1D condensate. In Sec. III we implement the barrier protocol and discuss the operation of an atomtronic dc-SQUID. In Sec. IV we discuss the emergence of phase-slips due to the barrier motion. In Sec. V we analyze the voltage-flux oscillations. We conclude in Sec. VI.

II. SIMULATION METHOD
We consider a lithium gas in a toroidal trap, in which the lithium atoms are in molecular form of two bound 6 Li atoms. The trapping potential is harmonic along the transverse directions with the trapping frequencies ω ρ and ω z , see below. The trapping energies are larger arXiv:2204.03000v1 [cond-mat.quant-gas] 6 Apr 2022 than the mean-field energy µ of the gas, i.e. ω ρ,z > µ. Therefore, the system is in a quasi-1D regime. We emphasize that we do not approximate the gas as a purely-1D system, via a single-band approximation, but rather include transverse motion by including several of the lowest bands, resulting in a quasi-1D description, as we expand on below. Spatially, the condensate is confined to several µm in the transverse directions. The underlying Hamiltonian of a 3D interacting BEC iŝ where g = 4πa s 2 /m D is the interaction strength proportional to the s-wave scattering length a s , and m D is the molecular mass, i.e. m D = 2m Li . The potential V (r) is a harmonic ring trapping potential V trap (r) = m D ω 2 ρ (ρ − R) 2 + ω 2 z z 2 /2. We set the ring circumference to 250 µm, so the radius is R ring ≈ 40 µm. The trapping frequencies are ω ρ = 2π × 4 kHz and ω z = 2π × 1.5 kHz. For our numerical calculation, we discretize the spatial motion along the azimuthal direction on a 1D lattice of 250 sites with discretization length l = 1 µm [41]. We work in cylindrical coordinates (ρ, θ, z). The transversal degrees of freedom ρ and z are treated by expanding the field operatorψ for each site i aŝ where φ are the eigenfunctions of the quantum harmonic oscillator. The quantum numbers m = (m 1 , m 2 ) are written as a tuple for compact notation. We restrict ourselves to 0 ≤ m 1 , m 2 ≤ 2, i.e. we include two excited states in each transverse direction. With this expansion in higher bands in the transverse directions, the system is modelled as a multi-band Hubbard model. For our numerical simulations we use the truncated Wigner approximation and therefore approximate the quantum dynamics by a semiclassical evolution [42,43]. This allows us to replace the operatorsψ by complex numbers ψ. We sample the initial state from a continuous Wigner function, which has the form of a Gaussian for a non-interacting system [43]. Each state of the initial ensemble propagates under the classical equations of motion where J = 2 /(2m D l 2 ) is the tunneling energy and V b,k,n,i correspond to the external barrier potential V b that we introduce in Sec. III. The coefficients for the onsite interaction and the barrier potential are defined as and where the basis functions φ j are taken to be real. These coefficients are given in Appendix A.
As described above, we initialize the dynamics by sampling the Wigner function of a non-interacting system at zero temperature. The ground state is populated by a large number of particles N = 40000. All excited states only contain quantum fluctuations leading formally on average to an increase of 1/2 particles per state or lattice site. The interaction is then slowly increased to the desired value of U/J = 0.01 over 3.8 s. This results in heating of the condensate to a temperature comparable to the mean-field energy. The resulting state is the initial state of the simulation of the physical system, for which we set the initial time to t = 0. We add two mobile barriers to create an atomtronic SQUID, as we describe in Sec. III. As our observables, we calculate the local density n(θ) = |ψ(θ)| 2 and the global phase winding Ω w = θ δφ 0 (θ), where the phase difference δφ 0 (θ) = φ 0 (θ + l) − φ 0 (θ) is between −π and π. φ 0 (θ) is the local phase corresponding to the lowest (condensate) mode. We calculate these observables for each sample and then average over the initial ensemble to obtain the average density and the average phase winding.

III. IMPLEMENTATION OF AN ATOMTRONIC SQUID
A solid state realization of a dc-SQUID consists of two Josephson junctions in a superconducting ring. In the following, we implement such a device as an atomtronic system in the so-called voltage state, where a constant bias FIG. 2. SQUID dynamics in a ring condensate. Column density n(ρ, θ) = |ψ(ρ, θ)| 2 at different times during the barrier protocol with the barrier strength V0/µ = 1.7. At earlier times (a-c), the two barriers move with the same velocity v φ by choosing v1,2 = 530 µm/s and v bias = 0, which models a magnetically induced screening current. The arrows indicate the barrier velocities and direction of motion. The stirring motion continues for 2.2 s, which allows the condensate to relax in the phase winding state of the lowest energy. At later times times (d-f), the barriers move at different velocities v1,2 = ±v bias /2+vΦ, with a bias velocity of v bias = 106 µm/s, and thereby start to approach each other. The strong barrier results in a resistive regime and a density imbalance is formed, in analogy to a voltage. current I bias is applied to the setup. Furthermore, we apply the equivalent of a magnetic flux Φ through the ring. As we demonstrate below, the voltage shows characteristic oscillations depending on the magnetic flux that is applied to the SQUID, in analogy to a solid state SQUID. The required current regime derives from the critical current I c of a Josephson junction, at which the junction switches from the non-resistive to the resistive regime. Given the two Josephson junctions in the SQUID, the minimal bias current is 2I c , for a nonzero voltage across the SQUID.
The existence of these voltage-flux oscillations is due to the fact that the flux can only pass through the ring in quantized multiples of the flux quantum Φ 0 . If the magnetic flux is not an integer multiple of Φ 0 , i.e. Φ = nΦ 0 (n = 0, 1, 2, ..), a screening current I s is formed in the ring to enhance or reduce the flux to an integer multiple of Φ 0 . Therefore, the total current at the two junctions is modified to I junction = I bias /2 ± I s . This leads to a lowered critical current for higher |I s | which in turn increases the voltage. So the voltage-flux curve oscillates with a periodicity of Φ 0 [44].
An atomic analogue of solid-state SQUID is experimentally realized by Ref. [39]. Here we propose to measure voltage-flux oscillations as follows. We use two barrier potentials, of the form: to model the Josephson junctions of width σ = 0.9 µm. θ b1 and θ b2 are the angular position of the two bar-riers. Unlike the electronic SQUID where the barriers are stationary, we consider mobile barriers, i.e. θ b1,b2 = θ b1,b2 (t). Instead of letting a current I bias flow through the junction, we move the barriers with a velocity ±v bias /2. Note that this scenario of equal velocity barriers has been experimentally realized by Ref. [38].
To model the magnetic flux we add the velocity v Φ , so that the two velocities are v 1,2 = ±v bias /2 + v Φ .
The barrier protocol is described in Fig. 1. We gradually ramp up the barrier strength V 0 over 75 ms. Inspired by the setup in Ref. [38] we set the initial distance between the barriers to x 0 = 100 µm, i.e. a value less than half the circumference. This provides more space for the barrier movement later on. Then we proceed to slowly accelerate both barriers to v Φ by choosing the same v 1,2 and v bias = 0, as depicted in Fig. 1. At this stage the barriers move in the same direction with equal velocity, which we show in the simulations in Figs. 2(a-c). The distance between the barriers remains constant. Due to this stirring motion, a phase winding is formed in the condensate, indicating rotary motion. As we elaborate on in the next section, the relaxation of the condensate to the phase winding state of lowest energy, without suppressing the phase coherence entirely, is slow for typical experimental time scales. Therefore, we continue the stirring motion for 2.2 s for the system to relax. We note that in a typical solid state SQUID, the relaxation time to the current state with the lowest energy is fast compared to the time scales of the operation. In an atomtronic SQUID these time scales are in general not strongly separated, which provides an intrinsic challenge of atomtronics to emulate solid state devices.
The final step is to slowly increase v bias to its final value over the course of 150 ms, as indicated in Fig. 1. At this point the barriers move with unequal velocities v = ±v bias /2+v Φ as can be seen in the simulation results in Figs. 2(d-f). We show the example v bias /2 < v Φ so that both barriers still move in the same direction, but with different velocities. As a result, a density difference between the two subsystems can emerge, depending on the regime of operation. Within the superfluid regime of the Josephson junctions, the atoms tunnel through the barriers and the densities are equal in the two subsystems. However, in the resistive regime the tunnel current is not sufficiently large to maintain a density balance. Hence, the condensate density increases on the side of the barriers on which the two barriers approach each other and a density difference emerges in the system [ Fig. 2(f)], analogous to a voltage in an electronic SQUID.

IV. PHASE SLIP DYNAMICS
In an electronic SQUID, the magnetic flux penetrates the ring in integer multiples of nΦ 0 , with n being an integer. The wave function in the ring acquires a phase winding of 2πn accordingly. As the SQUID relaxes to its state of lowest energy, phase slips occur at Φ = (n + 1/2)Φ 0 , because this minimises the flux provided by the screening current to ≤ Φ 0 /2. We demonstrate the analogous behaviour for the atomtronic SQUID. Here, the relative velocity between condensate and barriers is minimised by the rotation of the condensate. Every 2π of phase winding is associated with a rotation velocity v 0 = 2π /(m D R ring ). This velocity v 0 plays the role of a flux quantum in this system. So the phase slips occur at v Φ = (n + 1/2)v 0 , if the condensate relaxes to the lowest energy state. We investigate this phase slip dynamics for v bias = 0 and nonzero v Φ as in Figs. 2(a-c). We calculate the average phase winding Ω w as described in Sec. II and analyze its dependence on the stirring time t and barrier height V 0 . In Fig. 3(a) we show the average phase winding in the system at t = 1.9 s, for different values of V 0 . As a dashed line, we show the phase winding increasing in unit steps, at half-integers of the flux quantum, as described above. The average phase winding that has accumulated in the condensate at that time deviates from the idealized, equilibrium expectation (dashed line). For low barriers (V 0 /µ = 0.9), no phase slips are observed for flux velocities up to v Φ /v 0 ∼ 3.5. For v Φ /v 0 4, the dependence of Ω w on v Φ is approximately linear, rather than an approximate step-like behavior. So the condensate does not relax to the rotational state of lowest energy, and SQUID operation is not possible. The origin of the slow relaxation, and the condensate remaining in a metastable state, is that the critical velocity at the barrier is too high. To emulate an electronic SQUID, the condensate has to relax to the lowest energy phase winding, to be in a resistive state for phase slips to occur. For that, we modify the barrier height. Indeed for V 0 /µ = 1.3 the critical velocity is lower and the steps are more pronounced while not sufficiently so to be a SQUID. A good step-like behaviour is observed for V 0 /µ = 1.7, where the critical velocity is reduced to a value similar to v Φ /v 0 ∼ 0.5. Note that the relaxation is generally more effective near the center of a step, i.e. at whole integer values, where screening currents are relatively low. At half integer values there are two ground states, one with v s = v 0 /2 and another one with v s = −v 0 /2 plus a phase winding. These states have nearly equal energy at this rotation velocity, and the tendency of the system to relax is suppressed. We note that such barrier-height induced deviations for phase slips were also studied in a single-barrier based atomic SQUID [40] and the suppression of the phonon velocity due to the barrier height was pointed out in toroidal condensates [22], which is similar to the reduction of the critical velocities in stirred BECs [45,46].
To illustrate the relaxation process, we consider the optimal barrier height V 0 /µ = 1.7 and analyze the time evolution of Ω w (t) at different stirring times. The results are shown in Fig. 3(b), where we focus on the first phase slip at v Φ /v 0 = 0.5. The behaviour for different barrier heights is similar for larger v Φ as we saw in Fig. 3(a). At t = 0.38 s the system is not relaxed as the step-like behaviour is not visible. At t = 0.76 s, we see the emergence of a clear plateau in the phase. Finally, at t = 1.9 s the phase is reasonably converged and a clear plateau is visible at 2π.

V. VOLTAGE-FLUX RELATION
As described earlier, we propose to induce a density imbalance in the ring by the moving barriers. This density imbalance, which corresponds to a chemical potential difference, is analogous to the voltage that emerges in a SQUID within the resistive regime. We now discuss this imbalance in more detail and its dependence on the rotational flux for V 0 /µ = 1.7. We determine the density imbalance ∆n = n R − n L , where n R and n L are defined as the density of the initially larger and smaller segment of the ring, respectively, see Sec. II and Fig. 2(a). In Fig. 4(a) we show the time evolution of ∆n(t) at v Φ /v 0 = 0.5, 0.7 and 1. Early in the dynamics, small density oscillations are exited by the motion of the barrier. These are the equivalent of plasma oscillations [47,48] and create some undesirable noise in the system. We note that these oscillations are observed experimentally as well [38]. However, with a relative amplitude of ∼ 2.5 % they are more significant than the ones observed in the experiment, where ∆n/n ∼ 1 %. Over 2 s they damp out significantly, and to fluctuations with relative amplitudes of ∆n/n ∼ 1 %.
In the absence of a bias current, the average imbalance of the condensate is zero, as shown in the time evolution before t ∼ 2.4 s in Fig. 4(a). After the relaxation process at t ∼ 2.4 s, we slowly turn the bias velocity to a maximum value of v bias = 10.6 µm/s, which is a factor of 10 smaller than the value used in Fig. 2. This results in a linear growth for ∆n(t), visible in Fig. 4(a). The magnitude of ∆n(t) depends on the flux velocity v Φ /v 0 . For half integer multiples, i.e. v Φ /v 0 = (n + 1/2), we expect a high resistance, i.e. strong imbalance, because the critical velocity is minimal at this point. At whole integer values the resistance should be minimal as no screening currents are present in the system. This is consistent with the results shown in Fig. 4(a). We repeat the above protocol for multiple values of v Φ /v 0 and extract ∆n(t) at t = 3.8 s using a linear fit. We show these results in Fig. 4(b). As described above, we find a minimum of ∆n at full integer values of v Φ /v 0 and a maximum at half integer values of v Φ /v 0 . This dependence on v Φ /v 0 results in oscillations of the density imbalance, which is analogous to voltage-flux oscillations of a solid-state SQUID. The observed oscillations are direct evidence for quantum interference of currents. The voltage-flux relation of the dc-SQUID in the overdamped limit is [33] which assumes negligible screening and yields the time averaged junction voltage as a function of the applied flux Φ. I c is the critical current of a single junction, I is the total current, R is the resistance and Φ 0 is the magnetic flux quantum. We emphasize that the condensate realization of a SQUID is not in the overdamped regime, and that we merely use this expression as a fitting function, and more generally as a comparison, see also Appendix B for comparison to a dc-SQUID model with nonzero capacitance and inductance. We fit this expression to our results and show the fit as a continuous line in Fig. 4. We find a qualitatively good agreement with the analytic expression of Eq. 7. In particular, the periodic nature of the dependence is captured, which is the central, defining property of a SQUID. However, the fit seems to underestimate the amplitude of the oscillations, which is captured better by the nonzero capacitance and inductance dc-SQUID model as described in Appendix B. From the fit, we extract the parameters I c R = 3.1 and I/(2I c ) = 1.005. The current is more than twice the critical current, which is consistent with the SQUID being in the resistive regime.

VI. CONCLUSIONS
In conclusion, we have put forth a proposal of how to create an atomtronic SQUID in a toroidal quasi-1D condensate. For that purpose, we simulate the condensate dynamics via a numerical implementation of a multiband truncated Wigner approximation. We note that this method is ideally suited for the dynamics of fluctuating quasi-1D condensates, such as in thin-ring geometries [49]. We stir the condensate with two mobile barriers to induce rotation of the cloud, which is equivalent to a magnetic flux applied to a conventional SQUID. This allows us to study the phase-slip dynamics and its dependence on the stirring time and the barrier height. For long stirring times and strong barrier heights the average phase winding results in a step-like behavior as a function of the flux velocity. We then operate the SQUID in this regime at a constant bias velocity to create a density imbalance, which shows characteristic oscillations with a periodicity of the flux velocity quantum. This highlights the voltage regime for the operation of atomtronic SQUIDs and demonstrates the quantum interference of currents.

VII. ACKNOWLEDGEMENTS
We thank Kevin Wright, Luigi Amico and Giacomo Roati for insightful discussions. This work is supported by the Deutsche Forschungsgemeinschaft (DFG) in the framework of SFB 925 -project ID 170620586 and the excellence cluster 'Advanced Imaging of Matter' -EXC 2056 -project ID 390715994, and the Cluster of Excellence 'QuantumFrontiers' -EXC 2123 -project ID 390837967.

Appendix A: Numerical Coefficients
We define dimensionless coefficients c via where U = g/l 3 and l osc is in units of the discretization length l. These are given by an integral where φ n (x) = (π) −1/4 (2 n n!) −1/2 exp(− are dimensionless harmonic oscillator functions. The normalization of of the c coefficients is chosen such that c 0000 = 1. Note that due to the symmetry properties of the Hermite functions, c abcd is zero when a + b + c + d is odd. Also c abcd is invariant under permutation of the indices. We provide the coefficients used in our calculation in table.
The coefficients for the barrier potential are determined by This results in We numerically evaluate the integral in Eq. A6 to obtain the coefficients.

Appendix B: Voltage-flux relation of dc-SQUIDs
For practical SQUIDs the inductance of the circuit and the capacitance of Josephson junctions are taken into account, corresponding to the circulating screening current and the junction displacement current, respectively. This general case of the dc-SQUID circuit is described by the equations [50] with the integer constraint with Φ 0 = π e being the flux quantum and Φ being the applied flux. φ 1 and φ 2 are the phase differences across the junctions, I c is the critical current and C is the capacitance of single junction. I cir is the circulating screening current and L is the inductance. R N is the resistance. ∆n (µm A1. Simulated density imbalance ∆n as a function of vΦ/v0, which is the same as Fig. 4(b). The continuous line is the fit with the over-damped model of dc-SQUIDs in Eq. 7, while the dashed line corresponds to the nonzero capacitance and inductance SQUID model of Eq. B13.
We introducē and and φ cir =πβ L I cir .
For parameter β L we have n 1d is the average density projected onto the spatial degree of freedom along the condensate. In terms of rescaled units we have: We fulfill this constraint viā where R(x) rounds the real number x to the nearest integer number. The voltage is We numerically solve Eq. B13 to determine the voltageflux relation V (φ). The screening parameter β L and the Stewart-McCumber parameter β C both have to be smaller than unity to avoid hysteretic V (φ) curves. As an example, we therefore choose β C = 0.1 and β L = 0.1. I = 2.01 is chosen according to the fit result in Sec. V. We calculate the voltage-flux relation and scale its magnitude by the factor α = 2.04, which accounts for the resistance, i.e., α = R N I c . This result is shown in Fig.  A1, which describes the simulation result better than the over-damped model. This is also reflected by the value of the resistance being smaller than that of the over-damped fit.