Free Access
Volume 541, May 2012
Article Number A78
Number of page(s) 8
Section The Sun
Published online 01 May 2012

© 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, I0 = 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 I0 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 Nt = 2 above z = 0. The photospheric distributions of Bz and Log10   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.

thumbnail 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, r1, of coordinate (x1,y1) of the plane P1 to the footpoint r2, of coordinate (X2,Y2) of the plane P2 (see Fig. 2). The plane P1 will for example correspond to footpoint with a positive magnetic flux while P2 will correspond to a negative flux. Two mapping exist which associate a footpoint on one plane to the other: the mapping Π12 from P1 to P2: r1(x1,y1) → r2(X2,Y2); and the inverse mapping Π21 from P2 to P1: r2(X2,Y2) → r1(x1,y1).

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

Evaluating D21 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.

thumbnail Fig. 2

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

For example, when evaluating the squashing degree Q at r1, it is better to integrate field lines from fixed points distant of δx1 = δy1 = δ and centered on r1. The field lines, separated on P1 along the x1 direction by δx1( = δ), have a footpoint separation on P2 which component along the X2 and Y2 directions are respectively equal to dX2x1 and dY2x1. In a QSL, dX2x1 and dY2x1 will typically be very large or very small relatively to δ. Because of numerical errors, the quantities such as X2/x1 ≃ dX2x1/δ will be evaluated more consistently than the quantities such as x1/X2. The same is valid for field lines separated by δy1 on P1. We conclude that, when numerically evaluating Q and the mapping properties at r1, Eq. (1) should be used to compute D12, and D21 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 Bz,1(x1,y1) and BZ,2(X2,Y2), which are respectively the normal component to P1 and P2 of the magnetic field evaluated in r1 and r2: (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 N12 & N21 of the two Jacobian matrix are: the Squashing degree, Q, for this field line is given by: while in theory Q12 is equal to Q21 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 P1 and P2 stays fixed and the aim is to determine the value of Q associated to an arbitrary point rc of V. We define a plane Pc represented by a cartesian referential such as rc has the coordinates (xc,yc). The plane Pc can be defined in the neighborhood of rc, Pc(rc), implying than its orientation could depend on rc (for example by taking perpendicular to the local magnetic field, see the end of Sect. 3.4). In this case, the local plane Pc(rc) can be different for different rc positions. However, for many applications, it may be more suited to consider an unique plane Pc, 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 rc, one first needs to determine the location of the two footpoints of the field line passing by rc. 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 P1 and negative on P2. 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 rc. 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(rc).

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 D12 or D21 from the position of the footpoints on P1 and P2 only. This intuitive methodology implies to first compute field lines from two directions on Pc around the point of interest rc as illustrated in Fig. 2. Four field lines are plotted from Pc distant from a fixed position δ centered on rc (δ being the distance between field lines position on each side of rc). One obtain the difference of position of the footpoints of these field lines, dx1, dy1, dX2 and dY2 on P1 and P2. These distances provides and information on how the flux tube is squashed. If one focus on the footpoint r1, Q12 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 Pc converges to P1 (dx1 and dy1 being computed in orthogonal directions). An equivalent relation can be obtained from the footpoint r2. 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 P1 (field lines are computed from Pc), the Jacobian matrix D12 is not properly estimated. Starting the field line integration from two orthogonal directions on Pc, does not generally implies that the corresponding footpoints on P1 and P2 will form an orthogonal system. Indeed, in a QSL, a square on Pc is deformed to possibly very elongated non-orthogonal quadrangles on P1 and P2. 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. dx1 or dy1) 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 rc, using smaller and smaller values of δ, we do not necessarily converge on a more precise value of Q. Indeed the ratios, such as dX2/dx1, 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 Pc under the approximation done.

The overall issues related to this method are related to the fact that one tries to compute D12 (and D21) using Eq. (1) (resp. Eq. (2)) from a grid defined on Pc while Eq. (1) is meant to be properly evaluated from a neighborhood of r1 (resp. r2).

thumbnail 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 Pc, 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 rc as Q is invariant along a field line. Figure 3 illustrates this method using the footpoint on P1. The Jacobian matrix D12 is evaluated from a set of field lines integrated from a distance δ along two orthogonal directions, in the plane P1. The components of the field lines footpoint distance on P2, as shown in Fig. 3, dX2x1, dY2x1, dX2y1 and dY2y1, are used to derive Q in Eq. (8) approximated such as: (11)In this method, as Q is not directly computed in rc, 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 P1 (or P2) and not on Pc. It is not directly possible to converge on rc, to use a more precise grid over Pc.

thumbnail 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 rc, we will use field lines computed from Pc in a local neighborhood of rc. 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 Pc to P1 and P2. We can now consider the mappings from rc to the corresponding footpoint of the field line: Πc1rc → r1; Πc2rc → r2 and the inverse mapping Π1c & Π2c from P1 & P2 to Pc: Π1cr1 → rc ; Π2cr2 → rc (see Fig. 4).

