<figure><img src="image_1.png" /></figure>

HAL
open science

# A model for the geomagnetic field reversal rate and constraints on the heat flux variations at the core-mantle boundary

Vincenzo Carbone, Tommaso Alberti, Fabio Lepreti, Antonio Vecchio

► To cite this version:

Vincenzo Carbone, Tommaso Alberti, Fabio Lepreti, Antonio Vecchio. A model for the geomagnetic field reversal rate and constraints on the heat flux variations at the core-mantle boundary. Scientific Reports, 2020, 10, 10.1038/s41598-020-69916-w. insu-03717007

HAL Id: insu-03717007

https://insu.hal.science/insu-03717007

Submitted on 8 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

OPEN

# A model for the geomagnetic field reversal rate and constraints on the heat flux variations at the core-mantle boundary

Vincenzo Carbone<fnref n="1" />,²✉, Tommaso Alberti D³, Fabio Lepreti D<fnref n="1" />,² & Antonio Vecchio⁴,⁵

A striking feature of many natural magnetic fields generated by dynamo action is the occurrence of polarity reversals. Paleomagnetic measurements revealed that the Earth's magnetic field has been characterised by few hundred stochastic polarity switches during its history. The rate of reversals changes in time, maybe obeying some underlying regular pattern. While chaotic dynamical systems can describe the short-term behaviour of the switches of the Earth's magnetic polarity, modelling the long-term variations of the reversal rate is somewhat problematic, as they occur on timescales of tens to hundreds of millions of years, of the order of mantle convection timescales. By investigating data of geomagnetic reversal rates, we find the presence of cycles with variable frequency and show that the transition towards periods where reversals do not occur for tens of million years (superchrons) can be described by a second-order phase transition that we interpret to be driven by variations of the heat flux at the core-mantle boundary (CMB). The model allows us to extract from the reversal sequence quantitative information on the susceptibility of the reversal rate caused by changes in the CMB heat flux amplitude, thus providing direct information on the deep inner layers of the Earth.

The Earth's magnetic field is generated, through the dynamo effect, by the motion of liquid metal in the core<fnref n="1" />. Paleomagnetic measurements revealed that the dipole component of the magnetic field has undergone several polarity switches during the Earth's history, the periods spent in a single polarity being stochastically distributed²⁻⁴. The occurrence of polarity reversals is a common property of the dynamo phenomenon being observed also in the Sun⁵, even though there the magnetic field reverses quasi periodically every ~11 years⁶,⁷, and in laboratory excited dynamo⁸. The physical mechanism which drives the Earth's magnetic field to reverse is still controversial. Some simple dynamical systems can describe the reversal dynamics⁹,¹⁰ up to timescales of several Myr (longer than the core magnetic diffusion time of about 10⁵ yr), and reversals have been observed in 3D global direct numerical simulations¹¹⁻¹⁴. Although statistical inferences of Earth's magnetic field reversals are conditioned by their relatively small number, it has been argued that they can be described as a sequence of stochastic independent events¹⁵, thus following a Poisson statistics¹⁶,¹⁷. As a different point of view, it has been shown that the sequence of periods spent in a single magnetic polarity is characterised by the presence of long-range correlations¹⁸,¹⁹ and that it can be better described by a Lévy process²⁰. The statistical difference could be ascribed to the fact that, on long time scales, the rate of reversals is not constant¹⁷, showing transitions from periods of rapid polarity switches to very long intervals of fixed polarity, called superchrons and lasting for a few tens of Myr, such as occurred between the Middle Jurassic and Middle Cretaceous periods. Data analysis of polarity reversals revealed that different cycles seem to be present in the variations of the rate, whose periods are related to typical time scales of mantle convection processes²¹⁻²⁷. On this basis, it has been suggested that the changes of the reversal rate may have been triggered by fluctuations in CMB heat flow, either global or localised in equatorial regions²³,²⁸⁻³⁰. In several works based on geodynamo numerical simulations, the effects of CMB

