Photons coming from an opaque obstacle as a manifestation of heavy neutrino decays

Within the framework of physics beyond the standard model we study the possibility that mesons produced in the atmosphere by the cosmic ray flux, decay to heavy Majorana neutrino and these mostly to photons in the low mass region. We study the photon flux produced by sterile Majorana neutrinos ($N$) decaying after passing through a massive and opaque object such as a mountain. In order to model the production of $N$'s in the atmosphere and their decay to photons, we consider the interaction between the Majorana neutrinos and the standard matter as modeled by an effective theory. We then calculate the heavy neutrino flux originated by the decay of mesons in the atmosphere. The surviving photon flux, originated by $N$ decays, is calculated using transport equations that include the effects of Majorana neutrino production and decay.


I. INTRODUCTION
The neutrino sector has provided through the discovery of neutrino flavor oscillations the most compelling evidence for physics beyond the Standard Model (SM ). However, other mysteries related with the same sector are still open. In particular, the tiny ordinary neutrino masses problem, for which the seesaw mechanism stays as one of the most straightforward ideas for solving it [1][2][3][4][5][6]. This mechanism introduces right handed sterile neutrinos that, as they do not have distinct particle and antiparticle degrees of freedom, can have a Majorana mass term leading to the tiny known masses for the standard neutrinos, as long as the Yukawa couplings between the right handed Majorana neutrinos and the standard ones remain small. Even for the low masses range for N here considered, the simplest Type-I seesaw scenario leads to a negligible left-right neutrino mixing |U lN | 2 ∼ m ν /M N ∼ 10 −9 [7][8][9]. The mixing U lN weighs the couplings of N with the S.M. particles and in particular with charged leptons through the V − A interaction Thus, as suggested in [9], the detection of Majorana neutrinos (N ) would be a signal of physics beyond the minimal seesaw mechanism, and its interactions could be better described in a model independent approach based on an effective theory. We consider a simplified scenario with only one Majorana neutrino N and negligible mixing with the ν L . In addition, the effective operators here presented allow the N -decay to one neutrino plus one photon.
This decay channel could account, in a particular parameter region, for some neutrino related problems as the MiniBOONE and [10,11] and SHALON [12] anomaly. The Majorana neutrino effective phenomenology regarding the relevant N decay modes and interactions is treated in [13,14].
In the present work we study the possibility that mesons produced in the atmosphere by the cosmic ray, decay to heavy Majorana neutrino and these mostly to photons in the low mass region. We consider the scenario in which this radiative decay could be detected as a photon flux coming from an opaque obstacle such as mountain, which will stop the photons produced before the obstacle as well as any other photons originated by another mechanism, leaving an observable survival photon flux generated after the obstacle.
In Sec. II, we briefly describe the effective operator approach. In Sec.III, we present the production mechanism of heavy neutrinos (N ) by meson decay and the N decay to photons.
In Sec.IV, we discuss the bounds on the effective operators coming from different experiment as 0νββ decay, colliders results, meson decay, Super-K and astrophysical observations, in the mass range of tens of MeV. We calculate the number of events of photons to be observed by a Cherenkov telescope at different distances from the opaque obstacle. We show the results as a contour plot for the number of events and include the region allowed by the bounds.
We leave to the Appendix A the study of meson decay to N in the effective formalism and in the Appendix B we present the calculation for the N → νγ decay in the LAB frame.
Finally, in Sec.V, we present our conclusions.

