Definition of the Spatial Propagator and Implications for Magnetic Field Properties

We present a theoretical framework to analyze the 3D coronal vector magnetic-field structure. We assume that the vector magnetic field exists and is a priori smooth. We introduce a generalized connectivity phase space associated with the vector magnetic field in which the basic elements are the field line and its linearized variation: the Spatial Propagator. We provide a direct formulation of these elements in terms of the vector magnetic field and its spatial derivatives, constructed with respect to general curvilinear coordinates and the equivalence class of general affine parameterizations. The Spatial Propagator describes the geometric organization of the local bundle of field lines, equivalent to the kinematic deformation of a propagated volume tied to the bundle. The Spatial Propagator’s geometric properties are characterized by dilation, anisotropic stretch, and rotation. Extreme singular values of the Spatial Propagator describe quasi-separatrix layers (QSLs), while true separatrix surfaces and separator lines are identified by the vanishing of one and two singular values, respectively. Finally, we show that, among other possible applications, the squashing factor [\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$Q$\end{document}Q] is easily constructed from an analysis of particular sub-matrices of the Spatial Propagator.

The build-up, storage, transport, and subsequent release of magnetic energy in the low-β solar corona is widely accepted as the basic requirement for solar coronal heating, the origin and generation of the solar wind, as well as eruptive phenomenology and space-weather prediction (see, e.g., Klimchuk, 2006, for a review). It is the geometric and topological features of the coronal magnetic field that govern these dynamics in the low-β coronal plasma environment (see Longcope, 2005, for a review). In general, non-trivial geometric and topological features such as null-point structure (Lau and Finn, 1990;Parnell et al., 1996), separator lines (Longcope and Cowley, 1996), bald patches (Seehafer, 1986;Wolfson, 1989;Low, 1992;Titov, Priest, and Démoulin, 1993), separatrix surfaces (Low, 1987;Somov, 1992;Priest, Heyvaerts, and Title, 2002), and quasi-separatrix layers (QSLs: Priest and Démoulin et al., 1996Démoulin et al., , 1997Titov, 1999;Milano et al., 1999) appear ubiquitously in the coronal magnetic field. In fact, the current state of global solar-wind generation models (see, e.g., Abbo et al., 2016, for a recent review), differ in the extent and complexity of these geometric and topological structures within the coronal magnetic field (e.g. Wang and Sheeley, 1990;Fisk, Schwadron, and Zurbuchen, 1998;Arge and Pizzo, 2000;Fisk, 2003;Cranmer, van Ballegooijen, and Edgar, 2007;Antiochos et al., 2011;. Since the vector magnetic field is reference-frame dependent, so too is the magneticconnectivity description (Hornig and Schindler, 1996); in fact, the entire concept of a magnetic-field line is not Lorentz invariant (Hornig, 1997). However, in a fixed reference frame, the geometric and topological features of the magnetic field constrain the dynamics. In general, the geometric features and topological constraints of the magnetic field are primarily important to understand where and how energy is stored and released in lowβ plasma environments such as the solar corona. In the presence of resistive (non-ideal) physics, field-line connectivity topology is no longer preserved, although separatrix structures remain. Separatrix surfaces and QSLs are associated with electric-current-sheet formation and reconnection, and therefore it is the geometry of these structures that constrain the storage of magnetic energy throughout the coronal volume, and determine both the location and the ability to release free energy (see the references above). In particular, QSL locations are related to observations of sudden flare brightening in Hα for a wide variety of solarflare phenomenology such as circular ribbon flares, two ribbon flares, and twisted flux rope (sigmoidal active region) morphologies (see, e.g., Janvier, 2017, and the references therein). Furthermore, separatrix surfaces and QSLs are the boundaries dividing regions of disparate connectivity in complex, multi-polar coronal structures. Hence, QSLs feature prominently in the rapid reorganization of the vector magnetic field and subsequent energy release in CME initiation mechanisms (see, e.g., Aulanier, Démoulin, and Grappin, 2005;Effenberger et al., 2011;Janvier et al., 2014;Schmieder, Aulanier, and Vršnak, 2015;Lynch et al., 2016, and the references therein).
The features of the magnetic connectivity map such as helicity, separatrix and QSL structures, and their effect on the system dynamics find their most explicit representations in terms of the algebraic and geometric descriptions of the field line and the fieldline bundle. In other words, the constraints on the system dynamics become transparent when cast in terms of the geometric structure and topological invariants of the field-line bundle and its behavior. Numerical magnetohydrodynamic (MHD) methodologies allow one to explore the locus of system dynamics for various magnetic-field models, thermodynamic heating models, and boundary conditions. While the MHD method (where applicable) is correct, the MHD equations and numerical solutions are often opaque to these explicit geometric structures (e.g. QSLs, separatrices, and separators) and their behavior. Moreover, many numerical difference schemes employ linear interpolation, which has the potential to limit these codes as an accurate representation of the system, especially in the vicinity of sharp gradients and other small-scale structures. Effectively, linear interpolation shifts the problem of the small-scale dynamics from the physical quantities on to finer and more complex numerical-grid resolutions. This is an extremely popular approach to numerical/computational solar-and space-plasma physics, but arguably leads to spurious effects (see, e.g., Edmondson, 2012, for a discussion); we offer no judgment regarding the veracity of these approaches. This perspective, however, requires the development of new mathematical tools/description/framework to analyze the geometric organization of the magnetic-field connectivity map.
The purpose of this article is the introduction, development, and presentation of a generalized connectivity phase space of field-line geometry associated to the vector magnetic field, in which the geometric structure and topological constraints are made explicit. This formalism does not alter the physics of Maxwell's equations, or MHD, but it is a fundamentally different framework that describes and analyzes vector magnetic fields. The basic assumptions of this framework are: i) the vector magnetic field is the primary (observable) quantity that satisfies some standard physical evolutions (e.g. Maxwell's equations, MHD induction, etc.), and the connectivity map is the secondary (derived) quantity; and ii) while large gradients may exist, the vector magnetic field is a priori smooth everywhere.
Under these assumptions, we derive a generally covariant measure of the local field geometry, called the Spatial Propagator, from the linearized variation of a field line. We demonstrate that the Spatial Propagator characterizes the geometric organization of a local bundle of field lines. Moreover, we identify topological invariants derived from the Spatial Propagator, as well as demonstrate the proper limiting connection to QSLs, separatrix surfaces, and separator lines. Beyond the limiting cases, the inclusion and analysis of existing and/or generated singular structures within the vector field are outside the scope of the present work.
The roadmap for this article is as follows: In Section 2 we lay out a precise mathematical definition for the field lines of a vector field (Section 2.1). We introduce the Spatial Propagator (Section 2.2) as a generalized spatial variation of an entire field-line solution, which describes the local bundle of field lines. Furthermore, we derive a direct relation between the Spatial Propagator and the local gradient of the vector field (Section 2.3), and hence the precise mathematical formulation of the geometric phase space, consisting of the integral curve solutions and their associated Spatial Propagators. Finally, we discuss the various representations of the Spatial Propagator (Section 2.4): covariance with respect to local curvilinear  coordinates, and equivalence with respect to affine transformations of the coordinate defined along the field line.
In Section 3 we characterize the geometric organization (dilation, rotation, anisotropic stretch, and connectivity gradient) of a local bundle of field lines using a mathematically equivalent kinematic analysis of a volume propagated along and deformed by the bundle. We explore the vector-field geometry using the Spatial Propagator in terms of volumetric dilation (Section 3.1), from which we identify a topological invariant measure that reflects the divergence-free condition of physical magnetic fields. We demonstrate that the singular values and singular vectors of the Spatial Propagator characterize the anisotropic stretch and rigid-body rotation deformations of the geometry (Section 3.2). Lastly (Section 3.3), we identify quasi-separatrix structures directly from the Spatial Propagator as extreme kinematic deformations, and we demonstrate the construction of the Q-factor (see, e.g., Titov, Hornig, and Démoulin, 2002;Titov, 2007) by a simple example.
We close with a description of other potential applications of the Spatial Propagator in Section 4. The major objects and main equations of this framework are summarized in Table 1.

The Integral Curve Description of Vector Field Geometry
Let M ⊆ R 3 be a subset of three-dimensional space, possibly with boundary ∂M. Let B (r) be a vector field for all positions r ∈ M. This field satisfies some set of dynamic evolution equations for some known initial and boundary conditions. Written with respect to a Cartesian 1 basis {ê x ,ê y ,ê z }, the vector field B (r) has component functions, (x, y, z), (x, y, z), where r is the spatial position with coordinates {x, y, z}.
In general, the vector field B (r) is time-dependent, and therefore, strictly speaking, so too are the integral curves (Section 2.1), and the Spatial Propagator (Section 2.2), as well as all other objects constructed therefrom. For ease of notation, throughout this work we suppress the functional time-dependence of all quantities; this may be interpreted as analyzing the system at a fixed time, or for very-low frequency dynamics, f c/L where c is the characteristic speed of communication, and L is a characteristic length scale. The dynamical description of the integral curves (Section 2.1), the Spatial Propagator (Section 2.2), etc. require a treatment of the full four-dimensional electromagnetic-field tensor (Jackson, 1999, Section 11.9), which is outside the scope of this work.

The Integral Curves of a Vector Field
In differential equation theory, a general vector field B (r) is everywhere tangent to a set of integral curves r (λ, r 0 ) that satisfy the initial value problem The family of integral curve solutions r (λ, r 0 ) are parametrized by a real number λ ∈ R, referred to as the connectivity parameter, and a three-component vector r 0 ∈ M, referred to as the reference point. A single integral curve, as a particular solution to Equations 2 and 3, is identified by a single fixed reference point r 0 . The connectivity parameter λ denotes the distance per unit field strength along a particular solution curve issuing from a particular reference point. We note that the rate of change with respect to λ is simply the directional derivative along the vector field, which may be written with respect to the coordinate representation The integral curve solutions to Equations 2 and 3 represent position vectors r = r(λ, r 0 ), the components of which are functions of (λ, r 0 ); written with respect to a Cartesian basis {ê x ,ê y ,ê z } r(λ, r 0 ) ·ê x = (r| λ,r 0 ) x = x(λ, x 0 , y 0 , z 0 ), r(λ, r 0 ) ·ê y = (r| λ,r 0 ) y = y(λ, x 0 , y 0 , z 0 ), r(λ, r 0 ) ·ê z = (r| λ,r 0 ) z = z(λ, x 0 , y 0 , z 0 ).
(5) 1 We remark, unless specifically noted otherwise, throughout this article we work with respect to the standard Cartesian coordinates {x, y, z} and orthonormal Cartesian unit basis vectors {ê x ,ê y ,ê z }.

Figure 1
Illustration of a general congruence, viz. field lines and reference points. A particular field line r(λ, r 0 ) is associated with a single fixed reference point r 0 . A congruence is defined by a set of reference points in some neighborhood The integral curve solutions r(λ, r 0 ) are known by various names depending on the nature and interpretation of the differential equations and vector field. In physics, integral curves for electric and magnetic fields are known as field lines, whereas the integral curves for a velocity field are called streamlines. In general dynamical systems theory, the integral curves of the governing differential equation system are referred to as trajectories or orbits.
Since we find application to solar magnetic fields, throughout this article we use field line to denote a particular solution to Equations 2 and 3 with fixed r 0 . Moreover, a flux tube is a congruence of integral curves (or simply a congruence); that is a local bundle of field lines defined by a particular set of reference points in some neighborhood r 0 ∈ 0 ⊆ M (see Figure 1).
We make a few observations and define some nomenclature regarding the field lines in physical systems of interest. The reference point r 0 is a free parameter, typically chosen on the system boundary, or to coincide with some known initial state within the system interior; for example, in a magnetized plasma this choice typically coincides with a highly-conductive parcel of plasma material. The λ = 0 datum defined by the reference point r 0 = r(0, r 0 ) is typically referred to as the launch footpoint in magnetic systems. The point r = r(L, r 0 ) for a finite λ = L, corresponding to the final point of integration of Equations 2 and 3 is often referred to as the target footpoint in magnetic systems; typically, the target footpoint is where the integral curves crosses the system boundary, or encounters a singular structure in the vector field.
The congruence of integral curves of a smooth vector field represents an equivalence class under general affine transformations of the connectivity parameter of the form λ → is a smooth, positive definite, scalar-valued function, 2 and b ∈ R an arbitrary constant. 3 The re-parametrized flow r( , r b ) satisfies Equations 2 and 3 for vector field X(r) = B(r)/f (r) and reference condition r b ; that is, 2 The function f must be at least C 1 differentiable. The positivity condition preserves field-line orientation. Furthermore, a solenoidal vector field ∇ · B = 0 with X = B/f leads to a constraint ∇ · X + X · ∇f = 0; then choosing a function f such that the vector field X is everywhere tangent to the level surfaces of f (i.e. X · ∇f = 0) preserves the solenoidal property to the vector field X. Furthermore, including a temporal dependence on f imposes a constraint equation on the X-evolution to preserve consistency with B. 3 The value of the constant b represents a simple relabeling of the reference point on the particular field-line solution; i.e. r 0 → r b .
A simple application of this λ → (λ) re-parametrization is the analysis of the B field in a system with boundaries ∂M. One has the freedom to choose the function f (r) (and constant b = 0) in order that the values = 0 and = 1 coincide with the initial and final points of the field lines taken at the boundary surfaces, r 0 ∈ ∂M and r(1, r 0 ) ∈ ∂M; typically at the photosphere for solar coronal applications. With = 0 and = 1 as boundary values, the connectivity parameter is not the physical (dimensional) arc length, but rather a normalized dimensionless arc-length parameter.
As a second application, consider f (r) = |B(r)|; assuming |B(r)| = 0. The vector field X(r) = B(r) |B(r)| = b(r) is well defined and is identified with the unit magnetic-field direction. For simplicity and without loss of generality, we may set the constant b = 0. Under this transformation, the system of Equations 2 and 3 becomes In this example, the connectivity parameter represents the physical arc length along the magnetic-field line. We refer to the representation constructed from the unit magnetic-field direction vector field and in which the connectivity parameter represents the physical arc length along the field line issuing from the reference point, as the arc-length representation.

The Spatial Propagator: A Local Congruence of Integral Curves
The relative geometry of a local congruence (i.e. dilation, anisotropic stretch, and rotation of the bundle of curves) is completely described by examining a field line under a spatial variation of the reference point r 0 → r 0 + h. We denote the spatial variation of the reference point by the reference shift vector h.
Recall that a fixed reference point r 0 is equivalent to a single field line r(λ, r 0 ); hence, when considering solutions r(λ, r 0 ) to Equations 2 and 3 for a set of r 0 ∈ 0 , a variation in the reference point r 0 → r 0 +h is equivalent to a variation of the entire field line. To see this, consider two spatially neighboring field-line solutions within the congruence, respectively: r(λ, r 0 ) and r(λ, r 0 + h). We may relate the component functions of the neighboring field lines by a Taylor expansion, such that for all λ and |h| h m , where h m is a local characteristic length scale defined by comparing the magnitude of the first-order with the higher-order variational terms (see Appendix A). Equation 10 describes the local organization of all field lines with reference points within an initial volume 0 of characteristic size | 0 | 1/3 h m ; that is to say, the congruence is local with respect to reference points r 0 + h ∈ 0 with |h| h m .
We remark, there is an implicit assumption in Equation 10 that the magnetic vector field B(r) is described by smooth component functions everywhere within the domain. In future work, we will explore the consequences of relaxing this assumption.
We define this first-order variation of the field line r(λ, r 0 ) with respect to a spatial variation in the reference point r 0 to be the Spatial Propagator F(λ, r 0 ). Then for all λ, the variational derivative may be represented as a 3 × 3 matrix whose component functions are simply the derivatives of Equation 5 with respect to the reference point components, Note that, at λ = 0, the field line reduces to the reference point r(0, r 0 ) = r 0 , and similarly r(0, r 0 + h) = r 0 + h. Hence, Equation 10 reduces to the definition of the reference shift vector r(0, r 0 + h) − r(0, r 0 ) = h, and Equation 11 at λ = 0 is simply the 3 × 3 identity matrix, (F| 0,r 0 ) i j = δ i j ; where δ i j is the Kronecker delta. The Spatial Propagator F(λ, r 0 ) may be considered to be the generalized gradient of an entire field line r(λ, r 0 ) issuing from the reference point r 0 for all λ. The difference vector v(λ, r 0 ) represents the spatial shift at each λ along the particular field line r(λ, r 0 ) to the corresponding point at λ along the neighboring field line r(λ, r 0 + h) (see Figure 2); the components of v(λ, r 0 ) are given by For a given spatial variation h, Equation 12 may be interpreted as a generalized directional derivative of an entire field line r(λ, r 0 ) in the direction v(λ, r 0 ). In other words, the Spatial Propagator F "propagates" every spatial variation h of the reference point r 0 along the particular field line in 3D space. The matrix representation of the Spatial Propagator F(λ, r 0 ) contains all of the geometric information of the local congruence; that is, it carries the local geometric organization of the field lines within a neighborhood 0 h 3 m of a particular field-line solution r(λ, r 0 ). We give a precise meaning to this statement in Section 3.

