Josephson Critical Currents and Related Effects in Ultracold Atomic Superfluid Sytems

Abstract

The Josephson and Proximity effects play a pivotal role in the design of superconducting devices for the implementation of quantum technology, ranging from the standard $Al$ based to the more exotic twisted high-$T_c$ junctions. Josephson critical currents have been recently investigated also in ultracold atomic systems where a potential barrier acts as a weak link. The unifying feature of the above systems, apart from being superconducting/superfluid, is the presence of spatial inhomogeneity, a feature that has to be properly taken into account in any theoretical approach employed to investigate them. In this work, we review the novel (dubbed LPDA for Local Phase Density Approximation) approach based on a coarse graining of the Bogoliubov–de Gennes (BdG) equations. Non-local and local forms of this coarse graining were utilized when investigating Proximity and Josephson effects. Moreover, the LPDA approach was further developed to include pairing fluctuations at the level of the non-self-consistent t-matrix approximation. The resulting approach, dubbed mLPDA ($modified$ LPDA), can be used whenever inhomegeneity and fluctuations effects simultaneously play an important role.

1. Introduction

When investigating inhomogeneous superconducting/superfluid systems, the most natural tool which comes to one’s mind is probably the BdG equations [1] (or equivalently the Gorkov’s ones [2,3]). This approach, which allows for a straightforward description of the physical system at hand, in its most known form, consists in solving Schroedinger-like equations for two-component single-particle eigenfunctions. At the cost of a long evaluation time and the need for a huge memory space, these equations return a complete set of orthonormal eigenfunctions in terms of which all physical quantities of interest can be evaluated. For these reasons, over the years, when dealing with interfaces or inhomogeneous systems in general, alternative equations such as Usadel [4,5] and Eilenberger [5,6] equations have been solved in place of BdG, at the cost of limiting the validity of the obtained results in a restricted region of the coupling-temperature phase diagram.

In 2014, Simonucci and Strinati [7] developed a new approach using as a starting point the self-consistent equation for the order parameter $\Delta \left(\mathbf{r}\right)$ written in terms of the finite temperature Green functions (solving this equation is equivalent to solving BdG equations [8]). They applied a double-coarse-graining procedure to the equation, assuming that the magnitude of the order parameter $|\Delta (\mathbf{r}\left)\right|$ varied on a larger length scale than that where the gradient of its phase varied. Doing so, they obtained a highly non-linear integral equation for $\Delta \left(\mathbf{r}\right)$, which allows for the investigation of systems with interfaces with no a priori assumption on the boundary conditions (see for example [9]). Performing a further approximation on that equation, which was dubbed NLPDA (Non-Local Phase Density Approximation) [10], they obtained a highly non-linear second-order differential equation for $\Delta \left(\mathbf{r}\right)$ named LPDA, for it entails a Local Phase Density Approximation. Arriving eventually at an equation like this one has been a long standing problem [11,12,13]. This equation could in principle be applied in any region of the coupling-temperature phase diagram (in fact, it reduces to the Ginzburg–Landau [14] and Gross–Pitaevskii [15] equations in the appropriate limits), thereby allowing for faster results than BdG (a factor of ${10}^5$ was estimated in ref. [10] when investigating the isolated vortex problem) with no a priori limitation in that regard (as instead is for Usadel or Eilenberger equations). The only inherent limitation of the LPDA approach is linked to the coarse-graining procedure itself: in ref. [10], it was established that the granularity size is of the order of the Cooper pair size ${\xi }_{pair}$, therefore limiting the reliability of LPDA results at low temperature in the BCS (Bardeen–Cooper–Schrieffer) limit, where ${\xi }_{pair}$ is of the same order of the healing length ${\xi }_{phase}$ (which, by definition, corresponds to the length scale over which the order parameter varies).

