text:
# Equilibrium Simulations of 2D Weak Links in 𝑝-wave Superfluids

## I INTRODUCTION

Recent experimental observationsberkeley ; paris  of the Josephson effect in weak links of superfluid 3He left theorists with some interesting problems. One of these was related to the interpretation of the “$`\pi `$ states”, local minima of the Josephson energy at phase differences other than $`0`$ (mod $`2\pi `$). Inspired by these findings, some theoretical work was done.viljas ; yippi ; viljas2 ; nishida ; zhang ; hatakenaka  By now there exists a reasonably good understanding of how such $`\pi `$ states could be obtained in small “pinhole” contacts,yippi ; viljas2 ; nishida  in arrays of apertures via the anisotextural effect,viljas ; viljas2  or in single large apertures in the Ginzburg-Landau (GL) regime.viljas ; thuneberg88  However, large physically relevant parts of the parameter space remain unexplored. In this paper we study the GL regime in more detail, and attempt to bridge the gap between the pinhole and the large-aperture limits.

In Ref. viljas,  we studied a three-dimensional (3D) circular aperture between two bulk volumes of B phase 3He by solving the GL equations in a full 3D lattice. With this calculation we were able to demonstrate that a local minimum of energy exists at phase differences close to $`\pi `$, which is associated with a separate “$`\pi `$ branch” in the current-phase relation. This solution was never clearly found to be a global minimum of energy, and no jumps to it from the strongly hysteretic “$`0`$ branches” could be observed. Therefore it appeared to be experimentally inaccessible. This calculation was very restricted, however, since the bulk order parameters were fixed to be equal on the two sides. In the present calculation we are able to impose much more general boundary conditions on the soft degrees of freedom of the order parameters, and to use the A phase in addition to the B phase. On the other hand, we have not attempted to carry out the calculations in 3D, but rather we use a two-dimensional (2D) model of a long slit-shaped weak link in a planar wall.thuneberg88  Nevertheless, as our experience obtained with the previous calculation of Ref. viljas,  suggests, we may expect very similar effects to exist in all weak links, regardless of their shape. In the A phase this is not quite so, as will be discussed below, because the critical current can even vanish in some very symmetrical situations.

Since this article is to appear in the proceedings of a winter school, our approach is in many ways tutorial. Furthermore, since the involved work is largely numerical, special attention is paid to explaining the computational methods. Throughout, we shall only deal with thermodynamic equilibrium, and the numerics are thus mostly related to a minimization of the GL free-energy functional on a square lattice. Our relaxation methods are actually quite standard ones, but have never been presented in detail in the specific context of 3He.

In Section 2 we discuss the Josephson coupling on a formal level, and remind how non-sinusoidal energy-phase relations can be obtained via a perturbative approach in 3He and unconventional superconductors alike. From Section 3 onward we turn to the case of 3He, and first review some basic issues and notations. We relate the symmetries of the free energy functional to the conserved supercurrents, and review the GL theory briefly. In particular, we consider spin supercurrents on an equal footing with the mass supercurrents everywhere. Section 4 introduces our 2D weak link model, with the divisions into numerical and asymptotic regions, and discusses physical ways to control the boundary conditions. Finally, an analysis of the four different phase combinations is given in Section 5, and some examples of current-phase relations (CPR’s) are presented. All critical currents and phase diagrams are summarized in Section 6. In Section 7 we end with some conclusions and discussion. The asymptotic solutions far from the weak link are developed in Appendix A, and an analysis of the numerical method is presented in Appendix B.

## II DC JOSEPHSON EFFECT REVISITED

The superfluid state of a BCS superfluid/superconductor is described by a gap matrix, or a pairing potential. This $`2\times 2`$ matrix has the form

$$\mathrm{\Delta }_{\alpha \beta }(\widehat{𝐤})=[\mathrm{\Delta }_s\mathrm{i}\sigma _y+𝚫𝝈\mathrm{i}\sigma _y]_{\alpha \beta }.$$
(1)
Here $`\mathrm{\Delta }_s`$ and $`𝚫`$ are the singlet and triplet pairing amplitudes, respectively, and their argument $`\widehat{𝐤}`$ parametrizes points on the Fermi surface. From the Pauli principle $`\mathrm{\Delta }_{\alpha \beta }(\widehat{𝐤})=\mathrm{\Delta }_{\beta \alpha }(\widehat{𝐤})`$ it follows that $`\mathrm{\Delta }_s`$ and $`𝚫`$ have even parity ($`s,d,g,\mathrm{}`$) and odd parity ($`p,f,h,\mathrm{}`$) in $`\widehat{𝐤}`$, respectively. The gap matrix also serves as the order parameter of the superfluid.

The dc (or *equilibrium*) Josephson effect results from a coupling of two superfluids (left, $`L`$, and right, $`R`$) through a weak link. Due to this coupling, the free energy of the system is generally changed by an amount $`F_J`$ which we call the *Josephson coupling energy*. We assume $`F_J`$ to be a function of the bulk values $`\mathrm{\Delta }^{L,R}`$ of the order parameters. From the order parameters we separate out phase factors by writing $`\mathrm{\Delta }^{L,R}=\mathrm{\Delta }_0^{L,R}\mathrm{exp}(\mathrm{i}\varphi ^{L,R})`$. Because of global gauge invariance $`F_J`$ should independent of the global phase. The dependence of $`F_J`$ on the phase difference $`\mathrm{\Delta }\varphi =\varphi ^R\varphi ^L`$ is called the energy-phase relation (EPR). Since the phase factors are $`2\pi `$ periodic, so is the EPR: $`F_J(\mathrm{\Delta }\varphi +2\pi )=F_J(\mathrm{\Delta }\varphi )`$. As is well known, the derivative of EPR with respect to $`\mathrm{\Delta }\varphi `$ gives the current-phase relation (CPR) for mass supercurrent, $`J_s(\mathrm{\Delta }\varphi )`$.barone  In the triplet case there may also exist a spin current, if $`𝚫^L𝚫^R`$.

In non-magnetic situations it is also reasonable to assume a symmetry of $`F_J(\mathrm{\Delta }\varphi )`$ under time-reversal (TR). If $`\mathrm{\Delta }^{L,R}`$ are both “TR invariant”, meaning that $`\mathrm{\Delta }_0^{L,R}`$ can be chosen real so that TR only complex conjugates the phase factor, then $`F_J(\mathrm{\Delta }\varphi )=F_J(\mathrm{\Delta }\varphi )`$.yipgl  Similarly, the CPR then satisfies $`J_s(\mathrm{\Delta }\varphi )=J_s(\mathrm{\Delta }\varphi )`$. Let us apply these results to the case of sufficiently small junctions, where we can assume $`F_J(\mathrm{\Delta }\varphi )`$ to be a single-valued function of $`\mathrm{\Delta }\varphi `$. Making a Fourier expansion of $`F_J(\mathrm{\Delta }\varphi )`$ and dropping terms based on the TR symmetry we get

$$F_J(\mathrm{\Delta }\varphi )=E_J^{(0)}\underset{n=1}{\overset{\mathrm{}}{}}E_J^{(n)}\mathrm{cos}(n\mathrm{\Delta }\varphi ),$$
(2)
with a similar sine expansion for $`J_s(\mathrm{\Delta }\varphi )`$. However, there exist “chiral” states, which are *not* TR invariant. An example is the A phase of superfluid 3He, to which we return shortly. In such cases the TR symmetry no longer implies $`F_J(\mathrm{\Delta }\varphi )=F_J(\mathrm{\Delta }\varphi )`$ in general, but in special cases more general symmetries of the form $`F_J(\beta \mathrm{\Delta }\varphi )=F_J(\beta +\mathrm{\Delta }\varphi )`$ may still be valid,yipgl  see Section 5.

An alternative expansion is to develop $`F_J`$ in powers of the order parameters $`\mathrm{\Delta }^{L,R}`$ (similar to the Ginzburg-Landau expansion). Comparing this expansion with the Fourier expansion (2), we see that the coefficient $`E_J^{(n)}`$ is proportional to $`\mathrm{\Delta }^{2n}[1+O(\mathrm{\Delta }^2)]`$, where $`\mathrm{\Delta }`$ denotes the amplitude of the order parameter. For small junctions the effective expansion parameter is the amplitude of the order parameter divided by temperature, $`\mathrm{\Delta }/k_BT_c`$. Since this is small at temperatures close to $`T_c`$, the lowest order terms dominate in the series (2). Usually the first order term $`F_J^{(1)}=E_J^{(1)}\mathrm{cos}(\mathrm{\Delta }\varphi )`$ is most important at all $`T`$, except in special cases where its amplitude is suppressed due to symmetry reasons. These symmetries for 3He or unconventional superconductors have been studied in several papers.viljas ; viljas2 ; yipgl ; degennes ; larkin  All other degrees of freedom of the order parameters but $`\mathrm{\Delta }\varphi `$ (and those related to the junction itself) are now embedded in the coefficients $`E_J^{(n)}`$, and these may be used as tuning parameters. For example, in some casespisquid ; sfs  the sign of the (normally positive) $`E_J^{(1)}`$ can be reversed, so that one obtains a “$`\pi `$ junction”, where the only minimum of a sinusoidal EPR is at $`\mathrm{\Delta }\varphi =\pi `$ rather than at $`\mathrm{\Delta }\varphi =0`$. In 3He-B such a trick may be done by controlling the order parameter textures with magnetic fields.yippi

For vanishing $`E_J^{(1)}`$ the higher order terms may still give a finite critical current. In particular, a finite $`E_J^{(2)}`$ term gives a $`\pi `$-periodic EPR due to the $`\mathrm{cos}(2\mathrm{\Delta }\varphi )`$ term. If the suppression of $`E_J^{(1)}`$ is only partial, then some interesting mixtures of the $`2\pi `$ and $`\pi `$ (and shorter) periodic components may be observed (see, for example, Ref. tanaka,  or Ref. viljas2,  and references therein). In particular, it is exactly these higher-order terms of Eq. (2) which result in the “$`\pi `$ states” of the $`p`$-wave (pinhole) junctions discussed in Refs. yippi, ; viljas2, ; nishida, ; zhang, . These $`\pi `$ states become more pronounced for $`TT_c`$, since the higher-order $`E_J^{(n)}`$’s are larger there. They are also the reason for the slanting of the CPR of an $`s`$-wave pinhole contact at low temperatures.kulik  However, in this case there is no way suppress $`E_J^{(1)}`$ in order to single out the smaller $`\pi `$ periodic components.

For large junctions we may expect two kinds of complications to the description based on Eq. (2). First, there is the effect of “kinetic inductance”: as a result of current conservation, there must be a finite phase gradient carrying the current $`J_s`$ into, out of, and within the weak link.likharev  Thus the phase difference $`\mathrm{\Delta }\varphi `$ itself depends on $`J_s`$, which makes $`F_J(\mathrm{\Delta }\varphi )`$ and $`J_s(\mathrm{\Delta }\varphi )`$ hysteretic (multivalued) in general. As a result, the expansion in the cosines, Eq. (2), is no longer valid as such. The hysteresis can be modelled by introduction of the “slanting parameter” and a self-consistent set of equations as in Refs. likharev,  and avenel, . Second, in the presence of a multicomponent order parameter the transmission properties of the weak link (reflected by the signs and absolute values of the coefficients $`E_J^{(n)}`$ above) may depend on the order parameter configuration inside the aperture. This configuration may change as a function of $`\mathrm{\Delta }\varphi `$, which changes the CPR’s.monien ; thuneberg88 ; viljas  Thus the order parameter must be solved self-consistently in and around the hole, with boundary conditions specifying $`\mathrm{\Delta }^{L,R}`$ only somewhere far from the junction. The onset of these large-aperture effects in weak links of superfluid 3He is the main subject of the this paper.

These changes can also be seen in the expansion of $`F_J`$ in the order parameters. Instead of $`\mathrm{\Delta }/k_BT_c`$, the effective expansion parameter in a junction of linear dimension $`D`$ turns out to be $`\mathrm{\Delta }D/\mathrm{}v_F`$, the gap divided by the (ballistic) Thouless energy. Using standard relations the expansion parameter can also be expressed as $`D/\xi _{GL}(T)`$, where $`\xi _{GL}(T)`$ is the Ginzburg-Landau temperature dependent coherence length. We see that the expansion breaks down when $`D/\xi _{GL}(T)1`$. This is in agreement with our results below, which show multivalued EPR for $`D\xi _{GL}(T)`$.

## III SUPERFLUID 3He

Below we shall only consider the triplet $`p`$-wave case where the gap vector $`𝚫`$ can be written as $`\mathrm{\Delta }_\mu (𝐱,\widehat{𝐤})=A_{\mu i}(𝐱)\widehat{k}_i`$ with a proper choice of the spin and orbital basis functions.leggett ; vollhardt  In practice this means superfluid 3He, but a similar analysis may be valid in possible triplet superconductors. In 3He, in the absence of magnetic fields, there are two stable bulk phases, the A and the B phase. These are known to correspond to the ABM state $`A=\mathrm{\Delta }_A\widehat{𝐝}(\widehat{𝐦}+\mathrm{i}\widehat{𝐧})\mathrm{exp}(\mathrm{i}\varphi )`$ and the BW state $`A=\mathrm{\Delta }_BR\mathrm{exp}(\mathrm{i}\varphi )`$, respectively. Here $`\mathrm{\Delta }_{A,B}`$ are the bulk gaps of the A and B phases. In the B phase $`R(\widehat{𝝎},\theta )`$ is a rotation matrix, which may be parametrized by the rotation axis $`\widehat{𝝎}`$ and rotation angle $`\theta `$. In the A phase $`\widehat{𝐦}\widehat{𝐧}`$, and one usually defines a third unit vector $`\widehat{𝐥}=\widehat{𝐦}\times \widehat{𝐧}`$ which gives the direction of relative orbital angular momentum of all Cooper pairs. Note that reversing the sign of $`\widehat{𝐝}`$ is equivalent to a phase shift by $`\pi `$, and that any phase shift by $`\varphi `$ is equivalent to a rotation of $`\widehat{𝐦},\widehat{𝐧}`$ around $`\widehat{𝐥}`$ by angle $`\varphi `$. The separation of the phase factor is therefore not unique, and it is further complicated by the existence of textures in the $`\widehat{𝐥}`$ field. However, in what follows, we can assume $`\widehat{𝐥}`$ to be constant most of the time, and if the same definitions are used consistently, no problems should arise. In the B phase the gap $`|𝚫(\widehat{𝐤})|`$ is isotropic, whereas in the A phase it has point nodes in the direction of the vector $`\widehat{𝐥}`$. To maximize the condensation energy, $`\widehat{𝐥}`$ is therefore always rigidly oriented perpendicular to solid surfaces. The presence of a specific direction $`\widehat{𝐥}`$ for the Cooper pair orbital angular momentum means that the A phase is not time-reversal invariant. Consequently, in the context of weak links, some of the usually “obvious” symmetries must be reconsidered when A phase is involved. These symmetry properties are discussed below, in Section 5.

The ABM and BW states are well-defined only in the hydrodynamic regime, *i.e.*, on large length scales. Close to surfaces, for example, the order parameter will be modified on the coherence length scale $`\xi _0`$ due to scattering of quasiparticles. A weak link involves length scales on the order of $`\xi _0`$ and a locally suppressed order parameter by definition, and therefore a more general treatment is required. In what follows, we let all the components of the order parameter vary freely close to the weak link, and the ABM or BW forms are fixed only on boundaries in the bulk liquid.

### III.1 Symmetries and Conservation Laws

The free energy of a $`p`$-wave spin triplet Fermi superfluid can be expressed as a functional of the order parameter field $`A(𝐱)`$. In a region $`\mathrm{\Omega }`$ bounded by $`\mathrm{\Omega }`$, we assume the free-energy expression to be of the form

