Critical Josephson current in BCS-BEC crossover superfluids

TL;DR

This paper presents a microscopic model for Josephson currents in ultracold fermionic superfluids across the BCS-BEC crossover, revealing a maximum near unitarity and matching experimental and numerical results.

Contribution

It introduces an analytic model for the critical Josephson current that accounts for interaction strength and barrier height, extending previous bosonic models to fermionic superfluids.

Findings

Critical current peaks near unitarity.

Model agrees with experimental data.

Analytic expressions valid for various coupling regimes.

Abstract

We develop a microscopic model to describe the Josephson dynamics between two superfluid reservoirs of ultracold fermionic atoms which accounts for the dependence of the critical current on both the barrier height and the interaction strength along the crossover from BCS to BEC. Building on a previous study [F. Meier & W. Zwerger, Phys. Rev. A, 64 033610 (2001)] of weakly-interacting bosons, we derive analytic results for the Josephson critical current at zero temperature for homogeneous and trapped systems at arbitrary coupling. The critical current exhibits a maximum near the unitarity limit which arises from the competition between the increasing condensate fraction and a decrease of the chemical potential along the evolution from the BCS to the BEC limit. Our results agree quantitatively with numerical simulations and recent experimental data.

Equations

$$\hslash j_c=\frac{\hslash I_c}{A}=2B{t}_{cc}({\mu }_B)n_c$$

$$t_{cc}({\mu }_B)=\frac{|t|({\mu }_B)}{4k({\mu }_B)B}{\mu }_B.$$

$$\hslash j_c=\frac{{\mu }_B{n}_c}{2k({\mu }_B)}|t|({\mu }_B),$$

$$\hslash j_c^{AB}=\frac{\pi }{2}{\Delta }_0\frac{k_F^2}{16{\pi }^2}|t{|}^2({\mu }_F)$$

$${\lambda }_0^{\mathrm{BCS}}=\frac{3\pi }{8}\frac{{\Delta }_0}{{ϵ}_F}$$

$$\hslash j_c^{AB}=\frac{{\mu }_F{n}_c}{2k_F}|t{|}^2({\mu }_F)$$

$$\frac{j_c}{|t|({\mu }_B)}=\frac{j_c^{AB}}{|t{|}^2({\mu }_F)}=\frac{n{v}_F}{8}{\lambda }_0\sqrt{\tilde{\mu }}.$$

$$|j_1|=j_c\frac{|t|({\mu }_B)}{4},$$

$$\hslash j_c(x,y)={\int }_{-R_z}^{R_z}dz{n}_c(\text{r}){\mu }_B(\text{r})\frac{|t|({\mu }_B(\text{r}))}{4k({\mu }_B)(\text{r})R_z}$$

$$\hslash I_c=\int d^{3}r{n}_c(\text{r}){\mu }_B(\text{r})\frac{|t|({\mu }_B(\text{r}))}{4k({\mu }_B)(\text{r})R_z},$$

$$V(z)=\frac{V_0}{{\cosh }^2(z/d)}$$

$$I_c=-\frac{4}{{\hslash }^2}\sum _{kq}|t_{kq}{|}^2{\int }_0^{\infty }\frac{d{\omega }_1}{2\pi }{\int }_0^{\infty }\frac{d{\omega }_2}{2\pi }\frac{B(k,{\omega }_1)B(q,-{\omega }_2)}{{\omega }_1+{\omega }_2}$$

Full text

I INTRODUCTION

