Accurate Determination of Zr Isotopic Ratio in Zircons by Femtosecond Laser Ablation MC-ICP-MS with “Wet” Plasma Technique

This work evaluates the use of femtosecond laser ablation multiple collector inductively coupled plasma mass spectrometry (fs-LA-MC-ICP-MS) for Zr isotopic analysis in zircons. The mass fractionation caused by instrumental mass discrimination was corrected by a combination of internal correction using Sr as an internal standard (coming from a NIST SRM 987 standard solution) and external correction using a matrix-matched standard. Several important instrument parameters were investigated, such as the effect of the addition of N2 and “wet” plasma condition, the mass fractionation behaviors between Zr isotopes and Sr isotopes, the position effect in laser ablation cell and the effect of laser ablation parameters (laser spot size and energy density). The Zr isotope compositions of seven zircons (GJ-1, 91500, Plešovice, Rak-17, Paki, Aus and Mala) were determined by the developed fs-LA-MC-ICP-MS and thermal ionization mass spectrometry (TIMS). Our fs-LA-MC-ICP-MS results for Zr isotope compositions agreed with TIMS analyses within analytical uncertainties, indicating the presented method is a suitable tool to resolve isotopic zoning in natural zircons. The results also suggest that GJ-1, 91500, Plešovice, Paki, Aus and Mala had the homogenous Zr isotope composition and could be considered as the potential candidates for the Zr isotope analysis in zircons, except Rak-17 which presented the large Zr isotope variation.


INTRODUCTION
Zirconium is classified as a high field strength element (HFSE) in geochemistry, and has been widely used to trace the differentiation in silicate magmatic systems and to understand the co-evolution of the Earth's mantle and crust (Tang and Li, 2021;Niu, 2004;Bau, 1996;Brenan et al., 1994;Woodhead et al., 1993). Recently, there has been an increasing interest in studying mass-dependent isotope variations of Zr in geological samples (Méheut et al., 2021;Guo et al., 2020;Tian et al., 2020b;Ibañez-Mejia and Tissot, 2019;Inglis et al., 2019;Kirkpatrick et al., 2019). Inglis et al. (2018) firstly reported the high precision Zr stable isotopic measurement method and the significant Zr isotopic fractionation (up to 0.57‰) in a suite of co-genetic lavas from Hekla volcano in Iceland spanning bulk compositional variations from basalt to rhyolite. Zircon crystallization was considered as the main cause of zirconium isoto-pic fractionation in the magma system. Because the higher coordination of Zr ions and longer Zr-O bond lengths in zircon relative to the melt would preferentially incorporate lighter Zr isotope composition in zircons. This expectation based on the first principles of stable isotope fractionation theory has been proved by the results from Inglis et al. (2019). However, another Zr isotope research from the FC-1 anorthositic cumulate found zircon and baddeleyite were isotopically heavier to the co-existing melt, which was contrary to the expectations from bonding rule (Ibañez-Mejia and Tissot, 2019). Although there are contradictions in the fractionation mechanism of Zr isotopes, the significant Zr stable isotopic fractionation during the magma differentiation indicated that Zr isotopes had the potential to trace both magmatic differentiation and crustal evolution.
Zircons have been considered as an important material of the Zr isotope research (Guo et al., 2020;Ibañez-Mejia and Tissot, 2019;Inglis et al., 2019;Kirkpatrick et al., 2019;Zhang et al., 2019). Furthermore, the comprehensive information can be provided in zircons to decipher the evolution of the Earth's crust and mantle, including the concentration of trace elements, O-Hf isotope compositions and U-Pb dating. The high precision determination of Zr isotopes can be performed by powder digestion + chemical purification combined with thermal ionization mass spectrometry (TIMS) or multiple collector induc-tively coupled plasma mass spectrometry (MC-ICP-MS) (Feng et al., 2020;Tian et al., 2020a;Inglis et al., 2018). However, these methods generally only report the whole-grain Zr isotopic composition. The microanalysis method allows the interand intra-crystalline isotope variations of mineral grains to be identified and distinguished on the scale of micrometers, such as laser ablation MC-ICP-MS (LA-MC-ICP-MS) and the secondary ion mass spectroscopy (SIMS) (Bao et al., 2021;Tang et al., 2020;Zhang et al., 2018). Recently, the significant Zr isotope zoning has been observed in single zircon grains. The micro-analyses of Zr isotopes are the powerful tools to identify the characteristics of these Zr isotope zoning (Kirkpatrick et al., 2019;Zhang et al., 2019). Kirkpatrick et al. (2019) reported the first in-situ analysis of δ 94/90 Zr in zircon using SIMS. The SIMS provided the high spatial resolution (8 μm) and good analytical precision (0.1‰, 2SE) but long analytical time (4 min/ analysis). Zhang et al. (2019) determined the Zr isotopes in zircons using laser ablation MC-ICP-MS (LA-MC-ICP-MS). The LA-MC-ICP-MS analysis required the larger ablation spot size (>16 μm) to obtain the good within-run precision (0.08‰, 2SE) and intermediate precision (0.11‰, 2SD) for δ 94/90 Zr, but had the higher analytical efficiency (~1 min) compared with SIMS.
In this study, we optimized the fs-LA-MC-ICP-MS strategy for the determination of Zr isotopes in zircons. Femtosecond laser system was applied in the Zr isotope analysis for suppressing the potential isotopic mass fractionation during the reaction process between laser with zircons. A two-volume ablation cell was used to eliminate the position effect in the Zr isotope analysis. The Sr isotope standard solution of NIST SRM 987 was introduced into ICP to produce a "wet" plasma condition and correct the Zr isotope ratio shift. The new analysis scheme and data correction mode provided a more robust Zr isotope analysis, which was conducive to the large number of samples and the long-term analysis. In addition, we investigated the homogeneity of Zr isotope composition on four international zircon reference materials (GJ-1, 91500, Plešovice and Rak-17) in detail. Moreover, three natural zircons (Aus, Mala and Paki) were commended as the potential Zr isotope materials.

