A data-driven framework for paleomagnetic Euler pole analysis

Owing to the inherent axial symmetry of the Earth’s magnetic field, paleomagnetic data only directly record the latitudinal and azimuthal positions of crustal blocks in the past, but paleolongitude cannot be constrained. An ability to overcome this obstacle is fundamental to paleogeographic reconstruction. The paleomagnetic Euler pole (PEP) analysis presents a unique means to recover such information in deep-time. However, prior applications of the PEP method have invariably incorporated subjective decisions into its execution, undercutting its fidelity and rigor. Here we present a data-driven approach to PEP analysis that addresses some of these deficiencies—namely the objective identification of change-points and small-circle arcs that together approximate an apparent polar wander path. We elaborate on our novel methodology and conduct some experiments with synthetic data to demonstrate its performance. We furthermore present implementations of our methods both as adaptable, stand-alone scripts and as a streamlined interactive workflow that can be operated through a web browser.


Introduction
Following from Euler's rotation theorem (or fixed point theorem), the motions of lithospheric plates across the surface of the sphere can be precisely described in the form of finite rotations-commonly called Euler vectors (or Euler poles)-which can be concisely expressed by three parameters: the latitude and longitude of a pole of rotation (θ, ϕ) and an angular displacement ψ. In plate tectonic and paleogeographic research, great effort is dedicated toward unraveling and chronicling the history of plate kinematics, which is tantamount to the determination and collation of Euler vectors. Improvements and innovations in the methods of Euler vector inference are therefore of broad importance.
An assumption common to all methods of Euler vector inference is that the motion of a given plate with respect to some arbitrary reference can be treated as stable over some nominal interval, such that it can be expressed by a single 'stage' Euler vector (Gordon et al., 1984;Cox and Hart, 1986). From the perspective of the system of reference, this stage Euler vector will be fixed and any point belonging to the plate will move along a plane perpendicular of that axis of rotation, tracing a circle centered on it. If an object passes from the domain of the reference to the plate (becoming a passive occupant of the latter), it too will follow the path of a circle during the stage interval, with a displacement proportional to its residence time on the plate. A series of such objects, passed from the reference domain to the plate progressively in time during the stage interval, would therefore form an arcuate array decorating the trace of a circle about the Euler vector. If a sufficient number of such objects could be identified, their distribution could be inverted to retrieve the location of the Euler vector, and their corresponding residence times on the plate could be used to determine its amplitude ψ.
At least two kinds of observational records present such age-progressive arcuate segments that can be used to infer stage Euler vectors: hotspot tracks and paleomagnetic apparent polar wander paths (APWPs) (Figure 1). Hotspot tracks record the motion of a plate with respect to the underlying mantle, whereas APWPs record the motion of a plate relative to the Earth's magnetic field, which is aligned with the planetary rotation axis. The inversion of paleomagnetic data to retrieve stage Euler vectorsor paleomagnetic Euler pole (PEP) analysis-is of particular interest because it offers the possibility to recover full kinematic descriptions of past plate motion, in contrast to the conventional analysis of paleomagnetic data, which constrains only paleolatitude and paleoazimuth. However, despite being first conceptually introduced more than half a century ago (Francheteau and Sclater, 1969), PEP analysis has seen limited application (e.g. Smirnov and Tarduno, 2010;Wu et al., 2017;Swanson-Hysell et al., 2019;Beck and Housen, 2003;Gordon et al., 1984). This appears to be due, at least in part, to the fact that most prior approaches to PEP analysis have incorporated subjective decisions into its execution. In order to secure the reproducibility of scientific results, our objective here is to introduce an unsupervised and objective framework with which to conduct PEP analysis-although our methodology is applicable to stage Euler vector inference more generally. We furthermore seek to provide an accessible implementation of this framework for community use by the use of Jupyter notebooks (Kluyver et al., 2016) that can be easily deployed using Binder (Project Jupyter et al., 2018). the Euler poles (EP; black poles) describing the kinematics of the host plate (grey polygon with static boundaries for simplicity) which gave rise to them over some interval of time (t 0 -t 10 ). A sequence of paleomagnetic poles and/or hotspot volcanoes may be inverted to discover the location of the corresponding EP.

