text:
# 1 Introduction

## 1 Introduction

Though there is no analytic proof of confinement in non-supersymmetric nonabelian gauge theories lattice simulations provide ample numerical evidence that nonabelian gauge theories confine quarks and gluons . Considerable effort has been spent to understand the underlying physical mechanism behind confinement. The two most popular models, explaining the physical reasons behind confinement, are the dual superconductor model  and the $`Z_N`$ vortex condensation model    .

Continuum gauge theories require gauge fixing. On the lattice, gauge fixing is not required, but possible. In the dual superconductor model the nonabelian part of the gauge is fixed, while a Cartan subgroup is not. As gauge fixing leaves behind a theory with Abelian gauge symmetry, this is called Abelian gauge fixing. The gauge fixing procedure has defects, which can be interpreted as monopoles. The gauge field, projected to its abelian components, and monopoles form a Coulomb gas that confines charged objects .

In the $`Z_N`$ magnetic vortex model of $`SU(N)`$ gauge theories the gauge is fixed completely, except for the discrete center of the group, $`Z_N`$. This is called center gauge fixing. The defects of this gauge fixing are lines, representing magnetic fluxes, vortices, in three dimensional space. The condensation and percolation of vortices can be easily shown to lead to an area law of Wilson loops, which is tantamount to confinement.

Lattice studies of abelian projected theories   found good agreement between numerical values of gauge invariant quantities (string tension, chiral condensate, etc.) obtained with and without abelian projection. This was assumed to be an indication of the validity of the dual superconductor picture.

Center projected pure lattice gauge theories    are $`Z_N`$ gauge theories. $`Z_N`$ gauge theories are theories of interacting vortices. Averages of Wilson loops and of other gauge invariant quantities calculated in the projected $`Z_N`$ gauge theory were also in good agreement with numerical values of gauge invariant quantities obtained without gauge fixing. This agreement implies that the center vortex picture is also a viable candidate for explaining confinement.

The analytic study of monopoles and magnetic vortices in pure gauge theories is, however, extremely difficult. Among existing studies, Polyakov’s classical Coulomb gas model model , and recent works defining vortices in the continuum  are of this nature. Strong coupling expansion  gives credence to the vortex condensation picture, as well. The analytic study of monopoles and vortices in gauge fixed pure gauge theories is very important. It is conceivable that finding an analytic proof of confinement may require an analytic proof of the existence and condensation of these classical objects in pure gauge theories.

Topological solutions, such as monopoles and vortices, appear in another class of theories under more controlled circumstances. Monopoles were found in Higgs theories with a single adjoint representation Higgs boson   and vortices were found in Higgs theories with two or more Higgs bosons    , as topologically stable classical solutions.

In an earlier work we conjectured that monopoles and vortices found in gauge fixed pure gauge theories are connected with monopoles and vortices appearing as classical solitons in Higgs theories . In this paper we intend to investigate this relationship and establish a firm link among these objects. Using the link between gauge fixed pure gauge theories and of Higgs theories the properties of monopoles and vortices and their interactions can be studied analytically in Higgs theories and the conclusions can be carried over to gauge fixed pure gauge theories.

We will rely on lattice realizations of gauge theories and Higgs theories throughout this paper. We will assume that at sufficiently small gauge coupling they are close to continuum theories. For the sake of simplicity we will only consider $`SU(2)`$ gauge theories, but our conclusions can be generalized to other nonabelian gauge theories without much difficulty.

## 2 Abelian gauge fixing and Higgs theories

We will first consider Abelian gauge fixing procedures in pure lattice gauge theories. The most popular gauges are the Maximal Abelian gauge (MAG) and the Laplacian Abelian gauge (LAG). These are opposite limits of a continuous set of gauge fixing procedures that we will define below. All of these procedures use lattice gauge configurations generated without gauge fixing.

MAG is defined  by maximizing a functional, $`S_{\mathrm{gf}}^V`$, of the gauge field over gauge transformations, $`V(𝒓)`$,

$`F_{\mathrm{gf}}(U^V)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝒓,\mu }{}}\mathrm{Tr}\left[\sigma _3U_\mu ^V(𝒓)\sigma _3U_\mu ^{V,}(𝒓)\right],`$ (1)
where $`U_\mu ^V(𝒓)=V(𝒓)U_\mu (𝒓)V(𝒓+𝝁)`$. $`𝝁`$ represents a lattice vector in the direction of the $`\mu `$th axis.

For future purposes we rewrite $`F_{\mathrm{gf}}(U^V)`$ in terms of quantities defined in the adjoint representation. As the choice of gauge is arbitrary we will choose $`U_\mu (𝒓)`$ as a gauge field at the maximum of $`F_{\mathrm{gf}}(U^V)`$, at $`V(𝒓)=e^{i\gamma _3\sigma _3}`$. Then we write $`S_{\mathrm{gf}}^V`$ as

$`F_{\mathrm{gf}}(U^V)`$ $`=`$ $`{\displaystyle \underset{x,\mu }{}}g_{3a}^A(𝒓)U_\mu ^{A,ab}(𝒓)g_{3b}^A(𝒓+𝝁),`$ (2)
where we defined the adjoint representation gauge transformation as

$$g_{ab}^A(𝒓)=\frac{1}{2}\mathrm{Tr}\left[\sigma _aV(𝒓)\sigma _bV^{}(𝒓)\right]$$
(3)
and the adjoint representation gauge field as

$$U_\mu ^{A,ab}(𝒓)=\frac{1}{2}\mathrm{Tr}\left[\sigma _aU_\mu (𝒓)\sigma _bU_\mu ^{}(𝒓)\right].$$

$`g_{ab}^A`$ and $`U_\mu (𝒓)`$ are orthogonal matrices. Note that the Cartan (Abelian) subgroup of the gauge group leaves the vector $`g_{a3}^A`$ invariant. In $`SU(2)`$ the normalized vector $`g_{a3}^A(𝒓)`$ fixes only two of the parameters of the gauge group.

LAG relaxes the constraint $`_a[g_{a3}^A(𝒓)]^2=1`$ . It maximizes (2) by finding the eigenvector of largest eigenvalue of the quadratic form defined by the right hand side of (2). The local normalization of the real vectors $`g_{a3}(𝒓)`$ then determines the gauge transformation through (3). The conceptual advantage of LAG over MAG is that the trajectories of monopoles can be defined as singularities of the gauge fixing procedure, in agreement with the original ideas of ‘t Hooft . At points where $`g_{a3}(𝒓)`$ vanishes the gauge transformation becomes singular. These points define the world lines of monopoles in four dimensional space.

We introduce now a gauge, that interpolates between MAG and LAG, by defining a functional

$`Z_A(U^V)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}𝑑f(𝒓)f^2(𝒓)\mathrm{exp}\left\{\stackrel{~}{F}_{\mathrm{gf}}(U^V,f)\right\},`$ (4)
where $`\stackrel{~}{F}_{\mathrm{gf}}(U^V,f)`$ is defined as

$$\stackrel{~}{F}_{\mathrm{gf}}(U^V,f)=\underset{𝒓}{}\underset{\mu }{}\left\{f(𝒓)f(𝒓+𝝁)\frac{1}{2}\mathrm{Tr}\left[\sigma _3U_\mu ^V(𝒓)\sigma _3U_\mu ^V(𝒓)\right]f^2(𝒓)\right\}+\frac{\lambda }{\eta ^2}\underset{𝒓}{}\left[f^2(𝒓)\eta ^2\right]^2.$$
(5)
The prescription for finding the appropriate abelian gauge is to find the maximum of $`Z_A(U^V)`$ in terms of the $`V(𝒓)`$ (or rather $`g_{a3}^A`$) and $`f(𝒓)`$. The gauge defined by (4) is similar to the LAG gauge-fixing functional, in the sense that the optimal local vector $`f(𝒓)`$ vanishes at certain points (or rather, owing to three constraints, on curves in spacetime). At these points the gauge transformation becomes singular. Just like in LAG, these points correspond to monopoles.

Notice now that the $`\lambda \mathrm{}`$ limit is equivalent to MAG. Indeed, this limit freezes out the integration over $`f(𝒓)`$ and leaves one with the exponential of (2), multiplied by $`\eta ^2`$. Furthermore, the $`\lambda 0`$ limit is equivalent to LAG. The introduction of the field $`f(𝒓)`$ is equivalent to relaxing the constraint on the magnitude of the gauge group element in adjoint representation, due to the replacement $`g_{3a}^A(𝒓)f(𝒓)g_{3a}^A(𝒓)`$. Only the sum, $`_{𝒓a}f^2(𝒓)(g_{0a}^A)^2(𝒓)=_𝒓f^2(𝒓)`$ is restricted, to make the integral finite. The general case of finite $`\lambda `$ that interpolates between these two gauge fixing methods is also an appropriate Abelian gauge and will henceforth be called the general Abelian gauge (GAG).

Next we propose an alternative to GAG, a conventional gauge fixing term, $`S_{\mathrm{gf}}`$, inspired by the functional $`Z_A(U^V)`$ and defined as

$$e^{S_{\mathrm{gf}}(U^V)}=Z_A(U^V)e^{\stackrel{~}{S}},$$
(6)
where

$$e^{\stackrel{~}{S}}=𝑑V^{}(𝒓)Z_A(U^V^{}).$$
(7)
The integral of (6) over all gauge transformations is 1 and the modified action, $`S=S_G+S_{\mathrm{gf}}`$, leads to the same expectation value of all gauge invariant quantities as the original plaquette action, $`S_G`$. We call this a smeared abelian gauge (SAG). Unlike in GAG, gauge configurations are understood to be generated in the presence of the gauge fixing term.

It is comparatively easy to establish the relationship between SAG and GAG. What is more important this relationship extends to non-gauge invariant quantities, such as gauge field configurations. Consider the Boltzmann factor in SAG,

$$e^{S_G(U^V)S_{\mathrm{gf}}(U^V)}.$$

$`S_{\mathrm{gf}}(U^V)`$ has a sharp minimum at $`g_{a3}^A(𝒓)=\delta _{a3}`$. The integral in the denominator of (6) can be evaluated by the saddle point method in the nonabelian parameters of the gauge group. Up to a determinant, it gives $`Z_A(U)`$, the maximum of the numerator of (6), over $`V(𝒓)`$. Though the integral of the Boltzmann factor over $`U_\mu (𝒓)`$ is gauge invariant, $`Z_A(U^V)`$ will give an appreciable contribution to the integral over $`U_\mu (𝒓)`$ only if $`g_{a3}^A(𝒓)\delta _{a3}`$. Thus, when one generates an ensemble of gauge field configurations using the Boltzmann factor $`\mathrm{exp}\{S_GS_{\mathrm{gf}}\}`$ one obtains only configurations that are very close to GAG configurations. The equivalence between GAG and SAG would become exact at $`\eta \mathrm{}`$.111One can rescale function $`f`$ by $`\eta `$ and then it is easy to see that all values of $`f`$ and $`g^A`$ but those maximizing the gauge fixing term are frozen out. In other words, at large $`\eta `$ SAG is only slightly smeared around GAG.

