Modeling geomagnetic spikes: the Levantine Iron Age anomaly

The Levantine Iron Age anomaly (LIAA) is a regional short-decadal geomagnetic strength field variation located at the Levantine region characterized by high intensities with maximum virtual axial dipole moments around 190 ZAm2. It has been constrained by archeomagnetic data coming from Eastern Europe and Western Asia between 1050 and 700 BC. The LIAA can be related to a fast and spatially localized geomagnetic positive anomaly (spike) at the Earth’s surface. In this study, we model the LIAA by using a Fisher–von Mises function that fits the most recent archeomagnetic intensity database in the region. A spherical harmonic analysis is implemented for this spike function to perturb a base model in order to build a global reconstruction (perturbed-model) that reproduces the spatial and temporal characteristics of the LIAA. Our results show the importance of harmonic degrees from n = 3–4 to n = 20 to reconstruct the anomaly extension suggested by the database. Two maxima linked with the LIAA are reproduced by our global perturbed-model at the Levantine region at 950 BC and 750 BC. A third maxima in intensity around 500 BC is also observed, affecting the whole Europe.


Introduction
Fast variation events associated with high intensities have been observed in the past geomagnetic field and labeled as geomagnetic spikes (e.g., Ben-Yosef et al. 2009;Shaar et al. 2011Shaar et al. , 2016;;Hong et al. 2013;Cai et al. 2017;Rivero-Montero et al. 2021; among others).However, the lack of these features in the present geomagnetic field (Jackson et al. 2000;Alken et al. 2021) makes their description and understanding of their origin a challenge in the study of geodynamic processes occurring in the Earth's outer core (e.g., Livermore et al. 2014;Terra-Nova et al. 2016;Sánchez et al. 2016).The continuous updating of the paleomagnetic and archeomagnetic database, which is essential to better constrain these rapid intensity anomalies in time and space, allows us to currently consider making a mathematical model that reproduces its characteristics on the Earth's surface.
The first geomagnetic record of these "geomagnetic spikes" was observed in the Levantine region around 900 BC (Ben-Yosef et al. 2009;Shaar et al. 2011), and named as the Levantine Iron Age anomaly (LIAA, Shaar et al. 2013Shaar et al. , 2016)).This event was initially characterized by a local intensity maximum associated to extremely high values of the virtual axial dipole moment (VADM) up to 190 ZAm 2 (note that the global mean VADM for the same time period was around 110 ZAm 2 according to the GEOMAGIA database of Brown et al. 2021).High-intensity records were also observed from Turkey (Ertepinar et al. 2012) and Georgia (Shaar et al. 2013) that confirmed this short-lived high intensity anomaly affecting these regions.Subsequent studies in the Levantine region allowed Shaar et al. (2016) to characterize the LIAA as a period of generally high field lasting around 350 years (∼ 1050 BC to ∼ 700 BC) in which at least two short-lived impulses of the geomagnetic field occurred.
In addition, high VADM values (around 130 ZAm 2 ) prior to the LIAA have also been recorded in eastern Asia (as Korea, Hong et al. 2013;and China, Cai et al. 2017).In the American continent a geomagnetic spike might have been recorded in two sedimentary sequences in Texas, United States of America (Bourne et al. 2016) although absolute intensity data with low VADMs from Mexico do not confirm this spike (Hervé et al. 2019).
Efforts have been done to explain the origin of these geomagnetic spikes in terms of geomagnetic field models.Livermore et al. (2014) indicated that the secular variation reported in the Levantine region is not compatible with the commonly accepted core-flow dynamics.Davies and Constable (2017) suggested that the LIAA could be a narrow intense flux patch in the core-mantle boundary (CMB) that grew in place and migrated northwestwards.On the other hand, Korte and Constable (2018) concluded that the explanation of a drifting intense flux patch is not compatible with other field observations and suggested that the spikes are linked to individual flux patches growing and decaying in situ.Osete et al. (2020) from a global archeomagnetic reconstruction, the so-called SHAWQ-Iron Age model, links the LIAA with a normal polarity flux patch at the CMB located bellow the Arabian Peninsula around 1000 BC that expanded towards the north-west, while decreasing in intensity, affecting Europe around 700-500 BC and vanishing in situ at around 100 BC.
However, most of the current paleomagnetic global reconstructions have difficulties to reproduce the LIAA in detail.An extreme example is the SHA.DIF.14kglobal model (Pavón-Carrasco et al. 2014) that was constructed using an archeomagnetic and volcanic database previous to the LIAA data publications and thus shows a flat strength field structure during the LIAA event (see Fig. 1).Other paleomagnetic reconstructions such as CALS10k.2 (Constable et al. 2016) or SHAWQ-Iron Age (Osete et al. 2020) that were built using the new LIAA data show a wide maximum in the paleointensity predictions in the Levantine region around 950 BC.But they smooth the anomaly (see Fig. 1) and cannot appropriately represent the LIAA as a short wavelength peak as suggested by the data.
The limitation of paleomagnetic reconstructions for reproducing small spatial wavelength features, as a spike, is related with its mathematical nature, the small amount of paleomagnetic information and the dating uncertainties.In spherical harmonic modeling, the spatial wavelength of a geomagnetic feature at the Earth's surface depends on the spherical harmonic degree n by, e.g., Kono (2015): where is the wavelength, a the Earth's mean radius and n the harmonic degree related to the geomagnetic poten- tial.Current spherical harmonic paleomagnetic reconstructions implement the paleofield till degree n = 10 , thus the minimum wavelength they are able to reproduce (using Eq. 1) is around 35° or 4000 km in surface (that is equivalent to an intensity feature-half wavelength-of width around 17° or 2000 km extent).In addition, the high dispersion of paleomagnetic and archeomagnetic data forces the use of a regularization approach in the paleomagnetic reconstructions that smooths the spatial and temporal variability of the models to better reproduce the real variability of the geomagnetic field.In this context, recent works have shown how the existence of geomagnetic spikes is also highly dependent of the error budget of the database (Livermore et al. 2021).
To avoid the limitations of global harmonic reconstruction models, Davies and Constable (2017) approached the problem in an interesting and innovative way.They created a synthetic spike field by using a circular Fishervon Mises probability density function.Then, the geomagnetic field at the Earth surface (or at the CMB) is the spike field plus a global harmonic reconstruction field model.These authors analyzed the global archeomagnetic database up to 2017 and concluded that the LIAA is a very thin anomaly of longitudinal width around 20° on surface.In their results the strength field shows no spike-like feature when power spectrum is calculated till degree 14 on Earth's surface but on the other hand, there is a significant power at harmonic degrees over 20.It is worth noting that this extension represented by harmonic degrees n > 20 is not related with an outer- core origin feature.Also these authors proved that if the source of the LIAA is a spike of less than 1° width on the CMB it must span till 60° in longitude on the Earth`s surface, which is not consistent with the archeomagnetic data available till 2017.
However, recent archeomagnetic contributions (Hervé et al. 2017;Osete et al. 2020;Rivero-Montero et al. 2021) suggest a wider structure for the spike that could be related with lower harmonic degrees (see Eq. 1).This motivates us to revisit this issue implementing a paleomagnetic reconstruction based on the most recent archeomagnetic intensity database.In this study, following the work of Davies and Constable (2017), we analyze the evolution of the LIAA as a perturbation of the radial component of the geomagnetic field.To do that, we first develop a synthetic harmonic model for the spike using the most updated archeomagnetic database in Europe and Western Asia.Then, we implement the previous spike model as a perturbation into a paleomagnetic global reconstruction to generate a new paleo-reconstruction that allows analyzing the spatial and temporal behavior of LIAA.

Data selection and weighting scheme to constrain the LIAA
To constrain the spike associated to the LIAA, we use the archeointensity data from the GEOMAGIA50V3.4database (Brown et al. 2021, https:// geoma gia.gfz-potsd am.de/) and the new data published up to middle of 2021.A first spatio-temporal selection was applied for the archeomagnetic data representing the spike from the whole database.Data selected come from longitudes between 20° W and 80° E and latitude between 10° N and 60° N (see Fig. 2).To select the age range of the data, we choose moving temporal windows of 100-year width centered every 50 years from 1150 to 400 BC, considering all the data in each window as contemporary.After applying the spatial and temporal filter, we have a total of 646 archeointensity data in the region (see Fig. 2).Additional file 1: Table S1 contains these data with the original references.
The available archeointensity data are classified depending on their quality.Following the classification given by Campuzano et al. (2019), we consider quality data (Q) and non-quality data (no-Q).The criteria used to determine the quality is based on the number of specimens for an accurate intensity mean (at least three specimens for Q data) and the laboratory protocol used to obtain the intensity.Q data are based on Thellier-Thellier type methods (Thellier and Thellier 1959) with pTRM checks and TRM anisotropy correction when needed.About the 39% of the selected data base corresponds to Q data and thus 61% of the total number of data represents the no-Q data (see histogram in Fig. 2).This provides, in each 100 years interval, an average of 77 Q data and 98 no-Q data.Since the number of Q data is not enough for our modeling purpose, we use both Q and no-Q data for a better coverage, but associate different weights.The temporal distribution is nearly homogenous, although it shows a tendency of more available data for recent ages, with a maximum of 272 data in [600 BC, 500 BC].Spatially, the final data base is not homogenously distributed within each 100-year time window (see Fig. 2a-h), with higher data density in Europe and Western Asia.
The reliability and uncertainty in the age and intensity value of each archeomagnetic datum have to be also considered in our model.For that, we design a weighting scheme for our data that involves three weights, associated with age uncertainty W t i , for each datum i, meas- urement uncertainty ( W F i ) and data quality ( W Q i ).The total weight is therefore the square root of the sum of the squares of each weight.
The age of the data is given by an uncertainty interval.The mean value for the data age uncertainty is around 90 years for the whole period with a maximum uncertainty of 250 years, so for 100-year width moving windows some data belong to different windows.In order to consider the age uncertainty contribution, we follow the methodology detailed in Pavón-Carrasco et al. (2009) where a normalized weight W t i is calculated as the fraction of the age uncertainty interval that falls into the temporal window over the total length of the age interval.Also the age of the data is relocated in order to belong to the window (when the mean value of the datum age falls outside the window).
With regard to the measurement uncertainty σ F i , we also use a normalized weight relative to the minimum value of this parameter.Following previous works (e.g., Campuzano et al. 2019), we have fixed a minimum value of 4 µT for the intensity uncertainty ( σ F MIN ).Then, the weight associated with the measurement uncertainty is given by: In addition, an extra weight is considered based on the quality of the archeointensity data.According to the classification of Campuzano et al. (2019) in Q data and no-Q data, we fix a weight W Q i of 10 for Q data and 1 for no-Q data.With this weight, we try to highlight the quality as the main contribution for the data weighting. (2) Thus we have three different weights, W t i and W F i take values between 0 and 1, and W Q i is 1 or 10.Then, the final weight for each datum is given by: and falls into the range between 1 and √ 102.

Fitting the data to a spike function
Following Davies and Constable (2017), the LIAA can be mathematically defined on the Earth's surface (radius r = a ) by a spike function in the radial component of the geomagnetic field.This fact is supported because this geomagnetic feature is associated with local maxima in intensity and inclination (see e. (3) Fig. 2 Spatial (a-h) and temporal (i) distributions of the archeomagnetic database used in this work.Data in orange represent the high quality archeointensity data (Q data) and green correspond to lower quality data (no-Q data).Note that some data belong to more than one temporal window due to the age uncertainty.The mean age of some data might be outside the window, but a fraction of the age uncertainty interval can belong to the window where and A is the amplitude parameter of the spike, θ c and ϕ c are the colatitude and longitude of the center of the spike and k is the concentration parameter of Fisher-von Mises distribution, that relates with the angular standard deviation σ (measured in °) by k = 81 2 /σ 2 .The second term in Eq. 4 is added to avoid the creation of a monopole which violates the Maxwell equation constraint ∇ • B = 0 (see Davies and ConsTable 2017 for more details).
Equation 4 needs to be approximated when k reaches high values because sinhk → ∞ and thus provides infinity values of B r (e.g., k > 700, B r is not defined).More informa- tion about this issue is given in Additional file 2: Appendix A1, Fig. S1.To solve this computing problem, we propose an approximation of Fisher-von Mises probability distribution f (θ, ϕ) = Ak e kcosβ 4π sinhk for high values of k: where we consider the definition of sinh(k) = e k −e −k 2 , thus, where Eq. 7 provides a good approximation for a Fishervon Mises distribution for high k.
Assuming an insulating mantle and considering the previous function (Eq.7), the geomagnetic field at the Earth's surface can be described by the sum of a geomagnetic field global model that does not take into account the LIAA and the previous spike function.In this study, we only use the radial field for the spike function and thus, the total radial geomagnetic field can be given by: where B base r is the radial field of the global model and B spike r is the radial field given by Eq. 4. Thus, the intensity of the geomagnetic field can be calculated and approximated as: (4) (5)  Carrasco et al. 2014).This model was chosen because it was implemented without the new high-intensity data from the Levantine region published after 2014 (Shaar et al. 2016(Shaar et al. , 2017)), consequently, this global model does not reproduce the LIAA event as seen in Fig. 1.This flat shape in the Levantine region given by the model makes it appropriate for adding the perturbation B spike r of Eq. 8.
Here, it is important to note that we are assuming that the horizontal components of the spike field ( B spike θ and B spike ϕ ) are approximated to zero in Eq. 9.This approximation simplifies our approach but we have to check if it is suitable.To do that, we perform a test using the perturbed-model (see next section) that considers both vertical and horizontal components of the spike.Then, we compare the intensity provided by the two terms of Eq. 9. Results (see Additional file 2: Appendix A2 and Figs.S2, S3 for details) indicate that the horizontal components of B spike do not play an important role in our approach and therefore, they can be considered as negligible in Eq. 9.
In order to fit the selected archeointensity data base by Eq. 4 to get the optimal parameters for the spike function, we implement a forward approach based on the Monte Carlo method and using a total of 2 × 10 5 iterations for each 100-year width temporal window, centered every 50 year from 1150 to 400 BC.In this forward approach, the parameters A, k, c , ϕ c (where c and ϕ c represent the latitude and longitude, respectively, of the center of the spike) of Eq. 4 vary simultaneously following a homogenous distribution.For each iteration, we first estimate the geomagnetic field elements ( B base r , B base θ , B base ϕ ) from SHA.DIF.14kmodel, calculated at the same location and time of the archeointensity data base.Then, the radial field of the spike function B spike r is derived from the random parameters A, k, c , ϕ c in a previously fixed interval.A previous test indicates the following ranges for each parameter: A ∈ (−5000 nT, 0 nT), k ∈ (0.200) and the center of the spike in c ∈ (20 o S, 80 o N), ϕ c ∈ (10 o E, 50 o E) .With the sum of both B r elements (Eq.8), we estimate the intensity of the geomagnetic field (Eq.9) that is compared with the archeointensity values.For each time window, the minimum value of the weighted root mean square (RMS) error between model and real data provides the best set of random parameters of the spike function.The weighted RMS error is given by: where F i represents the measured intensity given by the data, F i is the intensity reproduced by the model in the location of the data and W final i is the final weight for each datum as detailed in "Data selection and weighting scheme to constrain the LIAA" Section.
It is worth noting that the high dispersion of archeointensity data (see Fig. 1) can hide the spike maxima.Since the objective of this work is to describe the spike as an extreme geomagnetic event characterized by an intense and spatially narrow peak (rather than provide a new paleomagnetic global paleoreconstruction), we minimize the previous RMS error of Eq. 10 only using the higher archeointensity data that records the spike structure.Consequently, to get the optimal parameters in the direct approach, we only use the intensities higher than the weighted mean intensity value for each temporal window.In other words, we force the appearance of a short wavelength spike that fits the highest intensity data values.In contrast, the contribution of the data with intensity lower than this mean value is not considered in this fitting.They are probably well represented by the SHA.DIF.14kbase model.

Spherical harmonics fitting for the spike function: the spike harmonic model
Once obtaining the best parameters that define the spike function, a spherical harmonic analysis is performed to express the spike function in the same frame used for the global paleoreconstructions. Since the geomagnetic potential field can be represented as an expansion of harmonic functions, we can easily obtain the expression of the radial field derived from the spherical harmonic expansion of the potential as follows: where ( r, θ, ϕ ) are the spherical coordinates, g m n and h m n are the Gauss coefficients and P m n (cosθ) the Schmidt quasi-normalized Legendre functions, n is the harmonic degree and m the harmonic order.
The optimal parameters for B spike r are distributed between 1150 and 400 BC every 50 years.To get better resolution for the inversion, we apply a spline interpolation of the optimal parameters of B spike r to 25 year.Then, in order to obtain the Gauss coefficients of B harmonic−spike r , (10) we implement a temporal inversion using cubic B-splines (De Boor 2001): where P is the number of splines defined by the knot points (since the splines are defined by cubic polynomial pieces, P = 4).
Using a matrix representation, where g is the column array whose elements are g m n,p , E the data array and M the matrix that let Eq. 13 satisfy Eq. 11.The least-squares solution, the one that minimizes the sum of the squares of the residuals, is given by where M T is the transpose of M.
The goal is to obtain a spherical harmonic expansion for B spike r defined by Eq. 4, the named spike harmonic model.The input data for the temporal inversion are the B spike r on a 500-point regular grid over the Earth's surface every 25 years.To get a better resolution in the spike area, a denser regular grid (three times denser than the global grid) is considered only in a spherical cap of 30° centered at 15° N, 40° E, that is, covering the LIAA region as well as most Europe and Western Asia (see Additional file 2: Fig. S4).
The B spike r values given by fitting the parameters of Eq. 4 to the data, as explained in "Fitting the data to a spike function" section, have high values only in the region sampled by the denser grid and are close to 0 nT outside the region of the spike.Since these values are distributed homogeneously in both space and time, we do not have imposed any kind of temporal and spatial regularization in the fitting procedure (Eq.14).To get the solution for the inversion approach, the series representation of the radial component of the magnetic field given by Eq. 11 has to be truncated for a certain harmonic degree.In order to define the maximum harmonic degree for the expansion given in Eq. 11, we have performed different trials with n max between 15 and 20 but they reproduce B spike r with similar results.However, since the harmonic degree n = 20 provides shorten spatial wavelengths that better fit the geometry of the spike, we prefer to keep this maximum expansion for our spike harmonic model.( 12) Consequently, the spike harmonic model is developed up to harmonic degree 20 with temporal knot points every 25 years.
Assuming that the uncertainty of the optimal parameters for B spike r generated in the previous section follows a Gaussian distribution, we implement an iterative bootstrap method in the inversion to obtain an ensemble of the spike Gauss coefficients.In each iteration a Gaussian distributed random value of A, k, λ c , ϕ c in the range of their uncertainties is used to estimate the B harmonic−spike r by the inversion approach of Eq. 11.We carry out this procedure 5000 times, changing randomly the parameters for each iteration.The final spike harmonic model is obtained by the mean of the ensemble of 5000 set of Gauss coefficients and their uncertainties are given by standard deviations.

The spike function: a forward approach
The Monte Carlo forward approach allows to obtain the optimal parameters for each temporal window of the spike function that better fits the archeointensity data (both Q and no-Q data with different weight).Our results in Fig. 3 show the evolution of the optimal parameters and its uncertainty given by the error bars.
The intensity of the spike function, given by the parameter A, remains constant around 5 µT (Fig. 3a).However, A does not give a full information about the anomaly intensity because, as shown in Eq. 4, the amplitude of the Fisher-von Mises distribution function depends on A and k.The minimum value of the parameter k (Fig. 3b) is k min ≈ 19 for the vanishing of the anomaly in 450 BC and the maximum value k max ≈ 52 is reached for 1000 BC.We note the presence of another smaller maximum (k ≈ 25) around 600 BC.Since the parameter k is inverse related to the standard deviation, we can estimate a maximum width of about 4150 km for the LIAA beginning and a minimum width around 2500 km for the main intensity maximum in 1000 BC and around 3600 km for the second maximum in 600 BC.The evolution of the concentration parameter k with time (Fig. 3b) shows that the LIAA starts out as a wide structure, reducing its width till 1000 BC, where the anomaly has a smaller spatial wavelength, and starting growing in width till 800 BC.From then, the width stays more stable with a small decrease at 600 BC which corresponds with the second maximum in the parameter k.
The optimal location of B spike r is given by the coordinates of its center (see Fig. 3c, d).The spike emerges for the age period considered in this study, from 1150 to 400 BC.The optimal parameters are calculated for each temporal window of 100 years moving every 50 years.The error bars indicate the uncertainty of the optimal parameters at 1σ of significance around 40° N, 40° E, and reaches its maximum intensity at 35° N, 35° E between 1050 and 1000 BC.During the LIAA period (that is between 1050 and 750 BC) the center of the spike does not experience important variations, perhaps a slight variation in longitude of about 4° westward.After the LIAA times our model shows a slightly drifting to lower longitudes of the B spike r 's center till 500 BC, reaching longitude values around 15° E.
The maximum intensity of the anomaly can be calculated in the center of the spike, where cos β = 1.The intensity of B spike r shows an expected result in Fig. 3e, with a fast increasing tendency reaching a maximum of 40 µT around 1000 BC, where the LIAA is more intense and thinner.After that the intensity of B spike r decreases and shows a smaller maximum in 600 BC of around 20 µT.A drop in intensity of 10 µT is observed between 800 and 750 BC.
To quantify how the spike function reproduces the archeointensity data, we plot a histogram with the residuals between input and modeled data (see Additional file 2: Fig. S5).Residuals seem to be normally distributed centered in 0 with a small deviation to negative values which indicates that the model fits better the high-intensity values.This asymmetry pointed out by the residual histogram was expected, because we have imposed this bias in forward approach using only the higher archeointensity data to highlight the spike spatial structure (details in "Fitting the data to a spike function" section).
The optimal parameters with its uncertainty for the Fisher-von Mises function are also detailed in Additional file 3: Table S2.
To better represent the spatial and temporal evolution of the spike, we plot in Fig. 4 different snapshot maps of the spike function at the Earth's surface.Results show the anomaly seems to appear in the north-east of the Anatolian Peninsula around 40° N and moves slightly towards southern latitudes till the Levant Region where it reaches its maximum intensity at 950 BC.From 950 up to 750 BC, the position of the LIAA maximum remains in the Levant region, while the intensity maximum expands and its intensity decreases.A slight reactivation in intensity is observed around 650 BC associated with a slight migration of the center of the spike towards the west, reaching the south part of the Italian Peninsula around 650 BC.The center of the spike moves to the South of Greece, where it practically vanished by 450 BC.Note that the lack of data coming from North of Africa (see data base in Fig. 2) can constrain the LIAA boundary in this area only poorly, implying uncertainty about the described extension and drift.
In terms of spike intensity, results show a notable intensity oscillation, with a first intensity maximum around 950 BC (Fig. 4c) that is related with the LIAA event; and a slight reactivation around 650 BC that could be related with the maximum intensity values observed in Europe and eastern Asia (Osete et al. 2020;Rivero-Montero et al. 2021).

Spherical harmonic model of spike function: the spike harmonic model
In order to analyze the features of the LIAA event, we implement an inverse approach to obtain a spherical harmonic model for the B spike r by using synthetic data from the spike function as input data (see details in "Spherical harmonics fitting for the spike function: the spike harmonic model" section).Since the inversion is carried out using synthetic data from the spike function, we can implement the expansion in spherical harmonics till the chosen maximum degree n = 20 without risk of unreal oscillations.The Gauss coefficients for the spike harmonic model, B harmonic−spike r , are provided in Additional file 4: Table S3.
In Fig. 5, we plot a latitudinal profile (along the constant parallel that crosses the center of the spike) of the spike radial component used as input data of our inverse model and that given by the radial component of the spike harmonic model (obtained from harmonic degree n = 1 to n = 20).A general view shows a good agreement between the spike harmonic model and synthetic spike function.
To evaluate, in terms of harmonic degrees, the average energy of the spike harmonic model at the Earth's surface we use the power spectra (PS, Lowes 1974), given by Figure 6 provides the calculated PS of the spike harmonic model every 100 years (in color lines) from 1150 to 450 BC and the time-averaged power spectra (dashed gray line).As it can be seen, for the epochs around the LIAA event (i.e., from 1050 to 850 BC) the PS related to harmonic degrees over n = 3 show higher values than all the other PS values (outside the previous temporal windows) and the averaged PS.From n = 7, this difference becomes more notable, showing the importance of the harmonics between 7 and 20 to define the LIAA as a small wavelength and high-intensity feature of the paleofield.The PS calculated in the CMB is provided in Additional file 2: In order to obtain a new global model that reproduces the evolution of the LIAA, we consider a perturbation of the base model, adding the coefficients of the spike harmonic model (given in the previous "Spherical harmonic model of spike function: the spike harmonic model" section) to the SHA.DIF.14kcoefficients (developed from the dipole to degree 10) for the time range where the spike harmonic model has been obtained.Consequently, the new global paleoreconstruction is expanded up to harmonic degree n = 20.The uncertainties of both sets of Gauss coefficients (spike harmonic and SHA.DIF.14kmodels) have been considered, using a bootstrap method, to get the final set of Gauss coefficients uncertainties of the perturbed-model (coefficients and uncertainties are provided in Additional file 5: Table S4).
In order to illustrate a more clear vision of the geomagnetic field anomaly created by the spike, we plot in Fig. 7 the archeointensity data and the intensity of the magnetic field given by the perturbed-model.The intensity model fits the high archeointensity data for the region and period of study.The figure clearly shows how the LIAA spike has a localized structure in space and time (between 1050 and 850 BC), while the second maximum (around 600 BC) is wider and of lower intensity.
In addition, using this new global paleoreconstruction we can generate paleosecular variation curves (PSVC) at any point on the Earth's surface.In Fig. 8 different PSVCs at three different locations along the Mediterranean (coordinates marked with the yellow stars) are plotted covering the first two millennia BC.
In the Levantine region (Fig. 8a), the perturbedmodel (blue line with error bands) indicates an abrupt change in the intensity of the geomagnetic field (an increase of 25 µT) from 1200 to 1050 BC, reaching 75 µT, and a second and more intense pulse at 950 BC, when the field reaches its maximum value over 85 µT.A second maximum is observed at 750 BC (∼ 85 µT), after these maxima the paleofield values decreased up to 600 BC when the field was ∼ 65 µT.After that, a third relative maximum reaching a field value of 75 µT around 500 BC is detected.In summary, during the analyzed time period (when the perturbation has been added to the base model), the perturbed-model provides three main maxima located around 950 BC, 750 BC and 500 BC.The first two maxima (950 BC and 750 BC) correlate with the two events that characterize the LIAA (Shaar et al. 2016), while the third maximum is in agreement with the maximum observed by Osete et al. (2020) and Rivero-Montero et al. (2021).over the Earth's surface every 100 years.In pink-red lines, the PS before (1150 BC) and during LIAA maximum (1050 BC-850 BC) are plotted.In green-yellow the PS after LIAA maximum (750 BC-450 BC).The grey dash line shows the mean power spectrum during the whole period As expected, the largest differences between the perturbed-model and the base model SHA.DIF.14k(black line in Fig. 8) in the Levantine region are restricted during the period when the spike function is defined, i.e., from 1200 to 400 BC.In this period, the SHA.DIF.14kmodel presents low strength field values (note that the base model did not use the LIAA records) and only suggests the maximum of 500 BC, but about 10 µT lower than the perturbed-model.The SHAWQ-Iron Age model (red line in Fig. 8) clearly detects the three main maxima in coherence with our perturbed-model but it reaches lower intensity values with a difference ∼ 10 µT for the maxima of 950 BC and 750 BC and ∼ 8 µT for 500 BC.In addition, the maximum of 950 BC is observed by SHAWQ-Iron Age a bit earlier than the perturbed-model, around 1000 BC.
In Central Europe (Fig. 8b) the perturbed-model shows an increasing trend for the intensity from 1700 to 500 BC with a fast decay afterwards.During this increasing period the perturbed-model shows a wide maximum from 1150 to 950 BC when field reached values ∼ 65 µT.After that, there is a field decay of 10 µT and subsequently another wide maximum is detected between 750 and 500 BC.The perturbed-model provides intensities ∼ 80 µT during these maxima.The highest intensity values predicted by the perturbed-model are for the 500 BC maximum when field values are over 80 µT.The SHAWQ-Iron Age model also points out the three maxima at around 1000 BC, 750 BC and 500 BC.Again, due to the lack of data during the construction of the SHA.DIF.14kbase model, it is not able to reproduce the maximum of 1000 BC in Central Europe but observed a maximum at 500 BC (with lower amplitude, ∼ 15 µT, than the perturbedmodel) and an increase in intensity around 750 BC.
Finally, in Western Europe (Fig. 8c), the perturbedmodel seems to follow the SHA.DIF.14kbase model, perfectly matching outside the period between 750 and 400 BC, and reproducing a maximum around 500 BC but ∼ 10 µT more intense.The SHAWQ-Iron Age suggests three maxima of increasing intensity at around 1000 BC, 750 BC and 500 BC in Western Europe that is in agreement with the most recent Iberian intensity PSVC of Osete et al. (2020).However, the sequence of the first two maxima linked with the LIAA (1000 BC and 750 BC) are not predicted in this area by the perturbed-model due to two main reasons: (a) the chosen spike function presents a radial symmetry, and this geometry has not enough resolution to appropriately represent the geometry of these From a global perspective, we can summarize that the perturbed-model allows to describe the evolution of the intensity paleofield in the Mediterranean region during the first half of the first millennium BC.Three main maxima have been reconstructed.The maximum with higher intensity can be noted in the Levant around 1000-950 BC, the VADM reached over 160 ZAm 2 .This maximum has small longitudinal extension, so it can be noted in Central Europe with a decrease of 30% in VADM values but does not reach Iberia due to its small width.Another maximum in intensity is reproduced by the perturbed-model around 750 BC.This maximum is also noted in the Levant (VADM of 160 ZAm 2 ) but has a wider extension that makes it be clearly distinguished in Central Mediterranean (VADM of 150 ZAm 2 ).The peak of 500 BC has maximum intensity in Central Mediterranean and seems to have a bigger spatial extension, affecting the whole Mediterranean, from Iberia to Levant.The wide extension of this last maximum (Additional file 2: Fig. S7), affecting the whole of Europe, suggests that this feature is not related with the LIAA but a continental-scale feature of the geomagnetic field.These results are in agreement with previous studies.The two LIAA-linked peaks occurred around 1000 BC and 750 BC in the Levantine and surrounding areas (Ertepinar et al. 2012;Shaar et al. 2017;Rivero-Montero et al. 2021) and the maximum around 500 BC is considered not LIAA-linked, and widely extended as a maximum in intensity for the whole Europe from Canary Islands to Turkey, affecting probably Central Asia (Rivero-Montero et al. 2021).

Fig. 1
Fig. 1 Latitudinal profile (at 35° N) of the geomagnetic intensity field for a temporal window of 100-year wide centered at 950 BC given by the paleomagnetic models SHA.DIF.14k(Pavón-Carrasco et al. 2014), CALS10k.2 (Constable et al. 2016) and SHAWQ-Iron Age (Osete et al. 2020).In grey dots are represented the archeointensity data for the period [− 1000, 900] relocated to 35° N. Error bars show the uncertainty of the intensity data values g., Fig.8inPavón-Carrasco et al. 2021) where the vertical component (i.e., the radial field) takes an important relevance.Then, the anomalous field linked to LIAA should be dominated mainly by the vertical component.This function can be represented by a Fishervon Mises probability distribution in spherical coordinates for the radial field of the spike: geomagnetic field element (north and east elements, respectively).These horizontal elements along with the B base r element are synthetized by the SHA.DIF.14kglobal model (Pavón-

Fig. 3
Fig. 3 Optimal parameters of B spike r

Fig. 4
Fig. 4 Time evolution of the spike function B spike r at the Earth's surface in a spherical cap of 35° centered in 30° N, 30° E

FigFig. 5
Figure6provides the calculated PS of the spike harmonic model every 100 years (in color lines) from 1150 to 450 BC and the time-averaged power spectra (dashed gray line).As it can be seen, for the epochs around the LIAA event (i.e., from 1050 to 850 BC) the PS related to harmonic degrees over n = 3 show higher values than all the other PS values (outside the previous temporal windows) and the averaged PS.From n = 7, this difference becomes more notable, showing the importance of the harmonics between 7 and 20 to define the LIAA as a small wavelength and high-intensity feature of the paleofield.The PS calculated in the CMB is provided in Additional file 2: Fig.S6.(15)PS n = (n + 1) n

Fig. 6
Fig. 6 Power spectra (PS) of the modeled B harmonic−spike r

Fig. 7
Fig. 7 Intensity of the magnetic field in the Earth's surface calculated with the perturbed-model.The archeointensity data (in grey dots) used in this study for each 100-year width window are also plotted Fig.8a-c Secular variation of the geomagnetic intensity for the two first millennia BC reconstructed by the model SHA.DIF.14k(black), SHAWQ-Iron Age (red) and the new perturbed-model (blue) with a 95% confidence uncertainty interval.Yellow stars in the central map show the coordinates where the secular variation is plotted.In dots we plot the relocated archeomagnetic in a circle of 1000 km (orange for Q data and green for no-Q data)