© ESO, 2013

1. Introduction

The problem of reconstructing the solar coronal magnetic field from photospheric boundary data has been a very active research topic during the past few years. This is mainly due to the arrival of high-resolution and low-noise vector magnetographs either on the ground, such as THEMIS and SOLIS, onboard present solar space missions, such as Hinode and SDO/HMI, and on future ones, such as EST and SOLAR-ORBITER. The various methods proposed so far to reconstruct the coronal field have been reviewed very recently (Aly & Amari 2007; Wiegelmann 2008; Schrijver et al. 2006). We just recall here that the most frequently used ones are (i) the optimization method (Wiegelmann 2004), which uses all photospheric data and tries to minimize a cost function measuring the difference between the computed transverse field and the measured one; (ii) the vertical integration method, which also uses all photospheric data (Song et al. 2006); (iii) the relaxation methods (Valori et al. 2005), in which the coronal field is obtained as the result of an MHD evolution driven by boundary conditions depending on the observed field; and (iv) the Grad-Rubin methods (Sakurai 1981; Amari et al. 1999, 2006; Wheatland 2007; Wheatland & Regnier 2009), which only use a part of the observational data. There is also a hybrid method called OGRM (Amari & Aly 2010) which combines a Grad-Rubin iterative scheme with an optimization scheme by selecting the most appropriate boundary values of the force-free function α at each step.

Up to now, these methods have been essentially applied to the reconstruction of the coronal field of an active region. On such a scale, the curvature of the solar surface does not play an important role, and the region can be represented by a half space. However, the recent observational progress recalled above allow us now to dispose of high-resolution composite vector magnetograms made of several active regions or embedded in a full disk vector magnetogram of synoptic vector magnetograms. Furthermore, vector magnetic fields are also currently measured by different technics on the whole surface of some stars (Donati et al. 2007), although at a lower resolution and with less accuracy. Reconstruction of the global coronal field of the Sun and the other stars thus becomes possible in principle, but this requires implementing the various methods in the spherical geometry context because surface sphericity can obviously no longer be neglected. This has already been done for the optimization method (Wiegelmann 2007; Tadesse et al. 2009) and for the vertical integration method (Song et al. 2007).

In the present paper we present a new numerical code, called XTRAPOLS, which allows reconstruction of the nonlinear force-free field outside of a spherical region. This code is based on an iterative Grad-Rubin scheme implemented in spherical coordinates. We favor this approach because it corresponds to a well-posed boundary value problem. Moreover, it has already proved to be quite efficient in the previously considered Cartesian context, and has allowed us to obtain significant quantitative results when applied to observational data (Bleybel et al. 2002; Régnier et al. 2002, 2005; Régnier & Amari 2004; DeRosa et al. 2009; Canou et al. 2009; Canou & Amari 2010, 2012; Petrie et al. 2011). XTRAPOLS is a massively parallel code, and as such will be able to handle the large amounts of data currently provided by the magnetographs, with fields being measured on a very large scale and with very high resolution.

The paper is organized as follows. In Sect. 2 we describe the basic elements on which XTRAPOLS rests. We formulate the Grad-Rubin iterative scheme in the context of spherical geometry, explain the particular gauge that is used in the vector potential formulation such that ∇·B = 0 (the construction of the gauge is detailed in Appendix A), and explain some particular technical points. Section 3 is devoted to testing our code. We first describe the semi-analytical force-free magnetic field (Low & Lou 1990) that is used as the test case, and thus present the results obtained in the numerical reconstruction of that field with various numerical resolutions. Finally Sect. 4 gathers our concluding remarks together.

2. Method

We present in this section our method for reconstructing the coronal field in the domain D extending outside the sphere of radius r₀. We denote as r the position vector, and use spherical coordinates (r,θ,ϕ) related in the standard way to the Cartesian coordinates (x,y,z).

2.1. The Grad-Rubin iterative scheme

