Disentangle Neutrino Electromagnetic Properties with Atomic Radiative Pair Emission

We elaborate the possibility of using the atomic radiative emission of neutrino pair (RENP) to probe the neutrino electromagnetic properties, including magnetic and electric dipole moments, charge radius, and anapole. With the typical O(eV) momentum transfer, the atomic RENP is sensitive to not just the tiny neutrino masses but also very light mediators to which the massless photon belongs. The neutrino EM properties introduce extra contribution besides the SM one induced by the heavy W/Z gauge bosons. Since the associated photon spectrum is divided into several sections whose boundaries are determined by the final-state neutrino masses, it is possible to identify the individual neutrino EM form factor elements. Most importantly, scanning the photon spectrum inside the particular section with deviation from the SM prediction once observed allows identification of the neutrino EM form factor type. The RENP provides an ultimate way of disentangling the neutrino EM properties to go beyond the current experimental searches or observations.


Introduction
The existence of neutrinos is indicated by the continuous beta decay spectrum and informally proposed by Pauli in his famous letter to the participants of the Tübingen conference on radioactivity on December 4, 1930 [1].In exactly the same letter, Pauli suggested that "the neutrino at rest is a magnetic dipole with a certain moment" as the possible interaction that neutrino may have.So, the neutrino EM properties are actually the very first studied.While Pauli only mentioned magnetic moment, the general neutrino (ν) coupling to photons (A) contain four independent form factor types at the vanishing momentum transfer limit (q 2 → 0) [2,3], H EM = ν −iσ µν q ν (µ ν + iϵ ν γ 5 ) + γ µ − q µ/ q q 2 q 2 ⟨r 2 ν ⟩ 6 + a ν γ 5 νA µ (q), ( where µ ν , ϵ ν , ⟨r 2 ν ⟩, and a ν are the neutrino magnetic dipole moment (MDM), electric dipole moment (EDM), charge radius, and anapole, respectively.Within the three-neutrino paradigm, all the neutrino EM properties are Hermitian matrices.In the Standard Model (SM) of particle physics, the tree-level coupling between neutrino and photon is zero [4] while the loop-induced electromagnetic couplings are highly suppressed.While the neutrino charge radius and anapole are suppressed by the Fermi constant [5,6], the neutrino magnetic moment and electric dipole are further suppressed by the neutrino mass [5,[7][8][9].
However, all existing searches presented above have intrinsic caveats.With neutrinos not being detected and especially difficult to disentangle their mass eigenstates, constraints are obtained as a combination of parameters involving both the neutrino mixing and form factor matrix elements [27] to introduce various degeneracies [28].A single or just a handful of observations cannot disentangle and identify individual form factor elements.One possible solution [29] is the radiative emission of neutrino pair (RENP) from excited atoms [30][31][32][33][34][35].Since atomic transition energies are typically O(eV), close to the neutrino mass values, RENP is an ideal place to test the light mediators that couples with neutrino and electrons [36], including the massless photon.
In this work we further show that RENP can disentangle all the neutrino EM properties (µ ν , ϵ ν , ⟨r 2  ν ⟩, and a ν ) with details.In Sec. 2, we describe the general electromagnetic vertex properties to establish the formalism of neutrino EM form factors and summarize the existing experimental/observational constraints.The selection rules and details of how neutrino EM properties affect the RENP photon spectrum are presented in Sec. 3. Based on the theoretical formulas, we illustrate a practical strategy of scanning the photon spectrum to identify each form factor elements and their type in Sec. 4. Our discussions and conclusion can be found in Sec. 5.

Electromagnetic Form Factors
The neutrino EM interaction is characterized as the neutrino current (j with the momentum transfer q ≡ p i − p j as difference between the initial (p i ) and final (p j ) neutrino momenta.Since the neutrino current vertex Λ ij µ (q) is a function of the momentum transfer q, it can be expanded as a linear combination of q µ , γ µ and σ µν q ν to match the current Lorentz structure.Each type can be associated with a γ 5 to double the form factors F ij A (q 2 ) with A = 1, . . ., 6, [37], Then the coupling with a massless photon takes the form as L ∋ A µ j (ν) µ .
Note that the quantum electrodynamics (QED) gauge symmetry U (1) em is not broken even in the presence of neutrino EM interactions.Consequently, the U (1) em gauge invariance implies the conservation of the neutrino current, q µ j (ν) µ = 0, which restricts two of the 6 form factors.Of the 6 form factors, the first two In other words, F ij 1 has only off-diagonal elements.For the diagonal case, only F ij 3 survives and the current conservation holds with vanishing ūi/ qu j = (m i − m j )ū i u j and i = j.Note that the current conservation does not require on-shell photon but on-shell fermions.Then Eq. (2.2) reduces to, To make the physics picture more transparent, we have redefined, , respectively.The hermiticity of the interaction Lagrangian further requires that ( The form factors at q 2 → 0 give the values of the magnetic moment , and anapole a ij ν ≡ f ij A (0).In the zero momentum transfer limit, these form factors have vanishing effect with the only exception of the electric charge while the other three are suppressed by either q 2 or q ν .Consequently, the neutrino charge is highly constrained, |q ν | < 3 × 10 −21 e [22] and the next order expansion of f ij Q provides the charge radius, ⟨r For Majorana neutrinos, the current in Eq. ( 2.1) contains an extra term.While Dirac neutrino fields are formed with two sets of operators a and b, ν (D) ∼ au + b † v, the Majorana case has only one a, ν with switched indices i and j.The charge conjugation transformation matrix C in the spinors, u = v T C −1 , is used to swap p i ↔ p j in the above equation to, Since the charge conjugate transformation of the Lorentz structures follows, the extra contribution for Majorana neutrinos does not mix Lorentz structures and introduces at most a minus sign.
The only difference between Dirac/Majorana neutrinos with EM interactions is the form factor matrix pattern, In other words, some form factors are forbidden for Majorana neutrinos and consequently the number of independent parameters is significantly reduced.Observing nonzero diagonal element q ii ν , µ ii ν , or ϵ ii ν is a direct evidence of the Dirac nature of neutrinos.In addition, the existence of electric dipole moment ϵ ij ν indicates possible CP violation.For comparison, anapole violates both parity and charge conjugations, but preserves CP [37].

Existing Terrestrial and Celestial Tests
The constraints on the neutrino EM form factors are obtained from a variety of terrestrial experiments and astrophysical observations [37,38].Nevertheless, these existing measurements can only constrain combinations of the relevant form factor elements.In other words, there is no way to distinguish and identify the individual ones.

Terrestrial Tests with Electron Scattering
The terrestrial experiments typically use neutrino scattering with electron or nuclei via a tchannel photon exchange to probe the neutrino EM form factors. Either reactor, accelerator, or solar neutrinos are possible sources for such experiments.Currently, the strongest constraint comes from the neutrino scattering with electron.Such process can already be induced by the SM weak interactions with differential cross section [39,40], where the upper (lower) sign is for the (anti-)neutrino scattering.The SM neutrino and electron EW interactions have been summarized into the axial and vector current couplings Both the charged and neutral currents mediated by the W ± and Z gauge bosons can contribute.The neutrino mixing matrix elements U * ei and U ej arise from the charged current contribution with the flavor index e for the involved electron in the initial and final states.Also, the neutrino external states are in their physical mass eigenstates.In other words, Eq. (2.8) considers the incoming neutrinos as incoherent states.The more general case with coherence, especially for neutrino oscillation experiments, is elaborated below.
Non-zero neutrino charge radius [41] and anapole [42] modify the effective vector current coupling v ji → v ji as, (2.9) Although the charge randius and anapole are associated with the neutrino vector and axial currents, respectively, their contributions can altogether combine into v ij but not a ij .Note that v ij and a ij are actually associated with the electron bilinears in the effective four-fermion operators, not the neutrino ones.
Since both charge radius and anapole appear in the same place of v ij , an unfortunate degeneracy between them exists.Making it worse, the incoming neutrinos are typically coherent linear combination of mass eigenstates and oscillate over distance while the final-states neutrinos should be treated as the physical mass eigenstates.Being illusive, all the final-state neutrinos ν j contribute to the electron scattering event rate as an inclusive cross section.The actual observable is a summation over the indices i and j in Eq. (2.8) with coefficients given by the neutrino mixing matrix elements [27].In addition, the tiny neutrino masses are negligible and different final-state neutrinos have no actual difference.These lead to an accidental coincidence that the final-state neutrinos can also be taken as flavor states.For short baseline, no oscillation occurs and the couplings are written in the flavor basis, v αβ ≡ ij U αi U * βj v ij .Since a αβ = −δ eα δ eβ + δ αβ /2 and v αβ = −a αβ + 2 sin θ w δ αβ are diagonal, the differential cross section in Eq. (2.8) becomes, ν . (2.10) Note that the diagonal elements ⟨r 2 ν ⟩ αα and (a ν ) αα appear in the first line while the off-diagonal ones in the second line.Although the diagonal elements have interference with the SM contribution, the off-diagonal ones are standalone as square term in the second line.
For convenience, we have defined the effective charge radius as, With degeneracy among the charge radius, anapole, and the neutrino mixing matrix elements U αi U * βj , it is very hard to disentangle and identify the individual neutrino EM form factors.
A similar situation occurs for the magnetic (µ ν ) and electric dipole (ϵ ν ) moments.The cross-section in Eq. (2.8) receives corrections as [43], with X denoting any possible neutrino flavor or mass eigenstates.The effective magnetic moment |µ eff α | 2 is defined as combination of the magnetic (µ ν ) and electric dipole (ϵ ν ) moments, When defining the effective MDM µ eff α , we have used the property that µ † ν = µ ν and ϵ † ν = ϵ ν are hermitian.Since (µ ν ) ij and (ϵ ν ) ij appear in similar ways, the neutrino scattering process cannot distinguish them.Multiple degeneracies exist not just for the charge radius and anapole, but also EDM and MDM.Although there are four form factors each with 6 elements, only few independent combinations can be measured.

Terrestrial Tests with Nuclei Scattering
With low momentum transfer (E ν < O(10) keV), neutrinos interact coherently with atomic nucleus and the cross section is enhanced by the number of nucleons squared.In addition, the differential cross section receives one more enhancement from the photon propagator 1/q 2 for small q 2 .With larger cross-section, the coherent elastic neutrino-nucleons scattering (CEνNS) is a promising probe for neutrino EM properties.The CEνNS cross-section with EM interactions is slightly different from the electron one in Eq. (2.10) and Eq.(2.12) with just neutral current in the SM [46], where Z (N ) is the number of protons (neutrons) and F p (F n ) the proton (neutron) form factor. Similar to Eq. (2.10),only the diagonal charge radius elements have interference with the proton/neutron vector ( The axial current contribution is proportional to the difference between nucleons with spin-up and spin-down which is typically small for the atoms in CEνNS experiments and hence can be neglected [47].Clearly, the CEνNS constraints suffer from the parameter degeneracy discussed before.The COHERENT experiment bounds [48] on the diagonal matrix elements ⟨r 2 ν ⟩ eff αα and µ eff α are typically 1∼2 orders of magnitude larger than those in Eq. (2.14) and Eq.(2.15) [42,46], .17) at 90% C.L. and for the transition matrix elements, (2.18)

Solar Neutrino Measurements with Electron Scattering
In addition to man-made sources, the solar neutrinos can also be used to measure the neutrino EM properties.With nontrivial transition probability the electron neutrino in the solar interior to incoherent mass eigenstates when arriving the solar surface, as arising from the MSW effect [49][50][51] with adiabatic evolution [52], the effective MDM µ ⊙ ν and the effective charge radius ⟨r 2 ν ⟩ ⊙ become dependent on the neutrino energy through the effective mixing matrix U ei at the production places inside the Sun [53], With non-uniform matter density inside the Sun, the effective mixing matrix U ei and hence the transition probability are radius dependent functions.Consequently, the effective mixing matrix should contain an integration over the solar radius weighted by the matter density.However, the matter effect either ceases for low-energy neutrinos (E ν ≲ 1 MeV), with θ 12 = θ 12 , or dominates for high-energy neutrinos (E ν ≳ 5 MeV), with θ 12 = π/2.It is then more convenient to use different values for the effective magnetic moment in individual regions.The Borexino experiment [54] is sensitive to the pp and 7 Be neutrino spectra at E ν ≲ 0.42 and 0.86 MeV, respectively, to give constraints [18,55], at 90% C.L. The Super-Kamiokande (SK) experiment has a 5 MeV energy threshold and measures the 8 B neutrino flux with E ν ∈ [5,20] MeV [17] to obtain, at 90% C.L [17,56].
Since the solar neutrino flux is much lower than the man-made one, larger detector (than the kilogram detector used in the CEνNS experiment) is needed.In addition, with the pp flux dominating at low energy (E ν < 400 keV), a detector with lower energy threshold is better.These two factors renders an opportunity for the dark matter direct detection experiments to probe neutrino EM properties.While the Xenon1T experiment claimed a likely signal of the magnetic moment µ ⊙ ν = 2.2 × 10 −11 µ B at 2.7 σ level [19] in the June of 2020, PandaX-II could not confirm but put a stringent bound < 4.9 × 10 −11 µ B at 90% C.L. [20].The bound is further improved by the recent XENONnT [57] and LUX-Zepelin [21] measurements, respectively.While the strongest constraints for the charge radius comes from PandaX-II and XENONnT [58], at 90% C.L.

Neutrino Electromagnetic Decay
With non-zero neutrino masses, the presence of EM interactions allows for neutrino visible decay into photon.For m i > m j , the reaction ν i → ν j + γ may occur.Since photons are on their mass-shell q 2 = 0, only µ ν and ϵ ν that are proportional to a single q can contribute while the charge radius and anapole with q 2 suppression has vanishing effect.The decay rate [22], is a sum of matrix elements squared.So the constraint on individual element is not relaxed from the combined one.However, the constraints from observing neutrino decay is not that strong.
The reactor experiments put upper bounds on the decay width from the absence of decay photons while the solar neutrino flux reduces in the presence of neutrino decay [22].For the effective magnetic moment appears in the decay process, the bounds in [22] are, at 90% C.L. for reactor and solar neutrinos respectively.
The decay photons can also distort the CMB blackbody spectrum which fits the COBE/FIBRAS data quite nicely with precision at the level of 10 −4 [59] or equivalently [23], for the normal (NO) and inverted (IO) mass orderings, respectively.

Celestial Tests with Stellar Cooling
The strongest constraints for neutrino EM properties are obtained from astrophysical observations.The stellar cooling caused by plasmon decay (γ * → νν) are sensitive to the combination [60], where m γ * ∼ 10 keV is the plasmon mass [60].The constraint from the red giant cooling is |µ ⃝ ⋆ ν | < 2.2 × 10 −12 µ B at 90% C.L. [26] while the white dwarf cooling provides |µ ⃝ ⋆ ν | < 2.9 × 10 −12 µ B at 90% C.L. [24,25].These numbers convert to 4) × 10 −32 cm 2 for the charge radius and anapole.Although the constraint on individual elements cannot be relaxed with sum of squared terms, these numbers are subject to large theoretical uncertainty for the complicated environments inside a compact star [61].

Cosmological Constraints
Cosmological observations can also constrain the electromagnetic moments for Dirac neutrinos.The right-handed neutrinos can be thermally produced in the presence of nonvanishing µ ν or ϵ ν in the early Universe and modify the effective relativistic degrees of freedom, N eff .The bounds range from [62].Since no right-handed neutrinos can be produced without changing the neutrino chirality, this constraint does not apply for the charge radius or anapole.

Caveats
All bounds available have caveats.The cosmological bounds are valid for Dirac neutrinos only while the stellar cooling bounds suffer from uncertainties in the star modeling [61] and can be evaded in the presence of BSM light particles [63].Finally, all existing constraints are sensitive to combinations of EM parameters and even the neutrino mixing matrix elements.In other words, they suffer from parameter degeneracies [28].As we elaborate below, the RENP process of atomic transitions can provide a unique way of probing the individual elements of the neutrino EM form factors. Once experimentally realized, RENP would become a powerful tool for measuring the neutrino properties.

RENP with Neutrino EM Form Factors
As elaborated in Appendix A, the total Hamiltonian characterizing the RENP process has three parts, where H 0 describes the atomic energy levels such that The whole process |e⟩ → |g⟩+γ+ν ν is stimulated using two back-to-back lasers at frequencies ω and ω ′ , with ω < ω ′ .While the final-state photon is emitted in the same direction as the laser with smaller frequency, E γ = ω, the neutrino pair is emitted in the opposite direction with combined energy, E ν + E ν = ω ′ .In addition to E γ = ω, stimulation requires the frequencies to be tuned to match the energy difference between states, ω + ω ′ = E e − E g according to the global energy conservation [64,65].

Selection Rules for the QED and Weak Transitions
The stimulation is only possible if the excited state is metastable with a lifetime longer than 1 ms [32].In comparison, a typical single-photon emission process |e⟩ → |g⟩ + γ has O(1) ns lifetime [32].So the single-photon emission needs to be forbidden at the first order (E1 or M1).The emission of a single photon changes the total angular momentum from the initial (J e ) to the final one (J g ) by at most one single unity, J g − J e = 0, ±1, with the J g = J e = 0 The selection rules for the E1 and M1 transitions between the initial i and final f states.In addition, no transition with J i = J f = 0 can occur.transition forbidden.At the same time, the magnetic quantum number should also change at most by one unity, M Je − M Jg = 0, ±1.The relative parity between the initial (π e ) and final (π g ) states determines if the transition is E1 type with parity flip (π g = −π e ) or M1 type with parity conserved (π g = π e ).This is summarized in Table 1.For the purpose of forbidding the first-order E1 and M1 transitions, one may choose atomic states with either (1) J g and J e = 0 or (2) |J g − J e | > 1.Two possible atomic candidates found in the literature are Yb with J g = J e = 0 and Xe with J g = 0 and J e = 2 [64].
As summarized in Appendix A.2, the photon emission Hamiltonian up to leading terms is, The first term contains the spatial position operator x.Being parity odd, the position operator controls the E1 transitions.The second term connects the photon magnetic field (B = ∇×A (γ) ) with the spin operator ŝ ≡ σ/2.Since ŝ is parity even, the second term controls the M1 transitions.
The E1 and M1 transitions in QED have quite different strength.The E1 type photon emission is regulated by the dipole transition operator d gv ≡ e Ev−Eg ω ⟨g|x|v⟩ which is of the order of the Bohr's radius a 0 = (m e α) −1 times the energy difference E v − E g .For photon emission, the E1 type transition is typically two orders (∼ 1/α where α is the fine structure constant) larger than its M1 counterpart as detailed in Appendix A.2. So, the photon emission |v⟩ → |g⟩ + γ is selected to be of the E1 type to increase statistics.
Since each emitted neutrino has spin 1/2, the neutrino pair emission (|e⟩ → |v⟩ + ν ν) also changes the total angular momentum by at most 1 unity.Consequently, the emission follows the same selection rules as described in Table 1.With partiy violating weak interactions, the electron coupling to the neutrino current contains both vector and an axial-vector components, µ .As derived in Appendix A.3, the Hamiltonian responsible for neutrino pair emission is, The two terms inside the parenthesis are associated with the dipole operator being proportional to the Bohr radius.Since the dipole operator has odd parity, vj (ν) contributes an E1 type transition while aj (ν) σ an E1×M1 type one.Without 1/m e suppression, the third term provides the leading contribution which is M1 type.For the neutrino pair emission, the E1 transition is ∼ 3 orders (E/m e α with the atomic energy E ∼ O(eV)) of magnitude smaller than the M1 counterpart.
With E1 dipole transition for the photon emission and M1 type for the neutrino pair emission being the optimal choices for increasing statistics, the whole RENP process should be E1×M1 type.In practice, the transition types is selected by choosing the target atom states and the laser frequency tuning.This is true for not just the SM RENP process with QED and EW interactions but also for new physics searches such as the neutrino EM properties.

Transition Matrix Elements
Let us start from the simplest electromagnetic interactions on the electron side which is conceptually clearer to set the stage.The electric dipole moment operator controls the E1 transition with a single-photon emission, |v⟩ → |g⟩ + γ, between the virtual (|v⟩) and ground (|g⟩) states.In the Coulomb gauge, the electron couples with the electric field (E) of the emitted photon through the electric dipole operator (d gv ), where k ≡ (ω, ωp γ ) is the photon momentum.
For the SM interactions and Dirac neutrinos, the neutrino pair emission matrix element is [29][30][31]64], where s ve ≡ ⟨v|σ|e⟩/2 is the spin transition element from the excited (|e⟩) to the virtual (|v⟩) states.For Majorana neutrino, its field ν i ∼ ua i + va † i contains only a single type of creation/annihilation operator a i .Both ν and ν of the neutrino current j (ν) µ = ν i γ µ (1 − γ 5 )ν j can contract with any of the two final-state neutrinos.In addition to the usual contractions ⟨ν i (p ν i )|ν ∼ u(p ν i ) and ⟨ν j (p ν j )|ν ∼ v(p ν j ) that contributes to the Dirac matrix element in Eq. (3.5), another two contractions are possible with ⟨ν i (p ν i )|ν ∼ v(p ν i ) and ⟨ν j (p ν j )|ν ∼ ū(p ν j ) that appears for Majorana neutrinos [31].A relative sign in between comes from switching the fermion fields (operators).Putting everything together, the matrix element for Majorana neutrinos is, where we have implemented the property that the axial vector coupling matrix is Hermitian, a ji = a * ij , according to its definition below Eq. (2.8).
In the presence of neutrino EM interactions, the neutrino pair emission receives extra contribution.The vertex in Eq. (2.1) connects with the electron QED interaction (e ψe γ µ ψ e ) via photon mediatior.The resulting effective Lagrangian is where the neutrino current in the momentum space is j (νEM ) µ ≡ u(p ν i )Λ ij µ (q)v(p νj )/q 2 which contributes to the vector current (J V ) µ → j (νEM ) µ as defined in Appendix A.2, in a similar way as photon in Eq. (3.2).With the neutrino pair emission being practically selected to be an M1 transition, only the second term can survive although suppressed by the electron mass.For the purpose of probing the neutrino EM properties, it is much more advantageous to select E1 transition since the first term in Eq. (3.8) is proportional to the dipole operator and typically α −1 ∼ 100 faster than the M1 transition.Moreover, selecting the E1 type transition can also reduce the neutrino pair emission background from the SM weak W/Z interactions by 2 orders.So an ideal RENP experiment for probing neutrino EM form factors should use the E1×E1 transitions.However, we still adopt the E1×M1 type configuration to make a conservative illustration.
The corresponding matrix element for the M1 transition is, If E1×E1 can be experimentally implemented, the signal strength can be further enhanced by 4 orders while the background is suppressed by 4 orders.

Kinematics and Phase Space
To enhance the decay rate, a laser stimulation is needed.As a consequence, the photon momentum magnitude |p γ | = ω < ω ′ and direction (p γ ) are fixed by the incoming laser.While the neutrino pair carries the energy dictated by the atomic level difference E eg ≡ E e − E g and the photon energy ω, the momentum conservation requires the neutrino pair momentum to be opposite of the emitted photon, The equation above imposes a relation between the trigger laser frequency and the neutrino pair invariant mass, m 2 νν ≡ (p Because of the minus sign in front of m 2 ν ν , a maximum frequency threshold occurs at the minimum neutrino pair invariant mass m 2 ν ν | min = (m i + m j ) 2 [31,64,65], with m i (m j ) being the (anti-)neutrino mass.
If the neutrino pair is replaced by a photon, the two back-to-back stimulating lasers have the same frequency, ω = ω ′ = E eg /2 with vanishing m 2 ν ν .Since the neutrino pair emission only occurs for ω < ω ′ , the two-photon emission background |e⟩ → |g⟩ + 2γ can be suppressed by choosing the stimulating lasers with two different frequencies.It is practically possible to use the E1×E1 configuration to enhance the signal rate from the neutrino EM properties and at the same time suppressing the two-photon background.Although being higher-order in α, other QED backgrounds exists.For example, replacing the neutrino pair by a photon pair or any final state with more photons (|e⟩ → |g⟩ + nγ, n ≥ 3) produces a higher number of events than the RENP process [66].Extra measures are needed to reduce the higher order QED backgrounds to a reasonable level.The use of crystal waveguides is a possible solution to achieve an n-photon background free experiment [67,68].In our work we consider that all QED backgrounds can be removed.
The three-body final-state phase space reduces to the two-body one for the two final-state neutrinos [32] and the differential decay rate is [31,32], where the total matrix element M e→g is obtained from the second-order perturbation theory.
The full element M e→g is a combination of M D in Eq. (3.4) and the neutrino pair emission matrix element M νν [32], with energy difference E vg ≡ E v − E g and the atom number density n a .The neutrino emission matrix contains the two terms, M νν ≡ M W + M νEM , from the EW and EM contributions in Eq. (3.5) and Eq.(3.9), respectively.

Magnetic Moment and Electric Dipole Contributions
When only neutrino magnetic and electric moments are present, Eq. ( 2.3) simplifies and M νEM becomes, The integration d 3 p νj in Eq. (3.13) removes the three-dimensional δ function for the momentum conservation while the remaining δ function for the energy conservation adapts the angular part of d 3 p ν i to correlate the opening angle θ 0 between the neutrino and photon momentums, where ∆m 2 ij ≡ m 2 i − m 2 j and q 2 = E 2 eg − 2E eg ω are obtained from Eq. (3.10).The geometric limits | cos θ 0 | ≤ 1 imply an allowed neutrino energy range, where, Putting everything together, the differential decay rate in Eq. (3.13) becomes a function of the neutrino energy E ν i , The spin-average square of M νν has three parts.First, the weak-only component [31,36,64,65] is, where the quantity δ M = 0 for Dirac neutrinos and is δ M = 1 for Majorana neutrinos.The extra term ∝ 3δ M Re[a 2 ij ]m i m j that appears only for the Majorana neutrino case effectively changes the distribution of emitted photons and can be used to distinguish from the Dirac case and even the Majorana phases [31,32].The averaging is performed over the initial-and final-state spin states.For convenience, the spin operators have been eliminated using the spherical symmetry property of atomic states [31], where J e and J v are the total spins of the excited (|e⟩) and virual (|v⟩) states, respectively.The coefficient C ev is a numerical factor that is different for each atom.For both Xe and Yb, (2J v + 1)C ev = 2 [64].
The interference term is a combination of the weak and EM components, 2Re[M W M * µ,ϵ ], for Dirac neutrinos.After applying the spin-operator identity in Eq. (3.21), the interference term becomes, where E is defined in Eq. (3.18).The only difference between the neutrino magnetic moment and electric dipole is a minus sign in front of m j from the presence of γ 5 in the (ϵ ν ) ij σ µν γ 5 vertex.
While M W M * µ,ϵ is proportional to E − E ν i , the range of E ν i in Eq. (3.17) is anti-symmetric around E. Consequently, the interference contribution to the total rate Γ is zero.Similar cancellation occurs for Majorana neutrinos.Although such interference term can appear in the differential rate dΓ/dE ν i , the neutrino energy cannot be practically reconstructed due to its illusiveness.

The new physics contribution to the total rate Γ comes from the EM-only term, |M
Using the relation in Eq. (3.21) and further noticing that there is no interference between the magnetic (µ ν ) and electric (ϵ ν ) moment contributions due to the presence of a γ 5 matrix in the electric dipole interaction, the matrix element squared becomes, Similar to the Weak/EM interference term, the difference between the electric and magnetic moment contributions is a minus sign in front of m j .
Putting this result back into the differential decay rate in Eq. (3.19), together with the definition of θ 0 in Eq. (3.16) and the integration range in Eq. (3.17), the eletromagnetic total decay rate becomes [29], where the Heaviside Θ-functon defines the ω max ij thresholds in Eq. (3.12).The rate can be charactrized by a benchmark value Γ 0 [64,65], with η describing the fraction of the target volume V that is macroscopically coherent.The spectral function is then defined as the relative size, I µ,ϵ ≡ Γ µ,ϵ /Γ 0 .
Comparing the result in Eq. (3.26) with the SM-only case [31,36,64,65], the relative size 2 scales with µ 2 ν for the magnetic moment and similarly for its electric dipole counterpart.While the SM contribution is mediated by the heavy W/Z bosons, the neutrino magnetic and electric moment contributions  are inversely proportional to 1/q 4 .With typical q ∼ O(eV) for atomic processes, the SM contribution is suppressed by a factor of (q 2 /m 2 W ) 2 ∼ 10 −44 .However the presence of the electron magnetic moment suppresses the EM contribution by ω/µ B ∼ 10 −6 eV.All in all, a |(µ ν ) ij | ∼ 10 −11 µ B provides similar contribution to that of the SM.
The total spectral function, I ≡ I W + I µ,ϵ , is presented in the upper left (right) panels of Fig. 1 with the representative value of In both cases, the curves significantly differ from the SM contribution (black).The kinks in the curves represents the location of all six frequency thresholds in Eq. (3.12).In all curves, the magnetic moment contribution is slightly larger than the electric dipole one due to the difference in the sign in m j .The difference in the shape of I over the frequency range can be used to experimentally distinguish both contributions.

Charge Radius and Anapole Contributions
The charge radius and anapole always appear together with similar Lorentz structure, The total matrix element is then M ν ν = M W + M r,a .Similarly to the previous section, the spin averaged matrix element squared contains three terms, the SM-only term, the interference term between SM and EM contributions, as well as the pure charge radius and anapole terms.

The interference between M (D)
W in Eq. (3.5) and M r,a is, The sum over electron spins is performed using the spin relation in Eq. (3.21) to obtain, 2Re with E defined in Eq. (3.18).As before, this interference term Re M (D) W M * r,a ∝ E ν i − E is an odd function around E and hence does not contribute to the total rate after integration over E ν i .
The only contribution from the neutrino charge radius and anapole is the square of Eq. (3.29), |M r,a | 2 =4µ 2 B 1 2(2J e + 1) spins (0, s ev × q) µ (0, s * ev × q) ν Tr γ µ − q µ / q q 2 (3.32) The presence of γ 5 in the anapole term ensures that there is no interference between ⟨r 2 ν ⟩ ij and (a ν ) ij .After using Eq.(3.21) to eliminate the spin operators, the matrix element squared becomes, (3.33) The only difference between the charge radius and anapole terms is the minus sign in front of m j , which comes from the γ 5 matrix in the anapole interaction.Using the decay rate definition in Eq. (3.19) and the kinematic limits in Eq. (3.17), the total decay rate becomes, where Γ 0 is defined in Eq. (3.27) and the step function Θ(ω −ω max ij ) defines the frequency thresholds of Eq. (3.12).The relative size between the SM contribution in Eq. (3.28) and the anapole term is Γ With the expectation that the sensitivity is projected for the signal to have roughly the same event rate as the SM background, the neutrino charge radius and anapole can be probed at the level of 10 −26 cm 2 .Notice that an experiment with higher trigger frequency ω provides larger sensitivity on the charge radius and anapole due to the ω 2 factor.
The total spectral function I ≡ I W + I r,a is presented in the bottom left (right) panels of Fig. 1 for the representative value of |⟨r 2 ν ⟩ ij |/6 and |(a ν ) ij | of 4 × 10 −26 cm 2 .In both cases, the curves significantly differ from the SM contribution specially for lower frequencies.In all curves, the charge radius contribution is slightly larger than the anapole one due to the difference in the sign in m j , which modified (m i ± m j ) 2 /q 2 in Eq. (3.34).Since q 2 = E 2 eg − 2E eg ω decreases with increasing frequency, the relative difference between the charge radius and anapole curves becomes larger for higher frequency.In the next section, we show that the difference in the shape of I r and I a can be used to experimentally distinguish them.

RENP Sensitivity to Neutrino EM Properties
The advantage of RENP is its capability of separating the individual form factor elements f ij M,E,Q,A of Eq. (2.3) as first noticed in [29].From the hermiticity condition, each form factor contains 6 elements that enter the spectral function as magnitudes squared |f ij M,E,Q,A | 2 while their phases do not appear.The kinematics in Sec.3.4 allows the photon energy up to certain threshold frequency ω ij max defined in Eq. (3.12) for each neutrino pair ν i νj .With the trigger laser having a precision of δω ∼ 10 −5 eV [31] which is much better than E eg /2−ω max ij ∼ 10 −5 −10 −3 eV, it is possible to scan the various thresholds ω max ij with different combinations of neutrino mass eigenvalues and the regions in between.Scanning over 6 frequencies is sufficient to separate individual matrix elements as proposed in [29].For illustration, we take ω i = (1.069,1.07, 1.0708, 1.0712, 1.0716, 1.07164) eV shown as the black points in Fig. 2.
Once nonzero deviation from the SM prediction is observed, the experiment is also capable of identifying which neutrino property (µ ν , ϵ ν , ⟨r 2 ν ⟩, or a ν ) is contributing.If one element with specific indices i and j is experimentally observed to be nonzero, the spectrum function below the corresponding ω max ij is affected and can be scanned to identify which type that element is.At least another 3 frequency points are necessary to reconstruct the photon spectrum and then identify the relevant neutrino property.Since the spectrum function has larger value for smaller ω, we choose the three frequencies ω i = (1.0672,1.0678, 1.0684) eV shown as red points in Fig. 2 to increase event rate.Another reason of choosing lower frequencies bellow all ω max ij is that all elements can contribute to this region.However, if the nonzero element has larger ω max ij , it is also advantageous to add more frequency points.Both the low and high frequency  regions are sensitive to distinguishing ⟨r 2 ν ⟩ 12 and (a ν ) 12 from (µ ν ) 12 and (ϵ ν ) 12 with apparent difference.The difference between (µ ν ) 12 and (ϵ ν ) 12 is also larger at the high frequency region.
Following the setup in [29], we estimate the number of RENP photons at each frequency ω as N i = T Γ 0 I(ω i ), where T is the exposure time that we take to be equal among all the nine ω i values.While V is the target volume, I(ω) is the total spectral function containing the SM and EM contributions.For a typical Yb experiment with T = 10 days, we expect N (ω i=1,...,9 ) = (146, 136, 120, 107, 87, 82, 79, 39, 2.5) SM background events, respectively, and in total 759 events.The sensitivity is evaluated using the Poisson χ 2 function, where N true i (N test i ) are the event numbers with SM only (SM and EM) contributions at each frequency ω i .Fig. 3  ν ⟩ ij (bottom left), and (a ν ) ij (bottom right) with the normal ordering (NO) and m 1 = 0.01 eV.The dotted lines show the current bounds on the effective parameters from reactor (dotted blue for µ eff e or ⟨r 2 ν ⟩ eff ) [16] and solar (dotted purple for µ ⊙ e or ⟨r 2 ν ⟩ ⊙ ) [18, 55] neutrinos, as well as plasmon decay in red giants (dotted red for RG) [26] and white dwarfs (dotted yellow for WD) [24,25].

Discussion & Conclusion
The EM properties are not just the first neutrino interaction that was originally proposed by Pauli in exactly the same letter where he pointed out the existence of neutrino, but can also have sizable values to appear in various BSM models.Hence, observing the neutrino EM properties (magnetic and electric moments, charge radius, and anapole) can serve as a window of probing possible new physics.Unique probe and identification of each form factor elements are necessary which is a serious challenge for the existing experimental measurements and observations.We show that the RENP process is an ideal probe.The typically eV scale momentum transfer of atomic transitions is perfect for probing light mediators including the massless photon.Most importantly, the RENP process do not suffer from degeneracies in EM form factors and hence has the possibility of individual identification as a unique feature.Both kinetic energy (T e ), but contains potential terms related to J V and J A , as will become evident below.Using the above approximation Eq. (A.6) to substitute ϕ back into the left-hand side of Eq. (A.5a) and keeping the leading order in m −1 e , we obtain the non-relativistic electron Schrödinger equation, In other words, the electron Hamiltonian with both QED and EW interactions is, The first term is a kinetic energy in the presence of vector/pseudo-vector fields while the last two terms are the potential energy.
For later convenience, we reformulate the Hamiltonian as, for describing the atomic energy levels, the atomic transitions, and the radiative emission of neutrino pair as we elaborate below.

A.2. Atomic Energy Levels and Transition with QED
The Hamiltonian in Eq. (A.9) can in principle describe any vector and axial-vector interactions.
For electromagnetic coupling with photon, J 0 V = −eA 0 = −eV (x) and J V = −eA.To be exact, J 0 V is the electric potential due to the electron charge and ∇ × J V = −e∇ × A is actually the magnetic field coupled with electron spin.
The photon field has two contributions, A µ ≡ A µ forms the atomic Hamiltonian, H 0 , that describes the atomic energy levels, H 0 |a⟩ = E a |a⟩ with a = g, e, and v for the ground, excited, and intermediate virtual states.To be more concrete, H 0 is a combination of the electron kinetic energy and the classical electromagnetic potential obtained by replacements J V → −eA (0) and J A → 0 from Eq. (A.9),The first is proportional to the space coordinate operator x and hence contributes as dipole operator, d if ∝ ⟨i|x|f ⟩.On the other hand, the second term contains the electron spin operator (ŝ = σ/2) and the magnetic field (B = ∇ × A (γ) ).Both terms are contributed by the threedimensional vector potential A (γ) .
The two terms in Eq. (A.13) have definite spin and parity such that the induced transitions are subject to handy selection rules for the single-photon emission, |i⟩ → |f ⟩ + γ.The total angular momentum from the initial (J i ) to the final one (J f ) changes by at most one single unity, J f − J i = 0, ±1, with the J i = J f = 0 transition also forbidden.At the same time, the magnetic quantum number should also change at most by one unity, M J i − M J f = 0, ±1.The relative parity between the initial (π e ) and final (π g ) states determines if the transition is E1 type with parity flip (π f = −π i ) or M1 type with parity conserved (π f = π i ).Since the dipole operator reverses sign under parity transformation, x → −x, it induces E1 type transition while the magnetic moment operator in the second term of Eq. (A.13) is M1 type.
The two terms in Eq. (A.13) have different typical sizes.The dipole term is in unit of the Bohr radius a 0 = (m e α) −1 while the magnetic moment operator is regulated by the Bohr magneton (µ B = e/2m e ) which is suppressed by the electron mass.Although the dipole term also contains H 0 that gives the atomic energy at O(eV), the ∇ operator in the magnetic moment operator gives the photon frequency which is also at eV scale to balance.Consequently, the E1 transition is typically α −1 ∼ 100 times larger than its M1 counterpart.

Fig. 1 :
Fig. 1: The spectral function I(ω) for the Yb atom RENP processes with weak (SM) and neutrino EM interactions as a function of the trigger laser frequency ω for the normal ordering (NO) of neutrino masses and the lightest neutrino mass m 1 = 0.01 eV.