EXPERIMENT 1.Instrumentation
In-situ Zr isotope analyses were performed on a NEP-TUNE Plus MC-ICP-MS instrument (Thermo Fisher Scientific, Bremen, Germany), which was connected to a femtosecond laser ablation system at the State Key Laboratory of Geological Processes and Mineral Resources (GPMR), China University of Geosciences, Wuhan. The NEPTUNE Plus, a double focusing MC-ICP-MS, was equipped with seven fixed electron multiplier ICs, and nine Faraday cups fitted with 10 11 Ω resistors. The combination of the high sensitivity X skimmer cone and JET-sample cone was employed, and the mass spectrometer was operated in the low mass resolution mode. The Faraday collector configuration of the mass spectrometer was composed of an array from L4 to H4 to monitor 86 Sr + , 87 Sr + , 88 Sr + , 89 Y + , 90 Zr + , 91 Zr + , 94 Zr + (Table 1). An NWR Femto UC femtosecond system (New Wave Research, Fremont, CA, USA), which consisted of a 300 fs Yb: KGW femtosecond laser amplifier (PHAROS, Light Conversion Ltd., Vilnius, Lithuania) with a wavelength of 257 nm. The laser ablation was equipped with a two-volume cell, which has the constant distance between the laser ablation position with the aerosol extraction position. Helium gas was filled into the ablation cell, while argon was mixed into the sample-out line downstream from the ablation chamber prior to entering the torch. A signal-smoothing device was used downstream from the sample cell, which approximately eliminated the short-term variability of the signal (Hu et al., 2015). Details of the instrumental operating conditions and measurement parameters are summarized in Table 1.

