{'fraction_non_alphanumeric': 0.028375574156543636, 'fraction_numerical': 0.007717596586537341, 'mean_word_length': 4.433928797668821, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulations of black hole air showers in cosmic ray detectors

We present a comprehensive study of TeV black hole events in Earth's atmosphere originated by cosmic rays of very high energy. An advanced fortran Monte Carlo code is developed and used to simulate black hole extensive air showers from ultrahigh-energy neutrino-nucleon interactions. We investigate the characteristics of these events, compare the black hole air showers to standard model air showers, and test different theoretical and phenomenological models of black hole formation and evolution. The main features of black hole air showers are found to be independent of the model considered. No significant differences between models are likely to be observed at fluorescence telescopes and/or ground arrays. We also discuss the tau ``double bang'' signature in black hole air showers. We find that the energy deposited in the second bang is too small to produce a detectable peak. Our results show that the theory of TeV-scale black holes in ultrahigh-energy cosmic rays leads to robust predictions, but the fine prints of new physics are hardly to be investigated through atmospheric black hole events in the near future.

Event observables would be the secondary products of the BH or brane after their formation, i.e. extensive air showers originated by field emission in the decay phase. Hawking evaporation provides an emission mechanism for BHs [14]. SM quanta are emitted in the visible three-brane and can be detected. Branes may or may not evaporate, depending on their properties. However, the decay spectrum of massive excitations in string theories has been shown to be thermal [15]. This suggests that BH and brane decay signatures may be similar.
The observational signatures of BH events in the semiclassical approximation have been investigated in a number of recent publications. In Ref. [9] the authors found that BH interactions generate different air showers from SM interactions. BH air showers tend to rise faster and have larger muon content. A BH event produces a hadronic air shower occurring at a much greater depth in the atmosphere, i.e., a very deeply penetrating hadronic air shower. However, the inability of realistic detectors to observe the first interaction point hides most of the difference between BH and SM air showers. Given that present observatories are not large enough to study a large number of neutrino events, discrimination of BH and SM events is likely not to be achieved in the near future. Another characteristic of BH air showers is τ generation. Although the rate for these events is low if the ultrahigh-energy neutrino flux is at the level of the expected cosmogenic neutrinos, unusual "double-bang" air showers could signal a departure from SM interactions.
Whereas the semiclassical picture seems reasonable, the actual physics of subnuclear BH formation could be very different. In the last year or two, significant advances in the understanding of microscopic BH formation and air shower evolution have appeared in the literature. It is thus timely and worthwile to re-examine the observational signatures of BH air showers. To this purpose, we developed a thorough fortran Monte Carlo (MC) code to simulate the air showers induced by BH formation in neutrino-air collisions, which includes these theoretical refinements [16]. The MC has the same structure of the MC used in Ref. [9]. The code generates observable secondaries from BH evaporation using the PYTHIA generator [17]. These secondaries are then injected into the AIRES simulator [18] as primaries for the final air shower.
The purpose of our study is threefold. Firstly, we want to confirm the main findings of the previous investigation. Secondly, we want to test various proposals of BH vs. SM air shower discrimination that have appeared in the literature, such as the τ "double bang" effect [10]. Thirdly, we want to look for new ways of discriminating between different models of BH formation and evolution.
Our analysis will show that the main characteristics of BH air showers are essentially independent of the details of BH evolution. Because of large uncertainties and statistical fluctuations in air shower detection, it is also practically impossible to discriminate between alternative models of BH formation and evaporation. For instance, we will show that there is no significant observational difference between a model of BH formation based on the semiclassical black disk and the trapped-surface model [13], or between a model of BH evaporation with final explosive decay and stable remnant [19]. These results limit significantly the use of BH air showers (if they exist) to probe details of "new physics". We will also show that newly proposed signatures do not help in the task of discriminating BH vs. SM air shower detection. No observational trace of the "double bang" signature can be extracted from a realistic detector in the near future.

II. BASICS OF BH FORMATION AND EVOLUTION
In this section we briefly review the basics of BH formation and evolution, focusing on recent theoretical advances that have been included in the MC code.

A. BH formation and cross section at parton level
Thorne's hoop conjecture [20] states that a horizon forms when a mass M is compacted into a region with circumference smaller than twice the Schwarzschild radius R(M ) in any direction. At subnuclear level, this can be achieved by scattering two partons (ij) on the brane with c.m. energy √ s ij > M and impact parameter b < R(M ). This event can formally be described by the process ij → BH + E(X), where E(X) denotes collisional energy that does not contribute to the BH mass. This energy includes a bulk component of gravitational radiation and perhaps non-SM gauge fields, and a brane component of SM fields. If E(X) is zero, the hoop conjecture implies that the cross section for BH production is independent of the impact parameter (as long as b < R(M )) and equal to the geometrical black disk σ BD (s ij , n). If E(X) = 0, the cross section depends on the impact parameter, and is expected to be smaller than the black disk cross section. It is worth stressing that this picture is correct only if the BH is larger than the Compton length of the colliding quanta. (For discussions on the effect of wave packet size on the BH formation process, see Ref. [21].) A precise calculation of the collisional energy loss is essential to understand BH formation.
Many papers have been devoted to improve or disprove the hoop conjecture. The most popular model is currently the trapped-surface model [22,23,24], although alternative techniques have been explored [13]. The trapped-surface approach gives an upper bound on the gravitational component of E(X) by modelling the incoming partons as two Aichelburg-Sexl shock waves [25]. The Aichelburg-Sexl wave is obtained by boosting the Schwarzschild solution to the speed of light at fixed energy. The resulting metric describes a plane-fronted gravitational shock wave corresponding to the Lorentz-contracted longitudinal gravitational field. The parton scattering is simulated by superposing two shock waves travelling in opposite directions. The union of these shock waves defines a closed trapped-surface that allows to set a lower bound on the BH mass. The collisional energy loss depends on the impact parameter and increases as the number of spacetime dimensions increases. The BH mass monotonically decreases with the impact parameter from a maximum of about 60-70% of the c.m. energy for head-on collisions.
The trapped-surface result is consistent within one order of magnitude with the hoop conjecture. However, the partons are assumed to be pointlike, massless, spinless, and electrically neutral. The pointlike assumption fails for directions transversal to the motion [26]. Colliding partons generally have spin and charge. While size and spin effects are expected to be mostly relevant around the Planck energy, charge effects could dominate at higher energy. It should also be kept in mind that the trapped-surface model provides only a lower bound on the BH mass. An accurate estimate of the gravitational collisional energy loss would require the use of the full non-linear Einstein equations in higher dimensions. Since this is a virtually impossible task, alternative approximated models have been investigated. The gravitational energy emission in a hard instantaneous collision can be computed in the linearized limit [27]. This approach suggests that the trapped surface method overestimates the gravitational energy emitted in the process. For head-on collisions, the instantaneous method predicts the gravitational energy loss to be only about 10% of the c.m. energy. This result is in agreement with a perturbative calculation modelling the parton-parton collision as a plunge of a relativistic test particle into a BH with mass equal to the c.m. energy [28].
In conclusion, known methods for the estimate of the gravitational loss in relativistic scattering at parton level give a BH mass ranging between 60% and 100% of the c.m. energy. Today, the trapped-surface value and the black disk value can be considered as the lower and upper bounds on the BH mass, respectively.

B. Cross section at nucleon level
The total cross section for a super-Planckian event involving a nucleon is obtained by integrating the above cross section over the parton distribution functions. BHs formed in a neutrino-nucleon collision may dominate over the SM processes and stand a fair chance of detection. On the contrary, the branching ratio of the BH cross section in a nucleon-nucleon collision is ∼ 10 −9 . Therefore, BH detection in nucleon-nucleon interactions cannot be achieved with current and next generation detectors due to the low flux of UHECRs.
If the BH mass depends on the impact parameter, the generally accepted formula for the total cross section of the neutrino-nucleon process is where q i (x, −Q 2 ) are the Parton Distribution Functions (PDFs) [29,30] with four-momentum transfer squared −Q 2 , and fraction of the nucleon's momentum carried by the i-th parton √ x. z is the impact parameter normalized to its maximum value and x m = M 2 min /(sy 2 (z)), where y(z) and M min are the fraction of c.m. energy trapped into the BH and the minimum-allowed mass of the gravitational object, respectively. F is a form factor. The total cross section for the black disk model is obtained by setting F = 1 and y 2 (z) = 1.
Different sets of PDFs are defined in the literature. The PDFs are not known at energies above the TeV and for values of momentum transfer expected in BH formation. Equation (1) is usually calculated by imposing a cut-off at these values. The PDFs also suffer from uncertainties at any momentum transfer (∼ 10%) [9] and from the ambiguity in the definition of Q [31]. The momentum transfer is usually set to the BH mass or the inverse of the Schwarzschild radius. Although recent literature inclines toward the latter, there are no definite arguments to prefer either one or to exclude alternative choices. The uncertainty due to the ambiguity in the definition of the momentum transfer is about ∼ 10 − 20% [8].
The form factor and the amount of trapped energy depend in principle on energy, gravitational scale, geometry and physical properties of the extra dimensions and gravitational object. The trapped-surface method gives numerical values of order unity for these quantities. (See Refs. [22,23] and discussion above). However, these results depend on the way the trapped-surface is identified. Other models [32] give values which are more or less consistent with the trapped-surface method. With the lack of further insight, it is common practise in the literature to either choose the trapped-surface result or the simple black disk model.
The lower cutoff on the fraction of the nucleon momentum carried by the partons is set by the minimum-allowed (formation) mass of the gravitational object, M min . This threshold is expected to roughly coincide with the mass for which the semiclassical description is valid. This conclusion is motivated for spherically symmetric BHs by the following argument [33]: For M min /M ⋆ few, the Hawking entropy of the BH should be large enough to neglect strong gravitational effects. The semiclassical results are then extrapolated for smaller values with the assumption that the BH or its Planckian progenitor decays on the brane. However, this argument is based on Hawking's semiclassical theory and may not be valid at energies equal to few times the Planck mass. For example, the existence of a minimum spacetime length l m implies the lower bound on the BH mass [19,34]: BHs with mass less than M ml do not exist, since their horizon radius would fall below the minimum-allowed length. At fixed M ⋆ , the minimum-allowed mass grows as a power of l n+1 m . For n = 6 or 7 and l m M ⋆ 1, it follows M ml ∼ M min ≫ M ⋆ .

C. BH evolution
It is believed that the decay of microscopic BHs happens in four distinct stages: I. radiation of excess multipole moments (balding phase); II. spin-down; III. Hawking evaporation; IV. final explosion or formation of a BH remnant. Phases I-III rely on semiclassical results, provided that the entropy is sufficiently large. Phase IV is in the realm of quantum gravity.
Although some progress has been made, the understanding of balding phase and spin-down phase is still fragmentary. For example, the emission of radiation from a (n + 4)-dimensional rotating BH on the brane is known only for spin-0 fields [35]. Due to these limitations, some of these theoretical results cannot be implemented in MC simulations at the present stage. Moreover, phase I is not expected to lead to a significant amount of energy loss and the observational uncertainties (see below) are likely to dominate the theoretical uncertainties in phase II.
Many papers have been devoted to the investigation of the Hawking phase. Although several analytical and numerical results have been obtained [36], from the viewpoint of numerical simulations the situation is similar to the balding and spin-down phases. For instance, greybody factors for the graviton are not fully known in (n+4)-dimensions even for the spherically symmetric BH. This precludes their use in numerical codes, where a consistent use of greybody factors is required. The field content at trans-Planckian energies is also not known. Onset of supersymmetry, for example, could lead to other evaporation channels for the BH and large emission of non-SM or undetectable quanta during the decay phase. Finally, quantum effects may also affect the emission of visible quanta on the brane.
Quantum corrections to the Hawking phase can be phenomenologically described by assuming the existence of a minimum length of the order of the Planck length [37]. The existence of a minimum scale is a common consequence of most (if not all) theories of quantum gravity such as string theory, non-commutative geometry and loop quantum gravity. The presence of a cutoff at small spacetime distances leads to a modification of the uncertainty principle at Planck scales. Since the Hawking thermodynamical quantities can be derived by applying the uncertainty principle to the BH, the existence of a minimum length leads to corrections in the thermodynamical quantities [19,34].
At the end of the Hawking phase, the BH is expected to either non-thermally decay in a number n p of hard quanta or leave a remnant. In either case we must content ourselves with a phenomenological description, due to the lack of a theory of quantum gravity. The final non-thermal decay is usually described by setting a cutoff on the BH mass of the order of the Planck mass, Q min ∼ M ⋆ , and then equally distributing the energy Q min to a number n p of quanta. Since the decay is non-thermal, and in absence of any guidance from a theory of quantum gravity, the quanta are democratically chosen among the SM degrees of freedom. Note that Q min does not necessarily coincide with M min . The former gives the threshold for the onset of quantum gravity effects, whereas the latter gives the minimum-allowed mass of the classical object. From the above definitions, it follows M min ≥ Q min . The existence of a minimum length gives a natural means to set Q min . In that case, the modified thermodynamical quantities determine the endpoint of Hawking evaporation when the mass of the BH reaches M ml . This mass can be identified with the mass of the BH remnant [19,34].

III. BH GENERATOR
In this section we list the main characteristics of the MC generator used in the simulations [16]. The physics of BH formation and decay is determined by the following set of external parameters and switches in the MC code: The above parameters are briefly explained below. A more detailed explanation can be found at the MC generator web site [16].

A. BH formation and parton cross section
The MC does not require any lower or upper bound on the Planck mass M ⋆ . However, experimental constraints exclude values of M ⋆ 1 TeV and large values of M ⋆ do not allow BH formation in the atmosphere. Therefore, M ⋆ must be chosen with caution. Since n = 1 and n = 2 are excluded experimentally, and most of the theoretical models are limited to n ≤ 7, the number of extra dimensions n ranges from 3 to 7.
The MC includes three models for BH formation and cross section: Black disk, Yoshino-Nambu (YN) trappedsurface model [22], and Yoshino-Rychkov (YR) improved trapped-surface model [23]. This allows a comparison between air showers based on the black disk model [9] and air showers generated by BHs with significant gravitational loss at formation. Observable differences between different models of BH formation can be investigated, as suggested in Ref. [23].
The minimum BH mass M min is set in units of M ⋆ or M ml (if a minimum length is present, see below). This parameter is always larger than one, i.e. M min ≥ Max(M ⋆ , M ml ).

B. Total and differential cross section
The distribution of the initial BH masses is given by the differential cross section dσ/dM BH , where M BH = √ xs. The MC uses the (stable) cteq5 PDF distribution [30,38]. Since the use of different PDF distributions produces an insignificant uncertainty in the total and differential cross sections, other PDF distributions are not implemented in the MC. The uncertainty due to the choice of the momentum transfer is generally larger. Therefore, a switch allows to choose between BH mass or inverse of the Schwarzschild radius as definition of momentum transfer.
The part of c.m. energy of the neutrino-nucleon system which is not trapped or lost in gravitational radiation at formation is attributed to the nucleon remnant. For sake of simplicity, only neutrino-proton collisions are implemented in the MC code. A neutrino-neutron collision does not produce significant statistical differences in the nucleon remnant compared to a neutrino-proton collision. The proton remnant is successfully fragmented according to QCD in mesons/baryons (see, for example, Ref. [17]) and then is decayed with the PYTHIA generator along with the quanta created in the BH evaporation process.

C. BH evaporation
Due to the lack of results for the balding and spin-down phases described above, energy losses in these stages are assumed to be either negligible or included in the energy loss during formation. This is a reasonable assumption since the trapped-surface model likely overestimates the actual energy loss. Balding and spin-down effects are also not expected to produce detectable differences in BH air showers, given experimental uncertainties and statistical fluctuations. Nevertheless, keeping an open mind, we plan to include balding and spin-down effects in updated versions of the code, as soon as theoretical results become available.
A similar conservative approach is used in the Hawking phase, where only thermally-averaged greybody factors in four dimensions are implemented in the MC. This is justified by consistency reasons in the code (the full greybody factors for all fields are not known). As the SM fields are emitted on the brane, and given the observational uncertainties, the difference between thermally-averaged and exact greybody factors is not expected to be detectable. The particle content at trans-Planckian energies is assumed to be the minimal SU (3) × SU (2) × U (1) SM with three families and a single Higgs boson. The degrees of freedom c i and the thermally-averaged greybody factors Γ i are listed in Table I. The decay multiplicities per species N i are assigned according to the prescription of Ref. [39]: where f i (m) = 1 or 1 − 2 1−m for bosons or fermions, and the total multiplicity N is where S is the initial entropy of the BH.
The presence of a minimum length affects the BH evolution in the Hawking phase. If no minimum length is present, the MC evaporates the BH according to the Hawking theory. Alternatively, the BH evolution proceeds according to the modified thermodynamics of Ref. [19,34]. In both cases the evaporation ends when the BH reaches the mass Q min . This is set in units of M ⋆ (M ml ) if the minimum length is zero (nonzero). Note that the BH minimum formation mass M min and the endpoint of Hawking evaporation Q min are independent parameters. Four-momentum is conserved at  each step in the evaporation process by taking into account the recoil of the BH on the brane due to the emission of the Hawking quanta. The initial energy of the BH is distributed democratically among all the Hawking quanta with a random smearing of ±10%. This smearing factor is introduced on a purely phenomenological basis to take into account quantum uncertainties in the emission of each quantum.

D. BH final decay
The MC code allows for two different choices of final BH decay: Final explosion in a number n p of quanta or BH remnant. If n p = 0, the BH settles down to a remnant with mass Q min . If n p = 1. . . 18, the BH decays in a number n p of quanta by a n-body process with total c.m. energy equal to Q min .
A switch controls conservation of EM charge in the decay process (Hawking evaporation + final decay). The purpose of this switch is to allow for the existence of a charged BH remnant. If the EM charge is not conserved and n p = 0, the BH remnant carries a charge Q R , where 1 − Q R = Q H + Q N is the sum of the EM charge of the Hawking quanta plus the charge of the nucleon remnant. If the EM charge is conserved and n p = 0, the BH remnant is assumed to be electrically neutral, i.e. Q R = 0: The absolute value of the total charge in the Hawking quanta is |Q H | ≤ 2e/3 and Q N = 1 − Q H . This is justified from the fact that the BH charge should have been shed earlier in the evaporation process. (See, however, Ref. [40] for a different viewpoint.) It should be stressed that the air shower phenomenology of a charge remnant is not known and it is not clear how to track it in the atmosphere in a meaningful way.

IV. BASICS OF NEUTRINO AIR SHOWERS
This section presents the essentials of the theory and phenomenology of UHECR neutrino air showers.

A. Physics of neutrino air showers
UHECRs are believed to be a composite of protons and heavier nuclei. Ultrahigh-energy neutrinos are created as these UHECRs interact with the cosmic microwave background through photopion production (protons) or the infrared background (iron nuclei). A cutoff in the energy spectrum is expected at the threshold energy of the photopion production, known as the Greisen-Zatsepin-Kuzmin (GZK) cutoff. These neutrinos are called cosmogenic or GZK neutrinos. Cosmogenic neutrinos are almost "guaranteed" to exist, though they have not been observed yet. They are the most likely source of neutrinos to produce BHs in the atmosphere. The peak of the cosmogenic neutrino flux is around 10 17−18 eV (c.m. energy ∼ 10 -50 TeV). The flux depends on a number of parameters of the UHECR source such as spatial distribution, injection spectrum, abundance, maximum energy and cosmological evolution. These factors can affect the flux even by a couple of orders of magnitude. The depth of the first interaction point, X 0 , depends on the total cross section of the process. The column depth of Earth's atmosphere in the horizontal direction is 3.6 × 10 4 g cm −2 . The interaction length of a neutrino with energy E ν = 10 9 TeV is λ CC ≃ 1.1 × 10 7 g cm −2 for charged current (CC) interactions. The largest possible cross sections of BH events give shorter interaction lengths, but still larger than the column depth of Earth's atmosphere. Therefore, neutrinos interacting in the atmosphere can induce air showers at any X 0 . In contrast, SM hadronic interactions have large cross sections with X 0 high in the atmosphere. Considering deeply penetrating horizontal air showers effectively filters out SM hadronic air showers, while giving the most likely chance of interaction. The background for detecting BH air showers is limited to SM neutrino air showers.
The SM interaction channels are the CC and neutral current (NC) for all three flavors. The energy of the leading lepton in the final state is given by (1 − y)E ν , where the mean value of the inelasticity y is ≈ 0.2. The leading lepton in the NC interaction is a neutrino that does not contribute to the air shower and the cross section is lower than the cross section of the CC interaction, σ N C (E) ≈ 0.4 σ CC (E). The ν µ -CC produces a high energy µ that does not decay before reaching ground. The ν τ -CC produces a τ that also generally does not decay before reaching ground. Therefore, the most relevant background for BH air showers is the ν e -CC channel.

B. Air shower detection
Extensive air showers can be detected with fluorescence telescopes and ground arrays. Fluorescence telescopes observe the fluorescence light produced by the interaction of atmospheric nitrogen molecules with the EM component of the developing air shower. The fluorescence method pioneered by the Fly's Eye detector [41] and currently operated by HiRes [42] and the Pierre Auger Observatory (PAO) [43] is able to reconstruct the longitudinal development of the (mainly) e + e − component of the air shower. This technique provides a good estimate of the energy of the primary particle that initiates the air shower, since most of the energy of the air shower goes into the observable EM channel. This method also enables reconstruction of the shower maximum X m , i.e. the depth at which the cascade contains the maximum number of e + e − pairs. The quantity X m − X 0 is sensitive to the type of primary particle, its energy, and the kind of interaction initiating the cascade. The duty cycle is approximately 10%, as clear moonless nights are required.
Ground arrays record the "footprint" of the air shower. Various methods are used to detect charged particles on the ground. Some examples are plastic scintillators [44,45] and water Cerenkov tanks [43,46]. These detectors are spread over areas ranging from a few km 2 to a few thousand km 2 . For example, the fully functioning southern PAO will have 1600 water Cerenkov tanks each with surface area of 10 m 2 , covering 3000 km 2 . Arrival time, composition, and pattern of the ground signals are used to analyze the properties of the air shower. Ground arrays can be operated full time.
The best method of BH air shower detection is a combination of fluorescence and ground detectors, such as the PAO. In view of this, we study both the air shower longitudinal development and the muon content at ground level. For inclined air showers, the geomagnetic field affects the distribution of particles on the ground, which is very sensitive to the zenith angle. Therefore, we simply count the number of particles rather than study their distribution.

V. AIR SHOWER GENERATOR
The BH generator output consists of a list of elementary SM quanta which are decayed with PYTHIA. The secondaries of the BH decay (PYTHIA output) are boosted to the laboratory frame and injected in the air shower generator AIRES to obtain the air shower. In this section we describe in detail the air shower part of the simulation.
The AIRES code requires to input primary cosmic ray properties and simulation conditions. The required physical parameters are: In our simulations, the zenith angle is set to 70 • and X 0 is set to an altitude of 10 km, corresponding to a slant depth of 780 g cm −2 . The BH interaction takes place at the injection altitude. The threshold energies for tracking particles in the air showers are 100 keV for gamma rays, electrons and positrons, 1 MeV for muons, 1.5 MeV for mesons, and 150 MeV for nucleons. The geomagnetic field is set to the Pierre Auger Observatory (El Nihuil site). The thinning level is 10 −6 with weight limitation of 0.2. A more detailed explanation of each parameter and other possible options can be found in the AIRES manual.

A. BH air showers
The steps to initiate a BH air shower are: 1. The BH is decayed in the c.m. frame. The unstable quanta are hadronized or decayed instantaneously by PYTHIA, with the exception of top quarks and τ leptons. PYTHIA does not handle top quarks. Therefore, they are instantaneously decayed as t → bW before being injected in PYTHIA. The τ leptons are produced directly from BH evaporation, or from hadronization or decay of other particles. Depending on their energy, the τ s may decay before reaching the ground. In that case, they are decayed with PYTHIA but their secondaries are injected into AIRES at different atmospheric depths, according to their boost and free path (see below).
2. All secondaries from PYTHIA are boosted to the laboratory frame. The particles are tightly beamed due to their very high boost.
3. All secondaries are injected into AIRES as primaries of the air shower.
Note that neutral pions generated in the hadronization process are immediately decayed by PYTHIA in the c.m. frame. Their average energy in the laboratory frame is smaller than the critical energy, making them more likely to decay than interact.

B. SM air showers
The CC and NC SM air shower simulations follow Ref. [9]: 1. The differential cross section is integrated over the fraction of the total nucleon momentum carried by the parton for all possible values of y.
2. y is sampled from the previous distribution. The energy of the leading lepton is (1 − y)E ν .
3. The leading lepton of the CC interaction is injected into AIRES. The leading neutrino of the NC interaction is not observable and is not injected. The τ lepton in the ν τ -CC interaction is treated separately as in the BH air showers. The τ is decayed with PYTHIA and its secondaries are injected at the corresponding atmospheric depth.

VI. SIMULATION RESULTS: SM VS. BH AIR SHOWER DETECTION IN FLUORESCENCE DETECTORS
Simulations show that the characteristics of BH air showers and SM ν e -CC air showers with identical first interaction point X 0 are quite different. However, these differences can vanish with a suitable shift of X 0 for either of the two processes. As X 0 is not a fixed parameter for these interactions, the properties of BH air showers are of limited practical use for detection purposes. The characteristics of BH air showers and SM ν e -CC air showers are summarized in the following The BH air showers are similar to hadronic air showers. The ν e -CC air showers are comparable to air showers generated by photons. The hadronic nature of the BH air showers is due to the prevalence of hadronic channels in the Hawking evaporation phase. Their rapidity is due to the large number of hadrons initiating the air showers. On the contrary, the main interaction channels in the SM events are pair production and bremsstrahlung. These processes produce a smaller number of secondaries than a hadronic interaction. The Landau-Pomeranchuk-Migdal effect [47] also contributes in slowing the air shower for primary energies 10 7 TeV [48]. The BH air showers exhibit smaller fluctuations in X m than the ν e -CC air showers. Although the BH mass varies from shower to shower, their development is more or less stable because of the large number of BH secondaries. The BH air showers can be viewed as a superposition of many air showers with less energy. The larger fluctuations from shower to shower in the SM process are due to the fluctuations in the energy (1 − y)E ν carried by the leading lepton. The large variations in the total energy of the BH air showers are due to the presence of invisible channels (gravitons, neutrinos, and non-decaying τ s). For example, the proton remnant may fragment into a top quark, which decays as t → b + W . If the W decay mode is leptonic, a consistent part of the initial proton remnant energy may be carried away by a neutrino. On the contrary, most of the c.m. energy in the ν e -CC interaction is observable. This leads to a stable air shower total energy.
The differences between SM and BH air showers can be quantified by choosing a benchmark model for the BH process and comparing this model to the SM process. This method also allows to differentiate the effects of various parameters and theoretical models from the stable characteristics of the BH air showers. With some guidance from the theory, a reasonable choice is: The left panel of Fig. 1 shows the ν e -CC air showers and the benchmark BH air showers (50 runs each, neutrino primary energy E ν = 10 7 TeV). The difference in the shower maxima is ∆X m = 214 ± 12 g cm −2 . Although the showers appear to be quite distinct, this difference is a consequence of the same choice of X 0 for both BH and SM air showers. The right panel of the figure shows ν e -CC air showers shifted so X m (BH) ≃ X m (ν e -CC). Therefore, BH and SM air showers can only be distinguished when X 0 − X m is clearly measured. Since present detectors cannot measure X 0 , ν e -CC air showers and BH air showers cannot be discriminated on an event-by-event basis [9]. We will see in Sect. VIII that this conclusion does not substantially change if the BH parameters are varied.

VII. SIMULATION RESULTS: SM VS. BH AIR SHOWER DETECTION IN GROUND AND HYBRID DETECTORS
Discrimination of BH air showers and ν e -CC air showers can be improved by the use of ground arrays. The best possible scenario for BH air shower detection is a technique that combines air fluorescence telescopes and a ground array. Since fluorescence telescopes are able to measure X m accurately, a good air shower discriminator is to fix X 0 and count the number of particles at various distances from X m . This is equivalent to fixing the detection level and varying the air shower first interaction point. The fluctuation due to the change of X 0 is negligible compared to fluctuations arising from other uncertainties. Figure 2 shows the number of muons at different atmospheric depths vs. the number of electrons at the air shower maximum for BH air showers (benchmark model) and ν e -CC air showers. The muons are measured from the ground array and the electrons are measured from the fluorescence telescopes. The BH air showers are characterized by a higher muon content than the ν e -CC air showers. Although the separation is not large enough to distinguish the air showers on an event-by-event basis, discrimination of BH and SM events is possible with enough statistics. The number of muons also depends on the detection level. Since different ∆Xs can be seen as different initial interaction points, the air showers in the lower right (upper left) panel of Fig. 2 can be  If only ground detection is possible, a good BH vs. SM discriminator is the number of electrons and muons at various atmospheric depths X m + ∆X (Fig. 3). As in the hybrid detection scenario, although the BH air showers show higher muon content than the ν e -CC air showers, their discrimination requires large statistics. The number of muons depend on the observation level and decreases as ∆X increases. In absence of an air fluorescence telescope to accurately measure X m , alternative techniques must be used to reconstruct the air shower maximum.

VIII. SIMULATION RESULTS: EFFECT OF BH PARAMETERS
The effects of various BH models can be studied by varying the parameters described in Sect. III and comparing the simulations to the benchmark model. The air shower longitudinal developments for different choices are shown in Figs. 4-8 and summarized in Table II. The depth of the air shower maximum is not significantly affected by changes in the BH parameters or the theoretical model. Fluctuations in the air showers are large; all X m values lie within one standard deviation from each other, with the possible exception of the simulations with nonzero minimum length. We conclude that the main characteristics of BH air showers described in the previous section are robust. This result can be qualitatively explained by noticing that most of the c.m. energy of the neutrino-nucleon collision is not trapped in the BH. Therefore, different choices of BH parameters do not produce large observable effects in the air shower development. The main factors determining the BH evolution (BH mass distribution, energy and spectrum of emitted quanta) are difficult to disentangle because their variations do not act to coherently increase or decrease the shower maximum. For instance, increasing the minimum-allowed BH mass from M min = 2 TeV to M min = 10 TeV increases the average BH mass in the air showers. This leads to a larger number of quanta. However, this property does not translate into a faster air shower   development; the average energy per quanta is smaller, and the two effects compensate each other. Two interesting facts are worth observing. Firstly, the benchmark case has the largest cross section and the largest X m . This result is mainly due to the choice of a relatively small fundamental Planck constant. Adding quantum effects or graviton loss at formation seems to decrease X m slightly. Secondly, the presence of a minimum length may possibly be the only BH physical signature distinguishable from the black disk model. However, in our simulations the choice of l min has been purposedly fine-tuned to the maximum-allowed value that allows BH formation with primary neutrino energy E ν = 10 7 TeV. Relaxing this choice leads to values of X m closer to the black disk result.

IX. SIMULATION RESULTS: BACKGROUND-FREE SIGNAL FROM τ DECAY
It has been suggested that τ leptons produced during BH events, either directly or through decay of other particles, can produce an observable signal similar to the double bang produced in ν τ -CC interactions [10]. The mechanism is the following. The BH evaporation initiates a first air shower, while a second air shower is initiated by the τ decay at a lower altitude. If the second bang is large enough to be observed, the double peak feature provides a background-free signature independent of the first interaction point.
The τ lepton in the ν τ -CC interaction carries on average about 80 % of the total c.m. energy. In a BH event, the fraction of the total BH mass going into the τ second bang can be estimated from Eq. (3) to be roughly Therefore, a τ from BH evaporation carries on average a few percent of the total BH mass. Since the energy trapped in the BH accounts generally only for a small fraction of the initial c.m. energy, the total energy in the τ channel is  on average less than 1% of the total air shower energy. This gives the τ a higher chance to decay before reaching ground, but a smaller energy deposit in the second bang, making the latter harder to detect. Moreover, τ s are not produced every time a BH is formed and the probability to observe an air shower with a τ decaying before reaching ground is relatively small.
The double bang signature can be studied by selecting only air showers containing at least one τ decaying in air. We simulated these air showers using the benchmark parameters, zenith angle 70 • and X 0 = 160 g cm −2 (altitude of the first interaction point = 20 km). The higher altitude gives a larger separation between the two bangs. Two possible scenarios were considered: i) τ s decaying at any depth in the atmosphere, and ii) at least one τ decaying at an altitude X > 0.75(X g − X 0 ), where X g is the slant depth of the ground ("low altitude τ s"). Case ii) represents the best possible scenario for double bang detection, as the second bang is expected to occur close to the detector. The longitudinal profiles of the low altitude τ air showers are shown in the left panel of Fig. 9. The second bang is visible around 2250 g cm −2 . The longitudinal development of the τ channel contribution to the air shower is plotted in the right panel. The peaks between 500 g cm −2 and 1500 g cm −2 are from events with multiple τ s, where at least one of these τ s decays at low altitude. As is expected, the τ component contributes only a minimal fraction to the overall air shower energy. Although the second bang is in principle detectable, it is within the fluctuation of the first bang's tail; current and next generation detectors cannot discern a few percent feature. Moreover, existing detectors such as the PAO have a limited field of view and cannot track the full profile.

X. CONCLUSIONS
We simulated extensive air showers initiated by TeV-scale BH events produced from neutrino interaction in Earth's atmosphere. These simulations were performed with an advanced fortran MC code [16] that includes most of the theoretical results of the recent literature. The BH air showers were compared to SM air showers and different models of BH formation and evolution were investigated. We also studied exotic signatures of BH events, such as the "τ double bang" signature. Our goal was to test various proposals of BH vs. SM air shower discrimination and look for new ways of differentiating models of BH formation and evolution.
Our results show that the main features of BH air showers are largely independent of the details of BH formation and evolution. Statistical fluctuations and limitations in detection techniques hinder the discrimination of alternative theoretical or phenomenological models. No difference between the black disk model and alternative models of BH formation, or between BHs with final explosive decay and stable remnant, can be detected with current UHECR observatories. Distinguishing SM and BH air showers with hybrid detectors is possible if enough statistics is gathered. The most promising way is to measure the air shower maximum with a fluorescence telescope and count muons at ground. The double bang signature cannot be observed by any realistic detector at the present stage. These results imply that the theory of TeV-scale BHs in UHECRs is robust, but BH air showers are hardly to probe details of "new physics" in the near future.