$$F_\mathrm{\Omega }[A]=_\mathrm{\Omega }\mathrm{d}^3xf(A,\mathbf{}A).$$
(3)
Neglecting in this functional any terms that couple the spin and orbital degrees of freedom or introduce preferred spin directions (nuclear dipole-dipole and dipole-field interactions), it must remain invariant under global gauge transformations \[$`U(1)`$\] and global spin rotations \[$`SO(3)^s`$\] of $`A`$. These are parametrized by the independent real parameters $`\varphi `$ and $`\theta _\alpha ,\alpha =x,y,z`$, respectively, which are independent of $`𝐱`$. Infinitesimally the transformations are written as

$`A_{\mu i}`$ $`A_{\mu i}^{}=e^{\mathrm{i}\delta \varphi }A_{\mu i}A_{\mu i}+\mathrm{i}A_{\mu i}\delta \varphi `$ (4)
$`A_{\mu i}`$ $`A_{\mu i}^{}=R_{\mu \nu }(\delta \theta _\alpha )A_{\nu i}A_{\mu i}ϵ_{\alpha \mu \nu }A_{\nu i}\delta \theta _\alpha ,`$ (5)
where $`R(\theta _\alpha )`$ is a (right-handed) rotation matrix. In this context Eq. (5) should not be considered as a *passive* rotation of the spin-coordinate system, but rather as an *active* transformation of the physical state. Suppose that originally the order parameter corresponds to a stationary point of $`F_\mathrm{\Omega }`$, *i.e.*, it is a solution of the Euler-Lagrange equations $`\delta F_\mathrm{\Omega }/\delta A=\delta F_\mathrm{\Omega }/\delta A^{}=0`$ with fixed boundary conditions $`\delta A|_\mathrm{\Omega }=0`$, for example. Then for variations $`\delta A`$ which leave $`f`$ invariant, one finds a conservation law $`\mathbf{}𝐣_N=0`$ for some generalized current $`𝐣_N`$. (This is a special case of *Noether’s theorem*.byron ) In particular, for the variations of Eqs. (4) and (5) the corresponding currents are the “mass supercurrent” and the three independent components of “spin supercurrent”

$`𝐣_s`$ $`={\displaystyle \frac{2m_3}{\mathrm{}}}(+\mathrm{i}A_{\mu i}{\displaystyle \frac{f}{\mathbf{}A_{\mu i}}}+c.c.)`$ (6)
$`𝐣_\alpha ^{\mathrm{spin}}`$ $`=+ϵ_{\alpha \mu \nu }A_{\nu i}{\displaystyle \frac{f}{\mathbf{}A_{\mu i}}}+c.c.,\alpha =x,y,z.`$ (7)
Orbital rotations do not in general keep the energy invariant and therefore no conserved “orbital supercurrent” exists. The physical interpretation of Eqs. (6) and (7) as “mass” and “spin” currents (and thus their normalization constants) cannot be seen from this phenomenological approach, but they can be verified from microscopic theory. Note that the same current expressions can be obtained by inserting into Eq. (3) the “gauge invariant derivative” prescriptions $`\delta _{\mu \nu }_i\delta _{\mu \nu }_i`$ $`+\delta _{\mu \nu }\mathrm{i}a_iϵ_{\mu \nu \alpha }b_{\alpha i}`$ and expanding to linear order in $`𝐚`$ and $`b_{\alpha i}`$.

The conservation laws $`\mathbf{}𝐣_s=0`$ and $`\mathbf{}𝐣_\alpha ^{\mathrm{spin}}=0`$ play an important role in the weak link problem since the current density distributions are closely related to the perturbations of the order parameter at the junction. Due to the conservation of currents, the perturbations have, in principle, infinite range. In reality the symmetry-breaking dipole interaction makes spin currents decay on the dipole length scale $`\xi _D`$. This is much larger than the scale considered in this paper and thus we can assume the conservation of both currents to be exact.

### III.2 Ginzburg-Landau Theory

In numerical calculations we use the Ginzburg-Landau (GL) expansion of $`F_\mathrm{\Omega }`$, which is valid at temperatures close to the transition temperature $`T_c`$, where the amplitude of $`A`$ is small. The GL theory and its applications have been thoroughly discussed in various articles.degennes ; leggett ; buchholtzfetter ; thuneberg87  Taking into account the bulk and gradient terms only, the GL energy density is

$$\begin{array}{cc}\hfill f(A,& A)=\hfill \\ & \alpha Tr(AA^\mathrm{T})+\beta _1|Tr(AA^\mathrm{T})|^2\hfill \\ & +\beta _2[Tr(AA^\mathrm{T})]^2+\beta _3Tr(AA^\mathrm{T}A^{}A^\mathrm{T})\hfill \\ & +\beta _4Tr(AA^\mathrm{T}AA^\mathrm{T})+\beta _5Tr(AA^\mathrm{T}A^{}A^\mathrm{T})\hfill \\ & +K_1_iA_{\mu i}^{}_jA_{\mu j}+K_2_iA_{\mu j}^{}_iA_{\mu j}\hfill \\ & +K_3_iA_{\mu j}^{}_jA_{\mu i}.\hfill \end{array}$$
(8)
This includes all the linearly independent terms up to fourth order (second order in gradients) which are invariant under global rotations of spin and orbital coordinates \[$`SO(3)^{s,o}`$\], global gauge transformations \[$`U(1)`$\], and under time reversal, *i.e.*, complex conjugation. These include the transformations of Eqs. (4) and (5). Other terms resulting from possible magnetic dipole-dipole and dipole-field interactions *etc.* could also be included, and they would not introduce major complications to the numerical calculation. However, these interactions only become important on length scales $`\xi _D10\mu `$m and $`\xi _H`$ which we assume to be much larger than the temperature-dependent coherence length $`\xi _{GL}(T)`$ as defined shortly. This restricts the validity of our results from above on the field and temperature scales, since $`\xi _HH^1`$ and $`\xi _{GL}(T)(1T/T_c)^{1/2}`$. In practice it may not be possible do accurate measurements so close to the critical temperature that the latter restriction would be a problem. Using the GL energy of Eq. (8) in Eqs. (6) and (7), the formulas for the currents given in Refs. buchholtzfetter,  and thuneberg87, , for example, are exactly reproduced.

The GL free-energy expansion was introduced above phenomenologically, with several parameters: $`\alpha ,\beta _i,K_i`$ and $`\gamma `$. Values for these can be calculated from quasiclassical theory.sr ; thuneberg87  All temperature-dependence of GL theory is in the coefficient $`\alpha (T)=N(0)(1T/T_c)/3`$, where $`N(0)=m^{}k_F/2\pi ^2\mathrm{}^2`$ is the normal-state density of states on the Fermi surface for one spin species. The gradient-energy parameters are $`\gamma =3`$, and $`K_1/(\gamma 2)=K_2=K_3=K`$$``$$`(7\zeta (3)/240\pi ^2)`$$`N(0)(\mathrm{}v_F/k_BT_c)^2`$. These are all weak-coupling (WC) results, but the strong-coupling (SC) corrections will not be used, since they are not very accurately known, and are probably small. Also, $`K`$ and $`\alpha `$ appear only in the natural length scale of GL theory, which we use as our unit of length. This is the temperature-dependent coherence length $`\xi _{GL}(T)=\sqrt{K/\alpha }=\xi (0)(1T/T_c)^{1/2}`$, where $`\xi (0)=\sqrt{21\zeta (3)/240\pi ^2}(\mathrm{}v_F/k_BT_c)`$ is one way to define the zero-temperature coherence length $`\xi _0`$. The WC values for the $`\beta `$ parameters are $`\beta _i^{WC}/\beta _{BCS}=(1,2,2,2,2)`$, for $`i=1,\mathrm{},5`$, where $`\beta _{BCS}(7\zeta (3)/240\pi ^2)N(0)/(k_BT_c)^2`$. Pressure-dependent SC corrections to these have been calculated in Ref. sauls, . The main effect of the SC corrections is to change the difference $`\mathrm{\Delta }f_{AB}=f_A^cf_B^c`$ of the A and B phase condensation energies, $`f_A^c=2\alpha \mathrm{\Delta }_A^2/2`$ and $`f_B^c=3\alpha \mathrm{\Delta }_B^2/2`$, so that the A phase can become stable at pressures above the polycritical one, $`p_0`$. The value of $`p_0`$ depends sensitively on the $`\beta `$ parameters, and our fit gives it at roughly $`p_0=28.7`$ bar, whereas experimentally $`p_021`$ bar. Pressures close to $`p_0`$ are important for studying weak links between the A and the B phases, since the phase boundary must remain pinned in the aperture. For a slit of width $`W`$ the condition for the stability of the boundary is $`|\mathrm{\Delta }f_{AB}|<2\sigma _{AB}/W`$, where $`\sigma _{AB}`$ is the surface tension of the A-B interface.viljasab

The boundary conditions for the order parameter on solid surfaces generally satisfy $`A_{\mu i}\widehat{s}_i=0`$, where $`\widehat{𝐬}`$ is the surface normal.degennes  Everywhere in this paper we use a more strict boundary condition $`A_{\mu i}=0`$. This corresponds to a microscopically rough surface, which scatters quasiparticles in a completely diffusive way. Most realistic surfaces are suspected to be of this type and, in addition, this is the simplest one to implement numerically.

## IV THE WEAK LINK PROBLEM

We now apply the above considerations to describe a weak link of the form shown in Fig. 1. It is a 2D approximation of the slit-type weak links used in experiments like those described in Refs. paris,  and avenel, . This “archetype” is here considered between two infinite volumes, but we could also insert it between two flow channels or other restricted 2D geometries.ullahfetter  The origin of coordinates is placed in the middle of the junction, with the $`x`$ axis running through the aperture. Throughout, we call the sides with negative and positive $`x`$ coordinate the left ($`L`$) and right ($`R`$) sides, respectively.

The region $`\mathrm{\Omega }`$ in Eq. (3) is now the one inside the outer circular arcs of radii $`R_{\mathrm{}}`$. In 3D the $`R_{\mathrm{}}`$ cutoff could be taken to infinity. In 2D an inconvenient finite value for $`R_{\mathrm{}}`$ is required due to the logarithmic $`\mathrm{ln}(R_{\mathrm{}}/r)`$ (and not $`1/r`$) asymptotic behavior of the phase fields \[see Eqs. (24) and (32) in Appendix A\]. This cutoff may be thought to describe some effective length scale, determined by the dimensions of the container, or a distance after which the decay turns three-dimensional ($`1/r`$), for example. The $`z`$ directional length of the assumed slit gives a good approximation for the latter.thuneberg88  In the A phase a more natural choice for the form of the cutoff arcs is actually ellipsoidal, with the perpendicular coordinate scaled down by $`\sqrt{(\gamma +1)/2}\sqrt{2}`$. This follows from the smaller superfluid density in the direction of the orbital angular momentum vector $`\widehat{𝐥}`$ — see Appendix A.

### IV.1 Coupling Energy and Soft Variables

The left and right boundary conditions are in general functions $`A_{\mu i}^L(s_L)`$, $`A_{\mu i}^R(s_R)`$ of $`s_{L,R}`$ which parametrize positions on the $`\mathrm{\Gamma }^{L,R}`$. The Josephson coupling energy is a functional $`F_J[A^L,A^R]`$$`=\mathrm{min}_AF_\mathrm{\Omega }[A]F_0`$. Here we explicitly subtract a term $`F_0`$, which is the energy in $`\mathrm{\Omega }`$ for the same $`A^{L,R}`$, but with the aperture closed. This is done to remove the bulk and surface contributions, but in 2D $`F_J`$ still depends on $`R_{\mathrm{}}`$ due to the currents, as explained in Appendix A. Here we have assumed that fixed boundary conditions are suitable, and that they define the minimum uniquely. Actually, there can be several local minima which are separated from each other by energy barriers. Therefore the current-phase relation may be multivalued, each CPR branch corresponding to a different type of minimum. Jumps (“phase slips”) between these branches are hysteretic.

As special cases we shall consider the A and B phases of 3He. There are then four main cases which can realized for a hole in a nonmagnetic planar wall. Using pairs of letters to denote the $`L`$ and $`R`$ phases, respectively, these are BB, AA with parallel ($``$), AA with antiparallel ($``$) $`\widehat{𝐥}`$’s, and the different configurations of AB. The asymptotic forms in all of these cases can be treated very similarly by writing $`A_{\mu i}=R_{\mu \nu }A_{\nu i}^{(0)}(x)e^{\mathrm{i}\varphi }`$, where the broken-symmetry variables, or “soft degrees of freedom” $`R(\widehat{𝝎},\theta )`$ and $`\varphi `$ are assumed to be constants on $`\mathrm{\Gamma }^{L,R}`$.

The function $`A^{(0)}(x)`$ is the order parameter calculated for a planar diffusive wall in the absence of a weak link, and thus includes the suppression at the wall. A one-dimensional minimization is required to determine it (see Appendix B). In the B phase $`A^{(0)}=\mathrm{diag}(\mathrm{\Delta }_{},\mathrm{\Delta }_{},\mathrm{\Delta }_{})`$, with $`A_{\nu i}^{(0)}(x\mathrm{})=\mathrm{\Delta }_B\delta _{\nu i}`$ in bulk. In the A phase $`A_{\nu i}^{(0)}=a(x)\delta _{\nu x}(m_i^{(0)}+\mathrm{i}n_i^{(0)})`$, where $`m_i^{(0)}=\delta _{iy}`$, $`n_i^{(0)}=\pm \delta _{iz}`$ and $`a(x\mathrm{})=\mathrm{\Delta }_A`$. Here the A and B phase bulk gaps $`\mathrm{\Delta }_A`$ and $`\mathrm{\Delta }_B`$ are defined in Appendix B. The positive (negative) sign of $`n_t^{(0)}`$ chooses the vector $`\widehat{𝐥}=\widehat{𝐦}\times \widehat{𝐧}`$ to point in the positive (negative) $`x`$ direction, which are two degenerate configurations. A more familiar form for the spin part of the A phase order parameter is obtained by defining the $`\widehat{𝐝}`$ vector with $`\widehat{d}_\mu =R_{\mu \nu }\delta _{\nu x}=R_{\mu x}`$.

Assuming $`A^{(0)}`$ to be fixed, it suffices to write $`F_J`$ only in terms of the bulk order parameters, or more precisely, the soft degrees of freedom $`\varphi `$ and $`R`$. If we require $`F_J[A^L,A^R]`$ to be invariant under global gauge transformations and spin rotations, it should be unchanged when $`A^{L,R}`$ are both multiplied by $`R_{}^{L}{}_{}{}^{1}e^{\mathrm{i}\varphi ^L}`$. Thus we have $`F_J=F_J[A^{(0)L}(x),R_{}^{L}{}_{}{}^{1}R^RA^{(0)R}(x)e^{\mathrm{i}\mathrm{\Delta }\varphi }]`$, where $`\mathrm{\Delta }\varphi =\varphi ^R\varphi ^L`$. Now, since $`A^{(0)L}`$ and $`A^{(0)R}`$ are fixed by assumption, we see that $`F_J`$ can be parametrized simply as

$$F_J=F_J(\mathrm{\Delta }\varphi ,R_{\mu i}^LR_{\mu j}^R),i,j=x,y,z$$
(9)
where an orbital-space rotation matrix $`\psi _{ij}R_{\mu i}^LR_{\mu j}^R`$, $`i,j=x,y,z`$ remains as an argument. In Ref. viljas2,  this was derived only for the BB case, but we now see that it is valid more generally. For the A phase only the three components $`\widehat{d}_\mu =R_{\mu x}`$ out of the nine $`R_{\mu i}`$ are relevant. This implies that for the AA case Eq. (9) can be simplified as

$$F_J=F_J(\mathrm{\Delta }\varphi ,\widehat{𝐝}^L\widehat{𝐝}^R)$$
(10)
and for AB

$$F_J=F_J(\mathrm{\Delta }\varphi ,\widehat{𝐎}),$$
(11)
where we have defined an orbital-space vector $`\widehat{O}_i=\psi _{xi}=\widehat{d}_\mu R_{\mu i}`$.yipab

Often it is also instructive to consider the following “tunneling” form of the coupling energy:degennes

$$F_J^{\mathrm{tun}}=Re[a\stackrel{~}{A}_{\mu x}^L\stackrel{~}{A}_{\mu x}^R+b\stackrel{~}{A}_{\mu y}^L\stackrel{~}{A}_{\mu y}^R+c\stackrel{~}{A}_{\mu z}^L\stackrel{~}{A}_{\mu z}^R].$$
(12)
This is the lowest order term ($`\mathrm{\Delta }^2`$) in the order-parameter expansion discussed in Sec. 2. Thus it is also the leading term in the Fourier expansion (2): $`F_J^{(1)}E_J^{(1)}\mathrm{cos}\mathrm{\Delta }\varphi =F_J^{\mathrm{tun}}+O(\mathrm{\Delta }^4)`$. There are a few more assumptions coming into Eq. (12). The suppression of the order parameter near walls causes an ambiguity what location $`𝐫`$ should be used for $`A(𝐫)`$ in Eq. (12). Here we have used $`\stackrel{~}{A}^{L,R}`$ in the bulk, and defined $`\stackrel{~}{A}=A/\mathrm{\Delta }_{A,B}`$ depending on the phase. Also, we are assuming orthorhombic symmetry $`[(2/m)(2/m)(2/m)]`$ as for the slit above. The coefficients $`a`$, $`b`$ and $`c`$ are different for the different phase combinations and temperatures, and $`c=b`$ if the symmetry of the junction is cylindrical $`[(\mathrm{}/m)(2/m)]`$ instead.viljas2

### IV.2 Josephson Currents

Assume now the situation of Fig. 1, either 2D or 3D. There can be no current flow through solid walls, and thus $`𝐣_s\widehat{𝐬}=𝐣_\alpha ^{\mathrm{spin}}\widehat{𝐬}=0`$ on surfaces, where $`\widehat{𝐬}`$ is the local surface normal. Therefore, the conserved DC mass and spin currents flowing through the junction are given by ($`\mathrm{d}𝐒\widehat{𝐫}>0`$)

$`J_s`$ $`={\displaystyle _{\mathrm{\Gamma }^R}}d𝐒𝐣_s={\displaystyle _{\mathrm{\Gamma }^L}}d𝐒𝐣_s`$ (13)
$`J_\alpha ^{\mathrm{spin}}`$ $`={\displaystyle _{\mathrm{\Gamma }^R}}d𝐒𝐣_\alpha ^{\mathrm{spin}}={\displaystyle _{\mathrm{\Gamma }^L}}d𝐒𝐣_\alpha ^{\mathrm{spin}}.`$ (14)
These inherit their symmetry properties from the Josephson coupling $`F_J`$, and the expressions revealing this dependence can be derived as follows.

Consider infinitesimal variations of $`F_\mathrm{\Omega }[A]`$ around the stationary point, which are of the form shown in Eqs. (4) and (5) but now with local variational parameters $`\delta \varphi (𝐱)`$ and $`\delta \theta _\alpha (𝐱)`$. We choose them so that $`\delta \varphi |_{\mathrm{\Gamma }^L}=\delta \theta _\alpha |_{\mathrm{\Gamma }^L}=0`$ and $`\delta \varphi |_{\mathrm{\Gamma }^R}=\delta \varphi ^R`$, $`\delta \theta _\alpha |_{\mathrm{\Gamma }^R}=\delta \theta _\alpha ^R`$. Due to the stationarity, only the boundary term contributes to linear order, and using Eqs. (6) and (7) we have

$$\delta F_\mathrm{\Omega }=(\mathrm{}/2m_3)J_s\delta \varphi ^RJ_\alpha ^{\mathrm{spin}}\delta \theta _\alpha ^R.$$
(15)
On the other hand, we can similarly transform the boundary values $`\mathrm{\Delta }\varphi `$ and $`\psi _{ij}`$, so that $`\delta \mathrm{\Delta }\varphi =\delta \varphi ^R`$ and $`\delta \psi _{ij}=`$ $`ϵ_{\alpha \beta \gamma }R_{\beta i}^LR_{\gamma j}^R\delta \theta _\alpha ^R`$. By expanding $`F_J(\mathrm{\Delta }\varphi +\delta \varphi ^R,\psi _{ij}+\delta \psi _{ij})`$ $`=F_J(\mathrm{\Delta }\varphi ,\psi _{ij})+\delta F_J`$ to first order in $`\delta \varphi ^R,\delta \theta _\alpha ^R`$, and equating $`\delta F_J`$ with Eq. (15), we find

$`J_s`$ $`={\displaystyle \frac{2m_3}{\mathrm{}}}{\displaystyle \frac{F_J}{\mathrm{\Delta }\varphi }}`$ (16)
$`J_\alpha ^{\mathrm{spin}}`$ $`=ϵ_{\alpha \beta \gamma }R_{\beta i}^LR_{\gamma j}^R{\displaystyle \frac{F_J}{(R_{\mu i}^LR_{\mu j}^R)}},\alpha =x,y,z.`$ (17)
For a BB junction these expressions were derived in Ref. viljas2, , but again we see that they are more general. Using the simplified forms of the coupling energy for AA \[Eq. (10)\] and AB \[Eq. (11)\] cases, Eq. (17) reduces for an AA junction to

$`J_\alpha ^{\mathrm{spin}}`$ $`=[\widehat{𝐝}^L\times \widehat{𝐝}^R]_\alpha {\displaystyle \frac{F_J}{(\widehat{𝐝}^L\widehat{𝐝}^R)}},\alpha =x,y,z,`$ (18)
and for an AB junction and to

$`J_\alpha ^{\mathrm{spin}}`$ $`=ϵ_{\alpha \beta \gamma }\widehat{d}_\beta R_{\gamma i}{\displaystyle \frac{F_J}{\widehat{O}_i}},\alpha =x,y,z.`$ (19)
Note that the expansions of $`F_J`$ leading to Eqs. (16) and (17) work at most locally, *i.e.*, on each branch of the solution-space separately. In some simple cases one might also use them in order to check the consistency of one’s numerics. To remove the ambiguity related to the A-phase phase factor, we shall *always* choose $`\widehat{𝐥}=\pm \widehat{𝐱}`$ with $`\widehat{𝐧}^{(0)}=\pm \widehat{𝐳}`$, as explained above. This is a natural choice since $`\widehat{𝐧}`$-reversed configurations are related by time reversal. When the phase factor is included, TR should be understood as a reversal of both $`\varphi `$ and $`\widehat{𝐧}^{(0)}`$. On the other hand, a reversal of $`\widehat{𝐧}^{(0)}`$ alone is equivalent to reversal of $`\widehat{𝐥}`$ at constant $`\varphi `$.

### IV.3 Controlling the Boundary Conditions

Above we assumed the order parameter at the $`R_{\mathrm{}}`$ cutoff to be fixed to the form $`A(x)=RA^{(0)}(x)e^{\mathrm{i}\varphi }`$. Here $`A^{(0)}`$ determines the bulk phase, including its modification at walls. In case of the A phase, $`A^{(0)}`$ must always be such that $`\widehat{𝐥}\pm \widehat{𝐬}`$, where $`\widehat{𝐬}=\pm \widehat{𝐱}`$ is the wall normal, to minimize loss of condensation energy. Besides the phase angle, this leaves only the rotation matrix $`R`$ to be determined by some hydrodynamic interactions. In the A phase this reduces to fixing $`\widehat{𝐝}`$ with a combination of the dipole-dipole $`(\widehat{𝐝}\widehat{𝐥})^2`$ and the dipole-field $`(\widehat{𝐝}𝐇)^2`$ bulk terms. Unfortunately, for $`𝐇\widehat{𝐥}`$ the configuration remains undetermined. Furthermore, it may be difficult to produce the most interesting situations where $`\widehat{𝐝}^L\times \widehat{𝐝}^R0`$ by any physical means. In principle, one way is to use magnetic fields of different directions on the two sides. (Some more exotic ways, such as an A-B interface, could be imagined.nishida ) In the B phase there is more freedom in controlling $`R`$. To start with, we always assume the rotation angle $`\theta `$ to be in the minimum $`\theta _00.58\pi `$ of the dipole-dipole $`(\mathrm{cos}\theta +1/4)^2`$ interaction far from the junction.viljas2  The remaining $`\widehat{𝝎}`$ vector is coupled to the wall (with normal $`\widehat{𝐬}`$) by a surface-dipole term $`(\widehat{𝝎}\widehat{𝐬})^2`$, and to magnetic field via a bulk term $`(\widehat{𝝎}𝐇)^2`$ and a surface term $`(𝐇R\widehat{𝐬})^2`$. As discussed in Ref. viljas2, , it is the surface-dipole and the surface-field terms which determine $`\widehat{𝝎}`$ close to surfaces at low and high fields, respectively. For the magnetic surface configurations we use the definitions A-D of Ref. yippi, , but refer to them with lowercase letters a-d.

In addition to the rotation matrix, also the phase difference $`\mathrm{\Delta }\varphi `$ must be specified. In a channel geometry this would be replaced by specifying the phase gradient, *i.e.*, the superfluid velocity.ullahfetter  We assume that the inverse Josephson frequency $`\omega _J^1=\mathrm{}/2\mathrm{\Delta }\mu `$ is much larger than all order parameter relaxation times (except that of $`\mathrm{\Delta }\varphi `$). In practice this should be well satisfied, since $`\omega _J`$ tends to be no higher than in the audio regime.berkeley  This makes the equilibrium concepts of energy-phase $`F_J(\mathrm{\Delta }\varphi )`$ and current-phase $`J_s(\mathrm{\Delta }\varphi )`$ relations sensible, and therefore warrants the present calculation. However, the situation in the A phase may again be more complicated due to the orbital viscosity phenomenon, which slows down the dynamics of the $`\widehat{𝐥}`$ vector.vollhardt  Therefore we cannot expect the calculation to properly describe the dynamics of phase slips, *i.e.*, jumps between branches of $`J_s(\mathrm{\Delta }\varphi )`$.

### IV.4 About the Numerical Implementation

The cutoff $`R_{\mathrm{}}`$ is arbitrary, but from the CPR’s for one choice, we may determine the CPR’s for any other one.thuneberg88  For example, we may do the numerical calculation for $`r<R_c`$ where $`R_c`$ is small (see Fig. 1), and then use Eq. (26) or Eq. (34) to correct the phase differences for some $`r=R_{\mathrm{}}R_c`$, so that the new CPR becomes a “slanted” version of the original. Whenever spin currents are present, the spin-rotation boundary conditions should also be modified using Eq. (27) or Eq. (35). However, in the B phase this will lead to rotation matrices which are not of the form $`R(\widehat{𝝎},\theta _0)`$ that we wish to have for large $`R_{\mathrm{}}`$ (although still $`R_{\mathrm{}}\xi _D`$).

Thus the calculation should in general be done for each (large) cutoff radius $`R_{\mathrm{}}`$ and each set of spin-orbit boundary conditions separately. The numerical minimization may still be done inside $`R_cR_{\mathrm{}}`$, and for $`r>R_c`$ the asymptotic solutions of Appendix A are used. In our calculations, the fitting of the solutions at $`R_c`$ is not done by comparing derivatives at any location on the boundary, but by the more physical requirement of *conservation of total mass and spin currents* over $`R_c`$. This leads to a somewhat cumbersome self-consistent iteration procedure, which is described in Appendix B with the rest of the numerics. Fortunately, these practical issues should not affect any of the above or the following analysis, except through the value of $`R_{\mathrm{}}`$, and thus a discussion of the solutions in the asymptotic region is postponed to Appendix A. If these “asymptotic corrections” are not used, then $`R_{\mathrm{}}=R_c`$.

## V CPR’S FOR AA, BB AND AB JUNCTIONS

In this section, we present a more careful analysis of the BB, AA, and AB junctions. In what follows, we always present results for BB at vapor pressure $`p=0`$ bar, for AA at melting pressure $`p=34.4`$ bar, and for AB at roughly the coexistence pressure of $`p=28.7`$ bar ($`\mathrm{\Delta }f_{AB}=0`$).

In numerical calculations of the CPR’s, the asymptotic corrections were usually taken into account, and we used the outer cutoff $`R_{\mathrm{}}/\xi _{GL}30`$. In this case the value for the inner cutoff was roughly $`R_c/\xi _{GL}=10`$, although in principle the results should be quite independent of it. As expected, the logarithmic decay of the phase corrections is slow, and as we varied $`R_c/\xi _{GL}`$ and $`R_{\mathrm{}}/\xi _{GL}`$ in the range $`10\mathrm{}60`$, only small differences in the CPR’s could be seen. Most of the time a lattice spacing of $`\mathrm{\Delta }x/\xi _{GL}=\mathrm{\Delta }y/\xi _{GL}=0.5`$ was used. Refinement did not lead to qualitative differences in the form of the CPR’s, and only to differences of at most a few percent in the critical currents. The minimization was usually carried out until the error (as measured by the norm of the gradient in the $`R_c`$ region, see Appendix B) was smaller than $`10^5`$. For given boundary conditions, an accuracy better by more than 10 orders was possible, but completely unnecessary, considering the uncertainties related to the handling of the bulk cutoffs. Current conservation at every lattice point was practically exact, except for points on the $`R_c`$ cutoff, where only an accuracy of $`|\mathbf{}𝐣_s|10^2`$ was achieved locally. However, the total current over the cutoff was required to be conserved down to $`10^4`$, if the asymptotic corrections were used. These values refer to the $`\rho _s`$ and $`\rho _{}`$ units used in Appendix A and B, which are used also in all the figures below.

### V.1 BB Junctions

The case of BB junctions was already considered in Ref. viljas2, , and we skip most of the analysis. Note that for BB time-reversal symmetry (complex conjugation of $`A`$) implies the simple relation $`F_J(\mathrm{\Delta }\varphi )=F_J(\mathrm{\Delta }\varphi )`$ in addition to $`2\pi `$ periodicity, and thus $`J_s(\mathrm{\Delta }\varphi )=J_s(\mathrm{\Delta }\varphi )`$, $`J_\alpha ^{\mathrm{spin}}(\mathrm{\Delta }\varphi )=J_\alpha ^{\mathrm{spin}}(\mathrm{\Delta }\varphi )`$. Inserting Eq. (12) with $`\stackrel{~}{A}=R\mathrm{exp}(\mathrm{i}\varphi )`$ into Eq. (17) gives the spin current as $`J_\alpha ^{\mathrm{spin}}=ϵ_{\alpha \mu \nu }[aR_{\mu x}^LR_{\nu x}^R`$ $`+bR_{\mu y}^LR_{\nu y}^R+cR_{\mu z}^LR_{\nu z}^R]\mathrm{cos}\mathrm{\Delta }\varphi `$. Notice that even if $`F_J^{(1)}`$ is identically zero, and thus the EPR is $`\pi `$ periodic, the spin current may remain $`2\pi `$ periodic.

Despite the apparent simplicity, the BB junction has a rich spin current structure. This is because, in addition to the Josephson spin currents induced by bulk boundary conditions, there exist spontaneous spin currents close to solid surfaces.zhangwall  An example of these is shown in Fig. 2, where small spin current loops are seen to be stabilized in different parts of the weak link. However, these are not “topological” vortices with a quantized spin velocity circulation driving them. The cores are usually associated with a small A-phase-like “orbital magnetization” $`ϵ_{ijk}Im(A_{\mu j}^{}A_{\mu k})/2`$ along the $`z`$ axis.