In addition, we recall that BdG equations rely on a mean-field decoupling which is expected not to be valid at finite temperature not only in the BEC (Bose–Einstein condensation) limit but also at unitarity, where the Cooper pair size is of the same size as the mean interparticle distance. As a consequence, in that region of the coupling-temperature phase diagram, the outcomes of the LPDA equation, as well as the NLPDA equation, are not quantitatively reliable. For this reason, in ref. [16], we developed the mLPDA approach, where pairing fluctuations at the level of the non-self-consistent t-matrix approximation are added directly to the LPDA equation in order not to lose its inherent numerical advantages with respect to the BdG. This inclusion allows for an improvement in both the pairing gap [17] and critical temperature [18] with respect to the mean-field estimations. The mLPDA approach was validated in ref. [19], where it was applied to investigate the experimental systems realized in refs. [20,21], obtaining good quantitative agreement with the experimental data for the Josephson critical current both at low temperature for various couplings across the BCS-BEC crossover and at unitarity for T varying from 0 to the critical temperature $T_c$.

In this work, we provide an overview of the NLPDA, LPDA, and mLPDA approaches and summarize the most significant results we obtained when investigating the Josephson effect. The usefulness of these methods, apart from the significant speedup in the calculations with respect to BdG equations, relies on the fact that the only underlying assumption is the smoothness of the gap profile, and for this reason, they can be applied to the whole coupling-temperature phase diagram in the broken-symmetry phase with no limitations. In particular, the mLPDA approach, allowing for the description of inhomogeneity and pairing fluctuations in a superfluid system within a single approach and with no need of external parameters, poses itself as a powerful method to investigate a plethora of physical mechanisms, such as the intrinsic depairing limit of the critical current in superconductors [22,23], the paraconductivity probed by the Josephson effect [24], the giant Proximity effect in high-$T_c$ superconductors [25,26], and the combined features of Proximity [27] and Josephson [28,29] effects in a SNS junction.

The paper is organized as follows. Section 2 introduces the NLPDA, LPDA, and mLPDA approaches, Section 3 presents the main numerical results obtained, and Section 4 provides our conclusions. Throughout the paper, we consider balanced spin populations and set $\hslash =1$.

2. Theoretical Approach

In this section, we briefly summarize the main features of the NLPDA, LPDA and mLPDA approaches. Both the specific form of the equations used when investigating Proximity and Josephson effects and the details of the numerical procedures can be found in the original works [9,16,30]. Here, instead, we intend to provide an overview of these approaches focusing on their building blocks, pointing out, at the same time, the main motivations for their development.

Throughout this work, we will span the BCS-BEC crossover using the dimensionless coupling parameter ${\left(k_F{a}_F\right)}^{-1}$, where $a_F$ is the scattering length of the two-fermion problem, and $k_F={\left(3{\pi }^{2}n\right)}^{1/3}$ is the Fermi wavevector associated with (uniform) particle density n. In the weak-coupling (BCS) regime, $a_F<0$ and the coupling parameter ${\left(k_F{a}_F\right)}^{-1}\le -1$, while in the strong-coupling (BEC) regime, $a_F>0$ and ${\left(k_F{a}_F\right)}^{-1}\ge +1$, across the unitary limit where $|a_F|$ diverges and ${\left(k_F{a}_F\right)}^{-1}=0$.

2.1. The NLPDA Approach

The NLPDA equation for the local order parameter, as derived in ref. [7], reads

$$-{\displaystyle \frac{m}{4\pi a_F}}\Delta \left(\mathbf{r}\right)=\int {\displaystyle \frac{\mathrm{d}\mathbf{Q}}{{\pi }^3}}{\mathrm{e}}^{2i\mathbf{Q}·\mathbf{r}}\Delta \left(\mathbf{Q}\right)K^{\mathbf{A}}\left(\mathbf{Q}\right|\mathbf{r}),$$

(1)

where $\Delta \left(\mathbf{Q}\right)$ is the Fourier transform of $\Delta \left(\mathbf{r}\right)$, and the kernel $K^{\mathbf{A}}\left(\mathbf{Q}\right|\mathbf{r})$ is defined by

