Revised error propagation of 40Ar/39Ar data,
including covariances

Pieter Vermeesch


Department of Earth Sciences, University College London, p.vermeesch [at] ucl.ac.uk

Abstract

The main advantage of the 40Ar/39Ar method over conventional K-Ar dating is that it does not depend on any absolute abundance or concentration measurements, but only uses the relative ratios between five isotopes of the same element –argon– which can be measured with great precision on a noble gas mass spectrometer. The relative abundances of the argon isotopes are subject to a constant sum constraint, which imposes a covariant structure on the data: the relative amount of any of the five isotopes can always be obtained from that of the other four. Thus, the 40Ar/39Ar method is a classic example of a ‘compositional data problem’. In addition to the constant sum constraint, covariances are introduced by a host of other processes, including data acquisition, blank correction, detector calibration, mass fractionation, decay correction, interference correction, atmospheric argon correction, interpolation of the irradiation parameter, and age calculation. The myriad of correlated errors arising during the data reduction are best handled by casting the 40Ar/39Ar data reduction protocol in a matrix form. The completely revised workflow presented in this paper is implemented in a new software platform, Ar-Ar_Redux, which takes raw mass spectrometer data as input and generates accurate 40Ar/39Ar ages and their (co-)variances as output. Ar-Ar_Redux accounts for all sources of analytical uncertainty, including those associated with decay constants and the air ratio. Knowing the covariance matrix of the ages removes the need to consider ‘internal’ and ‘external’ uncertainties separately when calculating (weighted) mean ages. Ar-Ar_Redux is built on the same principles as its sibling program in the U-Pb community (U-Pb_Redux), thus improving the intercomparability of the two methods with tangible benefits to the accuracy of the geologic time scale. The program can be downloaded free of charge from http://redux.london-geochron.com.

1 Introduction

Let z be a function f of two variables x and y:

z = f (x,y)
(1)

then standard error propagation of z by first order Taylor expansion yields:

     (   )      (   )
  2    ∂f-2  2    ∂f- 2 2   ∂f-∂f-
σ z =  ∂x   σx +  ∂y   σy + 2∂x ∂ycov(x,y)
(2)

where cov(x,y) is the ‘covariance of x and y’. Current practice in 40Ar/39Ar geochronology generally assumes that the third term of Equation 2 can be safely neglected. For example, consider the 40Ar/39Ar age equation:

T = -1-ln(1+ JR )
    λ40
(3)

with λ40 the decay constant of 40K, J the neutron irradiation parameter (see Section 11) and R the 40Ar/39ArK-ratio (where 40Ar is the radiogenic argon component and 39ArK is derived from neutron reactions on 39K). Then the age uncertainty is currently calculated as (Berger and York1970McDougall and Harrison1999Koppers2002):

 2   J2σ2R +-R2σ2J
σT =  λ240(1 + RJ)
(4)

which assumes that cov(R,J) = 0. This assumption cannot be correct because both R and J are calculated using the same mass fractionation corrections, detector calibrations, interference corrections and radioactive decay corrections. The analytical uncertainty associated with each of these factors results in correlated errors between R and J. Ignoring these error correlations affects both the precision and accuracy of the resulting 40Ar/39Ar ages.

The problem of correlated errors is not limited to R and J alone. It crops up literally everywhere in the 40Ar/39Ar method. In fact, a covariant structure is deeply engrained into the very DNA of the method, which is based on five isotopes (36-40) of a single element (Ar). This paper will show that, because the 40Ar/39Ar method is based on ratios rather than absolute abundances, it is subject to the peculiar mathematics of ‘compositional data’ (Section 2). Correlated errors are created during mass spectrometry, when the ion detector signals are extrapolated to ‘time zero’ and blank corrections are made (Sections 3 and 4). They occur as a result of mass fractionation corrections and detector inter-calibrations (Section 5). They arise when accounting for the effect of radioactive decay on 39Ar (from K), 36Ar (from Cl) and 37Ar (from Ca) (Section 7), or whenever an interference correction is made (Section 8). Error correlations occur when calculating J-factors (Section 11) and, as we have already seen at the beginning of this section, when applying the J-factor to solve the age equation (Section 12). Error correlations must also be taken into account when calculating the weighted mean of several 40Ar/39Ar age analyses (Section 13). Finally, the methods presented in this paper provide a simple and elegant way to account for the systematic biases that occur as a result of the uncertainty in the 40K decay constant and the atmospheric 40Ar/36Ar ratio (Section 12).

Thus, the existence of correlated errors affects every aspect of the 40Ar/39Ar method. The paper at hand presents an analytical solution to this problem as an alternative to the numerical approximations proposed elsewhere (Scaillet2000). A new computer code called Ar-Ar_Redux was developed with the aim to facilitate the adoption of the rigorous data reduction and error propagation methods presented herein (Section 14).

2 40Ar/39Ar as a compositional data problem

As mentioned in Section 1, the 40Ar/39Ar-age calculation is based on the 40Ar39ArK-ratio (R, see Equation 3), which can be calculated as follows:

    1 − a+ b+ c
R = ----------- − f
        d− e
(5)

with

a = [40  ]
 --Ar
 36Ara[36  ]
 --Ar
 40Arm (6)
b = [40Ar]
 36Ara[36Ar]
 37Arca[ 37Ar-]
  40Arm (7)
c = [40  ]
 36Ar
   Ara[36  ]
 38Ar
   Arcl[38   ]
 40Ar
   Arm (8)
d = [    ]
 39Ar
 40Arm (9)
e = [39Ar]
 37Arca[37Ar ]
 40Ar-m (10)
f = [40  ]
 --Ar
 39Ark (11)

in which ‘a’ stands for ‘air’, ‘ca’ for ‘Ca-salt’, ‘k’ for ‘K-glass’, and ‘cl’ for ‘Cl decay products’. The subscript ‘m’ stands for either ‘sample’ or ‘fluence monitor’. The meaning of this equation and the significance of the subscripts will be elaborated in later sections of this paper. The important point which needs to be made here is that Equations 6-11 only contain ratios, and do not depend on the absolute abundances of the different argon isotopes. In statistical terms, 40Ar/39Ar-measurements are said to be ‘compositional data’ and are subject to the peculiar mathematics of the compositional dataspace or ‘simplex’ (Aitchison1986). To illustrate the profound implications of this point, consider the simple situation of a K-bearing sample containing neither Ca nor Cl. In this case, terms b, c and e in Equation 5 disappear, which leaves us with a simple three component system comprised of 36Ar, 39Ar and 40Ar. Because we are only interested in the relative abundances of these three isotopes, they can be normalised to unity and plotted on a ternary diagram (Figure 1). It is well known that common summary statistics such as the arithmetic mean and standard deviation are unreliable in this data space. This is because the ternary diagram occupies a narrowly restricted subspace of the realm of real numbers. These restrictions cause problems because standard data reduction methods commonly assume that the data follow a Normal distribution, which requires support from -to +. The solution to this conundrum is to transform the data from the simplex to a Euclidean ‘logratio space’, in which standard Normal theory can be safely used (Aitchison1986Vermeesch2010).

In addition to opening compositional data to standard statistical analysis, the logratio transformation also simplifies the algebra of 40Ar/39Ar data reduction. This is because many of the calculations required for processing 40Ar/39Ar-data involve multiplication and exponentiation, which reduce to simple addition and multiplication after taking logs. The next sections of this paper will show how the raw mass spectrometric data can be cast into a logratio covariance structure for further processing, for both multi-collector (Section 3) and single collector (Section 4) instruments.


PIC
Figure 1: 40Ar/39Ar-data are compositional data, in which only the ratios between components matter, and not their absolute abundances. This is reflected in the fact that 40Ar-39Ar-36Ar data can be renormalised to unity and plotted on a ternary diagram (left). There is a one-to-one mapping between this so-called ‘simplex’ and Euclidean logratio space (right).


3 Multi-collector data

