Neutron stars with a toroidal magnetic field

# Equilibrium models of relativistic stars with a toroidal magnetic field

J. Frieben and L. Rezzolla
Max-Planck-Institut für Gravitationsphysik, Albert-Einstein-Institut, Am Mühlenberg 1, D-14476 Golm, Germany
Department of Physics, Louisiana State University, Baton Rouge, LA 70803, USA
###### Abstract

We have computed models of rotating relativistic stars with a toroidal magnetic field and investigated the combined effects of magnetic field and rotation on the apparent shape (i.e. the surface deformation), which could be relevant for the electromagnetic emission, and on the internal matter distribution (i.e. the quadrupole distortion), which could be relevant for the emission of gravitational waves. Using a sample of eight different cold nuclear physics equations of state, we have computed models of maximum field strength, as well as the distortion coefficients for the surface and the quadrupolar deformations. Surprisingly, we find that non-rotating models admit arbitrary levels of magnetization, accompanied by a growth of size and quadrupole distortion to which we could not find a limit. Rotating models, on the other hand, are subject to a mass-shedding limit at frequencies well below the corresponding ones for unmagnetized stars. Overall, the space of solutions can be split into three distinct classes for which the surface deformation and the quadrupole distortion are either: prolate and prolate, oblate and prolate, or oblate and oblate, respectively. We also derive a simple formula expressing the relativistic distortion coefficients, which allows one to compute the surface deformation and the quadrupole distortion up to significant levels of rotation and magnetization, essentially covering all known magnetars. Such a formula replaces Newtonian equivalent expressions that overestimate the magnetic quadrupole distortion by about a factor of 6 and are inadequate for strongly relativistic objects like neutron stars.

###### keywords:
gravitational waves – stars: magnetic field – stars: neutron.
pagerange: Equilibrium models of relativistic stars with a toroidal magnetic fieldApubyear: 2012

## 1 Introduction

In addition to rapid rotation, also strong magnetic fields can introduce significant deformations in neutron stars, as shown, for instance, by Bocquet et al. (1995) and Cardall, Prakash & Lattimer (2001), who have computed fully non-linear models of relativistic stars with a poloidal magnetic field. At the end of the collapse of the core of a massive star, differential rotation could create strong toroidal magnetic fields of the order of inside the hot proto-neutron star (Bonanno, Rezzolla & Urpin 2003; Naso et al. 2008; Bonazzola & Haensel, unpublished). As a result, realistic models of magnetized relativistic stars require the simultaneous presence of both poloidal and toroidal field components. As pointed out in Gourgoulhon & Bonazzola (1993), however, this requires a formalism capable of dealing with non-circular space–times (i.e. with convective currents in the meridional planes) like the one presented here, but which has so far been implemented only in a perturbative scheme (Ioka & Sasaki, 2003, 2004). Triggered also by the increasing interest in strongly magnetized neutron stars due to their relation with soft-gamma repeaters and anomalous X-ray pulsars (Duncan & Thompson, 1992; Thompson & Duncan, 1996), a growing number of studies adopting perturbative techniques have appeared investigating either the field geometry and neglecting the influence of the magnetic field on the matter distribution (Ciolfi et al., 2009), or solving the coupled Einstein–Maxwell–Euler system, from which the magnetic deformation can be calculated (Ioka & Sasaki 2004; Colaiuda et al. 2008; Ciolfi, Ferrari & Gualtieri 2010; Gualtieri, Ciolfi & Ferrari 2011; Yoshida, Kiuchi & Shibata 2012). Until recently, however, fully non-linear models of relativistic magnetized stars were restricted to purely poloidal magnetic fields (Bocquet et al., 1995; Cardall et al., 2001), for which the generated space–time is circular like in the unmagnetized case (Carter, 1973). Following the recent insight that also a magnetic field with only a toroidal component is compatible with the circularity of space–time (Oron, 2002), studies of relativistic models of rotating stars with a toroidal magnetic field have emerged (Frieben & Rezzolla 2007; Kiuchi & Yoshida 2008; Kiuchi, Kotake & Yoshida 2009; Yasutake, Kiuchi & Kotake 2010; Yasutake, Maruyama & Tatsumi 2011). These new studies have complemented earlier Newtonian investigations (Sinha, 1968; Sood & Trehan, 1972; Miketinac, 1973), which being simpler, allowed for the investigation of more complex field geometries, in particular of mixed poloidal and toroidal magnetic fields (Tomimura & Eriguchi 2005; Yoshida, Yoshida & Eriguchi 2006; Haskell et al. 2008; Lander & Jones 2009; Fujisawa, Yoshida & Eriguchi 2012).

In this work, we are mainly concerned with the deformation of neutron stars endowed with a toroidal magnetic field because they might be an important source of gravitational radiation due to the prolate deformation induced by the magnetic field and provided that the axes of symmetry and of rotation are different (Bonazzola & Gourgoulhon, 1996). Furthermore, Jones (1975) has pointed out that viscous processes can trigger a secular instability, which drives the axis of symmetry of the prolate star into the plane perpendicular to the angular momentum vector transforming it into a bar-shaped rotating source of gravitational radiation (Cutler, 2002; Stella et al., 2005). From the astrophysical point of view, both the deformation of the stellar surface as well as the distortion of the matter distribution are relevant and can be measured by appropriate quantities, namely the surface deformation (or apparent oblateness) and the quadrupole distortion. Previous studies of relativistic models of stars with a toroidal magnetic field have provided a broad survey of non-rotating and rotating models for varying masses, radii and magnetic fluxes. The focus of this work is a complementary one to earlier studies: in order to obtain a comprehensive picture of the impact of a toroidal magnetic field on relativistic stars, we exclusively study models of fixed baryon mass corresponding to a gravitational mass of in the unmagnetized and non-rotating case. Although we expect neutrons stars to come in a narrow but non-zero range of masses, restricting to a single value of the gravitational mass has the important advantage that we can explore with unprecedented precision both the deformations introduced by magnetic fields and those introduced by rapid rotation. As we will comment later on, our increased accuracy has allowed us also to discover novel results and equilibrium configurations.

Our neutron stars are modelled assuming the matter to be a perfect fluid at zero temperature well described by a single-parameter equation of state (EOS), and as having infinite conductivity, as required by the ideal magnetohydrodynamics (MHD) limit. We do not consider multifluid models, which would allow for the stratification of neutron star matter, nor do we treat the protons in the interior as a superconducting fluid, which would greatly alter the magnetic properties of associated equilibrium models. Nevertheless, important results obtained for extremely magnetized models with field strengths of the order of can be readily extended to configurations of much lower and more realistic field strengths of the order of . More specifically, in addition to a standard polytropic EOS, we have considered a sample of seven realistic EOSs resulting from calculations of cold catalysed dense matter, namely the APR EOS (Akmal, Pandharipande & Ravenhall, 1998), the BBB2 EOS (Baldo, Bombaci & Burgio, 1997), the BN1H1 EOS (Balberg & Gal, 1997), the BPAL12 EOS (Bombaci, 1996), the FPS EOS (Pandharipande & Ravenhall, 1989), the GNH3 EOS  (Glendenning, 1985) and the SLy4 EOS (Douchin & Haensel, 2001). Altogether, this set of EOSs spans a wide range of physical properties and should cover any realistic description of neutron stars. For the Pol2 EOS, we have explored systematically the space of solutions and computed the corresponding surface deformation and quadrupole distortion. In addition, we have also computed the distortion coefficients, which allow one to compute the deformation of neutron stars up to large magnetizations and rotation rates through a simple algebraic expression, following the procedure devised in Cutler (2002).

The plan of the paper is the following. In Section 2, we give an overview of the novel results obtained in this study, in Section 3, we discuss the theoretical framework on which our approach, whose numerical implementation and testing is discussed in Section 4, is based. Results for static magnetized models are presented in Section 5 and for rotating magnetized models in Section 6. In Section 7, we deal with the case of moderate magnetic field and rotation and derive empirical distortion coefficients before presenting our conclusions in Section 8.

## 2 General Overview

Given the complexity of the numerous results found and the risk that the most important ones may be lost in the details, in the following we briefly summarize what we believe are the most salient properties of (relativistic) stars with purely toroidal magnetic fields. We recall that we measure with the deformation of the surface shape, while we measure with the (quadrupolar) deviation of the matter distribution from a spherically symmetric one. Overall, equilibrium models of relativistic stars with a toroidal magnetic field exhibit the following properties:

1. Non-rotating and magnetized models exhibit a prolate surface deformation and a prolate quadrupolar deformation, i.e. , , both of which decrease as the magnetization parameter is increased.

2. Rotating and unmagnetized models exhibit an oblate surface deformation and an oblate quadrupolar deformation, i.e. , , both of which increase as the rotation frequency is increased.