Samples
Three international zircon reference materials (GJ-1, 91500, Plešovice) were used in this study. GJ-1 was used for tuning and as the bracketing standard for isotope measurements. These zircons are widely used as reference materials for the micro analysis of U-Pb, Hf and O isotopes. A new international zircon reference material Rak-17 was analyzed in this study. Rak-17 was published by the G-Chron programme, which was an international proficiency testing programme devoted to U-Pb dating and developed by the International Association of Geoanalysts (IAG). Rak-17 was the millimeter-sized grains and has been determined 206 Pb/ 238 U of 295.56 ± 0.21 Ma by the isotope dilution-TIMS.
In addition, three natural zircon megacrysts were analyzed in this study. Paki zircons are deep red megacrysts (119 g) from Pakistan with a single grain diameter of 10-20 mm. Mala zircons come from Malawi with dense structure, large size of 60 mm and whole weight of 608 g. Aus zircons are transparent megacrysts (128 g) from Australia with a single grain diameter of 30-40 mm. All three zircon crystals were purchased through jewelers, so there is no accurate sampling location. In order to assess their homogeneity of Zr isotope composition, approximately 60% of three zircon megacrysts have been crushed and sieved to 1-2 mm. Then individual fragments were randomly selected and mounted on epoxy mounts. More than 30 fragments in each zircon megacryst were analyzed by fs-LA-MC-ICP-MS randomly. A part of the fragments was milled to powder using agate grinding and was acid-digested for Zr isotope ratio measurements using solution-MC-ICP-MS.

Analytical Method 1.3.1 fs-LA-MC-ICP-MS analysis of Zr isotopes
The laser ablation system, MC-ICP-MS and spray chamber were combined as in Fig. 1. The solution introduction system (glass spray chamber and 100 μL·min -1 nebulizer) was used to introduce the Sr isotope solution (NIST SRM 987, 100 ng·g -1 ) for the isotope ratio shift correction. A high flow rate of argon gas (1.0 L·min -1 ) was applied for the stable Sr signal. It needs to be emphasized that the Sr isotope solution was diluted by Milli-Q water rather than 2% HNO 3 . In addition, a small amount of N 2 (10 mL·min -1 ) was added to the carrier gas flow behind the signal-smoothing device by a simple Y connector. A zircon reference material, GJ-1, was used to optimize the instrumental parameters, including the He and Ar gas flow rates, the torch position, the RF power setting, and the source lens settings. The routine data acquisition consisted of one block of 100 cycles (0.524 s integration time per cycle): the first 20 cycles for background collection (no laser ablation) and the remaining 80 cycles for signal collection. The Sr solution (NIST SRM 987, 100 ng·g -1 ) was aspirated during the entire analytical procedure. Typical laser operating parameters are listed in Table 1.
In this study, Sr element was used as doping element to correct the Zr isotopic fractionation behavior. Many different data reduction methods can achieve this aim. We selected the exponential law (Eq. 1) to describe the isotopic fractionation behaviors of target elements and doping elements.
where r and R denote the measured and isotopic-fractionationcorrected (or "true") ratios, respectively, m i and m j are the atomic masses of the isotopes of interest, and f is the mass frac-tionation factor. The mass fractionation factor of Sr (f Sr ) can be obtained by assuming 88 Sr/ 86 Sr of 8.375 209 (Zhang et al., 2018). Then f Sr was used to correct Zr isotope mass fractionation. Previous researches have proved that the mass fractionation factors of target elements and doping elements were not consistent (Thirlwall, 2002;Maréchal et al., 1999). The corrected Zr isotope ratios by doping element of Sr had a system deviation from the reference value. Therefore, standard-sample bracketing method (SSB) was used to solve the system deviation problem (Yang, 2009). The SSB method is a widely used calibration method in LA-MC-ICP-MS. A matrix matched standard material with a known Zr isotopic composition is measured before and/or after the sample. In this study, the analytical sequence proceeded by bracketing "3-4" samples with the measurement of "1" standard. The relative difference in the  isotope ratios between the sample and standard was calculated using δ notation, as shown in Eq. 2 (using 94 Zr/ 90 Zr as example), and was expressed in parts per thousand.
The zircon reference materials GJ-1 was employed as the matrix-matching standard. All δ values were reported relative to GJ-1. Two Zr isotope ratios were calculated in this study, including 94 Zr/ 90 Zr and 94 Zr/ 91 Zr. The within-run precision for one single spot analysis was calculated by the combination of the SE (standard error) from samples and standards using the following equation.
where SE sample is the standard error of the Zr isotope ratio of the sample, SE standard-1 and SE standard-2 are the standard error of the Zr isotope ratio of the former and the latter bracketing standard, respectively. All data reduction for the fs-LA-MC-ICP-MS analysis of Zr isotope ratios was conducted using "Iso-Compass" software . Data were reduced in the following order and based upon user-selected background and sample integration intervals. The average background was first subtracted from the raw Zr signals. A simple 89 Y 1 H + interference correction was implemented, using the measured 89 Y signal with a yield of 89 Y 1 H + of~0.000 01 obtained by measuring a xenotime (YPO 4 ). The concentration of Sr in zircons is low. Therefore, we did not consider the possible impact of Sr in the sample.

