Non-local equation for the superconducting gap parameter

# Non-local equation for the superconducting gap parameter

S. Simonucci and G. Calvanese Strinati Division of Physics, School of Science and Technology
Università di Camerino, 62032 Camerino (MC), Italy
and
INFN, Sezione di Perugia, 06123 Perugia (PG), Italy
###### Abstract

The properties are considered in detail of a non-local (integral) equation for the superconducting gap parameter, which is obtained by a coarse-graining procedure applied to the Bogoliubov-de Gennes (BdG) equations over the whole coupling-vs-temperature phase diagram associated with the superfluid phase. It is found that the limiting size of the coarse-graining procedure, which is dictated by the range of the kernel of this integral equation, corresponds to the size of the Cooper pairs over the whole coupling-vs-temperature phase diagram up to the critical temperature, even when Cooper pairs turn into composite bosons on the BEC side of the BCS-BEC crossover. A practical method is further implemented to solve numerically this integral equation in an efficient way, which is based on a novel algorithm for calculating the Fourier transforms. Application of this method to the case of an isolated vortex, throughout the BCS-BEC crossover and for all temperatures in the superfluid phase, helps clarifying the nature of the length scales associated with a single vortex and the kinds of details that are in practice disposed off by the coarse-graining procedure on the BdG equations.

###### pacs:
74.20.Fg,03.75.Ss,05.30.Jp,74.25.Uv

## I Introduction

Non-locality lies at the heart of the phenomenon of superconductivity and is related to the finite spatial size of Cooper pairs. The first recognition of this phenomenon came from the work by Pippard Pippard-1953 (); Pippard-1955 (), who realised that (in a clean system) the supercurrent at a spatial point is determined by a spatial average of the vector potential over a neighboring region with the size of a Cooper pair.

When solving for the gap equation in inhomogeneous situations to obtain the gap parameter , a second length scale (the healing length ) enters the problem to describe the typical distance over which tends to recover its bulk value in the presence of a localized perturbation. Both lengths and depend on the coupling value of the inter-particle interaction Pistolesi-1994 (); Pistolesi-1996 (), which gives rise to Cooper pairs and their condensation to begin with, as well as on temperature Palestini-2014 (). The relative behaviour of with respect to is expected to determine the relevance of the local vs non-local behaviour of with respect to its surrounding values.

Gor’kov Gorkov-1959 () first realized the importance of distinguishing between the two lengths and , in the process of deriving the Ginzburg-Landau (GL) equation from the coupled differential equations for the normal and anomalous single-particle temperature Green’s functions (or, alternatively, from the integral version of the Bogoliubov-de Gennes (BdG) equations deGennes-1966 ()). Although Gor’kov’s derivation applies to the weak-coupling limit whereby is much larger than the inter-particle distance (as specified by the inverse of the Fermi wave vector where is the particle density), the local (differential) GL equation could be retrieved close to the critical temperature where . In the opposite case of strong coupling, when Cooper pairs turn into composite bosons with size , the condition can be realized also at zero temperature. In this case, it is the local (differential) Gross-Pitaevskii (GP) equation for composite bosons that can be derived from the BdG equations, as shown in Ref. Pieri-2003 ().

The way and evolve with respect to each other at zero temperature characterises the BCS-BEC crossover, whereby in the BCS limit and in the BEC limit Pistolesi-1994 (); Pistolesi-1996 (). In practice, the crossover between these two limits can be obtained by varying the inter-particle attraction, in such a way that the system evolves from a BCS state where pairs of (opposite spin) fermions are described by Fermi statistics, to a BEC state where two-fermion dimers (or composite bosons) are described by Bose statistics. A substantial amount of theoretical work Eagles-1969 (); Leggett-1980 (); NSR-1985 (); Randeria-1989-90 () had preceded the explicit experimental realization of the BCS-BEC crossover with ultra-cold Fermi gases Jin-2003 (); Hulet-2003 (); Grimm-2003 (); Jin-2004 (). In these systems, the attraction between opposite-spin fermions can be taken to be of zero range in space and instantaneous in time, conditions that will be assumed to hold in the rest of this paper in line with the original Galitskii approach for a dilute Fermi system Galitskii-1958 (). In condensed matter, where the inter-particle interaction can have more complicated forms, on the other hand, consideration of the BCS-BEC crossover was originally suggested by the fact that in high- superconductors the product is of order unity, corresponding to the “unitary” regime which is intermediate between the BCS and BEC limits Pistolesi-1994 (). Growing evidence for the occurrence of this crossover has lately emerged also in two-band superconductors with iron-based materials Chubukov-2016 ().

