Southward Drift of Eastern Australian Hotspots in the Paleomagnetic Reference Frame Is Consistent With Global True Polar Wander Estimates

Eastern Australian hotspots produce the longest continental tracks on Earth ( > 2,000 km). These hotspots are generally assumed to be stationary with respect to one another, an observation reinforced by our analysis. If any motion between them occurred, it is within our 95 % uncertainties of 370 km of accumulated motion between hotspots. In contrast, motion between eastern Australian hotspots and the spin axis is indicated by changing paleolatitudes of hotspots through time. Reconstructing Australian hotspot paleolatitudes makes use of the global reference frame established by the geocentric axial dipole hypothesis. The classic paleomagnetic approach to establishing paleolatitude typically uses studies of discrete formations of known age, but the nearly continuous record of volcanism over 34 Ma permit using paleomagnetic data from multiple, different studies that overlap in magnetization age in conjunction. This improves the robustness and precision of our paleolatitude estimates for Eastern Australian hotspot volcanoes. A northward shift of hotspot paleolatitudes between 15–24 Ma of ∼300–1,010 km ( 95 % confidence limits) is observed for volcanoes along the Comboyne, Canobolas and Tasmantid tracks, and for interpolated positions along the Cosgrove and Lord Howe tracks but is significantly resolved only for volcanoes along Comboyne and Canobolas tracks, as well as for two ages from the Stradbroke and Britannia seamounts on the Tasmantid track. Before 24 Ma most data plot to the North of a fixed present day hotspot location in the paleomagnetic reference but 95 % confidence limits overlap with it. Notable exceptions are an age from Hillsborough Volcano (Cosgrove Hotspot), and an age from Horse Head (Lord Howe Hotspot) that lie a small distance to the North km, and kmNorth, respectively. Importantly, paleolatitude changes observed along eastern Australian tracks are all explainable by the net rotation of the solid Earth, or true polar wander, with hotspots embedded in the global moving hotspot reference frame.


