Calibrating 36Cl production rates

Calibrating ³⁶Cl production rates

ACE determines HLSL production rates by calibrating independently dated datasets. For ³⁶Cl, this procedure is complicated by the fact that this nuclide is produced by three independent processes:

Spallation of ⁴⁰Ca
Spallation of ³⁹K
Neutron capture of ³⁵Cl

As a result, to calibrate for ³⁶Cl a multiple regression must be performed to determine the relative importance of these processes. A further complication with ³⁶Cl is that state of the art models of neutron capture of ³⁵Cl rely on a formalism of epithermal (energies between 0.1 MeV and 0.5 eV) and thermal (around 0.025 eV) neutron fluxes, where epithermal neutrons are generated by the moderation of fast neutrons, and thermal neutrons are generated by the moderation of epithermal neutrons (eg Phillips et al., 2001):

The total production of ³⁶Cl from neutron capture is the sum of these equations. To calibrate via a least squares regression the above equations must be expressed in terms of ${P}_{f0}$ , the production rate of epithermal neutrons from fast neutrons at the air/ground interface, which is the source of epithermal neutrons and therefore also the original source of thermal neutrons. Here we explain how we find Pf0 using a multiple regression in ACE.

Calculations using the diffusion equation solution for neutron capture (detailed below) show that the total production rate of ³⁶Cl from neutron capture is a linear function of ${P}_{f0}$ , plus a constant term which is directly proportional to muon production rates:

$Pn = A \space {P}_{f0} + B \space P_{\mu n o}$ (1)

The dependence of neutron capture on muons results from the fact that fast muons can generate photodisintegration reactions which emit neutrons, which can then be captured by ³⁵Cl (Gosse and Phillips, 2001, Section 3.3.4). We verify this using Maxima, a free algebraic manipulation software program very much like Mathematica. Maxima can be downloaded here:

http://maxima.sourceforge.net/

(Macintosh users can install maxima automatically using fink.) When installed the linear dependence of neutron capture production rates on ${P}_{f0}$ can be verified using the batch file Diffusion Equation Linearity File. Maxima programs are very simple, so that the learning curve is very shallow. The first lines of the code add some libraries for later use, line 3 is the following:

Pn = S_th*ST*Q_th*f_th/LAMBDA_th_ss*(phistar_th_ss+(1+R_mu_p)
     *FscrDELTAphistar_eth+(1+R_mu_p*R_th)*FscrDELTAphistar_th_ss
     +R_mu_p*phistar_th_ss)+S_th*ST*Q_eth*(1-P_eth_ss)
     *f_eth/LAMBDA_eth_ss*(phistar_eth_ss+(1+R_mu*R_eth)
     *FDELTAphistar_eth+R_mu*phistar_eth_ss);

Which is the sum of equations 3.54 and 3.55 from Gosse and Phillips (2001) given above with some additional terms for thermal neutron scaling (S_th), sample shielding (ST), and sample depth rescaling (Q_th and Q_eth). All terms match Gosse and Phillips notation as closely as possible, eg in the above ${Q}_{th}$ = Q_th.

The maxima batch file then expands the original expression for production rates from neutron capture into an expression which explicity contains ${P}_{f0}$ and $P_{\mu n o}$ . For example line 4 of the maxima file takes the original expression for production from neutron capture (using the %, operator at the beginning of the line) and substitutes Q_th for more explicit terms (Equation 3.80, Gosse and Phillips, 2001):

Q_th = (phistar_th_ss*LAMBDA_f*(1-exp(-Zs/LAMBDA_f))+(1+R_mu_p)
       *FscrDELTAphistar_eth*L_eth_ss*(1-exp(-Zs/L_eth_ss))+
       (1+R_mu_p*R_th)*FscrDELTAphistar_th_ss*L_th_ss
       *(1-exp(-Zs/L_th_ss))+R_mu_p*phistar_th_ss*LAMBDA_mu
       *(1-exp(-Zs/LAMBDA_mu)))/(Zs*(phistar_th_ss+(1+R_mu_p)
       *FscrDELTAphistar_eth+(1+R_mu_p*R_th)*FscrDELTAphistar_th_ss
       +R_mu_p*phistar_th_ss));

This process of subsitution continues to line 22, which contains the expression for $R_{\mu}$ (Equation 3.52, Gosse and Phillips, 2001):

R_mu = S_mu * P_mu_n_0 / (S_th * Pf_0 * R_eth);