thumbnail 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: Π12 = Πc2°Π1c, hence the Jacobian matrix D12 can be expressed relatively to the Jacobian matrix D1c and D2c: (12)The jacobian matrix D1c is equal to: (13)with (xc,yc) the cartesian coordinates of rc in Pc. As one start computing field lines from the plane Pc this quantity can be difficult to evaluate. It is more precise to use the Jacobian matrix Dc1. As Π1c = (Πc1)-1 on have: (14)with Δc1 the determinant of the matrix Dc1, which is such as (15)with Bn,c(xc,yc) the value of the magnetic field component normal to Pc at rc.

Hence, on can derive a form for the Jacobian matrix D12 which only depends on derivative taken along the plane Pc: (16)with (17)The norm N12 is thus given by: (18)with the elements of D ∗ . The squashing degree is thus: One verify that if one consider that the plane Pc is equal to the plane P1 (resp. P2) 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 P1 and P2 from a set of field lines integrated from a distance δ along two orthogonal directions in Pc around rc. The components of the field lines footpoint distance on P1 (resp. P2) are as shown in Fig. 4: dX1xc, dY1xc, dX1yc and dY1yc (resp. dX2xc, dY2xc, dX2yc and dY2yc). They are used to derive the norm of the jacobian N12 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 Pc 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.

thumbnail 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 Dc2 and Dc1 (using the formulation introduced in this study). They first consider six points at a distance δ around rc, two along each cartesian direction, and they compute their footpoints on P1 and P2. The position differences of the points are projected on a plane, here called , taken locally perpendicular to the field line passing by rc (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 P1 (and P2). A least square fit of these relationships with a matrix having constant coefficients provides Dc1 (resp. Dc2). D1c 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 P1 and P2. As is locally orthogonal to the field in rc, 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 Pc. Since in Savcheva et al. (2012a) no convergence is performed around rc, 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 Pc. 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 Pc. 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 Pc. 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 rc, 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 Pc 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).

thumbnail 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 xc = v ≈  − 20, yc = 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 1022 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 rc, 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 rc. Since the spatial scales involved in the mapping gradients can change by many orders of magnitude depending on the position of rc, 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.


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.


  1. Antiochos, S. K., Mikić, Z., Titov, V. S., Lionello, R., & Linker, J. A. 2011, ApJ, 731, 112 [Google Scholar]
  2. Aulanier, G., Pariat, E., & Démoulin, P. 2005, A&A, 444, 961 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Aulanier, G., Pariat, E., Démoulin, P., & DeVore, C. R. 2006, Sol. Phys., 238, 347 [NASA ADS] [CrossRef] [Google Scholar]
  4. Aulanier, G., Golub, L., DeLuca, E. E., et al. 2007, Science, 318, 1588 [NASA ADS] [CrossRef] [Google Scholar]
  5. Aulanier, G., Török, T., Démoulin, P., & DeLuca, E. E. 2010, ApJ, 708, 314 [Google Scholar]
  6. Bagalá, L. G., Mandrini, C. H., Rovira, M. G., & Démoulin, P. 2000, A&A, 363, 779 [NASA ADS] [Google Scholar]
  7. Buechner, J. 2006, Space Sci. Rev., 122, 149 [NASA ADS] [CrossRef] [Google Scholar]
  8. Chandra, R., Schmieder, B., Mandrini, C. H., et al. 2011, Sol. Phys., 269, 83 [NASA ADS] [CrossRef] [Google Scholar]
  9. Cohen, O., Kashyap, V. L., Drake, J. J., Sokolov, I. V., & Gombosi, T. I. 2011, ApJ, 738, 166 [NASA ADS] [CrossRef] [Google Scholar]
  10. Démoulin, P. 2005, Proceedings of the International Scientific Conference on Chromospheric and Coronal Magnetic Fields, ESA SP-596, 30 August–2 September, 596, 22 [Google Scholar]
  11. Démoulin, P. 2006, Adv. Space Res., 37, 1269 [NASA ADS] [CrossRef] [Google Scholar]
  12. Démoulin, P., Priest, E. R., & Lonie, D. P. 1996a, J. Geophys. Res., 101, 7631 [NASA ADS] [CrossRef] [Google Scholar]
  13. Démoulin, P., Priest, E. R., & Mandrini, C. H. 1996b, A&A, 308, 643 [NASA ADS] [Google Scholar]
  14. Démoulin, P., Bagalá, L. G., Mandrini, C. H., Hénoux, J., & Rovira, M. G. 1997, A&A, 325, 305 [NASA ADS] [Google Scholar]
  15. Gaizauskas, V., Mandrini, C. H., Démoulin, P., Luoni, M. L., & Rovira, M. G. 1998, A&A, 332, 353 [NASA ADS] [Google Scholar]
  16. Galsgaard, K., Titov, V. S., & Neukirch, T. 2003, ApJ, 595, 506 [NASA ADS] [CrossRef] [Google Scholar]
  17. Gorbachev, V. S., & Somov, B. V. 1989, SvA, 33, 57 [Google Scholar]
  18. Longcope, D. W. 2005, Liv. Rev. Sol. Phys., 2, 7 [Google Scholar]
  19. Longcope, D. W., & Strauss, H. R. 1994, ApJ, 437, 851 [NASA ADS] [CrossRef] [Google Scholar]
  20. Lugaz, N., Downs, C., Shibata, K., et al. 2011, ApJ, 738, 127 [NASA ADS] [CrossRef] [Google Scholar]
  21. Mandrini, C. H., Démoulin, P., Rovira, M. G., & de La Beaujardière, J.-F. 1995, A&A, 303, 927 [NASA ADS] [Google Scholar]
  22. Mandrini, C. H., Démoulin, P., Bagalá, L. G., et al. 1997, Sol. Phys., 174, 229 [NASA ADS] [CrossRef] [Google Scholar]
  23. Mandrini, C. H., Démoulin, P., Schmieder, B., et al. 2006, Sol. Phys., 238, 293 [NASA ADS] [CrossRef] [Google Scholar]
  24. Masson, S., Pariat, E., Aulanier, G., & Schrijver, C. J. 2009, ApJ, 700, 559 [NASA ADS] [CrossRef] [Google Scholar]
  25. Masson, S., Aulanier, G., Pariat, E., & Klein, K.-L. 2011, Sol. Phys., 394 [Google Scholar]
  26. McKenzie, D. E., & Canfield, R. C. 2008, A&A, 481, L65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  27. Milano, L. J., Dmitruk, P., Mandrini, C. H., Gómez, D. O., & Démoulin, P. 1999, ApJ, 521, 889 [NASA ADS] [CrossRef] [Google Scholar]
  28. Pariat, E., Aulanier, G., & Démoulin, P. 2006, in SF2A-2006: Proceedings of the Annual meeting of the French Society of A&A, ed. D. Barret, 559 [Google Scholar]
  29. Roussev, I. I., Forbes, T. G., Gombosi, T. I., et al. 2003, ApJ, 588, L45 [NASA ADS] [CrossRef] [Google Scholar]
  30. Savcheva, A., & van Ballegooijen, A. A. 2009, ApJ, 703, 1766 [NASA ADS] [CrossRef] [Google Scholar]
  31. Savcheva, A. S., van Ballegooijen, A. A., & DeLuca, E. E. 2012a, ApJ, 744, 78 [NASA ADS] [CrossRef] [Google Scholar]
  32. Savcheva, A. S., Pariat, E., van Ballegooijen, A. A., Aulanier, G., & DeLuca, E. E. 2012b, ApJ, submitted [Google Scholar]
  33. Titov, V. S. 2007, ApJ, 660, 863 [NASA ADS] [CrossRef] [Google Scholar]
  34. Titov, V. S., & Démoulin, P. 1999, A&A, 351, 707 [NASA ADS] [Google Scholar]
  35. Titov, V. S., Hornig, G., & Démoulin, P. 2002, J. Geophys. Res., 107, 1164 [Google Scholar]
  36. Titov, V. S., Galsgaard, K., & Neukirch, T. 2003, ApJ, 582, 1172 [NASA ADS] [CrossRef] [Google Scholar]
  37. Titov, V. S., Mikić, Z., Linker, J. A., Lionello, R., & Antiochos, S. K. 2011, ApJ, 731, 111 [Google Scholar]
  38. Török, T., Kliem, B., & Titov, V. S. 2004, A&A, 413, L27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  39. Török, T., Aulanier, G., Schmieder, B., Reeves, K. K., & Golub, L. 2009, ApJ, 704, 485 [NASA ADS] [CrossRef] [Google Scholar]
  40. Török, T., Chandra, R., Pariat, E., et al. 2011a, ApJ, 728, 65 [NASA ADS] [CrossRef] [Google Scholar]
  41. Török, T., Panasenco, O., Titov, V. S., et al. 2011b, ApJ, 739, L63 [NASA ADS] [CrossRef] [Google Scholar]
  42. Valori, G., Kliem, B., Török, T., & Titov, V. S. 2010, A&A, 519, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  43. van Ballegooijen, A. A. 2004, ApJ, 612, 519 [NASA ADS] [CrossRef] [Google Scholar]
  44. Wilmot-Smith, A. L., Hornig, G., & Pontin, D. I. 2009, ApJ, 704, 1288 [NASA ADS] [CrossRef] [Google Scholar]
  45. Wilmot-Smith, A. L., Pontin, D. I., & Hornig, G. 2010, A&A, 516, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

All Figures

thumbnail 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
thumbnail Fig. 2

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

In the text
thumbnail Fig. 3

Cartoon illustrating the second method to compute Q at rc.

In the text
thumbnail Fig. 4

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

In the text
thumbnail 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
thumbnail 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
thumbnail 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

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.