text:
# Model of Inhomogeneous Impurity Distribution in Fermi Superfluids

## I Introduction

The standard treatment of impurities in a metal is based on averaging over the locations of the impurities.AGD  This means that the probability of a quasiparticle being scattered is independent of its location. To be definite, we call this the homogeneous scattering model (HSM). The purpose of this paper is to study the case where the impurity distribution varies in space. We define a model of inhomogeneous scattering, which is as simple as possible and represents a medium that is uniform and isotropic on a large scale. We call it the isotropic inhomogeneous scattering model (IISM).ErkkiModels  The model is motivated by experiments on superfluid 3He in a porous aerogel, and all our numerical computations concentrate on this case. However, the model is independent of the pairing symmetry and therefore can equally be applied, for example, to $`s`$\- or $`d`$-wave superconductors.

Liquid 3He is an ideal example of unconventional superfluid because it is naturally pure, it has a spherical Fermi surface, and its $`p`$-wave pairing state is well understood. Therefore it is of interest to study impurity effects in this superfluid. The addition of impurities to helium can be done by using porous aerogel so that typically 98% of the volume is occupied by 3He. It was found experimentally that the superfluid transition temperature in aerogel is reduced but remains sharp.Porto95  Also other superfluid properties such as the superfluid density and the NMR shifts were found to be modified by aerogel. Many experiments studying this system have been made during the last seven years.halprev

The HSM is an attractive model for 3He in aerogel because of its simplicity. Essentially all theoretical calculations for this system are based on it. ErkkiModels ; baramidze ; M ; ErkkiGL ; PhysicaB ; raineraero ; BK2 ; higashitani ; SharmaSauls ; mineevkrotkov ; Baramidze2 ; SharmaLT23 ; yamamoto  However, already the first comparisons to experiment showed that the HSM is insufficient quantitatively. In particular, it was found that the order parameter is more suppressed than the superfluid transition temperature $`T_\mathrm{c}`$. A natural explanation for this comes from the fact that the scattering in aerogel is not homogeneous, as has already been demonstrated using the IISM.ErkkiModels  Similar results have been reached in calculations on unconventional superconductors. FKBS

Preliminary results of the IISM have been reported in Refs. ErkkiModels,  and LT22, . In this paper we present the IISM in detail. Our studies are based on the quasiclassical theory of Fermi liquids. We discuss the assumptions of this theory (Sec. II), and how these are satisfied for 3He in aerogel (Sec. III). The inhomogeneous scattering model is introduced in Sec. IV. For the case of $`p`$-wave pairing we introduce the order parameters of A- and B-type phases (Sec. V). We calculate several quantities including the critical temperature, order parameter, and superfluid density (Sec. VI). The results are compared with experiments on 3He in aerogel (Sec. VII). The equations of quasiclassical theory and the details of calculations are discussed in the Appendix.

## II Quasiclassical theory

On a microscopic scale a pure system is described in terms of particles (conduction electrons or 3He atoms) and their interaction. Because the interactions are strong, this leads to a complicated many-body problem.

The characteristic length in superconductivity or superfluidity is the coherence length. We define this quantity as

$$\xi _0=\frac{\mathrm{}v_\mathrm{F}}{2\pi k_\mathrm{B}T_{\mathrm{c0}}},$$
(1)
where $`v_\mathrm{F}`$ is the Fermi velocity. To be precise, we have used the superfluid transition temperature $`T_{\mathrm{c0}}`$ of a pure system. The coherence length is typically much larger than the Fermi wave length $`\lambda _\mathrm{F}=2\pi /k_\mathrm{F}`$. In 3He $`\xi _0`$ depends on pressure and changes from 16 nm at the solidification pressure to 77 nm at the vapor pressure whereas $`\lambda _\mathrm{F}0.7`$ nm.

The theory that is designed to work on the scale $`\xi _0\lambda _\mathrm{F}`$ is the quasiclassical theory.Serene  It treats the system as a dilute gas of weakly interacting quasiparticles. In quasiclassical theory all the many-body physics that takes place on the microscopic scale $`\lambda _\mathrm{F}`$ is eliminated. It only appears through phenomenological parameters like the Fermi-surface, Landau Fermi-liquid parameters and transition temperature $`T_{\mathrm{c0}}`$.

Let us consider any external objects in the system. These objects are characterized by a strong potential, on the order of the Fermi energy. As a consequence the state of the system is modified in the vicinity of the object. A theoretical analysis of these atomic scale changes is again difficult because of the strong interactions between particles.

In the quasiclassical theory the effect of external objects is twofold.BR  First, the phenomenological parameters discussed above are changed. These parameters are determined by processes on the Fermi-energy scale and are therefore of short range. Assuming that the surface area of the external objects times the atomic length scale ($`\lambda _\mathrm{F}`$) is a small fraction of the total volume, this effect is small and is neglected in the following. Second, the objects affect the low energy processes (energy $`k_BT_{\mathrm{c0}}`$) directly via a scattering of quasiparticles. The range of this effect is long, on the order of the coherence length \[Eq. (1)\]. Therefore it leads to substantial modification of the superfluid properties.

An important length characterizing the scattering is the mean free path $`\mathrm{}`$. The simplest case is to consider the limit $`\lambda _\mathrm{F}/\mathrm{}0`$. This is the quasiclassical limit where the Fermi wave length effectively disappears from the theory. First order corrections in $`\lambda _\mathrm{F}/\mathrm{}`$ lead to effects like weak localization, which are neglected here. In quasiclassical theory, the scattering is represented by a collection of “scattering centers.” The main assumption is that the quantum interference between two scattering centers is neglected. Technically this can be achieved considering an ensemble average where the locations $`𝐫_i`$ of the scattering centers are uncertain by a distance on the order of $`\lambda _\mathrm{F}`$ or more. The size of one scattering center is limited by the condition that it has to be localized on the scale of $`\xi _0`$. The scattering properties of a center can be parametrized by scattering phase shifts $`\delta ^{(l)}`$, which are taken at the Fermi energy in the normal state. (For simplicity we label the different partial waves by a single index $`l`$, but there is no need to restrict to spherically symmetric scattering centers. The scattering could also be spin-dependent, but it is also neglected here for simplicity.) Thus a complete description of the scattering needed in the quasiclassical theory consists of distribution functions $`n_i(𝐫)`$ and scattering phase shifts $`\delta _i^{(l)}`$ of the scattering centers $`i=1,2,\mathrm{}`$.

Within each scattering center an exact quantum treatment is allowed in principle. However, because of uncertainty about the microscopic processes at surfaces, the phase shifts $`\delta _i^{(l)}`$ cannot be calculated from first principles. Instead, one has to use some models, for example hard spheres.Merzbacher  In some quantities the phase shifts only appear in certain combinations, for example, the transport cross section

$$\sigma =\frac{4\pi }{k_\mathrm{F}^2}\underset{l=0}{\overset{\mathrm{}}{}}(l+1)\mathrm{sin}^2(\delta ^{(l+1)}\delta ^{(l)}).$$
(2)
Let us try to clarify some consequences of the assumptions made above. The energy one usually is trying to calculate is on the order of $`f_{\mathrm{cond}}\xi _0^3`$. Here $`f_{\mathrm{cond}}k_\mathrm{B}T_\mathrm{c}/\lambda _\mathrm{F}^2\xi _0`$ is the superfluid condensation energy per volume, and a typical volume $`\xi _0^3`$. Above the individual scattering centers were required to be small, $`\sigma \xi _0^2`$. This implies that the energy $`f_{\mathrm{cond}}\sigma \xi _0`$ associated with a single impurityTKR  is small in comparison. The typical number of impurities in volume $`\xi _0^3`$ is $`N\xi _0^2/\sigma `$, which for random impurities implies a fluctuation $`\delta N\sqrt{N}`$. The corresponding fluctuation in energy is by factor $`\sqrt{\sigma }/\xi _0`$ smaller than $`f_{\mathrm{cond}}\xi _0^3`$. This has to be neglected since there are other neglected contributions that are on the same order of magnitude. Thus impurity averaging in the quasiclassical approximation implies a scattering medium where the fluctuations in the impurity density are neglected.

## III Aerogel

The structure of aerogel, as relevant for 3He experiments, is discussed in Ref. Porto99, . Here we repeat some main points. Aerogel consists of small SiO2 particles of diameter $`3`$ nm, which are coalesced together to form a self-supporting structure. Experiments with 3He typically use aerogels with open volume fraction 98% or more. According to small-angle x-ray scattering measurements, there is a “fractal” range in the particle cluster up to a “correlation length” $`\xi _a100`$ nm. Above this scale the structure looks homogeneous. Computer simulations give a picture of widely spaced aerogel strands implying a long mean free path $`>100`$ nm.

It seems reasonable that the quasiclassical description above can be applied to liquid 3He in aerogel. The atomic layer on the SiO2 surfaces occupies only one per cent of the volume and can be neglected for many purposes. (Magnetic properties make an exception because the susceptibility in this layer is much larger than in pure liquid.Sprague95 ) Since a major part of the liquid is within a coherence length from SiO2,Porto99  the scattering effect leads to a substantial modification of the superfluid properties. It also seems that the scattering from aerogel can be represented by incoherent scattering centers whose size is small compared to $`\xi _0`$ because random variations on the order of $`\lambda _\mathrm{F}`$ are likely to develop already at much smaller distances.

The smallest reasonable choice for a scattering center is a single SiO2 particle (diameter $`2R_a3`$ nm). This is large in comparison to $`\lambda _\mathrm{F}`$: $`k_\mathrm{F}R_a10`$. According to hard sphere phase shifts this means that only $`1\%`$ of the scattering takes place in the s-wave channel, and $`99\%`$ is left to higher partial waves in Eq. (2). This dominant contribution of higher partial waves has several important consequences.

Firstly, the phase shifts $`\delta _i^{(l)}`$ are random numbers. This is because even the phase shifts (modulo $`\pi `$) of a hard sphere with a fixed $`R_a`$ are pseudorandom numbers for $`k_\mathrm{F}R_a1`$. Adding to this the varying particle size, the surface roughness and varying orientations of touching neighbor particles, it is simply impossible that the result would be anything else but random. (The randomness, of course, is valid only for partial waves $`l`$ that contribute essentially to scattering, i.e., for $`l<k_\mathrm{F}R_a`$.)

A consequence of the random phase shifts is that only the number of scattering phase shifts, or equivalently, the cross section $`\sigma `$ \[Eq. (2)\] is important in describing a scattering center. Thus a sufficient description of the scattering is obtained by specifying only $`\sigma _i`$ and $`n_i(𝐫)`$.

A second important consequence of the large particle size is that both $`\sigma _i`$ and $`n_i(𝐫)`$ are independent of pressure. This is crucial for comparison with experiments, because when fitting is needed, it can be done at one pressure only, and the predictions of the model get fixed at all pressures. We argue as follows. It is reasonable to assume that the aerogel is independent of the hydrostatic pressure, implying that $`n_i(𝐫)`$ is also. The Fermi wave vector changes by 10% over the pressure range from the vapor pressure to the solidification pressure. This could induce some pressure dependence in the cross section $`\sigma _i`$. For example, hard-sphere $`\sigma `$ depends essentially on $`k_\mathrm{F}`$ in the region $`k_\mathrm{F}R_a1`$.Salomaa  However, in the physical region $`k_\mathrm{F}R_a10`$ the dependence of $`\sigma `$ on $`k_\mathrm{F}`$ is very weak. Thus we conclude that both $`\sigma _i`$ and $`n_i(𝐫)`$ are independent of pressure.

The large $`k_\mathrm{F}R_a`$ is potentially bad news for theory because calculations that take higher partial waves into account are very complicated, see Ref. PhysicaB, . The promising conclusion of these calculations is that, at least in some cases, the results including higher partial waves are almost identical to those including only s-wave under the following conditions: (i) one uses the same transport mean free path $`\mathrm{}`$ and (ii) one uses either random phase shift $`\delta ^{(0)}`$ or fixed $`\mathrm{sin}^2\delta ^{(0)}\frac{1}{2}`$ in the s-wave calculation. Here we assume that this correspondence holds more generally. Thus we calculate only s-wave and present results for $`\mathrm{sin}^2\delta ^{(0)}=\frac{1}{2}`$. Finally, the randomness of the phase shifts also simplifies the numerical calculations since it implies that some components of the propagator vanish (see the Appendix).

The use of several impurity species is important in studying anisotropic scattering where the preferred direction varies in space. In this paper we limit to isotropic scattering. In order to simplify the notation we select all scattering centers to have equal cross section, so that the only scattering parameters are $`\sigma `$ and the total impurity density $`n(𝐫)`$.

## IV Setting the model

The simplest possible impurity profile $`n(𝐫)`$ is a constant. It implies a location independent (transport) mean free path $`\mathrm{}=(n\sigma )^1`$. This homogenous scattering model (HSM) has successfully been used to model the impurities in superconductors.AGD ; general  For impure $`p`$-wave superfluid (3He) it has been used to calculate the critical temperature,larkin  the order parameter and superfluid densityErkkiModels ; PhysicaB ; raineraero ; higashitani , properties near the superfluid transition,M ; ErkkiGL  properties in magnetic field,baramidze ; BK2 ; SharmaSauls ; mineevkrotkov ; yamamoto  density of states,SharmaSauls ; mineevkrotkov  thermal conductivitySharmaLT23  and estimate strong-coupling effects.Baramidze2

The predictions of the HSM are compared to experiments in Refs. ErkkiModels, ; ErkkiGL, ; raineraero, ; higashitani, ; SharmaSauls, ; mineevkrotkov, ; SharmaLT23, ; Porto99, ; Hook,  and Lawes, , see also below. Compared to pure 3He, the HSM gives the right tendency and can work even quantitatively in some cases, but there can be differences up to a factor of 5 in the suppression factors (see below). In some papers the experiments are compared with data that is calculated in the unitary limit of s-wave scattering, which seems to give better agreement than other phase shifts. As explained above, we believe this is misguided, and random or intermediate phase shifts should be used instead. The observed differences can more naturally be explained by inhomogeneous scattering as will be shown below, at least for some quantities.

