Paleomagnetic Evidence for Inverse Correspondence between the Relative Contribution of the Axial Dipole Field and CMB Heat Flux for the Past 270 Myr

We provide an evaluation of the paleolatitudinal dependence of the paleosecular variation throughout the Paleozoic-Mesozoic transition – linked to the high geomagnetic reversal frequency interval Illawarra Hyperzone of Mixed Polarity (IHMP; ~266.7–228.7 Myr). Our findings were compared with those for intervals of distinctive geomagnetic reversal frequencies within the Phanerozoic. Our results for the IHMP were conducted through estimates of angular dispersion (SB) of virtual geomagnetic pole (VGP) data groups, taken from a high quality paleomagnetic database. Model G was fitted to these data, providing its shape parameters a and b (respectively related to the antisymmetric and symmetric harmonic terms for the time-average geomagnetic field). Results for the IHMP exhibited compatible patterns with two well-known intervals of higher reversal frequency – Jurassic and the last 5 Myr. A comparison of b/a ratio results – considered as an efficient indicator for the relative contribution of the axial dipole field – for the last 270 Myr, indicated an inverse correspondence with the relative core-mantle boundary (CMB) heat flux, according to recent discussions, clarifying the physical meaning of the Model G shape parameters a and b.

a solution for a field generated by a spherical geodynamo -in relation to the non-axial dipole contribution. Such conditions have been linked to 'superchrons' (~10 7 yr, single geomagnetic polarity periods), as discussed by Biggin et al. (ref. 25 ) for the CNS, and for the 262-318 Myr Permian-Carboniferous Reversed Superchron (PCRS; ref. 26 ). Conversely, a lower dipolar contribution was reported for intervals of higher reversal frequency, such as the Jurassic 25,27 and the last 5 Ma 28 .
Such information can be acquired by evaluations of the ancient geomagnetic field through analyses of paleosecular variation (PSV), related to the spatio-temporal variability in both direction and intensity of the EMF 8,22 . It provides an independent way of investigating the EMF evolution through geological time, hence it is adequate for assessing information on the time-averaged field, and its dipolar and non-dipolar contributors 4,25,29,30 . The PSV is commonly obtained by the angular dispersion (S) of virtual geomagnetic poles (VGPs) datasets, given by: where N and Δ i are, respectively, the number of VGPs and the angular difference between the ith VGP and the mean VGP. A phenomenological model that has been successfully employed for evaluation of S -which demonstrated a clear relation between reversal frequency and the latitudinal dependence of VGP dispersions 24 , was proposed by McFadden et al. (ref. 31 ). This approach (Model G) considers that the VGP angular dispersion results from the contribution of two independent "families" -dipole (S D ) and quadrupole (S Q ) families, which are respectively related to odd and even l-m spherical harmonic terms (i.e., asymmetric and symmetric around the equator region): where λ is the paleolatitude, and a and b are the Model G shape parameters (which are empirical constants that are respectively related to the quadrupole (symmetric) and dipole (antisymmetric) families of the field). From hemispherically averaged VGP dispersion datasets carried out from 0-5 Ma lava flows, McFadden et al. (ref. 27 ) reported a possible correspondence for the past 160 Myr between the reversal frequency and the ratio b/a -which provides an empirical evaluation of the relative contribution of antisymmetric (b) to symmetric (a) harmonics terms of the geodynamo. Furthermore, Coe and Glatzmaier (ref. 24 ) reported by means of modeling simulations of the geodynamo that the symmetry of the time-averaged field -which can also be indicated by the ratio b/a -can be a better predictor of reversal frequency in comparison to the intensity evaluations.
Nevertheless, some important questions are still far from being completely elucidated about the extension of the large-scale variations for the reversal frequency, and its connections to the CMB heat flux fluctuations (linked to the long-term mantle dynamics) throughout the Phanerozoic. For instance, there are no reported discussions so far for: (i) a possible lower contribution of the antisymmetric family for the high reversal rate interval known as Illawarra Hyperzone of Mixed Polarity (IHMP; ~266.7-228.7 Myr). The IHMP is characterized by a high mean geomagnetic reversal frequency (comprising tens of polarity reversal events from the end of PCRS (Late Permian) to the lowermost Triassic [32][33][34], and is possibly related to some of the prominent geodynamic events that took place during the Paleozoic-Mesozoic transition 35,36 ; (ii) the extension of the original evaluation by means of b/a ratio as a function of reversal frequency proposed by McFadden et al. (ref. 27 ) and Coe and Glatzmaier (ref. 24 ) for Pre-Jurassic times, to achieve a better description of such behavior throughout the Phanerozoic; (iii) comparisons about the mean CMB heat flux and the b/a ratio, in order to verify a possible correspondence between both factors.
In this work, we aim to address these points, in order to provide new information for the discussions that linked the long-term variations of the geomagnetic reversals, the geodynamo's stability and the geodynamic processes throughout the Phanerozoic.

Methods
IHMP: selection criteria for the paleomagnetic database. In order to assess of the paleolatitudinal dependence of the paleosecular variation for the IHMP interval (~266.7-228.7 Myr), we conducted a pre-selection of paleomagnetic studies available in literature for this time interval, comprising of 112 works published between 1990 and 2018 based on igneous rocks. Such preliminary database research was carried out by means of academic search engines (e.g., Web of Science (https://www.webofknowledge.com/) and Scopus (https://www. scopus.com/home.uri)) and the IAGA's Global Paleomagnetic Database (http://www.ngu.no/geodynamics/ gpmdb/). Regarding the scarcity of studies based on highly sensitive magnetometers, which were often associated to low accuracy rock magnetism investigations, we did not consider datasets prior to 1990, according to similar procedures adopted by De Oliveira et al. (ref. 26 ).
From the preliminary dataset, we built the "final" paleomagnetic database by means of the following selection criteria: (1) all works that did not provide directional, characteristic remanent magnetization (ChRM) data per site and site coordinates, as well as at least ten sampling sites (N < 10) were ruled out; (2) preference was given to the selection of works which provide high-quality paleomagnetic poles in accordance to the Van der Voo (ref. 37 ) quality criteria; (3) the selected studies shall be related to level ≥ 4 of the GPMDB Demagcode procedure protocol 38,39 as reliable analyses of VGP dispersion datasets can be prevented due to the employment of inadequate Scientific REPORTs | (2019) 9:282 | DOI:10.1038/s41598-018-36494-x demagnetization procedures 40 ; (4) only studies that succeeded in the recalculation of its paleomagnetic pole (s) and associated paleolatitude (s) by means of its ChRM directional data and site coordinates were considered. In order to remove spurious data that could be possibly related to eventual excursional fields or lightning occurrences that may influence the VGP angular dispersion, due to the size of the paleomagnetic datasets De Oliveira et al. (ref. 26 ), all selected paleomagnetic datasets were submitted to the Vandamme (ref. 41 ) iterative method. We  ; S: VGP dispersion; S B , S u and S l : respectively, between-site VGP dispersion and its associated upper and lower 95 per cent confidence limits (obtained by the bootstrap method); ΔS: difference between S and S B ; λ: paleolatitude; A Aikhal; B Udzha; C Abagalakh; D Delkan; E Group C; F Sytikanskaya; G Group B. Datasets 2, 8, 10 and 15 (in bold/italic) provided A95 values that fall out of the A95min/A95max, range, and hence were not considered for the paleomagnetic data processing and the subsequent Model G curve fitting for the IHMP (see section "IHMP: selection criteria for the paleomagnetic database"). The resulting paleomagnetic database from application of selection criteria #1 -#4 is constituted of 16 VGP datasets, provided by 12 paleomagnetic studies (which corresponds to ~14.3% of the pre-selected works), from igneous-based lithologies (Table 1; Supplementary Information Tables S1 and S2). However, as some of the datasets exhibit considerably high k-values (>200), we adopted an additional procedure to evaluate whether such corresponding VGP distributions represent adequate PSV samplings, by means of application of the Deenen et al. (ref. 43 ) criteria. It provides a N-dependent A95 envelope defined by a range of upper (A95 max ) and lower (A95 min ) limits, in which the observed A95 shall be within for a sufficient PSV sampling. As discussed by some authors (e.g., ref. 43,44 ), datasets that provide A95 > A95 max may contain additional scatter contributors, whereas A95 < A95 min could be considered as an indicator for an EMF spot-reading record. It was noticed that four of the select datasets (datasets # 2, 8, 10 and 15) provided A95 values that fall out of the A95 min /A95 max range, and hence were not considered for the paleomagnetic data processing and the subsequent Model G curve fitting for the IHMP.
IHMP paleomagnetic data processing. From the paleomagnetic database, all VGP angular dispersion data were calculated by means of Eq. (1). Upper and lower limits for S (S u and S l , respectively) were estimated as suggested by the bootstrap method. Obtaining angular dispersion data due to the PSV (S B ) can be done by minimization of sampling and measurement errors 25 by means of the following relationship: where n and S w are, respectively, the average number of samples per site and the within-site dispersion. The relation S n / w 2 is the correction factor for the within-site dispersion of a given VGP dataset, which is given by 42 : where α 95 is the mean value of α 95 for the VGP dataset. S B data are also displayed in Table 1. The mean difference between S and S B is quite small (~1.9°), which could be an indirect indicator for the adequacy of the selection criteria adopted in this work. For the evaluation of VGP dispersion data regarding the paleolatitude for the IHMP, we considered the S B (λ) data. Model G curve fitting. For evaluation of the paleolatitudinal dependence of the VGP dispersion data to the selected S B (λ) dataset for the IHMP, we performed a curve fitting based on the Model G (ref. 31 ) by means on the Levenberg-Marquardt method, which is an iterative regression method for solving nonlinear least square problems, by means of a stabilization parameter that assures the convergence of the goal function for a minimum value by choosing Steepest Descent or Gauss-Newton methods (ref. 45 ). It was done by means of the modulus "scipy.optimize.leastsq", available at the Python online repository ScyPy (https://scipy.github.io/devdocs/generated/scipy.optimize.least_squares.html). From the best Model G fitted curve, we carried out the shape parameters a and b for the S B (λ) dataset to the IHMP, which will be discussed later.

Results
Evaluation of the paleolatitudinal dependence of the VGP dispersion data for the IHMP. By the hemispheric representation of the selected database along with its corresponding paleolatitudes (Fig. 1), it was not possible to observe any evidence for an equatorial asymmetry between the S B dispersion datasets related to both Southern and Northern hemispheres (open and full circles, respectively) -which could be reasonably explained by the assumption of the GAD hypothesis, as previously discussed by Biggin et al. (ref. 25 ). The S B (λ) distribution, in association to its best fitted Model G curve (which resulted in shape parameters . . ), clearly exhibits a low paleolatitudinal dependence trending pattern (ranging from S B ~13.8° to ~17.0° at (paleo)latitudes = 0° and 90°, respectively). All the three S B (λ) curves exhibit similar shapes, which is compatible to a low (paleo)latitudinal dependence due to smaller antisymmetric contribution during high reversal rate intervals. Nevertheless, the IHMP interval (average reversal rate of ~5.9 Myr −1 ) exhibits higher S B values at low paleolatitudes in comparison to those reported for lower reversal frequency intervals, as the CNS 25 (~8.7°) and the PCRS 26 (~9.4°) -and similar to the observed to the 0-5 Myr and Jurassic intervals, of similar reversal frequency (4-5 Myr −1 and 4.6 Myr −1 , respectively).
Furthermore, in order to compare the observed paleolatitudinal trending pattern and shape of the VGP dispersion curve for the IHMP with other high mean reversal frequency intervals we also demonstrated in Fig. 1, the Model G best-fit curves respectively provided for Jurassic times from Group 1 of Biggin et al. (ref. 25 ) and for the last 5 Myr 28 .
All curves exhibit the same low paleolatitudinal trending patterns, which has been discussed in literature (e.g., ref. 25 and ref. 27 ) as being due to a major symmetric family contribution in comparison to the influence from the antisymmetric family. Such effect leads to higher (lower) values for the shape parameter a (b) in comparison to low reversal frequency intervals, as the CNS (Johnson & McFadden,ref. 4 ). The IHMP (red) and Jurassic (blue) curves evolved similarly within the 0-90° paleolatitudinal interval, although the IHMP S B (λ) curve exhibit lower S B at lower and higher paleolatitudes. The VGP dispersion curves for both Jurassic and 0-5 Myr intervals provided shape parameters that are compatible to those found for IHMP (Jurassic: = . . 0 014 . Additionally, the mean reversal rate for the Jurassic 25 (~4.6 Myr −1 ) and the 0-5 Myr 28 (~4-5 Myr −1 ) intervals are quite similar. We estimated the mean average reversal frequency for the IHMP (for more detail, see description in "Evaluation of the time evolution of the b/a ratio" section) as ~5.9 Myr −1 , for the ~266.7-228.7 Myr suggested for this period, which is higher than the previous two intervals. By comparison, the higher (lower) values of mean average reversal frequency (b/a ratio) found for IHMP in comparison to the last 5 Myr and Jurassic could indicate the inverse relationship between mean reversal rate and b/a ratio, as expected, and the even lower influence of the antisymmetric family for the IHMP.
Evaluation of the time evolution of the b/a ratio. As discussed by several authors 13,14,18,20 Table 2. Secondary (symmetric) (a l u ) and primary (antisymmetric) (b l u ) harmonic terms, and estimates for its relative contribution (b/a ratio), carried out from Model G fitting curves from this work and other studies based on VGP dispersion analyses for most of the Phanerozoic. u (l): upper (lower) limits for the shape parameters is itself comparable to the variations of the heat flux patterns over the CMB, as suggested by numerical modeling works of mantle convection 46,47 .
In order to contribute to this debate, we also conducted an evaluation aiming to track the time evolution of the relative contribution of dipole/non-dipole fields derived from paleomagnetic data -by means of b/a ratios -and its possible correspondence with time variations in relative CMB heat flux throughout most of the Phanerozoic. The results for b/a ratios were provided both by this work and other studies, which together comprise contiguous, million-year scale intervals that exhibited high and low mean reversal rates throughout the Phanerozoic: (1) PCRS 26 ; (2)  It is noticeable that the time evolution of the b/a ratio matches, in an inverse relationship, the smooth trending pattern for the relative CMB heat flux from the PCRS to the present times, as provided by Olson & Amit (ref. 9 ). It is important to highlight that the b/a ratios were carried out with PSV analyses from Model G fittings of VGP dispersion curves, which are not of straightforward interpretation in terms of physical processes, because their origins rely on a number of different factors 19 .
Nevertheless, our results point out that the relative contribution of equatorially antisymmetric to symmetric spherical harmonics terms, given by the Model G, could be inversely related to the CMB heat flux variations, indicating that higher axial (non-axial) dipole contributions may be expected for lower (higher) relative CMB heat flux intervals for the last 270 Myr. As discussed previously, high/low b/a ratios would be considered, for a given time interval, as a predictor of low/high reversal frequency states 24 -which in turn could reflects high/low CMB heat flow conditions, as discussed by some authors 47,49,50 .
Such observations would shed some light on the physical meaning of the Model G shape parameters a and b, what can partially explain the adequacy of this phenomenological model for most of the Phanerozoic. Surely new investigations aiming to extend back in time the b/a ratio coverage herein presented, and with more time resolution, are demanded to verify the hypothesis.