This expression contains ${P}_{f0}$ in the denominator, so that the dependence of neutron capture production rates is an inverse function of this term as well as direct functions of phistar_eth_a and phistar_eth_ss (Equation 3.26, Gosse and Phillips, 2001).

At this stage the maxima file has expressed the neutron capture production rate explicitly in terms of ${P}_{f0}$ . To see this expression click here. Line 25 of this file expands the entire expression into products of ${P}_{f0}$ :

isolate(expand(%),Pf_0);

Which makes for a very long expression. However, as a test we can insert numerical values into this expression calculated by other programs to test how $Pn$ varies with ${P}_{f0}$ . As for expressions, we can subsitute numerical values into Maxima. Line 29 of the Maxima script has the following:

%, f_th = 0.007997574;

which literally means “find all occurences of f_th in the previous expression and insert the value 0.007997574”. For actual numbers we use output from the 36Cl calibration sample named “9353” from the Phillips et al (1997) calibration database (detailed here), but in fact the numbers are irrelevant. As long as they are non-zero, they will decide the coefficients between $Pn$ and ${P}_{f0}$ , terms that we are not actually interested in at this stage.

The maxima file then subsitutes numerical values for all terms in the neutron capture production rate equation except for ${P}_{f0}$ and $P_{\mu n o}$ . At the end of the maxima procedure, the output is:

Pn = .01046760791954185 Pf_0 + .002679086726981858 P_mu_n_0

Note the correspondence of this expression to that of Equation (1) above, with actual values for $A$ and $B$ . This means that despite the complexity of the diffusion equation solution for neutron capture, the production rate of ³⁶Cl is a linear function of both ${P}_{f0}$ and $P_{\mu n o}$ (a result which could also be verified using the maxima coeff command). This allows us to use a multiple linear regression when determining the relative contribution of the three production mechanisms for ³⁶Cl. Furthermore, we have the coefficients for the linear regression in the maxima file, we just need to combine them nicely in coefficients of ${P}_{f0}$ and $P_{\mu n o}$ .

To find the coefficients needed for the multiple linear regression, we use another, very similar Maxima file that instead of subsituting numerical values expands the equations and collects like coefficients of ${P}_{f0}$ and $P_{\mu n o}$ . The file Diffusion Equation Refactor Terms does just this, although to make things simple you will need to restart maxima to use it properly. After substituting in all components of ${P}_{f0}$ and $P_{\mu n o}$ , it finds coefficients of ${P}_{f0}$ using the following maxima command (line 28):

coeff(%,Pf_0,1); /*Linear Term*/

Means “in the previous expression, find the linear coefficient of Pf_0“. For the sample 9353 above, this would be .01046760791954185. However we wish to express it in terms of variables so that all samples can be determined. The output from this maxima expression without numerical subsitution is:

(%o93) 0 = %t91 + %t90 + %t89 + %t88 + %t87 + %t86 + %t85 + %t84 + %t83 +
           %t82 + %t81 + %t80 + %t79 + %t78 + %t77 + %t76 + %t75 + %t74 +
           %t73 + %t72 + %t71 + %t70 + %t69 + %t68 + %t67 + %t66 + %t65 +
           %t64 + %t63 + %t62 + %t61 + %t60 + %t59 + %t58 + %t57 + %t56 +
           %t55 + %t54 + %t53 + %t52 + %t51 + %t50 + %t49 + %t48 + %t47 +
           %t46 + %t45 + %t44 + %t43 + %t42 + %t41 + %t40 + %t39 + %t38 +
           %t37 + %t36 + %t35 + %t34 + %t33 + %t32 + %t31 + %t30 + %t29 +
           %t28 + %t27 + %t26 + %t25 + %t24 + %t23 + %t22

Meaning that the coefficient of ${P}_{f0}$ is the sum of terms 22 through to 91 in the maxima output. This is a lof of expressions, but represents the neutron capture production rate in its most expanded and least compact form. Many of the terms are similar to each other, and by trial and error they can be added together to make fairly compact expressions. For example line 30 of the maxima file groups together the first 10 expressions, more than 10 makes for a very long output:

f90(factor(%t22 + %t23 + %t24 + %t25 + %t26 + %t27 + %t28 + %t29 + %t30 + %t31));

The terms are added together, factored as much as possible, and output in fortran 90 form so they can easily be used in ACE. The output is:

