Shapiro steps in strongly-interacting Fermi gases
Giulia Del Pace Diego Hernández-Rajkov Vijay Pal Singh Nicola Grani Marcia Frómeta Fernández Giulio Nesti Jorge Amin Seman Massimo Inguscio Luigi Amico Giacomo Roati









[E-mail: ] delpace@lens.unifi.it
Department of Physics, University of Florence, 50019 Sesto Fiorentino, Italy
European Laboratory for Nonlinear Spectroscopy (LENS), University of Florence, 50019 Sesto Fiorentino, Italy
Istituto Nazionale di Ottica del Consiglio Nazionale delle Ricerche (CNR-INO) c/o LENS, 50019 Sesto Fiorentino, Italy



European Laboratory for Nonlinear Spectroscopy (LENS), University of Florence, 50019 Sesto Fiorentino, Italy
Istituto Nazionale di Ottica del Consiglio Nazionale delle Ricerche (CNR-INO) c/o LENS, 50019 Sesto Fiorentino, Italy


Quantum Research Centre, Technology Innovation Institute, Abu Dhabi, UAE


Department of Physics, University of Florence, 50019 Sesto Fiorentino, Italy
European Laboratory for Nonlinear Spectroscopy (LENS), University of Florence, 50019 Sesto Fiorentino, Italy
Istituto Nazionale di Ottica del Consiglio Nazionale delle Ricerche (CNR-INO) c/o LENS, 50019 Sesto Fiorentino, Italy




European Laboratory for Nonlinear Spectroscopy (LENS), University of Florence, 50019 Sesto Fiorentino, Italy
Istituto Nazionale di Ottica del Consiglio Nazionale delle Ricerche (CNR-INO) c/o LENS, 50019 Sesto Fiorentino, Italy


European Laboratory for Nonlinear Spectroscopy (LENS), University of Florence, 50019 Sesto Fiorentino, Italy
Istituto Nazionale di Ottica del Consiglio Nazionale delle Ricerche (CNR-INO) c/o LENS, 50019 Sesto Fiorentino, Italy


Instituto de Física, Universidad Nacional Autonoma de Mexico, C.P. 04510 Ciudad de Mexico, Mexico


Department of Engineering, Campus Bio-Medico University of Rome, Rome, Italy
European Laboratory for Nonlinear Spectroscopy (LENS), University of Florence, 50019 Sesto Fiorentino, Italy
Istituto Nazionale di Ottica del Consiglio Nazionale delle Ricerche (CNR-INO) c/o LENS, 50019 Sesto Fiorentino, Italy


Quantum Research Centre, Technology Innovation Institute, Abu Dhabi, UAE
INFN-Sezione di Catania, Via S. Sofia 64, 95127 Catania, Italy
Dipartimento di Fisica e Astronomia, Università di Catania, Via S. Sofia 64, 95123 Catania, Italy


European Laboratory for Nonlinear Spectroscopy (LENS), University of Florence, 50019 Sesto Fiorentino, Italy
Istituto Nazionale di Ottica del Consiglio Nazionale delle Ricerche (CNR-INO) c/o LENS, 50019 Sesto Fiorentino, Italy




§ ABSTRACT 
We report the observation of Shapiro steps in a periodically driven Josephson junction between strongly-interacting Fermi superfluids of ultracold atoms. We observe quantized plateaus in the current-potential characteristics, the height and width of which mirror the external drive frequency and the junction nonlinear response.

Direct measurements of the current-phase relationship showcase how Shapiro steps arise from the synchronization between the relative phase of the two reservoirs and the external drive.
Such mechanism is further supported by the detection of periodic phase-slippage processes, 

in the form of vortex-antivortex pairs.
Our results are corroborated by a circuital model and numerical simulations, overall providing a  clear  understanding of Shapiro dynamics in atomic Fermi superfluids. Our work demonstrates phase-coherent and synchronization effects in driven strongly-interacting superfluids, opening
prospects for studying emergent non-equilibrium dynamics in quantum many-body systems  under external drives.


Shapiro steps in strongly-interacting Fermi gases
    G. Roati
    September 9, 2024
=================================================

Driven many-body systems exhibit rich and complex behaviors that far extend beyond their equilibrium properties 
<cit.>.
A paradigmatic example is provided by Josephson junctions (JJs) under driving currents. A Josephson junction comprises two superconducting islands separated by a thin insulating barrier. Josephson showed that a dissipationless current I can flow across the junction
, maintained solely by the phase difference ϕ between the two superconducting reservoirs <cit.>. 
Above a critical current I_c, the junction transitions into the so-called voltage-state: a finite voltage V develops across the junction, which is directly related to the 
phase difference dynamics as provided by 
the Josephson-Anderson relation V=ħ/2eϕ̇, ħ being the reduced Planck constant and and e the electron charge. 
In this regime, the  Josephson current time evolution is characterized by the frequency ω_J=2e V/ħ. The junction is thus 
an element that operates in cycles with a natural frequency controllably set by the extrernally injected  current <cit.>.


When the JJ is driven by a combined DC & AC bias with frequency ω, 
the junction dynamics results to  be synchronized with the external drive:  ω_J=n ω,  n being a positive integer. This synchronization results in quantized voltage plateaus, known as Shapiro steps, appearing at ⟨ V⟩_t = n ħω/2e in the time-averaged voltage-current characteristics <cit.>. Shapiro steps arise from a phase-locking process where 
ϕ advances by 2π n during each cycle of the driving field, establishing coherence between the junction phase and the driven field.





Since their observation in superconducting Josephson junctions (SJJ) <cit.>, Shapiro steps have been investigated in various nonlinear systems, ranging from charge density waves <cit.>, colloidal systems <cit.>, superconducting nanowires <cit.>, and ^3He weak-links <cit.>. Shapiro steps in SJJs play a crucial role in metrology, calibration within electronics, and in defining quantum standards <cit.>.

Recent years have witnessed significant progress in

developing and operating circuit architectures employing ultracold atoms <cit.>. 
These 
 atomtronics circuits feature enhanced control and flexibility of their working conditions, 

that may go beyond conventional electronics. 
One of the most promising atomic platforms for implementing such circuits consists of ultracold Fermi gases <cit.>. 
These systems provide the unique opportunity to realize different regimes of superfluidity, ranging from the Bardeen-Cooper-Schrieffer (BCS) limit of weakly-bound fermion pairs to a Bose-Einstein condensate (BEC) of tightly bound molecules, including the intermediate universal, strongly-correlated unitary regime <cit.>. 

In particular, unitary Fermi gases (UFG) serve as prototypical examples of strongly-interacting systems, sharing similarities with other strongly-interacting Fermi systems found in nature, spanning a broad range of energy and length scales, from nuclear matter to neutron stars <cit.>.
Research on quantum transport phenomena within these systems has explored a variety of configurations, including mesoscopic two-terminal setups <cit.>, ring geometries <cit.>, optical lattices <cit.>, and scenarios involving disorder <cit.> or periodic modulation of the harmonic trapping potential <cit.>.






Josephson effects have been observed in both three-dimensional <cit.> and two-dimensional <cit.> unitary superfluids, providing insights into the microscopic properties of these gases. However, despite various theoretical proposals <cit.>, the observation of Shapiro steps in atomic Josephson junctions has remained elusive. 
Recently, Shapiro steps have been predicted to occur by periodically oscillating the position of a thin tunneling barrier within a weakly interacting Bose-Einstein condensate <cit.>. While 
such protocol defines new avenues for exploring Shapiro steps in atomic junctions, in particular for Fermi superfluids, it remains an open question whether strong interactions could affect their dynamics. Interactions are expected to enhance phase-coherence across the junction and potentially reduce decoherence effects associated with the finite compressibility of atomic superfluids <cit.>, but their role on synchronization processes is yet to be understood. Investigating these effects is crucial for deepening our understanding of strongly-interacting systems and for revealing their microscopic behavior under external drives.











 







In this work, we report the observation of Shapiro steps in a periodically driven atomic Josephson junction connecting two strongly-interacting Fermi superfluids. Similar to SJJs, we find that the external driving frequency determines the height of Shapiro steps, while the driving 
amplitude affects their width. Moreover,

we access 
a detailed microscopic understanding of the physical mechanism behind the Shapiro steps by directly measuring the phase dynamics across the junction. By analyzing the current-phase relation under the external drive and observing the phase evolution across different Shapiro steps, we demonstrate 
the synchronization between the junction phase and the external drive. We observe that the Shapiro dynamics is accompanied by 
vortex-antivortex pairs emissions, in agreement with the results of Ref. <cit.>.



Our results showcase the response of a many-body strongly-interacting system to an external drive, resulting in synchronized dynamics.








§ CREATING AND DRIVING THE ATOMIC JUNCTION


We produce strongly-interacting Fermi superfluids of N ∼2e4 atom pairs  by cooling a balanced mixture of the first and third lowest hyperfine states of ^6Li to temperatures ∼ 0.3(1)  T_c, where T_c is the superfluid critical temperature <cit.>. Interparticle interactions are parametrized by 1/k_Fa, where k_F = √(2mE_F)/ħ is the Fermi wave vector, E_F the Fermi energy, and a the s-wave scattering length. We tune a in the vicinity of the Feshbach resonance at 690G <cit.>, enabling the exploration of various superfluid regimes across the BEC-BCS crossover. 
In this work, we focus on two different interaction regimes: a UFG 
with Fermi energy E_F/h = 11.8(4)kHz, and a 

molecular BEC at 1/k_Fa = 3.3(1) with E_F/h = 9.2(3)kHz. 
The gas is confined in the x-y plane using a box-like potential with dimensions L_x × L_y = 125 × 17.5  μm^2, and along the vertical z direction 
by a strong harmonic confinement <cit.>. 
The trap configuration establish a nearly homogeneous in-plane atomic density with the superfluid remaining in the three-dimensional regime. To engineer the atomic JJ, we superimpose a thin repulsive optical barrier onto the atomic sample, which can be displaced along the x-axis, see Fig. <ref> (A). 
The barrier intensity profile and position are dynamically controlled by a digital micromirror device (DMD), whose designed pattern is projected onto the atomic plane using a high-resolution microscope objective <cit.>. The barrier features a FWHM

of 0.9(1)μ m  along the x-direction <cit.>, 
comparable to the healing length of the investigated superfluids, ranging between 0.6-1.7 μm for typical gas parameters. The barrier height, V_0, is set to be larger than the gas chemical potential, μ_0, with V_0 / μ_0= 1.3(1) for the unitary, and 1.2(1) for the BEC JJs, ensuring that the junction operates always in the weak coupling or tunneling regime <cit.>.

We inject currents into our junction by precisely moving the tunneling barrier along specific trajectories controlled with

the DMD <cit.>. 

