1 Introduction

Spontaneous fission of naturally occurring ²³⁸U produces linear tracks of radiation damange in the crystal lattice of apatite, zircon and sphene at a rate of ~215 tracks per ng of U per million years. In a series of seminal papers published during the 1960s, Robert Fleischer, Paul B. Price and Robert Walker showed how these fission tracks can be revealed to the human eye (aided by an optical microscope) by acid etching, forming the basis of a unique geochronometer (Price and Walker, 1962, 1963; Fleischer et al., 1964, 1965, 1975). In its most basic form, the fundamental age equation of the fission track method can be written as follows:

(1)

where λ is the total decay constant of ²³⁸U (1.55125×10^-10a^-1; Jaffey et al., 1971), λ_f is the fission decay constant (with estimates ranging from 7.9 to 8.7×10^-17a^-1 and an IUPAC-IUGS-recommended value of 8.45 ± 0.10 ×10^-17a^-1; Holden and Hoffman, 2000), ρ_s is the density (tracks per unit area) of the spontaneous fission tracks on an internal crystal surface, [²³⁸U] is the current number of ²³⁸U atoms per unit volume, and L is etchable length (or, more precisely, the mean equivalent isotropic length; Galbraith and Laslett, 1988) of the fission tracks. During most of the 20^th century, the only way to reliably estimate [²³⁸U] was by proxy, via induced tracks produced by neutron irradiation. During the 1970s, the fission track community experimented with a number of alternative designs to turn induced track densities into ²³⁸U-concentrations (such as the subtraction, population, and re-etch methods; Gleadow, 1981; Hurford and Green, 1982; Galbraith, 1984; Wagner and Van den Haute, 1992), before ultimately settling on the external detector method (EDM), in which the induced tracks are recorded in a mica- (or plastic) detector attached to the polished grain mount during irradiation. This detector is subsequently etched and counted separately (Fleischer and Hart, 1972; Hurford and Green, 1983; Hurford, 1990). The EDM emerged as the dominant method because it offers significant advantages over the other approaches. First, the EDM was recognised as the first geochronological method capable of routinely producing single grain age estimates, leading to the development of detrital geochronology. Second, it is least affected by the presence of uranium zoning. Third, it allows the fission decay constant (which used to be poorly constrained) to be hidden in a calibration constant (ζ) by reformulating the age equation as follows:

(2)

where ρ_i is the surface density of the induced fission tracks in the mica detector and ρ_d is the surface density of the induced fission tracks in a dosimeter glass of known (and uniform) U concentration. The latter value is needed to ‘recycle’ the calibration constant from one irradiation batch to the next as neutron fluences might vary between irradiations, or within a sample stack if the nuclear reactor is not well thermalised. ρ_s, ρ_i and ρ_d are unknown but can be estimated¹ by counting the number of tracks N_* over a given area A_* (where * is either s for ‘spontaneous’, i for ‘induced’ or d for ‘dosimeter’):

(3)

It is customary for the spontaneous and induced fission tracks to be counted over the same area (i.e., A_s = A_i). This can be done using an automated microscope stage (Smith and Leigh-Jones, 1985; Dumitru, 1993) or by simply repositioning the mica detector on the grain mount after etching (Jonckheere et al., 2003). Using these measurements, the apparent fission track age ( ) is given by

(4)

where is obtained by applying Equation 4 to an age standard and rearranging. Equation 4 represents a classic case of a ‘matched pairs’ experimental design (Galbraith, 2010). By counting the spontaneous and induced tracks over exactly the same area, the age calculation reduces to a simple comparison of two Poisson-distributed variables (N_s and N_i). This enables an explicit maximum likelihood formulation, which greatly simplifies all subsequent statistical analyses (Galbraith, 2005). As a result, it is fair to say that the fission track method represents the gold standard among geochronometers in terms of statistical rigour. Unfortunately, the EDM also has a number of practical shortcomings:

1. It requires irradiation with thermal neutrons. This greatly increases sample turnaround times and poses administrative and safety headaches. These problems are only getting worse with time as research reactors are becoming increasingly rare.
2. Etching the external (mica) detectors requires handling hazardous hydrofluoric acid.
3. The method is tedious and time consuming, as it requires counting three different track densities (spontaneous, induced and dosimeter).

The late 1980s and 1990s saw the emergence and maturation of ion (SIMS) and laser (LA-ICP-MS) microprobe technology allowing geoscientists, for the first time, to determine ppm-level ²³⁸U-concentrations of solid samples with microscopic resolution directly, rather than by proxy via induced fission of ²³⁵U. By the early 2000s, LA-ICP-MS technology had matured enough to be used for fission track dating. A proof-of-concept study by Hasebe et al. (2004) laid the foundations for this new approach, which has gradually gained popularity since (Hasebe et al., 2009; Chew and Donelick, 2012; Soares et al., 2014; Gleadow et al., 2015; Abdullin et al., 2016). But despite this popularity, relatively little work has been done to develop the statistical analysis of this new type of fission track data (Galbraith, 2010). Important issues such as error propagation and the treatment of ‘zero track’ grains are largely unsolved for LA-ICP-MS-based fission track (LAFT) dating. This limits the applicability of the method to young samples, and potentially compromises data quality in compositionally zoned crystals. The present paper addresses these issues by defining four different strategies towards LAFT dating, using either an absolute dating (Section 2) or a zeta calibration (Section 3) approach and using either a single or multiple laser spots per grain to account for uranium heterogeneity (Section 4). Section 5 introduces a method to calculate meaningful ages and uncertainties for very young and uranium-poor samples lacking any spontaneous fission tracks. With the age equation in place, we can redefine the pooled age and visualise LAFT data on radial plots, allowing us to assess whether the single grain ages are consistent within the analytical uncertainties (Section 6). Using a logarithmic transformation, this qualitative assessment can be formalised with a Chi-square test for age homogeneity. Samples failing this test are said to be ‘overdispersed’ and may be modelled using continuous or finite mixtures (Section 6). The methods described in this paper were implemented in a software package called IsoplotR. Section 7 applies IsoplotR to two real fission track datasets, providing practical examples of the absolute and zeta calibration approach using single and multiple laser spots per grain, and including zero track grains and finite mixtures.

2 ‘absolute’ fission track dating

The development of LA-ICP-MS has allowed fission track geochronologists to routinely determine the ²³⁸U-concentration of sand-sized crystals of apatite, zircon and sphene directly, rather than by proxy through neutron-induced fission of ²³⁵U. This opens up the possibility to use Equation 1 to calculate fission track ages:

(5)

where q is an ‘efficiency factor’ (0 < q < 1, Iwano and Danhara, 1998; Enkelmann and Jonckheere, 2003; Jonckheere, 2003; Soares et al., 2013) that depends on the mineralogy, etching conditions and observer; and [²³⁸Û] is the measured ²³⁸U-concentration, expressed in atoms of ²³⁸U per unit volume. Please note that uranium sitting near the counting surface of the grain will contribute more to the spontaneous track budget than more deeply seated uranium. So in the presence of uranium zoning, the LA-ICP-MS measurements should be weighted inversely with depth. For the remainder of this paper, we will assume that this depth weighting has either been done, or that the grains are not significantly zoned perpendicular to the etched surface. The effect of lateral zoning is covered in Section 4. Error propagation of Equation 5 follows the usual first order Taylor expansion (neglecting the systematic uncertainties in A_s, L and q):

(6)

where s[a] stands for “the standard error of a”.

3 The ‘zeta calibration’ approach

The apparent age estimates produced by Equation 5 are at most as accurate as the value for [²³⁸Û]. Unfortunately, it is notoriously difficult to measure isotopic concentrations accurately by LA-ICP-MS. Such concentration measurements are typically done by monitoring the isotope of interest (i.e., ²³⁸U) relative to a stoichiometric isotope (^xX, where x = 43 and X = Ca for apatite, or x = 29 and X = Si for zircon) and comparing the measured signal ratios [²³⁸Û∕^x] with those of reference materials of known composition. Unfortunately, no such materials exist for apatite, zircon and sphene (although some natural samples come close to fulfilling the requirements; Soares et al., 2014), leading most laboratories to use NIST standard glass (SRM610 and SRM612; Pearce et al., 1997) instead (Longerich et al., 1996; Norman et al., 1998; Liu et al., 2008). However, glass interacts very differently with UV laser light than most minerals, and therefore the accuracy of LA-ICP-MS based U-concentration measurements often leaves much to be desired. To avoid these problems and remove the need to measure absolute concentrations, all the poorly constrained factors that relate the raw mass spectrometer measurements (i.e., the ²³⁸U/^xX-signal ratios) to the ²³⁸U-concentrations can be grouped into a ‘zeta’ calibration constant similar to that used for the EDM (Equation 4):