Direct Relation Between the Spatial Propagator and the Vector Field
Like the field lines of the congruence, the Spatial Propagator F(λ, r 0 ) is a function of both the connectivity parameter λ and reference point r 0 . In order to calculate the Spatial Propagator, one may integrate all field lines, and then construct the difference Equation 10 between any two neighboring field lines, taking care to evaluate the scale length h m for every field line. However, in typical physical systems of interest the governing equations describe the evolution of the vector field (e.g. Faraday's law, MHD induction, etc.), and the associated field lines are derived therefrom. Hence, we seek a formulation of the Spatial Propagator directly from the vector field, which allows the simultaneous calculation of the spatial behavior of all field lines within the h m neighborhood.
To illustrate, we give a simple, qualitative derivation as follows. The connectivity parameter λ and the reference point r 0 may be considered independent variables in the family of smooth integral curve solutions r(λ, r 0 ) to the system of Equations 2 and 3. From this perspective, the spatial variation with respect to the reference point may be constructed directly from Equation 2 by Since the vector field B(r) is assumed to have smooth components, then the family of integral curve solutions r(λ, r 0 ) is smooth in the variables λ and r 0 ; hence, mixed partial derivatives commute, and Equation 13 may be written ∂ ∂λ By Equation 11, Equation 14 is a matrix evolution equation for the (unknown) Spatial Propagator F(λ, r 0 ). The ∂B(r) ∂r term is simply the rate of change of the vector field with respect to the spatial position r = r(λ, r 0 ); represented as a 3 × 3 matrix with respect to Cartesian basis, where each partial derivative component is evaluated along the field-line component functions, Equation 5. Substituting Equations 11 and 15 into Equation 14, the Spatial Propagator F(λ, r 0 ) satisfies the following first-order matrix differential equation, in components with respect to a Cartesian basis {ê x ,ê y ,ê z }: where δ i j is the Kronecker delta. We remark that by Equation 4 the matrix ODE, Equation 16, is the Lie transport (see, e.g., Kobayashi and Nomizu, 1963, p. 29) of the propagated shift vector v along the vector field B, where the reference shift vector h is a constant vector independent of the coordinates; that is (B · ∇)h = 0. In general, Equations 2 and 16, along with their respective initial conditions given by by Equations 3 and 17, are generally referred to as the equations of variation (see, e.g., Arnold, 1992, pp. 223 -225), and they constitute a well-posed problem. Hence, the Spatial Propagator F(λ, r 0 ) solution exists and is unique for all λ and r 0 (see, e.g., Bernstein, 2018Bernstein, , pp. 1193Bernstein, -1195, for existence and uniqueness proof).
For any reference point r 0 and all λ, the congruence of scale h m is a particular solution to the equations of variation consisting of both the integral curve r(λ, r 0 ) and corresponding Spatial Propagator F(λ, r 0 ). The congruence may be interpreted as a local phase-space of scale h m for the field-line connectivity consisting of elements r(λ, r 0 ) and F(λ, r 0 ). With respect to the global Cartesian coordinates, these phase space elements are represented by a three-component vector and 3 × 3-component matrix, where {ê x ,ê y ,ê z } are the orthonormal basis vectors, and {ê x ,ê y ,ê z } the dual basis (see Appendix B.1).

General Representations of the Spatial Propagator
There are two fundamental representation categories of the Spatial Propagator F(λ, r 0 ): the matrix representation with respect to a particular basis setê i , and the representation with respect to the connectivity parameter λ. First, we consider the basis representation of the congruence; that is the field line r(λ, r 0 ) and associated Spatial Propagator F(λ, r 0 ) are, respectively, a three-component vector and 3 × 3 matrix representation constructed with respect to a particular basis set. Fundamentally, the matrix representation of the congruence follows from the vector field B(r) and its gradient matrix ∇B(r); that is constructing B and ∇B with respect to Cartesian basis vectors lead to the Cartesian representation of the congruence: Equation 19.
In the Cartesian representation, the basis vectors are independent of the spatial position and hence Equations 11 and 15 find a particularly simple form. However, general curvilinear basis vectors are position dependent (see Appendix B.1); for example, consider the spherical-polar orthonormal 4 basis vectors, e r (θ, φ) = sin θ cos φê x + sin θ sin φê y + cos θê ẑ e θ (θ, φ) = cos θ cos φê x + cos θ sin φê y − sin θê ẑ prevalent within solar and space physics. In order to incorporate the spatial position dependence of general curvilinear basis vectors in this formulation, we impose the condition that the Spatial Propagator F(λ, r 0 ) transform as a tensor under general coordinate transformations.
The tensor-condition requires the ∇B(r) object to be the matrix representation of the covariant differential of the vector field (see, e.g., Kobayashi and Nomizu, 1963, pp. 143 -144). The components of the matrix representation of the covariant differential of the vector field with respect to general local curvilinear coordinates q i is given by (see Appendix B.1 for derivation and discussion). The first term is simply the rate of change of the component functions with respect to the local coordinates. The ( | r ) i jk in the second term are coefficients that account for the spatial position dependence of the basis vectors q i (q j ). Equation 16 with Equation 21, and the initial condition Equation 17, generalizes the field-line deviation to general curvilinear coordinates (see, e.g., Tassev and Sevcheva, 2017, for application in spherical-polar orthonormal coordinates). In Appendices B.2 and B.3 we develop explicit matrix representations of the covariant differential of the vector field ∇B(r) with respect to the common spherical-polar bases used by the solar and space-physics community (see, e.g., Tassev and Sevcheva, 2017).
Second, we consider the connectivity parameter representation of the congruence. For a given vector field B(r), the congruence solution r(λ, r 0 ) and F(λ, r 0 ) to Equations 2 and 16 with respect to λ is said to be given in the natural representation. That is the connectivity parameter λ has units of arc length per vector field magnitude, and the congruence elements, r(λ, r 0 ) and F(λ, r 0 ), reflect this functional dependence.
Recall from the end of Section 2.1, the congruence represents an equivalence class under general affine transformations of the connectivity parameter, λ → (λ) = f (r)λ + b for smooth, positive-definite functions f (r) and scalars b. Hence, the matrix differential Equations 16 and 17 transform accordingly; with respect to Cartesian coordinates where X(r) = B(r) /f (r), and the positions are evaluated along r = r( , r b ). Equation 22 with initial condition 23, generalizes the field-line deviation to include affine transformations in the connectivity parameter representations (see, e.g., Tassev and Sevcheva, 2017; Scott, Pontin, and Hornig, 2017, for application). Explicitly, recall the arc-length representation example, Equations 8 and 9, with f (r) = |B(r)|, and |B(r)| = 0. The vector field X(r) = b(r) is well defined and identified with the unit magnetic-field direction. In this example, the arc-length representation of the Spatial Propagator with respect to Cartesian coordinates satisfies where the positions are evaluated along the arc-length representation of the field line, r = r( , r 0 ).
In the application to solar magnetic fields, constructing the congruence solution from data (e.g. Solar Dynamics Observatory/Atmospheric Imaging Assembly 171 Å images) requires the solutions to Equations 8 and 24, with initial conditions 7 and 25, respectively. From this perspective, we have decoupled the magnetic-field strength estimation from the field-line trajectory estimation in the construction of the Spatial Propagator.

Geometric Deformation of a Congruence
Geometry describes the measurable lengths and angles associated with the system configuration. The value of the congruence formalism is that all of the geometric behavior of a local bundle of field lines is contained within a single integration of Equations 2 and 16, along with their respective initial conditions 3 and 17. The field line r(λ, r 0 ) and corresponding Spatial Propagator F(λ, r 0 ) implicitly incorporate all geometric information of the local bundle, which follows naturally from the vector magnetic-field structure: B and ∇B. The generic geometric configuration of the congruence may be described by a combination of dilation (Section 3.1), stretch, rotation (Section 3.2), and gradients in the connectivity structure (Section 3.3).
We recall that by assumption the vector field components are smooth functions. Hence, by standard theorems of existence, uniqueness, and extension for ordinary differential equations (see, e.g., Hirsch and Smale, 1974;Arnold, 1992;Taylor, 1996), each field-line reference point r 0 ∈ 0 is mapped smoothly and uniquely to the point r = r (λ, r 0 ) ∈ λ (see Figure 1); that is the neighborhood volume 0 is mapped, smoothly and uniquely, to the neighborhood volume λ . Hence, the geometric deformation of the local congruence is reflected in the deformation of the initial volume 0 by propagation along the congruence into λ .

Congruence Dilation and the Determinant of the Spatial Propagator
The dilation of a congruence is the compression/expansion of the constituent field lines, and it is completely described by the determinant of the 3 × 3, non-singular (invertible) matrix representation of the Spatial Propagator F(λ, r 0 ). This type of geometric deformation is quantified by the dilation of a volume λ ⊆ R 3 propagated through each λ along the congruence. In this section, we make explicit the functional dependence of the deformed volume on both the connectivity parameter λ and the reference point r 0 , by denoting λ = (λ, r 0 ). The signed differential volume element in Cartesian coordinates is constructed from the coordinate differentials dx iê i , We may choose differential-reference-shift vectors to coincide with the Cartesian coordinate differentials dh i ≡ dx iê i , and hence write the local signed differential volume element at the reference point d (0, r 0 ) = d 3 x.
Each differential-reference-shift vector dh i may be propagated along the congruence under the action of F(λ, r 0 ) (see Figure 3), such that where for each i = x, y, z, the F(λ, r 0 ) ·ê i is a vector; expanded with respect to Cartesian basis (see Appendix C.1), Hence, the local signed differential volume element propagated along the congruence is completely determined by the action of the Spatial Propagator on the basis vectors: By a standard result from vector calculus in three dimensions, the ratio of propagated volume element to initial volume element is given by the determinant of the matrix representation of F(λ, r 0 ): In Appendix C.1, we calculate Equation 30 explicitly for a Cartesian basis. Furthermore, identity 30 is valid for arbitrary orthonormal basis sets (see, e.g., Nickerson, Spencer, and Steenrod, 1959, pp. 95 -97, for a standard treatment). From a geometric perspective, the determinant of the 3 × 3 matrix representation of F(λ, r 0 ) acts as the Jacobian in a change of basis for the signed volume element propagated along the congruence (see Figure 3). Since the differential-reference-shift vectors dh i are mapped under the action of the Spatial Propagator into the propagated differential vectors dv i (λ, r 0 ), the signed differential volume elements d (0, r 0 ) and d (λ, r 0 ) are centered on the respective points r 0 = r(0, r 0 ) and r = r(λ, r 0 ). Hence, 3D integration with respect to the volume measure d (λ, r 0 ) generates a total volume (λ, r 0 ), deformed with respect to the initial shape (0, r 0 ) and centered on the field line r(λ, r 0 ) at each value of λ. Moreover, by Equation 30 the deformed total volume at every λ may be cast as an integral with respect to the reference volume measure d (0, r 0 ) at the reference point r 0 , The dilation of the total volume (λ, r 0 ) by propagation along the congruence is therefore completely determined by the determinant of the Spatial Propagator det F(λ, r 0 ) at each λ. Geometrically, the total volume (λ, r 0 ) decreases (volumetric compression) for 0 < det F(λ, r 0 ) < 1, remains constant (isochoric propagation) for det F(λ, r 0 ) = 1, or increases (volumetric expansion) for det F(λ, r 0 ) > 1, respectively.
We remark, det F(λ, r 0 ) < 0 corresponds to a mapping from a proper right-handed basis set, to a left-handed basis set. Such an F(λ, r 0 ) is discontinuous and requires singularities in the vector field B, which we do not consider in this work.
Moreover, using Equation 31, we can quantify the rate of volumetric dilation by propagation along the congruence, where we have made use of the identity (see Appendix C.2 for derivation).
In Equations 32 and 33, the Tr(∇B(r) ) term is evaluated at each point r = r(λ, r 0 ). In particular, every physical magnetic vector field is everywhere divergence-free: Tr(∇B(r)) = ∇ · B(r) = 0 for all r ∈ M. Hence, by Equations 17 and 33, for all finite λ, This implies that the total volume is conserved under propagation along a congruence; that is (λ, r 0 ) = (0, r 0 ) is constant regardless of deformation. Hence, we identify det F(λ, r 0 ) = 1 for all λ and every r 0 as the simplest topological invariant of a congruence generated by a smooth magnetic vector field. This true topological invariant reflects the divergence-free condition.

Congruence Stretch and Rotation and the Singular Value Decomposition of the Spatial Propagator
The anisotropic stretch and rigid-body rotation of a congruence are equivalent to the kinematic description of the propagated volume undergoing similar deformation; that is, the λ-evolution of the propagated vectors v(λ, r 0 ) with respect to a local volume-centered orthogonal reference frame. The congruence stretch and rotation is described by a deformation scaled along mutually orthogonal, body-centered axes of the propagated volume. This basic geometric picture is illustrated in Figure 4. We seek a geometric description in which the orthogonality of the reference frame is preserved under rotation and the lengths of the direction vectors are scaled. This description initially suggests the λ-evolution of the eigen-decomposition of the Spatial Propagator matrix F, in which the eigenvalues and eigenvectors play the role of scale factors and basis Figure 4 Illustration of an anisotropic stretch and rigid-body rotation of an orthogonal frame centered about the central field-line trajectory r(λ, r 0 ) within the propagated volume. In contrast to Figure 3, a specific orthogonal framer α is shown centered on r(0, r 0 ). It is rotated along the central field line into the specific orthogonal framel α centered at r(λ, r 0 ) while the vector lengths are simultaneously scaled by the corresponding σ α as defined in Section 3.2.
directions. The eigenvalues for a general, non-singular, 3 × 3 matrix F are determined by a cubic polynomial with real coefficients, the roots of which may be real or complex-valued depending on the sign of the classical cubic discriminant. In the case of non-negative cubic discriminant, the eigenvalues are all real (perhaps with multiplicity > 1), the geometric interpretation is scaling along linearly independent eigenvectors (since F is a non-singular matrix). However, mutual orthogonality of the eigenvectors follows only in the special case that the Spatial Propagator matrix is symmetric, F = F T (where F T denotes the matrix transpose). Moreover, a negative cubic discriminant corresponds to complex-valued eigenvalues and eigenvectors, the geometric interpretation of which is not a stretch and rotation of the propagated volume. Hence, the eigen-decomposition of the Spatial Propagator F does not provide the correct stretch and rotation description.
To determine the correct description, consider that a general stretch is described by a change in the length of propagated vectors v (λ, r 0 ) within the propagated volume. The length of a vector is described by the norm Equation 35 is a real-valued function of λ that involves the symmetric matrix F T · F, as opposed to the Spatial Propagator F alone. By a standard theorem of linear algebra (see, e.g., Halmos, 1958, Section 79), the eigenvalues of the symmetric matrix F T · F are all real (possibly with multiplicity > 1), and the eigenvectors are everywhere mutually orthogonal. More generally, the symmetric matrices F T · F and F · F T have identical eigenvalues, as well as orthogonal, albeit different, eigenvector bases. Hence, for a general, non-singular Spatial Propagator F, the real eigenvalues and real orthogonal eigenvectors of the symmetric matrices F · F T , respectively F T · F, may be used to infer indirectly the geometric stretch and rotation interpretation of F.
Formally, for each r 0 and all λ, there exists a factorization of F(λ, r 0 ), called the singular value decomposition (SVD), given by where P(λ, r 0 ) is a 3 × 3 diagonal matrix, and the R l (λ, r 0 ) and R r (λ, r 0 ) are 3 × 3 orthogonal matrices (e.g. Bernstein, 2018, pp. 555 -558). The diagonal entries of the matrix P(λ, r 0 ) are non-negative functions σ α (λ, r 0 ) for α = 1, 2, 3, given by where q α (λ, r 0 ) are the eigenvalues of the matrices F T · F and F · F T . The columns of R l (λ, r 0 ), denoted byl α (λ, r 0 ), are the mutually orthogonal eigenvectors of the symmetric matrix F · F T ; respectively, the columns of R r (λ, r 0 ), denoted byr α (λ, r 0 ), are the mutually orthogonal eigenvectors of the symmetric matrix F T · F. For each reference point r 0 and fixed λ, the set of values σ α for α = 1, 2, 3 are called the singular values of F(λ, r 0 ), and the corresponding set of vectorsl α andr α are called the singular vectors of F(λ, r 0 ). The anisotropic stretch and rigid-body rotation of the propagated volume, equivalently the congruence geometry, is completely determined by the singular values and singular vectors of the Spatial Propagator F(λ, r 0 ). Geometrically, an orthogonal framer α in the reference volume 0 is rotated into the orthogonal framel α in the propagated volume λ while the vector lengths are simultaneously scaled by the corresponding σ α , as shown in Figure 4.
We may decompose this simultaneous stretch and rigid-body rotation effect through a polar decomposition of the Spatial Propagator. Since R l (λ, r 0 ) is an orthogonal matrix (R −1 = R T , where superscript "T " denotes the matrix transpose) for each r 0 and all λ, then R T l (λ, r 0 ) · R l (λ, r 0 ) = I, where I is the identity. Hence, we may rewrite Equation 36 as By the exact same arguments, making use of the orthogonality of R r (λ, r 0 ), we find an equivalent rewrite of Equation 36 as The equivalent decompositions in Equations 38 and 39 are referred to as the left-polar decomposition and right-polar decomposition, respectively, of the Spatial Propagator (see, e.g., Halmos, 1958, Section 83). Unlike the general Spatial Propagator F(λ, r 0 ) itself, the matrices V(λ, r 0 ) and U(λ, r 0 ), denoted, respectively, the left-stretch and right-stretch, are symmetric, positive definite, 3 × 3 matrices defined by Moreover, the matrix R(λ, r 0 ) is a 3 × 3, orthogonal matrix defined by The left-stretch V(λ, r 0 ), right-stretch U(λ, r 0 ), and rotation R(λ, r 0 ) matrices all naturally inherit their (global) coordinate representation from the matrix representation of the Spatial Propagator F(λ, r 0 ). Equation 40 represent two coordinate rotations by the orthogonal matrices R l (λ, r 0 ) and R r (λ, r 0 ) that diagonalize, respectively, the V(λ, r 0 ) and U(λ, r 0 ) matrices. The anisotropic stretch and rigid-body rotation geometric interpretation, however, becomes explicit by application of the standard spectral decomposition theorem of linear algebra (see, e.g., Halmos, 1958, Section 79). Since the matrices V(λ, r 0 ) and U(λ, r 0 ) are symmetric, there exist orthogonal eigen-basesl α (λ, r 0 ) andr α (λ, r 0 ), with respect to which the matrix representations are By Equations 40 and 41, the eigen-basesl α (λ, r 0 ) andr α (λ, r 0 ), of the left-stretch V(λ, r 0 ), respectively right-stretch U(λ, r 0 ), matrices are the singular vectors of the Spatial Propagator F(λ, r 0 ).

Figure 5
Illustration of the polar decomposition of the Spatial Propagator. The right-polar decomposition, F = R · U, stretches/compresses by a factor of σ α along the corresponding directionr α , followed by a rotation aligning withl α ; the left-polar decomposition, F = V · R, rotateŝ r α to align withl α , followed by a stretch/compression by a factor of σ α along the corresponding directionl α .
Furthermore, inverting Equation 43 yields, Using Equations 42 and 44 in the left-polar decomposition Equation 38, the Spatial Propagator F(λ, r 0 ) matrix may be represented with respect to the singular vector basisl α (λ, r 0 ) andr α (λ, r 0 ), Geometrically, for each λ, the action of the matrices V(λ, r 0 ) and U(λ, r 0 ) is to scale the propagated volume by the singular values of the Spatial Propagator F(λ, r 0 ), and the action of the rotation matrix R(λ, r 0 ) is to rotate the singular directions of the reference state 0 to align with the singular directions of the propagated state λ (see Figure 5). Thus, the symmetric matrices V(λ, r 0 ) and U(λ, r 0 ), and orthogonal matrix R(λ, r 0 ), respectively, characterize the smooth anisotropic stretch and rigid-body rotation of the propagated volume along the congruence.
The basisl α (λ, r 0 ) describes the singular axes of the stretched and rotated state λ . Similarly, for any λ the basisr α (λ, r 0 ) identify the singular axes of the stretched but unrotated reference state 0 . Hence, for any λ the basisr α (λ, r 0 ) are the singular directions of the reference state 0 . We note, however, a spherical reference state 0 has no preferred directions; this follows from the fact that the initial condition F(0, r 0 ) = I possesses a single degenerate singular value (eigenvalue) σ = 1 of multiplicity 3. The resolution of this is that a unique, non-degenerate, orthogonal basisr α (λ, r 0 ) in the reference state 0 only exists and follows from a Spatial Propagator F(λ, r 0 ) with non-degenerate singular values.
From a computational standpoint, since F is non-singular, the matrices V, U, and R may be determined directly from F, F · F T , and F T · F. Using Equations 38 and 39, and that R is an orthogonal matrix, V, U, and R follow immediately: Since the matrices F · F T and F T · F are diagonalizable with respect to their respective eigenbases, the square-root operation is well defined. It may be shown the full 3 × 3 matrix representation of the Spatial Propagator F(λ, r 0 ) is directly dependent on the current distribution. Consider a quasi-static magnetic field B(r) with a non-trivial current distribution J (r) that satisfies Ampere's Law: μ 0 J (r) = ∇ × B(r). The covariant differential of the vector field ∇B may be decomposed into symmetric and anti-symmetric parts, where the superscript T denotes the matrix transpose. The matrix J × (r), denoted "Jsupercross", is an antisymmetric, 3 × 3 matrix representation of the current distribution (with respect to a Cartesian basis), The Spatial Propagator governing matrix ODE Equation 16 may be written The initial condition Equation 17 remains unchanged. We note that dropping the quasisteady assumption leads to a similar ODE with the addition of a "supercross" displacementcurrent-density matrix term; we do not pursue such high-frequency dynamics here.
Since the left-stretch V(λ, r 0 ) and right-stretch U(λ, r 0 ) matrices are symmetric, through Equation 49, the rotation matrix R(λ, r 0 ) implicitly requires a non-trivial anti-symmetric part of the covariant differential ∇B(r), equivalently a non-trivial (parallel) current distribution J (r). Moreover, under reasonable conditions we identify a class of magnetic fields B(r) for which the rotation matrix R(λ, r 0 ) depends only on the parallel current distribution J (r). Geometrically, the congruence rotation describes the twist of all neighboring integral curves about a central axis field line, and is therefore a measure of the twist helicity. The relationship between the Spatial Propagator, rotation matrix, (twist-) magnetic helicity, and (parallel) current distribution is beyond the scope of this work.

Quasi-Separatrix Layers and the Q-Factor from the Spatial Propagator
We have shown in the previous Sections 3.1 and 3.2 that the congruence geometry is equivalent to the kinematics of a volume undergoing dilation, stretch deformation, and rigid-body rotation under the action of the Spatial Propagator F(λ, r 0 ). Hence, a separatrix surface reflects the total-collapse of the 3D volume into a 2D surface as it is propagated by the congruence. The quasi-separatrix layer is similarly identified by an "extreme" kinematic deformation of the 3D volume; "extreme" being a subjective term and requiring a (user-defined) quantified threshold.
We remark, recalling from Section 3.1 that the divergence-free condition for vector magnetic fields leads to the topological invariant det F(λ, r 0 ) = 1. Geometrically, the volumetric kinematics are incompressible; 5 that is, the 3D shape deforms under the action of F(λ, r 0 ) while the total volume remains constant. Hence, the separatrix surfaces in a vector magnetic field are described by infinite stretch deformation of the volume.
For a fixed reference point r 0 , the 3D stretch deformation of a congruence at each λ is quantified by the three non-negative singular values σ α ∈ R ≥0 of the Spatial Propagator F(λ, r 0 ). There is an inter-dependence between the singular values that is governed by the divergence-free condition. The relative magnitudes of the singular values effectively reflect the ∇B matrix eigenvalue structure; the full algebraic proof is beyond the scope of this work. Hence, from a geometric perspective, the relative strengths of the set of singular values may be used to identify 2D separatrix surfaces, 1D separator lines, and 0D null points within the field structure. These topological features are identified by the basic geometric interpretations of the singular values following from the following features.
i) If any one singular value becomes zero, σ i → 0, then the other two correspondingly diverge to infinity, σ j , σ k → ∞; geometrically, the 3D field-line bundle has collapsed into a 2D surface. ii) If any two of the singular values become zero, σ i , σ j → 0, then the remaining one correspondingly diverges to infinity, σ k → ∞; geometrically, the 3D field-line bundle has collapsed into a 1D line. iii) The divergence-free condition prevents all three singular values from approaching zero simultaneously; geometrically, the 3D incompressible field-line bundle cannot collapse to a 0D point.
For example, in the natural representation the fan surface and spine line(s) associated with a magnetic null are identified in the limit as |λ| → ∞ by the vanishing of one and two singular values, respectively; the sign of λ depends on the field-line polarity. More generally, however, the identification of other non-trivial topological features that exhibit the vanishing of one, or multiple, singular values at finite and/or infinite |λ| remains an open question.
Up to this point, we have investigated the geometric meaning, consequences, and relations between the singular values of the full-3D Spatial Propagator in various limiting cases. The QSL is often described in terms of the "squashing" of a flux tube (see, e.g., Titov, 2007;Tassev and Sevcheva, 2017;Scott, Pontin, and Hornig, 2017); essentially, a relational measure between the eccentricity of the flux-tube cross-sectional area at the system boundaries, and hence is inherently 2D. Moreover, the squashing of a flux tube is defined entirely independent of rotation. For the remainder of this section, we demonstrate the dimensional reduction of the full-3D Spatial Propagator to describe the 2D squashing of the congruence. In particular, we illustrate a simple construction of the popular squashing factor Q (Titov, Hornig, and Démoulin, 2002;Titov, 2007), from the Spatial Propagator. We note, our aim here is neither a reformulation, nor a more efficient computation of the squashing factor Q over that available in the open literature (see, e.g., Tassev and Sevcheva, 2017). Rather the purpose is simply to show exactly how the popular Q-value may be derived from, and therefore fits within, the general Spatial Propagator framework.
Throughout this article we have assumed that the vector magnetic field B(r) has smooth component functions, Equation 1, that may be described with respect to a single global Cartesian coordinate chart. Hence the congruence solution components, Equation 19, described with respect to the same global Cartesian coordinate chart, are also smooth, singlevalued functions of the connectivity parameter λ, and reference condition coordinates r 0 . Furthermore, by construction, this framework is valid for any global coordinate chart that covers the system (see Appendix B.1 for general coordinate chart formulations, and Appendices B.2 and B.3 for spherical-polar formulations).
Recall from Section 2.1 that we identify λ = 0 with the reference point r 0 , corresponding to the beginning of the congruence, and we identify the fixed, finite value λ = L with r = r(L, r 0 ), corresponding to the end of the congruence. In general, we take the reference point r 0 and end point r(L, r 0 ) on the system boundary. 6 From this perspective, each particular field-line solution with L > 0 provides a unique connectivity map between disjoint points on the system boundary, and the Spatial Propagator quantifies the local geometric organization of this connectivity map within its characteristic scale.
Consider for simplicity a smooth vector field B(r) in a bounded Cartesian box; that is, we let the coordinate domains be x ∈ [−L x , L x ], y ∈ [−L y , L y ], and z ∈ [−L z , L z ]. The (inward) unit normaln at the boundaries is parallel (up to a sign) to the respective Cartesian basis vectors {ê x ,ê y ,ê z }. Moreover, consider a congruence solution, r(λ, r 0 ) and F(λ, r 0 ) (Equation 19) with characteristic scale h m , consistent with the given vector field B(r) such that the reference point r 0 is chosen at the z = −L z boundary, and the end point r(L, r 0 ) follows on the z = L z boundary (see Figure 6).
At r 0 , the boundary surface z = −L z has a unit normaln(r 0 ) =ê z . We may construct two orthogonal reference shift vectors {h 1 , h 2 } tangent to the z = −L z boundary plane; that is, both {h 1 , h 2 } satisfy a tangency constraint, h α ·n(r 0 ) = 0 for each α = 1, 2. Without loss of generality, we may choose the {h 1 , h 2 } proportional to the basis vectors {ê x ,ê y }. Written with respect to the global {ê x ,ê y ,ê z } basis, the reference shift vectors are Moreover, the components of each reference shift vector h α are subject to the characteristic scale constraint; for α = 1, 2, where h m is defined in Appendix A. By Equation 12, these reference shift vectors {h 1 , h 2 } are propagated along the congruence to the target footpoint r(L, r 0 ). The propagated vectors at λ = L are also written with respect to the global {ê x ,ê y ,ê z } basis, Since by assumption the vector field B(r) is everywhere smooth, the Spatial Propagator F(λ, r 0 ) is non-singular for all λ, and hence all linearly independent {h 1 , h 2 } are mapped to linearly independent {v 1 (L, r 0 ), v 2 (L, r 0 )}. Furthermore, while we have the freedom to choose the reference shift vectors h α tangent to the boundary surface at the launch footpoint r 0 , in general the propagated vectors v α (L, r 0 ) at the target footpoint r(L, r 0 ) are not tangent to their respective boundary surface. However, we may project the v α (L, r 0 ) onto the tangent boundary plane (see Figure 6) with r = r(L, r 0 ) at the target footpoint. The superscript (t) denotes the projection tangent to the boundary. In our simple example, at the target footpoint, r = r(L, r 0 ), the relation n(r) = −ê z and Equation 53 yield We remark that the construction/projection procedure Equation 53 is valid for any surface with normaln(r) on the system boundary.
The local properties of the mapping of vectors h α tangent to the boundary plane at r 0 , into vectors v (t) α (L, r 0 ) tangent to the boundary plane at r = r(L, r 0 ) are described by the 2 × 2 matrix D; the components of which are The matrix D is the sub-matrix of F(L, r 0 ) obtained by removing the column corresponding to the coordinate direction parallel to the launch-surface normal and row corresponding to the coordinate direction parallel to the target-surface normal; in this example, remove the column (F| L,r 0 ) z j and remove the row (F| L,r 0 ) i z . This principle is the same for general curvilinear coordinate systems, where general boundary-surface unit normals are not parallel to coordinate basis directions at all positions on the boundary surface. Hence the tensorial nature of the Spatial Propagator allows one simply to rotate the proper components of the 3 × 3 matrix into coordinates parallel to the surface normals at the point on appropriate boundaries, and then remove the appropriate rows/columns to construct the proper 2 × 2 sub-matrix; we leave such analysis for future work.
The field lines define a map from h-coordinates with origin at r 0 , to v (t) -coordinates with origin at r = r(L, r 0 ). By Equations 16 and 17, the basis vectors h α are propagated, via Lie transport (Scott, Pontin, and Hornig, 2017), from r 0 to r = r(L, r 0 ). The matrix D may be interpreted as the Jacobian of this map, given by the corresponding components of the matrix representation of the Spatial Propagator (F| L,r 0 ) i j . The squashing factor Q, may be immediately constructed; in the notation of Titov, Hornig, and Démoulin (2002), where the norm N  is given by and is the determinant of the Jacobian matrix D, In general, Equation 58 is the minor of the sub-matrix corresponding to the (F| L,r 0 ) i j matrix component, where the i-index corresponds to the coordinate direction parallel to the launchsurface normal, and the j -index corresponds to the coordinate direction parallel to the targetsurface normal.
Moreover, the eigenvalues {d 1 , d 2 } of the 2 × 2 Jacobian matrix D, given by the roots of the characteristic equation which, after some algebra, reduces to We remark that the line-tying condition does not imply a fixed Spatial Propagator; only the proper sub-matrices and their combinatorics forming the squashing factor Q are preserved.
The Q-Map (Titov, 2007), identifies the separatrix and QSL field structures. However, on their own, these regions offer only the possible current-sheet formation sites in the field structure. Whether or not electromagnetic stresses accumulate resulting in subsequent current-sheet formation depends on the details of energization and stress injection; e.g. separatrix and QSL structures undergoing rigid-body motion will not develop a current sheet (Aulanier, Pariat, and Démoulin, 2005;Aulanier et al., 2010;Janvier et al., 2014); whereas more general motions, such as shearing or twisting motions, will develop a current sheet unstable to reconnection Effenberger et al., 2011;Janvier et al., 2013).

