Archaeomagnetic investigations in Bolgar (Tatarstan)

The objective of this study is to provide a well-dated point for a future palaeosecular variation (PSV) reference curve for western Russia. For this purpose archaeomagnetic and magnetic property analyses were carried out on a pottery kiln unearthed at the UNESCO World Heritage site of ancient Bolgar, having a rather precise age dating. The archaeological context provided an age between 1340 and 1360 C.E. The characteristic remanence vector was determined through alternating field demagnetisation and Thellier-Thellier palaeointensity experiments. Some innovations were introduced regarding palaeointensity. The check testing the equality of blocking and unblocking temperature was redefined. This allowed waiving the commonly used additional zero-field cooling steps during the Thellier-Thellier experiment. Another innovation concerns the calculation of archaeointensity at structure level. A Bayesian approach was introduced for averaging individual specimen archaeointensities using a prior probability distribution of unknown uncertainties. Next, an additional prior probability distribution was used to correct for cooling rate effects. This resulted in a lower uncertainty compared to common practice and in eluding time consuming cooling rate experiments. The complex magnetic mineralogy consists of maghaemite, multi-domain haematite and Al-substituted haematite. Some samples contained also some non-stoichiometric magnetite. The magnetic mineralogy was determined through hysteresis loops, backfield and remanence decay curves, measurements of the frequency dependence of magnetic susceptibility and through low temperature magnetisation curves. Accompanying high-temperature thermomagnetic analyses revealed an excellent thermo-chemical stability of the studied specimens. Directions obtained from alternating field demagnetisation and those extracted from archaeointensity experiments are congruent and have low uncertainties. The obtained archaeomagnetic results are fairly in agreement with global geomagnetic field models and contemporary PSV data of the wider area. The geomagnetic field vector obtained for ancient Bolgar is of high quality, deserving thus its inclusion in a future PSV reference curve for European Russia.

Heritage sites. The complex represents traces of an ancient city, dating back to the beginning of the Xth century and which is associated with the establishment of the Volga-Bulgarian state, according to historiography. The culmination of this development reaches the city of ancient Bolgar in the first half of the XIVth century, i.e. at the period when the territory was included in state of the Golden Horde. From second half of the ХIIIth to the beginning of the XVth century was Bolgar the capital of the Volga-Bulgarian ulus of the Golden Horde. This period is characterised by erected monumental buildings of white stones and bricks such as for instance mosques, baths, mausoleums, bazaar. Bolgar is considered as one of the key places in international trade along the Volga waterway during the Middle Ages. The geographical spread of imported ceramic ware confirms wide international ties, reaching from China and Central Asia to Spain. Various crafts underwent significant development on the city territory. Metallurgical, pottery, jewellery, bone-carving, glass-processing and other workshops were discovered in the Volga-Bulgarian settlement through systematic archaeological research in the complex and around, since 1938. Perennial research allowed to develop a stratigraphy of the occupation layer of Bolgar. It includes 7 layers: starting from the 7th layer, i.e. the pre-urban settlements of Imenkov culture and the early Bulgarian layer 6 indicating the stage of birth of Bolgar in the Xth century until layer 2 containing deposits related to the Russian village of the XVIIth  XIXth centuries and layer 1 -the modern village. Most significant and widespread indications of settlement were found in layer 4, i.e. the layer of the city's Golden Horde epoch. Following A.P. Smirnov's expeditions, it is divided into two horizons -an early Golden Horde and the later Golden Horde (Khlebnikova, 1987).
The object of this study are two pottery kilns from the territory of the Volga-Bulgarian settlement. The kilns were identified and investigated in the course of archaeological work in 1980 and 1982. Both were located in the outskirts of the South-Western town of the Golden Horde period and at the Northern edge of the seasonal reservoir which is called Lake Gollanskoe. The pottery kilns were part of a workshop for burning of ceramic products such as the large pots and one-handed jugs (page 134 in Vasilieva, 1988). The complex itself consisted of a near-kiln pit and five kilns around. All fire kiln facilities are of the same type -"two-tiered, recurring, with upward movement of gases" (Fig. 79 in Vasilieva, 1988). Most likely the excavated kilns did not function at the same time, but were successively replaced by each other for several decades. After the research in 1982, the kiln complex was preserved. Re-opening of the kilns for study purposes and subsequent museum conservation was carried out in 2012. The samples concerned in the present study were collected at that moment.
Initially, the researchers of the Bulgarian settlement believed that the pottery workshops in the area around ancient Bolgar is limited to 20 to early 60-ies of the XIVth century (Khlebnikova, 1987) based on at the coin and other findings. Further investigations narrowed the operating age of the kiln complex between the forties and sixties of the XIVth century (Vasilieva, 1988). Two silver coins were found in the layer adjoining the Kiln 1 (Table 1). According to the coin typology developed by S.А. Yanina, they are attributed to the end of the XIIIth century (Yanina, 1954) and would suggest an earlier time of settlement development in the area surrounding the kilns. However, the pottery complex itself cuts the layer in which these two silver coins were found. The kilns are consequently somewhat younger than the coins, i.e. later than the end of the XIIIth century. The kilns were definitively not anymore in use after the sixties of the XIVth century when the city was captured by one of the Horde princes Bulak-Timur. Perennial archaeological research evidenced that settlement activities in the area surrounding the city had completely stopped and were not restored after the devastating consequences of the events of 1361 (Khlebnikova, 1987). Concerning the kiln site of the present study, this is evidenced by the absence of any numismatic artefacts of the late XIVth to early XVth centuries. The latest numismatic finding which originates from the layer overlying the kilns, is the Hyzr Khan copper coin of Gulistan's 13611362 coinage (Fedorov-Davydov, 1987). Finally, the study of Vasilieva (1988) confined the period of kiln use between 1340 and 1360 C.E.
The kilns from which samples were selected, are situated in the southwestern part of the production complex. The two kilns, Kiln 1 and Kiln 2, which were presumed to be well-preserved, were re-discovered during a high-precision magnetic survey carried out in the ancient town of Bolgar (54.98°N, 49.05°E) in 2011. Unfortunately it turned out that only Kiln 1 was well-preserved and remained in-situ. Research proceeded in 2012 with archaeomagnetic investigations of several samples taken from pottery Kiln 1 ( Fig. 1) consisting of baked clay. Two independently oriented block samples A and B (~15  20  8 cm) were taken from the floor of Kiln 1 (Fig. 2). The sample's azimuth was determined with a precision compass which was placed 5 to 6 m away from the kiln. At this distance the compass needle was not influenced by the magnetic anomaly of the archaeological structure as revealed by the preliminary magnetic survey. The magnetic azimuth of the samples was traced by spanning a thread between compass needle and sample; the value was corrected for the local magnetic declination, i.e. 12.2° in Bolgar (http://www.geomag.bgs.ac.uk/cgi-bin/igrfsynth). The surface of the sample was scraped under bubble level control until a horizontal surface was obtained.
In the laboratory, both block samples A and B were cut into cubic specimens of 2 cm edge length using a diamond blade saw. In total, 59 specimens were prepared. The specimens were hardened with sodium silicate, because some were crumbly and had cracks. To do so, the sodium silicate was diluted with water in a proportion of 1:10. This mixture was then used to cover some specimens completely or for filling cracks of other ones. Material remaining from both block samples was used for magnetic property analyses. These analyses were performed at the Laboratory for Environmental Magnetism   Stud. Geophys. Geod., 64 (2020)

EXPERIMENTAL METHODS
Various rock magnetic analyses were carried out to characterise the magnetic mineralogy of the baked clay from Kiln 1. These include isothermal magnetisation curves with a coercivity spectrometer (Burov et al., 1986;Jasonov et al., 1998) in magnetising fields between +1500 and 1500 mT. For the measurement, about 3 g of powder was pressed into a paper box and mounted on the PVC sample holder of the rotation magnetometer. Additional hysteresis measurements were carried out on selected specimens at 20 and 300 K using a Quantum Design MPMS3 Evercool having a noise floor of less than 0.08 nAm 2 . About 25 mg of powder was pressed into a PVC capsule and clipped into a sample holder made of brass. Three types of measurements were carried out: 1. Zero field cooled (ZFC) hysteresis loops at 20, 150 and 300 K; 2. +5 T field cooled (FC) hysteresis loops at 20 and 150 K; 3. low-field magnetisation warming curves from 20 to 300 K. The experimental procedure regarding the latter was as follows. After the FC hysteresis loop at 20 K, the sample was magnetised in a 5 T field at 20 K, which was applied for 60 s. Then the field was lowered to 1 mT. Next, the magnetisation M 1 mT was measured during continuous warming from 20 to 300 K at a rate of 5 K/min. The small field of 1 mT remained applied until the end of the measurement. The curves of experiment 3 can be considered as FC warming curves.
The decay of the isothermal remanent magnetisation after application of 500 mT field (M r, 500 mT ), was monitored for 100 s using a coercivity spectrometer. This type of measurements allowed determining the decay viscosity coefficient (S d ), a parameter indicating remanence contribution of grains with sizes near the superparamagnetic (SP)/stable single domain (SSD) grain size threshold to the M r, 500 mT . Further, magnetic low-field susceptibility was measured with a Bartington dual frequency sensor operating 0.465 and 4.65 kHz ( LF and  HF , respectively) for determining the frequency-dependent magnetic susceptibility, i.e.

 
, expressed in %, and thus testing the presence of SP grains (Dearing et al., 1995). Finally, the thermo-chemical stability was checked by differential thermomagnetic analysis (DTMA) employing an in-house made astatic Curie balance having a noise floor of about 1030 nAm 2 . Few milligrams of powdered specimen material were heated in air from room temperature to 800°C; at a heating rate of 100°C/min and in a magnetising field of 200 mT. Heating was repeated under same conditions after the specimen cooled down to room temperature, which took about 10 min.
Regarding archaeomagnetic analyses, the remanent magnetisation vector components were measured with an upgraded three-axis, model 760, (2G Enterprises), cryogenic magnetometer with DC SQUIDS. The magnetometer is equipped with a Cryomech 4K, model PT405, pulse tube cryorefrigerator. An automated in-line sample handler and two alternating field (AF) demagnetisation coils allowed three-axial demagnetisation of the 2 cm cuboidal specimens and to measure them at different axial rotation positions (i.e. +x, +y, x, y). The stability of the direction of the natural remanent magnetisation (NRM) was tested by AF demagnetisation. Seven specimens from block sample A and 42 from block sample B, were stepwise demagnetised in ten steps by AF fields with increasing maximal field amplitudes from 0 to 100 mT; sometimes also up to 120 or 140 mT. The three orthogonal remanence components were measured four times in +z position, and for palaeointensity determination also in -z position. The characteristic component of NRM (ChRM) was determined by principal component analysis (PCA) as proposed by Kirschvink (1980) including the origin of the coordinate system.
Four specimens from sample A and six from sample B, were selected for palaeointensity determination following the double heating in-field protocol of Thellier and Thellier (1959). The determination was carried out in twelve steps covering a temperature range from 23 to 550°C. The specimens were heated to the desired temperature, then held for about 40 minutes at this temperature and afterwards cooled to room temperature. A field of +70 T or 70 T was applied during each heating/cooling cycle. A thermal demagnetiser MMTD80A from Magnetic Measurements Ltd., including a calibrated solenoid for thermoremanence (TRM) acquisition, was used for this purpose. Prior to thermal treatment, the applied field of the heating/cooling chamber was mapped with a three axial fluxgate magnetometer in order to determine the zone of uniform field distribution. Within this zone, the field precision varies between 0.1 and 1.5 ppm at an accuracy of about 0.08%. The ChRM direction was obtained by vector subtraction and PCA. Remanent magnetisation was measured with the aforementioned cryogenic magnetometer, and thermochemical alteration was assessed by four partial thermoremanence (pTRM) and residual NRM (rNRM) checks (cf. Fig. A1 in Appendix A for definition). Eventually, the archaeointensity was determined by fitting a weighted slope through the points of the Arai diagram (Nagata et al., 1963) taking into account the standard deviation of both, pTRM and rNRM (e.g. York, 1966). The standard deviations of pTRM, rNRM, pTRM checks and rNRM checks were calculated from the standard deviations of the measured entities following the law of propagation of uncertainty which can be traced back to Gauss (1821Gauss ( , 1823Gauss ( , 1826, see also Sheynin (1979) and references therein.

RESULTS
4 . 1 . R o c k m a g n e t i c a n a l y s e s The room temperature remanence magnetisation curves (Fig. 3a, b) of specimens from both blocks are similar. Isothermal remanent magnetization (IRM) acquisition and backfield curves show a significant change in slope between 60 and 100 mT; field strengths of 1500 mT are insufficient to saturate the magnetisation. There is thus evidence for at least two magnetic mineral populations with significantly different coercive forces, i.e. a dominating low coercivity phase with coercivities mainly between 5 and 60 mT, and a minor high coercivity phase with coercive forces mainly between 200 and 600 mT (cf. insets in Fig 3a,b). This is confirmed by hysteresis loops carried out in fields up to ±5 T (Fig. 3c,d) evidencing a wasp-waisted shape, which indicate also the presence of at least two magnetic mineral populations with strongly contrasting coercivities (cf. Roberts et al., 1995). The wasp-waisted behaviour is even more pronounced at low temperatures (blue curves in Fig. 3c,d). This is probably related to blocking processes at low temperature, affecting grains which have a magnetically viscous behaviour at room temperature. The loops are closed above 3 T, but do not show signs of saturation due to paramagnetic contributions. The field cooled warming curves (Fig. 3e,f) evidence several characteristic slope changes, around 60, 80 K (only B-121), between 160 and 175 K and around 265 K. These may be related to low-temperature magnetic transitions (cf. Section 5.1). The isothermal remanence (M r ) acquisition and backfield curves measured at room temperature with the rotation magnetometer (cf. text). The insets show the 1st derivative, i.e. the remanent coercivity distribution; c), d): hysteresis loops at 20 K (after zero-field cooling, ZFC) and at room temperature (not corrected for high-field contributions, see Appendix C, remarkable slope change around 80 K seen in Fig. 3e is probably caused by viscous grains which become unblocked during warming and which lose partly their magnetisation. These are probably haematite grains, because they seem unaffected by the 1 mT field applied during warming. Regarding the differential thermo-magnetic analysis (DTMA), the progression of the 2nd heating (red curves in Fig. 4a,b) follows largely the one of the first heating. This observation is consistent for both studied specimens, i.e. specimen 1-13 from block

Fig. 4. Top:
Differential thermomagnetic analysis -integral diagrams (upper part) and differential diagrams (lower part) of the first and second heating of а) specimen 1-13 from block sample A, b) specimen 2-24 from block sample B. Bottom: Decay of remanent magnetisation at field of 500 mT (M r, 500 mT ) of c) specimen 1-141 from block sample A, and d) specimen 2-211 from block sample B. The IRM is normalised to the IRM at 0.4 s. After rescaling the time axis logarithmically, the data was fitted with a linear model. The slope of the best fit represents the decay viscosity coefficient S d . Other magnetic properties measured on the same specimens are given for comparison: B cr -coercivity of remanence, B c -coercive force, M desc, 500 mT -induced magnetisation on the descending branch of hysteresis at 500 mT,  FD -frequency-dependent magnetic susceptibility. sample A and 2-24 from block sample B. A significant magnetisation increase in the 2nd heating curve is not observed, which would be the case if new minerals had formed during the 1st heating. The derivative plots of 1st and 2nd heating are rather similar, except for the interval 420600°C. Here, some minor change is observed in both specimens, but which does not change significantly the magnetisation at room temperature. The curve progression is slightly different when comparing both specimens and indicates a variable composition of magnetic minerals. In general, the DTMA reveals only very minor magneto-mineralogical alteration when heating the specimens up to 800°C.
Measurements of IRM decay with time reveal the presence of magnetically viscous grains. The data can be well fitted with the theoretical exponential decay model (cf. yellow curves in Fig. 4c,d), but some deviation from the mean slope is observed at short times 03 s. Probably grains with different relaxation times are present in the sample. A model with two relaxation times yields a better fit, but as the analysis of viscous effects falls beyond the scope of this paper, this will not be further discussed. The examples given in Fig. 4c,d show that the decay viscosity coefficient S d varies amongst the different specimens. Between 1 and 11 % of the M r, 500 mT is lost in 100 s.
The frequency-dependent magnetic low-field susceptibility was determined as mentioned previously and expressed as  FD (cf. Section 3). The average is 7 ± 2% covering a range of values between 1 and 11%. This indicates, in addition to the viscous fraction (cf. above), a variable presence of superparamagnetic grains with sizes close to small stable single domain gains.

. 2 . A r c h a e o m a g n e t i c d i r e c t i o n s
Alternating field demagnetisation is insufficient to completely demagnetise the NRM in most of the specimens (Fig. 5). In some specimens about 60% of the NRM still remains after exposition to AF peak fields of 120 mT. This proves the presence of highly coercive minerals partly carrying the NRM. Most specimens show a single NRM direction (cf. Zijderveld diagrams in Fig. 5,. It is considered as primary because the direction seen in Zijderveld diagrams (Zijderveld, 1967(Zijderveld, , 1975 can be well extrapolated towards the origin of the diagram. Moreover, the maximum angular deviation (MAD) of the regression line varies only between 0.1 and 0.6°. In individual cases, a small secondary overprint is evident and was removed by AF peak fields of 5 or 10 mT (e.g. Fig. 5, specimen 2-322). The characteristic component of the NRM (ChRM) was determined by line fitting, generally using the field steps between 5 or 10 and 100 mT; the origin of the Zijderveld diagram was taken as additional point.
Two approaches were followed for obtaining the mean ChRM direction for Kiln 1. Because both block samples A and B were oriented independently, the first, hierarchical approach consisted in calculating the Fisher mean for sample A and B separately. The average direction for Kiln 1 was calculated then from both sample means. Results are given in Table 2. The mean inclination of both samples is very close, but declination varies by about 4°. The combined kiln mean falls fairly between both block sample means, but the 95% significance level ( 95 ) is rather large with about 5°. Therefore, in a second approach, the kiln mean was calculated as Fisher mean of all specimens, consequently resulting in a much lower  95 . The application of Fisher statistics was justified by a probability quantile-quantile plot (not shown) following Fisher et al. (1987). The two test statistics for declinations and inclinations are M u = 1.028 and M e = 0.567, respectively. They are well below the critical values for rejecting the hypothesis of a Fisher distribution at a 95% significance level, i.e. 1.207 crit u M  and 1.094 crit e M  . The hierarchical approach was not considered for further interpretation of the kiln mean ChRM direction, because it is based only on two samples only.
Thermal demagnetisation was not performed in the usual way. Instead, this information was extracted from the archaeointensity experiments by vector subtraction (cf. Section 3.3, Table 3). The thermal treatment appears to be more successful compared to AF demagnetisation. Depending on the specimen, only between 5 and 25% of the NRM remain after heating to 550°C; (cf. blue curves in Fig. 6

. 3 . A r c h a e o i n t e n s i t y
Rock magnetic results and NRM AF demagnetisation curves suggest the presence of low and high coercive remanence carriers. Nevertheless, archaeointensity determination was attempted because the material appeared to be well-baked and in thermo-chemical equilibrium (cf. Fig. 4), and because the individual specimens recorded the geomagnetic field direction with high confidence (cf. Fig. 5). In order to identify possible thermochemical alteration, conventional pTRM checks were carried out. In contrast, the equality of blocking and unblocking temperatures, which is another assumption of the Thellier and Thellier (1959) method, was not checked as usual, i.e. by an additional zero field cooling step. Instead, this equality was tested by combining the magnetisation vectors for positive and negative fields ( Details of the derivation are given in Appendix A. A weighted slope was determined by fitting a linear model to the data using between 7 and 10 points. The coefficients of determination R 2 are ≥ 0.9927 for all ten specimens. All pTRM and rNRM checks fall almost always into the 99.99% confidence interval of the mean slope and can therefore be both considered as successful. It is hence concluded that the samples are thermally stable and do not alter during the archaeointensity thermal treatment. The equality condition of unblocking and blocking temperature is also satisfied. Further, the Coe quality parameters (Coe et al., 1978) are comparable to those of other studies dealing with massive baked clay (e.g. Gómez-Paccard et al., 2006). Examples of archaeointensity determinations are given in Fig. 6. The Fisher average of the residual NRM directions and its uncertainty is comparable to the independently determined ChRM directions discussed in Section 4.2 (cf . Tables 2 and 3). The pTRM was also properly acquired: its mean Fisher inclination deviates only 0.4 ± 0.8° from 90°. The archaeointensity data analysis demonstrates that NRM and TRM vectors could be separated correctly, but also that both remanences were respectively properly demagnetised and acquired. This was not expected given the complex magnetic mineralogy. Certainly, the latter has some influence on the uncertainty envelope of the archaeointensity, e.g. specimen 2-221 in Hysteresis and isothermal remanent magnetisation curves indicate both the presence of low and high coercivity components, both able to carry remanent magnetisation (cf. insets Fig. 3a,b). The studied material can be considered as well-heated under oxidising conditions during kiln operation. This is evidenced by its strong reddish colour, magnetic property and archaeointensity analyses.
In general, magnetite-bearing material from burned structures is blackish or blackish/greyish in colour, and its magnetic properties change markedly when heated above 150°C under oxidising conditions (Spassov and Hus, 2006). This is not the case for the specimens studied. They are reddish in colour and are thermally stable, even after heating to 800°C under oxidising conditions (cf. Fig. 4). Magnetite and iron sulphides can therefore be excluded as low coercivity remanence carriers. However, the DTMA experiments still indicate minor changes between 420 and 600°C, when comparing 1st and 2nd heating (cf. Fig. 4a,b derivative plot). This is interpreted as oxidation of small magnetite contributions embedded and a haematite matrix, rather than a maghaemite-tohaematite transformation. The haematite matrix had protected these magnetite grains from oxidation during kiln operation but not during laboratory heating up to 800°C. Indeed, the FC warming curves evidence a characteristic slope change at 80 K (Fig. 3e), which can be associated with the Verwey transition of non-stoichiometric magnetite (Shepherd et al., 1985). However, this phase was detected only in one of the four specimens analysed. In the remaining three there was no indication for a Verwey transition, neither around 80 K nor around the ordinary temperatures, i.e. between 110 and 125 K. The low coercivity component is therefore associated with maghaemtite in most of the specimens.
Haematite and goethite are potential candidates for the high-coercivity phase. Goethite dehydrates and is transformed via magnetite into haematite when heated under laboratory conditions (cf. Özdemir and Dunlop, 2000). It is thus unlikely that pristine goethite survived the multiple heating-cooling cycles during operation. In contrast, post-baking goethite may have formed in cracks or fissures, i.e. goethite formation after the kiln was abandoned. Supposing this was the case, a significant difference between 1st and 2nd DTMA heating would be expected. This however is not the case and excludes goethite completely as high coercivity remanence carrier. The FC warming curves show in several samples two slope changes: at 160175 K and at 265 K. The Morin transition (T M ) of haematite is not fixed at a certain temperature; albeit the observations of Morin (1950), who observed it at 250 K. Instead T M covers a broad temperature range depending on grain size. The T M of MD haematites ranges from 250 to 265 K, the one of submicron SD haematites from 238 to 260 K, while a T M between 238 and 245 K may reflect sizes below 100 nm (Morrish, 1994;Özdemir et al., 2008). The observed slope change in the FC warming curves (Fig. 3e,f) is thus interpreted as Morin transition and related to MD sized haematite grains. However, another slope change is observed around 160175 K, which is also interpreted as T M . In contrast to the previous one, it indicates a smaller haematite grain size fraction and/or the Al substitutions in the haematite crystal lattice. In unsubstituted haematites T M is either above 180 K or absent (Vandenberghe et al., 2001), but this seems not the case according to Fig. 3e,f and insets. It is therefore argued that indeed Al-substituted haematites are present in the samples. This argumentation is reinforced when taking into account that samples are mostly composed of clay which contains large amounts of Al. According to the T M range observed here, i.e. 160175 K, the smaller haematite phase may contain up to 8% Al (da Costa et al., 2002).
The thermomagnetic curves in Fig. 4, particularly in Fig. 4a indicate a significant slope change around 190°C. As this slope change occurs in both heating and cooling curves it may me associated with the Curie temperature (T C ) of -Fe 2 O 3 , a metastable iron oxide phase with reported T C between 207 and 237°C (e.g. Dézi and Coey, 1973;Popovici et al., 2004;Tuček et al., 2010).
-Fe 2 O 3 forms only under certain conditions: 1. Sufficiently high decomposition temperatures of -Fe 2 O 3 and 2. Restrictions on the diffusion of the conversion gas and the agglomeration of the nanoparticles formed (cf. Machala et al., 2011, and references therein). Both conditions may be fulfilled in baked clay from archaeological artefacts, because they are operate at high temperatures and the clay contains silica layers. The latter could prevent agglomeration or limit gas diffusion.
Archaeological baked clay is often heated above 800°C, i.e. above the -Fe 2 O 3 → → -Fe 2 O 3 transition temperature. Moreover, the heating rate in kilns etc. is rather slow, i.e. °C/h. In contrast, Brázda et al. (2014) applied a heating rate of 10°C/min when carrying out their differential thermal analysis. Even at this much faster rate they observed the  →  transition. Therefore one would not expect -Fe 2 O 3 to occur in baked clay from archaeological sites. If its presence is still proven, such as for instance by Raman or X-ray diffraction (e.g. López-Sánchez et al., 2017;Calvo-Rathert et al., 2019), this means that the archaeological artefact had not been heated above 750°C.
As seen in Fig. 4a,b, the samples were heated up to 800°C, i.e. above the  →  transition. This transitions occurs, depending on heating rate and -Fe 2 O 3 crystal size, between 750 and 825°C (Dézi and Coey, 1973;Brázda et al., 2014). A possible conversion to haematite can therefore not be excluded in our case. Figure 4a,b may thus not serve as evidence for -Fe 2 O 3 presence.
Hysteresis loops measured at different temperatures can also provide evidence for -Fe 2 O 3 presence. The coercive force should decrease from 300 to 80 K and then increase again at lower temperatures (e.g. López-Sánchez et al., 2016). Such behaviour was not observed. Instead, the coercive fore increases steadily from 300 to 150 to 20 K (cf. insets in Fig. 3c,d) in both samples. Gich et al. (2005) and López- Sánchez et al. (2016) observed that the coercive force at 150 and 20 K is either similar or significantly lower than that at 300 K. Regarding ultra-small particles, i.e. < 4 nm, the coercive force changes only very little between 4 and 300 K (Dubrovskiy et al., 2015). As all this is not the case for samples B1-121 and B2-H (insets in Fig. 3c,d), hysteresis data disproves an -Fe 2 O 3 presence in the samples studied here.
Although the specimens carry a stable NRM direction and are able to acquire laboratory remanences (pTRM, IRM), magnetic analyses revealed the presence of SP and magnetically viscous grains (Fig. 4c,d). The latter are responsible for the small viscous NRM overprint as discussed in Section 4.2. Magnetic viscosity was tested only for a few specimens, which neither allowed the amplitude of S d to be associated with a specific magnetic mineral, nor to identify differences between block sample A and B. However, the two examples show that S d and  FD factor are significantly higher in the sample with higher coercivity of remanence B cr .
In summary, the argumentation identifies maghaemite and haematite as the principal remanence carriers, while magnetite occurs in trace amounts. The haematite is partly Al-substituted according to low-temperature experiments. The AF demagnetisation shows that the NRM is mainly carried low and high coercivity magnetic minerals, i.e. maghaemite and haematite. Despite the presence of haematite, a large part of the NRM is already demagnetised at 550°C. This may be due to lower unblocking temperatures because of the small haematite grains and a lowering of the Néel temperature of the Al-substituted haematite phase (e.g. Jiang 2014Jiang , 2015. The remanence carriers are in thermo-chemical equilibrium, and therefore it is expected the pTRM acquired during the archaeointensity experiment is not affected by magneto-mineralogical changes. The pTRM checks should be positive, i.e. falling close to the mean slope in the Arai diagram. The presence of magnetically viscous grains justifies to exclude the first point in the Arai diagram from the mean slope fit.

. 2 . A r c h a e o i n t e n s i t y
The majority of the specimens complies well with the most commonly used parameters for assessing the quality of individual palaeointensity results, such as for instance n > 5, f ≥ 0.35,  ≤ 0.1, q ≥ 5, MAD anc ≤ 6, n pTRM ≥ 3 (e.g. Paterson et al., 2014, and references therein). A description of these parameters can be found at https://earthref.org/PmagPy/SPD/home.html.
Having a closer look to the results in Table 3 one recognises that the individual determinations can be divided roughly into two classes, one around 62 T and another around 72 T. Both classes have rather excellent quality indicators. This raises immediately the question which of the two classes yields the correct result. Certainly the variable contribution of low and high-coercive phases between specimens accounts most probably for differences in the unblocking temperature spectra and may probably be responsible for this classification. However, it would be subjective to exclude one or the other class.
In order to obtain the archaeointensity of a sample or a structure, a mean or a weighted mean (e.g. Coe et al., 1978;Chauvin et al., 2000) can be used. Using an arithmetic mean implies that the uncertainties of the individual determinations are neglected and that given a small specimen set, the mean would be strongly biased by extreme values. In contrast, a weighted mean reduces only the influence of extreme values associated with large uncertainties and is not affected by a single value that has a large uncertainty but is close to the mean. Summarising, in the case of the present archaeointensity results, neither mean nor weighted mean yield a meaningful average value. This holds also true when referring to the median, which is indeed less influenced by extreme values but fails averaging a multivariate distribution of values. Regarding a site/structure, palaeointensity values of individual samples/specimens are often grouped by more than one value. This may be due inhomogeneous baking temperature distributions resulting in magnetomineralogical and unblocking temperature differences. Therefore another approach based on Bayesian principles was developed as explained in the following. An archaeointensity experiment of a single specimen yields a value F and its associated measure of uncertainty, s. The latter is taken as the standard error of the slope in the Arai diagram. However, F is only a single value and other values of F, called f i , might also be probable.
Assuming that f i are distributed in some way around F, their true standard deviation  i is unknown. Because the Bayesian concept considers all unknown quantities in a probability model for observed data as random variables (e.g. Wakefield, 2013), it is justified to assume that f i  F and hence f i are normally distributed.
It is clear that the true unknown uncertainty,  i , can be hardly less than the measures of uncertainty, s i . The prior distribution   i P  , should tend to zero as  i gets large and should be zero for  i < s i . A way to express this mathematically is: following Jeffrey's (1932) initial concept of a prior proportional to 1  and Sivia's (1996) suggestion of a proportionality to 2 1  that simplifies the calculations.
The posterior probability of an archaeointensity f at specimen level can now be calculated, given that f i are normally distributed around F: Solving the integral yields then the probability distribution for a specimen intensity determination f In this way cogent information from the archaeointensity experiment is taken into account in the form of s without knowing the true individual measurement uncertainties  i .
A number of recent studies have stressed the importance of correcting archaeointensities for differences in the rate at which the studied artefacts cooled in the laboratory and in the environment in which they were previously heated. Theoretical considerations suggest a figure of a few percent for each order of magnitude difference in cooling rate in the case of single domain grains (Fox and Aitken, 1980;Dodson and McCelland-Brown, 1980). It is difficult to estimate how rapidly an artefact may have cooled in the past. A reasonable assumption for Thellier-Thellier experiments would be that a difference in cooling time between one and two orders of magnitude could be expected between the experiment, i.e. 30 to 60 min and the original heating, e.g. ~25 h to 72 h. Using these figures a cooling rate correction of between 2.5 and 20% may be necessary, according to Fig. 3 in Poletti et al. (2013) and references therein. The decade logarithmic relation between the cooling rate and the magnitude of correction suggests a simple distribution for correcting, such as for instance a probability proportional to the inverse of the size of the effect. Denoting the cooling rate correction as r, an appropriate prior probability distribution would be: The posterior probability of the cooling rate corrected intensity is calculated by multiplying Eq. (6) with Eq. (5) and by replacing the term f i  F by takes into account the difference between f i and the cooling rate corrected mean.
The marginalisation of r yields then the likelihood for the field strength at specimen level F The likelihood of function for the ancient field strength for a number k of archaeointensities from the same feature F is then Stud. Geophys. Geod., 64 (2020) Under the assumption that there is no prior knowledge of the field strength, the final probability can be directly equated to the probability of the field strength within the window of time represented by the specimens. Albeit prior information of the field strength could be implemented from other geographical locations or sources (e.g. Arneitz et al., 2017Arneitz et al., , 2019, this was disregarded here in order not to disturb the independence of the current observation. The probability functions of all 10 specimens calculated using Eq. (8) are shown in Fig. 7a, while the combined probability is given in Fig. 7b. According to the latter, the cooling-rate corrected archaeointensity F pal for Kiln 1 falls with a probability of 95% within the interval between 58.6 and 62.9 T.
Posterior probability distribution functions cannot be characterised by standard deviation and mean. In order to still compare the result of Fig. 7b with the weighted average given in Table 3 and with other published data, the length 95% confidence interval is used. Regarding the Bayesian approach this interval is 4.3 T while for the weighted mean approach it is 7.2 T.
The Bayesian approach used here has helped to overcome the initially mentioned dilemma which class to consider as correct archaeointensity. A narrowing of the confidence interval was also achieved.

. 3 . C o m p a r i s o n w i t h o t h e r r e g i o n a l P S V d a t a
There is a general interest to extend the geomagnetic field record into the past and to combine these results with theoretical field models or to use them for analyses of the geomagnetic field (e.g. Hongre et al., 1998;Jackson et al., 2000). Therefore the main objective of the present investigation is to provide an independently dated Fig. 7. a) Archaeointensity F pal probabilities for Kiln 1 at specimen level taking cooling rate constraints into account; calculations according to Eq. (8). b) Archaeointensity determination at structure level by multiplying the probabilities shown in a) according to Eq. (9). The archaeointensity falls with a probability of 95% into the interval between 58.6 and 62.9 T. archaeomagnetic field vector for Kiln 1, based on ChRM direction/archaeointensity and archaeological context dating. Such investigations are of particular importance, as archaeomagnetic results around the Bolgar region and its neighbourhood are rare or inexistent as initially mentioned in the introduction.
The comparison of the PSV results obtained for Bolgar with nearby data is almost impossible. An exception, however, are data from sediment cores of Lake Aslikul which is only 364 km away, but which yields only directional data (cf. Nourgaliev et al., 1996). The chronostratigraphy of Lake Aslikul sediments relies on pollen samples and was further improved by comparing the relative pollen palaeotemperature distribution with a 14 C dated regional palaeotemperature master curve (cf. Nourgaliev et al., 2000a and references therein). The Lake Aslikul inclination and declination record starts at 600 yr B.P., which corresponds exactly to 1350 A.D..
The sampling site in the Aral Sea is with about 1200 km already much more distant from Bolgar than Lake Aslikul. The advantage of using data from the Aral Sea for comparison is the rather high sedimentation rate of about 16 mm/yr between 1295 and 1470 A.D. This record was dated with 14 C on macrofossils present in the lake sediment. These datings were calibrated using the OxCal program (Bronk Ramsey, 1995) which allowed to establish a confident time scale for this sediment record . The value extracted from the Aral Sea PSV curve refers to 655 yrs BP, i.e. 1295 A.D.
The PSV records from lakes Naroch and Svir, which are about 1400 km away from Bolgar, were dated with 14 C allowing the construction of a timescale combining data from both lakes (Nourgaliev et al. 2000b. Because of the low sedimentation rate in both lakes, the age model of Nourgaliev et al. (2005) takes into account the phenomenon of delayed acquisition of natural remanent magnetisation (cf. Bleil and von Dobeneck, 1999).The PSV value taken from this record refers to 600 yrs BP, i.e. 1350 A.D.
Archaeointensity data for comparison were taken from recently published studies from baked materials of the Moscow and Novgorod region which are respectively about 780 and 1150 km away from Bolgar (Salnaia et al., 2017a,b). Baked clay brick fragments were from mediaeval churches were investigated in both studies and the archaeointensity determined with the Triaxe method (Le Goff and Gallet, 2004). The uncertainty in age is with 1 to 7 years negligible small, thanks to historical sources. As however non-oriented samples were available only, no directional information could be obtained.
Additionally, older published data was taken from the catalogue of Rusakov and Zagniy (1973) containing Moldavian and Ukrainian archaeointensities from burned archaeological structures such as kilns, furnaces and hearths. The quality of these data is inferior compared to those published by Salnaia et al. (2017a,b). This results on one hand from the much higher dating uncertainties and on the other from the low number of samples/specimens studied per site. Moreover, it is unclear whether these data were cooling rate corrected or not. Unfortunately, Rusakov and Zagniy (1973) did not give any units for the archaeointensity, but is suggested that their values are given in Gauss (Table 1, column 10 in Rusakov and Zagniy, 1973). Another problem was to determine the uncertainty of their values, because only the uncertainty of the average ratio of archaeointensity ( F  ) to present day intensity (F 0 ) is given. Checking some of their data reveals that they used slightly different present day field intensities. It is also not clear to Stud. Geophys. Geod., 64 (2020) which geographic location these intensities would refer to. In order to calculate the uncertainty of their archaeointensities, the present day field intensity was recalculated from their data dividing F  by the average 0 F F  . This procedure was still preferred instead of using a general International Geomagnetic Reference Field (IGRF) value for 1973 A.D. referring to the geographic centre of the Ukraine. The as obtained F 0 were then multiplied with the corresponding mean square errors of the average 0 F F  . The sites selected for comparison are between 1200 and 1450 km away from Bolgar.
All directional and intensity data discussed here above and plotted in Fig. 8 were relocated to the geographic coordinates of Bolgar. This was done following the 'conversion via pole' method (Tarling, 1988;Noel and Batt, 1990) for declination and inclination. Regarding archaeointensity relocation, a geocentric axial dipole model was employed.
As thoroughly discussed by Casas and Incoronato (2007), relocation of geomagnetic field vectors results in a systematic error which depend inter alia on the geographic region. While for most European studies relocation uncertainties are inferior compared the experimental ones and thus negligible (e.g. Tema and Kondopoulou, 2011;Tema et al., 2013), this is not the case in the present study. Values of 0.2° and 0.35 T per 100 km were extracted from Fig. 5 in Casas and Incoronato (2007) for the western part of Russia for direction and intensity, respectively. The most distant data relocated to Bolgar are roughly 1400 km away and yield thus maximal uncertainties up to 2.8° and 4.9 T which both exceed the experimental uncertainties. Accordingly, relocation uncertainty ellipses were calculated for depending on the actual distance of the considered sites to Bolgar and plotted together with the experimental uncertainties in Fig. 8.
There are different possibilities to compare the obtained PSV results with existing data. One approach would be to compare the results of Bolgar with well-established PSV curves, which exist for countries of geographic Europe like Bulgaria, France, UK and many others. Another would be to establish a regional PSV model. The far distance to Europe, i.e. > 2000 km, and the lack of regional data are neither favourable for relocating existing European PSV curves to Bolgar (cf. Casas and Incoronato, 2007) nor to construct a regional PSV curve such as for instance through spherical cap harmonic analysis (e.g. Pavón-Carrasco et al., 2009). Instead, a model performant for global PSV calculations was used for comparison; concretely model CALS10k.1b from Korte et al. (2011). Albeit the lower spatial and temporal coverage between 0 and 3 ka due to bootstrap averaging, CALS10k.1b was preferred against other models (e.g. CALS3k.4) because it emphasises the robust field structure at the core mantle boundary, which is a reasonable approach given the lack of PSV data around Bolgar. This approach was implemented via the GEOMAGIA50.v3 (Brown et al., 2015). In addition, the Bolgar data are compared with the recently developed PSV model BIGMUDI4k. 1 of Arneitz et al. (2019), a global self-consistent model that is based on Bayesian iterative principles. The difference between both PSV model curves is attributed to a) differences in input data, i.e. BIGMUDI4k.1 does not include continuous sedimentary records, which means that the region of interest is not supported by data over the time interval of interest, b) different time regularisations, whereby BIGMUDI4k.1 is less regularised, due to the higher data density after 1600 AD.  (Korte et al., 2011, blue curve, with orange uncertainty envelope) and BIGMUDI4k.1 (Arneitz et al., 2019, red curve, with green uncertainty envelope). Regarding a) and b), areal PSV directions (black) were extracted from sedimentary records of the Aral Sea , Lake Aslikul  and Lake Naroch/Svir (Nourgaliev et al., 2002b. The uncertainties of declination and inclination were derived from marginalised variances (cf. Appendix B). The archaeointensities in c) for the Moscow and Novgorod regions are from Salnaia et al. (2017a,b) and for Ukraine from Rusakov and Zagniy (1973). Black ellipses refer to the experimental uncertainties of the published data, while grey ellipses refer to the estimated relocation uncertainty (see text). The original, nonrelocated data are summarised in Appendix C, Tables C.2 and C.3.
Stud. Geophys. Geod., 64 (2020) The lack of sedimentary data and lesser time regularisation leads to larger variations of BIGMUDI4k.1 with respect to CALS10k.1b before the onset of historical field records, particularly for the declination (Arneitz A. and Egli R., Zentralanstalt für Meteorologie und Geodynamik, Vienna, personal communication, 2020).
The determined mean ChRM vector for Kiln 1 of D = 359.9°, I = 69.3° and F pal between 58.6 and 62.9 T is fairly in agreement with the CALS10K.1b and BIGMUDI4k.1 model curves calculated for Bolgar for the time interval of 13401360 A.D., i.e. D = 0.9°, I = 69.2°, F pal = 61.9 T, and with other PSV data from Belorussia, Kazakhstan, Russia and the Ukraine, cf. Fig. 8.

CONCLUSIONS
Successful archaeodirections and intensities were obtained from Kiln 1 unearthed in ancient Bolgar. Despite the fairly good quality of the determined archaeointensities, a difference of ~10 T was found between both studied block samples. This is due to a variable magnetic mineralogy consisting of haematite and maghaemite. The magnetic mineralogy, however, neither influenced the acquisition of thermoremanent magnetisation nor the directional stability of the natural remanence during demagnetisation. This is evidenced by congruent ChRM directions from both, thermal and alternating field demagnetisation and accurate TRM acquisition. The latter being almost parallel to the applied field, having small uncertainties. A present viscous remanence component, as revealed by IRM decay experiments, had only a minor influence on the NRM directional stability. Equally present superparamagnetic grains were in thermo-chemical equilibrium and did not provoke neoformation of remanence carriers during the rockmagnetic and archaeointensity experiments.
The newly introduced concept for testing the equality of blocking and unblocking temperature was successfully applied and demonstrates that the commonly used zero-field cooling step for testing this condition can be foregone.
It was further demonstrated that neither a mean nor a weighted mean yields a meaningful average intensity value at structure level. Instead, a new approach based on Bayesian principles was developed to determine a representative archaeointensity at structure level. This Bayesian approach is particularly valuable when only a few intensity data are available for an archaeological structure. According to this new approach, time consuming and alteration affected cooling rate experiments are not an imperative need for reliable archaeointensity determination.
The comparison of the ChRM vector determined for Bolgar with a global geomagnetic field model and with other PSV data from the wider area validate the obtained archaeomagnetic results. Given the small uncertainty in age and geomagnetic entities, the data of Bolgar can be classified as of high quality and shall be included in a future regional PSV curve for western Russia.
Eventually, this study underlines once more the importance of the systematic relocation error. In the present case it exceeds the random experimental uncertainty, when relocating published PSV data over large distances.
Instead, the vector sum of Rearranging Eq. (A.10) yields: where n is the number of specimens and f the 95% quantile of the Fisher distribution (cf. The uncertainties ellipses of declination and inclination shown in Fig. 8a,b were calculated according to Eq. (B.4). The uncertainty envelopes of the PSV model were calculated in a similar manner (cf. Brown et al., 2015).
APPENDIX C   [K] no indication T M [K] 175182 and 265267  Table C.2. In the alternatingfield (AF) method, D and I is declination and inclination, respectively, obtained by the corresponding method;  95 is semi-angle of the cone of confidence around the mean direction in which the true direction occurs with a confidence level of 95%; n is number of samples or specimens used for demagnetisation experiments; lower and upper F pal stand for the lowest and highest probable archaeointensity value, respectively; Meth. -method used for archaeointensity determination; TT is Thellier-Thellier method, TR is triaxial method; n is number of samples or specimens used; CR means cooling rate correction. References are the same as in Table C