$$K^{\mathbf{A}}\left(\mathbf{Q}\right|\mathbf{r})\equiv \int {\displaystyle \frac{\mathrm{d}\mathbf{k}}{{\left(2\pi \right)}^3}}\left[{\displaystyle \frac{1-2f_F\left(E_{+}^{\mathbf{A}}(\mathbf{k};\mathbf{Q}|\mathbf{r})\right)}{2E^{\mathbf{A}}(\mathbf{k};\mathbf{Q}|\mathbf{r})}}-{\displaystyle \frac{m}{{\mathbf{k}}^2}}\right],$$

(2)

where m is the fermion mass, $f_F\left(E\right)={(e^{E/\left(k_{B}T\right)}+1)}^{-1}$ the Fermi function at temperature T ($k_B$ being the Boltzmann constant), and

$$E_{\pm }^{\mathbf{A}}(\mathbf{k};\mathbf{Q}|\mathbf{r})\equiv \sqrt{{\left({\displaystyle \frac{{\mathbf{k}}^2}{2m}}+{\displaystyle \frac{{\mathbf{Q}}^2}{2m}}-\overline{\mu }\left(\mathbf{r}\right)-{\displaystyle \frac{\mathbf{A}\left(\mathbf{r}\right)}{m}}·\mathbf{Q}\right)}^2+{|\Delta \left(\mathbf{r}\right)|}^2}\pm {\displaystyle \frac{\mathbf{k}}{m}}·(\mathbf{Q}-\mathbf{A}\left(\mathbf{r}\right)),$$

(3)

$$2E^{\mathbf{A}}(\mathbf{k};\mathbf{Q}|\mathbf{r})\equiv E_{+}^{\mathbf{A}}(\mathbf{k};\mathbf{Q}|\mathbf{r})+E_{-}^{\mathbf{A}}(\mathbf{k};\mathbf{Q}|\mathbf{r}),$$

(4)

where $\overline{\mu }\left(\mathbf{r}\right)\equiv \mu -V\left(\mathbf{r}\right)-{\mathbf{A}}^2\left(\mathbf{r}\right)/\left(2m\right)$ is the local chemical potential (where V is a single-particle external potential, and $\mathbf{A}$ is the vector potential). The thermodynamical chemical potential $\mu$ is obtained by supplementing the NLPDA equation with the density equation

$$n\left(\mathbf{r}\right)=\int {\displaystyle \frac{d\mathbf{k}}{{\left(2\pi \right)}^3}}\left\{1-{\displaystyle \frac{{\xi }^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})}{E^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})}}\left[1-2f_F\left(E_{+}^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})\right)\right]\right\},$$

(5)

where

$$\begin{array}{cc}\hfill {\xi }^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})& ={\displaystyle \frac{{\mathbf{k}}^2}{2m}}-\mu \left(\mathbf{r}\right)+{\displaystyle \frac{1}{2m}}{\left(\nabla \phi \left(\mathbf{r}\right)-\mathbf{A}\left(\mathbf{r}\right)\right)}^2,\hfill \\ \hfill E^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})& =\sqrt{{\xi }^{\mathbf{A}}{\left(\mathbf{k}\right|\mathbf{r})}^2+{|\Delta \left(\mathbf{r}\right)|}^2},\hfill \\ \hfill E_{+}^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})& =E^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})+{\displaystyle \frac{\mathbf{k}}{m}}·\left(\nabla \phi \left(\mathbf{r}\right)-\mathbf{A}\left(\mathbf{r}\right)\right),\hfill \end{array}$$

(6)

where $\mu \left(\mathbf{r}\right)=\mu -V\left(\mathbf{r}\right)$ and $\phi \left(\mathbf{r}\right)$ is the phase of the order parameter ($\Delta \left(\mathbf{r}\right)=|\Delta \left(\mathbf{r}\right)|e^{i\phi \left(\mathbf{r}\right)}$).

