Issue 
A&A
Volume 533, September 2011



Article Number  A125  
Number of page(s)  16  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201117105  
Published online  13 September 2011 
Forcefree twisted magnetospheres of neutron stars
Departament de Física AplicadaUniversitat d’Alacant,
Ap. Correus 99,
03080
Alacant,
Spain
email: daniele.vigano@ua.es
Received:
19
April
2011
Accepted:
8
July
2011
Context. The Xray spectra observed in the persistent emission of magnetars are evidence for the existence of a magnetosphere. The highenergy part of the spectra is explained by resonant cyclotron upscattering of soft thermal photons in a twisted magnetosphere, which has motivated an increasing number of efforts to improve and generalize existing magnetosphere models.
Aims. We want to build more general configurations of twisted, forcefree magnetospheres as a first step to understanding the role played by the magnetic field geometry in the observed spectra.
Methods. First we reviewed and extended previous analytical works to assess the viability and limitations of semianalytical approaches. Second, we built a numerical code able to relax an initial configuration of a nonrotating magnetosphere to a forcefree geometry, provided any arbitrary form of the magnetic field at the star surface. The numerical code is based on a finitedifference timedomain, divergencefree, and conservative scheme, based of the magnetofrictional method used in other scenarios.
Results. We obtain new numerical configurations of twisted magnetospheres, with distributions of twist and currents that differ from previous analytical solutions. The range of global twist of the new family of solutions is similar to the existing semianalytical models (up to some radians), but the achieved geometry may be quite different.
Conclusions. The geometry of twisted, forcefree magnetospheres shows a wider variety of possibilities than previously considered. This has implications for the observed spectra and opens the possibility of implementing alternative models in simulations of radiative transfer aiming at providing spectra to be compared with observations.
Key words: stars: magnetic field / stars: neutron / magnetohydrodynamics (MHD) / scattering
© ESO, 2011
1. Introduction
Anomalous Xray pulsars (AXPs) and soft gammaray repeaters (SGRs) are a class of neutron stars (NSs) characterized by high Xray quiescent luminosities, short Xray bursts, and giant flares (for SGRs). They are believed to be magnetars (Duncan & Thompson 1992; Thompson & Duncan 1996): young (P/2Ṗ ~ 10^{2}−10^{5} yr) NSs with very strong magnetic fields ~10^{14}−10^{15} G. Magnetars are supposedly born with a very short spin period (few milliseconds), which causes the formation of an intense inner magnetic field via dynamo processes and a subsequent rapid loss of rotational energy by means of magnetic braking. As more and better data become available in the last decade^{1}, the separation between the two traditional classes (SGRs and AXPs) has become thinner, and they are no longer considered different classes. All SGRs and most of AXPs have shown sporadic bursts in the X and γ bands, and the recent discovery of transient magnetars has made the situation even less clear. Currently, our interpretation of the term magnetars is being reconsidered (see Mereghetti 2008; Rea & Esposito 2011 for observational reviews). In general, magnetars are characterized by long periods P ~ 2−12 s and high Xray persistent luminosities L_{X} ~ 2 × 10^{33}−2 × 10^{35} erg s^{1}, which point to the magnetic field decay as the main source of energy, unlike radio pulsars, fed by the loss of rotational energy. One key ingredient in those models that try to explain the variety of magnetar phenomenology and their spectra is the presence of toroidal fields, of the same order of magnitude or even larger than the poloidal component, both in the interior and in the external magnetosphere.
The object of this paper is the modeling of the magnetosphere of magnetars, which is crucial because it is one of the main ingredients in the formation of spectra. The difficulty in building fully consistent solutions of NS magnetospheres is such that, after more than 40 years, the aligned rotator model proposed by Goldreich & Julian (1969) still remains very popular for explaining the basic features of a rotating magnetosphere. We must point out that the aligned rotator model, with all its variants, is surely a naive description of pulsars that has been strongly criticized (e.g., Michel 1980), among other reasons, simply because it cannot explain why pulsars pulsate. A further refinement consists of the ideal MHD approximation for forcefree magnetospheres, which leads to the socalled pulsar equation. This is simply the balance of the electromagnetic forces around the rotating neutron star under the assumption of axial symmetry and it neglects the pressure, inertial, and gravitational terms. The rotationally induced electric field is required to be perpendicular to the magnetic field. The latter can in general have an azimuthal component that twists the poloidal field lines. This toroidal component has to be sustained by currents, whose shape is the source term of the pulsar equation: it is the fundamental parameter to be chosen in order to find the corresponding geometric configuration of the magnetic field.
The problem is analytically solvable only for a few simple, adhoc choices of the current, such as in Michel (1973b), but these solutions generally present unphysical features in some border regions. An alternative numerical approach by Contopoulos et al. (1999) is useful to construct consistent forcefree models, with a smooth matching across the light cylinder. In this model the corotating region only hosts an untwisted dipolar magnetic field, and the focus is on the bundle of open twisted lines that cross the light cylinder and the resulting output power. The model was refined by Gruzinov (2005) and extended to the timedependent, misaligned case by means of 3D simulations by Spitkovsky (2006) and Kalapotharakos & Contopoulos (2009).
Nowadays, Xray instruments allow us to obtain more and more information from the spectra of magnetars. An especially interesting feature in the persistent emission of all of the magnetar candidates is that their spectra can be well fitted with a thermal component (0.4−0.7 keV) plus a hard nonthermal tail, described by a power law with photon index β ~ 3−4 (Mereghetti 2008). This hard component is commonly attributed to resonant Compton scattering, which becomes important in the presence of a strong current (i.e., high charge density). For this reason, much attention has been given to the modeling of currents flowing in a twisted magnetosphere and obtaining estimates of the charge carrier density and the corresponding resonant optical depth. Lyutikov & Gavriil (2006) propose a semianalytical 1D model able to describe the basic properties of radiation transfer across a twisted magnetosphere. More accurate 3D Monte Carlo simulations have been performed for the nonrelativistic (Fernández & Thompson 2007; Nobili et al. 2008a) and relativistic cases (Nobili et al. 2008b). Fernández & Davis (2011) discuss the implications of the outgoing polarization. The implementation in XSpec and the systematic fits to observational data have been successfully performed by Rea et al. (2008) and Zane et al. (2009) for the 1D and 3D cases, respectively. All these very remarkable results rely on the same family of models: the selfsimilar twisted dipole proposed by Thompson et al. (2002), extended to higher multipoles by Pavan et al. (2009), which represents a semianalytical, easytoimplement solution to the pulsar equation. However, this model lacks generality, as it relies on a very particular choice of the form of the current.
We note that forcefree magnetosphere models rely on a purely macrophysical approach but a fully coherent microphysical description is still lacking: how are the currents and electric fields in a given model generated and sustained? Which particles are present (electrons, ions, pairs)? Where do they come from? Which velocity distribution can be assigned to each kind of particle? In the magnetar framework, there are some excellent attempts to give a microphysical description of coronae. Beloborodov & Thompson (2007) deal with the problem considering a deviation from forcefree equilibrium, which has to be introduced to explain the production of pairs, but many questions remain open. Given the complexity of the problem, one has to decide whether to focus on detailed microphysical processes of an approximate macrophysical solution or to improve the largescale electrodynamics solutions by gradually removing some of the simplifying hypotheses.
In this article, we present our efforts to obtain some more general configurations of a twisted magnetosphere and to estimate the dependence of the spectrum on the geometric configuration. In Sect. 2, we set the formalism needed to describe the problem and present some analytical and semianalytical generalizations of previous works. In Sect. 3 we describe the numerical code used to build forcefree configurations and the battery of tests performed to verify the code. We pay special attention to the effect of key input parameters and the convergence of the method. In Sect. 4 we present new results for magnetosphere models, and explore the sensitivity of the relevant output (current, charge distribution, and resonant optical depth) to different input parameters. In Sect. 5 we summarize our findings and outline future improvements.
2. Analytical and semianalytical solutions
Fig. 1 Relation between magnetic flux function Γ, enclosed current function I(Γ), and toroidal field. 
2.1. The pulsar equation
We consider an NS with a magnetic field B and angular velocity Ω. Assuming that the plasma is a perfect conductor in vicinity of the star, rotation induces an electric field given by (1)where c is the velocity of light. Neglecting centrifugal (i.e., inertial), collisional, and gravitational terms compared to electromagnetic forces, the force equilibrium equation is expressed by (2)where ρ_{e} is the charge density and J the current density. Assuming axial symmetry (aligned rotator) and working in spherical coordinates (r,θ,φ), we introduce the magnetic flux function Γ(r,θ), in terms of which the poloidal field is (3)where is the unit azimuthal vector. With the choice of gauge Γ = 0 on the axis, the function Γ(r,θ) is related to the φcomponent of the potential vector by (4)and it is constant along a field line (B_{p}·∇Γ = 0). Thus its value labels the axisymmetric surface S_{Γ} given by the azimuthal rotation of one field line. The magnetic flux flowing within the surface S_{Γ}, with value Φ_{Γ}(B) = 2πΓ, is conserved by definition. A toroidal component (i.e., φcomponent in axial symmetry) of the magnetic field is present if sustained by poloidal currents. The enclosed current I flowing within S_{Γ} (see blue arrows in Fig. 1 and note that the forcefree condition implies that currents flow parallel to the surfaces S_{Γ}) has to be a function of only Γ, with the only requirement being that I(0) = 0. Thus we have the freedom of choosing the enclosed current function I(Γ), and by Ampère’s law, the φcomponent of the magnetic field is related to I by (5)For a given magnetosphere model, it is useful to define the twist of a field line as the azimuthal displacement between its surface footprints (for a dipolarlike configuration, one footprint in each hemisphere). The integral of the local displacement δφ(r,θ) = B_{φ}δθ/B_{θ}sinθ along the magnetic field line, l_{Γ}, between its surface footprints at magnetic colatitudes θ_{1,2}(Γ) is the line twist (6)where the dependence r(θ,Γ) can be found by solving the field line equations. The function Δφ_{tw}(Γ) is used to characterize a magnetosphere.
We then return to the equilibrium Eq. (2). In axial symmetry, only poloidal components of the electric field are allowed. Considering the poloidal component of Eq. (2) in cylindrical coordinates (ρ,φ,z) and denoting the radius of the light cylinder by r_{l} = c/Ω, we obtain the socalled pulsar equation (Scharlemann & Wagoner 1973; Michel 1973b): (7)where we indicate the partial derivatives of Γ with a subscript. The terms (ρ/r_{l})^{2} arise from the Coulomb force ρ_{e}E. The pulsar equation is actually the GradShafranov equation for forcefree fields. In general, Eq. (7) has to be solved numerically. Even for trivial choices of untwisted configurations (I(Γ) = 0), the corotation velocity of the charged particles distorts the vacuum solutions, resulting in open field lines near the light cylinder ρ ~ O(r_{l}) (Michel 1973a). The only fully analytical solution including rotation and a smooth matching across the light cylinder is the split twisted monopole presented by Michel (1973b), in which the field lines are radial and the magnetic field changes sign between the northern and southern hemispheres. A current sheet at the equator is needed to ensure the divergencefree condition. This is a windlike solution inappropriate to describing the region with closed field lines. Later, the numerical solution by Contopoulos et al. (1999) for the first time provided a complete description of a magnetosphere, whose closed lines are untwisted, while the wind zone has a configuration that is very similar to the split monopole.
In the case of a nonrotating star (magnetars are slow rotators and this may be a reasonable approximation for some problems), no electric field is induced and the equilibrium equation is simply J × B = 0: the forcefree condition only requires that currents flow parallel to magnetic field lines. In this limit, we can reformulate the whole problem with the simple equation (8)where the function α(Γ) is related to the enclosed current function introduced in Eq. (5) by (9)The pulsar equation in spherical coordinates is now expressed as (10)Next we expand the magnetic flux function in Legendre polynomials, P_{l}(μ): (11)where μ = cosθ, a_{l}(r) is the dimensionless radial function, r_{ ⋆ } the radius of the star, and Γ_{0} the magnetic flux normalization. Denoting the strength at pole by B_{0}, we hereafter choose (12)This leads to the following expression for the poloidal magnetic field components: (13)We can obtain the governing differential equation from Eq. (10): (14)Finally, by using the orthogonality relations of Legendre polynomials, we have (15)The righthand side depends on the functional form of α(Γ) and is responsible for the coupling between different multipoles. For constant α the equation has analytical solutions (see Sect. 2.2), while the problem is more difficult for other choices. The rest of this section is devoted to reviewing some possible choices that make the problem (semi)analytically solvable in the nonrotating case, which applies to magnetars.
2.2. Constant α: Bessel solutions
The choice of constant α = k/r_{ ⋆ }, where k is a dimensionless parameter, leads to decoupled equations for each single multipole, so that for each l, we must independently solve the equation (16)The analytic solutions of this equation are the spherical Bessel functions of the first and second kinds: where x = kr/r_{ ⋆ }. The constant ratio r_{ ⋆ }/k gives the typical scalelength of the magnetic field. The physical solutions (from which the vacuum solution B_{p} ∝ r^{−(l + 1)} can be recovered in the limit k → 0) are represented by the functions of the second kind, explicitly: (19)where c_{l} is the weight of the lmultipole.
This choice of constant α has been followed in several studies of the solar corona (Chiu & Hilton 1977; Seehafer 1978) or applied to the region with open field lines of the pulsar magnetosphere (Scharlemann & Wagoner 1973). The Y_{l}(x) functions are oscillatory at large distances, and the magnetic field components change sign as r varies (at fixed parameter k). The problem with this family of analytical solutions is that, at large distances, all components (which have the same radial dependence for any l) decay too slowly: B_{r} → r^{2},B_{θ} → r^{1},B_{φ} → r^{1}. Thus this configuration cannot be a solution for the whole space, as it would imply infinite magnetic energy in an infinite volume. Also, these solutions cannot be continuously matched with vacuum, because it would require that, at the same radius r_{out}, B_{φ}(r_{out}) = 0 and B_{r}(r_{out}) ≠ 0, a condition that cannot be satisfied because those two components are both proportional to the Y_{l}(x) functions and have the same zeros.
2.3. Selfsimilar models
Low & Lou (1990), followed by Wolfson (1995) and other authors, studied a particular class of selfsimilar solutions to describe the opening of the solar coronal magnetic field lines due to their shear. Thompson et al. (2002) applied the same approach in the magnetar framework. Assuming a radial powerlaw form for the magnetic flux function and a radial dependence α ∝ 1/r, the enclosed current function is (22)where (23)is the current normalization, related to the amount of current, like the dimensionless parameter k_{ss}. Equation (10) becomes a nonlinear secondorder differential equation for the angular function F(μ): (24)where primes stand for derivatives of F(μ) with respect to μ. The magnetic field is given by (25)The three components of the magnetic field have the same radial dependence, so we are seeking for selfsimilar solutions. The maximum line twist, which is the one suffered by a polar line (Γ → 0), is called global twist, and it can be taken as the parameter that uniquely describes the lpole family of selfsimilar solutions. Due to selfsimilarity, Eq. (6) for Γ → 0 becomes very straightforward: (26)Equation (24) has to be numerically solved by imposing three physical requirements: F(1) = F(−1) = 0 (on the axis, only the radial field component is allowed), and F′(1) = −2 (to fix the normalization). Given a value of k_{ss}, one can find an infinite number of solutions characterized by an eigenvalue p. Each solution represents a different multipole (except the dipole for which there are two solutions). Fixing k_{ss} = 0 (i.e., I_{0} = 0, no current) implies an integer value of p, so that we recover the vacuum multipolar solution with F(μ) = (1 − μ^{2})dP_{p}(μ)/dμ where P_{p} is the pth Legendre polynomial. For each multipole, a unique relation k_{ss}(p) (or I_{0}(p) once the values of B_{0} and r_{ ⋆ } are fixed) defines the family of solutions. In Fig. 2 we show the curve of parameters we obtained for the twisted dipole family. Our results agree with the dipolar solutions of Thompson et al. (2002) and the higher multipole solutions of Pavan et al. (2009). We have used these models in the numerical code (section Sect. 3) for testing purposes and for the sake of comparison with other numerical solutions.
Fig. 2 Curve I_{0}(p) describing the family of selfsimilar twisted dipoles. 
However, all these solutions are of limited generality because they have a defined symmetry with respect to the equator, and a linear combination of solutions for different multipoles is not a new solution, due to the nonlinear character of the problem.
2.4. General axisymmetric solutions
The righthand side depends on the functional form of α(Γ) and is responsible for the coupling between different multipoles. We have explored several possible analytical forms of α(Γ). Aiming at more general solutions, we explore several possible analytical forms of α(Γ) in Eq. (15). In search of semianalytical solutions, we consider power laws of the form (27)with k the dimensionless parameter regulating the amount of currents, and Γ_{0} given by Eq. (12).
Obviously the trivial choice α = 0 gives the vacuum (i.e., untwisted) solutions. The choice of α constant has already been considered in Sect. 2.2. Some choices with q < 0, such as q = −1/2, present mathematical problems and are definitely not physical, and high values of q make the problem unsolvable with semianalytical methods. Therefore, only a few simple choices may be useful because of their simplicity. An interesting case is q = 1/2. Assuming that Γ > 0 in the domain of integration, this particular choice leads to (28)and the following ODE for each radial function f_{l} = a_{l}(r)r/r_{ ⋆ }: (29)where each dimensionless Gaunt factor g_{lmn} involves the integral of the product of three Legendre polynomials (see Appendix B for a detailed derivation). Numerically it is ~O(1). Therefore, the source term couples each multipole l with an infinite number of multipoles m and n. One can solve the problem by truncating the series and allowing only for a finite number of multipoles. This approach is more general than other simple cases, but it has many more degrees of freedom.
3. Numerical solutions
In the previous section we have discussed some analytical and semianalytical solutions that share the same drawback: the arbitrary choice of the enclosed current function I(Γ) or, equivalently, α(Γ). Some of these solutions are nonphysical, in the sense that they can neither be extended to infinity nor matched to vacuum solutions. All these limitations make the (semi) analytical approach insufficient for general purposes, because we have no physical argument for preferring one particular form of the current over another. The alternative is to find numerical solutions of the nonlinear, forcefree equations describing an NS magnetosphere. We have three main reasons for working in this direction. First, we expect these solutions to be more general and in some cases very different from the semianalytical ones. By studying new numerical equilibrium solutions, we hope to gain insight into the form of the most realistic choices of enclosed current functions. Second, a reliable numerical code could in principle be extended, if adding the rotationallyinduced electric field, to build a consistent corotating magnetosphere model from which the distribution of charge density can be computed directly. Third, the ability to build magnetospheric solutions for any prescribed form of the magnetic field at the surface is very useful for future studies of the dynamics of the magnetosphere coupled with the physics of the crust.
3.1. The magnetofrictional method
For simplicity, we begin with the twodimensional nonrotating case, assuming that the forcefree approximation is valid up to an outer radius r_{out}, which reduces the problem to finding solutions of Eq. (10). In the socalled magnetofrictional method (Yang et al. 1986; Roumeliotis et al. 1994), one begins with an initially nonforcefree configuration and defines a fictitious velocity field proportional to the Lorentz force, v_{f} = νJ × B/B^{2}, where .ν is a normalization constant that sets the time unit in the induction equation. Hereafter we use ν = 1. This results in a fictitious electric field: (30)This fictitious electric field enters into the induction equation as a frictional term that forces the solution to relax to a forcefree configuration.
In the original method, Roumeliotis et al. (1994) write the magnetic field as (31)where Γ is the magnetic flux function defined in Eq. (3), and β = φ − γ(r,θ), where γ(r,θ) is a scalar potential related with the freedom in the choice of the form of the enclosed current. The induction equation becomes a system of two advection equations for those two functions: (32)A static solution is achieved if and only if v_{f} = 0, since the velocity field, ∇Γ, and ∇β are all orthogonal to the magnetic field by definition.
We apply the same idea but, instead of evolving the functions Γ and β, we evolve the magnetic field components directly by solving the induction equation (33)where the fictitious electric field can be written as (34)where E_{f} is a measure of the deviation from the forcefree condition, because J ∥ B is accomplished if and only if E_{f} ≡ 0. The main reason for solving the induction equation instead of Eqs. (32) is to allow for future extensions of the code by considering a real, rotationallyinduced electric field. The disadvantage is that we have to be more careful when setting boundary conditions for the electric field, because we could converge to stationary solutions characterized by ∇ × E_{f} = 0, which are not necessarily forcefree.
3.2. Linear analysis of the magnetofrictional method
We now consider a background, uniform magnetic field B_{0} and a small perturbation δB ∝ e^{i(k·r − ωt)}. In the linear regime, the equations read as (35)Explicitly, the last equation can be written as (36)If the perturbed current is orthogonal to the background magnetic field (either longitudinal perturbations δB ∥ B_{0} or transverse perturbations with k ∥ B_{0}), the first term on the righthand side of Eq. (36) vanishes and the dispersion relation is purely dissipative: (37)Any perturbation of this type will be dissipated on a timescale ∝ k^{2}. In contrast, for transverse perturbations with both δB and k orthogonal to B_{0}, the current is parallel to the background magnetic field, and the two terms in the righthand side of Eq. (36) cancel out, so that the perturbation does not evolve (a neutral mode with ω = 0). Therefore, the magnetofrictional method is designed to dissipate all induced currents nonparallel to the magnetic field but allows for stationary solutions with currents parallel to the field. Since the largest lengthscale in our problem is set by the size of the numerical domain, λ_{max}πr_{out}, the typical diffusion timescale on which we expect to converge to a forcefree solution is .
3.3. The numerical method
We work in spherical coordinates (r,θ,φ) under the assumption of axial symmetry. We employ a fully explicit finite difference time domain (FDTD) method (Taflove & Brodwin 1975) with a numerical grid equally spaced in θ and logarithmic in the radial direction, unless the outer radius is very close to the surface of the star (located at r = 1), in which case we employ a linearly spaced radial grid. Our typical resolution varies between 30–200 points in the radial direction and 30–100 points in the angular direction. At each node (θ_{i},r_{j}), we define all components of B^{(i,j)} and consider the three surfaces centered on this node and normal to the unit vectors . The elementary areas are calculated by the following exact formulas: Since the righthand side of Eq. (33) contains ∇ × E_{f}, it is useful to apply Stokes’ theorem as follows. The magnetic flux across is approximated by (41)and the magnetic fluxes are advanced in time using (42)where the numerical circulation of the electric field along the edges of the surface S_{n} is approximated by using the values of E_{f} in the middle of the edges, whose lengths are Explicitly, the circulation of E_{f} along the edge of each surface is In Fig. 3 we show the location of the variables needed for the time advance of , (red), and (green).
We also make use of Stokes’ theorem to calculate the current components at each node:
Then, the values of J^{(i,j)} and B^{(i,j)} directly provide , as defined by Eq. (34). For each cell (i,j) the local divergence of B is defined as the net magnetic flux flowing across the surfaces divided by the cell volume (fluxes through the toroidal surfaces S_{φ} cancel due to axial symmetry): (49)with (50)The numerical method ensures that the local divergence is conserved to machine accuracy by construction. As a matter of fact, Eqs. (42) and (49) imply (51)According to Eq. (46), the last equation is written as a sum of toroidal elements E_{φ}l_{φ}. They are evaluated twice at each of the four surrounding edges (with center θ_{i ± 1},r_{j ± 1}), with opposite sign, so they cancel exactly. Therefore, provided there is an initial divergenceless field, it is guaranteed that the numerical divergence defined by Eq. (49) remains zero to roundoff level at each time step.
We must mention that we also tried a method with a staggered grid (Yee 1966), in which each ncomponent of the magnetic field is defined only at the center of the normal surface, , while the electric field components are defined in the middle of its delimiting edges. Methods based on staggered grid are wellsuited to solving Maxwell’s equations, because it provides a natural way to timeadvance one field by means of the circulation of the other one. However, we are not dealing with the true Maxwell’s equations, but rather with an artificial electric field, Eq. (34). Evaluating the dot product J·B requires the interpolation between two or four values of several of the six mutually displaced components. Considering the red components in Fig. 3, for instance, the calculation of also requires , which would be not defined at that location. The unavoidable interpolation errors prevent the code from completely relax to E_{f} = 0, except in the trivial case of untwisted configurations. For this reason, we decided to work with a standard grid.
Fig. 3 Location of the variables in the numerical grid. 
3.4. Boundary conditions
At the polar axis, we impose the vanishing of all angular components of magnetic field and currents: B_{θ} = B_{φ} = J_{θ} = J_{φ} = E_{θ} = E_{φ} = 0. At the surface, we have to fix the magnetic field components as provided by some interior solution. However, an arbitrary choice of poloidal and toroidal fields may not be compatible with a forcefree solution. We decided to impose E_{f} = 0, at the surface, which is equivalent to keeping the value of the radial component fixed at the surface, B_{r}(1,θ) and therefore to fixing the angular dependence of the magnetic flux function, Γ(1,θ). As a consequence, B_{θ} and B_{φ} are allowed to vary on the first radial grid point.
The external boundary is set at r = r_{out}. We have explored two different boundary conditions: E_{f}(r ≥ r_{out}) = 0 and the continuous matching to external vacuum solutions. The first choice is equivalent to fix the radial component B_{r}(r_{out},θ), while allowing for B_{θ} and B_{φ} to evolve. Coupling to vacuum solutions can be done using the spectral Legendre decomposition of the radial field at the outer surface, or more precisely, of the magnetic flux function, in terms of which the vacuum boundary condition for each multipole is easily imposed (e.g., as in Pons & Geppert 2007). Then we can reconstruct the angular dependence of B_{θ}(r_{out},θ). The vacuum region is characterized by B_{φ} = 0 and the absence of currents or fictitious electric fields. This implies that a current sheet J_{θ}(r_{out}) ≠ 0 is needed to ensure current conservation.
If we choose E_{f} = 0 as outer boundary condition, the code can actually converge to E_{f} ≡ 0 at a roundoff level, because mathematically this is the only solution compatible with ∇ × E_{f} = 0. The price to pay is a forced matching of the inner solution with a fixed value of B_{r}(r_{out}). In contrast, if we couple with vacuum, there is no guarantee that the final solution is E_{f} = 0 everywhere. We discuss below the influence of the different boundary conditions on the results.
3.5. Convergence criterion and tests
Since the magnetofrictional method is based on introducing an artificial, viscous electric field that drives an arbitrary initial configuration into a forcefree state, we need a convergence criterion to decide when our solutions are acceptable. For that purpose, we keep track during the run of the following quantities:

Volumeintegrated magnetic energy (total and contributionfrom the toroidal field)

Volumeintegrated energy stored in the fictitious electric field

Total volumeintegrated helicity (see Appendix A for a discussion about this definition).

Volumeintegrated absolute value of the divergence of both, B and J. These two quantities are expected to vanish at roundoff level by construction.

An average throughout the entire magnetosphere volume of the local angle between current and magnetic field^{2}(52)where the sum is performed over each node (i,j).

A last important crosscheck is the consistency of the functions I(Γ) and α(Γ). First, we check that, for each ncomponent, the three functions α_{n}(r,θ) = 4πJ_{n}/cB_{n}(r,θ) are the same. Second, at each radius or angle, the relation (9) has to be satisfied.
Hereafter, we show the electromagnetic quantities in units of B_{0}, r_{ ⋆ }, c, so the magnetic field B scales with B_{0}, the current density J with cB_{0}/r_{ ⋆ }, the enclosed current I with cB_{0}r_{ ⋆ }, and the magnetic flux function Γ with . To test our code and to fix our convergence criteria for the realistic models, we performed a battery of tests. In the first basic test we considered the analytical vacuum dipole with A_{φ} = B_{0}sinθ/2r^{2}, B_{φ} = 0 and checked the ability of the code to maintain this solution. Due to the discretization errors, a little numerical toroidal current appears, and consequently a nonvanishing toroidal fictitious electric field. These errors are small (Δℰ_{b}/ℰ_{b} ≲ 1%), but it is interesting to see how long it takes to dissipate the perturbation to obtain ℰ_{e} = 0 to machine accuracy. In Table 1 we show the results with different resolutions and values of outer radius, always keeping the timestep close to the maximum value allowed by the Courant condition. The expected behavior is obtained, with the constant of proportionality depending on the grid resolution, which affects the strength of the initial numerical current.
Fig. 4 Influence of the outer boundary condition and the location of the external boundary on the evolution of the ratio ℰ_{e}/ℰ_{bφ} (left) and (right). 
Time needed to dissipate the numerical currents of the vacuum dipole.
The second, less trivial test is provided by the analytical models discussed above: the Bessel and selfsimilar solutions (Sects. 2.2, 2.3). Again we began with an initial model consisting of a known solution, and let the system dissipate the currents that come from discretization errors. To obtain the initial models, the selfsimilar solutions require the numerical resolution of the nonlinear ODE (24), while the analytical Bessel solutions are directly implemented. We tried different parameters for the Bessel solutions (varying k) and selfsimilar models (varying the multipole index and the global twist). For every model tested with analytical solutions, we observe a very slight numerical readjustment of the configuration and the code rapidly reaches the relaxed state, with relative changes in ℰ_{b}, ℰ_{bφ}, ℋ less than ~1%. We had to pay special attention to work with sufficient radial resolution in the case of highly twisted Bessel models k ≳ 1, due to their oscillatory radial dependence. For low resolution, the code may find, after a large scale reconfiguration, a completely different solution with smaller k, which is numerically more stable. If the resolution is high enough, all analytical solutions are found to be stable.
We have also studied the evolution of vacuum dipolar solutions with an additional toroidal field for different values of r_{out}. In this case the initial currents are due not only to discretization errors, but also to an inconsistency in the initial model. Moreover, the mean angle defined in Eq. (52) initially has a finite value . In Fig. 4 we show how some convergence monitors evolve as a function of time for three cases with an initial toroidal field of the form B_{φ} = 0.1B_{0}sinθ/r^{3} (), but different external boundary conditions: matching with vacuum at r_{out} = 10 or imposing E_{f} = 0 at r_{out} = 10 or 100. For comparison, we also show results for two analytical solutions: a Bessel solution and a twisted selfsimilar model with the same helicity.
Finally, we tested a vacuum dipole perturbed by a weak toroidal field. This is a case of physical interest for quasi periodic oscillations of magnetars, as discussed in Timokhin et al. (2008) and Gabler et al. (2011). Given a background poloidal field described by Γ, we chose an arbitrary functional form I(Γ) and built the toroidal field according to Eq. (5). As expected, the perturbed configuration is stable and the stationary solution is rapidly reached after a small readjustment. Typically, for max(B_{φ}) = 0.1 B_{0}, we have and changes ΔB/B ~ 1%.
In general, the magnetic energy is not conserved, since the system has to dissipate part of the current to reach a forcefree configuration. This effect is more evident for initial configurations with high helicity. When the outer boundary condition E_{f} = 0 is imposed, the helicity is conserved within a few percent, as expected (see Appendix A for the helicity conservation theorem), and both electric field and mean angle eventually vanish (to machine accuracy). However, when r_{out} is large, or when vacuum boundary conditions are imposed, configurations with high initial helicity take a much longer time to relax (see Fig. 4). In all cases, the relaxation process is faster near the surface, where the configuration of the magnetosphere is more important for our purposes.
On the basis of these results, our convergence criteria for accepting that a configuration has reached a forcefree state are hereafter ℰ_{e}/ℰ_{bφ} < 10^{8} and , plus the requirement that both quantities are monotonically decreasing with time. Some short, initial relaxation phase, in which some largescale reconfiguration occurs is possible. We chose to compare the electric energy to the magnetic energy contribution from the toroidal field, which is much more restrictive than simply the ratio of electric to magnetic energy, especially for low helicity.
Fig. 5 Current density distribution J near the surface for solutions obtained with the same initial data and boundary conditions, E_{f} = 0, but varying r_{out} = 5,10,50,100 (left to right). 
4. Results
With the numerical code described above, we can obtain very general solutions of forcefree, twisted magnetospheres. We discuss separately the influence of the following relevant parameters:

the location of the outer radiusr_{out} in combination with the external boundary condition;

angular and radial dependence of the initial toroidal field;

initial twist and helicity, fixed by the functional form and the strength of the initial toroidal field;

the geometry of the initial poloidal field.
In Fig. 5 we compare the distribution of currents in a solution obtained by imposing E_{f} = 0 at r_{out} = 5,10,50, or 100. In all cases the initial poloidal component is a vacuum dipole solution with A_{φ} = B_{0}sinθ/r^{2} and a toroidal field of the form B_{φ} = 0.1B_{0}sinθ/r^{3}. We observe that, near the axis, the solutions are clearly affected by the location of the external boundary, if it is not far enough from the surface (r_{out} ≲ 10). In such a case, the influence of the external boundary is important, and it introduces artificial features, although the current distribution in the equatorial region is similar in all cases. The final configurations become almost indistinguishable when r_{out} = 50 or 100. Taking r_{out} ≳ 100, we ensure that the numerical noise caused by the interaction with the external boundary is negligible. By neglecting the contribution from the open field lines, for which the twist is illdefined, the global amount of twist is similar in all models (~1.2 rad).
Fig. 6 Enclosed current I(Γ) for different boundary conditions (E_{f} = 0 or coupled with vacuum) at r_{out} = 5,10,50, or 100. 
The function I(Γ) for the same four cases is shown in Fig. 6, together with two more cases corresponding to r_{out} = 5 and r_{out} = 100 but replacing the external boundary condition E_{f} = 0 by a smooth matching with vacuum solutions. The use of a different boundary condition affects the final solution only for low values of r_{out}. Matching with external vacuum implies that no currents can flow through the boundary, so that I(Γ) = α(Γ) = 0 along every field line crossing the outer boundary. As a consequence, a forcefree configuration, coupled with a vacuum, will be characterized by a plateaux I(Γ) = 0 for Γ < Γ_{c}, with Γ_{c} labeling the first closed field line. This means that a bundle of open field lines is potential. The fraction of open field lines (the length of the plateaux for low Γ) is large only for low values of r_{out}. For Γ > 0.2 (equatorial region), all curves coincide. Increasing r_{out} > 100 further has no visible effect on the models with this helicity. For models with higher helicity, the interaction with the boundary becomes more important, and r_{out} needs to be accordingly increased to minimize the boundary effects.
The next step is to explore the influence of the strength and form of the initial toroidal field. Table 2 summarizes the parameters of the initial models employed and the initial mean angle between B and J, Eq. (52). Models AF are obtained starting from a dipolar poloidal component with strength at pole B_{0} and a toroidal component given by the general form (53)The angular part is chosen to be of the form f_{in}(θ) = sin^{d}θ, with d being a positive integer. The radial dependence of models A to E is a power law, g_{in}(r) = r^{−s}, while in model F we use a rapidly decaying function: (54)The initial configuration of model G is asymmetric with respect to the equator: it is a superposition of dipolar and quadrupolar poloidal components plus a toroidal component, as follows: (55)We fix r_{out} = 100 and the external boundary condition to E_{f}(r_{out}) = 0 for all models, except for model C (highest helicity), where r_{out} = 500 to reduce the influence of boundary, as discussed before.
Table 3 shows the features of our final configurations: helicity, maximum line twist, maximum value of current density, and parameters of I(Γ) = I_{0}(Γ/Γ_{0})^{1 + 1/p}. We also include (model S1 and S2) two selfsimilar solutions (Sect. 2.3) of similar helicity, where all components have the same radial dependence B_{i} ~ r^{−(p + 2)}. In this case the current function is analytical, I = I_{0}(Γ/Γ_{0})^{1 + 1/p} with (p,I_{0}) belonging to the family of solutions in Fig. 2. In the other cases I_{0} and p are obtained with fits to the numerical function.
Comparison between two selfsimilar solutions (S1S2) and our numerical solutions (AG).
Fig. 7 Selfsimilar model S2. From left to right: poloidal magnetic field lines (white) and strength of the toroidal component (colored logarithmic scale); current density distribution J in the near region (r ≤ 5) (colored linear scale); angular profiles of the magnetic field components close to the star surface, and angular profiles of the current density. 
The final geometry of the magnetic field and currents for all these models is shown in Figs. 7 to 14. For k_{tor} ≲ 0.1 (models A and B), the initial poloidal field remains almost unaltered, and the behavior of the solutions is nearly linear. Toroidal field strength, helicity, current density J, enclosed current function I(Γ), and global twist scale linearly with k_{tor}. In contrast, for k_{tor} ≳ 0.1, the high initial helicity results in a larger twist angle, up to several radians, which in turn corresponds to a highly deformed poloidal field. A direct comparison of models A, B and C, which differ from each other only in the strength of the initial toroidal component, illustrates this effect: models A and B have the same shape with just a different scale factor, but model C is qualitatively different.
The general features in the low k_{tor} models do not differ much from selfsimilar solutions, because the initial conditions were close to a slightly twisted dipole. In selfsimilar models, two features are the absence of radial currents on the axis and a higher concentration of currents around the equatorial plane. Conversely, in our models, currents are spread more over the angular direction, and we allow for the existence of radial current on the axis. As a consequence, comparing numerical solutions A and B with selfsimilar models with comparable helicity, the former reaches lower maximum values of current density with a higher global twist. We also note that in the most extreme case (model C) the angular dependence of the toroidal field and radial currents becomes steeper^{3}. It is also interesting to compare models C and S2 (Figs. 10 and 7). Both have a similar helicity, but the global twist is higher in model C, while the maximum current density is higher in model S2.
Comparing models B and D, which only differ in the angular dependence of the initial data and the normalization (fixed to obtain the same helicity), we find that the converged solutions are very similar, except near the axis where model D has no radial currents. The effect of varying the initial radial dependence can be estimated by comparing model D to model E. It seems that in this case the final solution keeps memory of the initial model: the converged solution shows a toroidal field that decreases with distance faster in the model E than in model D.
In model E a tiny toroidal field (note the color scale in the figures) appears near the axis. This is likely a numerical artifact that can be partly ascribed to the (narrow) bundle of lines that depart from polar region and interact with the outer boundary. As a matter of fact, these structures are stronger for low values of r_{out}, as already shown in Fig. 5. Moreover, the numerical dissipation of the current is slower near the axis and longer runs are needed to reach more restrictive convergence criteria.
The different radial dependences of the magnetic field components in models B, E and F, together with S1, are shown in Fig. 15. The selfsimilar solution has the same radial dependence r^{−(p + 2)} for the three components, but in the numerical solutions the toroidal field can decrease faster (model E and F) or slower (model B) than the poloidal components. Furthermore, the radial behavior may depend on the magnetic colatitude θ, too. In addition, the radial dependence of the poloidal components is close to the power law r^{3} when the twist is low, but it may significantly deviate from a power law for the more complex models. The departure from selfsimilar solutions is likely to have a visible effect on the observed spectrum, as discussed in more detail in the next section.
A different solution that can only be found numerically is model G, the asymmetric configuration. Most of the current density is concentrated at a high latitude θ_{m}.
An interesting case (not shown) consists of an initial model with toroidal and poloidal fields of opposite parity (e.g., a dipolar poloidal field plus a quadrupolar toroidal field). In this case, the total initial helicity is zero and, since this is a conserved quantity, the current is dissipated and a potential solution is found by the numerical code. This is consistent with the fact that the numerical evolution converges towards the most trivial solution with the same helicity.
Fig. 15 Radial profiles of the magnetic field components for models S1, B (top), E and F (bottom). We show B_{r} on the axis and B_{θ},B_{φ} at the equator for r ∈ [1,30] . 
Finally, to better understand the differences among models, a comparison between the enclosed current function I(Γ) is very helpful, as shown in Fig. 16. In principle, if we know this function or an approximation fits the results of numerical simulations, the pulsar equation (Eq. (7)) can be solved and the magnetospheric structure can be determined. In our models A, B, D, E, G, the current I(Γ) is monotonic and can be well fitted by a single power law I_{0}(Γ/Γ_{0})^{1 + 1/p}, but the values obtained for the parameter p are not consistent with the value of I_{0} describing the selfsimilar family of solutions (Fig. 2), as expected. As a matter of fact, for models A, B, D, the values of p are greater than 1 (the selfsimilar dipolar family is described by p ∈ [0,1] ). In contrast, model E lies quite close to the selfsimilar solution. The enclosed current of models C and E varies more rapidly, describing the concentration of currents in a smaller angular region. For some models, a powerlaw fit to I(Γ) simply does not work. In model G, even if the symmetry with respect to the equatorial plane is broken, the resulting enclosed current I(Γ) is not very different from models A, B, D, and S1 (after rescaling accordingly to the factor k_{tor}). The difference is the mapping between Γ and the surface footprints. In all symmetric models, the maximum magnetic flux (proportional to Γ) is located at the equator, while it corresponds to θ_{m} ≠ π/2 for model G. Depending on the viewing angle, this has a strong imprint on the processed spectra.
Fig. 16 I(Γ) for different models. The curve of model A is magnified by a factor of 10. 
4.1. Resonant optical depth
A prime candidate mechanism for generating the hard tail component observed in magnetar Xray spectra is the resonant Compton upscattering of photons in an external magnetic field. We denote the cyclotron frequency by ω_{B} = ZeB/mc, where Ze and m are the charge and mass of the particle. For electrons, this corresponds to an energy ħω_{B} = 11.6(B/10^{12} G) keV. In the Xray band, the plasma frequency is much lower than the photon frequency, thus taking the polarization vector orthogonal to the propagation direction (semitransverse approximation) works well (see Canuto et al. 1971; and Ventura 1979). The normal modes of propagation are in general elliptically polarized and the degree of ellipticity depends on the ratio ω/ω_{B} and θ_{kB}, defined as the the angle between and B. It is common to introduce the ordinary (O) and extraordinary (X) modes. In the limit of propagation perpendicular to the magnetic field (θ_{kB} = π/2), the modes are linearly polarized, with the polarization vector parallel (Omode) or orthogonal (Xmode) to B. For parallel propagation (θ_{kB} = 0), we recover the circularly polarized modes.
Fig. 17 Estimated resonant optical depth (multiplied by the microphysical factor κ) for radial photons as a function of the polar angle θ. We show results for models S1, A, B, C (first row), D, E, F, G (second row) for different photon energies as indicated in the legend. We set B_{0} = 10^{13} G. 
The cross section of O and E modes depends in a nontrivial way on θ_{kB}, and on the ratio ω/ω_{B}. The cross section of the Omode is close to the Thomson cross section σ_{T} for θ_{kB} near to π/2, it scales with sin^{2}θ_{kB} if ω ≪ ω_{B}, and it does not depend dramatically on the frequency. Conversely, for low photon energies ħω ≪ ħω_{B}, the cross section of the Xmode is strongly suppressed due to the reduced mobility of charged particles across magnetic field lines. Moreover, at photon frequency ω = ω_{B}, the Xmode becomes resonant (Ventura 1979). If we neglect the thermal velocity of the particles (cold plasma approximation) and the natural width of resonance due to the radiationdamping force, the resonant cross section can be approximated by a delta function: (56)where the factor (1 + cos^{2}θ_{kB}) arises from the assumption of unpolarized light. The first attempts to consider this process quantitatively in magnetars are presented by Lyutikov & Gavriil (2006). They study a simplified, semianalytical 1D model by assuming that seed photons are radially emitted from the NS surface with a blackbody spectrum (BB) and the resonant Compton scattering occurs in a thin, plane parallel magnetospheric slab. By assuming a bulk velocity of electrons, but neglecting all effects of their recoil (Thomson limit ħω_{B} ≪ m_{e}c^{2}/γ, where γ is the Lorentz factor of electrons), it can be estimated how a BB spectra is affected. Lyutikov & Gavriil (2006) find an average upscattering for the transmitted radiation (forward scattered photons), while the mean energy of the reflected radiation remains the same. Later, Nobili et al. (2008a,b) predicted spectra obtained via a new 3D Monte Carlo code, accounting for polarization, QED effects, and relativistic cross sections. Although including our magnetospheric models in such sophisticated Monte Carlo simulations is beyond the scope of this paper, we can estimate some of the possible effects by looking at the behavior of relevant quantities.
In particular, we can estimate the resonant optical depth, given by the integration of Eq. (56) along a given line of sight of constant θ (Lyutikov & Gavriil 2006): (57)where n_{Z} is the density of scatterer particles, and all the quantities are evaluated at the resonant radius r_{res}. As the latter depends on the photon energy, the energy dependence of the optical depth is given basically by the local ratio (1 + cos^{2}θ_{kB})n_{Z}/dB/dr (provided that r_{res} lies above the star surface). If we assume charge separation, the particle density is proportional to the current density, J = κcn_{Z}Ze. The dimensionless factor κ depends locally on the plasma composition and on the (bulk and/or thermal) velocity of the particles. Following Thompson et al. (2002) and Nobili et al. (2008a), we can (very roughly) estimate the optical depth along a line of sight (assuming κ constant along it), and considering only radially directed photons: (58)If we consider selfsimilar magnetosphere models under these naive assumptions, the optical depth becomes independent of r_{res}, hence, independent of where the scattering happens. This is because the local ratio (1 + cos^{2}θ_{kB})J/dB/dr is the same for each angle, since for every component i, B_{i} ∝ r^{−p + 2}, and J_{i} ∝ r^{−p + 3}. Furthermore, the optical depth does not depend on the normalization of the magnetic field, B_{0}.
In contrast, the optical depth in our numerical models depends on the photon energy, because the components of J and B cannot be described by the same power law. In Fig. 17 we show the estimated resonant optical depth, Eq. (58), for our models, compared with the selfsimilar model S1. We plot κτ_{res} as a function of the polar angle θ, for different energies of the seed photons, ħω = 0.5,1,2,4,8 keV, taking B_{0} = 10^{13} G. If the magnetic field is predominantly poloidal (low helicity), the optical depth is roughly given by the ratio between toroidal field (which provides poloidal current) and poloidal field at the resonant radius. Increasing the photon energy, the resonant radius will be closer to the surface, modifying the estimate of J and dB/dr which are involved in Eq. (58). Thus, if the toroidal field decays slower than the poloidal component, a higher photon energy implies a lower ratio J/dB/dr at the resonant radius. Looking at the radial profiles (Fig. 15), we can understand why the optical depth increases with the photon energy for models E and F, while for other models the optical depth is greater for soft photons. A more precise prediction of how the spectrum is modified when using one or another magnetosphere model requires more elaborated calculations than those presented here. We also point out that, due to the linear relation between B and ω_{B}, increasing the energy is equivalent to decreasing the normalization B_{0} by the same factor.
5. Conclusions
The study of forcefree magnetosphere models from an analytical or semianalytical point of view have led to remarkable results in a few cases. In the literature, the twisted configurations usually implemented in the modeling of synthetic spectra are often based on the selfsimilar solutions by Thompson et al. (2002), which are restricted to a twisted dipole or to a single higher multipole (Pavan et al. 2009). Since semianalytical studies of a general combination of multipoles are more difficult due to the nonlinear character of the equations, we faced the problem by constructing and thoroughly testing a numerical code that relaxes an arbitrary initial model to a forcefree solution with given boundary conditions. With our numerical simulations we can find general solutions of twisted magnetospheres with complex geometries. Our only input is the radial field at the surface, which can be completely general (i.e., provided by numerical simulations of the internal evolution).
Our main conclusion is that the possible family of relaxed configurations depends on the initial data. As already pointed out by Contopoulos et al. (2011) in a different context, the solution is not unique, given the surface radial magnetic field and the outer boundary conditions. This reflects the freedom we have in choosing the enclosed current function I(Γ), for a given Γ(θ) at the surface. Starting with two different configurations with the same value of volumeintegrated helicity does not necessarily lead to the same final configuration. Thus we were able to find a variety of solutions, which can be qualitatively different from the selfsimilar models.
Some of the new magnetosphere models that we found are characterized by a high degree of nonlinearity, twists up to a few radians, and a nontrivial functional form of the threading current I(Γ). These models establish a new basis for generalizing the study of radiation transfer in neutron star magnetospheres, as already extensively done with the selfsimilar models (Rea et al. 2008; Fernández & Thompson 2007; Zane et al. 2009). In particular, our approach provides a more general framework for studying how the geometry of the magnetic field and currents affects the output Xray spectrum of a magnetar. If resonant Compton scattering is a crucial ingredient in forming the spectra, relaxing the selfsimilarity constraint in the models may have a significant effect, as our preliminary estimates of optical depths have shown. This opens the possibility of implementing alternative models in simulations of radiative transfer. In the future we aim at providing spectra to be compared with observations.
The updated catalog of these sources can be found at http://www.physics.mcgill.ca/~pulsar/magnetar/main.html
Acknowledgments
We are grateful to P. CerdáDurán, M. Gabler, and T. Akgün for interesting discussions and comments. This work was partly supported by CompStar, a Research Networking Program of the European Science Foundation and the Spanish grant AYA 201021097C0302. D. Viganò is supported by a fellowship from the Prometeo program for research groups of excellence of the Generalitat Valenciana.
References
 Beloborodov, A. M., & Thompson, C. 2007, ApJ, 657, 967 [NASA ADS] [CrossRef] [Google Scholar]
 Canuto, V., Lodenquai, J., & Ruderman, M. 1971, Phys. Rev. D, 3, 2303 [NASA ADS] [CrossRef] [Google Scholar]
 Chiu, Y. T., & Hilton, H. H. 1977, ApJ, 212, 873 [NASA ADS] [CrossRef] [Google Scholar]
 Contopoulos, I., Kazanas, D., & Fendt, C. 1999, ApJ, 511, 351 [Google Scholar]
 Contopoulos, I., Kalapotharakos, C., & Georgoulis, M. K. 2011, Sol. Phys., 269, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Duncan, R. C., & Thompson, C. 1992, ApJ, 392, L9 [NASA ADS] [CrossRef] [Google Scholar]
 Fernández, R., & Davis, S. W. 2011, ApJ, 730, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Fernández, R., & Thompson, C. 2007, ApJ, 660, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Gabler, M., Cerdá Durán, P., Font, J. A., Müller, E., & Stergioulas, N. 2011, MNRAS, 410, L37 [NASA ADS] [Google Scholar]
 Gaunt, J. A. 1930, Roy. Soc. London Philosoph. Trans. Ser. A, 229, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Goldreich, P., & Julian, W. H. 1969, ApJ, 157, 869 [NASA ADS] [CrossRef] [Google Scholar]
 Gruzinov, A. 2005, Phys. Rev. Lett., 94, 021101 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Kalapotharakos, C., & Contopoulos, I. 2009, A&A, 496, 495 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Low, B. C., & Lou, Y. Q. 1990, ApJ, 352, 343 [NASA ADS] [CrossRef] [Google Scholar]
 Lyutikov, M., & Gavriil, F. P. 2006, MNRAS, 368, 690 [NASA ADS] [CrossRef] [Google Scholar]
 Mereghetti, S. 2008, A&A Rev., 15, 225 [NASA ADS] [CrossRef] [Google Scholar]
 Michel, F. C. 1973a, ApJ, 180, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Michel, F. C. 1973b, ApJ, 180, L133 [NASA ADS] [CrossRef] [Google Scholar]
 Michel, F. C. 1980, Ap&SS, 72, 175 [NASA ADS] [CrossRef] [Google Scholar]
 Nobili, L., Turolla, R., & Zane, S. 2008a, MNRAS, 386, 1527 [NASA ADS] [CrossRef] [Google Scholar]
 Nobili, L., Turolla, R., & Zane, S. 2008b, MNRAS, 389, 989 [NASA ADS] [CrossRef] [Google Scholar]
 Pavan, L., Turolla, R., Zane, S., & Nobili, L. 2009, MNRAS, 395, 753 [NASA ADS] [CrossRef] [Google Scholar]
 Pons, J. A., & Geppert, U. 2007, A&A, 470, 303 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rea, N., & Esposito, P. 2011, in HighEnergy Emission from Pulsars and their Systems, ed. D. F. Torres, & N. Rea, Science Proc. (Springer), 247 [Google Scholar]
 Rea, N., Zane, S., Turolla, R., Lyutikov, M., & Götz, D. 2008, ApJ, 686, 1245 [NASA ADS] [CrossRef] [Google Scholar]
 Roumeliotis, G., Sturrock, P. A., & Antiochos, S. K. 1994, ApJ, 423, 847 [NASA ADS] [CrossRef] [Google Scholar]
 Scharlemann, E. T., & Wagoner, R. V. 1973, ApJ, 182, 951 [NASA ADS] [CrossRef] [Google Scholar]
 Seehafer, N. 1978, Sol. Phys., 58, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Spitkovsky, A. 2006, ApJ, 648, L51 [NASA ADS] [CrossRef] [Google Scholar]
 Taflove, A., & Brodwin, M. E. 1975, IEEE Trans. Microwave Theory and Techniques, 23, 623 [Google Scholar]
 Thompson, C., & Duncan, R. C. 1996, ApJ, 473, 322 [NASA ADS] [CrossRef] [Google Scholar]
 Thompson, C., Lyutikov, M., & Kulkarni, S. R. 2002, ApJ, 574, 332 [NASA ADS] [CrossRef] [Google Scholar]
 Timokhin, A. N., Eichler, D., & Lyubarsky, Y. 2008, ApJ, 680, 1398 [NASA ADS] [CrossRef] [Google Scholar]
 Ventura, J. 1979, Phys. Rev. D, 19, 1684 [NASA ADS] [CrossRef] [Google Scholar]
 Wolfson, R. 1995, ApJ, 443, 810 [NASA ADS] [CrossRef] [Google Scholar]
 Yang, W. H., Sturrock, P. A., & Antiochos, S. K. 1986, ApJ, 309, 383 [NASA ADS] [CrossRef] [Google Scholar]
 Yee, K. S. 1966, IEEE Trans. Antennas Propagat., 302 [Google Scholar]
 Zane, S., Rea, N., Turolla, R., & Nobili, L. 2009, MNRAS, 398, 1403 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Magnetic helicity
For a given magnetic field configuration, the usual definition of magnetic helicity (A.1)is not unique due to the gauge freedom in the vector potential, A → A + ∇Ψ, where Ψ is a scalar function. However, in axial symmetry (no φdependence), the following definition of helicity (A.2)is gaugeindependent because . We have employed this definition of helicity in numerical code. It can also be shown that this quantity is conserved during the evolution. Taking the time derivative (denoted by ∂_{t}) and using the induction equation, we have (A.3)If the poloidal electric field vanishes at the boundaries, the total helicity is conserved.
With the definition (A.2), the explicit expression of the magnetic helicity for selfsimilar models (Sect. 2.3) is (A.4)and for a dipolar Bessel model (Sect. 2.2) we have (A.5)
Appendix B: A nontrivial semianalytical configuration
Here we show the derivation of Eq. (29) in Sect. 2.4. Substituting in Eq. (15), with Γ given by Eq. (11), leads to (B.1)The integral at the righthand side can be written as where the third term comes from sgn′(x) = 2δ(x). The first and third terms in this equation are zero. Next, we assume by simplicity that Γ ≥ 0 in the whole domain, so that sgn(Γ) ≡ 1, and we define the dimensionless radial functions f_{m}(r) = a_{m}(r)r/r_{ ⋆ }. Hereafter we drop the dependences on μ, r, and the integration limits for conciseness. Using the Legendre equation we obtain from Eq. (11) (B.2)and we can express the integral ℐ in a more compact form (B.3)where we have also used the recurrence relations between Legendre polynomials. Finally, the ODE (15) for each f_{l}(r) is (B.4)where The integral of the product of three Legendre polynomials can be expressed by Wigner 3 − j symbols:
(B.5)Alternatively, we can express the factors above in terms of associated Legendre polynomials: (B.6)These integrals can be evaluated analytically using Gaunt’s formula (Gaunt 1930) (B.7)The integral ℐ_{3} is nonzero if and only if both the following relations are satisfied:

b − c ≤ a ≤ b + c, a triangular condition;

a + b + c is even.
For each fixed l in Eq. (B.6), the conditions above are satisfied for an infinite number of pairs (m,n), with g_{lmn} having numerical values of ~O(1). Thus Eq. (B.4) couples each lmultipole with an infinite number of other multipoles.
It should be noted that this is valid only if Γ ≥ 0 in its whole domain (r,μ) ∈ [1,r_{max}] × [−1,1] , otherwise we cannot use the formulas above (Eqs. (B.5), (B.7)). A non positive Γ would make the calculation of these factors much harder.
All Tables
Comparison between two selfsimilar solutions (S1S2) and our numerical solutions (AG).
All Figures
Fig. 1 Relation between magnetic flux function Γ, enclosed current function I(Γ), and toroidal field. 

In the text 
Fig. 2 Curve I_{0}(p) describing the family of selfsimilar twisted dipoles. 

In the text 
Fig. 3 Location of the variables in the numerical grid. 

In the text 
Fig. 4 Influence of the outer boundary condition and the location of the external boundary on the evolution of the ratio ℰ_{e}/ℰ_{bφ} (left) and (right). 

In the text 
Fig. 5 Current density distribution J near the surface for solutions obtained with the same initial data and boundary conditions, E_{f} = 0, but varying r_{out} = 5,10,50,100 (left to right). 

In the text 
Fig. 6 Enclosed current I(Γ) for different boundary conditions (E_{f} = 0 or coupled with vacuum) at r_{out} = 5,10,50, or 100. 

In the text 
Fig. 7 Selfsimilar model S2. From left to right: poloidal magnetic field lines (white) and strength of the toroidal component (colored logarithmic scale); current density distribution J in the near region (r ≤ 5) (colored linear scale); angular profiles of the magnetic field components close to the star surface, and angular profiles of the current density. 

In the text 
Fig. 8 Same as Fig. 7 for model A. 

In the text 
Fig. 9 Same as Fig. 7 for model B. 

In the text 
Fig. 10 Same as Fig. 7 for model C. 

In the text 
Fig. 11 Same as Fig. 7 for model D. 

In the text 
Fig. 12 Same as Fig. 7 for model E. 

In the text 
Fig. 13 Same as Fig. 7 for model F. 

In the text 
Fig. 14 Same as Fig. 7 for model G. 

In the text 
Fig. 15 Radial profiles of the magnetic field components for models S1, B (top), E and F (bottom). We show B_{r} on the axis and B_{θ},B_{φ} at the equator for r ∈ [1,30] . 

In the text 
Fig. 16 I(Γ) for different models. The curve of model A is magnified by a factor of 10. 

In the text 
Fig. 17 Estimated resonant optical depth (multiplied by the microphysical factor κ) for radial photons as a function of the polar angle θ. We show results for models S1, A, B, C (first row), D, E, F, G (second row) for different photon energies as indicated in the legend. We set B_{0} = 10^{13} G. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.