Conclusion and Future Applications
In this article we introduce a generalized field-line connectivity phase space associated with the vector magnetic field in which the geometric and topological features of the system are made explicit. The fundamental assumption is that the vector magnetic field is a priori smooth everywhere. The basic elements are the field line and its linearized variation, the Spatial Propagator: Equation 19. Equations 2 and 16, with initial conditions from Equations 3 and 17, provide a direct formulation of these phase-space elements in terms of the vector magnetic field and its spatial derivatives. Furthermore, the field line and Spatial Propagator are constructed with respect to general curvilinear coordinates and the equivalence class of general affine parameterizations.
The geometric interpretation is that the Spatial Propagator characterizes the organization of the local bundle of field lines. Since the vector field is everywhere smooth, so too are the field-line and Spatial Propagator solutions smooth and unique. The geometric organization of the local bundle is completely equivalent to a kinematic description of a volume centered on the particular field-line solution and undergoing deformation by transport along the field. This deformation kinematics is characterized by volumetric dilation, anisotropic stretch, and rotation. The volumetric dilation (expansion/compression) is completely described by the determinant of the Spatial Propagator: Equation 30. For the vector magnetic field, the determinant of the Spatial Propagator is a topological invariant everywhere equal to unity, Equation 34, which reflects the divergence-free condition. The anisotropic stretch and rotation kinematics are described by the singular values and singular vectors of the Spatial Propagator. The singular values, Equation 37, characterize the general anisotropic stretch of the congruence volume. The congruence rotation is a simple rigid-body rotation kinematics between the orthonormal principal-direction bases: Equation 43, from the reference volume 0 to the propagated volume λ .
Extreme singular values identify QSLs within the system; true separatrix surfaces and separator lines within the system are identified in the limiting cases of one, or two, zero singular values, respectively. Moreover, the Q-factor is simply constructed from analysis of the particular sub-matrix of the Spatial Propagator obtained by removing the column corresponding to the coordinate direction parallel to the launch-surface normal and row corresponding to the coordinate direction parallel to the target-surface normal.
This magnetic connectivity phase-space framework opens up extensive directions in geometric and topological analysis of vector magnetic fields. For example, in future work we will relax the a priori smooth vector magnetic-field assumption in order to analyze both existing singular structures such as current sheets and their formation. Moreover, the magnetic helicity may be decomposed into twist and writhe components (Moffatt and Ricca, 1992); in future efforts we will show that, accounting for relative shearing of the volume, twist helicity may be described with the congruence rotation, while the writhe helicity is related to the embedding of the field-line trajectory and propagation of the reference shift vector in 3D space.
In the present article the field-line bundle and Spatial Propagator are presented and analyzed from a geometric perspective. Since the field lines are curves (i.e. spatial positions) in a three-dimensional space whose tangent vector is everywhere parallel to the magnetic-field vector, evolutionary dynamics (ideal or otherwise) of the field-line bundle are often described by imposing the frozen-in condition (e.g. infinite plasma conductivity) somewhere within the system domain (typically at the system boundary). This allows one to describe and follow the dynamics, not of the field line itself, but rather of the parcel of plasma to which the field line is connected. In non-ideal systems; the difference between the ideal motion and the actual motion of the field is attributed to resistive slipping; this is the origin of "slip-running reconnection" and other resistive slip phenomenology. In follow-up analysis, we use this formalism to describe the dynamics of a field-line bundle and associated Spatial Propagator without reference to plasma conductivity, or material parcels, but rather with respect to pure electric and magnetic fields; such a description necessarily requires a treatment of the full four-dimensional electromagnetic-field tensor.
In particular, we will show that the magnetic field generalizes to the antisymmetric, dual electromagnetic-field tensor * F μν (t, x) (see, e.g., Jackson, 1999, Section 11.9). The field line generalizes to a 4-vector flow field φ μ (λ; t 0 , r 0 ), where the spatial components are identified with the canonical notion of a field line. In addition, the propagator generalizes to a second-order mixed tensor field F μ ν (t, x) where the spatial components are identified with the Spatial Propagator described in this article. Moreover, we will show that expanding these four-dimensional generalizations into explicit spatial and temporal components recovers the definition of a magnetic-field line as three spatial constraint equations, and three time-dynamic equations for the field-line components, where S is the Poynting vector (Jackson, 1999, pp. 608 -610). An explicit spatial and temporal decomposition of the propagator follows similarly. Finally, this framework allows for the analysis of extrinsic thermodynamic properties, such as the total mass or total magnetic energy, of the field-line bundle through Equations 31 and 32. In particular, we derive an analogous first law of thermodynamics applied to the pure magnetic-congruence geometry, where E and B are the electric and magnetic fields, T is the matrix representation of the electromagnetic stress tensor, ρ is the net charge density, and S is the Poynting vector (Jackson, 1999, pp. 608 -610); all quantities are evaluated along the field line r = r(λ, r 0 ). Equation 63 is the total electromagnetic energy of a congruence formulated in terms analogous to mechanical work and energy generation. We will derive Equation 63 and explore applications to coronal phenomenology (e.g. coronal heating, active-region stability, flares, and CME initiation) in future work.
Since the vector field B(r) is assumed smooth everywhere, it may be shown that the L 1,2 norm of the first-and second-order variation matrices F and M are non-singular; that is, there exists a finite F, M ∈ R such that Hence, Equation 66 becomes where we defined the characteristic scale h m ≡ 2F/M; that is, the characteristic scale is the minimum scale such that the linear term is at least of the same magnitude as the secondorder term. Then, for all |h| h m , the higher-order terms of Equation 69 are all negligible relative to the first-order term, and the first-order variation F(λ, r 0 ) describes the relative behavior of the congruence.
In the special case that the congruence solution r(λ, r 0 ) is linear in reference point r 0 , then all components of (M| λ,r 0 ) i jk = 0 and all higher-order derivatives vanish, and h m → ∞. Hence, for all h ∈ R 3 , Furthermore, for completeness, in cases of high symmetry the higher-order variations may vanish through some finite order n. In such examples, the characteristic scale is defined as h m ≡ ( (n+1)!F M (n+1) ) 1/n with n ≥ 2, where M (n+1) is the L 1,2 norm of lowest-order non-zero variation matrix. We do not pursue such cases here.
An extremely popular approach to computational solar and space-plasma physics is the linear interpolation methodologies employed in numerical algorithms for the computation of vector-field components between grid points (see, e.g., Tassev and Sevcheva, 2017). Such approaches simply shift the problem of investigating small-scale dynamics from the physical vector magnetic field on to finer and more complex numerical-grid resolutions. Depending on the size of h m relative to the characteristic scales of the system, this may be a fundamental limitation as an accurate representation of the vector-field structure. The formalism presented here is a fundamentally different approach in which the investigative burden remains on the physical quantities, as opposed to simply increasing the computational expense to handle finer and finer mesh structure.

Appendix B: Matrix Representation of the Covariant Differential of a Vector Field with Respect to General Basis
In this appendix, we briefly review some aspects of canonical differential geometry (e.g. general coordinate charts, basis and dual basis vectors, and the differential operator ∇). In particular, we construct the 3 × 3 matrix representation of the covariant differential ∇B(r) with respect to general, local curvilinear bases. The term covariant is standard nomenclature; in this context meaning the object transforms as a tensor under general coordinate transformations. We construct the matrix representation of the covariant differential ∇B(r) with respect to a general coordinate chart (Appendix B.1), followed by respective applications to the spherical-polar curvilinear coordinate basis (Appendix B.2), and the spherical-polar orthonormal basis (Appendix B.3) in common use by the solar and space-physics community (see, e.g., Tassev and Sevcheva, 2017). We note that this appendix is not intended to be a full exposition of differential geometry. We describe only those aspects necessary and sufficient for the construction of the covariant differential of the vector field in curvilinear coordinates; see Kobayashi and Nomizu (1963) for a full axiomatic approach.

B.1.1 General Curvilinear Coordinates, Basis and Dual Basis Vectors
Let M ⊆ R 3 be a subset of three-dimensional space. Consider a general curvilinear coordinate chart (M, q) such that the position of each point p ∈ M is described by smooth coordinate functions q : M → R 3 , r(p) = q 1 (p), q 2 (p), q 3 (p) . (71) For ease of notation, we drop the formal dependence on the point p ∈ M, labeling it with only the coordinates q i . The coordinate chart (M, q) admits a natural set of coordinate basis vectors, {q 1 , q 2 , q 3 }. The (not necessarily orthonormal) basis vector q i is the unique element of the tangent bundle that acts on the (smooth) position vector by taking its derivative with respect to the coordinate q i , such that the position vector Equation 71 has components written with respect to the {q 1 , q 2 , q 3 } basis, The coordinate basis vectors {q 1 , q 2 , q 3 } are often written explicitly with respect to Cartesian basis set {ê x ,ê y ,ê z }. Assume that the general curvilinear coordinates are smoothly related to Cartesian coordinates; that is, there exist smooth coordinate transformation functions, such that the position vector Equation 71 has components written with respect to the {ê x ,ê y ,ê z } basis Since the Cartesian basis vectorsê i are independent of position, the definition in Equation 72 leads to each basis vector q i expanded in components with respect to the standard Cartesian basis In general, the components of the Jacobian are functions of position: q i = q i (q j ). Moreover, for a given basis q i there exists a unique dual basis q i defined as those linear functionals that satisfy (see, e.g., Kobayashi and Nomizu, 1963, pp. 4 -6 for details). Since the coordinate transformation functions Equation 74 are assumed smooth, that is the Jacobian is locally invertible, the inverse coordinate transformation functions exist. Hence, similar to Equation 76, the dual basis q i may be written in components with respect to the standard Cartesian basis In general, the components of the inverse Jacobian are functions of position, q i = q i (q j ).