Real aerogel has voids where scattering is negligible. This can be modelled by an impurity density $`n(𝐫)`$ which depends on the location $`𝐫`$. It is in principle possible to use a realistic $`n(𝐫)`$ for aerogel. This has the drawback that the computational volume should be large in order to get a representative sample, and this implies heavy numerical effort. Here we prefer the opposite limit of a simple model $`n(𝐫)`$. The simplest possibility would be a plane-wave variation $`n(𝐫)=n_0+n_1\mathrm{cos}(𝐪𝐫)`$. A stronger version of this would be equally spaced scattering planes. In the limit of very strong scattering in the planes this leads to isolated slabs.ErkkiModels  The problem with all these models is that they are anisotropic. For example, the superfluid density would depend on the direction of the superfluid velocity $`𝐯_s`$. In comparison to experiments one should use some average over the directions, but this neglects the process how the averaging really takes place by a nonuniform current distribution.

The purpose of the isotropic inhomogeneous scattering model (IISM) (Ref. ErkkiModels, ) is to incorporate a non-constant $`n(𝐫)`$ with spherical symmetry. We take a spherical volume of radius $`R`$ and use an impurity density $`n(r)`$ that depends only on the radial coordinate $`r`$. We assume that these spheres fill all the space. This last point is not strictly possible, but represents an approximation that is similar to using spherical approximation for a Wigner-Seitz unit cell.AM  The calculation in a single unit cell can represent states where the superfluid order parameter $`\stackrel{}{A}(𝐫)`$ has the Bloch form

$$\stackrel{}{A}(𝐫)=\stackrel{}{A}^{(0)}(𝐫)\mathrm{exp}(i𝐪𝐫),$$
(3)
where $`\stackrel{}{A}^{(0)}(𝐫)`$ is a strictly periodic order parameter. In the present case the wave vector $`𝐪`$ is the imposed phase gradient that is related to the superfluid velocity $`𝐯_\mathrm{s}=(\mathrm{}/2m)𝐪`$ defined on a scale larger than $`R`$.

At the surface of the IISM sphere we impose the boundary condition that an exiting quasiparticle is effectively returned to the sphere at the diametrically opposite point, see Fig. 1. In the case of current-carrying states, there has to be a phase shift which corresponds to the phase factor in Eq. (3). Otherwise the state of the quasiparticle (momentum, spin) is unchanged. The boundary condition can be expressed mathematically for the Green’s function as

$`\stackrel{˘}{g}(\widehat{𝐤},R\widehat{𝐫},ϵ_m)=\mathrm{exp}\left(\mathrm{i}𝐪R\widehat{𝐫}\stackrel{˘}{\tau }_3\right)\stackrel{˘}{g}(\widehat{𝐤},R\widehat{𝐫},ϵ_m)\mathrm{exp}\left(\mathrm{i}𝐪R\widehat{𝐫}\stackrel{˘}{\tau }_3\right),`$ (4)
see the Appendix for notation.

The model reduces to the HSM in the limit that $`n(r)`$ is independent of $`r`$.

Because of the spherical approximation, it is worth while to properly define all large-scale ($`R`$) quantities. We define macroscopic quantities as averages of the corresponding local quantities over the sphere,

$$A\frac{3}{4\pi R^3}_{r<R}d^3rA(𝐫).$$
(5)
In particular, corresponding to the local mass current density of the superfluid, $`𝐣_s`$, we have the macroscopic $`𝐣_s`$. The superfluid density $`\stackrel{}{\rho }_s`$ is then defined by $`𝐣_s=\stackrel{}{\rho }_s𝐯_s+O(v_s^2)`$ with $`𝐯_\mathrm{s}=(\mathrm{}/2m)𝐪`$ as defined above. The average mean free path $`\mathrm{}_{\mathrm{ave}}`$ is defined by $`\mathrm{}_{\mathrm{ave}}^1=\sigma n(r)`$.

The parameters specifying the scattering are the radius of the sphere $`R`$, the average mean free path $`\mathrm{}_{\mathrm{ave}}`$, and the shape of the impurity density $`n(r)`$. For $`n(r)`$ we use two different analytic forms

$`n(r)`$ $`=`$ $`c\left[\left({\displaystyle \frac{r}{R}}\right)^j{\displaystyle \frac{j}{j+2}}\left({\displaystyle \frac{r}{R}}\right)^{j+2}\right],`$ (6)
$`n(r)`$ $`=`$ $`c^{}\left[\mathrm{cos}^j\left({\displaystyle \frac{\pi r}{2R}}\right)+b\right],j2.`$ (7)
with parameters $`j`$ and $`b`$. Here the prefactors $`c`$ and $`c^{}`$ are determined by the parameters $`\mathrm{}_{\mathrm{ave}}`$ and $`\sigma `$ ($`c`$, $`c^{}>0`$). The functions in Eqs. (6) and (7) are shown in the inset of Fig. 2. We call them void and cluster profiles, respectively, because the former has strongest scattering at the boundary and the latter at the center of the sphere. Both profiles have zero derivative at $`r=R`$ in order to have a smooth impurity density everywhere.

An attractive feature of the IISM is that the symmetries of a homogeneous system are preserved on a large scale $`R`$. Thus one can apply phenomenological large-scale theories such as Ginzburg-Landau and hydrodynamic theories as for a homogeneous medium.Volovik ; ErkkiGL ; fomin  The only change is that the parameters of these theories are modified by the inhomogeneity. Some of these parameters are calculated below.

## V Order parameter

We apply the IISM for superfluid 3He. The main assumption in addition to those already mentioned is the use of the weak-coupling approximation. The dipole-dipole interaction is neglected because it is unimportant on the scale of a few $`\xi _0`$, which we study here. The equations and some details of the numerical implementation are discussed in the Appendix.

Even with all simplifications, the computational effort in the IISM is quite substantial. For example, the computer code has five nested loops in addition to the one needed for the iteration of the order parameter, and there are several stages of initialization, interpolation, and data collection.

Based on calculations with the HSM, no new superfluid phases of 3He are expected in the presence of scattering.ErkkiModels  However, the order parameters of the A and B phases are modified by inhomogeneous scattering. The general forms can be deduced using symmetry arguments. The B phase order parameter for $`𝐯_s=0`$ has the form

$`\stackrel{}{A}(r,\varphi ,\theta )=e^{\mathrm{i}\chi }\stackrel{}{R}\left[\mathrm{\Delta }_r(r)\widehat{𝐫}\widehat{𝐫}+\mathrm{\Delta }_a(r)(\widehat{𝜽}\widehat{𝜽}+\widehat{\mathit{\varphi }}\widehat{\mathit{\varphi }})\right],`$ (8)
where spherical coordinates ($`r,\varphi ,\theta `$) are used. The phase $`\chi `$ and the rotation matrix $`\stackrel{}{R}`$ are arbitrary constants. The calculation determines the real-valued radial $`\mathrm{\Delta }_r(r)`$ and angular $`\mathrm{\Delta }_a(r)`$ functions. These functions are shown in Fig. 2 for three different scattering profiles $`n(r)`$.

We see that the order parameter is inhomogeneous, and is mostly suppressed in regions where the scattering is strong. For the void profile \[Eq. (6)\] the order parameter components have maxima $`\mathrm{\Delta }_{\mathrm{max}}=\mathrm{\Delta }_r(0)=\mathrm{\Delta }_a(0)`$ at the center and minima at the surface. The absolute minimum $`\mathrm{\Delta }_{\mathrm{min}}`$ is given by $`\mathrm{\Delta }_a(R)`$.

We define two different averages of the spatially varying order parameter. The usual average is given by

$$\mathrm{\Delta }_{\mathrm{ave}}^2=\frac{1}{3}\mathrm{\Delta }_r^2+2\mathrm{\Delta }_a^2.$$
(9)
An “NMR average” is defined by

$$\mathrm{\Delta }_{\mathrm{NMR}}^2=\frac{1}{15}2\mathrm{\Delta }_r^2+6\mathrm{\Delta }_r\mathrm{\Delta }_a+7\mathrm{\Delta }_a^2.$$
(10)
This form can be justified as follows. The frequency shifts of nuclear magnetic resonance in superfluid 3He are determined by the dipole-dipole interaction energy Leggett

$$f_\mathrm{d}=g_\mathrm{d}[|Tr\stackrel{}{A}|^2+Tr(\stackrel{}{A}^{}\stackrel{}{A})].$$
(11)
The coefficient $`g_\mathrm{d}`$ is a phenomenological parameter that depends on a cut-off energy, which is on the order of the Fermi energy.btex  According to the principles of Sec. II we assume it is not changed by the impurity. We calculate the dipole-dipole energy \[Eq. (11)\] for the B-phase order parameter (8) and find

$$f_\mathrm{d}=4g_\mathrm{d}\mathrm{\Delta }_{\mathrm{NMR}}^2\mathrm{cos}\vartheta (1+2\mathrm{cos}\vartheta ),$$
(12)
where $`\vartheta `$ is the rotation angle of $`\stackrel{}{R}`$. This has exactly the same form as in pure 3He-B, the only change being that the energy gap $`\mathrm{\Delta }`$ is replace by $`\mathrm{\Delta }_{\mathrm{NMR}}`$. Thus the NMR properties within IISM are the same as in bulk liquid except for this renormalization.

We point out that the order parameter in Eqs. (8)–(17) is defined using the off-diagonal part of the mean-field self energy (25). It should be noted that the order parameter is not simply related to the energy gap in the excitation spectrum. Only in pure, homogeneous superfluid the energy gap equals $`\mathrm{\Delta }_{\mathrm{ave}}=\mathrm{\Delta }_{\mathrm{NMR}}`$. The excitation spectrum in the HSM has been studied in Refs. SharmaSauls,  and mineevkrotkov, , but there are no studies yet for the IISM.

The A-phase order parameter is more complicated since it is anisotropic. We select the cylindrical coordinates ($`\rho ,\varphi ,z`$) so that $`z`$ is along the anisotropy axis $`\widehat{𝐥}`$. The order parameter for $`𝐯_s=0`$ can be written as

$`\stackrel{}{A}(\rho ,\varphi ,z)=e^{\mathrm{i}\varphi }\widehat{𝐝}\left[\mathrm{\Delta }_\rho (\rho ,z)\widehat{𝝆}+\mathrm{i}\mathrm{\Delta }_\varphi (\rho ,z)\widehat{\mathit{\varphi }}+\mathrm{\Delta }_z(\rho ,z)\widehat{𝐳}\right],`$ (13)
where $`\widehat{𝐝}`$ is an arbitrary constant unit vector. The functions $`\mathrm{\Delta }_i`$ are real but now they depend on two coordinates $`\rho `$ and $`z`$. They satisfy symmetry relations $`\mathrm{\Delta }_\rho (\rho ,z)=\mathrm{\Delta }_\rho (\rho ,z)`$, $`\mathrm{\Delta }_\varphi (\rho ,z)=\mathrm{\Delta }_\varphi (\rho ,z)`$ and $`\mathrm{\Delta }_z(\rho ,z)=\mathrm{\Delta }_z(\rho ,z)`$. For the A phase we define

$`\mathrm{\Delta }_{\mathrm{ave}}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Delta }_\rho ^2+\mathrm{\Delta }_\varphi ^2+\mathrm{\Delta }_z^2`$ (14)
$`\mathrm{\Delta }_{\mathrm{NMR}}^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{\Delta }_\rho ^2+\mathrm{\Delta }_\varphi ^22\mathrm{\Delta }_z^2.`$ (15)
In pure homogeneous superfluid $`\mathrm{\Delta }_{\mathrm{ave}}=\mathrm{\Delta }_{\mathrm{NMR}}`$ equals the maximum energy gap in the A phase. The dipole-dipole energy \[Eq. (11)\] in the A phase (13) is given by

$$f_\mathrm{d}=2g_\mathrm{d}\mathrm{\Delta }_{\mathrm{NMR}}^2(\widehat{𝐝}\widehat{𝐥})^2,$$
(16)
and the justification for $`\mathrm{\Delta }_{\mathrm{NMR}}`$ (15) is completely analogous to the case of the B phase.

In the case of a finite superfluid velocity (taken to be in the $`z`$ direction) we limit our calculations to the B phase where the order parameter takes the form

$`\stackrel{}{A}(\rho ,\varphi ,z)=e^{\mathrm{i}\chi }\stackrel{}{R}\left[\mathrm{\Delta }_{\rho \rho }\widehat{𝝆}\widehat{𝝆}+\mathrm{\Delta }_{\varphi \varphi }\widehat{\mathit{\varphi }}\widehat{\mathit{\varphi }}+\mathrm{\Delta }_{zz}\widehat{𝐳}\widehat{𝐳}+\mathrm{\Delta }_{\rho z}\widehat{𝝆}\widehat{𝐳}+\mathrm{\Delta }_{z\rho }\widehat{𝐳}\widehat{𝝆}\right].`$ (17)
Now $`\mathrm{\Delta }_{\mu i}`$ are complex and depend on $`\rho `$ and $`z`$. They satisfy symmetry relations $`\mathrm{\Delta }_{\rho z}(\rho ,z)=\mathrm{\Delta }_{\rho z}^{}(\rho ,z)`$, $`\mathrm{\Delta }_{z\rho }(\rho ,z)=\mathrm{\Delta }_{z\rho }^{}(\rho ,z)`$ and $`\mathrm{\Delta }_{ii}(\rho ,z)=\mathrm{\Delta }_{ii}^{}(\rho ,z)`$, where $`i=\rho ,\varphi ,z`$. The mass supercurrent $`𝐣_\mathrm{s}(𝐫)`$ has non-zero $`\rho `$ and $`z`$-components inside the sphere. One case is illustrated in Fig. 3.

It can be seen that the current density is smallest in regions of strong scattering. Naturally, the current is conserved, $`\mathbf{}𝐣_\mathrm{s}(𝐫)=0`$. In the void profile, the current has to go through the scattering regions at the cell boundary, whereas in the cluster profile (not shown) the transport current can flow past the scattering region situated at the center. Regardless of the profile, the averaged current is parallel to $`𝐯_s`$ and independent of the direction. Thus the superfluid density tensor $`\stackrel{}{\rho }_s`$ reduces to a scalar $`\rho _s`$.

## VI Results