To illustrate the calculations in the remainder of this paper, consider the following sequence of analyses: b1 (first blank), u1 (first sample), s1 (first age standard), u2 (second sample), b2 (second blank), s2 (second standard), s3 (third standard) and b3 (third blank). In a multicollector mass spectrometer, each of the five argon isotopes appearing in Equation 5 are monitored simultaneously through time (t) and can be cast into an [n × 5] matrix format, with n the number of integrations (i.e. t = {t1,t2,...,tn}):

         ⌊                                                       ⌋
           36Ar(x,t1) 37Ar(x,t1) 38Ar(x,t1) 39Ar(x,t1) 40Ar(x,t1)
         || 36Ar(x,t2) 37Ar(x,t2) 38Ar(x,t2) 39Ar(x,t2) 40Ar(x,t2)||
M (x,t) = |⌈     ..          ..          ..          ..          ..    |⌉
           36   .     37   .     38   .     39   .     40   .
             Ar(x,tn)   Ar(x,tn)   Ar(x,tn)   Ar(x,tn)   Ar(x,tn)
(12)

where ‘x’ stands for ‘blank’, ‘sample’ or ‘standard’. The same formulation can be used for the interference monitors (particularly Ca) but further discussion of these will be deferred to Section 8 and Appendix A. Because the measurements are done simultaneously on all five detectors, any random variation in, say, the filament voltage or trap current will simultaneously affect all signals, resulting in correlated residuals. The blank correction is made by subtracting the time-resolved signal of the nearest blank measurement (b) from that of the analysis (x), resulting in a new matrix B(x,b,t):

          ⌊ 36  b       37  b       38   b      39  b       40  b      ⌋
          | 36Arb(x,t1)  3A7r b(x,t1)  38Ar (bx,t1) 39Arb(x,t1)  40Arb(x,t1) |
B(x,b,t) = ||   Ar (.x,t2)   Ar .(x,t2)   Ar (.x,t2)   Ar (.x,t2)    Ar (.x,t2) ||
          ⌈      ..          ..           ..           ..           ..     ⌉
            36Arb(x,tn)  37Arb(x,tn)  38Arb(x,tn) 39Arb(x,tn)  40Arb(x,tn)
(13)

with

iArb(x,tj) = iAr(x,tj) − iAr(b,tj)
(14)

for i = {36, 37, 38, 39, 40} and j = {1, ..., n}. Our goal is to extract 4-element vectors of logratios from these [n×5] matrices of blank corrected mass spectrometer signals, taking into account any correlated errors. The easiest but by no means only way to achieve this is by forming the logratios prior to regression, yielding an [n×4] matrix for each analysis:

          ⌊  [36Arb(x,t)]   [37Arb(x,t )]   [38Arb(x,t )]   [ 39Arb(x,t )]⌋
          | l[40Arb(x,t11)]  l[40Arb(x,t11)]  l[40Arb(x,t11)]  l[ 40Arb(x,t11)]|
          || l 36Arb(x,t2)   l 37Arb(x,t2)  l 38Arb(x,t2)  l  39Arb(x,t2) ||
L (x,b,t) = ||  40Arb(x,t2)     40Arb(x,t2)    40Arb(x,t2)     40Arb(x,t2) ||
          |⌈      ...             ...            ...            ...     |⌉
            l[36Arb(x,tn)]  l[37Arb(x,tn)]  l[38Arb(x,tn)] l[39Arb(x,tn)]
              40Arb(x,tn)    40Arb(x,tn)    40Arb(x,tn)    40Arb(x,tn)
(15)

where ‘l’ stands for ‘natural log’ and 40Ar is used as a common denominator for all the ratios denoted by ‘m’ in Equation 5. We thus obtain five time-resolved logratio matrices, one for each run in the analysis sequence. These five matrices can be assembled into one [n × 20] matrix, which is naturally partitioned into three groups by the blanks.

G(t) = [L (u1,b1,t) L(s1,b1,t) | L(u2,b2,t) L(s2,b2) | L (s3,b3,t)] = [g1|g2|g3]
(16)

where the first group (g1) consists of sample u1 and standard s1, which share blank b1; the second group (g2) consists of sample u2 and standard s2, which share blank b2; and the third group consists of standard s3, which is the only analysis using blank b3. It is reasonable to expect the blank-corrected logratio signals to be correlated within each group, but uncorrelated between groups. We therefore extrapolate the logratio signals to t=0 (‘time zero’) in blocks, and concatenate the resulting logratio intercepts into a single 20-element vector:

X = [X(g1) X (g2) X (g3)]
(17)

with X(gi) the vector of logratio intercepts of the ith group, obtained by joint (non)linear regression. The [20 × 20] covariance matrix of X is given by:

      ⌊               ⌋
        Σg1  08,8  08,4
ΣX  = ⌈ 08,8  Σg2  08,4 ⌉
        04,8  04,8  Σg3
(18)

where Σgi is the covariance matrix of the ith group’s intercepts and 0i,j denotes a zero matrix of size [i×j]. One well known problem with the logratio transformation is the handling of zero or negative values. In the context of argon mass spectrometry, this occurs in one of two situations: (a) 36Ar (and 38Ar) in the atmospheric correction of extremely clean samples and (b) 37Ar in the Ca-interference correction of ‘expired’ samples. The zero value problem can be avoided by performing generalised linear regression of the ratios (using a logarithmic link function to ensure positive intercepts, Nelder and Wedderburn1972), or to cast the regression problem into a more sophisticated maximum likelihood form (Wood2015). A comprehensive discussion of these alternative methods falls outside the scope of the present paper and will be deferred to a future publication.

4 ‘Peak-hopping’ data

In single collector mass spectrometers, the various argon isotopes cannot be monitored simultaneously, but must be measured separately. This is achieved by separately scanning (‘hopping’) over the mass range of the argon isotopes by varying the field strength of the mass analyser. Thus, each mass has its own time scale ti, for i = 36, 37, 38, 39 and 40, resulting in a set of five time resolved data vectors M(x,i,ti) for each run x:

           ⌊          ⌋
             iAr(x,ti1)
           || iAr(x,ti2) ||
M (x,i,ti) = |⌈     ..    |⌉
             i   . i
              Ar(x,tn)
(19)

Because the five isotope signals are measured at different times, we can safely assume their residual noise to be uncorrelated. Again, blank correction is done in time-resolved mode, but separately for each isotope. This results in five (one for each run) times five (for each isotope) n-element ratio vectors:

            ⌊  [i     i   i     i ] ⌋
            | l[iAr(x,t1i)− iAr(b,t1i)] |
L(x,b,i,ti) = || l  Ar(x,t2)−.  Ar(b,t2)  ||
            ⌈  [        ..         ] ⌉
              l iAr (x,tin)− iAr(b,tin)
(20)

These vectors are assembled into five [n × 5] matrices, each of which is partitioned into three groups according to the shared blank corrections:

         [                                                          ]
G (i,ti) = L (u1,b1,i,ti) L(s1,b1,i,ti) | L (u2,b2,i,ti) L(s2,b2,i,ti) | L (s3,b3,i,ti) = [gi1|gi2|g3i]
(21)

Joint regression to t=0 yields a 5-element vector of log-intercepts for each isotope:

Z(i) = [Z(g1,i) Z (g2,i) Z (g3,i)]
(22)

with [5 × 5] covariance matrices

       ⌊  Σi  02,2 02,1⌋
Σ    = ⌈  02g,12 Σi   02,1⌉
  Z(i)     01,2 01g,22 Σi
                     g3
(23)

where Σgji is the covariance matrix of the jth group’s iAr intercepts. Next, we bring the ratio-intercept data for all five isotopes together into a single 25-element vector

Z = [Z (36) Z(37) Z (38) Z(39) Z(40)]
(24)

with [25 × 25] covariance matrix

     ⌊                                    ⌋
     | ΣZ (36)  05,5    05,5    05,5    05,5  |
     ||  05,5   ΣZ (37)   05,5    05,5    05,5  ||
ΣZ = |⌈  05,5    05,5  ΣZ (38)   05,5    05,5  |⌉
        05,5    05,5    05,5   ΣZ(39)   05,5
        05,5    05,5    05,5    05,5   ΣZ(40)
(25)

Finally, we form 20 logratios with the following matrix operation:

X = Z JX
(26)

The associated [20 × 20] covariance matrix is given by:

Σ   = J′ Σ  J
  X    X  Z  X
(27)

with JX the [25 × 20] Jacobian matrix of the subtraction operation and JXits transpose:

     ⌊                          ⌋
     | 15,5  05,5  05,5  05,5  − 15,5 |
J′X = |⌈ 05,5  15,5  05,5  05,5  − 15,5 |⌉
       05,5  05,5  15,5  05,5  − 15,5
       05,5  05,5  05,5  15,5  − 15,5
(28)

where 1i,i is an [i × i] identity matrix. We have now cast the raw mass spectrometer data in a common logratio format X (through either Equation 17 or 26) and associated covariance structure ΣX (Equation 18 or 27). From here on, multicollector and peak-hopping data can be treated on an equal footing.

5 Detector calibration

The different ion detectors in a multicollector mass spectrometer do not necessarily respond equally to ion beams of equal mass and size. The measured ratio of the beam intensities at t=0 will therefore not necessarily equal the true isotopic ratio. This issue obviously does not occur in single collector instruments. Although the latest generation of multicollector noble gas mass spectrometers quantify the relative sensitivities internally through an electronic detector intercalibration, this section describes a data reduction protocol for a conventional (‘analog’) detector calibration. Suppose that there are five detectors, one for each argon isotope, and denote these by d[36], d[37], d[38], d[39] and d[40]. The relative sensitivities of detectors d[36] and d[40] can be quantified by comparing the measured 40Ar/36Ar intensity ratio of an air shot with the known atmospheric ratio, as part of the mass fractionation correction (Section 6). The relative sensitivities of the remaining detectors, d[37]-d[40], on the other hand, are calibrated by steering a fixed 40Ar beam from an air tank across them. The resulting signals of this ‘peak hopping’ experiment are extrapolated to t=0 using the methods described in Section 4, resulting in four log-intercepts and their variances. No blank corrections are needed because we are only interested in the total amount of gas present in the mass spectrometer and not in the air composition itself. If the calibration experiment is repeated multiple times, then the measurements can be combined by taking the arithmetic mean of the logs (Section 13). To apply the detector calibration correction, we simply add the difference of the log-intercepts to the data, in matrix form. First, we append the log-intercepts of the calibration data to the sample vector.

 ∗
X  = [X  Z(d[37]) Z (d[38]) Z(d[39]) Z (d[40])]
(29)

with [24 × 24] covariance matrix ΣX:

     ⌊  ΣX      020,1        020,1        020,1        020,1    ⌋
     |  01,20 σ[Z(d[37])]2       0          0           0      |
Σ∗ = ||  01,20      0      σ[Z(d[38])]2      0           0      ||
 X   |⌈  01,20      0           0      σ[Z(d[39])]2      0      |⌉
        01,20      0           0          0       σ[Z(d[40])]2
(30)

where X is a 20-element vector of sample and standard measurements (Equation 17) and ΣX its covariance matrix (Equation 18), Z(d[i]) indicates the log intercept of 40Ar measured by detector d[i] at ‘time zero’, and σ[Z(d[i])] is its standard error. Then the detector calibrated data (C) and their [20 × 20] covariance matrix (ΣC) are obtained by:

C = X∗ J
        C
(31)

and

ΣC = J′ Σ∗ JC
      C  X
(32)

respectively, where JC is the [24 × 20] Jacobian matrix of the detector calibration and JCis its transpose:

     ⌊                             ⌋
       14,4 04,4 04,4  04,4  04,4  J∗C
  ′  || 04,4 14,4 04,4  04,4  04,4  J∗C ||
JC = || 04,4 04,4 14,4  04,4  04,4  J∗C ||
     ⌈ 04,4 04,4 04,4  14,4  04,4  J∗C ⌉
       04,4 04,4 04,4  04,4  14,4  J∗C
(33)

with

     ⌊                ⌋
        0   0   0   0
J∗C = || − 1  0   0   1 ||
     ⌈  0  − 1  0   1 ⌉
        0   0   − 1 1
(34)

Note that, if all the measurements (samples, age standards and interference monitors) use the same detector calibration, then the associated analytical uncertainties cancel out in the age calculation (Section 12) and we can set σ[Z(d[i])]2 = 0 i in Equation 30.

6 Mass fractionation

The five argon isotopes of interest span a mass range of 10%. The sensitivity of both single- and multicollector instruments varies with atomic mass, and significant errors can occur if the resulting ‘mass fractionation’ is uncorrected for. The mass fractionation factor can be quantified by comparing the measured signal ratios of an air shot with its known isotopic ratio (298.56 ± 0.31, Lee et al.2006). For multicollector instruments, each detector has its own mass fractionation correction factor. For detectors d[37], d[38] and d[39], these are obtained by peak hopping between masses 36 and 40. For d[40] and d[36], we can quantify the fractionation by directly monitoring the 36Ar/40Ar-ratio in multicollection mode. The exponential form of the kinetic isotope fractionation correction (Young et al.2002) conveniently reduces to a linear equation in a logratio context:

 [iAr ]   [ iAr |d[i]]    l[i]− l[j] (       [40Ar ] )
l jAr- = l jAr|d[j]  + l[40]−-l[36] A (j)+ l 36Ar
                                               a
(35)

where iAr|d[j] stands for the iAr signal measured on detector j and A(j) is the ‘time zero’ intercept of l[36Ar|d[j]]
 40Ar|d[j]a, except if j = 40 on a multicollector instrument, in which case A(j) is the ‘time zero’ intercept of l[36Ar|d[36]]
 40Ar|d[40]a. To apply Equation 35, we append the air shot data and the true air ratio to the calibration-corrected logratio intercepts:

     [         [40Ar ] ]
C ∗ = C A (40) l 36--
                 Ar  a
(36)

whose [22 × 22] covariance matrix ΣC can be written as:

     ⌊                      ⌋
 ∗   ⌈  ΣC     020,1 2  020,1 ⌉
ΣC =   01,20  σ[A(40)]   0
       01,20      0      0
(37)

Note that Equation 37 does not specify the analytical uncertainty of the atmospheric reference ratio. This is because any uncertainty resulting from an incorrect air-ratio at this point will cancel out during the atmospheric argon correction (Section 10). Recasting Equation 35 in matrix form, the fractionation correction of the sample and fluence measurements can be written as:

     ∗
F = C  JF
(38)

with [20 × 20] covariance matrix

Σ  = J′ Σ∗ J
 F    F  C  F
(39)

where JF is the [22 × 20] Jacobian matrix of the mass fractionation correction and JFis its transpose:

     ⌊                           ∗ ⌋
     | 14,4 04,4 04,4  04,4  04,4  JF∗ |
  ′  || 04,4 14,4 04,4  04,4  04,4  JF∗ ||
JF = |⌈ 04,4 04,4 14,4  04,4  04,4  JF∗ |⌉
       04,4 04,4 04,4  14,4  04,4  JF∗
       04,4 04,4 04,4  04,4  14,4  JF
(40)

with

     ⌊ − 1.000  − 1.000 ⌋
     | − 0.740  − 0.740 |
J∗F = |⌈ − 0.487  − 0.487 |⌉
       − 0.240  − 0.240
(41)

7 Decay corrections

Two of the five argon isotopes of interest are radioactive: 37Ar (t12 = 34.95 ± 0.08 days, Renne and Norman2001) and 39Ar (t12 = 269 ± 3 years, Stoenner et al.1965). A correction is required for the loss of these isotopes during the time elapsed between irradiation and analysis:

l[iAr ]∘ = l[iAr](τ)+ r(λi,τ)
(42)

where l[iAr] is the total amount of isotope i formed during irradiation, l[iAr](τ) is the amount remaining at a time τ after the end of the irradiation and r(λi) is the amount lost due to radioactivity when the decay constant is λi. Using a similar approach to Wijbrans and McDougall (1986), r(λi) can be calculated as:

          ⌊        ⌋    ⌊                           ⌋
           ∑             ∑  Pj (  1          1    )
r(λi,τ) = l⌈   PjΔtj⌉ − l⌈   λ-- eλiΔτj − eλi[Δτj+Δtj] ⌉
            j             j  i
(43)