In the context of the BdG and related equations, recently a method was devised to obtain a non-linear differential equation for the gap parameter by performing a suitable spatial coarse graining of the BdG equations, which deals with the smoothness of the spatial variations of the magnitude and phase of on a different footing for smoothing out short-range details of the gap parameter Simonucci-2014 (). This equation (referred to as a local phase density approximation (LPDA) to the BdG equations) was found to recover both the GL equation in weak coupling close and the GP equation in strong coupling at low temperature. In Ref. Simonucci-2014 () the LPDA equation was applied at any temperature in the superfluid phase throughout the BCS-BEC crossover, for the test case of an isolated vortex for which also an accurate numerical solution of the BdG equations is available to compare with Simonucci-2013 (). This test led to an extremely good agreement between the LPDA and BdG calculations essentially for all couplings and temperatures, with the exception of the BCS (weak-coupling) regime at low temperature where deviations between the two calculations have emerged. In Ref. Simonucci-2014 () the reason for this discrepancy was attributed to the fact that, in this regime, the vortex size (or healing length ) becomes comparable with the Cooper pair size Pistolesi-1994 (); Pistolesi-1996 (), thereby questioning the validity of a local (differential) approach like the LPDA equation in this restricted regime of coupling and temperature. Later, the LPDA approach was successfully applied to generate in a self-consistent way large arrays of vortices throughout the BCS-BEC crossover Simonucci-2015 (), which can be produced by setting an ultra-cold trapped Fermi gas into rotation. In this way, it was possible to account for the experimental data that had provided the first direct evidence of the superfluid phase in these systems Zwierlein-2005 ().

Although the end result of the coarse-grained derivation of Ref. Simonucci-2014 () was the local (differential) LPDA equation, a non-local (integral) version of the LPDA equation (that can accordingly be referred to as the NLPDA equation) was also reported in that reference, at an intermediate step between the original BdG equations and the LPDA equation. The (integral) NLPDA equation for the gap parameter contains a non-local kernel that depends on the superfluid gap itself in a highly non-linear way. However, the (non-local) NLPDA equation was not further examined in Ref. Simonucci-2014 (), where the attention was concentrated only on the (local) LPDA equation. Examination of the (non-local) NLPDA equation can also be of interest in itself, since it can give access to problems that are difficult to deal with using the (local) LPDA equation.

Purpose of this paper is to consider in detail the non-local (integral) NLPDA equation for the gap parameter and study the properties of its kernel, in order to determine its spatial range for all couplings and temperatures in the superfluid phase throughout the BCS-BEC crossover. This turns out be a non-trivial task, especially in the weak-coupling (BCS) regime at low temperature where the kernel shows rapid and slowly decaying spatial oscillations. This study eventually enables us to identify the spatial extent of the “granularity” associated with the coarse-graining procedure on which the NLPDA (and, as a consequence, the LPDA) equation rests, as well as to understand the reason for the failure of the LPDA equation in weak coupling at low temperature mentioned above. We will find that, for all couplings throughout the BCS-BEC crossover and temperatures in the superfluid phase, the spatial range of this kernel coincides with the Cooper pair size , a quantity which was independently determined in Ref. Palestini-2014 () by analyzing the pair correlation function for opposite-spin fermions. In this context, a method will also be implemented for solving numerically the (integral) NLPDA equation in an efficient way. This method will be explicitly utilized to study an isolated vortex and to compare the results with the solutions of both the LPDA Simonucci-2014 () and BdG Simonucci-2013 () calculations, for all couplings throughout the BCS-BEC crossover and temperatures in the superfluid phase. This test calculation will highlight the length scales associated with an isolated vortex by the three (BdG, LPDA, and NLPDA) calculations, a result that will be especially instructive in weak coupling at low temperature.

