On the Magnetic Squashing Factor and the Lie Transport of Tangents

The squashing factor (or squashing degree) of a vector field is a quantitative measure of the deformation of the field line mapping between two surfaces. In the context of solar magnetic fields, it is often used to identify gradients in the mapping of elementary magnetic flux tubes between various flux domains. Regions where these gradients in the mapping are large are referred to as quasi-separatrix layers (QSLs), and are a continuous extension of separators and separatrix surfaces. These QSLs are observed to be potential sites for the formation of strong electric currents, and are therefore important for the study of magnetic reconnection in three dimensions. Since the squashing factor, Q, is defined in terms of the Jacobian of the field line mapping, it is most often calculated by first determining the mapping between two surfaces (or some approximation of it) and then numerically differentiating. Tassev & Savcheva have introduced an alternative method, in which they parameterize the change in separation between adjacent field lines, and then integrate along individual field lines to get an estimate of the Jacobian without the need to numerically differentiate the mapping itself. But while their method offers certain computational advantages, it is formulated on a perturbative description of the field line trajectory, and the accuracy of this method is not entirely clear. Here we show, through an alternative derivation, that this integral formulation is, in principle, exact. We then demonstrate the result in the case of a linear, 3D magnetic null, which allows for an exact analytical description and direct comparison to numerical estimates.