-f_eth*LAMBDA_f*(P_eth_ss-1)*(Da*LAMBDA_f**2*L_eth_a*L_eth_ss**2* &
     SIGMA_eth_ss+Da*LAMBDA_f**3*L_eth_ss**2* &
     SIGMA_eth_ss-D_eth_ss*LAMBDA_f**2*L_eth_a*L_eth_ss**2*R_eth* &
     SIGMA_eth_a-Da*LAMBDA_f**3*L_eth_ss**2*R_eth*SIGMA_eth_a+Da* &
     LAMBDA_f**4*L_eth_ss*R_eth*SIGMA_eth_a+D_eth_ss*LAMBDA_f**4* &
     L_eth_a*R_eth*SIGMA_eth_a+D_eth_ss*Da*L_eth_a*L_eth_ss**2* &
     R_eth+Da**2*LAMBDA_f*L_eth_ss**2*R_eth-Da**2*LAMBDA_f**2* &
     L_eth_ss*R_eth-D_eth_ss*Da*LAMBDA_f**2*L_eth_a* &
     R_eth-D_eth_ss*Da*L_eth_a*L_eth_ss**2-D_eth_ss*Da*LAMBDA_f* &
     L_eth_ss**2)*S_th*ST/ &
     (LAMBDA_eth_ss*(Da*L_eth_ss+D_eth_ss*L_eth_a)*(LAMBDA_f**2* &
     SIGMA_eth_a-Da)*(LAMBDA_f**2*SIGMA_eth_ss-D_eth_ss)*Zs)

Unfortunately maxima has no Python output to use directly into ACE, but the codes are very similar. In fact the only changes are “\” instead of “&” for newline characters and “self.L_eth_ss” instead of “L_eth_ss” for compatibility with the rest of ACE. In the maxima file the term above is labelled term %t95 (for term 95). The equivalent code is then copied and pasted into the ACE Component SumCalibrationComponents.py as self.term95:

        self.term95 = \
          -self.f_eth*self.LAMBDA_f*(self.P_eth_ss-1)*(self.Da*self.LAMBDA_f**2*self.L_eth_a*self.L_eth_ss**2* \
          self.SIGMA_eth_ss+self.Da*self.LAMBDA_f**3*self.L_eth_ss**2* \
          self.SIGMA_eth_ss-self.D_eth_ss*self.LAMBDA_f**2*self.L_eth_a*self.L_eth_ss**2*self.R_eth* \
          self.SIGMA_eth_a-self.Da*self.LAMBDA_f**3*self.L_eth_ss**2*self.R_eth*self.SIGMA_eth_a+self.Da* \
          self.LAMBDA_f**4*self.L_eth_ss*self.R_eth*self.SIGMA_eth_a+self.D_eth_ss*self.LAMBDA_f**4* \
          self.L_eth_a*self.R_eth*self.SIGMA_eth_a+self.D_eth_ss*self.Da*self.L_eth_a*self.L_eth_ss**2* \
          self.R_eth+self.Da**2*self.LAMBDA_f*self.L_eth_ss**2*self.R_eth-self.Da**2*self.LAMBDA_f**2* \
          self.L_eth_ss*self.R_eth-self.D_eth_ss*self.Da*self.LAMBDA_f**2*self.L_eth_a* \
          self.R_eth-self.D_eth_ss*self.Da*self.L_eth_a*self.L_eth_ss**2-self.D_eth_ss*self.Da*self.LAMBDA_f* \
          self.L_eth_ss**2)/ \
          (self.LAMBDA_eth_ss*(self.Da*self.L_eth_ss+self.D_eth_ss*self.L_eth_a)*(self.LAMBDA_f**2* \
          self.SIGMA_eth_a-self.Da)*(self.LAMBDA_f**2*self.SIGMA_eth_ss-self.D_eth_ss))

For calibration purposes these terms are copied into ACE and added up to form the coefficient for ${P}_{f0}$ . When coefficients for all samples have been determined the multiple linear regression can be performed. The result is the calibrated production rates for ³⁶Cl which can be obtained for different scalings and calibration datasets using the Experiment Browser:

Click on image to expand

This entry was posted on Thursday, June 11th, 2009 at 1:57 pm and is filed under Software News. You can follow any responses to this entry through the RSS 2.0 feed. You can skip to the end and leave a response. Pinging is currently not allowed.

You must be logged in to post a comment.

ACE : Age Calculation Engine

Calibrating ³⁶Cl production rates

Leave a Reply