Relaxed Equilibrium Configurations to Model Fossil FieldsI – A first family

# Relaxed Equilibrium Configurations to Model Fossil Fields I – A first family

V. Duez CEA/DSM/IRFU/SAp, CE Saclay, F-91191 Gif-sur-Yvette Cedex, France; AIM, UMR 7158, CEA - CNRS - Université Paris 7, France
11email: vincent.duez@cea.fr,stephane.mathis@cea.fr
S. Mathis CEA/DSM/IRFU/SAp, CE Saclay, F-91191 Gif-sur-Yvette Cedex, France; AIM, UMR 7158, CEA - CNRS - Université Paris 7, France
11email: vincent.duez@cea.fr,stephane.mathis@cea.fr
Received 19 October 2009 / Accepted 06 March 2010
###### Key Words.:
Magnetohydrodynamics (MHD) – Plasmas – Magnetic fields – Sun: magnetic fields – Stars: magnetic fields
###### Abstract

Context:The understanding of fossil fields origin, topology and stability is one of the corner stones of the stellar magnetism theory. On one hand, since they survive over secular time-scales, they may modify the structure and the evolution of their host stars. On the other hand, they must have a complex stable structure since it has been demonstrated by Tayler and collaborators that simplest purely poloidal or toroidal fields are unstable on dynamical time-scales. In this context, the only stable configuration which has been found today is the one resulting of a numerical simulation by Braithwaite and collaborators who have studied the evolution of an initial stochastic magnetic field, which is found to relax on a mixed stable configuration (poloidal and toroidal) that seems to be in equilibrium and then diffuses.

Aims:In this work, we thus go on the track of such type of equilibrium field in a semi-analytical way.

Methods:In this first article, we study the barotropic magnetohydrostatic equilibrium states; the problem reduces to a Grad-Shafranov-like equation with arbitrary functions. Those latters are constrained by deriving the lowest-energy equilibrium states for given invariants of the considered axisymmetric problem and in particular for a given helicity which is known to be one of the main actor of such problems. Then, we obtain the generalization of the force-free Taylor’s relaxation states obtained in laboratory experiments (in spheromaks) that become non force-free in the self-gravitating stellar case. The case of general baroclinic equilibrium states will be studied in Paper II.

Results:Those theoretical results are applied to realistic stellar cases, namely to the solar radiative core and to the envelope of an Ap star, and discussed. In both cases we assume that the field is initially confined in the stellar radiation zone.

Conclusions:

## 1 Introduction

Spectropolarimetry is nowadays exploring the stellar magnetism across the whole Hertzsprung-Russel diagram (Donati et al. 1997, 2006; Neiner 2007; Landstreet et al. 2008; Petit et al. 2008). Furthermore, helioseismology and asteroseismology are providing new constraints on internal transport processes occuring in stellar interiors (Turck-Chièze & Talon 2008; Aerts et al. 2008). In this context, even if standard stellar models explain the main features of stellar evolution, it is now crucial to go beyond this modelling to introduce dynamical processes such as magnetic field and rotation to investigate their effects on stellar structure and secular evolution (Maeder & Meynet 2000; Talon 2008). To achieve this aim, secular MHD transport equations have been derived in order to be introduced in stellar evolution codes. They take into account in a coherent way the interaction between differential rotation, turbulence, meridional circulation, and magnetic field (Spruit 2002; Maeder & Meynet 2004; Mathis & Zahn 2005), while non-linear numerical simulations provide new insight on these mechanisms (Charbonneau & MacGregor 1993; Rudiger & Kitchatinov 1997; Garaud 2002; Brun & Zahn 2006). If we want to go further, the simplest modifications of static structural properties such as density, gravity, pressure, temperature, and luminosity induced by the magnetic field have also to be systematically quantified as a function of the field geometry and strength (Moss 1973; Mestel & Moss 1977; Lydon & Sofia 1995; Couvidat et al. 2003; Li et al. 2006; Duez et al. 2008; Li et al. 2009; Duez et al. 2010).

However, an infinity of possible magnetic configurations can be investigated since the different observation techniques only lead to indirect indications on the internal field topologies through the surface field properties they provide. Furthermore, since the simplest geometrical configurations like purely poloidal and purely toroidal fields are known to be unstable (Acheson 1978; Tayler 1973; Markey & Tayler 1973, 1974; Goossens & Veugelen 1978; Goossens & Tayler 1980; Goossens et al. 1981; van Assche et al. 1982; Spruit 1999; Braithwaite 2006, 2007), the best candidates for stable geometries are mixed poloidal-toroidal fields (Wright 1973; Markey & Tayler 1974; Tayler 1980; Braithwaite 2009).