The paper is organized as follows. In Section II the properties of the kernel of the NLPDA equation are studied in detail, both in wave-vector and real space, and the spatial extent of this kernel is determined for all couplings throughout the BCS-BEC crossover and temperatures in the superfluid phase. Knowledge of these properties will also enable us to set up an efficient strategy for the numerical solution of the NLPDA equation. This will be explicitly done in Section III for the case of an isolated vortex embedded in an infinite superfluid, for all couplings throughout the BCS-BEC crossover and temperatures in the superfluid phase. The range validity of the (local) LPDA equation will also be explicitly determined from an analysis of the (non-local) NLPDA equation from which it is derived. Section IV gives our conclusions. In Appendix A a model kernel is introduced in wave-vector space, whose analytic solution in real space will help us to identify the origin of the spatial oscillations of the kernel of the NLPDA equation that show up in weak coupling at low temperature. In Appendix B a method is implemented for solving numerically the NLPDA equation, by devising a novel approach for calculating the Fourier transforms from real to wave-vector space and vice versa.

## Ii The kernel of the Non-Local NLPDA gap equation

The following non-local (integral) gap equation

 −m4πaFΔ(r)=∫dRΔ(R)∫dQπ3e2iQ⋅(r−R)KA(Q|r) (1)

was obtained in Ref. Simonucci-2014 () at an intermediate step in the process of deriving the local (differential) LPDA equation by a coarse-graining procedure of the BdG equations. The kernel of this equation reads:

 KA(Q|r)=∫dk(2π)3{1−2fF(EA+(k;Q|r))2EA(k;Q|r)−mk2} (2)

where is the fermion mass, the scattering length of the two-fermion problem, the Fermi function at temperature ( being the Boltzmann constant),

 EA±(k;Q|r) =  ⎷(k22m+Q22m−¯μ(r)−A(r)m⋅Q)2+|Δ(r)|2 (3) ± km⋅(Q−A(r)),

and . In the above expressions, is the vector potential, the local chemical potential in the presence of an external (trapping) potential , and the magnitude of the local gap parameter.

Apart from its role as an intermediate step in deriving the local (differential) LPDA equation, no further consideration was given in Ref. Simonucci-2014 () to the non-local (integral) NLPDA equation (1), since all attention was concentrated on the LPDA equation. Here, our interest is to focus directly on the NLPDA equation (1) and study in detail the properties of its kernel as well as of its Fourier transform

 KA(R|r)=∫dQπ3e2iQ⋅RKA(Q|r), (4)

aiming at determining its spatial range in real space. To this end, it will be sufficient to consider the pivotal case with , , and , where is the uniform mean-field value of the gap parameter for the homogeneous system Physics-Reports-2017 (). Accordingly, we set to simplify the notation.

This analysis will be carried out for given temperature in the superfluid phase and coupling throughout the BCS-BEC crossover, whereby the coupling parameter ranges from in the weak-coupling (BCS) regime when , to in the strong-coupling (BEC) regime when , across the unitary limit (UL) when diverges. In practice, the “crossover region” of most interest is approximately limited by the interval .

A. Properties of the kernel

We begin by considering a number of properties of the kernel in -space, which is given by the expression (2) with the further account of the above provisions.

Spherical symmetry. When , the kernel depends only on . This is because the transformation where is a three-dimensional rotation can be compensated by an analogous rotation of the integration variable in Eq. (2). We thus write in place of (and, correspondingly, in place of where ).

Non analyticity at . At zero temperature, the Fermi function in Eq. (2) is non-vanishing as soon as its argument becomes negative. When this occurs, the Fermi function has a step singularity, which reflects itself in a “kink” in the kernel as a function of at a critical value . This kink, in turn, considerably affects the large- behaviour of the Fourier transform , to be discussed below. The value of is determined as follows. When , , and , the argument of the Fermi function in Eq. (2)

 E+(k;Q)= ⎷(k22m+Q22m−μ)2+|Δ|2+km⋅Q (5)

first approaches zero for , such that with . Setting , the condition for the expression (5) to vanish becomes

 x2−2x(εk+μ)+(εk−μ)2+|Δ|2=0 (6)

where , whose solutions

 x±(k)=εk+μ±√4εkμ−|Δ|2 (7)

are both acceptable provided for . However, only attains a minimum for , in correspondence to which . When , on the other hand, the argument of the Fermi function in Eq. (2) never vanishes and the kernel is a smooth function of . [Recall in this context that, at the mean-field level, the zero-temperature chemical potential changes its sign at about the coupling value .]

The critical value for the kernel differs, in general, from the critical value of the Landau criterion for superfluidity associated with pair-breaking excitations, which is given by the expression (cf., e.g. Section 4.6 of Ref. Spuntarelli-2010 ()):

 (QLc)2m=√μ2+|Δ|2−μ. (8)

The value of approaches only asymptotically in the weak-coupling (BCS) limit when and .

Figure 1 shows the dependence on the coupling of the two critical wave vectors and , which are obtained using the mean-field values of and at zero temperature.

At finite temperature, on the other hand, the Fermi function in Eq. (2) is a smooth function of its argument, resulting in a smooth dependence of the kernel on . This feature, in turn, will make the dependence of on less problematic than at zero temperature.

Small- behaviour. In the homogeneous case with a uniform gap parameter, the gap equation (1) reduces to the form:

 −m4πaF=K(Q=0)≡I0 (9)

with the notation introduced in Ref. Simonucci-2014 () (in which we now set , , and ). This condition identifies the value of at self-consistency. In particular, one sees that is positive (negative) on the BCS (BEC) side of the crossover where (), and vanishes at unitarity where . Near , decreases quadratically in with coefficient

 d2K(Q)dQ2∣∣∣Q=0=−2I1m (10)

where the notation was also introduced in Ref. Simonucci-2014 () (in which we again set , , and ). This expansion up to quadratic order in about was utilized in Ref. Simonucci-2014 () to derive the (differential) LPDA equation from the (integral) NLPDA equation.

It was also shown in Ref. Simonucci-2014 () that, at zero temperature, both and can be calculated analytically in terms of elliptic integrals according to the method of Ref. Marini-1998 (). A plot of as a function of the coupling (given in Fig. 7 of Ref. Simonucci-2014 ()) shows that is a monotonically decreasing function of coupling from the BCS to the BEC limits. These results were also used in Ref. Brand-2014 () to apply the LPDA equation to the study of the snake instability of dark solitons.

Large- behaviour. For large values of , the term in Eq .(5) can be neglected irrespective of the sign of , such that . In addition, becomes negative only for positive values of when . At zero temperature, the contribution to the kernel for large originating from the presence of the Fermi function then reads:

 ∫dk(2π)3fF(E+(k;Q))E(k;Q)≃mQ2∫dk(2π)3Θ(√2mμ−|k+Q|) (11) = m(2mμ)3/26π2Q2

where is the Heaviside unit step function of argument . This contribution is sub-leading with respect to the remaining part of the integral that defines the kernel for large , namely,

 ∫dk(2π)3{12E(k;Q)−mk2}≃∫dk(2π)3{mk2+Q2−mk2} (12) = −mQ22π2∫∞0dkk2+Q2=−m|Q|4π.

The expression (12) thus gives the leading contribution to for large , irrespective of coupling. This asymptotic result remains valid even at finite temperature, provided that .

For later convenience, we identify an asymptotic kernel defined for all values of by the expression (12), namely,

 K∞(Q)=−m|Q|4π (13)

irrespective of coupling and temperature.

Angular integration. At any temperature, the -integration that enters the definition of the kernel can be reduced to a one-dimensional numerical integration over , by performing analytically the integration over the angle . This is done with the use of the standard integral (where and ):

 ∫dx1aeλx+b=1b[x−1λln(aeλx+b)]. (14)

From Eq. (2) one then obtains the expression:

 K(Q)=∫∞0dkk2(2π)2 (15) ×{1E(k;Q)[mkBTkQln(e(E(k;Q)+kQ/m)/kBT+1e(E(k;Q)−kQ/m)/kBT+1)−1]−2mk2}

where as before.

Profile of for various couplings and temperatures. Figure 2 shows the wave-vector dependence of the kernel calculated numerically from the expression (15) in the range , for various couplings and temperatures from up to close to . These plots confirm the downward quadratic dependence of near for all couplings and temperatures, as well as the presence of a kink at at zero temperature for coupling values before the chemical potential changes its sign (i.e., for ). One notes further from this figure that, as soon that the Fermi function in Eq. (2) becomes smooth for increasing temperature, the kink singularity in is also smoothed out.

In addition, Fig. 3 shows the wave-vector dependence of the kernel over a more extended range of Q in the case of unitarity at zero temperature, evidencing how converges, in practice, for large enough to its asymptotic expression (12) (or, else, to the asymptotic kernel (13)).

B. Properties of the modified kernel

The asymptotic behaviour (12) of the kernel suggests us that, for the later purpose of calculating the Fourier transform , it is convenient to interpret the kernel in the sense of distributions Lighthill-1958 (). This is done by introducing a test function of the Gaussian form , that can be included in the kernel itself by the definition:

 Kσ(Q)=K(Q)e−Q2/σ2. (16)

It is further understood that the limit will be taken (at least formally) only after the Fourier transform

 Kσ(R)=∫dQπ3e2iQ⋅RKσ(Q) (17)

will be calculated. Although we shall consider values of up to , it will turn out that (or even less) will be sufficient for most purposes. Similarly, a related definition

 Kσ∞(Q)=K∞(Q)e−Q2/σ2 (18)

is also introduced for the asymptotic kernel (13).

Spatial oscillations on the BCS side at zero temperature and the effect of temperature. Figure 4 shows the spatial profile of the kernel (multiplied by ) for the same couplings and temperatures of Fig. 2 and with a common value of . Note that, at and for the couplings , the kernel presents regular oscillations of wave vector , which get quickly damped as soon as the temperature is increased. No oscillations are instead present at any temperature for the coupling , as expected.

Long-range tail of the kernel . To identify the decay rate of the amplitude of the oscillations of the kernel at large , Fig. 5 shows multiplied by over an extended range of for two couplings at . In this way, the amplitude of is found to behave asymptotically like for both couplings. In Appendix A, by studying a model function in -space whose Fourier transform in -space can be evaluated analytically, we will verify that the tail of the kernel stems from the fact that, in three dimensions, the kink singularity of the kernel extends over a sphere of finite radius. Consistently, for the asymptotic kernel for which the kink singularity reduces to the single point , the power-law dependence of the tail of its Fourier transform becomes and no oscillation occurs in this case (cf. also Eq. (22) below).

Behaviour of near the origin. The behaviour of the kernel near the origin depends on the large- behaviour of the kernel . This, in turn, depends on but not on coupling and temperature, as it can be seen from the expressions (12) and (16). Figure 6 shows typical profiles of the kernel in the restricted range , for , , and various values of . One notes from Fig. 6 that, by increasing the value of from one panel to the next by a factor of two, a large increase results in the value of when (which eventually leads to a divergence when ). One can also verify that the zero of closest to the origin occurs at , thereby approaching when . We have also verified that, in the restricted spatial range of Fig. 6, plots with given but different values of and can hardly be distinguished from each other.

Numerical checks on the overall shape of . The following identities hold for the kernel at any coupling and temperature below :

 4π∫∞0dRR2Kσ(R)=K(Q=0)=I0, 16π∫∞0dRR4Kσ(R)=−∇2(K(Q)e−Q2/σ2)Q=0 = −∇2K(Q)|Q=0+6σ2K(Q=0)=6I1m+6I0σ2, (20)

