Probing Active-Sterile Neutrino Transition Magnetic Moments with Photon Emission from CE$\nu$NS

In the presence of transition magnetic moments between active and sterile neutrinos, the search for a Primakoff upscattering process at coherent elastic neutrino-nucleus scattering (CE$\nu$NS) experiments can provide stringent constraints on the neutrino magnetic moment. We show that a radiative upscattering process with an emitted photon in the final state can induce a novel coincidence signal at CE$\nu$NS experiments that can also probe neutrino transition magnetic moments beyond existing limits. Furthermore, the differential distributions for such a radiative mode can also potentially be sensitive to the Dirac vs. Majorana nature of the sterile state mediating the process. This can provide valuable insights into the nature and mass generation mechanism of the light active neutrinos.


I. INTRODUCTION
Coherent elastic neutrino-nucleus scattering (CEνNS) [1] was first observed in 2017 by the COHERENT collaboration with a statistical significance of 6.7σ [2]. Future CEνNS experiments with extremely low nuclear recoil thresholds now aim to detect neutrinonucleus scattering events with O(eV) momentum transfers [3,4]. These next-generation experiments will not only test the Standard Model (SM) CEνNS rate more precisely but also constrain physics beyond the SM. One such new physics (NP) scenario is the socalled Primakoff upscattering of light active neutrinos to heavy sterile neutrinos via a transition magnetic moment. This dipole portal could be a promising way to probe both the existence of sterile neutrinos and possible NP generating the dipole coupling.
The focus of the present work is a closely-related process that can also produce distinct signatures at CEνNS experiments. This is the upscattering of an incoming active neutrino to a sterile neutrino, which subsequently decays to an active neutrino and photon, see Fig. 1. Searches for such a radiative upscattering process have been suggested for the DUNE [5], IceCube [6], and Super-Kamiokande [7] experiments. This process can be used to probe transition magnetic moments in CEνNS experiments as well as to distinguish sterile neutrinos of Dirac or Majorana nature, indicating the corresponding nature of the active neutrinos [8][9][10][11][12][13][14][15][16].

II. NEUTRINO TRANSITION MAGNETIC MOMENTS
A transition magnetic moment between the three active neutrinos and a sterile neutrino (or gauge-singlet fermion) can be described by the effective Lagrangian, where F µν = ∂ µ A ν − ∂ ν A µ is the electromagnetic field strength tensor, ν αL is an active neutrino field of flavour α = {e, µ, τ }, N is the sterile neutrino field, and µ α νN are the active-sterile dipole couplings. Here, we introduce N as an additional light (m N m N ) sterile state with a sterile-sterile transition magnetic moment to N (with a transition dipole coupling µ N N ). The Lagrangian in Eq. (1) is only valid at energies below the electroweak (EW) scale because it is not invariant under the SM gauge symmetry. Above the EW scale, it must be matched onto operators at dimension-six and above containing the SU(2) L and U(1) Y gauge fields, SM Higgs doublet H, and lepton doublet L α . Coherent scattering processes take place at energies well below the EW scale and hence Eq. (1) remains applicable.
The active and sterile neutrinos in Eq. (1) can either be Dirac or Majorana fermions.
For example, the active neutrinos ν αL can either be the left-handed Weyl components of Dirac (ν α = ν αL + ν αR ) or Majorana (ν α = ν αL + ν c αL ) fields, where in the former case it is necessary to introduce additional (sterile) Weyl fields ν αR . Similarly, the sterile neutrino N (and N ) can either be composed of independent left-and right-handed Weyl fields (N = N L + N R ), or a single right-handed field (N = N c R + N R ). In principle, there are four possible combinations of Dirac and Majorana active and sterile neutrinos. However, it can be shown that for energy scales much greater than the active neutrino masses, E ν m ν , the rates for processes involving Dirac and Majorana active neutrinos are approximately equal, in accordance with the Dirac-Majorana confusion theorem [17,18].
In this work, we will consider sterile neutrinos with masses similar to the energy scale of the process, E N ∼ m N . Consequently, the difference between the rates for Dirac or Majorana sterile neutrinos can be significant.
The active-sterile dipole couplings µ α νN in Eq. (1) give rise to the Primakoff upscattering process ν α A → N A; see, e.g., Refs. [13,19,20]. For relativistic ν α , it can be shown that the differential cross section in the nuclear recoil energy E R for this process is the same for outgoing Dirac and Majorana N . The non-observation of deviations from the SM CEνNS nuclear recoil rate at experiments can therefore be used to constrain the dipole couplings µ α νN as a function of the sterile neutrino mass m N . Constraints on the dipole couplings µ α νN have been set by a variety of experiments (c.f. Fig. 3) and are generally flavour dependent [5].
The active-sterile dipole couplings µ α νN can also induce the radiative upscattering process ν α A → ν β Aγ, depicted in Fig. 1. The rate of the process is proportional to |µ α νN µ β νN | 2 and is therefore suppressed with respect to the Primakoff upscattering. However, the outgoing photon serves as an additional signal to the nuclear recoil and provides kinematical information that can be used to discriminate between Dirac and Majorana N . We would like to remark that because the outgoing neutrino is not detected, the actual process is ν α A → XAγ, where X may be an active neutrino ν β or a sterile state N . The rate is then generically proportional to |µ α νN X µ XN | 2 , where the sum is over all active and sterile neutrinos that are coupled to N via a transition magnetic moment.
In this work we will examine as a case study the upcoming NUCLEUS experiment located at the Chooz reactor facility [3,4]. The planned experimental arrangement and future upgrade make it feasible to detect the energies and angles of outgoing photons [21].
The (non)-detection of photons at NUCLEUS can thus constrain the dipole couplings We give the detailed calculation of the differential cross section for ν α A → XAγ in Appendix B. Here, we outline the main results and the salient features relevant to our study. The amplitudes for the ν α A → XAγ process are given in the scenarios where the sterile neutrino N is Dirac or Majorana, respectively, by where F µνρσ , given in Appendix B, contains the denominator of the sterile neutrino propagator, the hadronic current of the nucleus A, the propagator of the exchanged photon, and the four-momentum and polarisation of the outgoing photon. We have taken the dipole couplings to be imaginary, (µ XN ) * = −µ XN , and thus CP conserving if X and N have opposite CP phases [22]. In the Dirac case, N is created by the hermitian conjugate of the first term and annihilated by the second term in Eq. (1). In the Majorana case, both the second term in Eq. (1) and its hermitian conjugate annihilate N ; it is then possible to make the replacement P R → (P R + P L ) = 1 at the decay vertex. The P R and P L project out the momentum / p N and mass m N terms from the N propagator, respectively. The three-body phase space of the final state XAγ is described by four variables; we choose the nuclear recoil energy E R , photon energy E γ , nuclear recoil angle θ R , and photon angle θ γ (both angles defined with respect to the incoming neutrino direction).
For a CEνNS experiment to detect an outgoing photon, and therefore the process ν α A → XAγ, a radiative sterile neutrino decay N → Xγ must take place within the detector. The probability for this to occur is given by the N → Xγ branching ratio B N →Xγ = Γ N →Xγ /Γ N multiplied by the probability for a decay to take place within the detector, P det where L det is the detector length, Γ N is the total decay width of N , and the boost factors are given by βγ = γ 2 − 1 and γ = E N /m N . It is also possible that N decays via invisible channels (for example, to a light dark sector), which contribute to the total width as Γ N = Γ N →Xγ + Γ inv N , giving as the probability of observing a radiative decay inside the detector. For example, for m N ∼ 1 MeV, the upper limit on the active-sterile mixing squared from beta decays is Assuming Γ N = Γ N →Xγ + Γ inv N to be small and therefore the total N decay length to be much longer than the detector length, N = βγτ N = βγ/Γ N L det , the exponential in the decay probability P det N can be Taylor expanded to give B N →Xγ P det N = Γ N →Xγ L det /(βγ), which is independent of the total decay width. For a detector length of L det = 5 cm, the rate of ν α A → XAγ events is independent of an invisible N decay width up to Γ inv N ∼ 10 −11 MeV. Above this value, P det N ≈ 1 and the probability that ν α A → XAγ occurs within the detector is suppressed by the small branching ratio.
We therefore consider three benchmark scenarios in this work. The first is to assume that N → ν β γ is the only N decay channel. The branching ratio is B N →Xγ = 1 and the decay probability, P det The branching ratio is now suppressed, B N →Xγ ≈ Γ N →Xγ L det /(βγ) 1, while the decay probability is large, P det N ≈ 0.63. However, the product B N →Xγ P det N is roughly the same as in the previous scenario. The third case is to include an additional light sterile state N (m N m N ) and a non-zero sterile-sterile dipole coupling µ N N , but no invisible modes. Interestingly, the rate for N → N γ satisfies Γ N →N γ βγ/L det for µ N N 10 −4 µ B and L det = 5 cm, while the upper limit on the sterile-sterile transition dipole coupling from invisible vector [23]. However, stronger bounds of order µ N N 10 −6 µ B can be derived from LEP [13]. For µ N N ∼ 10 −6 µ B , the decay probability is much larger than for µ β νN ∼ 10 −8 µ B , increasing the expected rate for the radiative signal. We want to emphasise that in this scenario the radiative upscattering ν α A → XAγ can compete with Primakoff upscattering ν α A → N A in constraining the active-sterile dipole coupling.
Regardless of the benchmark choice above, when the sterile neutrino total decay width Γ N is much smaller than the sterile neutrino mass m N , it is possible to use the narrow width approximation (NWA). In this limit, the ν α A → XAγ differential cross section in E R can be decomposed in the NWA as the ν α A → N A cross section multiplied by the N → Xγ branching ratio, where X = {ν β , N , . . .}. The differential cross section in E R therefore has the same shape as for ν α A → N A but is suppressed by an additional factor of |µ XN | 2 . To yield a differential rate, Eq. (5) must be multiplied by the incoming neutrino flux dφν α dEν and decay probability P det N and integrated over the incoming neutrino energy. Majorana nature of N , as an overall factor of two can absorbed into the measured value of µ XN .
The ν α A → XAγ differential cross sections in the photon energy E γ and angle θ γ can also be computed in the NWA, but cannot be factorised as in Eq. (5). The double differential cross section in these two variables can be written as where L γ,D(M) µν H µν is the contraction of leptonic and hadronic currents for the process and ∆ 4 is the 4 × 4 symmetric Gram determinant formed from any four of the five incoming or outgoing four-momenta, both defined in Appendix B. Here, the variable t 1 = q 2 = −2m A E R corresponds to the four-momentum squared of the exchanged photon, s 1 = p 2 N to the four-momentum squared of the sterile neutrino, and F(t 1 ) is the nuclear form factor. The limits of integration t ± 1 are found by solving ∆ 4 = 0 for t 1 . In general, the integral cannot be performed analytically due to the presence of F(t 1 ).
In Fig. 2, we plot the double differential cross section of the ν α A → ν α Aγ process (for a 73 Ge target) as a function of the outgoing photon energy E γ and angle θ γ in the Dirac (left) and Majorana (right) cases, setting F(t 1 ) = 1 for simplicity and choosing the benchmark values E ν = 3 MeV (approximately the peak of the reactor neutrino flux at Chooz) and m N = 1 MeV (so that N can be produced on-shell). In the plots, we choose the values µ α νN = 3 × 10 −8 µ B and Γ inv N ≈ Γ N = 10 −11 MeV. The difference between the Dirac and Majorana scenarios is striking; while both cases predict a considerable number of events at small photon energies, irrespective of the photon emission angle, there are more forward emissions of high energy photons in the Majorana case.
The single differential distributions in E γ and θ γ are also different in the Dirac and Majorana scenarios. In the Dirac case the differential cross section in E γ decreases linearly from the minimum to maximum photon energies E − γ and E + γ , respectively, while in the Majorana case the distribution is flat, i.e.
, and Θ(x) is the Heaviside step function. Furthermore, the angular distribution in the lab frame peaks at slightly lower angles in the Majorana case compared to the Dirac case.
The lab-frame distributions are in exact agreement with the angular distributions in the rest frame of N , which can be derived purely from arguments of rotational and charge, parity and time reversal (CPT) invariance. Due to the conservation of angular momentum, Dirac N (N ) can only decay to left-polarised (right-polarised) photons γ − (γ + ) with an angular distribution in cos ϑ γ proportional to (1 + cos ϑ γ ) in the rest frame [10,[24][25][26][27][28][29].
Majorana N can decay equally to both left and right-polarised photons with angular distributions in cos ϑ γ proportional to (1 + cos ϑ γ ) and (1 − cos ϑ γ ), respectively; the total distribution is thus isotropic. The distinctive Dirac and Majorana energy distributions in Eqs. (7) and (8), respectively, can be readily derived by boosting these rest-frame angular distributions to the lab frame. The photon circular polarisation provides an additional handle on the nature of the sterile neutrino and the CP properties of the transition dipole coupling (see, e.g., Ref. [27]); a detailed analysis of the complementarity of a polarisation study is beyond the scope of this work and will be addressed elsewhere [30].

IV. SENSITIVITY TO ACTIVE-STERILE NEUTRINO MAGNETIC MOMENTS
In order to study the feasibility of our proposal, we will now examine the NUCLEUS experiment, which aims to detect CEνNS (ν e A →ν e A) with nuclear recoil energies as low as E R ∼ 10 eV. Situated at the very-near-site (VNS) of the Chooz reactor site, the detector will receive an electron antineutrino flux of φν e ∼ 10 12ν e cm −2 s −1 . In the near future, phase I of the experiment will use a 10 g Al 2 O 3 /CaWO 4 detector of size L det ∼ 5 cm, while in the far future, phase II will upgrade to a 1 kg 73 Ge detector of size L det ∼ 25 cm [21].
Interestingly, the detector will be sensitive to nuclear recoils and potentially to outgoing photons in the 1 keV to 10 MeV energy range. These can be detected at the cryogenic outer veto with an ionisation resolution of 50 -100 keV for O(MeV) photons [21]. A coincidence study between the nuclear recoil and an outgoing photon signal which can potentially lead to excellent background rejection. In this work, we therefore neglect any secondary backgrounds to the radiative mode as a first approximation.
For the Chooz reactor antineutrino flux, the maximum number of events are expected for sterile neutrino masses m N ∼ 1 − 5 MeV. We find that the Dirac and Majorana cases show different differential rates in the photon energy, as expected from the differential cross sections shown in Fig. 2. The Majorana case presents a more symmetric distribution in the outgoing photon energy, while the Dirac case differential rate is shifted towards lower outgoing photon energies. For a specific sterile neutrino mass m N , the cross section for the radiative upscattering processν e A →ν e Aγ grows with increasing incoming neutrino energy, with a resonant energy around m N . Consequently, a different sterile neutrino mass can be probed for each neutrino energy in the Chooz spectrum. Average nuclear recoil and photon energies are increased for larger masses m N . The total number of events falls off at low energies due to a smaller cross section, and at high energies due to the decreased flux. Since the number of events peaks at E ν ∼ m N , the decay width can be considered to be roughly constant. In this work, we only consider reactor antineutrinos as a proof-of-concept, though in principle a similar study could be performed for solar, atmospheric, or beam-dump neutrinos.
In order to identify the reach of the NUCLEUS experiment within the dipole portal parameter space, we consider both the Primakoff and radiative upscattering processes at the NUCLEUS detector. For the former process, the experiment will observe N obs nuclear recoils over its run time T . If the expected number of recoil events is given by The nuclear recoil background has been estimated by the NUCLEUS collaboration to be N bkg = 100 keV −1 kg −1 day −1 [4]. The number of recoil events induced by SM CEνNS, N νA , and Primakoff upscattering, N N A , are found by integrating the cross sections of these processes over the incoming energies of the Choozν e flux, from the minimum to maximum nuclear recoil energies, and multiplying by the detector mass and run time.
We now assume that the experiment does not see an excess of recoil events over the SM CEνNS signal and background, i.e. N obs = N bkg + N νA , giving χ 2 = N 2 N A /N exp . For simplicity, we do not include a nuisance parameter in χ 2 to account for the presence of systematic errors. To set bounds at 90% C.L., we find the allowed values of µ e νN satisfying χ 2 < 2.71.
For the radiative upscattering process, the NUCLEUS experiment should in principle be able to detect N γ obs coincidence events of a nuclear recoil accompanied by an outgoing photon. We assume a negligible SM background and therefore take the expected number of coincidence events, N γ exp = N νAγ , to be Poisson distributed. Here, N νAγ is found by again integrating the radiative differential cross section in Eq. (5) over the Choozν e flux, from the minimum to maximum nuclear recoil energies, and multiplying by the detector mass and operation time. To take into account the probability of the decay occurring within the detector, we must also multiply by the factor P det . Assuming that NUCLEUS does not observe any coincidence events, N γ obs = 0, we set bounds at 90% C.L. by finding the allowed values of µ e νN satisfying N γ exp < 2.30. In Fig. 3, we summarise the current and future constraints on the electron-flavour active-sterile dipole coupling µ e νN as a function of the sterile neutrino mass m N . We show current bounds from terrestrial experiments [5-7, 13, 15, 31-34] and astrophysical processes [6,13,15] as solid lines (with excluded areas filled) and the expected sensitivities As discussed previously, the sterile-sterile dipole coupling µ N N is not subject to the same constraints as the active-sterile couplings µ α νN (for example, those on µ e νN in Fig. 3). Consequently, the rate for the radiative upscattering process can be increased with respect to the Primakoff upscattering, leading to an improvement in sensitivity. For m N 1 MeV, the near-and far-future sensitivities now constrain smaller values of µ e νN compared to the bounds from Primakoff upscattering.

V. CONCLUSIONS
For the radiative upscattering mode, a signal would consist of the coincidence of a nuclear recoil and an outgoing photon, separated by the decay length N of the sterile state. In this work, we have proposed a novel approach in which the final state photons are searched for in a separate detector to the CEνNS target, and in particular we have studied the detection prospects of the reactor-based NUCLEUS experiment. As seen in Fig. 3, the current limits on the transition dipole coupling µ e νN derived from Primakoff upscattering at the COHERENT experiment are stringent, almost coinciding with sensitivity of the NUCLEUS experiment 1 kg upgrade in the radiative upscattering mode.
Using Primakoff upscattering, the 10 g NUCLEUS experiment will improve the sensitivity (red thin dashed), extending the limits to the region excluded by astrophysical observations for sterile neutrino masses m N 10 MeV. If the 10 g NUCLEUS experiment is indeed able to detect the Primakoff upscattering, it will provide an exciting motivation to search for the radiative upscattering mode in the 1 kg NUCLEUS upgrade.
Even though the radiative mode is doubly suppressed by the dipole coupling as |µ e νN | 4 , its unique (potentially background-free) signature suggests that a future detection is not outside the realm of possibility. Such an observation would act as a smoking gun for our specific scenario and allow to differentiate it from other mechanisms.
The reach of the radiative upscattering mode can also be improved, as the intermediate This would have further implications on the neutrino mass generation mechanism, as well as theories of leptogenesis explaining the asymmetry between matter and antimatter in the universe. Specifically, if the sterile neutrino is found to be Majorana, the active neutrinos will also be of Majorana nature due to a mass term induced by the transition magnetic moment. We refer the interested reader to Appendix A, which provides a brief review of theoretical models in which transition magnetic moments are generated. observable radiative upscattering rate together with smallness of active neutrino masses will hint towards a neutrino mass model with enhanced symmetry, e.g. a horizontal symmetry reinforcing the Voloshin mechanism [8] or an inverse seesaw mechanism [36][37][38]. Below we provide a brief overview of the relevant constraints and implications for Dirac and Majorana neutrino mass models.
Dirac neutrinos: If the sterile state N R is a Weyl field with a coupling to the SM neutrino ν L via the dipole interaction in Eq. (1), then the active-sterile transition magnetic moments µ νN can give rise to mass terms of the form L ⊃ m νNνL N R +h.c.. This becomes more apparent when looked at from an effective field theory point of view as discussed in Ref. [11]. An effective Lagrangian can be constructed as where the d ≥ 4 denotes the operator dimension, j runs over all independent operators of a given dimension, µ is the renormalisation scale, and Λ corresponds to the new physics (NP) scale where new heavy degrees of freedom are integrated out.
In a SM gauge group invariant theory an effective transition magnetic dipole moment can be generated by gauge-invariant, dimension-six (d = 6) operators with couplings to the SU(2) L and U(1) Y gauge fields W a µ and B µ . Above the EW symmetry breaking scale, due to renormalisation group (RG) running these operators will mix with other d = 6 operators that contain the SM lepton doublet L = (ν L , e L ) T , N R , and the SM Higgs field H. One such operator actually leads to generation of a Dirac neutrino mass term m νN after the EW symmetry breaking, as can be seen by considering the basis of independent operators at d = 6 that are closed under renormalisation [11] O (6) are the U(1) Y and SU(2) L field strength tensors, respectively, g 1 and g 2 are the corresponding gauge couplings, andH = iσ 2 H * . Starting from the Wilson coefficients C such that C 1,2 (µ = Λ) [11]. After the EW symmetry breaking the combination C where µ = v corresponds to the EW symmetry breaking scale. On the other hand O yielding a relation between δm νN and µ νN given by which leads to the constraint for Λ = 1 TeV. This implies that for a realistic active neutrino Dirac mass m ν 1 eV, However, we note that this constraint is strictly applicable when N R is a Weyl field forming a Dirac pair with ν L . The above constraint does not hold true in a number of general circumstances. One example is the scenario in which N is a Dirac fermion containing two Weyl fields (i.e. N = N R + N L ) and has a Dirac mass m N m obs ν , which can a priori be completely decoupled from the generation of the active neutrino masses. This is the Dirac scenario we consider in this work, and therefore the constraint in Eq. (A6) can be safely neglected. Another example is the scenario in which the tree-level Dirac mass term m SM νN is of the same order but of opposite sign compared to the δm νN generated by µ νN . While some might consider this cancellation to be unnatural, it is by no means improbable. A third example is a Dirac seesaw scenario where N is a Dirac fermion However, one can always tune the tree-level Yukawa coupling L Y ⊃ y νLLH N R + h.c. such that the resulting tree-level contribution to the Dirac mass, y ν v/ √ 2, nearly cancels the loop-induced contribution. In this way, one can lower the right-handed neutrino masses m N to as low as MeV scales while keeping the active neutrino masses at m ν 1 eV.
Exceptions to the relation in Eq. (A7) can naturally occur when additional symmetries are present. An example is the so-called Voloshin mechanism, where an approximate global SU (2) H symmetry is introduced such that (ν c L , N R ) transforms as a doublet under the SU (2) H [8]. While this symmetry allows for a SU (2) H singlet transition magnetic moment term of the formN R σ µν ν L −ν c L σ µν N c R , it naturally forbids the SU (2) H triplet contribution to the neutrino mass termN R ν L +ν c L N c R . Instances of recent models employing the Voloshin mechanism in the context of neutrino magnetic dipole moments can be found in Refs. [15,16].
In the presence of a sizeable mass mixing between active and sterile states (as can be the case in a typical type I seesaw) an active-sterile transition magnetic dipole moment can also be induced through loop diagrams involving charged leptons [9,10]. Such an active-sterile mass mixing induced contribution to the transition magnetic dipole moment is given by [13] On the other hand, in the presence of a transition magnetic moment between ν L and N , a loop contribution to the light active neutrino masses is induced through Fig. 4, which is directly proportional to the Majorana mass of N , where Λ is the cutoff scale for the UV completion of the model. For m N ∼ 1 MeV, Λ ∼ 1 TeV and m ν 1 eV, Eq. (A9) leads to |µ νN | µ B < 10 −8 . While a larger Majorana mass of N would lead to relaxation of the tight constraint from Eq. (A8), it would make the constraint from Eq. (A9) more stringent. However, one can easily circumvent these constraints by considering a scenario with a quasi-Dirac N , such as in the case of the inverse seesaw mechanism [36][37][38], where an approximate lepton number conservation (due to a very small Majorana mass splitting of a predominantly Dirac N pair) makes the active neutrinos very light, while the N being pseudo-Dirac can be significantly heavier.

Appendix B: Calculation of the Radiative Upscattering Process
In this appendix we derive the differential cross section for the radiative upscattering process ν α A → XAγ, X = {ν β , N , . . .}. For completeness, we also derive the differential cross section for the Primakoff upscattering process ν α A → N A.
Firstly, we note that the active and sterile neutrinos ν α and N (N ) may either be Dirac or Majorana fermions. Rates for the Primakoff and radiative upscattering processes can be therefore computed for the four possible combinations of Dirac or Majorana ν α and N (N ). We will see that in the limit of massless neutrinos in the initial and final In both the Dirac and Majorana cases, the transition magnetic moment between the light active neutrinos ν α (or a sterile neutrino N ) and the sterile neutrino N is described by the effective Lagrangian where X = {ν αL , N , . . .} and µ XN = {µ α νN , µ N N , . . .}. In the case where both N and X are Majorana, the second term can be rewritten to give where we have defined (µ XN ) * = −µ XN and used the charge-conjugation properties CP T L C −1 = P L and Cσ T µν C −1 = −σ µν . If X is Majorana and N is Dirac (or vice versa) it is not possible to make such a simplification; it is nonetheless convenient to write in the first scenario (Majorana X, Dirac N ), while in the latter case (Dirac X, Majorana N ), In the following we will use these interaction terms to investigate the Primakoff and radiative upscattering processes for Dirac and Majorana ν α , N and N .
We finally note that if X and N are both Dirac, the following effective Lagrangian terms can also be written However, as we are considering purely left-handed (right-handed) incoming neutrinos (antineutrinos), these terms do not contribute to the Primakoff or radiative upscattering processes.
Primakoff upscattering: The active-sterile transition magnetic moment µ α νN induces the Primakoff upscattering process shown to the left in Fig. 5

. An incoming Dirac or
Majorana active neutrino ν α exchanges a photon with a target nucleus A and scatters to an outgoing Dirac or Majorana sterile neutrino N . As shown in the diagram, the ingoing nucleus and neutrino and outgoing nucleus and sterile neutrino have four-momenta k 1 , k 2 , p 1 and p N , respectively. The four-momentum exchanged by the photon is thus We first consider the scenario where N is a Dirac fermion, i.e. N = N L + N R . An incoming neutrino ν α , created by the SM charged current j µ W =ν αL γ µ αL , may be a Dirac or Majorana fermion. If ν α is Dirac, the neutrino is annihilated by the second term in Eq. (B1). This results in the following matrix element for ν α A → N A, where J A λ is the hadronic current of the nucleus. In the following we use the hadronic current J A λ = −ieZ(ū A γ λ u A )F(q 2 ), where F(q 2 ) is a nuclear form factor, which describes the electromagnetic interaction of a spin-1 2 nucleus. If ν α is Majorana, the neutrino created by the SM charged current j µ W =ν αL γ µ αL can instead be annihilated by the first or second term in Eq. (B3). The second term corresponds to the propagation of a negative helicity neutrino and is again described by the matrix element in Eq. (B6). The first term on the other hand corresponds to the propagation of a positive helicity neutrino; because the neutrino is ultrarelativistic, this process is helicity suppressed by the small ratio m ν /E ν .
If ν α is Dirac, the SM charged current j µ † W =¯ αL γ µ ν αL instead creates an antineutrinoν α , which is annihilated by the first term in Eq. (B1). The matrix element is If ν α is Majorana, the neutrino created by the charged current j µ † W =¯ αL γ µ ν αL can again be annihilated by the first or second term in Eq. (B3). The first term induces the process ν α A →N A with the matrix element in Eq. (B7). The second term induces the process ν α A → N A but is helicity suppressed by m ν /E ν .
We next consider the scenario where N is a Majorana fermion, i.e. N = N c R + N R . If the incoming neutrino ν α produced by the charged current j µ W =ν αL γ µ αL is Dirac, it can only be annihilated by the second term in Eq. (B4). This induces the process ν α A → N A with the matrix element identical to Eq. (B6). Similarly, an antineutrino created by the charged current j µ † W =¯ αL γ µ ν αL can only be annihilated by the first term in Eq. (B4). This induces the processν α A → N A with a matrix element identical to Eq. (B7). An incoming Majorana neutrino ν α created by j µ W =ν αL γ µ αL can be annihilated by both terms in (B2). However, the contribution to the process ν α A → N A from the first term is helicity suppressed. Likewise, an incoming Majorana neutrino ν α created by j µ † W =¯ αL γ µ ν αL can also be annihilated by both terms in Eq. (B2), but the contribution from the second term is suppressed. Consequently, the matrix elements for these processes are also given by Eqs. (B6) and (B7), respectively.
The differential cross section for the process can now be found by taking the absolute square of the scattering amplitude, averaging over the spin of the incoming nucleus and summing over the spins of the outgoing nucleus and sterile neutrino. Neglecting the mass of the incoming neutrino, the differential cross section is calculated as where s = (k 1 + k 2 ) 2 is a Mandelstam variable, m A is the mass of the nucleus, and dΦ 2 is the two-body phase space of the outgoing nucleus and sterile neutrino, .

(B9)
In the second equality we have used the Dirac delta function (enforcing the conservation of four-momentum) to integrate over four of the three-momentum components. We have also re-expressed the integral in terms of φ (the azimuthal angle defining rotations around the incoming neutrino direction) and the Mandelstam variable t = q 2 = (k 1 − p 1 ) 2 . The physical region of the phase space is defined by the inequality ∆ 3 < 0, where is the 3 × 3 symmetric Gram determinant constructed from any three of the four incoming or outgoing four-momenta [39]. The spin-averaged and summed squared matrix elements for the process ν α A → N A can be written, for both Dirac and Majorana N , as where the leptonic and hadronic components are given by As mentioned previously, the form of the hadronic current H µν in Eq. (B13) is technically only valid for spin-1 2 nuclei, while most targets considered for CEνNS are spin-0. However, for small recoil energies, which is the regime of interest for most CEνNS experiments, the spin of the nucleus has a negligible impact on the cross section. We can now insert Eqs. (B9) and (B11) into Eq. (B8) and express the matrix element squared in terms of the Mandelstam variables s and t. The matrix element squared does not depend on the azimuthal angle φ, which can therefore be integrated over. This gives the following differential cross section in the Lorentz-invariant variable t, where the function F (a, b, c been introduced for convenience.
We would now like to determine the cross section as a function of lab frame variables.
In the lab frame, where E ν is the incoming neutrino energy and E R is the nuclear recoil energy. In the lab frame the Mandelstam variables are s = m A (m A + 2E ν ) and t = −2m A E R . We transform the variable from t to E R by multiplying by the Jacobian factor ∂t/∂E R = −2m A , which gives the standard result In this work, we are interested in the limit in which the momentum exchanged by the photon (and therefore the nuclear recoil E R ) is much smaller than the mass of the nucleus m A . We also assume that the sterile neutrino mass m N is comparible to the incoming neutrino energy E ν , but much smaller than m A . In terms of the Mandelstam variables this corresponds to s ≈ m 2 A t, m 2 N , and the function in Eq. (B14) simplifies to F (s, t, m 2 . This is equivalent to dropping the terms proportional to 1/E ν and 1/m A in Eq. (B15).
Radiative upscattering: We now outline how to compute the matrix element squared and differential cross section of the radiative upscattering process ν α A → XAγ, i.e. the Primakoff upscattering ν α A → N A followed by the radiative decay N → Xγ, where X can either be an active neutrino ν β or another sterile state N .
Before doing so, it is useful to derive the decay rate for the process N → Xγ induced by the transition magnetic dipole moment µ XN . If N and X are Dirac fermions, sterile neutrinos N and antineutrinosN decay to X andX, respectively. The matrix element where p 3 and are the four-momentum and polarisation of the outgoing photon, respectively. The matrix element forN →Xγ is determined from the second term in Eq. (B1) and is given by To compute the decay rate, we multiply the spin-and polarisation-summed squared matrix element by the two-body phase space of the outgoing neutrino and photon, As we will later be considering the sterile neutrino N as an intermediate particle in the scattering process ν α A → XAγ, we do not average over its spin. It is straightforward to find (neglecting the mass of X) spins, pols where we note the factor of two difference between the Dirac and Majorana case. Writing we see that the form of the squared matrix elements allows to integrate over φ and cos θ, The Majorana decay rate via the transition magnetic moment µ XN is a factor of two larger than the corresponding Dirac decay rate.
We now return to the radiative upscattering process ν α A → XAγ. The Feynman diagram for this process is shown to the right of Fig. 5, where we define the momenta of the incoming and outgoing particles and the Mandelstam variables s, s 1 , t 1 , s 3 and t 2 .
The variables s and t 1 = q 2 are equivalent to s and t for the ν α A → N A process. In the lab frame, these are where E γ and θ γ are the outgoing photon energy and angle, respectively, and θ R is the outgoing nuclear recoil angle. The angles are defined to lie between the direction of the incoming neutrino and the outgoing states.
We can again compare the scenarios where the sterile neutrino N is a Dirac or Majorana fermion. For Dirac N , an incoming active Dirac neutrino (or Majorana neutrino with negative helicity) triggers the process ν α A → XAγ. Conversely, an active Dirac antineutrino (or Majorana neutrino with positive helicity) induces the processν α A →XAγ.
In each case we neglect the incoming 'wrong' helicity Majorana neutrino. The amplitude for the ν α A → XAγ process for both Dirac and Majorana ν α can thus be written as where p N = p 2 + p 3 , Γ N is the total width of N , and F µνρσ ≡ − For Majorana N , additional processes are possible. For example, if ν α and X are Dirac, the lepton number violating processes ν α A →XAγ is allowed. The situation is similar if ν α and X are Majorana; now the outgoing state X can have negative or positive helicity.
However, we emphasise that the outgoing state X is not measured. We then need to sum the matrix elements for the processes ν α A → XAγ and ν α A →XAγ if X is Dirac (or simply the matrix element for where the decay vertex of N now contains the sum of chirality projectors (P R + P L ) = 1.
The P R projector isolates the momentum term / p N in the sterile neutrino propagator, while the P L projector selects the mass term m N . The former corresponds to the process ν α A → XAγ and the latter to ν α A →XAγ. We emphasise that lepton number is not measured in the final state. Considering instead an incoming Diracν α (or a Majorana ν α of predominantly positive helicity), the matrix element for is given by Eq. (B28) with the The differential cross section for the ν α A → XAγ process can be found by taking the absolute square of the scattering amplitude, averaging over the possible spins of the incoming nucleus and summing over the spins of the outgoing nucleus. We also sum over the polarisations of the outgoing photon. Neglecting the mass of light neutrinos (and any sterile state N ) in the final state, the differential cross section is given by where dΦ 3 is the three-body phase space for the outgoing nucleus, light neutrino and photon, In the second equality, the integral has firstly been decomposed into a pair of two-body phase spaces (and a trivial integral over the azimuthal orientation of the system) and then written in terms of four Mandelstam variables s 1 , t 1 , s 3 and t 2 . The function is the 4 × 4 symmetric Gram determinant, with ∆ 4 < 0 defining the physical region of the phase space [39].
where H µν is given in Eq. (B13) and the leptonic parts are Inserting Eqs. (B30) and (B32) into the cross section formula of Eq. (B29) now gives the Lorentz-invariant differential cross section where we have integrated over the azimuthal angle φ.
Connection to observables: From Eq. (B35) we wish to compute the differential cross sections in the experimental observables of interest, in particular the nuclear recoil energy E R , the outgoing photon energy E γ and angle θ γ (between the incoming neutrino and outgoing photon). Because Eq. (B35) depends on four Mandelstam variables, we need an other variable in the lab frame, which we choose to be the angle θ R (between the incoming neutrino and outgoing recoiling nucleus). We have already given the Mandelstam variables in terms of these labe frame quantities in Eqs. (B21) to (B25). We see that E R and θ R only appear in s 1 and t 1 and E γ and θ γ in s 3 and t 2 . To determine the single differential cross section in E R , the first step is to therefore integrate Eq. (B35) over s 3 and t 2 , i.e.
The limits of integration are such that the complete physical region of the phase space, ∆ 4 < 0, is integrated over. The limits s ± 3 (s 1 , t 1 , t 2 ) are hence found by solving ∆ 4 = 0 for s 3 , while t ± 2 (s 1 , t 1 ) are found by solving s + 3 = s − 3 for t 2 . Performing this integration for the full differential cross section in Eq. (B35), we obtain for Dirac N while for Majorana N , Taking the ratio of these cross sections, we see that they vary by the factor The first term corresponds to the process ν α A → XAγ which possible for both Dirac and Majorana N . The second term instead corresponds to the process ν α A →XAγ which requires a helicity flip of N and is only possible for Majorana N .
To determine the differential cross section in the relevant lab frame quantities E R and The differential cross sections in the nuclear recoil energy E R and angle θ R can now be computed by integrating over the remaining variable, The allowed region is bounded by the upper limits .

(B43)
In Fig. 6 we depict the kinematically-allowed region as the grey shaded area in the main plot. The kinematically allowed region can be seen to be independent of the sterile neutrino mass m N . In the sub-plots above and to the right, we show the double differential cross section in Eq. (B40) for fixed θ R = 0.5 rad (above) and E R = 10 −5 MeV (right) and three different values of m N . For illustrative purposes, we set the values of the transition magnetic moment and total decay width of N to be µ α νN = 3 × 10 −8 µ B and Γ N = 10 −3 MeV, respectively. A total decay width of this size would require additional invisible decay modes of N . We observe that, even though the double differential cross sections are non-zero over the entire kinematically allowed region, they are dominated by sharply peaked regions. This is a consequence of the total decay width of N being much smaller than the mass of N , justifying the use of the narrow width approximation (NWA).
In the Γ N m N limit, the following replacement can be made in Eq. (B35), which sets the intermediate N to be on-shell, i.e. s 1 = p 2 N = m 2 N . Inserting the expression of s 1 in terms of the lab frame variables in Eq. (B23) into s 1 = m 2 N allows to find the following relationship between the nuclear recoil angle and energy, or equivalently, where c R = cos θ R . In Fig. 6 we plot the curves of allowed values in the (E R , θ R ) plane from the condition s 1 = m 2 N for three different values of m N . From the plots above and to the right, we see that the double differential cross section is sharply peaked at these values of E R and θ R . For each value of m N there is an maximum recoil angle, as well as minimum and maximum recoil energies given by Eq. (B46) with c R = ±1. We now take Eqs. (B37) and (B38), make the substitution in Eq. (B44) and integrate over s 1 by setting s 1 = m 2 N . To obtain the differential cross section in the nuclear recoil energy, we finally multiply by the Jacobian factor ∂t/∂E R = −2m A to obtain In the NWA, the differential rate in the nuclear recoil for the ν α A → XAγ process is therefore the Primakoff upscattering cross section multiplied by the branching ratio for the radiative decay N → Xγ.
To obtain the differential cross sections for ν α A → XAγ in the outgoing photon energy E γ and angle θ γ , we can instead integrate Eq. (B35) over s 1 and t 1 as where again the limits s ± 1 (t 1 , s 3 , t 2 ) are found by solving ∆ 4 = 0 for s 1 and t ± 1 (s 3 , t 2 ) by solving s + 1 = s − 1 for t 1 . In general both integrals are non-trivial for the differential cross section in Eq. (B35) because s 1 appears in the factor of [(s 1 − m 2 N ) + m 2 N Γ 2 N ] in the denominator and t 1 appears in the nuclear form-factor F(t 1 ); it is therefore necessary to perform this integration numerically. Once this is done, it is possible to transform to the lab frame by multiplying Eq. (B49) by the Jacobian ∂(s 3 ,t 2 ) ∂(Eγ ,θγ ) = 4m A E ν E γ sin θ γ . To obtain the single differential cross sections in the variables E γ and θ γ , the remaining  Fig. 2, the contours depict the size of the double differential cross section in E γ and θ γ . The side plots now depict the double differential cross sections for fixed values of E γ (right) and θ γ (above) as dashed and dotted lines. Also shown are the single differential cross sections found by integrating over all allowed values of the other variable (solid lines).
variables are integrated over as where the maximum allowed photon energy is However, the NWA can also be used to simplify the calculation above. The substitution in Eq. (B44) can be made in Eq. (B49) and the s 1 integration performed by setting s 1 = m 2 N . However, the t 1 integral must be still be performed numerically due to the non-trivial dependence of F(t 1 ). Multiplying the double differential cross section by the Jacobian ∂(s 3 ,t 2 ) ∂(Eγ ,θγ ) , we obtain the double differential cross section in E γ and θ γ , In the following, we set F(t 1 ) = 1 and perform the integral over t 1 analytically. In the NWA, the limits of integration t ± 1 can be found by solving ∆ 4 = 0 for t 1 with s 1 = m 2 N . These limits correspond to the minimum and maximum photon energies In Fig. 7, we plot the double differential cross sections in E γ and θ γ for Dirac (left) and Majorana (right) N . We choose the values m N = 1 MeV, Γ N = 10 −11 MeV and µ α νN = 10 −7 µ B . In the sub-plots above and to the right of the contours, we also plot the double differential cross section for fixed values of the photon energy E γ (right) and angle θ γ (above). Furthermore, we plot the single differential cross sections in E γ and θ γ by integrating over the other variable as in Eqs. (B50) and (B51). Examining the single differential cross sections in the photon energy E γ , we see a stark difference between the Dirac and Majorana cases. In the former case, the cross section decreases linearly with the energy, while in the latter the cross section is constant. The minimum and maximum photon energies in Eq. (B53) can clearly be seen. Looking at the single differential cross sections in the photon angle θ γ , the difference between the Dirac and Majorana cases is less prominent; the Majorana cross section peaks at slightly lower angles compared to the Dirac cross section.
Differential rates: With the differential cross section for the Primakoff upscattering in Eq. (B15), we can now calculate the differential rate of nuclear recoil events in a CEνNS experiment as where dφν α dEν is the flux of incoming neutrinos ν α per cm 2 per second and E min ν is the minimum incoming neutrino energy that can produce a sterile neutrino of mass m N and a nuclear recoil energy E R , In Eq. (B55) we have divided by the mass number A of the target isotope multiplied by the proton mass m p to determine the differential rate per unit mass of the target material.
Similarly, the differential rate for nuclear recoil events that are coincident with an outgoing photon in the detector is given (in the NWA) by Here, we simply integrate over the flux dφν α dEν multiplied by the radiative upscattering cross section and the probability P det N = 1 − exp(− L det Γ N βγ ) for the decay to take place inside the detector. In the NWA, E min ν (E R ) is again given by Eq. (B56). As an example, in Fig. 8 (left) we plot the differential Primakoff upscattering rate in the nuclear recoil energy for the NUCLEUS experiment situated at the very-near-site (VNS) of the Chooz reactor site, for m N = 1 MeV and µ e νN = 10 −10 µ B . The flux of electron antineutrinos induces the processν e A →N A. In the plot we compare the rates for three different target materials; 73 Ge (grey), Al 2 O 3 (blue) and CaWO 4 (orange). For Al 2 O 3 and CaWO 4 we average over the mass numbers of the constituent isotopes. We also compare the Primakoff upscattering rates to the SM CEνNS processν e A →ν e A (dotted lines). For µ e νN = 10 −10 µ B , the number of Primakoff upscattering events is comparible the number of CEνNS events. The horizontal black dotted line indicates the predicted nuclear recoil background of 100 keV −1 kg −1 day −1 in the NUCLEUS detector at the VNS. The grey shaded region shows the range of nuclear recoils that can be detected by the experiment (i.e. a nuclear recoil threshold of 10 eV).
In Fig. 8 (right) we plot the differential rate in the nuclear recoil energy for the radiative upscattering processν e A →ν e Aγ, again for the NUCLEUS experiment at the VNS and for m N = 1 MeV and µ e νN = 3 × 10 −8 µ B . We assume that the outgoing antineutrino is of electron flavour,ν e A →ν e Aγ (i.e. there are no additional decay modes of N ), and that the nuclear recoil is accompanied by an outgoing photon. We compare the cases where the intermediate sterile neutrino N is a Dirac (solid lines) or Majorana (dashed lines) fermion, again for three different target materials.
One may first notice that there are peaks in these distributions, which otherwise have the same shape as the Primakoff upscattering distributions. The origin of these peaks is as follows; as we assume that N can only decay radiatively, the branching ratio in the integrand of Eq. (B57) is unity. For µ e νN = 3 × 10 −8 µ B and L det = 5 cm, the total width satisfies Γ N = Γ D(M) N →νeγ βγ/L det and the decay probability can be Taylor expanded to give P det N ≈ L det Γ N /βγ. The peaks occur at values of the nuclear recoil energy that minimise the minimum incoming neutrino energy in Eq. (B56) and hence minimise βγ = γ 2 − 1 ≈ Eν m N 2 − 1. This in turn maximises P det N and the integrand in Eq. (B57). In other words, for these values of the nuclear recoil, the incoming neutrino energy is as close as possible to m N . The produced sterile state is non-relativistic and has a shorter decay length N = βγτ N = βγ/Γ N than average. If more decays occur inside the detector, we would expect to observe a peak in coincidence events.
It is also useful to obtain the differential rate for the radiative upscattering process in where X γ = {E γ , θ γ }. In the NWA, the minimum incoming neutrino energy that can produce sterile neutrino of mass m N and a photon of energy E γ is For a given incoming neutrino energy E ν and sterile neutrino mass m N , the outgoing photon can be emitted at any angle in the range θ γ ∈ [0, π]. In the NWA, the minimum incoming neutrino energy is that which can produce a sterile neutrino of mass m N , In Fig. 9 we show the differential rates for theν e A →ν e Aγ process in the outgoing photon energy E γ (left) and photon angle θ γ (right) for a detector situated at the VNS of the Chooz reactor site. We again present the distributions for three target materials: