Coiling and Squeezing: Properties of the Local Transverse Deviations of Magnetic Field Lines

TL;DR

This paper investigates the local transverse deviations of magnetic field lines, introducing new measures like the squeeze and coiling rates to better identify magnetic reconnection sites and flux rope structures.

Contribution

It develops a theoretical framework linking magnetic field line deviations to fluid flow concepts, proposing improved proxies for detecting QSLs and flux ropes.

Findings

The squeeze rate effectively identifies QSLs and current sheets.

The coiling number provides an unbiased estimate of field line twists.

The proposed measures outperform traditional metrics in locating magnetic reconnection regions.

Abstract

We study the properties of the local transverse deviations of magnetic field lines at a fixed moment in time. Those deviations "evolve" smoothly in a plane normal to the field-line direction as one moves that plane along the field line. Since the evolution can be described by a planar flow in the normal plane, we derive most of our results in the context of a toy model for planar fluid flow. We then generalize our results to include the effects of field-line curvature. We show that the type of flow is determined by the two non-zero eigenvalues of the gradient of the normalized magnetic field. The eigenvalue difference quantifies the local rate of squeezing or coiling of neighboring field lines, which we relate to standard notions of fluid vorticity and shear. The resulting squeezing rate can be used in the detection of null points, hyperbolic flux tubes and current sheets. Once integrated along field lines, that rate gives a squeeze factor, which is an approximation to the squashing factor, which is usually employed in locating quasi-separatrix layers (QSLs), which are possible sites for magnetic reconnection. Unlike the squeeze factor, the squashing factor can miss QSLs for which field lines are squeezed and then unsqueezed. In that regard, the squeeze factor is a better proxy for locating QSLs than the squashing factor. In another application of our analysis, we construct an approximation to the local rate of twist of neighboring field lines, which we refer to as the coiling rate. That rate can be integrated along a field line to give a coiling number, $\mathrm{N_c}$. We show that unlike the standard local twist number, $\mathrm{N_c}$ gives an unbiased approximation to the number of twists neighboring field lines make around one another. $\mathrm{N_c}$ can be useful for the study of flux rope instabilities, such as the kink instability, and can be used in the detection of flux ropes.

Equations

$$\displaystyle v_{i}({\mbox{\boldmath{$r$}}}(\tau))=\frac{dr_{i}(\tau)}{d\tau}=m_{ij}r_{j}(\tau)\ .$$

$$m=\Theta +\Psi =\Theta +\Sigma +\Omega ,$$

$$\Theta \equiv \mathrm{I}\frac{\mathrm{tr}(m)}{2},\Psi \equiv m-\Theta ,\Sigma \equiv \frac{\Psi +{\Psi }^{\mathrm{T}}}{2},\Omega \equiv \frac{m-m^{\mathrm{T}}}{2}=\frac{\Psi -{\Psi }^{\mathrm{T}}}{2},$$

$$r_i(\tau +\delta \tau )=(\mathrm{I}+\delta \tau m{)}_{ij}{r}_j(\tau ).$$

$$\delta Z(z)$$

$$\delta C(c)$$

$$\mathrm{I}+\Psi \delta \tau \approx (\mathrm{I}+\Sigma \delta \tau )(\mathrm{I}+\Omega \delta \tau )=\delta Z\delta R\approx \delta R\delta Z.$$

$${\tan }^{-1}\left(\frac{e_y^{\pm }}{e_x^{\pm }}\right)={\tan }^{-1}\left(\frac{r_{\mathrm{sm},y}}{r_{\mathrm{sm},x}}\right)\pm i{\tanh }^{-1}\left(\frac{S_i}{S_a}\right),$$

$${\omega }_z={\nabla }_x{v}_y-{\nabla }_y{v}_x=m_{yx}-m_{xy},$$

$${\mathrm{N}}_{\mathrm{t}}=\frac{1}{2\pi }\left(\frac{{\omega }_z}{2}\Delta \tau \right).$$

$${\mathrm{N}}_{\mathrm{t}}=\frac{1}{2\pi }{\int }_{L}d\tau \frac{\alpha (\tau )}{2}.$$

$$Z(z_f)$$

$$C(c_f)$$

$$\mathrm{I}+\Psi \delta \tau =RA\delta B{A}^{-1}{R}^{-1},$$

$${\tilde{r}}_i(\tau +n\delta \tau )$$

$$\mathrm{I}+\Psi \delta \tau =R(\theta )Z\left(\sqrt{\frac{S_a}{S_i}}\right)\delta R({\omega }_c)Z^{-1}\left(\sqrt{\frac{S_a}{S_i}}\right)R^{-1}(\theta ).$$

$${\lambda }_{+}-{\lambda }_{-}=2i{\omega }_c,$$

$${\mathrm{N}}_{\mathrm{c}}\equiv \frac{1}{2\pi }{\omega }_c\Delta \tau =\frac{1}{4\pi }\Im \{{\lambda }_{+}-{\lambda }_{-}\}\Delta \tau ,$$

$${\omega }_z=m_{yx}-m_{xy}={\Psi }_{yx}-{\Psi }_{xy}=\left(\frac{S_i}{S_a}+\frac{S_a}{S_i}\right){\omega }_c,\text{mbox}implying:|{\omega }_z|\ge 2|{\omega }_c|,$$

$$\alpha =\left(\frac{S_i}{S_a}+\frac{S_a}{S_i}\right){\omega }_c,$$

$${\mathrm{N}}_{\mathrm{c}}\equiv \frac{1}{2\pi }{\int }_{L}d\tau {\omega }_c(\tau ).$$

$$\displaystyle\mathrm{I}+\Psi\delta\tau=R(\theta)\,S\big{(}\cot(\vartheta)\big{)}\,\delta Z\left(\frac{\rho_{\mathcal{Z}}}{2}\right)\,S^{-1}\big{(}\cot(\vartheta)\big{)}\,R^{-1}(\theta)\ .$$

$$\frac{d}{d\tau }\ln \left(\frac{r_{+}(\tau )}{r_{-}(\tau )}\right)={\lambda }_{+}-{\lambda }_{-}={\rho }_{\mathcal{Z}}.$$

$${\omega }_z=m_{yx}-m_{xy}={\Psi }_{yx}-{\Psi }_{xy}={\rho }_{\mathcal{Z}}\cot (\vartheta ).$$

$$r_i(\tau +\Delta \tau )=J_{ij}{r}_i(\tau )=(\mathrm{I}+\delta \tau m{)}_{ij}^{\Delta \tau /\delta \tau }{r}_i(\tau ),$$

$$J=\exp \left(\frac{\mathrm{tr}(m)}{2}\Delta \tau \right)\times \left[\cosh \left(\frac{{\rho }_{\mathcal{Z}}}{2}\Delta \tau \right)\mathrm{I}+\Psi \frac{2\sinh \left(\frac{{\rho }_{\mathcal{Z}}}{2}\Delta \tau \right)}{{\rho }_{\mathcal{Z}}}\right],$$

$$\displaystyle Q=2\left[\csc^{2}(\vartheta)\,\cosh\big{(}\rho_{\mathcal{Z}}\Delta\tau\big{)}-\cot^{2}(\vartheta)\right]\ ,$$

$$\displaystyle\frac{d\ln Q}{d\tau}\sim\frac{d\ln\mathcal{Z}}{d\tau}=\rho_{\mathcal{Z}}=\lambda_{+}-\lambda_{-}\quad\mbox{for real eigenvalues of $m$},$$

$$\left|\frac{\Delta \ln Q}{\Delta \tau }\right|$$

$$\displaystyle\mbox{Option 1:}\quad\frac{d\ln\mathcal{Z}_{1}}{d\tau}\equiv\rho_{\mathcal{Z},1}=0\quad\mbox{for complex eigenvalues of $m$}$$

$$\displaystyle\mbox{Option 2:}\quad\frac{d\ln\mathcal{Z}_{2}}{d\tau}\equiv\rho_{\mathcal{Z},2}=\frac{2|\omega_{c}|}{\pi}\ln\left(\frac{\alpha^{2}}{2\omega_{c}^{2}}-1\right)\quad\mbox{for complex eigenvalues of $m$.}$$

Full text

Sun:magnetic fields — Sun: magnetic topology

1. 1 Introduction

2. 2 A toy model: linear steady-state planar flow

3. 2.1 Flow decomposition

4. 2.2 Coiling and squeezing rates, and the geometric interpretation of vorticity and the force-free parameter

5. 2.2.1 Elliptical flows as Squeeze – Infinitesimal rotation – Un-squeeze maps

6. 2.2.2 Saddle flows as Shear – Infinitesimal squeeze – Un-shear maps

7. 2.2.3 Squeezing in elliptical flows

8. 3 Field-line deviations

9. 3.1 Evolution of the field-line deviation vector in 3D

10. 3.2 Evolution of the field-line deviation vector in the normal plane

11. 3.3 Specifying the bent basis uniquely

12. 3.4 Invariants

13. 4 The squeeze factor

14. 4.1 Review of the standard squashing factor calculation

15. 4.2 The derivative of the squashing factor with respect to the footpoint locations

16. 4.3 Constructing an approximate local squashing factor rate

17. 4.4 The squeeze factor

18. 4.5 The symmetrized squeeze rate

19. 5 Coiling number

20. 5.1 The generalized force-free parameter

21. 5.2 The coiling rate

22. 5.3 The coiling number

23. 6 Axially-symmetric, force-free flux ropes

24. 6.1 Generic, axially-symmetric, force-free flux ropes

25. 6.2 The Gold-Hoyle and Lundquist flux ropes

26. 6.3 Further flux-rope examples

27. 6.3.1 Flux rope with constant coiling rate

28. 6.3.2 Flux rope with constant coiling number per unit flux-rope length

29. 6.4 Discussion

30. 7 Summary

31. A Results for curvilinear coordinates

32. A.1 The field-line deviation equation

33. A.2 Calculating the standard squashing factor

34. A.3 The magnetic field gradient

1 Introduction

The structure of the coronal magnetic field and its restructuring determines the energetics of solar flares and coronal mass ejections. Partitioning the magnetic field into flux domains is often the first step in studying the structure of magnetic fields. Such partitioning can be accomplished by reconstructing the magnetic skeleton, which is formed by features such as null points, separatrix surfaces and separators. That type of partitioning has been studied extensively in the past (e.g. Gorbachev & Somov, 1988; Longcope, 1996; Parnell et al., 2010). A broader type of partitioning of magnetic fields has been relatively recently employed using sheets of large, yet continuous, variation of the magnetic field, called Quasi-Separatrix Layers (QSLs) (e.g. Longcope & Strauss, 1994; Priest & Démoulin, 1995). QSLs have been used successfully in studying solar eruptions (e.g. Savcheva et al., 2012a; Janvier et al., 2013, 2014; Liu et al., 2014; Savcheva et al., 2015) as well as laboratory magnetoplasma (e.g. Lawrence & Gekelman, 2009), and have been associated with locations with exponentially large values of the dimensionless squashing factor value, $Q$ (e.g. Titov et al., 2002; Titov, 2007), where large current build-up can develop (e.g. Longcope & Strauss, 1994).

Two magnetic features that are of significant theoretical and observational interest are hyperbolic flux tubes (HFTs; Titov, 2007; Savcheva et al., 2012a, b) – formed at the intersection of QSLs – and flux ropes. The latter consist of magnetic field lines wrapping around a common flux rope axis. Flux ropes play a fundamental role in solar eruptions (e.g. Liu et al., 2016) by exhibiting instabilities, such as torus or kink instability (e.g. Aulanier et al., 2010; Török et al., 2004). Meanwhile, HFTs can generically form under flux ropes when the QSL wrapping around a flux rope, separating it from the overlying magnetic arcade, intersects with itself. Observable features associated with HFTs are flare ribbons whose locations and motions are reproduced by the photospheric traces of the corresponding HFT (Liu et al., 2014; Savcheva et al., 2015).

Without prior knowledge, locating magnetic features of interest – even as basic as flux ropes – in realistic magnetic field models is a non-trivial task since the domain decomposition obtained using QSLs can be extremely complex. As not all QSLs are associated with producing current sheets (indeed, QSLs are exhibited even in potential fields, carrying no current; see e.g. Titov (2007)), quite often one has no other option but to resort to picking by hand the QSLs of interest through ad-hoc thresholding in current density or magnetic field (e.g. Savcheva et al., 2015). Non-QSL-based methods for detecting flux ropes have been proposed (e.g. Yeates & Mackay, 2009; Yeates & Hornig, 2016; Lowder & Yeates, 2017), yet again they rely on ad-hoc thresholds: in the case of the method described by Lowder & Yeates (2017), that threshold is even time-dependent.

An important property of flux ropes used in their stability analyses is the twist number (e.g. Berger & Field, 1984), measuring the number of twists a field line in the flux rope makes around the flux rope core. As the twist number is a non-local quantity, requiring a proper choice of flux rope axis (e.g. Guo et al., 2017), one usually resorts to computing the local twist number111We avoid using the symbol $\mathcal{T}_{w}$ to denote the local twist number since Berger & Prior (2006) use that symbol to denote both the local and non-local twist (see their equations (12) and (16)). Meanwhile, Liu et al. (2016) use $\mathcal{T}_{w}$ to denote only the local approximation to the twist, calling the non-local twist $\mathcal{T}_{g}$; see their equations (7) and (13). Since we focus on the local deviations of field-lines, in most of this paper, we are concerned only with the local twist number and how it compares with the coiling number we introduce. However, in Section 6, we do compare those local quantities with the non-local twist number for simple flux-rope configurations., $\mathrm{N_{t}}$, instead, which is proportional to the component of the current parallel to the magnetic field (see Berger & Prior, 2006, eq. 16). For non-diverging current distributions, $\mathrm{N_{t}}$ approaches the exact twist number at the core of a flux rope (e.g. Liu et al., 2016). As twist is a dimensionless quantity, it is tempting to use the natural inequality, $\mathrm{N_{t}}\gtrsim 1$, as a robust flux-rope detection threshold. That threshold would seem to imply that nearby field lines exhibit at least $\mathcal{O}(1)$ twists around one another to qualify as being part of a flux rope. However, in this paper we analytically show that the twist number gives a biased estimate of the number of turns infinitesimally separated field lines make around one another. Indeed, in a follow-up paper (Savcheva, 2019), we apply the above threshold to realistic magnetic fields, and we find that along with flux ropes, it picks out many sheared, untwisted structures, which further corroborates that $\mathrm{N_{t}}$ gives a biased estimate of the local field-line twisting.

In this paper, we construct an unbiased local equivalent of the twist number, which we call the coiling number, $\mathrm{N_{c}}$. In (Savcheva, 2019), we show that the natural inequality $\mathrm{N_{c}}\gtrsim 1$ can be used successfully as a detection threshold for flux ropes. For each field line, $\mathrm{N_{c}}$ is given as an integral of a local coiling rate ($\omega_{c}$) over the field line length.

We also introduce a quantity, complementary to the coiling rate, which we dub the squeezing rate ($\rho_{\mathcal{Z}}$), measuring the local logarithmic rate of squeeze of neighboring field lines. The field-line integral over $\rho_{\mathcal{Z}}$ gives $\ln(\mathcal{Z})$, with $\mathcal{Z}$ defined as the squeeze factor. We show that $\mathcal{Z}\sim Q$ in a sense made precise in the paper; and therefore, $\rho_{\mathcal{Z}}$ gives an estimate of the local contributions to $Q$ – with the latter having no exact corresponding local squashing rate. The squeezing rate is large in the vicinity of HFTs and null points, and can be used for their detection.

To construct the quantities described above at a particular location, we pick the field line passing through that location – we call that field line the reference field line below – and explore the kinematics of the displacements of nearby magnetic field lines relative to that reference field line. We focus on the transverse part of those deviations, calculated in a plane normal to the reference field line (see Fig. 1). As the normal plane is moved along the length of the reference field line, the intersection points between that plane and the nearby field lines can be described by a planar “flow”. By construction, that flow exhibits a critical point at the location where the reference field line intersects the normal plane. As we focus on the kinematics of nearby field lines, that transverse flow is approximately linear. Therefore, the types of transverse field-line flow are given by the standard equilibrium solutions of linear systems, such as a saddle, cycle, spiral or node (see Fig. 2). The type of solution is unique for the normal plane that is only minimally rotated around its normal when moved along the reference field line (see Section 3.3). That solution is determined solely by the two eigenvalues of the gradient of the velocity of the nearby field lines in the normal plane. Here, velocity is defined as the time derivative of the transverse deviations, where “time” corresponds to the field-line length parameter along the reference field line. Even at this point, one can conjecture that critical points corresponding to saddles and nodes exhibiting large eigenvalue difference, can be used in locating HFTs and the vicinity of null points; whereas spirals and centers are the typical flow patterns in flux ropes. Thus, we construct our $\rho_{\mathcal{Z}}$ and $\omega_{c}$ using the properties of the transverse field-line flow.