INTRODUCTION
Hotspot tracks are linear arrays of age-progressive volcanoes that mark the surface terminations of upwelling, hot mantle plumes (Wilson, 1963;Morgan, 1971). At their oldest end, hotspots are sometimes associated with a large igneous province, or flood basalt, presumed to denote the arrival of a plume head (Courtillot et al., 2003). The age progressive track that follows is associated with a plume tail, a structure that anchors the hotspot, the plume tail's surface termination, to lower mantle convection (Steinberger, 2000).
Presently active hotspots (e.g., Hawaii, Iceland) have an active volcano and a topographic high (related to the buoyancy flux of material flow from the mantle) above them, as well as a low shear wave velocity (V s ) region beneath them (Courtillot et al., 2003). Plumes have a finite lifespan, and so these features can be expected to disappear as activity wanes.
Hotspot track geometry reflects the relative motion between the tectonic plate and the hotspot, modified to a lesser degree by lateral magma transport along the lithosphere-asthenosphere boundary (Sleep, 1996;Sleep, 1997). Plate motions can be described by rotations about an Euler pole relative to a reference frame (Cox and Hart, 1986). Euler rotations can be defined by fitting the geometry and age progressions of a pair (or more) of hotspot tracks underlying a given plate (e.g., Hawaii and Louisville on the Pacific Plate), defining plate motion relative to the mantle (absolute plate motion), anchoring relative rotation models to a common reference. Commonly Euler rotations are defined relative to a neighboring plate (this is the case when reconstructing seafloor isochrons) and, by combining rotations for multiple plate pairs, a hierarchy can be constructed to define a relative plate rotation model. The observation that hotspot tracks get older in the direction of plate motion led to the early hypothesis that hotspots were stationary in the mantle (Morgan, 1971), and as such, might provide a means for defining a mantle reference frame, so linking surface tectonic processes to the deep Earth. But fixed Pacific hotspot reference frames failed to predict the age progression and geometry of those in the Indo-Atlantic realms (and vice versa) (Molnar and Atwater, 1973;Duncan, 1981;Molnar and Stock, 1987;Müller et al., 1993;Wessel and Kroenke, 2008). Uncertainties in the relative rotations between plates (forming plate circuits), used to translate hotspot geometries from one plate to another were originally hypothesized to account for these misfits, with possible Eocene extension between East and West Antarctica being the supposed cause of misfit between Pacific and Indo-Atlantic fixed hotspot reference frames. Failing to account for all the degrees of freedom when fitting rotations to hotspot tracks also contributed (Andrews et al., 2006) Sharp changes in hotspot track geometry, such as the Hawaiian-Emperor bend have been reproduced by Indo-Atlantic fixed hotspot reference frames by taking this circuit (Torsvik et al., 2017), but relative motion between hotspots is still indicated (the predicted length of the Emperor bend being ∼800 km shorter than observed) (Torsvik et al., 2019).
Relative motion between hotspots, inferred from the misfit between different, fixed hotspot reference frames, indicates that a combination of hotspot motion, non-rigid plate behavior, and errors in plate motions is required to explain them (Wang et al., 2017). Although the rate of relative motion between hotspots is estimated to be quite slow ( 4.5±3.8 mm a −1 ) (Wang et al., 2017), a new generation of hotspot reference frames emerged, based upon plume conduits embedded within whole-mantle convection cells (Steinberger and O'Connell, 1998;Steinberger, 2000;Torsvik et al., 2008;Doubrovine et al., 2012). Using a plate model for surface boundary conditions, this innovation allowed for hotspot motion, accommodating the motion inferred between earlier, fixed reference frames, while also defining a mantle reference of global extent. Such models predicted that hotspots are anchored by their plume tails to lower mantle return flow, moving away from subduction zones, and toward mid-ocean ridges (Steinberger, 2000).
Paleomagnetism has also proven to be a useful tool for examining hotspot motion, yielding paleolatitude estimates for extinct volcanoes along hotspot tracks. Paleomagnetic data provide latitude constraints with respect to Earth's geocentric axial dipole (symmetry prevents reconstruction of paleolongitude). If hotspots are fixed in the mantle, and true polar wander (TPW), that is, the rotation of the entire solid Earth (both lithosphere and mantle) with respect to Earth's axis of rotation were negligible, then paleolatitudes for extinct hotspot volcanoes should be invariant, and coincide with the hotspot's present day position. On the other hand, if hotspots were moving in latitude, then the paleolatitudes along their tracks should reflect this history. Paleolatitudes for both Hawaii and Kerguelen indicate past motion of these hotspots with respect to the spin axis (Grommé and Vine, 1972;Kono, 1980;Gordon and Cape, 1981;Cox and Gordon, 1984;Tarduno and Cottrell, 1997;Antretter et al., 2002;Tarduno et al., 2003;Doubrovine and Tarduno, 2004;Tarduno, 2007;Tarduno et al., 2009). As such, paleolatitude estimates reflect the combination of any hotspot motion together with TPW. Correcting for TPW involves subtracting absolute plate motion (in the mantle reference frame) from a continent's apparent polar wander path (APWP), that describes its motion with respect to the rotation axis. Conventionally this is done by rotating paleomagnetic data from the continents to a common reference frame, constructing a global apparent polar wander path (e.g., Torsvik et al., 2012), and then subtracting absolute plate motion away. That said, TPW corrected paleolatitudes for Hawaii and Kerguelen are incompatible with being caused entirely by TPW (Antretter et al., 2002;Torsvik et al., 2017). As such, a component of hotspot drift (with respect to the mantle) has been inferred from these results. Ridge-plume interactions have also been suggested to account for latitude changes at the time these hotspots were passing over or nearby ridges (Site 1140 for Kerguelen and Detroit Seamount for Hawaii, Antretter et al., 2002;Tarduno et al., 2009), suggesting that lithosphereasthenosphere topography might play some role in augmenting hotspot volcano distribution.
Paleolatitudes for both of these oceanic hotspots were calculated from inclination-only drillcore of seamounts, by determining a virtual geomagnetic latitude (λ VGL ) for each flow mean inclination (I m ) (Doubrovine and Tarduno, 2004), and then solving for the maximum-likelihood estimate of virtual geomagnetic latitudes (McFadden and Reid, 1982;Arason and Levi, 2010). The precision of such an approach is limited by the number of flows sampled, and, as samples are all derived from the same location, also might be influenced by local, non-dipole contributions (Torsvik et al., 2017). In contrast, continental hotspot tracks offer several advantages over their oceanic counterparts. They are relatively easier to obtain samples from, while full orientation information affords complete treatment of magnetization vectors, in contrast to inclination only data. This last factor is a major advantage, as contemporaneous results from other rocks can be used in reconstruction, affording both greater precision, as well as being less prone to any local, non-dipole effects. Any position on a continent can be reconstructed to its ancient paleolatitude given a corresponding paleomagnetic pole, by computing the paleocolatitude (p), the angle between pole and site positions (Butler, 1992). In this case the uncertainty on the pole (A 95 ) can be used to compute the uncertainty on the paleolatitude. Given coordinates of an extinct hotspot volcano and a paleomagnetic pole corresponding to the age the volcano was active, its paleolatitude at that time can be found.
Intraplate volcanism is widespread in Eastern Australia during the Cenozoic (Figure 1) and provides an ideal natural laboratory to study hotspot motion. Three types of magmatic provinces are classified onshore: 1) central volcanoes are complexes with differentiated magmatic products dominated by early basaltic flows and late-stage, felsic plugs intruding central vents; 2) leucitite suite edifices are low-volume features formed by distinctive potassium-rich, silica-poor lavas; and 3) lava fields are characterized by thin, effusive flows, scoria cones and maars (Wellman and McDougall, 1974). The location and age of lava fields follow no obvious trends , and are thought to result from edge-driven convection and shear-induced decompression melting over localized regions of thinner lithosphere (Conrad et al., 2011;Davies and Rawlinson, 2014;Aivazpourporgou et al., 2015;Rawlinson et al., 2017). In contrast, the central volcanoes and leucitite magmatic provinces define clear age-progressive tracks that reflect the Cenozoic motion of Australia over several hotspot anomalies. These tracks provide an important link between hotspots in the Pacific and Indo-Atlantic, and, together with the fact that they have not been used to determine hotspot reference frames, mean that their paleolatitudes provide a useful test case for evaluating the global moving hotspot reference frame (GMHRF), global apparent polar wander pathl (APWP) and their associated TPW estimates Torsvik et al., 2012).
From the roughly North to South trajectory of eastern Australian hotspot tracks, it is generally suggested that the underlying hotspots moved little with respect to one another, FIGURE 1 | Distribution of eastern Australian Cenozoic Volcanism, hotspot tracks and lithospheric thickness. Central volcanic fields are indicated in black, lava field provinces in gray (Wellman and McDougall, 1974). Hotspot tracks are indicated by dashed gray lines, with ages (in Ma) indicated (Cohen et al., 2007;Cohen et al., 2008;Knesel et al., 2008;Sutherland et al., 2012;Cohen et al., 2013). The red circle indicates a 100 km diameter melt region for the Cosgrove hotspot above its present day position as inferred from tomography (Montelli et al., 2006 and that volcano termination accurately reflects the location of the underlying mantle plume (McDougall and Duncan, 1988;Knesel et al., 2008;Cohen et al., 2013;Crossingham et al., 2017;Mortimer et al., 2018). By applying linear, or segmented linear regression to fit volcano latitude with age along individual tracks, these authors derived rates for Australian plate motion from the slope, although strictly speaking, these estimates represent the relative motion between each hotspot and the Australian plate. Longitude is ignored in these analyses, and although multiple hotspots on the same plate occur, tests for relative motion between these hotspots have not been performed. If the fixed hotspot assumptions are valid, then hotspot tracks will be concentric small circles of one another, as would be expected for plate motion alone over a set of stationary hotspots (Cox and Hart, 1986). As such, the angular distance between pairs of hotspots should remain constant with time (Harada and Hamano, 2000;Wessel et al., 2006). Mantle modeling (dataset S1, Doubrovine et al., 2012) suggests eastern Australian hotspots moved relatively little with respect to one another and the GMHRF ( ≤ 2 + ). Over the ∼35 Ma of eastern Australian hotspot activity, this might manifest as much as ∼ 350 km of relative motion between hotspot pairs, but actual tests of relative motion between these hotspots have not heretofore been carried out.
The Cosgrove track comprises the longest ( ∼2,250 km) and westernmost of Australia's hotspot tracks, beginning with the formation of Hillsborough Volcano ( 34-33 Ma) in central Queensland, followed by Nebo, Peak Range, Springsure and Buckland Volcanoes along a ∼ 450 km SSW trend (Cohen et al., 2013;Crossingham et al., 2018b). Volcanism along this track resumes south of a ∼ 650 km magmatic gap with the eruption of leucitite magmas ( 18-10 Ma) along a younger, ∼ 650 km segment Davies et al., 2015). The magmatic gap is attributed to suppression of partial melting below thick lithosphere, while reduced volumes of magma production, and a shift in lava composition that characterize the southern segment are attributed to diminished melt fractions produced under intermediate thickness lithosphere (Davies et al., 2015). Volcanoes of the Newer Volcanic Province ( ∼2.5 Ma) lie to the South of the Cosgrove leucitite along the predicted path of the hotspot, although the nature of their source is less clear. The relatively large E-W extent of this region has conventionally been used to argue against a plume related source, although Rawlinson et al. (2017) suggest that entrainment of plume material in an existing edge/shear driven upwelling region might explain its extent, and temporal coincidence with the Cosgrove hotspot. The present day position of this hotspot is estimated at 41 + S, 146 + E with tomographic methods (Montelli et al., 2006) although a clear, deep plume tail is not observed (French and Romanowicz, 2015).
Along the coast of eastern Australia the Comboyne hotspot track defines a SSW trend bookended by Fraser Island ( ∼ 31 Ma) and Comboyne volcanoes ( ∼ 16 Ma) (Ashley et al., 1995;Cohen et al., 2007;Knesel et al., 2008). A peculiarity of this track is the relatively long-lived, and higher volume locus of volcanism at the Main Range-Tweed provinces, which is associated with eastward offset from the otherwise linear trend of volcanoes . It, as well as deflections in the Tasmantid and Lord Howe seamount tracks has (under the assumption of hotspot fixity) been used to suggest a transient slowdown of the Australian plate between 26-23 Ma. Linking of these hotspot track distributions with unequivocal collision on Australia's margin has proved difficult to reconcile with geologic evidence. Originally the slowdown was used to infer the timing of collision between the Ontong-Java Plateau and the Melanesian arc, halting Pacific Plate subduction under the Australian Plate, and leading to Australian Plate subduction along SW Solomon Island trench (subduction polarity reversal) . Such a model implies Oligocene shutdown of Melanesian Arc volcanism, together with arc deformation and terrane exchange of the Solomon Islands, from the Australian Plate, to the Pacific. Although reduced in volume, volcanism continued in the Solomon Islands through Oligocene to Miocene times (Hackman, 1980;Chivas, 1982), seismicity and swath mapping indicate that subduction along the Melanesian Arc continued through this interval (Petterson et al., 1999), paleomagnetic poles do not indicate terrane exchange of the Solomon Block from Australia to Pacific during the Oligo-Miocene (Musgrave, 1990;Musgrave, 2013), and deformation of the arc began in Pliocene times (Petterson et al., 1997;Petterson et al., 1999). Alternatively, slowdown may have been associated with collision in Papua New Guinea (Jones et al., 2017), although the collisional history of this region is complex, and making an unambiguous link therefore is challenging.
Two remaining segments of volcanoes also occur onshore, from Bunya-Toowoomba (24-22 Ma), and Nandewar-Canobolas (19-12 Ma), that have been grouped together as the Canobolas hotspot track (Davies et al., 2015), although a track joining these volcanoes would have a large westward offset between Toowoomba and Nandewar volcanos, not observed in either tracks onshore, or the two tracks offshore. The close temporal and spatial association between Bunya-Toowoomba and Main Range-Tweed volcanoes suggest these Bunya-Toowoomba volcanoes might be associated with the Comboyne hotspot, rather than a separate hotspot. The lack of younger volcanoes South of Comboyne and Canobolas volcanoes, together with decreasing estimated volumes of magma at progressively younger volcanoes (Jones and Verdel, 2015) is good evidence that these hotspots are now extinct, their plume tails presumably shearing away once their flux decreased below a critical value (Steinberger and O'Connell, 1998).
In the Tasman Sea are the Tasmantid ( ∼33.5−6.5 Ma) , and Lord Howe ( ∼ 28 − 4 Ma) (McDougall et al., 1981;Mortimer et al., 2018) seamount chains. Courtillot et al. (2003) estimate the present day position of Tasmantid hotspot to be (39 + S, 150 + E) and Lord Howe to be (33 + S, 159 + E) although these are not included in the Montelli et al. (2006) hotspot catalog. The youngest dated volcanoes along these chains are both ∼ 6.5 Ma, at Gascoyne Seamount (Tasmantid) and Lord Howe Island (Lord Howe), and lie to the North of the predicted present day positions, although several smaller seamounts exist to the South of each of these estimates, that might be younger, related volcanoes. An unnamed seamount, possibly part of the Tasmantid hotspot track occurs at 40.1 + S, 157.1 + E, and the Flinders and Heemskerk/Zeehaen seamounts at 34.6 + S, 159.8 + E and 36.3 + S, 159.8 + E, respectively, might be part of the Lord Howe hotspot track. Without dating of rocks at these locations it is unclear if these hotspots remain active today, or went extinct in the Pliocene.

Hotspot Geometries
We performed geometric analyses on hotspot tracks in eastern Australia to evaluate predicted present day positions for them, and to explore the possibility of relative motion between eastern Australian hotspots. To estimate present day positions of hotspots we determined motion paths based upon Doubrovine et al. (2012) plate model, anchored to well dated volcanoes along each track.
Available age data along each track are influenced by geological preservation, and so are not uniformly sampled with time. As such, it is required to approximate track geometry and for this we use spherical splines (Thompson and Clark, 1981). The smoothing parameter controls the trade-off between smoothness of the spline curve, and accuracy of the approximation. Smaller values will more accurately fit the data, but the spline will have greater curvature (smoothing of one will interpolate between points), while greater values will reduce the curvature, by allowing departure of the curve from data. We follow the criteria of Doubrovine et al. (2012) in evaluating appropriate smoothing parameter 1) track geometry should be accurately reproduced, 2) ages along curves should increase monotonically, and 3) the distance between predicted and observed track locations should not exceed 150 km. Angular distances of less than this could conceivably be the result of lateral magma transport between the hotspot melt region and its surface termination (Sleep, 1996;Sleep, 1997). The lack of ages for many of the seamounts along the Lord Howe track is problematic, and we interpolate ages for undated seamounts along this track using our linear latitude model for this track. To test for relative motion between hotspots, we compared angular distance between pairs of tracks using spline fits. If angular distances remain constant with time, we can be confident that particular pair of hotspots have remained stationary with respect to one another. Conversely, changes in this distance over time indicate relative motion between the two hotspots compared.

