Maxwell C Brown
email: mcbrown@umn.edu

Gwenaël Hervé

Monika Korte

Agnès Genevey

Global archaeomagnetic data: The state of the art and future challenges

Keywords: archaeomagnetism, GEOMAGIA50, global models 1

Archaeomagnetic data are fundamental for our understanding of the evolution of Earth"s magnetic field on centennial to millennial timescales. From the earliest studies of the Thelliers, Aitken, Nagata and others in the 1950s and 1960s, archaeomagnetic data have been vital for extending our knowledge of the field to times prior to observational measurements. Today, many thousands of archaeomagnetic data allow us to explore the geomagnetic field in more detail than ever before.

Both regional time series of archaeomagnetic data and the inclusion of archaeomagnetic data in time-varying global spherical harmonic field models have revealed a range of newly discovered field behaviour. More sophisticated approaches to developing regional curves and global models have allowed us to resolve the field in certain regions more robustly and with greater resolution than previously possible. In this review we give an overview of the widely used global

Introduction

Archaeomagnetism is the study of the past direction and intensity of Earth"s magnetic field recorded by any type of manmade artefact or fired material. It is dependent on archaeological discoveries and advances that lead to a better description and understanding of our history and heritage. Although it was recognized at the end of the 19th century that fired materials can record Earth"s magnetic field [START_REF] Folgheraiter | Sur les variations séculaires de l"inclinaison magnétique dans l"antiquité[END_REF], it was not until the pioneering work of Émile and Odette Thellier beginning in the 1930s that the physical principles, methods and instrumentation necessary to accurately obtain the past direction and intensity of the geomagnetic field recorded by archaeological materials were developed [START_REF] Thellier | Sur l"aimantation des terres cuites et ses applications geophysiques[END_REF][START_REF] Thellier | Sur les propriétés de l"aimantation thermorémanente des terres cuites[END_REF][START_REF] Thellier | Sur l"intensité du champ magnétique terrestre dans le passé historique et géologique[END_REF].

Reviews by [START_REF] Thellier | Early research on the intensity of the ancient geomagnetic field[END_REF], Le [START_REF] Goff | Emile Thellier (1904-1987), a pioneer in studies of the "fossil" Earth"s magnetic field[END_REF] and [START_REF] Dunlop | Physical basis of the Thellier-Thellier and related paleointensity methods[END_REF] give excellent English language overviews of the Thelliers" most important contributions to the subject.

Archaeomagnetism established itself as a research field through the 1950s and 1960s, with proponents of the subject obtaining data from fired materials from locations globally. Studies from these decades reported data from Europe (e.g., [START_REF] Burlatskaya | Archeomagnetic evidence for the Earth"s magnetic field near Tiflis during the past[END_REF][START_REF] Aitken | Magnetic dating: Some archaeomagnetic measurements in Britain[END_REF][START_REF] Belshé | Some archaeomagnetic results from Greece[END_REF][START_REF] Chelidze | Some results of a study of the past geomagnetic field in Georgian SSR by J o u r n a l P r e -p r o o f Journal Pre-proof the archeomagnetic method[END_REF][START_REF] Bucha | Intensity of the Earth magnetic field during archeological times in Czechoslovakia[END_REF][START_REF] Kovacheva | Inclination of the Earth"s magnetic field during the last 2000 years in Bulgaria[END_REF], Northern Africa [START_REF] Athavale | Intensity of the geomagnetic field in prehistoric Egypt[END_REF], India [START_REF] Athavale | Intensity of the geomagnetic field in India over the past 4000 years[END_REF], China [START_REF] Deng | The geomagnetic field in Peking region and its secular variation during the last 2000 years[END_REF], Japan [START_REF] Watanabe | Secular variation in the direction of geomagnetism as the standard scale for geomagnetochronology in Japan[END_REF][START_REF] Nagata | Secular varaition of the geomagnetic force during the last 5000 years[END_REF][START_REF] Sasajima | Geomagnetic secular variation revealed in the baked earths in West Japan. Part 2 Change of the field intensity[END_REF], North America (e.g., [START_REF] Watanabe | Some results of an archaeomagnetic study on the secular variation in the southwest of North America[END_REF][START_REF] Schwarz | Original remanent magnetization of Ontario potsherds[END_REF] and South America (e.g., [START_REF] Nagata | Archeomagnetic intensity studies of South and Central America[END_REF][START_REF] Kitazawa | Intensity variation of the geomagnetic field during the past 4000 years in South America[END_REF]. Research continued though the 1970s, but it was not until the 1980s that there was a general increase in the number of studies reporting new archaeomagnetic data each year (Fig. 1); a trend that continued through to the 2010s. This has resulted in a large compilation of global data that has greatly improved our understanding of how Earth"s magnetic field has varied spatially and temporally on centennial to millennial timescales. based on primarily archaeomagnetic data (but also including volcanic data) were not developed until the construction of ARCH3k.1 [START_REF] Korte | Geomagnetic field for 0-3 ka: 2. a new series of time-varying global models[END_REF] (a three thousand year model), which was recently updated to a 10,000 year model (Constable et al., 2016). Spherical harmonic cap approaches using archaeomagnetic data have also been used to create regional models (e.g., [START_REF] Pavón-Carrasco | Initial SCHA.DI.00 regional archaeomagnetic model for Europe for the last 2000 years[END_REF][START_REF] Pavón-Carrasco | A regional archeomagnetic model for Europe for the last 3000 years, SCHA.DIF.3K: Applications to archeomagnetic dating[END_REF]. Varying approaches to modelling the archaeomagnetic field have been applied since, including [START_REF] Licht | Ensembles of low degree archeomagnetic field models for the past three millennia[END_REF], [START_REF] Pavón-Carrasco | A geomagnetic field model for the Holocene based on archaeomagnetic and lava flow data[END_REF], [START_REF] Sanchez | Modelling the archaeomagnetic field under spatial constraints from dynamo simulations: a resolution analysis[END_REF], [START_REF] Hellio | Time-correlation-based regression of the geomagnetic field from archeological and sediment records[END_REF], [START_REF] Arneitz | A Bayesian iterative geomagnetic model with universal data input: Self-consistent spherical harmonic evolution for the geomagnetic field over the last 4000 years[END_REF] and [START_REF] Mauerberger | Correlation based snapshot models of the archeomagnetic field[END_REF]. Concurrent to regional compilations of data and the development of global models there have been continued efforts to create global databases of archaeomagnetic data. The first global archaeomagnetic databases were paper lists of results, the first likely being the historical and archaeointensity compilation of [START_REF] Smith | Ancient geomagnetic field intensities-I Historic and archaeological data: Sets H1-H9[END_REF]. With the development of digital database structures, archaeomagnetic data could be compiled and updated more easily. Early efforts included those of [START_REF] Burlatskaya | Archeomagnetic determinations of geomagnetic field elements of the USSR Academy of Sciences[END_REF], [START_REF] Liritzis | A global archaeomagnetic intensity data bank[END_REF] and [START_REF] Daly | An updated and homogeneous world secular variation data base. 1. Smoothing of the archaeomagnetic results[END_REF], although the data were not available in a digital form. The first digital archaeomagnetic database that was easily accessible was the Plymouth archaeomagnetic directional database (ARCHEO97 and ARCHEO00) compiled by Don Tarling and last released in 1999. This was one of seven International Association of Geomagnetism and Aeromony (IAGA) databases available online to download as stand-alone programs. Two major efforts to compile all global archaeomagnetic data have been the ArcheoInt database [START_REF] Genevey | ArcheoInt: An upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment[END_REF] and the GEOMAGIA50 database [START_REF] Donadini | Database for Holocene geomagnetic intensity information[END_REF][START_REF] Korhonen | GEOMAGIA50: an archeointensity database with PHP and MySQL[END_REF]Brown et al., 2015b). Although GEOMAGIA50 largely subsumes the data within ArcheoInt, ArcheoInt contains additional fields that place archaeomagnetic results in their archaeological context and provides greater descriptive information regarding the acquisition of the data sets. The databases can be viewed as complementary. In addition, there is the HISTMAG database of [START_REF] Arneitz | The HISTMAG database: combining historical, archaeomagnetic and volcanic data[END_REF], which combines historical and archaeomagnetic data. There are also numerous archaeomagnetic data in the MagIC database (described in part in [START_REF] Shaar | Large geomagnetic field anomalies revealed in Bronze to Iron Age archeomagnetic data from Tel Megiddo and Tel Hazor, Israel[END_REF]; however, GEOMAGIA50 is currently the primary database for archaeomagnetic data. Unlike MagIC, GEOMAGIA50 includes only average data and is not designed to include results at the specimen level or raw measurements.

The site level data from GEOMAGIA50 has been used in numerous studies. In addition to being used to construct secular variation curves and global and regional field models, it has been used to J o u r n a l P r e -p r o o f Journal Pre-proof understand solar activity during the Holocene [START_REF] Usoskin | Solar activity during the Holocene: the Hallstatt cycle and its consequence for grand minima and maxima[END_REF] and to calibrate cosmogenic nuclide production stacks through the use of intensity data (e.g., authigenic 10 Be/ 9 Be ratios, [START_REF] Simon | Authigenic 10 Be/ 9 Be ratio signatures of the cosmogenic nuclide production linked to geomagnetic dipole moment variation since the Brunhes/Matuyama J o u r n a l P r e -p r o o f Journal Pre-proof boundary[END_REF].

An important consideration when using archaeomagnetic data for any purpose is the reliability of the data. This includes chronological controls and archaeomagnetic components (direction and intensity), which are most commonly determined from a thermoremanent magnetization (TRM): a magnetization acquired on cooling from firing temperature to room temperature. Archaeomagnetic directions can be influenced by post-cooling displacement and magnetic refraction (section 3.1.1) and obtaining reliable archaeointensities requires that numerous factors are considered (section 3.1.2). These includes thermal alteration during palaeointensity experiments (section 3.1.3), the influence of non-ideal magnetic remanence carriers (e.g., multi-domain (MD) grains, section 3.1.4), remanence anisotropy (section 3.1.5) and differences between natural and experimental cooling rates (section 3.1.6). All chronological determinations have an associated uncertainty, whether an archaeological age, determined through physical measurements (e.g., by radiocarbon dating or luminescence methods), or by a combination of approaches (section 3.2). Documenting such uncertainties is a challenge (section 5.1.1) and uncertainties should be carefully considered in any study looking to investigate field behaviour.

In this review we cover the current status of the global archaeomagnetic database (GEOMAGIA50; section 2), provide an overview of archaeomagnetic procedures, data quality, uncertainties and chronological controls (section 3) and explore advances in regional secular variation curve construction and global archaeomagnetic field modelling (section 4). The review ends with a discussion on the future challenges of the subject (section 5).

Overview of the GEOMAGIA50 archaeomagnetic database

In the following sections we give a brief history of the GEOMAGIA50 database (section 2.1), cover the abundance of archaeomagnetic data within the most recent version of the database (GEOMAGIA50.v3.4) (section 2.2), discuss the spatial and temporal distribution of data (section 2.3 and section 2.4), and provide an overview of the archaeological materials used to obtain archaeomagnetic data (section 2.5). The methods used to obtain archaeomagnetic and age data, as well as their uncertainties, are discussed in section 3.

J o u r n a l P r e -p r o o f

Journal Pre-proof

In this review we consider purely archaeomagnetic data. Data from volcanic materials (lava, volcanic ashes, obsidian) and speleothems (i.e. [START_REF] Latham | The paleomagnetism and U-Th dating of Mexican stalagmite, DAS2[END_REF][START_REF] Trindade | Speleothem record of geomagnetic South Atlantic Anomaly recurrence[END_REF], although stored in GEOMAGIA50.v3.4, are neglected for the purpose of this study. We also restrict our analysis to materials dated between 8000 BCE and today, and we do not include materials that have been dated using archaeomagnetic dating. country-or region-specific compilations were also added (see Brown et al. (2015b) for a list of compilations). Further details of version 1 of the database can be found in [START_REF] Donadini | Database for Holocene geomagnetic intensity information[END_REF], [START_REF] Korhonen | GEOMAGIA50: an archeointensity database with PHP and MySQL[END_REF] and Brown et al. (2015b). After numerous updates since original publication, 2762 archaeomagnetic entries from 109 studies remain from version 1 in the most up-to-date version of the database.

History of GEOMAGIA50 and its most recent compilation

No publication accompanied version 2 of the database; however, the data compilation is described in [START_REF] Donadini | Geomagnetic field for 0-3 ka: 1. New data sets for global modeling[END_REF]. Around 100 archaeomagnetic entries from version 1 of the database were updated in version 2. Archaeomagnetic directional results were added independently of whether they accompanied intensity data. This greatly increased the amount of data in version 2 of the database, with 3072 data from 130 studies added at this time that remain in the most recent update of the database (5834 entries from 240 studies in total).

The current version of the database is version 3, which was initially published in 2015 (Brown et al., 2015b). It marked a change from being hosted at the Scripps Institution of Oceanography, University of California-San Diego, to GFZ Potsdam (https://geomagia.gfz-potsdam.de/). Sediment data were also added in version 3 (Brown et al., 2015a). 1006 entries from 100 archaeomagnetic studies were added to version 3.1 of the database;

498 entries from 220 studies were added to version 3.2 (released in 2017); and 1717 entries from 109 studies in version 3.3 (released in 2019). GEOMAGIA50.v3.2 also incorporated a number of legacy studies (studies published prior to the inception of the database in 2004) that were missing J o u r n a l P r e -p r o o f dating methods used. It also includes the number of samples/specimens and specimen types investigated, and the types of archaeological materials the data were obtained from (section 2.5).

Full details of the fields within GEOMAGIA50 are given in Brown et al. (2015b). There is a large disparity in the global distribution of archaeomagnetic data (Fig. 2). Data from Europe dominates the database (Fig. 3 and Fig. 4a-c): 59% of all entries (including Russia), 51% without Russian data. The UK (10% of entries), France (9%), Russia (8%) and Georgia (5%) contribute to a significant portion the European entries (Fig. 3 and Fig. 4b). Many European countries individually contribute between 2% and 4% of the total number of entries. The UK comprises the largest number of all entries (961), which are primarily directional data (905). See [START_REF] Batt | Advances in archaeomagnetic dating in Britain: New data, new approaches and a new calibration curve[END_REF] for further details on the UK contribution. France is the second largest contributor with 890 entries (770 with directions, 162 with intensity). A large amount of data was added (520 entries) following the publication of Le [START_REF] Goff | An updated archeomagnetic directional variation curve for France over the past two millennia, following 25 years of additional data acquisition[END_REF]. It is worth noting that although 72% of directional entries come from Europe (omitting Russia), this region covers only 1-2% of Earth"s surface (depending on the definition of Europe) (Fig. 4b). Data from regions adjacent to Europe are also dense, with the Levant (Israel, Syria, Jordan), Egypt and Iraq contributing significantly to the database. The distribution of data (directions and intensities) from Europe and these regions is shown in Fig. 3.

Spatial distribution of archaeomagnetic data

Outside of Europe the United States of America (7% of entries), China (6%), Japan (4%)

and Mexico (3%) are the main contributors. All other nations make up 29% of entries. Although the number of Japanese entries in GEOMAGIA50 totals 370, the Japanese archaeomagnetism database of T. Hatakeyama (Okayama University, Japan) (http://mag.center.ous.ac.jp/en) lists 744 directional data and 59 intensity data, placing it third in the list of country entries. We aim to integrate this significant contribution with GEOMAGIA50 in the future.

Although there have been recent efforts to improve the global distribution of data, the Southern Hemisphere is currently poorly represented, with only 400 entries or 4% of all archaeomagnetic entries. 76 entries contain a direction and 340 an intensity. The disparity in data distribution is stark when it is considered that Africa and South America, which cover 9% of Earth"s surface when combined (32% of the land area), provide only 7% of the entries in the database (Fig. 4a). However, the amount of Southern Hemisphere data continues to improve. In 2012), [START_REF] Tarduno | Antiquity of the South Atlantic Anomaly and evidence for top-down control on the geodynamo[END_REF] and Hare et al. (2018). Previously only one study had published data from this region [START_REF] Henthorn | An archaeomagnetic study of Mgungundlovu[END_REF] and this was not added to GEOMAGIA50 until version 3.3. A number of South American countries have garnered new data. In the first version of GEOMAGIA50, there were no entries from Brazil, Argentina, Uruguay and Chile. In the past 10 years data have been obtained from all four: Brazil, 49 intensity entries [START_REF] Hartmann | Archeointensity in Northeast Brazil over the past five centuries[END_REF][START_REF] Hartmann | New historical archeointensity data from Brazil: evidence for a large regional non-dipole field contribution over the past few centuries[END_REF][START_REF] Hartmann | New archeointensity data from South Brazil and the influence of the South Atlantic Anomaly in South America[END_REF][START_REF] Poletti | Archeomagnetism of Jesuit Missions in South Brazil (1657-1706 AD) and assessment of the South American database[END_REF]), Argentina, 44 entries (e.g., Gómez-Paccard et al., 2019;[START_REF] Goguitchaichvili | First archaeointensity reference paleosecular variation curve for South America and its implications for geomagnetism and archaeology[END_REF], Uruguay, 6 entries [START_REF] Capdepont | Archaeomagnetism applied to ceramics from coastal archaeological sites in Uruguay[END_REF] and Chile, 1 entry [START_REF] Roperch | Secular variation of the earth"s magnetic field and application to paleomagnetic dating of historical lava flows in chile[END_REF]. In addition, data have been obtained from other South American countries J o u r n a l P r e -p r o o f

Journal Pre-proof south of the Equator. By far the most number of entries come from Peru, with 191 (e.g., Gunn and Murray, 1980;[START_REF] Yang | A comparison of archaeointensity results from Chinese J o u r n a l P r e -p r o o f Journal Pre-proof ceramics using Thellier"s and Shaw"s palaeointensity methods[END_REF]. Smaller contributions come from Bolivia (13 entries) (e.g., [START_REF] Nagata | Archeomagnetic intensity studies of South and Central America[END_REF][START_REF] Kitazawa | Intensity variation of the geomagnetic field during the past 4000 years in South America[END_REF] and Ecuador (23 entries) [START_REF] Kitazawa | Intensity variation of the geomagnetic field during the past 4000 years in South America[END_REF][START_REF] Bowles | Archaeomagnetic intensity results from California and Ecuador: evaluation of regional data[END_REF][START_REF] Herrero-Bervera | First archaeointensity results from Ecuador with rock magnetic analyses and 14C dates to constrain the geomagnetic field evolution in South America: Enhancing the knowledge of geomagnetic field intensity[END_REF].

The area between the tropics fairs better than the Southern Hemisphere, with nearly 10% of all entries coming from this latitude band. This includes the large and growing data set from Mexico (see, Hervé et al., 2019b,c;[START_REF] Mahgoub | Refined 3600 years palaeointensity curve for Mexico[END_REF]. New studies from India [START_REF] Basavaiah | Last 3000 years of geomagnetic field intensity from India: New reference palaeointensity data from two east coast archaeological sites and archaeomagnetic dating insights[END_REF][START_REF] Deenadayalan | Rock-magnetic and archaeomagnetic investigations on archaeological artefacts from Maharashtra, India[END_REF], western Africa [START_REF] Kapper | Reconstructing the geomagnetic field in West Africa: First absolute intensity results from Burkina Faso[END_REF][START_REF] Kapper | Novel insights on the geomagnetic field in West Africa from a new reference curve (0-2000 AD)[END_REF] and eastern Africa [START_REF] Osete | First archaeomagnetic field intensity data from Ethiopia, Africa (1615  12 AD)[END_REF] have contributed important intensity data from areas that are isolated from others globally. As we move closer to the equator the amount of available data shrinks with <2% of database entries from between 10  latitude. Six studies have produced new data in this latitude band over the past 10 years, with the first archaeomagnetic data from Kenya [START_REF] Tchibinda Madingou | First archeomagnetic data from Kenya and Chad: Analysis of iron furnaces from Mount Kenya and Guéra Massif[END_REF] and the Ivory Coast [START_REF] Kapper | Novel insights on the geomagnetic field in West Africa from a new reference curve (0-2000 AD)[END_REF], and others building on small data sets from Ecuador [START_REF] Herrero-Bervera | First archaeointensity results from Ecuador with rock magnetic analyses and 14C dates to constrain the geomagnetic field evolution in South America: Enhancing the knowledge of geomagnetic field intensity[END_REF] and Colombia [START_REF] Cejudo | First archaeomagnetic results from Colombia (the Bogotá Savanna Pre-Hispanic sites): Implications for the Caribbean absolute geomagnetic intensity variation curve[END_REF].

The spatial distribution of directional and intensity data are distinctly different (Fig. 2, Fig. 3 and Fig. 4). Countries that produce numerous directional data do not always produce large amounts of intensity data and vice versa. As stated already, the UK has the most directional entries, but few intensity entries. Conversely, China has the most abundant intensity data by country, but does not make the top ten countries for directional data. Russia, east of the Black Sea has abundant directional data, but sparser intensity data. To a lesser extend the same is true for the Ukraine, which produces the 5th most directional data (5% of all directional entries), but far fewer intensity data than other countries. India and Brazil have no directional data, but numerous intensity data. This disparity can be crucial in areas with sparse data coverage, where full vector data are particularly important for constraining field models, e.g., sites in West Africa, where few directional data have been obtained (Burkina Faso; [START_REF] Donadini | Directional changes of the geomagnetic field in West Africa: Insights from the metallurgical site of Korsimoro[END_REF], whereas intensity data are more plentiful [START_REF] Mitra | Two thousand years of archeointensity from West Africa[END_REF][START_REF] Kapper | Reconstructing the geomagnetic field in West Africa: First absolute intensity results from Burkina Faso[END_REF][START_REF] Kapper | Novel insights on the geomagnetic field in West Africa from a new reference curve (0-2000 AD)[END_REF]. The greater abundance of intensity data can be related to the availability of material to study (see section 2.5).

Temporal distribution of archaeomagnetic data

There is a large variability in the temporal distribution of data in GEOMAGIA50.v3.4 over the past 10,000 years (Fig. 5). Both the number of archaeomagnetic directions and intensity in J o u r n a l P r e -p r o o f [START_REF] Goff | An updated archeomagnetic directional variation curve for France over the past two millennia, following 25 years of additional data acquisition[END_REF].

There is a peak in archaeointensity age entries during the first millennium BCE, where there has been concerted efforts to characterize the Levantine intensity spike (see section 4.1.2).

There are notable minor peaks in the number of intensity entries during the Neolithic, with notable studies from the Neolithic and Bronze age from China (207 entries) (see [START_REF] Cai | High-fidelity archeointensity results for the Late Neolithic Period from Central China[END_REF], Iraq (179 entries) [START_REF] Sakai | Variation of geomagnetic field deduced from analysis of remanent magnetization of archaeological objects[END_REF]Nachasova andBurakov, 1995, 1998;Yutsis-Akimova et al., 2018a,b) and the rest of the Middle East (148 entries) (e.g., [START_REF] Kawai | Archaeomagnetism in Iran[END_REF][START_REF] Gallet | Archaeomagnetism at Ebla (Tell Mardikh, Syria). new data on geomagnetic field intensity variations in the Near East during the Bronze Age[END_REF][START_REF] Stillinger | Refining the archaeomagnetic dating curve for the Near East: new intensity data from Bronze Age ceramics at Tell Mozan, Syria[END_REF][START_REF] Shaar | Large geomagnetic field anomalies revealed in Bronze to Iron Age archeomagnetic data from Tel Megiddo and Tel Hazor, Israel[END_REF][START_REF] Gallet | Analysis of geomagnetic field intensity variations in Mesopotamia during the third millennium BC with archeological implications[END_REF], Bulgaria (136 entries) (e.g., [START_REF] Kovacheva | Archaeomagnetic database from Bulgaria: the last 8000 years[END_REF]Kovacheva et al., 2009a[START_REF] Kovacheva | Extended and revised archaeomagnetic database and secular variation curves from Bulgaria for the last eight millennia[END_REF][START_REF] Kostadinova-Avramova | Archaeomagnetic knowledge of the Neolithic in Bulgaria with emphasis on intensity changes[END_REF] and Spain (79 entries) [START_REF] Nachasova | Geomagnetic field intensity variation in Spain[END_REF][START_REF] Nachasova | Archaeomagnetic study of ceramics from the Neolithic Los Castillejos multilayer monument (Montefrio, Spain)[END_REF][START_REF] Carrancho | First directional European palaeosecular variation curve for the Neolithic based on archaeomagnetic data[END_REF]. We note that there are very few BCE data from the Southern Hemisphere. There are only a few data per century back to 6000 BCE (Fig. 5d). In contrast the Northern Hemisphere (Fig. 5c)

has 10 times or more data per century.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Overview of archaeological materials

A wide range of archaeological materials and structures can be used to obtain directional and intensity information (Fig. 6). Almost all data from archaeological material ( 99%) were recovered from baked clays that acquired a TRM roughly parallel and proportional to the ambient geomagnetic field at the time of their firing. A few other archaeological materials can carry a remanent magnetization acquired through different processes. In mural paintings, e.g., frescoes, red pigments with hematite can acquire a so-called pictorial remanent magnetization when paint is sufficiently liquid to enable hematite grains to orientate parallel to the geomagnetic field (e.g., [START_REF] Chiari | Pictorial remanent magnetization as an indicator of secular variation of the Earth"s magnetic field[END_REF][START_REF] Zanella | Archaeomagnetic results from mural paintings and pyroclastic rocks in Pompeii and Herculaneum[END_REF]. Through a related process, lime-plasters (e.g., [START_REF] Hueda-Tanabe | Archaeomagnetic studies in central Mexico-dating of Mesoamerican lime-plasters[END_REF] and unburnt adobe bricks (e.g., [START_REF] Games | The magnitude of the paleomagnetic field: a new non-thermal, non detrital method using sun-dried bricks[END_REF] can also acquire a remanent magnetization, when the plaster or the clay is mixed with water. These materials are promising, even though experimental uncertainties are generally higher than for baked clays.

Fifty types of materials and structures are listed in the current version of GEOMAGIA50; however, there are some that have been sampled more frequently than others. In Fig. 6 we list the 8 most commonly used. In some cases (13% of entries) the type of material that was used is not given in the database. There are clear differences in the materials used for directional and intensity studies. For directional analysis, in-place oriented structures are necessary. Therefore kilns, ovens and hearths, bricks, and burnt structures are frequently used. For intensity the materials do not need to be in-situ, which allows a more diverse array of materials to be pooled from. Pottery and ceramics, owing to their abundance and ease of sampling are therefore the most common for intensity analysis. Over recent years copper slags have been used owing to their magnetically appropriate characteristics for intensity experiments [START_REF] Shaar | Testing the accuracy of absolute intensity estimates of the ancient geomagnetic field using copper slag material[END_REF].

Materials suitable for intensity are often easier to access, because the material has already been sampled and the collections they are from are well-studied. Sampling of these objects is also less invasive. For directional studies, it is necessary to be reactive to an archaeological excavation.

In-situ structures are uncovered and maybe destroyed when working on rescue excavations.

Sometimes kiln-type structures are preserved because of an obvious archaeological interest; however, sampling is invasive and possibly incompatible with heritage conservation. These issues may partially explain why the proportion of direction and intensity studies varies in different countries (section 2.3).

J o u r n a l P r e -p r o o f

Experimental considerations and data quality

In this section we outline the methods that have been used to obtain archaeomagnetic data and date archaeological materials. We discuss how experimental methods and practices affect the accuracy and precision of archaeomagnetic and chronological data and address how uncertainties are represented in the database. For intensity experiments we cover alteration during heating, the influence of multi-domain grains, remanence anisotropy and the effect of cooling rate.

Dating methods applied to archaeological materials are varied and we group them into two categories: those that directly or indirectly date a material. We discuss the nuances of these methods when applied to archaeological materials, how they can be combined to create a site chronology, and their age uncertainties.

Archaeomagnetic measurements

Directions

Three approaches have commonly been used to recover directional data from archaeological materials. The first two involve stepwise removal (demagnetization) of a TRM by either heating to increasing temperatures (thermal demagnetization) or by increasing the alternating current of a field coil (alternating field (AF) demagnetization). For some entries in the database both approaches have been used in conjunction. An alternative approach is to use viscosity cleaning. Developed by Émile Thellier (see, [START_REF] Thellier | Sur la direction du champ magnetique terrestre, en France, durant les deux derniers millenaires[END_REF], viscosity cleaning has proven to be as effective as a complete demagnetization in isolating directions, when a sample records a single TRM component. See Le [START_REF] Goff | An updated archeomagnetic directional variation curve for France over the past two millennia, following 25 years of additional data acquisition[END_REF] for an overview of this two-step method. Unfortunately, 56% of entries in the database do not report the demagnetization method used. Of entries that do list a demagnetization method, alternating field (AF) demagnetization is J o u r n a l P r e -p r o o f

Journal Pre-proof the most commonly used (33%), followed by viscosity cleaning (28%) (largely from entries from France, Le [START_REF] Goff | An updated archeomagnetic directional variation curve for France over the past two millennia, following 25 years of additional data acquisition[END_REF], a mixture of AF and thermal demagnetization (23%), and solely thermal demagnetization (16%).

There are various factors that are likely to interfere with the accurate recovery of past field directions. First is the precision of the sampling and sample orientation, which is critical in archaeomagnetism where one tries to recover small directional variations. Conservation of structures and mechanical problems, such as the inward or outward sagging of the walls or a slight tilting of the kiln sole, can influence the precision and reliability of the archaeomagnetic direction.

The direction recorded by a structure can be further perturbed by magnetic refraction, whereby the magnetization of a structure can distort the magnetic field recorded, in particular when the magnetization is strong (e.g., [START_REF] Aitken | Archaeomagnetism: evidence for magnetic refraction in kiln structures[END_REF][START_REF] Hus | On the suitability of refractory bricks from a mediaeval brass melting and working site near Dinant (Belgium) as geomagnetic field recorders[END_REF]. This can also result from differential cooling as, for instance, may occur in large structures [START_REF] Lanos | The effects of demagnetizing fields on thermoremanent magnetization acquired by parallel-sided baked clay blocks[END_REF]. Understanding magnetic refraction requires dense sampling across all parts of a structure. Too much localized sampling can lead to a precise but biased mean direction.

Another factor that may bias remanence directions is the anisotropy of TRM. For bricks or tiles used to mason all or parts of a kiln, this effect results in a recorded direction that may deviate from the ancient field. Taking this effect into account requires the determination of an anisotropy tensor (see details on the correction for anisotropy effects in section 3.1.5). For baked clay ovens or hearths, the degree of anisotropy is usually considered to be weak, and does not impact the remanence direction, e.g., Kovacheva et al. (2009b) and Le [START_REF] Goff | An updated archeomagnetic directional variation curve for France over the past two millennia, following 25 years of additional data acquisition[END_REF]. However, it should be noted that a significant shallowing of inclinations of up to 13 was recently documented for thin oven soles [START_REF] Palencia-Ortas | New archaeomagnetic directions from Portugal and evolution of the geomagnetic field in Iberia from Late Bronze Age to Roman Times[END_REF], 2021). We further note that the GEOMAGIA50 database does not yet make it possible to assess whether or not the anisotropy effect has been evaluated and taken into account in the directional studies.

On the whole the precision of directional data within GEOMAGIA50.v3.4 is variable, but is in general of statistically good quality, with 80% of entries having  values are notably high and some k are very low. The precision of directional data can be difficult to quantify for some J o u r n a l P r e -p r o o f

Journal Pre-proof entries as 95  is not specified for 6% of directional entries and k is given for only 50%. We also note that the method of calculating 95  is not always noted in publications. There are two forms of the 95  equation; the original equation in Fisher (1953) and an approximation for a large number of samples (see, e.g., [START_REF] Butler | Paleomagnetism: Magnetic domains to geologic terrains[END_REF]. These can result in different values of 95

 if the number of samples is less than approximately 10.

Less than 2% of entries are based on the successful analysis of only one or two samples and have no associated 95  or k . When the number of successfully measured samples is at least equal to 3, k is greater than 100 for 80% of the entries reporting k (40% of the all directional results). Any study wishing to use directional data should assess the uncertainty that they are comfortable in incorporating into the analysis.  or 93% of all directional entries. Only k values < 4000 are plotted, totalling 2769 values (91% of all entries with k; 47% of all directional entries). Whether 95  is calculated using the full equation of Fisher (1953) or an approximation (see, [START_REF] Butler | Paleomagnetism: Magnetic domains to geologic terrains[END_REF] is not noted in the database as it is commonly not stated. (c) Uncertainty on archaeointensity estimates expressed as a percentage of the archaeointensity value (bin size = 1%). Note that the uncertainties plotted here are those given by the author and result from different approaches to calculating uncertainty.

Archaeointensity determinations

The linearity at low fields ( < 150  T) between geomagnetic field strength and the intensity of a TRM acquired on cooling in this field is the physical basis for intensity estimates. A detailed description of the protocols is beyond the scope of this article as there are numerous approaches and derivatives that can be used (e.g., [START_REF] Dunlop | Physical basis of the Thellier-Thellier and related paleointensity methods[END_REF][START_REF] Tauxe | 13 -Paleointensities[END_REF][START_REF] Tauxe | Essentials of Paleomagnetism, 5th Web Edition[END_REF], but we give an overview of those used for archaeological entries in GEOMAGIA50 and review the different experimental strategies used to detect and/or possibly

J o u r n a l P r e -p r o o f

Journal Pre-proof mitigate various effects that influence the intensity measurements.

GEOMAGIA50.v3.4 lists 25 palaeointensity methods and variants; however, these can be primarily classed into five main types, as listed in Fig. 8a. Thellier-type approaches that use double heating steps to impart a laboratory induced TRM as proposed by [START_REF] Thellier | Sur l"intensité du champ magnétique terrestre dans le passé historique et géologique[END_REF] make up 87% of all intensity entries in the database. Of the Thellier-type approaches, the original [START_REF] Thellier | Sur l"intensité du champ magnétique terrestre dans le passé historique et géologique[END_REF] method has been used more than any other method, followed by the Coe-Thellier approach [START_REF] Coe | Paleo-intensities of the Earth"s magnetic field determined from Tertiary and Quaternary rocks[END_REF]. The IZZI protocol [START_REF] Yu | Toward an optimal geomagnetic field intensity determination technique[END_REF] has increasingly been used over recent years as the revised order of the in-field and zero-field steps during the experiment aids in the identification of non-ideal (MD) grains that can bias intensity estimates (section 3.1.4). It currently makes up 11% of Thellier-type entries, but we anticipate it will be used increasingly over coming years. Other Thellier-type variants, such as that of [START_REF] Aitken | Determination of the intensity of the Earth"s magnetic field during archaeological times: Reliability of the Thellier technique[END_REF]), MT4 of Leonhardt et al. (2004), and the two specimen approach of [START_REF] Domen | A single heating method of paleomagnetic field intensity determination applied to old roof tiles and rocks[END_REF], make up only a minor contribution to the database. method [START_REF] Thellier | Sur l"intensité du champ magnétique terrestre dans le passé historique et géologique[END_REF], Coe-Thellier (Coe, 1967), Aitken [START_REF] Aitken | Determination of the intensity of the Earth"s magnetic field during archaeological times: Reliability of the Thellier technique[END_REF], IZZI [START_REF] Yu | Toward an optimal geomagnetic field intensity determination technique[END_REF], MT4 method [START_REF] Leonhardt | Multidomain behavior during Thellier paleointensity experiments: a phenomenological model[END_REF], the two specimen approach of [START_REF] Domen | A single heating method of paleomagnetic field intensity determination applied to old roof tiles and rocks[END_REF] and other non-specific Thellier-based methods. Note that for three entries two Thellier-type methods were used for the mean intensity given in the entry, therefore the individual mentions of Thellier-type methods totals 4610. Shaw methods include the original procedure [START_REF] Shaw | A new method of determining the magnitude of the palaeomagnetic field. Application to five historic lavas and five archaeological samples[END_REF] and modified versions by [START_REF] Kono | Reliability of palaeointensity methods using alternating field demagnetization and anhysteretic remanence[END_REF], [START_REF] Rolph | A new method of paleofield magnitude correction for thermally altered samples and its application to Lower Carboniferous lavas[END_REF], [START_REF] Shaw | Archaeointensity variations for the past 7,500 years evaluated from ancient Chinese ceramics[END_REF], [START_REF] Tsunakawa | The Shaw method of palaeointensity determinations and its application to recent volcanic rocks[END_REF], [START_REF] Yamamoto | Palaeointensity study of the Hawaiian 1960 lava: implications for possible causes of erroneously high intensities[END_REF]. Triaxe method is that of Le [START_REF] Goff | A new three-axis vibrating sample magnetometer for continuous high-temperature magnetization measurements: applications to paleo-and archeo-intensity determinations[END_REF].

Microwave methods are based on versions of the Thellier-type approaches listed above (see, e.g., Hill andShaw, 1999, 2007). The multispecimen entries include both [START_REF] Dekkers | Reliable absolute palaeointensities independent of magnetic domain state[END_REF] and [START_REF] Fabian | Multiple-specimen absolute paleointensity determination: An optimal protocol including pTRM normalization, domain-state correction, and alteration test[END_REF] approaches. Other methods are the approach of [START_REF] Walton | Archaeomagnetic intensity measurements using a SQUID magnetometer[END_REF] and the calibrated pseudo-Thellier method (de [START_REF] De Groot | Rapid regional perturbations to the recent global geomagnetic decay revealed by a new Hawaiian record[END_REF].

The remaining 13% of palaeointensity estimates were determined by variants of the Shaw

J o u r n a l P r e -p r o o f

technique proposed by [START_REF] Wilson | Paleomagnetism in Northen Ireland, Part I, The thermal demagnetization of natural magnetic moments in rocks[END_REF]), microwave variants of Thellier-type protocols [START_REF] Shaw | A comparison of archaeointensity results from Chinese ceramics using microwave and convential Thellier"s and Shaw"s methods[END_REF]Hill andShaw, 1999, 2007;[START_REF] Stark | Establishing a first archaeointensity record for the SW Pacific[END_REF]) (2%) and the two variants of the multispecimen parallel differential partial TRM (pTRM) method [START_REF] Dekkers | Reliable absolute palaeointensities independent of magnetic domain state[END_REF][START_REF] Fabian | Multiple-specimen absolute paleointensity determination: An optimal protocol including pTRM normalization, domain-state correction, and alteration test[END_REF]) (1%). The calibrated pseudo-Thellier method (de [START_REF] De Groot | Rapid regional perturbations to the recent global geomagnetic decay revealed by a new Hawaiian record[END_REF] and the approach of [START_REF] Walton | Archaeomagnetic intensity measurements using a SQUID magnetometer[END_REF] contribute less than 1% of entries.

Checking and/or correcting for thermal alteration

As noted above, most archaeointensity data have been obtained using protocols derived from the original [START_REF] Thellier | Sur l"intensité du champ magnétique terrestre dans le passé historique et géologique[END_REF] method. Its principle is based on the stepwise thermal demagnetization of the natural remanent magnetization (NRM, assumed to be a TRM) and its progressive replacement by a new TRM acquired in a laboratory field whose direction and intensity are controlled. The ratio between the remaining NRM and the partial TRM acquired after each heating/cooling step, with data usually displayed on an Arai-Nagata diagram [START_REF] Nagata | Secular varaition of the geomagnetic force during the last 5000 years[END_REF], allows an estimation of the past geomagnetic field intensity. The comparison of NRM lost to TRM gained requires the magnetic mineralogy of the specimen to remain unchanged during the thermal treatment. In order to assess alteration, Thellier suggested as early as 1946 a partial-TRM check (a pTRM check) [START_REF] Thellier | Sur l"intensité du champ magnétique terrestre, en France, à l"époque gallo-romaine[END_REF]. During the stepwise heating-cooling cycle, additional pTRM acquisition steps are added. After a number of heating steps, a lower temperature step is repeated and the pTRMs compared. This is done multiple times throughout the experiment, e.g., after every three heating steps, the first step of the three will be repeated. This alteration test is now common and always required for modern intensity studies using the Thellier method and derivatives (i.e. [START_REF] Coe | Paleo-intensities of the Earth"s magnetic field determined from Tertiary and Quaternary rocks[END_REF][START_REF] Aitken | Determination of the intensity of the Earth"s magnetic field during archaeological times: Reliability of the Thellier technique[END_REF][START_REF] Yu | Toward an optimal geomagnetic field intensity determination technique[END_REF]. It is important to underline that different approaches have been used to calculate the degree of alteration at each pTRM check (commonly expressed as a percentage) and the associated cut-off values to accept or reject a check or an intensity determination. 44% of Thellier-type intensity entries are accompanied by a pTRM check; however, this number hides the variability in the statistical cutoffs used (see [START_REF] Genevey | ArcheoInt: An upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment[END_REF][START_REF] Paterson | On improving the selection of Thellier-type paleointensity data[END_REF].

Monitoring magnetic susceptibility during heating has been used to check for the stability of the magnetic mineralogy; however, it must be noted that slight changes in susceptibility may not relate to changes in remanence carrying minerals or the formation of new remanence carriers

J o u r n a l P r e -p r o o f

Journal Pre-proof (rather changes in the susceptibility of magnetic minerals that do not have the capacity to hold or acquire a remanence). This approach was for example used for the part of the Bulgarian data set acquired in the 70s and 80s [START_REF] Kovacheva | Extended and revised archaeomagnetic database and secular variation curves from Bulgaria for the last eight millennia[END_REF]. Susceptibility monitoring was only used for 1.5% of intensity entries in the database.

Instead of rejecting samples for which alteration is judged too strong, another possibility is to correct for this effect. This was proposed by [START_REF] Burakov | Correcting for chemical change during heating in archeomagnetic determinations of the ancient geomagnetic field intensity[END_REF], with a protocol that additionally takes into account anisotropy of TRM. Several sets of data were acquired using this protocol (26 studies spanning 1986 to the present day). This protocol which has not been used in other laboratories is viewed with caution.

For 30% of database entries listing the use of a Thellier-type protocol, no alteration test was performed to check or correct for alteration: the linearity of the data points in the Arai-Nagata diagram over a large proportion of the unblocking temperatures was judged sufficient to testify of the absence of this effect. This concerns mainly studies published before the 1990s.

In the Triaxe method (Le [START_REF] Goff | A new three-axis vibrating sample magnetometer for continuous high-temperature magnetization measurements: applications to paleo-and archeo-intensity determinations[END_REF] measurements are made continuously in temperature, through successive series of heating and cooling, in zero field or laboratory field. The stability of the magnetic mineralogy is assessed by checking the stability of the ratio between the demagnetized NRM fraction and the acquired TRM fraction at each increasing temperature step.

This approach corresponds in a similar way to testing the linearity in an Arai-Nagata diagram, but the steps are spaced only 5 C apart: the data are therefore numerous (e.g., 60 data for a 300 C temperature interval) and the linearity is thus finely checked and also assessed through specific linearity tests (Le [START_REF] Goff | A new three-axis vibrating sample magnetometer for continuous high-temperature magnetization measurements: applications to paleo-and archeo-intensity determinations[END_REF].

To mitigate the risk of magnetic alteration, alternative methods have been developed. From the oldest to the most recent: the Shaw technique and derivatives [START_REF] Shaw | A new method of determining the magnitude of the palaeomagnetic field. Application to five historic lavas and five archaeological samples[END_REF][START_REF] Tsunakawa | The Shaw method of palaeointensity determinations and its application to recent volcanic rocks[END_REF][START_REF] Yamamoto | Palaeointensity study of the Hawaiian 1960 lava: implications for possible causes of erroneously high intensities[END_REF], the microwave technique (e.g., [START_REF] Walton | Application of ferrimagnetic resonance heating to palaeointensity determinations[END_REF][START_REF] Hill | Palaeointensity results for historic lavas from Mt Etna using microwave demagnetization/remagnetization in a modified Thellier-type experiment[END_REF] and the multispecimen protocol and adaptations [START_REF] Dekkers | Reliable absolute palaeointensities independent of magnetic domain state[END_REF][START_REF] Fabian | Multiple-specimen absolute paleointensity determination: An optimal protocol including pTRM normalization, domain-state correction, and alteration test[END_REF].

Most data obtained with the Shaw technique were acquired between 1975 and 1995 (e.g., [START_REF] Liritzis | Palaeointensity and thermoluminescence measurements on Cretan kilns from 1300 to 2000 BC[END_REF][START_REF] Shaw | Archaeointensity variations for the past 7,500 years evaluated from ancient Chinese ceramics[END_REF], but the method has seen a revival in recent years [START_REF] Kitahara | Archeointensity estimates of a tenth-century kiln: first application of the Tsunakawa-Shaw paleointensity method to archeological relics[END_REF][START_REF] Kitahara | Construction of new archaeointensity reference curve for East Asia from 2000 CE to 1100 CE[END_REF] in the form of the modified Tsunakawa-Shaw approach [START_REF] Yamamoto | Palaeointensity study of the Hawaiian 1960 lava: implications for possible causes of erroneously high intensities[END_REF]. The Shaw method involves only one heating in which the sample is heated above its Curie temperature allowing the acquisition of a full TRM. Prior to heating the NRM is stepwise

J o u r n a l P r e -p r o o f

Journal Pre-proof demagnetized using increasing alternating field (AF) steps. After heating the sample is again demagnetized using the AF steps as for the NRM. The linear relationship of the demagnetized NRM to TRM is then used to calculate an estimate of palaeointensity. Alteration is assessed through a comparison of coercivity spectra. Changes in an AF demagnetized anhysteretic magnetization (ARM) given before and after heating are compared. Later modifications to the method incorporated corrections to take into account alteration to the pre-and post-heating ARM spectra [START_REF] Kono | Reliability of palaeointensity methods using alternating field demagnetization and anhysteretic remanence[END_REF][START_REF] Rolph | A new method of paleofield magnitude correction for thermally altered samples and its application to Lower Carboniferous lavas[END_REF].

The microwave method follows the protocols of Thellier and Thellier and modified variants, e.g., the perpendicular single heating method [START_REF] Hill | The use of the "Kono perpendicular applied field method" in microwave palaeointensity experiments[END_REF], but thermal demagnetization is replaced by microwave demagnetization. The rationale is that microwave power should limit the rise in temperature of the sample matrix and reduce the possibility of alteration. However, some conversion to thermal energy to heat the matrix is likely and pTRM-checks test are now integrated in the microwave technique. Recent studies have also included checks for evaluating the cooling rate effect (e.g., [START_REF] Poletti | The cooling-rate effect on microwave archeointensity estimates[END_REF][START_REF] Ertepinar | Extreme geomagnetic field variability indicated by Eastern Mediterranean full-vector archaeomagnetic records[END_REF].

The multispecimen parallel differential partial pTRM method [START_REF] Dekkers | Reliable absolute palaeointensities independent of magnetic domain state[END_REF] started life as essentially a very simple method. Multiple specimens from a site were heated at the same temperature (below the temperature of alteration, but high enough for an appreciable decrease in NRM), but with a different field for each specimen aligned with the specimens NRM.

However, shortcomings in the method were evident and the method was expanded upon by [START_REF] Fabian | Multiple-specimen absolute paleointensity determination: An optimal protocol including pTRM normalization, domain-state correction, and alteration test[END_REF]. It was elaborated upon to correct for differences in the fraction of the pTRM imparted in each specimen, a specimen"s domain state, and included a step to monitor alteration.

Checking or correcting for the presence of multi-domain grains

Another possible factor for the failure of intensity determinations is linked to the presence of MD grains for which the laws of reciprocity and additivity of the partial TRMs are not obeyed [START_REF] Néel | Théorie du traînage magnétique des ferromagnétiques en grains fins avec applications aux terre cuites[END_REF]. Although the influence of MD grains on volcanic palaeointensity estimates has been investigated in detail, it has received less attention in archaeomagnetic studies. This is primarily a result of the different grain size distributions found in archaeomagnetic materials compared with volcanic rocks: archaeomagnetic materials are commonly dominated by pseudo-single domain grains, which are not effected by pTRM tails, whereas volcanic rocks frequently contain a MD

J o u r n a l P r e -p r o o f

Journal Pre-proof fraction where pTRM tails are significant (where a pTRM-tail results from a non-reciprocity between the blocking and unblocking temperatures). The influence of MD grains can be recognized on Arai-Nagata diagrams as a concave-up curve, whose misinterpretation can lead to underestimates or overestimates of intensity depending on which portion of the curve was used to calculate palaeointensity (e.g., [START_REF] Levi | The effect of magnetite particle size on paleointensity determinations of the geomagnetic field[END_REF][START_REF] Dunlop | Physical basis of the Thellier-Thellier and related paleointensity methods[END_REF]. The linearity of the data in the Arai diagram was often considered as a sufficient criterion to, if not exclude, at least consider that the proportion of MD grains is too small to critically affect the intensity determination. The presence of MD grains is now more directly investigated with either rock magnetic measurements, such as hysteresis curves, backfield curves and first order reversal curves (see, e.g., [START_REF] Day | Hysteresis properties of titanomagnetites: grain size and compositional dependence[END_REF][START_REF] Dunlop | Theory and application of the Day plot ( / rs s MM versus / cr c HH ) 1. Theoretical curves and tests using titanomagnetite data[END_REF][START_REF] Roberts | Domain state diagnosis in rock magnetism: Evaluation of potential alternatives to the Day diagram[END_REF], or through additional tests implemented during Thellier-type methods and microwave protocols, such as pTRM-tail checks (aiming at testing the independence of pTRM; [START_REF] Riisager | Detecting multidomain magnetic grains in Thellier palaeointensity experiments[END_REF] and additivity checks [START_REF] Krása | Experimental procedure to detect multidomain remanence during Thellier-Thellier experiments[END_REF]. Only 5% of intensity entries in the database list an MD check.

The IZZI protocol [START_REF] Yu | Toward an optimal geomagnetic field intensity determination technique[END_REF], a variant on the Thellier method, was designed to accentuate the influence of MD tails, evident by pronounced zig-zagging in the Arai-Nagata plot.

However, this method is sensitive to the direction of the laboratory field relative to the orientation of the NRM leading to over-or under-estimation of the pTRM-tail and with the field aligned with the direction of the NRM, MD tails can be suppressed.

In comparison to other protocols, the MSP-DSC method of [START_REF] Fabian | Multiple-specimen absolute paleointensity determination: An optimal protocol including pTRM normalization, domain-state correction, and alteration test[END_REF] has the advantage to (partially) correct intensities for domain state effect. The Triaxe protocol (Le [START_REF] Goff | A new three-axis vibrating sample magnetometer for continuous high-temperature magnetization measurements: applications to paleo-and archeo-intensity determinations[END_REF] mitigates the spurious effect of large grains because the laboratory TRM is almost a full one, mimicking the acquisition of the original TRM. The Shaw derivative of [START_REF] Yamamoto | Palaeointensity study of the Hawaiian 1960 lava: implications for possible causes of erroneously high intensities[END_REF] aims to remove all MD contributions by incorporating a low-temperature demagnetization step after each remanence acquisition and prior to AF demagnetization.

TRM anisotropy

An important parameter that may affect intensity determinations when analyzing baked clay artefacts is anisotropy of TRM (already touched upon in section 3.1.1). This anisotropy arises from the stretching of clay during the process of shaping an object, resulting in a preferential alignment of magnetic grains in the clay matrix (e.g., [START_REF] Rogers | Magnetic anisotropy in ancient pottery[END_REF][START_REF] Aitken | Archaeomagnetic determination of J o u r n a l P r e -p r o o f Journal Pre-proof the past geomagnetic intensity using ancient ceramics: Allowance for anisotropy[END_REF]. This effect may be particularly intense for pottery fragments and thin tiles and to a lesser extent to thick J o u r n a l P r e -p r o o f

Journal Pre-proof bricks, with biases up to several dozens of micro Tesla (e.g., [START_REF] Genevey | ArcheoInt: An upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment[END_REF][START_REF] Hervé | Fast geomagnetic field intensity variations between 1400 and 400 BCE: New archaeointensity data from Germany[END_REF][START_REF] Gómez-Paccard | New archeointensity data from NW Argentina (1300-1500 CE)[END_REF]. Conversely, it has been observed that this effect is generally less critical when analysing fragments made of clay, which are coarsely assembled, as they are usually taken from in situ structures (e.g., Kovacheva et al., 2009b). applied parallel to specimen natural remanent magnetization (NRM) direction during the palaeointensity method; other corrections are generally approaches that were insufficiently defined in a publication). See more details in section 3.1.5 and Section 3.1.6, For 38% of archaeointensity entries in the database remanence anisotropy was not investigated (Fig. 9a). In some cases data were obtained from less anisotropic materials and no measure of anisotropy was pursued. In a small number of entries where anisotropy was estimated, a correction was not necessary. This is most likely as the anisotropy was not considered to be significant. Different approaches have been proposed to evaluate remanence anisotropy.

Determination of a TRM anisotropy tensor for each analysed sample allows to evaluate the importance of this effect and to accurately correct the raw intensity determinations [START_REF] Veitch | An investigation of the intensity of the geomagnetic field J o u r n a l P r e -p r o o f Journal Pre-proof during Roman times using magnetically anisotropic bricks and tiles[END_REF][START_REF] Selkin | Long-term varations in palaeointensity[END_REF]. This approach was used for 30% of entries considering anisotropy in the database. The drawback of this approach is the time-consuming multiple heating steps (usually six), which increases the risk of mineralogical alteration. Aligning the laboratory field direction with the original NRM (25% of entries considering anisotropy) is an adequate alternative, as long as the degree of anisotropy is not too strong to bias significantly the direction (e.g., [START_REF] Aitken | Archaeomagnetic determination of J o u r n a l P r e -p r o o f Journal Pre-proof the past geomagnetic intensity using ancient ceramics: Allowance for anisotropy[END_REF]. Ideally, the laboratory field direction should be aligned with the ancient ambient field. This is achieved with the Triaxe protocol and MSP protocols where the direction of the laboratory field is adjusted so a TRM is imparted parallel to the primary TRM (see, Le [START_REF] Goff | A new three-axis vibrating sample magnetometer for continuous high-temperature magnetization measurements: applications to paleo-and archeo-intensity determinations[END_REF]. To minimize the effect of TRM anisotropy, [START_REF] Morales | Magnetic properties and archeointensity determination on Pre-Columbian pottery from Chiapas, Mesoamerica[END_REF] proposed to average the intensity values obtained for 6 specimens from the same fragment: here the specimens are oriented in such a way that the TRM is acquired in 6 orthogonal directions relative to a fixed J o u r n a l P r e -p r o o f

Journal Pre-proof arbitrary orientation. However, [START_REF] Poletti | Archeomagnetism of Jesuit Missions in South Brazil (1657-1706 AD) and assessment of the South American database[END_REF] and Hervé et al. (2019b) demonstrated that this approach results in larger standard deviations and possibly significant inaccuracies as high as 10-15  T.

As an alternative to the full determination of the TRM anisotropy tensor, it has been suggested to use other tensors to evaluate and correct for anisotropy; namely tensors of magnetic susceptibility (AMS; 14% of anisotropy assessed entries), anhysteretic remanent magnetization (ARM; 8%) or isothermal remanent magnetization (IRM; < 1% ). These substitutes are often quicker and easier to implement and avoid the six additional heatings during the thermal protocol.

However, the respective ellipsoids significantly differ in their shape and anisotropy degree from TRM ellipsoids (e.g., [START_REF] Chauvin | Paleointensity of the geomagnetic field recovered on archaeomagnetic sites from France[END_REF]. AMS can underestimate TRM anisotropy by several dozens of percent [START_REF] Gómez-Paccard | New archeointensity data from NW Argentina (1300-1500 CE)[END_REF]. In 22% of entries other types of anisotropy corrections have been applied, but either a method was not listed in the database or the method was not described in the publication.

Cooling rate effect

Another possible biasing factor for intensity determinations is the cooling rate dependence of TRM intensity [START_REF] Fox | Cooling-rate dependency of thermoremanent magnetisation[END_REF]. Ideally, to avoid such systematic bias, the cooling duration used for the acquisition of the laboratory TRM should be chosen to be identical to the original one when the primary TRM was recorded by the archaeological object. This is rarely possible as the original cooling time is usually long, ranging typically from half a day to a few days (with the notable exception of the slags, [START_REF] Shaar | Testing the accuracy of absolute intensity estimates of the ancient geomagnetic field using copper slag material[END_REF], while the laboratory cooling time is faster, generally from 0.5 up to 2 hours, depending of the type of oven and the size of the specimens.

For Thellier-Thellier data, the cooling rate effect can be evaluated through a comparison of the TRM acquired with a rapid cooling time (the one used routinely during the experiment) and a slow cooling time chosen to be close to the original one (e.g., [START_REF] Chauvin | Paleointensity of the geomagnetic field recovered on archaeomagnetic sites from France[END_REF][START_REF] Leonhardt | Cooling rate correction of paleointensity determination for volcanic glasses by relaxation geospeedometry[END_REF][START_REF] Poletti | The cooling-rate effect on microwave archeointensity estimates[END_REF]. This is performed for 28% of the intensity entries in the database and comprises 80% of entries that used a cooling rate correction (Fig. 9b). Precisely evaluating the duration of the past cooling is the main difficulty of the correction protocol. Experimental archaeology has provided constraints on this issue (e.g., [START_REF] Morales | Are ceramics and bricks reliable absolute geomagnetic intensity carriers?[END_REF][START_REF] Calvo-Rathert | Reproducibility of archaeointensity determinations with a multimethod approach on archaeological material reproductions[END_REF][START_REF] Genevey | New archeointensity data from French Early Medieval pottery production (6-10th century AD). Tracing 1500 years of geomagnetic field intensity variations in Western Europe[END_REF][START_REF] Schnepp | Validity of archaeomagnetic field recording: an experimental pottery kiln at Coppengrave, Germany[END_REF]Hervé et al., 2019a;[START_REF] Jones | Archeointensity of the Four Corners Region of the American Southwest[END_REF].

J o u r n a l P r e -p r o o f

Journal Pre-proof Archaeological information concerning, e.g., the estimated size of kilns, their morphology, and the type of firing (open or closed), may also help to assess the original cooling conditions. Another approach is to measure the cooling rate effect on TRM acquisition with increasingly slow cooling duration (therefore exploring different conditions of cooling) and to infer from this the error that would be made by under or over estimating the original cooling rate [START_REF] Genevey | Eight thousand years of geomagnetic field intensity variations in the eastern Mediterranean[END_REF][START_REF] Hartmann | Archeointensity in Northeast Brazil over the past five centuries[END_REF].

A different possibility is to apply a fixed correction for all samples from the same archaeomagnetic site, usually 5% or 10%. This "educated guess" concerns 7% of intensity data or 20% of entries which applied a cooling rate correction (Fig. 9b). This approach is based on the assumption that all fragments of the same archaeological object show the same TRM intensity dependence on cooling rate. Experimental studies have, however, pointed out that this effect is variable from one sample to another and (as predicted by theory) that the TRM intensity increases following a logarithmic law as a function of the ratio between an increasingly slow cooling time and a fixed rapid one [START_REF] Genevey | ArcheoInt: An upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment[END_REF]Hervé et al., 2019a). To avoid applying an educated guess correction to all fragments, it has been suggested to estimate at least for part of the collection the cooling rate effect and to apply an average correction to the other fragments (Kostadinova-Avramova and Jordanova, 2019).

Another important question is at what temperature to estimate the effect of cooling rate. In particular, Hervé et al. (2019a) showed too high of a temperature could greatly overestimate this effect and therefore underestimate the intensity value. This appears to depend on the magnetic mineralogy of the material analysed (see also [START_REF] Kostadinova-Avramova | Study of cooling rate effect on baked clay materials and its importance for archaeointensity determinations[END_REF].

The cooling rate effect is a challenging parameter to estimate and many studies have not explored this question (over 70% of entries in the database). However, some of these data were obtained with a relatively slow cooling time as part of routine intensity experiments (for example for Bulgarian dataset; [START_REF] Kovacheva | Extended and revised archaeomagnetic database and secular variation curves from Bulgaria for the last eight millennia[END_REF]: the cooling rate effect is therefore expected to affect them less strongly.

Optimally, we would like to be able to dispense with the question of the cooling rate effect.

It has been observed experimentally that the Triaxe protocol accounts for cooling rate (Le [START_REF] Goff | A new three-axis vibrating sample magnetometer for continuous high-temperature magnetization measurements: applications to paleo-and archeo-intensity determinations[END_REF][START_REF] Genevey | Evidence for rapid geomagnetic field intensity variations in Western europe over the past 800 years from new French archeointensity data[END_REF][START_REF] Hartmann | Archeointensity in Northeast Brazil over the past five centuries[END_REF][START_REF] Hervé | Fast geomagnetic field intensity variations between 1400 and 400 BCE: New archaeointensity data from Germany[END_REF][START_REF] Salnaia | New archeointensity results on a baked-clay tile collection from the New Jerusalem Monastery (Moscow region, Russia)[END_REF].

The multispecimen parallel differential pTRM method also seems to be insensitive to cooling rate (e.g., [START_REF] Schnepp | Validity of archaeomagnetic field recording: an experimental pottery kiln at Coppengrave, Germany[END_REF][START_REF] Calvo-Rathert | Reproducibility of archaeointensity determinations with a multimethod approach on archaeological material reproductions[END_REF], possibly because in this technique all J o u r n a l P r e -p r o o f

pTRMs are acquired at medium temperatures. However, this question still needs to be further explored (Schnepp et al., 2020a).

Intensity uncertainties

On the whole archaeomagnetic data within GEOMAGIA50.v3.4 have reasonably well constrained uncertainties (Fig. 7c). The majority of estimates have an uncertainty of less than 10%

of the intensity estimate (60% of intensity entries that report uncertainties), i.e. a few  T on most measurements. Some intensity measurements, however, have high uncertainties, ranging up to 40  T. They require careful evaluation prior to their inclusion in reference curves or for field modelling.

A caveat to all intensity uncertainties in the database is that they have been calculated in a variety of ways. Uncertainties may be reported as standard deviations (to 1 or 2  ), standard errors or they could be weighted. The type of intensity uncertainty is not noted in the database.

Care must therefore be taken when using intensity uncertainties when constructing field models and reference curves and using this field as a selection criteria.

Dating methods

Dating and its accuracy and precision are key elements for any archaeomagnetic study. For archaeological artefacts, the dating methods used and listed in GEOMAGIA50 are based on archaeological or historical constraints, or chronometric methods involving mainly radioisotopic and physicochemical measurements (Fig. 10). See [START_REF] Aitken | Science-Based Dating in Archaeology[END_REF] for an overview of scientific dating methods. We briefly describe the most salient aspects of these methods and their caveats in section 3.2.1.

The archaeological approach remains the most common and concerns almost 60% of the database. Behind the term "archaeological dating" is often hidden the use of a relative chronology, optically stimulated luminescence (OSL; 9 entries) or rehydroxylation (5 entries) were used.

Direct dating of TRM

One of the main sources to directly date TRM are document archives. A well-known example is the eruption of Vesuvius first described by Pliny the Younger, which destroyed the city of Herculaneum and Pompeii in 79 CE [START_REF] Evans | An archaeointensity investigation of a kiln at Pompeii[END_REF][START_REF] Evans | Archaeomagnetic results from southern Italy and their bearing on geomagnetic secular variation[END_REF]. But more commonly, these archives are used, for example, to precisely date the edification of religious or civil buildings (e.g., [START_REF] Schnepp | Paleomagnetic directions and thermoluminescence dating from a bread oven-floor sequence in Lübeck (germany): A record of 450 years of geomagnetic secular variation[END_REF][START_REF] Osete | First archaeomagnetic field intensity data from Ethiopia, Africa (1615  12 AD)[END_REF][START_REF] Salnaia | New archeointensity results on a baked-clay tile collection from the New Jerusalem Monastery (Moscow region, Russia)[END_REF][START_REF] Genevey | Refining the high-fidelity archaeointensity curve for Western Europe over the past J o u r n a l P r e -p r o o f Journal Pre-proof millennium: analysis of Tuscan architectural bricks (Italy)[END_REF] or short periods of activity of ceramics workshops (e.g., [START_REF] Genevey | Evidence for rapid geomagnetic field intensity variations in Western europe over the past 800 years from new French archeointensity data[END_REF]. Other objects, such as some amphoras, can be precisely dated directly through the identification of stamps [START_REF] Ben-Yosef | Six centuries of geomagnetic intensity variations recorded by royal Judean stamped jar handles[END_REF].

Among the chronometric methods used in archaeology, thermoluminescence (TL) and optically stimulated luminescence (OSL) are directly associated to the TRM acquisition. A firing above 400 C is time-zero of the method as at above this temperature the electron traps in quartz or feldspar grains in baked clays are emptied [START_REF] Aitken | Thermoluminescence dating[END_REF]. From this moment, traps progressively fill again under irradiation from the surrounding environment (mainly related to 40 K, 238 U, 235 U and 232 Th radioactive isotopes). In spite of the advantage of dating the same instance as the TRM acquisition, luminescence methods constitute only 3% of entries in the database. However, this method has been used in recent studies (e.g., Gómez-Paccard et al., 2012a;[START_REF] Schnepp | Paleomagnetic directions and thermoluminescence dating from a bread oven-floor sequence in Lübeck (germany): A record of 450 years of geomagnetic secular variation[END_REF][START_REF] Kondopoulou | Archaeomagnetic study and thermoluminescence dating of Protobyzantine kilns (Megali Kypsa, North Greece)[END_REF][START_REF] Cai | New constraints on the variation of the geomagnetic field during the late Neolithic period: Archaeointensity results J o u r n a l P r e -p r o o f Journal Pre-proof from Sichuan, southwestern China[END_REF][START_REF] Aidona | Archaeomagnetism and luminescence on Late Byzantine kilns in Thessaloniki (N. Greece): Implications for geomagnetic field variations during the last two millennia[END_REF]. Accurate luminescence dating requires a careful reconstitution of the radioactive environment of the baked clay since the last firing. The resulting long measurement time limits the use of the techniques [START_REF] Roberts | Optical dating in archaeology: thirty years in retrospect and grand challenges for the future[END_REF].

Another caveat of luminescence methods are age uncertainties of  5-10% (1 ), corresponding J o u r n a l P r e -p r o o f

Journal Pre-proof to  100-200 years for a 0 CE baked clay for example. However, this can be reduced if multiple TL measurements are made. It is worth noting that OSL does not always provide a direct dating of the TRM acquisition because time-zero of this technique can also be the last exposure to sunlight, offering the possibility to date the deposit of sedimentary layers around the studied baked clays.

Another method to directly date baked clay artefacts was proposed by [START_REF] Wilson | Dating fired-clay ceramics using long-term power law rehydroxylation kinetics[END_REF].

It is based on the process of rehydroxylation (RHX) of fired-clay ceramics after production.

Similar to luminescence methods, the principle is to start from a zero point by heating a sample up to 500 C (dehydroxylation) and then monitor precisely the sample"s weight gain in known environmental conditions over several weeks (through rehydroxylation). This allows the kinetics of the rehydroxylation process to be determined. Although promising for archaeologists, and in turn for archaeomagnetists, the relationship between mass gain and time has proved more complex than initially thought, with kinetics that appear to depend on the nature and/or firing conditions of the ceramic (in addition to the environmental conditions), and the applicability of the RHX method appears clearly compromised [START_REF] Bowen | Rehydration/rehydroxylation kinetics of a reheated XIX-Century Davenport (Utah) ceramic[END_REF]Le Goff andGallet, 2014, 2015). So far it has only been applied to two archaeomagnetic studies, both on Spanish ceramics [START_REF] Nachasova | Variations in geomagnetic intensity and temperature in the second millennium B.C. in Spain[END_REF][START_REF] Burakov | Archaeomagnetic study and rehydroxylation dating of fired-clay ceramics[END_REF].

Indirect dating of TRM

As mentioned in the introduction to this section, the archaeological approach remains the most used method of indirect dating. Archaeological dating is however a very generic term that integrates many different elements. The first is stratigraphy, which is essential for building chronologies for ancient multi-layered sites in the Middle East [START_REF] Shaar | Geomagnetic field intensity: How high can it get? How fast can it change? Constraints from Iron Age copper slag[END_REF][START_REF] Gallet | Analysis of geomagnetic field intensity variations in Mesopotamia during the third millennium BC with archeological implications[END_REF] and Eastern Europe (e.g., [START_REF] Kostadinova-Avramova | Contribution of stratigraphic constraints of Bulgarian prehistoric multilevel tells and a comparison with archaeomagnetic observations[END_REF]. Elements such as coins, fragments of ceramics or metallic artefacts (e.g., swords and fibulae) are also key for dating, if the evolution of their typology is well known. Together, these elements make it possible to define a post quem and ante quem terminus terminus (lower and upper age limit respectively) for an archaeological level or artefact. It is also important to understand whether there has been any nixing of the layers in the stratigraphy, which can limit chronological control.The central question is to precisely understand how the object analysed for archaeomagnetism is reliably related to these chronological constraints. This question is far from trivial, e.g., for settlements occupied over a long period. For intensity determination, one way to overcome this issue is to work directly

J o u r n a l P r e -p r o o f

Journal Pre-proof on dated pottery fragments, i.e. those whose shape or decoration is recognized and can be linked to a known local/regional typo-chronology.

The relative chronology given by the stratigraphy is fixed to the calendar scale by historical events or chronometric methods. Their precision and reliability are mainly related to the state of the art of archaeological research in the region for a certain period. For example, in Western Europe, precise typo-chronologies are firm for the Roman period (0-500 CE), but are "floating"

for the Neolithic period (6000-2000 BCE). These typo-chronologies, and more generally archaeological dating, are also likely to evolve according to the progress of knowledge. This is not a weakness insofar as the archaeomagnetic results remain accurate. However, it is important that dates associated with archaeomagnetic measurements reflect revisions to archaeological ages. For some regions there have been recent revisions to GEOMAGIA50 to accommodate new age information, e.g., Bulgaria [START_REF] Kovacheva | Extended and revised archaeomagnetic database and secular variation curves from Bulgaria for the last eight millennia[END_REF], United Kingdom [START_REF] Batt | Advances in archaeomagnetic dating in Britain: New data, new approaches and a new calibration curve[END_REF], Greece [START_REF] De Marco | An updated catalogue of Greek archaeomagnetic data for the last 4500 years and a directional secular variation curve[END_REF], USA [START_REF] Bowles | Archaeomagnetic intensity results from California and Ecuador: evaluation of regional data[END_REF][START_REF] Jones | Archeointensity of the Four Corners Region of the American Southwest[END_REF] and France (Le [START_REF] Goff | An updated archeomagnetic directional variation curve for France over the past two millennia, following 25 years of additional data acquisition[END_REF]. It must also be recognized that there are likely ages within GEOMAGIA50 that do not reflect advances in archaeological age determinations for specific times, regions or sites. Work can continue on sites for years to decades and the archaeomagnetic aspect of the excavation/project may not be the primary objective; new ages may come to light after the final publication of the archaeomagnetic work.

Another common method used to indirectly date archaeological materials is radiocarbon dating. Approximately 15% of entries have used radiocarbon dating as the sole chronological control or in conjunction with other dating methods. Charcoals from carbonaceous or ashy layers that are related to the last use of a kiln/fireplace or located in different horizons of the stratigraphy have frequently been used for dating (e.g., [START_REF] Shaar | Decadal-scale variations in geomagnetic field intensity from ancient Cypriot slag mounds[END_REF], but other materials such as seeds and bones have also been used. In comparison to the typochronological approach, its advantage is to give a precise date bound by experimentally derived uncertainties. However, the significance of this date relative to the TRM acquisition is not guaranteed. For example, the date can be affected by an old carbon/wood effect. Radiocarbon dates the formation of the organic cell and dating charcoals from reused woods or central tree rings can result in earlier dates up to a few centuries. A preliminary anthracological study is useful to identify such samples and select, if possible, materials with a short lifetime as burnt twigs, grasses or seeds.

A limitation of the method is that the abundance of radiocarbon within a sample can not be

J o u r n a l P r e -p r o o f

Journal Pre-proof simply related to a specimen"s age, through comparison to a decay product, as for example in 40 K/ 39 Ar dating; nitrogen produced by the decay of 14 C is not captured by the majority of materials [START_REF] Reimer | The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0-55 Cal kBP)[END_REF] Finally, we underline that the best way to minimize the risk that the true date of the TRM acquisition is not included in the given interval of age is to combine several chronometric and/or archaeological dates. This is often done by an archaeologist who has an overarching understanding of the site and its positioning in the regional fabric. More recently, mathematical techniques such as Bayesian chronological modelling (e.g., Bronk Ramsey, 2009; Lanos and Philippe, 2018) have brought additional insights into developing archaeological chronologies, especially for sites with complex stratigraphies (e.g., [START_REF] Shaar | Geomagnetic field intensity: How high can it get? How fast can it change? Constraints from Iron Age copper slag[END_REF].

Age uncertainties in GEOMAGIA50

Age uncertainties (expressed here as age ranges to accommodate the multimodal age probability distributions of calibrated radiocarbon ages) vary widely within the database, ranging from 0 years for some historically and archaeologically dated entries (1.4% of data) to 2900 years for an archaeologically dated oven from Germany (Schnepp et al., 2020b) (Fig. 11).

Approximately 6% of data (627 entries) have an age range  10 years; 30% have an age range

J o u r n a l P r e -p r o o f

Journal Pre-proof  50 years; and 50% of entries have age ranges of 100 years or less. Nearly all age ranges are less than 500 years ( 90%). There are spikes in the age ranges, with ranges of 100, 200, 300, 400 and 500 being more populous than others (Fig. 11b). The majority of these ranges are from archaeological dated materials and are assignments to specific centuries or across multiple centuries. In general, there is no correlation between age and age range. It is important to note that age ranges can be reported at differing precisions (e.g., 1 or 2 standard deviations) and they do not have the same form. For example, some age distributions will follow a normal distribution (e.g., uncalibrated radiocarbon ages and luminescence techniques), some a multimodal distribution (calibrated radiocarbon ages) and others a uniform distribution (e.g., archaeological ages assigned to a specific archaeological period). For a specific age within an age range, this means there will be differing probabilities of this age depending on the dating method used.

Archaeomagnetic field reconstructions

As described in section 2.3 and section 2.4 archaeomagnetic data are inhomogeneous in space and time. Furthermore, little can be garnered about the large scale geomagnetic field from individual data. Regional or global compilations of data are therefore necessary to gain a greater understanding of the temporal and spatial evolution of the field. This section will give an overview of the two main approaches to reconstructing the geomagnetic field on centennial to millennial time scales: regional secular variation curves and global spherical harmonic models.

Regional secular variation curves

The potential to combine individual archaeomagnetic data from different locations into composite archaeomagnetic curves for dating purposes was recognized in the 1950s (e.g., [START_REF] Cook | Archaeomagnetism. A preliminary report on Britain[END_REF][START_REF] Watanabe | Secular variation in the direction of geomagnetism as the standard scale for geomagnetochronology in Japan[END_REF] and a variety of reference curves have been obtained for several parts of the world since then (see [START_REF] Korte | Refining Holocene geochronologies using palaeomagnetic records[END_REF] for a detailed review). Because the geomagnetic field cannot be considered purely dipolar, field variations at one location (or in one region) are not representative of the evolution of the field as a whole. Smaller-scale non-dipole

J o u r n a l P r e -p r o o f

Journal Pre-proof contributions lead to deviations from a dipolar geometry, resulting in variations in direction and intensity that can vary from one region to another. Combining all global archaeomagnetic data into composite curves will not fully capture the evolution of the field and may obscure regional field structures. Therefore it has been common to develop regional archaeomagnetic curves. It is generally assumed that data within a radius of several 100 to a few 1000 km reflect similar field variations and can be combined to form a reference curve for a region (e.g., [START_REF] Tarling | Geomagnetic secular variation in Britain during the last 2000 years[END_REF][START_REF] Tema | New Italian directional and intensity archaeomagnetic reference curves for the past 3000 years: Insights on secular variation and implications on dating[END_REF]. A review by [START_REF] Korte | Refining Holocene geochronologies using palaeomagnetic records[END_REF] included an investigation of the spatial correlation length of geomagnetic variations and the possible influence of the distance to a curve for dating accuracy. However, strict guidelines cannot be given owing to the complex spatial and temporal evolution of the geomagnetic field over short time scales.

Archaeomagnetic reference curves of field directions, intensity, or all three field components have been developed over a number of decades for several European countries (e.g., Kovacheva et al., 2009a;[START_REF] Tema | New Italian directional and intensity archaeomagnetic reference curves for the past 3000 years: Insights on secular variation and implications on dating[END_REF]Schnepp et al., 2020b,a), Japan (e.g., [START_REF] Watanabe | Secular variation in the direction of geomagnetism as the standard scale for geomagnetochronology in Japan[END_REF][START_REF] Nagata | Secular varaition of the geomagnetic force during the last 5000 years[END_REF][START_REF] Kitazawa | Intensity ,of the geomagnetic field in Japan for the past 10,000 years[END_REF][START_REF] Sakai | Archaeointensity determinations from Western Japan[END_REF], China (e.g., [START_REF] Wei | Intensity of the geomagnetic field near Loyang, China between 500 BC and AD[END_REF][START_REF] Wei | The total intensity of geomagnetic field in southern China for the period from 4500 B.C. to A.D. 1500[END_REF][START_REF] Batt | New archaeomagnetic studies near Xi'an, China[END_REF][START_REF] Yang | A comparison of archaeointensity results from Chinese J o u r n a l P r e -p r o o f Journal Pre-proof ceramics using Thellier"s and Shaw"s palaeointensity methods[END_REF][START_REF] Shaw | Archaeointensity variations for the past 7,500 years evaluated from ancient Chinese ceramics[END_REF], and the United States of America (e.g., [START_REF] Watanabe | Some results of an archaeomagnetic study on the secular variation in the southwest of North America[END_REF]Sternberg, 1989a;[START_REF] Hagstrum | Archeomagnetic dating in western North America: An updated reference curve based on paleomagnetic and archeomagnetic data sets[END_REF][START_REF] Jones | Archeointensity of the Four Corners Region of the American Southwest[END_REF] (see [START_REF] Constable | Centennial-to millennial-scale geomagnetic field variations[END_REF] for a more detailed list with comprehensive references).

Several curves have been frequently updated with new data as they become available, e.g., France (e.g., [START_REF] Thellier | Sur la direction du champ magnetique terrestre, en France, durant les deux derniers millenaires[END_REF][START_REF] Bucur | The direction of the terrestrial magnetic field in France, during the last 21 centuries. recent progress[END_REF][START_REF] Chauvin | Paleointensity of the geomagnetic field recovered on archaeomagnetic sites from France[END_REF][START_REF] Genevey | Intensity of the geomagnetic field in western Europe over the past 2000 years: New data from ancient French pottery[END_REF][START_REF] Genevey | Evidence for rapid geomagnetic field intensity variations in Western europe over the past 800 years from new French archeointensity data[END_REF][START_REF] Genevey | New archeointensity data from French Early Medieval pottery production (6-10th century AD). Tracing 1500 years of geomagnetic field intensity variations in Western Europe[END_REF][START_REF] Gallet | Three millennia of directional variation of the Earth"s magnetic field in western Europe as revealed by archeological artefacts[END_REF]Hervé et al., 2013a,b;[START_REF] Goff | An updated archeomagnetic directional variation curve for France over the past two millennia, following 25 years of additional data acquisition[END_REF]. With efforts to improve data coverage for other regions, there are now curves for China [START_REF] Cai | Recent advances in Chinese archeomagnetism[END_REF], the Near East [START_REF] Gallet | New late neolithic (c. 7000-5000 BC) archeointensity data from Syria. Reconstructing 9000 years of archeomagnetic field intensity variations in the Middle East[END_REF][START_REF] Stillinger | Refining the archaeomagnetic dating curve for the Near East: new intensity data from Bronze Age ceramics at Tell Mozan, Syria[END_REF][START_REF] Shaar | Synchronizing geomagnetic field intensity records in the Levant between the 23rd and J o u r n a l P r e -p r o o f Journal Pre-proof 15th Centuries BCE: Chronological and methodological implications[END_REF][START_REF] Livermore | Archeomagnetic intensity variations during the era of geomagnetic spikes in the Levant[END_REF]), Mexico (Soler Arechalde et al., 2019;[START_REF] Mahgoub | Refined 3600 years palaeointensity curve for Mexico[END_REF] and South America [START_REF] Goguitchaichvili | First archaeointensity reference paleosecular variation curve for South America and its implications for geomagnetism and archaeology[END_REF].

Approaches to curve construction

The first step to building a reference curve is to relocate the distributed data to a central location (also known a reference location) to eliminate differences that result in directions or intensity at different locations purely from a dipole field geometry. For directions, this is commonly done by using the conversion-via-pole (CVP) method, whereby a directional pair with one set of geographic coordinates is transformed to a virtual geomagnetic pole (VGP) and the

J o u r n a l P r e -p r o o f

Journal Pre-proof subsequent VGP is then transformed to a new directional pair using the geographic coordinates selected for the reference curve [START_REF] Shuey | Geographic correction of archeomagnetic data[END_REF][START_REF] Noel | A method for correcting geographically separated remanence directions for the purpose of archaeomagnetic dating[END_REF]. Alternative approaches have also been used, which assume an axial dipole (a dipole aligned with Earth"s rotation axis where the geographic and magnetic pole are coincident) [START_REF] Aitken | Magnetic dating-III: further archaeomagnetic measurements in Britain[END_REF][START_REF] Thellier | Sur la direction du champ magnetique terrestre, en France, durant les deux derniers millenaires[END_REF], though this method is no longer common. In addition, average curves have been calculated

using VGPs and not relocated directions (e.g., in the USA; Sternberg, 1989b;[START_REF] Lengyel | A revision to the U.S. Southwest archaeomagnetic master curve[END_REF]. It is important to note that the non-dipolar nature of the archaeomagnetic field means that all relocation methods have an associated uncertainty [START_REF] Shuey | Geographic correction of archeomagnetic data[END_REF][START_REF] Casas | Distribution analysis of errors due to relocation of geomagnetic data using the "Conversion via Pole" (CVP) method: Implications on archaeomagnetic data[END_REF].

For archaeointensities, curves are typically constructed using intensity transformed to either a virtual dipole moment (VDM) or a virtual axial dipole moment (VADM) (e.g., [START_REF] Daly | An updated and homogeneous world secular variation data base. 1. Smoothing of the archaeomagnetic results[END_REF][START_REF] Yang | Variations in the geomagnetic dipole moment over the last 12000 years[END_REF]. A VDM is analogous to a VGP, as it uses inclination assuming a tilted dipole. VADM assumes a geocentric axial dipole configuration and allows intensity data lacking inclination to be compared. As with relocated directions, there is an intrinsic uncertainty when calculating a V(A)DM assuming a dipole field configuration when the field can have noticeable non-dipolar components. This can result in a dispersion between sites that is a reflection of non-dipolar field behaviour and not related to issues with how estimates of archaeointensity were obtained.

Different data fitting and smoothing methods have been employed to derive regional secular variation curves. Some early curves relied on hand drawn fits through the data (see, [START_REF] Thellier | Sur la direction du champ magnetique terrestre, en France, durant les deux derniers millenaires[END_REF][START_REF] Clark | Developments in archaeomagnetic dating in Britain[END_REF]. However, through time increasingly sophisticated mathematical approaches have been used to construct curves, applying methods that not only derive single curves through time (or through inclination and declination), but calculate uncertainties. Simple interpolation of individual field components (or means across time interval bins) with or without the estimation of curve uncertainties has been common (e.g., Sternberg, 1989a;[START_REF] Yang | Variations in the geomagnetic dipole moment over the last 12000 years[END_REF], but over the past thirty years a variety of mathematical approaches have been taken. Methods such as bivariate extensions of Fisher Statistics (Le Goff et al., 1992) continue to be used to produce curves for Europe, e.g., France (Hervé et al., 2013a;[START_REF] Goff | An updated archeomagnetic directional variation curve for France over the past two millennia, following 25 years of additional data acquisition[END_REF]. Recently bootstrap or Bayesian methods have been used to obtain curves, uncertainty estimates, and/or probability distributions [START_REF] Thébault | A bootstrap algorithm for deriving the archeomagnetic field intensity variation curve in the Middle East over the past 4 millennia BC[END_REF][START_REF] Hellio | Stochastic modelling of regional archaeomagnetic series[END_REF][START_REF] Livermore | Transdimensional inference of archeomagnetic intensity change[END_REF]. The Bayesian method of [START_REF] Lanos | Bayesian Inference of Calibration Curves: Application to Archaeomagnetism[END_REF] and [START_REF] Lanos | Hierarchical modelling of archaeomagnetic data and curve estimation by moving average technique[END_REF] is notable as it produces consistent curves for all J o u r n a l P r e -p r o o f

Journal Pre-proof three field components and provides curve uncertainties that consider uncertainties on the archaeomagnetic data and ages, and the data distribution. This approach has been used in the creation of a number of archaeomagnetic curves across Europe, e.g., Austria and Germany [START_REF] Schnepp | Archaeomagnetic secular variation in Germany during the past 2500 years[END_REF], and Bulgaria [START_REF] Kovacheva | Extended and revised archaeomagnetic database and secular variation curves from Bulgaria for the last eight millennia[END_REF]. The method of [START_REF] Livermore | Transdimensional inference of archeomagnetic intensity change[END_REF], which is published alongside open source code, also produces curve uncertainties for intensity and the posterior sample age distributions as a direct output.

Examples of regional field variations

There are two features noticed in several regions that have received attention in the last decades: archaeomagnetic jerks and intensity spikes. Rapid changes in directional variations seen in Bauer plots (declination against inclination) associated with an increase in intensity in French and Middle East data have been named "archaeomagnetic jerks" [START_REF] Gallet | On the possible occurrence of archaeomagnetic jerks in the geomagnetic field over the past three millennia[END_REF][START_REF] Gallet | Does earth"s magnetic field secular variation control centennial climate change? Earth Planet[END_REF][START_REF] Gallet | Geomagnetic field hemispheric asymmetry and archeomagnetic jerks[END_REF]. [START_REF] Gallet | Does earth"s magnetic field secular variation control centennial climate change? Earth Planet[END_REF] suggested that if archaeomagnetic jerks are global features, they could be associated with episodes of a tilted and enhanced dipole. Later, [START_REF] Gallet | Geomagnetic field hemispheric asymmetry and archeomagnetic jerks[END_REF] noted that they may correspond to maximum geomagnetic field hemispheric asymmetry, leading to most-eccentric dipole events, related to the dynamics of flux patches at mid-to high-latitudes. If archaeomagnetic jerks are regional, they may result from a recurring non-dipole field structure that influences Western Europe. Using the global field model CALS7k.2 [START_REF] Korte | Continuous geomagnetic field models for the past 7 millenia: 2. CALS7[END_REF], [START_REF] Dumberry | Azimuthal flows in the Earth"s core and changes in length of day at millennial timescales[END_REF] inferred that archaeomagnetic jerks are associated with a change in the dominant azimuthal flow direction at the top of the outer core below Europe. It is important to note that there is a clear difference in timescales between archaeomagnetic and geomagnetic jerks (e.g., [START_REF] Mandea | Geomagnetic and archeomagnetic jerks: where do we stand? Eos[END_REF]. Archaeomagnetic jerks do not appear unusually rapid compared to what we know from the present field. An archaeomagnetic jerk may last 100-200 years, whereas a geomagnetic jerk lasts 1 year.

In contrast, the intensity variations during geomagnetic intensity spikes during the Iron Age derived from archaeological materials in the Levant (e.g., [START_REF] Ben-Yosef | Geomagnetic intensity spike recorded in high resolution slag deposit in Southern Jordan[END_REF][START_REF] Shaar | Geomagnetic field intensity: How high can it get? How fast can it change? Constraints from Iron Age copper slag[END_REF], are much faster than field changes observed for recent and historical times. The intensity of the field was also far greater than seen today, exceeding twice today"s field strength [START_REF] Shaar | Large geomagnetic field anomalies revealed in Bronze to Iron Age archeomagnetic data from Tel Megiddo and Tel Hazor, Israel[END_REF], 2018) (Fig 12). [START_REF] Livermore | Archeomagnetic intensity variations during the era of geomagnetic spikes in the Levant[END_REF] suggest that six intensity spikes are required by the Levant data sets. Increases in intensity were also associated with a directional anomaly (most notably in inclination) and the combined directional-intensity anomaly is referred to as the

J o u r n a l P r e -p r o o f

Journal Pre-proof

Levantine Iron Age Anomaly [START_REF] Shaar | The first catalog of archaeomagnetic directions from Israel with 4,000 years of geomagnetic secular variations[END_REF]. As with archaeomagnetic jerks, it is currently unclear whether the Iron Age anomaly is regional or global in extent. Data from areas surrounding the Levant (e.g., Georgia, Turkmenistan, Uzbekistan, Cyprus, Greece, Bulgaria and Egypt)

indicate there is some evidence of an increase in intensity at the same (or similar) times to the Levantine spikes; however, although intensity reaches twice the present day field"s at these locations, the increase spans a broader time (a few hundred years) (Fig. 12). This maybe a true representation of the field behaviour during the Iron Age or it may reflect dating inaccuracies.

There is some evidence globally of an intensity increase coincident with the Iron Age anomaly and of spikes at other times (see, [START_REF] Korte | Archeomagnetic intensity spikes: global or regional geomagnetic field features?[END_REF]. Their origin is under discussion [START_REF] Livermore | Core-flow constraints on extreme archeomagnetic intensity changes[END_REF][START_REF] Davies | Geomagnetic spikes on the core-mantle boundary[END_REF][START_REF] Korte | Archeomagnetic intensity spikes: global or regional geomagnetic field features?[END_REF][START_REF] Troyano | Constraining the axial dipole moment variations during the historical period from new archeointensity results obtained in Bukhara (Uzbekistan, Central Asia)[END_REF] and difficult to explain given our current knowledge of the geodynamo. To further understand the driving mechanisms that generate high intensity spikes more high-quality archaeomagnetic data from several regions are necessary to fully characterize their regional or global behaviour. et al. (2021) as GEOMAGIA50 includes all data and a data selection protocol is not applied. The reader is referred to the above compilations for rationales related to data selection.

Global archaeomagnetic field models

Given its source in Earth"s outer core, the geomagnetic field is a global phenomenon and any studies that aim to decipher its driving processes must consider the global evolution of the field. In addition, variations in global field strength, expressed as a dipole moment are also of interest, e.g., in the context of estimating geomagnetic shielding against solar wind, galactic cosmic ray production, atmospheric ionization and solar activity (e.g., [START_REF] Usoskin | Solar activity reconstructed over the last 7000 years: The influence of geomagnetic field changes[END_REF][START_REF] Usoskin | Role of centennial geomagnetic changes in local atmospheric ionization[END_REF][START_REF] Usoskin | Regional millennial trend in the cosmic ray induced ionization of the troposphere[END_REF][START_REF] Usoskin | Solar activity during the Holocene: the Hallstatt cycle and its consequence for grand minima and maxima[END_REF]. A range of global archaeomagnetic field models have been derived over the past decades, from which maps of the field can be generated for Earth"s surface (e.g., Fig. 13) and core-mantle boundary (CMB) (e.g., Fig. 14). In the following sections we give an overview of the history of global archaeomagnetic field models, how modelling approaches have evolved over the J o u r n a l P r e -p r o o f

Journal Pre-proof past twenty years, and the current state of the art. We discuss how data selection and data uncertainties influence global models and how these have been treated in the most recent models.

We describe some of the major findings that global modelling has facilitated, but also note caveats to the modelling approaches.

Dipole moment reconstructions

According to the International Geomagnetic Reference Field (IGRF 13th generation; [START_REF] Alken | International Geomagnetic Reference Field: the thirteenth generation[END_REF] the present core field is dominated to about 93% by a dipole centered in the middle of the Earth and tilted with respect to Earth"s rotation axis, and to about 91% percent by an axial dipole, i.e. a dipole aligned with the rotation axis. For a purely axially aligned dipole field, the global dipole moment can be determined from a single intensity value and the latitude of the observation. For the moment of a tilted dipole, the inclination at the observation site is additionally required (see, e.g., [START_REF] Merrill | The Magnetic Field of the Earth: Palaeomagnetism, the core, and the deep mantle[END_REF]. However, when non-dipole field contributions are present, any dipole moment values determined in this way are biased depending on the strength of the local non-dipole field. It is often assumed that non-dipole field contributions average out when enough individual VDMs or VADMs are averaged in space and/or time, so that such an average V(A)DM is considered a valid approximation of the actual dipole moment (e.g., [START_REF] Merrill | The Magnetic Field of the Earth: Palaeomagnetism, the core, and the deep mantle[END_REF].

However, the validity of this assumption in unclear and at least for short intervals or a strongly inhomogeneous global data distribution, the resulting averaged V(A)DM is likely biased.

Several V(A)DM reconstructions from archaeomagnetic data (which also in general include volcanic data) span the past 10 to 50 kyr [START_REF] Mcelhinny | Variations in the geomagnetic dipole: I. The past 50 000 years[END_REF][START_REF] Yang | Variations in the geomagnetic dipole moment over the last 12000 years[END_REF][START_REF] Genevey | ArcheoInt: An upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment[END_REF][START_REF] Knudsen | Variations in the geomagnetic dipole moment during the Holocene and the past 50 kyr[END_REF][START_REF] Valet | Secular variation of the geomagnetic dipole during the past 2000 years[END_REF][START_REF] Usoskin | Solar activity during the Holocene: the Hallstatt cycle and its consequence for grand minima and maxima[END_REF]. [START_REF] Genevey | ArcheoInt: An upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment[END_REF] showed that based on the ArcheoInt database, VADM or mixed VADM/VDM curves from Eurasia differ notably from curves for the rest of the world and they constructed global curves for the past 3 kyr using equally weighted regional curves to avoid biasing from a heterogenous data distribution. [START_REF] Knudsen | Variations in the geomagnetic dipole moment during the Holocene and the past 50 kyr[END_REF] used GEOMAGIA50 (version 1) for a VADM reconstruction over the past 50 kyr, in time windows increasing from 500 years for the past 4 kyr to 4 kyr prior to 24 ka, and noted that field strength through the Holocene is higher than during the preceding 40 kyr.

History of archaeomagnetic field models

J o u r n a l P r e -p r o o f

Journal Pre-proof

Although early attempts to construct global archaeomagnetic field models date back to the early 1970s (e.g., [START_REF] Márton | Secular variation of the geomagnetic virtual dipole field during the last 2000 yr as inferred from the spherical harmonic analysis of the available archeomagnetic data[END_REF][START_REF] Braginskiy | Spherical analysis of the geomagnetic field based on archaeomagnetic data[END_REF], they have only received considerable attention over the past 20 years, when the data basis had become large enough to allow for more spatial detail and temporally continuous reconstructions. The recent history of purely archaeomagnetic field models is closely linked to models including palaeomagnetic sediment records in addition to archaeomagnetic and also volcanic data. A surge of models spanning back to 2 to 12 ka followed the publication of [START_REF] Hongre | An analysis of the geomagnetic field over the past 2000 years[END_REF]. This includes a series of 100-year snapshot models for the past 3000 years by [START_REF] Constable | Global geomagnetic field models for the past 3000 years: transient or permanent flux lobes? Phil[END_REF], and its first continuous equivalent [START_REF] Korte | Continuous global geomagnetic field models for the past 3000 years[END_REF]. Several recent reviews include overviews of all these models [START_REF] Constable | Centennial-to millennial-scale geomagnetic field variations[END_REF][START_REF] Korte | Archeomagnetic intensity spikes: global or regional geomagnetic field features?[END_REF][START_REF] Korte | Refining Holocene geochronologies using palaeomagnetic records[END_REF], and we focus on archaeomagnetic models (including volcanic data, but no sediment records) in the following discussion.

The first such models were ARCH3k.1 and ARCH3k_cst.1 [START_REF] Korte | Geomagnetic field for 0-3 ka: 2. a new series of time-varying global models[END_REF] for the time interval 1000 BCE to 1990 CE. ARCH3k.1 was initially based on all available archaeomagnetic and volcanic data that the authors were aware of (9605 values), with iterative outlier rejection. Journal Pre-proof 2019), and a proof-of-concept model for the past 1000 years [START_REF] Mauerberger | Correlation based snapshot models of the archeomagnetic field[END_REF].

Range of modelling approaches

Most global geomagnetic field models, whether covering recent, historical, archaeo-or palaeomagnetic times, are based on series of spherical harmonic (SH) functions that are fit to the data by mathematical inversion techniques. The geomagnetic field is conveniently described by a series of coefficients that scale with field contributions that can be described by a (tilted) dipole, quadrupole, octupole and increasingly shorter wavelength parts. Moreover, when assuming that

Earth"s mantle is electrically insulating, the SH representation can be downward-continued to provide an image of the field morphology at the top of Earth"s outer core, the CMB. For continuous models over certain time intervals the coefficients are smoothly varying time-dependent functions, mostly based on cubic B-splines when constructing historical to millennial scale models (see, e.g., [START_REF] Korte | Continuous global geomagnetic field models for the past 3000 years[END_REF][START_REF] Korte | Geomagnetic field for 0-3 ka: 2. a new series of time-varying global models[END_REF].

As a result of uncertainties in data and age (see section 3) a model cannot and should not fit all data exactly, and some form of smoothing constraint is implemented in the modelling. The simplest form is a truncation of the SH expansion at low degrees to limit the spatial variability of the model and a temporal parameterisation allowing only slow temporal changes. However, most modellers prefer a more flexible form of regularization, where the model parameterisation allows for more variability than expected to be resolved by the data, and the fit to the data is traded off against additional smoothness constraints in space and time.

Methodological differences among most archaeo-and palaeomagnetic SH models mainly lie in the choice and strength of smoothing constraints, and the treatment of outlying data. [START_REF] Hellio | Time-correlation-based regression of the geomagnetic field from archeological and sediment records[END_REF] in a new approach used statistical information about geomagnetic field evolution from satellite and observatory observations in temporal cross-covariance functions as a constraint in a Bayesian modelling frame. The method results in an ensemble of models with statistically coherent errors on the parameters. [START_REF] Arneitz | A Bayesian iterative geomagnetic model with universal data input: Self-consistent spherical harmonic evolution for the geomagnetic field over the last 4000 years[END_REF] also used a Bayesian approach when directly combining archaeomagnetic data with historical observations.

Two recent studies investigate new methods for snapshots in time as a step towards improved continuous models. approach, assuming the geomagnetic potential to be a Gaussian process rather than using SH basis functions. The method provides realistic regional model uncertainties depending on data distribution.

New modelling approaches provide additional relevant information on model resolution and model uncertainties. Models in general agree for regions or parameters that are well constrained by data, as can be seen in the maps of Fig. 13 and Fig. 14, where the different models appear more similar (panels in a) and have smaller uncertainties (panels in b) in the Northern, than the Southern Hemisphere. This can also been seen in time series from Europe and South Africa, where there is a better agreement of the models for Paris (dense data) than South Africa (sparse data) (Fig. 15). Differences in the data basis, outlier treatment and how uncertainties are weighted have a stronger influence on the models than the method used [START_REF] Sanchez | Modelling the archaeomagnetic field under spatial constraints from dynamo simulations: a resolution analysis[END_REF][START_REF] Korte | Archeomagnetic intensity spikes: global or regional geomagnetic field features?[END_REF]. When creating a field model (and assessing site dependent output), the underlying data basis should be considered, especially how well a model is constrained for a certain region and time or for a certain purpose. Moreover, all available models are smoothed representations of the actual field variability in both space and time. The amplitudes of rapid field changes in field models are not fully resolved and are likely underestimated.

Influence of data selection and distribution on global models

Data selection, weighting and distribution have a significant influence on the output of global models. Data selection follows two philosophies. The first is to use all available data without any prior selection, hoping that the signal to noise ratio will increase with the number of available data. The second philosophy is to make a prior data selection by imposing a set of quality criteria. This makes sense when the quality of the data is well understood and the information is available in global databases, e.g., studies from France (Le [START_REF] Goff | An updated archeomagnetic directional variation curve for France over the past two millennia, following 25 years of additional data acquisition[END_REF] or the Levant [START_REF] Shaar | Large geomagnetic field anomalies revealed in Bronze to Iron Age archeomagnetic data from Tel Megiddo and Tel Hazor, Israel[END_REF][START_REF] Shaar | Synchronizing geomagnetic field intensity records in the Levant between the 23rd and J o u r n a l P r e -p r o o f Journal Pre-proof 15th Centuries BCE: Chronological and methodological implications[END_REF]. However, this is currently not the case in many other regions of the world, where results are sparse and/or many of the results have been obtained decades ago, before some of the modern laboratory methods and tests providing modern quality criteria existed. It is worth noting that [START_REF] Korte | Geomagnetic field for 0-3 ka: 2. a new series of time-varying global models[END_REF] performed a comparison of models with and without prior data selection based on data and dating uncertainties and found no notable improvement when data selection was imposed.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Well distributed global data are the most relevant ingredient for an overall good global model. Recent models providing improved uncertainty estimates [START_REF] Sanchez | Modelling the archaeomagnetic field under spatial constraints from dynamo simulations: a resolution analysis[END_REF][START_REF] Hellio | Time-correlation-based regression of the geomagnetic field from archeological and sediment records[END_REF][START_REF] Mauerberger | Correlation based snapshot models of the archeomagnetic field[END_REF] quantify what has been qualitatively stated before (Korte et al., 2009): with the presently available data distribution (more precisely the scarcity of Southern Hemisphere data) archaeomagnetic field models provide limited information about the Southern Hemisphere geomagnetic field.

In Fig. 13 and Fig. 14 we show intensity at Earth"s surface and the radial field at the CMB for one snapshot in time. Models based on archaeomagnetic (and volcanic) data (ARCH10k.1, SHAWQ2k, AmR and COV-ARCH) and archaeomagnetic, volcanic and sediment data (CALS10k.2; Constable et al., 2016) and (COV-LAKE; [START_REF] Hellio | Time-correlation-based regression of the geomagnetic field from archeological and sediment records[END_REF]) can produce models with some broad similarities in intensity at Earth"s surface and radial field for the Southern Hemisphere, e.g., lower intensity patches extending across the Indian Ocean and southern Atlantic ocean. However, the precise locations and morphologies of intensity and radial field patches are different. This can be seen in the model COV-LAKE, which incorporates sediment data and includes more Southern Hemisphere data than its counterpart archaeomagnetic model (COV-ARCH); the use of sediment data results in different global intensity and radial field morphologies. Using sediment data in addition to archaeomagnetic data, but applying the same modelling approach, (e.g., COV-ARCH and COV-LAKE), results in reduced uncertainties on the model output for the Southern Hemisphere (Fig. 13b).

Figure 13: (a) Maps of intensity at Earth"s surface at 900 CE from six global field models:

ARCH10k.1 (Constable et al., 2016), SHAWQ2K [START_REF] Campuzano | Emergence and evolution of the South Atlantic Anomaly revealed by the new paleomagnetic reconstruction SHAWQ2k[END_REF], CALS10k.2 (Constable et al., 2016), AmR [START_REF] Sanchez | Modelling the archaeomagnetic field under spatial constraints from dynamo simulations: a resolution analysis[END_REF], and COV-ARCH and COV-LAKE [START_REF] Hellio | Time-correlation-based regression of the geomagnetic field from archeological and sediment records[END_REF]. (b) Maps of intensity uncertainty for AmR, COV-ARCH and COV-LAKE. Journal Pre-proof coverage is good (a), whereas notable differences exist for South Africa, where there are limited data (b). The included models are the archaeomagnetic models ARCH10k.1 (brown), SHAWQ2k

(yellow), AmR (red) and COV-ARCH (black), and the two models additionally including sediment records CALS10k.2 (green) and COV-LAKE (blue). We note that owing to the regularization applied in global modelling, all curves are reduced in temporal resolution in comparison to regional curves, e.g., intensity curves for Paris [START_REF] Livermore | Transdimensional inference of archeomagnetic intensity change[END_REF].

The lack of Southern Hemisphere data was explicitly considered by earlier versions of the SHA.DIF.14k model, which were European models based on regional rather than global basis functions [START_REF] Pavón-Carrasco | Regional modeling of the geomagnetic field in Europe from 6000 to 1000 B[END_REF]. Few archaeomagnetic (and volcanic) data in the Southern Hemisphere highlight the importance of using sediment records from the Southern Hemisphere to constrain the field in this region. Several other models (not discussed in detail here) mitigate the problem by including high resolution (mainly lacustrine) sediment records (see [START_REF] Constable | Centennial-to millennial-scale geomagnetic field variations[END_REF] and [START_REF] Korte | Refining Holocene geochronologies using palaeomagnetic records[END_REF] for reviews of these).

Major findings

The main applications of purely archaeomagnetic models are for dating purposes and for the calibration of relative intensities obtained from sediments. The advantage of models over regional reference curves for archaeomagnetic dating lies in that models can generate directional and intensity curves for any location without the need for re-location of data. Pavón-Carrasco et al.

(2011) presented a convenient Matlab tool to obtain age probability density functions from any combination of declination, inclination and intensity data. Estimated age ranges tend to be smaller if more than one field component is available. A range of published field models and reference curves are implemented in the published version of the tool, and additional ones can be incorporated by the user.

As absolute field strength cannot be retrieved from sediments, archaeological materials and volcanic rocks are the only sources available for obtaining palaeointensity. Based on our current compilation of archaeomagnetic data, both dipole moment reconstructions and global models show that the dipole moment was high around 2 to 3 ka and greater than today"s field (e.g., [START_REF] Constable | Centennial-to millennial-scale geomagnetic field variations[END_REF]. The Holocene maximum seems high compared to the preceding 40 kyr, as noted by [START_REF] Knudsen | Variations in the geomagnetic dipole moment during the Holocene and the past 50 kyr[END_REF] and the long-term palaeomagnetic average [START_REF] Tauxe | Long-term trends in paleointensity: The contribution of DSDP/ODP submarine basaltic glass collections[END_REF] J o u r n a l P r e -p r o o f

Journal Pre-proof linked to strong dipole moment variations [START_REF] Korte | Archeomagnetic intensity spikes: global or regional geomagnetic field features?[END_REF][START_REF] Hervé | How did the dipole axis vary during the first millennium BCE? New data from West Europe and analysis of the directional global database[END_REF], but their origin is not fully understood.

Studies of global field characteristics, such as symmetry (e.g., Constable et al., 2016) or the field morphology at the CMB [START_REF] Dumberry | Eastward and westward drift of the Earth"s magnetic field for the last three millennia[END_REF][START_REF] Nilsson | Persistent westward drift of the geomagnetic field at the core-mantle boundary linked to recurrent high-latitude weak/reverse flux patches[END_REF], with relevance for the theoretical understanding of the geodynamo, are preferably based on models including sediment records, which provide improved data coverage, in particular for the Southern Hemisphere. Asymmetries seen in the modern field have been found to persist over at least 10 kyr:

the field is weaker, but more variable on average in the Southern Hemisphere compared with the Northern Hemisphere, and secular variation tends to be stronger in the Atlantic and Indian Oceans compared with the Pacific (Constable et al., 2016). Although the magnetic flux morphology at the CMB changes notably with time, there are preferred or recurrent long-term patterns evident in time-averaged models, in particular nearly symmetrical patches of intense flux at high latitudes in both hemispheres (for more details see, e.g., [START_REF] Amit | The time-dependence of intense archeomagnetic flux patches[END_REF] and the review by [START_REF] Constable | Centennial-to millennial-scale geomagnetic field variations[END_REF]). Terra-Nova et al. ( 2017) more recently found recurring positions, but no preferred direction of motion and some correlation of flux evolution with lower mantle heterogeneities, supporting hypotheses of mantle control on the geodynamo [START_REF] Bloxham | Thermal core-mantle interactions[END_REF][START_REF] Bloxham | Time-independent and time-dependent behaviour of high-latitude flux bundles at the core-mantle boundary[END_REF]. Although both westward and eastward azimuthal flow motions seem to occur in the core over archaeomagnetic times [START_REF] Dumberry | Eastward and westward drift of the Earth"s magnetic field for the last three millennia[END_REF][START_REF] Wardinski | The evolution of the core-surface flow over the last seven thousands years[END_REF], recent studies show a clear dominance of westward drift, with rates between 0.07 /yr [START_REF] Nilsson | Reconstructing Holocene geomagnetic field variation: new methods, models and implications[END_REF][START_REF] Nilsson | Persistent westward drift of the geomagnetic field at the core-mantle boundary linked to recurrent high-latitude weak/reverse flux patches[END_REF] and 0.25 /yr [START_REF] Hellio | Time-correlation-based regression of the geomagnetic field from archeological and sediment records[END_REF][START_REF] Nilsson | Persistent westward drift of the geomagnetic field at the core-mantle boundary linked to recurrent high-latitude weak/reverse flux patches[END_REF].

The South Atlantic Anomaly (SAA) is an area stretching from southern Africa over the Atlantic to South America where the geomagnetic field intensity is notably lower than at comparable latitudes. It is known to have deepened and moved westward from about 1700 onwards from historical data [START_REF] Mandea | The magnetic field changing over the southern African continent: a unique behaviour[END_REF][START_REF] Hartmann | Time evolution of the South Atlantic magnetic anomaly[END_REF]. It is linked to the growth of patches of reversed flux at the CMB in the Southern Hemisphere [START_REF] Gubbins | Morphology of the geomagnetic field and implications for the geodynamo[END_REF][START_REF] Terra-Nova | Using archaeomagnetic field models to constrain the physics of the core: robustness and preferred locations of reversed flux patches[END_REF] and has been discussed as a trigger for geomagnetic field reversals [START_REF] Gubbins | Morphology of the geomagnetic field and implications for the geodynamo[END_REF][START_REF] Tarduno | Antiquity of the South Atlantic Anomaly and evidence for top-down control on the geodynamo[END_REF]. A unique connection to reversals is unclear as recent modelling of the field for times prior to the Holocene [START_REF] Brown | Earth"s magnetic field is probably not reversing[END_REF][START_REF] Panovska | One hundred thousand years of geomagnetic field evolution[END_REF] suggests that features similar to the SAA maybe recurrent and do not necessarily lead to reversals. It is of great interest to know the longevity of the SAA, and whether it J o u r n a l P r e -p r o o f

Journal Pre-proof is a recurrent feature of the field, as it may be linked to structures at the CMB that influence core flow and hence geomagnetic field generation [START_REF] Tarduno | Antiquity of the South Atlantic Anomaly and evidence for top-down control on the geodynamo[END_REF][START_REF] Tarduno | Subterranean clues to the future of our planetary magnetic shield[END_REF]. SHAWQ2k

and other models indicate that reverse flux appeared in the Southern Hemisphere as early as 900 CE east of Africa and evolve into the SAA [START_REF] Campuzano | Emergence and evolution of the South Atlantic Anomaly revealed by the new paleomagnetic reconstruction SHAWQ2k[END_REF] (Fig. 14).

Future challenges

Challenges for the future include addressing current inadequacies in the GEOMAGIA50 database, how to improve the temporal and spatial distribution of archaeomagnetic data, and advances in geomagnetic field modelling.

GEOMAGIA50

In previous sections we outlined some of the issues in using data from GEOMAGIA50 for modelling purposes; especially for data selection. One challenge is to homogenize the definition of intensity uncertainties (section 3.1.7) or at least indicate how it was calculated in a new database field. Furthermore, a large number of directional entries are missing k values and some are missing 95

 . All entries lacking these data need to be reassessed and the data added if missing; however, it is likely, especially for k , that the values were not given in the original publications.

A major deficiency in GEOMAGIA50 is the treatment of chronological metadata and we outline these issues and possible solutions in section 5.1.1. The definition of numbers of samples and specimens also require greater clarification, because they can differ for displaced and in-situ archaeological materials, and for lava flows (section 5.1.2).

Keeping the database up-to-date and useful for the scientific community remains a challenge. Given limited resources, it is not feasible to release a new update of the database when each new archaeomagnetic study is published. Instead over coming years, we intend to release an update at the end of each year containing all the new studies published that year. We note, for example, that the recent studies of [START_REF] Shaar | Synchronizing geomagnetic field intensity records in the Levant between the 23rd and J o u r n a l P r e -p r o o f Journal Pre-proof 15th Centuries BCE: Chronological and methodological implications[END_REF][START_REF] Troyano | Constraining the axial dipole moment variations during the historical period from new archeointensity results obtained in Bukhara (Uzbekistan, Central Asia)[END_REF] have yet to be included in version 3.4. of the database, but will be available in the 2021 release.

Archival of chronological data

Section 3.2.3 highlighted how complex the estimation of age uncertainties can be. All

J o u r n a l P r e -p r o o f

Journal Pre-proof methods have caveats and the reliability of the age information is intimately linked to knowledge of the archaeological context. It is a difficult task to develop a hierarchy of dating methods and deduce/calculate dates, while preserving the complexity of the dating process in the database. This may only be partly achievable in future revisions to the database. It requires more metadata which de facto increases the complexity of the database and its searchability. Ideally, additional fields could be added to GEOMAGIA50 giving more specific age information that could be used to sieve data for model or curve construction. However, owing to the range and complexity of dating methods, this is impractical. An alternative would be to add an accompanying text field to each entry describing the chronological controls, e.g. the cultural name or period and a short description of the dating results (as is given in ArcheoInt database [START_REF] Genevey | ArcheoInt: An upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment[END_REF]). This approach would not allow automatic data selection and would require manual assessments of the quality of data prior to modelling. Such functionality is not currently available and would require the assessment of all articles in GEOMAGIA50. [START_REF] Reimer | IntCal13 and Marine13 radiocarbon age calibration age curves 0-50,000 years cal BP[END_REF][START_REF] Reimer | The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0-55 Cal kBP)[END_REF]) (section 3.2.2). Therefore, there are likely some differences in field variations that are not geomagnetic in origin, but rather they relate to changes in calibration curves.

Although there is minimal revision for the past 10,000 years between the generations of calibration curves since 1998 (boundaries of the 95% intervals of date are generally modified by 5-10 years at the most), calibrated ages can be shifted up to several centuries for older periods [START_REF] Reimer | The IntCal20 Northern Hemisphere radiocarbon age calibration curve (0-55 Cal kBP)[END_REF]. Even if changes are small, radiocarbon ages should be recalibrated to keep them mutually consistent and to remove differences in field variations that stem from non-geomagnetic origins.

At a minimum, the experimental radiocarbon ages should be reported accompanied by the calibration method.

The uncertainties on calibrated ages are not treated ideally in the database. Currently only  uncertainties are given. Although this allows asymmetric uncertainties, i.e. the maximum and minimum ages at two standard deviations resulting from calibration, there is no way to record the full multi-modal probability distribution of a calibrated age. This is a significant limitation of the database, as it is important to take into account the irregular shape of the probability density function of calibrated radiocarbon dates in geomagnetic modelling and regional secular variation curves (e.g., [START_REF] Hellio | Stochastic modelling of regional archaeomagnetic series[END_REF][START_REF] Lanos | Bayesian Inference of Calibration Curves: Application to Archaeomagnetism[END_REF][START_REF] Hervé | Improvements in archaeomagnetic dating in Western Europe from the Late Bronze to the Late Iron Ages: An alternative to the problem of the Hallstattian radiocarbon plateau[END_REF][START_REF] Tema | Geomagnetic field secular variation in Pacific Ocean: A Bayesian reference curve based on Holocene Hawaiian lava flows[END_REF]Yutsis-Akimova et al., 2018a). Systematically storing uncalibrated ages (when available) and updating all radiocarbon ages and uncertainties after calibration is a major undertaking, but is an aim for the future.

Clarification of site-sample-specimen hierarchy

Reporting of the number of samples and specimens used for directional and intensity analysis requires evaluation. The database aims to follow the standard palaeomagnetic hierarchy, whereby each entry in the database is considered to be a site, a group of data related to a geological unit or archaeological context that has a unique age. A sample is treated as part of the site that was removed for further analysis. A specimen is a subdivision of a sample and it is this that palaeomagnetic or archaeomagnetic measurements are made on. In palaeomagnetism, this hierarchy works well. A site would be, e.g., a lava flow; a sample, a palaeomagnetic core drilled out of the lava flow; and a specimen, the subdivision of the sample (core) that was measured and the directional or intensity data were obtained from. The number of samples and specimens maybe J o u r n a l P r e -p r o o f

Journal Pre-proof similar, e.g., if one specimen was taken from each sample or the number of specimens could be more if multiple specimens were measured from a sample. In this case, the specimen numbers are averaged to give a sample mean and it is always the sample mean and number of samples that are used to calculate the site mean direction.

In archaeointensity studies this may work differently. For example, a piece of pottery may be related to an instance in time, but might not belong to a context with other pieces of pottery of the same age (a context being a site, analogous, e.g., to a lava flow). In this case the piece of pottery could be treated as both a site or a sample. The piece of pottery can be further divided for measurement and these divisions could be considered to be either samples or specimens. If the piece of pottery were treated as if it were to belong to a context with many other pieces of pottery, then it is just one sample of possibly many, therefore the number of samples would be one and there would be multiple specimens. If the piece of pottery is treated as a site, then the number of samples would be multiple and equal to the number of specimens.

In archaeomagnetic directional studies that use "in situ" structures, i.e. ones that have remained in their original position since firing, there is obviously no ambiguity; a site can be equated with the structure it derives from. In this case, and for directional determination, one should emphasize the importance of working on several independently oriented samples for deriving a mean direction.

Both approaches have their advantages, depending on how the data are to be used or from a consistency point of view (rigid use of the site, sample, specimen hierarchy). For example, obtaining multiple measurements of intensity from a piece of pottery can provide an accurate mean intensity determination, valuable for secular variation or curve construction. This would require that the piece of pottery is treated as a site and that the number of specimens is treated as the number of samples. In the database, the number of archaeointensity measurements used for the mean value is therefore treated in the same way as the number of samples from, e.g., a lava flow: they would have equal value. However, this leads to a mismatch in how different archaeological entries may be treated in the database (e.g., an in-situ structure versus a piece of pottery) and between archaeological and palaeomagnetic hierarchies, e.g., a lava flow would no longer be equal to a context, if the piece of pottery becomes the site. Treating a single artefact as belonging to a single context, regardless of whether there are other artefacts maintains the logic of the hierarchy, but may result in data not being included in further analyses if the number of samples is listed as J o u r n a l P r e -p r o o f one, e.g., if the data are filtered by the number of samples.

The fact that there is no common system of hierarchy is a current weakness of the database.

In the future, we propose (as in ArcheoInt) that N corresponds to the number of thermal units, i.e.

the number of units that can be considered to have been magnetized at the same time. It is fixed to 1 for a lava flow, an archaeological in-situ structure and a single fragment of pottery, or is higher than 1 for a group of baked clay fragments. A second number, n, would be the number of individual values from which the average and its uncertainty are calculated. Because the averages are not calculated homogeneously between data entries, either at the sample level or the specimen level, this solution has the disadvantage to mix samples and specimens, but it would clearly make it easier for data selection.

Improvements in data distribution

As noted in Section 4.2.4 the most important factor in improving global models and our overall understanding of the global archaeomagnetic field is more data from regions with sparse data distributions. With 50% of all data coming from Europe (Section 2.3), there is significant room for improving data coverage. Africa is a clear target given its large area and rich archaeological history. Studies on Burkina Faso, Ivory Coast, Mali, Ethiopia, Kenya, Zimbabwe and South Africa over the past 10 years have made significant steps in improving data coverage.

However, an increased emphasis should be placed on developing archaeomagnetic research projects in Africa. Similarly, South America is well positioned to provide useful data. Data from both southern Africa and South America will be key to unraveling the long term evolution of the South Atlantic Anomaly. The Indian subcontinent has a rich archaeological history and yet no archaeomagnetic directional data have been produced. Given its unique position, new directional data could aid in understanding of westward drift of the field across the Indian Ocean. Australia, New Zealand and Pacific islands also have the potential to provide further archaeomagnetic data.

This would be especially valuable as it could improve our understanding of field variations in the Southern Hemisphere and the Pacific, which are greatly lacking in data. It would be interesting to investigate/confirm the persistence of lower field variability in the western Pacific (Constable et al., 2016). Finally, the large amount of Japanese directional data not in the database should be added. There are hundreds of entries in this data set and it will be of great interest to see how this influences our understanding of field evolution in eastern Asia.

J o u r n a l P r e -p r o o f

Journal Pre-proof

Acquiring data in regions with few data presents practical challenges. Access to in-situ structures is crucial for full vector studies. However, sampling must take place shortly after excavation and this is not always practical. In most countries archaeomagnetists are limited by several constraints (e.g., travel, time, funding, and export licences). This explains why the spatial distribution of directions is fairly close to palaeomagnetic laboratories, as is the case with Europe.

A solution would be to develop local laboratories and/or networks of researchers trained in archaeomagnetic sampling, as well as to collaborate with and train local archaeologists. However, such efforts take time to implement, so we may only see a gradual increase in the amount of data from poorly represented areas.

Another aspect of improving the data distribution is to extend the database further back in time through the Iron Age, Bronze Age and the Neolithic. There is a tendency to focus on more dramatic field changes, e.g., spikes in intensity during the transition from the Iron to Bronze age in the Levant and surrounding areas; however, all times (or "quiet times") are equally valuable to study, as all field variations relate to the underlying geodynamo process. Furthermore, more detailed descriptions of field evolution through time will allow for the development of more accurate archaeomagnetic dating curves.

A consideration of the types of archaeological materials used to extend our knowledge of field variations to older times is also required. Baked clays are less frequent back through time, as are in situ structures, which are more likely to suffer eventual post-displacements. The baking degree of older clay-based materials is also usually lower, resulting in a less stable mineralogy prone to alteration and therefore less favourable for obtaining archaeointensity results. Increasing the success rate of archaeointensity experiments on such materials is a major challenge that could be overcome in the next decade by new approaches, such as scanning magnetometry and computed tomography [START_REF] De Groot | Determining individual particle magnetizations in assemblages of micrograins[END_REF]. Furthermore, a better understanding of the magnetic mineralogy of archaeological materials will aid in this research (e.g., [START_REF] Lopez-Sanchez | Occurrence of epsilon iron oxide in archeological materials: Effect on the natural remanent magnetization[END_REF].

To recover intensity variations in the early Holocene or even in the Pleistocene prior to apparition of ceramic production, an alternative to baked clays is the study of heated rock artefacts, such as burnt cherts [START_REF] Kapper | New directional archaeomagnetic data of burned cave sediments from Switzerland and geomagnetic field variations in Central Europe[END_REF][START_REF] Zeigen | Archaeomagnetism of burnt cherts and hearths from Middle Palaeolithic Amud Cave, Israel: Tools for reconstructing site formation processes and occupation history[END_REF].

Beyond the acquisition of data for older periods and/or for regions that are still poorly documented, an important challenge in archaeomagnetism remains to better understand the issue of data dispersion. Dispersion is characteristic for many data sets and hinders our ability to finely J o u r n a l P r e -p r o o f

Journal Pre-proof trace geomagnetic field variations through time. Besides age uncertainty, the sources of dispersion are more numerous for intensity data (undetected alteration, MD effects, uncorrected TRM anisotropy and cooling rate) than for directions. One major issue is the cooling rate correction because its absence can potentially result in a systematic overestimation. One could think to apply a correction factor to uncorrected data. On average, the correction factor seems to be 5-10% [START_REF] Genevey | ArcheoInt: An upgraded compilation of geomagnetic field intensity data for the past ten millennia and its application to the recovery of the past dipole moment[END_REF]. However, defining a suitable rate of correction is difficult, because it depends on the specific rock magnetic properties of a specimen and the equipment and protocols used in each laboratory.

Furthermore, what is the precision with which the direction and/or intensity of the geomagnetic field can be retrieved? Does dispersion reflect our current limitations in the acquisition of data? To better constrain these questions, one can ask whether there would be an interest in revisiting regions with a good data coverage and an active archaeological research, in order to acquire new precisely dated data, for example from the past few centuries where chronological constraints can be extremely tight. By limiting the influence of age uncertainties, this could help solve the above issues, which are crucial in using archaeomagnetism for archaeological purposes and for refining our knowledge of the evolution of the geomagnetic field.

Global geomagnetic field modelling

From a field modelling point of view, the above-mentioned improvements to the underlying data basis are paramount for improving the temporal resolution and full global spatial reliability of models. Both these aspects are relevant when using field models to infer geodynamo processes in the core, or when using their predictions for regional reference curves. As noted in section 4.2.4, the lack of archaeomagnetic data in the Southern Hemisphere and equatorial areas can be partly compensated for by using sediment data. Similarly, a recently renewed interest in speleothems may produce new high resolution time series for the Holocene in the coming years from locations where it is not possible to obtain archaeological or sediment data. [START_REF] Lascu | Speleothem magnetism[END_REF] give a detailed overview of the potential for speleothems to recover detailed field variations, with the study of [START_REF] Trindade | Speleothem record of geomagnetic South Atlantic Anomaly recurrence[END_REF] providing a detailed Holocene record from

Brazil and other studies resolving other geomagnetically interesting times in great detail (e.g., [START_REF] Lascu | Age of the Laschamp excursion determined by U-Th dating of a speleothem geomagnetic record from North America[END_REF][START_REF] Chou | Multidecadally resolved polarity oscillations during a geomagnetic excursion[END_REF].

Continuing the efforts to improve the treatment of data and dating errors and translating J o u r n a l P r e -p r o o f

Journal Pre-proof them into realistic model errors through methodological developments is also of interest for both these cases. Improved global archaeomagnetic field models may contribute to answering open questions about, e.g., the maximum possible rate of geomagnetic field change and the influences of lowermost mantle structure on the geodynamo (which is reflected in magnetic field morphology). They will also likely contribute to improved predictions of future geomagnetic field evolution by assimilation of data-based models into numerical simulations (e.g., [START_REF] Fournier | An introduction to data assimilation and predictability in geomagnetism[END_REF][START_REF] Tangborn | Impact of archeomagnetic field model data on modern era geomagnetic forecasts[END_REF].

Figure 2 :

 2 Figure 2: Geographical distribution of archaeomagnetic sites given in the six versions of GEOMAGIA50 to date (date range from 8000 BCE to 2000 CE; archaeomagnetically dated sites are excluded). Colours denote when the data were added to the database. (a) sites with directional data; (b) sites with intensity data. Some version 1 (v1) sites were updated with directional information and are shown in (a) as belonging to v1, although they were updated after the initial release of v1. If sites were removed during revisions of subsequent versions, they are not shown on the figure.

Figure 3 :

 3 Figure 3: Geographical distribution of archaeomagnetic sites given in the six versions of GEOMAGIA50 to date for Europe and surrounding regions (date range from 8000 BCE to 2000 CE; archaeomagnetically dated sites are excluded). (a) sites with directional data; (b) sites with intensity data. See Fig 2 for legend and other details.

Figure 4 :

 4 Figure 4: Pie charts of the number of archaeomagnetic entries (in brackets) in GEOMAGIA50.v3.4 by region (a-c) and by country (d-f) by data type: (a,d) directional and

Fig. 2

 2 Fig. 2 we show the increase of Southern Hemisphere data with each new version of GEOMAGIA50. Notable studies have that have obtained data from southern Africa are Neukirch et al. (2012),Tarduno et al. (2015) and Hare et al. (2018). Previously only one study had published

Figure 5 :

 5 Figure 5: Archaeomagnetic entries in GEOMAGIA50.v3.4 by age in 100 year bins. (a) directions, (b) intensity, (c) Northern Hemisphere data and (d) Southern Hemisphere data.

Figure 6 :

 6 Figure 6: Pie charts of the number of mentions (in brackets) of the archaeomagnetic materials used to determine (a) direction and/or intensity, (b) directions, (c) intensity. The eight most used materials are shown in each subplot, remaining material types are grouped under "Other". Note in (b), most inclination data only come from displaced bricks, making the assumption that they were fired on one of their sides. The number of directional and intensity entries do not match the number of materials given in the plots, as numerous entries were determined from multiple materials.

  values of k (the precision parameter;[START_REF] Fisher | Dispersion on a sphere[END_REF] into the thousands. Conversely, some 95

Figure 7 :

 7 Figure 7: Measures of uncertainty and precision (Fisher, 1953) on archaeomagnetic directional and intensity entries within GEOMAGIA50.v3.4: (a) 95% cone of confidence ( 95  ) (bin size = 1 degree); (b) precision parameter ( k ) (bin size = 100 k ). Only 95  values < 20 are shown, corresponding to 5542 entries or 99% of all entries with an 95

Figure 8 :

 8 Figure 8: Pie charts of the number of entries (in brackets) within GEOMAGIA50.v3.4 associated with different palaeointensity methods. (a) All palaeointensity methods (note that the total number of entries exceeds 5231 as multiple palaeointensity entries were derived from measurements using one or more methods). Thellier methods by type are given in (b): Original Thellier-Thellier

Figure 9 :

 9 Figure 9: Pie charts of the intensity entries in the GEOMAGIA5.v3.4 database noting (a) remanence anisotropy corrections and (b) cooling rate corrections. Number of uses of an approach are given in brackets. In (a) TRM = thermoremanent magnetization; ARM = anhysteretic remanent magnetization; IRM = isothermal remanent magnetization; NRM lab B parallel = laboratory field

  which itself is constrained by elements of absolute dating. The different types of dating methods are clearly complementary and the quality of the two approaches cannot be simply ranked, i.e. scientific dating does not always outrank archaeological observations, it depends on the specific context and an understanding of an archaeological site. The importance of sampling in close collaboration with an archaeologist is paramount for selecting materials whose TRM acquisition can be dated with the maximum precision and confidence. Two categories of methods to date J o u r n a l P r e -p r o o f TRM acquisition are distinguished here, either direct, i.e. directly concerning the analyzed material itself, or indirect, i.e. the material is dated by association with another dating element.

Figure 10 :

 10 Figure10: Pie chart of the 8 most commonly used methods to date archaeological materials in GEOMAGIA50.v3.4. Note, the numbers in brackets do not sum to the total number of entries in the database, as numerous entries have been dated using multiple methods. "Other" methods include whether accelerator mass spectrometry was used to obtain radiocarbon ages (133 entries, with frequent overlap with the calibrated and uncalibrated radiocarbon age entries), and if

Figure 11 :

 11 Figure 11: Archaeomagnetic age ranges within GEOMAGIA50.v3.4, binned by (a) 100 year age ranges and (b) 10 year age ranges. (a) the full span of age ranges; (b) truncated to age ranges  500 years.

Figure 12 :

 12 Figure 12: Virtual axial dipole moment (VADM) for the Levant and surrounding regions calculated from data in GEOMAGIA50.v3.4. Note that the data shown in (a) do not reflect the most recent interpretations (compilations) of Shaar et al. (2016), Shaar et al. (2020) and Livermore

ARCH3k_cst. 1

 1 was based on a smaller data set (6211 values) with prior data selection, excluding data with directional uncertainty 95 et al. (2013) presented model A_FM based on 9660 data, spanning 1000 BCE to 2000 CE. Similar to the two ARCH3k models, this was part of a study comparing archaeomagnetic data only models to models including sediment records. A model derived with the main purpose of archaeomagnetic dating is SHA.DIF.14k (Pavón[START_REF] Pavón-Carrasco | A geomagnetic field model for the Holocene based on archaeomagnetic and lava flow data[END_REF], spanning nearly the past 14 kyr based on 12779 data and following from a series of regional European models by the same group (e.g.,[START_REF] Pavón-Carrasco | A regional archeomagnetic model for Europe for the last 3000 years, SCHA.DIF.3K: Applications to archeomagnetic dating[END_REF]. ARCH10k.1 was derived mainly as a starting model for a reconstruction including sediment records for 8000 BCE to 1990 CE(Constable et al., 2016). A model with somewhat improved Southern Hemisphere data coverage due to recent efforts to improve the global data coverage and a data weighting scheme according to archaeomagnetic quality criteria, named SHAWQ2k, was presented by[START_REF] Campuzano | Emergence and evolution of the South Atlantic Anomaly revealed by the new paleomagnetic reconstruction SHAWQ2k[END_REF]. New modelling methods were explored for models AmR, spanning 1200 BCE to 2000 CE in 40-year snapshots[START_REF] Sanchez | Modelling the archaeomagnetic field under spatial constraints from dynamo simulations: a resolution analysis[END_REF], COV-ARCH, a continuous model for the past 3 kyr[START_REF] Hellio | Time-correlation-based regression of the geomagnetic field from archeological and sediment records[END_REF], BIGMUDI4k.1, an iterative approach simultaneously inverting palaeomagnetic, archaeomagnetic and historical records for the past 4,000 years(Arneitz et al., J o u r n a l P r e -p r o o f

  [START_REF] Sanchez | Modelling the archaeomagnetic field under spatial constraints from dynamo simulations: a resolution analysis[END_REF] use statistics from a numerical dynamo simulation as mean and covariance background constraints, which avoids subjective choices of regularization parameters and provides an improved understanding what global spatial resolution J o u r n a l P r e -p r o o f Journal Pre-proof can be retrieved from the data. Mauerberger et al. (2020) implemented a Bayesian non-parametric

Figure 14 :

 14 Figure 14: Maps of (a) the radial field (Br) and (b) its uncertainty for the core-mantle boundary at 900 CE for different field models. See Fig. 13 for model references.

Figure 15 :

 15 Figure 15: Time series of magnetic declination (D, top panels), inclination (I, middle panels) and intensity (F, bottom panels) for the past 3000 years at two locations as predicted by six different global magnetic field models. All models agree closely most of the time for Paris, where data

Figure 1 Figure 2 Figure 3 Figure 4 Figure 5

 12345 Figure1

Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Figure 11 Figure 12 Figure 13 Figure 14

 67891011121314 Figure 6

Figure 15

 15 Figure15

  Version 1 of GEOMAGIA50 primarily focused on compiling palaeointensity data and contained data from both archaeological materials and lava flows. Directional data were added only if they accompanied intensity data. Version 1 integrated the ArcheoInt database of Genevey

	et	al.	(2008)	and	the	IAGA	ARCHEO00	database

(http://www.ngdc.noaa.gov/geomag/paleo.shtml) compiled by Don Tarling. Data from other

Journal Pre-proof general

  decrease with age. This is most stark for BCE data, with 35% of all entries from this time.The number of BCE directions is substantially less (20% of total directions) than for CE (Common Era) directions. The contrast is less abrupt for archaeointensity data. Although the number of BCE intensity entries per century is in general less than for CE entries, 54% of all intensity data span 8000 BCE to 1 BCE.There are notable spikes in the number of directional and intensity entries for certain time periods. For directional data there are peaks in the number of directional data between 100 CE and 300 CE, 700 CE and 900 CE, 1100 and 1400 BCE, and 1700 BCE and 2000 (Fig.5a). The most populous century for directional results is the 19th (410 entries from 31 studies). Some peaks can be attributed to certain cultural periods, e.g., the high number of entries between 100 CE and 300 CE are from the peak of the Roman Empire, with the data set dominated by entries from present

day England, France, Hungary, Bulgaria and Spain. Other peaks are associated with concerted research initiatives in specific countries (or by certain research groups with dedicated focuses), e.g., the 700 CE to 900 CE peak is dominated by data from France for the High Middle Ages (Le

  Calibration can result in significant shifts in age. It is therefore important that calibration is clearly documented. As the database stands, there are still uncalibrated ages that are used for the date of an entry. An initial goal will be to reexamine all the uncalibrated ages and, if there is enough information available, calibrate those ages. As noted in section 3.2, uncalibrated ages are just one of a number of methods that may have been used at a site, e.g., uncalibrated ages may have been combined with archaeological ages or relative stratigraphic ages. In these cases, the dates can not be simply recalibrated, as the radiocarbon ages may form only one aspect of the structural framework for the stratigraphy. This requires that every paper with radiocarbon ages is reexamined in detail to see what scope there is for age re-evaluation.Journal Pre-proof calibration methods since this time and new calibration curves are published every few years (e.g.

Although attempts at incorporating greater radiocarbon information were made in version 2 of the database, there are numerous complications to successful implementation. Radiocarbon dates are currently entered in accordance with information provided in published articles, but not consistently and sometimes ambiguously. Ages maybe calibrated or not and calibrated dates can be reported with symmetrical or asymmetrical bounds (reflecting the calibrated age"s non-gaussian distribution) and at 68 or 95% confidence. It is not always clear at what level uncertainty is reported. Documenting calibration of radiocarbon dates is also problematic. GEOMAGIA50 contains studies from the 1960s onwards and numerous improvements have been made in J o u r n a l P r e -p r o o f

Acknowledgments

We thank the archaeomagnetic community for providing data when requested, especially Elisabeth Schnepp, Patrick Arneitz and Mary Kovacheva. We acknowledge Fabio Donadini, Kimo Korhonen and Lauri Pesonen for the development and design of GEOMAGIA50. GH was supported by the Labex LsScArBx (Bordeaux Archaeological Sciences Labex, ANR-10-LABX-52). Development of GEOMAGIA50 was supported by DFG SPP-1488. This is Institute for Rock Magnetism contribution 2008.

Conflict of interest

I declare no conflicts of interests on behalf of all authors.