Figures 3 and 4 show the CPR’s for apertures of two different widths $`W`$, but the same “depth” $`D`$, *i.e.*, thickness of wall, using the same series of boundary conditions: $`\widehat{𝝎}^R=\widehat{𝐱}`$, and $`\widehat{𝝎}^L`$ is in the $`xy`$ plane with azimuthal angle $`\eta _L`$ from the positive $`x`$ direction. These figures are equivalent to Fig. 11 of Ref. viljas2,  and, while corresponding to rather imaginary experimental conditions, they show quite clearly the essential differences between the CPR’s of narrow and wide apertures. In the narrow case (Fig. 3) the behavior follows well the perturbative picture of Eq. (2): the energy-phase and current-phase (mass or spin) relations are non-hysteretic and close to $`2\pi `$ periodic (co)sine functions. Only in cases where $`F_J^{\mathrm{tun}}`$ $`aR_{\mu x}^LR_{\mu x}^R`$ $`+bR_{\mu y}^LR_{\mu y}^R+cR_{\mu z}^LR_{\mu z}^R`$ vanishes (as it obviously does for $`\eta _L0.4\pi `$ and $`\eta _L0.7\pi `$) a small $`\pi `$ periodic contribution due to the remaining higher-order terms is visible, but with very small critical current. These higher harmonics become stronger and produce more pronounced $`\pi `$ states at low temperatures.yippi  For $`\eta ^L0.6\pi `$ the CPR is sinusoidal but “$`\pi `$ shifted”, whereas for $`\eta ^L\pi `$ the $`\pi `$ shift again vanishes. This is different from the pinhole case of Ref. viljas2, , where a $`\pi `$ shift is obtained for $`\widehat{𝝎}^L=\widehat{𝝎}^R=\pm \widehat{𝐱}`$. This is apparently due to the difference in the symmetry of the aperture: unlike for the pinhole, here we have $`cb`$.

On the other hand, in the wider case (Fig. 4) the CPR’s are already clearly hysteretic. The leftmost panels have $`\widehat{𝝎}^L=\widehat{𝝎}^R=\widehat{𝐱}`$, which is the situation studied in Ref. viljas,  for the circular 3D aperture. Here the behavior is similar: in addition to the usual “$`0`$ branch” there exists also a “$`\pi `$ branch” around $`\mathrm{\Delta }\varphi =\pi `$. Figure 5 depicts the corresponding behavior of the order parameter inside the aperture at exactly $`\mathrm{\Delta }\varphi =\pi `$ (compare with Fig. 3 of Ref. viljas, ). Here the constant rotation matrix has been dropped (*i.e.*, $`\theta =0`$) so that $`A`$ is diagonal at the infinities. On the “0 branch” the behavior of $`𝚫(𝐱,\widehat{𝐤})`$ is singular, whereas on the “$`\pi `$ branch” $`𝚫`$ avoids a singularity by “escaping into the third dimension”. The rotation of $`𝚫`$ around $`\widehat{𝐲}`$ with growing $`x`$ means that there is a spontaneous Josephson current of $`y`$ directional spin flowing through the aperture. Adding the spin rotation by $`\theta =0.58\pi `$ around $`\widehat{𝐱}`$ gives the $`y`$ and $`z`$ spin current components visible in Fig. 4. As pointed out in Ref. viljas, , the order parameter on the $`\pi `$ branch has a similar structure as in a double-core vortex,thuneberg87  where the virtual vortex cores are inside the wall, one on each side of the channel.

For parallel $`\widehat{𝝎}`$ vectors there can never be a jump from the 0 to the $`\pi `$ branch, and so the conclusion in Ref. viljas,  was that the branch may not be experimentally achievable. But this is not so. Once the rotational symmetry of the bulk boundary conditions around the surface normal is broken by perturbing $`\widehat{𝝎}^L`$ even slightly, the hysteresis of the $`0`$ branch is reduced, and the $`\pi `$ branch becomes the branch of minimum energy for $`\mathrm{\Delta }\varphi \pi `$. Since $`J_s^{}(\mathrm{\Delta }\varphi =\pi )>0`$, this “$`\pi `$ state” is even a (local or global) minimum of EPR and thus stable against small perturbations of $`\mathrm{\Delta }\varphi `$. Notice how in Fig. 4 the order parameter transitions between $`0`$ and $`\pi `$ minima continue to be associated with changes in the spin currents, although a simple interpretation as with Fig. 5 is no longer possible for large tilting angles of $`\widehat{𝝎}^L`$. But for small angles a jump onto, or away from the $`\pi `$ branch can be interpreted as a *phase slip by a half-quantum vortex*.viljas  Note also that $`\mathrm{\Delta }\varphi =0`$ tends to remain at least a local minimum of the EPR, *i.e.*, $`J_s^{}(\mathrm{\Delta }\varphi =0)>0`$, so that a sinusoidal but apparently $`\pi `$ shifted EPR (“$`\pi `$-junction”) tends to be avoided. In the narrow-aperture limit $`W/\xi _{GL}3`$ the situation $`J_s^{}(0)<0`$ appears to be more easily obtained (compare Figs. 3 and 4). The last panels of Fig. 4 correspond to $`\widehat{𝝎}^L=\widehat{𝝎}^R=\widehat{𝐱}`$, for which a $`\pi `$ branch with $`J_s^{}(\mathrm{\Delta }\varphi =\pi )<0`$ was found already in Ref. thuneberg88, . Our simulations confirm this old result, and show that the $`J_s^{}(\mathrm{\Delta }\varphi =\pi )<0`$ branch of Ref. thuneberg88,  is related to the branches with $`J_s^{}(\mathrm{\Delta }\varphi =\pi )>0`$ (stable $`\pi `$ state) by a continuous variation of the boundary conditions.

The limit between the two types of behavior, Figs. 3 and 4, is rather clear-cut, and it only depends on the width $`W`$. The transition occurs roughly at the width $`W/\xi _{GL}=3`$, which coincides with the analytic critical value $`\pi `$ for destruction of superfluidity in an infinite slab.fetterullah ; kjaldman  However, in the weak link case the order parameter amplitudes and critical currents always remain finite, although small. In wide apertures the CPR’s with $`\pi `$ states can have relatively *large critical currents* — comparable to those of $`0`$ branches in general. In the narrow apertures, however, there tends to be almost an order-of-magnitude difference between $`E^{(1)}`$ and $`E^{(2)}`$ in the GL regime.

Figure 6 shows the CPR’s for a more physical control parameter, namely the angle $`\theta _H`$ of an applied magnetic field. Here we must assume that the field is strong enough for the texture to be determined by the magnetic surface interaction $`(𝐇R\widehat{𝐬})^2`$, but small enough for $`\xi _H\xi _{GL}`$, since we have neglected the magnetic terms from our GL free energy. We choose to show the “bd” configuration defined in Ref. yippi,  with $`𝐇`$ in $`zy`$ plane, since this is an example of a “$`\pi `$ junction” with $`W/\xi _{GL}>3`$. In fact, a rather similar figure as Fig. 6 exists also for the “ad”, “ac” and “bc” configurations with field in the plane of the wall. For the “ab” or “cd” configurations no $`\pi `$ states nor $`\pi `$ shifts were seen, although many magnetic field directions (both in-plane and off-plane) were checked.

Based on these results we conclude that a rich variety of $`\pi `$ states and $`\pi `$ shifts can exist in single BB apertures in the GL regime, although it may be difficult to realize the required $`\widehat{𝝎}`$ configurations experimentally. Approximately the behavior of Fig. 4, where the unreachable $`\pi `$ branch becomes visible, could be obtainable by slowly turning on $`𝐇`$ while keeping it in a suitably chosen direction. Without a doubt, a qualitatively similar division into small and large aperture behavior will be valid at low temperatures also. To verify this, a quasiclassical calculation for the general-size weak link would be required. However, this appears to be too demanding considering that there is no reason to expect any essentially new phenomena. The critical currents and a phase diagram for BB are summarized in the next section.

### V.2 AA$``$ and AA$``$ Junctions

In AA junctions the time reversal symmetry is no longer quite so simple as for BB, due to the fact that complex conjugation also reverses the direction of $`\widehat{𝐥}`$. In both the AA$``$ and AA$``$ configurations TR reduces only to $`F_J(\widehat{𝐥}^L,\widehat{𝐥}^R,\beta \mathrm{\Delta }\varphi )=F_J(\widehat{𝐥}^L,\widehat{𝐥}^R,\beta +\mathrm{\Delta }\varphi )`$, where $`\beta =0`$ or $`\pi `$. (Note again that reversing $`\widehat{𝐥}`$ is done by flipping the reference direction $`\widehat{𝐧}^{(0)}`$.) However, reversing the directions of *both* $`\widehat{𝐥}`$’s simultaneously has no effect in our orthorhombically symmetric junction, so that $`F_J(\widehat{𝐥}^L,\widehat{𝐥}^R,\mathrm{\Delta }\varphi )=F_J(\widehat{𝐥}^L,\widehat{𝐥}^R,\mathrm{\Delta }\varphi )`$, and therefore the TR symmetry reduces to that of the BB case. The current-phase relations therefore satisfy $`J_s(\beta \mathrm{\Delta }\varphi )=J_s(\beta +\mathrm{\Delta }\varphi )`$ and $`J_\alpha ^{\mathrm{spin}}(\beta \mathrm{\Delta }\varphi )=J_\alpha ^{\mathrm{spin}}(\beta +\mathrm{\Delta }\varphi )`$ as for BB. Since $`\pm \widehat{𝐥}^{L,R}`$ correspond to the same branch(es), the symmetry also implies that the minimum of EPR is at either $`\mathrm{\Delta }\varphi =0`$ or $`\pi `$, depending on the directions of the $`\widehat{𝐝}`$ vectors. In our case $`\mathrm{\Delta }\varphi =0`$ is the minimum if $`\widehat{𝐝}=\widehat{𝐥}`$ for both $`L`$ and $`R`$.

For AA$``$ the critical current will aways vanish if the aperture has full rotation symmetry around the wall normal $`\widehat{𝐬}=\pm \widehat{𝐱}`$. This is because an orbital rotation by $`\theta `$ around $`\widehat{𝐬}`$ adds a phase $`\theta `$ to the order parameter if $`\widehat{𝐥}=\pm \widehat{𝐬}`$, and therefore $`\mathrm{\Delta }\varphi =\varphi ^R\varphi ^L`$ is changed by $`\pm 2\theta `$. Since this rotation is a symmetry operation for the aperture, it should not change $`F_J`$ (10). Thus $`F_J`$ should have the periodicity $`F_J(\mathrm{\Delta }\varphi +2\theta )=F_J(\mathrm{\Delta }\varphi )`$ for *all* $`\theta `$. This can only be satisfied if the phase dependence of $`F_J`$ vanishes altogether. In general, if the aperture has $`n`$-fold rotation symmetry ($`n1`$), the $``$ EPR’s are $`2(2\pi /n)`$ periodic. For example, for the orthorhombic slit considered in this paper $`n=2`$ and only a rotation by $`\theta =\pi `$ must be a symmetry operation. This implies no additional restrictions on $`F_J`$, since it is already $`2\pi `$ periodic. The case AA$``$ has nontrivial $`2\pi `$ periodic $`F_J`$’s irrespective of the geometry, except for the special case $`\widehat{𝐝}^L\widehat{𝐝}^R=0`$ (see below).

Figure 7 represents the different configurations for $`\widehat{𝐥}`$ inside the aperture. The middle one for both $``$ (b) and $``$ (e) is symmetrical with respect to the $`x`$ axis. Case (b) is very singular and does not usually correspond to the minimum energy in large apertures, but may be metastable at least for $`\mathrm{\Delta }\varphi 0`$. Here a radial singularity is depicted, but a “hyperbolic” one is also possible. In some cases $`\widehat{𝐥}`$ can also escape from the plane. The $``$ configuration (e), on the other hand, is usually the minimum-energy configuration for $`WD`$ and $`\mathrm{\Delta }\varphi 0`$, as may be expected. The other configurations correspond to pairs of degenerate cases, where $`\widehat{𝐥}`$ bends asymmetrically, as shown in Fig. 8 in more detail. Actually, the quantity shown is $`ϵ_{ijk}Im(\stackrel{~}{A}_{\mu j}^{}\stackrel{~}{A}_{\mu k})/2`$, where $`\stackrel{~}{A}=A/\mathrm{\Delta }_A`$, which reduces to $`\widehat{l}_i`$ in the bulk. Figure 8 also shows the mass current density fields associated with a phase difference $`\mathrm{\Delta }\varphi =0.3\pi `$. Because the superfluid density is twice as large in directions perpendicular to $`\widehat{𝐥}`$ than for those in parallel with it, the bending of the $`\widehat{𝐥}`$ fields results in asymmetrical current distributions close to the junction. These effects result in a more complicated numerical problem than for the isotropic B phase (see Appendix A).

Figure 9 shows CPR’s for $``$ and an aperture of fixed aspect ratio, but different dimensions. No $`\pi `$ branches exist as for BB, and the corresponding figure for $``$ would be very similar. It was noticed that, both in $``$ and $``$ cases, when a jump to another branch of $`J_s(\mathrm{\Delta }\varphi )`$ occurred, the asymmetric texture shifted from one degenerate configuration (a or d) in Fig. 7 to the other (c or f). At least in the $``$ case this could be understood by remembering that a phase slip of $`\mathrm{\Delta }\varphi `$ by $`2\pi `$ is equivalent to a rotation of the orbital space around $`\widehat{𝐬}`$ by $`\pi `$. However, these results are due to the minimization algorithm, and do not necessarily correspond to physical dynamics.

As pointed out above, for AA$``$ at small $`D`$ and $`W/\xi _{GL}3`$ the minimum-energy configuration is that of Fig. 7e when $`\mathrm{\Delta }\varphi =0`$ and no current flows. But when $`\mathrm{\Delta }\varphi `$ is increased, a change to one of the asymmetric $`\widehat{𝐥}`$ configurations (d or f) will occur, which causes a sudden reduction of the current and thus a discontinuity in CPR. This is illustrated in Fig. 10 — see also Figs. 14 and 15 below. The larger $`D`$ is, the smaller is the critical $`\mathrm{\Delta }\varphi `$ for this event. An accurate calculation of this critical value is difficult, however, since it depends on whether or not the asymptotic correction for $`\widehat{𝐥}`$ is used, and on how far the $`R_c`$ cutoff is taken. This introduces some amount of additional uncertainty in the forms of the CPR’s. For AA$``$ one of the asymmetric configurations (a or c) is always the minimum one (except for small $`D`$ and $`W`$, see Fig. 15) and thus in Fig. 9 all the CPR’s are smooth and well-behaved.

For $`\widehat{𝐝}^L\widehat{𝐝}^R=0`$ the CPR’s are $`\pi `$ periodic in both AA configurations, as can be expected by noticing that $`F_J^{\mathrm{tun}}`$ \[Eq. (12)\] vanishes in this case. However, as discussed in Section 4, there appears to be no easy physical way in which a situation with $`\widehat{𝐝}^L\times \widehat{𝐝}^R0`$ could be achieved. According to Eq. (18) this more-or-less rules out spin currents, and thus seems to render the AA CPR’s rather un-interesting, apart from the effects of Fig. 10 arising from the $`\widehat{𝐥}`$ texture. Theoretically this is not a problem, of course, and if we force the configuration $`\widehat{𝐝}^L\widehat{𝐝}^R=0`$, then the following observations can be made. In the small-hole limit (CPR’s continuous) the critical currents are considerably suppressed, and $`J_s(\mathrm{\Delta }\varphi )J_c\mathrm{sin}(2\mathrm{\Delta }\varphi )`$. However, in the large-hole limit (CPR’s hysteretic), the critical currents are not significantly affected. The CPR’s look just like in Figs. 9 or 10, but with the $`\pi `$ branches being added in between each $`0`$ branch. As the $`\widehat{𝐝}`$’s are gradually changed from $`\widehat{𝐝}^L\widehat{𝐝}^R=0`$, either the $`0`$ or the $`\pi `$ branches begin to rise in energy, and the corresponding CPR branches get smaller (similar to metastable “$`\pi `$ states” in BB). Finally the branches become unreachable and $`2\pi `$ periodicity without any $`\pi `$ states is regained. Again, these observations hold for both AA$``$ and AA$``$ alike. The most notable difference between them is in the critical currents, which are summarized in the next section.