Introduction
Magnetic reconnection, which is thought to be the primary mechanism for the release of magnetic free energy in the solar corona, is a process by which magnetic field lines change their connectivity, in violation of the frozen-in condition. In twodimensional (2D) reconnection, this process is restricted to occur only at magnetic nulls, and so the study of magnetic reconnection has naturally been connected to magnetic topology, with topological features such as nulls, spines, and separatrix surfaces playing key roles. For a review, see Longcope (2005). It has been observed, however, that in three-dimensional (3D) reconnection, it is possible for field lines to change their connectivity continuously within non-ideal regions that may or may not contain nulls Schindler et al. 1988) and it is, therefore, not entirely clear what the role of magnetic topology is in this context. What is clear is that reconnection requires a parallel electric field and that, given a finite electrical resistivity, electrical currents must mediate the process.
Despite the fact that reconnection can occur in regions where the magnetic field does not vanish (Priest et al. 2003;Pontin et al. 2005), the volumes around topological features in the magnetic field appear to be preferential sites for current accumulation. This has led to the identification of so-called quasi-separatrix layers (QSLs), which are magnetic flux tubes for which the gradients of the field line mapping exhibit large (even singular) values (Priest & Démoulin 1995;Démoulin et al. 1996Démoulin et al. , 1997Titov 1999). Because QSLs are themselves magnetic flux tubes, they can be described locally. And yet, the property that defines a QSL is a global metric, which depends on the choice of boundaries as well as integrated properties of the field. It is therefore possible to have an entire class of field line configurations with the same mapping, 1 which makes it difficult to form a direct connection between QSLs and specific sites of current accumulation.
To better understand the connection between QSLs and reconnection, such authors as Titov et al. (2003), Démoulin (2006), and Aulanier et al. (2005) have studied current formation in the case of a hyperbolic flux tube (HFT), with further efforts being made by Wilmot-Smith et al. (2009) in the case of braided fields and Janvier et al. (2013) in the case of a torus-unstable flux rope. Other investigations into QSLs and their role in magnetic reconnection include scaling laws for singular current layers (Craig & Effenberger 2014;Effenberger & Craig 2016), as well as laboratory plasma experiments (Gekelman et al. 2012;Lorenzini et al. 2016) and even consideration of magnetospheric plasmas, both computational (Restante et al. 2013) and observational (Wendel & Adrian 2013). In every case, two things are clear. First, there is not a one-to-one correspondence between QSL structures and current accumulation. This is fundamentally related to the nature of the field line mapping, which is a non-local property of the magnetic field. And, second, despite the non-locality of the mapping, QSLs continue to provide insights into regions where currents are likely to develop, and where, should reconnection occur, the overall structure of the field is likely to be dramatically changed.
The typical metric for studying the field line mapping is the Squashing Factor (Titov et al. 2002;Titov 2007), which is a geometrical measure that, nonetheless, reproduces topological structures in the limit of a discontinuous mapping, and the identification of QSLs is predicated on our ability to accurately calculate the squashing factor, given an arbitrary magnetic field and a well defined boundary. But while the definition of Q can be stated fairly unambiguously, the magnetic field itself is often sampled discretely, and even where an analytical field is used it is usually not possible to describe the mapping analytically, so implementation of the definition depends on how the mapping is constructed and subsequent approximations to its derivatives. In practice, this usually means choosing a representative sample of field lines that are, in some sense, close to each other, and then determining where (and whether) they intersect the surfaces on which the mapping is to be described. These locations are then used to construct a finite difference approximation of the Jacobian of the mapping, which then leads directly to an estimate of Q.
This method is relatively straightforward in principle; however, the quantitative result can vary significantly with only slight variations in its implementation, and it is not always clear whether the calculated value is representative of the underlying structure of the field, or some artifact of the implementation. This is particularly the case when calculating Q for an intermediate point along a field line, and it was this concern that led Pariat & Démoulin (2012) to perform a survey of different implementations of the technique, in which they reviewed the adverse effects that arise from various seemingly benign assumptions before finally offering a "best practices" method, therein called "method 3." The idea behind this method is to represent the field line mapping as a composition of two separate mappings, each from a boundary surface to a shared, intermediate surface. The full Jacobian is then constructed algebraically from its constituent parts, and due to the symmetry of the method, it does seem to be an improvement over earlier procedures.
Given a particular method for calculating Q, it remains a computationally nontrivial task to choose field lines in the vicinity of each point of interest, trace each of these to its end points, and then numerically construct the local value of Q, a procedure that must be repeated for each point at which Q is to be calculated. In an attempt to alleviate this computationally intensive task, Tassev & Savcheva (2017) developed an alternate method, which is based on the construction in Pariat & Démoulin (2012), but differs in its approach to calculating the Jacobian elements, employing a transport technique similar to that of Longcope & Strauss (1994). This reduces the computational load by 40% from the Pariat method-where the latter must independently trace five 3D field line trajectories, for a total of 15 scalar elements, the former only traces one field line, plus a pair of additional tangent vectors for a total of nine scalar elements. Furthermore, Tassev & Savcheva (2017) claim that their method allows for significant relaxation of the precision requirements for field line integration, leading to further reduction in computing time.
Yet, while the method of Tassev & Savcheva (2017) is computationally efficient, its formulation is based on a perturbative representation of the magnetic field, so it is not entirely clear what we should expect from this method in terms of accuracy. On the one hand, because the field line mapping is never stated explicitly, their method avoids the certain complications that can arise in situations where adjacent field lines are traced to different boundary surfaces. And, as has been demonstrated through their benchmarking of the QSL Squasher software (Tassev & Savcheva 2017), the perturbative approximation seems to agree very closely with similar benchmarks from, e.g., Pariat & Démoulin (2012). However, it remains to be seen whether this approximation will remain valid in the vicinity of, e.g., a magnetic null, where the field lines are known to have self-similar solutions that span multiple scales, and where the field line mapping exhibits a discontinuity.
The aim of this study is to place the method of Tassev & Savcheva (2017) on a firm theoretical footing and to extend the applicability of the method. The organization of this paper is as follows. In Section 2, we develop a formalism for extracting information about the Jacobian from the behavior of local coordinate tangents, which are transported along the field. This leads directly to a procedure (which we summarize in Section 2.4) for calculating the squashing factor, Q, for arbitrary points within a computational domain, through the identification of each point with a field line and its associated mapping. Then, in Section 3, we consider the familiar example of a linear, 3D magnetic null, for which Q can be calculated analytically using the Lie transport of tangents, and we show that the result is identical to the result obtained by Pontin et al. (2016) through explicit differentiation of the mapping. We then go on to compare these results to numerical estimates obtained by QSL Squasher, before concluding in Section 4.

