Coiling and Squeezing

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

Svetlin Tassev Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA Braintree High School, 128 Town Street, Braintree, MA 02184, USA Antonia Savcheva Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA
###### 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, . We show that unlike the standard local twist number, gives an unbiased approximation to the number of twists neighboring field lines make around one another. 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.

Sun:magnetic fields — Sun: magnetic topology
\correspondingauthor

Svetlin Tassev

## 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, (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 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 to denote only the local approximation to the twist, calling the non-local twist ; 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., , 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, 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, , as a robust flux-rope detection threshold. That threshold would seem to imply that nearby field lines exhibit at least 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 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, . In (Savcheva, 2019), we show that the natural inequality can be used successfully as a detection threshold for flux ropes. For each field line, is given as an integral of a local coiling rate () over the field line length.

We also introduce a quantity, complementary to the coiling rate, which we dub the squeezing rate (), measuring the local logarithmic rate of squeeze of neighboring field lines. The field-line integral over gives , with defined as the squeeze factor. We show that in a sense made precise in the paper; and therefore, gives an estimate of the local contributions to – 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 and 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, , 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 gives (Section 4), while the imaginary part equals (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 () given below as a function of position, :

 vi(\boldmath{r}(τ))=dri(τ)dτ=mijrj(τ) . (1)

The constant matrix above equals the velocity gradient, . Clearly is a critical point, which depending on the eigenvalues of , 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 () between two neighboring streamlines as a function of time (). 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 ) as a function of the field line length parameter . In that case, we show that one must replace with the transverse gradient of the normalized magnetic field, (cf. eq. (49)).

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

 m=Θ+Ψ=Θ+Σ+Ω , (2)

with:

 Θ≡Itr(m)2,Ψ≡m−Θ,Σ≡Ψ+ΨT2,Ω≡m−mT2=Ψ−ΨT2 , (3)

where is the identity matrix. Clearly, is the diagonal, trace part of , usually called the rate-of-expansion tensor. Then, is the traceless part of . The latter is in turn usually decomposed into a symmetric traceless part (the rate-of-shear tensor; ) and an antisymmetric part (the rate-of-rotation tensor; ). 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 , that renders the critical point of the remaining traceless part of the velocity gradient () into either a saddle (for real eigenvalues of ) or a center (for complex eigenvalues of ).

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

 ri(τ+δτ)=(I+δτm)ijrj(τ) . (4)

The infinitesimal transformation entering above, can be an infinitesimal squeeze (), shear (), scaling () or rotation () 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):

 δZ(z) = (1+zδτ001−zδτ), δS(s)=(1sδτ01), δC(c) = (5)

Those four transformations depend on four parameters: , , and , 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 as , 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, , of the velocity gradient. Note that depends on three degrees of freedom since has four elements, but we imposed the condition . One of those degrees of freedom amounts to a trivial choice of orthonormal basis, which is defined up to an overall rotation. Thus, has only two non-trivial degrees of freedom, which means that it can be decomposed in various ways using the infinitesimal maps , and , given in eq. (2.1). Let us explore those decompositions below.

Since the rate of shear, , is symmetric, it has two orthogonal eigenvectors, and because it is traceless, its eigenvalues are opposite in sign. Thus, can be treated as an infinitesimal squeeze map, , along the eigenvectors of . As is antisymmetric, is an infinitesimal rotation map, , 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 ) the decomposition , can be written as:

 I+Ψδτ≈(I+Σδτ)(I+Ωδτ)=δZδR≈δRδZ . (6)

As eq. (4) is linear in , we kept only terms to linear order in 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 . This can be seen by using from eq. (6) in combining equations (1) and (2). Thus, in the deformation map language, the standard decomposition of 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 : a saddle and a center.

