Multi-Stage Magma Evolution in Intra-Plate Volcanoes: Insights From Combined in situ Li and Mg–Fe Chemical and Isotopic Diffusion Profiles in Olivine

Understanding the timescales of magma evolution and ascent is essential for interpreting geophysical monitoring signals from active volcanoes. In this study, we explore the potential of diffusion-driven Li concentration and isotope zoning profiles recorded by magmatic olivine crystals to unravel time scales of magma evolution processes. Lithium is a fast-diffusing element and may provide the opportunity to investigate changes in magma composition during magma ascent, shortly before eruption. Lithium chemical and isotopic profiles were determined in olivines from two localities in the Massif Central volcanic region (France) that have previously been investigated for their Fe–Mg isotope systematics. The combined investigation of isotopic and chemical profiles makes it possible to distinguish between crystal growth and diffusion events. Extremely low δ7Li-values down to −30.7‰ (relative to the commonly used Li isotope standard IRMM-16) in the crystal core regions and elevated values at crystal rims (δ7Li ∼8 to 10‰), along with increasing concentrations from cores (∼3 to 1 μg/g) toward rims (12 to 6 μg/g) were found. The shape and orientation of both the chemical and isotopic profiles indicate that they were dominantly generated by Li diffusion into and within the olivine grains during magmatic differentiation. While Mg–Fe isotope and major element profiles have been modeled by a single diffusion event (Oeser et al., 2015), concentration and isotope profiles of Li indicate that a second diffusion event took place, that was not recorded by the Mg–Fe exchange diffusion couple. The first diffusion event was interpreted as reflecting the residence of the olivine crystals in a magma chamber. As diffusion coefficients for Fe–Mg exchange diffusion are very well determined, the time scales of this event are likely best quantified by Mg–Fe isotopic exchange diffusion modeling (Oeser et al., 2015). This event probably also generated the low δ7Li observed in olivine cores. Comparing the length of the Mg–Fe and Li profiles could thus be used to determine the less well-known diffusion coefficients of Li in the studied olivine crystals. The findings of this study indicate that Li diffusion at low Li concentration levels, as typically observed in natural olivine, may be not as fast as previously thought. The second diffusion event might represent a short-lived event, such as degassing, related to the ascent of the magma and/or magma cooling after emplacement of the lava. Such a process would only affect Li, which, in contrast to the refractory elements Fe and Mg, is volatile during degassing. The findings of this study show that, according to their different diffusion rates and physiochemical properties, the combined use of spatially resolved Li and Mg–Fe chemical and isotopic diffusion profiles, is a powerful tool to model even multi-stage evolution processes in magmatic systems.


INTRODUCTION
In principle, chemically zoned crystals can be the result of crystal growth in an evolving melt (growth zoning) or of chemical diffusion, as the result of disequilibrium, between minerals and the surrounding melt (diffusion zoning) (Costa et al., 2008). In the case of Li, this can be caused either by an increase in Li concentration in the magma, typically either due to magma mixing or crystal fractionation, or, due to a decrease in Li concentration in the melt, e.g., during degassing (Vlastélic et al., 2011;Lynn et al., 2018). The resulting zoning can be used to reconstruct the migration of crystals through dynamic plumbing systems with several reservoirs of compositionally different melts (Kahl et al., 2013(Kahl et al., , 2015. Olivine and clinopyroxene, as early forming crystals, are the most likely minerals to accept Li in their structure due to similar ionic radii of Mg and Fe in the octahedrally coordinated sites (Shannon, 1976). Nevertheless, Li behaves incompatibly during magma differentiation, with typical concentration levels of 3 to 8 µg/g in basalts and of ∼20 µg/g in rhyolites (Ryan and Langmuir, 1987;Ryan and Kyle, 2004) and with olivine/melt distribution coefficients of 0.2-0.35 (Ryan and Langmuir, 1987;Brenan et al., 1998). As olivine is a very abundant mineral in primitive basalts and chemical diffusion in olivine is well characterized, e.g., for Fe-Mg exchange  or for Li (Dohmen et al., 2010;Richter et al., 2017), it frequently serves for diffusion studies (e.g., Lynn et al., 2018;Oeser et al., 2018). However, growth and diffusive origins of zoning cannot easily be distinguished by the investigation of chemical zoning alone and, notably, only the latter bears timescale information.
At magmatic temperatures, diffusion results in large isotope fractionation, e.g., of Li, Fe, and Mg, which may be recorded as isotopic zoning (e.g., Richter et al., 2003). Isotopic zoning coupled with chemical zoning is a strong indicator for a diffusive origin of the zoning because high-temperature equilibrium isotope fractionation is very limited, e.g., for Li (Tomascak et al., 1999b;Jeffcoate et al., 2007;Parkinson et al., 2007), Fe (Weyer and Ionov, 2007), and Mg (Teng et al., 2007;Liu et al., 2011). The combined information of chemical and isotopic zoning of Fe and Mg were successfully used in a number of recent studies to investigate complex magmatic evolution processes, such as magmatic differentiation and magma mixing, in a variety of different settings and to evaluate the effects of the end member processes of crystal growth and pure diffusion on the development of compositional zoning Sio et al., 2013;Oeser et al., 2015Oeser et al., , 2018Collinet et al., 2017). Lithium is the lightest alkali metal and has two stable isotopes with a relative mass difference of ∼17%. Accordingly, large isotopic fractionations are observed between different geochemical reservoirs (see compilation in Penniston-Dorland et al., 2017 and references therein). Large Li isotope fractionations have also been observed on the mineral scale, e.g., up to 29 in a clinopyroxene from the Solomon Islands volcanic rocks (Parkinson et al., 2007), 29 in San Carlos olivine and 38 for an orthopyroxene crystal (Jeffcoate et al., 2007), which is thought to be mainly caused by chemical diffusion. However, the origin of the low δ 7 Li-values in olivines remain debated. They may be directly caused by subduction-derived metasomatized material, such as isotopically light lower crust (Hamelin et al., 2009) or the ingress of light subduction-derived fluids in the mantle (Gu et al., 2016). On the other hand, they are assumed to be the result of Li diffusion in the mantle or during magmatic processes such as subduction or magma differentiation (Marschall et al., 2007a,b;Magna et al., 2008).
Here, we explore the suitability of using Li and Li isotope profiles recorded by chemically and isotopically zoned magmatic olivine to unravel magmatic processes and timescales of magmatic events. As Li is thought to diffuse faster than the Fe-Mg diffusion couple and also has specific physiochemical properties in melts such as volatility, it may record short magmatic processes which occur just before (or even after) eruption, and which are not otherwise recorded. The aim of this study is to investigate if the fast-diffusing Li isotope system is coupled with or decoupled from Mg-Fe exchange diffusion (Weyer and Seitz, 2012). We focus on chemically zoned magmatic olivine crystals from the Massif Central continental intra-plate volcanic system, which were previously investigated for Fe-Mg inter-diffusion processes and diffusion-driven Fe-Mg isotopic zoning of olivine , thus providing a framework for investigating Li in this study. Li concentrations and isotope ratios of olivine are analyzed by a newly developed matrix independent measurement method by femtosecond-laser ablation-MC-ICP-MS . By combining diffusion profiles of Fe-Mg and Li, we intend to decipher complex magma evolution scenarios, including cooling, magma mixing and degassing. Massif Central (France) were investigated in this study. Olivinephyric basanites from Roche Sauterre may originate from a large lava flow that merged into a small lava lake in a paleo-topographic low during the late Miocene to early Pliocene (Nehlig et al., 2001;Lorand et al., 2003;Richet, 2003). At the summit of the Banne D'Ordanche, an olivine-and clinopyroxene-bearing ∼710,000-year-old basanite from the north-west slope of a former strato-volcano was sampled (Richet, 2003). The samples were characterized regarding their major and trace element composition by Oeser et al. (2015). All investigated olivines are normally zoned with respect to Mg# (Mg# = [Mg]/([Mg] + [Fe])) with high Mg# (0.86 to 0.90) in their cores and low Mg# (0.74 to 0.8) in their rims. Nickel concentrations follow this trend. The width of Mg# zoning reaches up to 400 µm into the up-to-2 mm in diameter olivine crystals. Mg# zoning is accompanied by inversely correlated Mg and Fe isotopic signatures underlining a diffusive origin of the zoning. The inter-correlation indicates an inter-diffusion process with Mg diffusing out of the olivine and Fe diffusing into the olivine . Lithium concentrations and Mg# are summarized in Supplementary  Table S1 and backscattered electron images of the investigated crystals are displayed in Supplementary Figure S1.

Bulk Li Isotope Analyses With Solution Nebulization MC-ICP-MS
Two rock samples from Roche Sauterre (St6 and St3) were ground to a fine powder in order to perform whole rock δ 7 Li analyses. Sample dissolution and chromatographic Li purification was conducted as described in detail in Tomascak et al. (1999a), Bouman et al. (2004), and Seitz et al. (2004). Analyses of the purified Li fraction were conducted on a Thermo-Scientific Neptune Plus MC-ICP-MS for the simultaneous measurement of 6 Li and 7 Li. Following Seitz et al. (2004), a Cetac Aridus II desolvation unit equipped with a pneumatic nebulizer with an uptake rate of ca. 50 µl/min fitted into a PFA spray chamber was used. Sample analysis was performed sequentially, applying sample-standard-bracketing with the IRMM-16 Li reference solution. The reference material JB-2 (basalt powder, Geological survey of Japan) was measured relative to IRMM-16 with 3.2 ± 1.2 (2 SD) which agrees within uncertainties to the compiled literature values of 3.5-4.9 (Brant et al., 2012;Dellinger et al., 2014;Coogan et al., 2017).

Lithium Isotope Analyses With fs-LA-MC-ICP-MS
For in situ Li isotope analyses, a femtosecond laser ablation system (Spectra-Physics Solstice) was coupled to a multi collector inductively coupled plasma mass spectrometer (MC-ICP-MS, Thermo-Scientific Neptune Plus). The ablation beam has a pulse duration of ∼100 fs and a wavelength of 194 nm which is generated via frequency conversion from an infrared beam with 775 nm wavelength in a mirror and lens system and focused on the sample surface via a modified in-house built New Wave (ESI) stage combined with an optical microscope (Horn et al., 2006;Horn and von Blanckenburg, 2007). The laser spot size of ∼26 µm on the bracketing standard GOR132-G (Gorgona Island komatiite) allows for sufficient spatial resolution to resolve diffusion profiles. The protocol for Li isotope measurement was applied as described in detail by Steinmann et al. (2019). In brief, measurements were performed at relatively cool plasma conditions (900 W) in order to avoid matrix ionization in the plasma. In situ Li isotope ratio measurements were performed in static mode at low mass resolution which is sufficient to resolve atomic interferences. In order to keep background Li signals low, the measurements were performed under dry plasma conditions. The sample was mixed in a homogenization device to increase signal stability . For the detection of 7 Li a 10 13 amplifier coupled to a faraday cup was deployed, while a secondary ion multiplier was used for the detection of the less abundant 6 Li. Due to the slower signal response of the 10 13 amplifiers (as compared to 10 11 amplifiers) a tau correction (Kimura et al., 2016) was applied during data reduction. All measurements were performed using standardsample-bracketing with the komatiitic MPI DING reference glass GOR132-G as bracketing standard according to Eq. (1): Subsequently, all data were converted relative to the internationally used Li isotope standard IRMM-16 (with δ 7 Li IRMM−16 = δ 7 Li GOR132−G + 8.6; Jochum et al., 2006;Steinmann et al., 2019). Measurements of the bracketing standard GOR132-G were performed in raster ablation mode with a scan speed of 20 µm/s; the olivine profiles were measured in line ablation mode with lines arranged parallel to the crystal rim so that each line accounts for one measured δ 7 Li-value. Individual measurements consist of 180 cycles, each with an integration time of 1.049 s. The first ∼35 cycles were used for background correction, measuring only the gas blank without a laser ablation signal. This was followed by ∼130 cycles of sample ablation. As internal control for the accuracy of the method, the MPI DING reference glass T1-G [δ 7 Li = 1.6-2.4 , (Jochum et al., 2006;Le Roux, 2010;Xu et al., 2013)] was measured and yielded a long term value of 2.1 (δ 7 Li = 0.4 , 2 SD for n = 64 in 16 sessions throughout 22 months) in agreement with Steinmann et al. (2019). All individual data points along the profiles were conducted as ∼100 µm long line scans parallel to the crystal rim [similar to trace element and Fe-and Mg isotope analyses   (Figure 1)]. This was necessary to detect a sufficient number of Li ions in order to satisfy counting statistics due to the low concentration of Li in the samples, and in order to achieve a stable signal over the duration of the measurement. One crystal (St6-2b_ol4b, displayed in Figure 1) was analyzed across the whole crystal; all other crystals were analyzed from rim to core.

Diffusion Models
Diffusion has been modeled in several studies by applying an analytical or numerical approach to solve Fick's second law of diffusion with the aim of equilibrating an initially homogeneous olivine with new boundary conditions, e.g., by changing the conditions and concentration in the surrounding FIGURE 1 | Photomicrograph (reflected light) of sample St6-2b_ol4b showing the positions for major and trace element and isotope analyses. Black dots are from laser ablation for trace element analyses (indicated by the red arrow), which were conducted exactly along the profile previously analyzed for major elements with EMPA (not visible). Short black lines are from laser ablation for isotope analyses, which was conducted from both rims toward the core of the olivine grain (indicated by blue arrows).
melt (Costa et al., 2008). Most of the investigated olivine crystals have an elongated shape and we always analyzed the short distance to the core of those crystals (i.e., vertically to the long side, see Figure 1). For this purpose, a one-dimensional expression for the diffusion equation in a plane sheet (rather than a spherical) geometry in Eq.
(2) has been solved for diffusion modeling with analytical and numerical calculations (see Supplementary Material for equations).
This equation is valid for diffusive processes, which are not dependent on concentration, which was used for modeling Li diffusion in this study. The diffusivity of Li in olivine has been investigated in several studies following different approaches, either based on laboratory experiments at defined conditions or based on the analyses of natural samples.
(1) An experimental study performed by Dohmen et al. (2010) found two diffusion mechanisms for Li in olivine: one "fast" diffusion mechanism, operating at interstitial sites in the olivine structure and one "slow, " metal vacancy-controlled diffusion mechanism, operating at octahedral sites. The existence of these mechanisms was later confirmed by Richter et al. (2017). In natural systems, the slower mechanism, Li diffusion via octahedral sites, is assumed to dominate (Dohmen et al., 2010). From their experiments, conducted at 800 to 1200 • C and at controlled f O 2 (≈ WM buffer), Dohmen et al. (2010) have developed an equation from an approximated Arrhenius relation to determine the diffusivity of the "slow" diffusion mechanism in olivine as a function of temperature: log (D Li ) = −5.92 (±1.0) − 1.2847 × 10 4 /T(K) The resulting diffusion coefficients are still about an order of magnitude faster than those of most divalent cations, such as Mg, Fe, or Ni.
(2) Several other studies investigated the diffusivity of Li relative to that of other cations and observed an interdependence of the diffusivities of e.g., Fe, Mg, Ni, Mn, Li, and other trace elements (Qian et al., 2010;Spandler and O'Neill, 2010;Oeser-Rabe, 2015). The diffusivity of Li was determined relative to that of Mg-Fe exchange diffusion as given in Eq. (4). The factor denotes the product of diffusion coefficient and duration which equals the diffusive flux ( = D·t; e.g., Ganguly, 2002). For this purpose the concentration profiles of Li and Mg# (recalculated from MgO and FeO concentrations) are analyzed and Mg# is fitted with Eqs. (S.1) and (S.2). The concentration profile of Li is fitted by adjusting , which gives the relative diffusivity of Li depending on Mg-Fe. Eq. (4) is valid under the assumption that the diffusion time t is the same or very similar for all elements.
(3) We applied an approach similar to approach (2) for the olivine investigated in this study. As all olivine crystals, analyzed for their Li isotope composition in this study, were previously analyzed for their major and trace element compositions, including Fe, Mg, and Li, as well as for their Fe-Mg isotope compositions, we were able to determine Li diffusivity data relative to those of the well-determined Fe-Mg exchange diffusion. In detail, we used the timescales, as determined by Oeser et al. (2015), by modeling chemical and isotopic Mg and Fe inter-diffusion, which were interpreted to represent the residence of the olivine crystals in a magma between crystal formation for phenocrysts and the entrainment of xenocrysts in the host magma and eruption of the magma. Based on these time scales (fixed parameter t), the Li concentration and isotope profiles were fitted by adjusting D Li in order to determine the diffusivity of Li. Thus, the resulting Li diffusivity is based on modeling of two Li profiles, i.e., the concentration and isotope profiles. To be consistent with modeling by Oeser et al. (2015) for the same samples we applied constant temperature conditions for the system. We decided to assume a temperature of 1250 • C, as this was shown to be best suited for the investigated samples . We are aware that a rim Mg# of 0.8 indicates lower temperatures around 1200 • C (Roeder and Emslie, 1970). However, Costa et al. (2008) demonstrated that modeling with an integrated cooling rate does not affect the modeled time scales significantly as most of the diffusive flux occurs in the highest temperature regime of the system. In contrast to Fe-Mg diffusion, Li diffusion was generally modeled from both sides of the crystal, even if the profile was not measured as a complete transect. This was necessary because, due to the relatively fast diffusive character of Li, the crystal core composition was affected by diffusion from both of the opposite crystal rims, in several cases. For all crystals measured as half transects, boundary conditions for the unknown crystal rim are assumed to be in accordance with those of the measured rim. For St6-2b_ol4b, a complete transect was measured, and hence the measured values were applied in the model.
From the samples described above, olivine crystals with high Li concentrations and with the largest variations in Li concentrations from rim to core were selected for Li isotope analysis and modeling of profiles. For approach (1) the diffusivity is calculated according to Eq. (3), as the diffusivity depends on the temperature T. For approach (2) an analytical solution [Eq. (S.1) and Eq. (S.2) in Supplementary Table S1] of Eq. (2) is applied and for approach (3) a numerical solution by the method of finite differences of the diffusion equation in Eq.
The relation between the diffusion coefficients of the two isotopes of Li can be described by the empirical formula of Richter et al. (1999): Where D values are the diffusion coefficients of 6 Li and 7 Li, m values are the atomic masses of 6 Li and 7 Li in Dalton and the β-value describes the diffusion-driven isotope fractionation. β is an empirical constant which depends on the diffusion medium and is smaller than 0.5 (Richter et al., 1999). A β-value of 0.5 is the factor for a non-uniform gas (Chapman and Cowling, 1953), and the maximum value achievable. The β-value for Li has been determined experimentally and can be fitted by values of β = 0.4 ± 0.1 for the crystallographic a-and c-axis whereas the fractionation along the b-axis appears to be slightly lower (Richter et al., 2017). Transferring the relative difference in the diffusivities of the two Li isotopes, determined by Dohmen et al. (2010) to a β-value, a β Li of 0.19 may be assumed for the slow diffusion mechanism ( 6 Li diffuses 3% faster than 7 Li) and of 0.32 for the fast diffusion mechanism ( 6 Li is 5% faster than 7 Li). Compared to Richter et al. (2017), these values are slightly lower but comparable to Li diffusivity in molten oxides (Richter et al., 2003).

RESULTS
Four olivines from Banne D'Ordanche (BdOr-1_ol1, BdOr-1_ol2, BdOr-1_olxen1, and BdOr-1_olxen2) were analyzed regarding their Li isotope composition. Major-and trace element concentration data were taken from Oeser-Rabe (2015). The investigated olivine grains display Li concentrations varying from ∼8 to 9 µg/g at their rims to 1-2 µg/g in their core regions. The lava lake of Roche Sauterre was been sampled at two positions. Two olivines from one sample (St6-2b_ol4b and St6-2b_ol5) display concentrations varying from 7 µg/g at their rims to 3 µg/g in their cores. Two olivines in the other sample (St3-3a_ol1 and St3_3a_ol2) have somewhat lower Li concentrations ranging from 4 µg/g at their rims to 1 µg/g in their cores. All investigated grains display the lowest Li concentrations in their cores. Nevertheless, Li concentration profiles show different shapes (Figures 2A-C), i.e., some of them show a slight decrease of Li toward the very rim. Lithium concentrations and isotope compositions for all investigated samples are compiled in Table 1.
Lithium isotope profiles of all measured olivines display elevated δ 7 Li-values at the crystal rims and significantly   Oeser et al. (2015) and Oeser-Rabe (2015).
lower values in the core regions (Figure 2) with a Li isotope fractionation between rim and core ranging from 7 Li rim−core = 21.5 to 27.3 for Roche Sauterre and 7 Li rim−core = 13.5 to 38 for Banne D'Ordanche. In most cases 7 Li rim−core is equal-or close to the largest Li isotope fractionation 7 Li min−max observed for each grain, which indicates that the cores are strongly affected by diffusion and have not preserved their original Li isotopic composition. The δ 7 Li-values in the cores (−9.6 to −14.5 and −30.7 for one crystal) are clearly lower than δ 7 Li-values of ∼4 of volcanic bulk rock (Seitz et al., 2004;Magna et al., 2006;Jeffcoate et al., 2007). Similar low δ 7 Li-δ 7 Li values in natural olivines were observed for other olivines from the Massif Central (Gu et al., 2016) and for San Carlos olivine xenoliths (Jeffcoate et al., 2007). Lithium concentration and isotope profiles are broadly correlated, with high Li concentration and high δ 7 Li-values at the rim and decreasing Li concentrations and low δ 7 Li-values toward the core (Figure 2). The shape of the Li isotope profiles is consistent with the shape that would be expected from diffusiondriven isotope fractionation (Parkinson et al., 2007). Iron and magnesium isotope profiles, displaying diffusion-driven zoning of δ 56 Fe and δ 26 Mg are shown for comparison in Figures 2D-F (taken from Oeser et al., 2015). The width of Li isotopic zoning slightly exceeds that of Fe and Mg isotopic zoning, in most cases. All measured Li concentration and isotope ratios (this study) are listed in Supplementary Tables S2, S3 together with the Mg and Fe isotope compositions determined by Oeser et al. (2015).

DISCUSSION AND MODELING OF DIFFUSION EVENTS Boundary Conditions of the System
Coupled Li concentration and isotope profiles from the rims of the crystals toward their cores, and the extremely light isotopic composition of the olivine cores (δ 7 Li min -values down to −30.7 , Table 1) strongly indicate a diffusive origin of the zoning caused by faster diffusion of 6 Li into the crystal compared to 7 Li. The coupled zoning profiles for Li concentrations and isotope compositions are similar, though not identical in scale and magnitude, to previously investigated Fe concentration and isotope profiles (and anti-correlated Mg concentration and isotope profiles), which were generated by Fe-Mg exchange diffusion ; Figure 2). The Li concentration profile shapes observed for olivine from the Massif Central can be distinguished into three types (indicated by the green lines in Figures 2A-C, which represent best fits through the data), with always higher Li concentrations in the rims and lower concentrations in the cores. Type I (Figure 2A) shows decreasing Li concentrations toward the cores which is a typical profile for diffusion into the crystals. Type II ( Figure 2B) consists of two compositional plateaux with constant Li concentrations at outermost rims and cores and gradually decreasing Li concentrations in transition zones. Type III (Figure 2C) exhibits maximum Li concentrations at some distance from rims (∼100 µm) and then decreases of Li concentrations toward cores. These latter profiles may indicate two diffusion events with a first event of Li diffusion into the crystals and a second event causing Li diffusion out of the crystals. While the differences between the three types of Li concentration profiles are striking, the corresponding isotope profiles cannot be distinguished as clearly.
Before diffusion models can be established, the boundary conditions of the system at the start of diffusion have to be evaluated (Costa et al., 2008). These include the Li concentration (profile) before diffusion started and here also the initial isotope composition (or isotope profile). Here, we assume that the crystals were initially homogeneous and grew in equilibrium with the melt. Assuming that the melt initially had a similar Li concentration to that analyzed for whole rock samples [BdOr1 = 6.14 µg/g, St3-3a = 5.7 µg/g, and St6-2b = 5.78 µg/g ] and applying the Li olivine/melt partition coefficient determined by Brenan et al. (1998) for olivine with high Mg# (of ∼0.2), this results in equilibrium concentrations for olivine of 1.16 to 1.47 µg/g (see also Supplementary Figure S2). Lithium concentrations in the olivine cores range from 1 to 2 µg/g for the two olivine crystals from sample St3a from Roche Sauterre, as well as for three out of four Banne D'Ordanche olivines. The other two Roche Sauterre crystals (from St6) exhibit higher Li concentrations of 3.1 µg/g and 3.9 µg/g ( Table 1). Only one crystal from Banne D'Ordanche displays a Li core concentration of 8.7 µg/g, likely indicating significant Li diffusion into this core. Two of the crystals are found in a glomerocryst and hence are described as xenocrysts by Oeser et al. (2015). The initial concentration of such crystals is difficult to constrain. Notably, the rim Li concentration of most olivine grains is on the order of 4 to 8 µg/g and difficult to explain with olivine/melt partition coefficient on the order of 0.2 (Brenan et al., 1998), as this would indicate extremely high Li concentrations (on the order of 20 to 40 µg/g) in the melt at the time when diffusion ceased, i.e., at the end of fractional crystallization. Such high Li concentrations were not observed for any of the investigated basalts . However, although Brenan et al. (1998) suggested low olivine/melt partition coefficient on the order of ∼0.2, those were essentially determined for extremely forsterite-rich olivine (Mg# = 99 to 100), while a higher olivine/melt partition coefficient (of ∼0.35) was observed for an olivine with a slightly lower Mg# (of 92). Other authors have observed higher partition coefficients, between 0.21 and 0.56 (Taura et al., 1998;McDade et al., 2003;Ottolini et al., 2009), which are always in combination with lower Mg# than the forsterite-rich olivines form Brenan et al. (1998). Ottolini et al. (2009) have interpreted, that there might be an influence of the melt composition on the partition coefficient, and Taura et al. (1998) supposed a pressure dependence. Thus, Li partitioning into olivine potentially increases with decreasing Mg#. If so, the high Li contents, observed in the olivine rims, may primarily be driven by increasing Li partitioning into olivine during fractional crystallization, rather than by increasing Li contents in the melt, though Li concentrations in the melt at the end of fractional crystallization may also have been slightly higher than those observed in the sampled basalts.
Similar to the Li concentrations, the measured δ 7 Li bulk of the basalts may be assumed as the initial Li isotope composition of the olivine crystals before diffusion started. The isotopic composition of fresh basalts varies from +2 to +5 and does not appear to change significantly during fractional crystallization (Chan et al., 1992;Elliott et al., 2006;Tomascak et al., 2008). The degassing of a melt causes a transition of Li from the melt or crystal to a vapor phase and a decrease in δ 7 Li to lower values (Vlastélic et al., 2011). For Roche Sauterre St6 and St3 a δ 7 Li bulk of 3.3 ± 1.3 and 2.3 ± 2.0 were determined, respectively. For Banne D'Ordanche, no δ 7 Li bulk is available and a mean basaltic δ 7 Li value of 3.5 was assumed as the initial value for the modeling. Furthermore, the initial isotopic composition of the crystals is assumed to be homogeneous and the isotopic composition of the melt to be constant during fractional crystallization (e.g., Tomascak et al., 1999b). The initial Li isotope composition of the potential xenocrysts is uncertain and may as well be different from that of basalts. Mantle olivine may have heterogeneous δ 7 Li, which is ascribed to metasomatic overprinting (e.g., Nishio et al., 2004;Tang et al., 2012;Ackerman et al., 2013;Su et al., 2016) or diffusion in the mantle (e.g., Marschall et al., 2007a,b;Magna et al., 2008).

The First Diffusion Event: Constraints on Li Diffusion Coefficients and β Values
For the model applied here, an open system is assumed due to the quasi-infinite supply of Li from the surrounding melt. This quasi-infinite supply of Li also justifies the use of a fixed boundary concentration. As outlined above, a variety of Li diffusion coefficients, determined with different approaches, are available in the literature for modeling the observed chemical and isotopic zoning of Li. As will be shown below, the choice of Li diffusion coefficients is critical for the determination of diffusion timescales. We considered the Mg-Fe diffusion time determined by Oeser et al. (2015) as a fixed diffusion event to determine Li diffusivity. For this purpose we applied a combined approach of chemical and isotopic modeling of Mg-Fe exchange diffusion, because the diffusivity of Li in olivine is less well constrained than that of Mg-Fe . This event was interpreted as the residence time of the olivine crystals in a magma reservoir. Iron-Mg exchange diffusion was driven by the chemical evolution of the magma, as a result of fractional crystallization, which generated a chemical gradient between the melt and previously formed olivine crystals . This chemical evolution would have also generated a chemical gradient in Li between olivine and surrounding melt (and as described above, likely also a change in Li partitioning), forcing Li diffusion. This scenario is consistent with increasing Li concentrations toward olivine rims and decreasing δ 7 Li values toward the olivine cores, indicating diffusion of Li (and preferentially of isotopically light Li) into olivine, as would be expected during magma differentiation. We therefore assume that the light Li and Fe isotope-(and heavy Mg isotope) signatures were generated by the same event.
Based on this assumption, Li diffusion into olivine was modeled using the timescales given by Oeser et al. (2015) for the same samples (for BdOr1_ol2 = 239 days, for BdOr1_olxen2 = 252 days and for St3-3a_ol1 = 3.26 years) and fitting the measured Li concentration and isotope profiles. The initial concentration calculated above was assumed as the Li core concentration (C core ) and in some cases it was adapted to fit the measured profile. In the first modeled event a Li rim concentration (C rim1 ) of 3.5 to 21 µg/g was applied in order to match the slope of the measured profiles (Table 2). Notably, all crystals were modeled assuming equivalent diffusion and C rim from both crystal rims toward the core. The only exception was St6-2b_ol 4b, because the crystal was measured as a whole transect, and hence both measured rim concentrations were modeled. As mentioned above, the investigated olivines display three different types of Li zoning (Figure 2) and the shapes with constant or decreasing Li concentrations at the very rim (type II and III, Figures 2B,C) cannot have formed during a single diffusion event, but rather indicate a second diffusion event (with Li diffusion in opposite direction). For those olivines, including BdOr1_ol2, BdOr1_olxen1, St3-3a_ol1, and St6-2b_ol4b, the Li concentration and isotope composition measured at the very rim was excluded when modeling the olivine core composition and the slope of the core-rim transition during the first diffusion event. As will be shown below, the very rim composition of such olivines can be modeled with a second diffusion event.
Our results show that the initial concentration in crystal cores (C core ) changed during the course of diffusion for several grains, as a result of the faster diffusion of Li, as compared to that of Fe and Mg. Assuming a fixed diffusion timescale from Fe-Mg diffusion, the diffusion coefficient D Li was fitted to obtain congruent modeled and measured zoning profiles, for both, the Li concentration and Li isotope profiles, with a particular focus on Li core compositions. Another important variable for modeling Note that the duration of the first diffusion event was fix (taken from Oeser et al., 2015) and D Li was fitted, whereas for the second diffusion event, the duration, t step2 was a model outcome, β step1 is the modeled value from the first diffusion event and β step2 is the value determined from the second diffusion event. For comparison the experimentally determined β values of the slow diffusion process from Dohmen et al. (2010) and Richter et al. (2017) are 0.19 and 0.4 ± 0.1, respectively. * BdOr1-ol1 was modeled in a 3-step model with an additional step 0 before step 1 with C rim0 = 5 µg/g, β step0 = 0.25 and t step0 = 2760 days. + t step1 durations were taken from Oeser et al. (2015). # St6-2b_ol4b was measured from rim to rim, for modeling C rim was adjusted relative to C measured .
the isotopic profiles of Li is the β-value, which describes isotope fractionation during diffusion. It also has to be fitted in order to match the amplitude of the measured isotope profiles. β-values were determined in this study by fitting the modeled isotope profiles to the measured ones (as done for Fe and Mg in previous studies, e.g., Sio et al., 2013;Oeser et al., 2015 and others). With this approach, β-values determined in this study do not exceed 0.25 which is lower than the β-values of 0.4 ± 0.1 determined by Richter et al. (2017) for the crystallographic a-and c-axes, (for the b-axis somewhat lower β-values were determined), and also lower than that those that have been modeled by the same authors for a natural olivine (β = 0.30-0.36). The β-values obtained here indicate a 1.6% (β = 0.1) to 3.9% (β = 0.25) higher relative diffusivity of 6 Li as compared to 7 Li. This is similar to the results of Dohmen et al. (2010), who determined a ∼3% faster diffusivity of 6 Li relative to 7 Li in their experiments, which corresponds to β = 0.19. These findings imply, that β-values for Li may depend on intensive thermodynamic variables which have not yet been explored in detail. As the zoning profiles of several of the investigated olivines indicate a second diffusion event (see below), a slightly lower β value (than required for a perfect match of the core composition in only one diffusion step) was generally assumed for modeling of the first event, in order to match the measured isotopic zoning with a combined model of both diffusion events. Thus, the final β-value for the first diffusion event could only be determined after modeling both diffusion events. This is described in more detail below. Notably, assuming that some of the olivine grains (type I, Figure 2A) were unaffected by the second diffusion event would result in a slight underestimation of the β value for those samples. For example, olivine grain BdOr1-olxen2 was modeled with both one and two diffusion events. This revealed β-values of 0.16 for one diffusion event (Figure 3) and 0.13 for the first of two diffusion events (Supplementary Figure S3F), respectively. However, the assumption of a second diffusion event had no effect on the modeled values for D fit . All modeled D Li and β-values are listed in Table 2.
Diffusivities of Li relative to the diffusion couple Mg-Fe were obtained from the equations for the analytical solution in Supplementary Table S1 and yield D Li /D Mg−Fe ratios of ∼2 for BdOr1, which means the diffusivity of Li is twice as fast as the diffusivity of Mg-Fe in BdOr1. For St3-3a and St6-2b the diffusivity of Li is ∼1.5 times faster than the diffusivity of Mg-Fe (Figure 4). These results are similar to previous findings by Qian et al. (2010) and Spandler and O'Neill (2010). The parametrization for the determination of the diffusion coefficient of Mg-Fe diffusion in olivine by  depends on temperature (T), pressure (P), oxygen fugacity (f O 2 ), the mole fraction of the fayalite component (X Fe ), and the crystallographic orientation of the measured profile. With Eq. (4) and the relative diffusivities of Mg-Fe and Li and the absolute calculated diffusion coefficients for Mg-Fe (D Mg−Fe ) the diffusion coefficient for Li (D Li ) is obtained (Table 3).
The effect of using different Li diffusion coefficients was evaluated, by comparatively modeling Li diffusion with the experimentally determined diffusion coefficient (D exp ) (for this comparison, we used the slower diffusion mechanism of Dohmen et al., 2010) and the fitted diffusion coefficient from this study (D fit ). In Figure 3 the effect of both models on the shape of the chemical and isotopic diffusion profiles is compared (for sample BdOr1_olxen2). For the comparison in Figure 3, a crystal was selected, for which the effect of a potential second diffusion event was negligible. The experimentally determined diffusion coefficient for Li from Dohmen et al. (2010) was adjusted to a temperature of 1250 • C , according to Eq. (3), although, experimental diffusion coefficients were only determined in the temperature range between 800 • C and 1200 • C. This comparison reveals that applying the experimentally determined diffusivity of Li, the modeled Li concentrations and isotope compositions do not match the observed profiles when assuming a diffusion time from Oeser et al. (2015) (Figures 4A,B). In turn, applying D exp and fitting the modeled Li and δ 7 Li to the measured profile would indicate a significantly shorter diffusion time scale than that determined by modeling Mg-Fe exchange diffusion (Figures 3C,D). This would imply, however, that the development of Li concentration gradients was entirely decoupled from the development of Fe-Mg concentration gradients for all investigated samples, i.e., was generated by a much later event. Exchange diffusion between Fe and Mg was likely driven by the chemical differentiation of the  (Dohmen et al., 2010) and fitting the time scale to match the Li concentration and isotope profiles yields a diffusion time of t = 22 days (blue line); for comparison the profile is also modeled with the same time applying the fitted diffusivity from (A,B) (red line). All fits were modeled assuming a diffusivity D exp = 4.42*10 -15 m 2 /s according to a temperature of 1250 • C and D fit = 3.5*10 -16 m 2 /s. FIGURE 4 | Relative modeling of Mg# and Li (µg/g). The Mg# and Li profiles are fitted by adjusting Ã Mg-Fe , and Ã Li , respectively. The ratio of the Ã-values is calculated giving the relative diffusivity of Li relative to Mg-Fe, various curves for the Li profile show the variability and uncertainty in the modeled factor for relative diffusities. (A) BdOr1-olxen2 is modeled with an infinite analytical approach assuming a fixed rim concentration. The modeled profile for Li yields a factor for the Li diffusivity relative to that of the Fe-Mg exchange (indicated by Mg#) Ã Li /Ã Mg-Fe from 1 to 3, which may be considered as the uncertainty of this approach. (B) St3-3a_ol1 is modeled with a semi-infinite analytical approach assuming the equilibration of two chemically distinct layers of the crystal over time, yielding a relative diffusivity factor within an uncertainty of Ã Li /Ã Mg-Fe from 1 to 2. melt in the magma chamber. Thus, this scenario would imply that both Li concentration in the melt, as well as Li partitioning between olivine and melt, which itself is dependent on Mg# of the olivine, remained about constant for most of the time and then only changed at the end of magma differentiation. In contrast, assuming the time scales for the individual samples, given by modeling of Fe-Mg diffusion  and fitting the Li diffusivity to the measured profiles, a narrow range of D fit does acceptably fit to all analyzed Li concentration-and Li isotope profiles. We therefore conclude that diffusion timescales for Li were likely similar to those for Fe-Mg and that Li diffusivities in the investigated olivines are likely more similar to those suggested by Qian et al. (2010); Spandler and O'Neill (2010) and Oeser-Rabe (2015), rather than to those determined by Dohmen et al. (2010) (here labeled as D exp ). These findings demonstrate that additional experimental investigations on Li diffusivity at lower concentration levels are urgently demanded.
Differences between experimentally determined and fitted diffusion coefficients may arise from differences in the initial and boundary conditions between those imposed in experiments and those assumed for the investigated natural samples, e.g., Li concentration in the melt (or in the powder of the experiments) surrounding the olivine crystals. In the investigated samples, the melt Li concentrations were estimated (from bulk rocks and model fits), between 5.7 and 6.14 µg/g. In contrast, the powder mixture of 90% San Carlos olivine and Li 2 SiO 3 , used in the experimental setup of Dohmen et al. (2010), had significantly higher Li concentrations. Notably, it was observed by Dohmen et al. (2010) that Li diffusivity in olivine depends on the Li concentration and appears to decrease with decreasing Li concentrations.
Furthermore, other intensive thermodynamic variables of the system, including the chemical composition of the crystals and melt and f O 2 , which may not have yet been investigated in sufficient detail, may affect Li diffusion in olivine. In a recent experimental study, Richter et al. (2014) demonstrated a strong sensitivity of Li diffusion in clinopyroxene to f O 2 . Also the experiments of Dohmen et al. (2010) imply that the diffusion mechanism of Li in olivine depends on environmental variables, such as f O 2 and aSiO 2 (amongst other factors). These parameters affect the vacancy concentrations in olivine which control the diffusion mechanism (Dohmen et al., 2010). For a series of trace elements in olivine, including Al, Cr and Be, several studies observed a dependence of the diffusivity on factors, such as silica activity (aSiO 2 ) and oxygen fugacity (f O 2 ) (Jollands et al., 2016a,b;Zhukova et al., 2017). The diffusivity of Al in forsterite at low aSiO 2 for example is up to three orders of magnitude slower than at high aSiO 2 (Zhukova et al., 2017). These studies show that diffusion in olivine may depends on parameters involving the activity of SiO 2 and f O 2 .

The Second Diffusion Event
While our models broadly reproduce the low Li isotope compositions observed in olivine cores, modeled Li concentrations do not match the rim areas of several of the olivine grains investigated (e.g., BdOr1_ol2; BdOr1_olxen1; St3-3a_ol1). Similarly, the Li isotope composition does not match the observed δ 7 Li profiles for these olivine grains. Rather, the change in slope of the Li concentration profile in Figure 2C indicates that the flux of Li into the crystal (assumed in the first diffusion event) changed to a flux of Li out of the crystal during a second diffusion event at a later stage (Figure 5). The latter is not visible in all Li concentration profiles and potentially has not affected all olivine crystals to the same extent.
Such a second diffusion step was modeled using the Li diffusivity that was determined by modeling the first diffusion event, but assuming a different Li rim concentration (C rim2 ; Figure 5). C rim2 was estimated based on the bulk rock Li concentration of the samples [5.78 µg/g for St6-2b, 5.7 µg/g for St3-3a and 6.14 µg/g for BdOr-1; Oeser et al. (2015)] and slightly adapted to better fit the observed Li concentrations ( Table 2). This second diffusion event was modeled by fitting the Li isotope and concentration profiles to those observed at the crystal rims. As the diffusivity of Li was assumed to be given by the results of the first diffusion step, the modeled profiles were fitted by adjusting the duration of diffusion, which was unknown for this second event. Applying such a model, yields durations of 20 to 80 days for Banne D'Ordanche, 100 to 200 days for St3-3a from Roche Sauterre and 30 to 120 days for St6-2b from Roche Sauterre for this event (Figure 6). For the second modeling step the β-value was readjusted in order to fit the amplitude of the Li isotope diffusion profile. As described above, modeling of the second diffusion event also has an effect on the amplitude of the isotope profiles modeled in the first diffusion event. Best fit modeling of the olivine core-(first diffusion event) and rim (second diffusion event) isotope compositions furthermore, resulted in different β-values for the first and second diffusion event ( Table 2).
The second diffusion event of Li diffusion out of the olivine crystals indicates a decrease of the Li concentration in the melt, which may be associated to the ascent of the melt, or even have occurred after the emplacement of the lava. During the exsolution of H 2 O-rich vapor in the course of decompression Li is largely mobilized in magmatic systems (Sakuyama and Kushiro, 1979). For example, Neukampf et al. (2019) observed a significant decrease in Li concentrations in plagioclase, which they interpreted to be driven by a Li loss in the melt during degassing of the magma. The low Li concentrations in the melt caused Li diffusion out of the crystals. Similarly (though to a lower extent than observed by Neukampf et al., 2019), decreasing or constant Li concentrations at the olivine rims, observed in this study, may be related to degassing-driven Li loss in the melt during magma ascent or emplacement, resulting in Li diffusion out of the olivine crystals. Notably, the second event, which operated in an opposite direction to the first diffusion event, is not visible in the Mg-Fe chemical or isotopic zoning. This finding may support degassing as the driving force for diffusion during the final melt evolution, which would have affected fluid mobile elements, such as Li, but not refractory elements, such as Fe and Mg. Hence, though Li and Fe-Mg diffusion were likely coupled during the first diffusion event, they were apparently decoupled during the FIGURE 5 | (A1) and (B1) show the Li concentration and isotopic composition of the first diffusion event for BdOr1_ol2 with a duration of t = 239 days, C core = 1.1 µg/g, C rim1 = 12 µg/g and β = 0.25. (A2) and (B2) show the model for the second diffusion event (which is based on the first diffusion event) with t = 40 days, C rim2 = 6 µg/g and β = 0.2, the Li diffusion coefficient was determined by fitting the observed Li concentrations and isotope compositions to 5.5*10 -16 m 2 /s. Core Li concentrations and isotope compositions are marked by a dashed line. second event. This final decoupling was likely the result of the different physiochemical properties and diffusivity of Li, as compared to Fe-Mg. Accordingly, Li concentrations and isotopes may unravel processes which cannot be seen with Mg-Fe isotopes alone.
Some profiles are not matched perfectly by modeled curves (e.g., Figure 6D and Supplementary Figure S3B) which may be the result of cutting effects during the preparation of the sections. That is some sections probably do not sample the core of the crystals, which would result in an underestimation of the diffusive mass flux and isotope fractionation (Weyrauch et al., 2019). Furthermore, some sections were probably not cut perpendicularly to crystal rims, which may be indicated by asymmetric zoning (e.g., crystal St6-2b1_ol4) in the case of rim-to-rim profiles. This effect would result in artificially elongated Li diffusion profiles. For Mg and Fe diffusion modeling, the orientation of the crystallographic axes was determined by Oeser et al. (2015) by electron backscatter diffraction analyses. This is necessary due to the difference in diffusivities along the crystallographic axes . However, crystallographic orientation does not seem to have a large effect on Li diffusivity in olivine (Dohmen et al., 2010) and thus was not considered.
The time scales of the second diffusion event determined by Li diffusion modeling are 20-80 days for Banne D'Ordanche and 30-200 days for Roche Sauterre which is longer than the time scales obtained by Lynn et al. (2018) for the Keanakâko'i Tephra (Kîlauea Volcano, Hawai'i), also using Li diffusion modeling (only Li concentrations). The authors conclude a similar process of syn-eruptive degassing in their study with time scales of only hours to days which is 1-2 orders of magnitude faster than that estimated with the fitted D Li values of this study. The difference in the determined time scales of the degassing process is essentially the result of the different diffusion coefficients used in the studies. In this study, the same diffusion coefficient, as determined from the first diffusion event ( Table 3, derived from the comparison of Li-with Fe-Mg diffusion profiles) has been assumed for both diffusion events. However, as the diffusion coefficient may depend on temperature, f O 2 and the fluid content of the melt (Coogan et al., 2005; it may have changed during the second diffusion event, at which the melt may have become water saturated during degassing. Accordingly, it cannot be excluded that Li diffusion was faster during the second diffusion event and the modeled time scales overestimate the duration of a potential degassing process. Notably, however, the basalts investigated in this study, originate from massive flows, or in the case of Roche Sauterre, even from a small lava lake (Nehlig et al., 2001;Lorand et al., 2003;Richet, 2003). Thus degassing, as recorded by Li diffusion,  (1) relative to the diffusivity of Mg-Fe exchange diffusion (D Lirelative ), (2) by fitting D Li to the observed Li concentration and isotope profiles assuming the time scales from Fe-Mg exchange diffusion modeling , and (3) the experimentally determined D Li and D Mg−Fe (Dohmen et al., , 2010 calculated for 1250 • C.  *Calculated at P = 500 Pa, T = 1250 • C, fO 2 = NNO-2, crystallographic orientation, X Fe and concentration profiles from Oeser et al. (2015). (D Li /D Mg−Fe ) relative is calculated by multiplying the actual diffusivity of Mg-Fe at 1250 • C and recalculated to the crystallographic orientation of the profile, with Ã L i/Ã Mg−Fe dcalculated with approach (2) according to Eq. (4). # Experimentally determined D Li calculated at 1250 • C according to Eq (3) from Dohmen et al. (2010).
may have occurred mostly during cooling of the magma after lava emplacement and thus lasted much longer than timescales typically expected for magma ascent. Assuming the outcrop of Roche Sauterre samples represents a small lava lake, this scenario would explain the longer degassing timescales for samples from this locality as compared to those from Banne D'Ordanche. It may also explain that some samples appear to be more affected by degassing than others, as degassing rates would depend on the position of the samples in the lava flow/lake.

CONCLUSION
Combined in situ Li concentration and isotope analyses were performed on zoned olivines from the Massif Central volcanic region. A diffusive origin for Li concentration and isotope profiles was deducted from correlations with Mg-Fe isotopic inter-diffusion. However, in contrast to Fe and Mg, Li indicates two diffusion episodes. The first likely represents the residence time of olivine in a (deeper) part of the plumbing system, as implied by Mg-Fe isotope diffusion modeling . This diffusion event is likely driven by an increase in Li concentration in the melt and, likely also, by an increased partitioning of Li into olivine, driven by decreasing Mg# in olivine during fractional crystallization. Lithium diffusion into the olivines generated low δ 7 Li-values in the cores of the crystals. A second and shorter diffusion event, recorded by several of the investigated olivine grains resulted in Li diffusion out of olivine crystals, indicating a Li decrease in the melt. This second event is likely related to degassing of the melt during its ascent and after lava emplacement, resulting in the mobilization of Li into a vapor phase. Such a diffusion event can be deduced by modeling Li concentration and Li isotope profiles simultaneously. The second diffusion episode is not recorded in Mg-Fe profiles, which highlights that Li profiles can preserve information about magmatic evolution that is not recorded by Mg-Fe zoning, potentially due to the volatile character of Li in magmatic fluids. These findings demonstrate that combining information from chemical and isotopic zoning of Li with those from Fe-Mg is a very powerful tool for unraveling the timing of complex magma evolution scenarios, such as commonly observed in magma plumbing systems.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
LS contributed to the analyses of Li isotope profiles with femtosecond-laser ablation-MC-ICP-MS, solution nebulization isotope analyses with MC-ICP-MS and proceeding Li column chromatography, and data processing and interpretation, diffusion modeling, and writing the manuscript. MO contributed to the project design, analytical support with fs-LA-MC-ICP-MS, support with diffusion modeling, and data interpretation, support with manuscript writing. IH contributed to the analytical support with fs-LA-MC-ICP-MS and support with manuscript writing. SW contributed to the project design, support with data interpretation, and support with manuscript writing.

ACKNOWLEDGMENTS
This study was supported by the Deutsche Forschungsgemeinschaft (grant # OE 653/1-1) which is gratefully acknowledged. We are grateful for the reviews of TK and IP and editorial handling by DN, whose constructive comments significantly helped to improve this paper.