Transport Formulation
Because the notion of a field line is fundamental to the discussion of the squashing degree, we begin by examining field lines and their mappings. In particular, we consider certain properties of the Jacobian of the mapping and how these can be written in terms of elements of the tangent spaces of the mapped surfaces, without explicit reference to the mapping itself. We then describe how the mapping can be understood in terms of flows of coordinate trajectories, and how the basis vectors associated with these coordinate flows are transported along field lines. And since the coordinate flows are, themselves, tangent to the mapped surfaces, we may then associate the transport of these basis vectors along field lines with the transformation of elements of the tangent space under the field line mapping, and it is this association that is the underpinning of the method of Tassev & Savcheva (2017).
Consider a Cartesian coordinate system with a position vector, x x y z , , = ( ), and a generic vector field, . We shall denote a field line of B that passes through the point p as a smooth curve, X p s; B ( ), satisfying where s is a parameterized position along the curve, and is associated with the directional derivative (of a generic function, h) along B through We refer to the collection of all such curves as the flow of B, which we denote simply as X s B ( ). Since Equation (2) is an initial condition, we are free to set s to be zero at any point we choose; however, once the choice is made, the value of s at every other point along the same field line is then determined by Equation (1). Therefore, given a smooth, simply connected surface transversal to B, which we call S 0 , we may freely assign the value s=0 to every point within this surface. Then, provided that the domain of interest is sufficiently small so that field lines emanating from the surface S 0 do not return to it, we can define a mapping . 5 That is, for every point x in the domain, there is an associated value x s( ) and a point of origin ) .

The Jacobian
If the flow of B is smooth, then it sweeps out a volume extending from the surface S 0 , and this volume is foliated by isosurfaces of the monotonically increasing (differentiable) parameter x s( ), and each of these isosurfaces, which we refer to collectively as the surfaces S s , will be diffeomorphic to S 0 . Let us consider the mapping from S 0 to a particular value of , which defines the "forward" surface S f , and the "forward" mapping p X p S S s : ; .
Denoting the location of a given field line by local coordinates y i and Y j in the S 0 and S f surfaces, respectively, the dependence Y y j i ( ) is given explicitly by f f , whose Jacobian is then given in index notation by From the Jacobian, we can then form the Jacobian norm, which is related to the trace of the symmetrized Jacobian as which then leads to the squashing factor, as given by Titov et al. (2002), where  | | is the determinant of the Jacobian. While traditional methods for calculating Q rely on numerical estimates of the various terms in Equation (7), we wish to characterize the Jacobian without explicitly solving for the dependence in Equation (6), and we can do this by inspection of the push-forward, which relates vectors in the tangent spaces of S f and S 0 . Consider two vector fields, U 0 and V 0 , that define a locally orthonormal coordinate basis 3 on S 0 . Because U 0 and V 0 are elements of the tangent space of S 0 , they can be associated with corresponding elements of the tangent space of S f , and their transformation properties are then governed by the mapping f f , so that , 1 0 as seen in Figure 1. If we take the local coordinates in S 0 and S f to be Euclidean, we can compute the length of these vectors in the usual way, and this will reveal information about the deformation of the mapping. In matrix notation, the inner product of U f with itself is Similarly, taking V 0 to be another basis vector in S 0 , we find that From these relations, we can then construct the Jacobian norm, and from that Q, but before we do that let us consider the more general case, which consists of a composition of mappings, 4 one from S 0 to S f , denoted f f , and another from S 0 to a "backward" surface S b , denoted b f . We shall assume that each is invertible, and that their composition is a mapping from S b to S f , given by the transition map S S : .
The Jacobian of this transition map is then the matrix product of the Jacobian of each of the mappings, With these definitions, and a bit of algebra, the Norm of the Jacobian for the composite mapping takes the form which, by the properties of the trace, is then We can then make use of the fact that while the inverse of the symmetrized backward Jacobian is Combining Equations (18)- (20), we then find that and, from this expression, we can then construct the complete squashing factor for the mapping from S b to S f as where we have identified the determinant of each Jacobian as the ratio of the normal components of the magnetic field in S 0 , S b , and S f , 5 i.e.,