We now consider SAG from a different angle. Aside from $`\stackrel{~}{S}`$, the action is just that of a traceless adjoint Higgs theory. To prove this statement, we use the gauge invariance of the measure, $`dU_\mu (𝒓)`$, of $`S_G`$, and of $`\stackrel{~}{S}`$ to introduce an extra integration over gauge transformation $`V`$, and write the partition function as

$$Z=𝑑V(𝒓)𝑑U_\mu (𝒓)e^{S_G(U)S_{\mathrm{gf}}(U^V)}=𝑑U_\mu (𝒓)𝑑\mathrm{\Phi }(𝒓)e^{S_G(U)S_H(U,\mathrm{\Phi })\stackrel{~}{S}(U)},$$

where we defined the Higgs action as

$$S_H(U,\mathrm{\Phi })=\frac{1}{2}\underset{𝒓\mu }{}\mathrm{Tr}\left[\mathrm{\Phi }(𝒓)U_\mu (𝒓)\mathrm{\Phi }(𝒓+𝝁)U_\mu ^{}(𝒓)\mathrm{\Phi }^2(𝒓)\right]+\frac{\lambda }{\eta ^2}\underset{𝒓}{}\left\{\frac{1}{2}\mathrm{Tr}[\mathrm{\Phi }^2(𝒓)]\eta ^2\right\}^2.$$
(8)
Here we defined $`\mathrm{\Phi }(𝒓)=_af(𝒓)g_{a3}^A(𝒓)\sigma _a`$ and used the relationship

$$𝑑ff^2(𝒓)𝑑V(𝒓)𝑑\mathrm{\Phi }(𝒓)𝑑\psi (𝒓),$$
(9)
where $`\psi `$ is the Euler angle parameterizing the little group of $`\mathrm{\Phi }(𝒓)`$. Consequently, the integrand is independent of $`\psi `$. To prove (9) we use the representation of the Haar measure by Euler angles and (3).

Clearly, (8) is just a standard lattice Higgs action. To show that SAG is indeed equivalent to a Higgs theory we only need to investigate whether the presence of the term $`\stackrel{~}{S}`$, defined in (7), changes the nature of the action.

First of all, $`\stackrel{~}{S}(U)`$ is a gauge invariant function of the gauge field, $`U_\mu (𝒓)`$. Though it is nonlocal, as we will show later, it becomes local in the relevant, $`\eta \mathrm{}`$, $`g0`$ limit. It is also the generating functional of connected diagrams of an adjoint scalar field theory in a background gauge field, $`U_\mu (𝒓)`$. It must be composed of gauge invariant combinations of the gauge field, Wilson loops. $`\stackrel{~}{S}(U)`$ also has the property of local $`Z_2`$ invariance, i.e. unlike the standard plaquette action it is invariant under the local replacement $`U_\mu (𝒓)U_\mu (𝒓)`$. This is the property of a gauge action with an adjoint representation trace of the plaquette product of gauge fields.

To investigate the behavior of $`\stackrel{~}{S}`$ we first use a hopping parameter expansion . That implies expanding the integrand in the hopping term, integrating over the gauge group, and re-exponentiating $`𝑑g(𝒓)Z_H^g(\lambda )`$. The hopping parameter expansion assumes small $`\eta `$.

The leading nontrivial term of this expansion is

$$\stackrel{~}{S}(U)\frac{8(f_{\mathrm{av}}^2)^4}{81}\underset{𝒓\mu \nu }{}|U_{\mu \nu }^P(𝒓)|^2=\frac{8(f_{\mathrm{av}}^2)^4}{81}\underset{𝒓\mu \nu }{}\left[U_{\mu \nu }^{P,A}(𝒓)1\right]$$

where $`U_{\mu \nu }^{P,A}(𝒓)`$ is the plaquette in adjoint representation and

$$f_{\mathrm{av}}^2=\frac{_0^{\mathrm{}}𝑑ff^2e^{(\lambda /\eta ^2)(f^2\eta ^2)^2}}{_0^{\mathrm{}}𝑑fe^{(\lambda /\eta ^2)(f^2\eta ^2)^2}}\eta ^2.$$

This adds a negative adjoint representation term to the usual plaquette term of the gauge action,

$$\frac{1}{g^2(a)}U_{\mu \nu }^P(𝒓)\frac{1}{g^2(a)}U_{\mu \nu }^P(𝒓)\mathrm{\Delta }\beta _AU_{\mu \nu }^{P,A}(𝒓),$$

where in the leading order approximation $`\mathrm{\Delta }\beta _a=8(f_{\mathrm{av}}^2)^4/81.`$

Provided the hopping parameter expansion is convergent, higher order terms can only have one of three effects:

1. They can renormalize $`\mathrm{\Delta }\beta _A`$,
2. Loop terms larger than the plaquette (nonlocal terms) may appear in the adjoint representation,
3. Characters in higher order integer isospin representations may also appear.

Lattice gauge theories with actions mixing the characters of plaquettes in fundamental and adjoint representations have been studied thoroughly . The conclusion of those investigations is that though such a modification of the action leads to new phase transition lines in the space of couplings. None of those transitions is deconfining. In the continuum limit, when all couplings tend to zero, all of these theories are equivalent to a nonabelian gauge theory with a modified coupling.

In fact, there is no real need to evaluate all these terms, as we can predict the resulting effective action using simple considerations. The inclusion of the gauge fixing term does not change the expectation value of any gauge invariant quantity of the pure gauge theory. These quantities are dependent only on the gauge coupling, or alternatively, on the scale of the running coupling, $`g`$. The lattice coupling runs as

$$\frac{1}{g^2(a)}\frac{1}{g_0^2}2\beta _G\mathrm{log}[\mathrm{\Lambda }(a)/\mathrm{\Lambda }_0],$$

where $`\mathrm{\Lambda }(a)1/a`$ is the lattice scale and $`\beta _G=11/24\pi ^2`$; $`\mathrm{\Lambda }_0`$ and $`g_0`$ refer to another, fixed scale.

The action is also equivalent to an adjoint Higgs action with an effective gauge coupling (the modification coming from the contribution of $`\stackrel{~}{S}`$), $`g_{\mathrm{eff}}`$. The effective gauge coupling is such that the inclusion of the Higgs boson changes its running into that of the non-gauge fixed coupling, namely

$$\frac{1}{g_{\mathrm{eff}}^2}\frac{1}{g_0^2}2(\beta _G+\beta _H)\mathrm{log}(\mathrm{\Lambda }(a)/\mathrm{\Lambda }_0),$$

where $`\beta _H=1/24\pi ^2`$. Physical quantities, depending on the gauge field only, must be the same when the Higgs theory action $`S=S_G+\stackrel{~}{S}+S_H`$ is used as when the $`S_G`$ action is used alone. As, at small $`g(a)`$ these quantities are dependent only on the scale, $`\mathrm{\Lambda }(a)`$, the scales in the two calculation must coincide. This leads to a relationship between the effective coupling and the original gauge coupling

$$g_{\mathrm{eff}}^2(a)g^2(a)\frac{\beta _G}{\beta _G+\beta _H}=\frac{11}{12}g^2(a),$$
(10)
valid at sufficiently small $`g(a)`$.

Let us return now to the discussion of the non-locality of $`\stackrel{~}{S}`$.222The author is indebted to P. de Forcrand for calling his attention to this subject. We can show that the nonlocality goes away in the limit of $`\eta \mathrm{}`$ and $`g0`$ (continuum limit), not only in the limit of $`\eta 0`$. This is the very limit in which SAG goes over into GAG and the lattice gauge theory goes over into the continuum theory. First of all, after rescaling $`\mathrm{\Phi }\mathrm{\Phi }\eta `$, we have $`S_H\eta ^2\stackrel{~}{S}_H`$, where $`\stackrel{~}{S}_H`$ is independent of $`\eta `$. As

$$\stackrel{~}{S}=\mathrm{log}𝑑\mathrm{\Phi }e^{S_H},$$

the large $`\eta `$ limit is equivalent to the semiclassical limit as $`\eta ^2`$ plays the role of $`1/\mathrm{}`$. In this limit only single loop terms contribute to $`\stackrel{~}{S}`$. In particular, all reference to the self-coupling of the Higgs field, $`\lambda `$, is of $`O(\eta ^2)`$, as it appears in multi-loop terms only. As the external fields in $`\stackrel{~}{S}`$ are gauge fields and the internal lines are Higgs propagators (there are no internal gauge lines!) all the couplings of the one loop term are gauge couplings. Then the calculation of $`\stackrel{~}{S}`$ is very similar to the calculation of the running gauge coupling using the background gauge field method . Thus, among the one loop contributions to the logarithm of the functional determinant only diagrams with at most 4 background gauge fields are divergent. As it is shown on Ref.  the divergent contributions can be combined into the local term

$$\mathrm{log}det\mathrm{\Delta }=d^4x\left(\frac{1}{4}C(F_{\mu \nu }^a)^2+\text{ finite terms}\right),$$

where $`C`$ is a logarithmically divergent constant, depending on the transformation properties of the field in the loop. When the running coupling is calculated  one must consider the contribution of the dynamical gauge field, the fermion fields, and the ghost field, as well. We only need to calculate a scalar field loop, so the contribution is slightly different. We obtain

$$\stackrel{~}{S}=\frac{1}{2}\mathrm{log}(a\mathrm{\Lambda }_{\mathrm{lattice}})\frac{1}{24\pi ^2}d^4x(F_{\mu \nu }^a)^2+\text{ terms finite at }a\text{=0},$$

where the scale was chosen to be $`\mathrm{\Lambda }_{\mathrm{lattice}}`$. As in the scaling limit, the gauge coupling in the pure gauge term scales with the lattice parameter as

$$\frac{1}{g^2}=\beta _G\mathrm{log}\left(\frac{1}{a\mathrm{\Lambda }_{\mathrm{lattice}}}\right).$$

$`\stackrel{~}{S}`$ combines with $`S_G`$ to result in the replacement of $`g^2`$ by $`g_{\mathrm{eff}}^2`$, as indicated in (10). At large $`\eta `$ and in the continuum limit the role of $`\stackrel{~}{S}`$ is to renormalize the gauge coupling. Then, in the continuum limit, at large $`\eta `$, the gauge fixed pure gauge theory is indeed equivalent to a Higgs theory. It is clear that at smaller values of $`\eta `$ or/and at larger gauge coupling nonlocal effects are more pronounced and $`\stackrel{~}{S}`$ is dependent on the parameters of the Higgs potential. Nonlocal effects are of $`O(1/\mathrm{log}(a\mathrm{\Lambda })`$ and of $`O(\eta ^2)`$. We have also seen, however, that in the continuum limit at small $`\eta `$ these theories are also equivalent.

The equivalence of these lattice theories does not alone provide an obvious relationship between monopole configurations found in them, as the monopoles themselves are not gauge invariant objects. To prove such an equivalence we need to examine not only partition functions but also Boltzmann factors.

Assuming that the gauge coupling is small enough to apply continuum scaling, gauge configurations of the continuum theory and the lattice theory are essentially equivalent. To find classical solutions of the continuum theory we need to minimize the action in terms of the fields. The maximum of the Boltzmann factor should be calculated by minimizing $`S_G+\stackrel{~}{S}+S_H`$ in terms of the gauge fields and gauge transformations, as well. This minimization is equivalent to the minimization in terms of the Higgs field and the gauge field. In the continuum limit this minimization should lead to a combination of configurations containing a variable number of finite energy classical solutions, monopoles  . Monopoles are identified by the radial dependence of the magnitude of the Higgs field. The magnitude of the Higgs field itself is gauge invariant. The final step in identifying the monopoles with those of GAG is the diagonalization of the Higgs field with a gauge transformation.

The conclusion of the above considerations is that

1. The appearance of monopoles in gauge fixed pure gauge theories follows from their existence in Higgs theories.
2. The analytic form, interaction, condensation, etc. of monopoles should be very similar to those in Higgs theory.
3. The mechanism of confinement, provided monopoles are responsible for it, is the same in pure gauge and Higgs theories.

Before turning to vortices let us apply our method to calculate the distribution of a single monopole in Abelian projected GAG. To do that we take the the ’t Hooft-Polyakov monopole and diagonalize the Higgs boson by a gauge transformation. The gauge and Higgs fields are given by

$`𝑨(𝒓)`$ $`=`$ $`𝒓\times 𝝈F(r),`$
$`\mathrm{\Phi }(𝒓)`$ $`=`$ $`𝒓𝝈W(r),`$ (11)
where $`F(r)`$ and $`W(r)`$ are the radial wave functions.

Then the abelian gauge transformation is a rotation around the $`z`$ axis by an angle $`\varphi `$, followed by a rotation around the $`y`$ axis by an angle $`\theta `$, where $`\varphi `$ and $`\theta `$ are polar coordinates. Applying the same gauge transformation to the gauge field and projecting to the abelian component we obtain (we set the gauge charge equal to 1)

$$\frac{1}{2}\mathrm{Tr}[\sigma _z(UA_\mu U^{}+iU_\mu U^{})]=\frac{z}{r}_\mu \varphi .$$

The magnetic field obtained from this vector potential is that of a single magnetic charge at the origin

$$𝑩=\frac{1}{2}\frac{\widehat{𝒓}}{r^2}+\text{ Dirac string}.$$

It is worth noting that the abelian projection of gauge field (11) without abelian gauge fixing leads to

$$A_\mu (𝒓)=_\mu \varphi r^2F(r),$$

that is a more like a vector potential for a magnetic flux than for a monopole.

## 3  Central gauges and Higgs theories

Now we turn to center gauges and vortices. Lattice studies have shown the relevance of $`Z_2`$ vortices in pure lattice $`SU(2)`$ gauge theory. These studies use two alternative gauges, the Maximal Center gauge (MCG)  and the Laplacian Center gauge (LCG)  . The gauge transformations are followed by center projection.

MCG is defined through the maximization of the following functional (we use notations, identical to the ones used for the discussion of Abelian gauges) over gauge transformation $`V(x)`$:

$$F^C(U)=\frac{1}{4}\underset{x\mu }{}|\mathrm{Tr}U_\mu ^V(𝒓)|^2=\frac{1}{6}\underset{𝒓\mu }{}\mathrm{Tr}\left[\underset{i}{\overset{3}{}}\sigma _iU_\mu ^V(𝒓)\sigma _iU_\mu ^V(𝒓)\right].$$
(12)
$`F^C(U)`$ can be rewritten in terms of adjoint representation quantities as

$$F^C(U)=\frac{1}{3}\underset{𝒓\mu }{}\underset{iab}{\overset{3}{}}g_{ai}^A(𝒓)U_\mu ^{A,ab}(𝒓)g_{bi}^A(𝒓+𝝁).$$
(13)
If one writes the gauge field in the canonical form

$$U_\mu ^g(𝒓)=a_0+i𝒂\mathbf{}𝝈,$$

where $`a_0^2+𝒂^2=1`$, then $`F^C(U)=a_0^2`$. The maximization of $`a_0^2`$ results in a gauge field as close to the center of the gauge group ($`a_0=\pm 1`$) as possible.

In fact, it is sufficient to keep only two terms of the sum over $`i`$ in (12) and (13), as the adjoint representation gauge transformation is defined uniquely by two of the orthogonal vectors $`g_{a1}^A`$ and $`g_{a2}^A`$. This is used in LCG, which relaxes the normalization and orthogonality conditions on the vectors $`g_{ai}^A`$. Vortices are found at points where the two vectors are parallel. These points form lines in three dimensional space. At these points the gauge transformation is singular. In what follows we will restrict ourselves to such two term gauge fixing functionals, though everything what we do can easily be generalized to the three term case, such as the original MCG.

MCG and LCG can again be extended into a Generalized Center Gauge (GCG), which also fixes the gauge, up to the center of the $`SU(2)`$ gauge group, $`Z_2`$. A convenient form for the gauge fixing functional is

$$Z_C(U^V)=𝑑\mathrm{\Phi }^1(𝒓)_0^{\mathrm{}}𝑑f(𝒓)\rho ^2(x)\mathrm{exp}\left\{\stackrel{~}{F}_{\mathrm{gf}}(U^g,f)S_H(U,\mathrm{\Phi }^1)S_{\mathrm{mix}}(V,f,\mathrm{\Phi }^1)\right\}$$
(14)
where $`\stackrel{~}{F}_{\mathrm{gf}}(U^V,f)`$ and $`S_H(U,\mathrm{\Phi }^1)`$ were defined in (5) and (8), respectively, while the ‘mixing term’ is given by

$$S_{\mathrm{mix}}(V,f,\mathrm{\Phi }^1)=\frac{\stackrel{~}{\lambda }}{\eta ^2}\underset{𝒓}{}[f(𝒓)g_{a3}^A(𝒓)\mathrm{\Phi }_a^1(𝒓)c\eta ^2]^2,$$

where $`c`$ is the cosine of the asymptotic angle between $`g_{a3}^A`$ and $`\mathrm{\Phi }_a`$. When we maximize (14) in terms of $`V(𝒓)`$ and $`f(𝒓)`$ we obtain an abelian gauge. The exponent will be dominated by the configuration of the Higgs field, gauge transformation, and $`f(𝒓)`$ that maximizes the exponent. In contrast to the Abelian gauge fixing of the previous section this limit fixes the abelian gauge group, as well. The Abelian subgroup is fixed by a rotation around the 3rd axis (the little group of $`g_{a3}^A(𝒓)`$), such that $`\mathrm{\Phi }_2^1(𝒓)=0`$. Vortices appear on lines in the 3D space where $`\mathrm{\Phi }_a^1(𝒓)`$ is parallel with $`g_{a3}^A(𝒓)`$ and where such a gauge fixing is impossible. Alternatively, in a slightly more symmetric manner, we could also define the gauge transformation by rotating the sum of the normalized Higgses parallel to the $`z`$ axis and rotating the difference into the $`(xz)`$-plane.

When $`c=0`$ and $`\lambda ,\stackrel{~}{\lambda }\mathrm{}`$ GCG reduces to MCG (or rather an analogous 2 term gauge), as $`\mathrm{\Phi }_a^1(𝒓)`$ is kept normalized to $`\eta `$ and orthogonal to $`g_{a3}^A(𝒓)`$. When $`\lambda _1=\lambda _2=0`$ GCG reduces to LCG.

Next we introduce the smeared center gauge (SCG) such that the gauge fixing term in the Lagrangian is defined as

$$e^{S_{\mathrm{gf}}}=Z_C(U^g)e^{\overline{S}},$$

where

$$e^{\overline{S}}=𝑑V^{}(𝒓)Z_C(U^V^{})$$

In a manner similar to our discussion of Abelian gauges we can show that on large lattices SCG is very close to GCG and we can also show that SCG is completely equivalent to a lattice Higgs theory with two adjoint Higgses. This follows because after a gauge transformation one can introduce the second Higgs field as

$$\mathrm{\Phi }_a^2(𝒓)=V(𝒓)g_{a3}^A(𝒓)$$

and prove that the action can be written as $`S=S_G+S_{\mathrm{gf}}=S_G+\overline{S}+S_H^2`$, where $`S_H^2`$ is a Lagrangian for two adjoint Higgses. $`S_G+\overline{S}`$ is again a modified pure gauge action. The modified gauge coupling should be equal to the standard pure gauge coupling minus the contribution of two adjoint Higgs bosons, as $`S`$ is equivalent to a gauge fixed version of the pure gauge theory. Thus the relationship between the gauge couplings is $`g_{\mathrm{eff}}^2(a)(11/13)g^2(a).`$

The Boltzmann factor of the GCG must be dominated by vortex configurations, if the Boltzmann factor of the Higgs theory is. To make the connection between GCG and SCG distributions we need to fix the vortex gauge such that one of the Higgs bosons is rotated parallel to the 3rd axis in isospace while the second Higgs boson is rotated into the $`13`$ plane. Thus, the ‘thick’ vortices appearing in GCG have the same location and radial dependence as the vortices of the Higgs theory.

To illustrate the equivalence of GCG and SCG we calculate the shape of a single vortex in center gauge fixed pure gauge theory, by making use of a classical solution of the Higgs theory with two adjoint Higgs bosons . The solution is of the form

$`\mathrm{\Phi }^{1,2}(𝒓)`$ $`=`$ $`U(\varphi )[\chi _d(\rho )\sigma _3\pm \chi _{}(\rho )\sigma _1]U^{}(\varphi ),`$
$`A_\mu (𝒓)`$ $`=`$ $`a(\rho )_\mu \varphi \sigma _3,`$ (15)
where the boundary conditions require that $`\chi _{}(0)=a(0)=0`$ and $`a(\mathrm{})=1`$, the topological charge. $`\rho `$ and $`\varphi `$ are cylindrical coordinates, along with $`z`$ and

$$U(\varphi )=e^{i\varphi \sigma _3/2}.$$

The normalized components of the Higgs field also satisfy boundary conditions at infinity: $`\chi _{}(\mathrm{})=\sqrt{(1c)/2}`$, $`\chi _d(\mathrm{})=\sqrt{(1+c)/2}`$. $`c`$ is the angle between the two Higgs bosons at infinity. (15) satisfies the constraint that the two Higgs fields become parallel at the location of the vortex, along the $`z`$ axis. The $`z`$ axis is also the locus of the singularity of $`U(\varphi )`$.

To predict the form of a thick vortex on the lattice, in pure gauge theory, we need to bring the Higgs bosons to a standard form. This is the form in which the fields have only components in the $`(xz)`$ plane and the sum of the Higgses is diagonal. To get to that form we need to use the singular gauge transformation $`U^{}(\varphi )`$. As $`U`$ commutes with the gauge field of (15) the only change of the gauge field will be the addition of the inhomogeneous term, $`iU^{}_\mu U`$. This will result in the following change in the gauge field (the gauge charge has been set equal to unity):

$$A_\mu (𝒓)A_\mu (𝒓)+iU^{}_\mu U=[a(\rho )1]_\mu \varphi \sigma _3.$$

The gauge field remains diagonal after the gauge transformation. It is singular on the $`z`$-axis and it goes to zero exponentially at $`\rho =\sqrt{x^2+y^2}\mathrm{}`$.

To understand what happens on the lattice and after center projection we rewrite the Higgs fields, the vector potential, and the gauge transformation into a form appropriate for the lattice. The Higgs field (15) is already in such a form. The lattice gauge field is a unitary matrix. The gauge field in (15) can be conveniently written as

$$U_\mu (𝒓)=\mathrm{exp}\left\{\frac{i}{2}a(\rho (𝒓))[\varphi (𝒓+𝝁)\varphi (𝒓)]\sigma _3\right\},$$
(16)
where $`\rho (𝒓)`$ and $`\varphi (𝒓)`$ are the cylindrical coordinates at lattice points. The singular gauge transformation that brings the Higgs field into the $`(xz)`$ plane is

$$V(𝒓)=\mathrm{exp}\left\{\frac{i}{2}\varphi (𝒓)\right\}.$$

Thus, the gauge transformed gauge field is

$$U_\mu ^V(𝒓)=V(𝒓)U_\mu (𝒓)V^{}(𝒓+𝝁)=\mathrm{exp}\left\{\frac{i}{2}[a(\rho (𝒓))1][\varphi (𝒓+𝝁)\varphi (𝒓)]\sigma _3\right\}.$$

Let us examine the gauge field near the $`z`$ axis. At $`\rho 0`$ the radial component of the original gauge field is negligible, $`a(\rho )\rho ^2`$, and the new gauge field has the form of a pure gauge transformation. When one traverses around the $`z`$ axis one needs to define a surface in the three dimensional space, with a boundary at the $`z`$ axis, such that on this surface the value of $`\varphi (x)`$ jumps by $`2\pi `$. This surface can be chosen arbitrarily, so we choose the half-$`(xz)`$ plane containing the negative half of the $`x`$ axis. If the gauge coupling small, the argument of the gauge field $`U_\mu (x)`$ is small and at central projection projects to the identity matrix. When we intersect the half plane, however,

$$U_y(𝒓)e^{\frac{i}{2}[\pi (\pi )]}=1,$$

and the gauge field projects to the non-trivial element of the center, -1. This sheet of negative elements of the center ends in the $`z`$-axis, where a $`Z_2`$ vortex appears. Thus, we showed that in the center gauge, after center projection a $`Z_n`$ vortex appears wherever a vortex appears in the Higgs theory. At the same time we have given a representation for the thick vortex, obtained after transformation to center gauge, but before center projection.

Strictly speaking, a secondary vortex also appears at some distance from the $`z`$ axis. The location of that vortex is determined by the equation

$$a(\rho )1\frac{1}{2},$$

which ensures that the jump of the argument of $`U_\mu `$ is smaller than $`\pi /2`$. At large $`\rho `$ the radial component of the gauge field $`a(\rho )e^{m_A\rho }`$, where the auxiliary gauge mass, $`m_A=g(a)\eta `$. Then the secondary vortex will appear at $`\rho 1/m_A`$. If, in the continuum limit, $`g(a)`$ is sufficiently small then this vortex will be far away from the $`z`$ axis, not affecting the confinement picture. The location of the primary vortex is gauge invariant and of the secondary vortex depends on the gauge. On the contrary, if we perform the center projection without transforming (15) to center gauge, then the $`Z_2`$ vortex does not appear at the $`z`$ axis, we would only have the above described secondary vortex and the simple vortex condensation picture would not be recognizable.

## 4 Conclusion

We have shown that pure gauge theories in Abelian or Center gauges are equivalent to Higgs theories with one or more adjoint representation Higgs bosons. The equivalence extends to gauge configurations. Monopoles and vortices and their interactions, condensation, etc. in gauge fixed pure gauge theories can be studied analytically by investigating similar objects in Higgs theories, where they appear as classical solutions.

We encounter a problem when we apply our method of connecting pure gauge theories with Higgs theories in the limiting cases of Maximal Abelian and Maximal Center gauges. Then our construction leads to gauged nonlinear sigma models. As far as we know, the relevant gauged nonlinear sigma models do not have classical solutions in 3+1 dimensions. The reason is simple: Higgs fields in nonlinear sigma models cannot vanish. In gauged linear sigma models the location of the zeros marks the center of monopoles or vortices. This is, of course, in agreement with the fact that monopoles and vortices do not appear as singularities of the gauge transformation in MAG and MCG either. It would be interesting to investigate, using nonlinear sigma models, why the projected theories still show many of the characteristics found in General Abelian and General Center gauges.

The link between gauge fixed pure gauge theories and Higgs theories can be employed to investigate important questions concerning confinement. One example that comes to mind is the recently uncovered intriguing relationship between monopole and vortex condensations  . Simulations show that monopoles do not form a Coulomb gas as was originally assumed, but rather they line up to form vortices, thereby making vortices more fundamental for the mechanism of confinement.

Another important question is why multiple non-abelian vortex configurations appear at all. ‘t Hooft’s argued  that in a $`SU(2)/Z_2`$ theory the homotopy group has only one non-trivial class, so vortices should coalesce into a single vortex or no vortex at all. The situation is somewhat different in $`U(1)`$ gauge theory where the homotopy group is $`Z`$, and infinitely many kinds of vortices exist. Lattice evidence supports, however, the existence of multi-vortex configurations. It is conceivable that studying interactions of vortices in Higgs theories one will be able to answer this important question, as well.

The author thanks the U.S. Department of Energy for partial support through grant #DE FG02-84ER-40153.


format: text

text:
# Superfluid vortices in neutron stars

## 1 Introduction

A natural probe of the internal stucture and dynamics of neutron star crusts is represented by pulsar glitches, which are sudden accelerations of the star’s rotational frequency. The observation of more than seventy events in about thirty pulsars provides a general view of the phenomenon. At present the standard model to explain pulsar glitches focuses on the interaction of superfluid vortices with the crust of the neutron star (Pines & Alpar (1985)). Some variants include crust breaking, resulting in a starquake, or stellar plate tectonics due to the action of magnetic threading of proton vortex lines.

The standard model is built in analogy to type II superconductors, where the attractive interaction of vortex cores with the lattice defects of the metal prevents the motion of vortex lines through a metal sample. An important but unknown parameter is the pinning energy of a vortex to a nucleus, the difference in energy between the configurations in which a vortex and a nucleus are superimposed and the one in which they are far apart and non-interacting. Usually this quantity is calculated simply as the loss of condensation energy between the two configurations. Since the core of a vortex is composed of normal matter and the neutron superfluid which occupies the nuclear volume is, for sufficiently high densities, a weaker superfluid (i.e. has a smaller local pairing gap) than the superfluid outside the nuclei, it is energetically favourable for a vortex to remain attached to a nucleus. However, this analysis neglects important effects like the change in kinetic energy (Epstein & Baym 1988; Pizzochero, Viverit & Broglia 1997) and the presence of bound states (De Blasio & Elgarøy (1999)).

Furthermore, the vortex is depicted as a cylinder composed of uniform normal matter with a sharp radius equal to the BCS coherence length of the superfluid (the BCS coherence length is the mean square radius of Cooper pairs, see De Blasio et al. (1997) for the case of hadronic matter). This is a very rough model since the gap changes smoothly inside the vortex and not abruptly, as assumed. The vortex-nucleus interaction results from effects requiring a full quantum mechanical solution of the BCS equations in the cell. Another approach which has been attempted is the Ginzburg-Landau model (Epstein & Baym (1988)). This model provides a good treatment of the kinetic energy density, which is an important ingredient in the energy balance within the interacting volume. Unfortunately, the pairing is not well treated in the Ginzburg-Landau equations since the two main conditions for their applicability, that the temperature should be close to the critical one and that there should be only slow and smooth spatial variations in the order parameter, are not fulfilled in a neutron star crust.

The relevant quantity, the vortex-nucleus pinning force, can in principle be calculated from the variation of the pinning energy with respect to the vortex-nucleus distance. An approximate expression can be found by dividing the pinning energy by the minimal distance where the two objects can be considered as separate. The pinning force enters in two important ways in the standard model for pulsar glitches. Vortices are released from the crust when the frequency lag between the superfluid and the crust reaches a critical value. They unpin in clusters, decreasing suddenly the angular momentum of the superfluid. The glitch results from the conservation of angular momentum, causing the crust to increase its spinning rate. Therefore, the critical frequency is proportional to the value of the local pinning force.

The second issue where pinning force is relevant is postglitch relaxation. This is particularly interesting, because time analysis following a glitch is now available for some glitches and especially for the Vela pulsar. In the vortex creep model, vortices that have abandoned their pinning centers repin with other nuclei within characteristic times which also depend on the value of the pinning force. The pinning force is of the order $`E_p/\xi `$, where $`\xi `$ is an appropriate length scale of the vortex–nucleus interaction potential, roughly the size of the vortex core, and $`E_p`$ is the pinning energy. This length is usually taken to be the BCS coherence length, but recent calculations (De Blasio & Elgarøy (1999)) show that the vortex core can be significantly smaller than the BCS coherence length. This indicates that if the pinning energies were comparable, the pinning force should be larger. It is clear that a better understanding of the pinning model for pulsar glitches will be possible when more reliable calculations of the structure of a superfluid neutron vortex are available.

Again, the study of superfluid vortices in neutron stars is not exhausted by considering solely neutron vortices in the crust. Below the inner crust, neutrons couple in the $`{}_{}{}^{3}P_{2}^{}`$ channel to generate a triplet superfluid, while protons, having a much smaller density, couple in the singlet channel. The flux of protons around the vortex core generates a magnetic field parallel to the rotational axis. In addition, a large number ($`10^{19}`$ cm<sup>-3</sup>) of proton flux lines is generated by the very high magnetic field ($`10^{12}`$G) inside the core. The interaction of these flux lines with neutron vortices drifting towards the exterior of the star might be an efficient mechanism for magnetic field decay (see for example Ruderman (1997)).

The study of these physical situations may require a more accurate model for vortices in hadronic matter. A fully quantum mechanical model of vortices in type II superconductors has been extensively investigated in some publications (Gygi & Schlüter (1991); Hayashi et al. (1998)). It is possible to study the properties of neutron vortices parallelling the model for superconductors. Experiments with scanning-tunnelling microscopy and refined numerical calculations (Gygi & Schlüter (1991)) have confirmed the theoretical prediction Caroli, de Gennes & Matricon (1964); Bardeen et al. (1969) that bound states are formed in the center of the core of a vortex in type-II superconductors. The pairing gap decreases to zero at the center of the core and grows to an asymptotic value within distances of the order of or larger than the coherence length.

