Josephson junction dynamics in a two-dimensional ultracold Bose gas

We investigate the Berezinskii-Kosterlitz-Thouless (BKT) scaling of the critical current of Josephson junction dynamics across a barrier potential in a two-dimensional (2D) Bose gas, motivated by recent experiments by Luick \textit{et al.} arXiv:1908.09776. Using classical-field dynamics, we determine the dynamical regimes of this system, as a function of temperature and barrier height. As a central observable we determine the current-phase relation, as a defining property of these regimes. In addition to the ideal junction regime, we find a multimode regime, a second-harmonic regime, and an overdamped regime. For the ideal junction regime, we derive an analytical estimate for the critical current, which predicts the BKT scaling. We demonstrate this scaling behavior numerically for varying system sizes. The estimates of the critical current show excellent agreement with the numerical simulations and the experiments. Furthermore, we show the damping of the supercurrent due to phonon excitations in the bulk, and the nucleation of vortex-antivortex pairs in the junction.


I. INTRODUCTION
For a Josephson junction created by a superconductorinsulator-superconductor interface, the Josephson relation I s = I c sin φ relates the supercurrent I s to the phase difference φ of the order parameter across the junction [1].I c is the critical current of the junction, which is the maximal supercurrent across the junction.This connection is the defining functionality of quantum mechanical devices, such as superconducting quantum interference devices (SQUIDs).
Josephson junctions are utilized in phenomenological models of high-temperature superconductors, which describe these materials as stacks of two-dimensional (2D) systems coupled by JJs.Josephson junction arrays also serve as a model to describe transport phenomena in optically driven high-temperature superconductors [37].2D systems such as thin film superconductors [38] or Josephson junction arrays [39] undergo a Berezinskii-Kosterlitz-Thouless (BKT) transition within the XY universality class [40][41][42].The superfluid phase has quasi-long-range order and is characterized by a scale-invariant exponent τ of the single-particle correlation function g 1 that decays algebraically at large distances, g 1 (r) ∼ |r| −τ /4 .At the transition the exponent assumes the critical value τ c = 1 which is accompanied by a universal jump of the superfluid density.
Recently, Ref. [31] reported on a study on Josephson junction dynamics of an ultracold 2D gas of 6 Li atoms, which is realized by separating two uniform 2D clouds with a tunneling barrier.Using a strong barrier higher than the mean-field energy, the experiments measure the current-phase relation of an ideal junction, and the critical currents in the crossover from tightly bound molecules to weakly bound Cooper pairs.In this paper, we establish a connection between the quasi-order scaling of 2D superfluids and the critical current of the Josephson junction.Specifically, we demonstrate the BKT scaling of the critical current of a Josephson junction coupling two 2D Bose gases using classical field simulations.The Josephson junction is created by a tunneling barrier between two 2D clouds of 6 Li 2 molecules, motivated by the experiments of Ref. [31].We find an interplay of the bulk and junction dynamics that is influenced by the barrier height and the temperature.Depending on these parameters, we find multimode (MM) and second-harmonic (SH) contributions to the current-phase relation.For large barrier heights we find that the junction dynamics displays ideal Josephson junction (IJJ) behavior, i.e. it obeys the nonlinear current-phase relation I(φ) = I c sin(φ).We map out the MM, the SH, the IJJ, and an overdamped regime as a function of barrier height and temperature.We determine the critical current numerically based on the current and phase dynamics at the barrier.The numerically obtained critical current, and an analytical estimate that we derive, show excellent agreement with the experimental values of Ref. [31].We confirm the BKT scaling of the critical current by performing simulations for varying system sizes.The exponent of the critical current across the transition demonstrates agreement with the exponent of the corresponding equilibrium system, if the system is in the IJJ regime.Finally, we address the damping of the current and identify the damping mechanism which is due to phonon excitations in the bulk, and the nucleation of vortex-antivortex pairs in the junction.This paper is organized as follows.In Sec.II we describe our simulation method.In Sec.III we show the condensate dynamics and its dependence on the barrier height and the temperature.In Sec.IV we determine the critical current for an ideal junction and compare it to the simulation and the experiment.In Sec.V we show the power-law scaling of the critical current for varying system sizes.In Sec.VI we discuss the dissipation mechanism of the current, and in Sec.VII we conclude.