In our calculations the unbounded domain D is replaced by the bounded numerical box Ω_b. The latter is chosen to be the spherical shell Ω_b = { r: 0 < r₀ < r < r₁ } ⊂ R³ and it is thus bounded by the two spherical surfaces S₀ = { r = r₀ }  (with inner unit normal ) and S₁ = { r = r₁ }  (with inner unit normal ). The Grad-Rubin scheme in Ω_b is identical to the one in a Cartesian box previously described in details in Amari et al. (2006). It alternates between solving an hyperbolic problem and an elliptic one.

The hyperbolic problem is defined by and the elliptic one by Here h_n and λ are the given boundary values for and α, respectively, and is the part of ∂Ω_b where h_n > 0.

Equation (1) means that α^(k) keeps a constant value all along a magnetic line of B^(k). To determine α^(k), we have thus to fix its value λ at only one of the two footpoints of the latter. Here we have chosen the footpoint located on the part of ∂Ω_b of positive polarity. But we could have just as well chosen to fix its value at the footpoint located on the part of ∂Ω_b of negative polarity, where h_n < 0.

The iterative process is started by choosing B⁽⁰⁾ to coincide with the unique potential field B_π whose normal component on the boundary is h_n; i.e., B⁽⁰⁾ is the solution of

2.2. Vector potential formulation

To address the ∇·B = 0 constraint, XTRAPOLS introduces a vector potential A^(k) and ensures that the iterated numerical field B^(k) = ∇ × A^(k). To determine A^(k) in a unique way, we impose the gauge conditions where X_t denotes the component of X tangent to the boundary, and A_π is a specific vector potential of the current-free field B_π; i.e., (11)Here, A_π is chosen to be the “restricted DeVore gauge” (RDVG) described in detail in Appendix A. A_π is thus required to satisfy which implies that (14)To compute A_π on the boundary S₀ ∪ S₁:

• We introduce a scalar potential χ_π(r,θ,ϕ) such that (15)where ∇ _⊥  is the orthoradial component of the gradient. This representation is allowed by Eqs. (12), (13).

• We impose χ_π(r₀,θ,ϕ) as the unique solution with zero average value of the equation (16)Solving this equation determines A_π on S₀.

• We introduce the standard representation B_π = ∇ψ_π in Ω_b, and compute the scalar potential ψ_π by solving the boundary value problem (17)

• χ_πis thus given by (18)from which we get A_π on S₁ by Eq. (15).

2.3. Some specific technical points

In our vector potential formulation, α is determined by solving which is just a rewriting of Eqs. (1), (2). The boundary value of α is thus directly imposed at each iterative step, rather than being progressively increased to its nominal value through a second iterative loop as was done in Amari et al. (1999). In other words, we only keep the inner Grad-Rubin iteration loop in the two-level iteration procedure used in that paper.

Equations (19), (20) are solved by the method of characteristics. For the field B^(k), the characteristic passing through the point r is represented here by the function X^(k)(s,r) that solves (Amari et al. 1999, 2006) The parameter s running along is thus taken to vanish at r and to increase in the direction of B^(k), and it takes negative values on that part of comprised between the footpoint located on and r.

For any node r_h ∈ Ω_b at which α^(k) is defined, one has (23)where is the intersection of with . Because λ is only known at nodes that do not generally coincide with , we use an interpolation from the four nearest neighbors. We define the location of α^(k) on the cell vertices. We use either a high-order Adams-Bashford integration scheme with adaptive step size or a fourth-order Runge Kutta scheme, which also allows us to capture the ending point of the characteristics accurately at the limits of the computational box.

Computation of the elliptic part of each iterative step is done on a spherical grid. As for the hyperbolic part, we found that computing the characteristics on the spherical grid led to some awkward behavior near the poles (i.e., for θ near 0 and π). We thus used a somewhat more elaborate method that turned out to be robust, but slower than the previous one. It consists in writing the equation of the characteristics in Cartesian coordinates and computing a characteristic by performing at each point in Ω_b a mapping from spherical to Cartesian coordinates and a mapping back to spherical ones.