In the present paper we extend our previous work (De Blasio & Elgarøy (1999)), hereafter Paper I, on the microscopic structure of neutron vortex lines. In Paper I we made use of the Bogoliubov-de Gennes equations (de Gennes (1989); Ketterson & Song (1999)), that have been successfully developed in studies of vortices in type II superconductors (Gygi & Schlüter (1991); Hayashi et al. (1998)) and more in general non-homogeneous superconductivity.

In Paper I we were mostly interested in the vortex core size, since this quantity influences the pinning energy and force. We will here continue our study of this quantity, providing more complete results and extending our calculations towards a fully microscopic calculation of pinning energies. Furthermore, we will show results for quantities such as density distributions, pairing potentials and vortex tension. The calculations presented here are also more refined than those in Paper I, as they are fully three-dimensional and include effects like the Hartree mean field and a density-dependent pairing force. Furthermore, since our earlier paper was brief on the technicalities, we will give a more detailed account of the formalism and the numerical solution of the relevant equations.

## 2 The model

The $`{}_{}{}^{1}S_{0}^{}`$ neutron superfluid in the inner crust is spatially nonhomogeneous due to the presence of a nuclear lattice and the vortices induced in the superfluid by the rotation of the star. To describe this system one needs to go beyond the standard BCS formalism for fermion pairing. A formalism for this has existed in solid state physics for several years in the form of the Bogoliubov-de Gennes (BdG) equations (de Gennes (1989)). For readers not familiar with these, we will give a sketch of their derivation.

The theory is most conveniently formulated in terms of the neutron field operators $`\psi _\sigma (𝐫)`$ and $`\psi _\sigma ^{}(𝐫)`$ which respectively destroy and create a neutron with spin projection $`\sigma `$ in position $`𝐫`$. As a starting point, we choose a pure pairing Hamiltonian

$$H=H_0+H_{\mathrm{int}}=H_0|g|d^3r\psi _{}^{}(𝐫)\psi _{}^{}(𝐫)\psi _{}(𝐫)\psi _{}(𝐫),$$
(1)
where $`|g|`$ is the pairing strength (we use the notation $`|g|`$ to avoid confusion about the sign of the pairing strength, which in some works is defined to be negative, in others to be positive). This corresponds to a zero-range pairing interaction of the form

$$v(𝐫,𝐫^{})=|g|\delta (𝐫𝐫^{}).$$
(2)
and necessitates the use of a momentum cutoff confining the interaction to a narrow set of states near the Fermi surface. Note that in general the pairing strength is space dependent due to density variations induced by the mean field or by the vortex. We will return to this point later, for the time being we assume that it is a constant.

The one-body hamiltonian has the form

$$H_0=d^3r[\psi _{}^{}(𝐫)h_0\psi _{}(𝐫)+\psi _{}^{}(𝐫)h_0\psi _{}(𝐫)]$$
(3)
where

$$h_0=\widehat{t}+W(𝐫)+S(𝐫)=\widehat{t}+Y(𝐫)$$
(4)
and

$$\widehat{t}=\frac{\mathrm{}^2^2}{2m}\lambda .$$
(5)
In the above equations $`\lambda `$ is the chemical potential, $`W(𝐫)`$ is an external potential representing a nucleus, $`S(𝐫)`$ is a self-consistent potential generated by the neutrons themselves while $`Y(𝐫)`$ stands for the total field.

By selective averaging, or, more formally, through applying the saddle-point approximation to the grand canonical partition function which follows from this Hamiltonian, one can derive a mean field Hamiltonian

$$H_{\mathrm{MFA}}=H_0+d^3r[\psi _{}^{}(𝐫)\psi _{}^{}(𝐫)\mathrm{\Delta }(𝐫)+\mathrm{\Delta }^{}(𝐫)\psi _{}(𝐫)\psi _{}(𝐫)]+\zeta $$
(6)
with

$$\mathrm{\Delta }(𝐫)=|g|\psi _{}(𝐫)\psi _{}(𝐫)$$
(7)
being sometimes called “the complex order parameter”, “the pairing potential”, or “the pair function”. Since this is a mean-field approximation, quantum fluctuations in the order parameter is neglected. The quantity

$$\zeta =d^3r\frac{|\mathrm{\Delta }(𝐫)|^2}{|g|}$$
(8)
is a c-number which plays no role in the free energy minimization leading to the Bogoliubov-de Gennes equations, but which is important in the energy calculations. Here, and in the following $`\mathrm{}`$ denotes thermal averages. The mean field Hamiltonian (6) can be diagonalized by a canonical Bogoliubov transformation

$`\alpha _{i,}(t)`$ $`=`$ $`{\displaystyle d^3r[U_i^{}(𝐫)\psi _{}(𝐫,t)V_i^{}(𝐫)\psi _{}^{}(𝐫,t)]}`$ (9)
$`\alpha _{i,}^{}(t)`$ $`=`$ $`{\displaystyle d^3r[U_i(𝐫)\psi _{}^{}(𝐫,t)+V_i(𝐫)\psi _{}(𝐫,t)]},`$ (10)
with the inverse transformation being

$`\psi _{}(𝐫,t)`$ $`=`$ $`{\displaystyle \underset{i}{}}[U_i(𝐫)\alpha _{i,}(t)+V_i^{}(𝐫)\alpha _{i,}^{}(t)]`$ (11)
$`\psi _{}^{}(𝐫,t)`$ $`=`$ $`{\displaystyle \underset{i}{}}[U_i^{}(𝐫)\alpha _{i,}^{}(t)V_i(𝐫)\alpha _{i,}(t)].`$ (12)
In these equations, $`i`$ is a label for the single-particle eigenstates. The coefficients $`U`$ and $`V`$ of this Bogoliubov transformation must obey the canonical relations

$`{\displaystyle \underset{i}{}}[U_i(𝐫)U_i^{}(𝐫^{})+V_i(𝐫^{})V_i^{}(𝐫)]`$ $`=`$ $`\delta (𝐫𝐫^{})`$ (13)
$`{\displaystyle \underset{i}{}}[U_i(𝐫)V_i^{}(𝐫^{})U_i(𝐫^{})V_i^{}(𝐫)]`$ $`=`$ $`0`$ (14)
$`{\displaystyle d^3r[U_i(𝐫)U_i^{}^{}(𝐫)+V_i(𝐫)V_i^{}^{}(𝐫)]}`$ $`=`$ $`\delta _{ii^{}}`$ (15)
$`{\displaystyle d^3r[U_i(𝐫)V_i^{}(𝐫)U_i^{}(𝐫)V_i(𝐫)]}`$ $`=`$ $`0,`$ (16)
in order for the quasiparticle operators to follow Fermi statistics. Carrying out this transformation, the Hamiltonian takes the simple form

$$H_B=U_0+\underset{i}{}E_i(\alpha _{i,}^{}\alpha _{i,}+\alpha _{i,}^{}\alpha _{i,}).$$
(17)
The constant term $`U_0`$ is the ground state energy of the system. The Bogoliubov-de Gennes equations for the transformation coefficients which bring the Hamiltonian to the diagonal form $`H_B`$ follow from the Heisenberg equations of motion for the field operators, $`i\stackrel{.}{\psi }=[\psi ,H_B]`$, that is, in matrix form

$$i\frac{}{t}\left(\begin{array}{c}\psi _{}(𝐫,t)\\ \psi _{}(𝐫,t)\end{array}\right)=\left(\begin{array}{cc}\widehat{t}+Y(𝐫)& \mathrm{\Delta }(𝐫)\\ \mathrm{\Delta }^{}(𝐫)& \widehat{t}Y(𝐫)\end{array}\right)\left(\begin{array}{c}\psi _{}(𝐫,t)\\ \psi _{}(𝐫,t)\end{array}\right).$$
(18)
The $`\psi `$-operators are not eigenvectors of the Hamiltonian and thus have no definite frequency. We can, though, express them through the eigenoperators of the Hamiltonian, $`\alpha `$, $`\alpha ^{}`$:

$$\alpha _{i,}(t)=\alpha _{i,}e^{iE_qt};\alpha _{i,}^{}(t)=\alpha _{i,}^{}e^{iE_qt}.$$
(19)
By doing so, and gathering the terms with $`\alpha _{q,}`$ and so on, one obtains the BdG equations

$$\left(\begin{array}{cc}\widehat{t}+Y(𝐫)& \mathrm{\Delta }(𝐫)\\ \mathrm{\Delta }^{}(𝐫)& \widehat{t}Y(𝐫)\end{array}\right)\left(\begin{array}{c}U_i\\ V_i\end{array}\right)=E_i\left(\begin{array}{c}U_i\\ V_i\end{array}\right).$$
(20)
Solving these equations, one gets both the amplitudes $`U`$ and $`V`$ along with the excitation energies $`E_q`$ of the system. Once the amplitudes are known, the expectation value of any single-particle operator $`𝒪`$ can be obtained from

$$𝒪=2\underset{E_i>0}{}d^3r[U_i^{}𝒪U_if(E_i)+V_i𝒪V_i^{}(1f(E_i))],$$
(21)
where $`f(E_i)=(1+\mathrm{exp}(E_i/k_BT)^1`$ is the Fermi distribution function, $`k_B`$ is Boltzmann’s constant and $`T`$ is the temperature. More explicitly, we write the BdG equations as

$`\left({\displaystyle \frac{\mathrm{}^2^2}{2m}}+Y(𝐫)E_F\right)U_i(𝐫)+\mathrm{\Delta }(𝐫)V_i(𝐫)`$ $`=`$ $`E_iU_i(𝐫)`$ (22)
$`\left({\displaystyle \frac{\mathrm{}^2^2}{2m}}+Y(𝐫)E_F\right)V_i(𝐫)+\mathrm{\Delta }^{}(𝐫)U_i(𝐫)`$ $`=`$ $`E_iV_i(𝐫)`$ (23)
where $`U`$ and $`V`$ are the quasiparticle amplitudes, $`E_F`$ is the Fermi energy, $`\mathrm{\Delta }(𝐫)`$ is the (space-dependent) pairing potential and $`m`$ is the neutron mass. In the above equations we have replaced the chemical potential with the Fermi energy. The resulting violation in particle number conservation is found to be very small. In equations (22,23) the subscript $`i`$ represents all relevant quantum numbers. The pairing potential has to be calculated self-consistently as