II. EFFECTIVE MAJORANA INTERACTIONS
In this work, we study the observable effects of a heavy sterile Majorana neutrino N decaying to photons after passing through a massive and opaque object such as a mountain.
Thus, we need to model the interactions of N with ordinary matter in order to describe the N production by meson decay in the atmosphere and the subsequent N -decay to photons.
Being N a SM singlet, its only possible re-normalizable interactions with SM fields involve Yukawa couplings. But, as we discussed in the introduction, these couplings must be very small in order to accommodate the observed tiny ordinary ν masses. In this work, we take an alternative approach, considering that the sterile N interacts with the light neutrinos by higher dimension effective operators, and take this interaction to be dominant in comparison with the mixing through the Yukawa couplings. In this sense, we depart from the usual viewpoint in which the sterile neutrinos mixing with the standard neutrinos are assumed to govern the N production and decay mechanisms [15,16].
We parameterize the effects of new physics by a set of effective operators O constructed with the standard model and the Majorana neutrino fields, satisfying the SU (2) L ⊗ U (1) Y gauge symmetry [9]. The effect of these operators is suppressed by inverse powers of the new physics scale Λ, which is not necessarily related to the Majorana neutrino mass m N .
The total Lagrangian is organized as follows: grangian density for the interaction of right handed Majorana neutrinos N with bosons, leptons and quarks. We list the dimension 6 operators that can be generated at tree level or one-loop level in the unknown fundamental ultraviolet theory, and which are baryon-number conserving. The first subset includes operators with scalar and vector bosons (SVB), and a second subset includes the baryon-number conserving 4-fermion contact terms: where l i , u i , d i and L i , Q i denote, the right handed SU (2) singlets and the left-handed SU (2) doublets, respectively for the family i. The following one-loop level generated operators coefficients are naturally suppressed by a factor 1/16π 2 [9,17]: As we will show these operators contribute to N production by meson decay in the atmosphere and the subsequent N decay to photons.
In order to obtain the necessary interactions, we derive the relevant pieces of the effective Lagrangian terms involved in the calculations. We take the scalar doublet after spontaneous We have contributions to the effective Lagrangian coming from (3), related to the spontaneous symmetry breaking process: and the four-fermion interactions involving quarks and leptons from (4) where for simplicity in the notation we omit here the family index. The one-loop generated operators are suppressed by the 1/(16π 2 ) factor but, as we show in [13], these play a major role in the N -decay. In particular for the low m N range studied here and for the relative size of the coupling constants we will discus in the section IV, the dominant channel decay is N → νγ, which is produced by different terms coming from the operators in (5).
where P (a) is the 4-moment of the incoming a-particle. Moreover c W = cos(θ W ) and The constants α (i) L j , with j = 1, 4, are associated to the specific operators: The complete Lagrangian for the effective model is presented in an appendix in the recent work [14]. For completeness we include in Fig.1 a plot with the principal decay channel in the low mass region. A sum on particle and antiparticle is understood and the channel N → ν lep include the three body decays to leptons N → ν l i l + i l − i , ν l i l + j l − i , ν l i ν l iν l i In Sec. IV, we discuss the different bounds we consider on the coupling α (i) J and the strategy we follow to take it into account for the prediction on the fluxes.

III. PHOTON FLUX BY HEAVY NEUTRINO DECAY IN THE ATMOSPHERE.
We consider the transport equation for charged pions and kaons in the atmosphere including the interaction and decay terms [18]: where X is the slant path in g/cm 2 , ρ(X) is the density of the atmosphere, M represents the meson π or K and N a nucleon. The cosmic nucleon flux is parametrized as where γ = 1.7 and A 0 = 1.8 is the normalization constant for the initial cosmic flux. The constants c π and c K are m π /τ π and m K /τ K respectively, with τ π,K the mean lifetime of pion and kaon.
The constants Z i j are the spectrum-weighted moments for the inclusive cross section for a incident particle i colling with air nucleus and producing a outgoing particle j, with i, j = π, K, N . On the other hand the attenuation length constants Λ N , Λ π , Λ K are related With the usual approximations: (i) the hadron flux can be factorized φ N (E, X) = E −α φ N (X), (ii) the interaction length is independent of energy and (iii) the differential cross section is Feynman scaling, the solution for the meson flux is [18] The transport equation for the heavy neutrino N has as dominant contributions the absorption due to the N -decay and N -regeneration coming from the decay of mesons π ± and K ± .
We consider only the dominant meson decay process M → N µ, where M represents the meson π or K. The calculation of these decays in the effective theory we consider is shown in the Appendix A and the obtained result for the width is: where being q = d, s for π or K, respectively.
Continuing with the transport equations, we have for the Majorana neutrino N with an absorption term given by the N -decay and two source terms coming from the meson decay, and where with The branching ratios B r (π → N µ) and B r (K → N µ) can be written as Finally, the expressions for the decay distributions are and and for the integrals on the meson flux z π max z π min dzφ π (E/z, X) and insert the corresponding meson flux obtained in Eq.13 with M labeling the meson π and K. The integration in the z-variable is direct: where the function H M (E, u, X) reads With the above definitions, we have the heavy neutrino flux given by where We show in Fig.(2)  for the couplings intensity discussed in the next section. For high energy, the fluxes are independent of the value of m N , while for lower energy heavy neutrinos present a lower flux due to a shorter decay time.
As can be seen from Fig. 1 (see also [14]), for the masses of N considered (tens of MeVs), the dominant decay channel is N → γν. This channel decay was calculated in [14] Γ Thus, the total width for the low mass region is In order to study the production of photons by the heavy neutrino decays, we consider the coupled transport equations into the first equation and solving, we have the solution Thus, if we call l the distance between the obstacle and the detector and l i the traveled distance inside it, then we can write the photon flux arriving the detector as: where we have removed the photons produced inside the obstacle because they get absorbed.
In the next section we present our numerical results. We will consider the different bounds on the effective operators and the predictions for the photon flux, including the number of events to be detected.