The in-plane homogeneous density of the superfluid ensures that the injected current I(t) is independent of its longitudinal position, and it is determined only by the instantaneous barrier velocity v(t).
The instantaneous current reads I(t) = v(t) ñ L_y = v(t) N/L_x, with ñ=∫ n_3D(z) dz = N/L_x L_y the planar density. To introduce a biased alternating current, we move the barrier according to the trajectory x(t) = v_DC t + x_ACsin (ω t)  <cit.>, giving rise to the instantaneous injected current 
I(t) = I_DC + I_ACcos(ω t), with I_DC = v_DC N/L_x, and I_AC = ω x_AC N/L_x, (Fig <ref>B). We explore the junction operation under different driving conditions by tuning the trajectory parameters v_DC, x_AC, and ω, as illustrated in Fig. <ref> B for bias (i) and composite current injections (ii).
The junction dynamics can be portrayed with the well-known washboard potential analogy <cit.>, as sketched in Fig. <ref> C. The phase of the JJ follows the same equation of motion as that of a classical particle subject to a viscous force moving in a tilted washboard potential, where the tilt is given by the external current I. For a DC bias above I_c, the particle enters the running-phase regime, where ⟨ϕ̇⟩_t ≠ 0, and the junction operates resistively [Fig. <ref> C-(i)]. For an AC drive, in the overdamped regime, where friction overtakes inertia, the velocity locks to the driving tilt, resulting in a steady but synchronized descent, crossing n potential wells during each modulation cycle [Fig. <ref> C-(ii)]. This resonant motion locks the mean particle velocity to fixed values ⟨ϕ̇⟩_t = nω, leading to distinct, equally-spaced velocity steps, the Shapiro steps.






§ SHAPIRO STEPS IN A UNITARY JUNCTION





We probe the response of our junction by measuring the I - Δμ characteristic curves as a function of the bias current, I_DC, after three driving periods T <cit.>, where T= 2π / ω. Here, Δμ = μ_L - μ_R represents the chemical potential difference between the two reservoirs measured from in-situ density images. 

Curves for I_AC = 0 and I_AC>0 for typical driving parameters are shown in Fig. <ref> A-B. 



In the undriven scenario
, we recover the hallmark behavior of a current-biased JJ, featuring a zero-imbalance plateau for DC currents below I_c <cit.>. 
For the AC-driven JJ instead
, we observe the emergence of distinct plateaus, with number, widths, and heights strongly depending on the driving parameters. 

To characterize the plateau properties, we perform a phenomenological fit of each step in the I - Δμ curve using multiple, patched sigmoidal functions <cit.>, from which we get the height Δμ_n, and the width Δ I_n of each step (inset of Fig. <ref> A).

The step height is in agreement with Δμ_n = n  ħω for various driving frequencies 
independently on I_AC (Fig. <ref> C), as it 
is determined exclusively by the external driving frequency, while the driving strength I_AC affects only the position and width of each step (Fig. <ref> D).



To quantitatively describe the observed I - Δμ characteristics, we exploit the resistively and capacitively shunted junction (RCSJ) circuital model <cit.>. In this model, the junction is set in parallel with resistive and capacitive components, with the total circulating current being: 

    I(t) = I_c sinϕ + ħ Gϕ̇ + ħ C ϕ̈.

The supercurrent contribution is given by the sinusoidal current-phase relation I_c sinϕ, the resistive and capacitive contributions are G ϕ̇ and C ϕ̈, respectively.
Similarly as for SJJ, the Josephson-Anderson relation Δμ = -ħϕ̇ connects the phase dynamics to the chemical potential build-up <cit.>, with  Δμ playing the role of the junction voltage. 

The Stewart-McCumber parameter β_c = I_c C/ ħ G^2 determines the dynamic regime of the junction, discriminating between underdamped (β_c ≫ 1) and overdamped (β_c ≪ 1) 
<cit.>. 
In the latter, the undriven solution is analytically given by Δμ = G^-1√(I_DC^2-I_c^2), 


which we use to fit the DC measurements

to extract I_c and G. 
In the driven scenario instead, the I-Δμ characteristic is evaluated by numerically solving Eq.<ref>. As reported in Fig. <ref> A-B,

the overdamped solutions of the RCSJ model display the characteristic Shapiro steps matching the measured ones. 

We note that even though we 

estimate  
β_c≈ 14(2) for the unitary superfluid, the overdamped RCSJ numerical solutions well capture the observed steps.





In particular, the step widths obtained from the numerical solutions of the overdamped RCSJ model are in agreement with the experimentally measured values (Fig. <ref> D). Note that this result is analogous to the Bessel function behavior of the Shapiro steps found analytically under voltage-driven modulation <cit.>, Δ I_n/I_c = |J_n (V_ac/ħω)|, but quantitatively deviates from it at small driving amplitudes, as expected under current-driven modulation <cit.>.  









§ PHASE DYNAMICS AND SYNCHRONIZATION


To link the 
Shapiro steps with the phase dynamics, we tune our atomic junction to operate in the BEC regime at 1/k_Fa=3.3(1). Here, 

the reduced interactions 
allow to leverage matter-wave interference between the expanding reservoirs to directly measure the relative phase ϕ <cit.>, accessing 
both Δμ-I and ϕ-I characteristics <cit.>, as reported in Fig. <ref> A-B 
for ω = 2 π×175Hz, I_AC/I_c = 0, 
and I_AC/I_c =1.4(1). 
In the non-modulated case for I_DC<I_c, we recover the expected DC Josephson trend <cit.>, i.e. ϕ = arcsin(I_DC/I_c).
The measured ϕ shows larger errorbars only for currents near and above I_c, attributed to the running phase in the voltage state <cit.>.







The driven scenario instead displays two main features: a diminished step contrast  in the I-Δμ curve as compared to the unitary junction, and a well-defined phase ϕ 

within the plateau region even at finite Δμ, presenting significant fluctuations only in the transition regions between plateaus <cit.>.
To support our findings, 

we compare the experimental results 
with numerical simulation performed with a dynamical classical-field method under similar experimental conditions <cit.>.
The numerical simulations quantitatively agree with our results for Δμ-I and ϕ-I characteristics, capturing the smoothness of the Shapiro steps in the experimental data.

The agreement with the numerical simulation evidences that the observed lower contrast of the plateaus in the BEC is inherent to our junction, which nevertheless displays Shapiro steps dynamics in a broad range of modulation parameters

(Fig. <ref> of Supplementary Materials).








We demonstrate the dynamical synchronization between the junction phase and the external driving by monitoring ϕ as a function of time 

in the 0-th and 1-st Shapiro steps (Fig. <ref> C), additionally, the 2-nd Shapiro step is shown in Fig. <ref> of Supplementary Materials.
The interference fringes show an oscillating behavior in both cases, highlighted by 
the measured relative phase 
ϕ (Fig. <ref> D), which showcases the phase-locking and synchronized dynamics induced by the driving.









Note that 
phase locking effects have also been observed in  a superfluid system in a Floquet driven tilted double well potential <cit.>.
Whereas in the 0-th step, the phase oscillates in phase with the injected current without reaching 
the critical values for phase slips, in the 1-st step ϕ reaches this critical limit once every cycle, except for the transient behavior of the first modulation, leading to periodic 
phase slips, marked in Fig. <ref> C-D (ii) by the red arrows.

The phase-slippage process produces a clear spike in ϕ for both the dynamical classical-field simulation and overdamped RCSJ solution, associated to an increase of the errorbar on the measured phase <cit.>.

Depending on the 
junction dimensionality, phase slips manifest as different topological excitations, from a soliton in 1D

<cit.>, to vortex-antivortex pair in 2D <cit.>, and vortex rings in 3D <cit.>. 
The kind of excitations created is determined by which is energetically favorable <cit.>, or arises from 
a multi-step process such as soliton decay into different topological excitations
<cit.>.


In our BEC junction, we directly observe the presence of vortex-antivortex pairs as

localized density depletions in the superfluid after a short time of flight <cit.> (Fig. <ref> E). 
For the 1-st Shapiro step, we observe vortices generated on the near-left of the barrier (i-ii) that subsequently move towards the bulk of the reservoir (iii-iv). After an additional modulation period, a second vortex-antivortex pair is similarly emitted from the barrier region (v).
The dynamical classical-field simulations corroborate this mechanism (Fig. <ref>-<ref> of Supplementary Materials): 
at the end of the first cycle, the 
phase difference accumulated along the junction releases a solitonic excitation from the barrier, which 

then quickly decays into vortex-antivortex pairs, via snake-instability <cit.>.

The periodic emission of vortex-antivortex pairs in Shapiro steps is a proxy of the underlying synchronized phase dynamics of the junction, which undergo n phase slippage processes in the n-th step of each modulation cycle, as pictorially illustrated in the washboard potential analogy of Fig. <ref> C(ii) for n=1. 




By measuring the number of vortex pairs N_d as a function of the injected current (Fig. <ref>), we probe the synchronization also in the strongly-interacting regime.




In the non-modulated case
, vortices are observed to proliferate for I_DC>I_c, signaling a contribution to the junction conductance due to collective excitations, in agreement with previous observations in three-dimensional atomic JJs <cit.>.
For the AC drive instead
, the number of detected vortex pairs shows a step-like trend that closely resembles the Shapiro steps observed in the I- Δμ curve measured under the same driving conditions. 
In particular, in the 2-nd Shapiro step, N_d doubles with respect to the 1-st step, highlighting that, as a result of the underlying synchronized phase dynamics, the number of phase slips is proportional to the step number. 




The dissipative and resistive dynamics giving rise to the Shapiro steps in our atomic JJ is dominated by collective excitations, playing a similar role of quasiparticles produced by Cooper pair breaking in SJJ. In fact, in all explored configurations, the junction always works in the regime of  Δμ≪Δ, with Δ the superfluid gap, where broken pairs are energetically suppressed and the only accessible excitations are collective ones, namely phonons, solitons and vortices <cit.>. 








§ CONCLUSIONS



We have observed Shapiro steps in a current driven atomic JJ composed of two weakly-coupled strongly-interacting Fermi gases. 
Shapiro steps appear as quantized plateaus in the relative chemical potential, the analog of the voltage drop in SJJs, occurring at integer multiples of the external driving frequency. As the driving amplitude increases, the widths of these plateaus exhibit a non-linear dependence, arising 
from the 
interplay between the external drive and the non-linear characteristics  of the atomic Josephson junction.  


We probe the synchronization process at the heart of Shapiro steps by directly accessing the microscopic phase dynamics of the junction featuring periodic phase-slippage processes 
which manifest as vortex-antivortex pairs. 

At the more fundamental level, our results also disclose the interplay between quantum phase coherence and dissipation dynamics in a strongly correlated Fermi system.



Our work opens prospects for a deeper comprehension of how quasiparticles and quantum phase slips contribute to quantum transport, potentially leading to advancements in understanding the complex dynamics within driven superconducting networks <cit.>. 

Our driven atomtronic circuit 
can be indeed exploited

to investigate driven atomic JJ arrays, whether arranged linearly or in annular configurations <cit.>. In turn, this opens the way to simulate important paradigms in quantum synchronization, as provided by Kuramoto-like models <cit.>, by driven atomic JJ in the presence of strong interactions and with disordered couplings between the links <cit.>.
Finally, as Shapiro steps provide the voltage standards in quantum electronics <cit.>, we could leverage the quantized steps to access a direct measurement of chemical potential differences <cit.>. 
This approach would be specifically  valuable in the intermediate regimes of the crossover superfluids, where the equation of state is challenging  <cit.>.








We note that Shapiro steps have been very recently observed also in a driven atomic JJ with weakly interacting bosons. 


























§ ACKNOWLEDGMENTS

We thank Francesco Marino, Ludwig Mathey, and Klejdja Xhani for the discussions and Francesco Scazza for careful reading of the manuscript. G.R. and G.D.P. acknowledge financial support from the PNRR MUR project PE0000023-NQSTI. G.R. acknowledges funding from the Italian Ministry of University and Research under the PRIN2017 project CEnTraL and the Project CNR-FOE-LENS-2023. The Authors acknowledge support from the European Union - NextGenerationEU for the “Integrated Infrastructure initiative in Photonics and Quantum Sciences" - I-PHOQS [IR0000016, ID D2B8D520, CUP B53C22001750006]. This publication has received funding under the Horizon Europe program HORIZON-CL4-2022-QUANTUM-02-SGA via project 101113690 (PASQuanS2.1).
J.A.S. acknowledge financial support from CONAHCyT grant CF-2023-I-72; DGAPA-UNAM-PAPIIT grant IN105724; CIC-UNAM grant LANMAC-2024, and European Community's Horizon 2020 research and innovation program under grant agreement n° 871124.



