Cosmic Ray Muons in Laboratories Deep Underground

We provide comprehensive calculations of total muon fluxes, energy and angular spectra, and mean muon energies in deep underground laboratories - under flat overburdens and mountains and underwater - using our latest calculation code, MUTE v3. For precise modeling, we compiled rock densities and chemical compositions for various underground labs, as well as topographic map profiles of overburdens, and integrated them into our calculations. Our results show excellent agreement with available data for most underground sites when using the latest surface muon flux model, daemonflux. Moreover, since our calculations do not rely on underground measurements of muons or other secondaries, we can verify the consistency of measurements across different detectors at different sites. MUTE is an open-source, publicly available program, providing a solid framework for accurate muon flux predictions in various underground environments, essential for applications in cosmic ray physics and dark matter searches.


I. INTRODUCTION
For many decades, measurements of muons have been conducted underwater and underground to study the energy spectrum and composition of primary cosmic rays [1,2].Despite these long-standing efforts, large uncertainties surrounding atmospheric muon and highenergy neutrino fluxes persist, with sustained challenges to their reduction [3][4][5].Shielded by the Earth's natural overburden, large-volume underground detectors offer a unique opportunity to study high-energy muons without the technical challenges or high cost of single-purpose surface-based spectrometers capable of resolving multihundred TeV range energies.
Addressing these uncertainties is of paramount importance for several reasons.Precise estimation of the muon flux can significantly benefit fields like rare-event searches, including direct dark matter searches.The detectors for these searches are typically installed deep underground in mines or beneath mountains in order to limit exposure to cosmic-ray muons [6].However, cosmic-ray muons produce secondary neutrons by interacting with the rock or surrounding detector materials, which pose significant identification challenges, as they can mimic low-energy dark matter signals.Knowledge of the underground muon rate, energy spectrum, and angular distribution is necessary to estimate muoninduced backgrounds to design effective shielding, including muon and neutron vetos.For the past few decades, analytical calculations [7,8] and empirical parametric fits to vertical-equivalent muon intensity data [9,10] have served as inputs to Monte Carlo muon transport codes such as mum [11], music and musun [12], and Geant4 [13].However, these methods often demand significant computing resources and do not thoroughly address the treatment of errors.fluka [14], meanwhile, can be used to calculate muon fluxes starting from a cosmic ray flux parametrization; however, it also comes with computational cost, in particular for high energies.
In our previous work in Ref. [15], we concentrated on the development of the open-source mute code, designed to calculate underground muon fluxes and other related observables without using those fluxes as input.mute is flexible, efficient, and precise, and combines two modern computational tools, mceq [16,17] and proposal [18], for surface muon fluxes and underground transport respectively.We previously compared verticalequivalent underground intensities computed with mute using a flat overburden to experimental measurements.In the present study, we expand on this by scrutinizing the residual body of related data using an updated version of mute in combination with the latest and most precise surface flux model, the daemonflux model [5].daemonflux makes use of two data-driven models to achieve uncertainties on surface fluxes of less than 10% up to 1 TeV: Global Spline Fit (gsf) [19] for the primary cosmic ray flux, and the Data-Driven Model (ddm) [20] for hadronic interactions.In this work, we calculate total underground and underwater muon fluxes, as well as underground angular distributions, energy spectra, and mean energies.We use topographic maps for laboratories situated under mountains and we incorporate the density and composition of the rock above each lab into our calculations, which are significant sources of systematic uncertainty.We compare our mute predictions with experimental data from different underground and underwater sites worldwide.The new release of mute is publicly available so our results can be easily reproduced and used for further studies. 1

II. COMPUTATIONAL METHOD
For practical applications, the muon intensity at a certain depth underground has been historically parametrized using functions termed depth-intensity relations (DIR).D.-M.Mei and A. Hime [10] give the following DIR for depth ranges between 1 and 10 km.w.e.: where h is the vertical depth of the lab, and I 1 , I 2 , λ 1 , and λ 2 are determined by fitting to experimental data.Comparatively, many collaborations [21][22][23][24][25] give the following three-parameter relation: where A, h 0 , and α are fit from data.Parametric fits like these can provide reasonable estimates for underground labs, as demonstrated in Ref. [10].There are, however, several caveats.Firstly, the fits are done to verticalequivalent underground intensity data, which is known to only approximately correspond to the true vertical intensity at zenith angles below ∼ 30 • [1,7,15].In addition, various datasets are typically combined without compensating for systematic uncertainties through correction functions or nuisance parameters, and hence the fit results might be significantly biased towards experiments that have underestimated their uncertainties.They are also limited in the sense that they cannot provide additional physical information about the muons underground, such as their energy and angular distributions, which must be calculated using separate parametrizations.The computational scheme of mute (v1.0.1) used for labs under flat overburdens is described in Ref. [15].The central quantity in the computation scheme from which all other underground and underwater observables are calculated is the underground flux, Φ u .This is calculated as a convolution between the surface muon fluxes from mceq or daemonflux, Φ s , and a surfaceto-underground transfer tensor from proposal, U , and is given by where E u is the underground muon energy, X is the slant depth, θ is the zenith angle, E s is the surface muon energy, and * denotes a convolution over the surface muon energies.The slant depth (or grammage) X is the primary measure of distance in our calculations.Generally, it is defined as where the integration is performed along the line of sight ℓ ′ from the center of the detector or lab, situated at the distance ℓ from the intersection with the surface in a specific direction, which is defined by θ and the azimuthal angle ϕ.For homogeneous rock or for a given average density, this expression simplifies to The second equivalence is valid for flat overburdens located at a fixed vertical depth d in km.For convenience, it is common to convert distances and slant depths into kilometer water equivalent units (km.w.e.= 10 5 g cm −2 ): X km.w.e.= X/ρ water h km.w.e.= dρ/ρ water .
The muon flux at the surface is assumed to have azimuthal symmetry.Furthermore, the altitude dependence of the surface flux can be safely neglected in our calculations since muons with energies relevant for underground fluxes at h ≥ 0.5 km.w.e.(i.e.energies above ∼ 100 GeV) are produced in the upper atmosphere, around 10 km above sea level.For these reasons, we consider Φ s to be a function only of E s and θ.The conventions for definitions and notations used for "flux" and "intensity" in our work are explained in Appendix A.
In this followup study, we introduce mute v2 [26] and v3, capable of calculating muon observables for laboratories beneath mountains.For the case of a flat geometry above the lab (labs under flat overburdens or underwater), due to the assumed symmetry in the azimuthal angle, Eq. ( 5) provides a direct geometric correspondence between the slant depth and zenith angle for a given vertical depth.This relationship reduces the underground or underwater flux in Eq. ( 3) at the specified depth to a function of just two variables: Φ u (E u , θ).For mountains, however, the nonuniform shape of the overburden introduces dependence on the azimuthal angle to the FIG. 1. Slant depths of the Gran Sasso mountain, X(θ, ϕ), in terms of the zenith and azimuthal angles [24].Depths greater than 14 km.w.e. have been masked (shown in black), as 14 km.w.e. is the default maximum slant depth for calculations by mute.In all cases, muons traveling along θ = 0 • are defined as vertical down-going muons.
amount of rock a muon has to travel through, meaning the geometry of the mountain has to be taken into account in the calculations.To do this, topographic maps of the mountains in terms of the depth in spherical coordinates, X(θ, ϕ), are obtained from the labs, such as the one for Laboratori Nazionali del Gran Sasso (LNGS) in Italy, shown in Fig. 1.For this reason, the slant depth and the zenith angle are kept as separate variables in the calculation of underground fluxes for mountains: Φ u (E u , X(θ, ϕ), θ).The explicit dependence on θ is relevant because the surface fluxes depend on the zenith angle but not the azimuthal angle due to the symmetry at high energies.
For both overburden types (flat and mountainous), the differential underground muon intensity, I u , is given by integrating the underground flux over the underground energy: where E th is a choice for a threshold energy, typically defined by a specific detector.Since the typical energies of muons surviving the transit to deep underground labs end up in the GeV range underground, we set E th to zero by default for all labs for generality.For a flat overburden lab situated at a fixed vertical depth, Eq. ( 7) reduces to a function of one variable, I u (θ).It remains a doubledifferential quantity for mountains as in Eq. ( 7), with X(θ, ϕ) given by the topographical map.We implement the calculation of the double-differential intensity in mute v2 by precomputing an underground intensity matrix I u (X, θ) on a constant grid of depths X and zenith angles θ.We then interpolate this distribution to a specific overburden profile X(θ, ϕ) to obtain a matrix of lab-specific intensities I u (θ, ϕ).The elements of this matrix correspond to the bins of the mountain maps X(θ, ϕ) provided directly by the experiments.Mountain [29] Using these underground intensities, we calculate various physical observables for different underground sites worldwide.The underground sites used are listed in Table I, along with their average rock density and vertical depths in km as found in the literature.We present flat overburden calculations for the Waste Isolation Pilot Plant (WIPP), Soudan Underground Laboratory, Boulby Underground Laboratory, Stawell Underground Physics Laboratory (SUPL), Sanford Underground Research Facility (SURF) located in the former Homestake Gold Mine, and the Sudbury Neutrino Observatory Laboratory (SNOLAB).For laboratories under mountains, we have used topographic maps of the Yangyang Underground Laboratory (Y2L) centered around the COSINE-100 detector, the Kamioka Observatory under the Ikenoyama mountain centered around two locationsthe KamLAND [27] and Super-Kamiokande [27,28] detectors-LNGS under the Gran Sasso mountain centered around the Large Volume Detector (LVD) [24], Laboratoire Souterrain de Modane (LSM) centered around the Fréjus detector [21], and China Jinping Underground Laboratory (CJPL-I) centered around the Jinping Neutrino Experiment (JNE) [29].In the case of labs under mountains, it is necessary to center the map around a specific detector's location due to the calculations' dependence on the azimuthal angle.We selected these labs based on the availability of experimental data with which we can compare our results.
By default, mute provides surface-to-underground transfer tensors for muon energy loss in standard rock (Z = 11, A = 22, ρ = 2.65 g cm −3 [2,41]), fresh water (ρ = 1.000 g cm −3 ), and sea water (ρ = 1.03975 g cm −3 [18]) for a range of slant depths spanning 0.5 km.w.e. to 14 km.w.e.The latest release, mute v3, provides explicit transfer tensors for the specific rock density and composition at the locations of the labs listed in Table I, resulting in a more precise and reliable muon flux calculation.

III. MODELING ROCK COMPOSITION
The flux of muons underground and underwater is governed largely by the energy losses of the muons as they travel through rock and water.When traveling through matter, high-energy muons relevant to underground experiments lose energy via pair production, bremsstrahlung, photonuclear interactions, and ionization.The cross sections for these interactions depend on the medium's chemical composition, making accurate modeling of media a significant factor in flux calculations.The effects of each of these and how they are implemented into our calculations with mute v3 are briefly explained in the following subsections.

A. Rock density
All internal proposal calculations for energy loss are performed in units of grammage, meaning they are agnostic to changes in the density.The rock density, therefore, only enters the calculations when converting physical depths in km to grammage in km.w.e. via Eq.(6).
For the labs under flat earth in Table I, the rock density has been accounted for in the lab's vertical depth, h.This has sometimes led to different vertical depths being used in this work compared to those found in the literature.For the labs under mountains, some mountain maps were already available in units of km.w.e., in which case no conversions were made.However, when the slant depths in the mountain profile maps were given in units of km, they were converted into km.w.e. of laboratory rock by multiplying by the corresponding rock density in Table I.This simple linear scaling of the slant depths is an estimate of the impact of the average rock density alone, as chemical differences between different rock types are considered only for energy losses and not when calculating depths.This is in contrast to parametric polynomials given by some experiments, such as LVD [24] and SNO [25], for example.These experiments use underground muon intensity data to construct conversion formulas from standard rock to laboratory rock, which aim to account for all differences between the rock types when converting depths.
Although rock density is essential for calculating total muon fluxes underground, many labs have significant uncertainties on their rock density, or no uncertainties at all, as seen in Table I 2. Total underground muon fluxes as a function of rock density under the Gran Sasso mountain centered around LVD in Hall A of LNGS.The error bars represent the uncertainty from daemonflux.The horizontal lines and their error bands represent the measurements from MACRO (Hall B) [42], Borexino (Hall C) [43], and LVD (Hall A) [44].The vertical line and its error band represent the measured rock density of (2.72 ± 0.05) g cm −3 , given in Table I.
misjudge its rock density if too little information on the rock composition is available.Because of the relation between rock density and underground muon flux and the precision of the results produced by mute, there is potential for mute to constrain uncertainties on rock densities or even correct misjudged rock densities.We have done this for LNGS in Fig. 2, and compare the results to three experimental measurements from MACRO, Borexino, and-the most recent LNGS measurement-LVD.The good agreement with LVD suggests that it is likely possible to further constrain rock densities for other labs using mute; however, we note that the precision of the surface flux model is not yet good enough to resolve the slight differences between different experimental halls.

B. Rock composition
The cross sections for muon energy loss interactions are dependent on the average number of protons, ⟨Z⟩, and the average atomic mass, ⟨A⟩, of the propagation medium.The cross sections are proportional to Z 2 /A for pair production [45] and bremsstrahlung [46,47], to ⟨A⟩ for photonuclear interactions [48,49], and to ⟨Z⟩ for ionization [50].Although overburdens are typically made up of multiple types of rock, for the optimization of the program, homogeneous media defined by average ⟨Z/A⟩ and Z 2 /A values have been used in proposal.The definitions of these average values are given in Appendix B, and the values used for each lab are listed in Table II.
TABLE II.Rock composition definitions in terms of average numbers of protons and nucleons.In some cases, multiple contrasting values for ⟨Z⟩ and ⟨A⟩ are found in or derived from the literature, as is the case with Soudan (see Refs. [10,32,51]), Kamioka (see Refs. [27,52]), LSM (see Refs. [21,53,54]), and SNOLAB (see Refs. [25,55,56]).In these cases, the set of values for which the most detailed information was provided was used.⟨Z/A⟩ and Z 2 /A values are sometimes found directly in the literature, though, in some cases, only percent weights for minor component chemical compositions are found.In these latter cases, the average numbers of protons and nucleons were calculated from these compositions.All available chemical compositions are listed in Table V in Appendix B. Lastly, sometimes values are published as ⟨Z⟩ and ⟨A⟩ instead of ⟨Z/A⟩ and Z 2 /A .In general, ⟨Z/A⟩ ̸ = ⟨Z⟩ / ⟨A⟩, meaning algebraic conversion from one pair of values to the other is not possible; however, the approximation error for this inequality is sufficiently small (less than 1% in all cases) that it is neglected, and we approximate ⟨Z/A⟩ = ⟨Z⟩ / ⟨A⟩ where needed to fill Table II.

Laboratory
A detailed geological survey of the Stawell, Victoria area in Australia was published in Ref. [61], but no information on the average rock composition is provided.It is known that the rock above SUPL is mainly basalt [35,36,62], but the chemical composition of different types of basalt can vary widely, with the percent weight of SiO 2 lying between 45% and 52%, which can lead to significant variations in the energy loss and therefore total underground fluxes.Due to a lack of precise geochemical measurements, this work used standard rock for all SUPL calculations.
The effect of changing the ⟨Z/A⟩ and Z 2 /A values on the energy loss results from proposal can be significant.Figure 3 shows the ratio of underground muon intensities using the laboratory rock from Table II to those using standard rock, demonstrating that increasing the Z 2 /A value leads to a decrease in the muon intensity results, as physically expected.There is a deviation from standard rock of -27% in the underground intensity for CJPL-I at the vertical depth of the lab (approximately 6 km.w.e.), while the effect is minimal for Y2L, Kamioka, and LNGS at their respective vertical depths (less than -5%).However, the effect is amplified as slant depth increases, as the muons travel through higher quantities of rock.Therefore, muons reaching the labs at higher zenith angles are affected more greatly by changes to the chemical composition.Because of this depth-dependent energy loss effect, DIRs and mean underground energy parametrizations that are not corrected for rock type, like those in Eqs. ( 1) and (2), are less suitable for highprecision calculations, highlighting the need for detailed simulations in underground muon flux calculations.
Lastly, it should be noted that because the values presented in Table II are averages, and overburdens typically contain various types of rocks in different quantities, the ⟨Z/A⟩ and Z 2 /A values will contain some margin of error that is not taken into account.While these average values are assumed to be sufficient approximations for the definitions of the propagation media, more detailed information about rock compositions can be helpful for more precise results.

C. Sternheimer parameters
To describe ionization losses, proposal uses the Bethe-Bloch equation [2].This equation includes a density correction term, −δ/2, as described in Refs.[18,63], where δ is given by Here, x 0 and x 1 demarcate transition regions of the function form of the δ parameter [63].δ 0 , a, b, and c are additional empirical parameters, all together termed the "Sternheimer parameters."Values for these parameters are typically taken from [18,[64][65][66], and Table VI in Appendix B lists the values for the media used in this study.
Although important for energy loss, the only relevant rock type for which published Sternheimer parameters exist in the literature (in addition to standard rock) is Fréjus rock, the rock above LSM [18].For this reason, the calculations for LSM presented in this work were done using the full definition of Fréjus rock.For other laboratories, however, because values for the Sternheimer parameters do not exist in the literature, all calculations with mute use the Sternheimer parameters for standard rock as an approximation.This is with the exception of WIPP, which uses an overburden of salt.
The difference in energy loss between standard rock and Fréjus rock is shown in the bottom panel of Fig. 3.This panel shows the effect of using the Fréjus ⟨Z/A⟩ and Z 2 /A values with the dashed line, the impact of using the Fréjus Sternheimer parameters with the dotted line, and the effect of using both with the solid line.The Sternheimer parameters are shown to contribute a nearly constant 10% decrease in underground muon intensity across all depths.Therefore, we estimate at least a 10% systematic error on predictions for other laboratories due to the lack of knowledge of Sternheimer parameters.

IV. TOTAL MUON FLUX
The total muon flux is the main physical observable of interest for underground and underwater experiments for the purposes of muon-induced background studies.It is calculated by integrating the underground or underwater intensity from Eq. ( 7) over the solid angle: We have performed this calculation for labs underground and underwater and present the results below.

A. Total flux underground
The total underground flux has been calculated for the sites listed in Table I using daemonflux as the surface muon flux model and U.S. Standard Atmosphere [67] as the atmospheric density model.The results are given in Table III and are shown in Fig. 4 along with experimental measurements.
In Fig. 4, the total underground flux is plotted against "equivalent vertical depth."For labs under flat overburdens, this is the depth of the lab, comparable to depths listed in Table I.For labs under mountains, however, the equivalent vertical depth is defined, as in Ref. [10], as the depth a lab would be at under a flat overburden given the total underground muon flux seen by the lab.We allow the mute points for labs under mountains in Fig. 4 to float along the x axis and fit their depths so they lie directly on the mute curve, which has been calculated for standard rock.We report these inferred depths in Table III and label them h SR .The uncertainties on these values come from the error band on the mute curve.
The official data points shown in Fig. 4 for labs under mountains are seen to be offset from the mute result along the x axis, such as the CJPL-I lab under the Jinping mountain, whose depth we infer to be (6.20 ± 0.13) km.w.e.This is notably different from the value of 6.72 km.w.e.quoted in, for example, Refs.[29,74], TABLE III.Experimental and predicted total muon fluxes for underground laboratories.The predicted total flux values were used to infer equivalent vertical depths by allowing the mute calculation to float along the x axis and fitting it to the curve for standard rock in Fig. 4, resulting in a standard-rock equivalent vertical depth, hSR, for all labs.which is stated to be the result of a scaling factor F from cosmic ray leakage due to the topographic profile of the mountain.6.72 km.w.e.is, in fact, the straight vertical depth above the lab (termed the "engineering depth"), not the equivalent vertical depth as used in our work.This is also the reason for the offset between the experimental and calculated points for KamLAND and Super-Kamiokande, while a similar offset is present for LNGS and LSM caused by the use of average vertical depths in the literature.The inconsistencies in the literature in defining a meaningful vertical depth for laboratories under mountains make the results difficult to visually interpret.Therefore, Fig. 5 shows the ratio of the results in Fig. 4 to the mute predictions in Table III independent of lab depth.
The main uncertainties in mute come from the surface flux model, daemonflux, and the modeling of the overburden.The former are represented by the gray error bands in Figs. 4 and 5, and contribute an error of ∼7% for deep depths (higher muon energies) and around 2% for shallow depths.The latter are represented by horizontal gray bars and are entirely attributed to the range of allowed average rock densities listed in Table I.These two errors are shown separately to emphasize their independence from each other.The interaction model uncertainties are independent of the lab, in contrast to the knowledge of the rock density, which is either lim-ited by the precision of the geological surveys or by the nonuniformity of the overburden, making it difficult to compute a representative mean density.The interaction model errors have been discussed explicitly in Ref. [15], where the rock density variations were implicitly included in the errors of the experimental data.The magnitudes of the uncertainties in Table I are large in many cases, leading to significant variations in the total flux seen in Fig. 5. Error bars are not shown for labs where no uncertainty on the rock density was found in the literature.The mountain maps cited in Sec.II are assumed to be exact, so no error was considered for the slant depth values.Lastly, the proposal Monte Carlo simulations in mute were performed with N = 10 6 events per energy and zenith angle bin, making the statistical error across all calculations negligible.
In general, discrepancies between predictions and measurements are not unexpected because the measurement conditions are not always possible to fully reproduce.Many factors can have a significant effect on the results, including the composition of the rock, the location of the lab under the mountain, the location of the detector in the lab, the energy threshold of the detector, the angular acceptance of the detector, and analysis decisions like energy cuts, binning choices, whether the analysis is inclusive of muon bundles or just single muons, and the treatment of seasonal variations in the collection  III at the vertical depths of the labs quoted in the literature.These depths are sometimes given as minimum, average, or straight vertical depths ("engineering depths"), rather than the equivalent vertical depths used in this work.The curve from the parametric formula in Eq. ( 1) from Ref. [10] is also shown for comparison.and analysis of the data.More detailed information and analyses of these factors for a given detector may further improve the accuracy of mute results.Although these details were not always available for each laboratory for this study, mute is provided as an open-source tool for experiments to tailor to their specific use cases.
Despite the many potential factors contributing to the accuracy of the calculations, mute in combination with daemonflux provides an accurate description of the underground muon flux data within uncertainties for nearly every lab.The precision of the calculations of the total underground flux and the equivalent vertical depth in Table III are good enough, for example, to resolve the difference in location between the KamLAND and Super-Kamiokande detectors under the Ikenoyama mountain.This agreement between prediction and measurement is also indicated when comparing our results to measurements of the muon seasonal variation amplitudes in one of our previous works [75].As a result, the total underground muon flux may be one of the most robust observables to use to constrain primary cosmic ray flux models or hadronic yields because the results are modelindependent.However, because of the magnitudes of the systematic uncertainties coming from the rock density errors in Fig. 5, the constraints provided to flux model calibrations might not be as strong as they potentially could be [5,76] II  and VI).A dagger ( † ) indicates that our calculation used standard rock rather than being tuned to the rock above the given laboratory.A double dagger ( ‡ ) indicates that the point used is a prediction calculated from simulation and does not come from experimental data.Where multiple data points are shown, the left-to-right order of the points is the same as the order of the measurements in brackets in the labels.Similar ratios produced with other surface flux models can be found in Appendix C.
In the next subsections, we discuss in more detail the results shown in Fig. 5 for a subset of laboratories.

