Deriving An Equation, Using An Integral From A Scientific Paper

Mon. Not. R. Astron. Soc. 385, 531–542 (2008) doi:10.1111/j.1365-2966.2008.12861.x Modelling magnetically deformed neutron stars B. Haskell,1� L. Samuelsson,1,2 K. Glampedakis1,3 and N. Andersson1 1School of Mathematics, University of Southampton, Southampton SO17 1BJ 2NORDITA, Roslagstullsbacken 23, 106 91 Stockholm, Sweden 3SISSA, via Beirut 2-4, 34014 Trieste, Italy Accepted 2007 December 14. Received 2007 December 14; in original form 2007 June 6 ABSTRACT Rotating deformed neutron stars are important potential sources for ground-based gravitational wave interferometers such as LIGO, GEO600 and VIRGO. One mechanism that may lead to significant non-asymmetries is the internal magnetic field. It is well known that a magnetic star will not be spherical and, if the magnetic axis is not aligned with the spin axis, the deformation will lead to the emission of gravitational waves. The aim of this paper is to develop a formalism that would allow us to model magnetically deformed stars, using both realistic equations of state and field configurations. As a first step, we consider a set of simplified model problems. Focusing on dipolar fields, we determine the internal magnetic field which is consistent with a given neutron star model and calculate the associated deformation. We discuss the relevance of our results for current gravitational wave detectors and future prospects. Key words: gravitational waves – magnetic fields – stars: neutron. 1 IN T RO D U C T I O N It is well known that a non-axisymmetric deformation of a rotating neutron star will lead to a time-varying quadrupole moment and could provide a good source of gravitational waves. In fact, rapidly rotating neutron stars are important potential sources of continu- ous gravitational waves for interferometric detectors such as LIGO, GEO600 and VIRGO as well as the planned next generation inter- ferometers which should be able to target them specifically. There have even been suggestions of narrow-banding advanced LIGO in order to target the low-mass X-ray binaries (see e.g. Brady & Creighton 2000), as it is thought that gravitational waves may be playing a role in setting the spin equilibrium period of these systems (Bildsten 1998; Andersson et al 2005). There are a number of mechanisms that may lead to a neutron star being deformed away from symmetry. First of all, the crust of a neutron star is elastic and can support ‘mountains’. The size of the deformation that can be sustained depends on many factors, such as the equation of state and the evolutionary history of the crust. This problem has been studied by Ushomirsky, Cutler & Bildsten (2000) and Haskell, Jones & Andersson (2006). Another mechanism for producing asymmetries is an oscillation mode developing an insta- bility, driven by gravitational radiation reaction, such as the r-mode instability proposed by Andersson (1998) (see Andersson 2003, for a relatively recent review). Deformations can also be caused by the magnetic field. Neutron stars are known to have significant magnetic fields and it is well known that a magnetic star will not remain spherical (Chandrasekhar & Fermi 1953; Katz 1989). If the magnetic axis is not aligned with the rotation axis the deformation �E-mail: haskellb@soton.ac.uk will not be axisymmetric and this will lead to gravitational wave emission. In fact, Cutler (2002) has suggested that a strong toroidal field could force a precessing neutron star to become unstable and ‘flip’ to become an orthogonal rotator.1 This would be an optimal configuration for gravitational wave emission. Heyl (2000) has sug- gested that magnetized white dwarfs may be interesting sources of gravitational waves. Several numerical studies have been directed at understanding the gravitational wave emission of magnetically distorted neutron stars (see e.g. Bonazzola & Gourgoulhon 1996; Ioka & Sasaki 2004; Tomimura & Eriguchi 2005). Finally, Melatos & Payne (2005) have studied the closely related problem of mag- netically confined accretion, which may lead to large deformations. The aim of this paper is to investigate deformations due to the in- terior magnetic field in more detail. We want to develop a framework that would allow us to quantify the relevance of these asymmetries for any given stellar model and magnetic field configuration. As a starting point it makes sense to study some simple model problems. This will give us a better idea of the nature of the problem and what the key issues are. Moreover, it is a natural way to proceed given that the interior field structure is uncertain. The best current models (see e.g. Braithwaite & Spruit 2006) tell us that the field will tend to have a mixed poloidal and toroidal nature. Any analysis should be able consider such generic configurations. It is, of course, the case that magnetically deformed stellar models have been studied previously and we can benefit greatly from the existing literature. It is particularly important to appreciate that the magnetic field configuration may be constrained for a given stellar model. In effect, the magnetic field must be solved for together with 1 The dynamics of a precessing magnetic star has been discussed in great detail by Wasserman (2003). C© 2008 The Authors. Journal compilation C© 2008 RAS D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 532 B. Haskell et al. the fluid configuration. This statement is quite obvious, but it has not always been appreciated in discussions of gravitational waves from magnetically deformed stars. We pay special attention to this issue by discussing the general problem in Section 2, and then working out the details for the stellar models we consider in Appendix A (uniform density) and Appendix B (polytropes). In many respects the discussion in the appendices is an adaptation of already existing results. This discussion provides the main input for the deformation calculation, which is presented in Section 3. The potential impact of our results on gravitational wave observations is discussed in the concluding section. Throughout the paper, we consider only dipolar fields, but the developed formalism is general and can easily be extended to any magnetic field. 2 MAG NETIC FIELDS IN STELLAR I N T E R I O R S Chandrasekhar & Fermi (1953) were the first to realize that a star would not remain spherical in the presence of a strong magnetic field. They calculated the deformation of an incompressible star with a constant dipolar field by minimizing the energy of the config- uration. The case of a constant-density star with an internal poloidal field matched to an external dipole was later considered by Ferraro (1954) and Goosens (1972), by solving the Euler equations. In our analysis we shall use the latter approach. However, before consider- ing magnetic deformations it is important to understand which are the permissible field configurations, given the equation of state. We shall see that the range of permitted fields is quite restricted for the simple stellar models that we want to consider. 2.1 Formulating the problem Let us consider the equilibrium configuration of a magnetic star. We will assume that the magnetic energy is small compared to the gravitational energy and that magnetic effects can be treated as a perturbation of a spherical, non-magnetic, background. This should always be the case for realistic neutron star parameters. For simplicity, we will focus on non-rotating stars. Then the equations of hydromagnetic equilibrium are ∇p ρ + ∇� = (∇ × B) × B 4πρ = L 4πρ , (1) where B is the magnetic field, p the pressure, ρ the density and L defines the Lorentz force. The gravitational potential � obeys Poisson equation ∇2� = 4πGρ. (2) The magnetic field must also obey, from Maxwell’s equations, ∇ · B = 0 (3) Now, let us suppose that we want to consider a model described by a barotropic equation of state, ρ = ρ(p). Then, following Roxburgh (1966) and taking the curl of equation (1), we obtain an equation for the magnetic field: ∇ × [ B × (∇ × B) ρ ] = 0. (4) Since this equation contains the density it should be solved si- multaneously with equation (1). In other words, the magnetic field structure is constrained by the density profile. As we discuss in Appendices A and B, this constraint can be quite severe. The main problem arises at the stellar surface. It is clear that the requirement that (4) holds when ρ → 0 may have significant impact on the magnetic field configuration. It is quite easy to argue that the constraint (4) may not be that relevant for more realistic models. This was understood by Mestel (1956) a long time ago. A possible resolution would be to consider a more complex equation of state, e.g. including the temperature dependence. For ρ = ρ(p, T) we should replace (4) by( ∂ρ ∂T ) p ∇T × ∇p − ρ2∇ × [ B × (∇ × B) 4πρ ] = 0. (5) In this equation the first term can (at least in principle) balance the second in the surface region, thus relaxing the constraint on the magnetic field. However, this would involve non-spherical vari- ations in the temperature leading to meridional circulation. This obviously complicates the problem considerably. In fact, as far as we understand, there has not (yet) been much progress on solving this problem. Another potential way to avoid imposing (4) would be to allow the neutron star crust to supply the stresses needed to allow the currents required for a smooth matching to the exterior. One would also need to analyse the transition from the fluid to the exterior vacuum in detail. However, this is not simply a matter of joining the interior fluid magnetohydrodynamics equations to a set of exterior vacuum Maxwell equations. In reality, one would want to account for the pulsar magnetosphere and the presence of electron- positron pairs etcetera. This is (again) a much harder problem than we would want consider at this point. Accepting that the constraint (4) needs to be considered for our models, we next assume that the magnetic field only produces small deviations from a spherically symmetric background model (essen- tially, we assume that the ratio of magnetic to gravitational potential energy is small). This allows us to expand all our variables in the form ψ(r, θ ) = ψ0(r) + ψ1(r)Pl, (6) where Pl are the standard Legendre polynomials and ψ1 is a small perturbation of O(B2). We will concentrate on quadrupole (l = 2) deformations, simply because they are optimal from the gravita- tional wave emission point of view, However, the formalism applies to the general problem, in which case the perturbation is given by a sum of Legendre polynomials. We can first of all solve the structure equations in the absence of a magnetic field and obtain the background model. The result is then fed into equation (4) to determine B, which we will need in order to solve equations (1) and (2) to first perturbative order for the quantities ρ, p and �. Restricting ourselves to axisymmetry, the φ component of the magnetic force must be zero, as there is nothing to balance it in equation (1). Hence [B × (∇ × B)]φ = 0. (7) With this understanding, let us examine some general magnetic field solutions for various background models. If one splits the magnetic field into two components, a poloidal one Bp = (Br , Bθ , 0) and a toroidal one Bt = (0, 0, Bφ), and introduces a stream function S(r, θ ) such that Br = 1 r2 sin θ ∂S ∂θ , (8) Bθ = − 1 r sin θ ∂S ∂r . (9) Then equations (3) and (7) reduce to (Roxburgh 1966) Bp · ∇(r sin θ Bt) = 0 (10) C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 Modelling magnetically deformed neutron stars 533 which gives Bφ = β(S) r sin θ , (11) where β is some function of the stream function S. This means that the toroidal part of the field is a function of the poloidal part. Equation (4) then takes the form ∂ ∂r { 1 ρr sin θ ∂S ∂θ [ 1 r sin θ ∂2S ∂r2 + 1 r3 ∂ ∂θ ( 1 sin θ ∂S ∂θ )] + β ρr2 sin2 θ ∂β ∂θ } − ∂ ∂θ { 1 ρr sin θ ∂S ∂r [ 1 r sin θ ∂2S ∂r2 + 1 r3 ∂ ∂θ ( 1 sin θ ∂S ∂θ )] + β ρr2 sin2 θ ∂β ∂r } = 0. (12) In the following we shall focus on dipole solutions. We can then take a stream function of form S(r, θ ) = A(r) sin2 θ (13) or directly look for solutions of the form B = r̂{W (r) cos θ} + θ̂{X(r) sin θ} + φ̂{iZ(r) sin θ}, (14) where the phase of the φ̂ component is chosen in such a way as to produce real valued equations in the following. For this kind of field, equation (3) gives rW ′(r) + 2 [W (r) + X(r)] = 0, (15) where the prime indicates differentiation with respect to r. We now have a formulation of the background problem. We need to solve equation (4) (or equivalently equation 12), which depends on the density profile ρ(r). It is thus necessary to prescribe an equa- tion of state for the stellar matter. We have chosen to consider two simple configurations. The case of a uniform density fluid is con- sidered in Appendix A, while we analyse the problem for n = 1 polytropes in Appendix B. Our main reason for focusing on these simple models is that they permit the main part of the calculation to be done analytically. This has the advantage of leading to relatively simple expressions, where the dependence on the various parame- ters in the problem are explicit, for the magnetic deformation. 2.2 Boundary conditions To complete the solution for the magnetic field configuration, we need to implement boundary conditions at the surface of the star (at the centre it is sufficient to impose regularity, cf. the analysis in Appendices A and B). At the surface, the solenoidal nature of the field requires continuity of the radial component 〈Br〉 = 0, (16) where 〈Br〉 = Br ext − Br int indicates the discontinuity between the external and internal parts of the field at the surface. Furthermore, we shall require continuity across the boundary of the traction vec- tor t i = T ij n̂j , i.e. the projection of the total stress tensor along the normal unit vector.2 This is equivalent to requiring local force balance per unit area. Our condition is then 〈t i〉 = 0. (17) 2 While the main part of our discussion uses standard spherical coordinates, the discussion of the boundary conditions becomes clearer if we use a coordinate basis. This means that the Einstein convention of summation over repeated indices applies. The stress tensor is the sum of a fluid piece and a magnetic piece T ij = −pδij + 1 4π ( BiBj − 1 2 B2δij ) . (18) Let us consider the projection along the normal to the surface. The normal vector to the perturbed surface will have the form n̂S = r̂ + n̂1S, (19) where n̂1S indicates the correction of O(B 2). Projecting Tij along this vector, we obtain −p0δij n̂1S − δp δir + 1 4π ( BiBr − 1 2 B2δir ) , (20) where p0 is the background pressure and δp the first-order pertur- bation. Note that the magnetic term is already O(B2) and one can take n̂S = r̂ for this term. At the surface the background pressure p0 vanishes, so we must impose continuity for t r = −δp + 1 4π [ (Br )2 − 1 2 B2 ] , (21) t θ = 1 4π BrBθ , (22) tφ = 1 4π BrBφ. (23) As we already have 〈Br〉 = 0, the θ and φ components of equa- tions (23) demand that 〈Bθ 〉 = 〈Bφ〉 = 0, which then leads to 〈B2〉 = 0. As p = 0 in the exterior, it must be the case that δp = 0 at the surface [all quantities are now O(B2) so we can consider the surface of the unperturbed configuration]. The conclusion is that all components of the magnetic field must be continuous across the interface, i.e. that we have no surface currents. We should stress that we are in principle allowed to introduce surface currents at this point. This would lead to discontinuities in the θ and φ components of the field. However, we feel that unless there are physical arguments dictating the nature of such currents there is no reason to introduce them. If we allow for the presence of an arbitrary current at the surface we have too much speculative freedom and it is not clear to what extent various models make sense. Hence, we find it more natural to restrict ourselves to models where the above traction conditions apply. Note that the traction conditions are automatically satisfied if we have a purely toroidal field. The external magnetic field, in fact, then solves ∇ · B = ∇ × B = 0 (24) and the assumption of axisymmetry forces a vanishing toroidal component, Bφ = 0. In the case of a purely poloidal field it is sufficient to match to an external curl-free dipole. In the case of a mixed poloidal/toroidal field, however, one must have that at the surface Z(R) = 0. (25) In the uniform density case, the analysis in Appendix A led to Z = arW, cf. (A19). Hence, we must have W = 0 (provided that a �= 0). This means that, under axisymmetry, a mixed poloidal and toroidal dipolar internal field and a general multipolar external field are forced to obey Brint = Brext = 0 and Bφint = Bφext = 0 at the surface. It is thus immediately clear that we cannot match this kind of internal field to an external dipole field. C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 534 B. Haskell et al. One also has a condition on Bθ , as the exterior solution for a field which is regular at infinity is of the form Bext = ∑ l≥m { −(l + 1) Al rl+2 Y ml r̂ + Al rl+2 ∂θY m l θ̂ + imAl rl+2 sin θY ml φ̂ } . (26) It is clear that if the r̂ component must vanish at the surface, so must the θ̂ component. We thus see that, for a constant-density star, all components of the mixed magnetic field must vanish at the surface. The discussion in Appendix B shows that this is true also in the case of a polytropic equation of state. In conclusion one can therefore, for the barotropic equations of state that we consider, match an interior mixed poloidal/toroidal dipole field only to a vanishing external field. It is also clear that purely toroidal fields are not allowed, since (A17) is only compatible with the surface condition if c = 0. It is interesting to note that if we consider the simplest possible fluid configuration (uniform density) then the boundary conditions restrict us to a very limited class of magnetic fields: notably poloidal or mixed poloidal/toroidal fields that vanish outside the star. Moreover, the analysis of the n = 1 polytrope case in Appendix B shows that these limitations are not due to us choosing a pathological model for the equation of state. 3 MAG N E T I C D E F O R M AT I O N S So far we have discussed how the requirement that the magnetic field be consistent with the fluid configuration restricts the model parameter space. In particular, for polytropes it does not allow us to freely specify the toroidal component of the field in the mixed case. Also, for most of the models that we have considered, one can only consider fields which vanish in the exterior. This is clearly not physical, as one would expect the exterior field to be prevalently dipolar far from the star. This problem could easily be fixed by allowing surface currents (which one may expect as there are strong electric fields at the stellar surface). Nevertheless, given that we do not have a physical model for such currents, we will not introduce them here. We now turn to the main problem, how a given magnetic field deforms the shape of the star. As we want to be able to compare the results for different field configurations, it is worth thinking about how such a comparison would be carried out. After all, each model field will be naturally described by some set of parameters which may not be easily translatable from model to model (see Appen- dices A and B). Intuitively, one would expect the deformation of the star to depend on the ratio of magnetic to gravitational potential energy (see e.g. Cutler 2002). Hence, it is natural to use the magnetic energy as a measure. Given this, we will express our results in terms of the parameter B̄2 = 1 V ∫ B2 dV , (27) where V is the volume of the star. Having determined a set of magnetic field configurations which are consistent with a chosen stellar interior,3 we turn our attention 3 It should be noted that we are only requiring the field to be consistent here. In reality, one can constrain the model further by requiring that the field is also stable. For discussions of the stability problem, and references to the relevant literature, see Glampedakis & Andersson (2007) and Akgün & Wasserman (2008). to solving equation (1) in order to obtain the new equilibrium shape of the star. Inspired by the work of Saio (1982) and our own recent work on crustal deformations (Haskell et al. 2006), we shall define a new variable x(r, θ ) = r [1 + ε(r)Pl(cos θ )] , (28) where r is the standard radial variable and Pl is a Legendre poly- nomial representing the deformation. Note that we are considering only one polynomial here, but the formalism could easily be ex- tended to a sum of Legendre polynomials. The perturbed surface thus takes the form xs(R, θ) = R [1 + ε(R)Pl(cos θ )] , (29) where R is the radius of the unperturbed (spherical) star. We shall also assume that the pressure, density and gravitational potential take the form ψ(r, θ ) = ψ(r) + δψl(r)Pl, (30) where ψ is the background quantity and δψl is a small perturbation of O(B2). From now on we shall, unless there is a risk of confusion, write δψ instead of δψl . Equation (1) then leads to( dδp dr + ρ dδ� dr + δρ d� dr ) Pl = [(∇ × B) × B]r 4π = Lr 4π , (31) (δp + ρδ�) dPl dθ = r [(∇ × B) × B]θ 4π = rLθ 4π , (32) which must be solved together with the perturbed Poisson equa- tion (for l = 2), d2δ� dr2 + 2 r dδ� dr − 6 r2 δ� = 4πδρ. (33) Let us first of all tackle the case of an incompressible star. 3.1 Deformations of incompressible stars In the case of a uniform density star we consider the field in equa- tion (A14). This gives us for the Lorentz force Lr = 2 3 ( B R2 )2 ( R2r 3 − 2r 3 5 ) (1 − P2) , (34) Lθ = −1 3 ( B R2 )2 ( R2r 3 − r 3 5 ) dP2 dθ , (35) as well as, obviously, Lφ = 0. In the case of an incompressible star, there can be no l = 0 deformation, as this would not conserve the volume and therefore not conserve the mass of the star. We shall thus only consider the case of l = 2. For this case δ� inside the star takes the form δ� = −4 5 πGρε(R)r2. (36) Inserting this solution into (32) and evaluating at the surface gives δp(R) = 4 5 πGρ2ε(R)R2 − B 2 90π . (37) If we now remember that δp(r) = δp(x) − ε(r)r dp dr (r), (38) we have at the surface δp(R) = 0 − ε(R)R dp dr (R). (39) C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 Modelling magnetically deformed neutron stars 535 By using the background pressure p = 2πGρ2 1 3 (R2 − r2), (40) equation (37) gives us ε(R) = − 1 48 B2 π2Gρ2R2 (41) which agrees with the result of Goosens (1972). In order to assess the relevance of this deformation for gravita- tional wave emission we calculate the ellipticity, which is defined as � = Izz − Ixx I0 . (42) Here I0 is the moment of inertia of the spherical star, while Ijk is the inertia tensor: Ijk = ∫ V ρ(r) ( r2δjk − xjxk ) dV . (43) For the field in (A14) we have B̄ = 0.1B, and hence we find that � = −3 2 ε = 1 18 B2R4 GM2 ≈ 10−12 ( R 10 km )4 ( M 1.4 M� )−2 ( B̄ 1012G )2 . (44) The ellipticity is positive, so the star is oblate, as expected. As a sanity check of this result we can compare to the estimate used by Cutler (2002). For a field of 1015 G he assumes that � ≈ 1.6 × 10−6. Our more detailed calculation has led to a result that is almost a factor of 2 smaller. 3.2 n = 1 polytrope with a poloidal field Let us now consider the case of a star with an n = 1 polytropic equation of state. The magnetic field appropriate for this case is that in equation (B4), for which the Lorentz force has non-trivial l = 2 components4 Lr = − π 3B2s R 2(π2 − 6)y [2y 3 + 3(y2 − 2) (y cos y − sin y)] sin yP2(θ ), (45) Lθ = − π 3B2s R 2(π2 − 6)y [y 3 + 3(y2 − 2) sin y + 6y cos y] sin y dP2 dθ , (46) where y = πr/R and Bs represents the conventional dipole field strength at the surface. Having worked out the Lorentz force we can solve the perturbed hydrostatic equilibrium equations (31)–(32) and the Poisson equa- tion (33), with the condition that δ� and δ�′ be regular at the centre and match the exterior solution at the surface. From equation (32) 4 In the general case there will be both a radial and a quadrupole deformation. Since we are primarily interested in the gravitational wave aspects we ignore the l = 0 deformation in the discussion. This contribution would be important if one wanted to (say) study the oscillations of deformed magnetic stars. The required calculation is provided in Appendix C. we then obtain δρ = − π 5B2s 8(π2 − 6)2GR2ρcy × [ y3 + 3(y2 − 2) sin y + 6y cos y + 2(π 2 − 6)2ρc π4B2s yδ� ] , (47) where ρc is the central density. Inserting this into the right-hand side of equation (33) we can solve for δ�. We then get δ� = π 4B2s 4(π2 − 6)2ρcy3 × [−2y5 − 3(y4 + 4π2y2 − 12π2) sin y + y(y4 − 36π2) cos y]. (48) Finally, we can evaluate the distortion at the surface: ε = − δρ(R) R [ dρ dr ]−1 r=R = π 5 ( π2 − 24) Bs2 16R2ρc2G ( π2 − 6)2 (49) and obtain the ellipticity � = −3π 5 ( π2 − 12) R4Bs2( π2 − 6)3 GM2 . (50) Scaling to canonical neutron star values for the different parameters, and introducing the volume averaged field as in the constant-density case (here B̄ = 0.54Bs), we have � ≈ 2 × 10−10 ( R 10 km )4 ( M 1.4 M� )−2 ( B̄ 1012 G )2 . (51) In other words, for this model configuration the magnetic field induced ellipticity is two orders of magnitude larger than in our constant-density model. This is a useful indication of how much the result can vary for supposedly ‘similar’ field configurations if we change the equation of state. Of course, in making this statement one must keep in mind that the two field configurations really are different. 3.3 n = 1 polytrope with a toroidal field We next consider the case of a purely toroidal magnetic field in a star with an n = 1 polytropic equation of state. For the field given in (B13) the Lorentz force has l = 2 components, Lr = 2 3 B2 πRy [sin y + y cos y] sin yP2(θ ), (52) Lθ = 2 3 B2 πRy sin2 y dP2 dθ . (53) From equation (32) we can then write δρ = − 1 24R2ρcGπ2 [6π3ρc δ� − B2y sin y] (54) which, inserted on the right-hand side of Poisson equation (33), gives δ� = − B 2 36π3ρcy3 [(y5 − 15π2y) cos y + 5π(π − y2) sin y]. (55) We can then calculate the distortion at the surface ε = − 1 144 B2 ( π2 − 15) π2R2Gρc2 . (56) It follows that the ellipticity is given by � = 1 9π2 (π2 − 15) (π2 − 6) R4 B2 G M2 . (57) C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 536 B. Haskell et al. In terms of the average field which in this case corresponds to B̄ ≈ −0.17B, we have � ≈ −10−12 ( R 10 km )4 ( M 1.4 M� )−2 ( B̄ 1012 G )2 . (58) The deformation is now prolate, and notably about two orders of magnitude smaller than in the case of a poloidal field (for the same B̄). We can compare this result to the estimate used by Cutler (2002). For a 1015 G toroidal field he assumes that � ≈ −1.6 × 10−6 (the same magnitude as in the poloidal case). The deformation we have determined is slightly smaller than this. 3.4 General deformations Having considered some particular examples, let us now present the formalism for a more general field configuration and equa- tion of state. Assume that the magnetic field takes the form of equation (B6), for which the Lorentz force is Lr = −dA dr sin2 θ r4 [ (π2λ2r2 − 2)A + d 2A dr2 ] , (59) Lθ = −2A r5 cos θ sin θ [ (π2λ2r2 − 2) + d 2A dr2 ] , (60) where λ is a parameter that determines the ratio of the toroidal and poloidal field components (see Appendix B for details). Focusing our attention on the l = 2 components, we have (see Appendix C for a discussion of the associated radial deformation) Lr (l = 2) = dA dr 2P2 3r4 [ (π2λ2r2 − 2)A + d 2A dr2 ] , (61) Lθ (l = 2) = 2A 3r5 dP2 dθ [ (π2λ2r2 − 2)A + d 2A dr2 ] . (62) From the angular part of the perturbed hydrostatic equilibrium equa- tions (1) we obtain (for the l = 2 case that we are considering) δp = −ρδ� + rLθ = −ρδ� + 2 3 A r4 [ (π2λ2r2 − 2)A + d 2A dr2 ] (63) which substituted back into the radial part of the equations leads to δρ = − 1 3r5 ( d� dr )−1 [ −3 dρ dr δ�r5 + 2r dA dr A (π2λ2r2 − 2) −4A2 (π2λ2r2 − 4) + 2r3Ad 3A dr3 − 4r2Ad 2A dr2 ] . (64) This now allows us to compute the source term of the perturbed Poisson equation (33), which, for l = 2, reads d2δ� dr2 + 2 r dδ� dr − δ� [ 6 r2 + 4πG r5 ( d� dr )−1(dρ dr r5 )] =−8πG 3r5 ( d� dr )−1[ r dA dr A(π2λ2r2 − 2) − 2A2(π2λ2r2 − 4) + r3Ad 3A dr3 − 2r2Ad 2A dr2 ] . (65) The boundary conditions for this equation remain regularity at the centre and continuity of δ� and its derivative at the surface. For l = 2 the exterior gravitational potential has the form δ�ext ∝ 1 r3 . (66) Table 1. The ellipticity � for a sample of eigenvalues λ, representing dif- ferent mixtures of poloidal and toroidal fields. The results correspond to the parameters M = 1.4 M�, R = 10 km and we take the energy averaged amplitude of the poloidal part of the field to be 1012 G. The energy averaged toroidal field (B̄t) is given in each case. As can be seen, the star starts off with a small deformation, as the effects of the poloidal and toroidal compo- nents cancel each other, and becomes more and more prolate as the toroidal component of the field is increased. Note that we have a vanishing exterior field in all cases. λ � B̄t 2.362 −6.3 × 10−13 3.0 × 1012 3.408 −2.8 × 10−12 4.2 × 1012 7.459 −4.4 × 10−11 8.1 × 1012 25.488 −2.3 × 10−10 2.7 × 1013 33.491 −4.4 × 10−10 3.5 × 1013 58.495 −3.7 × 10−8 9.3 × 1013 318.499 −1.3 × 10−6 4.8 × 1014 Having computed δ� in the interior we can obtain δρ from equa- tion (64) and thus compute the ellipticity � and the deformation at the surface ε(R) = − δρ(R) R [ dρ dr ]−1 r=R . (67) As an example we carry out the calculation for a polytrope and a mixed poloidal/toroidal field. Then we can use the stream function from equation (B10) to determine the magnetic field and the Lorentz force. The result will obviously depend on which eigenvalue we choose for λ (recall that the field parameters are not continuous in this case). This way we can vary the relative strength of the toroidal part of the field. A sample of the obtained results is listed in Table 1. As expected, for low values of λ the effects of the poloidal and toroidal field work against each other and the deformation is comparatively small. It is possibly surprising that the cancellation is so efficient. If one were to, for example, add the deformations deduced for a pure poloidal field (51) to that for a purely toroidal field (58) using the appropriate values for B̄ then one would predict that the deformation should be � ∼ 10−10 for the first eigenvalue in the table. Furthermore, one would have expected the deformation to be oblate. We see that for the mixed field we are using the true deformation is two orders of magnitude smaller. It is also prolate. This illustrates that it is important to use both realistic field configurations and equations of state when estimating magnetic neutron star deformations. The results in Table 1 show that the star is always prolate, and as λ grows the toroidal field becomes dominant and the star becomes more and more prolate. 3.5 n = 1 polytrope, field confined to the crust Finally, let us consider the case when the field is confined to a region close to the surface. As discussed in Appendix B, we will call this the ‘crust’, even though we are neglecting its elastic properties. We consider the core to be unperturbed, and impose the condition (B14) at the inner boundary. Then the deformation is confined to the crust and, as the density is low in this region, the quadrupole moment, and thus the ellipticity, will in general be much smaller than if we were considering deformations of the whole star. However, it is possible that the field is confined to the crust after being expelled from a superconducting core. In reality, the flux expulsion takes such a long time that it will not be ‘completed’ in observed neutron C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 Modelling magnetically deformed neutron stars 537 Table 2. The ellipticity � for various eigenvalues λ when the field is confined to a region close to the surface (e.g. the crust). The results correspond to the parameters M = 1.4 M�, R = 10 km. We take the field to have the same energy as that of a field extended to the whole star with λ = 318.499. The energy-averaged values of the poloidal (B̄p) and toroidal (B̄t) magnetic field are also given. We should thus compare all results for the ellipticity in this table with that for λ = 318.499 in Table 1. The field is confined to a region close to the surface beginning at a radius rb = 9 × 105 cm, roughly corresponding to a 1 km thick crust. The core is taken to be spherical. As can be seen from the listed results, the deformations are larger than in the case of the field extending to the whole star, by up to an order of magnitude. λ � B̄p B̄t 32.98 −1.0 × 10−7 6.2 × 1014 1.2 × 1015 158.01 −5.5 × 10−7 5.8 × 1014 1.3 × 1015 321.00 −1.7 × 10−5 5.6 × 1014 1.4 × 1015 stars. Yet it is interesting to consider this problem as an extreme example. If the flux is completely expelled a large number of field lines are squeezed into a small region close to the surface of the star. This leads to a strong field and large deformations. This result can be seen in Table 2, where we have calculated the deformation of the star by assuming that the total magnetic energy in the crust is equal to that of the field extended to the whole star, as calculated in the previous section. The ellipticity is then larger than in the case of a field extending throughout the star, by up to more than an order of magnitude. This agrees with what was found by Bonazzola & Gourgoulhon (1996), who analysed the deformations due to a poloidal field confined to a thin shell close to the surface of a neutron star. In this context it is worth commenting on the very recent results of Akgün & Wasserman (2008), who consider a neutron star with a region of Type II superconductivity. Since the magnetic flux is then carried by flux tubes, this case is significantly different from our current problem. The main difference between the results is that in the superconducting problem one finds that the deformation scales with BHc1, where Hc1 ≈ 1015 G is the lower critical field, rather than B2. This immediately suggests that for B = 1012 G one would expect the deformation to be a factor of about 1000 larger in the su- perconducting case. The expectation is confirmed if one compares the results of Akgün & Wasserman (2008) to equation (58). The possibility that the core may sustain significantly larger deforma- tions than we predict in this paper provides strong motivation for more detailed analysis of the various issues concerning neutron star superconductivity. 4 D ISCUSSION We have presented a scheme for calculating the magnetic deforma- tions of a neutron star. We have considered the particular case of a dipolar field in a non-rotating star, but the formalism is readily ex- tended to other cases. The formalism can, for example, be adapted to the case of slow rotation if we take the magnetic axis to be aligned with the spin axis of the star (see Appendix D). We have considered the case of a purely poloidal field, the effect of which is to make the star oblate, that of a purely toroidal field which makes the star prolate, and the case of a mixed poloidal and toroidal field, which makes the star more and more prolate as the strength of the toroidal component is increased. We have seen that the condition that the magnetic field be consistent with the background equation of state can be quite restrictive. For most of the present set of model prob- lems, we are forced to consider only fields that vanish outside the star. This is clearly not physical as one would expect the field to be prevalently dipolar far from the star. This problem could be solved by introducing surface currents. Since we do not have a physical model for such currents we have not considered this possibility in detail. Nevertheless, our formalism is sufficiently general that it could be extended to any field configuration. It seems appropriate to conclude this investigation by discussing the relevance of our results for gravitational wave observations. A rotating magnetic neutron star becomes interesting from the gravi- tational wave point of view when the magnetic axis is not aligned with the spin axis. Then the deformation is no longer axisymmetric and we have a time-varying quadrupole moment. This is in contrast with the typically much larger (see Appendix D) rotational deforma- tion which is always axisymmetric and cannot radiate gravitational waves. The framework we have described can, in principle, be used for any neutron star equation of state or field configuration. It al- lows us to calculate the ellipticity, and thus the quadrupole, in order to make quantitative gravitational wave estimates. This is a useful exercise since we can compare our results to observational upper limits on the neutron star ellipticity, and perhaps constrain the pa- rameters of proposed theoretical models. To see how this works out in the present case, let us consider the recent upper limits on pulsar signals from the LIGO effort (Abbott et al. 2007). The strongest cur- rent constraint comes from PSR J2124−3358, a millisecond pulsar spinning at 202.8 Hz. For this system, LIGO obtains the constraint � < 7.1 × 10−7. This is impressive, because it restricts the maximal neutron star mountain to a fraction of a centimetre. The question is, how much does this constrain the magnetic field? From (51) and (58) we easily find that we must have B̄ < { 6 × 1013 G for a purely poloidal field, 8 × 1014 G for a purely toroidal field (68) (taking radius and mass to have the canonical values). Is this con- straint telling us anything interesting? The answer is probably no. The pulsars currently considered by LIGO are almost exclusively fast spinning, taken from the millisecond spin period sample. For these neutron stars the exterior dipole field inferred from the spin- down is weak. In the case of PSR J2124−3358 the standard estimate would suggest that the exterior field is about 3 × 108 G. It is, how- ever, thought that the interior field could be significantly stronger than that in the exterior. It could well be that our simple models provide useful representations of this interior field, and what is missing is a relatively small part that is irrelevant when it comes to deforming the star. Nevertheless, it does not seem reasonable to suggest that the interior field can be as much as six orders of magnitude stronger than the exterior one (for reasons of stability etcetera). Hence, our results indicate that LIGO has not yet reached the sensitivity, where it constrains any reasonable theory. Given this, one may ask to what extent future observations will improve on the current results. It is relatively straightforward to estimate how well one can hope do to if we recall that (using matched filtering) the effective gravitational wave amplitude increases as the square root of the number of observed wave cycles. For a spinning neutron star, where the spin frequency is essentially fixed during the observation, this means that the signal-to-noise ratio scales with the observation time as t1/2obs . The data used by Abbott et al. (2007) is based on ob- servations lasting the order of three weeks. If one could extend this to a full year, one would gain a factor of 4 or so in amplitude. This translates into an improvement of a factor of 2 in the constraint on the magnetic field, not enough to test a realistic theory model. In order to do much better one would have to improve the detector C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 538 B. Haskell et al. sensitivity. With the advanced LIGO upgrade, the sensitivity will increase by about one order of magnitude. With a one year obser- vation we then get a further √ 10 improvement in the constraint on the magnetic field. Advanced detectors may also have the ability to narrow-band and focus on a particular frequency. Suppose one were to do that and get (say) an additional order of magnitude improve- ment of sensitivity at the particular frequency for a given pulsar. This would still only improve the constraints in (68) by a factor of about 2×√10×√10 = 20. These rough estimates suggest that we should probably not expect to detect gravitational waves from mag- netically deformed millisecond pulsars. There are obvious caveats to this, in particular, associated with the expected superconductivity of the core. If all the magnetic flux is confined to the crust (which is unlikely but theoretically possible) then the deformation may be larger, cf. the results in Table 2. A region of Type II superconduc- tivity may also allow larger deformations, see Akgün & Wasserman (2008). Even though a rough estimate suggests that this will not affect the conclusions of the above discussion much, the problem is clearly worth further investigating. The situation is not quite as pessimistic if one considers slower spinning pulsars. Consider, for example, the current LIGO upper limit of � < 2 × 10−3 for the Crab pulsar (Abbott et al. 2007). For the models we have considered we must then have B̄ < { 3 × 1015 G for a purely poloidal field, 4 × 1016 G for a purely toroidal field. (69) At first sight, these limits may seem less interesting than those from the faster spinning system. However, we need to remember that young pulsars, like the Crab, are thought to have much stronger magnetic fields. From the spin-down one would estimate that the exterior dipole field is 4 × 1012 G for the Crab. Hence, the current observed upper limit is only a factor of about 1000 away from the expected field. One should be able to improve this by at least an order of magnitude with future detectors (the current LIGO detectors may be improved at lower frequencies, advanced LIGO will do better and VIRGO is expected to have good performance in the low-frequency regime) and a full year of observation. This suggests that gravitational wave observations may in the future be able to test theoretical models, where the interior magnetic field is more than two orders of magnitude stronger than the exterior field. The challenge now is to improve on the rough estimates we have presented in this paper. Future work needs to consider the defor- mation for realistic field configurations such as those worked out by Braithwaite & Spruit (2006). One must also understand the role of superconductivity better. There are also a number of interesting astrophysics questions, in particular concerning accreting systems. These involve the burial of the field by accreted matter, leading to a relatively weak exterior field, and potential accretion induced asymmetries (Melatos & Payne 2005). Finally, it may also be inter- esting to consider magnetars. After all, they are expected to have fields as strong as the constraints given in (68). Of course, they are also spinning very slowly which places them outside the bandwidth of any ground-based gravitational wave detector. Nevertheless, a newly born magnetar may occasionally pass through the detection window. Such an event should be observable from within the galaxy, and perhaps beyond. AC K N OW L E D G M E N T S This work was supported by PPARC through grant numbers PPA/G/S/2002/00038 and PP/E001025/1. BH would like to thank the Institut Henri Poincare–Centre Emile Borel for hospitality and support during part of this work. NA acknowledges support from PPARC via Senior Research Fellowship no PP/C505791/1. LS gratefully acknowledges support from a Marie Curie Intra-European Fellowship, contract number MEIF-CT-2005-009366. REFERENCES Abbott B. et al., 2007, Phys. Rev. D, 76, 042001 Akgün T., Wasserman I., 2008, MNRAS, 383, 1551 Andersson N., 1998, ApJ, 502, 708 Andersson N., 2003, Class. Quantum Gravity, 20, R105 Andersson N., Glampedakis K., Haskell B., Watts A. L., 2005, MNRAS, 361, 1153 Bildsten L., 1998, ApJ, 501, L89 Bonazzola S., Gourgoulhon E., 1996, A&A, 312, 675 Brady P. R., Creighton T., 2000, Phys. Rev. D, 61, 082001 Braithwaite J., Spruit H. C., 2006, A&A, 450, 1097 Chandrasekhar S., Fermi E., 1953, ApJ, 118, 116 Cutler C., 2002, Phys. Rev. D, 66, 084025 Ferraro V. C. A., 1954, ApJ, 119, 407 Glampedakis K., Andersson N., 2007, MNRAS, 77, 630 Goosens M., 1972, Ap&SS, 16, 286 Haskell B., Jones D. I., Andersson N., 2006, MNRAS, 373, 1423 Heyl J. S., 2000, MNRAS, 317, 310 Ioka K., 2001, MNRAS, 327, 639 Ioka K., Sasaki M., 2004, ApJ, 600, 296 Israel G. et al., 2005, ApJ, 628, L53 Katz J. I., 1989, MNRAS, 239, 751 Melatos A., Payne D. J. B., 2005, ApJ, 623, 1044 Mestel L., 1956, MNRAS, 116, 324 Monaghan J. J., 1966, MNRAS, 134, 275 Roxburgh I. W., 1966, MNRAS, 132, 347 Saio H., 1982, ApJ, 256, 717 Tomimura Y., Eriguchi Y., 2005, MNRAS, 359, 1117 Ushomirsky G., Cutler C., Bildsten L., 2000, MNRAS, 319, 902 Wasserman I., 2003, MNRAS, 341, 1020 APPENDI X A : UNI FORM DENSI TY STARS As discussed in Section 2, the magnetic field configuration cannot be freely specified for barotropic equations of state. In order to un- derstand how restricted our choice of magnetic field configurations may be it is useful to consider the simplest model, a constant-density star. Insert ρ = constant into (4) to get W (rZ) ′ + 2ZX = 0, (A1) W ( r2Z′′ − 2Z) + 2Z (rX′ − X) = 0, (A2) W [ r2X′′ − 4 (W + X)] + 2Z (rZ′ − Z) = 0. (A3) We will discuss the solution of these equations for three potential configurations, purely poloidal and toroidal fields as well as a (more realistic) mixed field. A1 Poloidal fields Let us first consider purely poloidal fields. In this case we can combine equation (A3) with equation (3) to get r2X′′′ + 4rX′′ − 4X′ = 0 (A4) which has the solution X = D + Cr2 (A5) C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 Modelling magnetically deformed neutron stars 539 with C and D constant and where the ∝ r−3 solution has been excluded to ensure regularity at the centre. This gives, for W, W = −D − 1 2 Cr2. (A6) Note that for C = 0 this corresponds to a uniform dipole field. The same solution can be obtained by solving for the stream function. It is then sufficient to take S(r, θ ) of the form in equation (13) and take β = 0 so that the magnetic field reduces to B = ( 2A cos θ r2 , −A′ sin θ r , 0 ) . (A7) We can solve equation (12) with the boundary conditions that the field must remain finite at the centre, i.e. A r2 , A′ r finite at r = 0 (A8) while at the surface the field must be continuous with an external curl free dipole field, so that we have A r + A′ = 0, at r = R. (A9) The solution to this problem was given by Ferraro (1954), who considered a constant-density star with a current density with φ component of the form J = B R2 r3 sin3 θ, (A10) where B is a constant parameter. In this case equation (12) reduces to ∇2S = B R2 r2 sin2 θ (A11) which yields for the stream function S = B r 2 2R2 ( r2 5 − R 2 3 ) sin2 θ. (A12) The magnetic field is then Br = −B [ 1 3 − 1 5 ( r R )2] cos θ, (A13) Bθ = B [ 1 3 − 2 5 ( r R )2] sin θ, (A14) which corresponds to taking C = −2 5 B R2 and D = B 3 (A15) in the previous solutions (A5) and (A6). A2 Toroidal field For a purely toroidal field the equations reduce to rZ ′ = Z (A16) which gives Z = −icr. (A17) This corresponds to a uniform current distribution inside the star, i.e. |∇ ×B| = 2c= constant, pointing along the z-axis. In fact, the only solution that is consistent with the surface boundary condition is c = 0, i.e. there is no non-trivial toroidal field solution. A3 Mixed poloidal/toroidal field Finally, let us consider the general case. If we formally solve equa- tion (15) for X, then equation (A1) becomes r2 [ W ( Z r )′ − ( Z r ) W ′ ] = r2W 2 ( Z rW )′ = 0. (A18) Hence, unless W = 0 or Z = 0 we have Z = arW, (A19) where a is a constant. Note that equation (A2) is implied by (A1) by the use of (15); thus the substitutions Z = arW and X = −W − r 2 W ′ (A20) solve equations (A1) and (A2). The last equation, (A3), then be- comes 1 2 rW [r2W ′′′ + 4rW ′′ − r(1 + a2r2)W ′] = 0. (A21) Thus we find that the solution can be written as W = 1 r3 [ C1(2ar − 1)e2ar + C2(2ar + 1)e−2ar ] + C3, (A22) leading to X = − 1 2r3 [ C1(4a 2r2 − 2ar + 1)e2ar − C2(4a2r2 + 2ar + 1)e−2ar ] − C3, (A23) where Ci are constants. In order to ensure regularity at the centre we must take C1 = C2, and the central values of the fields (which we shall denote with a subscript ‘c’) then become Wc = −Xc = 16 3 C1a 3 + C3 and Zc = 0. (A24) As we can see, the parameter 1/a is just a length-scale and can be absorbed in the other constants if we define a new dimensionless radial coordinate x = 2ar . This allows us to redefine the constant C1 → C1/(2a)3. (A25) The free parameters in the regular solution now correspond to the new C1 and C3 and are only two (if we exclude the trivial choice of scale or units given by a). Explicitly we have W = 2C1 cosh x x2 − 2C1 sinh x x3 + C3, (A26) X = C1 cosh x x2 − C1(x 2 + 1) sinh x x3 − C3, (A27) Z = C1 cosh x x − C1 sinh x x2 + 1 2 xC3. (A28) We can further interpret the parameters by noting that the central values of the fields are Wc = −Xc = 23 C1+C3(Zc = 0 still). We can thus use C1 as an overall scale for the field and define Ŵ = W/C1 and likewise for the other variables. Thus, using C3/C1 = Ŵc−2/3, we obtain Ŵ = 2 cosh x x2 − 2 sinh x x3 + Ŵc − 2/3, (A29) X̂ = cosh x x2 − (x 2 + 1) sinh x x3 − Ŵc + 2/3, (A30) Ẑ = cosh x x − sinh x x2 + 1 2 x(Ŵc − 2/3). (A31) C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 540 B. Haskell et al. We now have only one non-trivial parameter, Ŵc, which controls the ratio of poloidal to toroidal fields, while the other parameters a and C1 have been absorbed into the definition of our variables. They can simply be interpreted as scales of the problem (specifically a length-scale and the scale of the magnetic field). APPEN D IX B: POLYTRO PIC MODELS The uniform density results in Appendix A suggest that the magnetic field configuration is severely constrained. It is natural to ask if this restriction derives from having taken a somewhat pathological and simplistic stellar model. To establish that this is not the case, we want to allow for a non-uniform density distribution. It is natural to consider an n = 1 polytrope, as this model has the gross features of a ‘realistic’ neutron star equation of state. It also permits an analytic treatment. The allowed magnetic field must still satisfy equation (4), but the role of ρ is no longer trivial as we are considering a density profile ρ(y) = ρc y sin y, (B1) where we have introduced the dimensionless radius y = πr/R and ρc is the density at the centre of the star. B1 Poloidal field Following Monaghan (1996), we solve equation (12), imposing regularity at the centre and matching to an external dipole. The solution for the stream function is then S = 2 sin 2 θ 3y [y3 + 3(y2 − 2) sin y + 6y cos y] (B2) which leads to a field of form Br = Bs cos θ π(π2 − 6) [y 3 + 3(y2 − 2) sin y + 6y cos y], (B3) Bθ = Bs sin θ 2π(π2 − 6) [−2y 3 + 3(y2 − 2) (sin y − y cos y)], (B4) where Bs represents the conventional dipole field strength at the surface. B2 Mixed toroidal and poloidal fields Let us now assume a mixed poloidal and toroidal configuration. We will once again look for dipolar solutions and take a stream function of the form S(r, θ ) = A(r) sin2 θ . Furthermore, following Roxburgh (1966) we shall define β to be β = πλ R S. (B5) This is obviously a particular choice. However, from (11) we see that it represents the leading order behaviour for small β, i.e. weak fields. Moreover, this choice leads to a separable equation which means that the problem can be treated more or less analytically. By varying the parameter λ, which describes the relative strength of the toroidal part of the field, one can hope to get some insight into how the two field components interact in the mixed problem. The magnetic field thus takes the form B = ( 2A cos θ r2 , −A′ sin θ r , πλA sin θ rR ) . (B6) Then equation (12) can be written as A d dr [ 1 ρr2 ( 2A r2 − d 2A dr2 ) − π 2λ2A ρR2r2 ] = 0. (B7) We want to find a solution such that the field remains finite at the centre, cf. (A8), and that all the components of the magnetic field are continuous at the surface. As the toroidal field must vanish at the surface this forces the condition A = 0 at r = R. (B8) It must also be the case that A′ = 0 at r = R (B9) which derives from the condition that all components of the field must vanish at the surface. We can thus only consider fields that vanish outside the star. Solving the above equations for the density profile of an n = 1 polytrope and using (A8) and (B8), we find that the stream function takes the form A = BkR 2 (λ2 − 1)2y { 2π λy cos(λy) − sin(λy) πλ cos(πλ) − sin(πλ) + [(1 − λ2)y2 − 2] sin y + 2y cos y}, (B10) where Bk parametrizes the strength of the magnetic field and λ �= 1. The case λ = 1 is special, and it is easy to show that A cannot satisfy all the boundary conditions for this case. Imposing the remaining boundary condition (B9) we find an eigenvalue rela- tion for λ. Permissible fields must be such that πλ(λ2 − 1) cos(πλ) − (3λ2 − 1) sin(πλ) = 0 (λ �= 1). (B11) By taking higher values of λ we can increase the relative strength of the toroidal part of the field compared to the poloidal part. Solving the transcendental equation numerically we find that the first three eigenvalues are λ1 = 2.362, λ2 = 3.408, λ3 = 4.430. These results agree with the values given by Roxburgh (1966) and Ioka (2001). For large λ it is easy to show that the roots are well approximated by λn = n + 3/2. The key point is that the toroidal part of the field is not freely specifiable in this model. If we prescribe the strength of the poloidal field component, then we can find solutions with increasing toroidal fields, but these are discrete. B3 Purely toroidal fields In the case of a purely toroidal field equation (4) takes the form ∂ ∂r ( Bφ ρr sin θ ) ∂ ∂θ (Bφr sin θ ) − ∂ ∂θ ( Bφ ρr sin θ ) ∂ ∂r (Bφr sin θ ). (B12) The boundary conditions we have to impose are that the field vanish both at the centre and at the surface of the star, as before. We can take a solution of the simple form Bφ = B π sin y sin θ (B13) which will satisfy the boundary conditions since ρ = 0 at the surface (y = π) for the polytropic model we are considering. Note that, in contrast to the uniform density problem, the polytropic equation of C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 Modelling magnetically deformed neutron stars 541 state allows for a purely toroidal field. This is a natural consequence of the fact that the density vanishes at the surface of a polytropic star which makes it possible to satisfy the boundary conditions in this case. B4 Field confined to the neutron star ‘crust’ It has sometimes been assumed that, if the core of a neutron star is a Type I superconductor the magnetic field will be expelled from the core and confined to the crust, see e.g. Bonazzola & Gourgoulhon (1996). In reality this is not very likely to happen. Basically, the flux expulsion will not be completed on a time-scale relevant for observed neutron stars. The true structure will be more complex depending on the nucleation of macroscopic superconducting re- gions. The end result may well be that a Type I superconducting core is similar to a normal magnetized fluid. If this is the case, then our general analysis would apply. Still, it is interesting to con- sider the extreme possibility of complete flux expulsion. To discuss this situation one should in principle consider the full equations of hydrostatic equilibrium, including the elastic terms. However, in order to obtain some initial estimates for this problem, we will con- sider the case of a fluid with a magnetic field confined to a region close to the surface. This model would be relevant provided that the deformed shape represents the relaxed configuration of the crust. Then there is no strain, and thus no elastic terms in the equilibrium equations. We shall consider the same mixed field as in the previous section, i.e. of the form (B6). We shall, however, require the field to vanish inside a certain radius rb (which can be considered to be the base of the crust). The boundary conditions for the third-order dif- ferential equation (B7) at the surface thus remain A(R) = A ′ (R) = 0 together with a condition which comes from imposing continuity of the Br component of the magnetic field, i.e. imposing A(rb) = 0. (B14) We again get an eigenvalue problem for the parameter λ, allowing us to calculate the permitted ratios of toroidal to poloidal field strengths. It should be noted that we are not imposing continuity of the tractions at the inner boundary. In fact, there is a discontinuity in the Bθ component of the field at rb, which will lead to currents at the crust/core interface. In order to simplify the calculation we also take the core to be unperturbed, and simply impose continuity of the perturbation in the gravitational potential δ� and of its derivative δ�′. At best, this simple model can be seen as a rough representation of the problem with a superconducting core. There are a number of is- sues that one ought to worry about, and which we are not considering here. The main question concerns the true nature of a Type I super- conductor and whether one should expect that the magnetic flux is expelled from the entire core rather than (say) bunched up into much smaller, but still macroscopic, regions? Moreover, if the core forms a Type II superconductor (as is usually expected) then the magnetic flux is carried by flux tubes. This problem is significantly different from our idealized model and requires a separate analysis. We will not consider this problem here. The interested reader will find a useful discussion in Wasserman (2003) and Akgün & Wasserman (2008), who analysed the case of toroidal magnetic fields in a neu- tron star with a Type II superconducting region. A P P E N D I X C : R A D I A L D E F O R M AT I O N S As discussed in Section 3, it may be relevant to work out also the radial deformation due to the presence of the magnetic field. For general fluid configurations, any quadrupole deformation is likely to be accompanied by an l = 0 component. In order to determine the radial deformation we can assume that the new radial variable x defined in equation (28) labels the deformed gravitational equipo- tential surfaces. This means that we impose ∇x × ∇� = 0 (C1) which, if we write ε(x, θ) = D0(x) + D2(x)P2(θ ), (C2) leads to the condition D2 = − √ 4π 5 [ d� dr ]−1 δ� r . (C3) This expression allows us to calculate the quantity D2 throughout the star and can be of use, for example, when discussing oscilla- tions on a deformed background, cf. Saio (1982). Writing out the equations that need to be solved in the l = 0 case, the perturbed Euler equation and Poisson equation, we have d dr δ�0 = ψ0, (C4) d dr ψ0 = 4πGδρ0 − 2 r ψ0, (C5) d dr δρ0 = ρ P { L0(B) − ρψ0 − δρ0 [ d dr ( P ρ ) 1 + d dr � ]} , (C6) where ψ0 is introduced in order to obtain a first-order system. We have also assumed a linearized barotropic equation of state δP P = 1 δρ ρ (C7) and L0(B) is the l = 0 part of the Lorentz force, which can be cal- culated once the magnetic field has been specified. These equations can now be integrated, imposing regularity at the centre of the star and imposing that at the surface δρ0 = 0 and matching δ�0 and δ�′0 to an exterior solution of the form δ�ext0 ∝ 1/r. This will allow us to calculate D0 = − √ 4π [ d� dr ]−1 δ�0 r . (C8) APPENDI X D : G ENERAL DEFORMATI ONS : COMBI NI NG MAG NETI C FI ELDS A N D ROTAT I O N The framework discussed in the main body of the paper can easily be extended to the case, where the deformation is due not only to a magnetic field, but also to the star’s rotation. We can then write the equations of motion as ∇p ρ + ∇ψ = (∇ × B) × B 4πρ = L 4πρ , (D1) where ψ = � − 1/2�2r2 sin2 θ. (D2) � is the gravitational potential and � the constant rotation rate. As the new term due to rotation is written as the gradient of a scalar function, its curl will still vanish and the compatibility condition for the magnetic field in equation (4) remains the same. We can thus C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022 542 B. Haskell et al. still use the magnetic field from equation (B6). The equations of hydrostatic equilibrium, for the l = 2 case, now take the form dδp dr + ρ dδ� dr + δρ d� dr + 2 3 ρ�2r = Lr, (D3) δp + ρδ� + ρ� 2r2 3 = rLθ , (D4) where Lr and Lθ are given in equations (61) and (62). We can proceed as in the previous analysis and use equation (D4) to obtain δρ as a function of δ�, B2 and �. This allows us to write the perturbed Poisson equation (for l = 2) as d2δ� dr2 + 2 r dδ� dr − δ� { 6 r2 + 4πGdρ dr [ d� dr ]−1} =−4πG 3r5 [ d� dr ]−1[ 2r dA dr A(π2λ2r2 − 2) − 4A2(π2λ2r2 − 4) +2r3Ad 3A dr3 − 4r2Ad 2A dr2 − dρ dr �2r7 ] (D5) and thus obtain the surface deformation ε and the ellipticity � as defined in equations (66) and (67). As an example let us take an n = 1 polytrope with the purely poloidal field from equation (B4). The equations of hydrostatic equilibrium now include the rotation term as in equation (D4), where the Lorentz force is that of equation (46). We can thus apply exactly the same procedure as previously and solve for the perturbed quantities δ�, δρ and δp. This allows us to calculate the deformation at the surface ε = − π 5B2(24 − π2) 16R2ρ2c G(π 2 − 6)2 − 5 4π �2 Gρc . (D6) Working out the corresponding ellipticity we find � ≈ 8 × 10−11 ( R 10 km )2 ( M 1.4 M� )−4 ( B̄ 1012G )2 + 1.4 × 10−1 ( � �br )2 , (D7) where �br = (2/3) √ πGρ̄ is the break-up frequency (ρ̄ is the aver- age density of the non-rotating model). For our values this frequency would be f≈ 1250 Hz, which corresponds to a period P ≈ 0.8 ms. We see that, as expected, the rotational term completely dominates the magnetic term. Rotational deformations will, however, always be axisymmetric. The magnetic deformations can lead to a quadrupolar deformation and thus to gravitational wave emission if the magnetic axis is inclined with respect to the spin axis. This paper has been typeset from a TEX/LATEX file prepared by the author. C© 2008 The Authors. Journal compilation C© 2008 RAS, MNRAS 385, 531–542 D ow nloaded from https://academ ic.oup.com /m nras/article/385/1/531/1033831 by R adboud U niversity N ijm egen user on 15 June 2022