Marie Troyano
email: troyano@ipgp.fr

Yves Gallet

Agnès Genevey

Vladimir Pavlov

Alexandre Fournier

France Lagroix

Makhsuma Niyazova

Dzhamal Mirzaakhmedov

Analyzing the geomagnetic axial dipole field moment over the historical period from new archeointensity results at Bukhara (Uzbekistan, Central Asia)

Keywords: Archeomagnetism, Secular variation, Geomagnetic field intensity, Geomagnetic field modeling, Axial dipole moment

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Variations of Earth's magnetic dipole cover a wide range of timescales from a year or less to tens of millions of years. Three different frequency bands are evidenced by analyses of the dipole power spectrum from paleo-and geomagnetic data and simulations [START_REF] Constable | A paleomagnetic power spectrum[END_REF]Ziegler et al., 2011;[START_REF] Olson | From superchrons to secular variation: a broadband dynamo frequency spectrum for the geomagnetic dipole[END_REF][START_REF] Panovska | Observed periodicities and the spectrum of field variations in Holocene magnetic records[END_REF][START_REF] Bouligand | Frequency spectrum of the geomagnetic field harmonic coefficients from dynamo simulations[END_REF]Lesur et al., 2018): an ultra-low to low frequency band (UF), a transitional frequency band (TF), and a high frequency band (HF). The UF band comprises chrons and superchrons and is associated with the thermal evolution of the outer core. The TF band covers paleo-/archeomagnetic secular variations and is associated with geodynamo processes. Finally, the HF band contains the shortest periodicities of the axial dipole's variations (as observed from satellite data). These bands are separated by two cut-off frequencies T s (between UF and TF) and T f (between TF and HF), estimated by [START_REF] Hellio | Time-correlation-based regression of the geomagnetic field from archeological and sediment records[END_REF] from recent field statistics as T s = 100 kyr and T f = 60 yr, for the purpose of constructing the COV-ARCH model (more on global models below).

The axial dipole's power spectrum from numerical dynamo simulations corroborates these results [START_REF] Olson | From superchrons to secular variation: a broadband dynamo frequency spectrum for the geomagnetic dipole[END_REF][START_REF] Bouligand | Frequency spectrum of the geomagnetic field harmonic coefficients from dynamo simulations[END_REF], although the estimated characteristic timescale T f is longer (T f ∼ 10 2 -10 3 yrs), which is probably associated with the convective timescale in the outer core of order 150 yr. While secular variations recovered from global archeomagnetic models are representative of the low-frequency TF band, regional variation curves spanning the last few millennia based on high-quality archeomagnetic data could be associated with the high-frequency band, on time scales on the order of the convective turnover time (e.g. [START_REF] Genevey | New archeointensity data from French Early Medieval pottery production (6th-10th century AD). Tracing 1500 years of geomagnetic field intensity variations in Western Europe[END_REF][START_REF] Genevey | Refining the high-fidelity archaeointensity curve for Western Europe over the past millennium: analysis of Tuscan architectural bricks (Italy)[END_REF].

Studying past field variations requires the construction of time-dependant global field models from the compilation of direct (or instrumental) and/or indirect geomagnetic field measurements. One of the most widely used mod-els is the gufm1 model, which covers the past 400 years [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF] from 1590 to 1990, and which was constructed from a large set of direct geomagnetic measurements obtained in land-based observatories and by mariners during their voyages across the seas (e.g. [START_REF] Jonkers | Four centuries of geomagnetic data from historical records[END_REF], as well as from satellite data for the most recent period. However, our ability to instrumentally measure geomagnetic field intensities only dates back to the 1830s (Gauss, 1833). To overcome this lack of intensity data, [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF], following [START_REF] Barraclough | Spherical harmonic analyses of the geomagnetic field for eight epochs between 1600 and 1910[END_REF], impose a linear decay rate of 15 nT/yr to the axial dipole component between 1590 and 1840, i.e. a rate corresponding to a crude extrapolation back in time of the behavior observed since ∼ 1840. Since it is essential for the construction of the gufm1 model, and in general for our knowledge of geomagnetic field behavior during the historical period, this crude extrapolation has been tested against archeo-/paleointensity data (i.e. indirect measurements) provided by the study of the thermoremanent magnetization carried by archeological artifacts and volcanic deposits (e.g. [START_REF] Gubbins | Fall in Earth's magnetic field is erratic[END_REF][START_REF] Finlay | Historical variation of the geomagnetic axial dipole[END_REF]Genevey et al., 2009;[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] Suttie | Consistent treatment of errors in archaeointensity implies rapid decay of the dipole prior to 1840[END_REF][START_REF] Poletti | Continuous millennial decrease of the Earth's magnetic axial dipole[END_REF]. [START_REF] Hulot | Uniqueness of mainly dipolar magnetic fields recovered from directional data[END_REF] indeed establish that the geomagnetic field can be recovered from directional data alone, up to a constant multiplier (the uniqueness of the sought-after solution being guaranteed by the existence of two, and only two, poles at Earth's surface). The multiplicative constant is in practice provided by independent intensity measurements, each Gauss coefficient entering the mathematical description of the field being renormalized to account for the intensity measured at the specific location of interest. [START_REF] Gubbins | Fall in Earth's magnetic field is erratic[END_REF] follow this line of reasoning and this is the first study to use the set of indirect intensity data available between 1590 and 1840 to recalibrate the axial dipole component provided by gufm1 by the ratio of measured to predicted intensities at intensity determination sites.

Due to scattered data, they assume that a linear fit is indeed the most reasonable solution prior to 1840, but estimate that the axial dipole component between 1590 and 1840 had a rate of decrease of 2.28 ± 2.72 nT/yr, which is significantly lower than that proposed by [START_REF] Barraclough | Spherical harmonic analyses of the geomagnetic field for eight epochs between 1600 and 1910[END_REF] and used by [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF] (15.46 nT/yr and 15 nT/yr respectively).

Next, [START_REF] Finlay | Historical variation of the geomagnetic axial dipole[END_REF] combines both direct and indirect geomagnetic measurements to calculate a new geomagnetic field model between 1590 and 1840, without imposing a linear decrease in the axial dipole during this period (but with an artificial overweighting of the indirect records). He shows that this approach does not provide better results than those favoring no change in axial dipole during the 17th and 18th centuries. [START_REF] Suttie | Consistent treatment of errors in archaeointensity implies rapid decay of the dipole prior to 1840[END_REF] propose a radically different approach based on the statistical analysis of errors in the archeo-paleointensity data. In particular, the dataset available between 1840 and 1990 is used to estimate reasonable errors in the data, which are best assigned as fractions (∼ 15%) of the field intensity values expected from gufm1. When applied to data prior to 1840, and again assuming a linear evolutionary trend in axial dipole over this period, they find a rate of decay (∼ 11.9 nT/yr) close to what [START_REF] Barraclough | Spherical harmonic analyses of the geomagnetic field for eight epochs between 1600 and 1910[END_REF] found. They further show that if data errors are assigned as fractions of measured intensities, the decay rate is similar to that proposed by [START_REF] Gubbins | Fall in Earth's magnetic field is erratic[END_REF] and [START_REF] Finlay | Historical variation of the geomagnetic axial dipole[END_REF] (i.e., with either a slight change or no change at all in the axial dipole component over the 17th and 18th centuries). How-ever, this observation is the result of a bias toward lower field values, as their uncertainties are lower when given as a proportion of measured intensities.

For the different methods above, dispersion of paleo-archeomagnetic data is such that it prevents overcoming the assumption of a linear evolution of the axial dipole component over the historical period. In addition, [START_REF] Suttie | Consistent treatment of errors in archaeointensity implies rapid decay of the dipole prior to 1840[END_REF] demonstrate that the use of quality criteria on the dataset does not significantly change the conclusions. More recently, [START_REF] Poletti | Continuous millennial decrease of the Earth's magnetic axial dipole[END_REF] also use a selected global dataset with strict criteria covering the historical period . After converting intensity data into corresponding axial dipole moments and performing linear regression computations for datasets covering various time intervals, they reach a conclusion favoring a linear decreasing trend of the axial dipole over the historical period of ∼ 12.5 nT/yr, thus close to that advocated by [START_REF] Barraclough | Spherical harmonic analyses of the geomagnetic field for eight epochs between 1600 and 1910[END_REF] and [START_REF] Suttie | Consistent treatment of errors in archaeointensity implies rapid decay of the dipole prior to 1840[END_REF].

Given the dispersion observed in the global compilation of intensity data regardless of the selection criteria considered, Genevey et al. (2009) explore a different approach using a single consistent regional intensity dataset to recalibrate the g 0 1 coefficient of gufm1. The principle remains the same as above [START_REF] Hulot | Uniqueness of mainly dipolar magnetic fields recovered from directional data[END_REF], which assumes that the geometry of the geomagnetic field as provided by gufm1 is correct. While it potentially avoids the problem of global data scatter, and the almost insoluble issue of selecting only the most reliable data, it does raise the pending issue of which dataset is sufficiently reliable to be used to recalibrate the Gauss coefficients (an evaluation that will surely vary from one author to another). Genevey et al. (2009) use the set of accurate and precisely dated archeointensity results obtained in western Europe (700 km around Paris, France). Instead of a linear decrease of the axial dipole magnitude over the historical period, they find a significant decrease between ∼ 1590 and the second half of the 18th century, with a minimum magnitude during this period, followed by a moderate increase from ∼ 1800 to ∼ 1840 and then, the well-established linear decrease up to the present. As a follow-up to this first study, [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] analyse precisely dated architectural brick fragments from southern and northern Brazil. Despite a significant non-dipole field effect between these two regions associated with the South Atlantic Anomaly (SAA), the results obtained appear to support the evolution in dipole field moment proposed by Genevey et al. (2009). As a new development, the present study carried out in Central Asia (Bukhara, Uzbekistan) focusing on the 1590 to 1850 period aims to further constrain the accuracy of the non-linear dipole moment evolution deduced from the western European dataset.

Historical context and sampling

Situated on the Silk Road, Bukhara (39.8 • N, 64.5 • E, Fig. 1) has long been an important place for trade, Islamic education and religion in Central Asia, as evidenced by the many madrasas (or religious schools) and mosques still standing in the city's historic center. These old buildings were built throughout the medieval period from fired clay bricks. Their history, and more generally that of the city itself, is well known through abundant written testimonies preserved in the state archives of Uzbekistan.

[Figure 1 about here.]

For the period covered by our study, the sampled buildings were erected during three successive dynasties that ruled Bukhara from the mid-16th century to the beginning of the 20th century: the Shaybanid dynasty during the 16th century, the Djanid dynasty from the 17th to the mid-18th century, and the Manghit dynasty from the mid-18th century to the early 20th century.

The Shaybanid dynasty, which claimed to be descended from Genghis Khan, conquered Bukhara from the Timurids in the early 16th century and founded the khanate of Bukhara. Their domination for ∼ 100 years was interrupted by the Djanid dynasty (which also claimed to be descended from Genghis Khan), which then established its rule over Bukhara for about a century and a half. The Shah of Iran (Nader Shah) conquered the khanate around the mid-18th century, but the collapse of his empire a few years later led to the establishment of the Manghit dynasty. This dynasty was of Uzbek origin and ruled the Emirate of Bukhara until 1920 when Soviet Red Army troops invaded the city.

Throughout the above period, the city of Bukhara was divided into several small social units called guzars, whose history is well documented in the archives. Each guzar was led by a chief (aqsaqal) nominated by the elders and had its own mosque around which the community was structured [START_REF] Khalid | The Residential Quarter in Bukhara before the Revolution (The Work of O. A. Sukhareva)[END_REF]. [START_REF] Sukhareva | The neighborhood community of the late-feudal city of Bukhara[END_REF] (see [START_REF] Khalid | The Residential Quarter in Bukhara before the Revolution (The Work of O. A. Sukhareva)[END_REF] has done considerable work in compiling the oral and written testimonies on the guzars, providing extremely valuable information on the dating of even the smallest madrasas and mosques built in Bukhara over the past millennium.

Our archeomagnetic sampling was focused on several major and some minor buildings in and around Bukhara (Fig. 2, Table 1). Among the most important are three madrasas (with a group of fragments collected for each of them): two were built during the reign of Abdullah Khan (1583-1598), one of the most famous rulers of the Shaybanid dynasty. One was built for his own glory (Madrasa Abdullah Khan; BK03, Fig. 2a), and the second was to glorify his mother (Madrasa Modari Khan, built around 1561; BK04, Fig. 2b). The third madrasa was erected in ∼ 1651 -1652 by the Khan Abd al-Aziz of the Djanid dynasty (Madrasa Abd al-Aziz Khan; BK01). We also carried out a sampling in the Chor Bakr necropolis, built near Bukhara at the time of the Shaybanid dynasty, and at the location of older tombs dating from the 10th century. There, two groups of fragments associated with the tomb (or khazira) of Khwādja Saad, son of Khwādja Islām Juybārī, leader of the Sufi order (i.e., a mystical order of Islam), erected just before his death at the end of the 16th century, have been collected (Khwādja Saad tomb's wall and floor, BK05 and BK06 respectively; Fig. 2d; Table 1). The sampling also included the Ark citadel. This ancient fortress, which was last destroyed in 1920 and rebuilt several times during the history of the city, was inhabited by the rulers. It comprised several buildings surrounded by an imposing wall, among which we sampled the k ¯ānaqāh, which is a dwelling place for dervishes (adherents of Sufi orders) dating from the mid-18th century (BK14, Fig. 2f).

In addition, we sampled fragment groups from three minor mosques (Mosque Dostum Chor Oghasi, BK09; Mosque Magoki Kurpa, BK12; and Mosque Kemuhtagaron, BK13) and three madrasas from smaller neighborhoods (i.e.

the guzars; Madrasa Kunjak, BK07; Madrasa Rakhmanqul,BK08; Madrasa Rashid-al-Din, BK11, Table 1). It should also be mentioned that special care was taken to avoid restored wall segments and/or recycled bricks, which would result in an inaccurate dating.

[Figure 2 about here.]

In total, our archeomagnetic study is based on 13 groups of architectural brick fragments. The three above-mentioned dynasties are equally sampled :

five buildings belong to the Shaybanid dynasty, four buildings to the Djanid dynasty and four buildings to the Manghit dynasty. For each of the two fragment groups BK08 and BK11, the samples were collected in different rooms of the same building; in this case, each sub-subset has been identified but all fragments are considered to come from the same ensemble (e.g. BK08A

or BK08B; see supplementary Table S1). Particular care was taken with the available dating constraints and we selected those buildings that have age uncertainties of less than ±25 years, but in most cases these are less than 15 years. For some of these buildings, the construction is very well constrained by archives due to the social prominence of the people they were built for. In general, minor mosques and madrasas are first mentioned after construction in written documents, in particular for their inauguration or when a donor subsidized its use.

Sampling was carried out using an electrical driller with a water can.

From 10 to 18 cores, 2.5 cm in diameter and from 5 to 10 cm in length, were drilled per group (Fig. 3). A total number of 160 cores were analyzed in the present study.

[Figure 3 about here.]

3. Methods

Archeointensity determinations

All experiments were conducted in the paleomagnetic laboratory of the Institut de Physique du Globe de Paris (IPGP). The archeointensity determinations are based on the experimental protocol developed for the Triaxe magnetometer. This unique magnetometer allows continuous magnetization measurements (every ∼ 5 • C) of a small individual specimen (less than 1 cm 3 ) at high temperatures and under controlled field conditions, both in intensity and direction (Le Goff and Gallet, 2004).

The Triaxe procedure consists of five measurement series automatically performed between a low temperature T 1 , usually 150 • C, and a high temperature, T 2 , at which most of the magnetization carried by the specimen is erased:

-Step 1: After rapid heating from room temperature to T 1 , the specimen is heated in a zero field from T 1 to T 2 to demagnetize its natural remanent magnetization (NRM). The corresponding magnetization measurements give the M 1 series;

-Steps 2 and 3: The specimen is cooled from T 2 to T 1 (step 2) and next heated from T 1 to T 2 (step 3), both steps in a zero field, to allow characterization of the thermal variability of the NRM fraction that remains blocked at T 2 . The magnetization measurements give the M 2 and M 3 magnetization series, respectively;

-Step 4: The specimen is cooled from T 2 to T 1 in a laboratory field, the intensity of which is chosen close to the expected ancient field intensity, and its direction is automatically adjusted so that the direction of the newly acquired laboratory thermoremanent magnetization (TRM lab ) is parallel to the direction of the original TRM (i.e., NRM). This step therefore leads to the acquisition of a new TRM with unblocking temperatures between T 2 and T 1 (magnetization series M 4 );

-Step 5: The specimen is then heated again between T 1 and T 2 to demagnetize the TRM lab (magnetization series M 5 ).

The procedure ends with rapid cooling of the specimen to room temperature.

Intensity determinations are based on the ratio between the NRM and TRM lab fractions unblocked between T 1 and a temperature T i , varying from T 1 to T 2 . At any T i , these fractions are respectively determined by:

∆ 1 (T i ) = (M 1 (T 1 ) -M 1 (T i )) -(M 3 (T 1 ) -M 3 (T i )); (1) ∆ 5 (T i ) = (M 5 (T 1 ) -M 5 (T i )) -(M 3 (T 1 ) -M 3 (T i )). ( 2
)
And the intensity value at T i is given by:

R (T i ) = H lab ∆ 1 (T i ) ∆ 5 (T i ) . (3)
An intensity value is obtained for each specimen from the average of the R'(T i ) data derived from all temperatures T i between T 1 and T 2 (see more details and discussion in Le Goff and Gallet, 2004). Since intensity values should only be determined on the primary and single-vector magnetization acquired during the manufacture of the archeological artifacts, if a secondary magnetization is observed above T 1 up to T 1 ' (but below T 2 ) from the analysis of the NRM demagnetization data, then the reference temperature T 1 can be increased to T 1 '.

Compared to more conventional paleo/archeointensity methods that rely on stepwise demagnetizations and magnetization measurements carried out at room temperature, the Triaxe procedure has several advantages, including the fact that the TRM lab is acquired under thermal and field conditions relatively similar to those that led to the NRM acquisition. This helps to mitigate possible spurious effects that would result from the presence of multi-domain grains. The fact that the direction of the TRM lab is parallel to that of the original TRM eliminates the need for anisotropy correction on TRM acquisition (Le Goff and Gallet, 2004). In addition, experiments have shown that the use of R'(T i ) data allows to overcome the effect of the cooling rate on TRM acquisition (Le Goff and Gallet, 2004;Genevey et al., 2009;[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], see also a more general discussion on TRM anisotropy and cooling rate effects in Genevey et al., 2008).

The intensity data are then examined according to a set of quality criteria, which have remained the same as in previous studies in which archeointensity data obtained using the Triaxe procedure are reported (e.g. [START_REF] Gallet | Archaeomagnetism at Ebla (Tell Mardikh, Syria). New netic field intensity variations in Western Europe over the past 1100 years[END_REF][START_REF] Genevey | New archeointensity data from French Early Medieval pottery production (6th-10th century AD). Tracing 1500 years of geomagnetic field intensity variations in Western Europe[END_REF][START_REF] Genevey | Refining the high-fidelity archaeointensity curve for Western Europe over the past millennium: analysis of Tuscan architectural bricks (Italy)[END_REF]Gallet et al., 2014Gallet et al., , 2015Gallet et al., , 2020)). At the specimen level, the R'(T i ) data must involve at least 50% of the NRM still blocked at T 1 (or T 1 ') and the relative variations of R'(T i ) between T 1 (or T 1 ') and T 2 must not exceed 10%. A mean intensity value is determined at the fragment level when successful results are obtained from at least two specimens (note that in our study, we increase this number to three different specimens). In addition, the fragment-mean value is retained only if its standard deviation does not exceed 5% of the corresponding mean-intensity value. Finally, a mean intensity value is calculated at the level of a group of fragments when results are available from a minimum of three different fragments meeting the above criteria. The error is given as the standard deviation computed from the set of the retained intensity values at the fragment level. These criteria applied successively at the specimen, fragment, and fragment-group levels have proven to be extremely effective in constraining the quality and consistency of the intensity values obtained using the Triaxe method (e.g. [START_REF] Gallet | High-temperature archeointensity measurements from Mesopotamia[END_REF]Genevey et al., 2009;[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] Hervé | Fast geomagnetic field intensity variations between 1400 and 400 BCE: New archaeointensity data from Germany[END_REF].

Magnetic mineralogy characterizations

In addition to the archeointensity experiments, we also performed different analyses on the retained fragments to identify the magnetic minerals present in the samples and to further assess whether this magnetic mineralogy alters during heating. Analyses include, for all retained fragments, low-field susceptibility vs. temperature measurements (using a KLY3 kappabridge from Agico coupled with a CS3 furnace) and for at least two fragments from each retained group, the acquisition (using a Vibrating Sample Magnetometer Model 3900) of isothermal remanent magnetization (IRM)

and hysteresis loop measurements as well as the thermal demagnetization of three-axis IRM acquired (using a MMPM10 pulse magnetizer) in orthogonal fields of 1.5, 0.6 and 0.2 T [START_REF] Lowrie | Identification of ferromagnetic minerals in a rock by coercivity and unblocking temperature properties[END_REF]. Additionally, for a selection of representative samples, complementary low-temperature magnetization mea-surements are carried out using a magnetic property measurements system (MPMS XL-5 EverCool). The latter experiments include the following remanent magnetization measurements: 1) temperature cycling of an IRM acquired at room temperature (RT-SIRM) in a 2.5 T field down to 10 K and return to room temperature in a zero field (less than ± 500 nT), and 2) the thermal demagnetization from 10 K to 300 K of an IRM acquired at 10 K in a 2.5 T field following a zero-field cooling (ZFC) and 2.5 T-field cooling (FC)

pre-treatments from 300 K to 10 K. Both the RT-SIRM and ZFC-FC measurements were duplicated in a second series of experiments where the initial IRMs acquired in 2.5 T are demagnetized in a 300 mT using the MPMS's superconducting magnet in a field oscillation mode, a method introduced and validated in Lagroix and Guyodo (2017). The objective of the second series of experiment is to remove the contribution from low coercivity minerals to the total magnetization.

Archeointensity results

Magnetic mineralogy

IRM acquisition curves are reported in Fig. 4a. They show that saturation of the magnetization is often not completely achieved at 1T, but a clear inflexion in the magnetization curves is observed at ∼ 0.1 T. The thermal demagnetization of three-axis IRM further indicates that the magnetization is mostly carried by low-coercivity minerals with unblocking temperatures below 600 • C, which is consistent with the presence of (titano)magnetite (Fig. 4b-e). Fine grained hematite, lowering its unblocking temperature (e.g. [START_REF] Özdemir | Hysteresis and coercivity of hematite[END_REF], and/or epsilon iron oxide, a magnetic phase often observed in archeological artifacts (e.g. [START_REF] Genevey | New archeointensity data from French Early Medieval pottery production (6th-10th century AD). Tracing 1500 years of geomagnetic field intensity variations in Western Europe[END_REF][START_REF] López-Sánchez | Epsilon iron oxide: Origin of the high coercivity stable low Curie temperature magnetic phase found in heated archeological materials[END_REF][START_REF] Kostadinova-Avramova | Archaeomagnetic knowledge of Neolithic in Bulgaria with emphasis on intensity changes[END_REF] are mineral phases compatible with the above observations. A duality of low and high coercivity minerals is observed in a few hysteresis loops displaying slight constrictions (Fig. 5a). However, most of the hysteresis loops are not wasp-waisted, always exhibiting a monotonic decrease in loop opening with increasing field (Fig. 5b-c).

[Figure 4 Comparing the two provides information on the relative contribution of low or high coercivity minerals to the total remanence (Lagroix and Guyodo, 2017).

The lack of a Verwey transition in ZFC-FC data (Fig. 7) and primarily reversible RT-SIRM curves (Fig. 6a andc) or temperature suppressed Verwey transition (Fig. 6b) are compatible with titanomagnetite [START_REF] Kakol | Influence of titanium doping on the magnetocrystalline anisotropy of magnetite[END_REF][START_REF] Moskowitz | Low-temperature magnetic behavior of titanomagnetites[END_REF][START_REF] Muxworthy | Review of the low-temperature magnetic properties of magnetite from a rock magnetic perspective[END_REF]. Hematite is unambiguously identified in the demagnetized RT-SIRM data (right panels of Fig. 6a andc) from the observed Morin transition, which displays a remanence loss and partial recovery over a wide temperature range (∼ 235 K to 170 K) compatible with fine grained (0.1 to 1 µm) hematite [START_REF] Özdemir | Morin transition in hematite: Size dependence and thermal hysteresis[END_REF][START_REF] Özdemir | Hysteresis and coercivity of hematite[END_REF]. Another noteworthy observation is the kink seen at ∼ 70 K in both ZFC and FC curves which is also compatible with Ti-rich (50 to 60 % Ti; see [START_REF] Moskowitz | Low-temperature magnetic behavior of titanomagnetites[END_REF] titanomagnetite and the persistence of the kink after 300 mT demagnetization (except for BK04-10) finds an explanation in the significant increase in coercivity at low temperature of Ti-rich titanomagnetite (see for example fig. 15c in [START_REF] Almeida | Magnetic characterization of synthetic titanomagnetites: Quantifying the recording fidelity of ideal synthetic analogs[END_REF]. The 70 K kink could alternatively be related to the epsilon iron oxide phase (López-Sánchez et al., 2016, 2017). Lastly, behaviour suggestive of nanogoethite [START_REF] Guyodo | From nanodots to nanorods: Oriented aggregation and magnetic evolution of nanocrystalline goethite[END_REF], which would be of weathering origin, is occasionally observed (right panels of Fig. 6b and Fig. 7a).

[Figure 6 about here.]

[Figure 7 about here.]

New archeointensity data

We analyzed a total of 160 fragments (532 specimens) from 13 different archeological (historical) contexts. Most often, the magnetization of the samples is comprised between ∼ 30 and ∼ 140 × 10 -8 Am 2 (with a maximum of ∼ 500 × 10 -8 Am 2 ) and 42 of them are too weak to be measured with the Triaxe magnetometer (< 30 × 10 -8 Am 2 ), which has a measurement sensitivity on the order of ∼ 10 -8 Am 2 (Le Goff and Gallet, 2004). Among the 118 remaining fragments, 70 fragments are rejected due to non-linear or complex behavior compared to the nominal behavior described in Le Goff and Gallet (2004) (see also Genevey et al., 2009;[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] and because of scattered magnetization measurements. with 21 specimens; Table 1). Fig. 9 also illustrates the overall consistency of the data obtained for each group of fragments, resulting in small standard deviations. They range from 0.9 µT (BK08) to 2.0 µT (BK11), or between 1.7% (BK06) and 4.5% (BK11) of the mean intensity values.

[Figure 9 about here.]

The new archeointensity data obtained at Bukhara cover a time interval of ∼ 250 years, between the mid-16th century and ∼ 1800 (Table 1; Fig. 10).

A significant decrease in intensity values by ∼ 14 µT is observed from ∼ 1560 to ∼ 1725, leading to an average rate of change of ∼ -0.1 µT/yr. The second half of the 18th century is then marked by an increase of ∼ 4 µT until the early 19th, leading to a variation rate of ∼ 0.05 µT/yr. This rate of change is fairly comparable to that of the present-day field in Bukhara (∼ 0.08 µT/yr).

[Table 1 about here.]

Discussion

Comparison of the new archeointensity data with model predictions at

Bukhara

The new archeointensity data obtained in Bukhara are compared in Fig. 10 with the variations in intensities predicted by several geomagnetic field reconstructions (see Section 1). For the historical period, this is the gufm1 model [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF], and the models of Unlike the other models, BIGMUDI4k.1 is built from the simultaneous inversion of both direct and indirect data. The authors note a significant decrease in the dipole energy associated with an increase of the non-dipole energy around ∼ 1600. According to [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], this is mainly due to the large increase in the amount of data at the onset of the historical geomagnetic era, rather than a true geomagnetic feature. The subsequent increase in dipole energy is therefore probably artificial, as is the resulting intensity peak. The modeled field behavior during this period should therefore be considered with caution. The COV-ARCH model, integrating only archeomagnetic data, predicts an evolution close to those of [START_REF] Gubbins | Fall in Earth's magnetic field is erratic[END_REF] and [START_REF] Finlay | Historical variation of the geomagnetic axial dipole[END_REF]. This evolution shows a minimum intensity at the end of the 18th century, slightly later than the minimum intensity observed from our data. Interestingly, in this model, the high-frequency range of the axial dipole variations is constrained by a timescale on the order of the convective turnover time (recall Introduction).

However, regardless of the model, none of the expected intensity evolutions reach the low intensity values observed in the 18th century from the new Bukhara archeointensity data, with a minimum overestimate of ∼ 5 -6 µT.

Prior to this period, the Bukhara data show a rapid decrease in intensities between 1550 and the early 18th century with a rate quite similar to that of gufm1 (see also the model of [START_REF] Suttie | Consistent treatment of errors in archaeointensity implies rapid decay of the dipole prior to 1840[END_REF]. On the other hand, the new data require an increase in intensities during the first half of the 19th century, which is either absent or much more limited in model predictions.

Dispersion of archeointensity results in regional datasets

[Figure 10 about here.]

Comparison with other archeointensity results previously obtained in the Bukhara area, as well as elsewhere in western Eurasia, raises a critical problem related to the dispersion that generally characterizes the regional datasets. Here we distinguish four geographical areas within a 700 km-radius around the cities of Bukhara (Uzbekistan), Moscow (western Russia), Tbilisi (Georgia), and Thessaloniki (Balkans). In each zone, the data are reduced to the latitude of the corresponding city. Most of these data are rather old and were described in the ArcheoInt compilation (Genevey et al., 2008). In our study, they are selected using at first the same minimalist criteria as in Genevey et al. (2008) (referred to as G2008 below). A location map of the selected data is given in supplementary material (Fig. S1). These criteria were originally proposed to allow the discussion of old data acquired without all the quality criteria now considered necessary for any new study, whereas applying modern criteria would eliminate most (if not all) of them. They do not consider the intensity methods directly but instead rely on 1) the error (most often a standard deviation) on the average intensity, which must be known and less or equal to 15%; 2) the number of intensity determinations (Nint) used to derive an intensity mean. Nint is required to be greater than or equal to three when no pTRM-check was implemented or when this test does not apply. Otherwise, Nint must be greater than or equal to two. For objects recognized to be strongly anisotropic (such as pottery or tiles), Nint is required to be greater than or equal to three if anisotropy effects on TRM acquisition were not taken into account. In a second step, we consider stricter criteria requiring pTRM-check (when this test does apply) and intensity average derived from at least three independent fragments. In both cases, we select data whose age uncertainties are less or equal to ±50 years because we are interested in fairly rapid variations over a short time interval of ∼ 300 years. Further note that practically none of the available data have been corrected for the cooling rate effect on TRM acquisition.

Intensive work by Russian archeomagnetists in the 1970s and 1980s (S.

Burlatskaya, I. Nachasova and K. Burakov) resulted in three regional datasets, in Uzbekistan, around Moscow (Moscow, Gor'kiy, and Vologda) and in Georgia. These data share common features. They were acquired from analyzed baked bricks. The number of intensity determinations corresponds to the number of independent bricks studied. The method used is the original Thellier and Thellier (1959) protocol with the use of pTRM-check along the measurements cycle. This key element was, however, not specified in the articles (mainly published in Russian journals) but given by S. Burlatskaya, I. Nachasova and K. Burakov in a personal communication (2004). A single dataset was obtained both using the Thellier and Thellier (1959) protocol and an original method developed by [START_REF] Burakov | A method and results of studying the geomagnetic field of Khiva from the middle of the sixteenth century[END_REF], so-called thermal curves, derived from Wilson (1961)'s method (the TRM anisotropy being also taken into account). Both datasets are considered in our paper (identified by two different symbols), although we should point out that [START_REF] Nachasova | Geomagnetic field variations in Central Asia during the last 2000 years: An analysis of global data[END_REF] argue that the thermal-curves dataset is more reliable.

The results from Uzbekistan, more precisely obtained in the cities of Bukhara [START_REF] Burlatskaya | New determinations of the parameters of the ancient geomagnetic field for Mongolia, Soviet Central Asia, and Abkhazia. Geomagnetism and Aeronomy[END_REF](Burlatskaya et al., , 1986b)), Samarkand [START_REF] Burlatskaya | Some archaeomagnetic data indicative of the westward drift of the geomagnetic field[END_REF](Burlatskaya et al., , 1986b)), and Khiva [START_REF] Burakov | A method and results of studying the geomagnetic field of Khiva from the middle of the sixteenth century[END_REF] are of particular interest as they allow a direct comparison with the new archeointensity data reported in the present study (Fig. 11a). For Bukhara, the sampled sites are unfortunately not indicated in the original Russian publication (nor in the compilation of Burlatskaya et al., 1986b), but their estimation of the age of the buildings does not correlate with our new data, and the buildings sampled in the present study are therefore probably different. Despite some scatter, a fairly satisfactory agreement could be found for all results dating from ∼ 1700 to ∼ 1850, with values often lower than the intensities expected in Bukhara from the models. However, this satisfaction must be tempered by the fact that the older data (before ∼ 1700) appear systematically weaker than our own intensity values. The discrepancy would be even larger if a cooling rate correction (for instance 5% as suggested by Genevey et al., 2008) was applied to the data. This ambiguous information is not improved by considering stricter selection criteria as most of the Uzbek data meet those criteria (see Fig. S2 in supplementary data).

[Figure 11 about here.]

A fairly large scatter is observed for the data around Moscow, whether these data are selected using the G2008 set of criteria or stricter criteria (remember that only the number of fragments per site, Nint≥ 2 or ≥ 3, is the difference between the two selections, Fig. 11b and Fig. S2b). This dispersion questions the reliability of at least part of these data, as previously discussed in Salnaia et al. (2017a,b). The data seem consistent with a decreasing trend in intensities over the historical period. Nonetheless, as the cooling rate effect was not evaluated in the old Russian studies [START_REF] Nachasova | Magnetic field in the Moscow area from 1480 to 1840. Geomagnetism and Aeronomy[END_REF]Burlatskaya et al., 1986b), these data could also agree with lower than predicted values during the 18th century. Comparing the scant data obtained in Georgia with the model predictions leads to another contrast (see Fig. 11c). While fairly consistent agreement is observed for the data up to ∼ 1700, the results are more scattered from ∼ 1700 to ∼ 1850, and in general the results are less consistent with the expected intensity values. With an arbitrary correction of the cooling rate effect of 5%, some of the latter results would be too low, in particular those dating from ∼ 1800 to ∼ 1850.

Finally, the Balkan area incorporates results from Greece and Bulgaria (Aitken et al., 1989a;[START_REF] Spassov | Rock magnetic property and paleointensity determination on historical Santorini lava flows[END_REF][START_REF] Kovacheva | Updated archeomagnetic data set of the past 8 millennia from the Sofia laboratory, Bulgaria[END_REF]Kovacheva et al., , 2014)).

Contrary to the Russian datasets, whether pTRM checks are implemented or not is critical for distinguishing between the two selected datasets based on the G2008 versus the stricter set of criteria. (Fig. 11d and Fig. S2d). With the loose selection (Fig. 11d), the data available between ∼ 1550 and ∼ 1700 appear relatively scattered, with the Bulgarian data in particular generally higher than the values expected from the models. A limited decrease to account for the cooling rate effect would improve the agreement. However, when stricter criteria are applied, all but one data point from this time interval are eliminated (Fig. S2d). Still considering this selection, five data points remain for the 18th century (including two obtained by [START_REF] Spassov | Rock magnetic property and paleointensity determination on historical Santorini lava flows[END_REF] from the same volcanic event) and only one for the 19th century. We may note that two Bulgarian results from the second half of the 18th century do not seem to indicate lower values than those predicted by the models. In contrast, this would be the case for the only result dating from the early 19th century. At this stage, it is therefore difficult to draw a firm conclusion from these rare data as well as from the entire Balkan dataset.

Overall, Fig. 11 shows that the data available in each of the four areas discussed above are too scattered to show any consistent pattern of intensity variations, at least at the century scale. As pointed out by many authors, determining a set of selection criteria that allows for significant reduction in the regional data scatter is a challenge. In the present case, increasing the strictness of the criteria does not alter our conclusion (at most, it leads to the rejection of most data in the Balkans), as this was also previously observed and discussed for western Europe (Genevey et al., 2009[START_REF] Gallet | Archaeomagnetism at Ebla (Tell Mardikh, Syria). New netic field intensity variations in Western Europe over the past 1100 years[END_REF][START_REF] Genevey | Refining the high-fidelity archaeointensity curve for Western Europe over the past millennium: analysis of Tuscan architectural bricks (Italy)[END_REF].

Our purpose is not to analyze and discuss in detail all datasets currently available worldwide (see for instance [START_REF] Poletti | Continuous millennial decrease of the Earth's magnetic axial dipole[END_REF]. As observed in western Eurasia, data at the regional spatial scale are generally either too scant or too scattered to draw a clear evolution of intensities, which could lead to a meaningful comparison with the intensity values derived from the geomagnetic field models (see for instance Tema et al., 2017;[START_REF] Goguitchaichvili | Last three millennia Earth's magnetic field strength in Mesoamerica and southern United States: Implications in geomagnetism and archaeology[END_REF]Kapper et al., 2020, for Hawaii, Mesoamerica and West Africa, respectively). Due to the general dispersion of the archeointensity data at the regional scale, a linear evolution of the dipole moment as constrained over the entire historical period, notably between 1600 and 1800, by a recalibration of g 0 1 with these data, is a simple and reasonable approximation (e.g. [START_REF] Gubbins | Fall in Earth's magnetic field is erratic[END_REF]. However, this does not demonstrate that the axial dipole moment evolution is actually linear; the scatter and the small amount of data leave room for more complex, possibly hidden variations in axial dipole moment.

A non-linear evolution of the axial dipole moment over the historical period

One might consider two options for explaining the dispersion of the data, either the frequent presence of biased results masking the "true" regional field intensity evolution or an inherent limitation in archeointensity determinations. In other words, due to their lack of resolution (and/or underestimation of their uncertainties), archeointensity data could not reliably detect and describe century-scale intensity variations. It is worth pointing out that this (dull) option is in clear contradiction with the convincing detection in western Europe of century-scale intensity variations over the past ∼ 1500 years (Genevey et al., 2009[START_REF] Gallet | Archaeomagnetism at Ebla (Tell Mardikh, Syria). New netic field intensity variations in Western Europe over the past 1100 years[END_REF][START_REF] Genevey | New archeointensity data from French Early Medieval pottery production (6th-10th century AD). Tracing 1500 years of geomagnetic field intensity variations in Western Europe[END_REF][START_REF] Genevey | Refining the high-fidelity archaeointensity curve for Western Europe over the past millennium: analysis of Tuscan architectural bricks (Italy)[END_REF]. For the historical period, western Europe benefits from a fairly dense archeointensity dataset showing a smooth evolution, with reduced dispersion (see description of the data in the mentioned studies). This leads Genevey et al. (2009) to explore a different approach for the recalibration of gufm1 Gauss coefficients by using a limited but consistent regional dataset.

[Figure 12 about here.] Following [START_REF] Gubbins | Fall in Earth's magnetic field is erratic[END_REF] and Genevey et al. (2009), a ratio of measured to predicted intensity is determined for each data point of Genevey et al. (2009[START_REF] Gallet | Archaeomagnetism at Ebla (Tell Mardikh, Syria). New netic field intensity variations in Western Europe over the past 1100 years[END_REF][START_REF] Genevey | Refining the high-fidelity archaeointensity curve for Western Europe over the past millennium: analysis of Tuscan architectural bricks (Italy)[END_REF]. It is then used to recalibrate all Gauss coefficients from gufm1, allowing the computation of a new field intensity prediction at Bukhara. This procedure applied to all the western European Triaxe data, to which we add three Triaxe data obtained in Russia by Salnaia et al. (2017a,b), allows for the determination of a consistent dataset (Fig. 12a). The results

show a clear intensity decrease between ∼ 1600 and the first half of the 18th century, followed by an increase up to 1850. This intensity pattern is further evidenced by the computation of a mean intensity variation curve and its credible interval using the transdimensional Bayesian method recently developed by Livermore et al. (2018). Based on this consistency, we also use the same dataset to recalibrate the axial dipole component (g 0 1 ) given by gufm1 (Fig. 12b). While the g 0 1 values (provided in Table S2) are rather compatible with those of the models during most of the 17th century, significant differences are then observed with all models, with smaller recalibrated values, throughout the entire 18th century and the early 19th century. On the other hand, our study also shows that there is currently no dense regional archeointensity dataset in western Eurasia, as elsewhere considering also the dispersion of the data, that could clearly contradict this g 0 1 evolution. As previously suggested by Genevey et al. (2009) and now based on a larger collection of results, Fig. 12b strongly militates for a non linear evolution of the axial dipole field moment over the historical period, with a distinct minimum of |g 0 1 (t)| ≈ 29400 nT during the 18th century. The average rate of decrease of |g 0 1 (t)| during the 17th and the late 18th reaches ∼ -26 nT/yr, while the increase during the first half of the 19th reaches a rate of ∼ 34 nT/yr. These two variation rates are higher than the one observed over the past 150 years (∼ 15 nT/yr; e.g. [START_REF] Barraclough | Spherical harmonic analyses of the geomagnetic field for eight epochs between 1600 and 1910[END_REF][START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF].

The regional approach used above is based on the reliability and accuracy of the geomagnetic field geometry of the gufm1 model [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF]. The implication on the axial dipole field moment's evolution between 1590 and 1850, as mainly constrained by the available Triaxe archeointensity data, therefore depends on this reliability. However, several studies have highlighted discrepancies between archeomagnetic directional data and the directions predicted by the gufm1 model (see for instance Tanguy et al., 2011, for the western Indian Ocean). In France, Le Goff and Gallet (2017) have also shown that while satisfactory consistency is observed after ∼ 1675, the gufm1 predictions differ significantly from most direct directional measurements prior to this date.

The relatively low reliability of the gufm1 model during the 17th century and part of the 18th century should not be surprising, given that very few, if any direct inclination data are available before 1700 -1750, and more generally, given the poor spatial and temporal coverage of historical directional measurements between 1590 and ∼ 1700 (e.g. fig. 1a,e [START_REF] Jackson | Four centuries of geomagnetic secular variation from historical records[END_REF], in particular in Central Asia (see fig. 8 to 14 in [START_REF] Jonkers | Four centuries of geomagnetic data from historical records[END_REF]. This calls for caution when interpreting recalibrated variations of axial dipole moments as shown in Fig. 12b. However, at the scale of Western Eurasia (from Western Europe to Central Asia), the satisfactory modelling of the non-dipole effects in gufm1 is evidenced by the good consistency, upon recalibration, of the Triaxe archeointensity datasets shown in Fig. 12a. In addition, it should be noted that this consistency is poorer when the data from Western Europe and Russia are transferred to the latitude of Bukhara using a purely axial dipole field approximation (Fig. S3). Nevertheless, the clear non-linear dipole moment evolution deduced using the gufm1 model in a region across which significant non-dipole field effects are not expected to occur over the historical period (e.g. Pavón-Carrasco et al., 2014b) does not demonstrate its truly dipole origin. For this, we need a large set of reliable and geographically distributed archeointensity data from around the world.

Away from western Eurasia, [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] This feature could thus represent a true dipole feature contradicting a linear evolution of the g 0 1 term over the entire historical period. At this stage, however, we recognize that its global (dipole) nature has yet to be confirmed by the acquisition of new high quality archeointensity data.

Conclusions

The acquisition of nine new archeointensity data from Bukhara, Uzbek- Figure 8: Left panels: thermal demagnetization data obtained for three different specimens. Open (close) symbols refer to the inclinations (declinations). Right panels: corresponding R'(T i ) datasets obtained from the same specimens (see in Le Goff and Gallet, 2004)

  Fig. 4b-e also shows the presence of high-coercivity minerals whose unblocking temperatures do not exceed ∼ 550 • C, being sometimes as low as ∼ 200 • C or with an inflexion around this temperature (Fig 4b-d).

  about here.] Low-field susceptibility versus temperature (heating and cooling) curves yield two main observations (Fig.5d-i). First, heating and cooling curves are reversible or very nearly, which attests to the stability on heating of the magnetic mineralogy in the temperature range used for intensity determinations. Second, most susceptibility curves show a clear inflexion around 300 • C, arising from a range of susceptibility evolutions from rapid rates of change (Figs.5g,i) to slower monotonic ones (Figs.5d,e,f,h), in addition to a higher temperature inflexion above 500 • C. At this stage, we could propose that the inflexions indicate the presence of two families of (titano)magnetite differing by their grain size, their titanium content and/or their oxidation state (see below).[Figure5about here.] Low-temperature magnetization measurements bring additional insight into the magnetic mineralogy. Compared to the classic RT-SIRM and ZFC-FC experiments (left panels in Fig. 6 and 7 respectively), their 300 mT demagnetized counterparts highlight the temperature dependent behaviour of the high coervicity minerals (right panels in Fig. 6 and 7 respectively).

Fig. 8 (

 8 Fig. 8 (left panels), together with the corresponding R'(T i ) data (right panels). In general, the specimens are fully demagnetized at relatively low temperatures, below 450-500 • C. A single magnetization component is essentially isolated, even though a small secondary component probably of viscous origin is identified in most cases at low temperatures (below 150 • C) but also in some cases having slightly higher temperatures (around 200 • C).

  [START_REF] Gubbins | Fall in Earth's magnetic field is erratic[END_REF];[START_REF] Finlay | Historical variation of the geomagnetic axial dipole[END_REF];[START_REF] Suttie | Consistent treatment of errors in archaeointensity implies rapid decay of the dipole prior to 1840[END_REF] derived from gufm1 and calibrated for the 1590 -1840 time interval from a global compilation of archeointensity measurements. Over this time interval, they predict the same pattern of variation but with various amplitudes corresponding to the different rates of decay imposed on the axial dipole component (recall the related commentary in the introduction). The comparison is also extended to geomagnetic models covering longer time intervals (between 3000 and 14, 000 years): A FM (Licht et al., 2013), pfm9k.1[START_REF] Nilsson | Reconstructing Holocene geomagnetic field variation: new methods, models and implications[END_REF], SHA.DIF.14k (Pavón-Carrasco et al., 2014a), CALS10k.2 and ARCH10k.1b[START_REF] Constable | Persistent high paleosecular variation activity in southern hemisphere for at least 10 000 years[END_REF], COV-ARCH[START_REF] Hellio | Time-correlation-based regression of the geomagnetic field from archeological and sediment records[END_REF] and BIGMUDI4k.1[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] constructed using global archeomagnetic datasets. A number of these models are constrained by gufm1 over the historical period (i.e. CALS10k.2, ARCH10k.1b, SHA.DIF.14k) and the corresponding predictions fall within the range of gufm1-recalibrated models. The pfm9k.1 model predicts a similar evolution, although the predicted dipole moment is higher compared to the other models for this period. The authors interpret this overestimation as resulting from the introduction of directional sedimentary data[START_REF] Nilsson | Reconstructing Holocene geomagnetic field variation: new methods, models and implications[END_REF]. As this study is focused on the field variations over the his-torical period, the models constrained by gufm1 (CALS10k.2, ARCH10k.1b, SHA.DIF.14k) and those integrating sedimentary data (pfm9k.1) are not represented. A FM predicts an intensity evolution in Bukhara very close to the prediction from gufm1, with higher values than the observed intensities during the 18th and early 19th century. On the other hand, BIGMUDI4k.1 predicts a different intensity evolution, with a quasi-constant intensity during the 17th century and a well-marked intensity peak during the 18th century.

Fig. 12a shows

  Fig. 12a shows a direct comparison between the new data from Bukhara and the western European results recently upgraded and summarized in Genevey et al. (2019). Most of these results share the same (Triaxe) experimental methodology and obey the same set of selection criteria. This comparison takes into account the geomagnetic field geometry given by the gufm1 model.

  obtained in southern and northern Brazil coherent historical archeointensity data using the Triaxe protocol (so far this is the only Triaxe data obtained outside western Eurasia) thus sharing the same criteria as before. As pointed out by Hartmann et al. (2011), these results also show the gufm1 model's lack of reliability for the period before ∼ 1750. Note that based on an archeointensity result obtained in Ethiopia dated ∼ 1615, Osete et al. (2015) also arrive at the same conclusion. On the other hand, using the same recalibration method as previously used byGenevey et al. (2009), the Brazilian data are consistent with a minimum of the axial dipole magnitude |g 0 1 (t)| around the late 18th century, as shown by the western Eurasian Triaxe data (Fig.12).

  figure). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 11 Archeointensity results obtained in a 700-km radius from a) Bukhara, b) Moscow (Russia), c) Tbilisi (Georgia), d) Thessaloniki (Greece), reduced at the latitude of the corresponding location. The data are filtered using the G2008 set of criteria. Each panel also shows the predicted intensity evolution from various geomagnetic models at the corresponding location (continuous lines, errors are given as two standard deviations by shaded areas; see legend and text for details). . . . . 61

Figure 1 :Figure 2 :Figure 3 :Figure 4 Figure 5 :

 12345 Figure 1: General location map of Uzbekistan and Bukhara

Figure 6 :

 6 Figure6: Representative examples of RT-SIRM cycles of a 2.5 T field IRM (left panels) and a 2.5 T IRM partially demagnetized (right panels) with a 300 mT field generated by the superconducting magnet operating in oscillation mode. The blue (resp. orange) dots correspond to the cooling (resp. heating) step. The results are normalized to M n (corresponding to the initial RT-SIRM at 300 K).

Figure 7 :

 7 Figure7: ZFC-FC warming curves for the same fragments as in Fig.6of 2.5 T IRMs acquired at 10 K (left panels) and 2.5 T IRMs partially demagnetized at 10 K (right panels) with a 300 mT field generated by the superconducting magnet operating in oscillation mode. The blue (resp. orange) dots correspond to the ZFC (resp. FC) step. The results are normalized to M n (corresponding to the initial FC at 10 K).

Acknowledgements

We are grateful to Sirozh Mirzaakhmedov for his help during the sampling in Bukhara. We would like to thank Maxime Le Goff for insightful discussions and his constant and efficient technical assistance in the use of the Triaxe magnetometers. We also thank Nicolas Gillet for providing programs to compute the COV-ARCH predictions. We are grateful to Patrick Arneitz and one anonymous reviewer for their helpful comments on the manuscript. This study was financed by the Simone and Cino Del Duca Foundation of the French Academy of Science and by the INSU-CNRS program PNP. VP aknowledges support from the Russian foundation for basic research RFBR # 18-55-41005. VP was also partly supported by the Ministry of Science and High Education of the Russian Federation (grant no. 14.Y26.31.0029 in the framework of Resolution no.220 of the government of the Russian Federation). MPMS and VSM measurements were carried out on the IPGP-IMPMC Mineral Magnetism Analytical Platform financed by Region Île-de-France, IPGP, UPMC, CNRS-INSU and ANR. This is IPGP contribution no. xxxx. ular variation curves from Bulgaria for the last eight millennia. Physics of the Earth and Planetary Interiors, 236:79-94, 2014. F. Lagroix and Y. Guyodo. A new tool for separating the magnetic mineralogy of complex mineral assemblages from low temperature magnetic behavior. Frontiers in Earth Science, 5:61, 2017. M. Le Goff and Y. Gallet. A new three-axis vibrating sample magnetometer for continuous high-temperature magnetization measurements: applications to paleo-and archeo-intensity determinations. Earth and Planetary Science Letters, 229(1-2):31-43, 2004. M. Le Goff and Y. Gallet. A reappraisal of instrumental magnetic measurements made in Western Europe before AD 1750: confronting historical geomagnetism and archeomagnetism. Earth, Planets and Space, 69(1): 1-9, 2017. V. Lesur, I. Wardinski, J. Baerenzung, and M. Holschneider. On the frequency spectra of the core magnetic field Gauss coefficients. Physics of the Earth and Planetary Interiors, 276:145-158, 2018. A. Licht, G. Hulot, Y. Gallet, and E. Thébault. Ensembles of low degree archeomagnetic field models for the past three millennia. Physics of the Earth and Planetary Interiors, 224:38-67, 2013. P. W. Livermore, A. Fournier, Y. Gallet, and T. Bodin. Transdimensional inference of archeomagnetic intensity change. Geophysical Journal International, 215(3):2008-2034, 2018.

Fournaise: bearing on volcanic activity at La Réunion Island and geomagnetic secular variation in Southern Indian Ocean. Earth and Planetary Science Letters, 53.6 Figure 12: a) New intensity evolution in Bukhara predicted by gufm1 recalibrated with the Triaxe data from western Europe (blue dots, Genevey et al. (2009[START_REF] Gallet | Archaeomagnetism at Ebla (Tell Mardikh, Syria). New netic field intensity variations in Western Europe over the past 1100 years[END_REF][START_REF] Genevey | Refining the high-fidelity archaeointensity curve for Western Europe over the past millennium: analysis of Tuscan architectural bricks (Italy)[END_REF]) and Russia (grey dots, (Salnaia et al., 2017a,b)), with the mean intensity variations curve and its 95% credible interval (in blue). This curve is computed using the AH-RJMCMC algorithm from Livermore et al. ( 2018) using the following input parameters: σ move = 30 yrs, σ change = 5 yrs, σ birth = 5 yrs, K max = 150. The intensity priors are set to a minimum of 35 µT and a maximum of 60 µT, with a chain length of 100 million samples (see Livermore et al., 2018, for details on the parameters). To stabilize the prediction for the younger period, the prediction is tied to the intensity value predicted by gufm1 in 1860 (47.5 µT). b) Evolution of the axial dipole component g 0 1 over the past four centuries. Dots gives the recalibration of g 0 1 from gufm1 by the new archeointensity data from Bukhara (red dots), western Europe and the Russian datasets (blue and grey dots resp.), with the median variations curve and its 95% credible interval computed using the same parameters described above (except for the intensity priors set to -38 µT and -26 µT). The continuous lines give g 0 1 as provided by gufm1 and various derived models (see text for details and Table S2 for