© ESO, 2011

1. Introduction

Anomalous X-ray pulsars (AXPs) and soft gamma-ray repeaters (SGRs) are a class of neutron stars (NSs) characterized by high X-ray quiescent luminosities, short X-ray bursts, and giant flares (for SGRs). They are believed to be magnetars (Duncan & Thompson 1992; Thompson & Duncan 1996): young (P/2Ṗ ~ 10²−10⁵ yr) NSs with very strong magnetic fields  ~10¹⁴−10¹⁵ 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¹, 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 X-ray persistent luminosities L_X ~ 2 × 10³³−2 × 10³⁵ 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 force-free magnetospheres, which leads to the so-called 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, ad-hoc 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 force-free 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 time-dependent, misaligned case by means of 3D simulations by Spitkovsky (2006) and Kalapotharakos & Contopoulos (2009).

Nowadays, X-ray 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 semi-analytical 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 self-similar twisted dipole proposed by Thompson et al. (2002), extended to higher multipoles by Pavan et al. (2009), which represents a semi-analytical, easy-to-implement 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 force-free 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 force-free 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 large-scale 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 semi-analytical generalizations of previous works. In Sect. 3 we describe the numerical code used to build force-free 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 force-free 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 dipolar-like 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 so-called pulsar equation (Scharlemann & Wagoner 1973; Michel 1973b): (7)where we indicate the partial derivatives of Γ with a subscript. The terms (ρ/r_l)² arise from the Coulomb force ρ_eE. The pulsar equation is actually the Grad-Shafranov equation for force-free 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 divergence-free condition. This is a wind-like 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 force-free 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 Γ₀ the magnetic flux normalization. Denoting the strength at pole by B₀, 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 l-multipole.

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. Self-similar models

Low & Lou (1990), followed by Wolfson (1995) and other authors, studied a particular class of self-similar 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 power-law 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 second-order 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 self-similar 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 l-pole family of self-similar solutions. Due to self-similarity, 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, no current) implies an integer value of p, so that we recover the vacuum multipolar solution with F(μ) = (1 − μ²)dP_p(μ)/dμ where P_p is the pth Legendre polynomial. For each multipole, a unique relation k_ss(p) (or I₀(p) once the values of B₀ 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 I0(p) describing the family of self-similar 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 Γ₀ 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 semi-analytical 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 semi-analytical 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, force-free 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 semi-analytical 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 rotationally-induced 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 magneto-frictional method

For simplicity, we begin with the two-dimensional nonrotating case, assuming that the force-free approximation is valid up to an outer radius r_out, which reduces the problem to finding solutions of Eq. (10). In the so-called magneto-frictional method (Yang et al. 1986; Roumeliotis et al. 1994), one begins with an initially non-force-free configuration and defines a fictitious velocity field proportional to the Lorentz force, v_f = νJ × B/B², 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 force-free 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 force-free 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, rotationally-induced 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 force-free.

3.2. Linear analysis of the magneto-frictional method

We now consider a background, uniform magnetic field B₀ 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₀ or transverse perturbations with k ∥ B₀), 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₀, 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 magneto-frictional 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 force-free 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 divergence-less field, it is guaranteed that the numerical divergence defined by Eq. (49) remains zero to round-off level at each time step.

We must mention that we also tried a method with a staggered grid (Yee 1966), in which each n-component 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 well-suited to solving Maxwell’s equations, because it provides a natural way to time-advance 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 force-free 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 round-off 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 magneto-frictional method is based on introducing an artificial, viscous electric field that drives an arbitrary initial configuration into a force-free 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:

• Volume-integrated magnetic energy (total and contributionfrom the toroidal field)

• Volume-integrated energy stored in the fictitious electric field

• Total volume-integrated helicity (see Appendix A for a discussion about this definition).

• Volume-integrated absolute value of the divergence of both, B and J. These two quantities are expected to vanish at round-off level by construction.

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

• A last important cross-check is the consistency of the functions I(Γ) and α(Γ). First, we check that, for each n-component, 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₀, r _⋆ , c, so the magnetic field B scales with B₀, the current density J with cB₀/r _⋆ , the enclosed current I with cB₀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₀sinθ/2r², 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).

The second, less trivial test is provided by the analytical models discussed above: the Bessel and self-similar 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 self-similar 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 self-similar 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₀sinθ/r³ (), 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 self-similar 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₀, 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 force-free 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 force-free 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 large-scale 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, Ef = 0, but varying rout = 5,10,50,100 (left to right).

4. Results

With the numerical code described above, we can obtain very general solutions of force-free, 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₀sinθ/r² and a toroidal field of the form B_φ = 0.1B₀sinθ/r³. 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 ill-defined, the global amount of twist is similar in all models (~1.2 rad).

Fig. 6 Enclosed current I(Γ) for different boundary conditions (Ef = 0 or coupled with vacuum) at rout = 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 force-free 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 A-F are obtained starting from a dipolar poloidal component with strength at pole B₀ 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₀(Γ/Γ₀)^1 + 1/p. We also include (model S1 and S2) two self-similar 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₀(Γ/Γ₀)^1 + 1/p with (p,I₀) belonging to the family of solutions in Fig. 2. In the other cases I₀ and p are obtained with fits to the numerical function.

Fig. 7 Self-similar 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.

Fig. 8 Same as Fig. 7 for model A.

Fig. 9 Same as Fig. 7 for model B.

Fig. 10 Same as Fig. 7 for model C.

Fig. 11 Same as Fig. 7 for model D.

Fig. 12 Same as Fig. 7 for model E.

Fig. 13 Same as Fig. 7 for model F.

Fig. 14 Same as Fig. 7 for model G.

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 self-similar solutions, because the initial conditions were close to a slightly twisted dipole. In self-similar 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 self-similar 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³. 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 self-similar 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 self-similar 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 Br 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₀(Γ/Γ₀)^1 + 1/p, but the values obtained for the parameter p are not consistent with the value of I₀ describing the self-similar 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 self-similar dipolar family is described by p ∈  [0,1] ). In contrast, model E lies quite close to the self-similar 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 power-law 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 X-ray 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¹²   G) keV. In the X-ray 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 (O-mode) or orthogonal (X-mode) 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 B0 = 1013 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 O-mode is close to the Thomson cross section σ_T for θ_kB near to π/2, it scales with sin²θ_kB if ω ≪ ω_B, and it does not depend dramatically on the frequency. Conversely, for low photon energies ħω ≪ ħω_B, the cross section of the X-mode is strongly suppressed due to the reduced mobility of charged particles across magnetic field lines. Moreover, at photon frequency ω = ω_B, the X-mode 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 radiation-damping force, the resonant cross section can be approximated by a delta function: (56)where the factor (1 + cos²θ_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, semi-analytical 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_ec²/γ, 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²θ_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_ZZe. 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 self-similar 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²θ_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₀.

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 self-similar 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₀ = 10¹³ 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₀ by the same factor.

5. Conclusions

The study of force-free magnetosphere models from an analytical or semi-analytical 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 self-similar solutions by Thompson et al. (2002), which are restricted to a twisted dipole or to a single higher multipole (Pavan et al. 2009). Since semi-analytical 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 force-free 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 volume-integrated 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 self-similar 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 self-similar 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 X-ray spectrum of a magnetar. If resonant Compton scattering is a crucial ingredient in forming the spectra, relaxing the self-similarity 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.

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 2010-21097-C03-02. D. Viganò is supported by a fellowship from the Prometeo program for research groups of excellence of the Generalitat Valenciana.

References

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 gauge-independent 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 self-similar models (Sect. 2.3) is (A.4)and for a dipolar Bessel model (Sect. 2.2) we have (A.5)

Appendix B: A nontrivial semi-analytical 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 ℐ₃ is non-zero 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 l-multipole 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

All Figures

	Fig. 1 Relation between magnetic flux function Γ, enclosed current function I(Γ), and toroidal field.
In the text

	Fig. 2 Curve I0(p) describing the family of self-similar 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, Ef = 0, but varying rout = 5,10,50,100 (left to right).
In the text

	Fig. 6 Enclosed current I(Γ) for different boundary conditions (Ef = 0 or coupled with vacuum) at rout = 5,10,50, or 100.
In the text

	Fig. 7 Self-similar 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 Br 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 B0 = 1013 G.
In the text