### V.3 AB Junctions

The case of an AB junction is the most complicated one to analyze. In this case the requirement of time reversal symmetry implies $`F_J(\widehat{𝐥},\beta \mathrm{\Delta }\varphi )=F_J(\widehat{𝐥},\beta +\mathrm{\Delta }\varphi )`$, where $`\beta =0`$ or $`\pi `$ as in AA. This no longer simplifies to the intuitive BB form as it did for our AA junctions. Rather the EPR’s for $`\widehat{𝐥}=\pm \widehat{𝐱}`$ form completely separate, unsymmetrical branches, which are $`2\pi `$ periodic in general, and which can both still separate into hysteretic sub-branches as usual. The minimum of the EPR is also not restricted to $`\mathrm{\Delta }\varphi =0`$ or $`\pi `$ on either branch (see Fig. 11 below). Similarly, the currents only satisfy relations between the $`\pm \widehat{𝐥}`$ branches $`J_s(\widehat{𝐥},\beta \mathrm{\Delta }\varphi )`$ $`=J_s(\widehat{𝐥},\beta +\mathrm{\Delta }\varphi )`$ and $`J_\alpha ^{\mathrm{spin}}(\widehat{𝐥},\beta \mathrm{\Delta }\varphi )=J_\alpha ^{\mathrm{spin}}(\widehat{𝐥},\beta +\mathrm{\Delta }\varphi )`$. The roles of $`\mathrm{\Delta }\varphi =0`$ and $`\pi `$ can again be interchanged by flipping the $`\widehat{𝐝}`$ of the A phase.

To analyze the AB case further, we consider the tunneling term, Eq. (12). Assuming A phase on the left and B phase on the right, we again introduce the orbital space vector $`\widehat{O}_i=\widehat{d}_\mu R_{\mu i}`$ and find that $`F_J^{\mathrm{tun}}=Re\{[b\widehat{O}_y(\widehat{m}_y\pm \mathrm{i}\widehat{n}_y)`$ $`+c\widehat{O}_z(\widehat{m}_z\pm \mathrm{i}\widehat{n}_z)]\mathrm{exp}(\mathrm{i}\mathrm{\Delta }\varphi )\}`$. We see that if $`\widehat{d}_\mu =\pm R_{\mu i}\widehat{s}_i`$, *i.e.*, $`\widehat{𝐎}\widehat{𝐬}`$, then $`\widehat{𝐎}\widehat{𝐦},\widehat{𝐧}`$ and $`F_J^{\mathrm{tun}}`$ vanishes. This is true, for example, in the most symmetric (and, in the absence of magnetic fields, most probable) situation where both $`\widehat{𝐝}`$ and the B phase rotation axis $`\widehat{𝝎}`$ are perpendicular to the wall, but not necessarily if either of them deviates from this configuration. Incidentally, $`\widehat{d}_\mu =\pm R_{\mu i}\widehat{s}_i`$ is also the minimum-energy configuration for a free, planar AB interface, when $`\widehat{𝐬}`$ is the interface normal.cross77

More general conclusions can be based on the general form $`F_J`$ (11) in the symmetric case $`\widehat{𝐎}\widehat{𝐬}`$. Similarly as in the case of AA$``$, a rotation of the orbital space around $`\widehat{𝐬}`$ by angle $`\theta `$ is now equivalent to a shift of the phase difference $`\varphi `$ by $`\pm \theta `$. (The factor is different here since the B phase side is not affected by the rotation.) In the case of a circularly symmetrical aperture we must then have $`F_J(\mathrm{\Delta }\varphi +\theta )=F_J(\mathrm{\Delta }\varphi )`$ for all $`\theta `$ which, again, means that the phase dependence must vanish to all orders. For the slit junction this implies only $`F_J(\mathrm{\Delta }\varphi +\pi )=F_J(\mathrm{\Delta }\varphi )`$, *i.e.*, that the energy-phase relations are $`\pi `$ periodic. Generally, if the aperture has $`n`$-fold rotation symmetry, then the CPR’s must be $`2\pi /n`$ periodic. However, these simple conclusions no longer hold when $`\widehat{𝐎}\widehat{𝐬}`$, and in general the nontrivial $`2\pi `$ periodic behavior is obtained in any aperture geometry.

Figure 11 shows examples of CPR’s when $`\widehat{𝐝}`$ and $`\widehat{𝝎}`$ vectors are controlled by a field $`𝐇`$ strong enough so that $`\widehat{𝐝}`$ points in the direction corresponding to minimum possible dipole energy $`(\widehat{𝐝}\widehat{𝐥})^2`$ allowed by a strict condition $`\widehat{𝐝}𝐇=0`$. Here $`\widehat{𝝎}`$ is chosen to be in the “a” configuration.yippi  The failure of the phase-inversion symmetry $`F_J(\mathrm{\Delta }\varphi )=F_J(\mathrm{\Delta }\varphi )`$ is clearly visible; Figure 11 is for $`\widehat{𝐥}=\widehat{𝐱}`$, and the branch for $`\widehat{𝐥}=\widehat{𝐱}`$ is obtained by using the above symmetries for $`F_J`$. Figure 12, on the other hand, shows the CPR’s for a fixed aspect ratio $`2:1`$ (as in Fig.9 for AA$``$) but different dimensions in the $`\widehat{𝐎}\pm \widehat{𝐬}`$ case. The CPR’s are $`\pi `$ periodic as they should, and even show the presence of strong $`\pi /2`$-periodic admixtures, or separate “$`\pi /2`$ branches”.

The order parameter profile on the $`x`$ axis which goes through the AB interface, and the form of the $`\widehat{𝐥}`$ field there are shown in Fig. 13. Here we have $`\mathrm{\Delta }f_{AB}>0`$ ($`p>p_0`$) for illustration purposes, although the CPR’s were calculated for $`\mathrm{\Delta }_{AB}=0`$. The boundary of the $`\widehat{𝐥}`$ field, *i.e.*, the A-B interface is thus seen to bulge into the B phase. Even at $`\mathrm{\Delta }f_{AB}=0`$ the interface always settles into the B-phase end of the channel.viljasab  As in the AA cases, $`\widehat{𝐥}`$ tends to bend parallel to $`\widehat{𝐲}`$ inside the hole. This is the preferred configuration also for the A-B interface itself. The order parameter components shown on the right are the same as in Ref. schopohl,  or Ref. thuneberg91, , except for the components arising from the bending of $`\widehat{𝐥}`$ for $`x\mathrm{}`$. Similarly to the AA cases, a phase slip (jump from branch to another) was usually found to be associated with a transition of the $`\widehat{𝐥}`$ texture. Textures with mirror-symmetry with respect to the $`xz`$ plane and $`\widehat{𝐥}`$ pointing out of the plane in the $`\pm \widehat{𝐳}`$ direction at the interface were also often seen on some branches.

## VI SUMMARY AND PHASE DIAGRAMS

Figure 14 summarizes the critical mass currents $`J_c`$ for BB, AA, and AB junctions, plotted in the units defined in Appendix B. They have been calculated for several sizes in the range $`W/\xi _{GL},D/\xi _{GL}=0\mathrm{}10`$, which is certainly the most relevant one. For essentially larger dimensions the CPR’s become strongly multivalued, and the aperture cannot be considered as a weak link. Figure 15, on the other hand, is a phase diagram which relates the changeover between continuous and hysteretic CPR’s, and some of the $`\widehat{𝐥}`$ texture transitions to well-defined regions of $`D,W`$ plane. The results were calculated using a lattice spacing of $`\mathrm{\Delta }x/\xi _{GL}=\mathrm{\Delta }y/\xi _{GL}=0.5`$, and with a cutoff $`R_c/\xi _{GL}=20`$, but with no asymptotic corrections. In AA and AB, some amount of error will necessarily remain close to regions where $`\widehat{𝐥}`$ transitions occur, since $`\widehat{𝐥}`$ is not allowed to vary freely beyond $`R_c`$, which makes the texture somewhat too rigid. All of the results of Figs. 14 and 15 are for the most symmetric configurations only: in BB $`\widehat{𝝎}^L=\widehat{𝝎}^R`$, in AA $`\widehat{𝐝}^L=\pm \widehat{𝐝}^R`$, and in AB $`\widehat{d}_\mu =R_{\mu i}\widehat{s}_i`$. However, at least the narrow-aperture limit $`W/\xi _{GL}3`$ is independent of any the boundary conditions. The regime with $`W/\xi _{GL}3`$ (for small $`D`$) is where the description due to Eq. (2) is supposed to be valid, whereas for larger apertures hysteresis begins to set in. This transition in behavior is seen in Figs. 9 and 12.

In Fig. 14 the BB case is actually well-known. In Ref. ullahfetter,  it was calculated using a channel geometry, rather than the empty half-spaces on the two sides. The fact that the values in this reference are slightly larger may be just an indication of the restricted form of order parameter which was assumed there. Thus we confirm the expectation that the critical currents should only depend on the properties of the weak link itself, not the way in which the current is driven through it.

The other three panels of Fig. 14 show the interesting feature that if there is A phase on one of the sides, the critical currents will always approach the same values for large $`D`$. They are thus determined by the inside of channel, and not the complicated $`\widehat{𝐥}`$ structures or the A-B interface at the ends of the junction. The transition to the long-channel limit involves locking of $`\widehat{𝐥}`$ perpendicular to the channel walls ($`\widehat{𝐥}=\pm \widehat{𝐲}`$) at the middle of the junction, as in a parallel-plate geometry with plate separation $`W`$.fetterullah  For AA$``$ this happens smoothly, but in AA$``$ and AB there are rather sharp cusps in the $`J_c`$ curve beyond which $`J_c`$ is almost constant. In both cases this transition is due to a crossing over of the critical currents of two branches of CPR with different $`\widehat{𝐥}`$ configurations. For $`AA`$ the configurations are the symmetric (e) and the asymmetric (d or f) ones in Fig. 7, as discussed in the previous section. Some of the associated CPR’s were shown in Fig. 10, which corresponds to the $`W/\xi _{GL}=6`$ line of Fig. 14 — for small $`\mathrm{\Delta }\varphi `$ the symmetric $`\widehat{𝐥}`$ configuration is the one to determine $`J_c`$. The $`J_c`$-transition lines are shown in the phase diagram, Fig. 15, although it should be kept in mind that their positions may depend on the way in which the numerics was done. For small $`D`$ the behaviors of all four cases are quite different. For AA$``$ and AB $`J_c`$ tends to be suppressed, whereas for BB and AA$``$ it increases. In AB there are actually $`J_c`$ transitions at two different $`D`$ values. Both are due to switching of $`J_c`$ between the $`0`$ ($`+`$) or $`\pi /2`$ ($``$) CPR branches, which are seen in Fig. 12. The $`+/`$ regions are shown in Fig. 15. Notice that all these details are only valid for $`\widehat{d}_\mu =\pm R_{\mu i}\widehat{s}_i`$. For other configurations the $`F_J^{\mathrm{tun}}`$ term is also finite, and the critical currents will in fact tend to be larger (for small $`D`$) than presented in the figures above.

Also shown in Fig. 14 are the critical currents for some apertures of varying size, but fixed aspect ratio. These correspond essentially to Figs. 9 and 12 for $`AA`$ and AB, and to equivalent (but not shown) figures for BB and AA$``$. Small numerical differences are due to the different computational parameters. These curves are interesting because all the temperature dependence in the GL regime is in the length scale $`\xi _{GL}(T)(1T/T_c)^{1/2}`$ alone, and scaling $`D/\xi _{GL},W/\xi _{GL}`$ with factor $`l`$ is thus equivalent to changing the temperature: $`l(1T/T_c)^{1/2}`$. The aperture in the Paris experiments reported in Ref. paris,  was a slit of dimensions $`0.18\times 2.6\mu `$m2 in a $`0.1\mu `$m thick SiN membrane. The aspect ratio $`W:D`$ for this is close to $`2:1`$, which was therefore most frequently used in our calculations. The behavior at small $`l`$ appears to follow a power law $`J_c/J_{A,B}l^\delta `$. For all of the three BB cases the exponent is roughly the same, $`\delta =1.8`$. Similarly, for both AA$``$ and AA$``$ it is $`\delta =2.2`$, whereas for AB $`\delta =4.5`$. In BB the $`J_c`$ curves appear to be linear for large $`l`$, but this is only illusory. To point out another numerical accident, note that in the AA$``$ case the critical current transition appears to occur at exactly the aspect ratio $`W:D=2:1`$. There seems to be no obvious reason for this.

However, the striking similarity in the $`J_c`$ scaling behavior of the AA$``$ and AA$``$ configurations for small $`l`$ appears very surprising at a first glance. Since in the pinhole limit the critical current for AA$``$ vanishes, should it not at least vanish faster as the smaller dimensions are approached? Two things should be realized here. The first is due to the 2D nature of our model: as the slit is made very narrow, eventually all but the $`\widehat{𝐳}`$ directional component of the orbital part of the order parameter are suppressed. The remaining coupling through a *single* complex component is direction-independent, and thus the behavior for both AA cases becomes similar. Second, even for a 3D aperture, the “pinhole” behavior would not be exactly reached in this limit, because the quasiclassical pinhole calculation relies on a *nonlocal* coupling through the “$`f`$ functions”. The present GL calculation is *local*, and assumes that the *pairing potential* is everywhere nonzero inside the weak link — otherwise there could be no coupling at all. A different approach for a finite-size BB junction was recently adopted in Ref. endres, . These comments are also valid for the BB case, which explains the somewhat different result of Section 5 as compared with Ref. viljas2,

## VII CONCLUSIONS AND DISCUSSION

In the sections above we have presented a $`p`$-wave Ginzburg-Landau analysis of a two-dimensional weak-link model. We have studied the current-phase relations and mapped the critical currents of the weak link, when it is placed between two bulk volumes of either superfluid 3He-A or 3He-B. The BB case was studied before in several restricted calculations,monien ; ullahfetter ; thuneberg88 ; viljas  but here we were able to use much more general boundary conditions on the bulk phases. We showed that similar “$`\pi `$ states” as in BB pinholesyippi ; viljas2  exist also for small finite-size apertures in the GL regime. They are due to the second-order terms in Eq. (2) and are thus associated with small critical currents. On the other hand, we showed how this behavior gradually becomes more complicated in larger apertures, as the order parameter has more freedom to vary inside the hole. In this way, various kinds of $`\pi `$ states may exist together with either continuous or hysteretic CPR’s which have relatively large critical currents. This is consistent with the experimental findings of Ref. paris, , where, apparently due to the uncontrollability of the bulk textures, many kinds of CPR’s were seen. As in Ref. viljas2,  for the “anisotextural” effect in a large pinhole array, it was found that the hysteresis of CPR’s tends to be smaller when the bulk boundary conditions are less symmetric. This enabled jumps to the $`\pi `$ states which were unreachable in symmetric cases.viljas  Furthermore, we calculated the spin currents, and related their behavior to transitions between different order-parameter configurations in the aperture. We also studied AA and AB junctions, which have previously not been studied in any depth, at least not within the GL theory for large apertures. We analyzed their symmetries, and related transitions in the $`\widehat{𝐥}`$ field to changes in forms of the CPR’s and their critical currents. In AB junctions, the usual phase-inversion symmetry was shown to be broken, which results in more complicated CPR’s. In addition to all these, a hydrostatic theory of the asymptotic phase corrections is presented in Appendix A, and an analysis of the numerical relaxation method is given in Appendix B.

We have also carried out some tests with the channel geometry used in Ref. ullahfetter, , and again we found the critical currents to match with those in Fig. 14. Also a geometry with the weak link between two *parallel* flow channels was tested. It was found that the critical currents of a BB weak link are not affected by the parallel flow before the flow itself goes unstable by nucleating a vortex at some wall. However, these studies were not carried much further. What we have also not done above is an analysis of the effects of a strong magnetic field. The required magnetic GL terms would be relatively easy to include in the calculation, but they would add more dimensions to the already large parameter space. Actually there remain even some untested phase combinations which may be possible for junctions with small $`W/\xi _{GL}>3`$ and $`DW`$: a BAB or ABA heterostructure. Such situations may be obtainable at pressures close to $`p_0`$. This is because the A phase tends to be slightly more stable in restricted “parallel plate” geometries, and on the other hand because there are some hysteresis effects related to the equilibrium position of the A-B interface as the pressure is varied.viljasab  The BAB configuration was recently discussed as an explanation of some dissipation effects,dissipation  and this was, apparently, the initial reason for considering the BAAB structure in Ref. nishida, . Although it should be possible to prepare such configurations in our program in some pressure regimes, their stability in a real experiment is questionable. Of course, the A phase could be stabilized more strongly inside the aperture by using a localized magnetic field. Also, since a self-consistent, general-temperature calculation of a finite-size aperture (with strong-coupling effects included) is still missing, no final conclusions can be drawn at this stage. Some of these issues may be worth further studies at some later time.

## ACKNOWLEDGMENTS

The Center for Scientific Computing (CSC) is acknowledged for computer resources. A. Schakel, R. Hänninen, and J. Kopu are thanked for helpful discussions.

## Appendix A ASYMPTOTIC SOLUTIONS IN 2D

In this appendix we discuss the approximations used in the asymptotic regions between $`R_c`$ and $`R_{\mathrm{}}`$ of Fig. 1. The need for these derives from a requirement to keep the numerical computation at a minimum, but at the same time have the bulk boundary conditions imposed as far away from the junction as possible. In 3D the boundary conditions ($`R_{\mathrm{}}`$) could be taken all the way to infinity using this method, but in 2D some smaller value has to be chosen. Since in this case $`R_{\mathrm{}}`$ is rather arbitrary, this procedure is to some extent only cosmetic, enabling a smoother transition of the order parameter at the $`R_c`$ cutoff. Thus, whenever this was not needed (like in evaluating the critical currents), the asymptotic corrections were not used.

At $`R_{\mathrm{}}`$ the boundary conditions are given by order parameters of the form $`A_{\mu i}^{L,R}=R_{\mu \nu }^{L,R}A_{\nu i}^{(0)L,R}(x)e^{\mathrm{i}\varphi ^{L,R}}`$ as shown above. However, since $`R`$ and $`\varphi `$ are constant, these states carry no mass or spin currents, apart from the spontaneous spin currents parallel to B-phase surfaces.zhangwall  The order parameters must be corrected in the asymptotic regions in order for them to conserve the currents coming from the weak link aperture. In this respect, the A and B phases can again be treated in the same way, and the asymptotics on each side can be assumed to take the form:

$$A(𝐱)=\stackrel{~}{R}^s(𝐱)RA^{(0)}(x)[\stackrel{~}{R}^o(𝐱)]^Te^{\mathrm{i}(\varphi +\delta \varphi (𝐱))}.$$
(20)
Here $`\delta \varphi (𝐱)`$ is a small phase correction and $`\stackrel{~}{R}^s(𝐱)=R(\delta \theta _\alpha (𝐱))`$ a small spin rotation: $`R_{\mu \nu }^s=\delta _{\mu \nu }ϵ_{\mu \nu \alpha }\delta \theta _\alpha `$. In addition we have included a small orbital rotation $`\stackrel{~}{R}^o(𝐱)=R(\delta ϵ_i(𝐱))`$. This is only needed in the anisotropic A phase to include corrections to the $`\widehat{𝐥}`$ texture, which are present even in the absence of currents. In deriving the asymptotic corrections to $`\delta \varphi `$, $`\delta \theta _\alpha `$, $`\delta ϵ_i`$ and the energy, we shall make the hydrodynamic approximation that $`A^{(0)}`$ is everywhere in its bulk form. In the B phase this means that $`A_{\mu i}^{(0)}=\mathrm{\Delta }_B\delta _{\mu i}`$ and for the A phase $`A_{\mu i}^{(0)}=\mathrm{\Delta }_A\delta _{\mu x}(\delta _{iy}\pm \mathrm{i}\delta _{iz})`$. For convenience we shall denote below the small-angle fields $`\delta \varphi `$, $`\delta \theta _\alpha `$, $`\delta ϵ_i`$ with $`\varphi `$, $`\theta _\alpha `$, $`ϵ_i`$, respectively.

In this context we refer to all positions with respect to the origin $`(x,y)=(D/2,0)`$. As a first approximation, all the small-angle corrections will be assumed to depend only on the radial distance $`r=\sqrt{x^2+y^2}`$ from this point, possibly with some simple dependence on the azimuthal angle $`\phi =\mathrm{arctan}(y/x)`$ also.

### A.1 Asymptotics in BW State

We first deal with the case of B phase which is numerically simpler to handle than the A phase, because it is isotropic (meaning that it has an isotropic superfluid density tensor). There is no need to worry about similar symmetry breaking effects as with the A phase $`\widehat{𝐥}`$ vector (see below), and $`\stackrel{~}{R}^o=1`$. On the other hand, the “entanglement” of the spin and orbital parts of the order parameter makes their separation impossible, which complicates the analytic treatment. To facilitate the analysis of the gradient part of Eq. (8), we first change to a spin basis where the constant $`R`$ in Eq. (20) is replaced by a unit matrix. This is done by noting that $`\stackrel{~}{R}^sR=R\stackrel{~}{R}^s^{}`$, with $`\stackrel{~}{R}_{\mu \nu }^s^{}=\delta _{\mu \nu }ϵ_{\mu \nu \alpha }\theta _\alpha ^{}`$, and $`\theta _\alpha ^{}=R_{\alpha \beta }^T\theta _\beta `$, where the identity $`ϵ_{ijk}R_{jl}R_{km}=R_{in}ϵ_{nlm}`$ is useful. The spin current in this new basis is obtained with $`j_{\alpha i}^{\mathrm{spin}}{}_{}{}^{}=R_{\alpha \beta }^Tj_{\beta i}^{\mathrm{spin}}`$. However, to simplify notation we drop the primes below, remembering that they should appear on all quantities with Greek spin indices. With these definitions, the gradient energy (per unit length) can be writtencross

$$F_G^B=\frac{1}{2}\mathrm{d}^2x\left[\rho _s𝐯_s^2+\rho _{\alpha \beta ;ij}^{\mathrm{spin}}v_{\alpha i}^{\mathrm{spin}}v_{\beta j}^{\mathrm{spin}}\right],$$
(21)
and the mass and spin currents are

$`𝐣_s`$ $`=\rho _s𝐯_s`$ (22)
$`j_{\alpha i}^{\mathrm{spin}}`$ $`={\displaystyle \frac{\mathrm{}}{2m_3}}\rho _{\alpha \beta ;ij}^{\mathrm{spin}}v_{\beta j}^{\mathrm{spin}},`$ (23)
where $`\rho _s=(2m_3/\mathrm{})^22(\gamma +2)K\mathrm{\Delta }_B^2`$, $`\rho _{\alpha \beta ;ij}^{\mathrm{spin}}`$ $`=\rho _s[(\gamma +1)\delta _{ij}\delta _{\alpha \beta }`$ $`(\gamma 2)\delta _{\beta i}\delta _{\alpha j}`$ $`\delta _{\alpha i}\delta _{\beta j}]`$, $`𝐯_s=(\mathrm{}/2m_3)\mathbf{}\varphi `$ and $`𝐯_\alpha ^{\mathrm{spin}}=(\mathrm{}/2m_3)\mathbf{}\theta _\alpha `$. By making a variation in $`\varphi `$ and $`\theta _\alpha `$ and noting the symmetry $`\rho _{\beta \alpha ;ji}^{\mathrm{spin}}=\rho _{\alpha \beta ;ij}^{\mathrm{spin}}`$ we see that the extrema of $`F_G^B`$ are also found from the continuity equations $`\mathbf{}𝐣_s`$ $`=\rho _s\mathbf{}𝐯_s`$ $`=(\mathrm{}\rho _s/2m_3)\mathbf{}^2\varphi =0`$ and $`\mathbf{}𝐣_\alpha ^{\mathrm{spin}}=`$ $`_i\rho _{\alpha \beta ;ij}^{\mathrm{spin}}v_{\beta j}^{\mathrm{spin}}`$ $`=(\mathrm{}\rho _s/2m_3)[(\gamma +1)\mathbf{}^2\theta _\alpha `$ $`(\gamma 1)_\alpha _i\theta _i]=0`$. The latter set of equations can be directly diagonalized,vollhardt  but we assume the solutions to be radial right from the beginning: $`\varphi (r,\phi )=\stackrel{~}{\varphi }(r)`$ and $`\theta _\alpha (r,\phi )=\stackrel{~}{\theta }_\alpha (r)`$. Integration over angles in the energy \[Eq. (21)\] and minimization with respect to these independent phase fields yields, for example, the solution

$$\varphi (r,\phi )=\stackrel{~}{\varphi }(r)=\varphi ^c\frac{\mathrm{ln}(r/R_{\mathrm{}})}{\mathrm{ln}(R_c/R_{\mathrm{}})}.$$
(24)
These satisfy $`\stackrel{~}{\varphi }(R_c)=\varphi ^c`$ and $`\stackrel{~}{\varphi }(R_{\mathrm{}})=0`$. Similar expressions are obtained for the spin phases $`\theta _\alpha (r,\phi )`$. In the $`\rho _s`$ units defined in Appendix B, the energy corresponding to the solutions of the type of Eq. (24) is

$`\stackrel{~}{F}_G^B`$ $`={\displaystyle \frac{2\pi }{4(\gamma +2)}}{\displaystyle \frac{1}{\mathrm{ln}(R_{\mathrm{}}/R_c)}}\{(\gamma +2)(\varphi ^c)^2`$
$`+{\displaystyle \underset{\alpha }{}}(\gamma +1)[1{\displaystyle \frac{1}{2}}{\displaystyle \frac{\gamma 1}{\gamma +1}}(1\delta _{\alpha z})](\theta _\alpha ^c)^2\}`$ (25)
and the asymptotic currents are

$`\stackrel{~}{J}_s`$ $`=\pi \mathrm{ln}(R_{\mathrm{}}/R_c)^1\varphi ^c`$ (26)
$`\stackrel{~}{J}_\alpha ^{\mathrm{spin}}`$ $`=\pi (\gamma +1)\mathrm{ln}(R_{\mathrm{}}/R_c)^1`$
$`\times \left[1{\displaystyle \frac{1}{2}}{\displaystyle \frac{\gamma 1}{\gamma +1}}(1\delta _{\alpha z})\right]\theta _\alpha ^c,`$ (27)
where no summation over $`\alpha `$ is implied. By fitting the numerically calculated currents $`\stackrel{~}{J}_s`$, $`\stackrel{~}{J}_\alpha ^{\mathrm{spin}}`$ (transformed with $`R^T`$to the “primed” spin basis) to Eqs. (26) and (27), one can solve the parameters $`\varphi ^c,\theta _\alpha ^c`$. Then one can use Eq. (24) and the equivalent expression for $`\theta _\alpha `$ in Eq. (20) to update the the asymptotic order parameter, remembering first to transform $`\theta _\alpha ^c`$ back into the original “unprimed” spin coordinates. There is an asymmetry in the spin parts of Eqs. (25) and (27): the $`z`$ direction is slightly more rigid than the others. This results from the orbital-space (real-space) asymmetry of the 2D problem, which also affects the spin space in B phase. In the A phase this is not so, since there the spin and orbital parts of the order parameter are decoupled; see below.

### A.2 Asymptotics in ABM State

In the A phase we assume the following conditions for the vector $`\widehat{𝐥}=\widehat{𝐦}\times \widehat{𝐧}`$: (i) $`\widehat{𝐥}\widehat{𝐰}=0`$, where the orbital rotation axis $`\widehat{𝐰}`$ is constant and in $`yz`$ plane, so that (ii) $`\widehat{𝐰}\widehat{𝐳}1`$, and (iii) $`\widehat{𝐥}\widehat{𝐱}=\mathrm{cos}ϵ`$ where the orbital rotation angle $`ϵ`$ is small. ¿From condition (i) it follows that the phase $`\varphi `$ is well defined and acts as the superfluid velocity potential. Condition (ii) ensures that $`\widehat{𝐥}`$ is essentially in the $`xy`$ plane with $`𝐯_s`$, and (iii) simply means that the variations from $`\widehat{𝐥}=\pm \widehat{𝐱}`$ are small. Using a linearized approximation for $`v_{\alpha i}^{\mathrm{spin}}`$ the A phase gradient energy following from Eq. (8) can then be written

$$\begin{array}{cc}\hfill F_G^A& =\frac{1}{2}\mathrm{d}^2x[\rho _{ij}v_{\mathrm{s}i}v_{\mathrm{s}j}+\rho _{\alpha \beta ;ij}^{\mathrm{spin}}v_{\alpha i}^{\mathrm{spin}}v_{\beta j}^{\mathrm{spin}}\hfill \\ & +(\mathrm{}/2m_3)v_{\mathrm{s}i}C_{ij}(\mathbf{}\times \widehat{𝐥})_j+K_s(\mathbf{}\widehat{𝐥})^2\hfill \\ & +K_b|\widehat{𝐥}\times (\mathbf{}\times \widehat{𝐥})|^2+K_t(\widehat{𝐥}\mathbf{}\times \widehat{𝐥})^2]\hfill \end{array}$$
(28)
and the mass and spin currents are

$`j_{\mathrm{s}i}`$ $`=\rho _{ij}v_{\mathrm{s}j}+(\mathrm{}/2m_3)C_{ij}(\mathbf{}\times \widehat{𝐥})_j`$ (29)
$`j_{\alpha i}^{\mathrm{spin}}`$ $`={\displaystyle \frac{\mathrm{}}{2m_3}}\rho _{\alpha \beta ;ij}^{\mathrm{spin}}v_{\beta j}^{\mathrm{spin}},`$ (30)
where $`\rho _{ij}=\rho _{}\delta _{ij}(\rho _{}\rho _{})\widehat{l}_i\widehat{l}_j`$, $`\rho _{}\rho _{}=(\gamma 1)/(\gamma +1)`$, $`\rho _{}=(2m_3/\mathrm{})^22(\gamma +1)K\mathrm{\Delta }_A^2`$, $`\rho _{\alpha \beta ;ij}^{\mathrm{spin}}=\delta _{\alpha \beta }\rho _{ij}`$ $`𝐯_s=(\mathrm{}/2m_3)\mathbf{}\varphi `$, $`𝐯_\alpha ^{\mathrm{spin}}=`$ $`(\mathrm{}/2m_3)(\widehat{d}_\alpha ^{\mathrm{}}\widehat{d}_\beta ^{\mathrm{}}`$ $`\delta _{\alpha \beta })\mathbf{}\theta _\beta `$, $`K_s=K_t=(\mathrm{}/2m_3)^2\rho _{}/(\gamma +1)`$, $`K_b=(\mathrm{}/2m_3)^2\rho _{}\gamma /(\gamma +1)`$ and $`C_{ij}=C_{}\delta _{ij}(C_{}C_{})\widehat{l}_i\widehat{l}_j`$. Here $`\widehat{d}_\alpha ^{\mathrm{}}=R_{\alpha x}`$ is the bulk spin vector. The first three terms of Eq. (28) are just $`𝐯_\mathrm{s}𝐣_\mathrm{s}/2`$ and $`𝐯_\alpha ^{\mathrm{spin}}𝐣_\alpha ^{\mathrm{spin}}/2`$, but there are three additional terms resulting from the $`\widehat{𝐥}`$ texture alone. Here the actual value of the $`C`$ tensor will not be important to us, since with our assumptions the current component resulting from $`\mathbf{}\times \widehat{𝐥}`$ is either divergenceless ($`\widehat{𝐥}`$ in plane) or even vanishes ($`\widehat{𝐥}`$ constant). As a result, we shall neglect the $`C`$ terms also from the free energy, but keep all others initially.