As noted earlier, one can use other (non-standard) decompositions of , however. For saddles, one can write 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 . That decomposition can be achieved by choosing the (repeated) eigenvectors of to coincide with one of the eigenvectors of as well as with one of the eigenvectors of . That shear-squeeze decomposition is illustrated in the third row of Fig. 3. For centers, one can decompose as , 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 of the semi-major axis of elliptical flow around a center generated by a matrix is not readily available in the literature. Thus, we include that here for completeness. First, one needs to find the complex conjugate eigenvectors of . Then, after a bit of algebra, one can show that the following relationship holds between the and components of and : where , and is the (fixed) aspect ratio (the ratio of the semi-minor () to semi-major () 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 , 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, , where is the velocity field of the flow. A planar flow (in a plane spanned by the basis , ) has a vorticity vector which is normal (along the third basis vector ) to the flow and is given by

 ωz=∇xvy−∇yvx=myx−mxy , (7)

which has a magnitude equal to and is, therefore, a rotation invariant. Thus, the vorticity equals twice the rotational rate (; 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 . 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 ; and thus, the number of turns they make around is always less than one, despite the fact that the imaginary vane makes turns for a time interval :

 Nt=12π(ωz2Δτ) . (8)

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 (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 over .

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:

 Nt=12π∫Ldτα(τ)2 . (9)

Here, the integral is over a field line ; while is the generalization of the force-free parameter (which we define for any magnetic field as ). For magnetic fields, plays the role of the fluid vorticity above (compare eq. (7) with eq. (84)). In Section 5, we show that 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. 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 (), shear (), scaling () and rotation () maps, obtained by composing (infinitely many of) their corresponding infinitesimal maps (eq. (2.1)):

 Z(zf) = (zf001/zf), S(sf)=(1sf01), C(cf) = (cf00cf), R(θ)=(cos(θ)−sin(θ)sin(θ)cos(θ)) , (10)

with the finite map factors , , , defined through the equations above. Any rate of expansion () is uniquely captured by a scale map () 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, , of the velocity gradient.

In this section, we decompose 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:

 I+Ψδτ=RAδBA−1R−1 , (11)

with and 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 . The rotation maps do not affect the invariants of , 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 .

Note that the structure of the decomposition in eq. (11) is chosen to be such that the flow, , generated by the traceless part of (the constant) in eq. (4) can be solved to give the position of a fluid element after a finite interval (with ) as:

 ~ri(τ+nδτ) = (I+δτΨ)ij~rj(τ+(n−1)δτ)=(I+δτΨ)nij~rj(τ) (12) = (RAδBA−1R−1)nij~rj(τ)=(RA(δB)nA−1R−1)ij~rj(τ) = (RABA−1R−1)ij~rj(τ) .

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

Transformation can be either a shear or a squeeze (as any rotation can be absorbed into ). Therefore, the non-trivial combinations of the pairs are: (squeeze, shear), (shear, squeeze), (shear, rotation), (squeeze, rotation), where we eliminated repeated pairs as the composition of is still a transformation of the same type, , 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 , 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 ()=(squeeze, rotation) in eq. (11). The rotation matrices in that equation can be used to rotate the basis vectors by an angle so that they are aligned with the semi-axes of the elliptical flow around the center generated by . Then can be used to un-squeeze the elliptical flow and render it into a uniform circular flow, rotating with a constant angular frequency , which we will call the coiling rate. At each time-step (), that rotation is generated by . 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 (where and are the semi-minor and semi-major axes of one of those ellipses, respectively), the squeeze factor must equal as (eq. (2.2)) rescales the component of by and by . Thus, we can write:

 I+Ψδτ=R(θ)Z(√SaSi)δR(ωc)Z−1(√SaSi)R−1(θ) . (13)

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 (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 over the field line length, which in (Savcheva, 2019) we set at (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 in this decomposition are . As does not affect the difference of the eigenvalues of , we can write that the difference between the eigenvalues () of equals that of :

 λ+−λ−=2iωc , (14)

and therefore, the period with which fluid streamlines revolve around equals . 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 matches that of , which in our case is the composition of infinitesimal maps. Thus, we can define a coiling number for elliptical flows as:

 Nc≡12πωcΔτ=14πI{λ+−λ−}Δτ , (15)

with the coiling rate, , 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 . 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 in a finite depends on the initial position of the fluid element as well as on (cf. Section 5, as well as (Berger & Prior, 2006)). However, since the flow is periodic, when is an integer, all fluid elements would have made exactly turns around the origin, independent of their initial positions. Another way of thinking about is that the un-squeezed flow (middle step in the decomposition shown in the second row of Fig. 4) is uniformly rotating with ; 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)):

 ωz=myx−mxy=Ψyx−Ψxy=(SiSa+SaSi)ωc , implying: |ωz|≥2|ωc| , (16)

with the sign of chosen to match that of . 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, , 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 . At the very least, as discussed above, when is an integer, it matches the actual number of turns the fluid elements make, while 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, is the magnetic equivalent of (see the discussion around eq. (9)). In Section 5, we find that we can write eq. (16) for as (cf. eq. (95)):

 α=(SiSa+SaSi)ωc , (17)

with 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:

 Nc≡12π∫Ldτ ωc(τ) . (18)

Similar to the coiling number in the fluid context, 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 when 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 ()=(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 so that one of the basis vectors is aligned with one of the eigenvectors of . Then can be used to un-shear the flow and render the second eigenvector of orthogonal to the first. The flow would then be a saddle flow with orthogonal eigenvectors which, at each time-step (), can be generated by an infinitesimal squeeze map with a squeeze rate , where the factor of is introduced for convenience. After the infinitesimal squeeze, we need to shear back the flow with so that the directions of its asymptotes again match the eigenvectors of . If we denote the angle between the eigenvectors of as (see top-left panel of Fig. 4), then the required shear factor needs to equal , which one can check explicitly using eq. (2.2). Thus, we can write this decomposition as:

 (19)

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 in this decomposition are . As does not affect the difference in eigenvalues of , we can conclude that the difference between the eigenvalues of equals that of . Therefore, the (relative logarithmic) squeeze rate between the directions corresponding to the eigenvectors of , and hence (as does not affect the eigenvectors of either), is given by . Indeed, if are the two unit eigenvectors of corresponding to , then for two fluid parcels lying on the asymptotes of the flow (in the eigenvector directions relative to the origin), we can write and . Then, from eq. (1), one can check explicitly that , and therefore, we can summarize:

 ddτln(r+(τ)r−(τ))=λ+−λ−=ρZ . (20)

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

 ωz=myx−mxy=Ψyx−Ψxy=ρZcot(ϑ) . (21)

Thus, we can see that similar to the analysis of elliptical flows in the previous section, the vorticity (and hence, the force-free parameter for magnetic fields) for saddle flows has a clear dependence on the geometry of the flow. Through its dependence on the angle () 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 , 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 () (Titov, 2007) to quantify the divergence of neighboring field lines, separated by in the plane normal to one of them. Let us see how relates to . We review the calculation of 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 relating evaluated at two different locations along a field line: at and for some finite :

 ri(τ+Δτ)=Jijri(τ)=(I+δτm)Δτ/δτijri(τ) , (22)

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

 (23)

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:

 Q=2[csc2(ϑ)cosh(ρZΔτ)−cot2(ϑ)] , (24)

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

We are interested in large , which holds555In the limit , has repeated eigenvectors, and despite the eigenvalues being the same (implying ), can still grow large, as can be seen from eq. (24). However, one can check that in that case grows only as , and not exponentially, and therefore we disregard those special cases. for large real arguments of , in which case the logarithm of approaches . Therefore, we can approximate as an integral over the local squeezing rate, given by in eq. (20). We denote that approximation of with , 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)):

 dlnQdτ∼dlnZdτ=ρZ=λ+−λ−for real eigenvalues of m, (25)

with the squeezing rate, , defined as the squeezing rate of the flow in the decomposition666Note that 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 . 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 . 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, remains the same. This can be seen from eq. (24), which is valid even for complex eigenvalues of , in which case (compare equations (14) and (20)). Note that can periodically still grow large, however. It is maximized at , and every later along the field-line, in which case one can show that equals , where we used eq. (16) and (21) to eliminate 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 ) 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 , within a quarter of a period, changes by:

 ∣∣∣ΔlnQΔτ∣∣∣ = 2|ωc|πln(Qmax.2)=2|ωc|πln[12(S2aS2i+S2iS2a)] . (26)

So, at this point we have two different possibilities () for the logarithmic squeezing rate under elliptical flows. If one is not interested in periodic localized flux-tube squashing, one can use:

 Option 1:dlnZ1dτ≡ρZ,1=0for complex eigenvalues of m (27)

since 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:

 Option 2:dlnZ2dτ≡ρZ,2=2|ωc|πln(α22ω2c−1)for complex eigenvalues of m. (28)

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

The ambiguity in the definition of 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 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)), carries complementary information to that of : the former would depend only on the real part of the eigenvalue difference of , 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 () 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 , 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 and with being the field-line length parameter. We pick an initial condition such that when both field lines intersect the plane normal to the first field line at . The magnetic field lines and are calculated as the integral curves of the unit magnetic field vectors, . Thus, in Cartesian coordinates, we have:

 d\boldmath{x}(τ)dτ=^% \boldmath{B}(\boldmath{x}(τ)) , (29)