The TIMS analysis of Zr isotopes
All of zircons investigated in this study were analyzed by TIMS with double spike method relatively to the IPGP-Zr standard. These measurements were conducted at the State Key Laboratory of Geological Processes and Mineral Resources (GPMR), China University of Geosciences, Wuhan. Detailed descriptions of the TIMS procedure for Zr isotopic analyses are given in Feng et al. (2020).

RESULTS AND DISCUSSION 2.1 Effect of theAddition of N and "Wet" Plasma Condition
In this study, the doping element Sr was added by the traditional solution nebulization mode. A large quantity of water was inevitably introduced into ICP, producing a "wet" plasma condition. In addition, the yield of hydrides and oxides could increase when water was introduced in ICP and broken down to oxygen and hydrogen. Therefore, we investigated the effect of "wet" plasma on the signal intensity of Zr isotopes and the yield of 89 Y 1 H + through the analyses of GJ-1 and a xenotime. The results are shown in Fig. 2. The max signal intensity of 90 Zr increased from 3.5 to 5.9 V when the plasma condition changed from "dry" to "wet" (Fig. 2a). This result is contrary to previous reports. For example, Zheng et al. (2018) have reported a significant suppression of Fe isotopes signal if excess water was add-ed. Lin et al. (2019) found the Li intensity is reduced by a factor of~1.45 with the addition of the water. The Zr signal enhancement with the addition of water was contributed to the high first ionization energy of Zr (6.64 eV) and the improvement plasma temperature due to the introduction of hydrogen ion. In addition, if a small amount of N 2 (10 mL·min -1 ) was added into ICP under the "wet" plasma condition, the max signal intensity of 90 Zr can increase to 8.0 V. The mode of "wet" plasma combined with N 2 gas showed the good effect for the Zr isotope analysis in zircons.
The "wet" plasma could increase the yield of oxygen and hydrogen. Figure 2b shows that the yield of 89 Y 1 H + obtained from a xenotime increased from 1.0 × 10 -6 to 1.9 × 10 -6 under the max signal intensity condition, when the plasma condition changed from "dry" to "wet". If the mode of "wet" plasma combined with N 2 gas was used, the yield of 89 Y 1 H + would increase to 3.2 × 10 -6 . The level of 89 Y 1 H + yield does not affect the Zr isotope analysis. Therefore, we infer that the Zr isotope analysis in zircons would not be interfered by 89 Y 1 H + under the "wet" plasma condition.