In ref. [9], Equations (1) and (5) are simultaneously solved for an infinite homogeneous SN interface by setting $\mathbf{A}\left(\mathbf{r}\right)=0$ and letting the scattering length change across the interface. The coupling parameter varies from ${\left(k_F{a}_F\right)}_L^{-1}$ on the far left of the system to ${\left(k_F{a}_F\right)}_R^{-1}$ to the far right. The critical temperatures associated with those couplings satisfy the relation $T_c^L>T_c^R$ so that in the temperature interval $[0,T_c^R]$, the whole system is superfluid, while for T in $[T_c^R,T_c^L]$, the right region is normal. In that work, we obtain the penetration depth ${\xi }_N$, which describes the length scale over which the order parameter penetrates into the N region, at various temperatures in the interval $[0,T_c]$ and for different combinations of the couplings ${\left(k_F{a}_F\right)}_S^{-1}$ and ${\left(k_F{a}_F\right)}_N^{-1}$, extending the validity of the results obtained in ref. [31]. However, in the theoretical approach, an important feature is missing, namely, the flow of a superfluid current. In condensed matter experiments, indeed, a well-known procedure to obtain the penetration depth is to measure the critical current $I_c$ for various widths of the N region of an SNS junction (see, for example, ref. [28]). The physical constraint of the continuity equation is difficult to implement within the NLPDA approach from a numerical perspective, due to the non-locality of the equation, but it is needed when investigating superfluid currents. Accordingly, we resort to the LPDA approach.

2.2. The LPDA Approach

The LPDA equation is obtained expanding the kernel in Equation (2) in powers of $\mathbf{Q}$ up to the quadratic order [7] and reads

$$-{\displaystyle \frac{m}{4\pi a_F}}\Delta \left(\mathbf{r}\right)={\mathcal{I}}_0\left(\mathbf{r}\right)\Delta \left(\mathbf{r}\right)+{\mathcal{I}}_1\left(\mathbf{r}\right){\displaystyle \frac{{\nabla }^2}{4m}}\Delta \left(\mathbf{r}\right)-{\mathcal{I}}_1\left(\mathbf{r}\right)i{\displaystyle \frac{\mathbf{A}\left(\mathbf{r}\right)}{m}}·\nabla \Delta \left(\mathbf{r}\right),$$

(7)

with the notation

$${\mathcal{I}}_0\left(\mathbf{r}\right)\equiv \int {\displaystyle \frac{\mathrm{d}\mathbf{k}}{{\left(2\pi \right)}^3}}\left[{\displaystyle \frac{1-2f_F\left(E_{+}^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})\right)}{2E\left(\mathbf{k}\right|\mathbf{r})}}-{\displaystyle \frac{m}{{\mathbf{k}}^2}}\right],$$

(8)

$${\mathcal{I}}_1\left(\mathbf{r}\right)\equiv {\displaystyle \frac{1}{2}}\int {\displaystyle \frac{\mathrm{d}\mathbf{k}}{{\left(2\pi \right)}^3}}\left\{{\displaystyle \frac{\xi \left(\mathbf{k}\right|\mathbf{r})}{2E{\left(\mathbf{k}\right|\mathbf{r})}^3}}[1-2f_F\left(E_{+}^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})\right)]+{\displaystyle \frac{\xi \left(\mathbf{k}\right|\mathbf{r})}{2E{\left(\mathbf{k}\right|\mathbf{r})}^2}}{\displaystyle \frac{\partial f_F\left(E_{+}^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})\right)}{\partial E_{+}^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})}}-{\displaystyle \frac{\mathbf{k}·\mathbf{A}\left(\mathbf{r}\right)}{\mathbf{A}{\left(\mathbf{r}\right)}^2}}{\displaystyle \frac{1}{E\left(\mathbf{k}\right|\mathbf{r})}}{\displaystyle \frac{\partial f_F\left(E_{+}^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})\right)}{\partial E_{+}^{\mathbf{A}}\left(\mathbf{k}\right|\mathbf{r})}}\right\}.$$