3. Between these limiting cases, neutral lines and divide models into having prolate/oblate surface deformations and prolate/oblate internal deformations, respectively.

4. For non-rotating models no upper limit was found to the magnetization parameter , with stellar models that become increasingly prolate, but also increasingly extended as the magnetization is increased (see Figs 6 and 7).

5. The magnetic pressure associated with the toroidal magnetic field also causes an expansion in the outer layers of the star, in particular a growth of its equatorial radius. This effect is present also for non-rotating models, for which, however, the polar radius grows more rapidly than the equatorial one, yielding a prolate surface deformation, i.e. . However, as the rotation is increased and the magnetization decreased, models can be found which have a prolate interior deformation, i.e. , and an oblate surface deformation, i.e. .

6. For any given angular velocity, the model of maximum magnetization coincides with the mass-shedding model, i.e. the model for which centrifugal and gravitational forces are equal at the equatorial radius. This point is characterized by the appearance of a cusp at the equator.

7. For rotating models the magnetic pressure at the equator adds to the centrifugal force, favouring the loss of mass. As a result of the increase in the equatorial size, the mass-shedding frequency is systematically smaller than the corresponding one for unmagnetized models.

8. In the space of parameters considered, the average toroidal magnetic field strength is not a monotonic function of the intrinsic magnetization parameter , and after reaching a maximum value, it gradually decreases. This implies the existence of double solutions in certain parts of the space of parameters .

9. Although the space of parameters is potentially degenerate, the corresponding space of parameters is not. As a result, any stellar model considered is characterized uniquely by the values of the angular velocity and of the magnetization parameter .

Given these results, it is natural to divide our models into three classes: (1) models labelled PP for prolate–prolate, for which both apparent shape and matter distribution are prolate, i.e. , ; (2) models labelled PO for prolate–oblate, whose shape is oblate whereas their matter distribution is prolate, i.e. , ; (3) models labelled OO for oblate–oblate, which appear oblate and also exhibit an oblate matter distribution, i.e. , . The latter class of models had not been found by Kiuchi & Yoshida (2008). A schematic picture illustrating the three different classes for models below the maximum field strength limit is shown in Fig. 1. Finally, while different EOSs with their different stiffness introduce quantitative differences in the behaviour described above, they all follow the same qualitative behaviour.

## 3 Mathematical Setup

### 3.1 Basic assumptions

We assume that the space–time generated by the rotating star is stationary and axisymmetric, with Killing vector fields and associated with these symmetries. If, in addition, the space–time is asymptotically flat and there exists an axis where vanishes, then and commute (Carter, 1970). This enables us to choose coordinates with vector fields and . If furthermore the total stress–energy tensor satisfies the circularity conditions

 \mn@boldsymbolT⋅\mn@boldsymbole0=α\mn@boldsymbole0+β\mn@boldsymbole3, (1)
 \mn@boldsymbolT⋅\mn@boldsymbole3=λ\mn@boldsymbole0+μ\mn@boldsymbole3, (2)

convective currents in the meridional planes of constant are absent by construction. In this case, quasi-isotropic (QI) coordinates can be adopted and the line element reads

 ds2=gαβdxαdxβ=−N2dt2+Φ2r2sin2θ2(dϕ−Nϕdt)2+Ψ2(dr2+r2dθ2), (3)

where , , , and are functions of . We further introduce the Eulerian observer whose 4-velocity is the future-directed unit vector normal to hypersurfaces of constant . From equation (3), we infer .

The compatibility of the electromagnetic fields with the circularity condition was established long ago (Carter, 1973) for the case in which the electromagnetic field tensor is derived from a potential 1-form with components , whereas the contravariant components of the electric current vector read . This fact has been exploited for computing fully relativistic models of stars with a poloidal electromagnetic field (Bocquet et al., 1995; Cardall et al., 2001). However, it has been shown that the case of a toroidal magnetic field satisfies the circularity assumption too (Oron, 2002; Kiuchi & Yoshida, 2008). In the following, we adopt a vector potential which is orthogonal to the Killing vectors and , namely and . The covariant components of then become

 Aα=(0,Ar,Aθ,0). (4)

This particular form of ensures, by construction, the absence of any poloidal electric or magnetic field component. The only non-vanishing component of the antisymmetric Faraday tensor is then given by

 Frθ=∂Aθ∂r−∂Ar∂θ. (5)

For the electric field and the magnetic field as measured by observer , we then obtain

 Eα=Fαβnβ=(0,0,0,0), (6)
 Bα=−\normalsize12ηαβγδFγδnβ=ΦsinθΨ2(Nϕ[∂Ar∂θ−∂Aθ∂r],0,0,−[∂Ar∂θ−∂Aθ∂r]), (7)