Mass Fractionation Behaviors between Zr Isotopes and Sr Isotopes
Internal standardization by doping elements also has been applied to some isotope analyses, such as doping Ni to monitor the isotopic fractionation of Fe or Cu, doping Fe to monitor V, and doping Tl to monitor Pb (Schuth et al., 2017;Zhang et al., 2016;Lazarov and Horn, 2015;Oeser et al., 2014). The good correlation between the target isotopes and doping isotopes is essential for the successful internal standardization. Because target isotopes would be subject to two different mass fractionation processes, including laser ablation process and ICP-MS process. However, the doping isotopes only can monitor the mass fractionation processes in ICP-MS. We summarize the correlation between Zr isotope and Sr isotope in different analysis sessions over two months (Fig. 3). As shown in Fig. 3, the good correlation between 94 Zr/ 90 Zr and 94 Zr/ 91 Zr with 88 Sr/ 86 Sr (R 2 = 0.996 5 and 0.996 8, respectively) were observed under the long-term analysis time. The results indicate that there was a good coupling relationship between the 88 Sr/ 86 Sr added to the ICP through solution nebulization method and the Zr isotope ratios from the ablated aerosols. By monitoring the change of 88 Sr/ 86 Sr during the fs-LA-MC-ICP-MS analysis, the isotope ratio drift behavior of Zr isotopes can be accurately corrected. The analytical results of GJ-1 during a single session were shown in the Fig. S1 (supplementary information). The raw data of 94 Zr/ 90 Zr over 1 hour showed a bad reproducibility (2RSD = 0.41‰). After the correction by using Sr isotope as the internal standard, the reproducibility of 94 Zr/ 90 Zr can be improved to 0.08‰ (2RSD). The improvement of reproducibility confirmed the advantage of the doping element (Sr) correction method.
The common SSB method can also correct the isotope ratio drift by repeatedly measuring the reference materials. To properly apply the SSB technique, it is essential that the mass fractionation is the same for the sample and the standard, and the isotope ratio drift between bracketing standards should be small during the measurement session (Lu et al., 2020;Yang, 2009). However, in practice, it is difficult to maintain the ICP stability for a long time, especially for LA-ICP-MS analysis.
Accurate Determination of Zr Isotopic Ratio in Zircons by Femtosecond Laser Ablation MC-ICP-MS Therefore, when using the SSB method to carry out the highprecision Zr isotope analysis, it is often necessary to constantly adjust the instrumental parameters, which is time consuming. Internal standardization by doping elements can make up for the shortcomings of SSB method. Any characteristics that cause a change in isotope ratio in ICP can be recorded by the doping element, such as the change of plasma temperature, gas flow rate, mass loading or interface area of sample cone and skimmer cone. Therefore, using Sr solution to correct the Zr isotope ratios shift is conducive to long-term Zr isotope analysis in zircons, especially for detrital zircon samples, which are required a large number of analyses to adequately characterize the isotopic composition distribution.

Position Effect in Laser Ablation Cell
The position effect can lead to differences in the measured isotope composition for the same sample and has been reported by Xie et al. (2018). Previous studies have proved that the isotopic composition in the ablated sample aerosols is particle size dependent. In addition, the transport efficiency of aerosols with different sizes is variable at different sampling positions, which is caused by the change of the distance between the laser ablation position and the exit of the ablation cell. In this study, we used the two-volume cell (TV-2), which has the constant distance between the laser ablation position with the aerosol extraction position. The TV-2 cell is a 100 mm × 100 mm rectangle, which can place a large number of samples (Fig. 4). The zircon reference material GJ-1 was placed in the center of TV-2. Four natural zircons Paki were placed in the four corners of the TV-2. The analysis results of four Paki zircons show that the measured values of δ 94/90 Zr GJ-1 were consistent within uncertainty (Fig. 4). Other three zircons (91500, Plešovice, Mala) were also analyzed with the same method and the data were listed in the supplementary information (Table S1). All values of δ 94/90 Zr GJ-1 in the three zircons at different positions of the cell showed the consistent values within uncertainty, respectively. The results indicate that the TV-2 cell did not induce discernible spatial fractionation for the Zr isotope analysis in zircons.

