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 cm. 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 g |

consts.k_neg | summary yield for production by negative muon capture in quartz | atoms (stopped ) |

consts.sigma190 | measured cross-section at 190 GeV for production by fast muon reactions | cm |

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 cm sr s |

out.R_vert_slhl | Stopping rate of vertically traveling muons at the specified depths at sea level | muons g sr s |

out.phi_vert_site | Flux of vertically traveling muons at the specified depths at the site elevation | muons cm sr s |

out.R_vert_site | Stopping rate of vertically traveling muons at the specified depths at the site elevation | muons g sr s |

out.phi | Total flux of muons at the specified depths at the site elevation | muons cm yr |

out.R | Total stopping rate of negative muons at the specified depths at the site elevation | negative muons g yr |

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 g yr |

out.P_neg | Nuclide production rate by negative muon capture at the specified depths | atoms g yr |

out.H | Atmospheric depth between the site elevation and sea level | g cm |

out.LZ | Atmospheric attenuation lengths for vertically traveling muons stopping at the specified depths | g cm |

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 at sea level and high latitude is:

(42) |

for depths g cm. This is Equation (1) from H2002a, modified so that is in g cm. For greater depths, is given by Equation (2) of H2002a, similarly modified:

(43) |

The units of are muons cm s sr.

**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 at sea level and high latitude is the derivative of the flux of vertically traveling muons with respect to depth. It has units of muons g s. For depths g cm,

(44) |

where

(45) | ||||

(46) | ||||

(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,

(48) |

where

(49) | ||||

(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 cm at the surface will stop at a depth of 1000 g cm.

**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:

(51) |

where is the elevation in meters and hPa is the sea level pressure. Pressure can then be converted to the quantity of interest, that is, , the atmospheric depth in g cm between the site of interest and sea level, by .

If the atmospheric depth between sea level and the site of interest is , then the vertically traveling muon flux at the surface as a function of muon range at sea level can be scaled to the vertically traveling muon flux at the surface as a function of muon range at the site by:

where 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 and effective attenuation length 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:

(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 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, , which is equivalent to the muon stopping rate as a function of depth at the site of interest , 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 is the integral of the muon stopping rate as a function of depth at the site of interest , which we have just calculated, from infinite depth to depth . 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 is:

(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 g cm, 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:

where is the zenith angle and is given by Equation (4) in H2002a, modified so that is given in g cm:

(56) |

The total muon flux at a particular depth then consists of Equation 56 integrated over the entire upper hemisphere and has units of muons cm s. This is given by:

(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 is the derivative with respect to depth of the total muon flux as a function of depth , that is:

(58) | ||

(59) | ||

(60) | ||

(61) | ||

(62) |

which, as we have already calculated , , and , 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 (muons cm s ) and the total stopping rate of muons (muons g s ) at our site. We convert these to muons cm yr and muons g yr , respectively, by multiplying them by s yr.

Finally, we compute the stopping rate of negative muons (negative muons g yr ).

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

Following Equation (11) in H2002b, the production rate of nuclide (atoms g yr) from negative muon capture is:

(63) |

where

(64) | ||||

(65) | ||||

(66) |

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

Following Equation (14) in H2002a, the production rate of nuclide (atoms g yr) from fast muon interactions is:

where .

is the nominal zero-energy muon interaction cross-section for the reaction responsible for producing nuclide . Here is has units of cm. The muon interaction cross-section for a particular reaction is thought to depend on the muon energy as follows:

(68) |

is determined from the measured cross-sections at 190 GeV energy in Table 1 of H2002a using this equation. These values for Be and Al are:

(69) |

is the number density of atoms of the target element (atoms g). The values for O and Si relevant to Be and Al production respectively are:

(70) |

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

(71) |

and is the mean muon energy at depth and is given by Equation (11) of H2002a:

(72) |