Ambipolar diffusion in superfluid neutron stars

In this paper we reconsider the problem of magnetic field diffusion in neutron star cores. We model the star as consisting of a mixture of neutrons, protons and electrons, and allow for particle reactions and binary collisions between species. Our analysis is in much the same spirit as that of Goldreich&Reisenegger (1992), and we content ourselves with rough estimates of magnetic diffusion timescales, rather than solving accurately for some particular field geometry. However, our work improves upon previous treatments in one crucial respect: we allow for superfluidity in the neutron star matter. We find that the consequent mutual friction force, coupling the neutrons and charged particles, together with the suppression of particles collisions and reactions, drastically affect the ambipolar magnetic field diffusion timescale. In particular, the addition of superfluidity means that it is unlikely that there is ambipolar diffusion in magnetar cores on the timescale of the lifetimes of these objects, contradicting an assumption often made in the modelling of the flaring activity commonly observed in magnetars. Our work suggests that if a decaying magnetic field is indeed the cause of magnetar activity, the field evolution is likely to take place outside of the core, and might represent Hall/Ohmic diffusion in the stellar crust, or else that a mechanism other than standard ambipolar diffusion is active, e.g. flux expulsion due to the interaction between neutron vortices and magnetic fluxtubes.


INTRODUCTION
The structure and evolution of magnetic fields in compact stars is of considerable interest. In part this is due to the large amount of information available: magnetic field strengths have been measured, or more often inferred, for many neutron stars, allowing astronomers to build up a detailed picture of the distribution of fields strengths in them. Crucially, there seem to be pronounced systematic differences in field strength between different classes of neutron stars. The young neutron stars and most other 'normal' pulsars have fields of the order 10 12 G, while the older millisecond pulsars (MSPs) have much lower fields, of order 10 9 G. Similarly (relatively) low fields are believed to exist in the neutron stars located in accreting low-mass X-ray binary systems (LMXBs). In contrast, magnetars, as observed as anomalous X-ray pulsars (AXPs) and soft gamma repeaters (SGRs), seem to have exceptionally high field strengths, of order 10 15 G.
The way in which magnetic fields slowly evolve in neu-tron stars may well be crucial in understanding the relationship between these classes of objects. There are indications that millisecond pulsars acquire their relatively low field strength through a period of accelerated field decay while accreting in LMXBs, neatly linking the normal, LMXB and MSP populations.
Even more intriguingly, the structure, strength and, most of all, the evolution of the magnetic field may provide a unique insight into the state of the stellar interior. Field decay in the hot crust of an LMXB has been invoked to explain the reduction in field strength that accretion is believed to bring about, giving information on the conductivity of the crustal matter (see e.g. Cumming, Arras & Zweibel (2004)). More exotically, magnetic field evolution has been invoked to explain the restless behaviour of the magnetars. As was argued convincingly by Thompson & Duncan (1995, 1996, the existence of a 10 15 G strength field explains many of the observed features of the AXPs and the SGRs. In particular, Thompson and Duncan pointed out that the very slow rotation rates of the magnetars imply that the energy required to power their emission cannot be supplied by rotation, but a 10 15 G magnetic field, evolving on the inferred lifetime of the magnetars (of order 10 3 − 10 4 years), would provide a suitable energy reservoir.
In order for this model to be self-consistent, there must exist a mechanism capable of causing a 10 15 G magnetic field to evolve on a timescale of order 10 4 years. The issue of magnetic field evolution had been considered previously by Goldreich & Reisenegger (1992) (hereafter GR92), who showed that a variety of processes contribute to the field evolution in a neutron star. These processes are (i) Ohmic decay, in which the current and therefore the field decays through the frequent collisions between the different particle species; (ii) Hall drift, in which the magnetic field produced by the current acts back on the current itself, producing a (dissipation free) rearrangement of the magnetic flux; (iii) ambipolar drift, in which the charged particles (electrons and protons) and the magnetic field move relative to the neutral fluid (the neutrons). It was the last of these, ambipolar drift, that Thompson and Duncan were able to invoke to explain the activity of the magnetars. While not important for normal field strength stars, they noted that the quadratic scaling of the drift velocity with magnetic field strength led to field rearrangement timescales of order magnetar lifetimes for magnetar-like field strengths. In contrast, Ohmic decay is unlikely to be important in neutron star cores (Baym, Pethick & Pines 1969). The Hall drift, although not a dissipative process itself, might cause a rearrangement of the magnetic field on astrophysically relevant timescales, assuming a magnetar-strong magnetic field (GR92). However, as we discuss in this paper, this effect is likely to be suppressed in the presence of proton superconductivity. It should be noted, nevertheless, that the possibility remains that Hall and Ohmic processes, working in concert, might lead to significant field evolution in the stellar crust; see e.g. Pons, Miralles & Geppert (2009).
It is precisely this issue of magnetic field evolution in a neutron star core that we revisit in this paper. Given its apparent significance for the magnetar model, we pay particular attention to ambipolar diffusion. Our starting point is the classic work of Goldreich and Reisenegger (GR92), who modelled neutron stars as being composed of a neutronproton-electron plasma, with field evolution being driven by binary scattering processes between the species. The analysis of GR92 was essentially at the level of estimating the key timescales, rather than calculating explicit fully consistent solutions. Their analysis has recently been extended by Hoyos, Reisenegger & Valdivia (2008, in which the magnetic evolution problems were solved explicitly in a simplified geometry, to monitor the validity of the timescales estimated by GR92, and to understand the change in geometry brought about by the diffusive processes. Our aim in this paper is somewhat different. Rather than obtaining more explicit solutions to the original model, we instead will content ourselves with making simple estimates in the spirit of GR92, but with the addition of the crucial ingredient of neutron superfluidity. There is strong evidence, both from theory and observation, that mature neutron stars are likely to be sufficiently cold that the nuclear matter is in a superfluid state, see e.g. Sauls (1989). Our aim is to see what qualitative features result from this more realistic level of modelling, and, in particular, to see if timescales of order of the magnetar ages naturally emerge. That the timescales might change significantly is plausible, given that superfluidity is known to suppress particle collisions and reactions, but introduce a new mechanism for coupling the different fluids, namely mutual friction, in which the neutron and charged fluids are effectively coupled. This coupling might be brought about by the scattering of electrons off the vortices present in the rotating superfluid (Alpar, Langer & Sauls 1988). It is the changes in the field evolution timescales, and their implications for the magnetar model, that we seek to investigate.
The structure of this paper is as follows. In section 2 we set out the form of the equations to be solved, without specialising to any particular model for the fluids themselves. In section 3 we estimate diffusion timescales assuming that the neutron star matter is normal (i.e. not superfluid), essentially recovering part of the analysis of GR92. In section 4 we repeat the calculation assuming neutron and proton superfluidity. Finally, in section 5 we summarise our findings and discuss their implications for the magnetar model.

Ambipolar diffusion model and hydromagnetic equilibrium
Our hydrodynamical model for the neutron star core is a magnetised three-fluid system consisting of neutrons, protons and electrons (labelled by the index x = {n, p, e}), allowing for relative motion with velocities vx. Besides the proton-electron electromagnetic coupling, our model accounts for fluid coupling through binary particle scattering or vortex mutual friction (depending on whether some fluids are in a superfluid state or not). In addition, our model is 'transfusive' in the sense that it allows for particle conversion through reactions. Naturally, the model also includes several simplifications. In particular, in writing down our superfluid hydrodynamical equations of motion, we will neglect the effect of entrainment between neutrons and protons, which is present when both species are superfluid. This approximation is justified in view of the order-of-magnitude precision of our final results 1 .
We also need to clarify to what extent we account for the effects of the expected type II proton superconductivity in the magnetohydrodynamics of the outer neutron star core (for a recent discussion see Glampedakis, Andersson & Samuelsson (2010)). A key effect we do include in our analysis is the superconducting magnetic force appearing in the equations of motion in the place of the familiar Lorentz force, due to the presence of the quantised proton fluxtubes. An effect we do not include is the participation of the fluxtubes in the mechanics of mutual friction coupling between the fluids. This is an important omission in our model given that the dynamics of proton fluxtubes and their interaction with other stellar components (e.g. the neutron vortices) may play a central role in magnetic field evolution (see e.g. Ruderman, Zhu & Chen (1998)). The magnetic field properties (equilibrium, secular evolution) in a 'full' superconducting model will be discussed elsewhere.
The remaining simplifications adopted in our model have to do with the particular nature of ambipolar diffusion. By definition, this is a slow process (i.e. it acts on timescales much longer than any dynamical timescale) with the star evolving through a series of 'snapshots', each one of them representing a hydromagnetic quasi-equilibrium. In turn, and given that the magnetic field and particle reactions have a small overall impact on the bulk stellar structure, we can assume that at any given time the system represents a small deviation from a stationary, chemically equilibrated, non-magnetic star in uniform rotation.
At this point one has to be careful with the precise notion of hydromagnetic equilibrium. In the case of a single fluid star, this notion is unambiguous: it refers to a state where the fluid motion is that of rigid-body rotation and where there is a balance between the magnetic stresses, the fluid pressure and the centrifugal force (if rotation is significant). This force balance is the basis of all existing calcula- However, in more realistic multi-component stellar models, transfusion (particle reactions) can occur, and the nature of the hydromagnetic equilibrium may change due to the ability of particle reactions in driving fluid flow and viceversa (this will become apparent by inspection of the mass continuity equations (5) below). As will be discussed in detail, these flows themselves create frictional forces, from the drag of one species against another or (in the superfluid case) from the mechanism known as mutual friction. If the rate of particle scattering is high enough, the frictional forces connected with these diffusive flows might overcome the forces generated by pressure gradients and play a dominant role in balancing the magnetic stresses.
We are therefore led to make a distinction between two sorts of slowly evolving equilibria: friction-dominated equilibrium and transfusion-dominated equilibrium. Roughly speaking, the designation is as follows. At high temperatures, the rate of particle reactions is high, which tends to eliminate chemical potential imbalances and the corresponding pressure gradients, while the rate of particle scattering is also high, so the drag force induced by the fluid flow is strong. In this high temperature regime, the key force balance is between magnetic forces and the frictional drag forces; hence the designation friction-dominated equilibrium. Conversely, at low temperatures, the rate of particle reactions is low, so that chemical potential perturbations and the corresponding pressure gradients persist and are therefore important, while the rate of particle scattering is also low, so the drag force associated with the diffusive flow is weak. In this low temperature regime the key force balance is between the magnetic forces and the pressure gradients connected with the chemical potential imbalances; hence the designation transfusion-dominated equilibrium. In our discussion below we will encounter both types of quasi-equilibria.
Of course, no previous studies have calculated such multi-component equilibria. The closest things that exist are the multifluid but non-magnetic equilibria of Prix (2004), who computed stationary equilibria of rotating two-fluid stars. The two fluids were both rigidly rotating, and both particle reactions and particle scattering were 'switched off', so as to allow solutions that were both stationary and did not involve frictional forces, only pressure gradients, rotation and gravity. We will refer to such simplified equilibria (magnetised or not) as 'ideal' equilibria. They are of significance in as much as they are the sorts of equilibria that have appeared previously in the literature, and can be relatively easily calculated. Construction of an ideal magnetic equilibrium would be a useful first step to calculating the more physical equilibria. Note that the ideal equilibrium is relatively close to the transfusion-dominated equilibrium described above, in the sense that in both, it is pressure forces rather than frictional forces that are key in maintaining equilibrium.

Basic formalism
We can now assemble the equations that describe our system. According to the model defined above, when we formulate the relevant hydromagnetic equations we can (i) ignore time derivatives (effectively filtering out wave phenomena) and (ii) work at leading order with respect to the deviation from a non-magnetic background. These same approximations were also used by GR92.
First we have the Euler equations, on a per-volume basis, referred to a frame rotating with the star. We write down one equation to describe the neutrons, and one to describe the combined proton-electron charged fluid. We will take as our background solution a rigidly rotating non-magnetic star. Using standard cylindrical coordinates {̟, z, ϕ} with the angular velocity taking the form Ωẑ, where the hat indicates a unit vector, we have: where ρx are the mass densities (with nx = ρx/mx the corresponding number densities), µx are the chemical potentials and Φ is the gravitational potential. In the proton-electron equation we used charge neutrality np = ne and neglected the electron mass (hereafter we use mp = mn = m). The perturbed hydromagnetic system is described by the pair of Euler equations (adopting the so-called Cowling approximation which entails δΦ = 0), where F cpl and Fmag are, respectively, the coupling and magnetic forces per unit volume. Their particular form depends on the state of matter (normal or superfluid) and will be specified later. The Euler equations are supplemented by the mass con-tinuity equations ∇ · (nx vx) = Γx (5) which feature the particle reaction rates 2 Γx. Baryon and charge conservation dictate that Γn = −Γp and Γp = Γe. These properties are indeed satisfied by beta decay in neutron star cores, typically operating through the modified Urca reactions. The background system, as described by eqns. (1) and (2), is characterised by chemical equilibrium, i.e. µn = µp + µe. This is no longer true in the magnetically perturbed system. The departure from chemical equilibrium is quantified by the 'chemical imbalance' parameter: For a small β we can write Γn = mλβ (7) where the coefficient λ is a function of density and temperature.
Note that for the ideal equilibrium defined earlier, the Coriolis terms in the Euler equations would be zero (all fluids rotate rigidly), the coupling force F cpl would be set to zero, and particle reactions eliminated by setting λ = 0.
Returning to the Euler equations (3) and (4) we can combine them and produce an equivalent (but slightly more useful) set of equations. Taking the difference between the two equations and manipulating slightly gives where xp = ρp/ρ is the proton fraction (ρ = ρn + ρp is the total density) and We will refer to (8) as the 'difference Euler equation'. Taking instead the sum of the two equations gives where u = vn + xp w is the neutron-proton average velocity, weighted by density. We will refer to (10) as the 'total Euler equation'. By combining the continuity equations in a similar way we get the set, Strictly speaking, the two terms on the left hand side of equation (12) could have comparable magnitudes, unless we invoke a weak stratification (in the sense that xp and therefore also np/nn vary only slightly over length scales of interest) which would allow us to neglect the second term. However, and despite its common use by many authors (including GR92), this is not a well justified approximation in realistic neutron star models. Perhaps a better reason for neglecting the second term in (12) would be to use (11) which suggests a meridional neutron velocity component ϕ × vn ∼ xp(φ × w). Following this reasoning we will neglect the second term in (12), but it should be pointed out that none of our (order-of-magnitude precision) results and conclusions would change had we kept this term, as long as vn w. Hence, we write where we have defined Before proceeding to our discussion of ambipolar diffussion we can make some non-trivial observations based on the 'difference' and 'total' equations (8) and (10) in the limit of an ideal hydromagnetic equilibrium (i.e. vx = 0, λ = 0, F cpl = 0). Firstly, from eqn. (8), we see that a star endowed with a magnetic field cannot be in chemical equilibrium, as it is the chemical imbalance β that provides the 'pressure' gradient ∇β necessary to balance Fmag/ρp, the magnetic force per unit mass of charged plasma. At the same time, from eqn. (10) we see that only in the special case when xp is a function of β is Fmag/ρ, the magnetic force per unit total mass, equal to the gradient of a scalar. This is reminiscent of the case of a single-fluid barotropic star, where the notions of stratification and chemical equilibrium are not relevant, and Fmag/ρ is always balanced by the gradient of a scalar (e.g. Haskell et al. (2008)). Finally, from either Euler equation it follows that an axisymmetric system has a vanishing azimuthal magnetic component F ϕ mag = 0.

General strategy
Our strategy will now be to look for solutions to the above set of equations, where we regard the magnetic field B and therefore also Fmag as given, so that the magnetic force on the right hand sides of eqns. (8) and (10) play the role of source terms, driving departures from chemical equilibrium and ambipolar diffusion of the charged particles with respect to the neutrons. We will assume that the magnetic field itself is axisymmetric, and that the subsequent evolution preserves this axisymmetry. This analysis is very much in the spirit of GR92-we only seek to make rough estimates of timescales, so will not attempt to construct fully self-consistent solutions, neither for the background nor for the perturbations. The departure from chemical equilibrium is given by the chemical imbalance β. We will use the proton-neutron velocity difference w as our diagnostic of ambipolar diffusion, and refer to w as the ambipolar diffusion velocity. Of course, ambipolar diffusion refers to the combined motion of the charged species (protons plus electrons) relative to the neutrons, so for the use of the term 'ambipolar diffusion' to be strictly accurate, it is necessary that the electron-proton velocity difference wpe = vp − ve be smaller than the protonneutron velocity difference: We will check if this inequality is satisfied in our calculations below. Note that when it is not satisfied, we would still expect the charged particle flow to generate magnetic field evolution; it is simply that it is no longer appropriate to describe it as classical ambipolar diffusion.
In the calculations that follow, we will use the perturbation equations to make rough estimates of the ambipolar diffusion velocity w = | w|, which can then be immediately converted into timescales (τ ≡ L/w, where L is a characteristic length scale). Within the order-of-magnitude precision of this calculation the timescale τ is essentially equal to the evolution timescale ∼ L/vp derived by the magnetic induction equation (see, for example, GR92). We will make repeated use of eqn. (8), the difference Euler equation. As noted above, the magnetic force on the right hand of this equation is to be regarded as the driving term, sourcing ambipolar diffusion. The terms on the left hand side all depend linearly on w. For the Coriolis term this dependence is explicit. For the coupling term this dependence is a function of the assumed state of matter (normal verses superfluid), as will be described below. The term ∇β represents the difference in pressure forces on the neutron and charged fluids, and can be related to the ambipolar diffusion velocity via eqn. (13). In the next two sections we explore the ambipolar timescales for normal and superfluid neutrons, accounting for the relative importance of the different terms as a function of temperature.

AMBIPOLAR DIFFUSION: NORMAL MATTER
Before tackling the superfluid problem, we will look at the case of normal matter. We will assume that the particle reactions are those of modified Urca, with an efficiency factor (Sawyer 1989) where T8 = T /(10 8 K) and ρ14 = ρ/(10 14 g/cm 3 ). The magnetic force can be identified with the usual Lorentz force (per unit volume): The coupling forces arise from binary particle collisions, so that the force is a function of the particle collision frequencies Dxy: where w = vp − vn as defined previously, and wen = ve − vn. These collision frequencies have been calculated by Yakovlev & Shalybkov (1990), who found that the neutronproton scattering (mediated by the strong nuclear force) was much more frequent than the neutron-electron scattering (mediated by the weak nuclear force), so we will include only the former in our coupling force: where Dpn ≈ 6.6 × 10 16 T 2 8 ρ −1/3 14 The difference and total Euler equations then become where D = (ρ/ρn)Dpn. As we will now show, rotation has a minimal impact on ambipolar drift motion. Indeed, we can see that the ratio of the Coriolis force to the drag force is ∼ 2Ω/D ≪ 1. It turns out that the same ratio controls the magnitude of the azimuthal component of the ambipolar drift. To see this, note that from eqns. (21) and (22) our assumption of axisymmetry implies an azimuthal component F ϕ which implies a dominantly poloidal ambipolar flow (this conclusion assumes that v ̟ n is not exceedingly larger than w ̟ , which is a reasonable expectation). Note that if we were to consider a non-rotating star then the same argument would have led to vanishing azimuthal components F ϕ mag = w ϕ = 0. We now wish to use (21) to make an estimate of the ambipolar diffusion velocity. At this point GR92 carried out a decomposition on eqn. (21) itself, splitting the equation into its solenoidal (zero divergence) and irrotational (zero curl) parts. This decomposition has the attractive feature that the ∇β term is purely irrotational, so the solenoidal projection of the equation represents a balance purely between the drag and Lorentz forces. However, there is no simple way of decomposing the Lorentz force itself into its solenoidal and irrotational parts, so we will instead make use of a slightly different decomposition. We will decompose at the level of the ambipolar velocity w, splitting it into what we will loosely term its 'solenoidal' and 'non-solenoidal' parts where w s and w ns are defined to satisfy With these formal definitions, the non-solenoidal part of the velocity field is that part of w which is sensitive to particle transfusion, while the solenoidal part is not associated with perturbations in the chemical potentials/pressure. When inserted into (21), and dropping the Coriolis force, this gives To make progress we can exploit that fact that both the drag coefficient D and the reaction parameterλ ≈ mλ are temperature dependent, allowing us to compare the relative sizes of the drag and pressure gradient forces acting upon the non-solenoidal flow: where we used the normalisations L6 = L/(10 6 cm), x01 = xp/0.1. This ratio becomes unity at a temperature α = 1 ⇒ Tα ∼ 4 × 10 8 ρ 1/12 Given the steep temperature dependence of the ratio α, it makes sense to consider the low and high temperature cases separately. If we begin with the T > Tα regime, from eqn. (27) we see that we can neglect the β 'pressure' gradient term in (26) to give We see that at high temperatures, the drag forces dominate the pressure gradients and ambipolar diffusion represents a balance between the drag force and the Lorentz force.
We therefore see that for T > Tα we are in the frictiondominated regime of section 2.1, or, equivalently, the star is in a friction-dominated equilibrium. In this regime the decomposition (24) and (25) is no longer useful and we can estimate a common diffusion timescale yr. (30) where B15 = B/(10 15 G). For T < Tα, eqn. (27) shows we can instead neglect the effect of drag on the non-solenoidal flow, so that eqn. (26) becomes We now need to be careful in extracting a diffusion timescale. In fact we can calculate two distinct timescales, by following GR92 and identifying that part of the Lorentz force that represents a departure of Fmag/ρp from a perfect gradient: Recall that in the single fluid case, the force Fmag/ρ was purely irrotational. Unfortunately, as discussed in section 2.2, there does not yet exist a calculation of the magnetic field structure of a multifluid star, so we cannot be sure how Fmag compares in magnitude with the full Lorentz force. For a generic initial stellar field, it is conceivable that | Fmag| ∼ | Fmag|. Another possibility (once the system enters the T < Tα regime) is that the flow w s , and therefore Fmag, is quenched as a result of short-timescale, dynamical evolution; if viable, this process cannot be described by our model. Given our ignorance on this point, all we can say is that | Fmag| | Fmag| and we will parameterise accordingly in the equations that follow.
With this in mind, one timescale is associated with the non-solenoidal flow maintained by the non-zero β which is responsible for balancing the magnetic force. This is described by which leads to an ambipolar diffusion timescale yr.
The second timescale is associated with the solenoidal flow sourced by Fmag, Given that | Fmag| | Fmag| we obtain a lower bound on the associated diffusion timescale: yr.
The timescales given above are in agreement with the ones obtained by GR92, aside from our solenoidal timescale (36) being a lower bound rather than an estimate (GR92 implicitly assumed | Fmag| ∼ | Fmag|). This level of agreement is not surprising given that, so far, our analysis closely followed the one in GR92.
The following scenario then suggests itself: after the birth of a neutron star, once all the fast dynamical perturbations not considered here have died away, the hot (T > Tα) star cools, proceeding through a series of friction-dominated equilibria as described by eqn. (29). When the star cools below Tα, the evolution can then be dominated by the solenoidal flow of eqn. (35), provided that w s is not suppressed on a short dynamical timescale. Only on timescales longer than that of eqn. (36), after which the solenoidal flow is likely to die away, will the evolution proceed slowly through a series of transfusion-dominated equilibria as given by eqn. Finally, we can check if this motion represents ambipolar diffusion, in the sense of eqn. (15), i.e. we need to verify that the proton-electron relative speed wpe is much smaller than the ambipolar drift velocity w. The speed wpe can be easily estimated from the Ampère law to be wpe ∼ 8 × 10 −10 B15 L6x01ρ14 cm s −1 We therefore see that condition (15) which is comparable to the velocity estimate of eqn. (37). We therefore see that in this case the rapid motion of the protons with respect to the neutrons implies that this motion cannot be interpreted as ambipolar diffusion in the traditional sense. However, as described in section 2.3, we would nevertheless expect the magnetic field structure to evolve on the timescales estimated. The conceptual difficulty associated with the violation of the condition (15), although potentially important for ambipolar diffusion in normal matter, becomes irrelevant in the presence of proton superconductivity (as we discuss in section 4.2).

AMBIPOLAR DIFFUSION: SUPERFLUID MATTER
In the presence of superfluidity it is clear that the preceding ambipolar diffusion model needs revision. Once a given particle species is in a superfluid state, the rates of physical processes where these particles participate (e.g. binary collisions and chemical reactions) are exponentially suppressed. Hence, we need to make the following modifications to the normal matter model: where Rpn and R sf are superfluid suppression factors. These are exponentially small numbers for temperatures well below the critical temperatures Tcn and Tcp for neutron and proton superfluidity. In the T ≪ Tcn, Tcp regime, and as long as the critical temperatures are not too different, these factors can be approximated by the following fits to the available rigorous results (Haensel, Levenfish & Yakovlev 2000Baiko, Haensel & Yakovlev 2001), Rpn ≈ 0.01 y 2 n yp e −(yn+yp) , R sf ≈ 3 × 10 −4 y 5 p e −2yn (41) where yn = 1.188 Tcn/T and yp = 1.764 Tcp/T (the difference in the numerical factors reflects the different type of pairing expected in the outer core: 1 S0 pairing for the protons and 3 P2 pairing for the neutrons, see, for example, Sauls (1989)). Note that in general Rnp = R sf .
As a result of the combined superfluid suppression in particle reactions and collisions, the notion of a frictiondominated equilibrium (section 2.1) becomes irrelevant; this is demonstrated below in detail. Thus, a superfluid star is expected to be close to an ideal hydromagnetic equilibrium (i.e. in the transfusion-dominated regime) for any temperature below the critical ones.
Besides the suppression described by (40), superfluidity also adds a new coupling mechanism between the fluids through vortex mutual friction. Mutual friction refers to the force mediated by the neutron vortices (which are responsible for the star's bulk rotation) as a result of their interaction with the charged particles and the magnetic field. The standard expression for the resulting coupling force (per unit volume, acting on the neutrons) is given by (Hall & Vinen 1956;Andersson, Sidery & Comer 2006), The dimensionless coefficient R is associated with the effective drag force experienced by a moving vortex. For 'standard' mutual friction, i.e. scattering of electrons off the vortex's magnetic field, the drag value is estimated to be (Alpar, Langer & Sauls 1988;Andersson, Sidery & Comer 2006) In this case the coupling between the neutron fluid and the charged particles can be considered weak. The opposite limit of strong mutual friction corresponds to Strong drag could emerge as a result of interactions between vortices and fluxtubes (Sauls 1989;Ruderman, Zhu & Chen 1998;Link 2003). The extreme case of perfect vortex pinning corresponds to B ′ → 1 and B → 0. We will not consider the case of perfect pinning in this paper, as it would require a model which incorporates the detailed dynamics of neutron vortices and proton fluxtubes. These effects are beyond the scope of this paper and will be addressed elsewhere. The total coupling force is now the sum of the drag and mutual friction forces: where the drag force takes on the same functional form as in the normal fluid case, aside from the presence of the suppression factor Rpn: The drag force is strongly suppressed at low temperatures, so that in a sufficiently cold star the mutual friction force will dominate the coupling. However, it is possible that for higher temperatures (but still below the critical ones) the drag force could dominate the mutual friction, despite the suppression of particle collisions. We can quantify this argument by means of the ratio Depending on whether f ≪ 1 or f ≫ 1 we will refer to the 'cold superfluid' and 'hot superfluid' regimes, respectively. The transition between these two regimes occurs at a temperature T f at which f ∼ 1 and, due to the exponentially varying Rnp, it is very rapid.
We can estimate T f by considering the simple case where Tcn = Tcp = Tc and adopting the weak mutual friction limit (44). Using a typical value Tc = 5 × 10 9 K in (41) we find that the interval (say) 0.1 < f < 10 corresponds to T f ≈ 2.7 − 2.9 × 10 8 K, P = 10 s (49) T f ≈ 2.9 − 3.2 × 10 8 K, P = 0.1 s (50) illustrating the steep variation of f with temperature, and the weak dependence of T f on spin period. These estimates for T f can be slightly reduced if Tc is decreased or Tcn < Tcp (as predicted by some pairing models, see discussion in Andersson, Comer & Glampedakis (2005)). The transition from 'hot' to 'cold' superfluid hydrodynamics is therefore potentially relevant for relatively hot neutron stars. The existing state-of-art cooling calculations (Pons, Miralles & Geppert 2009) suggest that magnetar interior temperatures are expected to takes values of a few times 10 8 K, (at the observed spindown ages 10 3 − 10 4 yr). Thus, magnetars could be found in either superfluid regime and our modelling should be able to account for both possibilities. We will therefore look at each case in turn in the following two subsections.
A final important effect of superfluidity (and in particular of type II proton superconductivity) is related to the magnetic force Fmag. This force is not given by the Lorentz force (17) but, instead, is related to the smoothaveraged self-tension of the proton fluxtube array. The explicit form of the superconducting magnetic force is given in Glampedakis, Andersson & Samuelsson (2010). For the pur-poses of this paper it is sufficient to use the rough estimate where Hc1 ≈ 10 15 G is the critical magnetic field associated with the self-energy of a single fluxtube (Tilley & Tilley 1990), whose density dependence we will neglect.

The 'cold superfluid' regime
We first discuss the regime T < T f where 'cold' superfluid hydrodynamics applies, i.e. f ≪ 1, so that the coupling force can be approximated as being entirely due to mutual friction. In order to account for vortex mutual friction we obviously need to retain a non-zero rotation Ω in our equations. Also, we will continue to assume an axisymmetric system. To begin with, it is convenient to decompose the various velocities into poloidal and toroidal components using standard cylindrical coordinates. For example, Inserting these in the Euler eqns. (8) and (10), together with the force (42), we obtain the difference Euler equation and the total Euler: The toroidal and poloidal parts of the difference Euler equation are: Based on these general equations we can arrive at several non-trivial results. Firstly, the z component of (56) shows that the superfluid system does not allow the establishment of chemical equilibrium (β = 0) unless F z mag = 0. However, it is almost certain that a generic magnetic field configuration does not comply with this requirement. Secondly, we can show that the magnetic force is nearly poloidal, i.e. F ϕ mag ≪ F P mag , which is consistent with the proximity of the system (see eqn. (64) below) to an ideal axisymmetric hydromagnetic equilibrium (which has F ϕ mag = 0 exactly). Indeed, from (54) and (56) we have where we assumed that the bracket term is of order unity. Finally, by combining (55) and (57) we obtain, for both weak and strong mutual friction, In other words, superfluid ambipolar diffusion is dominantly toroidal. This is in contrast with the properties of ambipolar diffusion in a non-superfluid star, see eqn. (23).
In order to study ambipolar diffusion within the present superfluid model we will mainly rely on the components (55) and (56) of the difference Euler equation, as we did in the non-superfluid case. Moreover, we will continue to decompose the ambipolar velocity w into solenoidal and nonsolenoidal components, as in eqn. (24).
Using the property (59) we can approximate (56) aŝ where the effective 'drag' coefficient D sf takes the following form in the limits of weak and strong mutual friction: To make our strategy maximally transparent, we can eliminate β from eqn. (60) using the continuity equation making the decomposition into solenoidal and nonsolenoidal velocity terms explicit: This equation is the superfluid analogue of eqn. (26), and shows that the structure of the ambipolar diffusion problem in the superfluid case very closely mirrors that of the normal case: we have a driving term, build out of the Lorentz force, on the right hand side, while on the left we have terms linear in the ambipolar velocity.
We can now proceed much as in the normal case, described in section 3. The terms on the left hand side of eqn. (63) are strongly temperature dependent, so we will examine the ratio of the non-solenoidal mutual friction coupling force and pressure gradient term: where we used w ns ϕ ∼ −xpw ns ̟ /B for weak mutual friction (the strong mutual friction result differs only by a factor xp). Hence α sf ≪ 1 for any relevant temperature and so we can neglect the mutual friction term linear in w ns ϕ on the left hand side of eqn. (63) to givê As in the normal case, care must be taken in extracting diffusion timescales. For a flow with a significant nonsolenoidal component, the pressure gradient term will dominate, leading to We can now estimate the diffusion timescale for the nonsolenoidal flow, but now it is important to distinguish between the toroidal and poloidal parts: as a result of (59), we should expect that the shortest diffusion timescale is associated with the toroidal component w ns ϕ . Indeed, specialising briefly to the weak drag case, the ratio between the toroidal and poloidal ambipolar timescales is Then, combining eqns. (59) and (66) leads to (again assuming weak mutual friction) The resulting ambipolar timescale is  59), the corresponding timescale connected with poloidal flow is even longer, so we will not write it down.
Again in analogy with the normal case, we can consider the case of a solenoidal flow; for such a flow, eqn. (65) becomes where, as before, Fmag/ρp denotes that part of the Lorentz force Fmag/ρp that represents departure from a perfect gradient. Again assuming only that | Fmag| | Fmag| we can make a numerical estimate of the corresponding ambipolar timescale: The timescale connected with the poloidal ambipolar diffusion can be estimated as: where we assumed weak mutual friction. As we pointed out earlier, the strong drag result is simply a factor xp shorter. In either case, the toroidal timescale is independent of B.
The ambipolar diffusion timescales (69), (71), (72) are among the main results of this paper. Compared to the corresponding results of the normal matter model (section 3), these timescales are significantly different. This is not surprising given the different nature of fluid coupling (mutual friction) in the superfluid model. Indeed, these results can be understood in simple terms as follows. The nonsolenoidal timescale is much longer in the superfluid case, as non-solenoidal flows generate pressure and therefore compositional changes, which take place on a timescale set by the slow particle reactions. The solenoidal timescales are much shorter in the superfluid case, as the particle coupling forces which inhibit such motion are suppressed. This last point means that any solenoidal part of the Lorentz foce will be rapidly eliminated, leaving the star in a very slowly evolving transfusion-dominated equilibrium. This contrasts with the normal fluid case, where the (lower bound) on the solenoidal timescale (equation (36)) was of order 10 5 years, and was therefore potentially relevant to the currently observed magnetar population. The implications of the obtained timescales for magnetar evolution are discussed further in section 5.

The 'hot superfluidity' regime
We now consider the 'hot' regime T > T f where matter is superfluid but particle collisions still provide the main coupling mechanism through the force Thus, we would expect the analysis of section 3 to be directly applicable here, after imposing the suppression (40) and replacing the Lorentz force with the superconducting force (51). Recall, however, that in the case of normal matter, we were able to safely ignore the effects of rotation, as in the difference Euler eqn. (21) the Coriolis force was found to be much smaller than the drag force. However, in the present superfluid case the ratio of the Coriolis to drag force is increased by a factor Rpn: Given that the factor Rpn may be very small, it is therefore not immediately clear if rotation will be important in this 'hot superfluidity' regime. Inspection of (41) reveals that for most of the T > T f regime the above ratio remains small. Then, we can readily modify the results of section 3. To begin with, the ratio α and the temperature Tα are modified as From these expressions we can safely conclude that the 'hot' superfluid system will be always in the regime T ≪ Tα where both friction-dominated and transfusion-dominated quasiequilibria could be accessible. The normal matter ambipolar timescales (34)  where Rnp ∼ 10 −23 is a typical value for the collisional suppression factor for temperatures close to T f , as estimated using eqn. (41). As it was the case in the 'cold' superfluid regime, the evolution of the solenoidal and non-solenoidal flows in the 'hot' regime are associated, respectively, with a very short and very long timescale. Thus, the system is expected to be in a transfusion-dominated quasi-equilibrium. Finally we check if the condition w ≫ wpe is satisfied. It is straightforward to show that this is the case for both regimes of superfluidity. For type II proton superconductivity the magnetic field appearing in the Ampère law is the so-called London field bL (Cartler, Langlois & Prix 2000; Glampedakis, Andersson & Samuelsson 2010). The London field itself is given by the standard expression (see Tilley & Tilley (1990) for details), The intrinsic weakness of the London field implies a relative proton-electron velocity wpe many orders of magnitude smaller than the numerical estimate (37). Related to this last point, the enormous weakening of the electric current j = enp wpe in the presence of type II proton superconductivity also implies a negligible magnetic field evolution due to the Hall drift.

DISCUSSION: IMPLICATIONS FOR MAGNETARS
We now discuss the main results of this paper, consisting of the ambipolar diffusion timescales (69), (71), (72) for the 'cold' superfluid regime and (76), (77) for the 'hot' superfluid regime, in the context of magnetic field evolution in magnetars. Ambipolar diffusion could provide a viable mechanism for magnetic field evolution in the observed magnetar population provided the associated timescale is comparable to the estimated ages ∼ 10 3 − 10 4 yr for these objects. Let us first examine the non-solenoidal timescales of eqns. (69) and (76). As a consequence of the superfluid suppression of particle reactions (represented by the exponentially small factor R sf ) all non-solenoidal timescales become extremely long, much longer than the corresponding normal matter timescales (30), (34) and the estimated magnetar ages. An interesting exception could be the 'cold' superfluid timescale (69) if the mutual friction coupling could be reduced at the level of B ∼ 10 −4 R sf . However, such a very small value of B is unlikely to occur if mutual friction is weak, given that the rather well understood electron scattering by neutron vortices leads to a much stronger coupling B ∼ 10 −5 − 10 −4 (Alpar, Langer & Sauls 1988;Andersson, Sidery & Comer 2006). In principle, one could have a very small B in the opposite limit of nearly perfect vortex pinning. However, in the absence of input from microphysics, B would have to be arbitrarily chosen to be of the order ∼ 10 −19 in order for our model to produce ambipolar timescales of the order of 10 4 yr. Regardless of this issue of 'fine-tuning', and as we pointed out earlier, the regime of vortex pinning would require a more sophisticated superfluid model with detailed vortex/fluxtube dynamics (e.g. the superconducting model of Glampedakis, Andersson & Samuelsson (2010)).
A similar situation emerges when we examine the solenoidal timescales (71), (72) and (77). They clearly show that any solenoidal component of the Lorentz force will be eliminated on a very short dynamical-like timescale, much shorter than the estimated ages of the magnetar population. This means that the superfluid case is very much simpler than the normal fluid case, where the (lower bound) on the solenoidal timescale (equation (36)) was comparable to the estimated magnetar ages, so that, in the normal fluid regime, it was not clear if such flows were physically relevant.
Thus, we conclude that ambipolar diffusion is unlikely to be the driving mechanism for magnetic field evolution in the observed magnetars unless we invoke mutual friction very much weaker than current modelling suggests. This does not sit well with the standard magnetar model, where magnetic field diffusion is needed to account for the observed activity of the SGRs and AXPs. There do, however, exist mechanisms alternative to core ambipolar diffusion that might supply the required field evolution. One possibility is a combination of Hall and Ohmic processes in the stellar crust (Pons, Miralles & Geppert 2009). More speculatively, it is possible that vortex-fluxtube interactions, not included in our model, might provide the crucial evolution mechanism; such a process will be considered elsewhere. Yet another possibility is that the inner core is composed of exotic matter such as hyperons or quark condensates to which our model does not apply. Thus the poorly constrained physics of the inner core may allow for processes that lead to magnetic evolution on the required timescale.
It is possible to consider situations where only one of the proton and neutron fluids form condensates without resorting to further calculation. Let us first consider the case where the neutrons are superfluid but the protons normal. This would apply e.g. if the magnetic field were to exceed the 'second critical field' above which proton superconductivity is suppressed, i.e. if B > Hc2 ∼ 10 16 G (Tilley & Tilley 1990). In such a regime the absence of entrained protons around the neutron vortices would suppress the mutual friction coupling force (Feibelman 1971;Sauls, Stein & Serene 1982), so that the hot superfluid regime considered in section 4.2 would be relevant. The lack of proton superfluidity would also have the effect of decreasing the importance of the superfluid suppression on the particle collision and reaction rates, so the suppression factors in this case would lie somewhere between the values given in eqn.
(41) and unity (Haensel, Levenfish & Yakovlev 2000Baiko, Haensel & Yakovlev 2001). From eqns. (76) and (77) we see this would have the effect of pushing both the solenoidal and non-solenoidal diffusions timescales in the direction of physical interest. Now consider the case where the neutrons are normal but the protons superconducting. Such a model formed the basis of the analysis of Arras, Cumming & Thompson (2004), who studied the link between magnetic field decay and thermal emission. In terms of our analysis, the lack of neutron superfluidity eliminates the mutual friction force, so that we can simply apply the 'hot superfluidity' model of section 4.2, where the mutual friction force was not absent but was negligibly small. We can therefore make direct use of eqns. (76) and (77), with suitable values of the suppression factors (Baiko, Haensel & Yakovlev 2001;Haensel, Levenfish & Yakovlev 2001). We find that the qualitative picture is unchanged, with the non-solenoidal timescale remaining much too long to be of physical interest, while the (lower bound) on the solenoidal timescale remains extremely short.
Our model has a number of obvious weaknesses, the most serious of which is the neglect of interactions between the neutron vortices and the magnetic fluxtubes. The inclusion of such effects would provide a logical extension of the analysis of this paper. Our analysis could also be improved by considering particular realistic field configurations. This would allow explicit calculations of how the magnetic field evolves, giving more accurate timescales, and information on how the geometry of the field evolves in time. Indeed, such calculations have been carried out recently by Hoyos, Reisenegger & Valdivia (2008 for non-superfluid stars. Related to this last point, our work illustrates the need for the construction of accurate axisymmetric equilibria, to act as backgrounds about which to monitor the field evolution. However, as we have endeavored to point out, the construction of such solutions in complicated by the simultaneous existence of magnetic forces and particle reactions: the magnetic stresses imply the existence of a nonzero chemical imbalance, which in turn implies the existence of particle reactions, relative motion of the species, and particle scattering. Interestingly, in the case of normal matter, our analysis in section 3 makes clear that, if the Fmag/ρp component of the Lorentz force (the part which is not the gradient of a scalar) is significant, then the frictiondominated equilibrium that results has particle scattering as an important component of the force balance. Such forces have been absent from all magnetic equilibrium calculations to date, both perturbative (e.g. Haskell et al. (2008); Ciolfi et al. (2009);Lander & Jones (2009)) and numerical (e.g. Braithwaite & Nordlund (2006)), and so these normal fluid calculations might be significantly in error. Whether such a component of the Lorentz force is relevant, or whether it is removed on some fast dynamical timescale by processes not considered here, is certainly an issue affecting the normal fluid case. Clearly, a scheme for producing pseudo-stationary equilibria needs to be developed, for both normal matter and superfluid/superconducting stars, which would then provide a background for more detailed investigations of secular field evolution.