§ AUTHOR CONTRIBUTIONS



G. D. P., V. P. S., L. A. and G. R. conceived the experiment. 
G. D. P., D. H.-R., N. G., M. F. F. and G. N. performed the experimental work, acquired the data and carried out the data analysis. 
V. P. S. performed numerical simulations. 
L. A and G. R. supervised the work.
All authors contributed extensively to the discussion and interpretation of the results, and took part in writing the manuscript.

















Science


















§ SUPPLEMENTARY MATERIALS
 










 §.§ Experimental method



  §.§.§ Gas preparation

We initially realize a superfluid atomic Fermi gas by evaporatively cooling the first and third lowest hyperfine state of ^6Li atoms in a red-detuned, cigar-shaped harmonic trap at 690G. To obtain a superfluid in the BEC regime, during the last part of the evaporative cooling procedure we sweep the magnetic field at 633G, where the molecular scattering length is a_M = 0.6  a = 1029  a_0. After the production of the superfluid, we adiabatically ramp up in 100ms a TEM_0,1-like laser beam, that confines the atoms in the ẑ direction creating a harmonic potential, and the hard wall potential with a rectangular shape in the x-y plane created with the DMD. Both these two potentials are realized with a blue-detuned laser at 532nm. The vertical harmonic confinement has a trap frequency ω_z =2π×685(5)Hz at unitarity and ω_z =2π×416(4)Hz in the BEC regime. Subsequently, we adiabatically ramp down the initial cigar-shape potential in 100ms, and we let the system equilibrate for at least 50ms before ramping up the tunneling barrier. 
The superfluid produced in such a hybrid trap is homogeneous in the x-y plane since the residual harmonic confinement arising from the TEM_0,1-like and the curvature of the magnetic field corresponds to 2.5Hz, therefore negligible for the dynamics studied in this work. 
The Fermi energy of the system in the hybrid trap is computed directly from a Thomas-Fermi approximation yielding <cit.>: E_F = ħ^2 k_F^2/2m = 2 ħ√(ħπω_z N/m A), where ħ is the reduced Plank constant, m the mass of a ^6Li atom, and A = L_x × L_y, the area of the box trap. 
The chemical potential of the unitary superfluid is calculated as μ_0/h = ξ^3/4 E_F/h = 5.9 (2)kHz, where ξ≃ 0.37 is the Bertsch parameter; whereas, for the molecular BEC, we calculate it as μ_0/h = (3 π N ħ^2 ω_z a_M / 2 A √(m) )^2/3 / h = 1.20(5)kHz <cit.>.

The healing length in the BEC regime is ξ_L = 0.59(1)μ m, while in the UFG the typical length scale of the superfluid is given by 2π/k_F ≃1.67(3) μ m. 
To produce the atomic JJ, we ramp up the barrier in its initial position in 10ms and wait for 30ms before starting the movement. The homogeneity of the cloud before ramping up the barrier ensures the same density and chemical potential in the two reservoirs at t = 0ms.



The chemical potential difference Δμ 
is extracted by taking an in situ image of the cloud after the barrier movement and  

computing μ_L (μ_R) from the measured number of atoms N_L (N_R) in the left (right) reservoir. This way of computing Δμ is analogous to a measurement of the time average ⟨Δμ⟩_t. In fact, as it is shown in Fig. <ref> A, the instantaneous density accumulation and depletion, created in the vicinity of the barrier because of its motion, propagate in time away from the barrier position inside the two reservoirs, traveling at the speed of sound (measured as in Ref. <cit.>). 
Consequently, the time evolution of the instantaneous chemical potential difference at the barrier position, Δμ_b (t), is mapped in different positions in space. 
The black horizontal dashed line in the figure indicates the time at which the measurement is performed, corresponding to three modulation periods. At this time, the initial density modulations have traveled up to the edge of the cloud covering the full reservoirs. For this reason, the time average value of Δμ_b (t) is mapped in the measured value of the global chemical potential difference, Δμ, between the two reservoirs.
This reasoning is expected to be less valid for higher frequencies, for which the period is shorter and the density modulations do not reach the edges at the end of the three cycles, leading to a 
lower value of the measured Δμ. 
Given the lower speed of sound in the BEC regime, the decrease happens at lower frequencies with respect to the UFG regime. 
We note that all the measurements of Δμ reported in this work have been performed by imaging the cloud after three modulation periods. This value is the longest number of cycles allowed by the size of the junction, for the range of barrier velocity spanned. For each ω, we also acquired a corresponding DC curve with I_AC=0, moving the barrier at constant velocity for a time of 3  T.





  §.§.§ Barrier characterization

The barrier intensity, shape, and position are controlled with the DMD. The effective movement of the barrier is obtained by creating a sequence of different light patterns with the DMD, with the barrier position following the desired equation of motion 

and shining it on the atomic plane through the high-resolution microscope objective with a frame rate of 12.5kHz. 
Even though the DMD-made barrier movement is intrinsically discretized, the small size of a single DMD-pixel imaged on the atomic plane, 0.25  μm, ensures the smoothness of the barrier movement. 
We characterize the properties of the barrier movement by monitoring the position of the barrier as a function of time in the case of Fig. <ref> B. The barrier position is obtained by performing a Gaussian fit of the density depletion of the in situ images as a function of time. The extracted positions are indicated by the white points in Fig. <ref> A. We fit the trajectory with the equation of motion of the barrier x(t) = v_DC t + x_ACsin (ω t). The fitted parameters are compatible within 1% error for the value of v_DC and 0.5% error in the case of x_AC and ω with respect to the set values.
The Fourier transform (FT) of the observed trajectories respect the ones with the same parameters but x_AC = 0μ m, x̅_b(t), shows an almost single-frequency behavior, with a small contribution of the second and third harmonics, as it is shown in Fig. <ref> C, confirming that our protocol provides a clean and monochromatic AC drive. The measurements reported in this work have been performed with values of x_AC in the range 1-6  μm, limited by the finite resolution of the DMD-projecting setup and the frame rate of the DMD, respectively.