where Pj is the power and Δtj the duration of the jth irradiation interval and Δτj is the time elapsed between the end of the jth irradiation segment and τ. At this point it is important to merge the data reduction pathways for the samples and fluence monitors with those of any co-irradiated K-glass and Ca-salt. This is because they are all affected by the same decay constant uncertainties, resulting in correlated errors. However, in this Section we will, for the sake of simplicity, assume that [        ]
36Ar∕37Arca, [         ]
 39Ar∕37Arca and [         ]
 39Ar∕40Ark have been obtained from elsewhere and do not need to be corrected for radioactive decay. For completeness, further details about the joint analysis of co-irradiated interference monitors with the sample are given in Appendix A. To apply the decay correction to the samples and fluence monitors, we first concatenate all the decay corrections into one 5-element vector:

r(i) = [r(λi,τ[u1]) r(λi,τ[s1]) r(λi,τ[u2]) r(λi,τ[s2]) r(λi,τ[s3])]
(44)

The [5 × 5] covariance matrix of which is given by:

        ′       2
Σr (i) = Jr(i) σ(λi) Jr(i)
(45)

where σ(λi) is the standard error of the iAr decay constant, and Jr(i) is the Jacobian matrix:

    [                                                       ]
     ∂r(λi,τ[u1]) ∂r(λi,τ[s1]) ∂r(λi,τ[u2])-∂r(λi,τ[s2]) ∂r(λi,τ[s3])
Jr =     ∂λi        ∂λi        ∂ λi        ∂λi        ∂λi
(46)

with the partial derivatives given by:

                  [                               ]∕       [                      ]
∂r(λi,τ-[x])  ∑   Pj-1-+-λiΔτj[x]   1+-λi(Δ-τj[x]+-Δtj)-   ∑     ---1---   -----1------
   ∂λi    =     λi   eλiΔτj[x]  −   eλi(Δ τj[x]+Δtj)         Pj eλiΔτj[x] − eλi(Δτj[x]+Δtj)
             j                                         j
(47)

Next, we append the vector of 10 decay corrections to the 20 fractionation-corrected logratio intercepts:

 ∗
F  = [F r(37) r(39)]
(48)

with [30 × 30] covariance matrix

     ⌊  ΣF    020,5  020,5  ⌋
Σ∗ = ⌈ 05,20  Σr(37)    0   ⌉
 F     05,20    0    Σr(39)
(49)

The decay correction can then be cast into matrix form as

D = F∗ J
        D
(50)

yielding a 20-element vector with covariance matrix

       ′  ∗
ΣD  = JD ΣF JD
(51)

using the [30 × 20] Jacobian matrix JD and its transpose JD:

  ′         ∗     ∗
J D = [120,20 JD(37) JD(39)]
(52)

with

      ⌊                               ⌋
         J∗D∗(i)  04,1   04,1  04,1   04,1
      ||  04,1   J∗D∗(i)  04,1  04,1   04,1 ||
J∗  = ||  04,1   04,1   J∗∗   04,1   04,1 ||
 D(i)  |⌈  0     0     0D(i)  J∗∗    0   |⌉
         04,1   04,1   04,1  0D(i)  J4∗,∗1
          4,1    4,1    4,1   4,1   D(i)
(53)

where

        ⌊   ⌋             ⌊   ⌋
        | 0 |             | 0 |
J∗D∗(37) = |⌈ 1 |⌉ and J∗D∗(39) = |⌈ 0 |⌉
          0                 0
          0                 1
(54)

8 Interference corrections

The 40Ar/39Ar-method pairs the natural radioactive decay of 40K to 40Ar with the synthetic activation of 39K to 39Ar. Unfortunately, neutron activation produces not only 39Ar but a host of other Ar-isotopes as well. The most important reactions are (McDougall and Harrison1999):

K :   39K(n,p)39Ar
  40K(n,p)40Ar
Ca :   40Ca(n,nα)36Ar
  40Ca(n,α)37Ar
  42Ca(n,α)39Ar
Cl :   35Cl(n,γ)36Cl β−
−−→36Ar
  37Cl(n,γ)38Cl  −
−β−→38Ar

The first five of these reactions can be characterised by mass spectrometric analysis of K-glass (40Ar/39Ar ratio) and Ca-salt (36Ar/37Ar and 39Ar/37Ar ratios). These ratios are directly incorporated into Equation 5 (parameters a, b and f). The chlorine decay products, on the other hand, are generally calculated from the independently determined and reactor-specific 36Cl/38Cl-production ratio and will be discussed in Section 9. If the K- and Ca-interference corrections are based on externally determined values, then we compile these with the decay-corrected sample and fluence measurements for further processing in Section 10:

I = [D  D(ca) D (k )]
(55)

where, using the notation of Section 7, D(ca) is a 2-element vector containing the decay-corrected 36Ar/37Ar- and 39Ar/37Ar-logratios of neutron-activated Ca, and D(k) is the 40Ar/39Ar-logratio of neutron-activated K. The corresponding [23 × 23] covariance matrix is given by

     ⌊  ΣD    020,2   020,1  ⌋
Σ  = ⌈ 02,20  ΣD(ca)  02,1  ⌉
  I    01,20   01,2   σ2
                     D (k)
(56)

After which we can proceed to Section 9 of this paper. If, on the other hand, Ca and K interferences are quantified by co-irradiated Ca-salts and K-glass, then we can explicitly include the resulting mass spectrometer uncertainties into the error propagation. Further details of this are provided in Appendix A. In summary, the vector I, obtained from either Equation 55 or Appendix A, contains all the information required to solve Equation 5 except for factor ‘c’, which is discussed next.

9 Cl-decay

In contrast with the K- and Ca-interferences, which can be directly characterised by mass spectrometric analysis of co-irradiated materials, the Cl-interference on 36Ar is generally calculated from an independently determined and reactor-specific 36Cl/38Cl-production ratio (Foland et al.1993Renne et al.2008). Let G(x) be the logratio of the chlorine decay products (i.e., l[         ]
 36Ar∕38Ar) in sample (or fluence monitor) x. Using the approach of Wijbrans and McDougall (1986) to account for the radioactive decay of Cl to Ar, we obtain:

       [36  ]
G(x) = l 38Cl + g(τ[x])
          Cl
(57)

with

       [   ∑     ( −λ [Δτ+Δt ]   −λ Δτ )]
g(τ) = l 1+--j-Pj-e--36-∑-j--j-−-e--36--j-
                    λ36  jPjΔtj
(58)

where λ36 is the 36Cl decay constant and τ, Δτi and Δti are as defined in Section 7. The decay corrections can be compiled into a single five-element vector

G = [G (u1) G(s1) G (u2) G(s2) G (s3)]
(59)

whose [5 × 5] covariance matrix is given by:

       [  ( [36Cl])2       ]
ΣG = J′G  σ l 38Cl     0    JG
              0    σ (λ36)2
(60)

with

    [                                                                ]
J =        1            1            1           1            1
 G    ∂G (u1)∕∂λ36  ∂G (s1)∕∂λ36  ∂G(u2)∕∂λ36 ∂G (s2)∕∂λ36  ∂G (s3)∕∂λ36
(61)

where the partial derivatives are given by:

        ∑     [                                                       ]
∂G-(x)  --j-Pj-(1-+-λ36Δ-τj[x])e−-λ36Δτj −-(1+-λ36[Δ-τj[x]+-Δtj])e−λ36(Δ-τj[x]+Δtj)
 ∂λ36 =           λ36∑j Pj [λ36Δtj +e− λ36(Δ τj[x]+Δtj) − e− λ36Δτj[x]]
(62)

Note that the Cl-interference correction implemented in Equation 8 does not account for the presence of atmospheric 38Ar and the production of 38Ar from K. Doing so is straightforward but adds considerably more complexity to Equation 5 (Appendix B).

10 40Ar39ArK

After all the preprocessing discussed in the previous sections, we have finally gathered all the ratios required to solve Equation 5. To this end, we compile all the information obtained thus far into a single vector of logratios

    [   [40Ar]    ]
U =  I l 36Ar   G
              a
(63)

and its [29 × 29] covariance matrix

     ⌊  Σ        0         0    ⌋
     |    I    ( [2403,1] )2   23,5 |
ΣU = ⌈ 01,23  σ  l 36AArr-a    01,5 ⌉
       05,23       05,1       ΣG
(64)

To simplify the notation in the remainder of this Section, it is useful to permute U and ΣU so that the Cl-interference data (G) are interspersed with the samples and fluence monitors:

U∗ = U P and Σ ∗ = P U P
              U
(65)

where P is the [29 × 29] permutation matrix

    ⌊                                                       ⌋
      14,4 04,1  04,4  04,1  04,4  04,1  04,4  04,1  04,4  04,1  04,4
    || 04,4 04,1  14,4  04,1  04,4  04,1  04,4  04,1  04,4  04,1  04,4 ||
    || 04,4 04,1  04,4  04,1  14,4  04,1  04,4  04,1  04,4  04,1  04,4 ||
    || 04,4 04,1  04,4  04,1  04,4  04,1  14,4  04,1  04,4  04,1  04,4 ||
    || 04,4 04,1  04,4  04,1  04,4  04,1  04,4  04,1  14,4  04,1  04,4 ||
P = || 04,4 04,1  04,4  04,1  04,4  04,1  04,4  04,1  04,4  04,1  14,4 ||
    || 01,4  1   01,4   0   01,4   0   01,4   0   01,4   0   01,4 ||
    || 01,4  0   01,4   1   01,4   0   01,4   0   01,4   0   01,4 ||
    || 01,4  0   01,4   0   01,4   1   01,4   0   01,4   0   01,4 ||
    ⌈ 01,4  0   01,4   0   01,4   0   01,4   1   01,4   0   01,4 ⌉
      01,4  0   01,4   0   01,4   0   01,4   0   01,4   1   01,4
(66)

Next we convert the logratio vector U into a vector of 30 ratios

W  = [V (u1) V (s1) V (u2) V (s2) V (s3)]
(67)

where

V(x) = [a(x ) b(x ) c(x) d(x) e(x) f(x)]
(68)

with a-f as defined in Equations 6-11. f(x) is the same for all analyses in this example but may vary between samples when combining different irradiations. W is calculated in matrix form by

          ∗
W  = exp[U  JV]
(69)

with JV the [29 × 30] Jacobian matrix:

     ⌊                         ⌋
        J∗V   05,6 05,6 05,6 05,6
     ||  05,6  J∗V  05,6 05,6 05,6 ||
     ||  05,6  05,6  J∗V  05,6 05,6 ||
JV = ||  05,6  05,6 05,6  J∗V  05,6 ||
     ⌈  05,6  05,6 05,6 05,6  J∗V ⌉
        J∗V∗  J∗V∗  J∗V∗ J ∗∗V  JV∗∗
(70)

where

     ⌊                  ⌋          ⌊                 ⌋
     | 1  0  0  0  0  0 |            0  1  0  0  0  0
  ∗  || 0  1  0  0  1  0 ||      ∗∗  || 0  0  0  0  1  0||
JV = |⌈ 0  0  1  0  0  0 |⌉ and JV = ⌈ 0  0  0  0  0  1⌉
       0  0  0  1  0  0              1  1  1  0  0  0
       0  0  1  0  0  0
(71)

The [30 × 30] covariance matrix of W is obtained by

       ′  ∗
ΣW  = JW ΣU JW
(72)

where the [29 × 30] Jacobian JW is given by

      ⌊  ∗                     ⌋
      | JW   05,∗6  05,6  05,6  05,6|
      || 05,6  JW   05,∗6  05,6  05,6||
JW  = || 05,6  05,6  JW   05,∗6  05,6||
      |⌈ 05,6  05,6  05,6  JW   05,∗6|⌉
        05∗,6∗  05,∗6∗  05,∗6∗  05,∗6∗  JW∗∗
        JW   JW   JW   JW   JW
(73)

with

     ⌊  a 0  0  0  0  0 ⌋          ⌊                  ⌋
     |  0 b  0  0  e  0 |            0  b  0  0  0  0
 ∗   ||  0 0  c  0  0  0 ||      ∗∗  || 0  0  0  0  e  0 ||
JW = |⌈  0 0  0  d  0  0 |⌉ and JW = ⌈ 0  0  0  0  0  f ⌉
        0 0  c  0  0  0              a  b  c  0  0  0
(74)

The five element vector R of 40Ar39ArK-ratios is calculated with Equation 5:

R = [R (u1) R(s1) R (u2) R(s2) R(s3)]
(75)

and its [5 × 5] covariance matrix is obtained by

ΣR  = J′ ΣW JR
       R
(76)

where JR is the [30 × 5] Jacobian matrix and JRis its transpose

     ⌊                                      ⌋
       J∗R(u1)   01,6     01,6    01,6     01,6
 ′   ||  01,6   J∗R(s1)   01,6    01,6     01,6  ||
JR = ||  01,6    01,6   J∗R(u2)   01,6     01,6  ||
     ⌈  01,6    01,6     01,6   J∗R(s2)   01,6  ⌉
        01,6    01,6     01,6    01,6   J∗R (s3)
(77)

with

        [   − 1        1          1     a(x)− b(x)− c(x) − 1 1 − a(x)+ b(x)+ c(x)   ]
J∗R(x) =  -------------------- ------------------------2------------------2---- − 1
         d(x)− e(x) d(x)− e(x) d(x) − e(x)     [d(x)− e(x)]         [d(x) − e(x)]
(78)

11 J-factors

The parameter J quantifying the production of 39Ar from 39K in the age equation is determined by analysing the argon composition of a co-irradiated fluence monitor with accurately known K-Ar age (Ts). This composition may vary across the irradiation stack due to neutron flux gradients in the reactor, which can be quantified by analysing several fluence monitors interspersed with the samples at known positions. The most appropriate J-factor for each sample is then obtained by simple linear interpolation:

       λ40Ts
J(x) = e----−-1-
        R(s|x)
(79)

where R(s|x) denotes the 40Ar39ArK-ratio of the fluence monitors interpolated to the position of sample x (which is henceforth referred to as p[x]). Applying this procedure to our two sample - three monitor case study, we form a four-element vector of sample ratios and interpolated fluence monitor ratios:

Y = [R(u1) R (u2) R(s|u1) R(s|u2)] = R JY
(80)

with [4 × 4] covariance matrix

      ′
ΣY = JY ΣR JY
(81)

where R is the vector of 40Ar39ArK-ratios for the samples and fluence monitors (Equation 75), JY is the [5 × 4] Jacobian matrix and JY is its transpose. Suppose that sample u1 sits between monitors s1 and s2 in the irradiation stack, and u2 sits between monitors s2 and s3. Then

     ⌊ 1     0      0     0       0     ⌋
     | 0     0      1     0       0     |
J′Y = || 0  p[u1]−-p[s1]  0  p[s2]−p[u1]   0     ||
     ⌈    p[s2]− p[s1]     pp[s[u22]−]−pp[[s1s2]]p[s3]−p[u2] ⌉
       0     0      0  p[s3]−p[s2]p[s3]−p[s2]
(82)

Finally, we use Equation 79 to generate a five-element vector of sample 40Ar39ArK-ratios, their respective J-factors, and the 40K decay constant:

Q = [R (u1) R(u2) J (u1) J(u2) λ40]
(83)

with [5 × 5] covariance matrix

         ⌊                    ⌋
       ′ ⌈ ΣY    20       0   ⌉
ΣQ  = JQ    0   σ (λ40)   20     JQ
            0     0     σ (Ts)
(84)

where JQ is the [6 × 5] Jacobian matrix and JQis its transpose:

     ⌊ 1  0   0      0      0      0    ⌋
     | 0  1   0      0      0      0    |
 ′   ||     1−eλ40Ts       Tseλ40Ts-λ40eλ40Ts ||
JQ = || 0  0R (s|u1)2   0λ Ts R(sλ|u1)Ts  R(s|uλ1)Ts ||
     ⌈ 0  0   0   1R−e(s|u402)2 TsRe(s4|0u2)-λ40Re(s|u402) ⌉
       0  0   0      0      1      0
(85)

The decay constant λ40 is included into Equation 83 because this parameter appears in both the expression for J (Equation 79) and the age equation (Equation 3), resulting in correlated errors.

12 Solving the age equation

The 40Ar/39Ar-ages of samples u1 and u2 are calculated by plugging the relevant items of vector Q into Equation 3, resulting in a 2-element vector T

T = [T (u1) T(u2)]
(86)

with [2 × 2] covariance matrix

      ′
ΣT = JT ΣQ JT
(87)

where JT is the [5 × 2] Jacobian matrix:

     [ -----J(u1)----      0       ----R(u1)-----      0      − l[1+J-(u1)R-(u1)]-]
J′T =   λ40[1+J (u1)R (u1)]     J(u2)     λ40[1+J(u1)R(u1)]    R (u2)      l[1+J λ(u2420)R (u2)]
             0      λ40[1+J(u2)R(u2)]      0      λ40[1+J-(u2)R-(u2)]− ----λ240-----
(88)

13 (weighted) mean ages

Given a vector of N age measurements (T = [T(u1) T(u2) ... T(uN)]), we can calculate the arithmetic mean age Ta as:

T¯a = (T 1N,1)∕N
(89)

with standard error

 2
σ (¯Ta) = (11,N ΣT 1N,1)∕N
(90)

Alternatively, to calculate the error-weighted mean Tw, first calculate its variance:

 2 ¯    (     − 1    )− 1
σ (Tw ) = 11,N ΣT  1N,1
(91)

then

T¯w = σ2(¯Tw)(T Σ −T11N,1)
(92)

The MSWD (‘Mean Square of the Weighted Deviates’, also known as ‘reduced Chi-square statistic’ outside geology) is a measure of the ratio of the observed scatter of the data points (T[ui]) around the mean value (T) to the expected scatter from the assigned errors (ΣT):

            1
M SW  D = -----[T − ¯T] Σ−T 1[T − T¯]′
          N − 1
(93)

If MSWD>1, then the samples are said to be ‘overdispersed’ with respect to the analytical uncertainty. This commonly occurs in very precise datasets, which have sufficient power to resolve minute levels of sample heterogeneity. In this case, the geologically meaningful levels of heterogeneity can be quantified using a ‘mixed effects’ model with two sources of analytical uncertainty:

                   2    2
T[ui] ∼ N [T¯,σ(T[ui]) + ζ]
(94)

where N[a,b] stands for “the Normal distribution with mean a and variance b”, and ζ2 is the ‘overdispersion’ (Vermeesch2010). Equation 94 can be solved by the method of maximum likelihood, which simultaneously estimates the average, its standard error, and the overdispersion.

14 Ar-Ar_Redux

The revised data reduction procedure outlined in this paper revisits every aspect of the 40Ar/39Ar method. Unfortunately, the matrix format of the calculations is incompatible with existing data reduction platforms such as ArArCalc (Koppers2002). A new computer code named Ar-Ar_Redux was developed to solve this problem and facilitate the adoption of the methods described herein. A prototype version of Ar-Ar_Redux currently exists as a package in the R programming environment, which is an increasingly popular open source alternative to Matlab, available free of charge on any operating system at http://r-project.org. A standalone program with graphical user interface is in development for future release. ‘Ar-Ar_Redux’ derives its name from ‘U-Pb_Redux’, which is a similar program developed by the U-Pb dating community (McLean et al.2011Bowring et al.2011). Both programs use a similar matrix formulation and, although U-Pb_Redux currently does not employ a logratio transformation, future versions of it will. The R-version of Ar-Ar_Redux can be downloaded free of charge from the ‘Comprehensive R-Archive Network’ (CRAN, http://cran.r-project.org). Appendix C gives a brief introduction to Ar-Ar_Redux, with further details provided at http://redux.london-geochron.com. The latter website will also host the standalone version of the program when it is ready for public release. Currently, Ar-Ar_Redux accepts input files that are compatible with the ARGUS-VI multicollector instrument, but other input formats can easily be implemented as well. Ar-Ar_Redux is intended to be a community-driven software platform, which can evolve to accommodate the demands and expectations of 40Ar/39Ar practitioners, and the reader is invited to contact the author with any questions or requests. The program is bundled with a real dataset, which was kindly provided by Prof. David Phillips of the University of Melbourne.

15 Discussion and conclusions

One might wonder how much difference the revised data reduction workflow makes compared to currently used procedures. The answer to this question depends on the particular details of the sample of interest. For example:

  • Error correlations are stronger when several samples share the same blank than when each sample has its own blank.
  • Large interference corrections result in strong error correlations.
  • Multicollector data are more strongly correlated than ‘peak hopping’ data.
  • Analysing co-irradiated interference monitors yields stronger error correlations than using externally provided interference corrections.

Regarding the latter two examples, it is important to note that correlated errors should not necessarily be considered undesirable, as long as they are properly quantified. It is only when covariances are ignored that uncertainties are overestimated, potentially significant age differences are blurred out and geologically meaningful information is lost. Experience tells that the covariance terms can be very substantial. For the test data provided with Ar-Ar_Redux, error correlations (defined as ρ(x,y) = cov(x,y)/[σ(x)σ(y)]) between aliquots of the same sample are on the order of 0.9.

Renne et al. (1998) make the distinction between ‘internal’ and ‘external’ errors. ‘internal errors’ can be conceptually defined as the natural variability that would arise if the same sample were dated multiple times under the same experimental conditions. ‘external’ errors include the systematic effects of decay constant uncertainty, the K/Ar ratio of the age standard, the air ratio etc. Renne et al. (1998) point out that “comparison of two different 40Ar/39Ar dates based on the same standard may legitimately ignore uncertainties in K-Ar data, decay constants, as well as all intercalibration factors common to both dates”. However, when comparing a 40Ar/39Ar-age with, say, a zircon U/Pb age, “it is important to consider all sources of systematic error in data from both radioisotopic systems”. Thus, great care must be taken which sources of uncertainty should or should not be included in the error propagation. In practical terms, this results in different analytical forms of the error propagation depending on the situation. This added complexity disappears entirely when using the methods presented in this paper. By processing the data in matrix form and explicitly taking into account covariances, the internal and external errors are jointly considered, with the latter corresponding to the off-diagonal terms of the covariance matrix. Revisiting Renne et al. (1998)’s two scenarios, we find that the difference between two 40Ar/39Ar dates based on the same standard may appear to be statistically insignificant compared to their respective variances, but statistically significant when the covariance terms are considered (Figure 2).

This paper has revisited many but not all aspects of 40Ar/39Ar data reduction. For example, it has not discussed isochrons, in which linear regression is used to deconvolve the radiogenic and inherited argon components without the need to assume an atmospheric composition for the latter. Although the least squares algorithms which are currently used for this purpose do take into account error correlations between the x- and y-variables (e.g., York1969), they ignore the covariance between different samples. Similarly, thermal modelling is done by jointly considering multiple analyses and finding best-fitting (‘Arrhenius’) trends to them. Current fitting algorithms do not account for the significant error correlations that exist between subsequent heating steps in a diffusion experiment. The covariant structure of linear regression naturally follows from the covariant age structure represented by Equations 86 and 87, but a detailed discussion of this will be deferred to a forthcoming publication.


Figure 2: A synthetic yet realistic example of two replicate age estimates of the same sample (T1 = 99 Ma and T2 = 101 Ma) plotted against each other as an error ellipse. Ignoring the covariances, the two dates appear to agree within two standard errors. Taking into account the off-diagonal terms of the covariance matrix (ΣT), however, reveals that the two samples are overdispersed with respect to the analytical uncertainties.

In summary, this paper presented a fresh look at the 40Ar/39Ar method, by recasting every aspect of it into a matrix form and rigorously keeping track of all covariances. Thus, the methods outlined in this paper put the 40Ar/39Ar method on an equal footing with the U-Pb method (McLean et al.2011). Using the same data reduction framework for both methods will improve their intercomparability, which in turn will benefit the accuracy and precision of the geologic time scale (Min et al.2000Kuiper et al.2008).

Acknowledgments