Therefore, it is necessary to go on the track of possible stable magnetic configurations in stellar interiors to evaluate their effects on stellar structure and to use them as potential initial conditions to study secular internal transport processes.

In this work, we thus revisit the pioneer works by Ferraro (1954); Mestel (1956); Prendergast (1956) and Woltjer (1960). Ferraro (1954) studied the equilibrium configurations of an incompressible star with a purely poloidal field. Prendergast (1956) (see also Chandrasekhar 1956a; Chandrasekhar & Prendergast 1956; Chandrasekhar 1956b) then extended the model to take into account the toroidal field, by solving the magneto-hydrostatic equilibrium of incompressible spheres. The obtained configurations seem to be relevant in regard of the most recent numerical simulations that may explain fossil fields in early-type stars, white dwarfs or neutron stars (Braithwaite & Spruit 2004; Braithwaite & Nordlund 2006; Braithwaite 2006) and of theoretical studies of their helicity relaxation (Broderick & Narayan 2008; Mastrano & Melatos 2008). The main generalization of Prendergast’s work we achieve here consists in relaxing the incompressible hypothesis in order to take into account the star’s structure (see also Woltjer 1960), which differs as a function of its stellar type and of its evolution stage, and to derive the minimum energy equilibrium configuration for a given mass and helicity, which are then applied to realistic models of stellar interiors.

Assuming that the Lorentz volumetric force is a perturbation compared with the gravity, we derive the non force-free magnetohydrostatic equilibrium. In this first article, we then focus on the barotropic equilibrium states family 111Barotropic states are such that their density and pressure gradients are aligned. They can be convectively stable or not, depending on their entropy stratification. We will precisely introduce their definition in §2.2. , for which the possible field configurations and the stellar structure are explicitely coupled; these may correspond to the numerical experiments by Braithwaite and collaborators. In this case, the problem reduces to a Grad-Shafranov-like equation (Grad & Rubin 1958; Shafranov 1966; Kutvitskii & Solov’ev 1994), similar to the one intensively used in fusion plasma physics. We then focus on its minimum energy eigenmodes for a given mass and helicity, which are derived and applied in order to model relaxed stellar fossil magnetic fields which are found to be non force-free. Arguments in favor of the stability of the obtained configurations are finally discussed (Wright 1973; Tayler 1980; Braithwaite 2009; Reisenegger 2009) and we compare their properties with those of relaxed fields obtained in numerical simulations (Braithwaite 2008). The case of general baroclinic equilibrium states will be studied in Paper II (Wright 1969; Moss 1975).

## 2 The Non Force-Free Magneto-Hydrostatic Equilibrium

In this work, we focus on the magnetic equilibrium of a self-gravitating spherical shell to model fossil fields in stellar interiors. To achieve this goal, we start from

 {0}=−∇P−ρ∇V+FL,whereFL={j}×% {B}. (1)

Eq. (1) must be satisfied in the interior of an infinitely conducting mass of fluid in the presence of a large-scale field together with the Poisson equation, , and the Maxwell equations, (Maxwell flux) and (Maxwell-Ampère).
, and are respectively the pressure, the density and the gravitational potential of the considered plasma. B is the magnetic field and j is the associated current, which is given in the classical MHD approximation by the Maxwell-Ampère’s equation. is the magnetic permeability of the plasma and is the Lorentz force.

### 2.1 Magnetic field configuration and the magnetohydrostatic equilibrium

If we consider only the axisymmetric case, where all physical variables are independent of the azimuthal angle (), can be written in the form

 {B}=1rsinθ∇Ψ(r,θ)×^eφ+1rsinθF(r,θ)^eφ, (2)

which is divergenceless. and are respectively the poloidal flux function and the toroidal potential; are the usual spherical coordinates and their unit-vector basis. Finally, the poloidal component of the magnetic field () is such that , so it belongs to iso- surfaces. The magnetohydrostatic equilibrium (Eq. 1) implies that the poloidal part (in the meridional plane in the axisymmetric case) of the Lorentz force () balances the pressure gradient and the gravitational force, which are also purely poloidal vectors, while, in the absence of any other force, its toroidal component () vanishes. We thus have

 FL={F}LP+FLφ^eφ=% {F}LP. (3)

Using Eq. (2), we obtain

 {F}LP=−1μ0r2sin2θ{(F∂rF+∂rΨΔ∗Ψ)^er +1r(F∂θF+∂θΨΔ∗Ψ)^eθ}, (4)

