© ESO, 2012

1. Introduction

Magnetic reconnection in line-tied low-plasma-β environment such as the solar corona is linked to the formation of intense field-aligned current sheets. These currents develop preferentially in specific location of the magnetic field, in particular regions where the connectivity of the magnetic field is discontinuous. Magnetic topology studies aims at determining such regions (see review of Démoulin 2005). Null points, separatrices and separators are classical region of magnetic field discontinuity (see review by Longcope 2005). Through the identification of such regions, magnetic topology studies have been able to explain the location and shape of flares ribbons (e.g. Gorbachev & Somov 1989; Mandrini et al. 1995; Démoulin et al. 1997; Gaizauskas et al. 1998; Bagalá et al. 2000; Mandrini et al. 2006).

However all flares cannot be directly associated with such magnetic field discontinuity. A generalized topological structure, the quasi-separatrix layers (QSLs) have been introduced by Démoulin et al. (1996b, see also review of Démoulin 2006. QSLs have been defined as regions where the mapping of the field lines has strong gradient. It includes possible discontinuous mapping, hence the presence of separatrices, as a particular case. They were localized by the estimation of the norm, N of the Jacobian matrix of the mapping of the field line connectivity (see Sect. 2). Titov et al. (2002) improved the concept by defining the squashing degree, Q, which, unlike N, is invariant along a field line. QSLs are thus 3D magnetic volumes of high Q, in which the magnetic connectivity varies strongly.

Similarly to separatrices, flare ribbons could be associated with QSLs in a large number of events (e.g. Mandrini et al. 1997, 2006; Masson et al. 2009; Chandra et al. 2011). Indeed, QSLs are also preferential sites for electric current build-up (Milano et al. 1999; Galsgaard et al. 2003; Aulanier et al. 2005; Buechner 2006; Pariat et al. 2006; Masson et al. 2009; Wilmot-Smith et al. 2009, 2010). Magnetic reconnection develops within QSLs, with the particularity that field line continuously reconnect with their neighboring field lines, leading to an apparent slipping of the field lines (Aulanier et al. 2006, 2007; Török et al. 2009; Masson et al. 2009, 2011), while classical reconnection at a separatrice is realized in one step.

QSLs have thus received an increased attention from the community. As the spatial extension of a QSL is typically a very thin 3D volume with a complex shape, the representation of QSLs is typically not straightforward. The 3D representation by an iso-surface of Q can be used (Titov et al. 2002, 2003; Titov 2007). A complementary way to represent complex QSLs can be done by studying the distribution of Q on 2D sections of the studied 3D domain. As the computation of Q is relative to the choice of boundary for the field line mapping, QSLs are usually represented at these boundaries (e.g. Démoulin et al. 1997; Buechner 2006; Masson et al. 2011; Titov et al. 2011). However, it may be necessary to represent QSLs in particular sections of the 3D domain (as in Aulanier et al. 2006; Antiochos et al. 2011; Savcheva et al. 2012a). The aim of this paper is to analyze several possible methods to estimate Q inside the 3D domain and to provide a formulation for the most suitable method. We focus on computing Q on 2D cuts because it has important applications and allows us to focus on the core of the methods.

For illustration, each method will be applied to the example of the magnetic configuration of Titov & Démoulin (1999,TD hereafter). It is formed by a twisted flux rope (part of a torus) in equilibrium in an active-region like magnetic field. The normalized parameters are chosen such as: R = 2, a = 0.45, d = 1, |L| = 1, and |q| = 1, I₀ = 2, where R and a are the large and small radius of the torus, d and q are the depth below the photospheric plane (z = 0) and the magnitude of the magnetic sources, L is half the horizontal distance between the sources, and I₀ is the current intensity of the line current, respectively (see Fig. 2 and Sect. 2.1 in Titov & Démoulin 1999). With these parameters the flux rope has a twist of N_t = 2 above z = 0. The photospheric distributions of B_z and Log₁₀   Q are given in Fig. 1. The TD model has been chosen because its topology has been thoroughly studied (Titov & Démoulin 1999; Titov 2007) and because it defines a structure commonly used as an initial configuration in numerical simulation of eruptions (e.g. Roussev et al. 2003; Török et al. 2004, 2011a,b; Valori et al. 2010; Cohen et al. 2011; Lugaz et al. 2011).

We also apply the methods to a more complex magnetic field derived from observations (Savcheva et al. 2012a). This magnetic field is based on a non-linear force-free field constrained by both the photospheric magnetogram and the observed coronal loops. This magnetic field is associated to a long-standing coronal sigmoid (observed during 6 days) within an active region. The aim of studying such a field is to test the performance of the methods on a more structured magnetic field as typically observed on the Sun.

The paper is organized as follow: after recalling the method to compute Q at the boundary of the domain in Sect. 2, we will present in Sect. 3 several methods to derive Q maps within the domain and investigate how precise they are. Each method is illustrated by its application to the TD model. Then, we compare the best two methods using an observed magnetic configuration (Sect. 4). Finally, we conclude in Sect. 5.

Fig. 1 Left panel: photospheric distribution (grayscale) of the vertical component of the magnetic field, Bz, for the TD model. The continuous (resp. dotted) isolevels of Bz corresponds to intensities of  ±  [1,2,3,4] . Right panel: photospheric distribution (color scale) of Log10   Q. A color version is available online.

2. Squashing degree estimation at the planar boundary of a 3D domain

The practical definition of QSLs in the solar context is inherently relative to the mapping of magnetic field lines from one footpoint to the other. The concept is justified by the fact that field lines are typically line tied in the dense photosphere, while a more general definition of QSLs is possible by analyzing the divergence of neighbor field lines (see Sect. 2.3 in Démoulin et al. 1996b). The computation of QSLs is therefore relative to the choice of the boundaries of the studied domain V. The norm N of the mapping and the squashing degree, Q are defined for given boundaries of the studied domain (Démoulin et al. 1996b; Titov et al. 2002; Titov 2007). By construction the squashing degree is constant for a given field line (Titov et al. 2002), hence, once the boundaries of the domain of study is fixed, Q can be given for the whole domain, each field line having a single value of Q.

The quantities N and Q, are usually derived from the Jacobian matrix of the field-line mapping from the boundary of one footpoint to the other. In the following we will consider that the boundaries associated with each footpoint are planar and thus can be represented by two independent cartesian referentials. For generality we will consider that these planes are arbitrarily oriented regarding to each others. More commonly these plans could be two sides of a cubic 3D domain of a numerical simulation. In the particular case of closed magnetic fields such as in an active region field (and as in the TD model used here) these two planes are the same one: one compute the mapping between the positive footpoint to the negative one of a photospheric-like boundary.

Let us consider a field lines which links the footpoint, r₁, of coordinate (x₁,y₁) of the plane P₁ to the footpoint r₂, of coordinate (X₂,Y₂) of the plane P₂ (see Fig. 2). The plane P₁ will for example correspond to footpoint with a positive magnetic flux while P₂ will correspond to a negative flux. Two mapping exist which associate a footpoint on one plane to the other: the mapping Π₁₂ from P₁ to P₂: r₁(x₁,y₁) → r₂(X₂,Y₂); and the inverse mapping Π₂₁ from P₂ to P₁: r₂(X₂,Y₂) → r₁(x₁,y₁).

The Jacobian matrix, D₁₂ & D₂₁ associated to each mapping are: As Π₁₂ and Π₂₁ are inverse functions (Π₁₂°Π₂₁ = Π₂₁°Π₁₂ = Id), D₂₁ is related to D₁₂ by: (3)with Δ₁₂ the determinant of D₁₂.

Evaluating D₂₁ is analytically equivalent using Eq. (2) or (3), however numerically it is generally not the case. In order to compute the Jacobian matrix, one needs to compute in two orthogonal directions a set of field lines sufficiently close to resolve the gradient of connectivities. This is achieved by computing progressively field lines closer to the selected point (cf. Sect. 3 of Aulanier et al. 2005). It is better numerically to realize such convergence where the value of Q is needed both because of finite difference precision (selection of the finite difference scale) and of large connectivity dependence within QSLs. Such fine treatment is needed to compute Q because of the very severe distortion of the mapping as illustrated in Fig. 2 of Longcope & Strauss (1994), even for broad QSLs.

In presence of QSLs, the Jacobian matrix components are typically very large or very small. As they are evaluated by finite difference of the position of the footpoints of the field lines, it is numerically more accurate to have a fixed denominator for all components.

Fig. 2 Cartoon illustrating the first possible method, but incorrect, to compute Q at rc.

For example, when evaluating the squashing degree Q at r₁, it is better to integrate field lines from fixed points distant of δ_x1 = δ_y1 = δ and centered on r₁. The field lines, separated on P₁ along the x₁ direction by δ_x1( = δ), have a footpoint separation on P₂ which component along the X₂ and Y₂ directions are respectively equal to d_X2x1 and d_Y2x1. In a QSL, d_X2x1 and d_Y2x1 will typically be very large or very small relatively to δ. Because of numerical errors, the quantities such as ∂X₂/∂x₁ ≃ d_X2x1/δ will be evaluated more consistently than the quantities such as ∂x₁/∂X₂. The same is valid for field lines separated by δ_y1 on P₁. We conclude that, when numerically evaluating Q and the mapping properties at r₁, Eq. (1) should be used to compute D₁₂, and D₂₁ should be derived from Eq. (3) instead of Eq. (2).

Titov et al. (2002) demonstrated that the determinant of the Jacobian matrix, Δ₁₂ & Δ₂₁ can be related to the ratio of B_z,1(x₁,y₁) and B_Z,2(X₂,Y₂), which are respectively the normal component to P₁ and P₂ of the magnetic field evaluated in r₁ and r₂: (4)let us note that this expression should be used numerically because it is much more precise than the computation of the determinant of the Jacobian matrix which involves partial derivatives.

The norms N₁₂ & N₂₁ of the two Jacobian matrix are: the Squashing degree, Q, for this field line is given by: while in theory Q₁₂ is equal to Q₂₁ it may not be the case numerically. In addition of the errors related to the computation of the Jacobian matrix at each footpoint explained previously, field line integration errors will also introduce variations. Indeed for a magnetic field represented in a discrete mesh, integration errors of the field lines would lead to finite errors on the location of the footpoints of the field lines. Since QSLs are both associated to large Q values and large Q gradients, except just in the thin central part (as in Fig. 1, right panel), a small spatial shift in the integration of a field line can shift significantly the value of Q. Therefore, it is more precise to numerically evaluate Q at the location where it is needed rather than transferring Q along field lines.

3. Squashing degree estimation within a 3D domain

Let us now focus on the way to compute Q in the volume of the studied domain V. The boundaries P₁ and P₂ stays fixed and the aim is to determine the value of Q associated to an arbitrary point r_c of V. We define a plane P_c represented by a cartesian referential such as r_c has the coordinates (x_c,y_c). The plane P_c can be defined in the neighborhood of r_c, P_c(r_c), implying than its orientation could depend on r_c (for example by taking perpendicular to the local magnetic field, see the end of Sect. 3.4). In this case, the local plane P_c(r_c) can be different for different r_c positions. However, for many applications, it may be more suited to consider an unique plane P_c, e.g. the large scale 2D cut on which Q is to be represented. We follow this latter assumption in the present paper as it simplifies the description and permits to focus more on the key-points of Q computations.

In order to determine the value of Q at r_c, one first needs to determine the location of the two footpoints of the field line passing by r_c. The field integration is assumed to retain the direction of the magnetic field and thus the footpoints are supposed to be properly computed (to the numerical precision) to their respective plane, e.g. positive on P₁ and negative on P₂. As Q is constant for all the points of a given field line, it should be possible to compute Q at any point of the field line. If Q is know analytically at a given point and the connectivity can be computed, it is easy to derive Q at r_c. This is what is used in Wilmot-Smith et al. (2009). However for non analytical fields the computation is not as straightforward and several methods can be possibly used to derive Q(r_c).

3.1. Pseudo-method 1

For didactic purpose, we first analyze an example of a method, which could appear intuitive, but which is incorrectly computing the partial derivatives of the connectivity matrices. This illustrates how tricky can be the computation of Q, and the cautions that must be taken.

When one compute a Q map in a 2D cut of V, as field lines are computed at every mesh point of this map, it is tempting to directly use the position of the computed field lines on the already defined grid and directly use Eqs. (1) and (2) to derive D₁₂ or D₂₁ from the position of the footpoints on P₁ and P₂ only. This intuitive methodology implies to first compute field lines from two directions on P_c around the point of interest r_c as illustrated in Fig. 2. Four field lines are plotted from P_c distant from a fixed position δ centered on r_c (δ being the distance between field lines position on each side of r_c). One obtain the difference of position of the footpoints of these field lines, d_x1, d_y1, d_X2 and d_Y2 on P₁ and P₂. These distances provides and information on how the flux tube is squashed. If one focus on the footpoint r₁, Q₁₂ can eventually be approximated with: (10)at first sight, this estimation can look coherent since it converges to the typical formula to compute Q when P_c converges to P₁ (d_x1 and d_y1 being computed in orthogonal directions). An equivalent relation can be obtained from the footpoint r₂. The results are illustrated in the upper left panel of Fig. 5. One observes values of Q which are smaller than 2 while Q is theoretically greater than 2 (Titov et al. 2002).

However, since one does not control the direction and distance of the location on P₁ (field lines are computed from P_c), the Jacobian matrix D₁₂ is not properly estimated. Starting the field line integration from two orthogonal directions on P_c, does not generally implies that the corresponding footpoints on P₁ and P₂ will form an orthogonal system. Indeed, in a QSL, a square on P_c is deformed to possibly very elongated non-orthogonal quadrangles on P₁ and P₂. At best, only one of the quadrangle diagonals can be oriented along one of the axes of a local system of coordinates, while the other diagonal will always be tilted with respect to the second coordinate axis. The non-orthogonal diagonals of the quadrangle, or their projection to the appropriate axis, are not proper jacobian elements, because each differential of independent variables (e.g. d_x1 or d_y1) of partial derivatives must be taken at a constant value of the other variable (y1 or x1, respectively). As the mapping of the field line integration grid is strongly distorted in a QSL, this leads to important errors on the computation of Q.

Finally if one compute field lines closer and closer from r_c, using smaller and smaller values of δ, we do not necessarily converge on a more precise value of Q. Indeed the ratios, such as d_X2/d_x1, are not necessarily more precisely estimated and the above problem, of non orthogonal derivative directions, remains. The estimation of the Jacobian does not converge with a smaller mesh size on P_c under the approximation done.

The overall issues related to this method are related to the fact that one tries to compute D₁₂ (and D₂₁) using Eq. (1) (resp. Eq. (2)) from a grid defined on P_c while Eq. (1) is meant to be properly evaluated from a neighborhood of r₁ (resp. r₂).

Fig. 3 Cartoon illustrating the second method to compute Q at rc.

3.2. Method 2

A more suitable method consists in first integrate the field line passing through P_c, then choose one of its footpoint, and compute Q using Eq. (8) or (9) from a set of field lines originating from a neighborhood of this footpoint. The derived value of Q is then attributed to r_c as Q is invariant along a field line. Figure 3 illustrates this method using the footpoint on P₁. The Jacobian matrix D₁₂ is evaluated from a set of field lines integrated from a distance δ along two orthogonal directions, in the plane P₁. The components of the field lines footpoint distance on P₂, as shown in Fig. 3, d_X2x1, d_Y2x1, d_X2y1 and d_Y2y1, are used to derive Q in Eq. (8) approximated such as: (11)In this method, as Q is not directly computed in r_c, field line integration errors may lead to some miss-location of the actual values of Q. This will tend to broaden the width of the QSL in region where QSL are particularly thin, i.e. Q is especially large. This is illustrated in Fig. 5, lower panels where Q is estimated from each footpoint in the two panels. Compared to the correct value of Q (upper right panel, cf. following section), Q presents a broader structure at the highest values of Q.

One also note some differences between the computation done from each footpoint. This asymmetry also results from field integration differences between footpoints. Our field integration is based on the D02CJF routine of the NAG library for which we have used a value of TOL equal to 10^-8. It means that the location of the footpoints are known to better than the 8th decimal when computing field line from one footpoint to the other (this can be checked by computing the field line back to the initial footpoint). Our field line computation is therefore particularly robust. However, as in the core of the QSL field lines diverge particularly strongly, even small integration errors leads to differences when computing Q from one footpoint or the other.

Another issue with this method is that the precision depends on δ evaluated on P₁ (or P₂) and not on P_c. It is not directly possible to converge on r_c, to use a more precise grid over P_c.

Fig. 4 Cartoon illustrating the third and most proper method to compute Q at rc.

3.3. Method 3

In order to properly compute Q at r_c, we will use field lines computed from P_c in a local neighborhood of r_c. But instead of directly trying to evaluate Eqs. (8, 9), we will derive the expression of Q using the Jacobian matrices derived from the mapping of the field from P_c to P₁ and P₂. We can now consider the mappings from r_c to the corresponding footpoint of the field line: Π_c1r_c → r₁; Π_c2r_c → r₂ and the inverse mapping Π_1c & Π_2c from P₁ & P₂ to P_c: Π_1cr₁ → r_c ; Π_2cr₂ → r_c (see Fig. 4).

Fig. 5 Distribution of Log10   Q in the plane y = 0 for the TD model using the different methods computed with 5122 points. Upper left: pseudo-method 1. Upper right: method 3 (proper method). Lower left: method 2 computed from the first footpoints. Lower left: method 2 computed from the second footpoints. A color version is available online.

By composition of the function we have: Π₁₂ = Π_c2°Π_1c, hence the Jacobian matrix D₁₂ can be expressed relatively to the Jacobian matrix D_1c and D_2c: (12)The jacobian matrix D_1c is equal to: (13)with (x_c,y_c) the cartesian coordinates of r_c in P_c. As one start computing field lines from the plane P_c this quantity can be difficult to evaluate. It is more precise to use the Jacobian matrix D_c1. As Π_1c = (Π_c1)^-1 on have: (14)with Δ_c1 the determinant of the matrix D_c1, which is such as (15)with B_n,c(x_c,y_c) the value of the magnetic field component normal to P_c at r_c.

Hence, on can derive a form for the Jacobian matrix D₁₂ which only depends on derivative taken along the plane P_c: (16)with (17)The norm N₁₂ is thus given by: (18)with the elements of D ^∗ . The squashing degree is thus: One verify that if one consider that the plane P_c is equal to the plane P₁ (resp. P₂) then Eq. (16) consistently gives Eq. (1) (resp. computed by Eq. (2)) and the classical values of N and Q are found. Titov (2007) has derived a covariant form of Q for any shapes of boundaries and/or coordinate systems. The present equations can similarly be generalized following the concepts presented in that paper.

Practically, to estimate Q from Eq. (19) one needs to compute the footpoints on P₁ and P₂ from a set of field lines integrated from a distance δ along two orthogonal directions in P_c around r_c. The components of the field lines footpoint distance on P₁ (resp. P₂) are as shown in Fig. 4: d_X1xc, d_Y1xc, d_X1yc and d_Y1yc (resp. d_X2xc, d_Y2xc, d_X2yc and d_Y2yc). They are used to derive the norm of the jacobian N₁₂ in Eq. (18), then Q in Eq. (20) is approximated such as: (21)with (22)for the TD model the result of the estimation of Eq. (21) is displayed in the upper right panel of Fig. 5. Method 3 corrects the problems encountered previously with method 2 (broadening and asymmetry of the Q distribution). The hyperbolic flux tube (HFT, region of highest Q values Titov et al. 2002) is clearly identified (especially with the color version of the figure where the dynamic of the plot is larger). Figure 6 presents two other examples of 2D Q maps obtained in 2 distinct plane of the 3D studied domain. The location and shape of the HFT is also clearly shown in these 2D cuts.

It is also straightforward to obtain a precise computation of Q over P_c by using smaller and smaller values of δ. This iteration with smaller δ values is very important since the numerical grid used to compute Q needs to resolve its sharp variations and δ should be adapted to the local mapping gradient. This implies that the final δ used is a sharp function of the position: δ should be smaller where Q is larger (e.g. the QSL thickness was found to scale approximately as 1/N in various magnetic configurations, Démoulin et al. 1996a,b). Equation (21) and the corresponding methodology as been used in previously published papers (Aulanier et al. 2005, 2006; Pariat et al. 2006; Masson et al. 2011) but with only a very short summary of the procedure and no comparison with other methods.

Fig. 6 Distribution of Log10   Q in the planes x = 0 (top panel) z = 0.229 (bottom panel) and for the TD model, with the photospheric plane z = 0 being the reference boundary for the computation. A color version is available online.

3.4. Comparison of method 3 with a previous procedure

Recently, Savcheva et al. (2012a) presented a original method to derive Q in cross-section of the 3D domain which presents similarities to method 3. Instead of determining D^∗, as performed here, they directly compute D_c2 and D_c1 (using the formulation introduced in this study). They first consider six points at a distance δ around r_c, two along each cartesian direction, and they compute their footpoints on P₁ and P₂. The position differences of the points are projected on a plane, here called , taken locally perpendicular to the field line passing by r_c (this orthogonal condition was implicitly used in Savcheva et al. 2012a; while not described, Savcheva, priv. comm.). The use of 6 points ensures that at least 4 of the projected points are significantly different. The distances of the six projected points are then related to the relative positions of the six footpoints on the plane P₁ (and P₂). A least square fit of these relationships with a matrix having constant coefficients provides D_c1 (resp. D_c2). D_1c is then inverted as in Eq. (14). The squashing degree is then determined with Eqs. (5, 8). The use of 6 points, while overdetermining the system to be solved, allow to estimate the precision of the matrix coefficients.

However, let us note that the field lines passing by the six original points (outside ) are not exactly the same as the field lines passing by the projections of these points on . Therefore the matrices which are computed are not strictly the ones which correspond to the Jacobian matrices of the mapping from to P₁ and P₂. As is locally orthogonal to the field in r_c, the difference is relatively small. It can nonetheless be significant in region of high Q, i.e. where the connectivity change significantly. A possible way to improve this would be to recompute the location of the field line footpoints for the points projected on P_c. Since in Savcheva et al. (2012a) no convergence is performed around r_c, their method is in any case limited in its ability to compute high value of Q. They indeed limit the value of the computed Q to 100. The above effect of projection is thus likely to be insignificant in their study.

Our previous derivation of Sect. 3.3 does not apply for the particular points where the field lines are tangent to the plane P_c. Indeed when a field line is strictly tangent to the plane of cut, the mapping cannot be determined. Practically, on discrete mesh, this is an extremely seldom situation. However, there are still some cases where the magnetic field is nearly parallel to P_c. In such cases, the computation of Q can be less accurate. Savcheva et al. (2012a) solved this problem by computing Q on the plane locally perpendicular to the magnetic field (Savcheva, priv. comm.).

Method 3 can be modified by using a similar local plane, . This suppress the problem of field lines tangent to P_c. The neighboring points being defined after , this allows to avoid the need to introduce points outside , and the effect of projection on (as in Savcheva et al. 2012a). This procedure break down only in the close vicinity of a magnetic null point where the computation of Q is already not possible (since it diverges to infinity). Thus, the use of a local plane to compute Q at any point r_c, permits to have 2D cuts (not limited to a plane) as well as 3D computations of Q within the volume V. Still, in many applications, it will be sufficient, and indeed recommended, to select first the plane P_c crossing the QSL in order to better compare the Q distribution to plasma parameters (e.g. electric currents, plasma density).

4. Application to an observed configuration

The TD model is an ideal magnetic configuration with a smooth magnetic field and only one QSL. In solar applications, the number of photospheric magnetic and electric current polarities implies a much more complex coronal magnetic field even in force-free extrapolations. Such complexity can also be present in MHD simulations. We present below an example of such complex field which enhance the differences between method 2 and 3.

We have performed a QSL computation on a non-linear force-free (NLFFF) reconstruction of an observed sigmoid (see McKenzie & Canfield 2008; Savcheva & van Ballegooijen 2009; Savcheva et al. 2012a,for more detailed information about the observational properties of this region). A flux rope insertion method (van Ballegooijen 2004) has been used to obtain the 3D magnetic structure of this region (this is explained in details in Savcheva et al. 2012a). The time of the model, 06:41 UT on February 12th 2007, precedes a solar eruption by about one hour. Unlike the TD configuration, the NLFFF model presents a high degree of complexity. The NLFFF model is indeed based on an observed magnetogram in which the field has various structures in a broad range of spatial scales (cf. Figs. 1 and 2 of Savcheva et al. 2012a). Numerous very thin QSLs are therefore present in the domain (Savcheva et al. 2012a,Sect. 5).

Figure 7 displays a vertical 2D cut in the middle of the inserted flux rope, roughly perpendicular to its axis. While in the case of the TD model differences between method 2 and 3 were relatively small (Fig. 5), and mostly located at the HFT, we remark here multiple differences. The method 2 is strongly affected by the complexity of the magnetic field and it tends to broaden significantly the QSLs. These errors will significantly affect the precise comparison between the location and width of the QSLs with other quantities (e.g. current densities).

Fig. 7 Distribution of Log10   Q in a vertical cut in the center of a NLFFF reconstruction of a sigmoid region observed on February 12th 2007 at 06:41 UT. Left: method 2 computed from the positive footpoints. Right: method 3 (proper method). A color version is available online.

Moreover, with the method 2, we also note a spotted distribution of Q while the distribution is smooth with the method 3 (Fig. 7). This is present at many locations within QSLs, and it is better seen in the upper left part of the Q distribution (around x_c = v ≈  − 20, y_c = z ≈ 40, and better seen on a zoom of the color representation of Q). This is another effect of the errors due to the transport of Q values along field lines. This effect is stronger as the magnetic field is spatially more complex.

One also observed relatively important differences in the highest values of Q obtained in the QSLs between method 2 and 3 (Fig. 7). As noted in Fig. 8 of Savcheva et al. (2012a), an HFT is present in the domain. This region should have the highest value of Q. However similarly to the TD case, the method 2 also failed at attributing the correct values at the HFT. This is because the convergence steps needed to properly resolve the connectivity gradients is not consistently done with the method 2. With the method 3, the HFT is directly identified and Q reaches value of the order of 10²² there.

Finally, let us note that the cut presented here was done in the vicinity of the cut presented in Fig. 7 of Savcheva et al. (2012a). In the latter, the absence of convergence, the relatively low resolution and the maximum thresholding value of 100 for Q do not allow an in-depth and quantitative study of the QSLs, even though similar structures are globally recognizable. A more detailed study of the magnetic topology of Savcheva et al. (2012a), in comparison with the TD model and the evolution of an MHD numerical simulation (Aulanier et al. 2010), using the method 3 presented here, will be presented in a forthcoming paper (Savcheva et al. 2012b).

5. Conclusion

In the present study we have analyzed some issues related to the computation of the squashing degree Q inside a 3D domain. We have analyzed three methods.

The first one is only a pseudo-method explicitly shown as an example of the possible tricks encountered when computing Q. Indeed, because such computation involves a combination of field line mapping, partial derivatives, and a broad range of spatial scales (with much smaller scales than for the magnetic field), the computation of Q is not straightforward. This is especially true for QSLs (i.e. where Q is large) since a tiny mesh size is required there to resolve the huge spatial gradients of the field line mapping.

The second method is an application of the original Q computation at the domain boundaries with a transport of Q along field lines within the domain. At a given position, it provides two values of Q: one transported from each magnetic polarity of the boundary.

The third method is designed to avoid the limits of the second method, in particular to avoid the errors associated to the mapping of the boundary Q distribution within the domain. For any point inside the domain, located at r_c, the mapping is computed towards the boundaries in both directions along field lines. We have derived an analytical expression for Q, Eq. (19), which needs to be substituted to the initial definition of Eq. (8). Numerically, the partial derivatives of both mappings are computed by finite differences with four points located in a plane in the close vicinity of r_c. Since the spatial scales involved in the mapping gradients can change by many orders of magnitude depending on the position of r_c, a local convergence is needed to find a grid finely adapted to resolve the gradients. Method 3 is indeed designed to more efficiently incorporate this convergence than method 2. Thus, with method 3, it becomes straightforward to represent the distribution of Q in a 3D domain and in particular on planar cuts through the domain.

We have illustrated our approach with two examples. The first one is a Titov & Démoulin (TD 1999) configuration which has a flux rope in equilibrium within an active-region like magnetic field. It has a smooth magnetic field with a main QSL. Figures 5, 6 presents examples of 2D Q maps obtained in distinct planes of the 3D studied domain. The location and shape of the HFT is clearly shown in this particular configuration. The second example is derived for the modeled magnetic configuration of an observed sigmoidal region (Savcheva et al. 2012a). Since the observed magnetogram has many polarities of various sizes, the computed coronal field has much more structures than in the TD configuration. Accordingly, the QSL pattern is complex (Fig. 7). It also contains an HFT, similarly to the TD configuration.

The method 2 permits a relatively precise computation of Q, when a convergence to small scales is implemented, and allows to determine large Q values (e.g. compare our Fig. 7 to Fig. 8 of Savcheva et al. 2012a). Differences are nonetheless present even in a smooth analytical model: artificial broadenings of the QSL; unequal results depending on the footpoint used to compute Q (Fig. 5). Method 3 corrects the defaults detected in method 2 by computing Q where it is needed, both avoiding errors in the Q transport along field lines and permitting a proper convergence to small scales. The improvements are already visible for the TD model (Fig. 5), and more evidently for a complex solar magnetic field (Fig. 7). We conclude that method 3 provides a useful practical tool to compute and represent QSLs in numerical domains issued from analytical, numerical or observation data set. The results presented above, with Q computed on 2D planar cuts of a 3D domain, will be useful to precisely compare Q distribution to the distributions of other quantities (e.g. electric current density, plasma flows) and to investigate the role of QSLs in solar physics.

Acknowledgments

A. Savcheva is deeply thanked for providing the data of the NLFFF model, for detailed discussions, and for access to additional information about her method to derive the QSLs in her study. The authors thank the referee for helpful comments which improved the clarity of the paper. The QSL computations have been performed on the multi-processors TRU64 computer of the LESIA.

References

All Figures

	Fig. 1 Left panel: photospheric distribution (grayscale) of the vertical component of the magnetic field, Bz, for the TD model. The continuous (resp. dotted) isolevels of Bz corresponds to intensities of  ±  [1,2,3,4] . Right panel: photospheric distribution (color scale) of Log10   Q. A color version is available online.
In the text

	Fig. 2 Cartoon illustrating the first possible method, but incorrect, to compute Q at rc.
In the text

	Fig. 3 Cartoon illustrating the second method to compute Q at rc.
In the text

	Fig. 4 Cartoon illustrating the third and most proper method to compute Q at rc.
In the text

	Fig. 5 Distribution of Log10   Q in the plane y = 0 for the TD model using the different methods computed with 5122 points. Upper left: pseudo-method 1. Upper right: method 3 (proper method). Lower left: method 2 computed from the first footpoints. Lower left: method 2 computed from the second footpoints. A color version is available online.
In the text

	Fig. 6 Distribution of Log10   Q in the planes x = 0 (top panel) z = 0.229 (bottom panel) and for the TD model, with the photospheric plane z = 0 being the reference boundary for the computation. A color version is available online.
In the text

	Fig. 7 Distribution of Log10   Q in a vertical cut in the center of a NLFFF reconstruction of a sigmoid region observed on February 12th 2007 at 06:41 UT. Left: method 2 computed from the positive footpoints. Right: method 3 (proper method). A color version is available online.
In the text