<footnote n="1">¹Dipartimento di Fisica, Università Della Calabria, Ponte P. Bucci, Cubo 31C, 87036 Rende, CS, Italy. ²Istituto Nazionale di Astrofisica (INAF), Direzione Scientifica, Rome, Italy. ³INAF-IAPS Istituto di Astrofisica e Planetologia Spaziali, Via Fosso del Cavaliere 100, 00133 Rome, Italy. ⁴Radboud Radio Lab, Department of Astrophysics/IMAPP, Radboud University, P.O. Box 9010, 6500GL Nijmegen, The Netherlands. ⁵LESIA - Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 place Jules Janssen, 92195 Meudon, France. ⁶email: vincenzo.carbone@fis.unical.it</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Upper panel: geomagnetic polarity back to 375.3 Ma. The polarity is denoted by 1 for the present polarity and 0 for reversed. A line connects the values, therefore black regions correspond to a high reversal rate and white regions to intervals of constant value. Lower panel: reversal rate time series obtained by averaging over sliding windows of length $\Delta t=8$ Myr (see Supplementary Information).</figcaption></figure>

heat flux amplitude and morphology have been investigated and discussed in the context of observed reversal rate changes³¹⁻³⁶.

## Results

We analyse data²⁶ obtained from the latest geomagnetic polarity time scale (2012), integrated with data from the late Devonian³⁷. The time sequence of polarity reversals was reconstructed for the whole Phanerozoic eon²⁶. The time sequence of reversals and the time evolution of the rate of reversals $\gamma(t)$, obtained by averaging reversals over sliding windows of width 8 Myr, are reported in Fig. 1. We analysed the geomagnetic series $\gamma(t)$ in terms of empirical modes, through the Empirical Mode Decomposition (EMD, see “Methods”), a technique which is particularly suitable to process non-stationary time series³⁸ like the reversal rate series. Through this technique, a time series is decomposed into empirical modes, called Intrinsic Mode Functions (IMFs), characterized by different frequencies and therefore it is possible to extract relevant timescales involved in the non-stationary process under investigation. For the reversal rate time series $\gamma(t)$ we obtained a sequence of 11 IMFs $C_j(t)$ ($j=0,1, ..., 10$), which are shown in Fig. 2 along with the residue (see “Methods”). For each mode we can extract an instantaneous frequency $\omega_j(t)$, whose time variations describe the non-stationary processes underlying the observed variability of the reversal rate, and the time average of $\omega_j(t)$ allows us to define a typical period $T_j = 2\pi/\langle\omega_j\rangle$ (where the time average is denoted by angular brackets) associated to each j-th mode.

The probability density functions of the instantaneous periods $2\pi/\omega_j(t)$ (see Fig. S1 of Supplementary information) indicate a compatibility with previously reported cycles. Periods shorter than about 40 Myr²¹,²²,²⁴,²⁶,³⁹ are found in the modes with $j \le 5$. Such variability has been related to CMB heat flux changes and plume dynamics within the Earth's mantle. Modes characterized by longer periodicities, associated by our analysis to EMD modes with lower frequencies, are perhaps hidden in the low-frequency Fourier peaks commonly used to recover reversal rate periodicities. Previously reported periodicities also include time scales longer than 100 Myr, which can arise from subducting lithospheric slabs reaching the CMB²⁷,³⁰. These time scales are compatible with the instantaneous periods found in the modes $j=9,10$. Moreover, note that the residue of the EMD (see Fig. 2) is a decreasing function of age, thus indicating that the global rate of reversals is not constant, and the global average persistence of the geomagnetic field in a single polarity state decreased, on average, going from 400 Ma to the present.

EMD results suggest that the geomagnetic system, for timescales longer than the magnetic diffusion time, can be modelled through transitions between chrons induced by a continuous underlying stochastic process, different chrons being characterized by the average frequency of switches $T_j^{-1}$ between different polarity states of the magnetic field. To identify the number of states present in the system, we describe the reversal rate variations in terms of a stochastic dynamical system and we assume that a transition among states of different reversal rates is triggered by a stochastic forcing, namely the continuous change of heat flux at CMB. We use the Langevin equation $dx = -U'(x)dt + \sigma dw$ to describe the dynamics of the rate changes, where $x$ is a state variable, which in our case represents an IMF or a sum of IMFs, $U(x)$ is a given potential, $U'(x) = dU/dx$, $\sigma$ is the noise level and $dw$ is a Wiener process, i.e., a stochastic process with independent Gaussian increments, which describes the
<figure><img src="image_5.png" /><figcaption>**Figure 2.** IMFs obtained through the EMD technique applied to the reversal rate time series. The residue of the decomposition is shown in the lowest right panel, while the IMFs $C_j(t)$, for $j=0,1, ..., 10$, are reported in the other panels.</figcaption></figure>