Studying the transverse flow of field lines has the advantage of simplicity as one can apply intuition from the study of planar fluid flows. We use that to our advantage as in Section 2 we introduce most of our results in the context of a toy model of a steady-state linear planar fluid flow. Only then (Section 3) do we address issues arising from the fact that the field-line flow is “time”-dependent (i.e. the flow changes as one moves along the reference field line); that the normal plane is constantly changing orientation due to the curvature of the reference field line (see Fig 1); and that the type of flow in the normal plane depends on overall rotations of the basis spanning the normal plane around the reference field line. We specify that basis uniquely by the requirement that the rate of rotation of individual neighboring field lines around the reference field line is independent of whether that rate is calculated using the standard non-local twist rate (cf. eq. (82)), or in the basis spanning the normal plane. In the end, we demonstrate that for that unique normal-plane basis, the transverse flow type is entirely determined by the non-zero eigenvalues of the gradient of the normalized magnetic field, \mbox{\boldmath{\nabla}}{\hat{\mbox{\boldmath{B}}}}, which match those of the gradient of the transverse velocity of the flow. The real part of the difference of the non-zero eigenvalues of \mbox{\boldmath{\nabla}}{\hat{\mbox{\boldmath{B}}}} gives $\rho_{\mathcal{Z}}$ (Section 4), while the imaginary part equals $2\omega_{c}$ (Section 5). In Section 6 we write our results for generic, axially-symmetric, force-free flux ropes, and show numerical solutions for example flux rope configurations. In Section 7 we present a summary of our result. In Appendix A we write our results in curvilinear coordinates.

2 A toy model: linear steady-state planar flow

In this section we write our coiling and squeezing rates for a toy model of a linear steady-state 2-dimensional (2D) flow. We use this model to show the relationship between the quantities we introduce in this paper and standard fluid notions, such as vorticity and shear. Along the way, we relate the flow vorticity to the coiling and squeezing rates using properties of the global geometric deformation of the flow.

2.1 Flow decomposition

Let us consider the steady-state 2D linear flow with a velocity field ($v$) given below as a function of position, $r$:

[TABLE]

The constant matrix $m$ above equals the velocity gradient, $m_{ij}=\nabla_{j}v_{i}$. Clearly {\mbox{\boldmath{r}}}=0 is a critical point, which depending on the eigenvalues $\lambda_{\pm}$ of $m$, can be a right-handed (RH) or left-handed (LH) stable/unstable spiral (two complex conjugate eigenvalues), a node (real eigenvalues of the same sign), a saddle (real eigenvalues of opposite signs); or in certain special cases: a center (purely imaginary eigenvalues) and an (im)proper node (see Fig. 2 for example illustrations). Even if a flow does not exhibit a critical point, the equation above can always be written for the (linearized) relative displacement vector ($r$) between two neighboring streamlines as a function of time ($\tau$). Then it is easy to see how eq. (1) can be later reinterpreted as giving the transverse deviation rate between two neighboring magnetic field lines (separated by a transverse displacement $r$) as a function of the field line length parameter $\tau$. In that case, we show that one must replace $m$ with the transverse gradient of the normalized magnetic field, {\hat{\mbox{\boldmath{B}}}} (cf. eq. (49)).

The velocity gradient ($m$) is usually decomposed as follows (e.g. Thorne & Blandford, 2017):

[TABLE]

with:

[TABLE]

where $\mathrm{I}$ is the identity matrix. Clearly, $\Theta$ is the diagonal, trace part of $m$, usually called the rate-of-expansion tensor. Then, $\Psi$ is the traceless part of $m$. The latter is in turn usually decomposed into a symmetric traceless part (the rate-of-shear tensor; $\Sigma$) and an antisymmetric part (the rate-of-rotation tensor; $\Omega$). This split is manifestly possible for all types of critical points. The rate of expansion corresponds to uniform contraction/expansion of the area (in 3D, that would be the volume) of a fluid element. Once that uniform expansion flow is subtracted from $m$, that renders the critical point of the remaining traceless part of the velocity gradient ($\Psi$) into either a saddle (for real eigenvalues of $m$) or a center (for complex eigenvalues of $m$).

The standard decomposition of $m$ above, however, is not unique. To see that, note that for infinitesimal $\delta\tau$, eq. (1) can be written in finite difference form as:

[TABLE]

The infinitesimal transformation $(\mathrm{I}+\delta\tau m)$ entering above, can be an infinitesimal squeeze ($\delta Z$), shear ($\delta S$), scaling ($\delta C$) or rotation ($\delta R$) map, or a composition of those. With an appropriate choice of basis, those infinitesimal maps can be written as (for the corresponding finite transformations, see eq. (2.2) below):

[TABLE]

Those four transformations depend on four parameters: $z$, $s$, $c$ and $\varpi$, which are defined through the equations above.

Clearly, the rate-of-expansion tensor represents the rate of scale change, and thus can be written in terms of $\delta C$ as $(\mathrm{I}+\delta\tau\Theta)=\delta C$, which is proportional to the identity matrix. Since it commutes with the rest of the maps, we need to focus only on the traceless part, $\Psi$, of the velocity gradient. Note that $\Psi$ depends on three degrees of freedom since $m$ has four elements, but we imposed the condition $\mathrm{tr}(\Psi)=0$. One of those degrees of freedom amounts to a trivial choice of orthonormal basis, which is defined up to an overall rotation. Thus, $\Psi$ has only two non-trivial degrees of freedom, which means that it can be decomposed in various ways using the infinitesimal maps $\delta Z$, $\delta S$ and $\delta R$, given in eq. (2.1). Let us explore those decompositions below.

Since the rate of shear, $\Sigma$, is symmetric, it has two orthogonal eigenvectors, and because it is traceless, its eigenvalues are opposite in sign. Thus, $(\mathrm{I}+\Sigma\delta\tau)$ can be treated as an infinitesimal squeeze map, $\delta Z$, along the eigenvectors of $\Sigma$. As $\Omega$ is antisymmetric, $(\mathrm{I}+\Omega\delta\tau)$ is an infinitesimal rotation map, $\delta R$, independent of the choice of orthonormal basis vectors. Thus, (up to an overall rotation aligning one of the basis vectors with one the principle axes of $\delta Z$) the decomposition $\Psi=\Sigma+\Omega$, can be written as:

[TABLE]

As eq. (4) is linear in $\delta\tau$, we kept only terms to linear order in $\delta\tau$ in the above equation. Thus, the velocity fields generated by each of the infinitesimal maps above can be simply added (as in linear superposition) to obtain the velocity field of the composite infinitesimal map $\delta R\delta Z$. This can be seen by using $\Psi$ from eq. (6) in combining equations (1) and (2). Thus, in the deformation map language, the standard decomposition of $\Psi$ into rate-of-shear and rate-of-rotation tensors corresponds to the superposition of an infinitesimal squeeze flow and an infinitesimal rotation flow. An illustration of that superposition is shown in the first two rows of Fig. 3 for the two possible flows generated by $\Psi$: a saddle and a center.

As noted earlier, one can use other (non-standard) decompositions of $\Psi$, however. For saddles, one can write $(\mathrm{I}+\Psi\delta\tau)$ using the Schur decomposition as a rotated upper triangular matrix, which can in turn be decomposed 222In the context of fluid flows, see (Keylock, 2017) for an example usage of the Schur decomposition; and see (Kolar, 2007) for an example of using, what they call, shear-plus-residual decomposition. as $\delta Z\delta S$. That decomposition can be achieved by choosing the (repeated) eigenvectors of $\delta S$ to coincide with one of the eigenvectors of $\delta Z$ as well as with one of the eigenvectors of $\Psi$. That shear-squeeze decomposition is illustrated in the third row of Fig. 3. For centers, one can decompose $(\mathrm{I}+\Psi\delta\tau)$ as $\delta R\delta S$, with the shear aligned with one of the semi-axes of the concentric elliptical trajectories around the center333A closed-form expression which allows one to determine the direction {\mbox{\boldmath{r}}}_{\mathrm{sm}} of the semi-major axis of elliptical flow around a center generated by a matrix $\Psi$ is not readily available in the literature. Thus, we include that here for completeness. First, one needs to find the complex conjugate eigenvectors \mbox{\boldmath{e}}^{\pm} of $\Psi$. Then, after a bit of algebra, one can show that the following relationship holds between the $x$ and $y$ components of {\mbox{\boldmath{r}}}_{\mathrm{sm}} and \mbox{\boldmath{e}}^{\pm}:

$\tan^{-1}\left(\frac{e^{\pm}_{y}}{e^{\pm}_{x}}\right)=\tan^{-1}\left(\frac{r_{\mathrm{sm},y}}{r_{\mathrm{sm},x}}\right)\pm i\tanh^{-1}\left(\frac{S_{i}}{S_{a}}\right),$

where $i=\sqrt{-1}$, and $S_{i}/S_{a}$ is the (fixed) aspect ratio (the ratio of the semi-minor ($S_{i}$) to semi-major ($S_{a}$) axis) of one of the concentric ellipses around the center. Taking the real and imaginary parts of that expression allows one to find the direction of {\mbox{\boldmath{r}}}_{\mathrm{sm}}, as well as the aspect ratio of the elliptical flow.. That non-standard decomposition is illustrated in the fourth row of Fig. 3.

In the rotation-squeeze decomposition of a center flow, one can show that the rotation rate is different than the one obtained in the rotation-shear decomposition; and furthermore, the latter is not unique as it depends on whether the shear is aligned with the semi-minor or semi-major axis of the elliptical flow. Similarly, the squeeze rate of a saddle flow depends on whether the flow is decomposed into squeeze plus shear or into squeeze plus rotation flows (for the latter, see Section 4.5). Thus, defining local squeeze and rotation rates is prone to pitfalls if one is not careful in stating precisely what one is trying to quantify. Below we include an example which is a clear illustration of that.

One standard quantifier of flow rotation is vorticity, \mbox{\boldmath{\omega}}\equiv\mbox{\boldmath{\nabla}}\times{\mbox{\boldmath{v}}}, where $v$ is the velocity field of the flow. A planar flow (in a plane spanned by the basis \hat{\mbox{\boldmath{x}}}, \hat{\mbox{\boldmath{y}}}) has a vorticity vector which is normal (along the third basis vector \hat{\mbox{\boldmath{z}}}) to the flow and is given by

[TABLE]

which has a magnitude equal to $\sqrt{-2\mathrm{tr}(\Omega^{2})}$ and is, therefore, a rotation invariant. Thus, the vorticity equals twice the rotational rate ($\varpi$; see eq. (2.1)) one infers from the rotation-squeeze decomposition of the velocity gradient tensor. As the squeeze map is symmetric, it does not produce vorticity for that decomposition. Meanwhile, in the rotation-shear decomposition of a center flow, both the infinitesimal shear and rotation produce vorticity, and therefore, the rotational rate associated with the infinitesimal rotation map in that decomposition is not as straightforward to relate to vorticity.

A standard intuitive way (e.g. Thorne & Blandford, 2017) to understand vorticity is to imagine placing a vane with orthogonal fins in the flow. Such a vane will spin at an angular velocity given by half the vorticity as can be seen by averaging the angular velocity obtained from eq. (1) for any two orthogonal vectors $r$. Such a vane would also rotate in sheared flows, as well as in saddle flows with non-orthogonal eigenvectors. Yet, in the case of shear or saddle flow, fluid elements do not revolve around {\mbox{\boldmath{r}}}=0; and thus, the number of turns they make around {\mbox{\boldmath{r}}}=0 is always less than one, despite the fact that the imaginary vane makes $\mathrm{N_{t}}$ turns for a time interval $\Delta\tau$:

[TABLE]

Moreover, even for elliptical flows, below we show that the rate at which a vane would rotate fails to match the average rate at which fluid elements revolve around {\mbox{\boldmath{r}}}=0 (see Section 2.2.1). Thus, clearly vorticity is not the quantity we should be using to quantify the number of turns a fluid element makes around {\mbox{\boldmath{r}}}=0 over $\Delta\tau$.

At this point, the reader should compare the above equation with the equation for the standard local magnetic twist number (Berger & Prior (2006); see also our eq. (97)), which we reproduce here for convenience:

[TABLE]

Here, the integral is over a field line $L$; while $\alpha$ is the generalization of the force-free parameter (which we define for any magnetic field $B$ as \alpha\equiv(\mbox{\boldmath{\nabla}}\times\mbox{\boldmath{B}})\cdot\mbox{\boldmath{B}}/B^{2}). For magnetic fields, $\alpha$ plays the role of the fluid vorticity above (compare eq. (7) with eq. (84)). In Section 5, we show that $\alpha$ equals twice the average of the angular rate of motion of neighboring field lines (whether that angular motion is due to rotation or shear) around a reference field line (cf. eq. (86)). The average is taken over neighboring field lines, assuming they are uniformly distributed in angle, similar to the way the orthogonal fins of the vane, floating in the fluid flow described above, sample streamlines equally spaced in angle. Thus, similar to the fluids case, the standard local twist number fails to give the correct number of turns magnetic field lines make around a reference field line when they exhibit saddle-like transverse flow. $\mathrm{N_{t}}$ fails to return the correct number of turns even for elliptical transverse field-line flows (see Sections 2.2.1 and 5.2). We address that issue in the next section.

2.2 Coiling and squeezing rates, and the geometric interpretation of vorticity and the force-free parameter

In the previous section, we wrote the vorticity in the standard way, in terms of the rotation rate tensor, which we related to an infinitesimal rotation map. In this section we reinterpret the vorticity of a 2D flow in terms of finite (as opposed to infinitesimal) deformation maps of the flow. Along the way, we will be able to construct local measures of the rotation and squeeze rates of neighboring field lines.

Let us start by writing the finite squeeze ($Z$), shear ($S$), scaling ($C$) and rotation ($R$) maps, obtained by composing (infinitely many of) their corresponding infinitesimal maps (eq. (2.1)):

[TABLE]

with the finite map factors $z_{f}$, $s_{f}$, $c_{f}$, $\theta$ defined through the equations above. Any rate of expansion ($\Theta$) is uniquely captured by a scale map ($C$) which, being proportional to the identity matrix, commutes with all other maps. Therefore, similar to the previous section, here we focus solely on the traceless part, $\Psi$, of the velocity gradient.

In this section, we decompose $\Psi$ into a set of finite transformations (given by eq. (2.2)) and one infinitesimal transformation (eq. (2.1)). The role of the finite maps is to undo certain deformations to the flow, so that we can use only one infinitesimal map to describe its evolution. When combining infinitesimal maps, we showed that the velocity field produced by them can be treated as a linear superposition of the velocity field produced by each individual infinitesimal map (see the discussion around eq. (6)). However, under the decomposition we are about to write down, not all maps are infinitesimal, and therefore one can no longer use linear superposition to study the map combinations. Yet, the interpretation of the results in the end will be as intuitive.

Let us focus on transformations of the type:

[TABLE]

with $A$ and $\delta B$ being a finite and an infinitesimal transformation map, respectively. The finite rotation maps above are used to rotate the flow to the principle axes of the transformation $A$. The rotation maps do not affect the invariants of $\Psi$, and we will therefore omit them from most of the discussion below. The decomposition shown in the equation above has three parameters (one for each of the maps), matching the degrees of freedom of $\Psi$.

Note that the structure of the decomposition in eq. (11) is chosen to be such that the flow, \tilde{\mbox{\boldmath{r}}}(\tau), generated by the traceless part of (the constant) $m$ in eq. (4) can be solved to give the position of a fluid element after a finite interval $\Delta\tau=n\delta\tau$ (with $n\to\infty$) as:

[TABLE]

In the last equality we used the fact that the composition of infinitesimal maps $\delta B$ results in a finite map $B$. Thus, the flow produced by $\Psi$ after some finite $\Delta\tau$ matches the expression for the infinitesimal flow (eq. (11)) with the infinitesimal $\delta B$ replaced by a finite $B$. Thus, one can think of this decomposition as a global deformation of the flow done by $A$, with the $\tau$-evolution captured by $B$ (see e.g. Fig 4, discussed below). This will allow us to interpret our results more easily below.

Transformation $A$ can be either a shear or a squeeze (as any rotation can be absorbed into $R$). Therefore, the non-trivial combinations of the pairs $(A,\delta B)$ are: (squeeze, shear), (shear, squeeze), (shear, rotation), (squeeze, rotation), where we eliminated repeated pairs as the composition of $A\,\delta A$ is still a transformation of the same type, $A$, which is not sufficient to parametrize generic elliptic and saddle flows.

Of the pairs listed above, one can check explicitly that the (squeeze, shear) pair gives a composition $Z\,\delta S\,Z^{-1}$, which is again an infinitesimal shear map, and therefore, cannot describe generic elliptic and saddle flows. The pair (shear, squeeze) has two real eigenvalues and therefore can describe saddles. The (shear, rotation) and (squeeze, rotation) pairs have two complex conjugate eigenvalues and therefore can describe centers. We find the latter pair more intuitive when discussing elliptical flows, and therefore we focus on it below. We return to saddle flows in Section 2.2.2.

2.2.1 Elliptical flows as Squeeze – Infinitesimal rotation – Un-squeeze maps

Let us first focus on describing centers using the pair ($A,\delta B$)=(squeeze, rotation) in eq. (11). The rotation matrices in that equation can be used to rotate the basis vectors by an angle $\theta$ so that they are aligned with the semi-axes of the elliptical flow around the center generated by $\Psi$. Then $Z^{-1}(z_{f})$ can be used to un-squeeze the elliptical flow and render it into a uniform circular flow, rotating with a constant angular frequency $\omega_{c}$, which we will call the coiling rate. At each time-step ($\delta\tau$), that rotation is generated by $\delta R(\omega_{c})$. After that, we need to squeeze the circular flow back into an elliptical flow. For an elliptical flow described by concentric ellipses with a fixed aspect ratio $S_{i}/S_{a}$ (where $S_{i}$ and $S_{a}$ are the semi-minor and semi-major axes of one of those ellipses, respectively), the squeeze factor must equal $z_{f}=\sqrt{S_{a}/S_{i}}$ as $Z(z_{f})$ (eq. (2.2)) rescales the component $r_{x}$ of $r$ by $z_{f}$ and $r_{y}$ by $1/z_{f}$. Thus, we can write:

[TABLE]

That decomposition is illustrated in the second row of Fig. 4.

Using equations (2.1) and (2.2), we can write out the right-hand-side of eq. (13) explicitly as a matrix and then calculate its invariants. We find that the eigenvalues444In the context of fluid flows, a vortex detection criterion using the imaginary part of the complex eigenvalues (referred to as the swirling strength of a vortex) of the velocity gradient was used by Zhou et al. (1999). Even though that criterion can be directly ported to magnetism by using the gradient of the magnetic field, our method differs in several significant ways. (1) As we show in Section 5, we use the gradient of the normalized magnetic vector field to calculate $\omega_{c}$ (which is why we choose to call it by a different name: the coiling rate). For fluids, normalizing the velocity field would break Galilean invariance for the vortex detection method; however, in the case of magnetic fields in the solar corona, we do have a preferred reference frame. In Section 3, we show that the non-zero eigenvalues of the gradient of the normalized magnetic field in 3D match those of the 2D gradient of the rate-of-deviation of neighboring field lines in the plane normal to a reference field line. This allows for a simple interpretation of our results. (2) Zhou et al. (1999) use a vortex detection criterion by choosing an ad hoc threshold for the swirling strength. Our flux rope detection criterion is a threshold on the integral of $\omega_{c}$ over the field line length, which in (Savcheva, 2019) we set at $\mathrm{N_{c}}\gtrsim 1$ (see eq. 18). In other words, our criterion is that neighboring field lines should make at least one winding around one another to classify as part of a flux rope. of $\Psi$ in this decomposition are $(\pm i\omega_{c})$. As $\Theta$ does not affect the difference of the eigenvalues of $m$, we can write that the difference between the eigenvalues ($\lambda_{\pm}$) of $m$ equals that of $\Psi$:

[TABLE]

and therefore, the period with which fluid streamlines revolve around {\mbox{\boldmath{r}}}=0 equals $2\pi/\omega_{c}$. This is not surprising in light of eq. (12), from which one can see that for the (squeeze, rotation) decomposition, the periodicity of the flow generated by $\Psi$ matches that of $B$, which in our case is the composition of infinitesimal $\delta R(\omega_{c})$ maps. Thus, we can define a coiling number for elliptical flows as:

[TABLE]

with the coiling rate, $\omega_{c}$, defined as the rotational rate of the flow in the decomposition of eq. (13) (bottom row of Fig. 4).

It is important to note that the same result above (eq. (14)) is obtained in the (shear, rotation) decomposition, highlighting the robustness of $\omega_{c}$. However, we find the result for vorticity obtained in the (squeeze, rotation) decomposition to be more intuitive (see below), which is why we focus on that decomposition.

Since the flow is elliptical, the angular velocity of the streamlines is not constant. Therefore, the exact number of turns a particular fluid element makes around {\mbox{\boldmath{r}}}=0 in a finite $\Delta\tau$ depends on the initial position of the fluid element as well as on $\Delta\tau$ (cf. Section 5, as well as (Berger & Prior, 2006)). However, since the flow is periodic, when $\mathrm{N_{c}}$ is an integer, all fluid elements would have made exactly $\mathrm{N_{c}}$ turns around the origin, independent of their initial positions. Another way of thinking about $\mathrm{N_{c}}$ is that the un-squeezed flow (middle step in the decomposition shown in the second row of Fig. 4) is uniformly rotating with $\omega_{c}$; and therefore, we can treat its rotation rate as an average of the actual fluid flow rotation rate in the sense that the periodicities of the two flows match.

Similar to the above analysis that lead to eq. (14), we can use eq. (13) to calculate the vorticity of the flow (eq. (7)):

[TABLE]

with the sign of $\omega_{c}$ chosen to match that of $\omega_{z}$. Thus, for center flows, we can see that vorticity is proportional to the coiling rate as well as to a geometric factor, related to the aspect-ratio of the elliptical flow. Thus, the local twist number, eq. (8), is boosted relative to the coiling number, eq. (15), by terms dependent on that aspect ratio. From the discussion above, we can see that the coiling number, $\mathrm{N_{c}}$, is a much better quantifier than the local twist number for the average number of turns fluid elements make around the origin in a finite $\Delta\tau$. At the very least, as discussed above, when $\mathrm{N_{c}}$ is an integer, it matches the actual number of turns the fluid elements make, while $\mathrm{N_{t}}$ always overestimates that number for non-circular flows.

The above result is easy to understand intuitively. Going back to our submerged vane example, a vane floating at the origin would be predominantly affected by fluid elements at the co-vertex of their streamlines (the points of closest approach to the origin, lying on the tips of the minor axis). For ellipses with large aspect ratio, one can see that at the co-vertex, fluid elements sweep large angles (relative to the origin) over short amount of time even for fluid flows with large period. Thus, our imaginary vane would rotate at an angular speed which is boosted relative to the average angular speed of the fluid elements around the origin.

In the context of transverse magnetic field-line deviations, $\alpha$ is the magnetic equivalent of $\omega_{z}$ (see the discussion around eq. (9)). In Section 5, we find that we can write eq. (16) for $\alpha$ as (cf. eq. (95)):

[TABLE]

with $\omega_{c}$ given by eq. (14), where the eigenvalues correspond to to the non-zero eigenvalues of the 3D gradient of the normalized magnetic field (compare equations (14) and (93)). Thus, for magnetic fields we can define a coiling number as:

[TABLE]

Similar to the coiling number in the fluid context, $\mathrm{N_{c}}$ is a much more robust estimate of the number of turns neighboring magnetic field lines make around a reference field line, getting that number exactly right for constant $\omega_{c}$ when $\mathrm{N_{c}}$ is an integer. That is again in contrast with the biased magnetic local twist number, which is boosted by the aspect ratio of the transverse elliptical field-line flow.

2.2.2 Saddle flows as Shear – Infinitesimal squeeze – Un-shear maps

Now let us focus on representing saddle flows using the pair ($A,\delta B$)=(shear, squeeze) in eq. (11). Similar to the discussion of elliptical flows above, we can use the rotation matrices in that equation to rotate the basis vectors by an angle $\theta$ so that one of the basis vectors is aligned with one of the eigenvectors of $\Psi$. Then $S^{-1}(s_{f})$ can be used to un-shear the flow and render the second eigenvector of $\Psi$ orthogonal to the first. The flow would then be a saddle flow with orthogonal eigenvectors which, at each time-step ($\delta\tau$), can be generated by an infinitesimal squeeze map $\delta Z(\rho_{\mathcal{Z}}/2)$ with a squeeze rate $\rho_{\mathcal{Z}}$, where the factor of $1/2$ is introduced for convenience. After the infinitesimal squeeze, we need to shear back the flow with $S(s_{f})$ so that the directions of its asymptotes again match the eigenvectors of $\Psi$. If we denote the angle between the eigenvectors of $\Psi$ as $\vartheta$ (see top-left panel of Fig. 4), then the required shear factor needs to equal $s_{f}=\cot(\vartheta)$, which one can check explicitly using eq. (2.2). Thus, we can write this decomposition as:

[TABLE]

That decomposition is illustrated in the first row of Fig. 4.

Following the analysis in the previous section, using equations (2.1) and (2.2), we can write out the right-hand-side of eq. (19) explicitly as a matrix and then calculate its invariants. We find that the eigenvalues of $\Psi$ in this decomposition are $(\pm\rho_{\mathcal{Z}}/2)$. As $\Theta$ does not affect the difference in eigenvalues of $m$, we can conclude that the difference between the eigenvalues of $m$ equals that of $\Psi$. Therefore, the (relative logarithmic) squeeze rate between the directions corresponding to the eigenvectors of $\Psi$, and hence $m$ (as $\Theta$ does not affect the eigenvectors of $m$ either), is given by $\rho_{\mathcal{Z}}$. Indeed, if \hat{{\mbox{\boldmath{v}}}}_{\pm} are the two unit eigenvectors of $m$ corresponding to $\lambda_{\pm}$, then for two fluid parcels lying on the asymptotes of the flow (in the eigenvector directions relative to the origin), we can write {\mbox{\boldmath{r}}}_{+}(\tau)=r_{+}(\tau)\hat{{\mbox{\boldmath{v}}}}_{+} and {\mbox{\boldmath{r}}}_{-}(\tau)=r_{-}(\tau)\hat{{\mbox{\boldmath{v}}}}_{-}. Then, from eq. (1), one can check explicitly that $d\ln(r_{+}(\tau)/r_{-}(\tau))/d\tau=\lambda_{+}-\lambda_{-}$, and therefore, we can summarize:

[TABLE]

Using the matrix entering on the right-hand-side of eq. (19), we can express the vorticity of the flow (eq. (7)) as:

[TABLE]

Thus, we can see that similar to the analysis of elliptical flows in the previous section, the vorticity (and hence, the force-free parameter $\alpha$ for magnetic fields) for saddle flows has a clear dependence on the geometry of the flow. Through its dependence on the angle ($\vartheta$) between the flow eigenvectors, the vorticity quantifies the shear necessary to make the asymptotes of the flow orthogonal. For flows with orthogonal asymptotes (produced by an infinitesimal squeeze), we have $\vartheta=\pi/2$, which means that there is no overall shear in this decomposition, and no vorticity.

In the context of transverse magnetic field-line deviations, one can use the squashing factor ($Q$) (Titov, 2007) to quantify the divergence of neighboring field lines, separated by $r$ in the plane normal to one of them. Let us see how $\rho_{\mathcal{Z}}$ relates to $Q$. We review the calculation of $Q$ in Section 4.1. Here we highlight just the parts necessary for this section. The squashing factor is a function (see eq. (61)) of the transformation matrix $J$ relating $r$ evaluated at two different locations along a field line: at $\tau$ and $\tau+\Delta\tau$ for some finite $\Delta\tau$:

[TABLE]

where we used eq. (4). One can decompose $m$ into a vector times the Pauli matrices vector, then use the exponentiation of a Pauli vector formula to obtain in the limit of $\delta\tau\to 0$:

[TABLE]

where we used equations (3) and (20). This can in turn be used to calculate the squashing factor (using eq. (61)). After some algebra, we obtain:

[TABLE]

with $\vartheta$ given by eq. (21). The equation above is valid for any steady-state planar flow.

We are interested in large $Q\gg 1$, which holds555In the limit $\vartheta\to 0$, $m$ has repeated eigenvectors, and despite the eigenvalues being the same (implying $\rho_{\mathcal{Z}}\to 0$), $Q$ can still grow large, as can be seen from eq. (24). However, one can check that $Q$ in that case grows only as $\Delta\tau^{2}$, and not exponentially, and therefore we disregard those special cases. for large real arguments of $\cosh()$, in which case the logarithm of $Q$ approaches $\ln(Q)\sim\rho_{\mathcal{Z}}\Delta\tau$. Therefore, we can approximate $\ln(Q)$ as an integral over the local squeezing rate, given by $\rho_{\mathcal{Z}}$ in eq. (20). We denote that approximation of $Q$ with $\mathcal{Z}$, and we will refer to the latter as the squeeze factor. Therefore, the squeeze factor is the solution to the equation (compare with eq. (76) and (78)):

[TABLE]

with the squeezing rate, $\rho_{\mathcal{Z}}$, defined as the squeezing rate of the flow in the decomposition666Note that $\lambda_{+}-\lambda_{-}$ also equals the squeezing rate one would obtain from the infinitesimal squeeze plus infinitesimal shear decomposition of the flow (third row of Fig. 3). In that decomposition, the vorticity of the flow equals the negative of the shear rate entering in $\delta S$. In contrast, the decomposition used in this section allows one to interpret vorticity through geometric distortions of the flow (eq. (21)), bringing it in parallel with the discussion presented in Section 2.2.1 for elliptical flows. of eq. (19) (top row of Fig. 4).

From the above equation, we can see that we can treat the squeezing rate as a local approximation for the logarithmic rate of squashing for constant $m$. In Section 4.3 we generalize the above equation (using somewhat different arguments) to field lines in 3D with non-zero curvature and show that in that case, the eigenvalues entering above correspond to the two non-zero eigenvalues of the 3D gradient of the normalized magnetic field.

2.2.3 Squeezing in elliptical flows

Our result (eq. (25)) from the previous section applies to squeezing under saddle flows. In the context of magnetic fields, under elliptical transverse flows, a flux tube is periodically squashed and un-squashed; thus, on average, $Q$ remains the same. This can be seen from eq. (24), which is valid even for complex eigenvalues of $m$, in which case $\rho_{\mathcal{Z}}=2i\omega_{c}$ (compare equations (14) and (20)). Note that $Q$ can periodically still grow large, however. It is maximized at $\Delta\tau=\pi/(2\omega_{c})$, and every $\Delta\tau=\pi/\omega_{c}$ later along the field-line, in which case one can show that $Q$ equals $Q_{\mathrm{max.}}=S_{a}^{2}/S_{i}^{2}+S_{i}^{2}/S_{a}^{2}\geq 2$, where we used eq. (16) and (21) to eliminate $\vartheta$ from eq. (24). Comparing with eq. (67) which we write down in our review Section 4.1, a flux tube, with an initially circular cross-section, within a quarter of a period (by “period” we loosely refer to the field-line length corresponding to $2\pi/\omega_{c}$) is transversely squashed into an ellipse with an aspect ratio equal to the square of the aspect ratio of the integral lines of the transverse flow. In the next quarter of a period, the flux tube gets completely un-squashed back into a circle, and so on. Thus, one can define an average squashing rate (localized) within each quarter of a period: starting with an initial $Q=2$, within a quarter of a period, $\ln(Q)$ changes by:

[TABLE]

So, at this point we have two different possibilities ($\rho_{\mathcal{Z},1|2}$) for the logarithmic squeezing rate under elliptical flows. If one is not interested in periodic localized flux-tube squashing, one can use:

[TABLE]

since $Q$ periodically returns to 2 as discussed above. If localized periodic squashing is of interest, however, one can take that into account by defining the squeeze factor as:

[TABLE]

In the equation above, we used eq. (17) to express eq. (26) using the generalized force-free parameter. In Section 5, we express $\alpha$ and $\omega_{c}$ as scalars constructed using the gradient of the normalized magnetic field.

The ambiguity in the definition of $\mathcal{Z}$ above for elliptical flows is a consequence of the fact that the choice of squeezing rate is, in the end, application specific. In other words, our definition for $\mathcal{Z}$ offers fine-grained control over which effects causing squeezing are included and which – not. Note, however, that in defining the coiling number in eq. (15) we ignored any partial coiling caused by transverse saddle flows. So, under Option 1 above (eq. (27)), $\mathcal{Z}$ carries complementary information to that of $\mathrm{N_{c}}$: the former would depend only on the real part of the eigenvalue difference of $m$, while the latter depends only on the imaginary part777Indeed, an infinitesimal rotation map can be diagonalized to give an infinitesimal squeeze map with an imaginary squeeze factor; see eq. (2.1)..

3 Field-line deviations

In this section, we write an equation for the relative displacement of neighboring field-lines in the 2D subspace spanned by the plane normal to those field lines. Those transverse field-line deviation vectors are changing as one moves the normal plane down the field lines. Thus, the flow of the transverse deviations with the field-line length parameter ($\tau$) is planar, and takes the form of eq. (1), which was the starting point of our analysis in the previous section. However, in order to write that equation in terms of the magnetic field $B$, one has to address issues arising from the fact that generic magnetic field lines have curvature, and therefore the normal plane does not have a constant orientation in 3D. We address those complications in this section.

To simplify the analysis, we write our equations in a 3D Cartesian basis. We generalize the results to curvilinear coordinates in Appendix A.

3.1 Evolution of the field-line deviation vector in 3D

Let us pick two neighboring field lines with radius vectors {\mbox{\boldmath{x}}}(\tau) and {\mbox{\boldmath{y}}}(\tau) with $\tau$ being the field-line length parameter. We pick an initial condition such that $\tau=\tau_{0}$ when both field lines intersect the plane normal to the first field line at {\mbox{\boldmath{x}}}(\tau_{0}). The magnetic field lines {\mbox{\boldmath{x}}}(\tau) and {\mbox{\boldmath{y}}}(\tau) are calculated as the integral curves of the unit magnetic field vectors, {\hat{\mbox{\boldmath{B}}}}. Thus, in Cartesian coordinates, we have:

[TABLE]

and similarly for {\mbox{\boldmath{y}}}(\tau).

Then we can construct the field-line deviation vector \delta{\mbox{\boldmath{x}}}(\tau)\equiv{\mbox{\boldmath{y}}}(\tau)-{\mbox{\boldmath{x}}}(\tau). Using (29), the deviation vector \delta{\mbox{\boldmath{x}}}(\tau) is a solution to the linearized equation:

[TABLE]

Let us write this equation in component form, suppressing the $\tau$ dependence on the right hand side (we use the Einstein summation convention):

[TABLE]

where the last equality defines the (generally) non-symmetric matrix $M$. Index placement in this section is unimportant as we work in a Cartesian basis.

Let us additionally define the deviation vector {\mbox{\boldmath{r}}}(\tau)\equiv{\mbox{\boldmath{y}}}(\tau^{\prime}(\tau))-{\mbox{\boldmath{x}}}(\tau), with $\tau^{\prime}$ chosen such that {\mbox{\boldmath{r}}}(\tau) is orthogonal to the field line at {\mbox{\boldmath{x}}}(\tau):