We start by studying the critical temperature $`T_\mathrm{c}`$. We calculate the quantity $`T_\mathrm{c}/T_{\mathrm{c0}}`$, the critical temperature relative to the critical temperature in the absence of scattering. In the HSM this depends only on the parameter $`\xi _0/\mathrm{}`$. The dependence turns out to be the same as calculated by Abrikosov and Gorkov for $`s`$-wave superconductors in the presence of magnetic impurities.AG  This result was generalized to the nonmagnetic $`p`$-wave case in Ref. larkin, . In the IISM $`T_\mathrm{c}/T_{\mathrm{c0}}`$ depends only on the ratio $`\xi _0/\mathrm{}_{\mathrm{ave}}`$ and the impurity profile $`n(r)`$. In particular, it is independent of the phase (A or B) and the phase shifts $`\delta ^{(l)}`$. The relative $`T_\mathrm{c}`$ is plotted as a function of $`\xi _0/\mathrm{}_{\mathrm{ave}}`$ in Fig. 4 for different scattering profiles.

One sees that in all cases the HSM gives the lowest $`T_\mathrm{c}`$. This is natural since the inhomogeneity implies existence of regions where the scattering is less than the average, and in these regions the order parameter nucleates at a higher temperature. For example, for the steep void profile with $`R/\mathrm{}_{\mathrm{ave}}=2`$ we find $`T_\mathrm{c}=0.8T_{\mathrm{c0}}`$ at the average impurity where superfluidity in the HSM is completely suppressed ($`\xi _0/\mathrm{}_{\mathrm{ave}}=0.28`$). Generally one concludes that the larger the inhomogeneity (in amplitude and in length scale), the larger $`T_\mathrm{c}`$ is obtained at a given $`\xi _0/\mathrm{}_{\mathrm{ave}}`$. One notices that the cluster profile with no background scattering ($`b=0`$) differs qualitatively from the other curves in Fig. 4. This profile is exceptional because it has long quasiparticle trajectories where the scattering is negligible.

Let us now turn to the amplitude $`\mathrm{\Delta }`$ of the order parameter. In the HSM $`\mathrm{\Delta }/k_\mathrm{B}T_\mathrm{c}`$ depends mostly on the relative temperature $`T/T_\mathrm{c}`$ but also on the scattering rate $`\xi _0/\mathrm{}`$, the phase shifts $`\delta ^{(l)}`$, and the phase (A or B).PhysicaB  The dotted lines in Fig. 5 represent four different values of $`\xi _0/\mathrm{}`$ for B-phase with $`\mathrm{sin}^2\delta ^{(0)}=\frac{1}{2}`$.

In the IISM the parameter $`\xi _0/\mathrm{}`$ is replaced by $`\xi _0/\mathrm{}_{\mathrm{ave}}`$ and $`n(r)`$. The order parameter depends on the location $`𝐫`$ (8). In order to describe its temperature dependence we use parameters $`\mathrm{\Delta }_{\mathrm{max}}`$ and $`\mathrm{\Delta }_{\mathrm{min}}`$, and $`\mathrm{\Delta }_{\mathrm{ave}}`$ (9). The temperature dependence of the three characteristic numbers is shown by solid lines in Fig. 5.

We see that there can be strong variation since $`\mathrm{\Delta }_{\mathrm{min}}^2`$ and $`\mathrm{\Delta }_{\mathrm{max}}^2`$ are quite different. The average $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ is considerably below the HSM curve that gives the same $`T_\mathrm{c}`$. Therefore $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ is more suppressed than the critical temperature $`T_\mathrm{c}`$.

The temperature dependence of $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ varies. Let us first consider the case in the right-hand panel of Fig. 5, where the shape of $`\mathrm{\Delta }_{\mathrm{ave}}^2(T)`$ is concave near $`T_\mathrm{c}`$. This is typical for a small inhomogeneity amplitude. In our case it also means that the average scattering is small since the minimum scattering vanishes. We see that at low temperatures $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ agrees very well with the HSM curve that is calculated with the same average scattering rate. At higher temperatures $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ deviates from this since the true $`T_\mathrm{c}`$ is higher than the one given by the average scattering rate.

In the case of the left hand panel of Fig. 5, the temperature dependence of $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ is quite linear near $`T_\mathrm{c}`$. In fact, the linear range is wider (relative to $`T_\mathrm{c}`$) than in the pure case or HSM. This happens when the amplitude of the inhomogeneity is large. Large inhomogeneity means large average scattering, thus the HSM results based on the average scattering rate are more suppressed than in the case of concave $`\mathrm{\Delta }_{\mathrm{ave}}^2(T)`$. In the particular case of Fig. 5(a), the HSM result is quite small since $`\xi _0/\mathrm{}=0.275`$ is near the critical value of complete suppression. In this case $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ is nowhere nearly approximated by the HSM result. However, the linear $`\mathrm{\Delta }_{\mathrm{ave}}^2(T)`$ can also appear in cases where the HSM with average scattering rate is not that suppressed, and provides a good approximation for $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ at low temperatures. This takes place at small $`R/\xi _0`$, where the proximity coupling between different regions tends to average out the inhomogeneity. In this limit the range of scattering rates $`\xi _0/\mathrm{}_{\mathrm{ave}}`$ where $`\mathrm{\Delta }_{\mathrm{ave}}^2(T)`$ is concave near $`T_\mathrm{c}`$ seems to vanish.

The concave shape of $`\mathrm{\Delta }_{\mathrm{ave}}^2(T)`$ in Fig. 5(b) can be understood so that the different regions of an inhomogeneous sample have transitions more or less independently of each other: At $`T_\mathrm{c}`$ only $`\mathrm{\Delta }(r=0)=\mathrm{\Delta }_{\mathrm{max}}`$ starts to grow but $`\mathrm{\Delta }_a(r=R)=\mathrm{\Delta }_{\mathrm{min}}`$ and $`\mathrm{\Delta }_r(R)`$ (not shown) remain negligible until they start to grow at a lower temperature. In spite of the inhomogeneity, the onset of superfluidity indicated by $`\mathrm{\Delta }_{\mathrm{max}}`$ is very sharp giving a well defined $`T_\mathrm{c}`$.

The superfluid density $`\rho _\mathrm{s}`$ is conveniently expressed relative to the density of the liquid, $`\rho _\mathrm{s}/\rho _{\mathrm{tot}}`$. In addition to parameters discussed above, this depends on the Fermi-liquid parameter $`F_1^s`$ in both the HSM and IISM. (In the IISM also other parameters could contribute, but they are neglected here.) Since $`F_1^s10`$ is quite big, it has a strong effect on the results. The HSM results for $`\stackrel{}{\rho }_s`$ in both the A and B phases have been calculated in Refs. PhysicaB,  and higashitani, .

In the IISM we calculate $`\rho _\mathrm{s}`$ for the B phase. Some curves are plotted in Fig. 6.

