Our understanding of geomagnetic field intensity prior to the era of direct instrumental measurements relies on paleointensity analysis of rocks and archaeological materials that serve as magnetic recorders. Only in rare cases absolute paleointensity datasets are continuous over millennial timescales, in sub-centennial resolution, and directly dated using radiocarbon. As a result, fundamental properties of the geomagnetic field, such as its maximal intensity and change rate have remained a subject of lively discussion. Here, we place firm constraints on these two quantities using Bayesian modeling of well-dated archaeomagnetic intensity data from the Levant and Upper Mesopotamia. We report new data from 23 groups of pottery collected from 18 consecutive radiocarbon-dated archaeological strata from Tel Megiddo, Israel. In the Near East, the period of 1700–550 BCE is represented by 84 groups of archaeological artifacts, 55 of which were dated using radiocarbon or a direct link to clear historically-dated events, providing unprecedented sub-century resolution. Moreover, stratigraphic relationships between samples collected from multi-layered sites enable further refinement of the data ages. The Bayesian curve shows four geomagnetic spikes between 1050 and 600 BCE, with virtual axial dipole moment (VADM) reaching values of 155–162 ZAm2 – much higher than any prediction from geomagnetic field models. Rates of change associated with the four spikes are ~0.35–0.55 μT/year (~0.7–1.1 ZAm2/year), at least twice the maximum rate inferred from direct observations spanning the past 190 years. The increase from 1750 BCE to 1030 BCE (73 to 161 ZAm2) depicts the Holocene’s largest change in field intensity.

Observational records of rapidly varying magnetic fields strongly constrain our understanding of core flow dynamics and Earth’s dynamo. Archeomagnetic analyses of densely sampled artefacts from the Near-East have suggested that the intensity variation during the first millennium BC was punctuated with two geomagnetic spikes with rates of change of intensity exceeding 1 μT/y, whose extreme behaviour is challenging to explain from a geodynamo perspective. By applying a new transdimensional Bayesian method designed to capture variations on both long and short timescales, we show that the data considered only at the fragment (thermal-unit) level require a complex intensity variation with six spikes, each with a duration between ~30-100 years. However, the nature of the inferred intensity evolution and the number of spikes detected are fragile and highly dependent on the specific treatment of the archeomagnetic data. No spikes are observed when the data are considered only at the level of a group of fragments from the same archeological context, with a minimum of three different artefacts per context. Furthermore, the number of spikes decreases to zero when increasing the error budget for the intensity within reasonable levels of 3-6 μT and the data age uncertainty up to 50 years. Thus, depending on the choices made, the Near-Eastern data are compatible with a broad range of time-dependence, from six spikes at one extreme to zero spikes on the other, the latter associated with much more modest rates of change of ~0.2-0.3 μT/y, comparable to secular variation at other periods and in other regions.