The intensity of the profile created to produce the tunneling barrier is calibrated with a phase imprinting method <cit.>. 
We create a small JJ of dimension 37.5 × 17.5  μm^2 and imprint a phase on one of its reservoir

by applying a homogeneous optical potential U_0 for a short time interval Δ t. For Δ t < h/μ, we operate an 

almost pure phase excitation, with the phase imprinted by the light on the illuminated reservoir given by ħϕ = U_0 Δ t. To calibrate the potential U_0, we measure ϕ from the interference arising in the time-of-flight (TOF) expansion of the junction, 

using the phase of the non-imprinted reservoir as a reference. 
Figure <ref> A-B shows an example of an interference pattern, which we fit to extract ϕ with a sinusoidal function. 
By repeating this process for different values of the imprinting pulse power, we extract a calibration for the imprinted potential U_0 (Fig. <ref> C).




To characterize the barrier height and width we image the DMD-produced barrier intensity pattern on a service camera in the DMD-projecting optical system. We acquire images of the barrier in different positions throughout the area of the JJ and extract their height and size via a Gaussian fit along the x- direction. To account for the finite resolution of the imaging system after the service camera, before performing the fit, we convolve the images with a Gaussian of FWHM= 0.63  μm, well-approximating the Point Spread Function (PSF) of the microscope objective. 

The barrier properties are barely varying while changing their position, providing the average values reported in the main text for V_0 and FWHM. In particular, to obtain the value of the barrier height V_0 we compare the fitted barrier height from the image with the pattern used to create the phase imprinting and employ the calibration factor to convert it into energy units. 






  §.§.§ Phase and vortex detection

The relative phase between the two reservoirs is measured from the interference fringes arising after a TOF expansion of 2ms, after abruptly switching off all the confinements. The short TOF employed for this measurement ensures that interference fringes appear only close to the barrier region (Fig. <ref>), allowing for a precise measurement of the relative phase at the junction ϕ. In particular, to quantitatively extract ϕ we restrict to a region around the instantaneous position of the barrier (dashed red rectangle area in the figure) and fit the density profile integrated along the y-direction with a sinusoidal function. We constrain the wavelength of the fringes to the same value for all the acquired interferograms, as this quantity depends only on the TOF. We fix the best fringe wavelength as the average wavelength from a preliminary sinusoidal fit to all the data acquired under different conditions. In Fig. <ref> B-C we report the measured phase evolution in the second Shapiro step, under the same conditions as in Fig. <ref> of the main text. 
Both in the experimental data and in the overdamped RCSJ model numerical solutions, two phase-slips processes are evident in each modulation cycle. 


 


To detect vortices in the BEC regime we employ a TOF technique as well, but, to avoid the simultaneous presence in the images of interference fringes, we ramp down the barrier in 0.24ms, wait 4ms and then perform the TOF expansion. To increase the contrast of the vortices, we ramp down the DMD-made potential during the first 1ms of TOF expansion. 



To observe the presence of vortices in the UFG regime we add to the TOF expansion a fast sweep of the Feshbach magnetic field to map the system in a molecular BEC, where vortices are visible as clear holes in the density <cit.>.
In particular, at the end of the movement of the barrier
, we wait 1ms 
before rapidly removing the barrier from the system in 1.5ms 
and then we wait 1.5ms 
before performing the TOF procedure, 

starting the magnetic field ramp 
after the end of the barrier movement.
In both BEC and UFG regimes, the waiting time between the barrier ramp down and the TOF expansion is expected not to affect significantly the number of vortices in the system. 
We note that, since N_d of Fig. <ref> of the main text is measured after three modulation cycles, we would expect to see a higher number of vortex-antivortex pairs, corresponding to the periodic phase-slippage process. 
However, the finite size of the junction along the y-direction facilitates the vortices to escape from the bulk density, so that the measured N_d most likely corresponds to the number of vortex-antivortex pairs emitted during the last modulation cycle.
We note that with a similar protocol as for detecting vortices, employing a fast sweep during the TOF, the relative phase at the junction could be measurable also at unitarity from the interferograms <cit.>. However, this technique provides more noisy interferograms, and the lower signal/noise together with the presence of vortices makes the extraction of ϕ harder. For this reason, we decided to perform the phase measurement in the BEC regime. 





 §.§ Circuital model simulations





We study the expected behavior of Δμ and of ϕ in the framework of the RCSJ model, in which the phase evolution is described by Eq. (<ref>) of the main text. We solve numerically the equation in the overdamped regime (β_c ≪ 1), neglecting the capacitance, i.e. the second derivative term. 
In this limit, the I_AC=0 case solution Δμ = G^-1√(I_DC^2-I_c^2) is used to extract the values of G and I_c by fitting the observed value of Δμ as a function of I_DC. These values are then used in order to simulate the dynamic in the I_AC≠ 0 case. Figure <ref> displays the results for values of the parameters as the case plotted in Fig. <ref>A-B of the main text, using the Josephson-Anderson relation Δμ = -ħϕ̇. In particular, we extract the expected average value of Δμ from the time average value of the phase derivative ⟨ϕ̇⟩ for a total evolution time of 10periods. Figure <ref>C shows the expected behavior in the absence of modulation (I_AC=0, blue line) and the predicted shape of Shapiro steps for a I_AC=1.3  I_c. In Fig. <ref>D we plot the corresponding phase at the end of 3 complete periods for the same cases of I_AC=0 and I_AC=1.3  I_c. This plot shows a similar trend to the one shown in Fig. <ref>B of the main text, apart for the regions in the transition between the steps in the modulated case, where this quantity shows a not well-defined behavior. In particular, Fig. <ref>E, shows the value of the distance between the phase at a given value of I_DC resulting from different initial values of the phase ϕ_0: Δϕ^2 (I_DC) = ∑_i,j ( ϕ(I_DC,ϕ_0 = ϕ_0,i) - ϕ(I_DC,ϕ_0 = ϕ_0,j))^2. This value diverges in this transition region.