with , and

 Δ∗Ψ≡∂rrΨ+sinθr2∂θ(1sinθ∂θΨ) (5)

is the usual Grad-Shafranov operator in spherical coordinates. On the other hand, since , we get ; the non-trivial values for are therefore obtained by setting

 F(r,θ)=F(Ψ). (6)

Then, we obtain and , leading to the final expansion of the Lorentz force

 {F}L=A(r,θ)∇Ψ, (7)

where

 A(r,θ)=−1μ0r2sin2θ(F∂ΨF+Δ∗Ψ). (8)

Therefore, the poloidal component of is nonzero a priori, the field being thus non force-free in this case.

This point has here to be discussed. First, Reisenegger (2009) demonstrated that the magnetic field can not be force-free everywhere in stellar interiors (see the demonstration in the appendix A of his paper). Note in this context that the “force-free” configurations obtained by Broderick & Narayan (2008) verify this theorem because they have current sheets with a non-zero Lorentz force on the stellar surface. Moreover, Shulyak et al. (2007, 2009) showed how the atmosphere of a CP star can be the host of a non-zero Lorentz force. Therefore, from now on, we consider the non force-free equilibrium.

Now, if we take the curl of Eq. (1), we get the static vorticity equation

 (9)

that governs the balance between the baroclinic torque (left-hand side; see Rieutord (2006) for a detailed description) and the magnetic source term. Then, as has been emphasized by Mestel (1956), the different structural quantities such as the density, the gravitational potential, and the pressure relax in order to verify Eq. (1) for a given field configuration (see Sweet (1950); Moss (1975), and Mathis & Zahn (2005) §5.). Thus, the choice for is left free.

### 2.2 The barotropic equilibrium state family

Magnetic initial configurations are one of the crucial unanswered question for the modelling of MHD transport processes in stellar interiors. To examine this question, Braithwaite and collaborators (Braithwaite & Spruit 2004; Braithwaite & Nordlund 2006) studied the relaxation of an initially stochastic field in models of convectively stable stellar radiation zones. The field is found to relax, after several Alfvén times, to a mixed poloidal-toroidal equilibrium configuration, which then diffuses towards the exterior.

We choose here, using an analytical approach, to find such field geometries, which are governed at the beginning by the magnetohydrostatic equilibrium.

To achieve this aim, we focus in this first article on the particular barotropic equilibrium states (in the hydrodynamic meaning of the term) for which the field configuration is explicitely coupled with the stellar structure, since we have in this case

 (10)

Those are the generalization of the Prendergast’s equilibria which take into account the compressibility and which have been studied in polytropic cases by Woltjer (1960); Wentzel (1961); Roxburgh (1966); Monaghan (1976).

Let us first recall the definition of the barotropic states. In fluid mechanics, a fluid is said to be in a barotropic state if the following condition is satisfied (see Pedlosky (1998) in a geophysical context and Zahn (1992) in a stellar one):

 ∇ρ×∇P={0}; (11)

in other words, the baroclinic torque in the vorticity equation (Eq. 9) vanishes. Then, the surfaces of equal density coincide with the isobars since the density and the pressure gradients are aligned. This does not imply any question of equation of state, which in stellar interiors can take the most general form ( being the temperature). Moreover, this does not presume anything about the stratification of the fluid, which can be stably stratified or not. For example, a star in solid body rotation is in a barotropic state (as opposed to baroclinic) (see once again Zahn 1992). Let us illustrate this point with the simplest case of a non-rotating and non-magnetic stellar radiation zone. In this case, the hydrostatic balance is given by . If we take the curl of this equation, we obtain the stationary version of the thermal-wind equation:

 −∇ρ×∇Pρ2=−∇×[∇V]={0} (12)

and the star is thus in a barotropic state in the hydrodynamic meaning of the term.

This has to be distinguished from the point of view of thermodynamics where a barotropic equation of state is such that while a non-barotropic equation of state is such that .

Then, it is clear that a fluid with a barotropic equation of state is automatically in an hydrodynamical barotropic state. However, in the case of a fluid with a non-barotropic equation of state, the situation is more subtle: in the case where the curl of the volumetric perturbing force vanishes (i.e. ) the fluid is in an hydrodynamical barotropic state, while in the general case it is in a baroclinic situation. Then, a fluid with a non-barotropic equation of state can be in a barotropic state even if it is only for a specific form of the perturbing force. In this first work, we choose to examine the first equilibrium family in which the Lorentz force verify the barotropic balance described by Eq. (11) in a stably stratified radiation zone. The second general case (cf. Mestel 1956) will be studied in paper II.