B n B n B B 23
This expression for the squashing factor in terms of vectors that are transformed under the field line mapping is exactly the expression that appears in Equation(11) of Tassev & Savcheva (2017), after associating the basis vectors U and V with their "field line displacement" vectors a and b. However, where their expression was derived from consideration of the covariance of Q and explicit inspection of the Jacobian elements given in Pariat & Démoulin (2012), using the local coordinate system in S b and S f , this formulation makes no explicit reference to any particular coordinate system. And while their description was particular to displacements perpendicular to the field, our formulation places no requirement on the orientation of the surfaces S 0 and S f with respect to B, and can therefore be applied to mappings of arbitrary orientation.

Transport of Tangents
Given that the squashing factor can be constructed from the transformation properties of vectors in the tangent spaces of S 0 , S f , and S b , it remains to determine these transformation properties explicitly. In the previous section, we considered the vectors U 0 and V 0 , which form a locally orthonormal coordinate system in S 0 . We shall now consider the flows generated by these vectors, which we call X p Then, because U 0 and V 0 span a coordinate system, their flows necessarily commute, i.e., . 25 For convenience, we will refer to the composition of these flows as X p , ; We can generalize the definition of U 0 and V 0 (which exist only in S 0 ) to a pair of vector fields, U and V, defined continuously along the flow of B, so that U , and similarly with V. This is done by requiring that . 27 Equations (26) and (27), together with the commutation of flows of U and V in S 0 , and the Jacobi identity, guarantee all three flows commute, so that the composition of these flows,  =| | | |, and similarly for f  | |, a result obtained from consideration of the transformation properties of the area twoform associate with U V , which shows that the elements of the transformed tangent space encode all of the information necessary for calculating N bf 2 , and therefore Q.
is independent of the ordering of U, V, and B, as demonstrated in Figure 2. These properties imply that the Lie bracket for the various pairings of U, V, and B must all vanish, so that The foregoing conditions simultaneously guarantee that , m n { } are a viable coordinate pair in any isosurface of s, and that their value is preserved along any field line of B so that, in fact, the triplet of values s, , m n { } is a viable coordinate system for the entire volume swept out by the flows of U, V, and B. And from these conditions, we can develop a prescription for constructing U and V throughout the domain, by taking as initial conditions the values of U 0 and V 0 , given at a point p S 0 0 Î , and then integrating along the flow of B, with And, just as we did with the field line parameter s, we can map values of m and ν onto the entire domain, of which values are preserved along flows of each other's basis vectors so that Thus, it follows that for any isosurface of s, given by S s , U and V will locally be elements of the tangent space of S s , and must transform accordingly. This connection makes exact the transport condition described in Tassev & Savcheva (2017) and the connection between the transport of tangents and the construction of Q in terms of said tangents. Namely, since U and V are basis vectors for a coordinate system whose flows commute with the field line parameterization, they must evolve according to the transport equation along a given field line. And since they are also elements of the tangent space for isosurfaces of the field line parameter, they must transform according to the Jacobian of the field line mapping. The Jacobian is, therefore, fundamentally related to the transport equation, and its various properties can be understood by inspection of the transport of tangents along field lines.