Paleomagnetic Euler pole analysis
If a given plate's motion relative to the planetary spin-axis is constant for some interval of time, the paleomagnetic poles aquired by the plate over that interval will trace an arc across the surface of the sphere. The fundamental task in PEP analysis is to invert an assemblage of paleomagnetic data to retrieve the average instantaneous kinematics that were acting during such an interval (i.e. the stage Euler vector). However, because plate motions intermittently change according to force balance adjustments, APWPs comprise a discrete series of arcuate segments separated by change-points whose number and timing is unknown. Thus, PEP analysis also necessitates the division of paleomagnetic data into a number of temporal subsets-or groups-such all the paleomagnetic poles of a given group represent the same (unknown) stage Euler vector. We are thus presented with two tasks that we may treat as distinct optimization problems: 1. Grouping paleomagnetic data into subsets that manifest distinct stage Euler vectors. 2. Determination of the best-fit stage Euler vector for each paleomagnetic subset.
Here we develop an unsupervised search process to solve these problems in the reverse order by first determining the best-fit stage Euler vector to all possible APWP segments (all ordered subsets of the paleomagnetic data), and then identify the least-cost sequence of Euler vectors approximating the total APWP.

Methodology
We start by considering that we have a series of N paleomagnetic poles that we call p i , whit i = 1, 2, . . . , N. Each one of these paleomagnetic poles have an associated date t i and can be represented in Cartesian coordinates (x i , y i , z i )-with x 2 i + y 2 i + z 2 i = 1and also by their latitude and longitude, (θ i , ϕ i ). Without loss of generality, we assume that the paleomagnetic poles are time ordered, meaning t i < t j for i < j. We refer to P = {p 1 , p 2 , . . . , p N } as the full sequence of paleomagnetic poles that defines the APWP. Furthermore, given two indices i and j with i < j, we denote by P i− j = {p i , p i+1 , . . . , p j } the subsequence of paleomagnetic points contained in the interval of time between t i and t j . We use q to refer to Euler poles, where q i is the instantaneous Euler pole at time t i . Lastly, we denote with R(q, ω) the matrix describing a rigid counterclockwise rotation around the axis q by the angle ω, which thus represents the average instantaneous rotation over the time interval-or stage-i.e. a stage Euler vector.

Gradient-based Euler vector optimization from paleomagnetic data
In this section we consider the problem of fitting an arc (whether belonging to a small-or great-circle) to a subsequence of paleomagnetic poles. Given k and m with k < m, we want to find the coordinates of the Euler pole q and an angle φ ∈ [0, π/2] such that the points in P k−m are approximated by the trace of a circle with angular radius φ, centered on Euler pole q (Figure 2a).
We infer the parameters of the circle by minimizing the sum of the squares of the geodesic distance from each data point to the circle trace. The geodesic distance can be computed as the absolute difference of the angle between the Euler pole q and each with (x q , y q , z q ) the coordinates of the Euler pole q. Then, we consider an optimization problem that seeks to minimize the cost function given by where the optimization is perform with respect to the three parameters of the circle that we write as the three-dimensional vector Θ = (θ q , ϕ q , φ). In the noiseless case where all the points lie along a small-or great-circle and satisfy the equation xx q + yy q + zz q = cos φ for all p ∈ P k−m , the quantity C k−m (Θ) achieves its minimum at zero. Since each (x i , y i , z i ) is typically affected by random perturbations, there is no exact fit between the set of points and a small or great circle trace.
In contrast to the ordinary least squares method, this optimization problem has no analytical solutions. Although an optimal solution can be found numerically via grid searching (e.g. Gordon et al., 1984), such an approach converges slowly. This approach was also considered in Gray et al. (1980) and Fujiki and Akaho (2009), where the authors use a constrained optimization method to find the Cartesian coordinates of the pole. Here we employ the nonlinear conjugate gradient method of Polak and Ribiere (1969) to find an optimal solution more efficiently. We take advantage of the solver provided in the open-source SciPy package (Virtanen et al., 2020) to perform the computations. It is worth noting that even though we have the following constraints θ q ∈ [0, π], ϕ q ∈ [0, 2π) and φ ∈ [0, π), the periodicity of the cosine and sine functions allow us to treat this problem as an optimization problem without constraints.
The final update yields a single Euler pole which best describes the subsequence of paleomagnetic poles ( Figure 2b). Its cost function gives the sum of the residuals between the predicted and observed paleomagnetic poles, and thus, a measure of goodness of fit. Using the spherical law of cosines, from the geographic coordinates of the inverted Euler pole (θ q , ϕ q ) we can infer a stage rotation angle ψ, that is, the rotation angle needed to translate the youngest paleomagnetic pole p k to the oldest one p m , at an average angular velocity ω = ψ/(t m−1 − t k ). To illustrate the performance of our methodology, an example Jupyter notebook with synthetic data is provided via Binder (https://mybinder.org/v2/gh/LenGallo/PEPy/HEAD).

Identifying the least-cost sequence of Euler vectors approximating an APWP
Equipped with a method to invert a given set of paleomagnetic data to recover the best-fitting Euler vector describing them, we now present a method that seeks to identify the unknown change-points (i.e. plate-motion changes) that divide an APWP into discrete segments associated with distinct stage Euler vectors. This is accomplished by seeking that sequence of change-points and Euler vectors which maximizes the goodness of fit between the predicted and observed paleomagnetic poles.
Our goal is to identify a subsequence of ordered indices (i.e. change-points) j 1 , j 2 , . . . , j m+1 , with j 1 = 1 and j m+1 = N + 1, such that each one of the m segments P j k − j k+1 belongs to the same circle. This implies that q (i) = q j 1 = . . . = q j i+1 −1 and where we define by q (i) and φ (i) the Euler pole and radius of the circle in the ith segment-or stage-of the APWP. In order to do so, we seek to minimize the following cost function where the optimization is performed over all possible increasing paths ( j 1 , j 2 , . . . , j + m + 1) of unknown length m with j 1 = 1, j m+1 = N + 1; and with respect to the Euler pole q (k) and rotation angle φ (k) associated to each one of the m segments P j k − j k+1 . The last term λm in the cost function (3) corresponds to a regularization added to avoid overfitting, with λ ≥ 0 a hyperparameter of the model. Large values of λ will increase the value of the cost function, such that a good-fitting APWP with fewer segments will be favored. Note that equation (3) can be rewritten as with Θ (k) the parameters of the minimum circle for the segment P j k − j k+1 . Based on the method introduced in section 3.1, we can find the Euler pole and amplitude for each one of the possible segments P i− j for each possible value of i and j, and then compute the minimum inside the parentheses in equation (4).
Our approach to minimize the cost function of equation (4) proceeds with an exploratory data analysis to characterize all possible pathways P i− j , and then employs a graph representation of the APWP to discover the optimal sequence of change-points ( j 1 , j 2 , . . . , j m+1 ) and the number m which yield the minimum sum of squared misfits between the set of paleomagnetic poles and the modelled arc segments describing them.
Graphs are compact mathematical structures that amount to a set of connected objects, where objects are referred to as nodes and connections are referred to as arrows (or edges for undirected graphs). Arrows can have an associated weight that indicates the cost of traversing them. An APWP can then be considered as a network where the nodes correspond to paleomagnetic poles and pairs of nodes are connected by arrows, to which a numerical weight can be assigned (Figure 3). We may thus model an APWP P = {p 1 , p 2 , . . . , p N } where the arrow from p i to p j exists for all i, j such that t j > t i and at least 3 paleomagnetic poles (in order to define a plane, at least 3 non-colinear points are required). The weight of each arrow is given by the optimal value of the function C i− j plus λ. In order to be consistent with our previous notation, the contribution of the last paleomagnetic pole p N is included in the weight of each arrow with its tail in p N , that is, the weight of each arrow from p i to p N is given by the goodness of fit of the sequence P i−(N+1) = {p i , . . . , p N } plus λ. In this context, a candidate approximation of the APWP can be seen as the directed path (p j 1 , p j 2 , . . . , p j m+1 ) with p j 1 = p 1 and p j m+1 = p N , where each pair of consecutive nodes in the path (p j k , p j k+1 ) is connected by an arrow. The total cost of each alternative series of directed paths is given by the sum of the individual weights of the arrows in the path, which coincides with equations (3) and (4). The characterization of all possible pathways P i− j is accomplished by a moving window searching algorithm, where each possible P i− j is analyzed.
In graph or network theory, Dijkstra's algorithm allows identification of the shortest path between two given nodes (Dijkstra, 1959), such that the sum of the weight of its constituent edges is minimized. In an ideal (noiseless) case, the sum of the weight of its constituent edges connecting the youngest node and the oldest one, is null. We perform shortest path searches through the open-source package SciPy (scipy.sparse.csgraph.shortest_path, Virtanen et al., 2020).
The hyperparameter λ is adjusted to improve the accuracy of the model to avoid overfitting, namely, by preventing an overestimation of the number of segments, m. We apply the 'Elbow' method to find the λ value that yields the optimum value of m (Bholowalia and Kumar, 2014). This is achieved by plotting the total cost against the number of segments m and evaluating their tradeoffs: (too) small values of m will add large total costs but at some threshold value of m the cost will be observed to decrease sharply. Above this threshold, further increasing m will not yield significant further reductions to the total cost. The graphical 'elbow' in the plot thus points to this optimal value of m.
We summarize our algorithm to find the optimal sequence of Euler vectors to approximate a given sequence of paleomagnetic poles P as follows: 1. For each time-ordered pair of paleomagnetic poles p i and p j , find the circle in P i− j = {p i , . . . , p j } that minimizes the value of C i− j plus λ. 2. Construct a directed graph and select λ through the Elbow method to avoid overfitting. 3. Use Dijkstra's algorithm to find the shortest path between the first pole p 1 and the last pole p N . 4. From the optimal path, we recover the parameters of the stage Euler vectors for each APWP segment.

Results
To illustrate the application of our methodology, we execute a simple test involving idealized synthetic paleomagnetic data. We start by generating a stochastic model of the drift of a plate over 200 Myr by randomly generating stage Euler vectors-represented in the matrix form as R(q, ω)-which remain constant during a temporal stage of given duration. The geographic coordinates of each Euler vector (θ q , ϕ q ) are randomly selected from a uniform distribution on the unit sphere and their angular velocities ω are chosen from a uniform distribution constrained between 0.5 • /Myr and 2.5 • /Myr (which yields plate velocities up to 28 cm/yr). The corresponding duration of each Euler vector is randomly selected from a uniform distribution between 10 and 30 Myr, and Euler vectors are drawn until their summed durations reach or exceed 200 Myr. Using this stochastic kinematic model, we assemble a synthetic APWP populated by paleomagnetic poles every 1 Myr (Figure 4a). To replicate a more realistic scenario with noise, we perturb the geographic coordinates of the paleomagnetic poles by Gaussian random noise up to 0.1 • . Note that the latter step is only taken to demonstrate the capacity of our method to handle noise; a more comprehensive consideration of the feasibility of PEP given standard paleomagnetic noise has yet to be addressed.
Next, we proceed to invert the synthetic paleomagnetic data according to the methodology outlined above. First, for each candidate APWP segment (all ordered subsets of the paleomagnetic data) we compute the best-fitting Euler vector by minimizing equation (4). We then organize the computed costs as a graph in the form of an adjacency matrix, which can be rendered visually as a heatmap (Figure 4b). Visual inspection of this heatmap reveals strong contrasts in the magnitude of the cost function, which tends to increase sharply with the passage of a change-point. Using the regularization introduced in equation (4) and the Elbow method (Figure 4c), we identify the shortest path through this graph following Dijkstra's algorithm, which yields the number of segments m and a discrete set of change-points ( j 1 , j 2 , . . . , j m+1 ). From this sequence of change-points, we retrieve the parameters of the minimum circle Θ (k) for each segment P j k − j k+1 , enabling us to infer the complete kinematic history through the series of stage Euler vectorsR(q, ω). The workflow can be run interactively in the browser via Binder (https://mybinder.org/v2/gh/LenGallo/PEPy/HEAD) with a synthetic example in the Jupyter Notebook (APWP.ipynb).
To measure the performance of our graph-based PEP methodology in terms of accuracy, we generate 1000 stochastic models of the drift of plate and invert the paleomagnetic data from each for the series of underlying stage Euler vectorsR(q, ω). Comparisons between the pairs of true R(q, ω) and estimatedR(q, ω) stage Euler vector series are then assessed by a variety of metrics. Using the inverted kinematic model, we assemble a synthetic APWP and compute the time-dependent geodesic distance (i.e. residual) between it and the known APWP. In Figure 5a we show the ensemble median of this residual, along with the 0.25-75 and 0.16-0.84 percentile ranges representing standard confidence intervals. As could be expected, we observe that average misfits tend to accumulate backward in time, but we note that they never exceed 3 degrees. In 5b-c we further dissect the discrepancies: in panel b we show the ensemble median of the time-depedendent geodetic distance between Euler poles, while panel c shows the difference in angular velocities. These results reveals that the average misfit between the true and estimated stage Euler vectors of any given simulation is low. Furthermore, we compare the true kinematics and our inferences quantifying the discrepancy of the drift of a given point after rotation by the series R andR. We start by generating 5 nodes evenly distributed in latitude (in the central meridian) and we then calculate a time-dependant geodesic distance between the inferred drift and the true drift for each From this stochastic kinematic model, we assemble a synthetic APWP (dots) populated by paleomagnetic poles every 1 My. b) Heatmap visualization of the adjacency matrix illustrating the magnitede of the cost function for each candidate APWP segment. c) Scatter plot illustrating the Elbow method: the total cost (y axis) is plotted against the number of segments m (x axis), color coded by the hyperparameter λ. The graphical 'elbow' in the plot thus points to this optimal value of λ. point. Last, we compute the mean of these distances for each time in order to illustrate the time-dependant drift (i.e. total motion) misfit. Figure 5d shows that the median of the drift misfit never exceeds 4 degrees. Finally, we address the velocity surface misfit associated with our kinematic inferences. With that aim, the surface velocity at each of the former points is calculated through the cross-product between each one of the former points and the true and inferred Euler vectors.
Considering all these ensemble medians to be representative of the misfit associated with the workflow in the absence of significant noise, we conclude from these metrics that our methodology can robustly retrieve the kinematics of a plate from its APWP.

Discussion
Conceptually, in providing a means to retrieve the full kinematic histories of tectonic plates from paleomagnetic data, PEP analysis presents an enormously powerful tool. And yet, despite its foundations being laid more than half a century ago (Francheteau and Sclater, 1969), lingering doubts about the rigor and viability of PEP analysis have left it largely stranded in obscurity.
One of the most significant and persistent shortcomings in PEP analysis has been the subjective determination of which APWP segments to invert, and which inversions to retain. As noted by previous studies (e.g. Tarling and Abdeldayem, 1996;Wu and Kravchinsky, 2014), the optimal fits to real paleomagnetic data tend toward short segments associated with proximal circles of small-radii, which together often imply improbable rates of plate motion change, unrealistic angular velocities and/or geologically dubious kinematics. As a consequence, the selection of which segments to fit has been done a priori by visual inspection (e.g. Gordon et al., 1984;May and Butler, 1986;Wu and Kravchinsky, 2014), undercutting the rigor and reproducibility of the analysis. In some instances the best-fitting small-circle solutions have also been dismissed in favor of either subjective assessments of more feasible solutions (Tarling and Abdeldayem, 1996) or forced great-circle fits (Smirnov and Tarduno, 2010). Although Smirnov and Tarduno (2010) consider the use of great-circle fits to be a conservative measure, we note that approach confines all stage Euler vectors to the equatorial plane. Rose (2016) presented a significant step forward in introducing a Bayesian approach to PEP analysis (later applied by Swanson-Hysell et al. (2019)), which addressed some of these issues of subjectivity and presented a pathway toward a more comprehensive treatment of uncertainty in such analyses. However, in its present formulation, their approach still regularize the problem through the choice of the number of change-points.
Our methodology addresses several of these deficiencies, and notably the problem of objectively selecting change-points, by identifying the least-cost solution without supervision. The least-cost solution is a full approximation of the APWP, such that the location of the change-points and the parameters of the stage Euler vectors defined between them are fully and automatically determined by the algorithm. Although the method employs a regularization parameter to avoid overfitting, the Elbow method allows a more objective selection of this parameter. Other parameters that can be modified by the user (number of inversion seeds, maximum duration of a given stage)-while capable of affecting the outcome if set to too low values-are used to improve efficiency and so can be set to arbitrarily high values. Admittedly, our experiments here have been conducted on idealized synthetic data, and it can be anticipated that experiments with real (noisy) data may introduce circumstances in which the least cost solution can otherwise be shown to be unrealistic (e.g. requiring impossibly high rates of rotation). However, the flexibility of our methodology is such that minor adaptions to our optimization framework (e.g. specifying constraints to rotation rates in the optimization problem) allow such issues to also be tackled systematically.
Despite the advancements presented herein, there remain a number of significant challenges associated with the application of PEP analysis to real paleomagnetic data. The first is the feasibility of the methodology in light of the noise accompanying paleomagnetic data, arising from both intrinsic (geomagnetic secular variation) and extrinsic (e.g. age uncertainties, tectonic rotations, inclination shallowing) uncertainties. Although the cost function can be further constrained to avoid unrealistic solutions, there exists some level of noise above which the uncertainty in the solution space is so vast that the optimal solution is rendered meaningless. However, a systematic exploration of these noise thresholds, as well as the trade-offs between noise level and the rate of plate motion, have yet to be conducted.
A second challenge concerns the phenomenon of true polar wander (TPW). TPW is a rotation of the entire solid Earth which occurs in response to changes in the planetary moment of inertia following changes to Earth's internal distribution of mass (Goldreich and Toomre, 1969). Because TPW always acts to restore the largest inertial axis to the planetary axis of rotation, TPW always occurs about an equatorial axis. Thus, during a TPW event, the entire tectonic plate system (together with the mantle) will move relative to the spin axis by a common angular displacement. Paleomagnetic data, in recording the movement of the lithosphere relative to the magnetic field (which remains parallel to the spin-axis during TPW) therefore represents some composite signal of TPW and individual plate motion. Before PEP analysis can be applied to any given paleomagnetic dataset, the TPW signal needs to be removed from it. Unfortunately, time-dependent estimates of TPW are themselves uncertain and often controversial. More significantly, prior to Cretaceous time, TPW estimates are derived from plate kinematic models (Steinberger and Torsvik, 2008); any results from PEP analysis employing these TPW estimates would thus undermine their own foundation wherever they challenged the underlying kinematic model. Future progress on this front will likely necessitate a joint modelling approach.
We end with a reminder that the methodology and tools presented herein, while motivated by the desire to improve upon existing approaches in PEP analysis, are not restricted to paleomagnetic applications. They could, for example, be equally applied to hotspot tracks. There are also broader applications of the small-circle fitting operation (which is intentionally presented as a standalone tool), as arise, for example, in structural geology (supplementary materials figure?).

Conclusions
Owing to the inherent axial symmetry of the Earth's magnetic field, paleomagnetic data only directly record the latitudinal and azimuthal positions of crustal blocks in the past, but paleolongitude cannot be constrained. The paleomagnetic Euler pole concept presents a unique means by which to estimate motions in deep-time. We have presented a methodology that enables PEP analysis to be conducted in a fully data-driven way. Our methodology follows the solution of two distinct optimization problems. We first find the parameters of the best-fitting circle to a set of paleomagnetic poles through a gradient-based Euler vector optimization method. Then, through an exploratory data analysis, all the possible segments at every point of the APWP are assessed in order to transform the data into a weighted, undirected graph. This graph representation of the APWP allows the identification of that sequence of stage Euler vectors which together present the least cost (shortest path of the graph) description of an observed APWP. We demonstrated the performance of our algorithm with a set of simple experiments with synthetic data drawn from stochastic kinematic histories, from which our inversions recovered kinematic parameters that closely resembled the known values. We closed with a discussion of some challenges that remain for the application of PEP analysis to real paleomagnetic data, and a note that the tools presented herein have applications beyond PEP analysis.

Acknowledgments
This project has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 101025975 (L.C. Gallo) and from the Research Council of Norway through the Centres of Excellence funding scheme, project number 22327. The simulations were performed on resources provided by UNINETT Sigma2 -the National Infrastructure for High Performance Computing and Data Storage in Norway.