On magnetic equilibria in barotropic stars

# On magnetic equilibria in barotropic stars

Cristóbal Armaza11affiliation: Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile. , Andreas Reisenegger11affiliation: Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile. , and Juan Alejandro Valdivia22affiliation: Departamento de Física, Facultad de Ciencias, Universidad de Chile, Casilla 653, Santiago, Chile.
###### Abstract

Upper main sequence stars, white dwarfs and neutron stars are known to possess stable, large-scale magnetic fields. Numerical works have confirmed that stable MHD equilibria can exist in non-barotropic, stably stratified stars. On the other hand, it is unclear whether stable equilibria are possible in barotropic stars, although the existing evidence suggests that they are all unstable. This work aims to construct barotropic equilibria in order to study their properties, as a first step to test their stability. We have assumed that the star is a perfectly conducting, axially symmetric fluid, allowing for both poloidal and toroidal components of the magnetic field. In addition, we made the astrophysically justified assumption that the magnetic force has a negligible influence on the fluid structure, in which case the equilibrium is governed by the Grad–Shafranov equation, involving two arbitrary functions of the poloidal flux. We built a numerical code to solve this equation, allowing for an arbitrary prescription for these functions. Taking particularly simple, but physically reasonable choices for these functions with a couple of adjustable parameters, all of the equilibria found present only a small () fraction of the magnetic energy stored in the toroidal component, confirming previous results. We developed an analytical model in order to study in more detail the behavior of the magnetic energy over the full range of parameters. The model confirms that the toroidal fraction of the energy and the ratio of toroidal to poloidal flux are bounded from above for the whole range of parameters.

magnetohydrodynamics (MHD); stars: magnetic field; stars: neutron; stars: white dwarfs

## 1. Introduction

In contrast to the magnetic fields of low-mass stars (e.g. the Sun), magnetic fields in upper main sequence stars, white dwarfs, and neutron stars are quite stable on time scales comparable to their life times. However, not all of these stars are actually “magnetic,” in the sense that massive stars and white dwarfs follow a bimodal distribution in which only a fraction within % of each star type stands out for having strong magnetic fields when compared to the remaining (Liebert et al., 2003; Ferrario & Wickramasinghe, 2006; Aurière et al., 2007; Kawka et al., 2007; Donati & Landstreet, 2009; Wade et al., 2014). Throughout this work, we will be interested in the stable fields in upper main sequence stars, white dwarfs, and neutron stars.

From a theoretical point of view, the persistence of magnetic fields in massive stars and compact remnants motivates the interest in what physical conditions are involved in sustaining such equilibria. Both analytical and numerical studies have generally pointed out that stable equilibrium configurations require a poloidal component of the magnetic field to stabilize a toroidal one, and vice versa (Prendergast, 1956; Braithwaite & Spruit, 2004; Braithwaite & Nordlund, 2006; Akgün et al., 2013), while purely poloidal and purely toroidal equilibria undergo intrinsic instabilities related to their geometries (Markey & Tayler, 1973; Tayler, 1973; Wright, 1973; Flowers & Ruderman, 1977; Kiuchi et al., 2008; Ciolfi et al., 2011; Lasky et al., 2011; Marchant et al., 2011). The simulations of Braithwaite and collaborators showed that initially random fields often evolve naturally into nearly axisymmetric, toroidal-poloidal “twisted-torus” configurations on short (Alfvén) time scales.

In addition, an important ingredient to consider is that the matter within the stars considered is non-barotropic, that is, the pressure is a function of mass density and another thermodynamic quantity, which is conserved in a fluid element if the latter is perturbed from its equilibrium position (Reisenegger, 2009). This second thermodynamic variable accounts for stable stratification, i.e., stability against convective motion. On short timescales, massive stars and white dwarfs are stratified by entropy gradients, while in neutron stars this role is played by a varying chemical composition. Stable stratification could be a crucial ingredient in stabilizing these equilibria, inhibiting radial motions and hence suppressing fluid instabilities occurring within these stars, as supported by recent numerical simulations (Mitchell et al., 2014). As a matter of fact, it has recently been suggested that a plausible condition of dynamical stability for these magnetized stars can be expressed as

 2a(Emag/Egrav)

where is the magnetic-to-gravitational energy ratio, is the fraction of magnetic energy stored in the poloidal component, and is a pure number accounting for how stratified the star is (Braithwaite, 2009; Akgün et al., 2013). The additional fact that all of these objects are expected to contain a gravitational energy much larger than the magnetic energy (Reisenegger, 2009) would imply that the toroidal component can be stabilized by a much weaker poloidal field, but the converse is not true. This is of particular relevance in the context of energy released by magnetars, where it has been claimed that a much stronger internal toroidal field would account for observed outbursts emitted from these objects (Thompson & Duncan, 1996).

Stable stratification, although crucial on short time scales, could be eroded if some processes change either entropy or chemical composition. In massive stars and white dwarfs, heat conduction is too slow to significantly change entropy during their life time. For neutron stars, the stable stratification could be overcome by changing the chemical composition by means of -decay and/or ambipolar difusion (Goldreich & Reisenegger, 1992; Pethick, 1992; Reisenegger, 2009), making them effectively barotropic over long time scales.