[TABLE]

Above we used the fact that {\hat{\mbox{\boldmath{B}}}}({\mbox{\boldmath{x}}}) is the tangent vector for the field line passing through $x$. Choosing the field lines {\mbox{\boldmath{x}}}(\tau) and {\mbox{\boldmath{y}}}(\tau) infinitesimally close to each other, we can write $\tau^{\prime}(\tau)=\tau+\varepsilon(\tau)$ with $\varepsilon$ being small. Thus, we can linearize:

[TABLE]

where we used eq. (29). Setting \mbox{\boldmath{r}}(\tau)\cdot{\hat{\mbox{\boldmath{B}}}}({\mbox{\boldmath{x}}}(\tau))=0, we find \varepsilon(\tau)=-\delta{\mbox{\boldmath{x}}}(\tau)\cdot{\hat{\mbox{\boldmath{B}}}}({\mbox{\boldmath{x}}}(\tau)), which in turn implies888 The same result was obtained by Scott et al. (2017). (in matrix notation):

[TABLE]

where the projection operator is given in component form by: P_{\perp,ij}({\mbox{\boldmath{x}}}(\tau))=\mathrm{I}_{ij}-\hat{B}_{i}({\mbox{\boldmath{x}}}(\tau))\hat{B}_{j}({\mbox{\boldmath{x}}}(\tau)). From eq. (32) we can see that $r$ equals the projection of \delta{\mbox{\boldmath{x}}} onto the plane normal to the reference field line at location {\mbox{\boldmath{x}}}(\tau).

We can write an evolution equation for {\mbox{\boldmath{r}}}(\tau) similar to eq. (30) for \delta{\mbox{\boldmath{x}}}(\tau). To do that, we need to find the derivative of the projected deviation, $r$, which is given by eq. (32). Thus, we need to take the derivatives of {\hat{\mbox{\boldmath{B}}}} and then $P_{\perp}$ with respect to $\tau$:

[TABLE]

where $\mathrm{T}$ denotes matrix transpose. Let us also write down the following useful identities:

[TABLE]

The second line above follows from $2\hat{B}_{i}\nabla_{j}\hat{B}_{i}=\nabla_{j}\hat{B}^{2}=\nabla_{j}1=0$. These properties of {\hat{\mbox{\boldmath{B}}}} and its derivatives will lead to important simplifications below, which is the reason we wrote the field-line equation (eq. (29)) using the normalized magnetic field, and not the magnetic field itself.

Combining equations (31), (32), (33) and (34), we can write:

[TABLE]

where we used the fact that $P_{\perp}^{2}=P_{\perp}$. The term $(\mathrm{I}-P_{\perp})M^{\mathrm{T}}$ above takes into account the rotation of the plane orthogonal to {\hat{\mbox{\boldmath{B}}}} due to the changing direction of {\hat{\mbox{\boldmath{B}}}}, which in turn is due to the curvature of the field lines (cf. eq. (50) and the discussion around it).

3.2 Evolution of the field-line deviation vector in the normal plane

Let us remind the reader what the end goal is for this paper. It is obtaining local approximations to the squashing factor rate as well as to the field-line coiling rate. In order for those scalars to be locally defined around some field-line location {\mbox{\boldmath{x}}}(\tau), they can depend only on the first few derivatives of \delta{\mbox{\boldmath{x}}}(\tau), and {\mbox{\boldmath{r}}}(\tau) in particular. If all derivatives mattered, that would allow one to construct a Taylor expansion of \delta{\mbox{\boldmath{x}}}(\tau) (for analytic \delta{\mbox{\boldmath{x}}}(\tau)) along a non-infinitesimal field-line length, which will render the quantities non-local. As we will see later on, we in fact need only the first derivative of \delta{\mbox{\boldmath{x}}}(\tau) to be able to construct those quantities. Therefore, as the coiling and squashing rates are scalar quantities, we can conclude that they can only depend on invariant quantities constructed from the appropriate index contractions of the matrix $M_{ij}$ (entering in eq. (31)) with itself and with $\hat{B}_{i}$. That already includes possible contractions with $\tilde{m}$ (eq. (35)) as it is constructed out of $M$ and {\hat{\mbox{\boldmath{B}}}}. Thus, in what follows we are going to focus exclusively on the first derivative of the transverse deviation vector, {\mbox{\boldmath{r}}}(\tau).

That derivative does not necessarily lie in the plane normal to {\hat{\mbox{\boldmath{B}}}}(\tau). To see that, one can multiply both sides of eq. (35) by $P_{\perp}$ to find that generally:

[TABLE]

The reason for the inequality is that, by construction, {\mbox{\boldmath{r}}}(\tau) always lies in the plane normal to {\hat{\mbox{\boldmath{B}}}}, and that plane changes orientation with {\hat{\mbox{\boldmath{B}}}}.

However, in this paper we are only interested in transverse field-line deviations, disregarding the curvature of the field lines (cf. eq. (50)), which arises from the $\tau$-dependent {\hat{\mbox{\boldmath{B}}}} direction. So, to study those transverse deviations, we can pick a $\tau$-dependent 2D basis (subject to certain constraints; see below) which spans the changing normal plane, and write $r$ in that 2D transverse basis. Let us denote the components of $r$ in one of those (non-unique; see below) $\tau$-dependent 2D transverse orthonormal bases as: $r^{\diamondsuit}_{i}=(r_{x},r_{y})$. We can trivially extend that basis to 3D by including {\hat{\mbox{\boldmath{B}}}} as the third orthogonal unit vector, so that the trivial extension of $r^{\diamondsuit}_{i}$ to 3D is $(r_{x},r_{y},0)$ and the components of {\hat{\mbox{\boldmath{B}}}} in that basis are $(0,0,1)$.

By construction, the components $r^{\diamondsuit}_{i}$ must be related to the components in the fixed Cartesian basis by a $\tau$-dependent rotation matrix, $R_{ij}(\tau)$. Extracting only the first two non-zero components, we can write (with the sum written out explicitly to highlight the upper limit):

[TABLE]

The matrix $R$ is not unique since we can always make additional rotations in the normal plane, while leaving the components of {\hat{\mbox{\boldmath{B}}}} equal to (0,0,1) in the rotated frame. We address that issue in Section 3.3.

In order to make progress, we find it simpler to undo the rotation on the right-hand-side above up to some $\tau_{0}$ using the constant $R^{\mathrm{T}}(\tau_{0})$ matrix (which equals $R^{-1}(\tau_{0})$ as $R$ is orthogonal). Let us denote the components of vectors and operators expressed in that rotated coordinate system with a superscript ${(\tau_{0})}$. Then the components, $r^{{(\tau_{0})}}_{i}$, of $r$ in that 3D basis are given by:

[TABLE]

As we discussed in the beginning of this section, we will be focusing only on the first $\tau$ derivative of {\mbox{\boldmath{r}}}(\tau) around some $\tau_{0}$. For infinitesimal $\delta\tau$, we can write that derivative in finite difference form, which implies that we need to focus only on {\mbox{\boldmath{r}}}(\tau_{0}) and {\mbox{\boldmath{r}}}(\tau_{0}+\delta\tau). The former can be easily obtained in component form from the equation above: $r^{(\tau_{0})}_{i}(\tau_{0})=r_{i}(\tau_{0})$; while the latter can be written out as:

[TABLE]

where the last equation defines the infinitesimal rotation matrix $\delta R$ (assuming $R(\tau)$ is picked to vary smoothly). In what follows, we write an equation for $\delta R$.

The rotation matrix $\delta R$ infinitesimally rotates the basis vectors for the quantities evaluated at $(\tau_{0}+\delta\tau)$, but leaves those evaluated at $\tau_{0}$ unchanged since $\delta R(\delta\tau=0,\tau_{0})=\mathrm{I}$ as can be seen from the definition above. We loosely refer to this mixture of basis vectors at the two different locations, $\tau_{0}$ and $(\tau_{0}+\delta\tau)$, along the field line as the bent basis, as it “bends” from $\tau_{0}$ to $\tau_{0}+\delta\tau$ in such a way so as to leave the components of {\hat{\mbox{\boldmath{B}}}} unchanged (see the discussion above). In other words, $\delta R$ is chosen such that $\hat{B}_{i}^{(\tau_{0})}(\tau_{0}+\delta\tau)=\delta R_{ij}\hat{B}_{j}(\tau_{0}+\delta\tau)=\hat{B}_{i}(\tau_{0})$.

Using the Rodrigues’ rotation formula999One can use the Rodrigues’ rotation formula to find the rotation matrix ($\tilde{R}$) in the plane defined by two unit vectors, $a$ and $b$, that takes the components $b_{i}$ to the components $a_{i}$. That matrix is given by:

\tilde{R}_{ij}=\mathrm{I}_{ij}+2a_{i}b_{j}-\frac{(a_{i}+b_{i})(a_{j}+b_{j})}{1+{\mbox{\boldmath{a}}}\cdot{\mbox{\boldmath{b}}}}\ .

It is straightforward to check that indeed $\tilde{R}\tilde{R}^{\mathrm{T}}=\mathrm{I}$, and $\tilde{R}_{ij}b_{j}=a_{i}$. We use this result to obtain eq. (40)., after a bit of algebra, one can show that the sought for rotation matrix transforming vector components to the bent frame is given by:

[TABLE]

with all operators evaluated at $\tau_{0}$. Indeed, checking that $\delta R\delta R^{\mathrm{T}}=\mathrm{I}$ to linear order in $\delta\tau$ is a straightforward exercise, which guarantees that $\delta R$ is a rotation matrix. Then, using eq. (33), one can check that to linear order in $\delta\tau$:

[TABLE]

where we used eq. (34).

Note that similar to $R(\tau)$, $\delta R(\delta\tau,\tau_{0})$ is not unique as we can additionally make a rotation in the plane normal to {\hat{\mbox{\boldmath{B}}}}. We address that ambiguity in Section 3.3. However, let us first use the rotation matrix $\delta R$ to calculate $r^{(\tau_{0})}_{i}(\tau_{0}+\delta\tau)$. First, note that eq. (35) can be written with finite differences as:

[TABLE]

From that, we obtain:

[TABLE]

where we used equations (40) and (34). Therefore, in the bent basis, we can write the first derivative of the components of {\mbox{\boldmath{r}}}(\tau) at $\tau_{0}$ as:

[TABLE]

where we used the fact that $r^{(\tau_{0})}_{i}(\tau_{0})=r_{i}(\tau_{0})$ since in the bent frame, quantities evaluated at $\tau_{0}$ are not rotated by $\delta R$. We can see that $m^{(\tau_{0})}$, defined in the last equality above, is simply the transverse part of $M$, which is reassuringly reasonable given the non-trivial way we arrived at it.

Now we are ready to write an equation for the first derivative of the components of {\mbox{\boldmath{r}}}(\tau) at $\tau_{0}$ in the 2D transverse basis. Let us apply the constant rotation matrix $R(\tau_{0})$, which is independent of both $\tau$ and $\delta\tau$, to both sides of eq. (44):

[TABLE]

From equations (37) and (38), one can write the vector components $R_{ij}(\tau_{0})r_{j}^{(\tau_{0})}$, entering on both sides above, as:

[TABLE]

with the $i=3$ component vanishing. Moreover, since $m^{(\tau_{0})}$ is the transverse part of $M$ (eq. (44)), that means that $R_{ij}(\tau_{0})m^{(\tau_{0})}_{jk}\left(R^{\mathrm{T}}(\tau_{0})\right)_{kl}$ entering in eq. (45) has non-zero elements only in the upper-left $2\times 2$ block since {\hat{\mbox{\boldmath{B}}}} has components (0,0,1) in the normal-plane basis extended to 3D (see discussion above). We denote that non-zero block as the 2-by-2 matrix $m$:

[TABLE]

In other words, the trivial extension of $m$ (by padding it with zeros) to 3D equals $m^{(\tau_{0})}$ up to a 3D rotation rotating the components of {\hat{\mbox{\boldmath{B}}}}(\tau) to (0,0,1).

Furthermore, it is easy to see from eq. (44) that the $\tau$ derivative of the components of {\mbox{\boldmath{r}}}(\tau) evaluated at $\tau_{0}$ in the bent basis is entirely in the normal plane as well:

[TABLE]

where we used eq. (34). The above result is in contrast to what we obtained in the fixed Cartesian basis, eq. (36). Thus, once we rotate that derivative with $R(\tau_{0})$ (left-hand-side of eq. (45)), that derivative will have a vanishing third component, and its first two components can be written as $dr_{i}^{\diamondsuit}/d\tau$ with $i=1,2$.

To conclude, after rotating with $R(\tau_{0})$, the rotated $m^{(\tau_{0})}$, $r_{i}^{{(\tau_{0})}}$ and $dr_{i}^{(\tau_{0})}/d\tau$ are explicitly in the 2D subspace of the normal plane. From now on, equations involving quantities evaluated in that 2D transverse subspace will be denoted with a $\diamondsuit$ at the start of each line, instead of modifying the symbol for each such quantity. Collecting the results of the discussion above, in the transverse 2D subspace, eq. (45) can be written as:

[TABLE]

with $m$ given by eq. (47). The above equation gives the first derivative of the components of $r$ entirely in the 2D transverse subspace, which was what we set out to find.

As a side note, earlier we argued that the we are not interested in the changes in the orientation of the normal plane, which intuitively must be due to the curvature ($\kappa$) of the field line {\mbox{\boldmath{x}}}(\tau). Let us check that explicitly. The part of $M$ which is not captured by $m^{(\tau_{0})}$ is $M(\mathrm{I}-P_{\perp})$ since from eq. (34) we have $m^{(\tau_{0})}=P_{\perp}MP_{\perp}=MP_{\perp}$. We can write that as:

[TABLE]

where we used eq. (29) and the definition of $M$ (eq. (31)), as well as the Frenet-Serret formulas from differential geometry. In the above equation \hat{\mbox{\boldmath{N}}} is the normal unit vector to the field line, which along with the tangent vector (in our case, that is {\hat{\mbox{\boldmath{B}}}}, which runs contrary to standard differential geometry notation) and the binormal unit vector, form the basis for the Frenet-Serret frame of a curve. The torsion of the curve does enter101010See equation (3.140) of the manuscript “Fundamentals of Plasma Physics” by J. D. Callen, archived at https://web.archive.org/web/20170622161238/http://homepages.cae.wisc.edu/~callen/chap3.pdf in the antisymmetric part of $m^{(\tau_{0})}$, and therefore is part of the force-free parameter $\alpha$ as well (cf. eq. (94) and surrounding discussion). But it will not be of special significance to our analysis, since it is defined in terms of the changes to the field-line tangent and normal vectors, which we explicitly disregard once we focus on the field line motions in the normal plane (but see footnote 15). And indeed, torsion and curvature are properties of a single field line, while in this paper we explore the transverse behavior of an ensemble of nearby field lines.

3.3 Specifying the bent basis uniquely

Before we move on, there is one more issue to address. As mentioned above, the bent basis (and hence $\delta R$) is not unique as one can make an infinitesimal rotation in the normal plane and still keep the difference between the components $r_{i}^{(\tau_{0})}(\tau_{0}+\delta\tau)$ and $r_{i}(\tau_{0}+\delta\tau)$ small. To fix the bent basis uniquely, we require that the infinitesimal angle ($\delta\phi$) of rotation of \hat{\mbox{\boldmath{r}}} in the normal plane between $\tau$ and $(\tau+\delta\tau)$ is the same, whether it is calculated using the components of \hat{\mbox{\boldmath{r}}} in the fixed 3D Cartesian basis or using the components of \hat{\mbox{\boldmath{r}}} in the 2D transverse basis (and therefore in the bent basis as well) if we pretended that basis was held fixed (see below). In the 3D Cartesian basis, that angle can be found using (cf. eq. (82)):

[TABLE]

where $\epsilon_{ijk}$ is the Levi-Civita symbol. In the last equality above, we linearized $\hat{r}_{k}(\tau_{0}+\delta\tau)$ and applied the antisymmetry of $\epsilon_{ijk}$, combined with:

[TABLE]

where in the last step we used eq. (35). Plugging in the expression for $\tilde{m}$ (eq. (35)) in eq. (51) and using eq. (34), we find (cf. eq. (83)):

[TABLE]

We want the above expression to match the “naive” angle, $\delta\phi^{2d}$, that one would calculate from the components of $r$ in the 2D transverse basis assuming the basis vectors were held constant:

[TABLE]

Since (as discussed above) the 2D basis spans a plane in the bent basis orthogonal to {\hat{\mbox{\boldmath{B}}}}(\tau), one can rewrite the equation above in the bent basis:

[TABLE]

where we used the fact that the bent basis only changes the components of quantities evaluated at ($\tau_{0}+\delta\tau$). The equation above can be easily checked by writing it in a basis such that the components of {\hat{\mbox{\boldmath{B}}}} equal (0,0,1) (which is what the rotation matrix $R(\tau_{0})$ achieves in eq. (45)). Plugging the expression for $m^{(\tau_{0})}$ (eq. (44)) above, using P_{\perp}{\mbox{\boldmath{r}}}={\mbox{\boldmath{r}}}, and comparing with eq. (53), we find that indeed:

[TABLE]

Thus, the rotation matrix $\delta R$ is the unique infinitesimal rotation matrix that does not introduce any spurious rotations of the components of $r_{i}$ in the normal plane. Thus, the normal-plane basis is uniquely specified by the requirement that the rate of rotation of individual neighboring field lines around the reference field line is independent of whether that rate is calculated using the standard non-local twist rate (cf. equations (51) and (82)), or in the basis spanning the normal plane. We will refer to that property (eq. (56), combined with equations (53) and (55)) of the bent basis (defined by $\delta R$) as the basis being minimally rotated.

3.4 Invariants

As we discussed in the beginning of Section 3.2, we are interested in the local invariants one can build out of the matrices introduced so far. So, let us show that the following useful relationship between the invariants of $M$, $\tilde{m}$, $m^{(\tau_{0})}$ and $m$ (defined in equations (31), (35), (44) and (49), respectively) holds111111A standard result in linear algebra states that the eigenvalues, and hence the determinant, of a square matrix can be written as a function of the traces of the various powers of that matrix.:

[TABLE]

for integer powers $n\geq 1$. The last equality follows trivially from the fact that $m^{(\tau_{0})}$ is the 3D extension of $m$ (see the discussion around eq. (47)).

Below we show the first two equalities in eq. (57) for $n=1,2$. The result for the rest of the powers follows analogously. Indeed, we can use eq. (34) and the cyclic property of the trace to write:

[TABLE]

From eq. (57), we can conclude that two of the eigenvalues of $M$, $\tilde{m}$, $m^{(\tau_{0})}$ and $m$ match, while the third eigenvalue of $M$, $\tilde{m}$ and $m^{(\tau_{0})}$ is zero. Indeed, from the definition of the projection operator, it is easy to see that $m^{(\tau_{0})}=P_{\perp}MP_{\perp}$ has an eigenvalue of zero, corresponding to an eigenvector equal to {\hat{\mbox{\boldmath{B}}}}. Similarly for $\tilde{m}$:

[TABLE]

where we used eq. (34). Note, however, that even though $M$ must have a right eigenvector with a corresponding zero eigenvalue, that does not necessarily need to be {\hat{\mbox{\boldmath{B}}}} (which is, however, a left eigenvector of $M$ as can be seen from (34): {\hat{\mbox{\boldmath{B}}}}^{\mathrm{T}}M=0).

4 The squeeze factor

In this section we obtain an approximation to the squashing factor as an integral over purely local quantities. We call that approximation the squeeze factor in reference to the decomposition we used in Section 2.2.2. Indeed, the results of this section parallel those for the simple 2D planar flow presented in Section 2.2.2. Yet, our approach in this section will be markedly different; thus, offering an alternative way to understand our results.

To derive an expression for the squeezing factor, we start by writing down the general equations for the squashing factor, $Q$. In doing so, we schematically follow Titov (2007); however, we simplify the analysis by restricting it only to orthonormal basis vectors (see below).

4.1 Review of the standard squashing factor calculation

We start with a reference field line (dashed line in Fig. 5) for which we want to calculate $Q$. That quantity is constant along any field line once the footpoints of the field line have been chosen. Let us pick three points along the length of the reference field line at values of the field-line length parameter ($\tau$) equal to $\tau_{0}$, $\tau_{1}$ and $\tau_{2}$, with points $\tau_{1}$ and $\tau_{2}$ infinitesimally close to each other. In this section, we focus only on points $\tau_{0}$ and $\tau_{1}$, which we will treat as the footpoints of the reference field line.

We can launch a set of neighboring field lines, intersecting in a circle the plane normal to the reference field line at $\tau_{0}$. The cross-section of that flux tube of neighboring field lines is generally deformed into an ellipse as one follows the field lines (see Fig. 5). The amount of that deformation is quantified by $Q$ as discussed below.

Let us introduce the following notation: \mbox{\boldmath{r}}^{(i)}\equiv\mbox{\boldmath{r}}(\tau_{i}) (with $i=0,1,2$). From eq. (32), we know that each of those vectors lies in the plane orthogonal to the reference field line at the respective {\mbox{\boldmath{x}}}(\tau_{i}) (see Fig. 5). We will refer to each one of those normal planes as plane ($i$) from now on.

Let the map between the components of \mbox{\boldmath{r}}^{(1)} in an orthonormal basis spanning plane (1) and the components of \mbox{\boldmath{r}}^{(0)} in an orthonormal basis spanning plane (0) be given by:

[TABLE]

with $J$ being a $2\times 2$ matrix. As we did for eq. (49), we place a $\diamondsuit$ next to all equations written in the varying 2D basis of the normal planes as a reminder that we are not working in the fixed 3D Cartesian basis.

The squashing factor (Titov, 2007; Pariat & Démoulin, 2012), $Q$, between $\tau_{0}$ and $\tau_{1}$ can then be written as:

[TABLE]

The squashing factor is invariant under rotations in either plane (0) or (1) (Titov, 2007). For an easy way to see this, note that under a rotation in plane (0), we have $r^{\prime(0)}=R_{(0)}r^{(0)}$ and under a rotation in plane (1) we have $r^{\prime(1)}=R_{(1)}r^{(1)}$, with $R_{(0,1)}$ being rotation matrices. Substituting those expressions in $r^{(1)}=Jr^{(0)}$, we can write $r^{\prime(1)}=J^{\prime}r^{\prime(0)}$ with $J^{\prime}=R_{(1)}JR_{(0)}^{-1}$, which is the transformation equation for $J$ under rotations in both planes (0) and (1). Plugging $J^{\prime}$ in eq. (61), we obtain:

[TABLE]

where we used the cyclic property of the trace, as well as the fact that the rotation matrices are orthogonal. Thus, we can see that indeed Q is independent of both $R_{(0)}$ and $R_{(1)}$.

We find it useful to introduce yet another $2\times 2$ matrix:

[TABLE]

It is the matrix of the cross-sectional ellipse, which is formed at the intersection of the flux tube of neighboring field lines with plane (1). To see that, note that (up to a choice of units) one can write \left(\mbox{\boldmath{r}}^{(0)}\right)^{2}=1 for all tube field lines, since, by construction, we picked a flux tube that intersects plane (0) in a circle. Using eq. (60), we can then write:

[TABLE]

from which one can see that $E$ is indeed the matrix of the cross-sectional ellipse in plane (1). Note that the fact that the quadratic form (in {\mbox{\boldmath{r}}}^{(1)}) specified by eq. (64) is an ellipse can be seen from the inequality $\mathrm{det}\left(E\right)=1/\mathrm{det}\left(J\right)^{2}>0$. Comparing the above equation with the implicit equation for an ellipse written in canonical form, and using the rotation invariance of the trace and determinant of $E$ (which can be confirmed following the proof of the invariance of $Q$ above), one can write:

[TABLE]

where $s_{a}$ and $s_{i}$ are the semi-major and semi-minor axes of the cross-sectional ellipse, respectively.

Using the definition of $E$, eq. (63), we can express $Q$ from eq. (61) as:

[TABLE]

In writing the above equation, we used the fact that since $E$ is a $2\times 2$ matrix, we have $\mathrm{tr}(E^{-1})=\mathrm{tr}(E)/\mathrm{det}(E)$. Combining equations (65) and (66), we can write the squashing factor as the sum:

[TABLE]

where $s_{a}/s_{i}\geq 1$ is the aspect ratio of the cross-sectional ellipse in plane (1). That result makes the invariance of $Q$ under rotations in plane (1) explicit (Titov, 2007).

4.2 The derivative of the squashing factor with respect to the footpoint locations

In this section, we will take the derivative of the squashing factor with respect to the location of one of the footpoints of the reference field line for which $Q$ is calculated. To do that, let us move one of the footpoints an infinitesimal distance, $\delta\tau$, along the reference field line from $\tau_{1}$ to $\tau_{2}=\tau_{1}+\delta\tau$ (see Fig. 5) and use the results from the previous section.

Plugging eq. (60) in eq. (49) evaluated at $\tau_{1}$, and using the fact that \mbox{\boldmath{r}}^{(0)} was arbitrarily chosen, we can write:

[TABLE]

Taking the derivative of $E$ in eq. (63) with respect to $\tau$, and combining with eq. (68), we find

[TABLE]

Let us now take the derivative of $Q$ with respect to $\tau$. Using equations (66) and (69), we find:

[TABLE]

To obtain the above equation, we used the fact that $E$ is a symmetric matrix (as can be seen from eq. (63)), along with the identity $\mathrm{det}\left(E\right)=\exp\left[\mathrm{tr}\left(\ln(E)\right)\right]$, which allows us to write the familiar equation for the derivative of a determinant: $d(\mathrm{det}\left(E\right))/d\tau=\mathrm{det}\left(E\right)\mathrm{tr}\left(E^{-1}dE/d\tau\right)$.

Note that in eq. (70) we can replace $m$ with its symmetric part as any anti-symmetric part of $m$ has no contribution as $E$ is symmetric. However, eq. (70) can still depend implicitly on the anti-symmetric part of $m$ through $E$, given by eq. (69).

The system of equations (69) and (70) specifies the evolution of $Q(\tau)$ as long as $m(\tau)$ is known, and as long as we specify initial conditions for $E$ and $Q$. Initially, we can set $\tau=\tau_{0}$, for which the tube of neighboring field lines has a circular cross-section by construction. Thus, after setting $\tau_{1}=\tau_{0}$ in eq. (64), we can read off the initial condition for $E$: $E_{ij}(\tau_{0})=\mathrm{I}_{ij}$. From eq. (66), the initial condition for $Q$ is, therefore, given by $Q(\tau_{0})=2$. To close the system of equations, we would need an expression for $m(\tau)$. However, we will avoid writing an explicit expression for the $2\times 2$ matrix $m$ and instead in the next section we will use the results of Section 3.4 for its invariants.

The right-hand side of eq. (70) is the local logarithmic contribution to the squashing factor, which one may be tempted to interpret as the local (logarithmic) squashing factor rate. However, it depends on the cross-sectional ellipse matrix $E$, which in turn depends on the history of how the flux tube of neighboring field lines was deformed up to $\tau$. Thus, this rate is not a truly local quantity. In what follows, we show a way around that issue.

4.3 Constructing an approximate local squashing factor rate

In this section, we manipulate eq. (70) to obtain the squeeze factor ($\mathcal{Z}$): an estimate of the squashing factor as an integral over purely local quantities ($\rho_{\mathcal{Z}}$):

[TABLE]

We use the symbol $\mathcal{Z}$ to denote that approximation in anticipation of recovering the results of Section 2.2.2.

In order to construct the approximate $\mathcal{Z}$, we explored a variety of ad hoc modifications to eq. (70) which result in a purely local right-hand side. To decide which modification is most “reasonable”, we performed numerical experiments with realistic simulated data, the final results of which are presented in (Savcheva, 2019), where we show maps of $\mathcal{Z}$ and $Q$ which are morphologically comparable. To keep with the spirit of this paper, however, below we present an analytical motivation for our choice for $\mathcal{Z}$ (but see Section 4.5 for a discussion of an alternative), developed after the numerical experiments informed our intuition. The final result of this section matches that of Section 2.2.2 (eq. (25)). However, that result was obtained for constant $m$. Here we derive that result using kinematic arguments for tube-like flux surfaces undergoing deformations by slowly-varying $m$ in the presence of non-zero field-line curvature.

We will treat the following two cases separately: $m$ having two complex conjugate eigenvalues, and $m$ having two real eigenvalues. Let us start with the latter.

If $m(\tau)$ has two real eigenvalues and is changing only slowly, then the cross-sectional ellipse will be preferentially stretched in the direction of the (normalized) eigenvector \hat{\mbox{\boldmath{v}}}_{+} of $m$, which corresponds to its largest eigenvalue $\lambda_{+}$ (such that m\hat{\mbox{\boldmath{v}}}_{+}=\lambda_{+}\hat{\mbox{\boldmath{v}}}_{+}). Indeed, we are going to assume that the semi-major axis of the cross-sectional ellipse is collinear with \hat{\mbox{\boldmath{v}}}_{+}. The ellipse matrix, $E$, can then be approximated locally as $E_{\mathrm{loc.}}\sim E$:

[TABLE]

In the above equation, \hat{\mbox{\boldmath{v}}}_{\perp} is defined to be orthonormal to \hat{\mbox{\boldmath{v}}}_{+}, and is therefore not necessarily an eigenvector of $m$. With the construction above, \hat{\mbox{\boldmath{v}}}_{+}^{\mathrm{T}}E_{\mathrm{loc.}}\hat{\mbox{\boldmath{v}}}_{+}=1/s_{a}^{2}, \hat{\mbox{\boldmath{v}}}_{\perp}^{\mathrm{T}}E_{\mathrm{loc.}}\hat{\mbox{\boldmath{v}}}_{\perp}=1/s_{i}^{2} and \hat{\mbox{\boldmath{v}}}_{+}^{\mathrm{T}}E_{\mathrm{loc.}}\hat{\mbox{\boldmath{v}}}_{\perp}=0, which show that indeed the cross-sectional ellipse is oriented as described above. Note that $R_{+}$ is an orthogonal matrix as its columns are the orthonormal vectors \hat{\mbox{\boldmath{v}}}_{+} and \hat{\mbox{\boldmath{v}}}_{\perp}, and is, therefore, a rotation matrix.

Next we substitute $E_{\mathrm{loc.}}$ for $E$ in eq. (70). First, we need to calculate:

[TABLE]

where we used the cyclic property of the trace. Note that121212As $R_{+}^{\mathrm{T}}mR_{+}$ equals an upper triangular matrix, as seen in eq. (73), it is the real Schur form of $m$. The off-diagonal element can be interpreted as a finite shear due to the non-orthogonal eigenvectors of $m$, which causes a non-zero force-free parameter $\alpha$. See Section 2.2.2 for further discussion. $R_{+}^{\mathrm{T}}mR_{+}$, which appears above, is a rotation of $m$, and therefore its eigenvalues must match those of $m$. Thus, we can conclude that \hat{\mbox{\boldmath{v}}}_{\perp}^{\mathrm{T}}m\hat{\mbox{\boldmath{v}}}_{\perp}=\lambda_{-}, where $\lambda_{-}$ is the smaller of the two eigenvalues of $m$. We can then write:

[TABLE]

Substituting in eq. (70) and simplifying, we obtain:

[TABLE]

In order for the above expression to be purely local, we need to make it independent of the semi-major and semi-minor axis of the cross-sectional ellipse. A trivial way to do that is to maximize the above expression, which entails taking $s_{i}\ll s_{a}$. Indeed, this inequality holds for large $\tau$ when $m$ is constant.

Thus, we are assuming that the cross-sectional ellipse is infinitely squeezed, with a semi-major axis oriented along the eigenvector {\mbox{\boldmath{v}}}_{+} of $m$. With that assumption, we can write the squeeze factor $\mathcal{Z}$ – our approximation to the squashing factor – as:

[TABLE]

This expression is completely local as sought for. Analytically, it is rather neat as it is simply the difference in the logarithmic rates with which the field lines are being pulled along the eigenvectors of the transverse flow set up by $m$ (cf. eq. (20)).

Note that the eigenvalue difference above depends on the antisymmetric part of $m$. This is due to our assumption that $m$ is slowly varying; and therefore, the squeezed cross-sectional ellipse ends up oriented along one of the eigenvectors of $m$ and not, for example, its symmetrized counterpart (see Section 4.5). Thus, one may worry that the above expression is sensitive to physically irrelevant infinitesimal rotations between planes (1) and (2), resulting from poor choice of basis. Fortunately, $m$ has non-zero eigenvalues matching $m^{{(\tau_{0})}}$ (Section 3.4), which in turn was written in the minimally rotated bent basis (Section 3.3). This guarantees that the components of $r$ expressed in the normal plane rotate at the same rate as those expressed in the fixed 3D Cartesian basis, i.e. the antisymmetric part of $m$ captures only rotations that are indeed physical. For further discussion, see eq. (56) and the discussion around it.

Equation (76) is valid for real eigenvalues of $m$. For two complex conjugate eigenvalues, the expression would result in an imaginary contribution to $\ln(Q)$ and is thus invalid. To deal with that case, we can use the fact that the eigenvalues will then have equal real parts. Thus, the transverse flow set up by $m$ is a stable or unstable spiral (or center), which periodically rotates the cross-sectional ellipse, returning it to the same aspect ratio after a full rotational period (assuming constant $m$). Therefore, from eq. (67) we can see that $Q$ would on average remain unaffected by such a flow. Thus, a reasonable choice for $\rho_{\mathcal{Z}}$ would be to set it to zero in this case. That corresponds to eq. (27) in Section 2.2.3.

Another way of motivating setting $\rho_{\mathcal{Z}}$ to zero for complex eigenvalues of $m$ is to note that the initial orientation of the cross-sectional ellipse can be arbitrary as there is no preferred “stretching” direction set up by the flow of $m$. Thus, for $E_{\mathrm{loc.}}$ we can again use eq. (72) with $R_{+}$ this time being a rotation matrix rotating the cross-sectional ellipse at a random angle $\theta$. Plugging that ansatz into (70) and averaging over the angle $\theta$, after some algebra we again find:

[TABLE]

If the localized periodic squashing of the flux tube is of interest, however, one can use, eq. (28) for complex eigenvalues of $m$ (see Section 2.2.3).

Equations (76) and (77) are the final result of this section. Next, we write down the eigenvalue difference in eq. (76) using the local properties of the magnetic field.

4.4 The squeeze factor

In order to evaluate eq. (76), we need to calculate the eigenvalues of the $2\times 2$ matrix $m$. As we showed in Section 3.4, the two non-zero eigenvalues of $M$ must match the two eigenvalues of $m$. The difference between those two eigenvalues for the $2\times 2$ matrix $m$ can be written as $\lambda_{+}-\lambda_{-}=\sqrt{2\,\mathrm{tr}(m^{2})-\left(\mathrm{tr}(m)\right)^{2}}$, which can in turn be written as $\lambda_{+}-\lambda_{-}=\sqrt{2\,\mathrm{tr}(M^{2})-\left(\mathrm{tr}(M)\right)^{2}}$ as the third eigenvalue of $M$ is zero. In Appendix A (cf. eq. (A3)), we show that in curvilinear coordinates, $M_{\ j}^{i}$ is (unsurprisingly) given by the covariant derivative $\nabla_{j}\hat{B}^{i}$. Thus, in any coordinate system we can finally write equations (76) and (77) as131313It may be of interest to point out the relationship between $Q$ and the Lyapunov exponents of an attractor flow specified by some {\hat{\mbox{\boldmath{B}}}}(\tau). For a particular attractor, the Lyapunov exponents ($l_{i}$) are defined through the eigenvalues ($\Lambda_{i}(\tau)>0$) of the symmetric matrix $J^{\mathrm{T}}(\tau)J(\tau)$ as:

l_{i}=\lim\limits_{\tau\to\infty}\frac{\ln\big{(}\Lambda_{i}(\tau)\big{)}}{2\tau}\ .

Assuming that the footpoints are $\tau\to\infty$ apart along the field lines; and that the field is defined over a region, where taking that limit is meaningful; and that the field-line flow asymptotically is described by an attractor – all assumptions that clearly break down for magnetic fields in the solar corona – one can then use eq. (61) to find $Q=(\Lambda_{1}+\Lambda_{2})/\sqrt{\Lambda_{1}\Lambda_{2}}$. Therefore, in the limit $Q\gg 1$, we obtain $\ln Q(\tau)\sim\tau|\Delta l|$, where $\Delta l$ is the difference between the two Lyapunov exponents of the flow in the normal plane. Comparing with eq. (78), we can see that $|\Delta l|$ can then be approximated by the $\tau$-average of $\Re(\Delta\lambda)$. For small $\tau$, assuming that we still have $Q\gg 1$, one can then identify $\rho_{\mathcal{Z}}$ as the difference between the local Lyapunov exponents (Abarbanel et al., 1991). :

[TABLE]

with $\Re(x)$ denoting the real part of $x$. To set the minimum value of $\mathcal{Z}$ to that of ${Q}$, we use $\mathcal{Z}(\tau_{0})=2$ as an initial condition for eq. (78). Note that for an imaginary eigenvalue difference, if one is interested in taking into account the localized periodic squashing of the flux tube then, instead of setting $\rho_{\mathcal{Z}}$ to zero in that case, one can use eq. (28) (see Section 2.2.3).

Note that the equation above is invariant of the direction of integration along the reference field line as required for $Q$ (Titov et al., 2002; Titov, 2007). Moreover, since $\rho_{\mathcal{Z}}\geq 0$, the contributions to the squeeze factor, $\mathcal{Z}$, cannot be undone as one integrates over a field line, especially if one uses eq. (28) in the case of an imaginary eigenvalue difference. That is in contrast to $Q$, which can remain small in the presence of large localized squashing, as long as the flux tube around the field line in question is squashed and then un-squashed between its footpoints. As $Q$ is growing exponentially large in the vicinity of HFTs and nulls, a large $\rho_{\mathcal{Z}}$ can be used as an indicator for the presence of those features.

The first equality in eq. (78) is approximate in a sense that can be quantified by comparing equations (70) and (76), with the difference arising from using the cross-sectional ellipse matrix ($E$) in $Q$, while using its local approximation ($E_{\mathrm{loc.}}$) in $\mathcal{Z}$. Thus, one can think of $\rho_{\mathcal{Z}}$ integrated over a stretch of $\tau$ as the extra logarithmic squashing of an already infinitely squeezed ellipse (that is, the quadratic form, corresponding to $E_{\mathrm{loc.}}$) that is pointed along the eigenvectors of $m$. Given the choice of $E_{\mathrm{loc.}}$, described in Section 4.3, that orientation of the cross-sectional ellipse can be assumed to be caused by a slowly evolving $m$, which preferentially squashes flux tubes along its eigenvectors; or can be thought of as being produced by a deliberate choice of footpoint location ($\tau_{0}$), which produces a tube with that particular cross-sectional ellipse at some $\tau$.

For an alternative justification for approximating the squashing factor with $\mathcal{Z}$, compare the result above with the discussion in Section 2.2.2, and eq. (25) in particular. Equation (78) is the final result of this section. It is the generalization of eq. (25) for field lines with non-zero curvature.

4.5 The symmetrized squeeze rate

As we discussed previously, there is an inherent ambiguity of how one quantifies the squashing rate using local quantities. In the previous section, we presented a kinematic argument for our choice of $E_{\mathrm{loc.}}$, eq. (72). We assumed that $m$ varies slowly enough, so as to guarantee that the cross-sectional ellipse is stretched along the eigenvector of $m$ corresponding to the largest eigenvalue.

One may wonder how our results would change if we picked $R_{+}$ in eq. (72) to be a rotation matrix with a rotation angle which maximizes the resulting $d\ln\mathcal{Z}/d\tau$. After lengthy algebra, one can confirm that such a choice still recovers equations (75) and (76), but this time the eigenvalue difference ($\Delta\lambda_{s}$) that would enter in (76) would correspond to the symmetrized matrix $m_{s}\equiv(m+m^{\mathrm{T}})/2$. Thus, this corresponds to an ellipse aligned with the eigenvectors of $m_{s}$, and not of $m$. It should be clear, however, that such an orientation for the cross-sectional ellipse can only happen by chance as a transient in a regime when $m$ is changing fast.

In the language of Section 2.1, $\Delta\lambda_{s}$ equals the difference between the eigenvalues of the rate-of-shear tensor ($\Sigma$; see eq. (3)) as a non-zero rate-of-expansion does not affect $\Delta\lambda_{s}$. Therefore, $\Delta\lambda_{s}$ gives the relative logarithmic rate of squeeze in the standard rotation-squeeze decomposition (first and second row of Fig. 3). Thus, the symmetrized squeeze rate, $\Delta\lambda_{s}$, is non-zero for elliptical flows. One can show that for any $2\times 2$ matrix $m$, $\Delta\lambda_{s}^{2}=\Delta\lambda^{2}+\alpha^{2}$, with the generalized force-free parameter, $\alpha$, being the difference between the two off-diagonal elements of $m$ (e.g. eq. (84)). Combining with equations (14) and (17) we arrive at:

[TABLE]

which is clearly different from the physically motivated localized rate of squashing that we obtained earlier for periodic flows: equations (26) and (28) in Section 2.2.3.

Similarly, a dependence of $Q$ on the eigenvalues of $m_{s}$ can be obtained from eq. (24) for the steady-state 2D flow (discussed in Section 2) for small $\Delta\tau$, such that $Q$ is not exponentially boosted. Expanding $Q$ in that equation to second order in $\Delta\tau$, one can show that $Q\approx 2+\Delta\tau^{2}\Delta\lambda_{s}^{2}$. So, clearly $\Delta\lambda_{s}$ can play a role in regions of low $Q$. However, in this paper, we are only interested in $Q\gg 1$ which, as we argue in Section 2.2.2, depends on the eigenvalues of the unsymmetrized $m$.

Despite these arguments against using the symmetrized eigenvalue difference, one may still consider using $\Delta\lambda_{s}$ in (76), especially if interested in regions that do not necessarily contribute to an exponential boosting of $Q$. Therefore, we did modify QSL Squasher (Tassev & Savcheva, 2017) to allow for the calculation of $\rho_{\mathcal{Z}}$ using $\Delta\lambda_{s}$ in (76) as an option.

5 Coiling number

In this section we write down our local approximation for the coiling rate of neighboring field lines in the normal plane. We relate it to the generalized force-free parameter, and we interpret the latter in two different ways: as an angle average of the angular rate of motion of neighboring field lines, and as a geometrically boosted coiling rate. The results of this section parallel those for the simple 2D planar flow presented in Section 2.2.1.

5.1 The generalized force-free parameter

For the purposes of calculating the local twist number (cf. our eq. (97), as well as eq. (16) of (Berger & Prior, 2006)) the force-free parameter, $\alpha$, can be generalized for any magnetic field as:

[TABLE]

Note that $\alpha$ can also be written as:

[TABLE]

The last term above vanishes because it is proportional to the dot product of two orthogonal vectors. Therefore, we obtain:

[TABLE]

Let us now calculate the angular rate $\omega$ at which neighboring field lines are rotating around a reference field line. The rate is given by (e.g. Berger & Prior (2006); Liu et al. (2016); cf. our eq. (51)):

[TABLE]

To simplify the above expression, let us use eq. (52) and express $\tilde{m}$ from (35). We find that the specific combinations of dot and cross products allow for the elimination of many of the terms. In the end, we obtain (cf. eq. (18) of Liu et al. (2016), as well as our eq. (53)):

[TABLE]

To compare the expression for $\omega$ with that for $\alpha$, let us pick a basis in which {\hat{\mbox{\boldmath{B}}}} has components $(0,0,1)$. In that basis, we can write $r$ as $(r_{x},r_{y},0)$. Plugging those into (81) and (83), we obtain:

[TABLE]

Clearly, $\omega$ depends on the location of the neighboring field line, while $\alpha$ does not. So, let us write $r_{x}=r\cos(\phi)$ and $r_{y}=r\sin(\phi)$ in polar coordinates and average $\omega$ over $\phi$. We obtain:

[TABLE]

Comparing with (84), we can finally write:

[TABLE]

The expression above is coordinate independent and represents the final result of this section. It allows us to interpret $\alpha$ as the average angular rate of motion of neighboring field lines around the reference field line passing through the location at which $\alpha$ is calculated. The average is taken over the neighboring field lines, assuming they are initially uniformly distributed in angle.

The above result mirrors the discussion around eq. (8) for vorticity (the fluid flow equivalent of $\alpha$) as a measure of the angular rate of rotation of a floating vane uniformly “sampling” (in angle) nearby streamlines. Thus, the angular average in eq. (86) is non-zero even for sheared saddle transverse flow of neighboring field lines, when none of the field lines make full turns around one another. Thus, $\alpha$ is a poor measure of local magnetic field line twist. To improve on it, we introduce the coiling rate below.

5.2 The coiling rate

Below we obtain our expression for the coiling rate. We start with the unambiguously defined rate of rotation of individual neighboring field lines in the normal plane. Then we show how the coiling rate is related to a specific averaging of the rate of rotation of individual field lines in a way which is distinct from the averaging that lead to the generalized force-free parameter (see eq. (86)).

So, let us start by exploring the behavior of the vector {\mbox{\boldmath{r}}}=(r_{x},r_{y}) in the plane orthogonal to {\hat{\mbox{\boldmath{B}}}}. That vector obeys a first-order differential equation (DE) (cf. eq. (49)) with $(r_{x},r_{y})=(0,0)$ being a critical point. Let us modify the right-hand side of eq. (49) by taking the traceless part of $m$, and denote its solution vector with $r^{\Psi}$:

[TABLE]

with $\Psi$ defined in eq. (3). Using equations (44) and (47), we can write the 3D extension of the above equation in the bent basis as:

[TABLE]

with $\Psi^{(\tau_{0})}$ being the traceless-transverse part of $m^{(\tau_{0})}$. Above we used the fact that $\mathrm{tr}(P_{\perp})=2$. Note that $\Psi^{(\tau_{0})}$ is then simply the traceless-transverse part of $M$.

We can check that \hat{\mbox{\boldmath{r}}}^{\Psi} rotates at the same rate as \hat{\mbox{\boldmath{r}}} in the normal plane if they start off collinear. To do that let us initialize the two vectors in 3D such that \hat{\mbox{\boldmath{r}}}^{\Psi{(\tau_{0})}}(\tau_{0})=\hat{\mbox{\boldmath{r}}}(\tau_{0}). The angle \hat{\mbox{\boldmath{r}}}^{\Psi} sweeps in the normal plane up to $\tau_{0}+\delta\tau$ is by construction equal to the angle its 3D extension $\hat{r}^{\Psi{(\tau_{0})}}$ sweeps, which in turn can be calculated as:

[TABLE]

The last line above equals $\delta\phi^{(\tau_{0})}(\tau_{0})$ (cf. eq. (55)), and therefore $\delta\phi^{(\tau_{0})}_{\Psi}$ indeed matches $\delta\phi^{(\tau_{0})}$.

Therefore, in order to study how $r$ is rotating in the normal plane, it is sufficient to understand the rotation of {\mbox{\boldmath{r}}}^{\Psi}. So, from now on let us focus on eq. (87) (and where necessary, its 3D extension, eq. (88)). As we are interested in rotations in the normal plane (and not in saddle-like behavior), let us focus on the case where $\Psi$ has a positive determinant. For constant $\Psi$, the DE in eq. (87) has a solution which is a center at {\mbox{\boldmath{r}}}^{\Psi}=0. Then, the phase portrait of {\mbox{\boldmath{r}}}^{\Psi} is a set of concentric ellipses. Let the semi-major axis of one of them is given by $cS_{a}$, and the semi-minor axis – by $cS_{i}$ (with $c$ being a constant). We can then pick $c$ and its units such that $S_{i}S_{a}=\omega_{c}$, with $\omega_{c}$ being the angular rate of rotation of {\mbox{\boldmath{r}}}(\tau) which, in anticipation of recovering the results of Section 2.2.1, we will call the coiling rate. In the basis aligned with the semi-axes of the ellipse, the DE can then be written as:

[TABLE]

where the last equality allows us to read off the $2\times 2$ matrix $\Psi$ in this particular basis. As a check, let us write the general solution to the system above: $r_{x}(\tau)=cS_{a}\cos\left[\omega_{o}\tau+\varphi\right]$, $r_{y}(\tau)=cS_{i}\sin\left[\omega_{o}\tau+\varphi\right]$, with $c$ and $\varphi$ being constants of integration. This confirms that eq. (90) has as a solution the ellipse described above.

Note that the ellipse with semi-axes specified by $s_{a}$ and $s_{i}$ in eq. (65) is the cross-sectional ellipse formed by the intersection of a flux tube of neighboring field lines with a plane orthogonal to the reference field line. In contrast, the ellipse with semi-axes specified by $S_{a}$ and $S_{i}$ above is the ellipse which each one of those field lines would describe in the normal plane around the reference field line under the constant flow specified by eq. (90).

From eq. (90), we can write the following relationships between $S_{i}$, $S_{a}$, $\omega_{c}$ and the invariants of the matrix $\Psi$:

[TABLE]

where we used the fact that $\Psi$ is traceless.

Now, let us allow for a varying {\hat{\mbox{\boldmath{B}}}}(\tau). We can still define $\omega_{c}$, $S_{i}$ and $S_{a}$ at each instant through the equations above. We can then substitute the 3D extension $\Psi^{(\tau_{0})}$ for $\Psi$ from eq. (88) in eq. (91) after noting that $\Psi^{(\tau_{0})}$ has one zero eigenvalue (for the eigenvector equal to {\hat{\mbox{\boldmath{B}}}}). Thus, we can write the instantaneous coiling rate as:

[TABLE]

where we used (34). Note that this can be written as (cf. eq. (78)):

[TABLE]

with $\Im(x)$ giving the imaginary part of $x$. Thus, up to a sign, $\omega_{c}$ is the imaginary part of one of the two complex conjugate eigenvalues of $M$. The sign of $\omega_{c}$ is picked to be the same as the sign of $\alpha$ (or equivalently, the average $\langle\omega\rangle$): positive for right-handed twist, and negative for left-handed.

Following the same procedure we used to obtain eq. (92), using eq. (91) we can calculate the following combination for non-zero $\omega_{c}$:

[TABLE]

One can show that the right hand side above equals $\alpha^{2}$ from eq. (81). The easiest way to see that is to write the equation in a basis in which {\hat{\mbox{\boldmath{B}}}}(\tau) has components (0,0,1). Thus, we can conclude that:

[TABLE]

Therefore, the above equation is the generalization of eq. (16) for field lines with non-zero curvature. Similar to our results of Section 2.2.1, we again see that the generalized force-free parameter, which is the magnetic field equivalent of vorticity of planar fluid flows, is boosted relative to the actual coiling rate by geometric factors related to the finite squeeze deformation needed to produce an elliptical transverse flow of field lines (see Section 2.2.1 for further discussion).

Equations (86) and (95) lead to the following inequality:

[TABLE]

Thus, we can conclude that the magnitude of the generalized force-free parameter, $\alpha$, can be enhanced relative to twice the coiling rate, $2\omega_{c}$, for two different reasons. (1) When the orthogonal flow produced by $M$ is not a spiral or center, i.e. when it is a saddle or a node, $\alpha$ generally can still assign a sense of rotation to that flow as long as the transverse part of $M$ is not symmetric (cf. right-hand side of eq. (94)). (2) When the orthogonal flow caused by $\Psi$ is a center (i.e. $M$ has a pair of complex conjugate eigenvalues), but the center produces an elliptical flow pattern, $\alpha$ is still enhanced relative to $2\omega_{c}$. The reason is that for elliptical flows (esp. with large aspect ratios), field lines passing close to the reference field line can sweep large angles even when $\omega_{c}$ is low, thus enhancing the average $\langle\omega\rangle$, which results in an enhanced $\alpha$ (see also the discussion after eq. (86)).

5.3 The coiling number

The local twist number $\mathrm{N_{t}}$ is defined as the integral over a field line $L$ (see Berger & Prior, 2006, eq. (16)):

[TABLE]

where the last equality comes from eq. (86). However, the above average includes contributions from saddle-like behavior of the field-lines in the transverse direction to {\hat{\mbox{\boldmath{B}}}}. Thus, one can envision constraining the integral only to those parts of the reference field line for which the phase portrait of the transverse flow is a center/spiral (i.e. $\Im(\lambda_{+})\neq 0$). Thus, one can define an alternative local twist number:

[TABLE]

The difference between $\mathrm{N_{t}}$ and $\mathrm{N_{t,im}}$ can be used to quantify the saddle-flow contributions to the local twist number.

We define the coiling number as:

[TABLE]

The coiling number, $\mathrm{N_{c}}$, gives an unbiased estimate of the number of turns neighboring field lines make around one another. It is unbiased in the following strict sense: for steady-state transverse flows (implying constant $\omega_{c}$) and integer $\mathrm{N_{c}}$, $\mathrm{N_{c}}$ matches the actual number of turns neighboring field lines make around one another (see Sections 2.2.1 and 5.2). In contrast, for elliptical transverse flows, the local twist number is biased by a factor depending to the aspect-ratio of the flow. The difference between $\mathrm{N_{c}}$ and $\mathrm{N_{t,im}}$ can be used to quantify that bias; see eq. (95).

6 Axially-symmetric, force-free flux ropes

In this Section, we are going to study the properties of generic, axially-symmetric, force-free flux ropes. We are first going to write down the general equations governing the magnetic field and current in such flux ropes, as well as the coiling and local twist rates. Then we are going to apply them to several specific flux-rope examples, including ones not discussed previously in the literature. We are going to collect our results in Fig. 6 along the way.

6.1 Generic, axially-symmetric, force-free flux ropes

Let us calculate the local coiling rate for flux ropes with cylindrical symmetry. In the orthonormal cylindrical basis (\mbox{\boldmath{\hat{r}}},\,\mbox{\boldmath{\hat{\phi}}},\,\mbox{\boldmath{\hat{z}}}), we impose the following conditions for the non-normalized magnetic field components:

[TABLE]

with non-zero $B_{\hat{\phi}}(r)$ and $B_{\hat{z}}(r)$, which are functions solely of the cylindrical radial coordinate $r$. In the above equation, $\alpha_{0}$ is a constant to be specified below. The limit in the second equality above is well-defined for non-diverging current distributions that close to the flux-rope axis are linear in $r$ (cf. eq. (6.1)). Given the above conditions, from equations (81), (93), as well as eq. (A.3) from Appendix A for the gradient of the magnetic field in cylindrical coordinates, we can express the coiling rate and the generalized force-free parameters as:

[TABLE]

From the first equation above, we can see that the coiling rate is defined for increasing $\hat{B}_{\phi}^{2}$, and vanishes when $\hat{B}_{\phi}$ has an extremum, or when $\hat{B}_{\phi}$ vanishes. The first condition is always satisfied near the axis of the flux ropes we consider with finite current distribution as can be seen from eq. (100). This behavior can be seen by comparing the right panels on the second and third rows of Fig. 6 for the sample of flux ropes discussed below.

Let us also write down the coiling number, local twist number, and the true non-local twist number, all evaluated per unit flux-rope length measured along the rope axis ($L_{z}$). We are going to denote those quantities by $n_{c}$, $n_{t}$ and $n_{\mathrm{true}}$, respectively. Using eq. (99), one can calculate $n_{c}$ as:

[TABLE]

where $\tau_{z}$ is the field-line length parameter along the axis of the flux rope, and the integrals are taken over a fixed flux-rope length, $L_{z}$. Similarly, we can write (see also Liu et al., 2016):

[TABLE]

One can eliminate the radial derivative of $\hat{B}_{\hat{\phi}}$ from eq. (6.1), and combine with equations (102) and (103) to write a concise relationship between the rates above:

[TABLE]

To obtain the equation above, we used the normalization condition:

[TABLE]

For generic flux ropes with cylindrical symmetry, one can Taylor expand $B_{\hat{\phi}}(r)$ and $B_{\hat{z}}(r)$ and then substitute those series in eq. (6.1). As long as the conditions in eq. (100) are satisfied, one can then explicitly check that:

[TABLE]

at the flux rope axis, which is confirmed in all examples in this section (compare the panels in the second row of Fig. 6 together with the left panel on the third row of the same figure). The last quantity is the true non-local twist rate (equal to $2\pi n_{\mathrm{true}}$; cf. eq. (103)) of the axially symmetric rope, and therefore, $\omega_{c}$ and $\alpha/2$ are equally good measures of twist at the core of such flux ropes141414The fact that $\alpha/2$ gives the correct coiling rate at the core of axially-symmetric flux ropes with smooth current density, was shown by Liu et al. (2016). In Section 6, we show that that requirement can be relaxed to include non-diverging current distributions, which close the flux-rope axis become linear functions of radius.. As one moves away from the axis of the rope, both $\alpha/2$ and $\omega_{c}$ diverge from the true non-local rate of twisting of field lines around the axis of the rope (see Fig. 6). One can quantify that deviation by using the $c_{3}$ parameter introduced by Liu et al. (2016). However, predicting that non-local twisting rate from purely local properties of the field is not possible for generic 3D flux ropes, which is why in this paper we focused only on predicting as accurately as possible the local rate of coiling of infinitesimally separated field lines.

To apply our analysis to simple flux ropes, first let us construct a generic, force-free, axially-symmetric flux rope using the normalized components of the magnetic field, $\hat{B}_{\hat{z}}(r)$ and $\hat{B}_{\hat{\phi}}(r)$. We do that by imposing the force-free condition \mbox{\boldmath{j}}=\alpha{\mbox{\boldmath{B}}}, where we use units for current, such that the current density is given by \mbox{\boldmath{j}}=\mbox{\boldmath{\nabla}}\times{\mbox{\boldmath{B}}}. The force-free condition implies that $j_{\hat{\phi}}/B_{\hat{\phi}}=j_{\hat{z}}/B_{\hat{z}}$, in which we can substitute the components of the normalized magnetic field, multiplied by a radially-dependent magnitude of the magnetic field $B(r)$, to find:

[TABLE]

where $B_{0}\equiv B(r=0)$ is a constant, and we applied eq. (105). From the above equation, we can see that $B(r)$ is a monotonically decreasing function (top-right panel of Fig. 6). The integral above converges as long as we require that $\hat{B}_{\hat{\phi}}$ decays at the rope axis sufficiently fast, which holds when the boundary condition, eq. (100), is satisfied.

One can use equations (105) and (107), combined with eq. (6.1), to show that the current density for a generic, axially-symmetric, force-free flux rope is given by:

[TABLE]

where a prime denotes a derivative in $r$.

If one is solving eq. (108) for $B(r)$, one needs to solve it for $B^{\prime\prime}(r)$ first. The equation has two roots, only one of them being physical. We write it down here for completeness:

[TABLE]

As can be seen from eq. (107), the quantity under the square root above is guaranteed to be non-negative by the requirement $0\leq\hat{B}_{\hat{\phi}}^{2}\leq 1$, which serves as a constraint on the allowed $j(r)$ profiles. Similarly, $\omega_{c}$ can be expressed solely using $B(r)$:

[TABLE]

Equations (6.1, 105, 107, 108, 109, 110), combined with the boundary condition, eq. (100), allow one to construct a generic, force-free, axially-symmetric flux rope, as long as only one of the functions $\hat{B}_{\hat{\phi}}(r)$, $\hat{B}_{\hat{z}}(r)$, $\alpha(r)$, $\omega_{c}(r)$, $j(r)$ or $B(r)$ is specified. Using those equations, one can explicitly check (by Taylor expansion) that the boundary condition, eq. (100), translates to the following behavior with $r\to 0^{+}$:

[TABLE]

where $\alpha_{0}\equiv\alpha(r=0)$. The above set of equations can be used to set the boundary condition at $r=0$ for equations (6.1, 105, 107, 108, 109, 110). So, for example, from eq. (6.1), we can see that at $r=0$ we have $B^{\prime}(r=0)=0$ and $\hat{B}_{\hat{z}}^{\prime}(r=0)=0$. Next, we are going to use those expressions to study specific examples of force-free, axially-symmetric flux-ropes151515 For completeness, we give the expressions for the curvature ($\kappa$) and torsion ($\tau_{r}$) of field lines in axially-symmetric flux ropes. Using the Frenet-Serret formulas from differential geometry, as well as eq. (A.3), we obtain:

$\displaystyle\kappa=\frac{\hat{B}_{\hat{\phi}}^{2}}{r}\ ,\ \ \tau_{r}=\frac{\hat{B}_{\hat{\phi}}\hat{B}_{\hat{z}}}{r}\ .$

(112) At the flux rope axis, those quantities approach:

$\displaystyle\kappa\approx\frac{\alpha_{0}^{2}r}{4}+\mathcal{O}(\alpha^{3}_{0}r^{2})\ ,\ \ \tau_{r}\approx\frac{\alpha_{0}}{2}+\mathcal{O}(\alpha_{0}r)\ ,$

(113) where we used eq. (6.1). Note that similar to the coiling rate, the torsion approaches the true non-local twist rate at the flux rope axis (cf. eq. (106))..

6.2 The Gold-Hoyle and Lundquist flux ropes

In the Gold-Hoyle model of a flux rope (GH; Gold & Hoyle, 1960), we have $B_{\hat{\phi}}=B_{0}br/(1+b^{2}r^{2})$ and $B_{\hat{z}}=B_{0}/(1+b^{2}r^{2})$ (with $b$ and $B_{0}$ being constant), thus satisfying eq. (100). The GH model has a constant non-local twist rate given by $n_{\mathrm{true}}(r)=b/(2\pi)$ (see eq. (103)). In Fig. 6, we show the results of Section 6.1 applied to the GH model.

Equation (6.1), applied to the GH model, gives:

[TABLE]

Comparing the first equality above with eq. (95), we can see that, independent of $r$, the transverse field-line flow for this model corresponds to a center with purely circular trajectories around it. Thus, for this flux rope model, $\alpha/2$ does not overestimate the local field-line coiling rate.

However, we can immediately point to a counterexample: the Lundquist model (Lundquist, 1950), which again satisfies eq. (100). In the Lundquist model, $B_{\hat{\phi}}=B_{0}J_{1}(\alpha r)$ and $B_{\hat{z}}=B_{0}J_{0}(\alpha r)$, where $B_{0}$ and $\alpha(r)=\alpha_{0}$ are constants, the latter matching the force-free parameter; and $J_{n}$ are the Bessel functions of the first kind. From eq. (6.1), we obtain:

[TABLE]

Clearly, the coiling rate approaches $\alpha/2$ at the flux rope core; however, away from the rope axis, the two diverge, indicating that the transverse flow becomes elliptical (cf. eq. (95)).

We demonstrate the properties of the Lundquist model in Fig. 6. One can see that the true non-local twist rate becomes infinite when $\hat{B}_{\hat{z}}(r)=0$, which happens at the zeros of $J_{0}(\alpha_{0}r)$: $\alpha_{0}r/2\approx 1.20$ and $2.76$, and so on (third row, right panel of Fig. 6). Beyond the first zero at $\alpha_{0}r/2\approx 1.20$, the coiling rate becomes imaginary for a range of radii, implying that locally the field lines are no longer twisting around one another, and instead they form a saddle flow in the normal plane. Thus, one can argue that the first zero of $\hat{B}_{\hat{z}}(r)$ delineates the innermost boundary of the Lundquist flux rope. Indeed, beyond that radius, the sense of winding of the field lines around the flux-rope axis reverses (third row, left panel of Fig. 6).

6.3 Further flux-rope examples

6.3.1 Flux rope with constant coiling rate

One may be curious whether there are closed-form expressions for axially-symmetric flux ropes with constant $\omega_{c}$. Combining eq. (6.1) with eq. (100) one can show that for $\omega_{c}=$constant, we have $\hat{B}_{\hat{\phi}}=\omega_{c}r$, with $\hat{B}_{\hat{z}}=\sqrt{1-\hat{B}_{\hat{\phi}}^{2}}$; and therefore, such a flux-rope is defined only within $r\leq|\omega_{c}|^{-1}$. Knowing the normalized components of the magnetic field, one can construct a force-free flux rope using eq. (107):

[TABLE]

where $B_{0}$ and $\omega_{c}$ are constants. Combining the above equation with equations (102) and (103), one can easily show that $n_{c}(r)=n_{\mathrm{true}}(r)$ for this flux rope, i.e. the coiling number per unit flux-rope length matches the true non-local twist number.

The corresponding $\alpha$ is given by:

[TABLE]

The current distribution for this flux rope diverges at $r\to|\omega_{c}|^{-1}$, has a non-zero minimum at $r=|\omega_{c}|^{-1}\sqrt{2-\sqrt{2}}\approx 0.765|\omega_{c}|^{-1}$ and a finite maximum at $r=0$. See Fig. 6 in which we demonstrate the properties of this flux rope with a constant $\omega_{c}(r)=\alpha_{0}/2$. This exhibits a core in which the current density is smaller than at flux-rope periphery. Indeed, one is free to use the equations of Section 6.1 to construct hollow-core flux ropes as well (see Fig. 6, where we include such an example); see below for a discussion.

6.3.2 Flux rope with constant coiling number per unit flux-rope length

Another analytic example is a flux rope with a constant coiling number per unit flux-rope length: $n_{c}=$const., with $n_{c}$ given by eq. (102). Combining that with eq. (6.1), one can show that in this case $\hat{B}_{\hat{z}}(r)=\exp\left[-(2\pi rn_{c})^{2}/2\right]$, with $\hat{B}_{\hat{\phi}}=\pm\sqrt{1-\hat{B}_{\hat{z}}^{2}}$. Once again, one can construct a force-free flux rope by a suitable choice of $B(r)$. That flux-rope is given by:

[TABLE]

with $B(r)$ given by eq. (107) with $\hat{B}_{\hat{\phi}}=B_{\hat{\phi}}/B$, and $B_{\hat{\phi}}$ given by eq. (118). The corresponding current distribution is peaked around the $r=0$ axis, similar to the majority of the models shown in Fig. 6. For completeness, we give the corresponding analytical coiling rate and $\alpha$ below:

[TABLE]

6.4 Discussion

One can similarly construct a flux rope with a constant local twist number per unit flux-rope length ($n_{t}(r)$=const.), as well as a flux rope for which the local twist number per unit flux-rope length matches the true non-local twist number ($n_{t}(r)=n_{\mathrm{true}}(r)$). We do not provide analytical expressions for those two flux-rope configurations. Instead, in Fig. 6 we show our results for them by numerically integrating the equations derived in Section 6.1.

From Fig. 6, one can see that the flux ropes with constant coiling, local twist, and non-local twist rates per flux-rope length (i.e. those in the figure denoted with $n_{c}(r)=\alpha_{0}/(4\pi)$, $n_{t}(r)=\alpha_{0}/(4\pi)$, and GH, respectively), and especially the first two, show quite similar properties in terms of current density, magnetic field, $\alpha$, and $\omega_{c}$. However, the true non-local twist rate becomes drastically different beyond $\alpha_{0}r/2\sim 1$: compare the left panel on the third row with both panels on the last row of Fig. 6.

From this set of examples, one can see that as long as $\alpha_{0}r/2\lesssim 1$, the non-local twist number is well approximated by both the coiling and local twist numbers to within several tens of percent. This is not surprising in light of equations (106) and (6.1). This statement is valid independent of whether we are looking at hollow-core or centrally-peaked current distributions. Indeed, in Fig. 6 we include a hollow-core flux rope example (with a current density given by the more or less random $j(r)=\alpha_{0}B_{0}\exp\left[(r\alpha_{0}/2)^{4}-(r\alpha_{0}/2)^{40}\right]$, which produces such a current profile).

However, the numerical results quoted above show that any study relying on determining an accurate non-local twist number far from the flux-rope axis, cannot rely solely on the local twist or coiling numbers. Instead, one would have to resort to finding the flux-rope axis to calculate the non-local twist rate (e.g. Liu et al., 2016). In turn, for complex magnetic field configurations, especially those found in global models of the solar corona, the first step to finding the flux-rope axis is finding the flux ropes themselves. As argued in Section 5 and as shown in a follow-up paper (Savcheva, 2019), using the natural $\mathrm{N_{c}}\gtrsim 1$ threshold as a flux-rope detection threshold is superior to using $\mathrm{N_{t}}\gtrsim 1$ due to the latter giving many false-positives. As the results discussed above (and shown in Fig. 6) indicate, once the presence of a flux rope has already been established, both the coiling and local twist number are equally reliable estimates of the non-local twist sufficiently close to the flux rope axis.

7 Summary

We study the properties of the local transverse deviations of magnetic field lines. Those deviation vectors follow a planar flow in a plane normal to one of the field lines (the “reference” field line, which can be any field line of interest: e.g. the one passing through the voxel for which we are calculating the coiling and squeezing rates) as the plane is moved along the length of that field line. As long as we focus only on field lines with infinitesimal separations, the resulting planar field-line flow is linear with a critical point at the location where the reference field line intersects the normal plane. Therefore, the types of transverse field-line flows are given by the standard equilibrium solutions of linear systems, such as a saddle, cycle, spiral or node. The type of solution is unique for the normal plane that is minimally rotated around its normal when moved along the reference field line (see Section 3.3). We show that that solution is determined by the two non-zero eigenvalues of the gradient of the normalized magnetic field.

We argue that the eigenvalue difference, $\Delta\lambda$, can be used to quantify the degree of squeezing or coiling of neighboring field lines. Thus, when $\Delta\lambda$ is real, we define it as the squeezing rate ($\rho_{\mathcal{Z}}$) which, when integrated over a field line, gives $\ln(\mathcal{Z})$, with $\mathcal{Z}$ defined as the squeeze factor. We demonstrate that $\mathcal{Z}$ can be approximated by the squashing factor, $Q$. Thus, the squeeze factor offers an alternative, with which one can check the robustness of a QSL detection.

QSL maps have become largely synonymous with maps161616However, calculating maps of the squashing factor are computationally expensive, which lead Tassev & Savcheva (2017) to propose an alternative, computationally cheap, crude method for locating QSLs using so-called FLEDGE maps. of the squashing factor, $Q$, with QSLs identified as the locations of exponentially large $Q$ (e.g. Pariat & Démoulin, 2012). Yet, $Q$ has the following drawback: it can miss QSLs corresponding to regions of large local field-line squeezing, when that squeezing is undone before the field lines reach the footpoints between which $Q$ is calculated. As the squashing factor is an inherently globally-defined quantity, we show that one cannot construct an exact local squashing factor rate, which could account for such localized squeezing. However, we show that $\rho_{\mathcal{Z}}$ can be treated as a local approximation to the rate of squashing in a certain strict sense, described in the paper. As by definition, $\rho_{\mathcal{Z}}\geq 0$, the squeeze factor, $\mathcal{Z}$, does not suffer from the problem of missing unsqueezed regions described above. Thus, in that regard, $\mathcal{Z}$ is a better proxy for QSLs than the squashing factor. Moreover, given by just an integral over a field line, the squeeze factor is computationally cheaper than the squashing factor, which in addition, depends on computing the relative deviations of neighboring field lines (e.g. Tassev & Savcheva, 2017).

As $\mathcal{Z}$ grows exponentially large in the vicinity of HFTs and null points, one can use a threshold in the local $\rho_{\mathcal{Z}}$ to find the locations of nulls, HFTs and current sheets (see Savcheva (2019)). This is in contrast to the squashing factor, for which an equivalent, purely local squashing rate does not exist; and thus, one has to run a full non-local squashing factor calculation to obtain a similar HFT detection (e.g. Savcheva et al., 2016). One may argue that when calculated on a grid, $\rho_{\mathcal{Z}}$, lacks discriminatory power to distinguish a null point from a short HFT, for example. However, finding a null point from a discrete sampling of the magnetic field is inherently ambiguous as it depends on the interpolation scheme used to supersample the magnetic field inside each grid cell171717Whether a null is present inside a grid cell depends on whether the three surfaces defined by $B_{i}=0$ intersect at a point within the cell or not. As the shape of those surfaces is determined by the interpolation scheme used to interpolate $B$ within the cell, so does the presence of a null point. For an algorithm for finding null points using trilinear interpolation, see (Haynes & Parnell, 2007)..

When $\Delta\lambda$ is imaginary, we define it as twice the coiling rate ($\omega_{c}$) which, when integrated over a field line, gives the coiling number, $\mathrm{N_{c}}$. We show that unlike the standard local twist number (e.g. Berger & Prior, 2006, eq. 16), the coiling number gives an unbiased estimate of the local number of twists neighboring field lines make around one another. The local twist number is defined as a field-line integral over what-we-call the generalized force-free parameter ($\alpha$), which is proportional to the component of the current parallel to the magnetic field. As $\alpha$ is the equivalent of vorticity of the planar transverse field-line flow, it is non-zero for saddle flows which exhibit shear. Yet, field lines do not revolve around one another in such flows. Therefore, $\alpha$ clearly gives a biased estimate of local twist, resulting in $\mathrm{N_{t}}\geq 1$ even in regions exhibiting shear, well outside of flux ropes. We find that $\alpha$ is biased even for elliptical/spiral transverse field-line flows. We offer a geometric interpretation of that bias: $\alpha$ is enhanced relative to the average rate of rotation of neighboring field lines with factors depending on the aspect ratio of the transverse field-line flow. We show that $\omega_{c}$ is unbiased in that regard, as for steady-state flows, $2\pi/\omega_{c}$ exactly equals the period of rotation of infinitesimally-separated field lines around one another. We also show that for non-diverging current distribution, $\omega_{c}$ approaches $\alpha/2$ at the core of axially-symmetric flux ropes, and therefore, is an equally good measure of twist at the core of such flux ropes (Section 6). The coiling number can be useful in the study of flux rope instabilities, such as the kink instability (e.g. Liu et al., 2016). Moreover, one can use the natural inequality $\mathrm{N_{c}}\gtrsim 1$ for detecting flux ropes in complicated magnetic field configurations. We demonstrate that in (Savcheva, 2019).

We give closed-form expressions for $\rho_{\mathcal{Z}}$ and $\omega_{c}$ in terms of the gradient of the normalized magnetic field in both Cartesian and curvilinear coordinates. We implemented the calculation of the squeeze factor and coiling number, as well as their respective local rates, in the public code181818The code is available here: https://bitbucket.org/tassev/qsl_squasher/. The code allows one to output the squeeze rate as calculated using either of the two options proposed for handling elliptical flows: equations (27) or (28); or using the squeeze rate corresponding to the symmetrized $m$ (Section 4.5). The code also outputs the local type of critical point (e.g. spiral, node, saddle) – thus, offering another type of partitioning of magnetic fields – as well as the trace of the gradient of the normalized magnetic field. The code can be easily modified to output the aspect ratio of elliptical flows using eq. (95); the orientation of the eigenvectors of the flow (or for elliptical flows – the orientation of the semi-axes; see footnote 3); as well as non-transverse quantities such as the curvature of the field lines using eq. (50). We are going to explore those field-line properties elsewhere. QSL Squasher introduced by Tassev & Savcheva (2017). In a follow-up paper, we show numerical calculations of those quantities for realistic magnetic fields (Savcheva, 2019).

We would like to thank Edward DeLuca for useful discussions and comments on the manuscript.

Appendix A Results for curvilinear coordinates

In this section, we obtain the expressions for the coiling and squeezing rates in curvilinear coordinates. Below we first write our results for any coordinate basis, and then we focus on spherical and cylindrical coordinates. We also address an approximation introduced in the squashing factor calculation in spherical coordinates in (Tassev & Savcheva, 2017).

A.1 The field-line deviation equation

For two infinitesimally separated points in space at positions $x$ and {\mbox{\boldmath{x}}}+d{\mbox{\boldmath{x}}}, the difference between their radius-vectors can be written as d{\mbox{\boldmath{x}}}=dx^{i}\mathbf{e}_{i}({\mbox{\boldmath{x}}}) only in coordinate basis \mathbf{e}_{i}({\mbox{\boldmath{x}}}) (indeed, this can serve as a definition of a coordinate basis). The equation for the field-line radius vector {\mbox{\boldmath{x}}}(\tau) (eq. (29)) can then be written in any coordinate basis as:

[TABLE]

and similarly for a neighboring field line with radius vector {\mbox{\boldmath{y}}}(\tau).

The field line deviation equation (eq. (30)) for \delta{\mbox{\boldmath{x}}}(\tau)={\mbox{\boldmath{y}}}(\tau)-{\mbox{\boldmath{x}}}(\tau) then takes the form:

[TABLE]

where we linearized in \delta{\mbox{\boldmath{x}}} and used the fact that the Christoffel symbols, $\Gamma^{k}_{\ ij}$, give the spatial derivative of the basis vectors: $\partial_{j}\mathbf{e}_{i}=\Gamma^{k}_{\ ij}\mathbf{e}_{k}$. In writing the above equation, we used the standard coma notation for partial derivatives, and the semicolon notation for covariant derivatives.

The equality between the first and the last lines in eq. (A2) is hardly surprising. From that equation, one can read off the $\tau$ derivative of the components of \delta{\mbox{\boldmath{x}}} as $\hat{B}^{i}_{\ ;k}\delta x^{k}$. That is correct if the components of both \delta{\mbox{\boldmath{x}}}(\tau+\delta\tau) and \delta{\mbox{\boldmath{x}}}(\tau) (entering in d\delta{\mbox{\boldmath{x}}}(\tau) in the finite difference sense, with $\delta\tau$ being infinitesimal) are written in the common basis \mathbf{e}_{i}({\mbox{\boldmath{x}}}(\tau)). This is exactly what we want for the generalization of $M_{ij}$ in eq. (31) to curvilinear coordinates. Thus, we can conclude that in any coordinate basis, \mathbf{e}_{i}({\mbox{\boldmath{x}}}(\tau)), (being careful about index placement):

[TABLE]

which is a rather unsurprising result (but see Section A.2). We used that result to promote the derivatives, entering in the squeezing and coiling rates (see equations (78) and (93)), to covariant derivatives. Note that those rates are scalars, and therefore the covariant derivative entering in them can be written in non-coordinate basis as well.

In Section A.3, we write the above equation explicitly in the spherical and cylindrical orthonormal bases. However, let us first clear up a possible confusion, which may arise from eq. (A3), about how the standard squashing factor in curvilinear coordinates is calculated.

A.2 Calculating the standard squashing factor

Let us take a brief detour and write the equations for the field line deviation components in a form suitable for calculating the standard (non-local) squashing factor in spherical coordinates, as done in eq. (6) of (Tassev & Savcheva, 2017). To obtain eq. (A3) for $M^{i}_{\ j}$, we had to write all vectors and operators entering in eq. (A2) in the same common basis \mathbf{e}_{i}({\mbox{\boldmath{x}}}(\tau)), including both \delta{\mbox{\boldmath{x}}}(\tau+\delta\tau) and \delta{\mbox{\boldmath{x}}}(\tau), which enter in d\delta{\mbox{\boldmath{x}}}(\tau) in the finite difference sense (with $\delta\tau$ being infinitesimal). However, in the standard calculation of $Q$, we need to evaluate the components of \delta{\mbox{\boldmath{x}}}(\tau+\delta\tau) in its corresponding basis \mathbf{e}_{i}({\mbox{\boldmath{x}}}(\tau+\delta\tau)), and not in the basis \mathbf{e}_{i}({\mbox{\boldmath{x}}}(\tau)). The reason is that one then wants to project out the components of \delta{\mbox{\boldmath{x}}}(\tau+\delta\tau) that are parallel to {\hat{\mbox{\boldmath{B}}}}(\tau+\delta\tau), which in numerical codes is given in the basis \mathbf{e}_{i}({\mbox{\boldmath{x}}}(\tau+\delta\tau)) and not in the basis \mathbf{e}_{i}({\mbox{\boldmath{x}}}(\tau)). Thus, we write the left-hand side of eq. (A2) as:

[TABLE]

where we used eq. (A1) and linearized in $\delta\tau$. Note again (first line above), that $\delta x^{i}(\tau+\delta\tau)$ and $\delta x^{i}(\tau)$ are written in different coordinate bases.

Setting equations (A2) and (A4) equal, we can see that the terms containing the Christoffel symbols cancel out. To do that, we note that the Christoffel symbols are symmetric in their lower two indices, which is true for any coordinate basis. The end result is:

[TABLE]

with the partial derivative entering above, and not the covariant derivative which enters in eq. (A3). The above equation is valid in any coordinate basis, and therefore can be used to write the field line deviation equation in cylindrical or other coordinates, as needed. Below we focus our attention to spherical coordinates.

We use the right-handed coordinate basis \mbox{\boldmath{e}}_{\phi},\ \mbox{\boldmath{e}}_{\theta},\ \mathbf{e}_{r}, which is related to the right-handed non-coordinate orthonormal basis \hat{\mbox{\boldmath{\phi}}},\ \hat{\mbox{\boldmath{\theta}}},\ \hat{\mbox{\boldmath{r}}} through:

[TABLE]

Note that we use the convention used in (Tassev & Savcheva, 2017), where $\theta$ stands for latitude, and not co-latitude, so one needs to be careful when comparing with results in the right-handed orthonormal basis [radius, co-latitude, longitude].

In numerical datasets, the magnetic field in spherical coordinates is usually given in the non-coordinate orthonormal basis given above. We denote the components of {\hat{\mbox{\boldmath{B}}}} in that basis with hats over the indices, i.e. $\hat{B}^{\hat{i}}=[\hat{B}^{\hat{\phi}},\hat{B}^{\hat{\theta}},\hat{B}^{\hat{r}}]$. The components of {\hat{\mbox{\boldmath{B}}}} in the coordinate spherical basis, we denote by dropping the hats over the indices: $\hat{B}^{i}=[\hat{B}^{\phi},\hat{B}^{\theta},\hat{B}^{r}]$. Since {\hat{\mbox{\boldmath{B}}}}=\hat{B}^{i}\mathbf{e}_{i} holds in any basis, then using eq. (A6) we can write:

[TABLE]

Noting that the components of \delta{\mbox{\boldmath{x}}} in the coordinate spherical basis are given by $\delta x^{i}=[\delta\phi,\delta\theta,\delta r]$, and using eq. (A7) to write eq. (A5) using the components of {\hat{\mbox{\boldmath{B}}}} in the non-coordinate spherical orthonormal basis, after a bit of algebra we find:

[TABLE]

with the components of \delta{\mbox{\boldmath{x}}} in the non-coordinate orthonormal spherical basis given by:

[TABLE]

similar to eq. (A7) above, as well as to eq. (8) of (Tassev & Savcheva, 2017).

Comparing eq. (A.2) above with eq. (6) of (Tassev & Savcheva, 2017), we can see that the terms neglected in eq. (6) of that paper correspond to the terms $\left[\delta r/r-\delta\theta\tan(\theta)\right]B^{\hat{\phi}}$ and $(\delta r/r)\hat{B}^{\hat{\theta}}$ which are correspondingly subtracted on the right-hand sides of the first two equations in (A.2) above. In the above equations, those terms are suppressed relative to the rest by the ratio of the typical scale over which {\hat{\mbox{\boldmath{B}}}} varies and the radius of the Sun. We implemented those terms in the open-source QSL Squasher code, introduced by Tassev & Savcheva (2017). After performing numerical experiments with realistic data, we find that the additional terms affect the calculation of $Q$ only negligibly, given the overall qualitative nature of $Q$. Yet, we recommend one to implement those terms as they have the potential to be important in spherical coordinates for magnetic fields with small gradients.

A.3 The magnetic field gradient

The local squeezing ($\rho_{\mathcal{Z}}$) and coiling rates ($w_{0}$) are given by equations (78) and (93) which, as argued in Section A.1, are valid in any basis. Here, we would like to write those equations explicitly in spherical and cylindrical coordinates.

To implement those equations numerically, it is convenient to write the covariant derivative of {\hat{\mbox{\boldmath{B}}}} in the non-coordinate orthonormal basis given by eq. (A6). Calculating that derivative is a standard exercise, but we include it here for the sake of completeness. To do that, we need to start with the metric for spherical coordinates and calculate the Christoffel symbols, which allow us to find the covariant derivative $\hat{B}^{i}_{\ ;j}$ in the coordinate spherical basis. Then we substitute in the components of {\hat{\mbox{\boldmath{B}}}} in the orthonormal basis (eq. (A7)). Finally, we write the tensor components in the non-coordinate orthonormal basis (eq. (A6)). We perform those steps with the open-source mathematics software, SageMath191919https://www.sagemath.org/. The final result is:

[TABLE]

Given that the above expressions are written in an orthonormal basis without using the fact that {\hat{\mbox{\boldmath{B}}}} was normalized, the first three lines add up to the divergence (as their sum is simply the trace $\hat{B}^{\hat{i}}_{\ ;\hat{i}}$) of any vector field {\hat{\mbox{\boldmath{B}}}} expressed in spherical coordinates; line 4 minus line 5 (one of the three non-zero antisymmetric parts of $\hat{B}^{\hat{i}}_{\ ;\hat{j}}$) gives the longitudinal component of the curl of that vector field; line 6 minus line 7 gives the latitudinal component of the curl; line 8 minus line 9 gives the radial component of the curl. It is straightforward to check that these results match standard textbook results for the divergence and curl of a vector field in spherical coordinates with a few differences due to the fact that our $\theta$ stands for latitude, and not the usual co-latitude. Those differences are as follows: in the right-handed non-coordinate basis [radius, co-latitude, longitude], if $\theta$ above is to stand for the standard co-latitude, then all one has to do to find $\hat{B}^{\hat{i}}_{\ ;\hat{j}}$ in that basis is to replace $\cos(\theta)$ with $\sin(\theta)$, as well as $\tan(\theta)$ with $(-\cot(\theta))$, in eq. (A.3).

One can follow the same procedure we outlined above to express the covariant derivative of {\hat{\mbox{\boldmath{B}}}} in the cylindrical orthonormal basis (\hat{\mbox{\boldmath{r}}},\hat{\mbox{\boldmath{\phi}}},\hat{\mbox{\boldmath{z}}}), with $r$ being the cylindrical radial coordinate:

[TABLE]

Since we are using an orthonormal basis in writing out equations (A.3) and (A.3), index placement is unimportant, and the generalized force-free parameter $\alpha$ (eq. (81)) can be written as:

[TABLE]

Equations (78), (93), (A.3) and (A12) give the final result in spherical coordinates for the local squeezing and coiling rates. These equations were implemented in the public code QSL Squasher and were used in the calculations presented in the follow-up paper (Savcheva, 2019).