(7)

where _icp is determined by analysing a reference material of known age, precisely as is done for the EDM. But unlike the EDM, the ζ-calibration constant for LAFT dating must be determined anew for every analytical session, to account for any changes in ablation and plasma conditions which may occur over time. Therefore, _icp is known as a ‘session zeta’ (Hasebe et al., 2009). Single grain age uncertainties for the zeta-calibration approach are given by

(8)

4 Uranium zoning

The U-uncertainties (s[²³⁸Û] in Equation 6 and s[²³⁸Û∕^x] in Equation 8) account for both the analytical error of the LA-ICP-MS measurements and the uranium heterogeneity of the sample. The EDM is less sensitive to uranium zoning, thanks to the ‘matched pairs’ experimental design discussed in Section 1. As long as the spontaneous and induced track densities are counted over the same area, the presence of U-rich or U-poor zones has no effect on the resulting age. Unfortunately, things are not so straightforward for LAFT data. Suppose, for example, that the analyst has placed a round laser spot in the top half of the strongly zoned grain shown in Figure 1. This would result in a high uranium concentration (or isotopic ratio measurement) and a small analytical uncertainty, but would be completely unrepresentative of the average composition of the grain. Blindly combining such a single spot measurement with the number of spontaneous fission tracks counted over the entire crystal would produce a precise but grossly inaccurate age. For old and uranium-rich samples, it is often possible to detect and avoid the problems caused by uranium zoning by carefully observing the spatial distribution of the spontaneous tracks. But this is generally impossible for young and/or uranium-poor samples. Two analytical strategies may be used to address this issue and capture the naturally occurring uranium dispersion by LA-ICP-MS.

Figure 1: Example of a compositionally zoned apatite crystal using the External Detector Method (EDM). In this old sample, the effect of an order of magnitude difference in uranium concentration is visible in both the spontaneous (left) and induced (right) track distribution. For young and/or uranium poor samples, the uranium zoning would only be visible in the induced track distribution. For LAFT data, it would not be visible at all. Modified after Ksienzyk et al. (2014).

A first approach is to only count the spontaneous fission tracks contained within the area occupied by the laser ablation spot (Figure 2a). This effectively imposes a ‘matched pairs’ experimental design on the LAFT data which, by definition, should capture uranium zoning effects. One can then safely plug the analytical errors of the LA-ICP-MS measurements into Equation 6 or 8 to propagate the age uncertainty. However, the single spot approach only uses a small part of the grains (so as to fit within the smallest one of them), containing a small number of spontaneous tracks. Because the latter dominate the uncertainty budget of LAFT data, the resulting single grain age estimates are imprecise despite the high precision of the isotope ratio measurements.

The only way to avoid these problems and enlarge the spontaneous track counting area is to somehow quantify the naturally occurring uranium heterogeneity outside the ablation spot. This can be achieved by ablating some or all of the grains multiple times (Figure 2b) and making some form of parametric assumption for the uranium heterogeneity. A reasonable assumption would be that the ²³⁸U-concentrations or ²³⁸U/^xX-ratios within each grain follow a lognormal distribution with location parameter u_j and scale parameter (for 1 ≤ j ≤ n, where n is the total number of grains):

(9)

where (u_j,v_j) stands for “a Normal distribution with mean u_j and variance v_j”. These parameters may be estimated by analysing multiple laser spots per grain. Let û_jl be the logarithm of the l^th (out of m_j) ²³⁸U-concentration or ²³⁸U/^xX-ratio measurement of the j^th grain, and let s[û_jl] be its standard error (which approximately equals the relative measurement uncertainty). If all the LA-ICP-MS measurements employ the same spot size, then u_j can be estimated by taking the arithmetic mean of the single spot measurements û_jl (for l = 1...m_j):

(10)

and v_j is estimated by:

(11)

This formula clearly does not work for small grains that only fit a single laser ablation spot (m_j = 1). It is therefore necessary to make the simplifying assumption that v_j = v for all j. In that case, v may be estimated from the pooled measurements as follows:

(12)

where n degrees of freedom are subtracted to account for the n estimated parameters (û_j, for j = 1...n). The single grain ages can be estimated by substituting exp[û_j] for [²³⁸U] or [²³⁸U∕^xX] in Equation 5 or 7. In principle, we could obtain the age uncertainties by substituting exp[û_j] for s[²³⁸Û] or s[²³⁸Û∕^x] in Equation 6 or 8. However, doing so would assign exactly the same relative U-concentration or -ratio uncertainty to all grains, irrespective of whether the ablated area (m_jA_icp, where A_icp corresponds to the spot size) was small or large compared to the area over which the spontaneous fission tracks were counted (A_sj). This problem can be solved by using the following hybrid expression:

(13)

which reduces to the LA-ICP-MS measurement uncertainties if the laser spot size equals the fission track counting area (i.e., _j = ∑ _l=1^m_j(s[û_jl])² if A_icp = A_sj), and to Equation 12 if the laser spot size is much smaller than the fission track counting area ( _j = if A_icp  ≪  A_sj).

Figure 2: Two strategies to deal with the uranium zoning problem shown in Figure 1. a – enforce a matched pairs experimental design by only counting the spontaneous fission tracks located within a single ablation spot (A_sj, solid lines). b – capture the magnitude of the uranium heterogeneity by placing several laser spots within a larger counting area (i.e., A_icp, dashed lines).

5 The ‘zero tracks’ problem

The uncertainty propagation (Equations 6 and 8) shows that the standard error of the fission track age varies in inverse proportion to the spontaneous fission track count. Because the number of spontaneous fission tracks is usually quite low (especially for the single spot approach), the analytical uncertainty of fission track ages is generally dominated by the 1∕N_s term. For example, decreasing N_s from 100 to 4 inflates the relative age error associated with the spontaneous tracks from 10% to 50%. This further increases to 100% if N_s = 1. For a sample lacking any spontaneous fission tracks (N_s = 0), Equations 2-8 yield a zero age with infinite uncertainty. This effectively precludes any subsequent statistical analysis, such as averaging, mixture modelling, or thermal history analysis. For EDM data, the zero track problem can be solved by adding half a count to both the spontaneous and the induced tracks. This significantly improves the precision of the ages whilst only having a minor effect on their accuracy (Galbraith, 2005, p.80). Unfortunately this approach cannot be directly applied to the LAFT method, due to the absence of induced fission tracks. A pragmatic solution to this problem is to recast the LAFT age Equations (5 and 7) into an ‘EDM-like’ form:

(14)

where _icp = 1 for the absolute dating method, and r_j is a constant of proportionality that converts the (log)normally distributed ²³⁸U-concentration or ²³⁸U/^xX-ratio measurements of the j^th grain (û_j) into a Poisson distributed ‘equivalent induced track count’ (_ij) over an area equal to that counted for the fossil tracks:

(15)

r_j can be estimated using the property that the mean of a Poisson variable equals its variance:

(16)

Combining Equations 15 and 16, it is easy to show that

(17)

which can be plugged into Equation 14 to estimate the fission track age. Its standard error is given by

(18)

where the first term can be removed if the absolute dating method is used.

6 Radial plots and statistical models

The fission track method is a low precision technique. As explained in the previous section of this paper, single grain age uncertainties of 100% or more are not uncommon, especially in young and uranium-poor samples. To overcome this limitation and ‘beat down the noise’, it is generally necessary to analyse multiple grains from a sample and average them. There are several ways to do this. The simplest approach is to assume that all the grains have exactly the same age. This age can be estimated by pooling all the spontaneous track counts (N_s = ∑ _j=1ⁿN_sj) and uranium measurements (û) together, with