II. SIMULATION METHOD
Motivated by the experiments [31] we study 2D clouds of 6 Li 2 molecules confined in a box of dimensions L x ×L y .We simulate the dynamics using the c-field method of Ref. [43].The system is described by the Hamiltonian Ĥ0 = dr (1) ψ ( ψ † ) is the bosonic annihilation (creation) operator.The interaction g is given by g = g 2 /m, where g is the dimensionless interaction, and m the molecular mass.g is determined by g = g0 / 1 − g0 2π ln(2.09kF z ) , with g0 = √ 8πa s / z [44].a s is the molecular s-wave scattering length, z = /(mω z ) the harmonic oscillator length in the transverse direction, and k F the Fermi wavevector.We discretize space on a lattice of size N x × N y and a discretization length l = 0.5 µm.Within the c-field representation we replace the operators ψ in Eq. 1 and the equations of motion by complex numbers ψ.We sample the initial states in a grand canonical ensemble having chemical potential µ and temperature T via a classical Metropolis algorithm.We choose the system parameters, such as the density n, g, and T to be close to the experiments.In particular, we choose n ≈ 1.24 µm −2 , T /T 0 = 0.3, and L x × L y = 20 × 40 µm 2 , which are the same as in Ref. [31].The critical temperature T 0 is estimated by T 0 = 2πn 2 /(mk B D c ), where D c = ln(380/g) is the critical phase-space density [45].We vary g in the range 1 − 4 to cover the BEC regime of the experiments.For additional simulations, we use n ≈ 2.25 µm −2 and g = 1.8, while we vary T /T 0 and the box size.
To create the Josephson junction we add a barrier term H b (t) = dr V (r, t)n(r), where n(r) is the density at the location r = (x, y).The barrier potential V (r, t) is given by where V 0 (t) is the time-dependent strength and w the width.The potential is centered at x 0 = L x /2.We choose w in the range (0.85 − 2) µm and V 0 in the range V 0 /µ ≡ Ṽ0 = 0 − 6, where µ = gn is the mean-field energy.We ramp up V 0 linearly over 150 ms and then wait for 50 ms.This splits the system in x-direction into two uniform 2D clouds, which we refer to as the left and right reservoir.To create a phase difference across the junction, we imprint a phase φ 0 on the left reservoir, resulting in the phase difference φ 0 = φ L − φ R , where φ L (φ R ) is the mean phase of the left (right) reservoir.The sudden imprint of the phase and the barrier lead to the dynamics displayed in Fig. 1.We calculate the x component of the current density defined as: We calculate j 1 with r = (x 0 − lê x , y), and j 2 with r = (x 0 , y). j = j 1 + j 2 /2 gives the averaged current density at the barrier center.This current fulfills the continuity equation of the density imbalance between the two reservoirs.The time evolution of the current is determined by I(t) = j(t)L y , see Fig. 6(a).We fit I(t) with the function f (t) = I 0 e −Γt sin(ωt + θ) to determine the magnitude I 0 ≡ |I 0 |, the damping rate Γ, the frequency ω, and the phase shift θ.