Under a voltage modulation, the Shapiro step half-widths follow the well-known relation Δ I_n/I_c = |J_n(V_AC/ħω)|. However, quantitative deviations occur under current modulation, especially for the first few steps, as shown in Fig. <ref>.







 §.§ Phase Fluctuations


In our system, we refer to the phase fluctuations in terms of the fluctuations of the relative phase given small deviations in the initial conditions of ϕ, and not to stochastic noise terms in the RCSJ model.
These fluctuations in the initial conditions are natural in our experimental setup and modify the measurement of the phase we perform after TOF, as discussed in the previous section.


To quantify the stability of the phase measurements we make use of the overdamped RCSJ model. We compare the phase evolution of from 100 simulations of Eq. (<ref>) given different initial conditions ϕ(t=0) = ε, where ε follow a uniform distribution in the domain [-2π/10, 2π/10]. The ε range is chosen based on the fluctuations observed without any external current injected in the JJ. We compare our experimental data to the simulation results after three modulation cycles, see Fig. <ref>.
Let us note that the data points in Fig. <ref> A-B display a larger error bar near the transitions between the plateau regions (B), and near and above I_c (A), also visible in the measured ϕ as larger errorbars. 

This behavior is recovered from the numerical simulations shown in Fig. <ref> D-E, where the variance measured over the ensemble of different initial conditions is non-zero only in the transition regions while maintaining a well-defined value in the central region of the plateaus.



 §.§ Classical-field simulation method


To simulate the dynamics of ^6Li_2 molecules on the BEC side we use classical-field dynamics within the truncated Wigner approximation <cit.>.
The condensate is described by the Hamiltonian



    Ĥ = ∫d r[  ħ^2/2m∇ψ̂^†( r) ·∇ψ̂( r)  + V( r) ψ̂^†( r)ψ̂( r)   + g/2ψ̂^†( r)ψ̂^†( r)ψ̂( r)ψ̂( r)],



where ψ̂( r) and ψ̂^†( r) are the bosonic annihilation and creation field operators, respectively. 
The 3D interaction parameter is given by g=4π a_D ħ^2/m_D, where a_D is the molecular s-wave scattering length and m_D = 2 m is the molecular mass.
Following the experiments we choose: a_D = 0.6 a=1029  a_0, where a_0 is the Bohr radius, V( r) = V(z) = m_D ω_z^2 z^2/2 the harmonic trapping potential, with ω_z = 2 π× 416Hz being the trap frequency in the z direction. 



For our simulations, we consider the same dimensions of the experimental Josephson junction, by mapping the real space on a lattice system of  253 × 37 × 8 sites with the lattice discretization length l= 0.5 μmm.
We note that l  is chosen such that it is to be comparable or smaller than the healing length ξ_L = ħ/√(2mgn_0) and the de Broglie wavelength, where n_0 is the density <cit.>. 
In the classical-field representation, we replace the operators ψ̂ in Eq. <ref> and in the equations of motion by complex numbers ψ. 
We sample the initial states ψ(t=0) in a grand canonical ensemble of temperature T and chemical potential μ via a classical Metropolis algorithm. 
For all simulations, we use T= 40  and adjust μ such that the total atom number N is about 2 × 10^4.  
Each initial state is propagated using the classical equations of motion. 

To create a Josephson junction we add the term ℋ_ex = ∫d r  V(x,t) n( r, t), 
where n( r, t) = |ψ( r , t)|^2 is the local density and V(x, t) is the barrier potential of the form 


    V(x,t)  = V_0 (t) exp[-  2( x-x_0- x(t) )^2/w^2].


Here, V_0(t) is the height and w is the width of the barrier.  
The initial location of the barrier potential is fixed at x_0= 32, while  x(t) is the driving term. 
As the barrier size is at the limit of the experimental resolution, the experimental estimation of its height is particularly challenging. 
For this reason, we perform the numerical simulation by fixing the value of the barrier size to be compatible with the confidence range of the experimental one and scan V_0/μ_0 to obtain the best agreement with the experimental data of AC and DC drive reported in Fig. <ref>.
We choose w=0.7, and ramp up V_0(t) linearly to V_0/μ_0=1.45 over 100 and then wait for 50 to achieve equilibrium, 
where μ_0 is the chemical potential averaged over the z direction. 
This creates a weak link at x_0 by separating the cloud into two subsystems 
(referred to as the left and right reservoirs).  
The driving term is given by  


    x(t) = v_DC t + x_ACsin(ω t),


where v_DC is the DC barrier velocity, x_AC is the AC driving amplitude and ω is the AC driving frequency <cit.>. 
Similarly to the experiment, we write the current as 
I(t) = I_DC + I_ACcos(ω t), where I_DC = v_DC I_c/v_c and I_AC = x_ACω I_c/v_c, with I_c being the critical current and v_c the critical velocity.  
For various driving parameters, we analyze I-Δμ characteristic curves as a function of the bias current I_DC after three driving periods. 
Δμ= μ_R - μ_L is the chemical-potential difference between the left (μ_L) and right (μ_R) reservoir, which is determined using the total atom number in each reservoir. 
In Fig. <ref>A of the main text, we show the simulation results of I-Δμ curves for undriven (I_AC/I_c=0) and driven (I_AC/I_c=1.4 and ω=2 π× 175) junctions.
We perform simulations of the AC response of the junction as a function of the driving frequency ω and report these results in Fig. <ref>B, which are in excellent agreement with the experimental findings.
To determine the height of Shapiro steps we fit I-Δμ curves with sigmoid functions as in the experiments. 



