Issue 
A&A
Volume 541, May 2012



Article Number  A78  
Number of page(s)  8  
Section  The Sun  
DOI  https://doi.org/10.1051/00046361/201118515  
Published online  01 May 2012 
Estimation of the squashing degree within a threedimensional domain
LESIA, Observatoire de Paris, CNRS, UPMC, Université Paris Diderot, 92190 Meudon, France
email: etienne.pariat@obspm.fr
Received: 24 November 2011
Accepted: 4 March 2012
Context. The study of the magnetic topology of magnetic fields aims at determining the key sites for the development of magnetic reconnection. Quasiseparatrix layers (QSLs), regions of strong connectivity gradients, are topological structures where intenseelectric currents preferentially buildup, and where, later on, magnetic reconnection occurs.
Aims. QSLs are volumes of intense squashing degree, Q; the fieldline invariant quantifying the deformation of elementary flux tubes. QSL are complex and thin threedimensional (3D) structures difficult to visualize directly. Therefore Q maps, i.e. 2D cuts of the 3D magnetic domain, are a more and more common features used to study QSLs.
Methods. We analyze several methods to derive 2D Q maps and discuss their analytical and numerical properties. These methods can also be used to compute Q within the 3D domain.
Results. We demonstrate that while analytically equivalent, the numerical implementation of these methods can be significantly different. We derive the analytical formula and the best numerical methodology that should be used to compute Q inside the 3D domain. We illustrate this method with two twisted magnetic configurations: a theoretical case and a nonlinear force free configuration derived from observations.
Conclusions. The representation of QSL through 2D planar cuts is an efficient procedure to derive the geometry of these structures and to relate them with other quantities, e.g. electric currents and plasma flows. It will enforce a more direct comparison of the role of QSL in magnetic reconnection.
Key words: magnetic fields / magnetic reconnection / magnetohydrodynamics (MHD) / Sun: magnetic topology
© ESO, 2012
1. Introduction
Magnetic reconnection in linetied lowplasmaβ environment such as the solar corona is linked to the formation of intense fieldaligned 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 quasiseparatrix 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 buildup (Milano et al. 1999; Galsgaard et al. 2003; Aulanier et al. 2005; Buechner 2006; Pariat et al. 2006; Masson et al. 2009; WilmotSmith 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 isosurface 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 activeregion like magnetic field. The normalized parameters are chosen such as: R = 2, a = 0.45, d = 1, L = 1, and q = 1, I_{0} = 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_{0} 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_{10} 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 nonlinear forcefree field constrained by both the photospheric magnetogram and the observed coronal loops. This magnetic field is associated to a longstanding 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, B_{z}, for the TD model. The continuous (resp. dotted) isolevels of B_{z} corresponds to intensities of ± [1,2,3,4] . Right panel: photospheric distribution (color scale) of Log_{10} 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 fieldline 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 photosphericlike boundary.
Let us consider a field lines which links the footpoint, r_{1}, of coordinate (x_{1},y_{1}) of the plane P_{1} to the footpoint r_{2}, of coordinate (X_{2},Y_{2}) of the plane P_{2} (see Fig. 2). The plane P_{1} will for example correspond to footpoint with a positive magnetic flux while P_{2} will correspond to a negative flux. Two mapping exist which associate a footpoint on one plane to the other: the mapping Π_{12} from P_{1} to P_{2}: r_{1}(x_{1},y_{1}) → r_{2}(X_{2},Y_{2}); and the inverse mapping Π_{21} from P_{2} to P_{1}: r_{2}(X_{2},Y_{2}) → r_{1}(x_{1},y_{1}).
The Jacobian matrix, D_{12} & D_{21} associated to each mapping are: As Π_{12} and Π_{21} are inverse functions (Π_{12}°Π_{21} = Π_{21}°Π_{12} = Id), D_{21} is related to D_{12} by: (3)with Δ_{12} the determinant of D_{12}.
Evaluating D_{21} 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 r_{c}. 
For example, when evaluating the squashing degree Q at r_{1}, it is better to integrate field lines from fixed points distant of δ_{x1} = δ_{y1} = δ and centered on r_{1}. The field lines, separated on P_{1} along the x_{1} direction by δ_{x1}( = δ), have a footpoint separation on P_{2} which component along the X_{2} and Y_{2} 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_{2}/∂x_{1} ≃ d_{X2x1}/δ will be evaluated more consistently than the quantities such as ∂x_{1}/∂X_{2}. The same is valid for field lines separated by δ_{y1} on P_{1}. We conclude that, when numerically evaluating Q and the mapping properties at r_{1}, Eq. (1) should be used to compute D_{12}, and D_{21} should be derived from Eq. (3) instead of Eq. (2).
Titov et al. (2002) demonstrated that the determinant of the Jacobian matrix, Δ_{12} & Δ_{21} can be related to the ratio of B_{z,1}(x_{1},y_{1}) and B_{Z,2}(X_{2},Y_{2}), which are respectively the normal component to P_{1} and P_{2} of the magnetic field evaluated in r_{1} and r_{2}: (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_{12} & N_{21} of the two Jacobian matrix are: the Squashing degree, Q, for this field line is given by: while in theory Q_{12} is equal to Q_{21} 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_{1} and P_{2} 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 keypoints 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_{1} and negative on P_{2}. 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 WilmotSmith 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. Pseudomethod 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_{12} or D_{21} from the position of the footpoints on P_{1} and P_{2} 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_{1} and P_{2}. These distances provides and information on how the flux tube is squashed. If one focus on the footpoint r_{1}, Q_{12} 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_{1} (d_{x1} and d_{y1} being computed in orthogonal directions). An equivalent relation can be obtained from the footpoint r_{2}. 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_{1} (field lines are computed from P_{c}), the Jacobian matrix D_{12} 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_{1} and P_{2} will form an orthogonal system. Indeed, in a QSL, a square on P_{c} is deformed to possibly very elongated nonorthogonal quadrangles on P_{1} and P_{2}. 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 nonorthogonal 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_{12} (and D_{21}) 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_{1} (resp. r_{2}).
Fig. 3 Cartoon illustrating the second method to compute Q at r_{c}. 
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_{1}. The Jacobian matrix D_{12} is evaluated from a set of field lines integrated from a distance δ along two orthogonal directions, in the plane P_{1}. The components of the field lines footpoint distance on P_{2}, 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 misslocation 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_{1} (or P_{2}) 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 r_{c}. 
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_{1} and P_{2}. We can now consider the mappings from r_{c} to the corresponding footpoint of the field line: Π_{c1}r_{c} → r_{1}; Π_{c2}r_{c} → r_{2} and the inverse mapping Π_{1c} & Π_{2c} from P_{1} & P_{2} to P_{c}: Π_{1c}r_{1} → r_{c} ; Π_{2c}r_{2} → r_{c} (see Fig. 4).
Fig. 5 Distribution of Log_{10} Q in the plane y = 0 for the TD model using the different methods computed with 512^{2} points. Upper left: pseudomethod 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: Π_{12} = Π_{c2}°Π_{1c}, hence the Jacobian matrix D_{12} 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_{12} which only depends on derivative taken along the plane P_{c}: (16)with (17)The norm N_{12} 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_{1} (resp. P_{2}) 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_{1} and P_{2} 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_{1} (resp. P_{2}) 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_{12} 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 Log_{10} 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 crosssection 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_{1} and P_{2}. 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_{1} (and P_{2}). 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_{1} and P_{2}. 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 forcefree 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 nonlinear forcefree (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 Log_{10} 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^{22} 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 indepth 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 pseudomethod 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 activeregion 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 multiprocessors TRU64 computer of the LESIA.
References
 Antiochos, S. K., Mikić, Z., Titov, V. S., Lionello, R., & Linker, J. A. 2011, ApJ, 731, 112 [Google Scholar]
 Aulanier, G., Pariat, E., & Démoulin, P. 2005, A&A, 444, 961 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Aulanier, G., Pariat, E., Démoulin, P., & DeVore, C. R. 2006, Sol. Phys., 238, 347 [NASA ADS] [CrossRef] [Google Scholar]
 Aulanier, G., Golub, L., DeLuca, E. E., et al. 2007, Science, 318, 1588 [NASA ADS] [CrossRef] [Google Scholar]
 Aulanier, G., Török, T., Démoulin, P., & DeLuca, E. E. 2010, ApJ, 708, 314 [Google Scholar]
 Bagalá, L. G., Mandrini, C. H., Rovira, M. G., & Démoulin, P. 2000, A&A, 363, 779 [NASA ADS] [Google Scholar]
 Buechner, J. 2006, Space Sci. Rev., 122, 149 [NASA ADS] [CrossRef] [Google Scholar]
 Chandra, R., Schmieder, B., Mandrini, C. H., et al. 2011, Sol. Phys., 269, 83 [NASA ADS] [CrossRef] [Google Scholar]
 Cohen, O., Kashyap, V. L., Drake, J. J., Sokolov, I. V., & Gombosi, T. I. 2011, ApJ, 738, 166 [NASA ADS] [CrossRef] [Google Scholar]
 Démoulin, P. 2005, Proceedings of the International Scientific Conference on Chromospheric and Coronal Magnetic Fields, ESA SP596, 30 August–2 September, 596, 22 [Google Scholar]
 Démoulin, P. 2006, Adv. Space Res., 37, 1269 [NASA ADS] [CrossRef] [Google Scholar]
 Démoulin, P., Priest, E. R., & Lonie, D. P. 1996a, J. Geophys. Res., 101, 7631 [NASA ADS] [CrossRef] [Google Scholar]
 Démoulin, P., Priest, E. R., & Mandrini, C. H. 1996b, A&A, 308, 643 [NASA ADS] [Google Scholar]
 Démoulin, P., Bagalá, L. G., Mandrini, C. H., Hénoux, J., & Rovira, M. G. 1997, A&A, 325, 305 [NASA ADS] [Google Scholar]
 Gaizauskas, V., Mandrini, C. H., Démoulin, P., Luoni, M. L., & Rovira, M. G. 1998, A&A, 332, 353 [NASA ADS] [Google Scholar]
 Galsgaard, K., Titov, V. S., & Neukirch, T. 2003, ApJ, 595, 506 [NASA ADS] [CrossRef] [Google Scholar]
 Gorbachev, V. S., & Somov, B. V. 1989, SvA, 33, 57 [Google Scholar]
 Longcope, D. W. 2005, Liv. Rev. Sol. Phys., 2, 7 [Google Scholar]
 Longcope, D. W., & Strauss, H. R. 1994, ApJ, 437, 851 [NASA ADS] [CrossRef] [Google Scholar]
 Lugaz, N., Downs, C., Shibata, K., et al. 2011, ApJ, 738, 127 [NASA ADS] [CrossRef] [Google Scholar]
 Mandrini, C. H., Démoulin, P., Rovira, M. G., & de La Beaujardière, J.F. 1995, A&A, 303, 927 [NASA ADS] [Google Scholar]
 Mandrini, C. H., Démoulin, P., Bagalá, L. G., et al. 1997, Sol. Phys., 174, 229 [NASA ADS] [CrossRef] [Google Scholar]
 Mandrini, C. H., Démoulin, P., Schmieder, B., et al. 2006, Sol. Phys., 238, 293 [NASA ADS] [CrossRef] [Google Scholar]
 Masson, S., Pariat, E., Aulanier, G., & Schrijver, C. J. 2009, ApJ, 700, 559 [NASA ADS] [CrossRef] [Google Scholar]
 Masson, S., Aulanier, G., Pariat, E., & Klein, K.L. 2011, Sol. Phys., 394 [Google Scholar]
 McKenzie, D. E., & Canfield, R. C. 2008, A&A, 481, L65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Milano, L. J., Dmitruk, P., Mandrini, C. H., Gómez, D. O., & Démoulin, P. 1999, ApJ, 521, 889 [NASA ADS] [CrossRef] [Google Scholar]
 Pariat, E., Aulanier, G., & Démoulin, P. 2006, in SF2A2006: Proceedings of the Annual meeting of the French Society of A&A, ed. D. Barret, 559 [Google Scholar]
 Roussev, I. I., Forbes, T. G., Gombosi, T. I., et al. 2003, ApJ, 588, L45 [Google Scholar]
 Savcheva, A., & van Ballegooijen, A. A. 2009, ApJ, 703, 1766 [NASA ADS] [CrossRef] [Google Scholar]
 Savcheva, A. S., van Ballegooijen, A. A., & DeLuca, E. E. 2012a, ApJ, 744, 78 [NASA ADS] [CrossRef] [Google Scholar]
 Savcheva, A. S., Pariat, E., van Ballegooijen, A. A., Aulanier, G., & DeLuca, E. E. 2012b, ApJ, submitted [Google Scholar]
 Titov, V. S. 2007, ApJ, 660, 863 [NASA ADS] [CrossRef] [Google Scholar]
 Titov, V. S., & Démoulin, P. 1999, A&A, 351, 707 [NASA ADS] [Google Scholar]
 Titov, V. S., Hornig, G., & Démoulin, P. 2002, J. Geophys. Res., 107, 1164 [Google Scholar]
 Titov, V. S., Galsgaard, K., & Neukirch, T. 2003, ApJ, 582, 1172 [NASA ADS] [CrossRef] [Google Scholar]
 Titov, V. S., Mikić, Z., Linker, J. A., Lionello, R., & Antiochos, S. K. 2011, ApJ, 731, 111 [Google Scholar]
 Török, T., Kliem, B., & Titov, V. S. 2004, A&A, 413, L27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Török, T., Aulanier, G., Schmieder, B., Reeves, K. K., & Golub, L. 2009, ApJ, 704, 485 [NASA ADS] [CrossRef] [Google Scholar]
 Török, T., Chandra, R., Pariat, E., et al. 2011a, ApJ, 728, 65 [NASA ADS] [CrossRef] [Google Scholar]
 Török, T., Panasenco, O., Titov, V. S., et al. 2011b, ApJ, 739, L63 [NASA ADS] [CrossRef] [Google Scholar]
 Valori, G., Kliem, B., Török, T., & Titov, V. S. 2010, A&A, 519, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van Ballegooijen, A. A. 2004, ApJ, 612, 519 [NASA ADS] [CrossRef] [Google Scholar]
 WilmotSmith, A. L., Hornig, G., & Pontin, D. I. 2009, ApJ, 704, 1288 [NASA ADS] [CrossRef] [Google Scholar]
 WilmotSmith, A. L., Pontin, D. I., & Hornig, G. 2010, A&A, 516, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Figures
Fig. 1 Left panel: photospheric distribution (grayscale) of the vertical component of the magnetic field, B_{z}, for the TD model. The continuous (resp. dotted) isolevels of B_{z} corresponds to intensities of ± [1,2,3,4] . Right panel: photospheric distribution (color scale) of Log_{10} Q. A color version is available online. 

In the text 
Fig. 2 Cartoon illustrating the first possible method, but incorrect, to compute Q at r_{c}. 

In the text 
Fig. 3 Cartoon illustrating the second method to compute Q at r_{c}. 

In the text 
Fig. 4 Cartoon illustrating the third and most proper method to compute Q at r_{c}. 

In the text 
Fig. 5 Distribution of Log_{10} Q in the plane y = 0 for the TD model using the different methods computed with 512^{2} points. Upper left: pseudomethod 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 Log_{10} 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 Log_{10} 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.