P_mu_total.m

Syntax:

out = P_mu_total(z,h,consts,dflag)

This function calculates the nuclide production rate due to muons at a particular surface elevation and depth below the surface. The method is described in:

Heisinger, B., Lal, D., Jull, A.J.T., Kubik. P., Ivy-Ochs, S., Neumaier, S., Knie, K., Lazarev, V., and Nolte, E., 2002. Production of selected cosmogenic radionuclides by muons: 1. Fast muons. Earth and Planetary Science Letters 200, pp. 345-355. (henceforth, H2002a).

Heisinger, B., Lal, D., Jull, A.J.T., Kubik. P., Ivy-Ochs, S., Knie, K., and Nolte, E., 2002. Production of selected cosmogenic radionuclides by muons: 2. Capture of negative muons. Earth and Planetary Science Letters 200, pp. 357-369. (henceforth, H2002b).

The input argument z is depth below the surface in g $\cdot$ cm$^{-2}$. This argument can be a vector.

The input argument h is the atmospheric pressure at the surface in hPa.

The input argument consts is a structure containing nuclide-specific constants. The fields are as follows:

consts.Natoms	number density of target atoms in quartz	atoms $\cdot$ g$^{-1}$
consts.k_neg	summary yield for production by negative muon capture in quartz	atoms $\cdot$ (stopped $\mu_{-}$)$^{-1}$
consts.sigma190	measured cross-section at 190 GeV for production by fast muon reactions	cm$^{-2}$

The argument dflag is a string variable telling the function what to return. If dflag = `no,' the output is simply a vector of production rates of the same size as the input argument z. If dflag = `yes,' the output is a structure containing diagnostic information, as follows:

out.phi_vert_slhl	Flux of vertically traveling muons at the specified depths at sea level	muons $\cdot$ cm $^{-2} \cdot$ sr $^{-1} \cdot$ s$^{-1}$
out.R_vert_slhl	Stopping rate of vertically traveling muons at the specified depths at sea level	muons $\cdot$ g $^{-1} \cdot$ sr $^{-1} \cdot$ s$^{-1}$
out.phi_vert_site	Flux of vertically traveling muons at the specified depths at the site elevation	muons $\cdot$ cm $^{-2} \cdot$ sr $^{-1} \cdot$ s$^{-1}$
out.R_vert_site	Stopping rate of vertically traveling muons at the specified depths at the site elevation	muons $\cdot$ g $^{-1} \cdot$ sr $^{-1} \cdot$ s$^{-1}$
out.phi	Total flux of muons at the specified depths at the site elevation	muons $\cdot$ cm $^{-2} \cdot$ yr$^{-1}$
out.R	Total stopping rate of negative muons at the specified depths at the site elevation	negative muons $\cdot$ g $^{-1} \cdot$ yr$^{-1}$
out.Beta	Factor describing the energy dependence of fast muon reaction cross-sections	nondimensional
out.Ebar	Mean muon energy at the specified depths	GeV
out.P_fast	Nuclide production rate by fast muon reactions at the specified depths	atoms $\cdot$ g $^{-1} \cdot$ yr$^{-1}$
out.P_neg	Nuclide production rate by negative muon capture at the specified depths	atoms $\cdot$ g $^{-1} \cdot$ yr$^{-1}$
out.H	Atmospheric depth between the site elevation and sea level	g $\cdot$ cm$^{-2}$
out.LZ	Atmospheric attenuation lengths for vertically traveling muons stopping at the specified depths	g $\cdot$ cm$^{-2}$

The calculation goes as follows.

1. Calculate the flux of vertically traveling muons as a function of depth at sea level and high latitude.

This is accomplished by Equations (1) and (2) in H2002a. The flux of vertically traveling muons at a depth $z$ at sea level and high latitude $\Phi_{v,0}(z)$is:

$$\Phi_{v,0}(z) = \frac {5.401 \times 10^{7}} {\left( z + 21000 \ri... ...+1000 \right) ^{1.66} + 1.567 \times 10^{5} \right] } e^{-5.5 \times 10^{-6} z}$$ (42)

for depths $z < 200,000$ g $\cdot$ cm$^{-2}$. This is Equation (1) from H2002a, modified so that $z$ is in g $\cdot$ cm$^{-2}$. For greater depths, $\Phi_{v}(z)$ is given by Equation (2) of H2002a, similarly modified:

$$\Phi_{v,0}(z) = 1.82 \times 10^{-6} \left[ \frac{121100}{z} \right] ^{2} e^{- \frac{z}{121100}} + 2.84 \times 10^{-13}$$ (43)

The units of $\Phi_{v}(z)$ are muons $\cdot$ cm $^{-2}\cdot$ s $^{-1}\cdot$ sr$^{-1}$.

2. Calculate the stopping rate of vertically traveling muons as a function of depth at sea level and high latitude, which is equivalent to the range spectrum of vertically traveling muons at the surface.

The stopping rate of vertically traveling muons at a depth $z$ at sea level and high latitude $R_{v,0}(z)$ is the derivative of the flux of vertically traveling muons with respect to depth. It has units of muons $\cdot$ g $^{-1}\cdot$ s$^{-1}$. For depths $z < 200,000$ g $\cdot$ cm$^{-2}$,

$$R_{v,0}(z) = \frac{d}{dz} \left( \Phi_{v,0}(z) \right) = -5.401\t... ...}{dz} - a \left( \frac{db}{dz}c + \frac{dc}{dz}b \right) } {b^{2}c^{2}} \right]$$ (44)

where

$$a = e^{-5.5 \times 10^{-6} z} \frac{da}{dz} = -5.5 \times 10^{-6}e^{-5.5 \times 10^{-6} z}$$ (45)

$$b = \left( z + 21000 \right) \frac{db}{dz} = 1$$ (46)

$$c = \left( z + 1000 \right) ^{1.66} + 1.567 \times 10^{5} \frac{dc}{dz} = 1.66 \left( z + 1000 \right) ^{0.66}$$ (47)

The negative sign is added because the muon flux decreases with increasing depth, so its derivative ought properly to be negative, but we would like a positive value for the stopping rate. For greater depths,

$$R_{v,0}(z) = \frac{d}{dz} \left( \Phi_{v}(z) \right) = -1.82 \times 10^{-6} \left[ \frac{df}{dz} g + \frac{dg}{dz} f \right]$$ (48)

where

$$f = \left[ \frac{121100}{z} \right]^{2} \frac{df}{dz} = \frac {-2 \left( 121100 \right)^{2} } {z^{3}}$$ (49)

$$g = e^{-\frac{z}{121100}} \frac{dg}{dz} = - \frac {e^{-\frac{z}{121100}}}{121100}$$ (50)

The stopping rate of vertically traveling muons as a function of depth is equivalent to the range spectrum of vertically traveling muons at the surface. That is, a muon that had a range of 1000 g $\cdot$ cm$^{-2}$ at the surface will stop at a depth of 1000 g $\cdot$ cm$^{-2}$.

3. Adjust the range spectrum of vertically traveling muons to a different elevation.

First, calculate the difference in atmospheric pressure between sea level and the elevation of interest. We use the standard atmosphere approximation to convert elevation to atmospheric pressure. The atmospheric pressure in hPa as a function of elevation is:

$$P_{atm}(z) = P_{atm,0} exp \left( -\frac{0.03417}{0.0065} \left[ \ln{288.15} -\ln{\left( 288.15 - 0.0065 h \right)} \right] \right)$$ (51)

where $h$ is the elevation in meters and $P_{atm,0} = 1013.25$ hPa is the sea level pressure. Pressure can then be converted to the quantity of interest, that is, $\delta z$, the atmospheric depth in g $\cdot$ cm$^{-2}$ between the site of interest and sea level, by $\delta z=1.019716 \left( P_{atm,0} - P_{atm}(z) \right)$.

If the atmospheric depth between sea level and the site of interest is $\delta z$, then the vertically traveling muon flux at the surface as a function of muon range $Z$ at sea level $R_{v,0}(Z)$ can be scaled to the vertically traveling muon flux at the surface as a function of muon range at the site $R_{v}(Z)$ by:

$$R_{v}(Z) = R_{v,0}(Z)e^{\frac{\delta z}{\Lambda_{\mu}(Z)}}$$ (52)

where $\Lambda_{\mu}(Z)$ is a range-dependent, that is, energy-dependent effective atmospheric attenuation length. These attenuation lengths are measured in studies of atmospheric muon fluxes. We follow H2002b and use a relation between muon momentum $P$ and effective attenuation length $\Lambda_{mu}(P)$ derived from:

Boezio, M, and 33 co-authors, 2000. Measurement of the flux of atmospheric muons with the CAPRICE94 apparatus. Physical Review D, 62, 032007.

The relation is:

$$\Lambda_{mu}(P) = 263 + 150P$$ (53)

In order to use this relation, we convert muon range to momentum using tabulated values in:

Groom, D.E., Mokhov, N.V., and Striganov, S.I., 2001. Muon stopping power and range tables 10 MeV - 100 TeV. Atomic data and nuclear data tables 78, pp. 183-356.

Having obtained $\Lambda_{mu}$ for muons stopping at the depths of interest, we can then obtain the range spectrum of vertically traveling muons at the surface at the site of interest, $R_{v}(Z)$, which is equivalent to the muon stopping rate as a function of depth at the site of interest $R_{v}(z)$, by applying Equation 53.

4. Latitudinal variability in the range spectrum.

Although the muon range spectrum is expected to change with latitude as well as elevation due to geomagnetic effects, this effect is expected to be small. We follow H2002b and ignore it.

5. Calculate the flux of vertically traveling muons at the site of interest as a function of depth.

The flux of vertically traveling muons as a function of depth at the site of interest $\phi_{v}(z)$ is the integral of the muon stopping rate as a function of depth at the site of interest $R_{v}(z)$, which we have just calculated, from infinite depth to depth $z$. That is, the muon flux at a particular depth is composed of all the muons which stop below that depth. Thus, the flux of vertically traveling muons as a function of depth at the site of interest $\phi_{v}(z)$ is:

$$\phi_{v}(z) = \int_{z}^{\infty} R_{v}(x) dx$$ (54)

We actually do this integral numerically, although it would be possible to obtain an analytical expression. We take the limit of integration to be $2 \times 10^{5}$ g $\cdot$ cm$^{-2}$, where the muon flux is essentially negligible for our purposes.

6. Calculate the total muon flux as a function of depth at the site of interest.

Following equations (3) - (5) in H2002a, the zenith angle dependence of the muon flux is:

$$\phi(z,\theta) = \phi_{v}(z)cos^{n(z)}\theta$$ (55)

where $\theta$ is the zenith angle and $n$ is given by Equation (4) in H2002a, modified so that $z$ is given in g $\cdot$ cm$^{-2}$:

$$n(z) = 3.21 - 0.297 \ln{ \left( \frac{z}{100} + 42 \right) } + 1.21 \times 10^{-5} z$$ (56)

The total muon flux at a particular depth $\phi(z)$ then consists of Equation 56 integrated over the entire upper hemisphere and has units of muons $\cdot$ cm$^-2 \cdot$ s$^{-1}$. This is given by:

$$\phi(z) = \frac{2\pi}{n(z+\delta z)+1} \phi_{v}(z)$$ (57)

7. Calculate the total muon stopping rate as a function of depth at the site of interest.

The total muon stopping rate as a function of depth $R{z}$ is the derivative with respect to depth of the total muon flux as a function of depth $\phi(z)$, that is:

$$R(z) = \frac{d}{dz} \left(\phi(z) \right)$$ (58)

$$= \frac{d}{dz} \left( \frac{2\pi}{n(z+\delta z)+1} \phi_{v}(z) \right)$$ (59)

$$= \frac{2\pi}{n(z+\delta z)+1} \frac{d}{dz}\left( \phi_{v}(z) \right) - \phi_{v}(z) \frac{d}{dz}\left( \frac{2\pi}{n(z+\delta z)+1} \right)$$ (60)

$$= \frac{2\pi}{n(z+\delta z)+1} R_{v}(z) - \phi_{v}(z) \left(-2\pi... ... \left( n(z+\delta z)+1 \right)^{2} \frac{d}{dz}\left( n(z+\delta z) + 1\right)$$ (61)

$$= \frac{2\pi}{n(z+\delta z)+1} R_{v}(z) - \phi_{v}(z) \left(-2\pi... ...-0.297 \times 10^{-2}}{\frac{z+\delta z}{100}+42} + 1.21 \times 10^{-5} \right]$$ (62)

which, as we have already calculated $n(z)$, $R_{v}(z)$, and $\phi_{v}(z)$, we can calculate easily. Again, a factor of -1 is added to obtain a positive number of stopped muons.

To summarize, we have now calculated the total muon flux $\phi(z)$ (muons $\cdot$ cm $^{-2} \cdot$ s $^{-1}$) and the total stopping rate of muons $R(z)$ (muons $\cdot$ g $^{-1} \cdot$ s $^{-1}$) at our site. We convert these to muons $\cdot$ cm $^{-2} \cdot$ yr $^{-1}$ and muons $\cdot$ g $^{-1} \cdot$ yr $^{-1}$, respectively, by multiplying them by $3.154 \times 10^{7}$ s $\cdot$ yr$^{-1}$.

Finally, we compute the stopping rate of negative muons $R^{-}(z)=0.44R(z)$ (negative muons $\cdot$ g $^{-1} \cdot$ yr $^{-1}$).

8. Calculate the nuclide production rate due to negative muon capture.

Following Equation (11) in H2002b, the production rate of nuclide $i$ (atoms $\cdot$ g $^{-1} \cdot$ yr$^{-1}$) from negative muon capture $P_{i,\mu -}(z)$ is:

$$P_{i,\mu -}(z) = R^{-}(z) f_{i,C}f_{i,D}f^{*}_{i}$$ (63)

where

$$f_{10,C} = 0.704 f_{26,C} =0.296$$ (64)

$$f_{10,D} = 0.1828 f_{26,D} =0.6559$$ (65)

$$f^{*}_{10} = 0.0043 f^{*}_{26} =0.022$$ (66)

9. Calculate the nuclide production rate due to fast muon reactions.

Following Equation (14) in H2002a, the production rate of nuclide $i$ (atoms $\cdot$ g $^{-1} \cdot$ yr$^{-1}$) from fast muon interactions $P_{i,\mu fast}(z)$ is:

$$P_{i,\mu fast}(z) = \phi(z) \sigma_{0,i} \beta(z) \left(\bar{E}(z) \right)^{\alpha} N_{t,i}$$ (67)

where $\alpha=0.75$.

$\sigma_{0,i}$ is the nominal zero-energy muon interaction cross-section for the reaction responsible for producing nuclide $i$. Here is has units of cm$^{-2}$. The muon interaction cross-section for a particular reaction is thought to depend on the muon energy as follows:

$$\sigma_{i}(E) = \sigma_{0,i}E^{\alpha}$$ (68)

$\sigma_{0,i}$ is determined from the measured cross-sections $\sigma_{i}(E)$ at 190 GeV energy in Table 1 of H2002a using this equation. These values for $^{10}$Be and $^{26}$Al are:

$$\sigma_{0,10} = \frac{0.094 \times 10^{-27}}{190^{\alpha}} \sigma_{0,26} = \frac{1.41 \times 10^{-27}}{190^{\alpha}}$$ (69)

$N_{t,i}$ is the number density of atoms of the target element (atoms $\cdot$ g$^{-1}$). The values for O and Si relevant to $^{10}$Be and $^{26}$Al production respectively are:

$$N_{t,10} = 2.006 \times 10^{22} N_{t,26} = 1.003 \times 10^{22}$$ (70)

where $\beta$ is a function of depth and is approximated by Equation (16) of H2002a:

$$\beta(z) = 0.846 = 0.015 \ln{ \left(\frac{z}{100} +1\right) } + 0.003139 \left[ \ln{ \left( \frac{z}{100} \right) } \right]$$ (71)

and $\bar{E}(z)$ is the mean muon energy at depth $z$ and is given by Equation (11) of H2002a:

$$\bar{E}(z) = 7.6 + 321.7 \left( 1 - e^{-8.059 \times 10^{-6} z} \right) + 50.7 (1 - e^{-5.05 \times 10^{-7} z})$$ (72)