The algorithm is parallel in the following sense. The elliptic solver is performed by full domain decomposition with our iterative preconditioned conjugate gradient algorithm to get A and B. The global field as well as the global mesh are therefore known for the hyperbolic part and for each process (each spherical subdomain Ω_i). But each process can consider the calculation of alpha independently for each set of starting point in Ω_i. Although it is not full parallelism in the sense that the magnetic field B^(k) is global, it clearly speeds up the calculation and does not overload the memory too much since only one global field is needed by each process.

Finally we point out that our convergence criterion for the sequence of fields B^(k) is chosen to be (24)where ϵ = 10^-6 and we have introduced the L²(Ω_b)-norm of a vector B, defined by (25)As shown in Sect. 3.2, a total number of iterations of about 30–50 is sufficient to achieve convergence in most cases, and the stopping criteria is robust with respect to changes in the data put as boundary conditions.

3. The test case

3.1. The Low and Lou solutions

For testing the quality of our reconstructing method, we used the semi-analytic Low & Lou (1990) force-free magnetic field, B_LL, which has now become a standard benchmark. To describe that particular field, we first set up a new Cartesian coordinate system (x′,y′,z′) centered on some point O′ and define associated spherical coordinates (r′,θ′,ϕ′) in the standard way. Expressed in these coordinates, B_LL is axisymmetric with respect to the axis O′z′ and mirror-symmetric with respect to the equatorial plane  { z′ = 0 } , so it admits the usual representation (26)in terms of the poloidal flux function A(r′,θ′) and the current function F(r′,θ′). Both A and F are taken to have the form respectively, where a and ν are real constants, μ′ = cosθ′, and P is a solution of a nonlinear ordinary differential equation satisfying P(± 1) = 0. We take here a solution that generates a quadrupolar field and has the “stretching parameter” ν = 1.

Next we fix the position of the auxiliary coordinate system (O′;x′,y′,z′) with respect to our original coordinate system (O;x,y,z). This is done by imposing the simple relations x′ = x, y′ = y − d, and z′ = z, with d > 0 a constant, from which we get r′ = (r² + d² − 2drsinθsinϕ)^1/2, cosθ′ = (r/r′)cosθ, and tanϕ′ = (rsinθsinϕ − d)/rsinθcosϕ. Straightforward but somewhat heavy algebra thus allows computing the components of B_LL in both the Cartesian coordinates (x,y,z) and the spherical ones (r,θ,ϕ). Clearly, the latter depend on ϕ.

A similar shifting procedure applied to the Low-Lou field has already been used in Wiegelmann (2007), and we choose for d the value selected in that paper. It will thus be possible to compare our results with those reported in the latter.

3.2. Results

Table 1 shows some of the diagnostics, often called figures of merit, already used for Cartesian coordinates calculations in Amari et al. (1999, 2006) and in Schrijver et al. (2006). If u and v are two vectors, and N is the number of computational nodes, we introduce

• the vector correlation (29)

• the Cauchy-Schwartz correlation (30)

• the normalized vector error (31)

• the mean vector error (32)

• the normalized magnetic energy (with respect to the energy of the exact solution) (33)

The various quantities above are evaluated for u equal to the exact solution B_LL and v equal to a numerical solution obtained with XTRAPOLS. They have to be compared to the reference values We also found important to report the values of the L^∞(Ω_b)-norm of ∇·B, (36)as well as the CPU time, both figures providing additional measures of the quality of the solution.

To compare with the Cartesian case, we show only the cases – corresponding to “Case 3” in Wiegelmann (2007) – where the surfaces S₀ and S₁, on which we impose boundary conditions corresponding to the Low-Lou solution, are fixed by setting r₀ = 1 and r₁ = 2.57, respectively, where S₁ may be considered as some kind of source-surface, although we do not require the field to be radial on it.

Figure 1 shows the convergence properties of XTRAPOLS. We have drawn the rate of convergence on a logarithmic scale as functions of the number of iterations, as well as the L²(Ω_b)-norm of the iterates expressed in nondimensional units. It should be noted that

• The convergence rate appears to exhibit a regular decrease fromthe very beginning of the computation.

• The number of iterations necessary to reach a relative error of 10^-5, say, decreases when the resolution increases. We needed 40 iterations at the lowest resolution, 32 at the intermediate one, and only 27 at the highest one.

• There is a strong variation in the L² norm of the solution (therefore of the total magnetic energy) at the first iteration, at which point the potential state is left over. This might be related to the fact that no progressive injection of α at the boundary (no two-level iteration loop) has been achieved. After a very short “adaptation” step, this wide behavior is damped, and ||B^(k)||_L²(Ωb) exhibits a smooth increase similar to the one obtained in Amari et al. (2006).

• Our convergence criterion ϵ = 10^-6 appears to be relatively strong since ||B^(k)||_L²(Ωb) has already reached its asymptotic value, on the order of 0.01, at about iteration 20.

• These results, which are complementary to the one above, show the benefit of increasing resolution on the value of the L²-norm of the solution.

As shown in Fig. 2, the field lines of the numerically computed nonlinear force-free configuration are very similar to those of the exact solution. They are quite far from those of the potential (current-free) solution.

4. Discussion

In this paper, we have presented a method for reconstructing the solar coronal magnetic field on the full Sun scale. This method is an adaptation to spherical coordinates of the Grad-Rubin algorithm, based on a well-posed formulation, which we have developed over the years in Cartesian coordinates. We implemented this algorithm in the massively parallel code XTRAPOLS, and tested its efficiency by applying it to reconstructing of a nonaxisymmetric version of the Low-Lou force-free field in spherical coordinates. This reconstruction has been shown to be quite successful. In particular, our results are comparable to those obtained by Wiegelmann (2007) for his “type 3” boundary conditions. We did not consider here the other types of boundary conditions studied in that paper.

Our results show that the code has good convergence properties, as confirmed by the various diagnostics and by the fact that the L²-norm of the solutions, and then the energy, approaches an asymptotic value quite rapidly. But there is not much improvement in some of the diagnostics when the numerical resolution becomes higher. This does not imply, however, that there is no interest in increasing the resolution. For instance, one may want the reconstruction to catch not only the gross structure of the magnetic field, but also the small-scale features that are associated with the active regions and are expected to have a negligible influence on the large-scale field. Unfortunately, the Low-Lou solution, which currently appears as the only analytical solution useful for testing a code, does not exhibit such multi-scale magnetic field and electric current density. The efficiency of the code in such a situation will then be tested by applying it to actual data. This is the next step in our program, and the results will be presented in forthcoming papers.

Our method allows us to deal with the very large amount of data provided nowadays by various vector magnetographs on the ground or onboard satellites, such as SOLIS and HMI. These instruments are currently used for studying transequatorial loops connecting active regions and for building up synoptic maps and composite maps made of full disk maps embedded in synoptics ones. It is important to note that the reliability of a reconstruction depends as much on the quality of these magnetograms as on the accuracy of the numerical algorithm. This quality depends in particular on the method that is used to invert Stokes parameters and then on the assumption made on the value of the filling factor that has to be introduced to deal with both strong (active) and weak (quiet) field regions.

Acknowledgments

We acknowledge support from NASA’s Heliophysics Theory Program, the STEREO/SECCHI Consortium, and the Centre National d’Études Spatiales. The numerical simulations performed in this paper have been done on the IBM Power 6 supercomputer of the Institute I.D.R.I.S of the Centre National de la Recherche Scientifique. We also thank the anonymous referee for his/her helpful remarks that helped improve the clarity of the paper.

References

Appendix A: DeVore and restricted DeVore gauges in spherical geometry

The gauge used in the main text for representing the potential field is the spherical analog of the gauge introduced by DeVore (2000) in the Cartesian context. We explain in this appendix how to construct this gauge for an arbitrary magnetic field B occupying the domain Ω_b, with the value r₁ = ∞ being allowed (in which case Ω_b = D). We also give formulae valid when B is axisymmetric. Although they are not used explicitly in the present paper, they have proved to be useful in some auxiliary computations we did on the Low & Lou (1990) force-free magnetic field, and we think that they may also be useful to some readers.

A.1. DeVore gauge

As is well known, the conditions satisfied by the magnetic field B, allow the introduction of a vector potential A₀ such that (A.3)and all the other gauges have the form A = A₀ + ∇f. We first define the “DeVore gauge” (DVG) A₁ by the condition (A.4)A₁can thus be obtained from A₀ by requiring the arbitrary function f to satisfy (A.5)Note that f is determined by the latter equation only up to an additive function of (θ,ϕ).

In DVG we have (A.6)where the indice  ⊥  indicates an orthoradial component. As ∇ _⊥  × A₁ is collinear to , we can write (A.7)and we get by a mere integration (A.8)This expression leaves the surface function A₁(r₀,θ,ϕ) undetermined.

A.2. Restricted DVG

Taking the divergence of both members of Eq. (A.8) leads to (A.9)Using the gauge arbitrariness still in A₁, we set A = A₁ + ∇ _⊥ f(θ,ϕ), and choose f such that (A.10)We thus get a particular form A of the DVG for which (A.11)and (A.12)We refer to A as the restricted DeVore gauge (RDVG).

Note that (A.13)in the case where ∇ _⊥  × B _⊥  = 0 in Ω_b (no radial current).

Condition (A.11) allows a new potential χ(θ,ϕ) to be introduced such that (A.14)hence (A.15)This equation determines χ(θ,ϕ) from B_r(r₀,θ,ϕ) up to an additive constant whose value may be fixed, e.g., by requiring that (A.16)

A.3. The case of a potential field

Consider the case of a potential magnetic field B_π in Ω_b (∇ × B_π = 0). Then we can introduce the alternative representation (A.17)in terms of the scalar potential ψ_π. The latter satisfies Laplace equation (A.18)and can then be determined (up to a constant) by solving this equation along with the boundary condition (A.19)On the other hand, the RDVG A_π for B_π does satisfy ∇·A_π = 0 (see Eq. (A.13)), and we can introduce a scalar potential χ_π(r,θ,ϕ), defined up to an additive function f(r), such that (A.20)while χ_π(r₀,θ,ϕ) is determined by Eqs. (A.15), (A.16). Using Eqs. (A.20) and (A.17) in Eq. (A.8) for A_π, (A.21)we obtain (A.22)whence after integration (A.23)No constant of integration needs to be introduced here as is immediately seen by setting r = r₀.

A.4. The case of an axisymmetric magnetic field

Let us consider an arbitrary magnetic field B that is axisymmetric in Ω_b with respect to the z-axis – that is to say B is left invariant by the rotations about that axis. Then B has the standard representation (A.24)in terms of the two scalar functions A(r,θ) – the poloidal flux function – and F(r,θ) – the poloidal current function. A(r,θ) is defined up to an additive constant, which we may fix by imposing (A.25)where the first equality is a mere consequence of the fact that B has zero flux through any sphere of radius r, r₀ ≤ r ≤ r₁. Moreover, the regularity of B at the two poles of the sphere requires that To construct a vector potential for B satisfying the DVG condition, we first consider the poloidal part B_p. Noticing that , we can write (A.28)with (A.29)Owing to Eq. (A.26), this quantity is well defined even when θ = 0,π, with (A.30)Clearly we have in the whole D, and then A_P even satisfies the RDVG conditions.

Next we consider the toroidal part of the field, for which we have to compute a vector A_T such that Let us try to find a solution to these equations of the form (A.35)Injecting this form into Eq. (A.33), we get (A.36)hence (A.37)Taking for instance the particular solution (A.38)we obtain (A.39)which is a well defined quantity for θ = 0,π owing to Eq. (A.27), with (A.40)Reciprocally, the A_T we have constructed gives back the toroidal field B_ϕ. It is then a solution to our problem, which actually turns out to also satisfy the RDVG conditions. We indeed have (A.41)

whence (A.42)The field B thus admits the vector potential (A.43)such that and is a RDVG vector potential.

As an application, consider the axisymmetric force-free field B_LL (Low & Lou 1990) centered on O (i.e., without the shift in y introduced in Sect. 3.1). For this field, where μ = cosθ, and P(± 1) = 0, and from the formulae above we immediately get (A.50)

All Tables

All Figures