We now seek to minimize the energy \[Eq. (28)\] with respect to the small-angle fields $`\varphi `$, $`\theta _\alpha `$ and $`ϵ`$. In principle these should be optimized simultaneously, but we shall separate the problems so that in optimizing $`\varphi `$ or $`\theta _\alpha `$ we assume $`\widehat{𝐥}=\pm \widehat{𝐱}`$, and in correcting $`\widehat{𝐥}`$ with $`ϵ`$ we assume that the velocities vanish. Assume, then, that we have $`\widehat{𝐥}=\pm \widehat{𝐱}`$, in which case the last three terms in Eq. (28) vanish. By making a variation in $`\varphi `$ and $`\theta _\alpha `$ and noting the symmetry $`\rho _{\beta \alpha ;ji}^{\mathrm{spin}}=\rho _{\alpha \beta ;ij}^{\mathrm{spin}}`$ we see that the extrema of the first two terms of $`F_G^A`$ are found simply from the continuity equations $`\mathbf{}𝐣_s=_i\rho _{ij}v_{sj}=0`$ and $`\mathbf{}𝐣_\alpha ^{\mathrm{spin}}=`$ $`_i\rho _{\alpha \beta ;ij}^{\mathrm{spin}}v_{\beta j}^{\mathrm{spin}}=0`$. Now if we change to new coordinates $`x^{}=x\sqrt{(\gamma +1)/2}`$, $`y^{}=y`$, then the superfluid density tensor $`\rho _{ij}^{}`$ will appear isotropic. In these coordinates the continuity equations are simply the Laplace equations $`\mathbf{}_{}^{\mathbf{}}{}_{}{}^{2}\varphi =\mathbf{}_{}^{\mathbf{}}{}_{}{}^{2}\theta _\alpha =0`$. We attempt to solve these in polar coordinates $`r^{},\phi ^{}`$ of the primed system, where the general solution is of the form

$`\varphi (r^{},\phi ^{})`$ $`=A\phi ^{}+B\mathrm{ln}(r^{}/a)`$
$`+`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(C_nr_{}^{}{}_{}{}^{n}+D_nr_{}^{}{}_{}{}^{n})\mathrm{sin}(n\phi ^{}\alpha _n).`$ (31)
The solution we are looking for is determined by setting boundary conditions on the $`R_c`$ cutoff and at infinity. We assume that the currents coming from the junction flow radially as $`r^{}\mathrm{}`$. This is described by the solution proportional to $`\mathrm{ln}r^{}`$. However, because the $`\widehat{𝐥}`$ texture tends to align perpendicular to surfaces, it usually becomes nonsymmetrical close to the weak link (see Section 5). Since the superfluid density is larger perpendicular to $`\widehat{𝐥}`$, most of the current will flow perpendicular to it. Therefore, the currents do not flow in and out of the junction mirror-symmetrically with respect to the $`xz`$ plane.

Such a flow field can be thought to consist of two components: one which is symmetrical and radial, and another which circulates in front of the opening as illustrated in Fig. 16. For the circulating part of the mass current we test two boundary conditions: $`_r^{}\varphi (R_c)\mathrm{sin}\phi ^{}`$ and $`_r^{}\varphi (R_c)\mathrm{sin}2\phi ^{}`$. Both of these describe a flow which is outward above the $`x`$ axis when it is inward below it and vice versa, but the distributions differ. The corresponding solutions are $`\varphi =\mathrm{sin}\phi ^{}/r^{}`$ and $`\varphi =\mathrm{sin}2\phi ^{}/r_{}^{}{}_{}{}^{2}`$ and the total velocity potential is thus of the form

$$\varphi (r^{},\phi ^{})=\varphi ^c\frac{\mathrm{ln}(r^{}/R_{\mathrm{}})}{\mathrm{ln}(R_c/R_{\mathrm{}})}+\varphi ^{c,a}\frac{R_c^n}{r_{}^{}{}_{}{}^{n}}\mathrm{sin}(n\phi ^{}),n=1,2.$$
(32)
Here the coefficient $`\varphi ^c`$ is again determined by equating the current obtained from Eq. (29) by integration over $`R_c`$ (surface EDF in Fig. 16) with the total current $`J_s`$ calculated numerically at the junction (surface AB). The coefficient $`\varphi ^{c,a}`$, on the other hand, is trickier. It could be determined by equating the current resulting from the second term of Eq. (32) over the upper part of the cutoff (surface DE) with the $`y`$-directional current $`J_s^a`$ calculated numerically over the surface CD. Note that in the case $`n=1`$ the circulating parts of the currents over surfaces DE, D$`\mathrm{}`$ and DF are equal, and hence no net $`y`$ directional current is introduced. This is not so for $`n=2`$, and, in fact, in this case the phase field gives an unphysical-looking phase gradient perpendicular to the wall. However, our approximation does not take into account the fact that the order parameter is suppressed at the wall and (at least for a diffusive surface) the superfluid density vanishes there. Therefore no *current* should flow through it. The phase correction of Eq. (32) should only be seen as a variational ansatz, and in practice we find that $`n=2`$ gives a better overall fit with lower total energies than $`n=1`$. Similar expressions and conclusions exist for the spin velocity potentials $`\theta _\alpha `$.

We now give the results for the current vs. phase correction relationships. Let us remind that in deriving these we have assumed that $`\widehat{𝐥}=\pm \widehat{𝐱}`$ in the asymptotic region. In the $`\rho _{}`$ units defined in Appendix B, the gradient energy of the asymptotic solution is

$`\stackrel{~}{F}_G^A=`$ $`{\displaystyle \frac{1}{2}}\pi \sqrt{{\displaystyle \frac{2}{\gamma +1}}}\{{\displaystyle \frac{(\varphi ^c)^2}{\mathrm{ln}(R_{\mathrm{}}/R_c)}}+{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{(\theta _\alpha ^c)^2}{\mathrm{ln}(R_{\mathrm{}}/R_c)}}`$
$`+{\displaystyle \frac{1}{2n}}[1(R_c/R_{\mathrm{}})^{2n}](\varphi ^{c,a})^2`$
$`+{\displaystyle \underset{\alpha }{}}{\displaystyle \frac{1}{2n}}[1(R_c/R_{\mathrm{}})^{2n}](\theta _\alpha ^{c,a})^2\}`$ (33)
while the currents are

$`\stackrel{~}{J}_s`$ $`=\pi \sqrt{2/(\gamma +1)}\mathrm{ln}(R_{\mathrm{}}/R_c)^1\varphi ^c`$ (34)
$`\stackrel{~}{J}_\alpha ^{\mathrm{spin}}`$ $`=\pi \sqrt{2/(\gamma +1)}\mathrm{ln}(R_{\mathrm{}}/R_c)^1\theta _\alpha ^c`$ (35)
$`\stackrel{~}{J}_s^a`$ $`=n\sqrt{2/(\gamma +1)}\varphi ^{c,a}`$ (36)
$`\stackrel{~}{J}_\alpha ^{\mathrm{spin},a}`$ $`=n\sqrt{2/(\gamma +1)}\theta _\alpha ^{c,a}.`$ (37)
Again, by fitting the numerically calculated currents to these expressions, one can solve the parameters $`\varphi ^c`$, $`\theta _\alpha ^c`$. The above methods of estimating the “asymmetric” currents are very unreliable, and in practice the corrections $`\varphi ^{c,a}`$ and $`\theta _\alpha ^{c,a}`$ were not used. They are presented here for completeness.

In addition there is the small correction to the direction of $`\widehat{𝐥}`$, which is mostly due to the complicated surface structure. There is also a small tendency for the $`\widehat{𝐥}`$ to align with $`\pm 𝐯_s`$ at high superflow velocities, but here we assume $`𝐯_s=𝐯_\alpha ^{\mathrm{spin}}=0`$. We do not require $`\widehat{𝐥}`$ to vary only in $`xy`$ plane, as long as it is rotates only around the fixed axis $`\widehat{𝐰}`$ which is in the $`yz`$ plane. The minimization of the last two terms of Eq. (28) is done by assuming that $`\widehat{𝐥}(r^{},\phi ^{})=\pm R(\widehat{𝐰},ϵ(r^{},\phi ^{}))\widehat{𝐱}`$, so that $`\widehat{𝐥}=\pm \widehat{𝐱}`$ for $`r^{}\mathrm{}`$ but also at the walls ($`\phi ^{}=\pm \pi /2`$). This can be described with the variational ansatz $`ϵ(r^{},\vartheta ^{})=\stackrel{~}{ϵ}(r^{})\mathrm{cos}\phi ^{}`$. If we further assume that $`(K_{s,t}+K_b)/2K_b(\gamma +1)/2`$ in order to drop some terms with $`y`$ and $`z`$ derivatives, the energy density resulting from the bending of $`\widehat{𝐥}`$ is just $`2f_g^l=K_b[(_xϵ_y)^2+(_xϵ_z)^2]`$, where $`\mathit{ϵ}=ϵ\widehat{𝐰}`$. Minimization of this with the above ansatz gives $`\stackrel{~}{ϵ}(r^{})\stackrel{~}{ϵ}^c(R_c/r^{})`$, and the energy

$$\mathrm{\Delta }\stackrel{~}{F}_G^A=\frac{3\pi \gamma }{16}\sqrt{\frac{2}{\gamma +1}}[1(R_c/R_{\mathrm{}})^2](\stackrel{~}{ϵ}^c)^2.$$
(38)
Here $`(\stackrel{~}{ϵ}^c)^2=(\stackrel{~}{ϵ}_y^c)^2+(\stackrel{~}{ϵ}_z^c)^2`$, and $`\stackrel{~}{ϵ}_{y,z}^c\widehat{l}_{y,z}^c(x^{}=R_c,y^{}=0)`$, These cannot be obtained from any conservation law, but instead we extract them from the numerical solution by requiring continuity of $`\widehat{l}_y`$ and $`\widehat{l}_z`$ and their $`x`$ derivatives along the $`x`$ axis. The $`\widehat{𝐥}`$ vector can be obtained numerically from the $`A`$ matrix inside the $`R_c`$ region with $`\widehat{l}_i=ϵ_{ijk}Im(\stackrel{~}{A}_{\mu j}^{}\stackrel{~}{A}_{\mu k})/2`$, as shown in Fig. 8.

## Appendix B NUMERICAL METHODS

This appendix considers in some detail the implementation of our numerical relaxation methods. We skip the easier part of calculating the 1D order parameter $`A^{(0)}(x)`$ for a planar wall, since it is only needed as a boundary condition for the 2D calculation. Of course, the same techniques can be used for that also, but in practice the number of variables in the 1D minimization is so small that perfecting the method is not of vital importance. At all solid surfaces, we only consider the diffuse-scattering boundary condition $`A_{\mu i}=0`$, see Section 3. Below we assume that the boundary conditions on the minimization region are fixed, and comment on the addition of the asymptotic corrections (Appendix A) in the last subsection. For the purpose of numerics, we must consider all quantities to be dimensionless, and we first discuss the unit reductions.

### B.1 Units

A natural way to scale any quantity $`L`$ with the units of length is $`\stackrel{~}{L}=L/\xi _{GL}`$, and a dimensionless order parameter is obtained with $`\stackrel{~}{A}=A/\mathrm{\Delta }_{A,B}`$, where $`\mathrm{\Delta }_{A,B}`$ are the bulk gaps of A and B phases, respectively. In the A phase $`\mathrm{\Delta }_A^2=\alpha /[2(2\beta _{245})]`$, where $`\beta _{245}=\beta _2+\beta _4+\beta _5`$, and for the B phase $`\mathrm{\Delta }_B^2=\alpha /[2(3\beta _{12}+\beta _{345})]`$, where $`\beta _{12}=\beta _1+\beta _2`$, $`\beta _{345}=\beta _3+\beta _4+\beta _5`$. Physically motivated units for the energies and currents are based on the the bulk superfluid densities $`\rho _{}=(2m_3/\mathrm{})^22(\gamma +1)K\mathrm{\Delta }_A^2`$ and $`\rho _s=(2m_3/\mathrm{})^22(\gamma +2)K\mathrm{\Delta }_B^2`$ (see appendix A). Using these we define the following units for the currents and energy (per unit length). In A phase $`J_A=(\mathrm{}/2m_3)\rho _{}`$, $`J_A^{\mathrm{spin}}=(\mathrm{}/2m_3)^2\rho _{}`$, and $`E_A=(\mathrm{}/2m_3)^2\rho _{}`$. For B phase $`J_B=(\mathrm{}/2m_3)\rho _s`$, $`J_B^{\mathrm{spin}}=(\mathrm{}/2m_3)^2\rho _s`$, and $`E_B=(\mathrm{}/2m_3)^2\rho _s`$. It is these units which are used in all the figures of this paper, and sometimes we denote $`\stackrel{~}{J}_s=J_s/J_{A,B}`$, $`\stackrel{~}{J}_\alpha ^{\mathrm{spin}}=J^{\mathrm{spin}}/J_{A,B}^{\mathrm{spin}}`$, and $`\stackrel{~}{F}=F/E_{A,B}`$. Note that for weak coupling (see Section 3), the A and B phase superfluid densities $`\rho _s`$ and $`\rho _{}`$ are equal, and so are the units. For cases where both A and B phase are involved simultaneously, we choose to use the B phase units.

However, for numerics, a more convenient reduction of the energy is given by $`\overline{f}=f/[\alpha \mathrm{\Delta }_{A,B}^2/2]`$, and then the reduced GL parameters are given by $`\overline{\beta }_i=(2\mathrm{\Delta }_{A,B}^2/\alpha )\beta _i`$, $`\overline{\alpha }=2`$ and $`\overline{K}=2`$. Using the $`\mathrm{\Delta }_{A,B}`$ we thus have in the A phase $`\overline{\beta }_i=\beta _i/(2\beta _{245})`$ and in the B phase $`\overline{\beta }_i=\beta _i/(3\beta _{12}+\beta _{245})`$. In the formulas below, these unit reductions are assumed to have been performed, but we do not write any tildes or overlines on the symbols.

### B.2 The Minimization Problem

The 2D numerical scheme is based on an optimization of the free energy $`F_\mathrm{\Omega }`$ \[Eq. (3)\] or, which is more-or-less equivalent, solution of the corresponding nonlinear Euler-Lagrange field equations. In fact, the method to be discussed can be applied (at least approximately) for solving a more general class of nonlinear equations which have the form $`G_{\mu i}(A(𝐱))=0`$, where $`\mu ,i=x,y,z`$. Nevertheless, we shall consider the the problem as that of optimizing (minimizing) the GL energy functional $`F_\mathrm{\Omega }[A]`$ in the region $`\mathrm{\Omega }`$ of Fig. 1. More mathematically speaking, for the given fixed boundary conditions $`A^{L,R}`$ \[of the form of Eq. (20)\] on $`\mathrm{\Gamma }^{L,R}`$, we wish to find $`F_J[A^L,A^R]`$ $`=\mathrm{min}_AF_\mathrm{\Omega }[A]`$. As discussed in Appendix A, the full numerical minimization is done only inside the region bounded by the $`R_c`$ cutoffs, but the computational region can be effectively increased by making use of the asymptotic corrections between $`R_c`$ and $`R_{\mathrm{}}`$.

In the following, all entities $`A`$ with the $`3\times 3`$ matrix structure $`A_{\mu i}`$ where $`\mu ,i=x,y,z`$ will be referred to as *spin-orbit matrices*, or *spin-orbit fields*. For simplicity of notation, let us define the following rotational invariants which are functions of the spin-orbit matrices $`A,B,C,D`$:

$`f_\alpha (A,B)`$ $`=ReTr(AB^\mathrm{T}),`$ (39)
$`f_1(A,B,C,D)`$ $`=ReTr(A^{}B^\mathrm{T})Tr(CD^\mathrm{T}),`$ (40)
$`f_2(A,B,C,D)`$ $`=ReTr(AB^\mathrm{T})Tr(CD^\mathrm{T}),`$ (41)
$`f_3(A,B,C,D)`$ $`=ReTr(AB^\mathrm{T}C^{}D^\mathrm{T}),`$ (42)
$`f_4(A,B,C,D)`$ $`=ReTr(AB^\mathrm{T}CD^\mathrm{T}),`$ (43)
$`f_5(A,B,C,D)`$ $`=ReTr(AB^\mathrm{T}C^{}D^\mathrm{T}),`$ (44)
$`g_1(A,B)`$ $`=Re_iA_{\mu i}^{}_jB_{\mu j},`$
$`g_2(A,B)`$ $`=Re_iA_{\mu j}^{}_iB_{\mu j},`$
$`g_3(A,B)`$ $`=Re_iA_{\mu j}^{}_jB_{\mu i}.`$ (45)
The invariants have several symmetries, some of which are listed below:

$`f_\alpha (A,B)`$ $`=f_\alpha (B,A)=f_\alpha (A^{},B^{}),`$ (46)
$`f_{1,2}(A,B,C,D)`$ $`=f_{1,2}(A,B,D,C),`$ (47)
$`f_{3,4,5}(A,B,C,D)`$ $`=f_{3,4,5}(B,A,D,C),`$ (48)
$`f_{15}(A,B,C,D)`$ $`=f_{15}(C,D,A,B),`$ (49)
$`f_{15}(A^{},B^{},C,D)`$ $`=f_{15}(A,B,C^{},D^{}),`$ (50)
$`g_{1,2,3}(A,B)=g_{1,2,3}`$ $`(B,A)=g_{1,2,3}(A^{},B^{}).`$ (51)
It is worth noting here that symmetry of Eq. (51) of the gradient terms is not necessarily valid in a discretized form, unless symmetric difference approximations are used for all derivatives. Using the definitions of Eqs. (39-45), the GL energy functional \[Eqs. (3) and (8)\] can now be rewritten as follows

$`F[A,A^{}]=`$ $`{\displaystyle }\mathrm{d}^3x[\alpha f_\alpha (A,A)+{\displaystyle \underset{i=1}{\overset{3}{}}}K_ig_i(A,A)`$
$`+{\displaystyle \underset{i=1}{\overset{5}{}}}\beta _if_i(A,A,A,A)].`$ (52)
For the purpose of numerics, it is very helpful to consider $`F`$ as a functional of both $`A`$ and $`A^{}`$ independently, as shown here explicitly. Note that the third gradient term in (45) can be integrated by parts and added to the first term, if a separate surface energy term is introduced. Then, if the surface and the associated boundary conditions are symmetric enough, the surface term may vanish altogether, thereby reducing the numerical effort somewhat. Since this is not always the case, we shall not assume it anywhere.

### B.3 Line Searching in Minimization

The basic component of any efficient minimization algorithm is a method for doing *line searches*, *i.e.*, doing one-dimensional minimizations of the function in some given direction. More precisely, one wishes to minimize

$$h(\lambda )=F[A+\lambda D,A^{}+\lambda D^{}]$$
(53)
with respect to the real parameter $`\lambda `$, given the “starting point” $`A`$ and a “search direction” $`D`$. The extrema of $`h(\lambda )`$ are found from the zeros of its derivatives, *i.e.*, by solving

$$h^{}(\lambda )=2Re(G(A+\lambda D),D)=0,$$
(54)
where $`G(A(𝐱))=\delta F[A,A^{}]/\delta A^{}(𝐱)`$, and we have defined a *scalar product* of two spin-orbit fields $`A`$ and $`B`$ as follows:

$$(A,B)=\mathrm{d}^3xTr(A(𝐱)B^\mathrm{T}(𝐱)).$$
(55)
It is easy to check that this is indeed bilinear, satisfies $`(B,A)=(A,B)^{}`$, and that $`(A,A)0`$, with $`(A,A)=0`$ only for $`A=0`$. The scalar product also gives a natural definition of a *norm* for the spin-orbit field: $`A=\sqrt{(A,A)}`$. Note that, with respect to the matrix indices, this is the so-called *Frobenius norm*.

The interesting point is that, since the GL free energy is only of fourth order (quartic) in the order parameter, the function $`h(\lambda )`$ is simply a fourth order polynomial with real coefficients

$$h(\lambda )=a\lambda ^4+b\lambda ^3+c\lambda ^2+d\lambda +e.$$
(56)
The fourth-order coefficient is clearly always positive, so that the graph of the function must look like the curve in Fig. 17.

The number of minima is 1 or 2, but only those with $`\lambda >0`$ will be of interest here. If one is now willing to take the small trouble of calculating the coefficients $`a`$ to $`e`$, then the *exact* minima of $`h(\lambda )`$ are immediately found by solving the real roots of a third-order equation $`h^{}(\lambda )=0`$. For the conjugate gradient method the accuracy of line searches is important, at least in theory, and we list below the formulas for $`a`$ to $`e`$ for the GL free energy discussed above:

$$a=\mathrm{d}^3x[\underset{i=1}{\overset{5}{}}\beta _if_i(D,D,D,D)]$$
(57)
$$b=\mathrm{d}^3x[4\underset{i=1}{\overset{5}{}}\beta _if_i(D,D,D,A)]$$
(58)
$`c`$ $`={\displaystyle }\mathrm{d}^3x\{2{\displaystyle \underset{i=1}{\overset{5}{}}}\beta _i[f_i(D,D,A,A)+f_i(D,A,D,A)`$
$`+f_i(D,A,A,D)]\alpha f_\alpha (D,D)+{\displaystyle \underset{i=1}{\overset{3}{}}}K_ig_i(D,D)\}`$ (59)
$`d=`$ $`{\displaystyle }\mathrm{d}^3x\{4{\displaystyle \underset{i=1}{\overset{5}{}}}\beta _if_i(D,A,A,A)2\alpha f_\alpha (D,A)`$
$`+{\displaystyle \underset{i=1}{\overset{3}{}}}K_i[g_i(D,A)+g_i(A,D)]\}`$ (60)
$$e=F[A,A^{}]$$
(61)
Here the last term of $`d`$ can now be simplified with Eq. (51), but only with the proviso that the difference approximations are symmetric. The coefficient $`e`$ (free energy at position $`A`$) is just a constant and unnecessary for the line searching procedure.

The numerical calculation can be further simplified by assuming the coefficients $`a`$ and $`b`$ to be zero. This is equivalent to a local quadratic approximation of the energy functional, which is reasonable if the center point $`A`$ is already close to the minimum of $`F`$ (see Fig. 17). Note that Eqs. (57-B.3) for the coefficients $`a`$ to $`d`$ can also be obtained directly from Eq. (54), given the GL equations $`G(A)=0`$. This is indeed the only way, if no such functional exists whose Euler-Lagrange equation this $`G(A)=0`$ would be. Of course, in that case there is no guarantee that the CG iteration will converge at all.

### B.4 The Conjugate Gradient Method

Perhaps the simplest approach for the minimization of a nonlinear function(al) is the *steepest descent* (SD) method. In this method, on each ($`k`$th) step of iteration, the line search direction is simply chosen to be that of the negative gradient: $`D_k=G_k`$, where $`G_k=G(A_k)`$. The line search condition \[Eq. (54)\] then implies that $`(D_k,D_{k1})=0`$, *i.e.*, successive search directions are perpendicular to each other with respect to the appropriate scalar product \[Eq. (55)\]. This leads to the well-known zigzag-type trajectory, which in many cases results in extremely slow convergence.bazaraa  The situation can usually be helped somewhat by skipping the costly line search procedure altogether, by choosing $`A_{k+1}=A_kcG_k`$, where $`c`$ is a given small positive constant. This procedure is in fact equivalent to a simple (“forward Euler”) time-discretization of the time-dependent GL (or TDGL) equations, where $`c`$ is proportional to the time step and the strength of dissipation.thuneberg87  By following the progress of this iteration one may thus get a rough idea of the dynamics of the order parameter when, for example, a phase slip in a weak link occurs.

At least in the case of the GL equations, a significant speedup of convergence can be achieved by taking advantage of the knowledge of previous search directions. In the standard *conjugate gradient* (CG) methodsbazaraa  the iteration step $`k`$ is taken as follows:

> Conjugate gradient iteration step
>
> * $`D_k=G_k+\beta _kD_{k1}`$
> * Do a line search of $`h(\lambda )=F[A_k+\lambda D_k,A_k^{}+\lambda D_k^{}]`$ to find its (smallest) positive minimum point $`\stackrel{~}{\lambda }`$
> * Set $`A_{k+1}=A_k+\stackrel{~}{\lambda }D_k`$

The constant $`\beta _k`$ can be calculated in one of several ways. Two of the most commonly used are the *Fletcher-Reeves* choice $`\beta _k=(G_k,G_k)/(G_{k1},G_{k1})`$, and the *Polak-Ribière* choice $`\beta _k=(G_kG_{k1},G_k)/(G_{k1},G_{k1})`$. The latter can be slightly more efficient, but it has the memory-economical disadvantage that the gradient $`G_{k1}`$ of the previous round must also be stored. The method is usually started with $`D_0=G_0`$ which is guaranteed to be a direction of descent. During the iteration the process can also be restarted every once in a while. This is because it is possible that inefficient search directions will begin to be generated after a few rounds of minimization, if the function(al) deviates very strongly from a quadratic form. Restarting should be done at least if the line search fails to find a positive $`\stackrel{~}{\lambda }`$, and possibly on every $`N`$th step or so, where $`N`$ is the total number of real variables in the discretized order-parameter field.bazaraa

It seems that in practice neither the accuracy of the line search, nor the choice for the $`\beta _k`$’s is very important, at least if the method is sometimes restarted. The most important part is really the basic philosophy of using the previous search direction as described above.

One of the benefits of exact line searching is that it provides the small *change* in energy between steps very naturally through $`\mathrm{\Delta }F_k`$ $`=a\stackrel{~}{\lambda }^4+b\stackrel{~}{\lambda }^3+c\stackrel{~}{\lambda }^2+d\stackrel{~}{\lambda }`$, so that $`e`$, the energy itself, does not need to be calculated on every round. This appears to be a more accurate way to obtain the change than application of $`\mathrm{\Delta }F_k=F[A_k,A_k^{}]F[A_{k1},A_{k1}^{}]`$ after completion of each iteration step. This is an important point when the limit of available floating point accuracy is eventually approached: one should always have $`\mathrm{\Delta }F_k<0`$ for a *descent* direction $`D_k`$, and it is wise to check that this is indeed the case. However, a mistake would be done by assuming that $`F[A_n,A_n^{}]=_{k=0}^{n1}\mathrm{\Delta }F_k+F[A_0,A_0^{}]`$, since this will involve summation of numbers whose orders of magnitude are very different from each other. Sufficiently close to the minimum an $`F`$ calculated in this way may no longer change at all during the iteration, even if the the gradient is still nonzero. The change in energy should therefore be used as a judge of convergence only with care.

### B.5 Comparison of Methods

Here we compare briefly the convergence of the different methods. We use an example calculation in the geometry of Fig. 1, with $`R_c/\xi _{GL}=25`$, $`W/\xi _{GL}=10`$, $`D/\xi _{GL}=6`$, with no asymptotic corrections ($`R_{\mathrm{}}=R_c`$), and B phase with $`\widehat{𝝎}^L=\widehat{𝝎}^R`$ on both sides. The discretization lattice spacing was chosen to be $`\mathrm{\Delta }x/\xi _{GL}=\mathrm{\Delta }y/\xi _{GL}=1`$. Figure 18 shows plots of the length of the gradient on the $`k`$th iteration round, $`G_k=G(A_k)`$, versus the round number $`k`$. Both the Fletcher-Reeves and Polak-Ribière variants of the CG method are compared to the SD method. The superiority of either form of CG is evident from this figure: while the CG methods reach the level $`G_k10^4`$ for $`k=200\mathrm{}300`$, the SD method requires almost $`k=4000`$ for it. The real time needed by SD is thus more than tenfold compared to CG. This does not reveal anything rigorous about the convergence classes, but for the CG methods the convergence must be superlinear, and probably close to quadratic.bazaraa  The scaling properties (behavior of convergence with increasing problem size) were also not tested. As an interesting detail, it was noticed that the A phase tended to converge faster than B phase.

Another possibility for obtaining faster convergence compared to the SD algorithm would be to use *quasi-Newton* (secant) methods, which utilize knowledge of the Hessian matrix or some approximation to it.bazaraa  However, in our example problem the CG method has several advantages over these methods also: (i) While analytic calculation of the second derivatives of $`F`$ is possible, the Hessian matrix can easily become huge and take up a lot of memory. This is not so severe in approximate schemes (such as BFGS) which do not store the whole Hessian matrix. (ii) The Newton method has an inconvenient property of being unstable, if the initial guess is not chosen close enough to the minimum. To avoid problems, special stability checks are needed. (iii) Most inconveniently of all, the construction of the Hessian matrix and solving the related linear system of equations would require the programmer to *order* the field variables in some arbitrary way to form a *vector*. In the case of a computational region looking like the one in Fig. 1 and not a rectangle, there is no obvious way to do such ordering. Any choice would require one to program obscure and lengthy routines for transformations between the “$`xy`$” and vector orderings. The CG method requires no such thing, and all field variables can be iterated in any order on each iteration step, independently of each other. (iv) The explicit use of $`A`$ and $`A^{}`$ as the independent variables instead of $`ReA`$ and $`ImA`$ is also simplest in CG type methods. The practical necessity of this is clear from the complicated structure of the free energy functional. Points (iii) and (iv) are also why an application of most ready-made library routines is not very attractive for problems in general $`p`$-wave Ginzburg-Landau theory.

### B.6 Updating the Asymptotic Values

Finally we consider the task of fitting the asymptotic solutions of Appendix A to the numerical solutions obtained by the above methods for $`A^{L,R}`$ fixed on $`r=R_c`$. This could be done in several ways, but here we choose a very simple one, and describe it only in therms of the B phase. The basic idea is iterative. Start with $`\varphi ^c=\theta _\alpha ^c=0`$ in Eq. (24), an thus no phase corrections in the the asymptotic order parameter \[Eq. (20)\]. Then (i) minimize $`F`$ inside the $`R_c`$ cutoffs, (ii) calculate the currents $`\stackrel{~}{J}_s`$ and $`\stackrel{~}{J}_\alpha ^{\mathrm{spin}}`$ in the junction, (iii) insert them in Eqs. (26) and (27) to get new $`\varphi ^c`$ and $`\theta _\alpha ^c`$, (iv) update Eq. (20) accordingly, and then repeat from (i) until satisfactory conservation of the currents over the $`R_c`$ cutoff is achieved. If we denote by “out” the current in the asymptotic region, and $`\mathrm{\Delta }\stackrel{~}{J}_s`$ $`\stackrel{~}{J}_s^{out}\stackrel{~}{J}_s`$ then the requirement was usually $`|\mathrm{\Delta }\stackrel{~}{J}_s|10^4`$, and similarly for the spin currents.

Note that $`\mathrm{\Delta }\stackrel{~}{J}_s`$ $`\stackrel{~}{F}_\mathrm{\Omega }/\varphi ^c`$, where $`\stackrel{~}{F}_\mathrm{\Omega }`$ is the energy in the whole $`\mathrm{\Omega }`$ for given $`\varphi ^c,\theta _\alpha ^c`$, and thus the process (i-iv) is seeking the minimum of $`\stackrel{~}{F}_\mathrm{\Omega }`$. Since this is obviously a very inefficient method, the benefit of fast convergence by the conjugate gradient method in the $`r<R_c`$ part of $`\mathrm{\Omega }`$ is lost if this iteration is continued very far. Fortunately, the reason for taking the corrections into account is mostly cosmetic in practice, and only a few (less than 10) iteration rounds may already lead to a reasonably continuous order parameter and current across the $`R_c`$ cutoff. In the case of A phase, also the correction of $`\widehat{𝐥}`$, and possibly the asymmetric phase corrections, should be updated on each round. See Appendix A for details.