Rotation of the Tangent Plane
So far we have shown how the coordinate unit vectors U and V transform between surfaces S 0 and S f , and how their transformation properties relate both to the Jacobian Norm, and the continuous transport of tangents along field lines of B. For these results to be applicable to mappings between arbitrary surfaces of interest, we must be able to specify the orientation of the tangent plane at the point p f (or p b in the case of the backward mapping) so as to be consistent with the orientation of the surface that is being mapped. However, because the surface S f is necessarily an isosurface of s, its surface normal is defined by x s  ( ). We must therefore consider alternate parameterizations of the flows of B, which will allow for arbitrary orientation of the tangent plane, and whose elements will be related through a simple transformation rule.
Due to the linearity of Equations (1) and (3), the choice of parameterization for the field lines of B is arbitrary, up to a smooth, non-zero scalar field. That is, if we consider the field lines of B x x B x l = ( ) ( ) ( ), these are given by hich are identical to field lines of B, but for the (now more general) field line parameter t, which is given implicitly as X p s t s; . 38 x X p B s; Then, since t(s) increases monotonically with s along X p s; B ( ), the flows of B and B are related by a simple substitution rule, which we can use to generalize the results of the previous section. Because t(s) depends on the spatial variation of x l ( ), isosurfaces of t will themselves depend on x l ( ), and will generally be transversal to isosurfaces of s. And since the field lines are independent of the parameterization, the mappings generated by flows of B are no less valid as descriptors of the magnetic field than those generated by B. Therefore, if we wish to characterize the mapping to a surface with an arbitrary orientation, n f , at the point p f , it suffices to consider the parameterization given by a particular choice of x l ( ) such that, at the point p f , the prescribed surface normal, n f , is parallel to t  . This will ensure that the point p f is an element of an isosurface of t, called S f , whose tangent plane is orthogonal to n f at p f . But how are we to choose x l ( ) without prior knowledge of t  ? As it happens, we do not need to; it is enough to know that such a choice could be made, and transformation properties under the mapping defined by flows of B can be recovered from the mapping defined by flows of B.
Consider the generalized coordinate vectors Ũ and Ṽ , which will be equal to U and V in S 0 , and will commute with B. Since the flows of B are tangent to the flows of B, which are themselves perpendicular to contours of μ and ν, we find that the coordinates , m n { } are unaffected by the new parameterization, 6 so that ust as before. We wish to evaluate Ũ explicitly at the point p f , which will lie at the intersection of the isosurfaces S f and S f , each being defined by the values p s s f = ( ) and . 46 (46) is nontrivial to evaluate as it depends on the spatial variation of x l ( ), which is left general. Fortunately, we do not need to evaluate it directly; it is enough to know that it enters Equation (45) only as a coefficient for a term parallel to B, indicating that the component of Ũ that lies perpendicular to B is independent of x l ( ). Therefore, for every choice of x l ( ) (and therefore every choice of n f ), we can describe U f and, by extension, V f , as

The derivative in Equation
are given such that U f and V f are elements of the tangent space of S f , i.e.,

n U n V
0. 49 In practice, we can forego the explicit extraction of the perpendicular component, and transform directly between tangents for various parameterizations (see Figure 3) by Since U lies in an isosurface of s, while Ũ lies in S, which is an isosurface of some other parameterization, the latter can be constructed from the former through the addition of a component that is purely parallel to B. substitution of the above definitions into each, from which we see that if, for example, we calculate U f and V f by direct integration, we can then construct =^, while also satisfying Equation (49). And, by this construction, the mapped surface, S f , is guaranteed to be transversal to B, since no finite value of U G could make U f parallel to B unless U B 0 0 0 = , which is prohibited by our initial choice of U 0 .

Summary and Application
In the foregoing analysis, we have developed a formalism for initiating a pair of coordinate tangent vectors, evolving them along field lines, forming the projection of these tangents onto surfaces of arbitrary orientation, and then constructing the squashing factor for the mapping between those surfaces. We now include a short summary of these results in the form of an explicit procedure for calculating Q.
1. Choose a point, p 0 , at which Q is to be evaluated, and initialize a pair of orthonormal unit vectors (here called U 0 and V 0 ). These may be chosen to be orthogonal to B; however, the requirement is only that they not be parallel to B. 2. Trace a field line of B, forward and backward from p 0 to its end points at p f and p b , and simultaneously integrate Equations (32) and (33) for U and V, taking U 0 and V 0 as initial conditions and using whatever parameterization is convenient. 3. At the points p f and p b , choose local surface normals n f and n b , which define the local tangent space for the mapping to be considered. These may be taken to coincide with the criteria for the field line end points, e.g., the normal to the boundary, or the magnetic field direction, or any other local mapping surface of interest. 4. Reproject the integrated values of U and V at p f and p b using Equations (50) and (51) to recover Ũ and Ṽ , which are elements of the tangent space for the mapping into the prescribed surfaces, noting that the components of Ũ and Ṽ that lie perpendicular to B are the same as those of U and V. 5. Form the determinant and Norm of the Jacobian from the projected tangents (Ũ and Ṽ ), or from locally sampled values of B in the case of the determinant, and then algebraically construct Q for field lines in the neighborhood of p 0 , using Equation (22).

3D Null-An Example
In order to demonstrate the method that we have described, we now calculate the squashing factor for a linear null, and we compare our result with the exact expression derived by Pontin et al. (2016). We first align our coordinate system to the eigenvectors of the magnetic field Jacobian at the location of the null, and define the (nondimensional) magnetic field to be The field lines of B originate and terminate in the z=a and x=b surfaces, respectively, and we take these to be the boundaries of our domain. In order to avoid ambiguity, we employ local coordinates , c j ( ), corresponding to x y , ( ) in the z=a coordinate plane. This field configuration will be the same for both the analytical derivation as well as the numerical estimate that follows.

Exact Solution
For consistency with Pontin et al. (2016), we now consider the case where the mapped surfaces, S f and S b , are tangent to the boundaries of the domain (x = b and z = a) with surface normals n x f =ˆand n z b = -ˆ, respectively. With this convention for the orientation of the mapped surfaces, the associated squashing factor is referred to as Q ∂ . Since we are interested in Q ¶ within the z=a plane, which is cospatial with the "backward" surface (S b ), one of the composite mappings is the identity and this simplifies the calculation somewhat; however, as the general construction allows for such cases, no modifications to the described method are required.
We begin with the field line equation for the flows of B, which can be solved analytically to give and, therefore, Note that the value of s f depends explicitly on χ/b, and therefore changes within the x=b surface, which is, itself, not an isosurface of s. Therefore, any vector that is an element of the tangent space of s at z=a will not be mapped directly into the x=b surface, but will need to be transformed as described in Section 2.3 before it can be used to characterize the mapping between z=a and x=b. We now select an orthonormal pair of tangent vectors in the z=a plane, which we will integrate along field lines of B. In the general case, we would use the magnetic field to choose a preferred direction, and then construct these tangents to span a plane locally perpendicular to B, so that B U V 1 0 0 0 = |ˆ· ( )| . This is a perfectly valid approach, and we have carried out the calculation (which is quite lengthy) and confirmed that it produces the desired result. However, in the case of the linear null, it is convenient to choose our tangents to be initially parallel to the z=a plane, which is also the backward mapping surface. Then, since the backward mapping is just the identity, U 0 and V 0 are automatically equal to U b and V b .
To proceed with the calculation, let U x, 6 2 0 =ˆ( ) V y. 6 3 0 =ˆ( ) The transport equations for U and V are simplified greatly by the linearity of B, so that Equations (32) and (33) become where U i and V i are the Cartesian components of U and V, parameterized along X p s; B ( ), and k ij is the matrix of coefficients of B, from Equation (53). Since k ij is diagonal, the components of U and V evolve in exactly the same way as the components of the field lines themselves, with and similarly for V. From this, we find that, at s s f = , We will also need these quantities in the backward mapping surface. But, of course, the backward mapping is just the identity, and the S 0 surface is the same as the S b surface, so The associated squashing factor, as given in Equation (22), is then  Recalling that c and j are the field line locations in the z=a plane, we see that this expression is identical to the expression in Pontin et al. (2016), but has been derived here with no explicit differentiation of the field line mapping.

Numerical Estimate
Having validated the transport formulation analytically, we wished to test the numerical implementation on the same magnetic field configuration. And, where Tassev & Savcheva (2017) verified the QSL Squasher code for a sample magnetic field with large but finite gradients in the field line mapping, we wished to inspect the numerical result in the vicinity of a genuine separatrix structure, at which the value of Q diverges. We therefore used the QSL Squasher code to compute explicitly the perpendicular squashing factor for the same linear, 3D magnetic null described in the previous section.
The magnetic field was given as in Equation (52), with k=0.5, making the null symmetric in x y , . The computational domain was bounded by x y z 1 , , 1 -< < ( ) , corresponding to a b 1 = = , and although the QSL Squasher utility allows for the source magnetic field to be sampled on any rectilinear grid, the linearity of the field in this case negates any advantage that might be offered by increased resolution, so a coarse grid of 128, 128, 128 ( )was used for the source field. We performed the calculation both for a 3D volume rendering (see Figure 5), with a grid resolution of 256, 256, 256 ( ) , and for 2D surface renderings in the z=a and x=b planes, both with grid resolutions of 512, 512 ( ). For field line integration, we chose the default adaptive stepper (see Tassev & Savcheva 2017), with an error tolerance of 10 −4 . We also experimented with Eulerian stepping, and found similar results with a (surprisingly relaxed) step size of 10 −3 , but with a significant increase in computation time. Additionally, QSL Squasher allows for adaptive mesh refinement; however, in this case, we are interested in comparing the result as calculated at grid centers, so this feature was not employed. The results of these calculations are presented below, with relevant coordinates normalized to the value of a b = ( ). Since QSL Squasher only calculates the perpendicular squashing factor, a direct comparison of the numerical results to Equation (77) is not possible. However, the methods employed in the derivation of Q ¶ also allow for the calculation of the squashing degree for mappings with local coordinate tangents that are perpendicular to B (see Appendix for details of the calculation). For the following, we refer to the perpendicular squashing degree as Q^, while the numerical estimate of Q^from QSL Squasher is referred to as Q*. As in Section 3.1 these will be parameterized by the location of field lines that originate in the z=a plane (with local coordinates x y z a , , = ( ) given as , c j ( )) and terminate in the x=b plane (with x b y z , , = ( )given as (Y,Z)). In Figure 6, Q ¶ , Q^, and Q* are plotted for 0 j = , with 0 1 c < < . Note that while Q ¶ exhibits a local minimum at 2 2 c~, Q^has no local minimum, and decreases monotonically with c. Comparing Q* to Q^, the fidelity is extremely good all the way down to the sampling resolution of 10 3 c~-, which confirms the validity of the transport formulation as a numerical technique, as well as an analytical description. From inspection, we see that for 0 j = and small c the squashing degree goes as Q 2 ĉ -, which is identical to the limit of Q ¶ in Equation (77). This is to be expected since, for field lines that originate at small values of c, the magnetic field is nearly orthogonal to the boundaries at x=b and z=a, so the distinction between Q ¶ and Q^is immaterial in this limit.
In Figures 7 and 8 ) is fractionally larger than the theoretical minimum value of Q 2  . The left and right sides of Figure 7 indicate Q ∂ and Q^, as given in Equations (77) and (80), again with k=0.5 and a b 1 = = . The local minimum in Q ∂ at 2 2 c~is visible in the "C" shaped contour that twice intersects the 0 j = axis. Throughout the domain, the values of Q ∂ are generally higher than Figure 5. Perpendicular squashing factor, as calculated using QSL Squasher and rendered using Paraview, is shown for a 3D liner null. Q* is represented by the log 10 cubehelix color scale (Green 2011). Representative field lines are shown for the null spine and separatrix fan (purple), as well as self-similar field lines in the x=y and x y =diagonal planes (red and blue). those of Q^at the same locations, owing to the projection effect when mapping between surfaces that are not perpendicular to B; however, along the 0 j = axis the values of Q ∂ and Q^appear to converge for small c, as expected.
In Figure 8, we again see Q^in the right side of the figure, while the left side now shows Q*, calculated with QSL Squasher. For this calculation, the numerical domain was bounded by y 1  | | , so field lines that would otherwise intersect the x=b boundary at Y 1 > | | , instead intersect the y=b boundary. As a result (and because k=0.5 for this calculation), we see that Q* is symmetric about c j =  .
Despite this, the agreement between Q^and Q* can be readily seen by comparing their contours, which are nearly indistinguishable in the region c j > | |. Indeed, the only discernible difference appears to be in the slightly noisier quality of Q*, which is to be expected from a numerical result.
Since both the analytical expression for Q^and the (Y,Z) coordinates of a given field line in the x=b plane depend on the coordinates , c j ( ) in the z=a plane, these can be easily inverted to give Q Y Z ,( ), which can then be compared to the calculated value of Q* on the x=b surface. This calculation is depicted in Figure 9, which shows Q^, Q ¶ , and Q* as parameterized by Z with x b Y , 0 = = . As expected, the squashing degree diverges as Z goes to zero, which is consistent with the existence of a separatrix fan in the z=0 plane. The behavior of Q in this slice is qualitatively similar to that of Q near the spine, except for a weaker overall dependence on Z/a compared to χ/b-since Q goes as 2 cnear the spine, and Z 2 c , we find that Q∼Z −1 near the fan plane. The same local minimum in Q ¶ is in evidence, and the fidelity of Q* is extremely good.
A 2D representation of Q in the x=b surface is given in Figure 10, this time comparing only Q^and Q*. Again, we see excellent agreement between the analytical and numerical results. The slight upward tilt of the contours as Y | | increases is consistent with the same behavior for increasing j in Figures 7 and 8-field lines at larger values of Y are generally longer and their cross-sections more distorted than those that intersect the boundary near Y=0.

Conclusion
While the squashing degree (Q) of a vector field is a useful measure for characterizing certain topological and geometrical properties of the field, numerical estimates of Q are prone to systematic uncertainty and limited in resolution due to the computational intensity of the calculation. The method described and developed by Tassev & Savcheva (2017) goes some way to alleviating this through the use of GPU computing and also through the implementation of a novel formulation for calculating Q in terms of transported displacement vectors. But while the method has definite computational advantages, the theoretical validity of its formulation has not been previously described in detail.
Here we have confirmed the validity of the method of Tassev & Savcheva (2017) and shown that the squashing factor, Q, can be    derived exactly from the transport of coordinate tangents along field lines. This derivation is based on a differential geometry approach to the field line mapping, and is independent of any perturbative representation of the field. Furthermore, we have shown that the technique can be used for mappings of arbitrary orientation, and we have demonstrated this explicitly in the case of a linear null, enclosed in a cubic domain.
We have also shown that the numerical result, obtained by application of the QSL Squasher code on the same linear null, is in excellent agreement with the analytical solution. In particular, since Q is estimated without finite differencing of integrated values, the precision of Q is directly related to the numerical precision of the field line integrator, and the code clearly benefits from this scaling, as evidenced by the high fidelity of the numerical result, even for field lines that pass extremely close to the null.
We therefore conclude that the method described by Tassev & Savcheva (2017) is theoretically sound, and should not be expected to give erroneous results when applied to generic field configurations, although the numerical accuracy of the integration scheme must be monitored, as in all cases. Furthermore, we suggest that future implementations and refinements of this technique could be developed with the inclusion of user-defined definitions for the local surface normal of the mapping, as well as options for specifying the initial configuration for the integrated tangent vectors, which need not be orthogonal to the magnetic field, and might offer improvements to numerical accuracy if selected through some a priori knowledge of the magnetic field. + + -+ + - This lengthy expression can be simplified by considering the behavior along the x axis ( 0 j = ) for which we see that which has exactly the same limiting form as Q ∂ for small χ/b, namely,