Relativistic models of magnetars

# Relativistic models of magnetars: the twisted-torus magnetic field configuration

R. Ciolfi, V. Ferrari, L. Gualtieri, J.A. Pons
Dipartimento di Fisica “G.Marconi”, Sapienza Università di Roma and Sezione INFN ROMA1, 00185 Roma, Italy
Departament de Física Aplicada, Universitat d’Alacant, 03080 Alacant, Spain
###### Abstract

We find general relativistic solutions of equilibrium magnetic field configurations in magnetars, extending previous results of Colaiuda et al. [Colaiuda et al. 2008]. Our method is based on the solution of the relativistic Grad-Shafranov equation, to which Maxwell’s equations can be reduced. We obtain equilibrium solutions with the toroidal magnetic field component confined into a finite region inside the star, and the poloidal component extending to the exterior. These so-called twisted-torus configurations have been found to be the final outcome of dynamical simulations in the framework of Newtonian gravity, and appear to be more stable than other configurations. The solutions include higher order multipoles, which are coupled to the dominant dipolar field. We use arguments of minimal energy to constrain the ratio of the toroidal to the poloidal field.

###### keywords:
stars:neutron, stars:magnetic fields

## 1 Introduction

In this paper we construct models of non rotating, strongly magnetized neutron stars, or magnetars [Duncan & Thompson 1992], in general relativity. We extend our previous work (Colaiuda et al. [Colaiuda et al. 2008]) based on a formalism developed in Konno et al. [Konno et al. 1999] and in Ioka & Sasaki [Ioka & Sasaki 2004], including the toroidal field in a twisted-torus configuration. The extension to this field geometry is accomplished with an appropriate choice of the function which determines, point by point, the ratio between the toroidal and poloidal components of the magnetic field. The non-linear relation among the functions defining the toroidal and poloidal fields, naturally leads to couplings between different multipoles, thus making inadequate the one multipole solution which is usually assumed.

A motivation to find consistent equilibrium solutions in general relativity with this particular geometry comes from recent progresses in numerical magneto-hydrodynamic (MHD) simulations, that have made possible to study the dynamics and the stability of magnetic stars. By following the time evolution of generic initial configurations, in the framework of Newtonian gravity and using polytropic equations of state, Braithwaite & Spruit [Braithwaite & Spruit 2004] and Braithwaite & Nordlund [Braithwaite & Nordlund 2006] (see also Braithwaite & Spruit [Braithwaite & Spruit 2006]) have found magnetic field configurations, which are stable on timescales much longer than the Alfvèn time: they decay only due to finite resistivity. These configurations are roughly axisymmetric; the poloidal field extends throughout the entire star and to the exterior, while the toroidal field is confined in a torus-shaped region inside the star, where the field lines are closed. These configurations were named twisted-torus. Furthermore, Yoshida et al. [Yoshida et al. 2006] have shown that such configurations are not significantly affected by rotation; Geppert & Rheinhardt [Geppert & Rheinhardt 2006] studied the dependence of magnetostatic equilibrium configurations on the rotational velocity and on the initial angle between rotation and magnetic axis, finding hints for the existence of a unique stable dipolar magnetostatic configuration, independent of the initial field geometry.

We must remark that this particular field geometry resulting from dynamical simulations is obtained assuming that outside the star there is vacuum; consequently, outside the star electric currents are forbidden and the magnetic field can only be poloidal. This implies that the toroidal field cannot extend to the exterior and that the field lines which cross the surface are purely poloidal, whereas the field lines confined inside the star can maintain a mixed (poloidal and toroidal) structure. This configurations appear to be stable on dynamical timescales, probably due to magnetic helicity conservation, which requires the persistence of a toroidal component of the field. Notice, however, that different solutions including a magnetosphere may be possible. In this case, the toroidal field could also extend to the external region leading to a twisted magnetosphere [Lyutikov 2006, Pavan et al. 2009].

In our perturbative approach, there is a free parameter which represents the ratio between the toroidal and the poloidal components of the magnetic field. We estimate the value of this parameter, by identifying the configuration of minimal energy at fixed magnetic helicity. We also mention that in our configurations the external field has mainly a dipole structure, with small corrections from higher multipoles. Furthermore, confirming a previous suggestion by Prendergast [Prendergast1956], the toroidal and poloidal fields have amplitudes of the same order of magnitude, whereas the energy associated to the toroidal field is an order of magnitude smaller than that of the poloidal field, since the former is confined in a relatively small region. Similar configurations have been found in Newtonian models including rotation, in Yoshida & Eriguchi [Yoshida & Eriguchi 2006], Yoshida et al. [Yoshida et al. 2006].

In this paper we consider non-rotating stars because observed magnetars have a very slow rotation rate, although high rotation rates may occur in the early stages of their evolution.

The structure of the paper is the following. The model is presented in Section 2. In Section 3 we discuss a configuration with purely dipolar magnetic field, neglecting the couplings with higher multipoles. In Section 4 we include the and field components; Section 5 accounts for the general case, including all multipoles and their couplings. In Section 6 we compute the total energy and the magnetic helicity, and estimate the parameter which determines the ratio between toroidal and poloidal fields by energy minimization; we also compute the magnetic energy, and compare the contributions of the poloidal and toroidal fields for different values of . In Section 7 we discuss the results and draw the conclusions.

## 2 Basic equations and formalism

We assume that the (non-rotating) magnetized star is stationary and axisymmetric. We further assume that the magnetic field acts as a perturbation of a spherically symmetric background describing a spherical star. The magnetized fluid is described within the framework of ideal MHD, in which the effects of finite electrical conductivity are neglected. Rigorously speaking, this approximation is only valid while the crust is still completely liquid and while the core matter has not yet performed the phase transition to the superfluid state, which is expected to occur at most a few hours after birth (see e.g. section 5.1 in [Aguilera et al. 2008] and references therein). The onset of superfluidity and/or crystallization limits the period during which magnetostatic equilibrium can be established. Both the melting temperature and the critical temperature of transition to the superfluid state, are between to K, and a typical neutron star quickly cools down below this temperature in a few hours. However, since the characteristic Alfvén time is on the order of s, depending on the background field strength, there is ample time for the magnetized fluid to reach a stable state, as shown in Braithwaite & Spruit (2006), while the state of matter is still liquid. After the crust is formed, the magnetic field is frozen in, and it only evolves on a much longer timescale due to Ohmic dissipation or, in some case, due to the Hall drift [Pons & Geppert 2007]. Therefore, it is reasonable to expect that the MHD equilibrium configurations set within the first day after formation, will fix the magnetic field geometry for a long time.

Here we first summarize the basic equations of ideal MHD in the framework of general relativity and then introduce the perturbative approach. Next, we obtain the form of the electromagnetic potential in the case of twisted-torus configurations, and derive the relativistic Grad-Shafranov equation. We use spherical coordinates, , where . A stationary axisymmetric space-time admits two killing vectors, and , and with our coordinate choice all quantities (including the components of the metric tensor ) are independent of and .

### 2.1 Equations of ideal MHD in General Relativity

According to a comoving observer with four-velocity , the stress-energy tensor of a perfect fluid with an electromagnetic field is

 Tμν=Tμνfluid+Tμνem, (1)

where

 Tμνfluid=(ρ+P)uμuν+Pgμν, (2)
 Tμνem=14π[(uμuν+12gμν)B2−BμBν]. (3)

As usual, Euler’s equations are found by projecting the equation orthogonally to

 (ρ+P)aμ+P,μ+uμuνP,ν−fμ=0, (4)

where is the Lorentz force and . Here, is the Maxwell tensor, in terms of which the electric and magnetic fields can be defined as

 Eμ≡Fμνuν,     Bα≡12ϵαβγδuβFγδ. (5)

The basic equations of ideal, general relativistic MHD are, then: (i) the continuity equation ; (ii) Maxwell’s equations ; (iii) the condition of a vanishing electric field in the comoving frame and (iv) Euler’s equations (4).

### 2.2 The perturbative approach and the form of the electromagnetic potential

We treat the magnetic field as an axisymmetric perturbation of a spherically symmetric background and seek for stationary solutions. The background metric is

 ds2=−eν(r)dt2+eλ(r)dr2+r2(dθ2+sin2θdϕ2), (6)

where , are solution of the unperturbed Einstein equations (the TOV equations) for assigned equations of state. The unperturbed 4-velocity is . To model the unperturbed neutron star we use the equation of state of Akmal, Pandharipande and Ravenhall called APR2 [Akmal at al. 1998], with a standard equation of state for the stellar crust (see Benhar et al. [Benhar, Ferrari & Gualtieri 2004]), which results in a neutron star of mass and radius km. We remark that our EOS accounts for the density-pressure relation in the crustal region, but not for its elastic properties. Our equations apply to a star where the solid crust has not formed yet, or to configurations with a relaxed crust where elasticity is irrelevant.

It can be shown (see for instance Colaiuda et al. [Colaiuda et al. 2008]) that are of order , and the perturbations are of order . Therefore, at first order in the magnetic field is coupled only to the unperturbed background metric (6), whereas the deformation of the stellar structure induced by the magnetic field, which we do not consider in this paper, appears at order . Furthermore, and . Note that in the Grad-Shafranov equation, which we solve to order , the metric perturbations do not appear; thus, to find the magnetic field configurations we do not need to solve Einstein’s equations. In Section 6 and in Appendix B, we will solve some components of Einstein’s equations, in order to evaluate the total energy of the system.

With these assumptions, the potential , at , has the form . With an appropriate gauge choice we can impose and write the potential as

 Aμ=(0,eλ−ν2Σ,0,ψ), (7)

where the components of are expressed in terms of two unknown functions, and .

A further simplification of is possible by exploiting the fact that . Using Maxwell’s equations and neglecting higher order terms, we find

 ~ψ,θψ,r=~ψ,rψ,θ, (8)

where . This result implies and allows us to write

 sinθΣ,θ=ζ(ψ)ψ, (9)

where is a function of of order .

The function represents the ratio between the toroidal and poloidal components of the magnetic field; different choices of this function lead to qualitatively different field configurations. The simplest case is [Ioka & Sasaki 2004, Colaiuda et al. 2008, Haskell et al. 2008], but with this choice (like with other simple choices of ) the field lines of the toroidal field reach the exterior of the star, where there is vacuum. However, the magnetic field in vacuum can only be poloidal (see, for instance, Colaiuda et al. [Colaiuda et al. 2008]), thus this solution presents an inconsistency. To avoid this inconsistency, one should consider a non-vacuum exterior, i.e. a magnetosphere, but modelling a neutron star magnetosphere is a quite difficult task, especially in general relativity. An alternative choice is to assume that the magnetic field is entirely confined inside the star [Ioka & Sasaki 2004, Haskell et al. 2008], but in this way the parameter must assume particular values; or, one can instead accept that the toroidal field has a discontinuity at the stellar surface, vanishing in the exterior [Colaiuda et al. 2008]; in this way the entire range of can be studied, but the discontinuity in the field will, for consistence, imply the existence of surface currents.

A different choice is made in this paper, where we assume that the toroidal field is confined in a toroidal region inside the neutron star, such that its field lines never cross the stellar surface, as in the twisted-torus configuration. As mentioned in Section 1, Newtonian numerical simulations [Braithwaite & Spruit 2004, Braithwaite & Nordlund 2006, Braithwaite & Spruit 2006] suggest that these configurations are indeed a quite generic outcome of the evolution of strongly magnetized stars.

The twisted-torus configuration can be obtained by choosing the following form of the function

 ζ(ψ)=ζ0[|ψ/¯ψ|−1] Θ(|ψ/¯ψ|−1). (10)

A similar choice has been made, in a Newtonian framework, in Yoshida et al. [Yoshida et al. 2006]. In equation (10), is a constant of order ; is a constant of order : it is the value of at the boundary of the toroidal region where the toroidal field is confined (this boundary is tangent to the stellar surface); finally, is the usual Heaviside function. With this choice, the function vanishes at the stellar surface, where , and the magnetic field

 Bμ =( 0,e−λ2r2sinθψ,θ,−e−λ2r2sinθψ,r, (11) −e−ν2ζ0ψ(|ψ/¯ψ|−1)r2sin2θΘ(|ψ/¯ψ|−1))

has no discontinuities.

### 2.3 The relativistic Grad-Shafranov equation

The Grad-Shafranov (GS) equation, which allows to determine the magnetic field configuration, can be derived from the -component of Maxwell’s equations

 Jϕ=−e−λ4π[ψ,rr+ν,r−λ,r2ψ,r]−14πr2[ψ,θθ−cotθψ,θ] (12)

and from the -components of Euler’s equations (4), as follows. Euler’s equations give

 fa = (ρ+P)aa+P,a+uauνP,ν (13) = (ρ+P)(ν2−eν2δut),a+P,a+O(B4).

For barotropic equations of state , the first principle of thermodynamics allows to write

 P,a=(ρ+P)(lnρ+Pn),a, (14)

then (13) yields

 fa=(ρ+P)χ,a, (15)

where . On the other hand, the -components of the Lorentz force can be written as (see Colaiuda et al. [Colaiuda et al. 2008])

 fa=ψ,ar2sin2θ~Jϕ, (16)

where, in the present case,

 ~Jϕ=Jϕ−e−νζ204π[ψ−3ψ|ψ/¯ψ|+2ψ3/¯ψ2] Θ(|ψ/¯ψ|−1).

Therefore,

 ψ,ar2sin2θ~Jϕ=(ρ+P)χ,a. (17)

From it follows that

 ψ,r(~Jϕ(ρ+P)r2sin2θ),θ−ψ,θ(~Jϕ(ρ+P)r2sin2θ),r=0,

which implies

 (~Jϕ(ρ+P)r2sin2θ)=F(ψ)=c0+c1ψ+O(B2), (18)

with , constants of order , respectively. Hence, turns out to be

 Jϕ = e−νζ204π[ψ−3ψ|ψ/¯ψ|+2ψ3/¯ψ2] Θ(|ψ/¯ψ|−1) (19) +(ρ+P)r2sin2θ[c0+c1ψ].

From Eqns. (12), (19) the relativistic GS equation at first order in is finally obtained:

 −e−λ4π[ψ,rr+ν,r−λ,r2ψ,r]−14πr2[ψ,θθ−cotθψ,θ] −e−νζ204π[ψ−3ψ|ψ/¯ψ|+2ψ3/¯ψ2] Θ(|ψ/¯ψ|−1) =(ρ+P)r2sin2θ[c0+c1ψ]. (20)

If we now define and expand the function in Legendre polynomials

 a(r,θ)=∞∑l=1al(r)Pl(cosθ); (21)

the GS equation rewrites as

 −sinθ4π∞∑l=1Pl,θ(e−λa′′l+e−λν′−λ′2a′l−l(l+1)r2al) −e−ν4πΘ(∣∣ ∣∣1¯ψ∞∑l=1alPl,θsinθ∣∣ ∣∣−1)ζ20[∞∑l=1alPl,θsinθ −3|¯ψ|⎛⎝∞∑l,l′=1alPl,θsinθ|al′Pl′,θsinθ|⎞⎠ +2¯ψ2⎛⎝∞∑l,l′,l′′=1alal′al′′Pl,θPl′,θPl′′,θsin3θ⎞⎠] =(ρ+P)r2sin2θ[c0+c1∞∑l=1alPl,θsinθ]. (22)

Here and in the following we denote with primes the differentiation with respect to .

Finally, projecting Eq. (22) onto the different harmonic components, we obtain a system of coupled ordinary differential equations for the functions . The projection is performed using the property

 2l′+12l′(l′+1)∫π0Pl,θPl′,θsinθdθ=δll′. (23)

If we consider the contribution of different harmonics, we need to solve a system of coupled equations, obtained from (22), for the functions .

### 2.4 Boundary conditions

The functions must have a regular behaviour at the origin; by taking the limit of the GS equation one can find

 al(r→0)=αlrl+1, (24)

where are arbitrary constants to be fixed.

Outside the star, where there is vacuum and the field is purely poloidal, equations (22) decouple, and can be solved analytically for each value of . The solution can be expressed in terms of the generalized hypergeometric functions (), also known as Barnes’ extended hypergeometric functions, as follows:

 al = A1 r−l F([l,l+2],[2+2l],z) (25) +A2 rl+1 F([1−l,−1−l],[−2l],z) .

where and and are arbitrary integration constants, which must be fixed according to the values of the magnetic multipole moments. Regularity of the external solution at requires for all multipoles. For example, for we have

 a1 ∝ r2[ln(1−z)+z+z22] a2 ∝ r3[(4−3z)ln(1−z)+4z−z2−z36] a3 ∝ r4[(15−20z+6z2)ln(1−z)+15z (26) −25z22+z3+z412].

At the stellar surface we require the field to be continuous. This condition is satisfied if and are continuous. For practical purposes, the boundary conditions at can be written as

 a′l=−lRflal (27)

where is a relativistic factor which only depends on the star compactness (in the Newtonian limit all ), and can be numerically evaluated with the help of any algebraic manipulator. For our model (), the values of for the first five multipoles are 1.338, 1.339, 1.315, 1.301, and 1.292, respectively.

In general, there are arbitrary constants to be fixed: the constants , associated to the condition (24), and . Thus, we need to impose constraints, of which are determined by the boundary conditions. conditions are provided by Eq. (27), i.e. by imposing continuity in of the ratios . The overall normalization of the field gives another condition, which is fixed by imposing that the value of the contribution at the pole is G (this corresponds to set km). The reason for this choice is that the surface value of the magnetic field is usually inferred from observations by applying the spin-down formula, and assuming a purely dipolar external field; for magnetars, the value of estimated in this way is G. The remaining condition will be imposed as follows.

In the case of a purely dipolar field (), we shall assume . In the general multipolar case, we choose to impose that the external contribution of all the harmonics, i.e. , is minimum.

## 3 The case of purely dipolar field

We begin discussing the simplest case of a purely dipolar configuration, in which all couplings with higher order multipoles are neglected in Eq. (22) (). In this case, for any assigned value of there exists an infinite set of solutions, each corresponding to a value of ; these solutions describe qualitatively similar magnetic field configurations.

However, when higher order harmonics are taken into account, as we will see in the next Section, the picture changes. For instance, when and the harmonic components are included, the equations for and decouple: the equation for is the same as in the purely dipolar case, but a solution for satisfying the appropriate boundary conditions exists only for a unique value of . Therefore, in the general case is not a truly free parameter (this is true also for ), and the fact that in the purely dipolar case it looks as such, is an artifact of the truncation of the multipoles. In order to provide a mathematically simple example, which will be useful to understand the structure of the twisted-torus configurations, in this Section we shall consider the simplest case .

By projecting Eq. (22) onto the harmonic, and neglecting all contributions from terms we find

 14π(e−λa′′1+e−λν′−λ′2a′1−2r2a1) −e−ν4π∫π0(3/4)Θ(∣∣∣−a1sin2θ¯ψ∣∣∣−1) =(3/4)∫π0c0(ρ+P)r2sin3θdθ=c0(ρ+P)r2. (28)

The tetrad components of the magnetic field (i.e. the components measured by a locally inertial observer) are:

 B(r) = ψ,θr2sinθ, B(θ) = −e−λ2rsinθψ,r, B(ϕ) = −e−ν2ζ0ψ(|ψ/¯ψ|−1)rsinθ⋅Θ(|ψ/¯ψ|−1), (29)

where .

The profiles of the tetrad components of the field inside the star, are plotted in Fig. 1 for increasing values of ; is evaluated in and , in . In Fig. 2 we show the projection of the field lines in the meridional plane, for km. Figs. 1 and 2 show that the toroidal field is confined within a torus-shaped region; its amplitude ranges from zero, at the border of the region, to a maximum, close to its center. At the stellar surface and in the exterior vanishes, showing that there is no discontinuity in the toroidal field. The panels of Fig. 1 show the field profiles for different values of : larger values of correspond to a toroidal field with increasing amplitude, confined in an increasingly narrow region close to the stellar surface, while the amplitude of the poloidal components (, ) decreases. We remark that this implies that inside the star we cannot have a twisted-torus configuration where the toroidal component dominates with respect to the poloidal one: if becomes larger with respect to and , the domain where it is non vanishing shrinks.

## 4 The case with l=1 and l=2 multipoles

We now proceed with our investigation considering the and contributions, and setting . The projection of the GS equation (22) onto the harmonics and gives the following coupled equations

 14π(e−λa′′1+e−λν′−λ′2a′1−2r2a1) −e−ν4π∫π0(3/4)Θ(∣∣∣−a1−3a2cosθ¯ψsin2θ∣∣∣−1) ⋅ζ20[−a1−3a2cosθ+3(a1+3a2cosθ) ⋅∣∣∣−a1−3a2cosθ¯ψsin2θ∣∣∣+2sin4θ(−a31−9a21a2cosθ =(ρ+P)r2(c0−45c1a1), (30)
 14π(e−λa′′2+e−λν′−λ′2a′2−6r2a2) +e−ν4π∫π0(5/12)Θ(∣∣∣−a1−3a2cosθ¯ψsin2θ∣∣∣−1) ⋅ζ20[−a1−3a2cosθ+3(a1+3a2cosθ) ⋅∣∣∣−a1−3a2cosθ¯ψsin2θ∣∣∣+2sin4θ(−a31−9a21a2cosθ =−47(ρ+P)r2c1a2. (31)

We integrate this system by imposing the boundary conditions discussed above, i.e. a regular behaviour at the origin (equation (24)) and continuity at the surface of with the analytical external solutions given by (26).

Let us first consider the simple case . Eqns. (30), (31) decouple, and become

 e−λa′′1+e−λν′−λ′2a′1−2r2a1 =4π(ρ+p)r2[c0−45c1a1] e−λa′′2+e−λν′−λ′2a′2−6r2a2 =−16π7(ρ+p)r2c1a2. (32)

There are four constants to fix () and three conditions: km (normalization) and the ratios and from the matching with the exterior solutions; thus, we need an additional requirement. We remark that we cannot impose as in the purely dipolar case, because the ratio depends only on , and the matching with the exterior solution is possible only for a particular value of , i.e. km.

If we impose that is minimum, we find that this condition yields the trivial solution (with non-vanishing ). Indeed, from Eqns. (32) it is straightforward to see that is a solution of the system. When , equations (30), (31) are coupled, but they still allow the trivial solution , which minimizes , with non-vanishing . The existence of this solution is a remarkable property of this system, and it is due to the fact that the integral in on the left-hand side of Eq. (31) vanishes for (the integrand becomes odd for parity transformations ). Hence, if we look for a solution which minimizes the contributions from the components at the stellar surface, we have to choose the trivial solution .

If, instead, we do not require that is minimum, and we assign a finite value to the ratio , we find a non-trivial field configuration which is non symmetric with respect to the equatorial plane. This feature is shown in Fig. 3, where the projection of the field lines in the meridional plane are plotted for and equal to , and , respectively.

## 5 The general case

When all harmonics are taken into account, there exist two distinct classes of solutions: those symmetric (with respect to the equatorial plane), with vanishing even order components (), and the antisymmetric solutions, with vanishing odd order components (). Both solutions satisfy the GS equation (22). Let us consider the symmetric class. When the integrals arising when equation (22) is projected onto the even harmonics, which couple odd and even terms, vanish since the integrands change sign under parity transformations. Therefore, the symmetric solutions can be found by setting , projecting Eq. (22) onto the odd harmonics and solving the resulting equations for . Similarly, the integrals in equation (22) projected onto the odd harmonics, vanish when ; thus, we can consistently set , and find the antisymmetric solutions using the same procedure.

In Section 4, we set the value of at the surface to be km and we minimized the contribution. It is clear that, since the and multipoles belong to different families, any attempt to minimize the relative weight of one with respect to the other leads to the trivial solution. The properties of equation (22) discussed above, tell us that if we cannot consistently set to zero the remaining odd order components . However, we have the freedom of setting to zero all even terms . Therefore, since we have chosen to minimize the contributions of the harmonics outside the star, we shall focus on the symmetric family of solutions (); we will briefly discuss an example belonging to the antisymmetric family in subsection 5.4.

### 5.1 The case with l=1 and l=3

We now consider the system of equations including only the and components. The projected system is

 14π(e−λa′′1+e−λν′−λ′2a′1−2r2a1) −e−ν4π∫π0(3/4)ζ20(ψ−3ψ|ψ/¯ψ|+2ψ3/¯ψ2) ⋅Θ(|ψ/¯ψ|−1)sinθdθ =[c0−45c1(a1−37a3)](ρ+P)r2, (33)
 14π(e−λa′′3+e−λν′−λ′2a′3−12r2a3) +e−ν4π∫π0(7/48)ζ20(ψ−3ψ|ψ/¯ψ|+2ψ3/¯ψ2) ⋅Θ(|ψ/¯ψ|−1)(3−15cos2θ)sinθdθ =215c1(ρ+P)r2(a1−4a3), (34)

where

 ψ=[−a1+a3(3−15cos2θ)2]sin2θ. (35)

We again impose regularity at the origin (Eq. (24)), continuity in of with the vacuum solutions for , given by Eq. (26), and we fix km by normalization. For the remaining constraint we choose the solution that minimizes the absolute value of . We find that there is a discrete series of local minima of , and we select among them the absolute minimum.

Fig. 4 shows the profiles of the tetrad field components (see Eq. (29)) obtained by numerically integrating Eqns. (33), (34), for different values of . is evaluated at , while , are evaluated at . As increases, the magnitude of the toroidal field becomes larger, but the region where it is confined shrinks, as already found in section 3. The projection of the field lines in the meridional plane is shown in Fig. 5 for the same values of . It shows that, for km, the magnetic field lines lie in disconnected regions, separated by dashed lines in the figure. Inside these regions, the function has opposite sign and no toroidal field is present. A similar phenomenon has been discussed in Colaiuda et al. [Colaiuda et al. 2008]. As we will see in the next Section, the occurrence of these regions is an artifact of the truncation in the harmonic expansion, and disappears as higher order harmonics are included.

For completeness we also mention that the solutions corresponding to the local minima of different from the absolute minimum, correspond to very peculiar field configurations (see Fig. 6). The function has nodes on the equatorial plane, therefore the field lines lie in disconnected regions; for a fixed value of , the number of nodes increases as increases. These peculiar solutions exist for any value of , and appear also when higher order harmonic components are considered. Thus, they are not artifacts of the -truncation.

### 5.2 The case with l=1,3,5

We now include the contribution. The three equations obtained by projecting the GS equation (22) onto are given in the Appendix A. The boundary conditions are essentially the same as in the previous Section; in particular, we look for the absolute minimum of , with fixed km.

In Fig. 7 the profiles of the tetrad components of the magnetic field are plotted for values of in the range km. Fig. 8 shows the projections of the field lines in the meridional plane corresponding to the same values of . Comparing the results with the case we see that the presence of the harmonic modifies the magnetic field shape, but most of the features discussed in the previous Section are still present.

An interesting difference is the following. While in the case for km we find field configurations which exhibit two disconnected regions where the function has opposite sign and the magnetic field lines are confined (regions within dashed lines in Fig. 5), this does not occur when the component is taken into account. This shows that the above feature has to be considered as an artifact of the truncation in the harmonic expansion.

### 5.3 Higher order multipoles

Up to now we have included components with , neglecting the contribution from . In order to test the accuracy of this approximation, we have studied the convergence of the harmonic expansion. To this purpose, we have solved the GS equation (22) including odd harmonic components up to , for and km, and we have computed the quantities

 Δ(5)(r,θ) = ∣∣∣ψl≤5(r,θ)−ψl≤3(r,θ)¯ψ∣∣∣, Δ(7)(r,θ) = ∣∣∣ψl≤7(r,θ)−ψl≤5(r,θ)¯ψ∣∣∣. (36)

These functions are shown in Fig. 9. They are plotted only inside the star since outside they are much smaller. Fig. 9 shows that the error in neglecting , quantified by the function , is for and for km. Furthermore, a comparison of and shows that the harmonic expansion converges.

### 5.4 An example of antisymmetric solution

Here we show an example of a solution belonging to the antisymmetric family corresponding to . In Fig. 10 we plot the field lines projected on the meridional plane, for km and km. We remark that the field lines are antisymmetric with respect to the equatorial plane; as a consequence, the total magnetic helicity is zero (see Section 6). Similar zero-helicity configurations have been considered in Braithwaite [Braithwaite 2008a].

## 6 Magnetic helicity and energy

The stationary configurations of magnetized neutron stars which we have found, depend on the value of the free parameter , i.e. on the ratio between the toroidal and the poloidal components of the magnetic field. In this Section, we provide an argument to assign a value to . Furthermore, we compute the magnetic energy of the system to compare the contributions from poloidal and toroidal fields.

The total energy of the system (the star, the magnetic field and the gravitational field) can be determined by looking at the far field limit () of the spacetime metric [Misner et al. 1973, Thorne 1980]. Following Colaiuda et al. [Colaiuda et al. 2008], Ioka & Sasaki [Ioka & Sasaki 2004], we write the perturbed metric as

 ds2 = −eν[1+2h(r,θ)]dt2+eλ[1+2eλrm(r,θ)]dr2 (37) +r2[1+2k(r,θ)](dθ2+sin2θdϕ2) +2i(r,θ)dtdr+2v(r,θ)dtdϕ+2w(r,θ)drdϕ

where, in particular, . The total mass-energy of the system is

 E=M+δM (38)

where is the gravitational mass of the unperturbed star and

 δM=limr→∞m0(r). (39)

In Appendix B, we discuss the equations which allow to determine . We remark that includes different contributions, due to magnetic energy, deformation energy, etc.

In order to evaluate the magnetic contribution to , it is convenient to use the Komar-Tolman formula (see for instance Straumann [Straumann2004], Chap. 4) for the total energy:

 E=2∫V(Tμν−