$$\mathrm{\Delta }(𝐫)=|g|\underset{i;0<|E_i|<\mathrm{}\mathrm{\Omega }}{}U_i(𝐫)V_i^{}(𝐫)(12f(E_i))$$
(24)
where the sum is over quasiparticle states having energy eigvenvalues with absolute values $`|E_i|`$ smaller than a cutoff $`\mathrm{}\mathrm{\Omega }`$, in our calculations taken to be 10 MeV. In general the set of quantum numbers $`i`$ includes the principal quantum number, the projection of the angular momentum $`\mu =\mathrm{},3/2,1/2,1/2,3/2,\mathrm{}`$ and the wavenumber $`k_z`$ parallel to the $`z`$-axis. As we want to describe neutron vortices with cylindrical symmetry, we write the quasiparticle states as

$`U_i(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{L}}}u_{n\mu k_z}(\rho )\mathrm{exp}[i(\mu 1/2)\theta ]\mathrm{exp}[i(k_zz)]`$ (25)
$`V_i(𝐫)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{L}}}v_{n\mu k_z}(\rho )\mathrm{exp}[i(\mu +1/2)\theta ]\mathrm{exp}[i(k_zz)]`$ (26)
where $`\rho `$, $`\theta `$, $`z`$ are the cylindrical coordinates, $`L`$ is the length of the cylinder, $`n`$ is a radial quantum number and $`\mu `$ is half an odd integer.

The contributions to the pairing potential are calculated in each subspace of fixed $`\mu `$ and $`k_z`$ and are then summed to give the total pairing potential. In Paper I we considered the simplified case of a pairing potential calculated at $`k_z=0`$. This approximation is the one currently used in self-consistent calculations of vortices in type-II superconductors, where the electron effective mass is anisotropic and increases strongly along directions perpendicular to the vortex axis (Gygi & Schlüter (1991)). In the case of vortices in neutron matter, variations in the effective mass are expected to be much smaller. In the present study we shall therefore go beyond this approximation and allow for different values of the wavenumber along the vortex axis. Since a direct comparison of vortices in neutron matter and in superconductors is important, we shall also refer to the case with $`k_z=0`$ whenever such a comparison can be revealing.

The angular dependence in equations (25,26) of the quasiparticle states follows from imposing a pairing gap of the form

$$\mathrm{\Delta }(𝐫)=\mathrm{\Delta }(\rho )\mathrm{exp}[i\theta ]$$
(27)
representing a vortex with one quantum of circulation. We are not interested in larger winding numbers, since they are expected at much higher energies. Following (Gygi & Schlüter (1991)) we expand the quasiparticle states in terms of cylindrically symmetric Bessel functions and impose the boundary condition $`\mathrm{\Delta }(R)=0`$ at the edges of a cylinder of radius $`R`$, where $`R`$ is typically $`50`$-$`100\mathrm{fm}`$. More specifically, the basis functions are chosen as

$$\varphi _{jm}(\rho )=\frac{\sqrt{2}}{RJ_{m+1}(\alpha _{jm})}J_m\left(\alpha _{jm}\frac{\rho }{R}\right),j=1,\mathrm{},N,$$
(28)
where $`m=\mu \pm \frac{1}{2}`$ is an integer and $`\alpha _{jm}`$ is the $`j`$th zero of the Bessel function $`J_m(x)`$. The dimension $`N`$ of the basis is chosen large enough to ensure convergence and stability of the quantities of interest. The quasiparticle amplitudes are expanded as ($`i=\{n\mu k_z\}`$)

$`u_i(\rho )`$ $`=`$ $`{\displaystyle \underset{j}{}}c_{ij}\varphi _{j\mu \frac{1}{2}}(\rho )`$ (29)
$`v_i(\rho )`$ $`=`$ $`{\displaystyle \underset{j}{}}d_{ij}\varphi _{j\mu +\frac{1}{2}}(\rho ).`$ (30)
We consider first the case without mean field, which is the one usually addressed in the calculations for type-II superconductors. For given values of $`\mu `$ and $`k_z`$ equations (22,23) can be written as a $`2N\times 2N`$ matrix eigenvalue problem

$$\left(\begin{array}{cc}T^{}& \mathrm{\Delta }\\ \mathrm{\Delta }^T& T^+\end{array}\right)\mathrm{\Psi }_n=E_n\mathrm{\Psi }_n,$$
(31)
where the superscript $`T`$ denotes the transpose of a matrix, $`\mathrm{\Psi }_n^T=(c_{n1},\mathrm{},c_{nN},d_{n1},\mathrm{},d_{nN})`$,

$$T^\pm =\frac{\mathrm{}^2}{2m}\left(\frac{\alpha _{j\mu \pm 1/2}^2}{R^2}+k_z^2k_F^2\right)\delta _{jj^{}},$$
(32)
and the matrix $`\mathrm{\Delta }`$ is given by

$$\mathrm{\Delta }_{jj^{}}=_0^R\varphi _{j\mu 1/2}(\rho )\mathrm{\Delta }(\rho )\varphi _{j^{}\mu +1/2}(\rho )\rho 𝑑\rho .$$
(33)
In cases where a mean field is present the matrix elements in equation (32) get an additional term. In our calculations we assume that the mean field of a nucleus has a Woods-Saxon shape,

$$W(\rho )=\frac{W_0}{1+\mathrm{exp}\left(\frac{\rho R_N}{a}\right)}$$
(34)
where typical values of the constants are $`W_030`$ to $`40\mathrm{MeV}`$, $`R_N5\mathrm{fm}`$, and $`a0.6\mathrm{fm}`$. Note that we only consider cylindrical nuclei in the present paper. The case of a spherical nucleus interacting with a cylindrical vortex line is quite complicated since the two objects have different symmetries. However, it is worth noting that several calculations (Lorenz, Ravenhall & Pethick (1993); Pethick & Ravenhall (1995)) predict that nuclei take on cylindrical shapes in parts of the inner crust of a neutron star, so our results have in fact a direct bearing on the actual problem.

Even without a nucleus present, the neutrons move in a self-consistent Hartree field of the form

$$S(\rho )=|g|n(\rho )$$
(35)
where $`n(\rho )`$ is the neutron density distribution within the cylinder, given by

$$n(\rho )=2L^1\underset{n,\mu ,k_z;0<E_{n,\mu ,k_z}<\mathrm{}\mathrm{\Omega }}{}|v_{n,\mu ,k_z}(\rho )|^2.$$
(36)
In the presence of a mean field with cylindrical symmetry, $`Y(\rho )=W(\rho )+S(\rho )`$, the matrix elements of the one-body operator in equation (32) acquire an additional term of the form

$$W{}_{jj^{}}{}^{\pm }=_0^Rd\rho \rho \varphi _{j\mu \pm 1/2}(\rho )Y(\rho )\varphi _{j^{}\mu \pm 1/2}(\rho ).$$
(37)
For generality the expressions have been written at finite temperature, making use of the Fermi occupation factors for quasiparticle excitations. Although the case at finite temperature presents no additional difficulty, in applications to neutron star superfluids the temperature is essentially zero compared to the Fermi temperature. Thus, all numerical calculations in this paper are carried out at $`T=0`$. One can use the symmetries of the BdG equations under the transformation

$$(U,V,E)(V^{},U^{},E),$$
(38)
i.e., if $`(U,V,E)`$ is a solution, then so is $`(V^{},U^{},E)`$, to rewrite equation (24) simply as

$$\mathrm{\Delta }(𝐫)=2|g|\underset{i;0<E_i<\mathrm{}\mathrm{\Omega }}{}U_i(𝐫)V_i^{}(𝐫)=2\frac{|g|}{L}e^{i\theta }\underset{n,\mu ,k_z;0<E_{n,\mu ,k_z}<\mathrm{}\mathrm{\Omega }}{}u_{n,\mu ,k_z}(\rho )v_{n,\mu ,k_z}(\rho )$$
(39)
where the sum now extends over positive eigenstates of the BdG equations only. This symmetry can be in principle exploited to solve the system (31) for only positive (or only negative) eigenstates.

A point-like pairing interaction $`v(𝐫,𝐫^{})=|g|\delta (𝐫𝐫^{})`$ has been used in previous work on superconductors (de Gennes (1989); Gygi & Schlüter (1991); Hayashi et al. (1998)). In the case of neutron matter this approximation should work whenever the range of the inter-particle interaction (which is of the order $`1\mathrm{fm}`$) is at least comparable to or smaller than the nucleus and the vortex. Exactly as for a system of electrons, a cutoff energy $`\mathrm{}\mathrm{\Omega }`$ has to be introduced for the sum in equation (39) to converge.

## 3 Results

### 3.1 Calculations for different gaps

One of the parameters still largely unknown in the physics of neutron stars is the value of the neutron pairing gap as a function of the density in the crust and in the interior. The reason is that the gap is very sensitive to the value of the parameters defining the neutron-neutron interaction at the Fermi surface, as can be seen from the weak-coupling formula for the gap

$$\mathrm{\Delta }2\mathrm{}\mathrm{\Omega }\mathrm{exp}\left(\frac{1}{N(0)\stackrel{~}{V}(k_F,k_F)}\right)$$
(40)
where $`\mathrm{}\mathrm{\Omega }`$ is an energy cutoff (often taken to be the Fermi energy), $`N(0)`$ is the density of states for a fixed spin projection at the Fermi surface and $`\stackrel{~}{V}`$ is the effective interaction between two Cooper-paired neutrons. In many studies (Baldo et al. 1990; Elgarøy & Hjorth-Jensen 1998; Khodel, Khodel & Clark 1996) the last quantity is approximated by the free neutron-neutron interaction, taken from various nucleon interaction models fitted to proton-proton and neutron-proton scattering data. The weak-coupling approximation is poor in neutron matter due to the strong momentum-dependence of the free neutron-neutron interaction, and the BCS gap equation has to be solved in its full complexity (Baldo et al. (1990)). While sometimes numerically tricky, this task can now be carried out accurately, and the results turn out to be insensitive to the model chosen for the nucleon-nucleon interaction as long as it fits the scattering data accurately (Khodel et al. (1996); Elgarøy & Hjorth-Jensen (1998)). However, the interaction between two neutrons is modified when other neutrons are present, and an important contribution comes from the so-called induced interaction, corresponding to the exchange of density and spin-density fluctuations between two neutrons (Migdal (1967)). The effect of this has been calculated by Wambach, Ainsworth, & Pines(1993) and Schulze et al. (1996) with somewhat different results. Wambach et al. (1993), calculating the induced interaction accurately, but using the weak-coupling approximation for the gap, predict a decrease by a factor of $`3`$ at all densities compared with the calculations using free nucleon-nucleon interactions. Schulze et al. (1996), solving the full gap equation, but using a somewhat rougher treatment of the induced interaction, predict a similar reduction of the gap at low densities, but an increase at higher densities. We know of no estimates for contributions beyond the induced interaction; some of these may turn out to be important. Given the uncertainties in the predictions for gaps in uniform neutron matter, we will in some of our calculations let this be a free parameter and study how other quantities depend on the size of the gap. In other calculations we will use the results of Wambach et al. (1993), since these are commonly employed in other neutron star studies, and have so far been consistent with the glitch data. Actually, we find the nice feature that some important properties of the vortex are not too sensitive to the gap. In a first group of calculations we shall thus keep the value of the gap free, using several values of the pairing force $`|g|`$ for a fixed Fermi wave number.