(9)

When considering a superfluid flow, the vector potential $\mathbf{A}\left(\mathbf{r}\right)$ entering Equation (7) formally plays the role of an “effective” vector potential, and its physical meaning differs from that in classical electrodynamics. For this reason, we identify $\mathbf{A}\left(\mathbf{r}\right)$ with a constant wavevector $-\mathbf{q}$. As a consequence, in the presence of an inhomogeneity, such as, for example, a barrier, the phase of the order parameter solution of the LPDA equation can be written as $\phi \left(\mathbf{r}\right)=2\mathbf{q}·\mathbf{r}+2\varphi \left(\mathbf{r}\right)$, where the term $2\varphi \left(\mathbf{r}\right)$ develops because of the inhomogeneity in the system. Within the LPDA approach, the expression for the local current density is

$$\mathbf{j}\left(\mathbf{r}\right)={\displaystyle \frac{1}{m}}\left(\mathbf{q}+\nabla \varphi \left(\mathbf{r}\right)\right)n\left(\mathbf{r}\right)+2\sum _k{\displaystyle \frac{\mathbf{k}}{m}}{\mathcal{G}}_{11}^{\mathrm{mf}}\left(k;\mathbf{q}\mid \mathbf{r}\right),$$

(10)

where $n\left(\mathbf{r}\right)$ is the local density (5), and ${\mathcal{G}}_{11}^{\mathrm{mf}}\left(k;\mathbf{q}\mid \mathbf{r}\right)$ is the single-particle finite-temperature Green’s function in the mean-field approximation (see ref. [16] for further details).

In ref. [30], we use the LPDA approach to investigate the effects of a potential barrier embedded by an otherwise homogeneous superfluid system. In particular, we develop a numerical procedure to efficiently solve the real part of Equation (7), superimposing at the same time the current conservation, i.e., the continuity equation, using the expression (10) (see the Appendix of ref. [30]). In that work, we determine the Josephson critical current $J_c$ for barriers of different geometries, for various couplings across the BCS-BEC crossover, and at different temperatures. The results obtained in that work, the Josephson characteristics, the order parameter profiles, and the temperature dependence of the Josephson critical current, are in good qualitative agreement with the results of previous theoretical [29,32,33,34,35] and experimental [28,36,37,38] works on Josephson junctions. However, those results are expected not to be quantitatively reliable in the region of the coupling-temperature phase diagram, where the mean-field decoupling, on which BdG (and as a consequence, the LPDA approach) relies, is not valid. For this reason, in the subsequent work [16], we develop the mLPDA approach.

2.3. The mLPDA Approach

We develop the mLPDA approach, including the effect of pairing fluctuations at the level of the non-self-consistent t-matrix approximation in a minimal albeit essential way directly on top of the LPDA approach. Accordingly, we retain the self-consistent equation for the gap parameter (7) and replace the expressions for the local density (5) and for the local current density (10) with the following ones

$$n\left(\mathbf{r}\right)=2k_{B}T\sum _n{e}^{i{\omega }_n\eta }\int {\displaystyle \frac{d\mathbf{k}}{{\left(2\pi \right)}^3}}{\mathcal{G}}_{11}^{\mathrm{pf}}\left(k,{\omega }_n;\mathbf{q}\mid \mathbf{r}\right),$$

(11)

$$\mathbf{j}\left(\mathbf{r}\right)={\displaystyle \frac{1}{m}}\left(\mathbf{q}+\nabla \varphi \left(\mathbf{r}\right)\right)n\left(\mathbf{r}\right)+2k_{B}T\sum _n{e}^{i{\omega }_n\eta }\int {\displaystyle \frac{d\mathbf{k}}{{\left(2\pi \right)}^3}}{\displaystyle \frac{\mathbf{k}}{m}}{\mathcal{G}}_{11}^{\mathrm{pf}}\left(k,{\omega }_n;\mathbf{q}\mid \mathbf{r}\right),$$

(12)

where ${\omega }_n=(2n+1)\pi k_{B}T$ (n integer) is a fermionic Matsubara frequency and $\eta$ a positive infinitesimal. The Green functions ${\mathcal{G}}_{11}^{\mathrm{pf}}$ entering these expressions are evaluated by means of the replacement $\mu \left(\mathbf{r}\right)=\mu -V\left(\mathbf{r}\right)$ in the spirit of a local density approach. Their explicit form can be found in ref. [16], where they are obtained by generalizing the expression of the t-matrix propagator in the presence of a stationary flow of momentum $\mathbf{q}$.

This approach has been used to investigate the role of pairing fluctuations in the Josephson effect and has been validated against the experimental data for the Josephson critical current of refs. [20,21]. The main results of this investigation are briefly presented in the next section.

3. Results

In the first part of this section, we briefly present the results we obtained when investigating the superfluid flow with associated momentum $\mathbf{q}=q\widehat{\mathbf{x}}$ in a SsS slab geometry where the superfluid is homogeneous along the y and z directions orthogonal to the superfluid flow and extends to infinity on both sides of the Gaussian barrier

$$V\left(\mathbf{r}\right)=V_{0}exp\left(-{\displaystyle \frac{x^2}{2{\sigma }^2}}\right)$$

(13)

with $V_0/E_F=0.1$ ($E_F$ being the Fermi energy) and $k_F\sigma =2.5$.

The second part of this section is instead devoted to the investigation of the experimental systems realized in refs. [20,21], where two superfluid reservoirs of ultracold trapped Li atoms are coupled through a repulsive barrier of the form

$$V\left(\mathbf{r}\right)=V_0\left(z\right)exp\left(-2{\displaystyle \frac{x^2}{w{\left(z\right)}^2}}\right),$$

(14)

with

$$V_0\left(z\right)=V_0/\sqrt{1+{\left({\displaystyle \frac{z}{z_R}}\right)}^2}\mathrm{and}w\left(z\right)=w\sqrt{1+{\left({\displaystyle \frac{z}{z_R}}\right)}^2},$$

where $z_R=\pi w^2/\lambda$ is the Rayleigh range of the laser ($\lambda$ being its wavelength) used to generate the repulsive barrier.

3.1. Case of Study: Homogeneous Superfluid

When investigating the Josephson effect, we are mainly interested in determining the Josephson critical current $J_c$. In order to do that, we solve the mLPDA equation for different values of q, and for each of them, we obtain the profile of the order parameter (both its phase and magnitude) which preserves the current density. In this way, evaluating the asymptotic phase difference $\delta \varphi$ across the barrier (13), we are able to draw the Josephson characteristic $J\left(\delta \varphi \right)$ [29,39] and identify its maximum with the critical current $J_c$. Typical outcomes of this procedure are shown in Figure 1, where we report, for selected points of the Josephson characteristic, the corresponding order parameter and density profiles (matching colors). These profiles refer to the region $x>0$ since, on physical grounds, we expect the profile of the magnitude of the order parameter to be symmetric, and that of the phase to be antisymmetric.

We apply this procedure to superfluids at different temperatures and for coupling spanning the BCS-BEC crossover. In this way, we investigate the dependence of the critical current $J_c$ on coupling and temperature separately but also determine the role of pairing fluctuations by comparing the Josephson critical currents obtained solving the LPDA and mLPDA equations when considering the same physical system. In Figure 2, we show this comparison both at low temperature for coupling spanning the BCS-BEC crossover (panel (a)) and at representative couplings for T ranging from 0 to $T_c$ (panels (b–d)). Looking at panel (a), we note that LPDA and mLPDA return at low temperature the same qualitative dependence on coupling for the Josephson critical current; however, the inclusion of pairing fluctuations lowers the values of $J_c$. This feature reverses when increasing the temperature as it can be seen in panels (b–d) of Figure 2. Moreover, we note that, while the LPDA approach predicts the same qualitative dependence of the Josephson critical current on temperature, the mLPDA returns a curve which passes from being convex to concave going from the BCS to the BEC limit of the crossover, across an essentially linear behavior at unitarity. This temperature dependence is shared by the critical velocity of superfluid ⁴He (see Figure 12 of ref. [40]), extending de facto the parallelism noted over the years between these two systems [41,42].