where the notation of Eqs. (9) and (10) has been used. These identities can be used as a check on the overall shape of the kernel obtained numerical for a given value of . These checks are important especially on the BCS side of unitarity at , where the kernel has rapid oscillations with a slowly decaying amplitude. In this context, special care requires the -integral on the left-hand side of Eq. (20), which is understood to contain also a Gaussian weight of the form (such that at the end of the calculation). Otherwise, this integral would not converge owing the tail of . Accordingly, we shall interpret the integral

 ∫dRπ3e2iQ⋅Re−R2/σ2R=(σR√π)3e−Q2σ2R⟶δ(Q) (21)

as approaching the Dirac delta function in the sense of distributions when Lighthill-1958 ().

Figure 7 shows the ratio of the left-hand side (lhs) to the right-hand side (rhs) of both Eq. (II) (upper panel) and Eq. (20) (lower panel) for different values of when and (for this calculation the value has been used in accordance with the above argument). The steady convergence of these results gives us confidence about the stability of our numerical calculations of the kernel for increasing (notwithstanding the divergence of this kernel at for increasing - cf. Eq. (24) below). In this context, we have also verified numerically that, for growing , the profile of tends uniformly toward an asymptotic profile, with a convergence rate that becomes slower as gets close to (where the -dependent singular behaviour shown in Fig. 6 occurs).

The kernel . The asymptotic kernel (13) in -space is a particular case of the model function studied analytically in Appendix A, where one sets and in Eq. (31) therein. The analysis carried out in Appendix A for the Fourier transform then yields for the function in the limit :

 Kσ∞(R)=im√πσ38π3R (22) × d2dy2[e−y24(erfc(iy2)−erfc(−iy2))]y=2σR

for all values of . Knowledge of the overall profile of given by Eq. (22) will be useful in what follows.

The profile of the kernel is shown in Fig. 8 (middle panel) over an extended range of . For comparison, the profile of the model function studied analytically in Appendix A (from which is obtained in the limit ) is also shown in Fig. 8 (upper panel), when equals the critical value of at unitarity taken from Fig. 1. In both panels, .

Knowledge of the limiting behaviours of will also be useful in what follows. When we obtain from Eq. (40) with and :

 Kσ∞(R→∞)≃m8π3R4 (23)

which shows no oscillatory behaviour, contrary to the expression (40) with . When we obtain from Eq. (42), again with and :

 Kσ∞(R→0)≃−mσ42π3 (24)

which, apart from a minor correction, diverges like the expression (42) in the limit .

Recovering the exponential behaviour of the kernel . It is apparent from Fig. 2 that the kink singularity of the kernel , which occurs at on the weak-coupling (BCS) side of the crossover, disappears either by increasing the temperature toward or by moving to the BEC side of the crossover even at . In both cases, the kernel becomes a smooth function of . To mimic this behaviour, we consider the following simple model function:

 G(Q)=−α√P20+Q2 (25)

with , which (apart from an overall constant shift) has the same kind of small- and large- behaviour of the full kernel . We then multiply this function by as it was done in Eqs. (16) and (18), and calculate numerically the spatial Fourier transform of the ensuing function in three dimensions. The result is shown in Fig. 8 (lower panel ) for and equal to the value of taken from Fig. 1 at unitarity. From this plot one concludes that the large- behaviour of the function has the form where the length is of the order of . It is interesting to note that this behaviour coincides with that of the kernel yielding the linear terms in the Gor’kov derivation of the Ginzburg-Landau equation Gorkov-1959 (), which is valid on the BCS side of the crossover and close to FW-1971 ().

These results imply that when passing, from a function like that given by Eq. (31) with a kink singularity on the real -axis, to a function like that given by Eq. (25) where the branch-cut singularity resides in the complex plane at , the large- behaviour of the corresponding Fourier transforms exhibits a drastic change, from the behaviour of Eq. (40) to the behaviour of .

C. The regularized kernel

From Eqs. (24) and (42) the kernels and are seen to tend to a common value when , provided is large enough. We have verified numerically that this result remains true over a finite (albeit small) range of , where and are seen to essentially coincide with each other. An example is shown in Fig. 9 for a given value of when and . Since the small- profile common to and is associated with the common large- behaviour of the corresponding kernels in -space, which is independent of coupling and temperature, the results of Fig. 9 suggest us to adopt the following procedure which allows us to concentrate on the large- behaviour of the kernel , that instead depends on both coupling and temperature.

For given coupling and temperature, we move from outwards and search for the point at which the two kernels and start to deviate from each other, say, within . This can be done for a set of values of , thus monitoring the convergence of the results for . The results of this procedure are shown in Fig. 10 throughout the BCS-BEC crossover, for several temperatures and three different values of . In all cases, the values of are not larger than , which represents a small length scale compared with the overall spatial extent of the kernel footnote-1 (). Note also that depends weakly on for all couplings and temperatures. This result is remarkable, in light of the fact that in the small- region (compared to ) both kernels and change instead considerably by varying . For instance, from Fig. 6 the positions of both the first zero and the first maximum in are seen to decrease by a factor of two from to , while from Fig. 10 is seen corrispondigly to change only by a few percents.

Once the point is identified by this procedure, for given we define a regularized kernel as follows:

 Kσreg(R)={0(R

In this way, a “hole” about is effectively introduced in the original kernel , thereby avoiding its strong divergence for but at the same time not affecting the determination of its spatial range, for which the behaviour when is irrelevant.

D. Spatial range of the kernel as a function of coupling and temperature

To determine the spatial range of the kernel (through its regularized version ) as a function of coupling and temperature, it is convenient to distinguish two cases when the kernel for large has:

(i) An oscillatory behaviour, like at from the BCS to the unitary regime;

(ii) An exponential behaviour, like at in the BEC regime or when approaching even in the BCS and unitary regimes.

For case (i), the spatial range of the kernel is determined by considering the behaviour of the function:

 Fσ(R)=∫RR0dR′R′2Kσ(R′)=∫R0dR′R′2Kσreg(R′). (27)

For given , this function converges asymptotically to a finite value in the limit , and this is so even at when its integrand has characteristic oscillations of wave vector with amplitude decaying like for large . Typical examples of the behaviour of vs are shown in Fig. 11 for , , and . At , one sees that has essentially converged to its asymptotic value as soon as it reaches the first maximum at , past which shows only a minor oscillatory behaviour around . At and , on the other hand, shows no oscillatory behaviour and reaches monotonically the asymptotic value at about the same value of identified at . This implies that, apart from the presence or absence of minor oscillations, for both couplings the overall shape of remains essentially the same upon varying the temperature. Accordingly, for coupling values we identify the value of obtained by the above procedure with the spatial range of the kernel at (expecting further that this range should only slightly depend on temperature).

For case (ii), the product is found to behave like for large , such that an exponential fit can be made directly on this product to extract the characteristic length . We have performed this fit, at for and at across the whole BCS-BEC crossover, again with the value . To connect with continuity the values of obtained here at for with the values of obtained previously for , we have set and determined the constant in such a way that upon approaching the BEC limit. By carrying out these calculations up to , we have obtained the value .

Figure 12 shows the coupling dependence throughout the whole BCS-BEC crossover of the range obtained at by the above “mixed” procedure (dots), together with the coupling dependence of the range obtained by the above exponential fit at (stars). In addition, the inset of Fig. 12 combines these data in the coupling dependence of the ratio (full line), which turns out to be about for all couplings in the range . Such a weak temperature dependence from up to , that we have obtained for the range of the kernel irrespective of coupling, is in line with the behaviour of the size of the Cooper pairs obtained in Ref. Palestini-2014 (). The inset of Fig. 12 shows also the coupling dependence of the quantity for (dashed line), which in the BCS limit and close to is expected to equal unity according to an analytic result due to Gor’kov FW-1971 (). Remarkably, our numerical calculations approximately reproduce this result not only in the BCS limit but also across unitarity.

E. Comparison of the length scales and

The coupling behaviour of reported in Fig. 12 is reminiscent of the coupling behaviour of the pair coherence length throughout the BCS-BEC crossover, which was obtained originally at in Ref. Pistolesi-1994 () in terms of the pair correlation function of opposite-spin fermions at the mean-field level. For the present purposes, it can be useful to relate at also with the approximate expression , where is the Landau critical wave vector given by Eq. (8). In the weak-coupling (BCS) limit, this expression yields , which coincides with the limiting value of obtained in Ref. Pistolesi-1994 () at . In the strong-coupling (BEC) limit, on the other hand, equals and vanishes in the relevant limit . In addition, in the intermediate coupling region the expression approximates reasonably well the values of at obtained in Ref. Pistolesi-1994 ().

Figure 13 compares the coupling dependence of at , as obtained in Ref. Pistolesi-1994 () (full line) and by the approximate expression (dashed line), with the coupling dependence of for (diamonds) taken from Fig. 12. Although the two quantities and have been obtained through quite different working procedures, the resemblance here between their coupling dependence appears to be rather remarkable. This resemblance persists also at finite temperatures up to , as shown by the results presented in Fig. 12.

To summarize, the above results are all consistent with one’s physical expectation that the range of the kernel of the non-local gap equation (1) should (at any temperature) be directly related to the size of the Cooper pairs, which represents the fundamental length scale of the BCS pairing theory for fermionic superfluidity deGennes-1966 ().

## Iii Numerical solution of the NLPDA equation for an isolated vortex

The validity of the (differential) LPDA equation was tested in Ref. Simonucci-2014 () for an isolated vortex embedded in an infinite superfluid, for which an accurate solution of the BdG equations is available to compare with across the whole BCS-BEC crossover for all Simonucci-2013 (). In Ref. Simonucci-2014 () deviations between the LPDA and BdG calculations were found in the BCS regime at low temperature, and their origin was attributed to the finiteness of the spatial range of the kernel of the (integral) NLPDA equation from which the LPDA equation was obtained in Ref. Simonucci-2014 () at a final step. No detailed analysis, however, was made in Ref. Simonucci-2014 () about the consequences of the finiteness of this spatial range, when solving for the gap parameter with a nontrivial spatial profile.

Here, we consider again the case study of an isolated vortex embedded in an infinite superfluid, for which the results of the NLPDA equation (1) can be tested against the results of the BdG equations and also compared with the results of the LPDA equation, over an extended region of the coupling-vs-temperature phase diagram. To this end, a strategy needs to be implemented to solve numerically the NLPDA equation in an efficient way.

A. Vortex solution

It was shown in Section II (cf. Figs. 2 and 3) that the kernel of the NLPDA equation (1) has a rather simple form in -space. This kernel starts at small as an inverted parabola with coefficients given by Eqs. (9) and (10) and ends up at large with the linear behaviour (12), and in between has a kink singularity at the critical wave vector (cf. Fig. 1) for and . Due to this kink singularity, the Fourier transform of this kernel in -space has an oscillatory behaviour with a slowly decaying tail for large , while the large- behaviour (12) results in a strong singularity at .

To avoid those features that can cause problems in the numerical solution of the integral equation, out of the two versions in which this equation can be written according to the definition (4), namely,

 −m4πaFΔ(r) = ∫dRΔ(R)K(r−R|r) (28) = ∫dQπ3e2iQ⋅rΔ(Q)K(Q|r)

it appears convenient to use the second version in -space. [According to the arguments of Section II, it is understood that the modified kernel of Eq. (16) enters Eq. (28), although this step is not strictly necessary when solving the gap equation where limits the relevant range of to values of much smaller than . In addition, the suffix has been dropped from Eq. (28) since we consider here the case with .]

To achieve self-consistency of the solution, the choice of the -version of Eq. (28) requires us to transform the profile of the gap parameter back and forth from - to -space as many times as needed. To this end, an efficient method is required not to loose numerical accuracy in the course of the repeated transformations. A method to fulfil this purpose is described in detail in Appendix B.

This method is especially suited when the symmetry of the problem reduces the - and -integration to one dimension. This is the case of an isolated vortex with cylindrical symmetry embedded in an infinite superfluid, for which the gap parameter takes the form:

 Δ(R)=Δ(ρ,φ,Rz)=Δ(ρ)eiφ (29)

where . Its Fourier transform then reads:

 Δ(Q)=∫dRe−2iQ⋅RΔ(R)=πδ(Qz)Δ(Q)eiφQ (30)

where and are the radial and azimuthal coordinates in the