On the other hand, and despite their lack of realism, barotropic equations of state have often been assumed to describe the matter within upper main sequence stars, white dwarfs and neutron stars (Yoshida & Eriguchi, 2006; Haskell et al., 2008; Lander & Jones, 2009; Ciolfi et al., 2009; Fujisawa et al., 2012; Pili et al., 2014; Bucciantini et al., 2015), because of the expected simplicity as a very first approximation to the realistic state of matter. Yet, it is important to identify correctly when one is actually modeling a barotropic fluid. The local adiabatic index of the fluid, , where is the specific entropy (or another quantity conserved in fast fluid displacements), will not generally be equal to the analogous index describing the equilibrium, . If , we say that the fluid is barotropic. A single-species, non-interacting degenerate Fermi gas, for example, does possess a barotropic equation of state, which in limiting cases of nonrelativistic and ultrarelativistic particles reduces to polytropes with and , respectively.

Barotropy strongly restricts the allowed range of possible equilibrium configurations and, as discussed, does not strictly represent the realistic stably stratified matter within these objects on short (Alfvén) time scales. Moreover, the question whether magnetic equilibria in barotropic stars can be stable or not is being answered negatively by recent studies (Lander & Jones, 2012; Mitchell et al., 2014).

Several authors (Ciolfi et al., 2009; Lander & Jones, 2009; Fujisawa et al., 2012; Gourgouliatos et al., 2013; Pili et al., 2014; Bucciantini et al., 2015) have explored the possible axially symmetric equilibria in barotropic stars, generally finding that the fraction of the total magnetic energy corresponding to the toroidal component, , is at most a few %, even in cases of comparable poloidal and toroidal magnetic field strength. Since stable stratification, which is absent in the barotropic case, is expected to be a key piece in the stability of the stars considered here, and given such a small fraction of toroidal energy apparently inherent to barotropic configurations, it is likely that all of these equilibria are dynamically unstable, as supported by recent simulations (Lander & Jones, 2012; Mitchell et al., 2014). More recent works (Ciolfi & Rezzolla, 2013; Fujisawa & Eriguchi, 2013) have shown, however, that higher fractions are actually possible, making a more extensive survey of these equilibria relevant. Studying properties of barotropic equilibria could also be relevant considering the scenario in which neutron stars would reach an effectively barotropic state after overcoming stable stratification by means of dissipative processes, as already discussed.

The present work is focused on obtaining a wide range of barotropic equilibria, paying attention to their main properties. In addition, these results may be considered as a starting point to study in more detail whether magnetic fields in barotropic stars can be stable or not. This paper is organized as follows. §2 presents the formalism used to construct barotropic equilibria, in order to introduce the notation used throughout this paper. In §3 we solve numerically the resulting equilibrium equation. We summarize the tests carried out to check our code and present the main results obtained. §4 expands the analysis of barotropic equilibria, introducing an asymptotic, analytical model to explore equilibrium configurations beyond numerical simulations. Finally, §5 presents the conclusions of this paper.

## 2. Magnetic equilibria in barotropic stars

We assumed that the star is a perfectly conducting fluid. In the ideal MHD approximation, the dynamical equilibrium of a star is described by the Euler equation,

 0=−∇P−ρ∇Ψ+1cJ×B, (2)

where , , , , and are the fluid pressure, mass density, gravitational potential, magnetic field, electric current density and speed of light, respectively, so that the last term in Eq. (2) is the Lorentz force per unit volume. Throughout this work, a spherical coordinate system with origin at the center of the star is used to describe all quantities. Also, we assumed that the current density drops toward the stellar surface and vanishes outside, as expected since the mass density and the charged-particle density do so. For simplicity, we considered axially symmetric magnetic stars. The magnetic field may then be expressed as the sum of a poloidal and a toroidal component, each determined by a single scalar function (Chandrasekhar & Prendergast, 1956),

 B=Bpol+Btor=∇α(r,θ)×∇ϕ+β(r,θ)∇ϕ, (3)

where , being the orthonormal basis of the coordinate system, with . By Ampère’s law, the electric current density reads

 4πcJ=∇β×∇ϕ−Δ∗α∇ϕ, (4)

 Δ∗ ≡r2sin2θ∇⋅[1r2sin2θ∇] =∂2∂r2+sinθr2∂∂θ(1sinθ∂∂θ) (5)

was introduced. For an axisymmetric equilibrium, the azimuthal component of the magnetic force must vanish, which is imposed by demanding that be a function of only, . This implies that both and are constant along poloidal field lines. Also, notice that, since the poloidal component is perpendicular to the toroidal one by construction, one can always split the magnetic energy in terms of the energy stored in each component,

 Emag=Epol+Etor=∫B2pol8πdV+∫B2tor8πdV, (6)

where the integration is carried out over all space.

### 2.1. Magnetic field outside the star.

Since we demanded vacuum outside the star, has to be a constant there (see Eq. (4)). This constant must be zero in order to have a finite magnetic field along the axis. Also, it is needed that

 Δ∗α=0, (7)

which is a linear, homogeneous partial differential equation, whose general solution after separating variables is

 α(r,θ)=∞∑ℓ=1sinθ[aℓrℓ+1+bℓr−ℓ]P1ℓ(cosθ), (8)

where is the associated Legendre polynomial of order with azimuthal index . Since the magnetic field must vanish far from the star, for all values of , and

 α =∞∑ℓ=1bℓsinθrℓP1ℓ(cosθ), (9) Br =1r2sinθ∂α∂θ=−∞∑ℓ=1ℓ(ℓ+1)bℓrℓ+2P0ℓ(cosθ), (10) Bθ =−1rsinθ∂α∂r=∞∑ℓ=1ℓbℓrℓ+2P1ℓ(cosθ) (11)

outside the star.

### 2.2. Magnetic field inside the star: the GS equation.

Throughout this work, we considered a barotropic equation of state, , in which case if we divide Eq. (2) by and take the curl of the result, the Lorentz force per unit mass must be the gradient of a certain axisymmetric scalar function , , where the factor is introduced for convenience. Replacing the decomposition of Eq. (3) in the Lorentz force, it turns out that must be a function of as well. Hence , and the so-called Grad–Shafranov equation (hereafter GS, Grad & Rubin 1958, Shafranov 1966),

 Δ∗α+β(α)β′(α)+ρr2sin2θχ′(α)=0, (12)

is obtained, where is the Grad–Shafranov operator given by Eq. (5). In the latter, a prime indicates derivative with respect to the argument. In this way, for given functions and , the GS equation may be solved for the unknown . This can be done self-consistently with the Poisson equation (relating the gravitational potential and the density profile) and the Euler equation, provided an equation of state. Imposing appropriate boundary conditions at the surface, this procedure gives , which determines the magnetic field inside the star. This approach has been successfully used to get (numerical) barotropic equilibria (Lander & Jones, 2009; Fujisawa et al., 2012), but in what follows we will briefly discuss an additional, useful approximation which allowed us to obtain suitable equilibria with less calculations.

It is clear that the magnetic field can deform the star, so that is in principle affected by . Nevertheless, the fact that the magnetic energy stored in the stars considered is much smaller than their gravitational energy (Reisenegger, 2009) suggests that this deformation is correspondingly small, so the density profile appearing in Eq. (12) can be taken to an excellent approximation to be as in the non-magnetic case. This approximation has been already considered (Ciolfi et al., 2009; Gourgouliatos et al., 2013). One advantage of this approach is that the results can be scaled to any (small) field strength, which is not true for the “self-consistent” one. Thus, in this work we looked for barotropic equilibria by solving the GS equation for a given density profile , having been particularly interested in comparing these two approaches.

### 2.3. Variational principle

Alternatively, the GS equation can also be obtained by extremizing the functional with respect to the scalar function , where the Lagrangian density is

 L=18πr2sin2θ[|∇α|2−β2(α)−2ρr2sin2θχ(α)], (13)

subject to the condition that must vanish at the stellar surface, i.e., homogeneous Dirichlet boundary conditions. This expression generalizes that given by Monaghan (1976).

### 2.4. Boundary conditions.

Maxwell’s equation imposes that the radial component of the magnetic field be continuous across the surface. Also, the and components of the magnetic field must be continuous in order to avoid surface currents. In terms of , continuity of demands that be continuous, which in turn implies that itself must be continuous,

 αin(R,θ)=αout(R,θ)=∞∑ℓ=1bℓRℓsinθP1ℓ(cosθ), (14)

where we have introduced the expansion for outside, while is a solution of the GS equation inside the star. Using the completeness of the set , one can extract the coefficients appearing in Eq. (14) in terms of the value of at the surface, leading to

 bℓ=Rℓ2ℓ+12ℓ(ℓ+1)Iℓ, (15)

where stands for the integral

 Iℓ≡∫π0P1ℓ(cosθ)α(R,θ)dθ. (16)

This choice of the coefficients ensures continuous. Continuity implies that

 ∂αin∂r∣∣∣r=R =∂αout∂r∣∣∣r=R=−1R∞∑ℓ=1ℓbℓRℓsinθP1ℓ(cosθ) =−1R∞∑ℓ=12ℓ+12ℓ+2sinθP1ℓ(cosθ)Iℓ, (17)

where in the last equality we have replaced the expression for found in Eq. (15). Finally, in order to have continuous, we demanded that at the surface, as explained in more detail in the next subsection.

### 2.5. Particular choices for the magnetic functions.