Y2L
The clearest outlier in Fig. 5 is the measurement for Y2L.Multiple measurements have been published for the COSINE-100 detector in the A5 hall of Y2L [31,68,77], all of which are in close agreement with each other but in disagreement with the mute prediction by about 20%, according to Fig. 5.An additional measurement from the KIMS experiment has been published in Refs.[78,79], but KIMS is located in the A6 hall of Y2L, which is displaced 200 m away from the A5 hall [77].As seen with the difference in mute results for the KamLAND and Super-Kamiokande detectors, which are spaced 150 m apart in the Kamioka laboratory [80], mute calculations are sensitive enough to capture these differences in the positions of detectors.Therefore, because the mountain map we received for Y2L is centered around the COSINE-100 detector in the A5 hall, we compare only to the most recent COSINE-100 measurement from Ref. [68].
A possible explanation for the discrepancy between our prediction and the COSINE-100 result is muon bundles.The flux calculated by mute corresponds to that of single muons; however, the COSINE-100 detector and analysis make no distinction between single and multiple muons.
According to Table III, the equivalent vertical depth of Y2L is 1.58 km.w.e.Using mute in combination with mceq, we approximate that, at this depth, 15% of events are expected to contain two or more muons per event.An excess of events in Fig. 9 in Ref. [77] suggests an approximate contribution of 10% due to muon bundles captured by the detector within detector acceptance.This implies a total muon flux that would be 15-20% higher than reported.We have computed a correction for this in Table III and have indicated this corrected value by a pale pink point for Y2L in Fig. 5.
It should be noted that the agreement we find is highly dependent on the rock density and vertical depth used for the lab.There is, however, wide variation in the depths of the lab quoted in the literature; see Refs.[32,51,69,[82][83][84], all of which give different vertical depths for the same level of the laboratory, ranging from 0.690 km to 0.780 km, and Refs.[51,69,[82][83][84][85][86], which give vertical depths ranging from 2.07 km.w.e. to 2.10 km.w.e.These differences in depth lead to a maximum possible discrepancy of over 20% between experimental measurements and our mute prediction in some cases.Therefore, for consistency, in Table I and in our calculations, we used the Soudan 2 density and vertical depth in km given in Ref. [32] to calculate a vertical depth in km.w.e.We then used that vertical depth in km.w.e. to calculate the total flux for a flat overburden and we compare that result to the measurement from Ref. [69].
Furthermore, as both Ref. [69] and Ref. [32] note, muon intensity changes as location in the laboratory, altitude, and rock density change, and there is significant variation in the rock density and composition above Soudan.For this reason, more accurate information on the overburden is needed for a detailed comparison.Although Fig. 10 of Ref. [69] shows a variation in the mountain profile above the lab of around 0.1 km, which is enough to have a great effect on mute results, we performed the calculation for a flat overburden because we had no access to this topographic map.Despite this, we still see overall very good agreement with the Ref. [69] measurement.

SUPL
As stated in Sec.III, aside from the rock density, no additional relevant quantitative information exists about the rock above SUPL.For this reason, standard rock was used in the proposal simulations.Despite not considering the rock composition, the mute result agrees well with the SABRE measurement from Ref. [70].We note that the result may change with more information about the rock density and composition.

SURF
Three total underground muon flux measurements have been published for SURF, from Homestake [72], MAJORANA [38], and LUX [73], all taken at the 4850level Davis laboratory location.
Detailed information on the rock's chemical composition above SURF is available in the literature [39,59,87,88].However, because of the significant variation in different rock types, with densities ranging from 2.43 g cm −3 to 3.26 g cm −3 , the uncertainty on the rock density remains high.With the uncertainties from the daemonflux model, the rock density, and the experimental systematics all taken into account, we find the best agreement with Homestake and the worst with MAJORANA.
It can be seen in Fig. 5 that these three measurements are not in good agreement with each other.Possible reasons for these disagreements are discussed in Ref. [73].For the case of the Homestake measurement, Ref. [72] finds a single-muon intensity of (4.91 ± 0.06) × 10 −9 cm −2 s −1 sr −1 for zenith angles up to 18 • where the uncertainty is solely statistical.Reference [73], however, corrects this to a value of (4.14±0.05)×10−9 cm −2 s −1 by considering both multiple muons and an angular range of 0 • ≤ θ ≤ 90 • .We estimate the fraction of muon bundles at the depth of SURF to be around 10%.This is in contrast to the Homestake-measured fraction of 4.4% [72].If all calculations were corrected consistently for muon bundles, we may see better agreement.
Additionally, despite all measurements being taken at the 4850-level of the mine, another factor is the variation in elevation of nearly 0.57 km.w.e. in the overburden above SURF (see Fig. 5 in Ref. [38]), as well as the lateral separation in locations between the Homestake and MAJORANA detectors, which, as discussed, can significantly affect mute results.

SNOLAB
Although agreement can be achieved between the mute prediction and the SNO data point from Ref. [25] when all systematic uncertainties are considered, it A curve is also included for the theoretical calculation of Bugaev, et al. [7].KM3NeT data is taken from Ref. [89].
should be noted that this data point is not well understood, and, like for other laboratories, agreement is dependent on the depth we consider for the lab.Despite some ambiguity about the correct depth to use for SNO-LAB, in our calculations, we take h = 5.920 km.w.e. based on the depth d = 2.092 km quoted in Ref. [25].
In general, mute in combination with daemonflux seems to overpredict the total muon flux at the depth of SNOLAB, which might be consistent with what we observe for the vertical intensity in Ref. [15] using ddm and also what we observe for CJPL-I in Fig. 5 at a similar depth.Because deep depths are associated with high surface energies, this may suggest that daemonflux does not calculate muon fluxes at high energies with high accuracy.Agreement improves when considering other hadronic interaction models, as shown in Fig. 11 in Appendix C.
Because a description of SNO rock in terms of Sternheimer parameters is not found in the literature for a full simulation of the energy loss, and because Ref. [25] provides the only published measurement of the total underground muon flux for SNOLAB, further investigation is not possible at the moment.However, future measurements of the total muon flux at SNOLAB may help resolve this issue.

B. Total flux underwater
Total muon flux measurements in water can be performed by water Cherenkov detectors, such as the KM3NeT detectors [89].The comparison between mute, a typical reference calculation [7], and an early KM3NeT measurement using a single Detection Unit is shown in Fig. 6.
When performing computations for laboratories under-water as opposed to under a flat overburden, the only change in the mute calculation is to the medium used for the surface-to-underground transfer tensor generation with proposal, from rock to water.Thus, the accuracy of the mute predictions in Fig. 5 seen for rock should apply to water as well.We note that KM3NeT's data is a factor of 2 below our prediction for sea water.A more recent measurement of the event rate as a function of the zenith angle from a more complete version of the KM3NeT detector [90] indicates a 30-40% discrepancy to a calculation made using a similar technique and the sibyll-2.3dinteraction model.Their result is approximately in agreement with our findings illustrated in Fig. 11 (see Appendix C) and also confirms that the previous results from Ref. [89] shown in Fig. 6 were significantly underestimated.

V. ANGULAR DISTRIBUTIONS
We have used data from the LVD experiment located at LNGS to check the consistency of mute predictions for the angular distribution of muons underground.The underground intensities from the LVD data, I u LVD , were calculated according to where K is a normalization constant representing the product of the detector lifetime and effective area to convert the I u LVD values to physical units, n(θ, ϕ) is the number of raw counts, and ε(θ, ϕ) is the acceptance of the LVD detector [24].K was calculated by requiring the total integrated underground fluxes from mute, Φ u tot , and from the LVD data to be equal: The intensities were similarly computed using mute with the method described in Sec.II, and the results are shown in Fig. 7 (top).In order to compare the experimental and predicted intensities, the residuals for each (θ, ϕ)-bin of the data were calculated as where δI u LVD are the statistical errors in each angular bin.The result is shown in Fig. 7 (bottom).Since, by construction from Eq. ( 11), the normalizations of the mute result and the LVD data match, only relative trends can be observed.The calculation is slightly below the data at near-vertical directions and above it for intermediate zenith angles, except for a small patch in the SW direction (∼ 225 • azimuth), where the situation is inverted.We calculate the χ 2 /ndf to be 1.34,where ndf is assumed to be the number of nonzero bins (ndf = 31756).
We first tested if we could obtain a better description of the data by allowing for a variable E th in our computation of the intensities with Eq. ( 7).However, we found that the prediction and the data were independent of small variations in E th .We therefore investigated whether additional systematic uncertainty would be needed to accommodate this minor disagreement by exploring the expected surface azimuthal symmetry.We binned the observed intensities in slices of equal zenith angle and slant depth and compared the widths of the distributions between data and simulation.We expected to see the same flux within these bins with minor deviations coming from uncertainties in the data and calculations, which would result in narrow distributions.We found, however, that the distributions from the data were always wider than those from the simulation.We con- clude, therefore, that the observed χ 2 /ndf is an effect of this additional systematic rather than an indication of a mismatch between model and data.Given this, the disagreement between mute and LVD is not concerning for predicting the total flux.
The one-dimensional projections of the angular distributions in the zenith and azimuthal directions have also been calculated using Once the distribution in Fig. 7 (top) is projected onto the zenith and azimuth axes (Fig. 8 top and bottom, respectively), we observe good compatibility of the result from mute within the small errors of the data, in particular for the zenith projection at near-horizontal directions.The worse agreement in the patch of positive residuals in the (θ ∼ 40 • , ϕ ∼ 225 • ) region in the bottom panel of Fig. 7 is also visible in the bottom panel of Fig. 8.This is an azimuthal range which contributes greatly to the flux, therefore affecting the χ 2 /ndf, but the overall agreement we find between the mute result and the LVD data is otherwise good.In general, comparing the experimental and calculated angular distributions can reveal critical information about the accuracy of data analyses of underground experiments.The overall agreement obtained with the LVD experimental data suggests that mute can be a powerful and helpful tool to inform on potential sources of errors and misreconstructions in data analyses for underground experiments, mainly due to the absence of statistical errors.

VI. UNDERGROUND SPECTRA AND MEAN ENERGIES
The muon energy spectrum for a specific underground lab, Φ u Ω , is defined by: In the case of labs under mountains, the zenith and azimuthal angles are those provided by the grid of the mountain map, whereas, for flat overburdens, Eq. ( 5) is used.Energy spectra calculated with mute for labs under flat earth and under mountains are shown in Fig. 9.The curves' vertical ordering corresponds to the labs' depths, with some influence from the individual labs' rock compositions.In general, however, as expected, the energy spectrum decreases as depth increases.The spectra were verified by integrating the curves and comparing the results to the total underground muon flux for the given lab found with Eq. ( 9), which is calculated via a different computational method, and the results matched within 1%.
Little data in physical units is available in the literature to compare our calculated energy spectra to, though there is a simulation result for Kamioka in Ref. [91], shown by the thick gray line in Fig. 9, which closely matches the mute curve (the difference is less than 10% for the majority of the displayed energy range).
The mean underground muon energy is also relevant to underground detectors, particularly in stopping vs through-going muons studies.Historically, it has been calculated by means of a parametrized energy spectrum, where A is a normalization constant with respect to the differential muon intensity at a given depth h [2,10,92].
Values for the parameters ϵ µ , b, and γ µ for standard rock are listed in Ref. [10] for the parametrizations from both  15).The energy spectrum given in Ref. [91] for Kamioka from a music/musun simulation is shown by the thick gray curve.
The mean energy is given by the first raw moment of the underground muon energy spectrum: where the integral in the denominator is the underground intensity, I u (X, θ), from Eq. ( 7).We have done this calculation using the energy spectra from both Eqs. ( 15) and ( 16) for the underground sites listed in Table I, with the results given in Table IV and plotted in Fig. 10.We have used the Lipari and Groom parameter sets with the h SR depths inferred from mute given in Table III.Because these parameters are for standard rock, we also provide mute results for standard rock at the h SR depths for comparison.Lastly, full mute predictions are provided, using the rock compositions listed in Table II and the labs' depths from Table I for flat overburden labs and the topographic maps for mountain labs.Overall, within the widths of the distributions, the results are compatible with each other.Additionally, we find agreement between our calculation for LNGS using the laboratory rock and the MACRO measurement from Ref. [94].

VII. CONCLUSION
We have presented a comprehensive compilation of the latest available references related to underground and underwater muon flux studies, and have compared to results from the muon intensity code mute to explore the physics of cosmic-ray muons.The latest version of our open-source code, mute v3, offers the ability to compute  IV.The mute (Full) points are plotted at the mute-inferred hSR depths given in Table III to be comparable to the mute (SR) curve.The error bars on the mute (Full) points and the error band on the mute (SR) curve represent the uncertainty coming from the daemonflux model.The MACRO point is a measurement for LNGS of (270 ± 18) GeV from Ref. [94].
underground energy spectra and angular distributions as well as the ability to set custom rock compositions.It furthermore expands beyond calculations for labs under flat overburdens by now offering all calculations for sea level and labs under mountains as well.As it is both fast and flexible, mute can be a valuable tool for laboratories to efficiently study their overburdens and data with little computational strain.
Using the data-driven daemonflux model, we achieve a remarkable accuracy, particularly at depths below 5 km.w.e., confirming consistency between surface and underground flux measurements.Moreover, our findings indicate that hadronic interaction models like sibyll-2.3dunderpredict muon fluxes in the TeV surface energy range, something that has also been observed by underwater measurements from the KM3NeT detector.
Our findings highlight the importance of both the density and the chemical composition of the rock overburden above underground labs.Detailed knowledge of these is crucial in order to achieve an accuracy of significantly below 15% in underground muon flux studies, and it is insufficient to use standard rock or parametric formulas.Therefore, we encourage experiments to study and publish details on the rock above their labs, namely average rock densities with uncertainties, ⟨Z/A⟩ and Z 2 /A values (or minor or major component geochemical measurements), and Sternheimer parameters.Throughout our work, we also found that many laboratories have varying reported vertical depths in the literature.To achieve high-accuracy results, we find that it is useful to have access to published topographic maps of overburdens in terms of X(θ, ϕ) in km.w.e., rather than only an average density, even for those labs under relatively flat earth.Additionally, it is helpful to publish intensity and flux data that has not been corrected to standard rock, allowing for external analyses to apply their own corrections in a consistent and transparent way.Lastly, we encourage laboratories to publish more detailed muon flux measurements, systematic uncertainties, and angular distributions whenever possible, even for decommissioned detectors.
This work has demonstrated that mute is a suitable tool for interpreting underground muon measurements.It also made clear the importance of accurate modeling of the characteristics of each laboratory's overburden in order to achieve sufficiently high precision for indirect studies of muon and neutrino flux uncertainties in the TeV to PeV range at the surface.

ACKNOWLEDGMENTS
We acknowledge the help of Marco Selvi and the LVD collaboration, who provided us with the topographic map of the Gran Sasso mountain and the LVD data for the muon angular distribution.We acknowledge, as well, the help of Michel Zampaolo and Luigi Mosca for the data of the Fréjus detector, and Holger Kluck and Wolfgang Rhode for a map of the Fréjus mountain.We also acknowledge Eunju Jeon and Hyun Su Lee, who provided us with a topographic map of the Yangyang Mountain and discussed with us muon bundles in the COSINE-100 detector, Itaru Shimizu and Shigetaka Moriyama for the maps of the Kamioka Mountain, and Shaomin Chen for the map of the Jinping Mountain.We want to thank Joseph Formaggio, Christopher Kyba, and Chris TABLE V. Percent weights (%) of minor component rock compositions for laboratories considered in this study, where available.Rock compositions were unavailable for WIPP, Boulby, and SUPL.In some cases, percent weights might not sum to 100%, as some minor components that do not contribute significantly to the makeup of the rock have been excluded from reporting.a Chemical compositions are available as reported here, but the ⟨Z/A⟩ and Z 2 /A values from Ref. [32] for Soudan and Ref. [21] for LSM were used instead of those computed using this table.b Multiple rock samples from Ikenoyama Mountain above the Kamioka lab are listed in Ref. [52].Here we have used the average of the two JR-1 and JA-3 igneous samples, which are widely distributed around the Kamioka area, as given in Ref. [58].This result is consistent with the recent underwater measurement by KM3NeT [90], where a simulation based on gsf and sibyll-2.3dis 30-40% lower compared to the observed event rate.The impact of using a hadronic interaction model other than sibyll is shown for the vertical equivalent intensities in Fig. 11 of Ref. [15] and can be expected to be similar in this work.

Laboratory
Media in proposal are modeled by three sets of values: A. The rock density, B. The rock composition, C. The Sternheimer parameters.

FIG. 4 .
FIG.4.Total underground muon flux vs equivalent vertical depth underground.The mute curve was calculated for a flat overburden of standard rock, and the points with thick black outlines were calculated with mute using topographic maps and rock compositions of the mountain.The latter have had their depth values fitted to lie on the mute curve, giving the laboratory's equivalent vertical depth.Data points are shown for the most recent experimental flux values listed in TableIIIat the vertical depths of the labs quoted in the literature.These depths are sometimes given as minimum, average, or straight vertical depths ("engineering depths"), rather than the equivalent vertical depths used in this work.The curve from the parametric formula in Eq. (1) from Ref.[10] is also shown for comparison.

FIG. 6 .
FIG.6.Total underwater muon flux vs slant depth.The mute curve was calculated for a flat overburden using sea water (ANTARES water) as defined in Table VI in Appendix B. A curve is also included for the theoretical calculation of Bugaev, et al.[7].KM3NeT data is taken from Ref.[89].

FIG. 7 .
FIG.7.Underground double-differential muon intensities for the Gran Sasso mountain as calculated by mute (top) and the residuals with respect to the LVD data (bottom) as a function of the zenith and azimuthal angles.In the figure of the residuals, black indicates masked bins for which the slant depth is greater than 14 km.w.e.Additionally, because the binning on the map of the Gran Sasso mountain is very fine, with large bin-to-bin fluctuations in the residuals, a Gaussian filter (σ = 1.5) has been applied to help identify any characteristic trends.

FIG. 8 .
FIG.8.One-dimensional projections of the zenith (top) and azimuthal (bottom) angular distributions for the Gran Sasso mountain as calculated by mute, compared to data from the LVD detector[24].The uncertainty bands on the mute curves come from the error in the daemonflux model.The uncertainties on the LVD points are derived from the acceptance and raw counts data.The total underground muon flux has normalized the spectra to mitigate uncertainties from the hadronic and primary models.

TABLE I .
Summary of underground sites used, their average rock density, ρ, their vertical depth in km, d, and their vertical depth in km.w.e., h.The depth in km.w.e. is calculated as h = ρd.The sites are sorted by increasing equivalent vertical depth.
[57]cause WIPP is a salt mine, its overburden consists primarily of NaCl.Therefore, the full definition of salt, including its Sternheimer parameters in Table VI in Appendix B, was used for WIPP.The ⟨Z⟩ and ⟨A⟩ values listed here are slightly different from those found in Ref.[57](where ⟨Z⟩ = 14.64 and ⟨A⟩ = 30.00)due to the presence of other elements in the rock samples in Ref.[57].b SR .The top panel shows the effects of changing solely the ⟨Z/A⟩ and Z 2 /A values in the proposal simulations, whereas the bottom panel shows the effect of changing both the ⟨Z/A⟩ and Z 2 /A values and the Sternheimer parameters.Error bands in both panels represent the uncertainty from the daemonflux model.Lab depths and mountain profiles are not taken into account.
FIG. 3. Deviation of underground muon intensities calculated using laboratory rock (LR) types from TableIIfrom intensities calculated using standard rock (SR) vs slant depth for four labs under mountains (top) and LSM (bottom).In both panels, intensities calculated using standard rock ⟨Z/A⟩ and Z 2 /A and Sternheimer parameters are represented by I u .FIG.5.The ratio of experimental total underground muon flux measurements to mute, using the results in TableIII.The error band represents the uncertainty in the model from daemonflux, while the horizontal bars represent uncertainty from the rock density values given in TableI.Data points are shown for the experimental flux values listed in TableIII.An asterisk ( * ) indicates the calculation uses Sternheimer parameters other than those for standard rock (see Tables

TABLE IV .
Mean underground muon energies in GeV.The first three columns are comparable to each other, as they show results calculated with standard rock using the hSR values from Table III as depths.The last column provides our full calculation for the expected mean underground muon energy for each lab site, and uses the laboratory rocks as defined in TableIIas well as the depths h from TableIfor flat overburdens and the topographic maps for mountain labs.The uncertainties come from the error in the daemonflux model.

TABLE VI .
[18]nheimer parameters for media used in this work.I is the ionization energy of the medium, and the δ0 parameter from Eq. (8) is 0 for all media.Values are taken from[18].