Stellar interiors, except just under the surface, are in a regime where , being the plasma’s magnetic pressure. On the other hand, in the domain of fields amplitudes relevant for classical stars (i.e. the non-compact objects), the ratio of the volumetric Lorentz force by the gravity is very weak. Therefore, the stellar structure modifications induced by the field can be considered as perturbations only from a spherically symmetric background (Haskell et al. 2008). Then, we can write , where and are respectively the mean density on an isobar, which is given at the first order by the standard non-magnetic radial density profile of the considered star, and its magnetic-induced perturbation on the isobar (with ). Thus to the first order, Eq. (10) on an isobar becomes

 −∇˜ρ×{g}eff¯¯¯ρ=∇×({F}L¯¯¯ρ)={0}, (13)

where the effective gravity , such that , has been introduced. This gives, using Eq. (7)

 ∇(A¯¯¯ρ)×∇Ψ={0}, (14)

which projects only along as

 ∂r(A¯¯¯ρ)∂θΨ−∂θ(A¯¯¯ρ)∂rΨ=0, (15)

so that there exists a function of such that

 A¯¯¯ρ=G(Ψ). (16)

Then, Eq. (8) leads to the following one ruling

 (17)

This equation is similar to the well-known Grad-Shafranov equation222The usual Grad-Shafranov equation is given by:
,
where the pressure is prescribed in function of . These describes the equilibrium between the magnetic force and the pressure gradient only.
which is used to find equilibria in magnetically confined plasmas such as those in tokamaks or in spheromaks (Grad & Rubin 1958; Shafranov 1966). However, here the source term is different and is directly related to the internal structure of the star through its density profile () (the general form of the Grad-Shafranov equation in an astrophysical context is discussed for example in Heinemann & Olbert (1978) and Ogilvie (1997)). Moreover, since the field has to be non force-free in stellar interiors (see the previous discussion in §2.1. and Eqs. 16, 7 and 8).

It is only applicable to the case of the barotropic state family. The equations for the general case will be studied in paper II.

### 2.3 F and G expansion

Let us now focus on the respective expansion of and as a function of .

First, since is a regular function, we can expand it in power series in :

 F(Ψ)=∞∑i=0λiRΨi, (18)

the being the expansion coefficients that have to be determined and a characteristic radius which will be identified below. On the other hand, must be regular at the center of the sphere; the first term () of the previous expansion is then excluded (cf. Eq. 2), the above expansion thus reducing to .

In the same way, can be expanded as

 G(Ψ)=∞∑j=0βjΨj. (19)

Then, Eq. (17) becomes:

 Δ∗Ψ+∑k>0ΛkR2Ψk=−μ0r2sin2θ¯¯¯ρ∞∑j=0βjΨj, (20)

where , being the usual Kronecker symbol. This is the generalization of the Grad-Shafranov-type equation obtained by Prendergast (1956) for the barotropic compressible states.

Thus, having assumed the non force-free barotropic magneto-hydrostatic equilibrium state leads to undetermined arbitrary functions ( and ) that must be constrained. To achieve this aim, we follow the method given in the axisymmetric case by Chandrasekhar & Prendergast (1958) and Woltjer (1959b) that allows to find the equilibrium state of lowest energy compatible with the constancy of given invariants for the studied axisymmetric system.

## 3 Self-Gravitating Relaxation States

### 3.1 Definitions and axisymmetric invariants

We first introduce the cylindrical coordinates where and . Then, B given in Eq. (2) becomes

 {B}(s,z)=1s∇Ψ(s,z)×^eφ+1sF(s,z)^eφ. (21)

Then, we define the potential vector such that and we get

 {B}=∇×{A}+Fs^eφwhereAφ(s,z)=Ψs. (22)

Next, we insert the expansion for the magnetic field B used by Chandrasekhar & Prendergast (1958) and Woltjer (1959a, b):

 {B}=−s∂zΦ(s,z)^es+1s∂s[s2Φ(s,z)]^ez+sT(s,z)^eφ (23)

where is the cylindrical unit-vector basis and where we identify using Eq. (21)

 Ψ=s2ΦandF=s2T. (24)

The Grad-Shafranov operator applied to can then be expressed as follow

 Δ∗Ψ=s2∇⋅(∇Ψs2)=[∂ss−1s∂s+∂zz]Ψ=s2Δ5Φ, (25)

where

We now introduce the two general families of invariants of the barotropic axisymmetric magneto-hydrostatic equilibrium states, which have been introduced by Woltjer (1959b) for the compressible case (see also Wentzel 1960):

 (26)
 III;q=∫VNq(s2Φ)TdV=∫V(s2Φ)qTdV, (27)

where and are arbitrary functions that have to be specified. Those latters are conserved as long as

 {B}⋅ˆ{e}r=0(\it i.e. Φ=T=0) (28)

on the boundaries.

### 3.2 Fossil fields barotropic relaxation states

Let us first concentrate on and relevant for fossil fields relaxation. First, if we set , we obtain

 III;0 = ∫VTdV=2π∫SBφdsdz (29) = 2π∫SBφdSφ=2πFφ

that corresponds to the conservation of the flux of the azimuthal field across the meridional plane of the star () in perfect axisymmetric MHD equilibria.

Then, if we set , we get

 (30)

where we thus identify the magnetic helicity (; see §5.1.) of the field configuration which is a global quantity integrated over the volume of the studied radiation zone.

Let us briefly discuss the peculiar role of this quantity in the search of stable equilibria. As emphasized by Spruit (2008), the magnetic helicity is a conserved quantity in a perfectly conducting fluid with fixed boundary conditions. However, in realistic conditions, rapid reconnection can take place even at very high conductivity, especially when the field is dynamically evolving, for example during its initial relaxation phase. Nevertheless, in laboratory experiments, like for example in spheromaks, the helicity is often observed to be approximately conserved, which leads to the existence of stable equilibrium configurations. In fact, if the helicity is conserved a dynamical or unstable field with a finite initial helicity () cannot decay completely, the helicity of vanishing field being zero. This is precisely what has been observed in the numerical experiment performed by Braithwaite & Spruit (2004) and Braithwaite & Nordlund (2006), where an initial stochastic field with a finite helicity decays initially but relaxes into a stable equilibrium.

In the context of laboratory low- plasmas, this process has been identified by Taylor (1974) and is thus called the Taylor’s relaxation.

For this reason, we now search as Chandrasekhar & Prendergast (1958) the final state of equilibrium, which is the state of lowest energy that the compressible star, preserving its axisymmetry, can attain while conserving the invariants , and in barotropic states.

To achieve this, we thus introduce the total energy of the system

 E = 12∫V{{B}2μ0+ρ[V+2U]}dV = 12∫V{1μ0[−s2ΦΔ5Φ+s2T2]+ρ[V+2U]}dV,

where is the specific internal energy per unit mass (Woltjer 1958, 1959b; Broderick & Narayan 2008). To obtain the minimal energy equilibrium state for the given invariants and a given helicity and azimuthal flux, we thus minimize with respect to , and . This leads, introducing the associated Lagrangian multipliers , to the following condition for a stationary energy:

 δE+∑naI;nδII;n+1∑q=0aII;qδIII;q=0. (32)

Following the method described in Chandrasekhar & Prendergast (1958) and Woltjer (1959b), we express and in function of , and . Since these variations are independent and arbitrary, their coefficients in the integrand of Eq. (32) must separately vanish, which gives

 1μ0Δ5Φ = ¯¯¯ρ∑naI;ndMn(s2Φ)d(s2Φ)+1∑q=0aII;qTdNq(s2Φ)d(s2Φ) (33) = ¯¯¯ρ∑naI;ndMn(s2Φ)d(s2Φ)+aII;1T,
 1μ0s2T=−1∑q=0aII;qNq(s2Φ)=−aII;0−aII;1s2Φ. (34)

These equations thus describe the minimal non force-free energy equilibrium states for a given helicity and azimuthal flux. From Eq. (17), we identify that is now constrained while required to ensure the non force-free character of the field is still arbitrary.

Let us now consider the first invariants family () given in Eq. (26) which are thus necessary to constrain . First, the non-magnetic global quantity, which is an invariant of the considered equilibrium, is the total mass of the stellar radiation zone . We thus set , leading naturally to consider the mass

 II;0=∫V¯¯¯ρdV=MRZ. (35)

However, since , we have thus to consider the highest-order invariant

 II;1=∫V(s2Φ)¯¯¯ρdV (36)

because of the non force-free behaviour of the field. This last invariant has been introduced by Prendergast (1956) in the axisymmetric non force-free magneto-hydrostatic incompressible equilibrium and it corresponds to the mass conservation in each flux tube described by the closed magnetic surface .

Furthermore, the considered radiation zone is stably stratified. Since in stellar interiors the magnetic pressure is very much less than the thermal one, the Lorentz force have only a negligible effect on the gas pressure (). Moreover, energy is required to move fluid elements in the radial direction because work has to be done against the buoyant restoring force which is thus very strong compared to the magnetic one. Therefore, the radial component of the displacement () which takes place during the adjustment to equilibrium is inhibited and due to the anelastic approximation justified in stellar radiation regions. Therefore, as emphasized by Braithwaite (2008), the mass transport in the radial direction is frozen (no matter can leave or enter in the flux tube) and can be used as a supplementary constrain in our variational method.

From Eqs. (33-34), we thus get the following equations describing the barotropic axisymmetric equilibrium state of lowest energy that the compressible star can reach while conserving its radiation zone mass (), the mass in each flux tube (), the flux of the toroidal field (), and a given helicity ():

 1μ0Δ5Φ=aI;1¯¯¯ρ+aII;1T, (37)
 1μ0T=−aII;1Φ−aII;0s2. (38)

Since the azimuthal field is regular at the origin, we get from Eq. (23) . Eliminating between Eqs. (37-38), we obtain

 Δ5Φ+[μ0aII;1]2Φ=μ0aI;1¯¯¯ρ (39)

that becomes, multiplying it by and using Eq. (24) & (25)

 Δ∗Ψ+[μ0aII;1]2Ψ=μ0aI;1¯¯¯ρr2sin2θ. (40)

We thus identify in Eq. (20)

 ⎧⎪⎨⎪⎩k=1j=0=and⎧⎪ ⎪⎨⎪ ⎪⎩aI;1=−β0aII;1=−1μ0λ1/R=, (41)

where we have constrained the initial arbitrary functions of the magnetohydrostatic equilibrium

 F(Ψ)=−μ0aII;1ΨandG(Ψ)=−aI;1. (42)

So, it reduces to

 Δ∗Ψ+λ21R2Ψ=−μ0¯¯¯ρr2sin2θβ0, (43)

the values of the real coefficients and being thus controled by the helicity () and the mass conservation in each axisymmetric flux tube defined by because of the non force-free stably stratified behaviour of the reached equilibrium.

This corresponds, as has already been emphasized, to the lowest energy equilibrium state for a given helicity (Bellan 2000; Broderick & Narayan 2008). The equilibrium state ruled by Eq. (43) are thus the generalization of the Taylor relaxation states in a self-gravitating star where the field is not force-free (i.e. ). Note also that some non force-free relaxed states have been identified in plasma physics (Montgomery & Phillips 1988, 1989; Dasgupta et al. 2002; Shaikh et al. 2008) that should be studied in a stellar context in a near future.

Let us note that in the case where is not considered () we recover the Chandrasekhar (1956a) force-free limit (see also Marsh 1992, for a generalization of the solutions) and the usual Taylor’s states for low- plasmas. The Prendergast’s model is recovered assuming a constant density profile (incompressible).

### 3.3 Green’s function solution

We are now ready to solve Eq. (43). If we introduce and if we set , it is recast in

 Lλ1Ψ=S, (44)

where

 Lλ1≡[∂rr+1−x2r2∂xx+λ21R2]. (45)

Using Green’s function method (Morse & Feshbach 1953) , we then obtain the particular solution associated with :

 Ψp(r,x)=−μ0β0∑lλl1Rsup[2l+32(l+1)(l+2)]× (46) {rjl+1(λl1rRsup)∫Rsupr[ξyl+1(λl1ξRsup)Jl(ξ)]dξ +ryl+1(λl1rRsup)∫rRinf[ξjl+1(λl1ξRsup)Jl(ξ)]dξ} ×(1−x2)C3/2l(x),

where

 Jl(ξ)=∫1−1S(ξ,x′)C3/2l(x′)dx′; (47)

and are respectively the spherical Bessel functions of the first and the second kind (also called Neumann functions) while are the Gegenbauer polynomials (Abramowitz & Stegun 1972). and , which are respectively the bottom and the top radius of the considered radiation zone, are introduced and we identify .

These functions (, , and ) are respectively the radial and the latitudinal eigenfunctions of the homogeneous equation associated with Eq. (44) :

 Lλ1Ψh=0. (48)

Then, if we express the solutions of this equation as , we get respectively

 (1−x2)d2gldx2+(l+1)(l+2)gl=0 (49)

and

 d2fldr2+⎡⎣(λl1Rsup)2−(l+1)(l+2)r2⎤⎦fl=0, (50)

giving

 gl=(1−x2)C3/2l(x) (51)

and

 fl = Kl1λl1rRsupjl+1(λl1rRsup) (52) +Kl2λl1rRsupyl+1(λl1rRsup),

and being real constants. are eigenvalues that allow to verify the boundary conditions discussed hereafter. One has to notice that has to vanish in order to preserve the regularity of the solution at the center. Applying this to Eq. (43), we finally obtain the general solution

 Ψ(r,θ)=Ψh+Ψp = sin2θ×{∞∑l=0Kl1λl,i1Rsuprjl+1(λl,i1rRsup)C3/2l(cosθ) − μ0β0λ0,i1Rsuprj1(λ0,i1rRsup)∫Rsupr[y1(λ0,i1ξRsup)¯¯¯ρξ3]dξ − μ0β0λ0,i1Rsupry1(λ0,i1rRsup)∫rRinf[j1(λ0,i1ξRsup)¯¯¯ρξ3]dξ}.

Notice that the particular solution for the poloidal flux function presents a dipolar geometry, owing to its angular dependence that follows the one from the source term .
Neglecting the density, we end up with the linear homogeneous equation, whose solutions are Chandrasekhar-Kendall functions (Chandrasekhar & Kendall 1957). Moreover, the source term being constituted only of a dipolar component, all the non-dipolar contributions are, according to this model, force-free.

The magnetic field is then given for by

 {B}=1r2sinθ∂θΨ^er−1rsinθ∂rΨ^eθ+λ0,i1RsupΨrsinθ^eφ. (54)

After a few manipulations, we can then express the current density as

 {j}P =1μ0∇×{B}T =λ0,i1μ0R{B}Pforce−free, (55) {j}T =1μ0∇×{B}P =λ0,i1μ0R{B}Tforce−free+β0Ê¯¯¯ρrsinθ^eφnonforce−free, (56)

where we recognize in the first term of the right hand side the force-free contributions and in the second the non force-free one, fully contained in the toroidal component.

The Lorentz force can as a matter of fact be written in the very simple form:

 {F}L={F}LP=β0Ê¯¯¯ρ∇Ψ. (57)

### 3.4 Configurations

The boundary conditions for which determine possible values for and have now to be discussed. Two major types of geometry are relevant for large-scale fossil magnetic fields in stellar interiors: initially confined and open configurations.

#### 3.4.1 Initially confined configurations

Let us first concentrate on the simplest mathematical solution in the case of a central radiation zone that initially cancels both at the center () and at a given confinement radius ().

Then, if we choose to cancel the coefficients for every , the condition is verified since . However, if we look at the magnetic field radial component behaviour at the center, it is easily shown, using Eq. (54), that if it is given by which does not cancel so that (where ). Therefore, this solution is multivaluated, thus physically inadmissible and .

Then, we consider the general case of a field initially confined between two radii and , owing to the presence of both a convective core and a convective envelope (as it is the case e.g. in A-type stars). We impose and , which gives the two independent equations for

 K01=μ0β0∫Rc2Rc1[y1(λ0,i1ξRc2)¯¯¯ρξ3]dξ (58)

and

 K01j1(λ0,i1)=μ0β0y1(λ0,i1)∫Rc2Rc1[j1(λ0,i1ξRc2)¯¯¯ρξ3]dξ; (59)

we here focus on the dipolar mode which is known to be the lowest energy per helicity ratio state (cf. Broderick & Narayan 2008). These can be formulated so that one first determines the value of according to

 j1(λ0,i1)∫Rc2Rc1[y1(λ0,i1ξRc2)¯ρξ3]dξ −y1(λ0,i1)∫Rc2Rc1[j1(λ0,i1ξRc2)¯ρξ3]dξ=0 (60)

and next computes following (59).

In the case where there is no convective core, as for example in central radiation zones of late-type stars such as the Sun, Eqs. (59) & (60) must be applied setting .

#### 3.4.2 Open configurations

This corresponds to the case of fields that match at the stellar surface (at , being the star’s radius) with a potential field as observed now in some cases of early-type stars such as Ap stars. Then, we have , being the associated potential.

In the case studied here, we focus on the first configuration (initially confined) since the search of relaxed solutions for given , , , and assumes that on the stellar radiation zones boundaries. This initial confined configuration will then become open one through Ohmic diffusion as in the Braithwaite and collaborators scenario.

## 4 Application to Realistic Stellar Interiors

To illustrate our purpose, we apply our analytical results (i) to model an initial fossil field buried below the convective envelope of the young Sun on the ZAMS, and then (ii) to model an initial field present in the radiation zone of a ZAMS magnetic Ap-star, whose lower and upper radiation-convection interfaces are respectively located at and at . In the first case, the parameter is determined such that the maximum field strength reaches the amplitude of , which is one of the upper limits given by Friedland & Gruzinov (2004) for the present Sun’s radiative core. In the second case, it is obtained such that it reaches the arbitrary value of . This value is approximatively of the same order of magnitude that the mean surface amplitude observed using spectropolarimetry for magnetic Ap-star which exhibits strong external dipolar magnetic behaviour (such as HD12288 (Wade et al. 2000)). We thus assume implicitly that such an initial confined internal field is a potential prelude to the multipolar one observed now at the surface, the latter state being acheived after a diffusive process that will be studied in a forthcoming paper.

### 4.1 Fossil fields buried in late-type stars radiative cores

The young Sun model used as a reference is a Cesam non-rotating standard one (Morel 1997), following inputs from the work of Couvidat et al. (2003) and Turck-Chièze et al. (2004).

In Fig. 2, three possible configurations for are given. We choose those corresponding to the first, the third and the fifth eigenvalues (given in table 1). Those are the generalization of the well-known Grad-Shafranov equation linear eigenmodes obtained in the force-free case (cf. Marsh 1992). The field is then of mixed-type ( is given for ), both poloidal and toroidal, and non force-free, properties already obtained by Prendergast (1956) in the incompressible case. The respective amplitudes ratio between the poloidal and the toroidal components will be described in §5., where the possible stability of such configurations will be discussed.

### 4.2 Fossil fields in early-type stars

Respective corresponding possible configurations in the case of an Ap star are given in Fig. 3.

The model is typical of an A2p-type star, with an initial mass . The solar metallicity is chosen as the initial one and the model is taken on the ZAMS, its luminosity being .

Obtained configurations are then mixed poloidal-toroidal (twisted) fields which may be stable in stellar radiation zones (cf. Braithwaite & Spruit 2004; Braithwaite & Nordlund 2006). Their configurations are thus given in both cases by concentric torus, the neutral poloidal points (where so that ) positions being function of the internal density profile of the star.

Let us emphasize here that the original approach of this work first consists in deriving the Grad-Shafranov-like equation adapted to treat the barotropic magnetohydrostatic equilibrium states for realistic models of stellar interiors. Such approach has already been applied to investigate the internal magnetic configurations in polytropes and in compact objects such that white dwarfs or neutron stars (see e.g. Monaghan 1976; Payne & Melatos 2004; Tomimura & Eriguchi 2005; Yoshida et al. 2006; Haskell et al. 2008; Akgün & Wasserman 2008; Kiuchi & Kotake 2008).

Then, the obtained arbitrary functions are constrained with deriving minimal energy equilibrium configurations for a given conserved mass, azimuthal flux and helicity that generalizes the relaxation Taylor’s states to the self-gravitating case where the field is non-force free.

## 5 Links between the field’s helicity, topology and energy

### 5.1 Helicity vs. mixity

Let us express the magnetic field in terms of magnetic stream functions (for the poloidal component of the field) and (for its toroidal part):

 B=∇×[∇×[ξP(r,θ)ˆer]+ξT(r,θ)ˆer]. (61)

Next, the vector potential A is given by the relation . Knowing that in the confined case the gauge choice is inconsequential, we can identify without further ado

 {A}=∇×[ξP(r,θ)ˆer]+ξT(r,θ)ˆer. (62)

The magnetic stream functions are then projected on the spherical harmonics

 ξP(r,θ) = ∑ℓ>0ξℓ0(r)Y0ℓ(θ), (63) ξT(r,θ) = ∑ℓ>0χℓ0(r)Y0ℓ(θ). (64)

From now on, we use Einstein summation convention where and the vectorial spherical harmonics basis such that any axisymmetric vector field can be expanded as

 \@vecu(r,θ)=uℓ0(r)R0ℓ(θ)+vl0(r)S0ℓ(θ)+wℓ0(r)T0ℓ(θ), (65)

where the vectorial spherical harmonics , , and are defined by:

 R0ℓ(θ)=Y0ℓ(θ)ˆer, S0ℓ(θ)=∇SY0ℓ(θ), T0ℓ(θ)=∇S×R0ℓ(θ), (66)

the horizontal gradient being defined as . (cf. Rieutord (1987)). Since

 ∇×(ξPˆer)=∇×(ξℓ0Y0ℓˆer)=∇×(ξℓ0R0ℓ)=ξℓ0rT0ℓ, (67)

we get from Eq. (62)

 {A}=χℓ0R0ℓ+ξℓ0rT0ℓ. (68)

On the other hand, we have

 B=ℓ(ℓ+1)r2ξℓ0R0ℓ+1r∂rξℓ0S0ℓ+χℓ0rT0ℓ, (69)

from which we finally obtain the expression for the helicity

 H = ∫R∗0∫Ω{[χℓ0R0ℓ+ξℓ0rT0ℓ]⋅ (70) [ℓ′(ℓ′+1)