As mentioned, is an arbitrary function of only. However, its shape is physically constrained as follows. Since outside and is constant along poloidal field lines, it has to be zero along poloidal field lines closing outside the star. That is, the toroidal field may lie only in regions where the poloidal field lines close inside the star. A particularly simple choice to achieve this is to set

 β(α)={s(α−αsurf)λifα≥αsurf0ifα<αsurf, (18)

where and are free parameters, while is the value of at the surface of the star on the equator, which is also the value of along the longest poloidal field line closing inside the star. This choice has already been used by previous authors (Ciolfi et al., 2009; Lander & Jones, 2009; Fujisawa et al., 2012; Lander & Jones, 2012; Akgün et al., 2013; Gourgouliatos et al., 2013). It is fundamental to notice that geometrically one does not know a priori where to locate the boundary between the volume with and without toroidal field. Also, is not a given parameter, but it has to be determined while solving the GS equation. We set in order to compare with previous works. It is expected that higher values of result in smaller contributions of the term in the GS equation, leading to a smaller toroidal field strength, although it turns out that there are no significant differences in the configurations when using larger . The parameter , instead, plays a significant role in the configurations obtained, so a separate subsection will be devoted to this dependence. Throughout this work, we also chose

 χ′(α)=k(constant). (19)

## 3. Numerical tests and results

We developed a finite-difference code to solve the GS equation for arbitrary choices of the functions , , and . Our code solves for on a polar grid of points in the -direction and points in the -direction, corresponding to the inside of the star, and it matches the solution smoothly to the multipolar expansion of Eq. (9). Since the infinite sum defining the multipolar expansion outside cannot be performed numerically, we truncate it to a maximum multipole , defined at the beginning of the method. The details of our code are summarized in the Appendix.

We solved the GS equation for different density profiles, considering the particular function given in Eq. (18) and normalized our results so that . We also normalized distances to stellar radius and density profiles to . The results can be easily scaled to other normalizations.

### 3.1. Tests performed to our code.

After obtaining numerical solutions, we tested whether they are actually solutions by finding that, for resolutions and , the inequality

 |Δ∗α+β′β+ρχ′r2sin2θ|≤ε∣∣αmax∣∣ (20)

held for on each point on the grid. Here is the maximum value of reached by the corresponding equilibrium over the space. This inequality held for all equilibria obtained with this resolution, while better accuracy was obtained for higher resolutions, as expected. Regarding boundary conditions, we confirmed that the coefficients that define outside the star were consistently calculated by getting that the inequality

 αin/αout≤1+ε (21)

held on each point on the surface, with for , . We also tested that the radial derivative of remained continuous by obtaining that the inequality

 ∣∣ ∣∣(∂αin∂r)r=1−(∂αout∂r)r=1∣∣ ∣∣≤ε∣∣αmax∣∣ (22)

held for on each point on the surface, again, for and . For the two latter tests, better resolutions improved the accuracy, but large () values of in the prescription of yielded values of , in both cases.

Regarding the number of multipoles outside the star, one could think that, since the multipolar expansions extend to , the larger , the more accuracy one gets. However, for a fixed resolution, we found that the accuracy when calculating the multipolar coefficients drops off after passing a certain critical value of . Since we expect that the equilibria we are looking for have higher multipoles, we are interested in finding an optimal . After repeating the previous tests for different resolutions, we found that gives the best results, so we fixed this parameter to that value.

We tested our code against known analytical solutions of the GS equation, corresponding to purely poloidal fields (for which identically). Here we present a test using the purely poloidal field with (Gourgouliatos et al., 2013)

 α(r,θ)BpR2=⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩[35x216−21x48+15x616]sin2θ,x<1sin2θ2x,x>1, (23)

where and is the magnetic field at the pole. This configuration is obtained by choosing , where is the central density. For this case, we found that

 ∣∣αnum−αanal∣∣≤ε|αmax| (24)

with on each point on the grid, for and . This excellent agreement increased the confidence in our numerical scheme.

We compared some of our equilibria with those obtained by Gourgouliatos et al. (2013) in the context of Hall equilibria in neutron star crusts. In terms of physics, the derivation and the physical interpretation of some quantities involved are quite different, but the mathematical form of the equation is exactly the same. Gourgouliatos’ code is based on a Gauss–Seidel method on a cylindrical grid extending outside the star as well, and assuming as exterior boundary condition a dipolar field at the boundary of the grid. Table 1 summarizes this comparison, in which we display the toroidal-to-poloidal energy ratio for different values of . As seen, both codes agree for all equilibria tested.

### 3.2. Magnetic polytropes

Keeping the same form of given in Eq. (18), we also found configurations using density profiles derived from non-magnetic equilibria with a polytropic relation between pressure and density,

 P0(ρ0(r))=K[ρ0(r)]1+1/n, (25)

where is a constant depending on physical parameters of the star, and is the polytropic index. As known, polytropes of index are not of physical interest as the star extends to infinity. Moreover, non-magnetic polytropes with are known to be unstable against radial perturbations (see for example Shapiro & Teukolsky 1983). Thus we restricted ourselves to obtain configurations up to , although, as said before, the code accepts any density profile.

As illustrated in Fig. 1 for the particular case of polytropes, for small and moderate values of , we found that the strength of the poloidal component is almost one order of magnitude higher than the toroidal one, while the poloidal field itself is around 10 times stronger near the center of the star compared to the field at the surface. For the latter range of , the configurations are mostly dipolar, with almost negligible relative amplitude of higher multipoles. A higher value of produces a configuration with a more substantial contribution of higher multipoles, although the dipole still dominates. As expected, increasing also increases the relative strength of the toroidal field with respect to the poloidal one, so, eventually, the amplitude of the toroidal field becomes comparable to the poloidal one. However, as Fig. 2 illustrates, the volume where the toroidal component is present shrinks in all cases, forming a thin ring of nearly circular cross section tangent to the surface of the star, a result already reported in the literature (Lander & Jones, 2009; Fujisawa et al., 2012; Gourgouliatos et al., 2013; Pili et al., 2014). Mathematically, this shrinkage can be understood as follows: in the GS equation, the term is an explicit function of the position but not of when constant. In that case, and considering that increases when choosing a larger , the term must become more negative because remains the same regardless of . In turn, the term is similar to the Laplacian of (recall Eq. (5)), that is, it depends on the curvature of , so the larger is, the more compact the poloidal lines are and the steeper the maximum in will be. The circular shape of the cross section occurs because, in the limit of , the toroidal current flowing along the ring is locally a straight line around which, since there is no matter to deform them, poloidal lines must form circular rings. This result is in agreement with the numerical results of Lander & Jones (2009), Fujisawa et al. (2012), Gourgouliatos et al. (2013), Pili et al. (2014) and Bucciantini et al. (2015), although Ciolfi et al. (2009) find a different behavior.

### 3.3. Fraction of the energy in the toroidal component

Table 2 shows the fraction of toroidal energy for the magnetic polytropes described before. In all cases, this fraction is only a few % of the total energy. Moreover, Fig. 3 suggests that is bounded by a maximum value when plotted as a function of the parameter , regardless of the density profile (see also Ciolfi et al. 2009, Lander & Jones 2009, Fujisawa et al. 2012, Gourgouliatos et al. 2013, and Pili et al. 2014). Table 3 shows the values of at which the maxima occur.

The existence of a maximum value for can be understood in terms of two competing effects: for low values of , the poloidal field is hardly affected by the weak toroidal field. Thus, the volume containing the latter does not change much with , , and , so

 EtorEmag=EtorEpol+Etor≈EtorEpol∝s2, (26)

i.e., the toroidal ratio increases quadratically as a function of . As increases more, so does the toroidal field, and the dependence with becomes less trivial. Since the region where this component lies shrinks, the volume integral defining will stop increasing at some and the toroidal fraction will start to decrease. In other words, the interplay between increasing while shrinking the volume where it is present produces an upper bound for . These results are qualitatively still true for other equations of state (Fujisawa et al., 2012) and when relativistic effects are included (Ciolfi et al., 2009; Pili et al., 2014).

### 3.4. On the validity of “weak” magnetic field

All our calculations are based on the assumption of a weak magnetic field, in the sense that the magnetic forces do not exert a significant effect on the stellar structure, i.e., to first order, those forces do not deform the primarily spherical shape of the star. Making use of this assumption, we took the simpler approach of obtaining magnetic equilibria by solving the GS equation only, assuming a fixed spherical density profile coming from a non-magnetic background star, instead of solving self-consistently the whole system of equations (Euler equation + Poisson equation + Maxwell equations) provided an equation of state , i.e., obtaining not only but also the fluid quantities, as functions of position. It is natural to wonder about the accuracy of the simpler approach used throughout this work compared to the self-consistent scheme. The latter method was already used by Lander & Jones (2009, 2012) (L&J) in the same astrophysical context. Table 4 exhibits a comparison of the fraction of toroidal energy for equilibria obtained by L&J versus our work, for two different polytropic indices. With these values, we can estimate the expected discrepancy in neglecting the effect of the magnetic field on the fluid as , which is in acceptable agreement with the relative error shown in Table 4. Since the magnetic field strength needed to generate a minimum distortion on the fluid is at least one order of magnitude above the fields for the objects considered here, then the approximation discussed should hold for the whole relevant range of magnetic fields. This has the implication of enormously simplifying the process of finding barotropic equilibria.

## 4. Asymptotic, analytic solutions

Numerical instabilities arise when increasing beyond a value (for fixed ) and our code has convergence problems. This occurs because more points on the grid are needed to resolve the toroidal volume, which becomes smaller. In this section we overcome this difficulty by introducing an asymptotic model for large , which allows us to study the global properties of the equilibria in the limit .

In order to motivate the model, we note that the term in the GS equation approaches zero as one approaches the stellar surface, because it is proportional to , which is expected to decrease monotonically until it vanishes at the surface. Thus, it is also expected that, as the volume occupied by the toroidal field shrinks and gets close to the surface, this term becomes significantly smaller in magnitude than the other two inside the volume, so that and roughly cancel each other. In fact, from our simulations we confirmed that

 r2sin2θρ|Δ∗α|∼r2sin2θρββ′≲10−3 (27)

for (again, we normalize ). In this way, one can expect that the relation

 Δ∗α+ββ′≈0 (28)

will become more and more accurate inside the toroidal volume when increasing . This is the central idea of the field model we are introducing: for large , we assume that Eq. (28) is valid, and solve for under that approximation. Of course, when is determined, the problem is solved and then one can obtain expressions for the toroidal flux, energy, and so on. The problem is even simpler if we recall that the volume where the toroidal field lies not only shrinks but also becomes a toroid regardless of the density profile. Because of this, we introduce a new coordinate system on the meridional half plane of the star, see Fig. LABEL:newcoords. From our simulations, it turns out that the smaller the volume, the more independent is with respect to the angle , so that for a sufficiently large , we can assume that does not depend on , but only on , therefore the GS operator reads (recall Eq. (5))

 Δ∗α≈∇2α=1ϵddϵ(ϵdαdϵ), (29)

where we have used that . Using the function of Equation (18), Equation (28) inside the toroidal volume can be written as

 1uddu(udηdu)=−η2λ−1, (30)

where we have defined the dimensionless function so that , with being the maximum value of (reached at , or ), and the dimensionless radius

 u≡λ1/2(αmax−αsurf)λ−1sϵ. (31)

Eq. (30) is a second-order ordinary differential equation for , hence two boundary conditions are needed to solve it. Firstly, and by definition, it is clear that . Also, , where the prime stands for the derivative with respect to the argument. This is because (and thus ) has a smooth maximum at . Eq. (30) must be integrated until the boundary of the toroidal volume, corresponding to , or , where is the first root of , a number to be determined. For there is no toroidal field, so this equation (and its solution) is valid for only. Eq. (30) accepts analytical solutions for a few exponents . For , the unique solution satisfying the boundary conditions is , with . For . the equation becomes a Bessel equation of order zero, with unique solution , being the first root. For all other values of , the equation is nonlinear, with no obvious analytical solution, so we wrote a fourth-order Runge-Kutta code to solve it numerically. Of course, is a pure number, depending on the value of and on nothing else. \@floatfigure\end@float

### 4.1. Toroidal quantities

Under the above assumptions, formulæ for the toroidal flux and energy are readily obtained. By definition, the toroidal flux through a meridional plane is

 Φtor ≈∫ϵsurf0β(ϵ)rmax2πϵdϵ =2πqλλrmaxs(α% max−αsurf)2−λ, (32)

where we defined

 qν≡∫usurf0uηνdu, (33)

a pure number. The toroidal energy can be computed as

 Etor ≈∫ϵsurf0β2(ϵ)8πr2max4π2rmaxϵdϵ =πq2λ2λrmax(α% max−αsurf)2. (34)

In principle, we are interested in obtaining analytical expressions as functions of the parameter . In Eqs. (4.1) and (4.1), , , and (which to first approximation equals the radius ) depend implicitly on , so it is desirable to obtain these dependences. This can be done as follows. First, let us consider the hypothetical case in which the toroidal field is effectively concentrated in an infinitesimally thin loop at the surface along the equator of the star. This would be the “” case. This limit of an infinitesimally thin loop is analogous to the problem of a thin, circular current loop. Outside this loop, is given by the solution of the GS equation for a purely poloidal field, namely,

 α=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∞∑ℓ=1[aℓ(rR)ℓ+1−f(r)δℓ1]sinθP1ℓ(cosθ)rR, (35)

where is the Kronecker delta and

 f(r)=k3r∫r0ξρ(ξ)(ξ3−r3)dξ, (36)

where (a constant) regulates the amplitude of the magnetic field (Gourgouliatos et al., 2013). The coefficients and must be fixed using boundary conditions. Continuity of at implies

 aℓ−f(R)δℓ1=bℓ. (37)

This also ensures continuity of the radial component of the magnetic field. This time, however, the tangential component of the magnetic field is not continuous at the boundary , but it has a discontinuity due to the presence of a thin loop carrying a current , with current density

 J=IRδ(cosθ)δ(r−R)eϕ. (38)

Here is the Dirac delta function. This current comes from a singular toroidal field, which produces an azimuthal current density

 4πcJϕ=(∇×B)|ϕ=−Δ∗αrsinθ=ββ′rsinθ, (39)

where in the last equality we have used Eq. (28). The second boundary condition to impose is

 ∂α∂r∣∣∣R+−∂α∂r∣∣∣R−=−4πIcsinθδ(cosθ). (40)

By direct replacement of the respective expressions for for and for , one can extract the coefficients by standard mathematical tools in order to finally get

 aℓ=f(R)+Rf′(R)3δℓ1+2πIRcP1ℓ(0)ℓ(ℓ+1), (41)

from which can be obtained from Eq. (37). Thus, if is known, the problem in this case is completely determined.

The motivation behind pointing out this hypothetical case is what follows. The general solution in Eq. (35) is of course also valid outside the toroidal region when the latter is a thin (but not infinitesimally thin) loop. In that case, however, the coefficients are not readily obtained as we would have to impose boundary conditions in a cumbersome geometry, assuming that the solution to the GS equation inside the toroidal volume is precisely that satisfying Eq. (30). Nevertheless, as the parameter increases, the coefficients should eventually converge to those in Eq. (41). Therefore, in the limit of a thin toroidal volume with finite cross section, it is expected that the solution to the GS equation will be that in Eq. (35) with coefficients and given approximately by those in (41), where is the current along the thin toroidal region. If this is the case, the current flowing can be easily obtained using the asymptotic model developed in this section. By definition, and recalling Eq. (39), we get

 I ≈c4π∫ϵsurf0λs2(αmax−αsurf)2λ−1r%maxη2λ−12πϵdϵ =cq2λ−12R(αmax−α% surf), (42)

where we have replaced . Eq. (42) implies that in this approximation, , and hence the coefficients, are implicit functions of . The asymptotic solution of the GS equation for large , outside the toroidal volume, is that given in Eq. (35), where the coefficients are approximately

 aℓ≈f(R)+Rf′(R)3δℓ1+q2λ−1(αmax−αsurf)πP1ℓ(0)ℓ(ℓ+1) (43)

and

 bℓ≈Rf′(R)−2f(R)3δℓ1+q2λ−1(αmax−αsurf)πP1ℓ(0)ℓ(ℓ+1). (44)

So, with this procedure, we can write an approximate solution for the GS equation in terms of its maximum value and its value at the surface between the region with and without toroidal field. In addition, and its derivatives must be continuous across the latter surface, i.e.,

 αsurf=α(R−ϵsurf,π/2) (45)

and

 αmax−αsurfϵsurfq2λ−1=∂α∂r(R−ϵsurf,π/2), (46)

respectively. Introducing the explicit form of the coefficients in the expansion for (Equation (35)), Equations (45) and (46) read

{widetext}
 αsurf=kρcR4[f(1−x)−f(1)+f′(1)3(1−x)2]+πq2λ−1(αmax−αsurf)F(1−x) (47)

and

 αmax−αsurfx=kρcR4q2λ−1[f′(1−x)−2(f(1)+f′(1))3(1−x)]+π(αmax−αsurf)F′(1−x), (48)

respectively, where we have redefined and then dropped the tilde. Also, we have introduced the parameter and the function

 F(z) ≡∞∑ℓ=1[P1ℓ(0)]2ℓ(ℓ+1)zℓ+1=z2∞∑ν=02ν+12ν+2[P2ν(0)]2z2ν =z28[42F1(12,12,2,z2)+z22F1(32,32,3,z2)]. (49)

Here is the Gauss hypergeometric function, so that converges for , while means derivative evaluated at . The interesting point with the previous formulæ is that Eqs. (47) and (48) allow to write both and in terms of the parameter , which is an implicit function of . Once one solves for and , along with Eq. (31) evaluated at ,

 usurf=λ1/2(αmax−αsurf)λ−1sϵsurf, (50)

where is a pure number, one gets the equations needed to write and as a function of . Since by construction the poloidal flux equals , we can define the quantity , which is the maximum poloidal flux reached by the configuration, and take the ratio . Using (4.1), we get

 ΦtorΦpol=qλλ1/2usurfx1−xαmax(x)−αsurf(x)αmax(x), (51)

which turns out to be an explicit function of . This equation can be combined with

 s=usurfλ1/2Rx(αmax(x)−αsurf(x))λ−1, (52)

(see Eq. (50)) which is also an explicit function of , to obtain a parametric relation between and . Eq. (51) allows to compare this model with the simulations.

### 4.2. The model versus the simulations

As discussed, the limit of having an infinitesimally thin toroidal volume is reached when , or . Fig. 4 displays vs. obtained from eqs. (51) and (52), as well as from the numerical simulations. In all cases, as increases, the asymptotic model approaches the numerically obtained curve, in good agreement with the approximations assumed within the model.

Having tested the model with the simulations, the question whether is bounded over the full range of can be answered. In this case, there is no trivial expression for as one should in principle integrate the asymptotic solution over all space outside the thin toroidal volume, and add the contribution inside the torus. This would demand non-trivial integration limits for the former integral, although the latter can be easily computed with the model of this section. In order to estimate the poloidal energy, we can use the known result for a circular loop of radius having a circular cross section of radius , with (e.g. Jackson 1975)

 Epol=12LI2;L≈4πRc2[ln(8Rϵsurf)−74], (53)

which diverges in the limit . The toroidal energy, on the other hand, goes to zero as , so the fraction also vanishes in that limit. The conclusion is that , as a function of , does possess a maximum value, and then decreases monotonically to zero. This would mean that all the configurations with as given in Eq. (18) and are restricted to have a small () fraction of the magnetic energy stored in the toroidal component. This choice has been the most used in the literature, so if one were interested in obtaining higher fractions, one should explore other prescriptions for (and perhaps for as well).

## 5. Summary and conclusions

We have developed a numerical code to model axially symmetric magnetic stars in ideal MHD equilibrium, allowing for both poloidal and toroidal component of the magnetic field, in a barotropic fluid. We presented several tests to check its accuracy, including comparisons with previous works in this area, finding a good agreement with them. From these comparisons, we showed that, for magnetic field strengths up to G, and assuming polytropic equations of state, the simpler approach of solving the resulting GS equation describing the equilibrium, and imposing a density profile derived from a non-magnetic equilibrium, gives essentially the same results as when obtaining equilibria with the self-consistent scheme of solving not only for the magnetic field, but also for the fluid quantities provided an equation of state describing the fluid.

Using our code, we found a relatively wide range of barotropic equilibria whose magnetic fields combine an internal toroidal field and a poloidal field extending to the exterior of the star as well. Fixing , numerical equilibria described by the piecewise function in Eq. (18) are poloidal-dominated, in the sense that the fraction of magnetic energy stored in the toroidal component of the magnetic field is , even for configurations with comparable poloidal and toroidal magnetic field strength. Our numerical simulations also confirm previous results showing a maximum value of this fraction as a function of the parameter appearing in . These properties do not depend significantly on the density profile assumed, which in this work varied between several very different polytropes.

Numerical simulations break down for large values of (). In order to obtain results beyond this limitation, we developed an analytical model which reproduces the behavior of the numerical solutions for large toroidal field. This model was able to mimic the main global properties in this regime, finding that both the fractions of magnetic fluxes and energies do have a maximum when plotted against , and decay to zero when goes to infinity.

It is likely that functions providing a larger region of poloidal field lines closing inside the star give larger toroidal energy, as proposed by Ciolfi & Rezzolla (2013). This is something we can readily explore as our code accepts arbitrary functions and . Also, having obtained barotropic equilibria allows to study their stability. This can be done either through time-dependent MHD simulations using the configurations we have obtained as initial condition, or by means of analyzing the change of the energy in the system when slightly perturbing the fluid inside the star, in an analogous treatment to that for non-barotropic, stably stratified stars followed by Akgün et al. (2013).

## Acknowledgments

This work was supported by CONICYT International Collaboration Grant DFG-06, FONDECYT Regular Grants 1110213, 1110135, 1150411 and 1150718, the Basal Center for Astrophysics and Associated Technologies (CATA; PFB-06), and a CONICYT Master’s Fellowship. We thank K. N. Gourgouliatos, S. K. Lander, K. Fujisawa, and P. Marchant for having provided part of their results, which yielded very stimulating discussions, and the anonymous referee for useful comments on the first draft of this paper.

## Appendix A Details on our numerical code

We are interested in obtaining the function that solves the generalized GS equation,

 Δ∗α=−F(r,θ,α(r,θ)),F(r,θ,α(r,θ))≡β(α)β′(α)+r2sin2θρ0(r)χ′(α). (A1)

We do this by writing a discrete (finite-difference) version of this equation and solving the resulting system of algebraic equations. We discretize the space inside the star into a polar grid of regular intervals in the radial direction and in the angular direction, of length and each, respectively. On this grid, each point is labeled by a pair , with coordinates and . Points along the axis are those with (northern hemisphere) and (southern hemisphere) and points on the surface of the star correspond to . Any quantity evaluated on a point is labeled as . The discretized GS equation reads

 αi+1,j−2αi,j+αi−1,j(Δr)2+1r2iαi,j+1−2αi,j+αi,j−1(Δθ)2−cotθjr2iαi,j+1−αi,j−12Δθ=−F(ri,θj,αi,j), (A2)

where we use central difference for the derivatives. This algebraic equation is valid for , . The value of at points lying on the axis (, and ) must remain fixed equal to zero, providing a boundary condition. Thus Eq. (A2) actually gives algebraic equations for the unknown variables , , , so () additional equations are required in order to close the system, which are obtained from the boundary conditions at the stellar surface. Demanding continuity of the radial derivative of gives the required number of equations,

 αNr+1,j−αNr−1,j2Δr=−1Rℓmax∑ℓ=12ℓ+12ℓ+2sinθjP1ℓ(cosθj)Iℓ≡gj (A3)

(see Eq. (17)). Variables , standing for the value of at points just outside the star, can be calculated by assuming that outside the star is written by the multipolar expansion in Eq. (9) (with ),

 αNr+1,j=αout(R+Δr,θj), (A4)

again, for , where we have introduced the parameter because we cannot perform the sum to infinity. The integral defined in Eq. (16), and involved in expansions and , can be calculated through Simpson’s rule and using the points at the surface, . This gives

 Iℓ≈(Nθ−2)∑j=0jevenΔθ3[P1ℓ(cosθj)αNr,j+4P1ℓ(cosθj+1)αNr,j+1+P1ℓ(cosθj+2)αNr,j+2], (A5)

so must be chosen as an even natural number. Of course, using this form of the expansions, we are explicitly assuming continuity of and its derivative with respect to .

When introducing the explicit form of in terms of the ’s, Eq. (A3) gives new independent equations relating part of the unknowns on the grid, so a consistent solution of the algebraic equations may be carried out. The final system of equations to solve is

 Gi,j≡αi+1,j−2αi,j+αi−1,j(Δr)2+1r2iαi,j+1−2αi,j+αi,j−1(Δθ)2−cotθjr2iαi,j+1−αi,j−12Δθ+F(ri,θj,αi,j)=0, (A6)

for , , and

 GNr,j≡αNr+1,j−αNr−1,j2Δr−gj=0, (A7)

for , where

 αNr+1,j=Δθ6sinθj(Nθ−2)∑~ȷ=0~ȷeven[αNr,~ȷTj,~ȷ+4αNr,~ȷ+1Tj,~ȷ+1+αNr,~ȷ+2Tj,~ȷ+2] (A8)

and

 gj=−Δθ6Rsinθj(Nθ−2)∑~ȷ=0~ȷeven[αNr,~ȷUj,