III. BULK VERSUS JUNCTION DYNAMICS
To characterize the junction dynamics we analyze the time evolution of the density and the phase of the system.As an illustration, we choose n = 2.25 µm −2 , T /T 0 = 0.3, and g = 1.8.We imprint a phase of π/4 on the left reservoir at t = 0, for Ṽ0 = 0, 0.8, and 2. We use w = 0.85 µm, which in terms of the healing length ξ is w/ξ = 2.4, where ξ = / √ 2µm = 0.35 µm.In Fig. 1(a) we show the time evolution of the density δn(x, t) = n(x, t)−n(x), which is averaged over the y direction and the thermal ensemble.n(x) is the averaged equilibrium density profile.For no barrier, Ṽ0 = 0, the phase imprint creates two density pulses that are visible as density increase and decrease in the left and right reservoir, respectively.The density pulses propagate in opposite direction with the sound velocity and are reflected by the box edges.For Ṽ0 = 0.8, the barrier confines the density wave partially in each of the reservoirs.For Ṽ0 = 2, the density waves are well confined within the reservoirs as flow between the reservoirs is obstructed by the barrier.Instead, the density waves tunnel across the barrier, resulting in coherent Josephson oscillations between the reservoirs.In Fig. 1(b) we show the time evolution of the phase δφ(x, t) = φ(x, t) − φ m (t) of a single trajectory, for the same Ṽ0 as in Fig. 1(a).φ m is the mean global phase.The sudden imprint of phase adds the mean phase difference φ 0 = φ L − φ R between the left and right reservoir at t = 0.During the time evolution a phase gradient develops within the reservoirs, which results in φ L/R being different from the phases close to the junction.For Ṽ0 = 0, the phase evolves linearly with the distance.As the barrier height is increased, the phase gradient within the reservoirs decreases.For Ṽ0 = 2, the phase gradient within the reservoir almost vanishes, resulting in φ L/R being the same as the phases in direct vicinity of the junction.This corresponds to IJJ dynamics.We now examine the current-phase relation (CPR) of the junction.We determine the current I 0 by fitting the time evolution of the current I(t) to a damped sinusoidal function as described in Sec.II.To obtain the CPR we calculate I 0 as a function of φ 0 .In Fig. 1(c) we show I 0 (φ 0 ) for various values of Ṽ0 at T /T 0 = 0.1, and in Fig. 1(d) at T /T 0 = 0.3.We analyze these CPR curves by fitting them with a multi-harmonic fitting function I(φ 0 ) = nmax n=1 I n sin(nφ 0 ), where we choose n max = 5.We find that the result of these fits can be grouped into three regimes, as a function of the barrier height Ṽ0 and the temperature.If the coefficient I 1 is dominant, the CPR reduces to the form of an IJJ, I(φ 0 ) = I 1 sin(φ 0 ).If both I 1 and I 2 are non-negligible, we refer to the regime as second harmonic (SH) regime.In Fig. 1(e) we indicate the crossover from IJJ to SH by depicting the ratio I 2 /I 1 .For lower temperatures and barriers, higher harmonical contributions become important.We indicate the multimode (MM) regime if n>2 (I n /I 1 ) 2 > 0.02, see also Appendix A. We note that CPR deviations were pointed out for Ṽ0 1 by Refs.[20,46,47] .Furthermore, we indicate the overdamped (OD) regime based on the analysis shown in Sec.VI.
The results depicted in Fig. 1(e) demonstrate that the IJJ regime is strongly sensitive to the temperature and the barrier height.This derives from the properties of the dynamical evolution shown in Figs.1(a) and (b).The initial phase imprint creates phonon pulses in the two reservoirs.For low temperatures and small barrier heights these pulses are weakly damped, which leads to the multimode regime.For increasing temperature and barrier height, fewer and fewer of the phonon modes of the system contribute.The increasing barrier height leads to a long tunneling time, which exceeds the damping time of more and more phonon modes, until the phase dynamics reduces to the dynamics of two global phases for each reservoir, as visible in Figs.1(a) and (b).However, if the barrier is increased further, eventually the dynamics become overdamped.Here, the two reservoirs dephase on the timescale of the tunneling rate.We note that the temperature dependence of the current-phase relation was measured in Nb/InAs/Nb junctions by Ref. [48] and in a weak link of 4 He superfluids by Ref. [49].