Josephson weak links [1, 2] represent a paradigmatic tool to investigate phase coherence between coupled superfluids [3, 4], and they have been realized in a variety of different setups: From solid state systems [5], where the Josephson effect has played a crucial role, e.g. to unveil the $d$-wave nature of pairing in high temperature superconductors [6, 7], to neutral superfluids of liquid Helium [8, 9] and weakly interacting Bose condensates of ultracold atoms [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. More recently, both theoretical [22, 23, 24, 25, 26, 27] and experimental [28, 29, 30, 31] efforts have been targeted to investigate the Josephson dynamics in the context of ultracold superfluids of strongly interacting fermion pairs. The so-called BCS-BEC crossover [32, 33], which relies on exploiting Feshbach resonances between two different hyperfine states of ultracold fermions [34], offers the unique possibility to investigate, within a single physical system, Josephson tunneling from the standard BCS limit of weakly bound fermion pairs to a Bose Einstein condensate of tightly bound molecules, including the scale invariant unitary Fermi gas at infinite scattering length as an intermediate phase.

In contrast with their bosonic counterparts, crossover superfluids are a challenging system to explore. On the experimental side, the strong interactions require, for coherent Josephson dynamics to be unveiled, the imprinting of thin barrier potentials on the sample, with thicknesses on the order of very few interparticle spacings [30, 31]. On the theoretical side, the crossover regime poses a difficult many-body problem for which many observables are accessible only by approximate methods [33]. In particular, microscopic calculations of the Josephson currents based on the solution of Bogoliubov de Gennes equations (BdGE) have been reported in Refs. [22, 26, 27]. Although this approach provides a qualitatively correct description of the entire crossover region, it is quantitatively reliable only in the BCS limit. An alternative strategy relies on the zero-temperature extended Thomas-Fermi models (ETFM) [23, 24, 25, 27]. These methods are based on a generalized Gross-Pitaevskii equation for the complex order parameter $\psi(\textbf{r})$ of bosonic pairs with mass $M=2m$ which tunnel across a barrier of effective height $V_{0}=2V_{0f}$. Here $m$ and $V_{0f}$ denote the mass and barrier height of a single fermion. The pair density is linked to the total fermion density $n=n_{\uparrow}+n_{\downarrow}$ by $|\psi(\textbf{r})|^{2}=n(\textbf{r})/2$ (for more details, see e.g. Ref. [27]). A key ingredient of the EFTM non-linear equation is the local chemical potential of pairs $\mu_{B}$, which is related to the chemical potential $\mu(n,a)$ of a uniform Fermi gas with total density $n$ and interspecies scattering length $a$ by $\mu_{B}(n,a)=2\mu(n,a)+|\epsilon_{b}|\equiv 2\mu_{F}(n,a)$. Here, $\epsilon_{b}$ is the binding energy of one dimer in the vacuum ($\epsilon_{b}=0$ for $a<0$). The advantage of this framework, compared with the mean field BdGE approach, is that $\mu(n,a)$ is an independent, purely thermodynamic input parameter for the calculation of the Josephson current. Hence, the correct chemical potential throughout the crossover can be accounted for by including the corresponding results obtained from independent Quantum Monte-Carlo calculations, see e.g. Ref. [35]. The EFTM method then properly describes the Josephson effect in the BEC limit. It fails even qualitatively, however, when extended to the unitary and BCS regimes [27, 30]. In particular, for fixed barrier heights, ETFM calculations yield a monotonically increasing Josephson current when moving from the BEC to the BCS limit. By contrast, the critical current is found to reach a maximum value near unitarity both experimentally [30, 31] and also in the numerical BdGE calculations [22, 26, 27]. Therefore, developing a unified model which provides a quantitatively correct description of the Josephson effect throughout the BCS-BEC crossover represents a relevant and timely problem, which is addressed in the present study.

Our work builds upon a previous microscopic description of the Josephson dynamics between two weakly interacting Bose Einstein condensates by one of us [12]. In this framework, the Josephson coupling between two superfluids connected by a high potential barrier is evaluated perturbatively in terms of a tunneling Hamiltonian for bosonic particles. To leading order, this approach gives rise to a critical current which is linearly proportional to the tunneling amplitude of condensed bosons, with a prefactor that characterizes the thermodynamics of the superfluid reservoirs. For lower barriers, the phase-dependent term in their coupling energy includes also a second order correction that involves the excitation of Bogoliubov phonon modes above the condensate, and which increases the maximum Josephson current across the barrier. Our aim in the following is to show that a straightforward extension of this model to neutral fermionic superfluids, which exhibit gapless Bogoliubov modes all along the BCS-BEC crossover, quantitatively reproduces both the results of numerical simulations and also the experimental data for the critical current over a wide range of parameters. Surprisingly, this holds even beyond the weak tunneling limit considered in the underlying microscopic model. We emphasize that by construction our framework does not account for the effects of real pair-breaking processes, which are irrelevant in the deep tunneling limit since the emerging Josephson coupling energy is then always much smaller than the gap for creating a single fermion excitation. Moreover, we stress that the present work does not intend to provide a description of the full current-voltage characteristic for a finite chemical potential difference $\Delta\mu\neq 0$ between the reservoirs. Instead, our sole aim here is to determine the critical current which sets the width of the ”zero-voltage” branch for arbitrary coupling strength and zero temperature.

Our results suggest that, throughout the BCS-BEC crossover, the critical Josephson current across a tunnel junction can be evaluated analytically solely on the basis of (i) the single-particle transmission amplitude for a pair at the relevant chemical potential, and (ii) bulk properties of the superfluid at equilibrium, specifically, the condensate density and the pair chemical potential. Most importantly, our approach provides a straightforward explanation for the experimentally observed maximum of the critical current density near unitarity which arises from two competing effects: the increase of the condensate fraction and the concurrent decrease of the chemical potential as a function of the dimensionless coupling strength $1/k_{F}a$. The identification of the unitary gas as the most stable fermionic superfluid [36, 37, 22, 38, 39] thus appears from a viewpoint which applies uniformly over the whole regime from weak to strong coupling, with no need to invoke different dissipation mechanisms, such as pair breaking on the BCS and phonon-like excitations on the BEC side [36, 22, 27].

II JOSEPHSON CURRENT FROM BCS TO BEC

II.1 First order Josephson current: BEC versus BCS results

In the following, we will argue that the results obtained in Ref. [12] for the Josephson current between two weakly interacting BECs can be extended in a straightforward manner to fermionic superfluids. This will provide an analytical result which holds throughout the full crossover from BCS to BEC. We start by considering a generic setup of a Josephson junction with the geometry depicted in Fig. 1a, where a rectangular barrier connects two homogeneous reservoirs. The results obtained for this idealized case will then be generalized to the harmonically trapped situation and Gaussian barriers sketched in Fig. 1b, relevant to compare our model predictions with experimental data [30, 31].

For barrier heights $V_{0}$ which considerably exceed the boson chemical potential $\mu_{B}$, it has been shown in Ref. [12] that the Josephson currents can be evaluated analytically in a systematic expansion in powers of an effective transfer Hamiltonian $\hat{H}_{T}$. In leading order, this gives rise to a non-dissipative Josephson current between the two superfluid reservoirs with a sinusoidal current-phase relation $I(\varphi)=I_{c}\sin{\varphi}$, where $\varphi$ is the difference between the condensate phases across the junction [40]. Assuming a homogeneous situation with transverse area $A=L^{2}$, the associated critical current density

[TABLE]

can be expressed as the product of the transfer matrix element $t_{cc}(\mu_{B})$, associated with coherent tunneling of bosons at an energy equal to their equilibrium chemical potential $\mu_{B}$, times the associated condensate density $n_{c}$. Formally, Eq. (1) also contains the total longitudinal size $2B$ of the system which however drops out, as expected. Indeed, for the junction geometry of Fig. 1a, the matrix element $t_{cc}$ has been explicitly evaluated in Ref. [12], and we conveniently recast it here as

[TABLE]

Here $k(\mu_{B})=\sqrt{2M\mu_{B}}/\hbar$ is the wave vector of a boson with mass $M$ and chemical potential $\mu_{B}$, while $|t|(\mu_{B})$ is the associated single-particle transmission amplitude. The critical current density can be therefore expressed in the simple form

[TABLE]

which shows that inter-particle interactions affect $j_{c}$ only through bulk properties of the superfluid. Specifically, they involve the chemical potential $\mu_{B}$ and the condensate density $n_{c}$ which, to lowest order in the interactions between the pairs, coincides with the full Bose density. It is important to notice that the microscopic tunneling amplitude $|t|(\mu_{B})$ is the one for a single boson. It depends on the interaction strength only to the extent that this determines the relevant energy at which the tunneling process occurs. As discussed in Ref. [12], the separation between single particle and many-body properties within Eq. (3) originates from the underlying assumption that the barrier height $V_{0}$ greatly exceeds the chemical potential $\mu_{B}$. Under these conditions, our major claim in the following is that Eq. (3) correctly describes the Josephson current at any coupling strength.

The fact that Eq. (3) may hold even for fermionic superfluids throughout the crossover from the BCS to the BEC limit is supported by the quite remarkable observation that, upon a proper identification of the variables, this result essentially coincides with the well known expression for the critical current between two BCS superconductors. Within BCS theory, the critical current of a Josephson junction in the tunneling limit is known to be equal to the current flowing in the normal state at a finite applied voltage $eV=\pi\Delta_{0}/2$, where $\Delta_{0}$ is the energy gap at zero temperature [1, 41, 3]. By evaluating the normal state conductance per unit area, which scales linearly with the number of discrete transverse modes, the Ambegaokar-Baratoff (AB) formula for the pair current density between two BCS superfluids coupled by a thin tunneling barrier can be expressed in the form [26]

[TABLE]

where $|t|^{2}(\mu_{F})$ is the transmission probability of a single fermion at the Fermi energy $\mu_{F}\to\epsilon_{F}$. While at a first glance the expressions in Eqs. (3) and (4) appear totally disconnected, it turns out that they are essentially equivalent. Indeed, in the BCS limit, the condensate density $n_{c}=\lambda_{0}\cdot n/2$ is linked to the gap by the relation (see e.g. Refs. [42, 43])

[TABLE]

where $n=k_{F}^{3}/(3\pi^{2})$ is the total fermion density and $\lambda_{0}$ denotes the condensate fraction. Hence, Eq. (4) can be recast in the form

[TABLE]

where the similarity with Eq. (3) is now apparent. Indeed, noting that $2k_{F}\sim 2\sqrt{2m\mu_{F}}/\hbar\to k(\mu_{B})$ in terms of the chemical potential $\mu_{B}=2\mu_{F}$ of pairs, the two expressions are perfectly equivalent if the transmission probability $|t|^{2}(\mu_{F})$ of a single fermion is replaced by the transmission amplitude $|t|(\mu_{B})$ of a pair.

More specifically, it is instructive to consider the thermodynamic prefactor, common to the BCS result Eq. (6) and the BEC one Eq. (3) for the critical current density, for a fermionic superfluid at arbitrary coupling. This quantity, which encodes the many-body properties of the system, can be suitably expressed in terms of the condensate fraction $\lambda_{0}=2n_{c}/n$ and the dimensionless chemical potential $\tilde{\mu}=\mu_{F}/\epsilon_{F}$ of one fermion, which enters into the BEC result Eq. (3) through the standard mapping $\mu_{B}(n,a)=2\mu(n,a)+|\epsilon_{b}|\equiv 2\mu_{F}(n,a)$ between the Bose and Fermi chemical potential discussed in the previous section. With $v_{F}=\sqrt{2\epsilon_{F}/m}$ as the Fermi velocity of the non-interacting Fermi gas one obtains

[TABLE]

Eq. (7) shows that the critical current density, in units of the natural scale set by $(n\,v_{F})$, is determined by the product of the condensate fraction $\lambda_{0}$ with the square root of the (normalized) fermion chemical potential $\sqrt{\tilde{\mu}}$.

The dependence of these two thermodynamic quantities on the coupling strength is shown in Fig. 2, which presents the results for $\lambda_{0}$ (green line) and $\tilde{\mu}$ (blue line) obtained within a Luttinger-Ward approach to the BCS-BEC crossover problem developed by Haussmann et al. [44]. Apparently, these two bulk properties exhibit an opposite behavior as a function of the coupling strength: while $\tilde{\mu}$ monotonically decreases from one towards zero upon moving from the BCS to the BEC limit, the condensate fraction $\lambda_{0}$ monotonically increases from exponentially small values in the BCS regime towards unity for a weakly interacting Bose superfluid. A qualitatively similar behavior is also obtained within a mean field BdGE approach, see e.g. Ref. [45] and references therein. However, the results of the two theories differ substantially on a quantitative level. In particular, for the gas at unitarity, mean field calculations yield a condensate fraction $\lambda_{0}\simeq 0.7$ instead of $0.51$, and a chemical potential $\tilde{\mu}\simeq 0.56$ instead of the value $0.36$ obtained within the Luttinger-Ward approach, which agrees quite well with the experimental value $\tilde{\mu}=\xi_{s}=0.37$ for the associated Bertsch parameter, see Ref. [46].

As a consequence of the competing trends of $\lambda_{0}$ and $\tilde{\mu}$, the thermodynamic pre-factor in Eq. (7) reaches a maximum near unitarity (see red line in Fig. 2), in qualitative agreement with both experimental [30] and numerical [22, 27] findings. Note, however, that the origin of the non-monotonic behavior in the present theory is quite different from the interpretation given by Spuntarelli et al. [22] for the case of weak barriers. There, following a Landau criterion, the rise and fall of the maximum Josephson current was ascribed to the existence of two distinct types of critical velocities associated with two different excitation branches: pair-breaking on the BCS and phonons on the BEC side, respectively. By contrast, our result Eq. (7) shows that in the tunneling regime this peculiar non-monotonicity emerges from the competition between $\lambda_{0}$ and $\tilde{\mu}$ uniformly throughout the crossover region. This conclusion is further supported by recent experimental studies [30, 31], which unambiguously identified vortex rings and phonons, rather than ”fermionic” pair-breaking excitations, as the microscopic mechanisms responsible for the breakdown of the dissipationless flow, even on the BCS side of the crossover.

In the following, we focus on the difference between the BEC result in Eq. (3) and the BCS one of Eq. (6) that, as highlighted in Eq. (7), amounts to the ratio between the tunneling amplitude of one pair and the transmission probability of one fermion, evaluated at their relevant energies. Naively, such a quantity could be expected to exhibit an exponential dependence upon the barrier thickness $d$ and the normalized height $V_{0}/\mu_{B}$, hence causing a large quantitative mismatch between the BCS and BEC results. However, employing the mapping of the ”fermionic” BdGE approach onto the ”bosonic” ETFM one [27], we have verified that this is not the case: The ratio of these two transmission terms is of order unity for a wide range of barrier geometries, featuring only a very weak dependence upon the barrier properties and the chemical potential. Indeed, applying Eq. (3) to the BCS limit of the crossover yields a critical current density that matches the BCS result Eq. (6) within a few $10\%$ accuracy. For instance, over an extremely wide range $k_{F}d\in[10^{-3},10^{2}]$ of dimensionless barrier thicknesses, the two terms match one another within $\pm 50\%$, both for rectangular and Eckart barriers of relative heights $V_{0}/\mu_{B}\in[1,10^{4}]$. Within the experimentally relevant regime of barrier widths $1\lesssim k_{F}d\lesssim 4$ and moderate heights, $V_{0}/\mu_{B}\lesssim 2$ [30, 31], the deviations are in fact smaller than $20\%$. In view of these results, in the following we will thus use the extension of the bosonic expression Eq. (3) to arbitrary couplings to describe the critical current of a Josephson junction along the entire BCS-BEC crossover.

II.2 Second order contributions and extension to inhomogeneous samples

In order to benchmark our analytic theory against numerical and experimental results, it is necessary to also take into account higher harmonics that contribute to the critical current. In fact, as already mentioned in the previous section, the tunneling limit of exponentially small values of $|t|(\mu_{B})$ is hard to achieve in practice and realistic barrier heights $V_{0}$ are close to, or only slightly larger than, the chemical potential [30, 31]. Therefore, higher order terms in the expansion in powers of the transfer Hamiltonian are non-negligible. In particular, in second order, these lead to an additional non-dissipative contribution $I_{1}\sin{2\varphi}$ to the coherent flow. It has half of the period $\varphi\to\varphi+2\pi$ of the fundamental Josephson current and, within our framework, it arises from condensate to non-condensate tunneling processes, involving the excitation of a Bogoliubov phonon. By suitably recasting the results of Ref. [12], the magnitude of the (negative) second harmonic current density can be linked to the dominant first order term via

[TABLE]

which still exhibits the separation between single-particle and many-body properties discussed above for the leading contribution given by Eq. (3). It is important to remark that Ref. [12], at second order in the transfer Hamiltonian, also predicts normal and phase-dependent dissipative currents. These terms are absent, however, at a vanishing value $\Delta\mu=0$ of the difference in chemical potential between the two superfluid reservoirs. As a result, these contributions are irrelevant for the determination of the maximum Josephson current, which is the aim of our present study.

Moreover, in order to test the predictions of our model against recent experimental data [31] and BEC results of EFTM numerical simulations [47] available in the literature, it is necessary to generalize the results obtained in Eq. (3) to the experimentally relevant case [30, 31] of inhomogeneous samples and Gaussian barriers, see Fig. 1b. In particular, the result in Eq. (3) for the (first order) critical current density can be easily extended to account for an inhomogeneous condensate density profile $n_{c}(\textbf{r})$ as

[TABLE]

Further integration of $j_{c}$ along the transverse directions $(x,y)$ then yields the total current

[TABLE]

which obviously recovers the result obtained in the homogeneous case [12], where $I_{c}$ scales linearly with the junction area $A$. In the generic inhomogeneous case, Eq. (10) may be regarded as a sum of transmitted particles over all infinitely small volumes centered at position $(r,z)$, each one containing a number $n_{c}(r,z)$ of condensed bosons, at energy $\mu_{B}(r,z)$.

By following a similar procedure and recalling Eq. (8), one can straightforwardly obtain an analytic expression also for the second order contribution to the Josephson current for inhomogeneous samples. While this latter correction is negligible in the deep tunneling limit $V_{0}/\mu_{B}\gg 1$, it may cause an up to 15% increase of the maximum supercurrent $I_{Max}$ in the experimentally relevant regime $V_{0}/\mu\sim 1$. This can be easily verified on the basis of the analytic results obtained by Goldobin et al. [48], that theoretically investigated how an arbitrary second harmonic term in the current-phase relation of a generic junction affects the maximum Josephson current.

III COMPARISON WITH BdGE NUMERICAL RESULTS FOR THE CROSSOVER REGION

In the following, we compare the predictions of our model for a fixed barrier geometry with numerical results based on the solution of time-dependent mean field BdGE, reported in Ref. [27]. To this end, we consider the specific tunnel junction setup of Zou and Dalfovo [27], analogous to the one sketched in Fig. 1a. In Fig. 3 the dashed blue line represents the critical fermion current ($2\,I_{Max}$) predicted by our analytic model for the barrier and reservoir parameters detailed in the figure legend, and employing the mean field continuum results for chemical potential and condensate density, respectively (see Ref. [45] and references therein).

As already emphasized in the previous section, we obtain a non-monotonic critical current as the interaction is changed from the BCS to the BEC regime. Compared with the behavior of $\lambda_{0}\,\sqrt{\tilde{\mu}}$ presented in Fig. 2, the shape of $I_{Max}$ as a function of $1/k_{F}a$ is quantitatively modified owing to the additional energy dependence of $|t|(\mu_{B})$, which monotonically decreases when moving from the BCS towards the BEC limit. In particular, this causes an overall suppression of the maximum Josephson current, more pronounced in the strong attraction limit of tightly bound pairs, where the chemical potential most strongly varies with the coupling strength.

Notably, our model appears in reasonable quantitative agreement with the results of BdGE simulations [27], see red dots in Fig. 3. In particular, the mismatch between our theory and the numerical results is largely ascribable to finite size effects connected with the small box employed in the simulation, strongly affecting the bulk superfluid properties, relative to the continuum ones. For instance, the agreement significantly improves when we employ the BdGE chemical potential values obtained within the actual tightly confining box [49], systematically higher than the continuum result, see blue solid line in Fig. 3. In this case, our theory matches the numerical data with deviations not exceeding 30$\%$ throughout the crossover region. Additionally, the residual mismatch can be largely attributed to effects not taken into account when evaluating our model predictions: First, the small box size may also cause an increased condensate fraction, relative to the continuum case, especially on the BCS side of the crossover. Second, the determination of the critical current by solving the time dependent BdGE was obtained by imparting a sizable initial excitation to the system [27], yielding an additional kinetic energy contribution to the overall particle energy, neglected in our model.

The ability of our analytic theory to yield quantitatively accurate Josephson currents in the crossover region is further demonstrated by the comparison presented in Fig. 4. There, we contrast our prediction based on Eq. (3) for the fermion current density ($2\hbar j_{c}$) of unitary superfluids with the numerical results obtained by Zou and Dalfovo [27] solving stationary BdGE, both as a function of the barrier height for a few values of the dimensionless barrier thickness (Fig. 4a), and vice-versa for different barrier heights as a function of the thickness (Fig. 4b). Following the convention of Ref. [27], the barrier height is given in the legends for one fermion, $V_{0f}=V_{0}/2$. Apparently, our analytic predictions quantitatively reproduce the numerical results remarkably well over a wide range of barrier geometries. Indeed, as long as the barrier height exceeds the superfluid chemical potential, sizable deviations are observed only for the case of anomalously thin barriers (see black squares in Fig. 4a), for which our theoretical model based on Eq. (3) would converge to the numerical data only for exceedingly high $V_{0}$ values. On the other hand and not surprisingly, our model fails even at the qualitative level for $V_{0}/\mu<1$, i.e. out of the tunneling regime, see e.g. blue line and triangles in Fig. 4b (BdGE calculations in the present box geometry yield at unitarity $\mu_{F}/\epsilon_{F}\simeq 0.7$ [49]).

IV COMPARISON WITH EXPERIMENTS

Finally, we compare our model predictions with experimental data that were obtained in the Lithium lab at LENS [31], and also with ETFM numerical simulations carried out in the BEC limit [47]. To this end, we evaluate the maximum Josephson current $I_{Max}$ accounting for both first and second order contributions [48], in the case of harmonically trapped samples and Gaussian barriers (see Fig. 1b), exploiting Eqs. (10) and (8), respectively. For simplicity, the Gaussian barrier is approximated with an Eckart potential of the form

[TABLE]

for which the single-particle transmission probability can be evaluated analytically, see e.g. Ref. [50]. In particular, it can be verified that a Gaussian profile with $e^{-2}$ waist $w$ can be excellently approximated by an Eckart function with $d\sim 0.6\,w$. Inserting the transmission amplitude of the Eckart barrier both in Eq. (10) for the first order current and in the analogous expression for the second harmonic based on Eq. (8), and accounting for the relevant superfluid atom number and trap frequencies employed by Burchianti et al. [31], we obtain our zero-temperature model predictions for the resulting maximum current $I_{Max}$ [48]. For this comparison, we exploited the chemical potential and condensate fraction obtained within the Luttinger-Ward approach [44]. An important point to emphasize in this context is that the approach provides a quantitatively precise description of the thermodynamic properties not only near unitarity, as shown in Fig. 2, but also in the BCS and the deep BEC limit. Indeed in the regime $1/k_{F}a>2$, the normalized chemical potential $\tilde{\mu}\to k_{F}a_{dd}/(3\pi)$ approaches the result associated with a dilute gas of Bosons with a dimer-dimer scattering length $a_{dd}=0.59\,a$ which is very close to the exact value $a_{dd}=0.6\,a$ [51].

In Fig. 5 we compare our theory predictions (solid lines) for the maximum currents $I_{Max}/(\frac{N}{2}\omega_{z})$, plotted as a function of the barrier height $V_{0}/\mu_{0}$, with experimental data [31] for BEC (green circles), crossover (red squares) and BCS (blue diamonds) superfluids, respectively. In the BEC limit, we also benchmark our model against results of numerical simulations [47] (empty green triangles). Here $\mu_{0}$ denotes the pair chemical potential at the center of the cloud evaluated for the superfluid gas in the absence of the repulsive barrier. The different regimes are marked in the legend in terms of the peak interaction parameter $1/k_{F}a$. We emphasize that, given the relatively small values of $V_{0}/\mu_{0}$ accessible in the experiments, second order effects are explicitly included in the evaluation of $I_{Max}$ following Ref. [48]. For each interaction regime, dashed lines represent our model predictions once a $\pm 5\%$ variation of the peak chemical potential is assumed when evaluating the transmission amplitude.

The comparison shown in Fig. 5 clearly demonstrates that our model is able to reproduce the results obtained in the literature so far with remarkable accuracy for all interactions and barrier heights. In particular, we emphasize that our analytic theory quantitatively matches the ETFM numerical results (green triangles) [47] obtained in the BEC limit even for barrier heights as low as $V_{0}/\mu_{0}\sim 0.7$, i.e. even when a sizable part of the trapped sample is out of the tunneling regime. There, inclusion of second order contributions to $I_{Max}$ is essential, as it yields corrections to the first order term up to 15%. Regarding the crossover and BCS regimes, our model correctly reproduces the observed trends, at least within the large uncertainties of the experimental data and of the zero-temperature approximation done in the theoretical analysis. In particular, we remark that inclusion of the non-perturbative results for both chemical potential and condensate fraction [52] is essential to obtain quantitative agreement with the data, which would not match our theory if mean field results were employed. As a consequence, although more experimental studies are required to fully validate our model within the crossover and BCS regimes, we stress that the measure of the critical Josephson current may provide a valuable experimental tool to determine the condensate density in crossover superfluids, a quantity which is hard to access by other means.

V CONCLUSIONS AND OUTLOOK

Building on a microscopic model for the Josephson dynamics between two weakly interacting Bose Einstein condensates [12], we have developed a simple theoretical framework which captures the corresponding critical current of strongly interacting fermionic superfluids throughout the BCS-BEC crossover. As testified by the comparison with available numerical [27, 47] and experimental [31] data, our model provides a quantitative description of the critical currents in generic junction geometries, solely based on bulk properties of the superfluid state and the knowledge of the single boson transmission amplitude. An important feature of our theory is that it provides a consistent explanation for the non-monotonic trend for the maximum Josephson current observed near unitarity. Its origin is connected with the interplay between a decrease of the chemical potential and a concomitant increase of the condensate fraction if the interactions are tuned from the BCS to the BEC limit. By contrast, EFTM models, to which our theory naturally connects deep in the BEC limit, fail when applied to intermediate and weak couplings even at the qualitative level [22, 27]. Indeed, in these approaches the condensate density by construction coincides with the superfluid one at all couplings. In the future, therefore, it will be interesting to test whether a modified EFTM formalism, which explicitly accounts for the reduction of the condensate fraction upon going from the BEC towards the unitary gas and eventually the BCS limit, yields results consistent with the ones obtained here or through BdGE approaches. We also anticipate that our model could be extended to describe dissipative normal currents [31], and to account for finite temperature effects. Moreover, starting from the present study, it might be possible to investigate how the critical current evolves from the tunneling regime considered here to the case of weak barriers studied in Refs. [22, 26, 37], in which pair-breaking excitations are expected to become relevant already at intermediate couplings close to unitarity [26].

On a rather fundamental level, it remains an open challenge to justify our use of the result in Eq. (3), which is based on viewing the Josephson effect throughout the BCS-BEC crossover in terms of the coherent tunneling of pairs, by developing a truly microscopic approach which involves the constituent fermions. Within such a description and in the tunneling limit, the Josephson effect is second order in the single fermion tunneling amplitudes $t_{kq}$. Quite generally, the exact critical current at zero temperature

[TABLE]

can be expressed in terms of the spectral functions $B(k,\omega)$ associated with the anomalous propagator of fermionic superfluids. In the BCS-limit, $B^{\rm BCS}(k,\omega)=2\pi\,u_{k}v_{k}\,[\delta(\omega-E_{k}/\hbar)-\delta(\omega+E_{k}/\hbar)]$ is determined by the gap parameter $\Delta_{k}=u_{k}v_{k}\cdot 2E_{k}$ and the excitation energies $E_{k}=\sqrt{\xi_{k}^{2}+|\Delta_{k}|^{2}}$ of the fermionic quasiparticles in the standard manner [53]. The expression Eq. (12) then reduces to the Ambegaokar-Baratoff result Eq. (4) which involves the tunneling probability $|t|^{2}(\mu_{F})$ evaluated right at the Fermi energy. In order to verify the non-monotonic dependence of the critical current on the dimensionless coupling constant $1/k_{F}a$ along the BCS-BEC crossover within such a fermionic approach one needs to determine the anomalous spectral functions $B(k,\omega)$ that enter Eq. (12) for arbitrary coupling, a problem which does not seem to have been discussed so far. Note Added: After submission of this work, our model has been successfully employed to analyze new experimental data by W. J. Kwon et al., arXiv:1908.09696, and also - in an extension to a two-dimensional setup - by N. Luick et al., arXiv:1908.09776.

Acknowledgements.

We thank F. Dalfovo and P. Zou for providing the chemical potential values relevant to analyze the data of Ref. [27], K. Xhani for sharing the ETFM results reported in Ref. [47] for trapped BECs, and B. Frank for the non-perturbative results shown in Figure 2 and employed to analyze the experimental data of Figure 4. We acknowledge fruitful discussions with F. Dalfovo, E. Neri, A. Recati, G. Roati, F. Scazza and K. Xhani. This work was supported by the ERC through grant no. 637738 PoLiChroM and in part also by the Nano-Initiative Munich (NIM).