and similarly for .

Then we can construct the field-line deviation vector . Using (29), the deviation vector is a solution to the linearized equation:

 dδ\boldmath{x}(τ)dτ=^% \boldmath{B}(\boldmath{y}(τ))−^% \boldmath{B}(\boldmath{x}(τ))≈(δ\boldmath{x}(τ)⋅\boldmath{∇}% \boldmath{x})^\boldmath{B}(\boldmath{% x}(τ)) . (30)

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

 dδxi(τ)dτ=(∇j^Bi)δxj≡Mijδxj , (31)

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

Let us additionally define the deviation vector , with chosen such that is orthogonal to the field line at :

 \boldmath{r}(τ)⋅^\boldmath{B}(\boldmath% {x}(τ))=0 .

Above we used the fact that is the tangent vector for the field line passing through . Choosing the field lines and infinitesimally close to each other, we can write with being small. Thus, we can linearize:

 \boldmath{r}(τ)≈\boldmath{y}(τ)−% \boldmath{x}(τ)+εd\boldmath{y}(τ)dτ=δ\boldmath{x}(τ)+ε^\boldmath{B}(\boldmath{y}(τ))≈δ\boldmath{x}(τ)+ε^\boldmath{B}(\boldmath{x}(τ)) ,

where we used eq. (29). Setting , we find , which in turn implies888 The same result was obtained by Scott et al. (2017). (in matrix notation):

 \boldmath{r}(τ)=P⊥δ\boldmath{x% }(τ) , (32)