To determine the ϕ-I characteristic of the junction we calculate the phase difference in the vicinity of the barrier 


    ϕ = ϕ(x-2l) - ϕ(x+2l),


where x(t) is the dynamic location of the barrier at time t. 
In Fig. <ref> we show the time evolution ϕ(t), averaged over many samples, at varying I_DC/I_c for the AC-driven system. From the phase change at the end of driving we obtain the ϕ-I curves, shown in Fig. <ref>B of the main text. 














  §.§.§ Vortex dipole creation: phase slip process


As mentioned in the main text, depending on the junction dimensionality, phase slips can manifest as different topological excitations <cit.>, in particular for our experimental conditions as vortex-antivortex pair in two dimensions, and arise from a multi-step process as soliton decays into a different topological excitation <cit.>.
Our classical-field simulations corroborate this mechanism as shown in Fig. <ref>-<ref>, where a phase difference is accumulated along the extended junction at the end of the second (Fig. <ref>) and third cycle (Fig. <ref>). After the solitonic excitation is released from the barrier, the soliton decays into multiple vortex-antivortex pairs possibly via a snake-instability <cit.>. It is important to note that the depinning of the soliton and the subsequent decay process happens on a timescale faster than the modulation period. The critical wavevector at which the snake instability occurs near unitarity <cit.>, and at unitarity reaches k_c∼ 0.93  k_F and decreases in the BEC side of the crossover according to the formula k_c ∼√(2)ξ_L^-1. 
The critical wavelength is therefore λ_c=2π/k_c = √(2)πξ_L≈ 2.5(2) μm in the BEC, and λ_c=2π/k_c= 1.8(2) μm in the UFG.
In particular, our junction transverse size is large compared to the critical length L_y/ λ_c ≈ 7, allowing multiple "snake" oscillations along the junction transverse direction.




 §.§ Shapiro steps in the BEC regime







In this section, we report the characterization of the Shapiro step properties for the BEC junction (1/k_Fa = 3.3±0.1), together with a detailed description of the phenomenological fitting procedure, employed to analyze all the I - Δμ curves, independently from the interaction regime.

The DC current-voltage data are fitted with the analytical solution of the overdamped RCSJ model, namely with Δμ = √( I_DC^2 - I_c^2)/G (dark blue dashed line of Fig. <ref> A), keeping both I_c and G as free parameter. 
On the other hand, to quantitatively characterize the Shapiro step height and width, we fit the AC I - Δμ curve with as many independent sigmoid functions as the number of visible steps (solid lines in Fig. <ref> A). In particular, we fit the n-th step with the function:

    Δμ = Δμ_n^rel/(1+ exp( -I_DC-I_n/Γ)))+ Δμ_n-1^rel,

where Δμ_n^rel is the relative height of the n-th step and I_n its position. From these parameters we compute the n-th step height as Δμ_n = ∑_i=0^nΔμ_i^rel and its width as Δ I_n = I_n+1-I_n. The step height characterization as a function of the modulation frequency in the BEC regime is reported in Fig. <ref>. Here, the experimental data (filled symbols) are compared with the results of classical-field numerical simulations (empty symbols), performed as described in the dedicated section, and analyzed with the same fitting procedure as described above. The numerical simulations are in very good agreement with the experimental results, both confirming the presence of Shapiro steps at Δμ = n ħω in the BEC regime in the explored range of frequency. In particular, in both cases, we observe a reduction of the measured step height for increasing modulation frequency, which could be due to the coupling of the drive with trap excitation along the z-direction. In fact, for AC drive frequencies approaching the vertical trap frequency, we expect the Shapiro dynamics to mix with trap excitations, possibly redistributing the injected energy between the two contributions. We observe a similar reduction also for the junction at unitarity at ω≃ 2 π× 500Hz, compatible with this scenario given the higher trap frequency in this regime. As already mentioned above, a reduction of the step height at higher modulation frequency could also arise from the finite speed of propagation of the generated Δμ at the barrier, given by the speed of sound. For high ω, the space propagated by the Δμ pulses during three modulation cycles can be shorter than the junction size, so that the measured Δμ can be lower than the stationary value. 



In Fig. <ref> A, we report also the overdamped RCSJ numerical solutions for the AC I - Δμ (red and blue dashed lines), obtained by using the value of I_c and G as extracted from the fit of the DC curve. The numerical solutions qualitatively represent the measured steps, despite a smoother trend visible both in the experimental data and in numerical simulations, visible also in the data reported in Fig. <ref> A of the main text. Such a discrepancy is most likely given by larger conductance at BEC (h C = 6.5(6)s), with respect to unitarity (h C = 1.8(2)s).


The conductance is here computed as the inverse of the charging energy E_c = 1/C = ∂μ_L/N_L + ∂μ_R/N_R≃ 4 ∂μ_0/N, and is therefore proportional to the compressibility of the gas. 
In the BEC regime, the larger compressibility, with respect to unitarity, allows for density excitations at higher contrast to propagate in the junction as a consequence of the AC drive, making the amplitude of the order parameters in the two reservoirs non-spatially homogeneous. As a result, the dynamic of the junction is no longer well describable in terms of uniquely its relative phase ϕ, and the RCSJ model becomes less accurate. We note that density excitations are observable also in the unitary junction as a consequence of the AC drive (see Fig. <ref>), but the strong interactions reduce their contrast and thus keep the overdamped RCSJ model a good approximation. 