Paleomagnetic Poles From VGP Level Data
Reconstruction of paleolatitudes using the pole approach requires a robust set of poles over the time interval concerned. For reference, we compare paleomagnetic poles presented in Hansma and Tohver (2019), as this work presents a significant number of new Australian paleomagnetic results, together with a comprehensive review of existing Cenozoic paleomagnetic data. While these poles do justice in reconciling paleomagnetic data with the speading history between Antarctica and Australia, we take advantage of the commonplace overlapping magnetization ages between Australian paleomagnetic studies to systematically determine combined paleomagnetic poles for the Cenozoic interval, accounting for positional uncertainty (A 95 ) as well as uncertainties in the age of paleomagnetic data by randomizing VGP ages within their bounds and sampling VGPs with replacement to generate bootstrap pole samples (Tauxe et al., 1991). This allows for a more robust characterization of Australia's APWP during the Cenozoic.
Our reasoning for this approach is justified because although PSV timesacles range from 10 4 to 10 5 years (approaching 5 Ma at extremes), limited geological preservation routinely makes adequate sampling impractical, if not impossible at this temporal resolution, with ∼ 100 sites required for confident sampling (Tauxe et al., 2003). APWPs are commonly fit to study-based paleomagnetic poles, so two averaging steps are applied, firstly VGPs within each study are averaged to their study based pole position, then these poles are averaged again. Multiple averaging is known to introduce bias, and unlike averaging VGPs, where PSV is being sampled, there is not a clear phenomonological basis for averaging between poles. The averaging intervals commonly applied during this stage are 10 Ma or larger in duration, larger than PSV timescales and, usually, also larger than the intervals of time represented by VGPs from each study. A slightly better approach often used is to weight each pole by the number of sites, but transitional VGPs are often dealt with inconsistently between studies. Others have grouped VGPs from multiple studies into combined poles (a concept dating to at least as early as McElhinny et al., 1974), and the notion that grouping of VGPs based on geology is arbitrary, is not new. This approach becomes problematic for the case of three studies (A, B, and C) each a collection of VGPs: pole A overlaps (in time) with B, and B with C, but A and C do not overlap. How should VGPs from these three studies be combined? A possible reason why, in spite of the growing number of temporal overlaps between paleomagnetic studies, that the issue of double averaging in paleomagnetism has received less attention, is that the identification of suspect data (Van der Voo, 1990), and improvement in sampling and laboratory methods (e.g., Kirschvink, 1980) was of larger concern in the effort to define reliable APWPs deeper in Paleozoic time.
We provide a synthesis of the extant body of site-level paleomagnetic data for Cenozoic Australia (557 sites total, excluding undated sites from laterites) from studies where basic demagnetization techniques have been performed. Magnetization ages have been updated using modern geochronoloic and biostratigraphic constraints (Figure 2; Supplementary Table S1). Multiple studies of coeval rocks from different formations allow us to both validate the consistency of results between studies, and also develop a more robust technique for the determination of paleomagnetic poles, compared to traditional grouping of data by individual geologic formation, and subsequent averaging of these formationbased poles. Our new approach returns the calculation of a paleomagnetic pole to site-level statistics, using VGPs of equivalent age, incorporated from potentially multiple coeval geologic formations. We argue that this approach is superior to existing methods because it provides 1) poles for specific age intervals based on the largest possible dataset, 2) a more robust averaging of PSV during intervals constrained by multiple Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 544496 studies, 3) consistent identification of transitional VGPs, and 4) systematic treatment of VGP age uncertainties. A further advantage is that this method reduces to a traditional, study based average in the case of paleomagnetic data that do not overlap in time.
We compile the directional vector of individual VGPs for each site, together with upper and lower age bounds for the magnetization. Where multiple, comprehensive studies of a province have been conducted, e.g., replication studies of Barrington Volcano (Wellman et al., 1969;Hill et al., 2002), and the Newer Volcanic Province (Rahman, 1971;Opdyke and Musgrave, 2004), we exclude the legacy study to maintain VGP independence, and prevent duplication. We then perform common mean tests between contemporaneous VGPs from different studies formations to ascertain whether the results across coeval intervals are self-consistent. This validation process is described in detail in Hansma and Tohver (2018), Hansma and Tohver (2019), but we note here that blanket demagnetized sites from the Older Volcanics (Mumme, 1963) reject common mean tests with multiple coeval studies, and are for this reason excluded. In some cases common mean testing is rejected marginally [e.g., Oligocene results described in Section 3.5.5 of Hansma and Tohver (2019)], but these all fall within one circular standard deviation of each other. It is plausible that VGPs from these studies, taken individually, each represent a slightly biased sampling of PSV, manifesting itself as slight differences between the means of each study. Rather than selecting a subset here, it is more robust to include all of these results in analysis. Of the 557 VPGs that pass selection, 220 are blanket demagnetized (Q ≥ 4 for these studies), although most of these magnetization directions are supported by more modern sites that sampled rock from the same region and employed stepwise demagnetization. Only 90 VGPs, specifically from Main Range, Liverpool Volcano and the Glasshouse Mountains are not directly corroborated by more modern demagnetization on the same rocks. Coeval studies employing stepwise demagnetization for all these time intervals exist, and as described above, common mean tests do not identify these results as suspect.
By first performing a common means test for contemporaneous VGP populations, we provide a better a priori test for suspect results, giving context to "quality filters" that are generally applied to individual formation-based poles (Van der Voo, 1990;Besse and Courtillot, 2002). Filtering based on a fixed quality score threshold is known to result in false positives and negatives, a caution that was both originally stressed (Van der Voo, 1990) and later reiterated (Van der Voo, 1993). Contextualization of results in the manner described above, followed by examination of their quality scores where differences are identified, does better justice to existing data. It can help minimize the false exclusion of lower score, but otherwise apparently reliable data, while also preventing the inclusion of any suspect data that quality schema fail to recognize.
To calculate poles, we apply a sliding-window approach after determining an optimum window size for the dataset. Two processes contribute to the observed distribution of VGP data with time: 1) short term variation as a consequence of PSV, and 2) FIGURE 2 | Magnetization age bounds for Australian Cenozoic paleomagnetic data. Number of sites indicated in parenthesis. Studies are italicized where limited sampling indicate that secular variation has not been averaged within the region itself.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 544496 long term variation caused by APW. There is a natural trade-off between averaging PSV, which is more reliably achieved with larger intervals, and minimizing the averaging of APW, where the smaller intervals are desirable. To identify the optimum interval we systematically vary the window size from a narrow interval, to progressively larger and larger ones (here from 1 to 20 Ma). At narrow intervals PSV is unlikely to be adequately sampled; as such the variance in VGP scatter across the dataset is expected to be large. As window size is increased, this variance will stabilize as PSV is sampled. Similarly, the median VGP scatter across all intervals is examined for different window sizes and is expected to be unstable when these are narrow, stabilize as PSV is sampled, and then gradually increase as the interval size grows beyond the ideal interval, due to the averaging of increasing amounts of APW together with PSV. For our dataset, window sizes below 10 Ma produce poles with large ranges of S, indicating that for some intervals of the Cenozoic Australian VGP dataset, such narrow window sizes do not adequately sample PSV (Figure 3). Window sizes ∼ 10 Ma and greater have less variation in VGP scatter, but the median VGP scatter gradually increases, a consequence of averaging significant amounts of APW due to plate motion, together with PSV. As such, we applied a 10 Ma sliding window, with 1 Ma steps over the data. Our results here suggest that it is likely that the nature of paleomagnetic sampling through time, dictated in large part by incomplete geological preservation, plays a larger role in limiting paleomagnetic resolution than shorter timescale PSV.
For reference, the procedure for determining poles is depicted in Figure 4, which we will describe in detail here. We calculate a Fisher mean pole where ≥ 50 sites are present, and where VGPs within the interval have an age range of > 3.5 Ma, and less than the window size (10 Ma). We identify transitional sites systematically with a variable VGP cut-off angle (Vandamme, 1994). VGP ages within each window are not usually normally distributed, so we assign the median age, rather than the center of the interval, or the mean VGP age to the pole calculated. We identify age outliers where VGP age is greater than 1.5 times the inter-quartile range, but this only excludes one datum, the 9.8 Ma Cosgrove leucitite (Hansma and Tohver, 2018), that occurs isolated in time from other Australian VGPs. Without any thresholding, for either minimum number of sites or minimum intervals of time, edge-effects associated with gaps and sparse intervals in the VGP record adversely affect the result near their boundaries. In our implementation, minimum time and site thresholding can be applied to the result after the fact (by setting these values to zero) but will increase runtime as the VGPcutoff angle (Vandamme, 1994) is determined iteratively and is called more often. The effect of gaps in the VGP record are managed in three ways: 1) if a window is incremented in time, but no new VGP data are found, then this pole is suppressed, because it is equivalent to narrowing the window size, 2) a minimum number of sites prevents pole calculation in sparsely sampled intervals, and 3) a minimum range for VGPs within each interval prevents pole calculation in intervals where, although enough sites are present, they are concentrated into a narrow time range (i.e., the effective PSV sampling interval is not met). Where the sliding window transitions over well represented, to poorly represented data almost completely, this last threshold is triggered.
For the Australian Cenozoic dataset, the early Eocene represents a sparsely sampled interval, with a few sites distributed within it. The minimum sites threshold prevents poles from being calculated here. The minimum time interval threshold was chosen to prevent the calculation of poles where the ages of VGPs are concentrated into a narrow range, far below than the 10 Ma suggested by VGP scatter analysis described earlier. In our dataset, this affects one search interval, between 34 and 44 Ma, where VGPs are concentrated in the 2 Ma between 34 and 36 Ma, and there are few VGPs in the remaining 8 Ma (36-44 Ma). Such narrow intervals of VGP averaging are unsatisfactory given VGP scatter analysis (Figure 3) indicating these poles are not robust. Interquartile ranges for sites within this interval were all below 1.5 Ma, with ranges of < 3.5 Ma. The pole longitudes for this interval lie at ∼100 + E, 20 + away from poles in the previous interval and the implied rotation for Australia is implausibly rapid. Examining VGP scatter for the poles found in this interval reveals they have significantly lower scatter S 13.5 + than poles in the interval before (33-43 Ma), as well as for older poles (although the older poles are Paleocene in age). This magnitude of VGP scatter is also lower than expected values based on geomagnetic field models (Tauxe and Kent, 2004), given that Australia during the early Oligocene was at middle to high latitudes >35 + S. There is uncertainty associated with the individual age of every VGP, regardless of accuracy and precision of the geochronological assay of the rock unit. This introduces a source of uncertainty in the set of poles that lie within a given interval, as VGPs with magnetization ages that straddle its boundary may, or may not fall within it. Although this is conventionally ignored in paleomagnetic analysis, we account for this by randomizing the set of VGP ages within their age bounds, and sampling with replacement (1,000 iterations, Supplementary Table S2). When reconstructing a feature over a given interval of time, the 95% confidence region surrounding bootstrapped poles with median ages within the interval is used.

Paleolatitude Reconstruction
Paleolatitudes were calculated using distance between corresponding paleomagnetic poles and sites for both for hotspot volcanoes (where hotspot magmas terminated at the surface and have been dated with geochronologic methods), as well as for the estimated position of each hotspot over time, using aforementioned spline fits (Section 2.1). This was done to quantify hotspot position with respect to the spin axis through time for each hotspot. Individual volcano estimates are based upon available geochronologic data presented in existing literature, following preferred ages of original authors (McDougall et al., 1981;McDougall and Duncan, 1988;Ashley et al., 1995;Cohen et al., 2007;McQueen et al., 2007;Cohen et al., FIGURE 4 | Control flow diagram (left) describing method used to generate bootstrapped poles and an accompanying example (right) for the 20-30 Ma search window. t start and t stop represent the age range desired for the entire search (the Cenozoic in this case), t max being the window size. Thresholding rules n min (minimum number of sites) and t min (minimum age range for VGPs in a given interval) can be specified. Thresholding can, of course, be applied to pole results after the fact (by setting n min and t min to zero), but results in longer runtime for reasons described in text. Age distributions within each interval are evaluated first, and thresholding rules applied (in the example histogram these fall beyond the search domain as VGPs are well represented over the interval). A bootstrap paleomagnetic pole is calculated, both including and excluding transitional sites (the example showing the variable VGP cutoff-angle in gray) (Vandamme, 1994). The window is then incremented, and the process repeated for multiple trials until the search is complete, leading to a collection of bootrapped poles from which the pole position and uncertainty can be estimated.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 544496 2008; Mortimer et al., 2010;Sutherland et al., 2012;Cohen et al., 2013;Crossingham et al., 2017;Mortimer et al., 2018;Crossingham et al., 2018a;Crossingham et al., 2018b). 40 Ar/ 39 Ar results, where available, were preferred, but K-Ar ages for Lord Howe Island (McDougall et al., 1981) were included following Mortimer et al. (2018). Several (20 total) age analyses were present in figures, but absent in tabulated data, and were digitized from figures: Cohen et al. (2008, Figure 7) includes ages for Nandewar, Warrumbungle, Canobolas, and Macedon volcanoes, while Mortimer et al. (2018, Figure 7) shows an age for Gifford's Guyot. The 2σ values in these cases were estimated from data point size, this value generally being larger than quoted errors for tabulated age data within these papers, and is 0.6 Ma for those from Cohen et al. (2008) and 1.0 Ma for Gifford's Guyot. For each sample, paleomagnetic poles with median ages within the 2σ uncertainties of age analyses were used for reconstruction. Hotspot locations that might not precisely coincide with volcano terminations were estimated using aforementioned spline fitting.
There are two processes that can affect observed paleolatitudes for hotspots, either their independent motion in the mantle, or the changing moment of inertia of the solid Earth (relative rotation between the solid Earth and the spin axis). Paleolatitudes for the present day position of hotspots were corrected for TPW using the preferred synthetic global apparent polar wander path (and its uncertainties) presented in Torsvik et al. (2012) together with the plate rotations relative to the GMHRF presented in Doubrovine et al. (2012), following the procedure described in Torsvik et al. (2017).

Geometric Analysis
The present day position for each hotspot was predicted by defining a motion path anchored to a well dated volcano at the age of activity along each track (Table 1) using the plate model of Doubrovine et al. (2012). The motion paths for each hotspot are depicted in ( Figure 5). The estimate for Cosgrove lies close to the position inferred by Montelli et al. (2006), being half a degree farther North in latitude and the same longitude. The estimate for Tasmantid hotspot is remarkably close to that estimated by Courtillot et al. (2003), but the position given by these authors for the latitude of Lord Howe hotspot corresponds to Lord Howe Island, where volcanic rocks are dated at ∼ 6 Ma (McDougall et al., 1981). Our estimate is farther South with a comparable longitude (35.7 + S, 159.4 + E), and lies close to the undated Heemskerk/Zeehaen seamounts.
Spline fits ( Figure 5) were used to interpolate along hotspot tracks without invoking parametric assumptions inherent to linear models. The spline fit provides the best estimate for the position of the hotspot at any given time because individual volcanoes may not terminate precisely above their hotspot, and volcanism along each track is not always continuous. With one exception, misfits between spline fits and volcanoes are less than our cut-off distance. Only one datum from Comboyne volcano (sample Q83 from Cohen et al., 2007) has a larger misfit (178 km), but other analyses from this same area are suitably fit. The paucity of ages along the Lord Howe hotspot track meant that several seamount ages were interpolated using a linear age-latitude model and the spline fit using these seamounts. The angular distance between predicted hotspot location at any given time and their volcano termination is normally distributed away from the predicted hotspot center, with 2σ 90 km ( Figure 6).
The angular distance between spline fits was used to test for relative motion between hotspot pairs during time intervals where both are active (Figure 7). Cosgrove and Comboyne hotspots appear most constant, but most hotspot pairs accumulate changes in angular distance over time. The cumulative effect of this is generally small, and within the TABLE 1 | Estimated present day positions of eastern Australian hotspots derived from motion paths , and from hotspot catalogs (Courtillot et al., 2003;Montelli et al., 2006)  FIGURE 5 | Spline fits to eastern Australian hotspots based on geochronologic data (squares, Ashley et al., 1995;McDougall et al., 1981;McDougall and Duncan, 1988;Cohen et al., 2007;McQueen et al., 2007;Cohen et al., 2008;Knesel et al., 2008;Mortimer et al., 2010;Cohen et al., 2013;Crossingham et al., 2017;Mortimer et al., 2018;Crossingham et al., 2018a;Crossingham et al., 2018b;Sutherland et al., 2012). Present day positions are shown in red, based upon motion paths (in gray) derived from Doubrovine et al. (2012) (Table 1).
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 544496 uncertainty of spline misfits for both hotspots (±368 km). The largest angular distance change was between Tasmantid and Lord Howe hotspot tracks (382 km), but ages were interpolated along the track and little significance is placed on this result. The second largest was between Cosgrove and Tasmantid tracks (361 km), but this was still within the uncertainties of the spline fit. While consistent with estimates of the relative motion of eastern Australian hotspots based upon mantle flow modeling , auxiliary dataset 1) it is not necessary to invoke movement between Australia hotspots to explain the distribution of volcanoes along each of their tracks.
Age vs. latitude trends of Eastern Australian hotspot volcanoes have been used to infer changes in Australian plate speeds during the Oligo-Miocene Cohen et al., 2013), but this interpretation is complicated by misfits between eastern Australian hotspot track geometries (Figure 7). Plate rate changes implied in these studies occur over narrow time intervals and in our analysis equate to only small misfits between volcano termination and predicted position ( Figure 6). While they may possibly represent changes in plate speed, they could also reasonably be explained by small amounts ( ∼ 90 km) of lateral magma transport away from the hotspot center.

Cenozoic Poles for Australia
Our combined paleomagnetic poles are in general agreement with APWPs inferred from sea-floor spreading between Australia and Antarctica during the Cenozoic, as well as conventional APWP analysis (Hansma and Tohver, 2019) (Figure 8). Our new approach provides a more continuous record of Australian paleomagnetic poles for intervals that are well-represented by VGPs from penecontemporaneous rock formations. Confidence limits on our poles are more conservative, due to the inclusion of VGP age uncertainties. Pole latitudes change with time at a faster rate after Eocene times, consistent with an increase in Australia-Antarctica spreading rate, from slower NE-SW spreading to faster N-S directed spreading during the Eocene (there is a paucity of Eocene VGP data and more detailed comparison is not possible) (Whittaker et al., 2007;Whittaker et al., 2013). The change in spreading direction, from NE-SW to N-S, cannot currently be   resolved from Australian Paleocene paleomagnetic data. For convenience, Table 2 gives a summary of poles within each search window, while Supplementary Table S2 gives results for all bootstrap samples. In the 0-10 Ma interval, paleomagnetic data mostly come from the Newer Volcanic Province (Opdyke and Musgrave, 2004), and here our pole calculation reduces to the pole position for these VGPs. This demonstrates that where VGPs from different studies do not overlap continuously in time our method reduces to a conventional pole calculation.
VGP scatter for bootstrapped poles show a decrease from Paleocene times to present day (Figure 9). Australia has moved from high latitudes to lower ones throughout the Cenozoic, and a decrease in VGP scatter is consistent with geomagnetic field models and recent lavas (Tauxe and Kent, 2004). Including transitional VGPs increases the scatter because scatter is sensitive to the presence of outliers. Recent lavas from low latitudes have VGPs that are less scattered than those from higher latitudes, and it is remarkable that such an observation can be observed in aggregate through time for Australia.
The proportion of transitional sites identified by variable VGP cut-off is possibly, but not clearly related to reversal frequency (Ogg, 2012) (Figure 10). Such a correlation would be expected if transitional sites are chance recordings of magnetic field reversals.
Correlation of their normalized values between 15 and 28 Ma, as well as between 55 and 60 Ma is encouraging, but transitional directions are recovered less often than would be expected before 5 Ma, and more often between 28 and 34 Ma. There are numerous reasons for imperfect correlation between these quantities. Transitional VGPs could be over-or under-represented by sampling. Transitional VGPs could reflect magnetic excursions. Other factors that might produce sites with VGPs apparently distant from their means may also play a role (e.g., unrecognized local rotation, local lightning remagnetization, or for the case of blanket demagnetized data, incompletely removed overprints). The interval compared does not contain superchrons, and does not sample intervals of low reversal rate, and as such does not sample the entire range observed for the polarity history of Earth. In spite of this, the correlations between reversal rate and proportion of transitional sites recovered are exciting and warrant further study.

Paleolatitudes for Australian Hotspots
Despite the relatively large uncertainties inherent in paleomagnetic data and the small effect size expected given TPW manifest since 35 Ma, paleolatitude estimates for eastern Australian hotspots are incompatible with hotspots that are fixed with respect to Earth's spin axis, but within uncertainty are all compatible with  Table 2). Discontinuities in poles through time occur over intervals where VGP data are insufficient for pole calculation (e.g., between 2-12 Ma). (B) Poles from A. plotted with their A 95 ellipses (in gray, darker gray where hotspot volcanoes are active) plotted together with synthetic APWP paths, indicated by squares (Seton et al., 2012) and triangles (Matthews et al., 2016) as well as conventional APWP construction (ages labeled in Ma) (Hansma and Tohver, 2019 hotspots that are embedded in the convecting deep mantle ( Figure 11). Estimates plotted in Figure 11 are tabulated in Supplementary Material. Contemporary estimates of TPW presented here Torsvik et al., 2012) have a resolution of 10 Ma, and even at our finer resolution, our results agree remarkably well with these estimates. Within uncertainties, differentiating between TPW models is not possible. Eastern Australian hotspots lie between the Pacific (Hawaii and Louisville) and Atlantic (New England, Reunion and Tristan) hotspots that constitute the basis for the GMHRF. Paleolatitude estimates from eastern Australian hotspots provide important validation of this reference frame for Oligocene and younger times. A northward shift of hotspot paleolatitudes between 15 and 24 Ma is observed for both volcanoes and spline interpolated positions along the Comboyne, Canobolas, and Tasmantid tracks, and for interpolated positions along the Cosgrove and Lord Howe tracks where volcano or seamount age data are absent ( Table 2). To illustrate the magnitude of this, at 23.2 Ma the Comboyne hotspot paleolatitude is −36.7 + ± 3.2 + (300-1,010 km North of the present day position, were it fixed to the spin axis). While this trend is significant for all interpolated positions along hotspot tracks, it is only significantly resolved for age data from volcanoes along Comboyne and Canobolas tracks (where there are numerous ages), as well as for two ages from the Stradbroke and Britannia seamounts on the Tasmantid track.
Before 24 Ma most data plot to the North of a fixed present day hotspot location in the paleomagnetic reference but 95% confidence limits overlap with it. Exceptions are found for an age from Hillsborough Volcano (Sample BC-98; Cohen et al.,  (Vandamme, 1994), with medians given by blue markers. Gray (unfiltered values) are also shown. Marker transparency of 1% per iteration was used. A trend of decreasing VGP scatter with time is consistent with northward motion for Australia, and falls within expected values based on field models (Tauxe and Kent, 2004) over this time interval.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 544496 FIGURE 10 | Number of transitional sites identified from Australian continental paleomagnetic data (top), and these values normalized, plotted with normalized reversal rate per Ma for 10 Ma intervals (below) (Ogg, 2012).
FIGURE 11 | Paleolatitudes for the (A) Cosgrove, (B) Comboyne, Canobolas, (C) Tasmantid, and (D) Lord Howe hotspots. Red (and green) data indicate paleolatitude of volcanic rock localities with geochronologic constraint (with 2σ age errors, and 95% confidence on paleolatitudes) along each hotspot track. The blue line represents paleolatitudes for spline fits with shaded 95% confidence region. Tabulated values are given in Supplementary Material. The black dashed line is a fixed point, in the paleomagnetic reference frame, at the present day position of each hotspot. The gray line is the TPW corrected position of this point through time (that is, fixed in the mantle reference frame) with associated shaded 95% confidence limits Torsvik et al., 2012).
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 544496 2013) along the Cosgrove track. It has a paleolatitude estimate of −37.0 + ± 3.2 + , 22-743 km North of the present day estimated position. If the (Montelli et al., 2006) hotspot estimate is used, this distance is 33 km greater. Paleolatitudes for the northern Peak Range (29-30 Ma) also lie North of these positions when using conventional paleomagnetic poles for Australia (Hansma and Tohver, 2019), but overlap in our analysis as pole uncertainties are more conservative, and so larger. Although there is a paucity of age data along the Lord Howe hotspot track ( Figure 11D), analytical ages as well as our interpolated ages follow a similar trend to the other hotspots. While most analytical ages along the Lord Howe hotspot overlap at 95% confidence with both paleomagnetic fixed, and TPW corrected values, an age from Horse Head (Mortimer et al., 2018; Sample DR15) along the Lord Howe track gives a paleolatitude of −30.5 + ± 3.7 + , 167 − 988 km North of the present day estimated position were it fixed with respect to the spin axis. As with other eastern Australian hotspots, within uncertainty, this motion is compatible with a hotspot fixed in the GMHRF.

CONCLUSION
The distribution of hotspot volcanoes along tracks in eastern Australia, combined with Australia's Cenozoic paleomagnetic data indicate that hotspots moved little with respect to one another and were not fixed to the geocentric axial dipole. That the changes in hotspot paleolatitude with respect to the paleomagnetic reference frame are consistent across all of the eastern Australian hotspots is remarkable, especially so as the motion observed corresponds with the predicted amount of true polar wander over the interval for which these hotspots are active . The results from eastern Australian hotspots, lie between the Pacific and Indo-Atlantic hotspots used to define this prediction, and so provide important validation for the GMHRF and TPW estimate for Oligocene and younger times.

DATA AVAILABILITY STATEMENT
The datasets analyzed for this study are included in the manuscript/Supplementary Material.

AUTHOR CONTRIBUTIONS
ET and JH both worked on the original idea of reconstructing eastern Australian hotspot paleolatitudes. JH developed the combined pole approach and wrote the Python implementation. ET and JH both wrote and edited the manuscript. We wish to thank both reviewers for their insightful comments that significantly improved the final version of this manuscript.

FUNDING
The University of Western Australia supported JH with a Robert and Maude Gledden Scholarship and was supported by an Australian Government Research Training Program (RTP) Scholarship. ET acknowledges Visiting Professor Fellowship from the Fundação CAPES.