where the projection operator is given in component form by: . From eq. (32) we can see that equals the projection of onto the plane normal to the reference field line at location .

We can write an evolution equation for similar to eq. (30) for . To do that, we need to find the derivative of the projected deviation, , which is given by eq. (32). Thus, we need to take the derivatives of and then with respect to :

 d^Bi(\boldmath{x}(τ))dτ=dxj(τ)dτ∇j^Bi(\boldmath{x}(τ))=^Bj(\boldmath{x}(τ))∇j^Bi(\boldmath{x}(τ))=Mij^Bj , dP⊥,ij(\boldmath{x}(τ))dτ=−d[^Bi(\boldmath{x}(τ))^Bj(\boldmath% {x}(τ))]dτ=−(M(I−P⊥)+(I−P⊥)MT)ij , (33)

where denotes matrix transpose. Let us also write down the following useful identities:

 Pn⊥=P⊥ , P⊥(I−P⊥)=0 , ^BiMij=0 ,P⊥M = M ,MTP⊥=MT . (34)

The second line above follows from . These properties of 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:

 dridτ=(MP⊥−(I−P⊥)MT)ijrj≡~mijrj , (35)

where we used the fact that . The term above takes into account the rotation of the plane orthogonal to due to the changing direction of , 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 , they can depend only on the first few derivatives of , and in particular. If all derivatives mattered, that would allow one to construct a Taylor expansion of (for analytic ) 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 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 (entering in eq. (31)) with itself and with . That already includes possible contractions with (eq. (35)) as it is constructed out of and . Thus, in what follows we are going to focus exclusively on the first derivative of the transverse deviation vector, .