As above, we compare the IISM with the HSM calculated with the same average scattering rate $`\xi _0/\mathrm{}_{\mathrm{ave}}`$. The results of both models are close to each other at low temperatures when the scattering rate is small. At higher temperatures $`\rho _\mathrm{s}`$ in IISM is larger because of its higher $`T_\mathrm{c}`$. Thus $`\rho _\mathrm{s}`$ is relatively more suppressed than $`T_\mathrm{c}`$. At larger scattering rates, $`\rho _\mathrm{s}`$ in the HSM becomes smaller and is completely suppressed at $`\xi _0/\mathrm{}=0.28`$, while $`\rho _\mathrm{s}`$ in the IISM stays finite much beyond that. All these results are qualitatively similar as discussed above for $`\mathrm{\Delta }_{\mathrm{ave}}^2`$. The discussion of concave or linear temperature dependence in the case of $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ cannot directly be applied to $`\rho _\mathrm{s}`$ since this behavior is largely masked by the effect of $`F_1^s`$. The superfluid density $`\rho _\mathrm{s}`$ is also different from $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ because of current conservation, which forces quite different current pattern in the void and cluster impurity profiles.

In order to reduce the effect of $`F_1^s`$, we define a bare superfluid density $`\rho _\mathrm{s}^\mathrm{b}`$ by

$`{\displaystyle \frac{\rho _\mathrm{s}^\mathrm{b}}{\rho _{\mathrm{tot}}}}={\displaystyle \frac{\left(1+\frac{1}{3}F_1^s\right)(\rho _\mathrm{s}/\rho _{\mathrm{tot}})}{1+\frac{1}{3}F_1^s(\rho _\mathrm{s}/\rho _{\mathrm{tot}})}}.`$ (18)
In the case of the HSM this is independent of $`F_1^s`$ and in the case of the IISM the dependence is small, only a few percent over the whole pressure range $`0\mathrm{}3.4`$ MPa. We calculate $`\rho _\mathrm{s}^\mathrm{b}`$ at an intermediate pressure of 1.56 MPa where $`F_1^s=10`$.

The order parameter and superfluid density can be expressed compactly by defining suppression factors

$`S_{\mathrm{\Delta }^2}(t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }^2(tT_\mathrm{c})}{\mathrm{\Delta }_0^2(tT_{\mathrm{c0}})}}`$
$`S_{\rho _\mathrm{s}^\mathrm{b}}(t)`$ $`=`$ $`{\displaystyle \frac{\rho _\mathrm{s}^\mathrm{b}(tT_\mathrm{c})}{\rho _{\mathrm{s0}}^\mathrm{b}(tT_{\mathrm{c0}})}}.`$ (19)
Here $`\rho _{\mathrm{s0}}^\mathrm{b}`$ is the bare superfluid density and $`\mathrm{\Delta }_0^2`$ the order parameter, both in pure superfluid. The parameter $`t`$ is the temperature relative to the transition temperature. Additionally, we wish to eliminate the parameter $`\mathrm{}_{\mathrm{ave}}`$, which is not directly measurable. This is achieved in plotting the suppression factors as a function of $`(T_\mathrm{c}/T_{\mathrm{c0}})^2`$. One such plot is shown in Fig. 7.

The important conclusion from Fig. 7 is that the suppression factors for the A and B phases are almost identical, although there is a clear difference in the values of $`\mathrm{\Delta }_{\mathrm{ave}}`$. This generalizes the result found previously within the HSM.PhysicaB  Also, the difference in $`S_{\mathrm{\Delta }^2}`$ between $`\mathrm{\Delta }_{\mathrm{ave}}`$ and $`\mathrm{\Delta }_{\mathrm{NMR}}`$ is small. The differences increase with increasing inhomogeneity, so that the suppression factors for $`\mathrm{\Delta }_{\mathrm{ave}}^2`$ and $`\mathrm{\Delta }_{\mathrm{NMR}}^2`$ differ by $`5`$% in our extreme case $`R=2\mathrm{}_{\mathrm{ave}}`$. All these differences are rather small and, in order to simplify the plots, we present below $`S_{\mathrm{\Delta }^2}`$ only for the B phase $`\mathrm{\Delta }_{\mathrm{ave}}`$.

Suppression factors for both $`\mathrm{\Delta }^2`$ and $`\rho _\mathrm{s}^\mathrm{b}`$ as functions of $`(T_\mathrm{c}/T_{\mathrm{c0}})^2`$ are plotted in Fig. 8.

The HSM results are shown for comparison. They are all above the diagonal. It is clearly visible that the effect of inhomogeneity is to bend these curves down. This means that $`\rho _\mathrm{s}^\mathrm{b}`$ and $`\mathrm{\Delta }^2`$ are more strongly suppressed than $`T_\mathrm{c}`$. With a steep profile and a large $`R/\mathrm{}_{\mathrm{ave}}`$ one obtains a strong suppression and a wide flat region in the suppression curve. Superfluid density is clearly more sensitive to the inhomogeneity than the order parameter.

## VII Comparison with experiments

A comparison of the calculated transition temperatures with experimental results is given in Fig. 9.

The pressure-independent $`\mathrm{}_{\mathrm{ave}}`$ is determined so that the theoretical curves and the experimental data intersect at the point where $`T_\mathrm{c}/T_{\mathrm{c0}}=0.7`$. For the samples measured in Manchester \[Fig. 9b\] the HSM gives a fairly good fit but a better correspondence is obtained with IISM by using $`R=\mathrm{}_{\mathrm{ave}}`$. For samples measured in Cornell and Northwestern \[Fig. 9a\] the closest curve is the IISM with $`R=2\mathrm{}_{\mathrm{ave}}`$ and $`j=8`$ in Eq. (6). This represents roughly the best fit obtained by the IISM.

Particularly interesting are the samples A and C measured in Cornell \[Fig. 9(a)\], which also have been studied using small-angle x-ray scattering. The values measured for the aerogel correlation scale $`\xi _\mathrm{a}`$ are 130 and 84 nm.Porto99  These are slightly smaller but on the same order of magnitude as our closest curve values $`R=210`$ and 140 nm, respectively. This supports the view that the better fit obtained by the IISM as compared to the HSM is not only due to more fitting parameters, but a more realistic modelling of the structure of the aerogel.

Some experimental points for both the order parameter and superfluid density are shown in Fig. 8. Each marker ($`\mathrm{}`$, $`\mathrm{}`$, $``$ etc.) identifies a data set that corresponds to one aerogel sample at a fixed reduced temperature $`t`$. The different points in each data set are obtained from measurements at different pressures. The superfluid density is measured by torsional oscillator and the interpretation of the experiments is rather straightforward. The order parameter is measured by NMR. This is possible because the frequency shift depends on the dipole-dipole energy, Eqs. (12) or (16). For interpretation one has to know whether the phase is A or B. The Stanford point ($`\mathrm{}`$) is based on seeing similar suppression in both phases.Barkerprl  The Northwestern data is drawn here by assuming A phase and $`\widehat{𝐝}`$ and $`\widehat{𝐥}`$ perpendicular to magnetic field, as they are in bulk 3He-A.Sprague95 ; Sprague96  However, it is not clear if this original assumption is correct.halprev  In case of the B phase the analysis would be more complicated because it would require an analysis of the texture.Alles ; Barkerprl ; btex

One can notice that the experimental data for the order parameter and superfluid density are qualitatively similar. They all fall clearly below the diagonal, and appear to extrapolate to $`S=0`$ already at a finite $`T_\mathrm{c}/T_{\mathrm{c0}}`$. This is in clear disagreement with the HSM results that are all above the diagonal. It can be seen that much better agreement is achieved with the IISM. Ideally each experimental data set should fall on some of the theoretical lines. As one can see, this is not quite the case in the IISM. A reasonable over-all fit to suppression factors and critical temperatures is achieved with the void impurity profile (6) with $`j=2`$ and $`R1.5\mathrm{}_{\mathrm{ave}}`$, but the choice depends on the properties one wishes to emphasize. The agreement could be improved by allowing for a pressure-dependent mean free path $`\mathrm{}`$. However, the most obvious reason for the remaining differences is that the scattering profile in the IISM has only one length scale whereas real aerogel has a wide distribution of length scales.

## VIII Conclusions

We have presented the isotropic inhomogeneous scattering model. We claim it is the simplest model of inhomogeneous scattering that is consistent with large-scale isotropy. Unfortunately, the computations needed are much more demanding than in the homogeneous scattering model. The model itself is independent of the pairing symmetry and thus could be used to study the effect of inhomogeneous impurity distributions in superconductors.

When applied to 3He in aerogel, the IISM gives better agreement with experiments than the HSM. We emphasize that this is not solely due to the IISM having more free parameters than the HSM. On the contrary, the inhomogeneity of aerogel is the most natural if not the only way to explain the differences between measurements and the HSM. The fitted parameters of IISM are in reasonable agreement with measurements and simulations on the structure of aerogel. The fit given by IISM is not perfect, though. The main problem with the IISM is that it contains only one length scale whereas real aerogel must have voids of various sizes.

###### Acknowledgements.

We wish to thank V. Dmitriev, W. Halperin, J. Hook, G. Kharadze, D. Osheroff, J. Parpia, D. Rainer, J. Sauls and J. Viljas for useful discussions. We also thank the Academy of Finland for financial support and CSC for providing computer resources.

## Appendix A

Here we present the relevant equations of the quasiclassical theory and discuss their numerical solution. The central quantity is the quasiclassical propagator $`\stackrel{˘}{g}`$.Serene  It is a $`4\times 4`$ matrix, whose components can be represented using the Pauli spin matrices $`\underset{¯}{\sigma }_i`$ as

$`\stackrel{˘}{g}=\left(\begin{array}{cc}g+𝐠\underset{¯}{𝝈}& (f+𝐟\underset{¯}{𝝈})\mathrm{i}\underset{¯}{\sigma }_2\\ \mathrm{i}\underset{¯}{\sigma }_2(\stackrel{~}{f}+\stackrel{~}{𝐟}\underset{¯}{𝝈})& \stackrel{~}{g}\underset{¯}{\sigma }_2\stackrel{~}{𝐠}\underset{¯}{𝝈}\underset{¯}{\sigma }_2\end{array}\right).`$ (22)
In equilibrium the propagator $`\stackrel{˘}{g}(\widehat{𝐤},𝐫,ϵ_m)`$ depends on the direction of momentum $`\widehat{𝐤}`$, on the location $`𝐫`$, and on the Matsubara energy $`ϵ_m=\pi k_\mathrm{B}T(2m1)`$, where $`m`$ is an integer. The propagator is determined by the Eilenberger equations,

$`[\mathrm{i}ϵ_m\stackrel{˘}{\tau }_3\stackrel{˘}{\sigma },\stackrel{˘}{g}]+\mathrm{i}\mathrm{}v_\mathrm{F}\widehat{𝐤}_𝐫\stackrel{˘}{g}`$ $`=`$ $`0,`$ (23)
$`\stackrel{˘}{g}\stackrel{˘}{g}`$ $`=`$ $`1.`$ (24)
Here $`[A,B]=ABBA`$ denotes a commutator and $`\stackrel{˘}{\tau }_i`$ denote Pauli matrices in the Nambu space; $`\stackrel{˘}{\tau }_i=\underset{¯}{\sigma }_i1`$. The self energy matrix $`\stackrel{˘}{\sigma }=\stackrel{˘}{\sigma }_{\mathrm{mf}}+\stackrel{˘}{\sigma }_{\mathrm{imp}}`$ consists of mean field and impurity contributions. For spin-triplet pairing the former part has the form

$$\stackrel{˘}{\sigma }_{\mathrm{mf}}(\widehat{𝐤},𝐫)=\left(\begin{array}{cc}\nu (\widehat{𝐤},𝐫)+𝝂(\widehat{𝐤},𝐫)\underset{¯}{𝝈}& 𝚫(\widehat{𝐤},𝐫)\underset{¯}{𝝈}\mathrm{i}\underset{¯}{\sigma }_2\\ \mathrm{i}\underset{¯}{\sigma }_2𝚫^{}(\widehat{𝐤},𝐫)\underset{¯}{𝝈}& \nu (\widehat{𝐤},𝐫)\underset{¯}{\sigma }_2𝝂(\widehat{𝐤},𝐫)\underset{¯}{𝝈}\underset{¯}{\sigma }_2\end{array}\right).$$
(25)
The real-valued functions $`\nu `$ and $`𝝂`$ depend on the Fermi liquid parameters $`F_l^{s,a}`$. Omitting $`F_3^s`$, $`F_1^a`$ and higher coefficients we have $`𝝂=0`$ and

$`\nu (\widehat{𝐤},𝐫)`$ $`=`$ $`{\displaystyle \frac{\pi k_\mathrm{B}TF_1^s}{1+\frac{1}{3}F_1^s}}{\displaystyle \underset{ϵ_m}{}}(\widehat{𝐤}\widehat{𝐤}^{})g(\widehat{𝐤}^{},𝐫,ϵ_m)_{\widehat{𝐤}^{}},`$ (26)
where $`\mathrm{}_{\widehat{𝐤}}`$ denotes angular average. The off-diagonal part of $`\stackrel{˘}{\sigma }_{\mathrm{mf}}`$ defines the order parameter, which for $`p`$-wave pairing is $`𝚫=\stackrel{}{A}\widehat{𝐤}`$. It is determined by the weak-coupling self-consistency equation

$`𝚫(\widehat{𝐤},𝐫)\mathrm{ln}{\displaystyle \frac{T}{T_{\mathrm{c0}}}}+\pi k_\mathrm{B}T{\displaystyle \underset{ϵ_m}{}}\left[{\displaystyle \frac{𝚫(\widehat{𝐤},𝐫)}{|ϵ_m|}}3(\widehat{𝐤}\widehat{𝐤}^{})𝐟(\widehat{𝐤}^{},𝐫,ϵ_m)_{\widehat{𝐤}^{}}\right]=0.`$ (27)
The impurity self energy $`\stackrel{˘}{\sigma }_{\mathrm{imp}}`$ is related to the forward-scattering $`t`$ matrix via

$$\stackrel{˘}{\sigma }_{\mathrm{imp}}(\widehat{𝐤},𝐫,ϵ_m)=n(𝐫)\stackrel{˘}{t}(\widehat{𝐤},\widehat{𝐤},𝐫,ϵ_m),$$
(28)
where $`n(𝐫)`$ is the impurity density. The usual derivation of this formula assumes a uniform impurity density $`n=\mathrm{const}`$.AGD  However, it is valid for an arbitrary function $`n(𝐫)`$ as long as it allows the impurity positions to be uncertain on the scale of the Fermi wave length $`\lambda _\mathrm{F}`$. We assume that only $`s`$-wave scattering is important. This allows us to write the $`t`$-matrix equation directly for the impurity self-energy

$$\stackrel{˘}{\sigma }_{\mathrm{imp}}(𝐫,ϵ_m)=n(𝐫)v\stackrel{˘}{1}+\pi N_\mathrm{F}v\stackrel{˘}{g}(\widehat{𝐤},𝐫,ϵ_m)_{\widehat{𝐤}}\stackrel{˘}{\sigma }_{\mathrm{imp}}(𝐫,ϵ_m).$$
(29)
Here $`2N_\mathrm{F}=m^{}k_\mathrm{F}/\pi ^2\mathrm{}^2`$ is the total density of states at the Fermi surface. The effective mass $`m^{}`$ is related to the atomic mass $`m_3`$ by $`m^{}=(1+\frac{1}{3}F_1^s)m_3`$. The scattering potential $`v`$ is related to the scattering phase shift $`\delta ^{(0)}`$ by $`v=(1/\pi N_\mathrm{F})\mathrm{tan}\delta ^{(0)}`$.

In the calculation we use the following symmetry of the propagator and self energy

$$[\stackrel{˘}{u}(\widehat{𝐤},𝐫,ϵ_m)]^T=\stackrel{˘}{\tau }_2\stackrel{˘}{u}(\widehat{𝐤},𝐫,ϵ_m)\stackrel{˘}{\tau }_2.$$
(30)
This limits the nonzero components of $`\stackrel{˘}{g}`$ to 10, and allows one to solve the Eilenberger equation with the “multiplication trick” of Ref. Zhang,  using only five components. In the current-free case all equations are immediately compatible with the symmetry in Eq. (30). In the case of nonzero current a possible problem could arise from the fact that in iterating the impurity equation (29) the product $`\stackrel{˘}{g}_{\widehat{𝐤}}\stackrel{˘}{g}_{\widehat{𝐤}}`$ might not be proportional to a unit matrix $`\stackrel{˘}{1}`$. The inconvenient terms always appear with odd powers of $`v`$. Considering that the phase shifts in aerogel are random, such terms average out and we can use the symmetry of Eq. (30) also for nonzero current.

By taking into account all the symmetries the impurity self energy in the IISM has the form

$`\stackrel{˘}{\sigma }_{\mathrm{imp}}(𝐫,ϵ_m)=\left(\begin{array}{cc}\mathrm{i}\delta +\mathrm{i}𝜹\underset{¯}{𝝈}& 𝐑\underset{¯}{𝝈}\mathrm{i}\underset{¯}{\sigma }_2\\ \mathrm{i}\underset{¯}{\sigma }_2𝐑^{}\underset{¯}{𝝈}& \mathrm{i}\delta +\underset{¯}{\sigma }_2\mathrm{i}𝜹\underset{¯}{𝝈}\underset{¯}{\sigma }_2\end{array}\right)+w\stackrel{˘}{1},`$ (33)
where $`𝐑`$, $`𝜹`$, and $`\delta `$ can be expressed analytically as functions of the propagator $`\stackrel{˘}{g}_{\widehat{𝐤}}`$.diplomit  The term $`w`$ is of no interest since it drops out in the Eilenberger equation (23). In cylindrical coordinates ($`\rho `$, $`\varphi `$, $`z`$) the only $`\varphi `$ dependence appears in the vector $`𝐑`$ in the $`A`$ phase, where $`𝐑=R_0(\rho ,z)e^{\mathrm{i}\varphi }\widehat{𝐳}`$ and $`R_0`$ is real. The terms $`\delta `$ and $`𝜹`$ are always real. In the A phase $`𝜹=0`$ and in the $`B`$ phase $`𝜹=\delta _\varphi \widehat{\mathit{\varphi }}`$ and $`𝐑=R_\rho \widehat{𝝆}+R_z\widehat{𝐳}`$.

From the propagator $`\stackrel{˘}{g}`$ (22) one can calculate the supercurrent density $`𝐣_\mathrm{s}`$ using

$$𝐣_\mathrm{s}(𝐫)=2\pi m_3v_\mathrm{F}N_\mathrm{F}k_\mathrm{B}T\underset{ϵ_m}{}\widehat{𝐤}g(\widehat{𝐤},𝐫,ϵ_m)_{\widehat{𝐤}}.$$
(34)
The superfluid density $`\rho _s`$ is obtained by calculating the averaged current at small values of $`𝐯_\mathrm{s}`$.

The first step in the calculation is to give initial values for the fields $`\stackrel{}{A}`$, $`\stackrel{˘}{\sigma }_{\mathrm{imp}}`$, and $`\nu `$. Normally we used values obtained from the HSM. The iteration of these fields was started by transforming the data from the cylindrical/spherical coordinates to “trajectory coordinates.” These are Cartesian coordinates with one coordinate, $`u`$, along the trajectory and two others ($`b`$ and $`t`$) specifying the position of the trajectory in the unit cell. Next the Eilenberger equation was solved along the trajectories. This was done by first calculating an unphysical exponentially growing propagator along each trajectory using the fourth order Runge-Kutta method. The end values at the cell boundary were saved to be used as starting values for next step in iteration. An exponentially decreasing solution was deduced using symmetries. The bounded physical propagator was obtained as a commutator of these two unphysical solutions as explained in Ref. Zhang, . The solution of the Eilenberger equation was repeated for all Matsubara frequencies $`|ϵ_m|ϵ_N`$. This was repeated for all trajectories. In the general case it represents loops for coordinates $`b`$, $`t`$ and the angle $`\alpha `$ between the trajectory direction and the $`z`$-axis. For the B phase without a current only the impact parameter $`b`$ with respect to the center of the sphere was actually needed because of symmetry.

In the next step better estimates for $`\stackrel{}{A}`$, $`\stackrel{˘}{\sigma }_{\mathrm{imp}}`$, and $`\nu `$ were calculated. For that the propagator was converted from trajectory coordinates back to cylindrical/spherical coordinates and the required angular averages were calculated. The contribution from higher Matsubara frequencies $`|ϵ_m|>ϵ_N`$ was approximated by a Ginzburg-Landau form ErkkiGL  with Matsubara sums evaluated using the Euler-MacLaurin formula. Using the updated values of the fields, the process was started from the beginning. The boundary condition (4) was used to calculate the value of the exponential propagator at the initial point of the trajectory from the value at the final point stored on the previous round. This loop was repeated until the fields converged.

In the numerical algorithm the number of discretized points in cylindrical coordinates $`\rho _i`$ and $`z_j`$ was approximately 40 for both in the range from $`0`$ to $`R`$. The discretization step in the Cartesian trajectory coordinates was approximately the same. Simple interpolation formulas were used for the transformations between coordinate systems. The angle $`\alpha `$ was typically discretized by eight Gaussian points in the range $`0<\alpha <\pi /2`$. The number of positive Matsubara frequencies $`ϵ_m`$ used was typically less than 20 for temperatures above $`0.1T_{\mathrm{c0}}`$. We made also more accurate test calculations. The qualitative behavior remains the same, but there are inaccuracies on the order of 2% in the results presented here.