The author would like to thank David Phillips and Erin Matchan (Melbourne) for providing pilot data which have greatly accelerated the development of Ar-Ar_Redux. The manuscript benefited from discussions with James Schwanethal and Noah McLean, and constructive reviews by Chris Hall and three anonymous reviewers. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement n259504 (‘KArSD’).

Appendix A: calculation of interference corrections by mass spectrometric analysis of co-irradiated monitor materials

Neutron reactions on Ca produce interferences on 36Ar and 39Ar, which can be corrected for by monitoring the 36Ar/37Ar- and 39Ar/37Ar-ratios of co-irradiated Ca-salts (Section 8). In this Section, we will use the same simplified regression methods as in Sections 3 and 4. If the three Ar-isotopes of interest are measured in multicollector mode, then their time resolved and blank corrected signals can be cast into the following [n × 2] logratio matrix:

           ⌊  [36Arb(ca,t1)]   [39Arb(ca,t1)] ⌋
           | l[337A6rbb(ca,t1)]  l[3739Arbb(ca,t1)] |
           || l 3A7Arrb((caca,,t2t))   l 37AArrb(c(caa,t,t2)) ||
L (ca,b,t) = ||       .  2         .   2  ||
           |⌈  [    ..    ]   [   ..     ]|⌉
             l 336A7Arrbb((caca,,tnt))   l 3937AArrbb(c(caa,t,tn))
                      n             n
(95)

resulting in a vector of logratio intercepts X(ca) and covariance matrix ΣX(ca). For the detector calibration, we replace Equation 29 with:

  ∗
X  (ca) = [X (ca) Z(d[37]) Z(d[39] Z (d[40])]
(96)

with covariance matrix ΣX(ca):

        ⌊                                           ⌋
        | ΣX (ca)     02,1   2     02,1         02,1    |
Σ∗X(ca) = |⌈ 01,2   σ [Z (d[37])]       0   2       0     |⌉
           01,2        0      σ[Z(d[39])]       0   2
           01,2        0           0      σ[Z(d[40])]
(97)

where, for the sake of notational simplicity, we have assumed that only a single Ca-salt measurement was made (accommodating duplicate analyses is trivial). Note that Equations 96 and 97 use Z(d[40]) instead of Z(d[36]), implying equal sensitivities of detectors d[36] and d[40]. This assumption is valid because the sensitivity difference between said detectors is accounted for by the mass fractionation correction. Equations 31 and 32 remain the same but use the following Jacobian matrix:

        [ 1 0  1   0   − 1 ]
J′C(ca) =   0 1  1  − 1  0
(98)

We thus obtain a two-element vector of sensitivity-corrected logratio intercepts C(ca) and its covariance matrix ΣC(ca). For the mass fractionation correction, we first append the air shot data to the calibration-corrected logratio intercepts:

        [            [    ] ]
 ∗                    40Ar
C (ca) = C (ca) A (37) l 36Ar a
(99)

with [4 × 4] covariance matrix ΣC(ca):

        ⌊                      ⌋
          ΣC (ca)    02,1    02,1
Σ ∗C(ca) = ⌈  01,2   σ[A (37)]2   0 ⌉
            01,2      0       0
(100)

Recasting in matrix form, the fractionation-corrected Ca-salt measurements and their covariances are given by:

F (ca) = C∗(ca) JF(ca)
(101)

and

ΣF (ca) = J′   Σ∗    JF(ca)
         F(ca)  C(ca)
(102)

respectively, where JF(ca) is the Jacobian matrix of the mass fractionation calibration and JF(ca)its transpose:

        [                     ]
J′    =   1  0 − 0.240  − 0.240
 F(ca)     0  1  0.487   0.487
(103)

For ‘peak hopping’ data, Equation 95 can be replaced with three vectors containing the logs of the time-resolved 36Ar, 37Ar and 39Ar signals, which may be processed as in Section 4 to calculate the logratio intercepts. Since detector calibration does not apply to single collector instruments, Equations 96-103 can be safely skipped. Next, we apply the decay correction which, as explained in Section 7, affects both 37Ar and 39Ar. At this point the data reduction of the Ca and K-interference monitors is merged with that of the samples and fluence monitors. This is achieved by collating their respective decay corrections:

r(i) = [r(λi,τ[u1]) r(λi,τ [s1]) r(λi,τ[u2]) r(λi,τ[s2]) r(λi,τ[s3]) r(λi,τ[ca]) r(λi,τ[k])]
(104)

the covariance matrices of which are given by Equation 45 with

      [ ∂r(λ ,τ[u]) ∂r(λ,τ[s ]) ∂r(λ,τ[u ]) ∂r(λ ,τ [s ]) ∂r(λ ,τ[s ]) ∂r(λ ,τ[ca])∂r(λ ,τ[k]&
Jr(i) =  ---i----1------i---1- ----i---2------i---2- ----i---3- ----i----------i-----
           ∂λi        ∂λi         ∂λi        ∂λi        ∂λi        ∂λi        ∂λi
(105)

To apply these decay corrections, we append them to the fractionation-corrected logratios:

 ∗
F  = [F F (ca) F(k) r(37) r(39)]
(106)

with [37 × 37] covariance matrix

     ⌊                                  ⌋
        ΣF    020,2   020,1   020,7   020,7
     || 02,20  ΣF (ca)  02,1    02,7    02,7  ||
Σ∗F = || 01,20   01,2   ΣF (k)   01,7    01,7  ||
     ⌈ 07,20   07,2    07,1   Σr(37)   07,7  ⌉
       07,20   07,2    07,1    07,7   Σr(39)
(107)

These values are then simply plugged into Equations 50 and 51:

     ∗
I = F  JD
(108)

      ′  ∗
ΣI = JD ΣF JD
(109)

where JD is the [37 × 23] Jacobian matrix and JDits transpose:

J ′D = [123,23 J∗D(37) J∗D(39)]
(110)

with

        ⌊                                             ⌋
          J∗D∗(37)   04,1    04,1    04,1    04,1   04,1  04,1
        ||  04,1   J∗D∗(37)   04,1    04,1    04,1   04,1  04,1 ||
        ||  04,1    04,1   J∗∗     04,1    04,1   04,1  04,1 ||
        ||  0      0      D0(37)  J∗∗     0     0    0   ||
J∗D(37) = ||   4,1     4,1     4,1    D(37)   ∗4∗,1    4,1   4,1 ||
        ||  04,1    04,1    04,1    04,1   JD(37)  04,1  04,1 ||
        |⌈   0      0      0      0      0    − 1   0  |⌉
            0      0      0      0      0    − 1   0
            0      0      0      0      0     0    0
(111)

and

        ⌊  ∗∗                                         ⌋
        | JD(39)   0∗4∗,1    04,1    04,1    04,1   04,1  04,1 |
        ||  04,1   JD(39)   0∗4∗,1    04,1    04,1   04,1  04,1 ||
        ||  04,1    04,1   JD(39)   04,1    04,1   04,1  04,1 ||
J∗    = ||  04,1    04,1    04,1   J∗D∗(39)   04,1   04,1  04,1 ||
 D(39)   ||  04,1    04,1    04,1    04,1   J∗D∗(39)  04,1  04,1 ||
        ||   0      0      0      0      0     0    0  ||
        ⌈   0      0      0      0      0     1    0  ⌉
            0      0      0      0      0     0    1
(112)

with JD(37)∗∗ and JD(39)∗∗ as in Equation 54. This completes the Ca-interference correction. The K-interference on 40Ar and (as discussed in Appendix B) 38Ar, can be corrected in a very similar manner by monitoring 40Ar/39Ar and 38Ar/39Ar in K-glass.

Appendix B: Cl-interference correction accounting for all sources of 38Ar

As mentioned at the end of Section 9, the Cl-intereference correction on 36Ar implemented in Equation 5 does not account for the presence of atmospheric 38Ar or the production of 38Ar from K. Doing so is straightforward but requires a reformulation of Equation 5:

     1− a+ b+ c − g− h+ i
R =  ---------------------− f
         d − e− j + k
(113)

with a-f as defined in Equations 6-11 and

g = [38   ]
 --Ar
 36Ara[36   ]
 --Ar
 38Arcl (114)
h = [40Ar ]
 36Ara[36Ar ]
 38Arcl[ 38Ar-]
  39Ark[ 39Ar-]
  40Arm (115)
i = [40   ]
 36Ar
   Ara[36   ]
 38Ar
   Arcl[ 38   ]
  39Ar-
   Ark[ 39   ]
  37Ar-
   Arca[37   ]
 40Ar
   Arm (116)
j = [     ]
 38Ar
 36Ara[     ]
 36Ar
 38Arcl[     ]
  39Ar-
  40Arm (117)
k = [38Ar ]
 36Ara[36Ar ]
 38Arcl[ 39Ar ]
  37Ar-ca[37Ar ]
 40Arm (118)

This formulation requires adjustment of Sections 10 and 11 and the addition of the [38Ar]
 39Ark to Section 8, which is omitted here for brevity.

Appendix C: A brief introduction to Ar-Ar_Redux

In its present form, Ar-Ar_Redux exists as a package in a statistical programming environment called R. After installing R from http://r-project.org, Ar-Ar_Redux can be installed by typing

install.packages(’ArArRedux’)

Once installed, the package can be loaded by typing

library(ArArRedux)

The first step in the data reduction procedure is to load the time resolved mass spectrometer signals and turn them into a vector of logratio intercepts with associated covariance matrix. The read function groups the calculations listed in Sections 3, 4 and 5:

X <- read(xfile="Samples.csv", masses=c("Ar37","Ar38","Ar39","Ar40","Ar36"),  
          blabel="BLANK#", Jpos=c(3,15), kfile="K-glass.csv", cafile="Ca-salt.csv",  
          dfile="Calibration.csv", dlabels=c("H1","AX","L1","L2"))

where xfile is the name of a file containing the time resolved mass spectrometer data of all the samples, fluence monitors and blanks; masses is a vector specifying the order in which the argon isotopes are listed within xfile; blabel is the prefix of the blanks listed in xfile; Jpos is a vector with the positions of the fluence monitors within the irradiation stack; kfile is the name of a file containing the time resolved mass spectrometer signals of co-irradiated K-bearing monitor glass, formatted in the same way as xfile; cafile contains the same information for the co-irradiated Ca-bearing salts; dfile contains the detector intercalibration data and dlabels is a list specifying the order in which the detectors are listed within dfile. Next, we form a list of two fractionation corrections, one for each denominator isotope used in Equation 5 (i.e. 37Ar, 39Ar and 40Ar):

fract <- list(fractionation("AirL2.csv",detector="L2",PH=TRUE),  
              fractionation("AirAX.csv",detector="AX",PH=TRUE),  
              fractionation("AirH1.csv",detector="H1",PH=FALSE))

where the fractionation function performs the calculations outlined in Section 6 and Appendix A; detector specifies the name of the detector of interest; and PH is a boolean flag indicating whether the data are collected in multicollector or ‘peak hopping’ mode. The last file that needs to be loaded contains the neutron irradiation schedule:

irr <- loadirradiations("irradiations.csv")

The process function carries out the fractionation, decay and interference corrections (Sections 6, 7, 8 and 9), interpolates the J-factors and calculates the ages (Sections 11 and 12):

ages <- process(X,irr,fract)

The following three lines are used to tabulate the results, view the covariance structure as a coloured correlation matrix, and calculate the weighted mean age of a subset (in this example samples S1-5) of the data, respectively:

summary(ages)  
corrplot(ages)  
weightedmean(ages,c("S1","S2","S3","S4","S5"))

Ar-Ar_Redux is very flexible. For example, all but the first four arguments to the read function are optional. If, for instance, no co-irradiated K-glass or Ca-salt were analysed, then it is possible to specify the interference corrections explicitly. A comprehensive overview of all the options falls outside the scope if this short Appendix. A more extensive tutorial is provided on http://redux.london-geochron.com. Contextual help within the R environment can be obtained from Ar-Ar_Redux’s built-in documentation. For example, to learn more about the read function, it suffices to type ?read at the command prompt.

References

   Aitchison, J., 1986. The statistical analysis of compositional data. London, Chapman and Hall.

   Berger, G.W., York, D., 1970. Precision of the 40Ar/39Ar dating technique. Earth and Planetary Science Letters 9, 39–44.

   Bowring, J., McLean, N., Bowring, S., 2011. Engineering cyber infrastructure for U-Pb geochronology: Tripoli and U-Pb_Redux. Geochemistry, Geophysics, Geosystems 12.

   Foland, K., Fleming, T., Heimann, A., Elliot, D., 1993. Potassium-argon dating of fine-grained basalts with massive Ar loss: Application of the 40Ar/39Ar technique to plagioclase and glass from the Kirkpatrick Basalt, Antarctica. Chemical Geology 107, 173–190.

   Koppers, A.A., 2002. ArArCALC – software for 40Ar/39Ar age calculations. Computers & Geosciences 28, 605–619.

   Kuiper, K.F., Deino, A., Hilgen, F.J., Krijgsman, W., Renne, P.R., Wijbrans, J.R., 2008. Synchronizing Rock Clocks of Earth History. Science 320, 500–504. doi:10.1126/science.1154339.

   Lee, J.Y., Marti, K., Severinghaus, J.P., Kawamura, K., Yoo, H.S., Lee, J.B., Kim, J.S., 2006. A redetermination of the isotopic abundances of atmospheric Ar. Geochimica et Cosmochimica Acta 70, 4507–4512.

   McDougall, I., Harrison, T.M., 1999. Geochronology and Thermochronology by the 40Ar/39Ar method. Oxford University Press, New York.

   McLean, N., Bowring, J., Bowring, S., 2011. An algorithm for U-Pb isotope dilution data reduction and uncertainty propagation. Geochemistry, Geophysics, Geosystems 12.

   Min, K., Mundil, R., Renne, P.R., Ludwig, K.R., 2000. A test for systematic errors in 40Ar/39Ar geochronology through comparison with U/Pb analysis of a 1.1-Ga rhyolite. Geochimica et Cosmochimica Acta 64, 73–98.

   Nelder, J.A., Wedderburn, R.W.M., 1972. Generalized linear models. Journal of the Royal Statistical Society. Series A (General) 135, pp. 370–384.

   Renne, P.R., Norman, E.B., 2001. Determination of the half-life of 37Ar by mass spectrometry. Physical Review C 63, 047302.

   Renne, P.R., Sharp, Z.D., Heizler, M.T., 2008. Cl-derived argon isotope production in the CLICIT facility of OSTR reactor and the effects of the Cl-correction in 40Ar/39Ar geochronology. Chemical Geology 255, 463–466.

   Renne, P.R., Swisher, C.C., Deino, A.L., Karner, D.B., Owens, T.L., DePaolo, D.J., 1998. Intercalibration of standards, absolute ages and uncertainties in 40Ar/39Ar dating. Chemical Geology 145, 117–152.

   Scaillet, S., 2000. Numerical error analysis in 40Ar/39Ar dating. Chemical Geology 162, 269–298.

   Stoenner, R., Schaeffer, O., Katcoff, S., 1965. Half-lives of argon-37, argon-39, and argon-42. Science 148, 1325–1328.

   Vermeesch, P., 2010. HelioPlot, and the treatment of overdispersed (U-Th-Sm)/He data. Chemical Geology 271, 108 – 111. doi:DOI: 10.1016/j.chemgeo.2010.01.002.

   Wijbrans, J.R., McDougall, I., 1986. 40Ar/39Ar dating of white micas from an Alpine high-pressure metamorphic belt on Naxos (Greece): the resetting of the argon isotopic system. Contributions to Mineralogy and Petrology 93, 187–194.

   Wood, S.N., 2015. Core Statistics. 1 ed., Cambridge University Press.

   York, D., 1969. Least squares fitting of a straight line with correlated errors. Earth and Planetary Science Letters 5, 320–324.

   Young, E.D., Galy, A., Nagahara, H., 2002. Kinetic and equilibrium mass-dependent isotope fractionation laws in nature and their geochemical and cosmochemical significance. Geochimica et Cosmochimica Acta 66, 1095–1104.