IV. CRITICAL CURRENT
To determine the Josephson critical current, we calculate the junction phase φ j ≡ φ j = φ 2,m − φ 1,m , where φ 1/2,m = φ x0∓2lêx are the mean phases calculated by taking an average of the fields over the y direction.We use the same density n as above, and T /T 0 = 0.2.We calculate the time evolution of the current I(t) and the phase φ j (t).In Fig. 2(a) we show the time evolution of I(φ j ) for Ṽ0 = 1, 2, and 3.In this small φ j regime a linear behavior is observed.The width of the distribution increases with increasing Ṽ0 due to increased phase fluctuations across show Ic( Ṽ0) for T /T0 = 0.1, 0.2, and 0.3, respectively, for w/ξ = 2.9 (squares) and 4.3 (circles).The continuous lines are the estimates of Eq. 5.
the barrier.We determine the critical current I c by the slope I/ sin φ j .Within the IJJ regime, this coincides with I 1 .I c decreases with increasing Ṽ0 , with I c = 83.8,15, and 2.9 ms −1 for Ṽ0 = 1, 2, and 3, respectively.In Figs.We derive an analytical estimate of the critical current by solving the mean-field equation and by considering single-particle tunneling across a rectangular barrier of height V > µ and width d, see Appendix B. We obtain the current density with and n 0 is the condensate density.κ is the damping parameter of the exponential wavefunction inside the barrier, which is determined variationally, and includes the mean-field repulsion under the barrier.We interpret this result for j c as the product of the density n 0 µ/ V + V 2 − µ 2 /2 at the barrier boundary and the velocity κ/m at the barrier center.Alternatively, it is instructive to rewrite Eq. 5 in terms of the bulk current density cn 0 , and the tunneling amplitude t 0 ( Ṽ , d) across the barrier as with c = µ/m is the sound velocity, ξ the healing length, and Ṽ = V /µ the scaled barrier height.We note that j c is described in terms of the bulk current and the tunneling amplitude for a 3D condensate in Refs.[28,29].We determine the estimate I c = j c Ly by using V = V 0 and d ≈ 1.2w, and by determining c(T ) and n 0 numerically.V 0 and w are of the Gaussian barrier used in the simulation.This value for d is set by fitting the simulated critical currents in the IJJ regime, which is close to our assumption d ≈ w used in Ref. [31].In Figs.2(b)-(d) we show the estimates as a function of Ṽ0 for T /T 0 = 0.1, 0.2, and 0.3.The estimates agree with simulated critical currents at all Ṽ0 > 1, for all T /T 0 .The agreement is particularly good for the IJJ regime.This suggests that the barrier reduction of I c is due to the tunneling amplitude t 0 ( Ṽ0 , w) that decreases with increasing Ṽ0 and w.We note that the barrier width reduction follows an exponential behavior given by Eq. 5.
In Fig. 3(a) we show I c as a function of Ṽ0 for w/ξ = 2.9 and various T /T 0 .As T /T 0 increases, I c ( Ṽ0 ) decreases with a rather sudden jump to very low values for T /T 0 > 0.6.I c and the bulk current are connected according to Eq. 7 via I B = I c ( Ṽ0 )/t 0 ( Ṽ0 , w).We thus divide I c ( Ṽ0 ) by t 0 ( Ṽ0 , w), see inset of Fig. 3(b).The results are almost independent of Ṽ0 as expected.By taking an average over the range Ṽ0 = 1 − 2, we obtain the mean value of I B which is shown in Fig. 3(b).The mean I B decreases with increasing T /T 0 and becomes small due to increased thermal fluctuations for T /T 0 ≥ 0.6.We compare this result to the actual bulk value I B = c(T )n 0 L y by determining n 0 and c(T ) numerically for various T /T 0 .This bulk result agrees at intermediate temperatures, where the system is near the IJJ regime.At low temperatures, the system is in the MM regime, where the above estimate is not valid.The bulk current is linked to the BKT scaling exponent that we determine in Sec.V. Finally, we compare the simulations to the measurements that are performed at several interaction strengths in the BEC regime [31].We use the range g = 1 − 4, and determine the critical current I c as described above.In ate from the experimental result, because the system enters the strongly interacting regime, where the c-field approach starts to deviate systematically.The shaded areas reflect the error bars of theory when assuming a 15% uncertainty of Ṽ0 as in experiment.We calculate the estimate of Eq. 5 by determining the condensate density numerically for all interactions.The condensate fraction n 0 /n has weak interaction dependence and is about 62%.We show the analytical estimates in Fig. 4. The estimates, including the 15% uncertainty of Ṽ0 , are consistent with both simulation results and the experimental results.

V. BEREZINSKII-KOSTERLITZ-THOULESS SCALING
We have shown above that the critical current depends on the condensate density.We use that dependence to extract the scaling exponent of the quasi-condensate.The condensate density scales algebraically with the system size as where τ (T ) is the temperature-dependent exponent, L the system length, and r 0 the short-range cutoff of the order of ξ.Employing this scaling, we first determine τ (T ) for a square-shaped system in equilibrium.We choose n ≈ 2.25 µm −2 , g = 1.8, and L in the range (2 − 128) µm.We calculate n 0 as a function of L for various T /T 0 for a system with periodic boundary conditions.In Fig. 5(a) we show the condensate fraction n 0 /n as a function of The continuous lines are the algebraic fits.The exponential fit (dashed line) is a good fit for T /T0 = 0.9 in panel (a).(c) Extracted exponents of the equilibrium system (squares) and jc (circles).The black cross corresponds to the measurements of jc in Ref. [31].
L for various T /T 0 .At low and intermediate T /T 0 , n 0 shows a power-law behavior as it decreases linearly with L on a log-log scale.At high T /T 0 , n 0 deviates from this power-law scaling and instead shows an exponential behavior which is characteristic of the thermal phase [50].
To confirm the power-law scaling, we fit n 0 /n to Eq. 9, with τ and r 0 as fitting parameters.We show the fits in Fig. 5(a).The algebraic fits describe the behavior very well for T /T 0 ≤ 0.6, whereas they fail to capture the dependence for T /T 0 > 0.6.For T /T 0 = 0.9 the dependence is captured by the exponential fit, while for T /T 0 = 0.7 the algebraic fit is better than the exponential fit.Next, we determine the scaling exponent of the critical current density.We choose n and g as above, and L in the range (40 − 128) µm.For the barrier, we use w/ξ = 2.9 and Ṽ0 in the range Ṽ0 = 1.2 − 2. We calculate the critical current density j c as a function of L for various T /T 0 , where j c is determined by the critical current described in Sec.IV.We obtain the condensate density n 0 by normalizing j c with the sound velocity c(T ) and the tunneling amplitude t 0 ( Ṽ0 , w), see Eq. 7. We then aver- age n 0 over the barrier heights employed.We expand on this normalization in Appendix C. In Fig. 5(b) we show the condensate fraction n 0 /n that is determined from the critical current density as a function of L for various T /T 0 .n 0 /n shows a power-law behavior for T /T 0 ≤ 0.6, which is confirmed by the algebraic fits shown in Fig. 5(b).In Fig. 5(c) we show the temperature dependence of τ for the equilibrium system and τ determined from the critical current density.The equilibrium value of τ increases linearly at low and intermediate T /T 0 , while it deviates from linear behavior at high T /T 0 in the crossover regime.We estimate the transition temperature with the BKT critical value τ c = 1, which gives T /T 0 ≈ 0.6.We note that this value of the critical temperature is renormalized to a lower value in thermodynamic limit [51].
We also note that this estimate is below T 0 of a weakly interacting 2D Bose gas [45].Above the transition, τ increases rapidly with the temperature.This temperature dependence of τ across the transition is described by the RG equations [52].Studies of BKT scaling in ultracold gases were reported in Refs.[53][54][55].In addition to the equilibrium value of τ we show the value of τ based on the critical current scaling.The results show excellent agreement with the exponents of the equilibrium system for the temperatures 0.2 < T /T 0 < 0.55, and follow the qualitative behavior outside of this temperature range.
The deviations below T /T 0 = 0.2 are due to multimode dynamics that influence the results of the critical current density, while the deviations for T /T 0 ≥ 0.55 are due to thermal excitations at the barrier, and the onset of overdamped dynamics.Furthermore, we determine the exponent τ from the measurements of the critical currents shown in Fig. 4. Making use of Eq. 7, these measurements yield the condensate fraction n 0 /n = 0.72 (8).For the system size L in the experiment and r 0 = ξ, the scaling of Eq. 9 results in an exponent of τ = 0.32 (12).We show this value of τ for T /T 0 = 0.3 in Fig. 5(c), which agrees with the exponents of the simulated critical current density and the equilibrium system.

VI. CURRENT DAMPING AND DISSIPATION MECHANISM
Here we analyze the damping of the supercurrent and identify the associated dissipation mechanism.As an illustration, we choose w/ξ = 2.9 and Ṽ0 = 2.0, and calculate the time evolution of the current I(t) as described above.In Fig. 6(a) we show I(t) for T /T 0 = 0.1, 0.3, and 0.5.The current oscillations are underdamped at T /T 0 = 0.1 and 0.3.The damping increases with increasing T /T 0 .For T /T 0 = 0.5, the current undergoes an overdamped motion.To quantify this observation we determine the oscillation frequency ω and the damping rate Γ as described in Sec.II.In Fig. 6(b) we show ω/ω 0 determined as a function of Ṽ0 for T /T 0 = 0.1, 0.3, and 0.5.ω 0 is the oscillation frequency for Ṽ0 = 0, which we refer to as the sound frequency.As Ṽ0 increases, ω/ω 0 decreases.This decrease is more pronounced for higher temperatures.We note that the dependence of ω/ω 0 on the barrier height differs qualitatively between the low ( Ṽ0 < 1) and high ( Ṽ0 > 1) barrier regimes.In Fig. 6(c) we show the results of Γ/ω as a function of Ṽ0 .For T /T 0 = 0.1, Γ/ω is small at all Ṽ0 , confirming underdamped motion.For T /T 0 = 0.3, Γ/ω increases at high Ṽ0 and is generally below 0.5 that we use as the definition of the temperature-induced overdamped limit.For T /T 0 = 0.5, Γ/ω increases rapidly with Ṽ0 and reaches the overdamped limit at Ṽ0 ≥ 1.5.
To identify the origin of the damping we examine the phase dynamics of a single trajectory of the ensemble.We calculate the phase φ(x, y) = φ(x, y) − φ m for the same parameters as in Fig. 6(a), where φ m is the mean global phase.In Fig. 7(a) we show φ(x, y) at t = 3.6 ms, after the phase imprint, for T /T 0 = 0.1, 0.3, and 0.5.The phase imprint develops a phase difference between the two reservoirs and a corresponding phase gradient across the barrier.At low temperature the reservoir phase is weakly fluctuating as demonstrated by φ(x, y) at T /T 0 = 0.1.The fluctuations of the phase increase with increasing T /T 0 .The reservoir phase is moderately and strongly fluctuating for T /T 0 = 0.3 and 0.5, respectively.This results in the creation of vortices, which is confirmed by the calculation of the phase winding around the plaquettes of the numerically introduced lattice.We calculate the phase winding around the lattice plaquette of size l × l using δφ(x, y) = δ x φ(x, y) + δ y φ(x + l, y) + δ x φ(x + l, y + l) + δ y φ(x, y + l), where the phase differences between sites are taken to be δ x/y φ(x, y) ∈ (−π, π].We show the calculated phase windings in Fig. 7(a).We identify a vortex and an antivortex by a phase winding of 2π and −2π, respectively.
For T /T 0 = 0.1, we observe only one vortex-antivortex pair inside the barrier and no vortices in the bulk.This scenario changes due to increased thermal fluctuations at high temperatures.For T /T 0 = 0.3, there is nucleation of multiple vortex pairs inside the barrier, and a few vortex pairs near the box edges.For T /T 0 = 0.5, we observe proliferating vortices in the regions of low densities around the barrier, and in the bulk.
To understand the role of vortex fluctuations at the barrier, we calculate the total number of vortices N v and average it over the thermal ensemble.In Fig. 7(b) we plot N v as a function of Ṽ0 for the same values of T /T 0 as in Fig. 7(a).At T /T 0 = 0.1, N v remains close to zero for barrier heights below a threshold value and beyond this N v increases with increasing Ṽ0 .At high temperatures the system features thermal vortices even in the absence of the barrier and the onset of barrier-induced vortices occurs at a Ṽ0 lower than that at low temperature.To determine this threshold Ṽc we calculate the differential vortex number N v ( Ṽ0 ) = N v − N v,0 , where N v,0 = N v ( Ṽ0 = 0).We show N v ( Ṽ0 ) in Fig. 7(c).We define Ṽc for which N v ( Ṽ0 ) approaches 1.This gives Ṽc = 1.47, 0.58, and 0.29 for T /T 0 = 0.1, 0.3, and 0.5, respectively.This onset of vortices is associated with the damping of the oscillation shown above.

VII. CONCLUSIONS
In this paper, we have established a direct connection between the Josephson critical current and BKT scaling in an ultracold 2D Bose gas using classical field simulations.For this, we have examined the dynamics across a Josephson junction created by a tunnel barrier between two uniform 2D clouds of 6 Li 2 molecules, which is motivated by the experiments of Ref. [31].Based on the current-phase relation, we have mapped out the multimode, the second-harmonic (SH), the ideal junction (IJJ), and the overdamped regime as a function of the barrier height and the temperature.For the IJJ regime, we have derived an analytical estimate of the critical current, which is in good agreement with the simulations and the experiments [31].We have demonstrated the BKT scaling of the critical current numerically by varying the system size.The scaling exponents of the critical current are in agreement with the exponents of the corresponding equilibrium system.Finally, we have addressed the damping of the current, which is due to phononic excitations in the bulk, and the nucleation of vortex pairs in the junction.
In conclusion, we have discussed the dynamics of atomic clouds in 2D, coupled via a Josephson junction, which results in the hybridization of the bulk and tunneling dynamics.As such, it combines and relates two foundational effects of quantum physics, in particular condensation and Josephson oscillations.Our results demonstrate a method to measure a static property of manybody order, in particular the condensate density, via a dynamical oscillatory process, in particular Josephson oscillations.Both the principle of this method, as well as the presented discussion of dynamical regimes of this system, can be applied to a wide range of quantum gas systems, to gain insight into their dynamical and static properties.To determine the algebraic scaling exponent of the quasi-condensate, we calculate j c as a function of Ṽ0 for varying system sizes with simulations of the squareshaped box.We use w/ξ = 2.9, and Ṽ0 in the range 1.2 − 2. The parameters n, g, and L are the same as in the main text in Sec.V. We determine j c as described in the main text.To obtain n 0 we divide j c by the sound velocity c(T ) and the tunneling amplitude t 0 ( Ṽ0 , w), see Eq. C3.In Fig. 9 we show the condensate fraction n 0 /n as a function of L for T /T 0 = 0.2, 0.3, and 0.4.As expected, the results of n 0 /n are almost independent of Ṽ0 and demonstrate a decreasing behavior with increasing L and T /T 0 .To determine the scaling exponent τ , we fit n 0 /n to the function n 0 /n = (L/r 0 ) −τ /4 , with τ and r 0 as fitting parameters.We show the determined values of τ for T /T 0 = 0.2, 0.3, and 0.4 in Fig. 10.The results are in agreement within the error bars for the barrier heights employed.For comparison, we average n 0 /n over the barrier heights employed and then determine τ from this averaged condensate fraction as we do in the main text in Sec.V.This result is also shown in Fig. 10, where it agrees with the determined exponents for Ṽ0 in the range 1.2 − 2.

FIG. 1 .
FIG. 1. Dynamical regimes.(a) Time evolution of the density δn(x, t) = n(x, t) − n(x), which is averaged over the y direction and the ensemble, for Ṽ0 = 0.0, 0.8, and 2.0.n(x) is the equilibrium density.Panel (b) shows the corresponding phase evolution δφ(x, t) = φ(x, t) − φm(t) of a single trajectory of the ensemble.φm is the mean phase.We imprint a phase of φ0 = π/4 on the left reservoir at t = 0 for n = 2.25 µm −2 and T /T0 = 0.3.The barrier width w is denoted by the two vertical dotted lines.Panels (c, d) show the current I0 as a function of φ0 for various Ṽ0 at T /T0 = 0.1 and 0.3, respectively.The continuous lines are fits with the fitting function I(φ0) = I1 sin(φ0) + I2 sin(2φ0).(e) The dynamical regimes are multimode (MM), second-harmonic (SH), ideal Josephson junction (IJJ), and overdamped (OD), see text.
2(b)-(d) we show I c ( Ṽ0 ) for T /T 0 = 0.1, 0.2, and 0.3, and the barrier widths w/ξ = 2.9 and 4.3.As expected, I c is smaller for wider barriers.I c also decreases with the temperature as we show below.

FIG. 3 .
FIG. 3. Temperature dependence.(a) Ic( Ṽ0) for w/ξ = 2.9 and various T /T0.(b) Mean bulk current IB (circles) determined from the normalized critical currents shown in the inset, see text.The error bar denotes the standard deviation.The bulk values of the current determined using the condensate density and phonon velocity are shown by the diamonds.

Fig. 4 FIG. 4 .
FIG. 4.Comparison to the experiments.The measurements of the critical current (crosses) are compared to the simulations (squares) and the analytical estimate (continuous line) for various ln(kFa2D) in the BEC regime.kF is the Fermi wavevector and a2D is the 2D scattering length.We use n ≈ 1.24 µm −2 , T /T0 ≈ 0.3, w = 0.81 µm, and Ṽ0 = 1.4,which are the same as the experiments.The shaded areas include a 15% uncertainty of Ṽ0 as in experiment.The measurement data is from Ref.[31].

FIG. 5 .
FIG.5.Determining the scaling exponent.(a) Condensate fraction n0/n of the equilibrium system as a function of system length L on a log-log scale for various T /T0.(b) Averaged n0/n determined from the critical current density jc via Eq. 7. The continuous lines are the algebraic fits.The exponential fit (dashed line) is a good fit for T /T0 = 0.9 in panel (a).(c) Extracted exponents of the equilibrium system (squares) and jc (circles).The black cross corresponds to the measurements of jc in Ref.[31].