That derivative does not necessarily lie in the plane normal to . To see that, one can multiply both sides of eq. (35) by to find that generally:

 P⊥d\boldmath{r}dτ≠d%\boldmath$r$dτ . (36)

The reason for the inequality is that, by construction, always lies in the plane normal to , and that plane changes orientation with .

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 -dependent direction. So, to study those transverse deviations, we can pick a -dependent 2D basis (subject to certain constraints; see below) which spans the changing normal plane, and write in that 2D transverse basis. Let us denote the components of in one of those (non-unique; see below) -dependent 2D transverse orthonormal bases as: . We can trivially extend that basis to 3D by including as the third orthogonal unit vector, so that the trivial extension of to 3D is and the components of in that basis are .

By construction, the components must be related to the components in the fixed Cartesian basis by a -dependent rotation matrix, . Extracting only the first two non-zero components, we can write (with the sum written out explicitly to highlight the upper limit):

 r♢i=3∑j=1Rij(τ)rj(τ), with i=1,2 . (37)

The matrix is not unique since we can always make additional rotations in the normal plane, while leaving the components of 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 using the constant matrix (which equals as is orthogonal). Let us denote the components of vectors and operators expressed in that rotated coordinate system with a superscript . Then the components, , of in that 3D basis are given by:

 r(τ0)i(τ)=RTij(τ0)Rjk(τ)rk(τ) . (38)

As we discussed in the beginning of this section, we will be focusing only on the first derivative of around some . For infinitesimal , we can write that derivative in finite difference form, which implies that we need to focus only on and . The former can be easily obtained in component form from the equation above: ; while the latter can be written out as:

 r(τ0)i(τ0+δτ)=RTij(τ0)Rjk(τ0+δτ)rk(τ0+δτ)≡δRij(δτ,τ0)rj(τ0+δτ) , (39)

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

The rotation matrix infinitesimally rotates the basis vectors for the quantities evaluated at , but leaves those evaluated at unchanged since as can be seen from the definition above. We loosely refer to this mixture of basis vectors at the two different locations, and , along the field line as the bent basis, as it “bends” from to in such a way so as to leave the components of unchanged (see the discussion above). In other words, is chosen such that .

Using the Rodrigues’ rotation formula999One can use the Rodrigues’ rotation formula to find the rotation matrix () in the plane defined by two unit vectors, and , that takes the components to the components . That matrix is given by: It is straightforward to check that indeed , and . 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:

 δR=I+δτ[(I−P⊥)MT−M(I−P⊥)] , (40)

with all operators evaluated at . Indeed, checking that to linear order in is a straightforward exercise, which guarantees that is a rotation matrix. Then, using eq. (33), one can check that to linear order in :

 ^B(τ0)i(τ0+δτ)=δRij^Bj(τ0+δτ)=δRij(δτ,τ0)(^Bj(τ0)+δτMjk^Bk(τ0))=^Bi(τ0) , (41)

where we used eq. (34).

Note that similar to , is not unique as we can additionally make a rotation in the plane normal to . We address that ambiguity in Section 3.3. However, let us first use the rotation matrix to calculate . First, note that eq. (35) can be written with finite differences as:

 \boldmath{r}(τ+δτ)=[I+δτ~m]\boldmath{r}(τ) . (42)

From that, we obtain:

 r(τ0)i(τ0+δτ) = δRijrj(τ0+δτ)=δRij[I+δτ~m]jkrk(τ0) (43) = [I+δτP⊥MP⊥]ijrj(τ0) ,

where we used equations (40) and (34). Therefore, in the bent basis, we can write the first derivative of the components of at as:

 dr(τ0)i(τ)dτ∣∣ ∣∣τ=τ0=r(τ0)i(τ0+δτ)−ri(τ0)δτ=(P⊥MP⊥)ijr(τ0)j≡m(τ0)ijr(τ0)j , (44)

where we used the fact that since in the bent frame, quantities evaluated at are not rotated by . We can see that , defined in the last equality above, is simply the transverse part of , 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 at in the 2D transverse basis. Let us apply the constant rotation matrix , which is independent of both and , to both sides of eq. (44):

 (45)

From equations (37) and (38), one can write the vector components , entering on both sides above, as: