text:
# Super-Arrhenius slowing down in the metastable liquid of hard spheres

## Abstract

It is demonstrated that a straightforward extension of the Arrhenius law accurately describes diffusion in the thermodynamically stable liquid of hard spheres. A sharp negative deviation from this behaviour is observed as the liquid is compressed beyond its stability limit, which can be compared with super-Arrhenius slowing down in supercooled liquids. It is suggested that the observed dynamical transition is caused by extra entropy barriers arising in the free-energy landscape. This interpretation is corroborated by the observation of characteristic anomalies in the geometry of a particle trajectory.

Computer simulations of liquids using hard spheres (HS) remain a key source of information about the most fundamental aspects of the liquid state. The main intellectual attraction of these models is in their purely geometric nature: all the observable properties of a HS liquid can be reduced to a single thermodynamic quantity - entropy, a measure of finite-range structural correlations hansen ; mountain . Thus, mapping a phenomenon exhibited by liquids with continuous interactions on a HS model yields its interpretation in terms of the statistical geometry of sphere packing.

An important aspect of liquid behaviour addressed using HS models is supercooled liquid dynamics yip ; gotze ; cohen ; speedy . This notion refers to the complex of dynamical anomalies observed in a liquid which remains in a metastable equilibrium below its melting point. For the one-component HS system, the stable liquid domain is bounded by the critical value of the packing fraction $`\eta _c=0.493`$ speedy  ($`\eta =\pi \rho \sigma ^3/6`$, where $`\rho `$ is the density and $`\sigma `$ is the HS diameter). Having been compressed beyond $`\eta =0.54`$, the HS liquid inevitably crystallizes speedy .

The most prominent indicator of the supercooled liquid regime, and, in fact, its defining feature ediger ; angell , is the characteristic behaviour of the relaxation dynamics which slows down with cooling progressively faster than it can be inferred from the Arrhenius prescription universally describing the temperature variation of the transport coefficients in stable liquids.

In contrast to atomistic models with continuous interaction, the HS systems lack the energy scale. Therefore, in order to compare the non-Arrhenius liquid dynamics as observed in the former with the respective behaviour of the metastable HS liquid, we have to resolve the fundamental problem of finding a description of the stable HS liquid dynamics that would be an adequate analog of the Arrhenius law. In this capacity, Hildebrant-Batchinski relation batch  connecting the diffusion coefficient with the free volume was suggested woodcock . The HS liquid diffusion woodcock  follows this relation in the stable domain and deviates from it for $`\eta >\eta _c`$ angwood . However, this deviation is positive and, therefore, cannot be regarded as an appropriate analog of the super-Arrhenius slowing down in supercooled liquids.

In this Letter, it is demonstrated that an earlier suggested relation between the diffusion coefficient and the thermodynamic entropy nat  can be interpreted as a straightforward extension of the Arrhenius law. Having been tested by the molecular dynamics, it is found to accurately describe the diffusion in the stable one-component HS liquid. Moreover, a significant slowing down as compared with this relation is observed in the metastable HS liquid domain, a direct analogue of the super-Arrhenius slowing down in supercooled liquids. This transition in the HS dynamics is discussed in terms of the free energy landscape.

The Arrhenius law asserts that the diffusion coefficient, $`D`$, scales with the temperature $`T`$ as

$$D=D_0e^{A/k_BT}$$
(1)
where $`A`$, the so-called activation energy, is interpreted as the average height of the free-energy barrier, per particle, that has to be crossed in order to perform an elementary step in the process of diffusion. It can be presented as $`A=E_ak_BTS_a`$, where $`E_a`$ and $`S_a`$ represent, respectively, the heights of the energy barrier and the entropy barrier involved. For the HS liquid, $`A/k_BT=S_a`$; $`e^{S_a}`$ can be understood as the inverse average number of attempts per a successful diffusive transition of a particle. Two possible types of the free-energy landscape can be discussed in the context:

(i) The available configurations are abundantly connected; the connectivity is facilitated by independent motions of individual particles, and each particle has an unrestricted immediate access to the entire space volume assigned to it by the structural constraints. Then the inverse average number of attempts needed for a particle to move from its current position scales as $`e^{A/k_BT}=e^s`$ where $`s`$ is the excess entropy, per particle, - the difference between the system’s entropy and that of the perfect gas at the same thermodynamic conditions.

(ii) The configuration-space connectivity is restricted by additional constraints, complementary to the ensemble-averaged structural constraints dominating type (i) landscape. In this case, the diffusive motions of an individual particle are coupled, in a presumably hierarchical manner, with the respective motions of other particles within a certain range pwa . These dynamical constraints give rise to extra high entropy barriers which divide the configuration space into a set of components palmer . These barriers dominate the relaxation dynamics and, in this way, control the diffusion rate. In terms of Eq. (1) this means that $`e^{A/k_BT}<e^s`$ It has to be emphasized that this type of free-energy landscape exists on a limited time-scale and, therefore, is quite distinct from the equilibrium distribution. However, there exists a conjecture ag  expressing the cooperativity range in terms of equilibrium entropy; this conjecture is not discussed here.

The described landscape regimes can be understood by keeping in mind an important topological peculiarity of the multidimensional configuration space ma  where each point is adjacent to almost all other points. This can be illustrated by the “small world” topology of abundantly connected multidimensional networks sw . On the other hand, if an entropy barrier characteristic of type (ii) landscape arises, it will be immediately close to almost all points of the configuration space with a profound effect on the general rate of dynamics.

The central postulate we adopt here is that the free-energy landscape that controls the relaxation dynamics in the stable HS liquid is of type (i). Having assumed this, Arrhenius law (1) can be transformed into the following simple form relating $`s`$ and $`D`$:

$$D=D_0e^s$$
(2)
If the time-scale is measured in terms of the Enskog collision frequency $`\mathrm{\Gamma }`$ chap , this relation becomes equivalent to the earlier suggested scaling law for atomic diffusion nat  which was successfully tested on a number of simple liquids, including the HS liquid using a two-body approximation for $`s`$ mountain . In the full $`s`$ version, the scaling law was tested on some metallic liquids babak , but for the HS liquid, such tests have not been done so far. In this context, a conjecture should be mentioned where $`De^{Bs}`$ with $`B1`$ rosen .

In this study, relation (2) was tested using a molecular-dynamics simulation of a HS system comprising 6912 identical particles. This simulation was necessitated by the fact that, although the HS diffusion has been explored in a number of simulations, a considerable uncertainty remains concerning its size-dependence at different densities speedy ; erpenbeck .

The assumption that the dynamics in a liquid is controlled by the excess entropy, a measure of the structural correlations, implies that the length should be measured in terms of the characteristic length of the structure which manifests itself in the position of the main peak of the structure factor. Fig. 1 shows that the latter scales with the density $`\rho `$ as $`\rho ^{1/3}`$; therefore, $`\rho ^{1/3}`$ is a convenient unit of length. Fig. 1 also shows that this unit of length is consistent with the equilibrium separation of the nearest neighbours in the Lennard-Jones liquid of which the structure factor can be reproduced by that of the HS liquid. The use of $`\sigma `$ as the unit of length that was adopted in the earlier test of the scaling law nat  on the HS liquid is apparently inappropriate at lower densities.

The HS excess entropy was calculated here by the thermodynamic integration using Boublik-Nezbeda approximation boublik  for the equation of state. A minor correction was introduced for $`\eta >0.47`$ speedy .

The simulation results are shown in Fig. 2. Each value of $`D`$ was obtained from the mean-square displacement averaged over $`10^4`$ collisions per particle. These results are compared with another set of the HS diffusion data which represent an extrapolation to the infinite system sizeerpenbeck ; alder . It is clear that with $`D_0=0.079`$, relation (2) accurately describes the diffusion in the one-component HS system within the range $`0.36<\eta <0.49`$ which covers the entire stable liquid domain.

$`D_0`$ corresponds to the diffusion rate in a hypothetical situation where the discussed liquid diffusion mechanism is still valid but $`s=0`$. The observed value of $`D_0`$ can be rationalized using the following arguments based on a simple toy model in the spirit of the lattice gas models. Consider a tessellation of the space into equal-size cells. Each cell contains one stationary particle which stays within the cell until its status is changed; besides, there is a number of mobile particles distributed in space. At each time-step, a mobile particle collides with the stationary particle that is in the same cell. As a result, the former becomes stationary, whereas the latter becomes mobile and traverse the cell boundary into a randomly chosen adjacent cell (if several mobile particles simultaneously occur in one cell, an arbitrary collision sequence is assumed). Thus, each particle performs a random walk with the step size $`\rho ^{1/3}`$. Assuming that the number of mobile particles is small, and keeping in mind that there is one particle jump per two collisions, the diffusion coefficient can be estimated as $`D\mathrm{\Gamma }\rho ^{2/3}/12`$, which is in a good agreement with the above value of $`D_0`$.

The results presented in fig. 2 demonstrate that in the metastable liquid domain $`\eta >\eta _c`$, the HS diffusion exhibits a rapidly increasing slowing down as compared with the Arrhenius-like behaviour conjectured by relation (2). Apparently, this behaviour can be regarded as a direct analog of the super-Arrhenius slowing down in the supercooled liquids with continuous interactions. Using the above arguments, the onset of this new dynamical regime can be interpreted as a clear signal that the height of entropy barriers exceeds the ensemble-averaged excess entropy, and the free-energy landscape that controls the HS liquid dynamics transforms from type (i) to type (ii).

The long-time configuration-space constraints conjectured to dominate the HS liquid dynamics in the metastable domain can be detected from their impact on a particle trajectory. The latter represents a 3D real space projection of the multidimensional trajectory of the system in its configurational space. The relaxation dynamics that unfolds in type (i) landscape represents a random walk in the configurational space, which must also be true for a single particle trajectory. The long-time constraints caused by the development of extra barriers characteristic of the type (ii) landscape give rise to distinct anomalies in the geometry of a particle trajectory. A quantity that can be conveniently employed as an indicator of a possible deviation from the random walk geometry in a particle trajectory is the maximum absolute value of the displacement from the original position within the time interval $`(0,t)`$:

$$Z(t)max\{R(\tau ),0<\tau <t\}$$
(3)
The ratio of its second moment $`Z^2(t)`$, the mean square maximum displacement, to the mean square displacement $`R^2(t)`$ for a random walk in 3D space is equal to $`1.4024\mathrm{}`$ seshadri . The evolution of this ratio as a function of $`R^2(t)`$ for various values of $`\eta `$ is presented in Fig. 3. For $`\eta =0.48`$, in the stable liquid domain, the random-walk limit is attained after a short initial period of ballistic behaviour, as soon as a particle leaves the cage of its nearest neighbours. The geometry of a particle trajectory changes significantly as the HS liquid is compressed beyond its stability limit: a new regime of diffusion develops where $`Z^2(t)/R^2(t)`$ exceeds the random walk value. This indicates that a particle trajectory has a higher chance than random to return to an earlier covered region. The extent of this apparent confinement effect in a particle trajectory amounts to the diffusive displacement of several particle diameters, and it increases rapidly, as well as the magnitude of the effect, with the increase of $`\eta `$.

According to the arguments presented above, the height of an entropy barrier can be regarded as a measure of the real-space extent of the local configurational transformation that the system has to perform to cross the barrier. Extra high barriers in type (ii) landscape correspond to extended correlations, both positional and dynamical. The rapid increase in the (time-limited) correlation length in the metastable HS liquid that has been detected here from the geometry of a particle trajectory represents a generic feature of the supercooled liquid behaviour fischer . Ergodicity restoring relaxation whereby the extra high entropy barriers are created and destroyed is facilitated by highly collective activated hopping. This dynamics was discerned in the metastable HS liquid using dynamical density-functional theory analysis kawasaki .

The profound similarity between the pattern of slowing down in the metastable HS liquid and that in conventional supercooled liquids indicates the entropic nature of this phenomenon ritort , and suggests a unifying scenario for its development in terms of the free-energy landscape transformation. The change in the packing geometry under cooling/compression beyond the liquid phase stability limit, presumably caused by reduction of the free volume, leads to formation of long-lived extra high entropy barriers and the loss of the short time-scale configuration-space connectivity. This results in changing the relation between between $`s`$ and $`D`$ which can be detected as a deviation from the scaling law (2).

This study was supported by the Swedish Research Council.


format: text

text:
# Sum rules and electrodynamics of high-𝑇_𝑐 cuprates in the pseudogap state.

## Abstract

We explore connections between the electronic density of states (DOS) in a conducting system and the frequency dependence of the scattering rate $`1/\tau (\omega )`$ inferred from infrared spectroscopy. We show that changes in the DOS upon the development of energy gaps can be reliably tracked through the examination of the $`1/\tau (\omega )`$ spectra using the sum rules discussed in the text. Applying this analysis to the charge dynamics in high-$`T_c`$ cuprates we found radically different trends in the evolution of the DOS in the pseudogap state and in the superconducting state.

One of the most enigmatic properties of cuprate high-temperature superconductors is the pseudogap in the spectrum of low-energy excitations developing primarily in weakly doped materials at a temperature T\* well above the superconducting transition $`T_c`$. First discovered through experiments probing spin-related behavior the pseudogap also leads to a characteristic modification of a variety of properties in the charge sector. Some aspects of the pseudogap have lead to an interpretation of this phenomenon in terms of a ”precursor” of the superconducting gap whereas another school of thought is proposing scenarios based on spin- or charge-density wave fluctuations. Using the novel sum rule analysis of the optical constants of cuprates measured for E vector along the CuO2 planes we are able to put constraints on the possible microscopic origins of the pseudogap.

It is generally agreed that the electromagnetic response of the conducting CuO2 planes is dominated by the nodal segments of the Fermi surface which remain nearly unaffected by the pseudogap developing close to the $`(\pi ;0)`$ and $`(0;\pi )`$ regions. Nevertheless, there are profound consequences of the pseudogap state for the in-plane carrier dynamics in many cuprates. The spectroscopic signatures of the pseudogap are best resolved in the spectra of the frequency dependent scattering rate:

$$\frac{1}{\tau (\omega )}=\frac{\omega _p^2}{4\pi }Re\left(\frac{1}{\sigma (\omega )}\right),$$
(1)
where $`\omega _p`$ is the plasma frequency and $`\sigma (\omega )=\sigma _1(\omega )+i\sigma _2(\omega )`$ is the complex conductivity. In the underdoped compounds (left panels in Fig. 1) the spectra of $`1/\tau (\omega )`$ are depressed at $`T<T^{}`$ and reveal a characteristic threshold structure at $`\mathrm{\Theta }500600`$ cm-1. At $`\omega <\mathrm{\Theta }`$ the scattering rate measured at $`T<T^{}`$ increases faster than linear with a crossover to a less rapid (nearly linear) increase at frequencies above the $`\mathrm{\Theta }`$ threshold. This behavior is universal and is observed in a large variety of materials.  As doping progresses to optimal and overdoped similar depression of $`1/\tau (\omega )`$ occurs only at $`TT_c`$(right panels in Fig. 1).

Conventional interpretation of these results is based on phase space arguments: a partial gap reduces the amount of final states available for scattering of the nodal quasiparticles therefore leading to a depression of $`1/\tau (\omega )`$ for energies below the pseudogap. The above argument is corroborated by a close correspondence between the spectra of the scattering rate extracted from IR optics and the form of the electronic self-energy obtained from the angular resolved photoemission spectroscopy (ARPES). In this paper we propose a new analysis suggesting a direct relationship between the behavior of $`1/\tau (\omega )`$ and the prominent features of the electronic density of states (DOS) $`D(\omega )`$. To elaborate on a connection between the spectra of $`1/\tau (\omega )`$ (Eq. 1) and the essential characteristics of the DOS in the studied materials we first turn to the BCS model. We calculated the conductivity of an $`s`$-wave superconductor both in the clean and dirty limits and determined $`1/\tau (\omega )`$ from Eq.1. A striking result is that the form of $`1/\tau (\omega )`$ reproduces the key features of $`D(\omega )`$ in a BCS superconductor including the gap at $`2\mathrm{\Delta }`$ followed by a sharp peak (Fig. 2). The divergence of the actual DOS at the gap edge is not found in $`1/\tau (\omega )`$. Also, the peak in $`1/\tau (\omega )`$ may occur somewhat above the gap energy and its location depends on the strength of (impurity) scattering. Nevertheless, the overall character of the single-particle density of states is reproduced in the $`1/\tau (\omega )`$ spectra. Similar connection also holds for $`d`$-wave superconductors.

An intriguing attribute of the model spectra displayed in Fig. 2 is revealed by the integration of $`1/\tau (\omega )`$. It appears that the area removed from the intra-gap region is fully recovered due to the overshoot of the spectra at $`\omega >2\mathrm{\Delta }`$. The reason for the above ”sum rule” can be clarified by the following considerations. Expressing $`\sigma (\omega )`$ Eq. 1 through the dielectric function $`ϵ_1(\omega )+iϵ_2(\omega )`$ one finds $`1/\tau (\omega )=[\omega _p^2/\omega ][ϵ_2/(ϵ_1^2+ϵ_2^2+12ϵ_1)]`$. The term on the RHS in the latter expression is to be compared with the integrand in the well-known sum rule:

$$_0^{\mathrm{}}𝑑\omega \frac{1}{\omega }\frac{ϵ_2}{ϵ_1^2+ϵ_2^2}=\frac{\pi }{2}$$
(2)
The difference between the two terms arising from $`12ϵ_1`$ in the denominator of the $`1/\tau (\omega )`$ expression results in less than 1$`\%`$ correction at frequencies $`\omega <0.5\omega _p`$. Therefore the balance between the areas associated with the intra-gap region and the overshoot in the $`1/\tau (\omega )`$ spectra that can be expressed as:

$$𝑑\omega [1/\tau ^A(\omega )1/\tau ^B(\omega )]0$$
(3)
is in fact expected from the sum rule arguments. In Eq. 3 indexes $`A`$ and $`B`$ refer to different states of the studied system (e.g.: normal, pseudogap, superconducting). If the $`\frac{\pi }{2}`$ sum rule (Eq. 2) and Eq.3 are applied to the results plotted in Fig. 2 this analysis provides an additional support for the idea of the direct correspondence between the $`1/\tau (\omega )`$ spectra and the $`D(\omega )`$. Indeed, the BCS model suggests that the states removed from the intra-gap region are recovered at energies above the gap. This fundamental conservation is also manifested in the behavior of $`1/\tau (\omega )`$ (Fig. 2). Moreover, in those situations when $`[12ϵ_1][ϵ_1^2+ϵ_2^2]`$ the DOS conservation can be quantitatively verified using Eq. 3.

In Fig. 3 we show an experimental example confirming connections between the structure seen in $`1/\tau (\omega )`$ spectra and the features of the DOS. We studied the response of the single crystals of Cr which is a spin- density wave (SDW) antiferromagnet with the Neel temperature $`T_N=312`$ K. The results of the scattering rate analysis are reported in Fig. 3 for the first time. We find that at $`T>T_N`$ the absolute value of $`1/\tau (\omega )`$ increases as $`\omega ^2`$ in accord with the Fermi-liquid (FL) theory. At 10 K the SDW gap is fully developed giving rise to a non- trivial form of the $`1/\tau (\omega )`$ spectra. The scattering rate is suppressed at $`\omega <500`$ cm-1 but then overshoots the 320 K spectrum with a maximum at $`\omega 900`$ cm-1. This behavior is similar to the results for the optimally doped cuprates plotted in Fig. 1. The main difference is that cuprates show a linear ”background” in $`1/\tau (\omega )`$ as opposed to the $`\omega ^2`$ background seen in the data for Cr. It is appropriate to compare the overshoot in the $`1/\tau (\omega )`$ spectra with the calculations for a BCS superconductor since the BCS theory is believed to produce an accurate representation of the DOS in a SDW system. While unmistakable similarities are revealed by such a comparison, it is important to keep in mind that in Cr only a part of the Fermi surface is affected by the SDW state. Therefore the gap in the DOS is incomplete which may account for a more gradual increase of $`1/\tau (\omega )`$ in the vicinity of the gap energy compared to theoretical prescriptions in Fig. 2. As pointed out above, the conservation of the DOS may be reflected through the balance of the ”intra-gap” and ”overshoot” areas (Eq. 3) provided the $`12ϵ_1`$ correction is insignificant. We found that in Cr the latter correction amounts to about 0.1 $`\%`$ of the $`ϵ_1^2+ϵ_2^2`$ at $`\omega <10,000`$ cm-1. Thus it is hardly surprising that Eq.3 is fulfilled for Cr with the accuracy of 10$`\%`$ with integration limited to $`\omega _c=1500`$ cm-1. Similar behavior of $`1/\tau (\omega )`$ spectra is also found in charge-density wave materials as well as in the heavy fermion materials revealing hybridization gap.

It is apparent from Fig.2 and Fig. 3 that the spectra of $`1/\tau (\omega )`$ capture the gross characteristics of the density of states in a conducting system especially in those situations when the DOS is (partially) gapped. We now return to the data for cuprates focusing on the implications of Eqs. 2-3 for the understanding of the charge response across the phase diagram. As pointed out above, the dominant features of the $`1/\tau (\omega )`$ spectra obtained for underdoped crystals in the pseudogap state are similar to those seen at $`TT_c`$ in the optimally doped samples. On a closer examination one finds differences as well. Suppression of $`1/\tau (\omega )`$ at $`\omega <\mathrm{\Theta }`$ in the pseudogap state does not lead to the development of the overshoot between the spectra measured at $`TT^{}`$ and $`T>T^{}`$. The data for the optimally doped crystals taken at $`TT_c`$ and $`TT_c`$ does reveal an overshoot at $`\omega 800`$ cm-1. In the latter case, the area under the overshoot is balanced out by the reduction of $`1/\tau (\omega )`$ at low energies with the accuracy better than 10 $`\%`$ in Tl2Ba2CuO6 crystal and about 15 $`\%`$ in the YBa2Cu3O6.95 sample. This balance is expected from Eqs. 2-3 because the contribution of $`12ϵ_1`$ term is about 0.1$`\%`$ of $`[ϵ_1^2+ϵ_2^2]`$ for either of the above superconductors. Thus the form of the $`1/\tau (\omega )`$ spectra for the optimally doped system is in accord with the notion of the transfer of the states from the intra-gap region to a peak at $`\omega >\mathrm{\Theta }`$. This behavior can be naturally attributed to the opening of the superconducting energy gap. Although, the absolute values of the optical constants for the underdoped compounds are also in the regime when the Eq. 3 ought to be satisfied, data gives no indications even for a partial recovery of the area associated with the far IR depression of $`1/\tau (\omega )`$ is the pseudogap state.

The fact that in cuprates the spectra taken above and below $`T^{}`$ ”violate” Eqs. 2-3 within the frequency range of Fig. 1 signifies the differences in the changes of the low energy DOS in the pseudogap state from the changes associated with superconductivity. Essentially, the analysis described above indicates that the states removed from $`\omega <\mathrm{\Theta }`$ at $`T<T^{}`$ disappear from the energy interval reliably sampled in our experiments (0.5 eV). This result is in accord with the tunneling studies of the temperature dependence of the electronic DOS in underdoped cuprates. Indeed, data obtained at $`TT_c`$ reveals a transfer of the electronic states from the intragap region in the quasiparticle peaks just above the gap. On the contrary, the formation of the pseudogap only leads to a depletion of DOS without visible traces of recovery of $`D(\omega )`$ at least at $`\mathrm{}\omega <`$ 0.1 eV. Specific heat measurements are also in accord with our findings since these experiments suggest a reduction of the DOS in the pseudogap state within the energy window accessible to thermodynamic probes( $`500`$ K). Yet, another indication of an energy scale extending beyond 0.1-0.5 eV that is involved into the pseudogap state response is provided by the IR studies of the interlayer $`c`$-axis conductivity $`\sigma _c(\omega )`$. Non-conservation of the low-energy spectral weight $`N_{eff}(\omega )=_0^\omega 𝑑\omega ^{}\sigma _c(\omega ^{})`$ for $`\omega <0.20.4`$ eV is inferred from the oscillator strength sum rule; these $`c`$-axis results parallel our findings deduced from the examination of the in-plane scattering rate using Eqs. 2-3. Thus one can conclude that both spectroscopic and thermodynamic studies of cuprates reveal anomalous changes of $`D(\omega )`$ in the pseudogap state distinct from the DOS effects associated with superconductivity. Hence, these results argue against a common origin of the pseudogap state and superconducting state.

Further support for the distinct genesis of the pseudogap and of superconducting gap is provided by the data for La2-xSrxCuO4 (La214) and Nd2-xCexCuO4 (Nd214) materials. The scattering rate analysis of the in-plane response of these systems reveals a characteristic threshold structure at $`\mathrm{\Theta }`$ 500-600 cm-1. In the double- or triple-layered materials the development of superconducting gap occurs approximately at the same energy as the pseudo-gap structure. However, in La214 and Nd214 compounds the signatures attributable to superconductivity are confined to energies below 40-50 cm-1 whereas the pseudogap feature is essentially identical to what is seen in the double- layered compounds with $`\mathrm{\Theta }500600`$ cm-1. It is difficult to account for the difference by more than one order of magnitude between the energy scales associated with the pseudogap and with superconductivity if both effects result from the same microscopics. The formation of the pseudogap is commonly discussed in the context of the density wave ideas. Our analysis in conjunction with the data for a canonical SDW system such as Cr argues against this point of view as well. Data shown in Fig. 3 is a direct testimony to a conventional transfer of the intragap states to the energy region right above the gap. Similar overshoot is also found in the quasiparticle relaxation of a charge-density wave system. Therefore, both spin and charge density waves are unable to account for the experimental situation in the pseudogap state in high-$`T_c`$ cuprates.

In conclusion, the sum rule analysis of the optical constants based on Eqs.2-3 argues against the decisive role of both superconducting and density wave fluctuations in the formation of the pseudogap. We emphasize that the scale associated with the spectroscopic signatures of the pseudogap dramatically exceeds the characteristic temperature $`T^{}`$. Such a disproportion between $`\omega `$\- and $`T`$-scales may signal relevance of the many-body effects to the formation of the pseudogap state. Authors are grateful A.V.Chubukov, J.P.Carbotte, D.Pines, T.Timusk, and especially to D.van der Marel for useful discussions. This work is supported by the DoE.


text:
# Infrared Emission from the Radio Supernebula in NGC 5253: A Proto-Globular Cluster?

## 1 Introduction

The dwarf I0/S0 galaxy NGC 5253 is actively forming stars, with abundant $`H\alpha `$ emission, WR stars, and many optically visible star clusters within the central 100 pc (Beck et al., 1996; Gorjian, 1996; Calzetti et al., 1997). More than half the radio emission from this galaxy is found to come from an intense source of diameter 1.3 pc at 2 and 1.3 cm (for a distance of 4 Mpc), in the star formation region but not identified with any optical counterpart (Turner et al., 1998). Turner et al. (2000) called this source the radio supernebula; it has a brightness temperature of $`10^4\mathrm{K}`$ at 2 cm and an optically thick spectrum that is either flat or slightly rising. Extra-galactic radio sources of such high $`T_b`$ and luminosity are usually non-thermal synchrotron sources with steeply falling spectra or AGN with flat spectra. Thermal radio sources, generally gaseous nebulae surrounding young stars (HII regions) have almost flat spectra, but are typically two to three orders of magnitude less luminous than this radio source. If this were an HII region it would require an extraordinarily high number of stars to excite the gas. If it were a non-thermal source it would also be extraordinary for its kind. From the radio emission alone, we cannot tell if this is an HII region, a supernova remnant (SNR), or an Active Galactic Nucleus (AGN).

We have therefore obtained high resolution images of the radio source in the mid-infrared, since an HII region will be a strong infrared emitter with a characteristic spectral energy distribution, different from a SNR or AGN. The observations and results are discussed in the next section, followed by our deductions as to the nature of the source and its properties.

## 2 Observations

The observations were made during the re-commissioning of the Long Wavelength Spectrometer (LWS) (Jones & Puetter, 1993) on the Keck I Telescope on Mauna Kea in March 1999. Two images were obtained: one with a 1 $`\mu `$m wide filter centered at 11.7 $`\mu `$m and one with a 0.5 $`\mu `$m wide filter centered at 18.7 $`\mu `$m. Fluxes were calibrated using the standard stars $`\beta `$ Gem and $`\alpha `$ Boo and are estimated to have uncertainties of 10%. In Figure 1 we show the 11.7 $`\mu `$m image in color superimposed on an optical photograph (Calzetti et al., 1997) from the Hubble Space Telescope (HST), and in Figure 2 the individual 11.7 and 18.7 $`\mu `$m images. The registration uncertainty is limited by the HST image and is estimated to be $``$1″.

The infrared source is barely resolved at a diameter of 0.58″ (FWHM) at 11.7 $`\mu `$m which is about twice the diffraction limit of the Keck. The source is 0.62″, only 25% larger than the diffraction limit, at 18.7$`\mu `$m. These sizes may be upper limits to the true source size since at the time of the observations the telescope settle times had not been accurately determined and inadequate settle times may have broadened the images. There is no significant difference between the appearances of the source at the two wavelengths, and allowing for the differences in beam size the radio and infrared images agree. The source has 2.2 Jy total flux at 11.7$`\mu `$m with S/N of 17 at the peak pixel and 2.9 Jy at 18.7$`\mu `$m, S/N of 12 at the peak. These agree with ground-based small-aperture 10.8$`\mu `$m fluxes ($``$ 1.4 Jy, Lebofsky & Rieke, 1979; Telesco, Dressel, & Wolstencroft, 1993) especially when we take into account that the source is brighter at longer wavelengths. While the ISO-SWS spectrum shows emission and absorption features within the 11.7$`\mu `$m band, our continuum flux is in agreement with the overall mid-infrared continuum (Crowther et al., 1999). IRAS fluxes at 12 and 25 $`\mu `$m are 2.6 and 12 Jy, respectively. The IRAS satellite beam included the entire galaxy; the fact that the IRAS fluxes are so close to what we measure indicates that this is a very bright and a very small infrared source.

## 3 The Infrared Source

### 3.1 What is It?

There are, as we stated above, three possible sources of the rising-spectrum radio emission from NGC 5253: an HII region, a SNR, or an AGN. The discovery that the radio source is also a strong mid-infrared emitter narrows down the possibilities. First, it rules out the presence of a bright radio supernova; supernovae are not significant mid-infrared sources. Most Galactic SNR are not detected at all in the IRAS bands (Arendt, 1989). A supernova explosion inside a dense molecular cloud will be very bright in the infrared (Shull, 1980), but that is calculated to be a short-lived stage, with the emission declining significantly on a time scale of a few years (Weiler et al., 1986). The agreement of our measurement with the IRAS fluxes from 1983 and the ground-based work from 1972 and 1986 (Rieke & Low, 1972; Lebofsky & Rieke, 1979; Telesco, Dressel, & Wolstencroft, 1993) argues against this model.

The possibility that the source is an AGN cannot be excluded so decisively. AGN are not themselves expected to be strong infrared emitters, although their spectra will depend on the amount and distribution of the dust nearby. Much of the observed range in infrared flux is probably due to the presence of both non-stellar central engines and starbursts in the same galaxy, and often in the same nuclear region. NGC 1068 is a nearby example of a galaxy with both an AGN and a powerful starburst in close quarters; at larger distances, such sources can be difficult if not impossible to separate (Heckman et al., 1995). The high spatial resolution of our observations and the closeness of NGC 5253 will however challenge any model which tries to combine an AGN and a starburst in the radio-infrared source. The conclusive test for the presence of an AGN is spectroscopic, not photometric, and depends on the detection of fast moving gas near the central engine. That the gas studied optically in NGC 5253 does not have broad lines is not relevant to the gas in the radio source, which must be heavily obscured as it is 70 times brighter than one would predict based on the $`H\alpha `$ line (Calzetti et al., 1997). The only spectroscopic information currently available on the gas inside the radio and infrared source is the $`10.5\mu `$m \[SIV\] line (Beck et al., 1996), which is no more than 60 $`\mathrm{km}\mathrm{s}^1`$ wide. Those observations cannot completely rule out a very weak plateau, as from a weak AGN. Another regime in which the influence of an AGN is usually strong is the X-ray: the X-ray spectrum and luminosity of NGC 5253 is perfectly consistent with a starburst in a clumpy medium (Martin & Kennicutt, 1995). In short, while a small and quite unusual AGN could perhaps be hidden in NGC 5253, there is nothing in the data that needs or suggests one.

The final possibility, which the data strongly favor, is that the source of the infrared and radio emission is an HII region ionized by thousands of young stars. We presented this model in Turner et al. (2000) based on radio data alone. The high mid-infrared fluxes of the NGC 5253 source, which exclude SNR and go far to exclude AGN, are perfectly consistent with the spectral energy distribution of an optically obscured HII region. The ratio of $`11.7\mu `$m flux to the thermal radio flux is 150, a factor of 5 higher than would be expected from pure $`Ly\alpha `$ heating (Ho et al., 1990) and in agreement with results from many other star formation regions, Galactic and extra-Galactic. The observed ratio in starburst galaxies is usually in the 120-250 range.

The model that the NGC 5253 source is an optically obscured HII region fits all the data at hand: the radio spectrum, which is optically thick at 6 and probably 2 cm (Turner et al., 1998, 2000), as are dense and compact HII regions, the infrared flux and its ratio to the radio emission, the infrared lines which agree with the radio continuum (Beck et al., 1996), and the great excess of infrared and radio strengths over that predicted from the optical. The NGC 5253 source could be a compact HII region such as W51 or K3-50 writ large, except that Galactic compact HII regions typically contain 1 bright star in 0.01 pc diameter; the NGC 5253 source must hold not one OB star but a populous cluster. The ionization requirements of the super nebula are immense, $`4\times 10^{52}sec^1`$. The stellar population and density which are calculated in Turner et al. (2000) for a range of IMFs are accordingly high: $`2\times 10^5`$ to $`2\times 10^6`$ stars and $`5\times 10^5`$ to $`10^6`$ $`M_{}`$ within a 1-2 pc region. The source requires not just a cluster of stars, but a large cluster with the stellar density, the size, and luminosity of a bright globular cluster.

### 3.2 On the Spectral Dominance and Spatial Concentration of the Source

The NGC 5253 super nebula is not only remarkable in itself; it dominates the galaxy’s total spectrum to an astonishing degree. 80-90% of the total galactic flux in the 12 $`\mu `$m bands comes from this 1-2 pc region. NGC 5253 as a whole is rather blue in the infrared: it is stronger at 60$`\mu `$m than at 100$`\mu `$m and its $`\nu F_\nu `$ peaks at 25$`\mu `$m; a significant fraction of the flux at 25$`\mu `$m may be due to a silicate emission feature (Crowther et al., 1999). The supernebula has been imaged in only two bands, but also appears extremely blue, with a deduced dust temperature (for $`\beta =1.5`$) of 180 K, consistent with the ISO mid-infrared spectrum (Crowther et al., 1999). We do not know what the supernebula contributes at 60 and 100$`\mu `$m, but even in the most extreme case that it contributes nothing, it would still provide $``$20%, or $`4\times 10^8\mathrm{L}_{}`$, of the total fluxes at wavelengths below 100$`\mu `$m, which constitutes at least 75% of the total bolometric luminosity ($`2.4\times 10^9\mathrm{L}_{}`$), radio to X-ray, of the galaxy. We note that the population of O stars required for the radio emission from the supernebula (Turner et al., 2000) will have a total luminosity of $`0.81.2\times 10^9L_{}`$, compared to $`1.7\times 10^9L_{}`$, the total observed infrared luminosity. This argues that all the infrared is distributed at the shorter wavelengths and that the super nebula gives at least half of the total IR luminosity, and potentially as much as two-thirds, in an infrared-dominated galaxy. So this giant HII region, less than 2 pc in diameter, may well dominate its host galaxy as thoroughly as does a bright AGN. NGC 5253 is the best example so far of the trend (Ho et al., 1990) that the bulk of the infrared emission in starburst galaxies is usually found in a much smaller volume than the IRAS beams, and implied by the finding (Bryant & Scoville, 1999; Downes & Solomon, 1998) that the molecular clouds which fuel star formation are concentrated in small active regions.

## 4 A Very Young Star Cluster and its Problems

If, as we have argued, the radio-infrared source is a super nebula or giant HII region ionized by thousands of young O stars, conditions in the source must be extreme. If we adopt a total population of $`10^6`$ stars (Turner et al., 1998, 2000) and a diameter of 1.5 pc, then the stellar separation is $``$0.02 pc. These cluster conditions, while unusual, are actually what would be expected from a just-formed or still forming super star cluster or globular cluster. At present the nebular density is high, $`4\times 10^4\mathrm{cm}^3`$, as deduced from the high radio optical depth (Turner et al., 2000) and there is enough dust in the region to obscure the optical emission. As the stellar activity continues the gas and dust may be expected to disperse, weakening the radio and infrared emission and permitting the cluster to be seen in the visible, which will then resemble super star clusters which are commonly seen in starburst regions.

We argue that the source is excited by a very young super star cluster, large enough to become a globular cluster; just how young? Since we do not see the stars we cannot find their ages directly. The age of the nebula is limited from its dynamical lifetime, which must be short. The gas in the nebula, which is at $`10^4`$ K, will be overpressured relative to the cooler, less dense ISM and will expand at the sound speed, $`10\mathrm{km}\mathrm{s}^1`$, so the observed supernebula should double its size in $`10^5`$ years.

The problem of lifetimes, reviewed in Kurtz et al. (2000), haunts the study of the Galactic ultra-compact HII regions that these young star clusters so resemble. Galactic ultra-compact HII regions are thought to be older than their dynamical ages, since the dense and clumpy molecular cloud environment of these sources slows their expansions and extends their lifetimes. In NGC 5253 there is no evidence for confinement by molecular gas: in fact, no CO is detected within 200 pc of the nebula (Turner, Beck, & Hurt, 1997; Meier, Turner, & Beck, 2001). Even at the low metallicity of NGC 5253, CO should be easily detectable based on $`\mathrm{L}_{\mathrm{IR}}`$. Another way to confine the nebula is with the $`B^2/8\pi `$ pressure of the magnetic field, which will become comparable to the thermal pressure of the NGC 5253 nebulae for B fields of a few $`\mu `$G. The B field strength in selected locations of NGC 5253 was estimated from minimum-energy arguments to be 50-57 $`\mu `$G (Turner et al., 1998), high for a galaxy but typical of a dense star-forming molecular cloud, and large enough that it may be important in the pressure balance of the nebula. Finally, the cluster is so massive and compact that it is not inconceivable that the nebula is actually gravitationally bound, in spite of the large pressure gradient.

The question of the lifetimes of the earliest and most obscured stages of star formation will become increasingly important as more of these very young sources are observed on scales from the Galactic Ultra-compact HII regions to the proto-globular cluster in NGC 5253. Luminous obscured infrared bright sources have been seen in the Antennae, although larger (50 pc) than the supernebula (Mirabel et al., 1998); and compact luminous IR sources similar to the supernebula have been seen in NGC 253 (Pina et al., 1992; Keto et al., 1993). Both of these galaxies also have large numbers of young super star clusters (Watson et al., 1996; Whitmore et al., 1999). High resolution radio and CO images of galaxies (Downes & Solomon, 1998; Kobulnicky & Johnson, 1999; Tarchi et al., 2000; Beck et al., 2000) are revealing other possible candidates, in different environments and perhaps at different stages of evolution. There actually may be many of these young sources yet to be discovered by high resolution infrared observations.

We are grateful to the staff of the Keck Observatory, and the special efforts of Randy Campbell, Fred Chaffee, and Barbara Schaefer in difficult circumstances. We also thank Peter Conti, Bruce Elmegreen, Tim Heckman, Kelsey Johnson, Ellen Zweibel, and Dick McCray for helpful discussions. JLT is supported by NSF grant AST 0071276 and VG thanks the National Research Council. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by JPL, California Institute of Technology, under contract with NASA.


text:
# Distributions of secondary muons at sea level from cosmic gamma rays below 10 TeV

## 1 Introduction

Muons detected at ground level arise mainly as decay products of charged mesons. These mesons are created abundantly in hadronic showers, most of which are caused by primary cosmic ray protons and nuclei interacting inelastically with the nuclei of the atmosphere. A small fraction of muons, however, have their first origin in photonuclear reactions induced by primary or secondary cosmic gamma rays. Several authors have performed analytical or Monte Carlo calculations of the muon flux produced in gamma showers (see references \[10-21\] of ). Most of these calculations were one dimensional and all of them referred to gamma energies much larger than 10 TeV.

Despite their relative scarcity compared to the large background of muons from hadron-generated showers, muons originating from primary gamma ray interactions are important for ground-based high statistics cosmic ray experiments which are sensitive to energies $`10`$ TeV such as MILAGRO  and GRAND . Cosmic gamma rays, unlike charged hadrons, are unaffected by the Earth’s and galactic magnetic fields and their direction points directly to the location of their source. Therefore, an accumulation of muon directions around a particular angle (and also at a particular time in the case of pulsed sources or gamma ray bursts) carries direct information about the location of the source and the neutrality of the primary causing the excess. The ability to distinguish these muons from background muons depends upon statistics, the strength and energy spectrum of the emitting source, the angular resolution of the experiment, and the degradation of the angular resolution due to all the physical processes which occur between the primary and the muon at detection level. This last point is the focus of this calculation. Knowledge of this angular resolution is important in order to determine a detection window which is large enough to include a good fraction of the gamma signal and, at the same time, small enough to minimize the large uniform background of hadronic origin.

The information about the primary gamma direction is degraded by multiple physical processes so the final measured direction of the muon is no longer collinear with the primary which originated the muon’s ancestors. In particular, the direction resolution is affected by the angle of production of the mesons in the primary interaction (and then in subsequent hadronic interactions of these mesons to create additional mesons), the decay angle of the muon relative to its parent meson, and the lengths of the charged particle paths including that of the detected muon at ground level. Along these path lengths there is Coulomb scattering and deflection in the Earth’s magnetic field for the trajectories of the mesons and the muon (usually most important for the muon which typically has the longest flight path and the lowest momentum.) In addition, these effects depend upon the energy of the particles involved and on the varying properties of the atmosphere (pressure, temperature.) The muons are produced by mesons of second, third, or further generations. The number of generations typically increases with increasing primary gamma energy and each new generation involves an additional production angle which tends to increase the angular disparity between muon and primary because of more steps in a random walk type process. An opposing effect at higher energies is that all angular distributions are more forward-peaked thus decreasing the angular disparity; the resulting effect is a combination of the two.

The purpose of the present study is to estimate the influence of the various effects on the ultimate angular resolution which can be obtained utilizing ground level muons. The complexity of the problem requires detailed analyses which can only be obtained by a full three-dimensional Monte Carlo (MC) calculation. In addition, since only a small fraction of all primary photons gives rise to an event of interest, it is desirable to resort to statistical variance reduction techniques to bring the computing task within manageable limits.

The Monte Carlo program FLUKA  provides an accurate and well-tested description of the hadronic and electromagnetic interactions with all the relevant physical effects, and a set of effective statistical tools to accelerate convergence of the results. The program has been used to calculate the geometrical properties of secondary muons at sea level arising from primary gamma cosmic rays with energies from 1 GeV to 10 TeV interacting in the Earth’s atmosphere. In this energy range the cross section for direct muon pair production is much smaller than that for hadroproduction and has been neglected. This range has been shown in  to be the interval which produces the maximum secondary muon flux at ground level for primary differential gamma spectra with spectral indices larger than 2.4 (a smaller index would yield a harder spectrum and give a maximum secondary muon flux arising from primary energies beyond those considered here.)

In this paper, only incident cosmic gamma rays at perpendicular (vertical) incidence on top of the atmosphere are considered. Ground-based experiments often choose angles close to normal for reasons of simplicity and because the detection rate is a maximum, the rate dropping rather rapidly as the angle deviates markedly from normal. For small angles from normal, the angular resolution results will probably not differ significantly from those presented in this paper.

In addition to the angular correlation of secondary muons at ground level, this study provides also additional information such as the distribution of heights at which the muons originate, the muon’s kinetic energy distribution at the ground and production level, the radial distribution of the muons, and the number of generations which preceded that of the muon.

## 2 Air shower calculations with FLUKA

The 2000 version of the FLUKA Monte Carlo code was used to simulate the electromagnetic and hadronic particle cascades induced by primary gamma rays in the atmosphere. The powerful biasing capabilities of the code allowed a simulation of the full three-dimensional shower in a single run from the top of the atmosphere down to ground level.

The physical models implemented in FLUKA which are relevant for the present study as well as various variance reduction techniques were already discussed in . There it was emphasized that these models were validated by comprehensive comparisons of FLUKA results with experimental data obtained mostly at accelerators but also in recent cosmic ray studies . Nevertheless, the simulation of an atmospheric shower and, in particular, of its muon component is a relatively new area of application of FLUKA. Additional examples showing that the code allows remarkably accurate predictions in this area are given in Section 2.1. Section 2.2 contains some details of the present gamma ray shower simulations which supplement those discussed in ; this reference contains a more complete description.

### 2.1 Muon production in air showers induced by hadronic primaries and comparison to experimental data

No data exist for the muon component of gamma ray-induced showers to which FLUKA results could be directly compared. On the other hand, data on muon flux and energy spectra have become available recently for atmospheric showers induced by cosmic ray protons and nuclei. These data were mostly obtained with spectrometers in balloon-borne experiments and cover the whole range of atmospheric depths from the top of the atmosphere down to sea level.

Fortunately, the mechanisms and models for muon production and transport in gamma ray-induced showers are largely the same as those in pure hadronic air showers. In both cases, the production of the pions which eventually decay into muons is described by the Dual Parton Model (DPM) . The only difference is that in the former case the photon is assumed to fluctuate first into a quark-antiquark state which then interacts hadronically, identical to a meson. This so-called Vector Meson Dominance model (VMD) for the hadronic photon fluctuation is a well-established concept (see  and references therein) and its application together with the DPM to hadronic interactions of photons has been proven to be very successful (see, for example, ). A good description of muon production from proton- and nuclei-induced air showers should therefore give confidence that the corresponding results obtained with FLUKA for photon-induced air showers are reliable as well. In the following a few examples are given.

A critical aspect of muon production in the atmosphere is the dependence of the total muon flux on the height above sea level. This dependence is commonly expressed as a function of depth in the atmosphere (in g/cm2 of air) which is a measure of the amount of air penetrated by the cascade. It has been measured by the CAPRICE experiment which was able to distinguish between positively and negatively charged muons. Results obtained with FLUKA for negative muons are compared with CAPRICE data  in Fig. 1. The different curves correspond to the muon flux calculated for different intervals of the muon momentum. To correspond to the average geometrical response of the CAPRICE apparatus, only muons with a polar angle of less than $`9^{}`$ with respect to the vertical (i.e., vertical muons) were scored. The simulations are based on sophisticated models and methods for the sampling of the identity, energy, and angle of the hadronic primary from the distributions of galactic cosmic rays on top of the atmosphere at the time of the measurements. More details can be found in . As can be seen from the figure, FLUKA gives a good description of the negative muon flux for all depths. Similar agreement exists for positive muons (not shown here).

Energy spectra of muons were recorded by the same experiment at ground (i.e., at an atmospheric depth of 886 g/cm2. FLUKA results for positive and negative muons are compared to these data in Fig. 2a. Again, good agreement between simulation and data is found.

Finally, the ratio of the fluxes of positively and negatively charged muons is compared to CAPRICE data  in Fig. 2b. Data and calculations are again for vertical muons at a depth of 886 g/cm2. The ratios are plotted as a function of the muon energy.

The above comparisons of FLUKA results and experimental data together with those published in earlier studies  indicate that the application of FLUKA to cosmic gamma ray showers should give reliable predictions as well.

### 2.2 The Monte Carlo simulation of cosmic gamma ray showers

In order to study the dependence of the shower properties on the primary photon energy, the showers were calculated for monoenergetic photons impinging vertically on top of the atmosphere (taken to be 80 km above sea level). In total 9 sets of simulations were performed for the following primary energies: 1, 3, 10, 30, 100, 300, 1000, 3000, and 10000 GeV. As in the earlier study  the atmosphere was approximated by 50 layers with a constant density in each layer and with layer-densities decreasing exponentially with increasing altitude. It has been verified  that this approximation does not affect the conclusions drawn in this study. The geometry is described in a right-handed orthogonal system with its origin at the intersection of the shower axis with sea-level elevation, $`z`$ points to the center of the Earth and $`x`$ and $`y`$ point North and East, respectively. Most of the results discussed throughout the paper were obtained without the effect of the magnetic field of the Earth included. However, for two primary energies (10 GeV and 1 TeV) additional simulations were performed which included the magnetic field of three different geographical locations (see Sec. 3.4).

The use of several variance reduction (biasing) techniques was essential to obtain results with reasonable statistical significance. They included leading particle biasing at each electromagnetic interaction, biasing of the photon mean-free-path with respect to photonuclear interactions, biasing of the decay length of charged mesons, and particle splitting at the boundaries of different air layers. Each of these biasing techniques alters the statistical weight of a particle in the cascade. A muon as produced and transported in the Monte Carlo simulation and which carries a certain statistical weight (hereafter called “MC muon”) contributes to distributions, yields, etc., of actual, i.e., measurable muons with a probability equal to its weight. The use of biasing techniques would not be appropriate to study fluctuations within the same shower which, however, is not the aim of this study. More details on the biasing techniques can be found in .

For each MC muon reaching sea level (detection level) the following information was recorded in a file for later analysis:

1. Information on the muon at detection level: muon charge, lateral coordinates $`(x,y)`$ with respect to shower axis, direction cosines with respect to the $`x`$ and $`y`$ axes, kinetic energy, statistical weight of the muon, and number of the event (i.e., primary cosmic ray photon) which produced this muon.
2. Information on the muon at its production vertex (i.e., meson decay vertex): lateral coordinates $`(x,y)`$ with respect to shower axis, direction cosines with respect to the $`x`$ and $`y`$ axes, height $`z`$, kinetic energy, statistical weight of the muon, and identity of the decaying meson.
3. Information on the grandparent at the parent production vertex: particle identity, lateral coordinates $`(x,y)`$ with respect to shower axis, direction cosines with respect to the $`x`$ and $`y`$ axes, height $`z`$, kinetic energy, statistical weight of the particle, and generation number.
4. Information on the photonuclear interaction vertex preceding the hadronic cascade in which the muon at detection level has been created: lateral coordinates $`(x,y)`$ with respect to shower axis, direction cosines with respect to the $`x`$ and $`y`$ axes, height $`z`$, energy, statistical weight of the photon, and generation number of the photon

The generation of a particle increases with each sampled discrete interaction (electromagnetic or hadronic), i.e., the primary photon is generation “1,” the generation of the electron and positron after the first pair production process would be “2,” etc. Delta ray production and Coulomb scattering are not considered to increase the generation number.

A summary of the simulated particles cascades for each primary gamma ray energy is given in Table 1. In the following, all results refer to the sum of positive and negative muons. The number of histories (number of primary gamma rays), $`N_\gamma `$, calculated for each primary energy decreases with energy in order to obtain uniform statistical significance on the muons which reach ground level for the different primary energies. At the lowest energy (1 GeV) the muon flux per gamma at sea level is very small so that even with a significant amount of biasing in the simulations it is difficult to reach adequate statistics. Less emphasis was therefore put on the simulation at 1 GeV. The column $`N_\mu `$ gives the total number of muons reaching detection level for the given number of primary photons $`N_\gamma `$. The predictions for the muon multiplicity at sea level per incident primary gamma ray $`N_\mu /N_\gamma `$ are listed in the last column. Similar values were already reported in , however with smaller statistical significance than those in Table 1.

## 3 Angular correlations

### 3.1 Angular distribution of muons at sea level

In this section the effect of the Earth’s magnetic field is not included; Section 3.4 has a separate discussion of this added effect which depends on the geographical location of interest. Fig. 3 summarizes the deviations of the angle of the muons which reach sea level from the primary gamma direction (vertical angle). It shows the distribution of the angle $`\mathrm{\Theta }_{xz}`$ which is the muon’s angle with respect to the $`z`$ axis projected onto the $`xz`$ plane (North-down plane). This angle is defined as $`\mathrm{\Theta }_{xz}=\mathrm{arctan}(\mathrm{cos}\theta _x/\mathrm{cos}\theta _z)`$ where $`\mathrm{cos}\theta _x`$ and $`\mathrm{cos}\theta _z`$ are the direction cosines of the detected muon. Negative projections have been reflected upon the positive values for $`\mathrm{\Theta }_{xz}`$ since positive and negative angles would be symmetric in the absence of the magnetic field. Also, the distributions in the $`yz`$ plane (East-down plane) are similar to those shown in Fig. 3.

Since the angle of the primary gamma with respect to the $`z`$ axis is zero (vertically incident primaries), any deviation of the muon direction from zero degrees is due to the various interaction and transport processes occurring in the cascade between the primary gamma ray and the muon which reaches sea level. Thus the results shown in Fig. 3 represent the correlation between the direction of the detected muon and the direction of the primary gamma ray.

All distributions in Fig. 3 are normalized to unit area to compare their shapes. The shape of the distributions for primary energies above about 30 GeV varies rather slowly. Therefore, only two histograms (100 GeV and 10 TeV) are shown for this energy range. At lower energies the shape of the distributions changes significantly as can be seen from the histograms for 3 and 10 GeV. The distributions become wider with decreasing primary energy. At low primary gamma ray energies, the muons have, on average, lower energies which cause larger angular deviations from the primary due to kinematics and scattering effects which are larger for the lower energies. This is illustrated in Fig. 4 for primary energies of 10 GeV and 10 TeV. An arbitrary cut of 2 GeV was applied to the angular distributions so that contributions from muons with energies below and above this cut could be seen separately along with the total.

Two measures of the widths of these distributions as a function of the energy of the incident primary cosmic gamma ray are given in Table 2. The width parameters are: the half-width at

half-maximum-height (HWHM) and the angular half-width containing 68% of the muons, i.e., the same number of events that a 1$`\sigma `$ cut would include if the distributions had a Gaussian shape (which they don’t have, because of the long tails). Both of these parameters provide a measure of the distribution widths which minimizes the distorting effect of a tail. In addition, the values for various cuts on the kinetic energy of the detected muons (1, 2, and 4 GeV) are presented. The cutoff energy results for primary energies of less than 10 GeV were omitted due to lack of statistics. For the HWHM calculation, the maximum height is assumed to be the value in the first bin of the corresponding histogram, i.e., the average value in the interval $`|\mathrm{\Theta }_{xz}|<0.25`$ degrees. The $`x`$ and $`y`$ rectangular coordinate system was chosen to represent the results of the muon’s angular distribution as the corrections for the additional deflection due to the Earth’s magnetic field differ in these two directions. In addition, $`(x,y)`$ is the natural coordinate system for GRAND. For some experiments a space angle is more natural. If we define $`R=\sqrt{x^2+y^2}`$, then the space angle resolution ($`\delta \mathrm{\Theta }_R`$(HWHM)) can be estimated directly from the numbers in Table 2 by multiplying the $`\delta \mathrm{\Theta }_{xz}`$(HWHM) values by $`\sqrt{2}`$ due to the symmetry of $`x`$ and $`y`$ in the absence of a magnetic field.

As can be seen in this table: (i) as the primary energy decreases below 10 GeV, the width of the angular distribution increases dramatically, (ii) as the energy rises above 300 GeV, there is only a small, gradual increase in the width of the angular distribution, and (iii) the width of the distribution narrows as the muon’s cutoff energy is raised; i.e., the angular resolution becomes better by eliminating the lower energy muons which have, on average, poorer angular resolution.

### 3.2 Correlation of muon angle with distance from the shower axis

In the preceding section, angular distributions were given for all muons regardless of their lateral distance from the shower axis. However, as the distance from the shower axis increases, the angle is systematically biased away from the shower axis. Fig. 5 shows the correlation between the projected angle $`\mathrm{\Theta }_{xz}`$ and the $`x`$ coordinate (i.e., the northward distance from the shower core) of the muon position at sea level for primary gamma ray energies of 10 GeV and 1 TeV. The angular distributions are given for different intervals in $`x`$ and each distribution is normalized to unit area. Here, negative values of $`\mathrm{\Theta }_{xz}`$ have not been reflected onto positive values and are explicitly shown. As expected, the peak of the $`\mathrm{\Theta }_{xz}`$ distributions are shifted toward positive values as $`x`$ becomes more positive and the distributions become wider as $`x`$ increases. Both effects depend on the energy of the primary gamma ray and are more pronounced at higher energies.

Hence, the angular distributions narrow significantly if they are limited to muons within a certain interval around the shower axis. Information on this correlation between the angle and the lateral distance of the muon from the shower axis allows experiments to improve the angular resolution by applying cuts to the data, e.g., in $`x`$. If the available statistics demand using all available muons, their angle at large $`x`$’s ($`y`$’s) could be corrected by the most probable angle at that $`x`$ (or $`y`$) to obtain a better estimate of the primary’s direction.

This correlation can also be used to infer the height from which the muon originated, a parameter difficult to obtain experimentally. If one extrapolates the muon direction backward (upward), the intersection point with the centroid of the shower (which can be obtained from the more numerous electrons in the same shower), is an estimator for the height where the muon was created. However, deviations of the muon production vertex from the centroid and the physical processes which alter the muon track between the production vertex and sea level obscure this information somewhat. For this reason it is not possible to study this question on a muon-by-muon basis. For example, tracks with projected angles close to zero which exist at almost all $`x`$-distances from the origin would yield, upon extrapolation, an infinite height-of-origin. This problem can be circumvented by finding the mean $`x`$ within an interval of $`x`$ and dividing this mean by the corresponding mean of $`\mathrm{tan}\mathrm{\Theta }_{xz}`$. The averaging serves another purpose as well: the randomness of the various physical processes, such as scattering, tend to cancel in the averaging.

### 3.3 Angular resolution for a spectrum of primary gamma ray energies

The results of the preceding sections can be combined to obtain the expected angular resolution for a spectrum of primary gamma rays. Various steps in the calculation are shown in Fig. 6a. Here, a differential energy spectrum $`d\mathrm{\Phi }_\gamma /dE_\gamma E_\gamma ^\alpha `$ with a spectral index of $`\alpha =2.41`$ is assumed; this $`\alpha `$ corresponds to the average of the spectral indices reported in the Third EGRET Catalog . Folding $`d\mathrm{\Phi }_\gamma /dE_\gamma `$ with the number of muons reaching sea level per primary photon, $`N_\mu /N_\gamma `$, yields the number of muons at sea level as function of the primary gamma ray energy, $`d\mathrm{\Phi }_\mu /dE_\gamma `$, with $`\mathrm{\Phi }_\mu =N_\mu /N_\gamma \times \mathrm{\Phi }_\gamma `$. Note that the differential spectra in Fig. 6a are multiplied by $`E_\gamma `$.

The flatness of the differential muon flux per gamma at sea level ($`E_\gamma \times d\mathrm{\Phi }_\mu /dE_\gamma `$) above 10 GeV signifies that the muons originate rather uniformly from a broad range of primary energies. Below 10 GeV, the muon flux decreases steeply as $`N_\mu /N_\gamma `$ rapidly approaches zero. Harder spectral indices ($`\alpha <2.41`$) would enhance the muon flux from higher primary energies and, conversely, softer indices would enhance lower energies.

Furthermore, in Fig. 6a the half-widths at half-maximum-height $`\delta \mathrm{\Theta }_{xz}`$ (HWHM) (given in Table 2 in the “$`E_\mu >0`$” column) are plotted. These widths are multiplied by the differential muon flux which is also shown in Fig. 6a. As can be seen, the larger angular widths at the lower energies ($`E_\gamma <3`$ GeV) do not contribute due to the small sea level muon flux at these energies.

The average angular resolution (i.e., average angular width) in the North-down plane, $`\delta \mathrm{\Theta }_{xz}(\mathrm{HWHM})`$, for a differential gamma ray spectrum with an index $`\alpha `$ is obtained from

$$\delta \mathrm{\Theta }_{xz}(\mathrm{HWHM})=\frac{\delta \mathrm{\Theta }_{xz}\times N_\mu /N_\gamma \times E_\gamma ^\alpha 𝑑E_\gamma }{N_\mu /N_\gamma \times E_\gamma ^\alpha 𝑑E_\gamma }$$
(1)
and shown in Fig. 6b. As can be seen, the resolution is about constant ($`1^{}`$) up to a spectral index of 2 and then rises for larger indices due to the increasing importance of the wider angular widths at low energy.

### 3.4 Effect of the magnetic field

The magnetic field of the Earth deflects positively (negatively) charged particles primarily eastward (westward) causing the angular distributions in the projected angle $`\mathrm{\Theta }_{yz}`$ (i.e., projected onto the East-down plane) to become wider. The size of the effect depends on the strength and direction of the field and, therefore, varies with geographic location. To obtain a measure of the magnitude of the added angular widths and how it depends on geographic location, three different locations are considered and summarized in Table 3. In order to estimate the maximum possible deflection, a location of maximum northward field component was chosen. Due to the offset of the magnetic field axis from the Earth’s center, this location is at about $`9^{}`$N, $`100^{}`$E (below referred to as “maximum”). In addition, two locations with intermediate magnetic fields were chosen: $`42^{}`$N, $`86^{}`$W (GRAND) and $`36^{}`$N, $`106^{}`$W (MILAGRO). The magnetic field components for these locations were obtained from the latest revision of the International Geomagnetic Reference Field (IGRF)  and are listed in Table 3. Here, the components of the magnetic field, B, are given in nanotesla (nT) and the components $`B_x,B_y,B_z`$ are directed North, East, and toward the center of the Earth, respectively. For uniformity of comparison, all values are for sea level elevation.

FLUKA is able to transport particles in arbitrary magnetic fields. However, since the variation of the field in the atmosphere with height above sea level and with time is only minor (about 4%), for simplicity the magnetic field used in the calculation for a certain location was assumed to be constant throughout the whole geometry. In addition, although the MILAGRO experiment is located at a height of 2170 m above sea level, all results reported for that location refer to sea level as this allows a direct comparison to the results at the other locations.

The angular distributions at sea level for the three magnetic field conditions and for primary gamma ray energies of 10 GeV and 1 TeV are shown in Fig. 7. In addition, the distributions obtained without magnetic field are given. As can be seen, for a 10 GeV gamma primary, the effect of the Earth’s magnetic field is considerable and changes with geographic location; as the northward component of the magnetic field increases, the width increases. However, at 1 TeV, there is only a slight widening of the distributions with increasing the North component of the field. To within the statistical uncertainties of the calculations, the total number of muons which reach sea level per primary gamma ray is independent of the effect of the Earth’s magnetic field for the two energies studied in this section (10 GeV and 1 TeV).

The half-widths of the angular distributions in the North-down ($`\mathrm{\Theta }_{xz}`$) and East-down planes ($`\mathrm{\Theta }_{yz}`$) containing 68% of the muons are presented for the different locations in Table 4. Again, the

widths are given for distributions containing only muons above certain kinetic energy thresholds (including zero threshold, or no cut). In addition, the corresponding values from the calculations without the effect of the magnetic field are given. The widths in the North-down plane are much less affected than in the East-down plane which contains the dominant effect of the magnetic deflection. This effect increases with increasing values of $`B_x`$ or decreasing values of muon momentum. In addition to the widths of the muon’s sea level angular distributions, Table 4 also shows the values at their production vertex (birth, E-d.(b)). Interestingly, these values are almost constant for the different field conditions. The effect of the magnetic field on the angular distributions is therefore mainly the deflection of the muon over its longer path and not the deflection of the parent pions or other charged particles preceding in the cascade due to their shorter path lengths and higher momenta as compared to the muons. In addition, the dominant factor in the final width (except for the highest magnetic field value) is the muon’s angular distribution at its production (birth) point.

## 4 The effect of pressure and temperature variations

Air pressure and temperature changes cause small variations in the density profile of the atmosphere and thus affect the muon flux at sea level. In order to obtain a quantitative estimate of these two effects, the following cases were studied: (i) the air density in each layer was increased by 3% while all other variables (geometry, cutoffs, etc.) were kept constant and (ii) the height above sea level of each slab boundary was raised by 5% and the corresponding density decreased in the same proportion such that the total thickness of air (in g/cm2) above the surface of the Earth remained constant. The former modification simulates an air pressure increase ($`p`$) whereas the latter simulates an overall 5% increase in absolute temperature ($`T`$) of the atmosphere. Since muons are produced at atmospheric heights far above those which determine the weather at the surface of the Earth, the mean $`T`$ which is involved in this calculation has little correlation with the temperature at the Earth’s surface. Thus the temperature effect calculated here is not expected to correlate with the Earth’s surface temperature but would require temperature measurements of the Earth’s atmosphere averaged from 0 to 20 km above the Earth’s surface as obtained, for example, in weather balloon measurements. The results presented in this section allow an investigation of, or corrections for, the muon flux due to small variations in the atmospheric pressure and height-averaged temperature.

The results for the total muon yield per primary gamma ray are summarized for two energies (10 GeV and 1 TeV) in Table 5. In all

cases the muon yields decreased by 3% to 6%, the decrease being more pronounced at the lower energy. As for increased pressure, a higher density increases the interaction probability of the pions and thus decreases the probability of decay, i.e., fewer muons are produced. As well, the added thickness of air which the muons must traverse is a barrier which some will fail to overcome. On the other hand, raising the boundaries of the layers and reducing their density causes two counterbalancing effects: (i) the pions propagate in a less dense medium and thus have a greater probability to decay rather than interact thus producing more muons, and (ii) the muons must propagate longer distances resulting in a greater probability that they will decay before reaching ground level. The calculated temperature dependence of the muon multiplicity ratio suggests that the enhanced muon decay probability dominates at these energies and is the more dominant at the lower energy.

## 5 General properties of muon production in gamma ray showers

Information on the muons at sea level and on their ancestors (i.e., the parent, grandparent, etc., see Sec. 2.2) allows a more detailed study of the properties of these ancestors.

### 5.1 The ancestors of the muons

Table 6 shows which parent contributes through its decay to the detected muons for the different primary gamma ray energies.

At the lowest energy (1 GeV) about 89% of the muons are produced in decays of negative pions and 11% in decays of positive pions. This asymmetry can be explained by the different interaction cross sections at low energy of pions of either charge and by the fact that these pions are mainly produced in secondary interactions of neutrons (see below). Correspondingly, the probability for positive pions to interact instead of decay is larger than for negative pions. Above 3 GeV this asymmetry disappears and decaying kaons begin to contribute to the sea level muon flux. The latter contribution increases with energy and amounts to about 5% at 10 TeV.

It is interesting to note that at all energies FLUKA predicts the relative contribution of a positive kaon parent to be always larger than that of a negative kaon by about 20%. This effect is due to the properties of the DPM describing inelastic hadronic interactions within FLUKA. In particular, it is a feature of the Reggeon contribution which describes particle production by one quark-diquark string stretched between a valence quark of the fluctuating photon and a diquark of a target nucleon. In this picture, kaon production involves the creation of a $`s\overline{s}`$ quark-antiquark pair and, in case of negative kaons, also the creation of an $`u\overline{u}`$ pair. On the other hand, positive kaons can be readily formed also by a $`u`$-quark of the fluctuating photon leading to the observed asymmetry.

The fractional contribution to the muon’s grandparent is given in Table 7. The calculations for 1 GeV show distinctly different features. About 97% of the muons at sea level originate from mesons produced in interactions of neutrons in the close vicinity of the detector. At all other energies photoproduction dominates the picture with a contribution decreasing with energy from 99% at 3 GeV to 60% at the highest energy; the remaining fraction is mainly from pions.

### 5.2 Distributions of the number of generations

The distributions of the generation number of the grandparent of the sea level muon is shown in Fig. 8. The lower the energy of the primary gamma ray the smaller is the atmospheric shower resulting in a relatively narrow distribution. This distribution is peaked at generation one for energies below a few hundred GeV where the parent is mainly produced in a photoproduction process of the primary photon (see also Table 7). At higher energies secondary hadron interactions contributes significantly to the production of the parent. These secondary hadrons are mainly of third or fourth generation and cause a shift of the peak of the total distribution at TeV energies. The large tail at these energies which extends up to more than 100 generations reflects photoproduction processes of secondary photons in the large electromagnetic shower.

### 5.3 The muon production heights

The distributions of heights above sea level at which the detected muons were produced are shown for four primary energies in Fig. 9. As expected, at high energy there is a decrease in height with increasing energy of the primary gamma ray because of the logarithmic increase of the shower length. However, as the energy of the primary gamma ray decreases below 10 GeV, the energy of the muon parent is close to the minimum energy required for the muon to penetrate the blanket of air between its production and sea level. Hence, the most probable height of muon production begins to decrease with energies decreasing below 10 GeV. The blip-up near zero height, most pronounced in the 3 and 10 GeV distributions, represents the excess contribution from very low energy muons.

### 5.4 Radial distributions

The radial distributions of the muons at their production vertices and at sea level are shown in Fig. 10 for five different gamma ray energies. The radial distance, $`R`$, is defined with respect to the shower axis ($`z`$ axis). The quantity $`dN/dA`$ denotes the number of muons per unit area and per primary gamma ray. The radial distributions of the production vertices are labeled “production.”

All distributions extend to more than 10 km. Whereas the sea level distributions have a relatively flat shape below $`R=2`$ km the production vertex distributions are increasing toward the shower axis and exhibit a change in slope or discontinuity at about 500 m which is most pronounced at low primary energy. This discontinuity indicates that two components with different shapes contribute. In order to further investigate this feature the muon production vertices were scored separately for those muons which originated from mesons produced in photoproduction processes and those from mesons produced by interactions of nucleons or other hadrons. The contributions from these components are shown for a primary energy of 10 GeV in Fig. 11. As can be clearly seen in the production vertex distributions (Fig. 11a) the photoproduction component dominates at small radii and the component due to interacting nucleons constitutes the tail at large $`R`$. The difference is smaller at sea level (Fig. 11b) due to multiple scattering of the muons in air.

### 5.5 Energy distribution of the muons

The kinetic energy distributions of muons at sea level and at the production vertex of these muons are shown in Fig. 12. Note, that the distributions are multiplied by the muon energy $`E_\mu `$ in order to enhance possible spectral structures at the higher energies. Whereas the sea level spectra are a smooth function with energy, the distributions of the muon energies at their birth exhibit a two-component structure. This two-component structure is again due to the superposition of muons from mesons generated in photoproduction processes and those from interacting hadrons. The two contributions are plotted separately in Fig. 13 for 10 GeV primary photon energy. It is clear that muons from photoproduction processes have a higher average energy.

## 6 Conclusions

Air showers caused by cosmic gamma rays with energies below 10 TeV were simulated using the Monte Carlo code FLUKA. The primary gamma rays were assumed to enter the atmosphere vertically. The reliability of FLUKA predictions in this energy region was explored by comparing its predictions for proton- and nuclei-induced air showers with experimental data; good agreement was found.

Many general properties of muon production in atmospheric gamma ray showers were studied. As expected, most of the muons which reach sea level are decay products of charged pions; the contribution from kaons increases as the primary energy increases. Below 100 GeV primary energy, the decaying mesons are produced mainly in photoproduction processes. At higher energies, hadronic interactions of nucleons and pions contribute significantly. Variations in atmospheric pressure and temperature were investigated: A 3.0% increase in atmospheric pressure decreases the muon yield at sea level by 4.0% for 10 GeV primaries (3.2% at 1 TeV). Similarly, a 5.0% increase in the height-averaged temperature decreases the muon yield by 6.0% (3.6% at 1 TeV).

The main goal of the study was to determine the angular resolution for primary cosmic gamma rays from measurements of secondary muon angles at ground level. The angular distribution of muons at sea level narrows as the muon energy increases since low energy muons constitute the tails of the angular distributions. Thus, the angular resolution can be improved by eliminating the lower energy muons. The half-widths at half-maximum-height of the muon’s projected angles are $`1.1^{}`$ above 100 GeV but rise to a value of $`3.2^{}`$ as the primary gamma ray energy is lowered to 10 GeV in the absence of the Earth’s magnetic field deflection. Primary energies below 10 GeV essentially do not contribute muons at ground level. The effects of magnetic bending were considered separately as they depend on several factors; for example, they are more important in the East-down plane, for low primary energies, low muon energies, and larger North-components of the Earth’s magnetic field.

The width of the angular distribution can be narrowed if the correlation of $`\mathrm{\Theta }_{xz}`$ (angle projected onto the North-down plane) with $`x`$ is taken into consideration (assuming the location of the core of the shower is measured). This narrowing can be accomplished either by eliminating large $`x`$ and $`y`$ values from the data or by compensating the measured angles knowing their expected mean values versus $`x`$ and $`y`$ as calculated from FLUKA.

Finally, the average muon angular correlation was calculated for a spectrum of primary gamma rays characterized by a spectral index. For a soft spectrum, e.g. a differential spectral index $`\alpha =3.0`$, the average projected angular width (HWHM) is $`2.6^{}`$; for a mean value, $`\alpha =2.41`$, the width is $`1.6^{}`$. If the spectrum is harder, e.g. $`\alpha 2.0`$, the angular width improves to a constant value of $`1.0^{}`$. Corresponding values for a space angle resolution ($`\delta \mathrm{\Theta }_R`$(HWHM)) are $`\sqrt{2}`$ larger. As mentioned above, if data with larger $`x`$ and $`y`$ values and/or the lower muon energies are removed, this angular correlation can be improved. This precise calculation of the angular correlations provides experiments which study the angular location of primary gamma ray sources by measuring muon angles with information on the angular resolution which is difficult to obtain experimentally.

## Acknowledgments

Part of this work was supported by the Department of Energy under contract DE-AC03-76SF00515. Project GRAND is funded through grants from the University of Notre Dame and private donations.

## Figure Captions

1. Dependence of the negative muon flux produced by hadronic primaries on the depth in the atmosphere shown for different intervals of the muon momentum. FLUKA results (solid lines) are compared to data obtained by the CAPRICE experiment  (points). The comparison is shown for different intervals in muon momentum. The curves and data points for the first five momentum intervals were shifted by constant factors as indicated.
2. a) Energy spectra of positive and negative muons at 886 g/cm2. FLUKA results (histograms) are compared to data obtained by the CAPRICE experiment . $`p_{\mathrm{Lab}}`$ is the detected muon momentum. b) Muon charge ratio as a function of kinetic energy measured by the CAPRICE experiment  (points) and calculated with FLUKA (solid line).
3. Angular distribution of muons at sea level. $`\mathrm{\Theta }_{xz}`$ is the angle of the detected muon projected onto the North-down or $`xz`$ plane. The distributions are given for different primary photon energies and are normalized to unit area. Effects of the Earth’s magnetic field are not included in this figure. Units of angle are in degrees.
4. (a) Angular distribution of muons at sea level for 10 GeV and (b) 10 TeV primary gamma ray energy. The distributions labeled “total” are the same as shown in Fig. 3. In addition, the contributions from muons with kinetic energies below 2 GeV and above 2 GeV are presented. Units of angle are in degrees.
5. Angular distributions of muons at sea level given for different intervals in the $`x`$ coordinate. Results are shown for primary gamma ray energies of (a) 10 GeV and (b) 1 TeV. All distributions are normalized to unit area. Units of angle are in degrees.
6. (a) The figure shows the following quantities as a function of the primary gamma ray energy: the energy spectrum of primary gamma rays ($`d\mathrm{\Phi }_\gamma /dE_\gamma `$) with a shape given by a differential spectral index $`\alpha =2.41`$ and arbitrary normalization, the multiplicity of muons at sea level per primary gamma ray ($`N_\mu /N_\gamma `$), the result of folding the multiplicity with the primary spectrum ($`d\mathrm{\Phi }_\mu /dE_\gamma =N_\mu /N_\gamma \times d\mathrm{\Phi }_\gamma /dE_\gamma `$; thin solid line with circles), the half-widths at half-maximum-height for all muon energies ($`\delta \mathrm{\Theta }_{xz}`$(HWHM)) in degrees, and the result of folding these widths with $`d\mathrm{\Phi }_\mu /dE_\gamma `$ (dot-dashed lines with diamond points). The calculated values are joined by lines to guide the eye. Note that all differential fluxes are multiplied by $`E_\gamma `$. (b) Average half-width at half-maximum-height for the angle projected onto the North-down plane as function of the spectral index of the primary spectrum; the ordinate is in units of degrees. Values for the space angle resolution ($`\delta \mathrm{\Theta }_R`$(HWHM)) would be $`\sqrt{2}`$ times larger.
7. Muon angular distribution at sea level including the effect of the Earth’s magnetic field for (a) 10 GeV and (b) 1 TeV primary gamma ray energy. The distributions are shown as a function of the muon’s angle projected onto the $`yz`$ (East-down) plane, $`\mathrm{\Theta }_{yz}`$, and for three different geographic locations: (1) $`42^{}`$N, $`86^{}`$W (GRAND), (2) $`36^{}`$N, $`106^{}`$W (MILAGRO), and (3) $`9^{}`$N, $`100^{}`$E (maximum). The magnetic deflection for $`\mathrm{\Theta }_{xz}`$ (North-down) is minimal and not shown. In addition, the distribution obtained without magnetic field (labeled “no field”) is given. Units of angle are in degrees.
8. Distribution in the number of generations in the shower (see text) carried by the particles creating the parent-meson. FLUKA results are shown for different energies of the primary photon.
9. Distribution of the muon production heights. Results are given for primary photon energies of $`3,10,100`$, and $`1000`$ GeV, respectively, and are normalized per primary photon. The top of the atmosphere is at 80 km in the calculations.
10. Radial distribution of the muon production vertices (i.e. meson decay vertices, dotted histograms) and of the muons at sea level (solid histograms). Results are given for primary photon energies of $`3,10,100,10^3,`$ and $`10^4`$ GeV, respectively (from the bottom curve to the top). Values have been normalized per primary photon.
11. Radial distribution of the muon production vertices for 10 GeV primary photons at normal incidence. The distributions are given for (a) the muon production vertex and (b) for sea level. In addition to the total distribution, the contributions from different “parents” of the decaying meson are given. For example, the histograms labeled “photons” show the radial distributions of muons from mesons which were produced by photonuclear interactions. Distributions are normalized per primary photon.
12. Kinetic energy spectra of muons at their production vertices (dotted histograms) and at sea level (solid histograms). Values are normalized per primary photon.
13. Kinetic energy spectra of muons from 10 GeV primary photons at normal incidence. The spectra are shown for (a) the production vertex and (b) for sea level. In addition to the total spectra which are identical to those shown in Fig. 12, the contributions from different “parents” of the decaying meson are given (similar to Fig. 11).