stochastic CMB heat flux. The potential function $U(x)$ can be evaluated by means of the Fokker-Planck (FP) equation associated with the Langevin model, describing the time evolution of the probability density function (pdf) $p(x,t)$ (see "Methods"). Moreover, if the potential function $U(x)$ can be written as a polynomial function of even order $k$ and positive leading coefficient, i.e., $U(x) = u_0 + u_1 x + u_2 x^2 + ... + u_k x^k$ and $u_k > 0$, its order is related to the number of available states for the reversal rates $n$, i.e., $n=k/2$$^{40,41}$. In this way, from the pdfs of the IMFs $C_j$ we can calculate the potentials for all the EMD modes (see Fig. S3 of Supplementary information) by assuming that the stochastic term, i.e., the Wiener process, is representative of processes occurring at timescales which are shorter than the mean timescale $T_j$ of each empirical mode $C_j$, and that the amplitude of the noise corresponds to the standard deviation of each $C_j$$^{40,41}$. Two types of potential shapes are present in the dataset, thus reflecting the number of possible states for the reversal rate at different time scales $T_j$. Namely, single-well potentials related to high-frequency CMB changes, for the set of modes $H = \{0 \le j \le 4\}$, and double-well potentials, related to low-frequency CMB changes, for the set of modes $L = \{5 \le j \le 10\}$. Some of the EMD modes corresponding to single-well potentials, namely the modes $j=2,3,4$, present periodicities (between 14 and 30 Myr, see Table T1 of Supplementary Information) that are close to cycles already identified using different techniques$^{21,22,24,26,39}$, as already mentioned above. The period of the EMD mode $j=10$ with a double-well potential is close to what has already been observed as a superchron period$^{25}$. In addition, the EMD analysis suggests the presence of characteristic intermediate time scales, in the range 50–100 Myr (see Table T1 of Supplementary Information). These periods, corresponding to asymmetric double-well potentials, are perhaps hidden in the large width of low-frequency Fourier modes reported in previous analyses$^{25,26}$.

The approach based on stochastic Langevin models has been proven successful in reproducing the dipole field variability observed both in paleomagnetic data covering the last 2 Myr$^{42,43}$ and in numerical geodynamo simulations$^{43}$. It has been shown that such models provide a good description of the axial dipole field dynamics and a reliable prediction, through the stationary solution of the FP equation, of its probability distribution. Here, we follow a similar approach and we assess the significance of our Langevin model by comparing, firstly, the partial reconstruction of the geomagnetic reversal rate signal obtained by summing the IMFs of the set of modes $L = \{5 \le j \le 10\}$ to a realisation obtained from the stochastic Langevin model (see Fig. S5 of Supplementary Information), and, then, the stationary solution of the FP equation to the pdfs of the partial reconstruction and of the Langevin model (see Fig. S6 of Supplementary Information). A quite good agreement is found between the pdfs, thus confirming the validity of our approach.

The dynamics of the obtained EMD modes allows us to interpret the transition from high-frequency chronos towards low-frequency superchrons as a kind of phase transition$^{44}$, that we assume to be driven by stochastic
<figure><img src="image_6.png" /><figcaption>Figure 3. The left panel shows the potential manifold U(r,Γ), the right panel shows the manifold projected on the (r,Γ) plane. The manifold has been drawn by connecting the potentials obtained for each EMD mode. The parameters u and h are obtained through best fits on the potentials obtained from the empirical modes.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4. Correlation time τ (black circles), estimated from the two-times correlation coefficient of the observed reversal rates for empirical modes from j=6 to j=10, as a function of the scaling parameter r. The superposed red curve corresponds to a fit obtained through an inverse power law function, reported on the figure along with the best fit parameters.</figcaption></figure>

fluctuations of the heat flux at the CMB. High-frequency chronos correspond to disordered states characterized by periods of rapid polarity reversals and stronger CMB activity. On the contrary, low-frequency rates correspond to more organized states, characterized by stable long residence times in a single magnetic polarity, with smaller CMB heat flux variations and weaker mantle plume activity. In the framework of mean-field approximation$^{44}$, let us consider a continuous order parameter from the set of standardised EMD modes $\Gamma = \{C_j^\sigma(t)\}$ (see Supplementary Information) and let us write up the potentials in terms of the manifold

$$U(\Gamma) = r\Gamma^2 + u\Gamma^4 - h\Gamma \quad (1)$$

shown in Fig. 3. The transition from the single-well potential to the asymmetric double-well one, happens when the parameter r changes sign. Then we define r in terms of CMB heat flow as $r=(Q-Q_c)/Q_c$, where Q_c represents a critical threshold. In other words, according to the mean field approach of phase-transitions$^{44}$, we assume that the transition to superchrons happens when the CMB heat flow Q becomes smaller than the critical value Q_c. The correlation time τ, estimated from the two-times correlation coefficient G(t₁, t₂) of the observed reversal rates for each mode

$$G(t_1, t_2) = \Gamma(t_1)\Gamma(t_2) \approx \exp\left[\frac{-|t_1 - t_2|}{\tau}\right], \quad (2)$$

(where the angular brackets denote time averaging) shows a power law dependence $\tau \approx 1/r^\alpha$ and thus diverges when Q = Q_c (see Fig. 4), where $\alpha = 0.48 \pm 0.03$, in close agreement with the scaling exponent $\alpha = 1/2$ required by the mean-field approximation of second-order phase transitions$^{44}$. This confirms that a kind of second-order phase transition is at work within the complex geodynamo system.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Estimated excess of CMB heat flow (black circles) as a function of the expected relative increase of reversals. The superposed red curve corresponds to an exponential function, obtained through a fit on the data, reported on the figure along with the best fit parameters.</figcaption></figure>

According to the mean field approximation of second-order phase transitions, the susceptibility $\chi$ is given by $\chi^{-1} = dh/d\Gamma$. The equilibrium solution in the mean field approach is obtained by minimising the potential $U$ (Eq. (1)) with respect to $\Gamma$. This gives $h = 2r\Gamma + 4u\Gamma^3$, from which $\chi^{-1} = 2r + 12u\Gamma^2$. In the mean field theory, for the single-well minimum, which corresponds to the disordered phase when the geomagnetic field reverses at a high rate, we have $\Gamma^2 = 0$ when $r > 0$, therefore $\chi^{-1} = 2r$. This result is very interesting in our case, as it allows us to infer the number of reversals induced, starting from a reference value $n_0$, by amplitude variations of the heat flux at CMB. In fact, since the heat flux variations due to a variation of the order parameter must be roughly proportional to the variations of the $h$ field of the model, namely $\Delta Q \approx \Delta h$, we get $\Delta Q \approx \chi^{-1} \Delta \Gamma \approx 2r \Delta \Gamma$. By estimating $\Delta \Gamma$ as $\Delta \Gamma \approx 1/(T_j - T_{j+1})$, using the characteristic periods $T_j$ of EMD modes $j = 6-10$, we can use the predictive property of the susceptibility to directly infer about the heat flow fluctuations at CMB, relative to a reference value $Q_0$, required to increase the reversal rate from the reference value $n_0$, as shown in Fig. 5. An empirical relation can be obtained through a fit on the data with the following exponential function

$$ \frac{\Delta Q}{Q_0} = A \exp \left[ \frac{\Delta n}{n_0 \nu} \right] \qquad (3) $$

where $\Delta n$ is the variation in the reversal rate and the best fit parameters result $\nu = 2.08 \pm 0.03$ and $A = 0.07 \pm 0.01$. This means that fluctuations of the order of about $\Delta Q \approx 0.18 Q_0$ of CMB heat flow should be enough to double the reversal rate in the geodynamo system.

## Discussion

It is worth remarking that our model represents a novel conceptual way of approaching geodynamo reversals over the last 400 Myr, by using the existing reversal dataset. This is very interesting because it allows us to provide constraints, directly from the observed geomagnetic reversal sequence, on the expected CMB heat flux variations. Using Eq. (3) with $n_0 = 1.43 \text{ Myr}^{-1}$, which is the median value of the reversal frequency over the interval covered by the considered dataset, and assuming $\Delta n = \gamma(t)$, in Fig. 6 we report the expected evolution of the heat flux variations required to account for the observed reversal rate. The spike between 160 and 170 Ma corresponds to the peak of the reversal rate in the same interval. It is worthwhile to point out that the number of reversals occurring in this interval is likely to be overestimated (by a factor up to about 2) because some anomalies caused by geomagnetic field fluctuations could be incorrectly interpreted as reversals$^{45,46}$. This could significantly reduce the amplitude of the spike in the reversal rate and, consequently, also the corresponding spike in $\Delta Q/Q_0$, however without making them to be suppressed.