where is the totally antisymmetric tensor associated with the metric . According to equations (6) and (7), the combination of a vanishing electric field and a toroidal magnetic field holds for any observer whose 4-velocity is a linear combination of the two Killing vectors and . The electromagnetic contribution to the stress–energy tensor is given by

 \mn@boldsymbolTαβ=14π(FακFκβ−14FκλFκλgαβ). (8)

Following the procedure adopted by Bonazzola et al. (1993), is split up into the total electromagnetic energy density , the Poynting 3-vector , and the electromagnetic stress 3-tensor as measured by the Eulerian observer . With the projection tensor , these quantities are obtained as projections of on to and orthogonal to , namely , , and . Specialized to the present case of no electric field and no poloidal magnetic field components, the electromagnetic contribution to the (3+1) matter variables [for an introduction to the (3+1) formalism of general relativity, cf.  Smarr & York 1978; Gourgoulhon 2012], namely the total energy density , the momentum density 3-vector , and the stress 3-tensor reads

 E=18π(BϕΦrsinθ)2, (9)
 Ji=0, (10)
 Srr=E,Sθθ=E,Sϕϕ=−E. (11)

In particular, all non-diagonal components of are zero, and . The circularity assumption requires that for the Poynting vector

 Jr=0,Jθ=0, (12)

and that the electromagnetic stress tensor satisfies

 Srϕ=0,Sθϕ=0. (13)

### 3.2 Einstein equations

The Einstein equations for the metric tensor defined by equation (3) and a general stress–energy tensor decomposed into the (3+1) quantities , and become

 Δν=4πΨ2(E+S)+Φ2r2sin2θ2N2(∂Nϕ)2−∂ν∂(ν+β), (14)
 ~Δ3(Nϕrsinθ)=−16πNΨ2Φ2Jϕrsinθ−rsinθ∂Nϕ∂(3β−ν), (15)
 Δ2[(NΦ−1)rsinθ]=8πNΨ2Φ(Srr+Sθθ)rsinθ, (16)
 Δ2ζ=8πΨ2Sϕϕ+3Φ2r2sin2θ4N2(∂Nϕ)2−(∂ν)2, (17)

where the following abbreviations have been used:

 ν≡lnN,ζ≡ln(NΨ),β≡lnΦ, (18)
 Δ2≡∂2∂r2+1r∂∂r+1r2∂2∂θ2, (19)
 Δ3≡∂2∂r2+2r∂∂r+1r2∂2∂θ2+1r2tanθ∂∂θ, (20)
 ~Δ3≡Δ3−1r2sin2θ, (21)
 ∂a∂b≡∂a∂r∂b∂r+1r2∂a∂θ∂b∂θ. (22)

The resulting system of four non-linear elliptic equations for the metric variables , , and can be solved iteratively once suitable boundary conditions of asymptotic flatness have been adopted. Additional details can be found in Bonazzola et al. (1993).

### 3.3 Maxwell equations

Since the electromagnetic field tensor is derived from a potential 1-form , the homogeneous Maxwell equations

 Fαβ;γ+Fβγ;α+Fγα;β=0 (23)

are satisfied by construction. The inhomogeneous Maxwell equations allow us to express the electric current 4-vector in terms of the Faraday tensor , where the alternative expression

 4πjα=1√−g(√−gFαβ),β (24)

with is used. The electromagnetic field tensor has only one non-vanishing contra- and covariant component, which can now be expressed in terms of the azimuthal component of the magnetic field ,

 Frθ=Ψ2ΦsinθBϕ, (25)
 Frθ=BϕΨ2Φr2sinθ. (26)

Combining equations (24) and (25), the poloidal components of the electric current 4-vector can be written as

 4πjr=1Ψ2ΦNr2sinθ∂(NBϕ)∂θ, (27)
 4πjθ=−1Ψ2ΦNr2sinθ∂(NBϕ)∂r. (28)

The remaining components and are zero, as expected. Note that the circularity condition, equations (1) and (2), forbids any meridional convective current contributing to , which is ensured by assuming the charge carriers to be massless. By taking the divergence of equation (24), it follows that thanks to the antisymmetry of , and this continuity equation for the electric current leads to a dependence between and , namely

 \mn@boldsymbol∇⋅\mn@boldsymbolj=1√−g[∂∂r(√−gjr)+∂∂θ(√−gjθ)]=0, (29)

which however is trivially fulfilled by virtue of equations (27) and (28) which express and as functions of a single quantity, namely , without any restriction from equation (29). A useful consequence of equations (27) and (28) is that the flow lines of the electric current coincide with the isocontours of , which allows one to visualize the current distribution without actually computing .

### 3.4 Equation of motion and condition on the Lorentz force

In the case of a perfect fluid the stress–energy tensor takes the form

 \mn@boldsymbolT=(e+p)\mn@boldsymbolu⊗\mn@boldsymbolu+p\mn@boldsymbolg, (30)

where is the energy density and the pressure as measured by the fluid comoving observer with 4-velocity . Since we assume the absence of meridional currents, can be written as a linear combination of the Killing vectors and , namely

 \mn@boldsymbolu=ut\mn@boldsymbole0+uϕ\mn@boldsymbole3. (31)

Introducing the Lorentz factor linking observers and and the fluid coordinate angular velocity , the physical fluid velocity and the Lorentz factor can be expressed, respectively, as

 U=ΦrsinθN(Ω−Nϕ), (32)
 Γ=(1−U2)1/2. (33)

Taking into account the contribution of the magnetic field to the (3+1) matter variables according to equations (9)–(11), the following expressions are obtained for a perfect fluid endowed with a toroidal magnetic field:

 E=Γ2(e+p)−p+18π(BϕΦrsinθ)2, (34)
 Jϕ=Γ2(e+p)ΦrsinθU, (35)
 Srr=Sθθ=p+18π(BϕΦrsinθ)2, (36)
 Sϕϕ=p+Γ2(e+p)U2−18π(BϕΦrsinθ)2. (37)

All other components are zero, and . The projection of on to spatial hypersurfaces of constant provides the equation of motion

 1e+p∂p∂xi+∂ν∂xi−∂∂xi(lnΓ)−1e+pFiβjβ=−F∂Ω∂xi, (38)

where is defined as

 F≡uϕut=Γ2ΦNUrsinθ. (39)

Assuming a single-parameter EOS, and , with the baryon number density, and zero temperature, we introduce the fluid log-enthalpy

 H≡ln(e+pnmb), (40)

with a mean baryon rest mass of . Using equations (25)–(28), (40), and discarding the case of differential rotation, equation (38) becomes

 ∂∂xi(H+ν−lnΓ)+(1e+p)(18πΦ2N2r2sin2θ)×∂(NBϕ)2∂xi=0. (41)

Equation (41) can only be satisfied if the Lorentz force term is derived from a scalar function as well, thus we require

 (1e+p)(18πΦ2N2r2sin2θ)∂(NBϕ)2∂xi=∂~M∂xi. (42)

Using the Schwartz theorem, the integrability condition of equation (42) can be written as

 ∂G∂θ∂(NBϕ)∂r−∂G∂r∂(NBϕ)∂θ=0, (43)

where . In other words, equation (43) states that the Jacobian of the coordinate transform is zero, and hence there exists a scalar function such that . Two different cases can be distinguished. (i) If , then does not depend on , and the relation implies that is constant, which is equivalent to the absence of any matter and . (ii) If , then the implicit-function theorem enables us to write . Discarding the case without matter, we retain possibility (ii) and conclude that

 NBϕ=g((e+p)Φ2N2r2sin2θ), (44)

with being an arbitrary scalar function. However, the regularity properties of an axisymmetric vector field in the case of spatial spherical coordinates impose some further restrictions on . For the covariant component of an azimuthal vector field , it can be shown (Bardeen & Piran, 1983) that

 Uϕ(r,θ)=r2sin2θm(r,θ), (45)

where is an arbitrary axisymmetric scalar function. Application of equation (45) to allows us to conclude that can be written as , with an arbitrary regular scalar function . The resulting expression for after application of equation (45) to equation (44) reads

 NBϕ=(e+p)Φ2N2r2sin2θf((e+p)Φ2N2r2sin2θ). (46)

By inserting equation (46) into equation (42), the gradient of the magnetic potential reads

 ∂~M∂xi=f4π∂∂xi[(e+p)Φ2N2r2sin2θf]=Ri. (47)

In general, cannot be determined in closed form. However, inspection of equation (47) reveals that for the sub-case of a monomial function (not to be meant as a summation over repeated indices), a solution can be written down immediately. In this case, equation (46) simplifies to

 NBϕ=λm((e+p)Φ2N2r2sin2θ)m+1, (48)

and the solution to equation (47) is given by an algebraic expression, namely

 ~M=λ2m4π(m+12m+1)((e+p)Φ2N2r2sin2θ)2m+1. (49)

The simplest function is obtained for and corresponds to a constant function of value . In this case, equation (46) reduces to

 NBϕ=λ0(e+p)Φ2N2r2sin2θ, (50)

and thus

 Bϕ=λ0(e+p)Φ2Nr2sin2θ, (51)

supplemented by equation (49), which yields the simplified expression

 ~M=λ204π(e+p)Φ2N2r2sin2θ. (52)

For all models computed in this study, the magnetic field component is given by equation (51) and the magnetic potential by equation (52). Because vanishes on the axis of symmetry, , and the integral of equation (41) reads

 H+ν−lnΓ+~M=Hc+νc, (53)

where and are the central values of and , respectively. Note that a general prescription for the determination of the magnetic field and the magnetic potential can be found in Kiuchi & Yoshida (2008) and Gourgoulhon et al. (2011).

### 3.5 Perfect-conductor relation

According to Ohm’s law and assuming an infinite conductivity for the neutron star matter, the electric field as measured by the fluid comoving observer has to vanish, namely

 E′α=Fαβuβ=(0,0,0,0). (54)

In the case of a purely poloidal magnetic field (Bocquet et al., 1995), the perfect-conductor relation equation (54) is non-trivial as it establishes a dependence between and and forces a dependence of the angular velocity on , i.e. . In the case of a purely toroidal magnetic field, on the other hand, the perfect-conductor condition is trivially satisfied and no condition is set on the rotation law, so that differentially rotating models can be built as in the unmagnetized case.

It is thus possible to consider differentially rotating configurations following the procedure for the unmagnetized case. The magnetic field as measured by the fluid comoving observer now becomes

 B′α=−\normalsize12ηαβγδFγδuβ=ΓΦsinθΨ2(Ω[∂Aθ∂r−∂Ar∂θ],0,0,−[∂Ar∂θ−∂Aθ∂r]). (55)

Note that and that .

### 3.6 Equation of state

To solve the system of equations introduced in the previous sections, we need a prescription relating the energy density and pressure to the baryon number density . The simplest of such relations is offered by the polytropic EOS, in which case the following identities hold:

 e(n)=mbn+κnγγ−1, (56)
 p(n)=κnγ, (57)

where is the polytropic constant and the adiabatic index. Combining equation (40), equations (56) and (57), the log-enthalpy and can be expressed as functions of each other, namely, as

 H(n)=ln(1+κmbγγ−1nγ−1), (58)
 (59)

Hereafter, we will assume and refer to this EOS as to Pol2. In addition to the Pol2 EOS, we have considered a sample of seven single-constituent one-parameter EOSs treating the neutron star matter as a perfect fluid and derived from very different models of the ground state of cold (zero temperature) dense matter in this work as they are provided by the lorene library.

More specifically, we have considered the APR EOS (Akmal et al., 1998), the BBB2 EOS (Baldo et al., 1997), the BN1H1 EOS (Balberg & Gal, 1997), the BPAL12 EOS (Bombaci, 1996), the FPS EOS (Pandharipande & Ravenhall, 1989), the Sly4 EOS (Douchin & Haensel, 2001), and the GNH3 EOS (Glendenning, 1985). All of these realistic EOSs are provided in tabulated form, which requires the interpolation of listed thermodynamical quantities , , and between contiguous sampling points. Thermodynamical consistency of the interpolated values is ensured by an interpolation procedure based on Hermite polynomials introduced by Swesty (1996), which is crucial for the accuracy of the resulting numerical models (Salgado et al., 1994).

The properties of the spherical non-rotating and unmagnetized reference models with a gravitational mass of can be found in Table 5. The softest EOS of this sample is the BPAL12 EOS, which yields the smallest circumferential radius of , whereas the stiffest one is the GNH3 EOS, which yields the largest circumferential radius of .

### 3.7 Global quantities

A number of global quantities which characterize the numerical models presented in this study can be computed and will be needed in the course of this investigation. These are given as follows. The gravitational mass:

 M≡∫NΨ2Φ(E+S+2NNϕJϕ)r2sinθdrdθdϕ, (60)

the total angular momentum:

 J≡∫Ψ2ΦJϕr2sinθdrdθdϕ, (61)

the total magnetic energy:

 M≡∫Ψ2ΦEr2sinθdrdθdϕ, (62)

the total kinetic energy:

 T≡\normalsize12ΩJ (63)

and the gravitational binding energy:

 W≡M−T−Mp−M, (64)

where the total proper mass is defined as

 Mp≡∫Ψ2ΦΓer2sinθdrdθdϕ, (65)

while the total baryon mass of the star is given by

 Mb≡∫Ψ2ΦΓnr2sinθdrdθdϕ. (66)

Other important quantities that are more closely related to the deformation of the star are the circumferential radius , which is defined through the stellar equatorial circumference as measured by the observer and is related to the equatorial coordinate radius according to

 Rcirc≡Φ(R,π/2)R, (67)

the surface deformation (or apparent oblateness)

 ϵs≡re/rp−1, (68)

where , are the equatorial and polar coordinate radii, respectively, and the quadrupole distortion of the star (Bonazzola & Gourgoulhon, 1996):

 ϵ≡−32IzzI, (69)

where is the quadrupole moment measured in some asymptotically Cartesian mass-centred coordinate system (Thorne, 1980), while the moment of inertia is defined as . We stress that no moment of inertia other than the latter can be defined in a meaningful way for axisymmetric rotating bodies and that the rotation must be rigid and about the axis of symmetry. Furthermore, we note that differs from the standard quadrupole moment and that in the case of QI coordinates, the latter can be read off from the asymptotic expansion of according to

 lnN=−Mr+M212r3+Qr3P2(cosθ)+⋯. (70)

Thorne’s quadrupole moment is the relevant quantity when the gravitational-wave emission from a distorted star has to be determined, and the corresponding values are extracted from the asymptotic expansion of certain components of the metric tensor in QI coordinates following the procedure presented in Bonazzola & Gourgoulhon (1996).

## 4 Numerical Implementation and Tests

### 4.1 Numerical scheme

The analytic scheme presented in Section 3 has been implemented by extending an existing multidomain, surface-adaptive spectral method for unmagnetized relativistic stars presented in Bonazzola, Gourgoulhon & Marck (1998),Gourgoulhon et al. (1999) and part of the lorene C++ class library for numerical relativity. We refer to Bonazzola et al. (1998) and Gourgoulhon et al. (1999) for a detailed description of the numerical method and of the tests carried out. The numerical solution is computed by iteration, starting from crude initial conditions of a parabolic log-enthalpy profile, no rotation, flat space and no magnetic field. Convergence is monitored by computing where is the sum of the absolute values of the log-enthalpy over all collocation points at iteration step . A convergent solution is assumed to be found when the normalized log-enthalpy residual decreases below a prescribed threshold, which is usually chosen to be of the order of or smaller. The computational domain is divided into a spherical nucleus containing the star and two shells covering the exterior, where the last one maps on to a finite interval the whole exterior space, from a certain radius up to spatial infinity. The default number of collocation points used in our study has been in the angular direction and in the radial direction, where the different numbers refer to the nucleus and the two shells, respectively. However, depending on the circumstances, the numerical resolution has been increased when necessary, i.e. up to and in the case of huge models with a matter distribution strongly deviating from a spherical one (cf. Table 3 and Fig. 6).

### 4.2 Comparison with Newtonian results

The results of the numerical scheme in the Newtonian limit have been compared with those of an earlier investigation in the Newtonian regime. Sood & Trehan (1972) have in fact computed linear perturbations of polytropic stars induced by a toroidal magnetic field corresponding to the Newtonian limit of equation (51), namely

 Bϕrsinθ=λ0ρrsinθ, (71)

where denotes the Newtonian mass density. In Table 1, we show normalized values of the resulting oblateness for different adiabatic exponents produced with the present code and the corresponding results from Sood & Trehan (1972). We note that  Sood & Trehan (1972) measured the deformation as where is the radius of the unperturbed model. Values computed according to this definition coincide with within the rounding error of tabulated values. Clearly, for the perturbation parameter , the respective values of agree to within , except for , for which the difference is about . However, for , certain functions intervening in the solution of the perturbed Lane-Emden equation become singular at the surface of the star (Das & Tandon, 1977), suggesting a larger numerical error for the case as computed by Sood & Trehan (1972).

### 4.3 Virial identities

In order to monitor the global error of the numerical models computed with the new code we have calculated the two general-relativistic virial identities GRV2 (Bonazzola, 1973; Bonazzola & Gourgoulhon, 1994) and GRV3 (Gourgoulhon & Bonazzola, 1994). We recall that these identities have to be fulfilled by any solution to the Einstein equations (14)–(17) and are not enforced during the iterative procedure that leads to the solution. The GRV2 identity, in particular, has been shown to be directly related to the global error of a numerical solution (Bonazzola et al., 1993), while the GRV3 identity is an extension of the virial theorem of Newtonian physics to general relativity [see Nozawa et al. (1998) for details of the practical implementation]. For strongly magnetized non-rotating models based on a polytropic EOS with and the moderate number of collocation points specified in Section 4.1, the violations to the identities GRV2 and GRV3 have been found to be of the order of or better. For the huge models listed in Table 3, corresponding values are of the order of which required the number of collocation points to be multiplied by a factor of 8 compared to their default values. In the case of rapid rotation, on the other hand, the numerical grid is not perfectly adapted to the stellar surface, and the spectral approximation suffers from Gibbs phenomena, as reported in Bonazzola et al. (1993). However, the corresponding values of the identities GRV2 and GRV3 are still of the order of . For realistic EOSs, values of GRV2 and GRV3 are of the order of or better.

## 5 Non-rotating Magnetized Models

Non-rotating models of neutron stars with a toroidal magnetic field always generate static space–times for which the time Killing vector is hypersurface-orthogonal. Note that this assumption does not always hold in the case of a poloidal magnetic field because the additional presence of an electric field gives rise to an azimuthal Poynting vector [cf.  Bonazzola et al. (1993)], which introduces a non-zero angular momentum even when the star does not rotate. As a consequence, the shift vector component does not vanish.

Although we are interested in the construction of rapidly rotating models (that we will present in Section 6), non-rotating configurations allow us to investigate the effects of the magnetic field without the influence of rotationally induced effects. As a representative example, Fig. 2 shows isocontours of the magnetic field strength and baryon number density for the non-rotating model PP-Pol2 built with the Pol2 EOS and for the maximum field strength (the physical properties of model PP-Pol2 are listed in Table 4). The unmagnetized reference model with a circumferential radius of is symbolized by a grey disc. In agreement with the simple structure function defined in equation (51), the toroidal magnetic field shown in the left-hand panel vanishes on the axis of symmetry, reaches a maximum value of in the equatorial plane deep inside the star and decreases towards the surface of the star, where it vanishes, so that the magnetic field is fully contained inside the star. The forces exerted by the magnetic field can be visualized through the magnetic potential , which shows a distribution similar to that of the magnetic field strength and that, for this reason, we do not report here. It should be noted that the magnetic potential is repulsive, since it is positive everywhere by construction. This is to be contrasted with the magnetic potential associated with the purely poloidal magnetic field adopted in Bocquet et al. (1995), that was instead attractive. As a result, the corresponding forces are the opposite, despite the isocontours appear to be very similar. Note also that in the present case of a purely toroidal magnetic field defined according to equation (51), the isocontours of coincide with the flow lines of the electric current , which are nested loops in the meridional plane. Finally, we note that vanishes at the surface of the star, implying that the latter coincides with an isosurface of the gravitational potential , as required by equation (53) when and .

The forces that can be derived from vanish on the axis of symmetry and reach their maximum in the equatorial plane. Between the centre of the star and the maximum of , they are directed inwards, inducing an approximately cylindrical compression of the central region of the star, which responds by a prolate deformation visible in the right-hand panel of Fig. 2. In the outer layers of the star, on the other hand, the magnetic forces are directed radially outward, pushing them away from the centre, which results in a growth of the dimensions of the star whose circumferential radius has increased now to a value of .

Both the apparent shape and the matter distribution are prolate, corresponding to negative values of the surface deformation and of the quadrupole distortion, namely and , respectively. The value of is significantly smaller than that of , and this due to the fact that terms of higher order in the multipole expansion of the gravitational potential fall off rapidly while, at the same time, the surface of the star must coincide with an isosurface of . This interpretation is supported by the fact that this difference is significantly smaller for small perturbations around this model, as confirmed by comparing values of respective distortion coefficients and reported in Section 7.

The strength of the magnetic field is controlled by the magnetization parameter (cf. equation 52), and the physical quantities associated with our models can be parametrized accordingly, in particular the maximum magnetic field strength . In Fig. 3, we compare our results with those reported in fig. 6 of Kiuchi & Yoshida (2008),222The data from fig. 6 of Kiuchi & Yoshida (2008) have been read off with the utility g3data see https://github.com/pn2200/g3data. considering in particular the variation of the gravitational mass , of the circumferential radius , of the central baryon number density