B.1.2 The Metric and the Inner-Product
The standard inner-product between two Cartesian vectors u and u is a binary operation defined by u · u = i u i u i , where u i and u i are the components with respect to the standard Cartesian basis. Moreover, the inner-product operation is symmetric [u · u = u · u], non-degenerate [u · u = u · u ] if and only if u = u , and positive definite [u · u ≥ 0] where the equality holds if and only if u = 0, More generally, the inner-product operation between two vectors with components given with respect to a general (curvilinear) basis is defined by a bilinear form g known as the metric (see, e.g., Kobayashi and Nomizu, 1963, p. 24), g(u, u ) ≡ u · u . The metric g is symmetric [g(u, u ) = g(u , u)], non-degenerate [g(u, u  Any vector u may be represented as a linear combination of smooth components with respect to the basis q i : Hence, the inner-product between two vectors with components given with respect to the basis q i is The scalar values g ij = g(q i q j ), are simply the inner-products between basis vectors q i . The g ij may be calculated using Equation 76, By similar arguments, we may construct (inverse) scalars g ij = g(q i , q j ) using Equation 78, In general, the scalars g ij = g ij (q k ) and g ij = g ij (q k ), are functions of position. Furthermore, Equations 81 and 82 satisfy k g ik g kj = k g ik g kj = δ i j .
The scalars g ij and g ij may be represented as 3 × 3 invertible matrices, We remark that it is conventional to define the g ij to be the components of the metric g and the g ij the components of the inverse metric g −1 .
By the non-degeneracy assumption, Equations 76, 78, 81, and 82, the metric and inverse metric provide a natural isomorphism between the q i and q i . By Equations 82 and 76 where in the last step we used the definition in Equation 78. A similar calculation with Equations 81, 78, and 76 yields j g ij q j = q i .

B.1.3 Curvilinear Coordinates and Orthonormal Bases
On the other hand, one may always construct an orthonormal basisq i , respectively dual basisq i , from linear combinations of the curvilinear coordinate basis q i and dual q i basis, The components of the transformation matrices are smooth functions of the coordinates M i j = M i j (q k ), respectively N i j (q k ). In particular, the matrix components M i j and N i j are determined by orthonormality conditions, Here the g kl and g kl are given by Equations 81 and 82. Furthermore, orthogonality between the orthonormal basis and dualq i ·q j = δ i j implies that the transformation matrices are related by N = M −1 .
We remark that by Equations 85 and 86, the coordinate basis q i and dual basis q i are strictly proportional when these basis sets are mutually orthogonal; that is when the metric is diagonal: q i = g ii q i and q i = g ii q i . In particular, for every orthonormal basis set the basis and dual basis vectors are identified,q i =q i ; specifically, the Cartesian basisê i =ê i .

B.1.4 The Covariant Differential Operator ∇
The covariant differential ∇ is a differential operator that satisfies the standard axioms (see, e.g., Kobayashi and Nomizu, 1963, pp. 123 -124). In curvilinear coordinates, ∇ is represented with respect to the dual basis q i via its action on smooth functions of position q i .
Let f : R 3 → R be a differentiable scalar function, then the covariant differential of f is (see, e.g., Kobayashi and Nomizu, 1963, p. 143); in Cartesian coordinates the basis and dual are identified,ê i =ê i , and Equation 89 reduces to the standard gradient. Furthermore, the standard differential is given by where ds = i dq i q i , and by Equation 77, q i · q j ≡ δ i j . Since the local bases q i are themselves smooth functions of position q i , the action of the operator ∇ on the basis vectors is a linear operation, where the linear connection coefficients are functions of position: k ji = k ji (q l ) (see, e.g., Kobayashi and Nomizu, 1963, p. 143).
A vector field B is a linear combination of smooth component functions B i = B i (q j ) with respect to the basis q i : The covariant differential ∇B of the vector field B is then where we used Equations 89 and 91 in the second step, and we relabeled dummy summation indices in the third step. Hence, the i, j th component of the matrix representation of ∇B(r) with respect to the local basis q i ⊗ q i is Equation 21: (see, e.g., Kobayashi and Nomizu, 1963, p. 144). The matrix representation of the covariant differential with respect to general curvilinear bases explicitly includes terms that describe both the rate of change of the vector components, as well as the manner in which the bases vectors change from point-to-point in the system.
The determination of the linear connection coefficients depends on the nature of the basis vectors. In the case of a curvilinear coordinate basis, the linear connection coefficients i jk may be written explicitly in terms of coordinate derivatives of the metric (see, e.g., Kobayashi and Nomizu, 1963, p. 160).
In the case of an orthonormal basis, the linear connection coefficients i jk are given by (see, e.g., Kobayashi and Nomizu, 1963, p. 78 Since, in general, the components and their derivatives of the transformation matrices are functions of position, so also are the structure constants functions of position: c ij k = c ij k (q l ). Finally, in the case of the Cartesian basis, the linear connection coefficients are all trivial, = 0, since both the basis and dual basis are independent of position. Therefore, for vectors with components constructed with respect to a Cartesian basis, the covariant differential reduces to the standard partial derivatives with respect to the Cartesian coordinates (see, e.g. Equation 15).

B.2 Spherical-Polar Coordinate Basis
The position of each point p ∈ M is written in terms of the spherical-polar curvilinear coordinate functions, with r ∈ [0, ∞), θ ∈ [0, π], and φ ∈ [0, 2π). The Cartesian spherical-polar coordinate transformation functions are x = r sin θ cos φ, The position vector is written, Hence, using the definition of the coordinate basis vectors of Equation 72, the sphericalpolar coordinate basis vectors are written in components with respect to the standard Cartesian basis: e r = sin θ cos φê x + sin θ sin φê y + cos θê z , e θ = r cos θ cos φê x + r cos θ sin φê y − r sin θê z , e φ = −r sin θ sin φê x + r sin θ cos φê y .
Using Equation 81 and the fact that the Cartesian basis are orthonormal, the metric components are immediately identified: e r · e r = g rr = 1, e θ · e θ = g θθ = r 2 , e φ · e φ = g φφ = r 2 sin 2 θ.
All other scalar products are zero. Hence, the spherical-polar coordinate basis vectors {e r , e θ , e φ } are orthogonal, but not orthonormal. Furthermore, the inverse metric components follow immediately: Hence, the only non-trivial linear connection components are r θθ = −r, Thus, with respect to the spherical-polar curvilinear basis the vector field is B = B r (r, θ, φ)e r + B θ (r, θ, φ)e θ + B φ (r, θ, φ)e φ ,

B.3 Spherical-Polar Orthonormal Basis
By Equation 87, the spherical-polar coordinate basis Equation 100, and the dual basis Equation 103, may be orthonormalized using the transformation matrices M and N: The spherical-polar orthonromal basis, respectively dual basis, are defined by linear combinations of the spherical-polar coordinate basis {e r , e θ , e φ } Equation 100, and the dual basis {e r , e θ , e φ } Equation 103, e r = e r ,ê θ = re θ ,ê φ = r sin θ e φ .
We note that Equation 116 is independent of the basis in which the matrix F is represented. Furthermore, the index position in which the pair (j, σ i (j )) appears is irrelevant.

C.1 Calculating the Determinant of a Non-Singular, 3 × 3 Matrix F by the Identity Equation 30 in an Arbitrary Basis
Let (M, h i ), with M ⊆ R 3 , be a coordinate chart with an arbitrary basis h i . The coordinate differential basis vectors are dh i = dh i h i . Written with respect to this basis, the signed differential volume element is d (h i ) = dh 3 · dh 1 × dh 2 = dh 3 h 3 · dh 1 h 1 × dh 2 h 2 = h 3 · (h 1 × h 2 ) dh 1 dh 2 dh 3 .
Suppose the linear map F may be represented as a non-singular, 3 × 3 matrix with respect to this basis h i and dual basis h i (see Appendix B.1) Since F is non-singular, its action is completely described by the transformation of the basis vectors. The propagated coordinate differentials follow immediately: The signed differential volume element with respect to the v i basis is then The triple vector product h k · (h i × h j ) reduces the number of terms in the i, j , k summation expansion to six; namely, only those involving the permutations of (i, j, k) = (1, 2, 3) with i = j = k. Hence, we renumerate the reduced i, j , k summation according to elements of the S 3 permutation group, Equation 115 (see, e.g., Kobayashi and Nomizu, 1963, p. 283): where we used Equation 118 in the second step.
(123) Furthermore, the differential volume element Equation 29 propagated along the congruence by the 3 × 3 matrix F is d λ = (F · dz) · (F · dx) × (F · dy) = (dz F ·ê z ) · (dx F ·ê x ) × (dy F ·ê y ) = (F ·ê z ) · (F ·ê x ) × (F ·ê y ) dx dy dz = (F ·ê z ) · (F ·ê x ) × (F ·ê y ) d 0 . (124) The general 3 × 3 matrix components of F with respect to the (Cartesian) basisê i ⊗ê j are The action of the 3 × 3 matrix F on the basis vectors is Then Hence, the triple vector product of the orthonormal Cartesian basis e i under the action of the linear map F is the determinant of the 3 × 3 matrix of F with the usual meaning, Finally, combining Equations 124 and 128, we recover Equation 30.

C.2 Differentiation of the Determinant of a Smooth, Non-Singular, Linear Map F(λ)
In this appendix we demonstrate Equation 33 for the determinant of a square 3 × 3 matrixvalued C n (n ≥ 1) function of a single real parameter: F(λ). Let F, M : R → R 3×3 be smooth, non-singular, linear, 3 × 3 matrix-valued functions of λ ∈ R; that is, the components F ij (λ) and M ij (λ) with respect to any basis representation are at least C 1 -functions. Furthermore, assume F(λ) and M(λ) satisfy the linear matrix ordinary differential equations for all λ ∈ R; in components with respect to any basis Equations 129 and 130 are a generalization of Equations 16 and 17. In particular, M ij (λ) = (∇B| r ) i j with r = r(λ, r 0 ) and initial condition (F 0 ) ij = δ i j . However, the following results are general, independent of the specifics of the basis representation and physical vector fields.