Keeping comparable values of the gap at infinity and changing $`k_F`$ we find that $`|g|`$ is approximately inversely proportional to the Fermi wave number. This is in rough agreement with the weak coupling formula which predicts (with $`\stackrel{~}{V}(k_F,k_F)=|g|`$)

$$|g|k_F=\frac{2\pi ^2\mathrm{}^2}{m}\frac{1}{\mathrm{ln}\left(\frac{\mathrm{\Delta }}{2\mathrm{}\mathrm{\Omega }}\right)},$$
(41)
i.e., $`|g|k_F\mathrm{constant}`$ modulo a slowly varying logarithmic factor.

Figure 1 shows the pairing gap as a function of the distance from the vortex core. As analyzed in Paper I, the gap increases from zero to an asymptotic value $`\mathrm{\Delta }_{\mathrm{}}`$, which is the value in infinite homogeneous matter. From the variation of the gap we can extract the appropriate size of the vortex. As in Paper I, we take the intersection of the tangent line at the origin $`\rho =0`$ with the line $`\mathrm{\Delta }_{\mathrm{}}`$ to define the length $`\xi _2`$. The quantities $`\xi _{50}`$ and $`\xi _{90}`$ are the distances from the axis where the gap reaches $`50\%`$ and $`90\%`$ of the value at infinity, respectively. The length $`\xi _0`$ is simply the BCS coherence length, i.e. the size of a Cooper pair. In Paper I we showed that with one single $`k_z=0`$ state the size of a vortex tends to change smoothly as a function of the pairing force $`|g|`$ and the Fermi wave number $`k_F`$. Table 1 of the present paper collects data for the vortex size for several different Fermi wave numbers and asymptotic gaps. The data show that the vortex size, although certainly not constant, has a limited range of variation, while the BCS coherence length changes dramatically. For a fixed $`k_F`$ and increasing $`\mathrm{\Delta }_{\mathrm{}}`$ the lengths generally tend to decrease. The reason for the different scaling of the vortex core compared with the BCS coherence length will be examined in Sect. 3.2 while consequences of astrophysical relevance will be discussed in Sect. 3.6.

Figure 2 shows the eigenvalues of the vortex state for the case of one single $`k_z`$ mode. Similar to that found in other systems exhibiting fermionic superfluidity, such as superconductors (Ketterson & Song (1999)) or finite nuclei (Schuck& Taruishi (1996)), a window of width $`2\mathrm{\Delta }_{\mathrm{}}`$ opens up where only bound states are present. These are visible as a branch that for large angular momenta $`\mu `$ approaches the energy of continuum states, very close to what one finds for vortices in type II superconductors (Gygi & Schlüter (1991)). The smallest eigenvalue has an energy of the order $`\mathrm{\Delta }_{\mathrm{}}^2/E_F`$. For more $`k_z`$ states, as shown in figure 3, there are in general more bound states for a given angular momentum $`\mu `$ and the energy is roughly $`E\mu \mathrm{\Delta }_{\mathrm{}}^2/E_F[1(k_z/k_F)^2]^{1/2}`$ (Caroli, de Gennes & Matricon (1964)). The eigenfunctions of the bound states are confined in the region of the vortex core and decrease exponentially far from the core. The wavefunctions of the continuum oscillate like scattering states with a radial wave number $`k_\rho `$ determined by the condition $`E_n^2=\mathrm{\Delta }_{\mathrm{}}^2+(\mathrm{}^2/2m)^2(k_\rho ^2+k_z^2k_F^2)^2`$.

### 3.2 Calculations with fixed gap at infinity

In a second group of calculations we assume that the pairing gap of neutron matter is given by the results of (Wambach et al. (1993)). In short, once the neutron Fermi wavenumber is fixed, we choose the value of the pairing strength $`g`$ which reproduces the given value of the pairing gap. This procedure is necessary when a nucleus is present in the center of the vortex core, because a nucleus modifies the local density. This in turn has strong influence on the local gap and condensation energy density, which are both strongly density–dependent.

We solve the BdG equations for a homogenous system at various densities, and require that the energy gaps thus obtained be equal to the ones of Wambach et al. This provides us with a density dependent pairing strength $`|g|=|g(n(\rho ))|`$. When solving the BdG equations for a vortex line, we naturally find that the neutron density varies within the cylinder. Thus, $`g`$ will also vary throughout the vortex, and this must be taken into account. The pairing potential is then modified to

$$\mathrm{\Delta }(𝐫)=2\underset{i;0<E_i<\mathrm{}\mathrm{\Omega }}{}|g(𝐫)|U_i(𝐫)V_i^{}(𝐫),$$
(42)
and similarly for other quantities involving $`|g|`$.

In the numerical solution of the BdG equations, we use $`50100`$ Gaussian mesh points for the radial coordinate, and some $`50100`$ angular momentum states and $`810`$ plane wave states (for the $`z`$ direction) in the expansion of the $`U`$ and $`V`$ amplitudes. Starting from initial approximations to $`\mathrm{\Delta }(\rho )`$ and $`n(\rho )`$ the BdG equations are solved by diagonalizing the resulting eigenvalue problem. This gives us approximations to $`U`$ and $`V`$, and from these new approximations to $`\mathrm{\Delta }(\rho )`$ and $`n(\rho )`$ are obtained. We iterate this procedure until the variations in these quantities are small from one iteration to the next. Usually this procedure converges after 5-10 iterations, and it is adequate to use $`\mathrm{\Delta }(\rho )=\mathrm{constant}`$ and $`n(\rho )=\mathrm{constant}`$ as initial approximations. Figures 4 and 5 show the pairing potential $`\mathrm{\Delta }(\rho )`$ and the neutron number density $`n(\rho )`$ respectively as functions of the distance from the axis for $`k_F=0.8\mathrm{fm}^1`$ in four different configurations: an isolated vortex at $`\rho =0`$, an isolated nucleus, vortex and nucleus both present, and uniform matter. When only the vortex is present the gap rises linearly to an asymptotic value within lengths specified by the parameters $`\xi `$ defined in Section 3.1. The fact that the pairing strength is now space-dependent due to the decrease of the density in the core visible in figure 5, changes only slightly the values of the vortex core radius. The shape of the gap is also only slightly altered by the variation of the pairing strength. The case where only the mean field is present shows a decrease of the pairing gap in the region occupied by the nucleus, qualitatively consistent with semiclassical approximations (Broglia et al. (1994)). When both a nucleus and vortex are present the picture becomes more complex due to the strong distortion of the velocity field induced by the nucleus. To understand the difference in behavior of the gap in the presence of a vortex, with or without the nuclear mean field, we analyze the problem using the Ginzburg-Landau (GL) equation

$$\frac{\mathrm{}^2^2}{4m}\psi (𝐫)+W(𝐫)\psi (𝐫)A\psi (𝐫)B|\psi (𝐫)|^2\psi (𝐫)=0$$
(43)
where $`\psi (𝐫)`$ is the order parameter and the self-consistent mean field $`S(\rho )`$ is neglected compared to $`W(\rho )`$, which is a very reasonable approximation. The parameter $`A`$ is negative for $`T<T_C`$, where $`T_C`$ is the critical temperature of the superfluid, while $`B`$ is always positive. Using cylindrical coordinates, taking $`\psi (𝐫)=e^{i\theta }f(\rho )`$, $`W(𝐫)=W(\rho )`$, and using units where $`\mathrm{}^2/4m=1`$, equation (43) can be written as

$$f^{\prime \prime }(\rho )+\frac{1}{\rho }f^{}(\rho )+\left[W(\rho )AB|f(\rho )|^2\frac{1}{\rho ^2}\right]f(\rho )=0.$$
(44)
The boundary condition at $`\rho =0`$ is $`f(0)=0`$. Looking at small $`\rho `$, we neglect the non-linear term, and distinguish between two cases. First, for the case of no nuclear mean field, equation (44) becomes

$$f^{\prime \prime }(\rho )+\frac{1}{\rho }f^{}(\rho )\frac{1}{\rho ^2}f(\rho )=0,$$
(45)
since $`A`$ is supposed to be small near $`T_C`$ and thus negligible compared to $`1/\rho ^2`$ for small $`\rho `$. This equation is easily seen to have the solution $`f(\rho )\rho `$. Thus, we expect the pairing potential to show a linear increase for small $`\rho `$ in the case of an isolated vortex, and this behaviour is also seen in our microscopic calculations. Looking at the second case, an attractive mean field $`W(\rho )W_0`$, assumed to be constant for the limited range in $`\rho `$ we are looking at now, equation (44) is reduced to

$$f^{\prime \prime }(\rho )+\frac{1}{\rho }f^{}(\rho )+\left(W_0\frac{1}{\rho ^2}\right)f(\rho )=0,$$
(46)
where we have absorbed the constant $`A`$ in the mean field $`W_0`$. Introducing a new variable through $`\rho =\frac{1}{\sqrt{W_0}}x`$, we obtain

$$x^2f^{\prime \prime }(x)+xf^{}(x)+(x^21)f(x)=0,$$
(47)
which the reader may recognize as the Bessel equation of order 1. The solution satisfying $`f(0)=0`$ is $`f(x)J_1(x)`$. Thus, in the case of a vortex interacting with a nucleus, we may expect the pairing potential to show oscillations at small $`\rho `$.

Figure 6 shows the relevant lengths for the case of polarized gap. The lengths remain limited ($`<15`$ fm) for essentially all densities and tend to decrease at higher density. This behavior is different from the BCS coherence length $`\xi _0`$ which scales like $`k_F/\mathrm{\Delta }_{\mathrm{}}`$. In the realistic case where the pairing strength is fitted to the gap in uniform matter, the coherence length $`\xi _2`$ decreases as a function of $`k_F`$. Thus, the density dependence of the pairing strength $`|g|`$ is also very important, and the vortex size appears to depend in a quite complex way on both $`|g|`$ and the Fermi wavenumber. In the figure we also show the analytic estimate of the vortex core size made by Kramer & Pesch (1974) whose value in our units (energies in MeV and lengths in fm) is

$$\xi _1=\frac{661.3}{k_F^2|g|}.$$
(48)
Since the pinning force between a nucleus and a vortex behaves like $`\mathrm{\Delta }^2k_F/\xi _2`$, an increase of the pinning force of the order $`\xi _0/\xi _2`$ can be expected. In addition we can confirm what we found in Paper I, that since the diameter of the vortex is always smaller than the lattice spacing, the vortex envelops at most one nucleus along a plane perpendicular to its axis.

Another interesting parameter is represented by the coherent flux of neutrons in the presence of a vortex. The flux is given by

$$J(𝐫)=\frac{\mathrm{}}{2mi}\left(\underset{i}{}\left[V_i(𝐫)V_i^{}(𝐫)V_i^{}(𝐫)V_i(𝐫)\right](1f(E_i))+\left[U_i(𝐫)U_i^{}(𝐫)U_i^{}(𝐫)U_i(𝐫)\right](1f(E_i))\right)$$
(49)
which in the present case becomes

$$J(\rho )=\frac{\mathrm{}}{2m\rho }[\underset{E_i<0}{}(\mu 1/2)(|u_i(\rho )|^2\underset{E_i>0}{}(\mu +1/2)|v_i(\rho )|^2].$$
(50)
Figure 7 shows the flux of neutrons with and without a nucleus present in the middle of the cell for the case of $`k_F=0.4`$ fm<sup>-1</sup>. The main effect of the nucleus is to decrease the flux. The coherent velocity field can be obtained dividing the flux by the neutron density. In this case the velocity is even more depleted in the centre due to the much higher density in the nuclear region than outside, a situation typical for low values of $`k_F`$.

### 3.3 Variations in $`\mathrm{\Delta }`$ due to a center of scattering

Pinning or antipinning between a vortex and a nucleus is due to local variations in both the kinetic energy density and the condensation energy. Although pinning will be examined in detail in a next section, it is interesting to study how the pairing gap changes in the cell when a nucleus alone is present. Some curves have already been presented in Figure 1, but here we shall be more systematic discussing of the applications to pinning calculations. The variation of the pairing gap as a function of the distance from the axis of a cylindrical nucleus when the vortex is absent can be seen in Figure 8 for some selected densities. Note that the presence of a mean field decreases the value of the gap. This is mainly due to the increase of the density in the nuclear region, since the pairing force $`|g|`$ depends on the local density. Note also that the variation of the gap extends beyond the range of the mean field, indicated with a dashed line. This is due in part to the long tail of the density distribution beyond the radius of the potential, an effect due to the presence of scattering states. Secondly, there is a proximity effect between the nuclear region and the free neutron gas due to Andreev scattering (that is, scattering due to spatial variations the pairing potential) of quasi-particles.

### 3.4 Pinning energy and vortex self-energy

The vortex self-energy or tension $`T`$, which is the energy carried by the vortex per unit length, is a relevant parameter for vortex-nucleus pinning. Let us consider the behavior of a single vortex line in interaction with the whole lattice. If the energy necessary for the vortex to reach the pinning centers is very high compared to the energy gained by pinning, the vortex responds stiffly to local deformations. The importance of this parameter can be grasped in a limiting case: for $`T\mathrm{}`$ the vortex behaves like a rigid cylinder. For a vortex moving perpendicularly to its axis the total pinning force with the lattice arises from stochastic summation from all the pinning centers and is proportional to the square root of the vortex length (see for example Tinkham (1996) for the case of superconductors). Thus, the pinning force cannot balance the Magnus force, which increases linearly with the vortex length. If this was the case, the pinning mechanism would not be effective for storing the energy released during a glitch (Anderson & Itoh (1975)). For a finite stiffness, the situation is more complicated. Solving the equation of motion for a vortex line passing through many centers of pinning would give the value of the total pinning energy which is able to be stored, but this is a difficult task due to the difficulty in finding the appropriate boundary conditions in the presence of many centers of pinning. It has been argued that although finite, the tension is too large to account for large glitches (Jones (1998)). If, on the other hand, the pinning model is correct, the value of the vortex tension is relevant for the dynamics of unpinning (Link & Epstein (1995)) and for the value of the effective pinning force per unit length.

The vortex tension is the energy difference per unit length between the configuration with and without a vortex, $`T=\left[E(\mathrm{vortex})E(\mathrm{uniform}\mathrm{matter})\right]/L.`$ To calculate the energy of the superfluid in each configuration one can take averages of Eq. (6). The result is

$$E=d^3r\frac{|\mathrm{\Delta }(𝐫)|^2}{|g|}+2\underset{E_n<0}{}E_n2\underset{E_n}{}E_nd^3r|V_n(𝐫)|^2.$$
(51)
Table (2) shows the calculated points as a function of the Fermi wave number compared to the classical calculation where the kinetic energy density is integrated from the coherence length $`\xi `$ up to an upper limit $`R_0`$:

$$T=4.38k_F^3\mathrm{ln}\left(\frac{R_0}{\xi }\right).$$
(52)
Evidently the presence of bound states does not change the vortex tension radically. Furthermore, constraining the vortex in a small cylinder can overestimate the role of bound states in the calculation of self-energy, because the vortex is macroscopic along the $`\rho `$ direction. We also expect that uncertainties in the pairing gap or in the cutoff length $`R_0`$ (that should be of the order of the vortex-vortex separation) may be much greater than the discrepancy we find between our microscopic calculations and the estimate (52). We conclude that equation (52) is reliable when applied to vortex calculations in a neutron star. Assuming a vortex-vortex separation of the order $`R_010^2`$ cm, this formula implies a tension of the order $`T10^9`$ erg cm<sup>-1</sup> at $`k_F=0.8`$ fm<sup>-1</sup>. At such values a vortex can deform only at length scales of the order of several hundreds of fm, nearly $`20`$ times larger than the lattice spacing, unless the pinning energy assumes unrealistically large values. How a vortex line with such small deformability can adapt its shape to gain a net pinning energy is still an unsolved problem (Jones (1997)). We show in the following that an answer may be provided by the high value of the elementary pinning force.

To calculate the pinning energy we need to consider four different configurations: a vortex superimposed on a nucleus, an isolated nucleus, an isolated vortex, and finally uniform matter. For each configuration, the energy needs to be calculated with Eq. (51). The pinning energy per unit length can be calculated as

$$E_p=\left[E(\mathrm{nucleus}\mathrm{and}\mathrm{vortex})E(\mathrm{nucleus})+E(\mathrm{uniform}\mathrm{matter})E(\mathrm{vortex})\right]/L.$$
(53)
Table (3) collects some pinning data calculated with Eq. (49) and (51). A negative sign corresponds to a pinning situation in which it is energetically favorable for a vortex to remain attached to a nucleus, while a positive sign indicates an antipinning situation, where vortices are repelled by nuclei. We find an antipinning-pinning transition at $`k_F=0.300.40`$ fm<sup>-1</sup>. A similar transition was found in a Ginzburg-Landau approach (Epstein & Baym (1988)). The pinning energy reaches values of the order of several MeV per unit length, corresponding to a very high value if compared to previous calculations performed solely on the basis of the condensation energy loss. As is also found in other models where the kinetic energy is accounted for (Epstein & Baym (1988)), part of these large values can be attributed to kinetic energy effects. For such large values of the pinning energy a vortex is likely to deform sensitively in order to be able to catch the pinning centers. In addition, the ions in the lattice will shift to the vortex in response to the high pinning force, a situation referred to as superstrong pinning. Since nuclei in most of the crust are spherical, one might question the applicability of our calculations, where nuclei are considered to be cylindrical. Although the proper use of spherical geometry would, of course, quantitatively change the results, we do not expect dramatic differences.

## 4 Conclusions

The results presented in this paper could be a step towards an understanding of vortex structure and interaction in neutron star crusts. We have found differences between the predictions of microscopic theory and those of more macroscopic models. In particular, we showed that the BCS coherence length cannot be used as a reliable estimate of the vortex core size, especially not at low and intermediate densities. By solving the self-consistent problem of a vortex in neutron matter, we have shown that the vortex core has a much more complicated structure than that which is usually assumed in simplified treatments of the problem, e.g.: 1) the gap changes from zero at the axis of the core up to an asymptotic value, and 2): there are bound states present.

We found that vortices in neutron star matter follow a behavior similar to that described by Kramer & Pesch (1974) for superconductors, namely, at zero temperature the vortex core size scales as the inverse of the Fermi wave number. In a superconductor in clean limit this effect is associated with a linear increase of the vortex core size with temperature. We have not investigated the temperature dependence of our results, because this is not very relevant in neutron star physics due to the low temperatures compared with the Fermi energy. Unfortunately, experimental data for clean type II superconductors are limited to the quite complex new cuprate systems. Scanning tunnelling microscopy has, however, confirmed the presence of bound states of quasiparticles in the vortex core, in good agreement with the Bogoljubov-de Gennes theory. In view of the small range of the nuclear force and of the simple structure of the Fermi surface, the BdG model could be a very good approximation for neutron star matter.

High pinning energies and small lengths imply a large pinning (or anti-pinning) force. In addition, the small value of the coherence length means that a vortex envelopes no more than one nucleus in a plane perpendicular to the vortex axis. The possibility for a vortex to store pinning energy is largely the result of a competition between the energy gain with pinning and the energy cost to deform. To be able to deform and catch a nucleus along its path, a vortex has to be relatively free to perform local deformations. Assuming that the energy of the vortex is solely due to kinetic effects, the tension turns out to be of the order $`T100k_F^3`$ MeV fm<sup>-1</sup> and is thus very high at densities corresponding to the deep inner crust. For a total pinning energy per nucleus $`E_p`$, a vortex can efficiently catch the pinning sites only if they are at least a distance $`a^2T/2`$ apart. For example, at $`k_F=1`$ fm<sup>-1</sup> a vortex deformation of the order of the lattice spacing can occur every $`8`$ or so lattice spacings, while at $`k_F=0.8`$ fm<sup>-1</sup> the pinning efficiency would be much better. Calculations of the geometry of vortex deformation in a one-dimensional crystalline structure have been performed by (Link & Epstein (1995)).

Finally, there are other unsolved problems in neutron star dynamics where a microscopic description of vortex states might prove helpful. The scattering of electrons off neutron vortices results in a coupling between protons and neutrons in the interior of the star. This is a possible coupling mode between the superfluid neutrons and the rest of the star, and depends on the microscopic structure of vortex lines (Sauls (1989)).

The proton vortices in the interior of the star induced by the rotation will probably be similar to the neutron vortices examined in this paper, but additionally they will generate a magnetic field parallel to the rotational axis. How these magnetic vortices interact with the (much more numerous) flux tubes induced by the high magnetic field penetrating into the star is only one of the problems connected to the dynamics of superfluid vortices below the star’s crust. Despite their importance in the evolution of the pulsar’s magnetic field, the properties of proton flux tubes and their possible interaction with neutron vortices in the interior, where neutrons couple in a triplet phase, are not well understood (Ruderman (1997)). A microscopic description of proton flux lines might also be applied to the study uf ultra-high magnetic fields and the evolution of magnetic field in magnetars.

###### Acknowledgements.

We would like to thank M. Hjorth Jensen, E. Osnes, L. Engvik and T. Engeland for illuminating discussions on the neutron star equation of state and on superfluidity in both neutron stars and finite nuclei, A. Sedrakian and C. Pethick for useful comments on the present model and M. Colpi for an interesting discussion about strong magnetic fields in collapsed stars. F. V. De Blasio was supported by the Marie Curie Research Training Grant under Contract No. ERB4001GT96383G.