We analyse data^{26} obtained from the latest geomagnetic polarity time scale (2012), integrated with data from the late Devonian^{37}. The time sequence of polarity reversals was reconstructed for the whole Phanerozoic eon^{26}. The time sequence of reversals and the time evolution of the rate of reversals *γ*(*t*), obtained by averaging reversals over sliding windows of width 8 Myr, are reported in Fig. 1. We analysed the geomagnetic series *γ*(*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^{38} 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 *γ*(*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 *ω*_{j}(*t*), whose time variations describe the non-stationary processes underlying the observed variability of the reversal rate, and the time average of *ω*_{j}(*t*) allows us to define a typical period *T*_{j} = 2π/〈*ω*_{j}〉 (where the time average is denoted by angular brackets) associated to each *j*-th mode.

The probability density functions of the instantaneous periods 2π/*ω*_{j}(*t*) (see Fig. S1 of Supplementary information) indicate a compatibility with previously reported cycles. Periods shorter than about 40 Myr^{21,22,24,26,39} are found in the modes with *j ≤ *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^{27,30}. 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* + *σ 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*, *σ* is the noise level and *dw* is a Wiener process, i.e., a stochastic process with independent Gaussian increments, which describes the 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*_{o} + *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 ≤ *j* ≤ 4}, and double-well potentials, related to low-frequency CMB changes, for the set of modes *L* = {5 ≤ *j* ≤ 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 ≤ *j* ≤ 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 chrons towards low-frequency superchrons as a kind of phase transition^{44}, that we assume to be driven by stochastic fluctuations of the heat flux at the CMB. High-frequency chrons 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 = left{ {C_{j}^{sigma } left( t right)} right}) (see Supplementary Information) and let us write up the potentials in terms of the manifold

$$ Uleft( Gamma right) = rGamma^{2} + uGamma^{4} – hGamma $$

(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*_{1}, *t*_{2}) of the observed reversal rates for each mode

$$ Gleft( {t_{1} ,t_{2} } right) = Gamma left( {t_{1} } right)Gamma left( {t_{2} } right) approx exp left[ {frac{{ – left| {t_{1} – t_{2} } right|}}{tau }} right], $$

(2)

(where the angular brackets denote time averaging) shows a power law dependence *τ* ≈ 1/*r*^{α} and thus diverges when *Q* ≈ *Q*_{c} (see Fig. 4), where *α* = 0.48 ± 0.03, in close agreement with the scaling exponent *α* = 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.

According to the mean field approximation of second-order phase transitions, the susceptibility *χ* is given by *χ*^{−1} = *dh*/*dΓ*. The equilibrium solution in the mean field approach is obtained by minimising the potential *U* (Eq. (1)) with respect to *Γ.* This gives *h* = 2*rΓ* + 4*uΓ*^{3}, from which *χ*^{−1} = 2*r* + 12*uΓ*^{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 *Γ*^{2} = 0 when *r* > 0, therefore *χ*^{−1} = 2*r.* 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 *ΔQ ≈ Δh*, we get *ΔQ ≈ χ*^{−1}*ΔΓ ≈* 2*r ΔΓ*. By estimating *ΔΓ* as *ΔΓ ≈* 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} }} = Aexp left[ {frac{Delta n}{{n_{0} nu }}} right] $$

(3)

where *Δn* is the variation in the reversal rate and the best fit parameters result *ν* = 2.08 ± 0.03 and *A* = 0.07 ± 0.01. This means that fluctuations of the order of about *ΔQ ≈* 0.18 *Q*_{0} of CMB heat flow should be enough to double the reversal rate in the geodynamo system.