By using Eq. (3) to make a comparison with the heat flux variation at CMB of numerical simulations, we can estimate $\Delta n/n_0 \approx 7-8$ when $\Delta Q/Q_0 \approx 2-3$. This is in agreement with direct estimates of the CMB heat flux evolution obtained using numerical geodynamo models with core evolution$^{30,47}$, which show that an increase of about 10 reversals, starting from few reversals, is achieved if the CMB heat flux varies by a factor of 2 or 3. Also simpler purely chemically-driven dynamos without core evolution seem to be compatible with the results presented here. For example, the simulations of Driscoll and Olson$^{35}$ show that variations in the thermo-chemical CMB buoyancy flux of a factor 2 can lead to a change from minor reversal activity (with reversal frequency less than ~ 1 Myr$^{-1}$) to frequently reversing activity with frequency ~ 4.5 Myr$^{-1}$. Olson et al.$^{34}$ considered a tomographic forcing of the CMB heat flux, finding that variations in its mean amplitude of a factor 2 lead to a relative increase of reversal frequency of a factor 6-10. It has also been found that the rate of reversals scales with the buoyancy-flux-based Rayleigh number $Ra^{48}$. Since the Rayleigh number is directly related to the heat flux at CMB, our relation (3) could be used to estimate values of $Ra$ on varying the fluctuations of the number of reversals and provide constraints on this dimensionless parameter. It is also worth to mention here the following still open research problem in this context. Mantle convection models predict temporal variations in total CMB heat flux of the order of few tens of percent over the last 400 Myr$^{30,49}$. According to the dynamo simulation results
<figure><img src="image_9.png" /><figcaption>Figure 6. Estimated relative CMB heat flux variations over the last 375 Myr obtained from the observed reversal rate variations by using Eq. (3). The dashed red line corresponds to the average CMB flux.</figcaption></figure>

discussed above and the stochastic model presented here, such variations are unlikely to produce significant changes in reversal frequency.

As a conclusion, we introduced a conceptual model for the reversal rate of the geomagnetic field, based on a second-order phase transition driven by stochastic fluctuations, which allows us to recover the CMB heat flux time variations which drove the reversals over the past 400 Myr. The CMB heat flux fluctuations are derived from the available reversal dataset and can be used to infer a constraint on the variations required to increase the number of geomagnetic reversals by a given amount.

## Methods

**Empirical mode decomposition.** The time evolution of the reversal rate $\gamma(t)$ has been analysed through the Empirical Mode Decomposition (EMD)$^{38}$. This decomposition is orthogonal and complete, and partial summations of EMD modes capture typical processes happening over different time scale ranges. In addition, EMD is able to avoid some limitations present in other decomposition analysis techniques. Differently from Fourier or wavelet analysis, EMD does not require any priori assumption on the functional form of the basis of the decomposition. This allows us to extract local nonstationarity and nonlinearity features from each time series.

Through this technique the reversal rate time series has been decomposed as

$$ \gamma(t) = \sum_{j=0}^{m} C_j(t) + r_m(t) \quad (4) $$

where each Intrinsic Mode Function (IMF) $C_j(t)$ represents an oscillating function with both amplitude and phase modulation in time $C_j(t) = A_j(t) \cos[\phi_j(t)]$, while $r_m(t)$ is the residual of the decomposition, obtained after all the oscillating functions have been extracted. The various IMFs are extracted recursively using a sifting procedure$^{38}$. As a first step the sets of local maxima and minima of $\gamma(t)$ are identified and interpolated by cubic splines, thus obtaining an upper and a lower envelope, respectively, as well as their local mean $m_1(t)$. The difference $H_1(t) = \gamma(t) - m_1(t)$ represents the first IMF if it satisfies two conditions: (1) the number of extrema and zero-crossings does not differ by more than 1; (2) the mean value of the envelopes obtained from the local maxima and minima is zero at any point. If these conditions are not fulfilled, the procedure is repeated on the $H_1(t)$ time series, and a new difference $H_{11}(t) = H_1(t) - m_{11}(t)$ is obtained, where $m_{11}(t)$ is the mean of the new envelopes. The process is iterated until the s-th iteration $H_{1s}$ fulfils the IMF properties. To avoid loss of information about amplitude and frequency modulations, a stopping criterion between two consecutive iterations is introduced through the parameter

$$ \eta = \sum_{t=0}^{T_{\max}} \left[ \frac{[H_{1(s-1)}(t) - H_{1s}(t)]^2}{H_{1(s-1)}^2(t)} \right] \quad (5) $$