Effect of Laser Ablation Parameters
The isotope fractionation induced by the laser ablation process is controlled by sample properties and laser ablation parameters, such as color, surface reflectivity, thermal property, crystal structure, and defect for samples and laser wavelength, energy density, pulse width, and repetition rate of the laser system. In this study, analytical standards and samples were zircons, indicating the weak matrix difference, especially the chemical composition. Therefore, we pay more attention to the assessment of the effect of laser ablation parameters, such as the laser spot size and energy density. Figure 5 shows the results of δ 94/90 Zr GJ-1 in GJ-1 as a function of laser spot size and laser energy. In Fig. 5a, GJ-1 was repeatedly measured at different spot sizes from 10 to 30 μm with the fixed laser energy of 50% and ablation frequency of 2 Hz. The analysis of GJ-1 with the spot size of 20 μm was used as the normalized standard. The results show that when the laser spot size ranged from 15 to 30 μm, the signal intensity of 90 Zr varied from 0.5 to 2 times, but δ 94/90 Zr GJ-1 in GJ-1 remained consistent within uncertainty. However, δ 94/90 Zr GJ-1 in GJ-1 at the spot size of 10 μm were sig-   nificantly low (-0.34‰--0.43‰). In Fig. 5b, GJ-1 also were repeatedly measured at different laser energy from 20% (1.7 J·cm -2 ) to 70% (3.1 J·cm -2 ) with the fixed spot size of 20 μm and ablation frequency of 2 Hz. The analysis of GJ-1 with the laser energy of 70% was used as the normalized standard. The results show that when the laser energy changed from 20% to 70%, the signal intensity of 90 Zr varied from 0.18 to 1.3 times, but δ 94/90 Zr GJ-1 in GJ-1 remained consistent within uncertainty.
Combining the results of the two experiments, we demonstrate that changing laser ablation spot size, laser energy and the signal intensity of Zr (0.18 to 2 times) within a large range did not cause the distinguishable Zr isotope fractionation. One exception was the GJ-1 analysis at the spot size of 10 μm, in which a significantly low value of δ 94/90 Zr GJ-1 was observed (-0.34‰--0.43‰). The possible reason was that the energy distribution of the femtosecond laser was Gaussian distribution, and was easy to produce a deep ablation hole under a small ablation spot. This may lead to potential down-hole fractionation of Zr isotopes. Therefore, we recommend that the laser spot size should be larger than 10 μm when using a femtosecond laser to analyze the Zr isotope ratios in zircons.