3.2. Case of Study: Trapped Ultracold Atoms

A main result of the present work, with the ensuing validation of the mLPDA approach, is the quite good agreement obtained with the recent experimental measurements of the critical current both as a function of coupling [20] and temperature [21]. In both the experiments, the physical system consists of N Li atoms confined by means of a cigar-shaped harmonic potential. The trapping frequency ${\omega }_x$ is considerably smaller than ${\omega }_y$ and ${\omega }_z$ so that the atomic cloud acquires an ellipsoidal shape. The latter is further affected by the superimposition of two hard-wall potentials (homogeneous along the y and z directions) at $\pm x_w$, which ensure the stability of the system and reduce the atoms from N to $N_w$. The potential barrier (14) is moved with different velocities along the x axis, starting from the center of the atomic cloud, in order to infer its critical current.

When evaluating the critical current to compare with the experimental data, we adopt some approximations. We treat the trapping potential within the local density approximation, neglect ${\omega }_x$ since the barrier spans only the central region of the atomic cloud, and assume the current density lines to be straight. As a consequence, upon dividing the atomic cloud in tubular filaments, we evaluate the critical current $I_c$ as the maximum of the Josephson characteristic $I\left(\delta \varphi \right)$ of the whole system, obtained by integrating the local characteristics $J\left(\delta \varphi \right)$ of the filaments.

In Figure 3, we show the mLPDA results and experimental data for the Josephson critical current. Panels (a–c) report $I_c$ as a function of coupling for three different heights of the barrier (14), while panel (d) shows the temperature dependence of the critical current at unitarity. The shaded green area in panels (a–c) correspond to the experimental uncertainty on the physical quantities involved in the simulation. The orange squares, representing the results obtained solving LPDA equations, are also reported for comparison. The first three panels of Figure 3 show that the qualitative behavior of the experimental data for $I_c$ at low temperature is well reproduced even at the LPDA level, but the inclusion of pairing fluctuations is fundamental to obtain a quantitative agreement. Furthermore, the good agreement between the mLPDA predictions and experimental data for the Josephson critical current at unitarity (see panel (d)) validates the mLPDA approach even at higher temperatures.

4. Discussion

In this work, we reviewed three different approaches which can be employed when investigating inhomogeneous superfluid systems. The NLPDA approach is a convenient choice when the boundary conditions due to the presence of interfaces in the system are not known. On the other hand, at the cost of a further approximation, the LPDA is an easier equation to handle and returns reliable results at the mean-field level as long as the size of the inhomogeneity is larger than ${\xi }_{pair}$. To conclude, the mLPDA approach poses itself as a good candidate when one plans to investigate systems with non-trivial geometrical constraints and to consider at the same time the effect of pairing fluctuations.

In the near future, we envision investigating the Josephson effect in an SNS slab geometry, allowing for the coupling to change with the spatial variable. In this regard, the inclusion of pairing fluctuations could play a key role in reproducing the features of the temperature dependence of the experimental data for the critical current (see Figure 1 of ref. [28]). Indeed, the identification between the Josephson critical current with the critical intrinsic current at the center of the Gaussian barrier (13) in ref. [44] suggests that the different qualitative temperature dependencies of the critical current along the BCS-BEC crossover (see Figure 2), obtained within the mLPDA approach, could be the key aspect to reproduce the change in slope in Figure 1 of ref. [28].