where $T_{\max}$ denotes the last time instant of the dataset. When $\eta < \eta_{cr}$, with the threshold fixed at $\eta_{cr} = 0.3$ in our case, the process is stopped. When the first IMF is calculated, the first residue $r_1(t) = \gamma(t) - C_1(t)$ is processed again, thus obtaining $C_2(t)$ and $r_2(t)$. The whole procedure is carried on until $C_m(t)$ or $r_m(t)$ are almost constant, or $r_m(t)$ is monotonic. From our time series we extracted 11 IMFs ($m=10$) along with the residual $r_m(t)$. It is worthwhile to remark that for each IMF we can extract the instantaneous frequencies $\omega_j(t) = d\phi_j/dt$, from which we can obtain the instantaneous periods $2\pi/\omega_j(t)$. The time averages of the frequencies $\omega_j(t)$ allows us to define the typical periods $T_j = 2\pi/\langle\omega_j\rangle$ associated to each mode.
**Langevin stochastic model.** The statistical properties of IMFs can be described by using a simple stochastic model based on a Langevin process

$$dx = -U'(x)dt + \sigma dw \quad (6)$$

where $x$ is the state variable (in our case an IMF or a sum of IMFs), $U$ is the potential function, acting as a drift term, $\sigma$ is the noise level, and $dw$ is a stochastic process. Specifically, the stochastic process mimicking the dynamics on timescales shorter than the characteristic period of each IMF, is a Wiener process, defined as a stochastic process with independent Gaussian increments, i.e., $dw$ is normally distributed with zero mean and variance $t$ (also called Brownian motion due to its historical connection with the related physical process). A very similar approach has been employed to describe the temporal behaviour of the dipole field both in paleomagnetic data and numerical geodynamo models$^{42,43}$. The potential function $U(x)$ can be derived from the stationary solution of the associated Fokker-Planck (FP) equation to the Langevin model which reads