Results of Reference Materials and Natural Zircons Using fs-LA-MC-ICP-MS
The precision of a single spot analysis using LA-MC-ICP-MS is called as repeatability or within-run precision, which is defined as the double standard error of the mean (2SE) and added the uncertainties from the normalized standards (Eq. 3). In this study, the within-run precision was mainly influenced by the Poisson noise (counting statistics) originating from the discrete nature of the particles in the ion beam. It means that the within-run precision would be strongly dependent on the signal intensity. Figure 6 shows the relationship between the 90 Zr signal intensity and the within-run precisions of δ 94/90 Zr GJ-1 (2SE) obtained from fs-LA-MC-ICP-MS in this study and ns-LA-MC-ICP-MS (Zhang et al., 2019). When the 90 Zr signal intensity was >9 V, both of fs-LA-MC-ICP-MS and ns-LA-MC-ICP-MS provided the same within-run precision of 0.05‰-0.13‰ (average values = 0.08‰). However, at the low signal intensity of <9 V, fs-LA-MC-ICP-MS showed the better within-run precision compared to ns-laser at similar signal intensities.
The intermediate precision, which is defined as twice the standard deviation (2SD), was obtained by the repeated analyses of GJ-1 over one year, and using GJ-1 itself as the external standard. This method can completely avoid the potential matrix effect. The results show that the average value of δ 94/90 Zr GJ-1 was 0.00‰ ± 0.08‰ (2SD, n = 539). The intermediate precision of 0.08‰ not only reflected the long-term repeatable range of our Zr isotope analysis method, but also showed the good homogeneity of Zr isotope composition in GJ-1. Other six zircon samples (91500, Plešovice, Rak-17, Paki, Aus and Mala) were also analyzed using the presented in-situ Zr isotopic method during one year period. GJ-1 was used as the external standard. The common zircon references of 91500 and Plešovice showed the significantly stable values of δ 94/90 Zr GJ-1 (-0.03‰ ± 0.11‰ (2SD, n = 242) and 0.13‰ ± 0.11‰ (2SD, n = 251), respectively) (Fig. 7), suggesting that they have the homogenous composition of Zr isotope and can be used as the reference materials for the Zr isotopic microanalysis. The results of three natural zircons obtained by repeatedly analyzing the different grains appeared to lack significant intergrain variations of Zr isotopic ratios (Fig. 7), indicating that they can be considered as the potential candidates for the Zr isotope analysis in zircons. The zircon Rak-17, which was used as the U-Pb dating reference by IAG, has a large intermediate precision of 0.22‰ (Table 2). Therefore, we did not suggest Rak-17 as an eligible reference for the Zr isotopic microanalysis, even though it has the homogenous U-Pb dating age.
The TIMS values of seven zircons (GJ-1, 91500, Plešovice, Rak-17, Paki, Aus and Mala) have been listed in Table 2. In this study, Zr standard reference NIST SRM 3169 was used as the external standard. Forlaboratory comparisons, the measured values of δ 94/90 Zr SRM3169 were converted to IPGP-Zr and GJ-1. The corresponding uncertainty (2SD) was calculated according to the error propagation formula. As a result, the δ 94/90 Zr GJ-1 values obtained by fs-LA-MC-ICP-MS are consistent with the results from TIMS for all of zircons, confirming the  accuracy of our method.

CONCLUSIONS
In this study, the Sr solution was spiked to ICP-MS by the nebulization method to monitor and correct the isotopic ratio drift in ICP. This new Zr isotope correction method is more suitable for long-term Zr isotope analysis in zircons, especially for detrital zircon samples, because of the less effects from the fluctuation of plasma. The adding N 2 and the "wet" ICP techniques were simultaneously applied to provide more robust plasma conditions. A 1.7-fold increase in Zr signal was observed under the "wet" plasma condition, while the yield of 89 Y 1 H + also remained at a very low level (3.2 × 10 -6 ) and did not affect the Zr isotope analysis. We investigated the position effect of the laser ablation cell and the influence of the femtosecond laser ablation parameters. The results show that there was no obvious Zr isotope position effect in the two-volume cell (TV-2). Changing laser ablation spot size, laser energy and the signal intensity of Zr (0.18 to 2 times) within a large range did not cause the distinguishable Zr isotope fractionation. One exception was the GJ-1 analysis at the spot size of 10 μm, in which a significant low value of δ 94/90 Zr GJ-1 was observed (-0.34‰--0.43‰). Therefore, the femtosecond lasers Zr isotope analysis should be carried out cautiously under the high spatial resolution mode (< 10 μm). fs-LA-MC-ICP-MS provided the better within-run precision compared to ns-laser at similar signal intensities. Under these optimized experimental conditions, a long-term external reproducibility of~0.1‰ (2SD) for δ 94/90 Zr is routinely obtained for zircons. Long-term repeated analyses of seven zircons (GJ-1, 91500, Plešovice, Rak-17, Paki, Aus and Mala) over one year showed they had the homogenous Zr isotope composition and could be considered as the potential candidates for the Zr isotope analysis in zircons, except Rak-17 which presented the large Zr isotope variation. The analytical results using fs-LA-MC-ICP-MS are in good agreement with the TIMS values, confirming the reliability of the proposed method.