(19)

The assumption of a single underlying age component can be tested by calculating the following statistic (Galbraith, 2010):

(20)

where

(21)

and

(22)

in which s[ _j] indicates the standard error of the j^th age measurement propagating the internal sources of analytical uncertainty only. This means that the uncertainties associated with _icp or with the fractionation correction of the LA-ICP-MS measurements should be omitted. Fortunately, those uncertainties are usually much smaller than 1∕N_s. The probability of observing a value χ_stat² under a Chi-square distribution with n - 1 degrees of freedom is called the ‘p-value’. If the latter falls below a pre-defined cutoff value of 0.05, say, then we have justifiable grounds to conclude that the data are ‘overdispersed’, rendering the pooled age inappropriate. The dispersion of fission track data can also be assessed visually using a radial plot (Galbraith, 1988, 1990). This is an (x_j,y_j)-scatterplot that sets the standarised estimates of the transformed variables z_j (Equation 21) out against their respective analytical precision:

(23)

where z_∘ is some conveniently chosen central value. This scatterplot may be decorated with a 2-sigma error bar about the origin and a radial scale drawn at some appropriate distance away from it. The power of the radial plot derives from its unique ability to simultaneously display the value and analytical uncertainty of measured quantities. The former is given by the slope of a line connecting the origin with any given sample point. The latter is given by its horizontal distance from the origin. In order for a fission track age population to be compatible with a single age component, the sample points should approximately plot within a band defined by two parallel lines going through the upper and lower limits of the 2-sigma error bar. Samples that fulfil this requirement generally pass the Chi-square test and can safely be represented by their pooled age. Samples that do not fit within a 2-sigma band and subsequently fail the Chi-square test may do so for several reasons. One possiblity is that their component grains exibit a range of annealing temperatures due to uranium heterogeneity, which may result in a continuous distribution of single grain fission track ages when subjected to a protracted history of slow cooling (Gleadow et al., 1986). In this case, the pooled age model is clearly inappropriate and is better replaced by a ‘random effects model’ in which the apparent fission track ages are assumed to follow a (log)normal distibution with two sources of variance:

(24)

whose unknown parameters (μ and σ) can be estimated from the data using the method of maximum likelihood (Galbraith, 2005). Let and be those estimates, then exp[] is known as the ‘central age’ and as the ‘(over)dispersion’. For the EDM, it can be shown that, when = 0, the central age equals the pooled age. This is not exactly the case for LAFT data, although the difference is usually small (Galbraith, 2010).

A second scenario in which fission track data might fail the Chi-square test is when the true age distribution consists of a finite mixture of discrete age components. If exp[μ_k] is the true age of the k^th (out of K) such components, and π_k is the proportion of grains belonging to it (with π_K = 1 -∑ _k=1^K-1π_k), then these 2K-1 unknown parameters can be estimated by numerically maximising the following log-likelihood function:

(25)

in which f(μ_k|z_j,s_j) equals the probability of observing a value μ_k under a Normal distribution with mean z_j and standard deviation s_j. Numerical recipes for solving Equations 24 and 25 are provided by Galbraith (2005).

a. b.

7 Examples

The methods introduced in this paper have been implemented in a software package called IsoplotR that exists within the statistical programming environment R. IsoplotR can be called from the command line by installing R from http://r-project.org and typing

at the command prompt. The package can be then be loaded using

Alternatively, most functions can also be accessed via a user-friendly graphical user interface (GUI) that can be accessed from http://isoplotr.london-geochron.com. This Section will illustrate some of IsoplotR’s functionality with two fission track datasets. The first of these was kindly provided by Dr. Noriko Hasebe (Kanazawa University) and applies the absolute dating method with multiple laser spots per grain to a sample of Fish Canyon apatite (28.8 ± 0.8 Ma, Gleadow et al., 2015). The second dataset was provided by Dr. David Chew (Trinity College Dublin) using the zeta-calibration method with a single laser spot per grain. It combines an artificial mixture of 30 measurements on Durango apatite (31.02 ± 0.22 Ma, McDowell et al., 2005) with a 29 further measurements of Mount Dromedary apatite (98.5 ± 0.5, McDougall and Wellman, 2011). The raw input data can either be copied and pasted into the GUI, or read from a .csv input file:

where the format argument to the read.data function records the input format of the fission track data. format=1 indicates the EDM, format=2 refers to LAFT dating using the zeta calibration method, and format=3 marks LAFT data using the absolute dating method. The files Hasebe.csv, DUR.csv, and Chew.csv are provided in the Supplementary Information. To calculate the fission track ages of the first dataset:

Note that the Fish Canyon dataset contains a zero track grain (FCT6 in Table 1) that has automatically been dealt with using the methods of Section 5. To plot the first dataset on a radial plot (Figure 3a):

which also calculates the central age of the dataset (32.6 ± 0.11 Ma with 7% of dispersion on a 1σ level). To set the ζ calibration factor for the second dataset:

where tst is a two-element vector containing the standard age and its analytical uncertainty. This yields a value (DUR$zeta) of 342 ± 23 a cm² (1σ). To fit a two component mixture to the data and visualise the results as a radial plot:

which produces two peaks at 31.3 ± 2.3 Ma and 104.4 ± 8.1 Ma, in good agreement with the known ages of Durango² and Mt. Dromedary apatite (Figure 3.b).

8 Discussion and conclusions

This paper has introduced four different analytical approaches towards LAFT dating (Table 3).

1. The simplest approach is to apply the absolute dating method (Equation 5) using the analytical uncertainty of a single LA-ICP-MS measurement per grain (Equation 6). Though straightforward in theory, this approach may be inaccurate in practice due to the lack of matrix-matched uranium concentration standards for LA-ICP-MS (Section 2). The single spot approach is also imprecise for most common samples due to the small number of spontaneous tracks counted within the area occupied by the ablation pit.
2. The accuracy can be improved by tying the sample measurements to a reference material of known age, analysed under the same conditions as the sample, using the ζ-calibration constant and Equation 7. However, the precision may still be low if a single laser spot is used, for reasons given under analytical approach 1.
3. The number of spontaneous track counts per grain may be increased by placing multiple ablation pits on some or all of the grains. When applied to the absolute dating method, this yields precise but potentially inaccurate results, for reasons given under approach 1.
4. Potentially the most accurate and precise approach to LAFT dating combines the ζ-calibration approach with multiple laser spots. The only caveat of this method is that it requires considerably more work than the single spot absolute dating approach, because it involves multiple measurements of not only the samples, but the reference material as well.

Using these four approaches, we have seen how LAFT dating can overcome the problems of uranium zoning and zero track grains. Nevertheless, it is fair to say that the EDM is still superior to LAFT dating in some ways. As we saw in Section 1, EDM data are based on the ratios of two independent Poisson variables. This produces statistical models that are convincing on physical grounds and well verified empirically. In contrast with this, LAFT data are based on mixed ratios of Poisson and (log)normally distributed variables that are not easily captured in parametric statistical models. So instead such models use the ages and errors for further inference, rather than the raw fission track data. These models generally assume (log)normality of the analytical uncertainties, which is only approximately correct. When a LAFT dataset fails the chi-square test for age homogeneity, say, then it might not do so for geological reasons. Instead, it is possible that the dataset fails the test simply because its parametric assumptions are not met. With this caveat in mind, LAFT offers clear practical advantages over EDM-based fission track data including higher throughput and the ease of double dating with the U-Pb method (Chew and Donelick, 2012). It seems very likely, then, that LAFT will continue to gain popularity at the expense of EDM-based fission track dating. IsoplotR aims to facilitate this transition, and to help extract geologically meaningful information from future LAFT datasets with the greatest possible accuracy and easy of use.

9 Acknowledgments

The author would like to thank Noriko Hasebe and David Chew for sharing the datasets in Tables 1 and 2. Raymond Jonckheere and Christoph Glotzbach are gratefully acknowledged for their prompt and meticulous reviews. The author was financially supported by ERC Starting Grant 259504 (‘KArSD’) and NERC Standard Grants NE/I008837/1 and NE/K003232/1.

References