IV. NUMERICAL RESULTS
The heavy Majorana neutrino couples to the three flavors families with couplings proportional to α (i) J /Λ 2 . For the case of the operator O (6),i N eφ these couplings can be related to the mixing angle between light and heavy neutrinos U l i N in Eq.1 [9] In analogy we can use the combination α (i) J v 2 /(2Λ 2 ) to represent the coupling intensity for all the operators. As it was discussed [13,14,19,20], the most restrictive bound on the operators O we have |U | 2 ≤ 8.8 × 10 −5 for the BELLE bound.
It is clear the different size between the contributions of both kind of operators. We maintain this hierarchy throughout the work decoupling the operators contributing to 0νββ.
for the first family. These curves were obtained for generic couplings compatible with 0νββ and the BELLE bound. The idea here is to show the relative behavior of these fluxes for different masses and distances to the obstacle. As we will discuss, there are most restrictive bounds that we will take into account when we show the final flux of photons.
In the considered mass range of tens of MeVs, the main sources of experimental bounds on the effective coupling α (i) J are the pion decay [7], the beam dump experiments [22], astrophysical observations, and the non-observation in Super-Kamiokande [23,24] of an excess of events coming from the decay of heavy neutrinos produced in the atmosphere. Another independent constraint on the effective operator coupling can be set based on the non-observation of heavy neutrino decays by the Super-K experiment. The heavy neutrinos produced by meson decays into the atmosphere would generate an excess of events in the detector. In order to estimate the importance of such effect, we calculate the fraction of neutrinos N that arrive at the Earth's surface and could decay inside the detector, which is located one kilometer deep (L deep ) and has a forty meters long edge (L e ). The flight distance is a function of the coupling and the N mass, L decay = L decay (m N , U 2 ). We calculate the fraction η(m N , U 2 ) as the ratio between the number of heavy neutrinos decaying inside the detector and the number of heavy neutrinos arriving at the Earth surface: where N (N ) and N (N ) We have considered the data reported in [23] and the discussion in [24]. We have found that a factor η = 10 −3 is a conservative value to impose the expected decay rates inside the detector does not exceed the rate of events detected by Super-Kamiokande [23] experiment.
In the plots of Figs.6a and 6b, we show the curve for which η = 10 −3 , which is a strong suppression factor. As we will see shortly there are regions of the parameters space where we still have an appreciable number of events and less one in a thousand heavy neutrinos N arriving the Earth decay inside the detector. In the same figure we include the upper limit for the coupling as obtained from the pion decay [7].
One further comment is in order at this point. In the Appendix A, we show the expression for the meson decay in the context of the effective theory we are studying. In this expression, we can see a strong contribution from scalar operators due to the light quarks masses in the corresponding denominators. In order to simplify the discussion, we will consider all the constants α J to be equal, but we have to take into account this important factor that determines the relative importance between the scalar and vectorial operators. For the pion decay, the coupling of scalar operators α scalar → α S 2 , α S 3 are accompanied by the big mass ratio mπ (mu+m d ) . Thus, is convenient to use the combination (m π /(m u + m d )) × α scalar to compare with the experimental bound. If we call α bound the corresponding bound, the value for α scalar to use in the production of N by meson decay is α scalar = ((m u +m d )/m π )×α bound .
In the case of N production by pion decay, we replace (m π /(m u + m d )) × α scalar → α bound .
But in the case of N production by K-decay, the replacement is We present the results for number of events that could be detected by a Cherenkov telescope through the observation of the electromagnetic showers originated by N decays after they have traversed an obstacle. In Figs.(6a) and (6b), we show the results as a contour plot for different number of events in the plane (m N , U 2) for different distances to the obstacle. We consider a generic detector like SHALON [12] with an effective area of 10 In the same figures, we show the bounds coming from pion decay and BELLE experiment [21]. Moreover we include the curves for which the N lifetime is the decay of mesons π and K, generically M decay. From the different Lagrangian presented in this section, we consider here the relevant pieces for the considered decay. Thus, we have We deduce the decay amplitude with q = d, s for π and K decay, respectively. In the last term, we need to rearrange the field operators in order to put together quarks fields in a sandwich and the lepton fields in another. In order to do that, we make a Fierz transformation to the last term taking into account a minus sign from the permutation of fermions, and then we have The calculation of the leptonic matrix element is straightforward, In order to calculate the hadronic matrix element, we have to rely on the symmetries. The matrix element 0|ūγ ν γ 5 q| M is a Lorentz 4-vector because the meson M is pseudoscalar andūγ ν γ 5 q is a pseudo 4-vector. The meson state is described by its four momentum q µ and nothing else, since the pion has spin zero. Therefore, q µ is the only 4-moment on which the matrix element depends and it must be proportional to q µ . Thus, we can write On the other hand, for the same reason, the matrix element of the 4-vector is In the case of the matrix element of the scalar or pseudo-scalar, we have to use the equation of motion where m q = m d , m s for the decay of π and K, respectively.
Putting it all together and integrating over the 2-body phase space, we obtain with M given in Eq.A3.
The result is We follow the development shown for µ-decay in the book of T. K.Gaisser [25], but in our case for the N → γν decay (we adapt the calculations presented in the Appendix of the recent work [26]). First, we obtain the N decay width in its rest frame, and then boost the result to the Laboratory frame. In the N rest frame, we have the following expression: being θ ν the direction of motion of the final ν taken from the Majorana neutrino N moving direction, and P = cos θ P where θ P is the angle between the Majorana neutrino spin direction in its rest frame, and its moving direction as seen from the Laboratory frame. The variable x represents the quotient between the final neutrino energy in the rest frame of the N and the mass of the Majorana neutrino: x = k 0 /m N . The functions F 0 (x) and F 1 (x) are To obtain the corresponding expression in the laboratory frame, we make the appropriate Lorentz transformations. Denoting by E ν and E N the Laboratory energies of the final neutrino and the Majorana neutrino, respectively, we have with z = E ν /E N and β N = 1 − m N 2 /E 2 N 1.
We implement the Lorentz transformation with the help of the δ-function, yielding We first integrate over θ ν and then we integrate over x in the interval (x min , x max ) with x min = z/(1 + β N ) and x max = min(1, z/(1 − β N )), obtaining For the low mass range considered in this work, the clearly dominant decay channel is the neutrino plus photon mode, and Γ tot LAB (E) = i=e,µ,τ Γ N →ν i γ LAB (E) + Γ N →ν i γ LAB (E) . Then we consider the γ decay channel, leading to the final γ photon distribution in the laboratory frame: Thus, after the indicated integrations in the evolution equations, the useful expression that we obtain is where z = 1 − y, x 0 = 1/2 and P = +1 for the right-handed Majorana neutrinos.