$$\frac{\partial p(x, t)}{\partial t} = \frac{-\partial}{\partial x} [U'(x)p(x, t)] + \frac{\sigma^2}{2} \frac{\partial^2}{\partial x^2} [p(x, t)] \quad (7)$$

where $p(x,t)$ is the probability density function (pdf). By searching for a stationary solution of the FP equation, we obtain $-\partial/\partial x \{[U'(x) p(x,t)] + \sigma^2/2 \partial/\partial x [p(x,t)]\} = 0$, whose solution is

$$p(x) = p_0 \exp \left[ \frac{-2U(x)}{\sigma^2} \right], \quad (8)$$

where $p_0$ is a normalization factor. Thus, by deriving the pdf of each IMF we can get the corresponding potential function as

$$U_j(x) = \frac{-\sigma^2}{2} \log \left[ \frac{p_j(x)}{p_0} \right] \quad (9)$$

being $p_j(x)$ the empirical pdf of each IMF and choosing $p_0$ such that $p_j(x)$ has a unit integral. Moreover, if the potential function $U(x)$ can be written as a polynomial function of even order $k$ and positive leading coefficient, i.e., $U(x) = u_0 + u_1 x + u_2 x^2 + ... + u_k x^k$ and $u_k > 0$, its order is related to the number of available states for the reversal rates $n$, i.e., $n = k/2^{40,41}$.

## Data availability

The data analysed in this work are available in the Supplementary materials of Melott, A.L., Pivarunas, A., Meert, J.G. & Lieberman, B.S. Does the planetary dynamo go cycling on? Re-examining the evidence for cycles in magnetic reversal rate. *International Journal of Astrobiology*, https://doi.org/10.1017/S1473550417000040 (2017).

Received: 3 September 2019; Accepted: 20 July 2020

Published online: 03 August 2020

## References

1. Merril, R. T. & McElhinny, M. W. *The Earth's Magnetic Field: Its history, origin and planetary perspective* (Academic Press, Cambridge, 1983).

2. Jacobs, J. A. *Reversals of the Earth's Magnetic Field* (Cambridge University Press, Cambridge, 1994).

3. Cox, A. Geomagnetic reversals. *Science* **163**, 237–245 (1969).

4. Valet, J.-P. & Fournier, A. Deciphering records of geomagnetic reversals. *Rev. Geophys.* **54**, 410–446 (2016).

5. Charbonneau, P. Dynamo models of the solar cycle. *Living Rev. Solar Phys.* **7**, 3–91 (2010).

6. Hathaway, D. H. The solar cycle. *Living Rev. Solar Phys.* **7**, 1–65 (2010).

7. Vecchio, A., Laurenza, M., Meduri, D., Carbone, V. & Storini, M. The dynamics of the solar magnetic field: polarity reversals, butterfly diagram and quasi-biennial oscillations. *Astrophys. J. Lett.* **749**, 27 (2012).

8. Monchaux, R. et al. Generation of a magnetic field by dynamo action in a turbulent flow of liquid sodium. *Phys. Rev. Lett.* **98**, 044502 (2007).

9. Bullard, E. C. The stability of a homeopolar dynamo. *Proc. Camb. Philos. Soc.* **51**, 744–760 (1955).

10. Pétrélis, F., Fauve, S., Dormy, E. & Valet, J.-P. Simple mechanism for reversals of earth's magnetic field. *Phys. Rev. Lett.* **102**, 144503 (2009).

11. Olson, P. L., Glatzmaier, G. A. & Coe, R. S. Complex polarity reversals in a geodynamo model. *Earth Planet. Sci. Lett.* **304**, 168–179 (2011).

12. Glatzmaier, G. A. & Roberts, P. H. A three-dimensional self-consistent computer simulation of a geomagnetic field reversal. *Nature* **377**, 203–209 (1995).

13. Amit, H., Leonhardt, R. & Wicht, J. Polarity reversals from paleomagnetic observations and numerical dynamo simulations. *Space Sci Rev* **155**, 293–335 (2010).

14. Glatzmaier, G. A. & Coe, R. S. Magnetic polarity reversals in the core. In *Treatise on Geophysics* 2nd edn, Vol. 8 (ed. Olson, P.) 283–298 (Elsevier, New York, 2015).

15. McFadden, P. L. Statistical tools for the analysis of geomagnetic reversals. *J. Geophys. Res.* **89**, 3363–3372 (1984).

16. McFadden, P. L. & Merril, R. T. Asymmetry in the reversal rate before and after the Cretaceous Normal Polarity Superchron. *Earth Planet. Sci. Lett.* **149**, 43–47 (1997).

17. Constable, C. On rates of occurrence of geomagnetic reversals. *Phys. Earth Planet. Inter.* **118**, 181–193 (2000).

18. Jonkers, A. R. T. Long-range dependence in the Cenozoic reversal record. *Phys. Earth Planet. Inter.* **135**, 253–266 (2003).

19. Ryan, D. A. & Sarson, G. R. Are geomagnetic field reversals controlled by turbulence within the Earth's core?. *Geophys. Res. Lett.* **34**, L02307 (2007).

20. Carbone, V. et al. Clustering of polarity reversals of the geomagnetic field. *Phys. Rev. Lett.* **96**, 128501 (2006).

21. Mazaud, A., Laj, C., de Seze, L. & Verosub, K. L. 15-Myr periodicity in the frequency of geomagnetic reversals since 100 Ma. *Nature* **304**, 328–330 (1983).
22. Raup, D. M. Magnetic reversals and mass extinctions. *Nature* **314**, 341–343 (1985).

23. McFadden, P. L. & Merrill, R. T. Lower mantle convection and geomagnetism. *J. Geophys. Res.* **89**, 3354–3362 (1984).

24. Stothers, R. B. Periodicity of the earth's magnetic reversals. *Nature* **322**, 444–446 (1986).

25. Driscoll, P. E. & Evans, D. A. D. Frequency of Proterozoic geomagnetic superchrons. *Earth Planet. Sci. Lett.* **437**, 9–14 (2016).

26. Melott, A. L., Pivarunas, A., Meert, J. G. & Lieberman, B. S. Does the planetary dynamo go cycling on? Re-examining the evidence for cycles in magnetic reversal rate. *Int. J. Astrobiol.* https://doi.org/10.1017/S1473550417000040 (2017).

27. Hounslow, M. W., Domeier, M. & Biggin, A. J. Subduction flux modulates the geomagnetic polarity reversal rate. *Tectonophysics* **742**, 34 (2018).

28. Loper, D. E. & McCarteny, L. Mantle plumes and the periodicity of magnetic field reversals. *Geophys. Res. Lett.* **13**, 1525–1528 (1986).

29. Lay, T., Williams, Q. & Garnero, E. J. The core–mantle boundary layer and deep Earth dynamics. *Nature* **392**, 461–468 (1998).

30. Biggin, A. J. et al. Possible links between long-term geomagnetic variations and whole-mantle convection processes. *Nat. Geosci.* **5**, 526–533 (2012).

31. Lay, T., Hernlund, J. & Buffett, B. A. Core-mantle boundary heat flow. *Nat. Geosci.* **1**, 25–32 (2008).

32. Valet, J.-P., Meynadier, L. & Guyodo, Y. Geomagnetic dipole strength and reversal rate over the past two million years. *Nature* **439**, 802–805 (2005).

33. Glatzmaier, G. A., Coe, R. S., Hongre, L. & Roberts, P. H. The role of the Earth's mantle in controlling the frequency of geomagnetic reversals. *Nature* **401**, 885–890 (1999).

34. Olson, P. L., Coe, R. S., Driscoll, P. E., Glatzmaier, G. A. & Roberts, P. H. Geodynamo reversal frequency and heterogeneous core–mantle boundary heat flow. *Phys. Earth Planet. Inter.* **180**, 66–79 (2010).

35. Driscoll, P. & Olson, P. Superchrons cycles driven by variable core heat flow. *Geophys. Res. Lett.* **38**, L09304 (2011).

36. Olson, P. & Amit, H. Mantle superplumes induce geomagnetic superchrons. *Front. Earth Sci.* **3**, 38 (2015).

37. Hansma, J. et al. Late Devonian carbonate magnet stratigraphy from the Oscar and horse spring ranges. Leonard Shelf, Canning Basin. *Western Aust. Earth Planet. Sci. Lett.* **40**, 232–242 (2015).

38. Huang, N. E. et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. *Proc. R. Soc. Lond. A* **454**, 903–995 (1998).

39. Mazaud, A. & Laj, C. The 15 m.y. geomagnetic reversal periodicity: a quantitative test. *Earth Planet. Sci. Lett.* **107**, 689–696 (1991).

40. Livina, V. N., Kwasniok, F. & Lenton, T. M. Potential analysis reveals changing number of climate states during the last 60 kyr. *Clim. Past* **6**, 77 (2010).

41. Alberti, T. et al. Natural periodicities and Northern Hemisphere-Southern Hemisphere connection of fast temperature changes during the last glacial period: EPICA and NGRIP revisited. *Clim. Past* **10**, 1751–1762 (2014).

42. Buffett, B. A., Ziegler, L. & Constable, C. G. A stochastic model for paleomagnetic field variations. *Geophys. J. Int.* **195**, 86–97 (2013).

43. Meduri, D. & Wicht, J. A simple stochastic model for dipole moment fluctuations in numerical dynamo simulations. *Front. Earth Sci.* **4**, 38 (2016).

44. Hohenberg, P. C. & Krekhov, A. P. An introduction to the Ginzburg-Landau theory of phase transitions and nonequilibrium patterns. *Phys. Rep.* **572**, 1–42 (2015).

45. Tominaga, M., Sager, W. W., Tivey, M. A. & Lee, S.-M. Deep-tow magnetic anomaly study of the Pacific Jurassic Quiet Zone and implications for the geomagnetic polarity reversal timescale and geomagnetic field behavior. *J. Geophys. Res.* **113**, B07110 (2008).

46. Ogg, J. G. Geomagnetic polarity time scale. In *The Geologic Time Scale* (eds Gradstein, F. M. et al.) 85–113 (Elsevier, New York, 2012).

47. Aubert, J., Labrosse, S. & Poitou, C. Modelling the paleo-evolution of the geodynamo. *Geophys. J. Int.* **179**, 1414–1428 (2009).

48. Driscoll, P. & Olson, P. Effects of buoyancy and rotation on the polarity reversal frequency of gravitationally driven dynamos. *Geophys. J. Int.* **178**, 1337–1350 (2009).

49. Zhang, N. & Zhong, S. Heat fluxes at the Earth's surface and core–mantle boundary since Pangeaformation and their implications for the geomagnetic superchrons. *Earth Planet. Sci. Lett.* **306**, 205–206 (2011).

## Acknowledgements

We are very grateful to the anonymous Reviewer for helpful and constructive comments and suggestions on the paper. We acknowledge F. Carbone for useful discussions.

## Author contributions

V.C. and F.L. contributed to the analysis of the observational data, scientific discussions and text writing. T.A. contributed to the analysis of the data and scientific discussions. A.V. contributed to software developments for data analysis and scientific discussions.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** is available for this paper at https://doi.org/10.1038/s41598-020-69916-w.

**Correspondence** and requests for materials should be addressed to V.C.

**Reprints and permissions information** is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_10.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International

License, which permits use, sharing, adaptation, distribution and reproduction in any medium or
format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the
Creative Commons license, and indicate if changes were made. The images or other third party material in this
article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the
material. If material is not included in the article's Creative Commons license and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2020