I$_{2}$ molecule for neutrino mass spectroscopy: ab initio calculation of spectral rate

It has recently been argued that atoms and molecules may become good targets of determining neutrino parameters still undetermined, if atomic/molecular process is enhanced by a new kind of coherence. We compute photon energy spectrum rate arising from coherent radiative neutrino pair emission processes of metastable excited states of I$_2$ and its iso-valent molecules, $|Av \rangle \rightarrow |Xv' \rangle + \gamma + \nu_i\nu_j$ and $|A'v \rangle \rightarrow |Xv' \rangle + \gamma + \nu_i\nu_j$ with $\gamma$ an IR photon and $\nu_{i(j)}$ $i(j)-$th neutrino mass eigenstates, and show how fundamental neutrino parameters may be determined. Energies of electronically excited states of I$_2$, including the effect of spin-orbit couplings were calculated by the multiconfigurational second order perturbation (CASPT2) method. Summation over many vibrational levels of intermediate states is fully incorporated. Unlike atomic candidate of a much larger energy difference such as Xe, I$_2$ transitions from a vibrational level $A(v=0)$ to $X(v' = 24)$ give an opportunity of determination of the mass type (Majorana vs Dirac distinction) and determination of Majorana CPV (charge-conjugation parity violating) phases, although the rate is much smaller.

It has recently been argued that atoms and molecules may become good targets of determining neutrino parameters still undetermined, if atomic/molecular process is enhanced by a new kind of coherence. We compute photon energy spectrum rate arising from coherent radiative neutrino pair emission processes of metastable excited states of I 2 and its iso-valent molecules, |Av → |Xv ′ + γ + ν i ν j and |A ′ v → |Xv ′ + γ + ν i ν j with γ an IR photon and ν i(j) i(j)−th neutrino mass eigenstates, and show how fundamental neutrino parameters may be determined. Energies of electronically excited states of I 2 , including the effect of spin-orbit couplings were calculated by the multiconfigurational second order perturbation (CASPT2) method. Summation over many vibrational levels of intermediate states is fully incorporated. Unlike atomic candidate of a much larger energy difference such as Xe, I 2 transitions from a vibrational level A(v = 0) to X(v ′ = 24) give an opportunity of determination of the mass type (Majorana vs Dirac distinction) and determination of Majorana CPV (charge-conjugation parity violating) phases, although the rate is much smaller.

Introduction
Neutrinos are most common particles next to 3 K photons in the present universe, yet their properties eluded comprehensive experimental determination since these neutral particles have only weak interaction. The conventional target for exploration of neutrino properties has been nuclei; stable nuclei in neutrino oscillation experiments and unstable nuclei in other neutrino experiments. Measured quantities derived so far by neutrino oscillation experiments are mass squared differences (∆m 2 ij ≡ m 2 i − m 2 j ) and mixing angles summarized [1] by ∆m 2 21 ∼ 7.5 × 10 −5 eV 2 , |∆m 2 31 | ∼ 2.5 × 10 −3 eV 2 , sin 2 θ 12 ∼ 0.31 , sin 2 θ 23 ∼ 0.42 , sin 2 θ 13 ∼ 0.024 .
It has been assumed in this analysis that there exist only three kinds of neutrino, which we also follow throughout this work. The other hint on the absolute mass scale is derived from cosmological arguments, giving i m i < O(0.5) eV [2]. The important quantities undetermined in oscillation experiments are (1) absolute neutrino masses (impossible to determine in oscillation experiments), (2) CP symmetry (C=charge conjugation, P = parity) violating phase (CPV phase in short), and (3) whether the massive neutrino belongs to the Majorana type [3] or not. There are two different kinds of ongoing experiments using unstable nuclei for measurements of still undetermined neutrino parameters; (1) the end point spectrum of beta decay of nuclei such as tritium for the measurement of an averaged absolute neutrino mass value, and (2) the search of neutrino-less double beta decay for verification of lepton number violation related to a finite Majorana type of masses. So far negative limits have been set in these experiments. The reason for the use of unstable nuclei is that the weak decay rate such as the nuclear beta decay increases with a high positive power of the released energy, usually 5th power, and the available nuclear energy of a few MeV gives a detectable rate. But with the expected small neutrino mass of a fraction of eV the energy mismatch becomes a serious problem, and the determination of the absolute neutrino mass value and exploration of undetermined important neutrino properties are more and more difficult.
With the advent of remarkable technological innovations manipulation of atoms and molecules may contribute greatly to fundamental physics. Neutrino physics may also be one of these areas. Atoms and molecules are target candidates of precision neutrino mass spectroscopy, as recently emphasized in Ref. [4], due to closeness of the energy released in their transition to expected neutrino masses. The relevant process of our interest is cooperative (and coherent, called macro-coherent subsequently) atomic de-excitation; |e → |g + γ + ν i ν j where ν i(j) , i(j) = 1, 2, 3 is one of neutrino mass eigenstates. γ in the present work refers to a photon in the visible to the infrared region. The initial state |e must be metastable, say its lifetime O(1) msec. The process, as shown in Fig. 1, exists as a combined effect of second order in Quantum Electro-Dynamics (QED) and weak interaction of the standard electroweak theory [5].
To obtain a measurable rate of the process, it is crucial to develop the macro-coherence [4,6], a new kind of coherence. The macro-coherent emission of radiative neutrino pairs is stimulated by two trigger irradiation of frequencies ω, ω ′ constrained by ω + ω ′ = ǫ eg / , ω < ω ′ , with ǫ eg = ǫ e − ǫ g the energy difference of initial and final states. The measured photon energy in the de-excitation is given by the smaller frequency ω. The macro-coherence assures that the three-body process, |e → |g + γ + ν i ν j , conserves both the energy and the momentum. Assuming that atoms in the states |e and |g can be taken infinitely heavy and the atomic recoil may be ignored, there exist threshold photon energies [7] at Each time the emitted photon energy decreases below a fixed threshold energy of ω ij , a new continuous spectrum is opened, hence there are six energy thresholds ω ij , i, j = 1, 2, 3 ((i, j) threshold for brevity) separated by finite photon energies. Determination of the threshold location given by Eq. (3), hence decomposition into six mass thresholds, is made possible by precision of irradiated laser frequencies at ω ≈ ω ij and not by resolution of detected photon energy.

2/24
The macro-coherently amplified radiative emission of neutrino pair has been called RENP (radiative emission of neutrino pair) [4], which is the core idea of our neutrino mass spectroscopy that may determine all unknown neutrino parameters. It has been argued that this method using atoms and molecules is ultimately capable of determining the nature of neutrino masses, Dirac vs Majorana distinction and measuring the new Majorana source of CPV phases [7,8]. It is crucial for explanation of the matter-antimatter asymmetry of the universe to verify the Majorana nature of neutrinos and determine CPV phases related to the Majorana case [9,10].
There is an important constraint on possible target atoms/molecules to obtain reasonable rates for realistic experiments. The atomic operator involved in RENP has the character of M1×E1 where the M1 operator (actually electron spin S in subsequent RENP formulas) governs the weak interaction Hamiltonian of neutrino pair emission and the E1 operator denotes the usual dipole interaction of QED. The total angular momentum change via virtual intermediate state requires that the LS coupling scheme should be broken [4], hence candidates should be sought in heavy atoms/molecules. This way we found Xe atom as a good candidate.
Xe atom is excellent, having a great discovery potential of RENP process, and further expected to determine the absolute neutrino mass scale and distinguish the normal mass vs the inverted mass hierarchy (denoted by NH and IH respectively) [4]. On the other hand, distinction of the mass type (Dirac vs Majorana) requires a smaller released energy, of order a fraction of 1 eV, as shown in Ref. [8]. In the present work we shall study I 2 and its isovalent molecules for this purpose.
Molecules have a number of merits for RENP: (1) homo-nuclear diatomic molecules such as I 2 have vibrational states among which the usual E1 transition (its Hamiltonian ∝ d · E) is forbidden, (2) the richness of vibrational and rotational levels makes it ideal to perform a systematic search of neutrino mass thresholds. Moreover, I 2 molecule has a number of metastable states with energies 1 eV above the electronically ground vibrational excited levels.
Electronic excited states of I 2 molecule have been intensively investigated by quantum chemical calculations [11][12][13]. The potential energy curves (PECs) and spectroscopic constants as well as transition properties of this system have been well examined, while the spin operator element which is relevant for the present neutrino mass spectroscopy has not been focused so much. We shall present in this work RENP spectral rate using molecular wave functions based on the first principle calculation.
We show in the present work that the proposed I 2 de-excitation scheme gives a chance of measurements, the Majorana vs Dirac distinction and determination of CPV phases, which are, in the quasi-degenerate case of neutrino masses, easier than NH vs IH distinction at low photon energies where rates are largest. A good sensitivity of the spectral shape to determination of the smallest neutrino mass of order meV is also shown for this target molecule.
Throughout this work we assume that the macro-coherent mechanism as proposed in Ref. [4]. In most of the rest except Sec. 3 where the atomic unit is used following the standard practice of quantum chemistry, we use natural units of = 1, c = 1.

RENP amplitude and rate formula
We shall recapitulate from Ref. [4] main features of the spectral formulas.
The amplitude corresponding to the Feynman-like diagram of Fig. 1 is given by the electroweak theory and reads as where U ei (their relation to necessary mixing angles and CPV phases given in Eq. (18)) is the matrix element of neutrino mixing (expressions in terms of measurable quantities given later), and ν i ( p, h) the neutrino plane wave function of momentum p and helicity h of mass m i , E the electric field of irradiated trigger laser of frequency ω < ǫ eg /2, S the electron spin operator, and d the electric dipole operator. Quantum numbers of states should include all of electronic, vibrational, and rotational ones. The sum over all vibrational modes of intermediate states, |p and |q , is particularly important for molecules. For simplicity we ignore Hönl-London factor [14] of order unity assuming that rotational degrees of freedom are frozen. Two terms in bracket of the right hand side of Eq. (4) correspond to two different vertices, weak M1 type of neutrino pair emission (its Hamiltonian ∼ G F S · ν † j σν i ) and QED E1 transition vertex in the second order perturbation theory as depicted in Fig. 1: the first term of molecular state change of |e → |p → |g shall be designated by E1×M1 and called (a) in the figure, and the second term of molecular state change of |e → |q → |g shall be designated by M1×E1 and called (b). Both of states |p , |q are summed over, since they are virtual intermediate states. Fig. 1 Feynman-like diagrams for RENP from Λ-type atom/molecule, |e → |g + γ + ν i ν j , with γ a photon and ν i(j) a neutrino mass eigenstate. Virtual intermediate states |p , |q should be summed over. Two-photon transition |e → |g + γ + γ may also occur via weak M1× E1 or E1×M1 couplings to |p or |q .
For isotropic medium without directional alignment by a magnetic field, this amplitude squared gives the basic rate formula (omitting the second contribution C b (ω) from |q in 4/24 Eq. (4)), In the overall rate factor Γ dm for diatomic molecules the directionality of trigger correlated with the molecular axis is taken into account by an extra 1/3 reduction factor.ǭ is a reference energy, to make both C a (ω) and I(ω) dimensionless. In the rest of this work we takeǭ = ǫ eg . The function θ(x) = 0 for x < 0, = 1 for x > 0 is the step function, giving rise to six mass thresholds in Eq. (9). δ M = 1 for the case of Majorana neutrinos and δ M = 0 for the Dirac neutrino. This term ∝ δ M exhibits the quantum mechanical interference intrinsic to identical fermions of Majorana particles [7], which arises from the anti-symmetric wave functions of two identical fermions. For the Dirac neutrino case, emitted particles (except the photon) are a neutrino and an anti-neutrino, two distinguishable particles, and hence the interference terms are absent.
For stored field strength | E| 2 we shall take its maximal value ǫ eg n with n the number density of excited state targets. More generally, this rate should be multiplied by the dynamical factor η ω (t) whose calculation requires solution of the master equation given in Ref. [4]. It is important to keep in mind that the medium polarization R i , i = 1, 2 between |g and |p states under trigger irradiation is contained in the dynamical factor given by an integrated quantity over the entire target at 0 ≤ x ≤ L along the trigger irradiation. Thus, a large macroscopic polarization is required for a large RENP rate. The overall maximal rate in the unit of 1/time is The spectral information is given by I(ω) which is calculated using neutrino parameters experimentally determined, Eq. (2). Calculation of the molecular factor C(ω) is the main subject of the next section.

Molecular factor
In this section, we investigate the possibility of the RENP experiment using I 2 molecule, with the electronically ground state X 0 + g as the |g state and the (metastable) electronically excited states A 1 u and A ′ 2 u as the |e state. An advantage of the A ′ 2 u state over the A 1 u state is its longer natural lifetime. In the simplified model for the RENP process |A ′ v → |Xv ′ + γ + ν i ν j in Ref. [4], only the I 2 A 1 u state is considered as the intermediate |p state. However, other electronic states may have similar or larger contribution to the molecular factor. We first calculate the molecular factor including several I 2 excited electronic states in the intermediate state summation, without considering vibrational states. There, we examine (1) which intermediate electronic state contributes the most significantly to the molecular factor, (2) relative weight of E1×M1 and M1×E1 amplitude to the RENP process. Second, the effect of nuclear motion (vibrational states) on the molecular factor is investigated considering a single intermediate electronic state, in which we will examine how the molecular factor depends on ( , a pair of the vibrational states on the |e and |g states, which may be useful in designing experiments in future. Finally, the magnitudes of the molecular factor are compared for isovalent molecules, I 2 , Br 2 and Cl 2 , to study a dependence of the RENP rate on the spin-orbit couplings.
In order to simulate realistic experimental situation, the I 2 energies including the spinorbit effects are accurately calculated within the framework of the Born-Oppenheimer approximation.

Detail of the calculation
The electronic excited states of I 2 molecule were calculated by the multiconfigurational second-order perturbation method (CASPT2) [15,16], based on the reference wavefunctions obtained by the state-averaged complete active space self-consistent field (CASSCF) method [17,18] with the atomic natural orbital-relativistic correlation consistent (ANO-RCC) all-electron basis set [19]. The CASSCF method, using a linear combination of configuration state functions to describe electronic wavefunction, is adequate for calculation of electronically excited states with relatively small computational cost. More accurate energies, including dynamic electronic correlation, can be obtained by the CASPT2 method using the second order perturbation to the zero-th order CASSCF wavefunctions. The scalar relativistic effects were introduced at the CASSCF level by keeping the spin-independent (scalar) terms of the two-component reduced Hamiltonian using the 4th order Douglas-Kroll-Hess method [20][21][22]. At this point, the spin-orbit effects are not included in the CASSCF wavefunctions and the CASPT2 energies. The spin-orbit effects were considered by the state-interaction method, where the spin-orbit matrix elements were evaluated between the CASSCF states using the Breit-Pauli operator [23], while the unperturbed diagonal 6/24 energies were replaced by the CASPT2 energies. Final energies are obtained by diagonalization of this spin-orbit Hamiltonian. Our procedure is very similar to that employed in the calculation of the I 2 potential energy curves by Malmquist et al. [24]. The D 2h symmetry was employed in our calculations. In the CASSCF calculation, total of 28 orbitals, 1-7a g (i.e. 1a g , 2a g , · · · , 7a g ), 1-3b 3u ,1-3b 2u ,1b 1g ,1-7b 1u ,1-3b 2g ,1-3b 3g and 1a u , were kept doubly occupied. In addition, the 1a g (1σ g ) and 1b 1u (1σ u ) orbitals were frozen during the calculation. The remaining 50 electrons were distributed over 26 active orbitals: 8-13a g ,4-6b 3u ,4-6b 2u ,2b 1g ,8-13b 1u ,4-6b 2g ,4-6b 3g and 2a u . The state-averaging was performed over 18 electronic states: the three lowest A g , one B 3u , one B 2u , one B 1g , one B 2g , one B 3g and one A u states with singlet spin multiplicity, and one B 3u , one B 2u , one B 1g , three B 1u , one B 2g , one B 3g and one A u states with triplet spin multiplicity. Note that these electronic states correlate with two iodine atoms in the 2 P state at the dissociation limit. By diagonalization of the spin-orbit Hamiltonian, total of 36 spin-orbit eigenstates were obtained. All these calculations were performed using the MOLPRO suite of programs [25]. Although the electric transition dipole moments between the spin-orbit eigenstates were obtained as part of the spin-orbit calculation in the MOLPRO programs, the matrix elements of the electronic spin operator S were not provided. Thus, we explicitly evaluated the spin matrix elements between the spin-orbit eigenstates using the output of the eigenvectors. Excitation energy (eV) I( 2 P 1/2 ) + I( 2 P 1/2 ) I( 2 P 1/2 ) + I( 2 P 3/2 ) I( 2 P 3/2 ) + I( 2 P 3/2 ) The potential energy curves of the calculated I 2 electronic states are shown in Fig. 2. In the figure, the spin-orbit eigenstates are correlated with pairs of the atomic states I( 2 P 3/2 )+I( 2 P 3/2 ), I( 2 P 1/2 )+I( 2 P 3/2 ) or I( 2 P 1/2 )+I( 2 P 1/2 ), as the inter-nuclear distance becomes large. The dissociation and excitation energies of the electronic states relevant to this work, as well as the vibrational energies on these states are summarized in Table 1. For the purpose of this work, our results agree well with the experimental values. Table 1 Calculated vertical and adiabatic (T e ) excitation energies, dissociation energies (D e ) and vibrational energies (ω e ) for the I 2 X, A ′ , A, B ′′ and B states. Experimental values are indicated in parentheses. a Reference [26]. b Reference [27]. c Reference [28]. d Reference [29]. e Unbound in the calculated result. f Reference [30]. g Reference [31]. h Reference [32].

State
Vertical excitation energies (eV)

Molecular factor with the fixed-nuclei approximation
In this work, we consider the I 2 A 1 u ( 3 Π u ) and A ′ 2 u ( 3 Π u ) states as |e state and the X0 + g ( 1 Σ + g ) ground state as |g state, where electric dipole transition between |e and |g is forbidden or only weakly allowed. In order to inspect which I 2 electronic states contribute to the RENP rate as the intermediate |p state, we calculated C au (ω) at the fixed internuclear distance of 2.9Å, which is located between the equilibrium points of the X and A/A ′ states, without considering rovibrational energies. All calculated spin-orbit eigenstates are included as intermediate |p states in evaluating molecular factor C au (ω). In Fig. 3, calculated molecular factor C a au (ω) corresponding to Fig. 1 (a) is shown along with the similar molecular factor C b au (ω) corresponding to Fig. 1 (b), which has expression similar to Eq. (8) but with spin and dipole operators being swapped as We only need energies around 0-0.5 eV in evaluating the RENP rate, although the molecular factors are plotted up to 5 eV in Fig. 3 in order to inspect relative contribution of different intermediate electronic states. At low-energy region below 0.5 eV, the molecular factor for the E1×M1 process is larger than that for the M1×E1 process in both the A ′ → X and A → X cases. In other word, the first term in the bracket of Eq. (4) dominates over the second term, which suggests that only the first term should be retained for the RENP using I 2 molecule. In the plots, we can observe several spike-like structures which represent contributions of the intermediate electronic states as labeled in the figure. When the A 1 u ( 3 Π u ) state is used as the |e state, the B 0 + u ( 3 Π u ) state predominates among contributions from the other states. On the other hand, when the A ′ 2 u ( 3 Π u ) state is used as the |e state, the A 1 u ( 3 Π u ) and B ′′ 1 u ( 1 Π u ) electronic states are especially important at low energy region relevant to this work, while the other states have negligible contributions.
From here on, only the molecular factor in Eq. (8), which corresponds to the E1×M1 process of Fig. 1 (a), will be evaluated. For the A → X process, we will consider only the B electronic state as the intermediate electronic state. For the A ′ → X process, the A electronic state will be mainly considered as the intermediate electronic state, since the contribution of the A state looks larger than the B ′′ state in the low energy region. Validity of this approximation, inclusion of only the A electronic state in the A ′ → X process, will be inspected in the Appendix by comparing amplitudes evaluated by the A  to estimate the RENP rate. In this section, we evaluate the molecular factors considering vibrational levels on the A ′ , A and X electronic states for the A ′ → X process, and on the A, B and X electronic states for the A → X process. The ab initio energies of the X, A ′ , A and B electronic states, the same as shown in Fig. 2, were fitted using the functional form taken from Ref. [33]. The electric dipole transition moment between the X and A states and between the X and B states, the spin transition moment between the A ′ and A states and between the A and B states, shown in Fig. 4, were fitted by polynomial functions. The vibrational energies and wavefunctions on each electronic state were obtained by direct 10/24 diagonalization of the vibrational Hamiltonian using the pointwise coordinate representation of the wavefunctions, the discrete variable representation (DVR) method [34], with R min = 2.2Å, R max = 7.0Å and ∆R = 0.006Å. Summation of the intermediate vibrational states, on the A electronic state in the A ′ → X process, or on the B electronic state in the A → X process, was taken up to v = 40. These parameters for the DVR basis and the number of vibrational states in the summation were sufficient to obtain converged result.
In Fig. 5, the calculated molecular factors at ω = 0.4 eV for the A ′ → A → X process and for the the A → B → X process are shown as functions of the vibrational levels on the initial (A ′ or A) and the final (X) electronic states. In the A ′ → A → X process shown in the left panel in Fig. 5, the molecular factor is the largest around the point (A ′ (v = 0), X(v ′ = 20)). Starting from this point, the higher intensity region extends to the upper left direction. In the A → B → X process shown in the right panel in Fig. 5, the molecular factor is the largest around the point (A(v = 0), X(v ′ = 24)). In this case, the higher intensity region extends to the upper right direction. The structure of these higher intensity regions reflects the fact that the equilibrium points of the X and A ′ (or A) states are separated, and one of or both of the vibrational states on these electronic states should be sufficiently excited in order to achieve favorable overlap.
In Fig. 6, the molecular factors C a au (ω) for the A ′ → A → X and the A → B → X processes are shown as a function of photon energy w. The vibrational level on the initial state, the A ′ or A electronic state, was selected to be v = 0, and the vibrational level on the final X electronic state was selected to be v = 20 in the A ′ → A → X case, and v = 24 in the A → B → X case. As shown in the figure, the molecular factor for the A → B → X process is about 50 times larger than that for the A ′ → A → X process. Excitation energy (eV) Excitation energy (eV) In order to inspect the effect of atomic weight, or the strength of the spin-orbit couplings, we calculated the molecular factor C a au for the E1×M1 process for Cl 2 and Br 2 molecules with the fixed-nuclei approximation. The arrangements of the potential energy curves of Cl 2 11/24 and Br 2 are very similar to that of I 2 , as shown in Fig. 7. In both Cl 2 and Br 2 , the ground, first and second lowest excited states correspond to the X0 + g ( 1 Σ + g ), A ′ 2 u ( 3 Π u ) and A1 u ( 3 Π u ) states, respectively, as in the I 2 case. Also, the location of the B 0 + u ( 3 Π u ) state is similar to that of the B state in I 2 molecule. The procedure of the calculation is the same as in the I 2 case described in Sec. 3.1 and 3.2. The molecular factors were evaluated for two different pathways as in the I 2 molecule. In the first case, we took the A ′ 2 u ( 3 Π u ) state as the initial state and the X0 + g ( 1 Σ + g ) as the final state. In the second case, the A1 u ( 3 Π u ) state was taken as the initial state and the X0 + g ( 1 Σ + g ) was taken as the final state. All other 34 electronic states were considered in the summation of the intermediate state. We selected R = 2.2 and 2.5Å for Cl 2 and Br 2 , respectively, which are the middle points between the equilibrium points of the initial and the final states. In Fig. 8, the molecular factors for I 2 , Br 2 and Cl 2 are compared. The molecular factors for the A → X process are about 10 times larger than those for the A ′ → X process in all cases of I 2 , Br 2 and Cl 2 molecules. The magnitude of the molecular factor is the largest for I 2 , the smallest for Cl 2 in both the A ′ → X and A → X processes. These differences of magnitudes reflect the differences in strength of the spin-orbit couplings in these three molecules. We also investigated the lighter molecule, O 2 . In this case, the metastable c 1 Σ − u state was selected as the initial |e state, and the X 3 Σ − g state as the final |g state. Using the same procedure as we used for the I 2 , Br 2 and Cl 2 , the molecular factors for the E1×M1 and M1×E1 are calculated with the fixed-nuclei approximation. The molecular factor of O 2 for the E1×M1 process is approximately 10 −11 ∼ 10 −10 in the energy range ω = 1 ∼ 2 eV, where the A 3 Σ + u state has the dominant contribution to the virtual state summation. For the M1×E1 process, the molecular factor is about 10 −7 in the energy range ω = 1 ∼ 2 eV, where the 1 1 Π g state has the dominant contribution to the virtual state summation. As expected from its small spin-orbit couplings, the molecular factor for O 2 molecule is smaller 12/24 than the other I 2 , Br 2 and Cl 2 molecules: the magnitudes of the molecular factors in the A → X process are ∼ 10 −2 ,∼ 10 −3 and ∼ 10 −4 , respectively.

RENP spectral rate
A large macroscopic polarization necessary for significant RENP rates is developed by two trigger laser irradiation of frequencies, ω and ω ′ with ω + ω ′ = ǫ eg . RENP amplitude is proportional to the polarization g|(R 1 − iR 2 )/2|p averaged over intermediate states |p . Hence transition to a final single vibrational state X(v ′ ) for |g is selected out for the macro-coherently amplified RENP at each experimental setup.
It is appropriate before detailed presentation of numerical results to explain how experimentally neutrino properties and parameters may be determined. Suppose that two continuous wave trigger lasers of frequencies, ω + ω ′ = ǫ eg , are irradiated in counter-propagating directions and two exciting pulse lasers of frequencies, ω P − ω S = ǫ eg (in order to induce a Raman-type excitation |g → |e ) are suddenly switched on. Only during this pulse irradiation RENP occurs giving a unique signal; asymmetric directional increase of light output of lower frequency ω. If statistics allows, one may hope to measure parity-violating quantities such as emergence of circular polarization from linear polarization. Parity violating quantities that indicate unambiguous signal of involved weak interaction appear with smaller rates, at least by 100, than parity conserving quantities. Observations at different combinations of (ω, ω ′ ) provide spectrum rates at different ω's, thus covering some range of frequencies that gives the experimentally observed spectrum. After this spectrum determination, one compares the data with theoretical prediction computed assuming some values of undetermined neutrino parameters and properties, and finally the best fit with theoretical calculation determines these parameters.
The experimental method sketched here is by no means unique, and one can think of other schemes. Moreover, many simulations have to be done to determine the signal level once the best method of background rejection is found. It is suggested that the most serious background of the two photon process may be rejected by formation of condensed solitons [4], which are a target state of coherent macroscopic target entangled by static field condensates.
The quantity C a (ω)I(ω)/C 0 = C a au (ω)I(ω) in atomic units is plotted in the following figures except in the last three figures where absolute rates are illustrated for the target parameters of n = 10 21 cm −3 and V = 10 2 cm 3 . We used for calculation of the spectral rate numerical values of the mixing angles θ 12 , θ 13 , as determined in neutrino oscillation experiments Eq. (2), thus giving numerical weights of six thresholds in Table 2. (As usual, the conventional definition of oscillation angle factors c ij = cos θ ij , s ij = sin θ ij is used in this table.) We also used mass constraints given by oscillation data Eq. (2). The smallest mass defined by m 0 differs in the NH case where m 0 = m 1 (< m 2 < m 3 ) and in the IH case where m 0 = m 3 (< m 1 < m 2 ). For large m 0 , NH and IH differences are relatively small: for instance, in the NH case m 0 = 100 meV gives three neutrino masses, m 1 = 100, m 2 = 100.37, and m 3 = 111.8 meV's, while in the IH case m 0 = 100 meV gives m 3 = 100, m 1 = 111.8, and m 2 = 112.1 meV's, their mass range within ∼ 10 %. The mass pattern for a large m 0 has been called the quasi-degenerate case.
In all spectral rate figures we take as the initial state the I 2 A(v = 0) state and as the final state various vibrational states of X(v ′ ). It turns out that numerically dominant 13/24 = 20). However, the absolute rate is ∼ 50 times smaller than in the case of A(v = 0) → X(v ′ = 24). In Fig. 9 a global photon spectrum in the entire energy region is shown for the transition A(v = 0) → X(v ′ = 24), where four different cases are plotted, NH Dirac in black solid, NH Majorana of (α, β − δ) = (0, 0) (the CP conserving case) in red solid, IH Dirac in black dashed, and IH Majorana in red dashed (all with the smallest neutrino mass of 40 meV) [35]. All these cases appear nearly degenerate in the plot.
On the other hand, the enlarged spectrum in the threshold region is shown for smaller mass values of m 0 in Fig. 10. One can clearly observe three kinks in the NH case which are identified as photon energy thresholds of (11), (12), and (22). The other three threshold kinks, (13), (23), and (33), in the NH case are further to the left in this figure. This figure suggests a good chance of determining the smallest neutrino mass at the precision level of 1 meV, if one has a large statistics data in the threshold region.
We shall next examine the possibility of Majorana-Dirac distinction along with determination of CPV phases. The relevant CPV phases α, β, δ in the Majorana case appears in the 14/24 In the Dirac case there is no CPV phase dependence of the photon energy spectrum rate to this approximation. The parameter δ alone is accessible independently in neutrino oscillation experiments. These phases appear multiplicatively, in the rate formula around three thresholds of (12), (13) and (23), as cos 2α , cos 2(β − δ) , cos 2(α − β + δ) , which are further multiplied by the weight factor of Table 2 [4] times a product of two masses m i m j , i = j. There is thus no doubt that the Majorana-Dirac distinction is easier for larger neutrino masses. We shall introduce a new notation of CPV phase β ′ ≡ β − δ to simplify the formulas. The case (α, β ′ ) = (0, 0) corresponds to CP conserving (CPC) Majorana neutrino pair emission. As is evident in Fig. 11, the Majorana-Dirac distinction in the proposed de-excitation of I 2 appears easier than NH-IH distinction at low enough photon energies where rates are largest. Numerically, the I 2 CPC Majorana rate for m 0 = 40 meV is significantly different from the Dirac rate by ∼ 0.07 at a photon energy 0.37 eV, while this difference for Xe J = 2 transition [4] never exceeds 0.2 % at all photon energies. Study of the sensitivity to the CPV phase α, β ′ however requires a more careful analysis, because the CPC case of (α, β ′ ) = (0, 0) gives the largest destructive interference due to the effect of identical Majorana fermions, hence the Majorana-Dirac distinction is easiest in this case. It is thus necessary to vary the CPV phases (α, β ′ ) in their allowed range. It turns out that with the given numerical weight factors of Table 2 dependence of spectral rates on β ′ is much weaker than on α. The most prominent kinks are at the (12) and (33) thresholds and the highest sensitivity to α arises from the (12) threshold. Under this circumstance we may introduce a useful concept of Majorana α−band which is defined by the region of spectral rates, bounded by the largest Majorana rate at α = π/2 and the smallest rate at α = 0. We illustrate this band structure in Fig. 12. Strictly, the α−band is defined for β ′ = 0, but dependence of the band shape on β ′ is small and we may approximately use the α−band terminology for any value of β ′ .
On the other hand, variation of β ′ gives much smaller band widths, and moreover these bands are separated as α values vary, as shown in Fig. 13. It implies that experimental determination of β ′ is much more difficult than α in the proposed I 2 de-excitation scheme. We shall concentrate on determination of the CPV parameter α in the following.
The spectral rate for a large value of the smallest neutrino mass, 250 meV (a quasidegenerate mass not excluded by the cosmological bound [2]) is shown in Fig. 14, which shows a great sensitivity to the Majorana-Dirac distinction for this quasi-degenerate case. Another example of spectral rate is shown for m 0 = 100 meV in Fig. 15 where the Dirac case and the Majorana case of α = π/2 become nearly degenerate in the plot. Thus, for relatively smaller mass values of m 0 one needs more detailed study, both by theoretical means and 16/24 24). NH in solid and IH in dashed. Different colors correspond to Dirac and a few Majorana cases; Dirac in black, Majorana of (α, β) = (0, 0) in red, Majorana of (α, β) = (π/2, 0) in blue, and Majorana of (α, β) = (π/4, 0) in green. The smallest neutrino mass of 250 meV is assumed. The transition A(v = 0) → X(v ′ = 25) gives similar, but slightly smaller rate curves. numerical simulations, for determination of the α parameter, which is beyond the scope of the present work.
It is important to note that comparison of two de-excitation to X(v ′ = 24 and 25) offers a good tool of unambiguous experimental identification of RENP (this can be done by a change of one of two trigger lasers). This is an advantage of molecular targets, since one may compare two spectra of different transitions without much relying on the absolute rate scale. The zero photon energy limit of absolute spectrum gives an overall magnitude of the RENP rate, which then indicates relative difference of these two different transitions. We plot in Fig. 16 16 Zero energy limits of NH absolute spectrum rates (in the unit of 10 −11 Hz) are shown for two transitions to X(v ′ = 24) (in solid) and X(v ′ = 25) (in dashed). The Dirac cases are given in black, while the Majorana cases of (α, β) = (0, 0) are given in red and the (π/2, 0) cases are in green. The target parameters are n = 10 21 cm −3 and V = 10 2 cm 3 . ǫ eg = 0.787 eV was used for X(v ′ = 25).
is an interesting phenomenon of crossover of rates occurring at two points of the smallest mass m 0 ∼ 140 meV and ∼ 390 meV. For the range of α parameter between these two values Dirac and Majorana (π/2, 0) rates to X(v ′ = 25) of smaller atomic energy separation overtake Majorana (0, 0) rate to X(v ′ = 24) of larger separation. The simultaneous use of different vibrational transitions of molecules may be effective for RENP identification and exploration of the neutrino mass range as well.
The actual rate at each photon energy is equal to the plotted values C au (ω)I(ω) multiplied by the factor Γ dm C 0 ∼ 1.0 × 10 −7 Hz n 10 21 cm −3 Note that we took η ω (t) = 1 which is however expected to be smaller than 1.

18/24
Let us compare I 2 RENP rate with Xe rate. The I 2 rate is calculated as ∼ 5 nHz at ω = 0, assuming targets with n = 10 21 cm −3 , V = 10 2 cm 3 , and η ω = 1. Under the same experimental condition of targets, the overall RENP rate for Xe atom is much larger, by ∼ 10 5 . The origin of this large rate difference is explained as follows. Dependence of the RENP rate (for definiteness at ω = 0) is roughly given by product of factors as where we used d and S to mean the magnitudes of dipole and the spin factor in the RENP amplitude. Relation ǫ eg ≈ ǫ pg approximately holds both for Xe and I 2 . An extra n-dependence in addition to the basic macro-coherence dependence ∝ n 2 arises by our assumption of taking the maximal RENP rate (η ω = 1), with the stored field energy density of | E| 2 = ǫ eg n. Difference of ǫ eg between Xe and I 2 is ∼ 10 and the spin factor is of the same order of unity for two targets. The Franck-Condon like suppression factor [36] intrinsic to molecules is only 0.03 for I 2 . The rest of rate difference by ∼ 300 is attributed to difference in the transition dipole moment of Xe and I 2 . Unlike the simple dipole transition 6s → 5p in Xe atom the I 2 molecular dipole arises from B( 3 Π u ) → X( 1 Σ + g ) involving different configurations and is much suppressed. We definitely need solid environment of the target number density of O(10 23 ) cm −3 for realistic I 2 RENP experiments, which is a challenge left to experimentalists. A molecule closer to Xe with respect RENP appears to be CsI.

Summary
An attractive feature of using molecules is that there are many vibrational states available as a final molecular state for RENP which makes easier experimental identification of the process. We computed the RENP spectrum rate for homo-nuclear diatomic molecules of isovalent series, Cl 2 , Br 2 and I 2 and found that the rate becomes larger as the atomic number increases, as is expected. Even the largest I 2 rate is, however, much smaller than Xe atom rate. The distinction of Dirac-Majorana neutrinos and the determination of some range of the new CPV phase α are possible at photon energies where rates are largest for the quasi-degenerate neutrino masses. The NH-IH distinction along with the smallest mass measurement can be done around the threshold region. state has, when the effect of nuclear motion is considered. Here, we evaluate the amplitude of the first term in the bracket of Eq. (4), including nuclear motion, separately for the A and B ′′ states, and estimate their contributions to the molecular factor. The reason for not treating molecular factor C a au directly, but evaluating the first term in the bracket of Eq. (4) is that the B ′′ state is repulsive, see Fig. 2, and it is difficult to treat continuum nuclear wavefunctions in the expression of Eq. (8). The first term in the bracket of Eq. (4), on the other hand, can be easily evaluated using the following expression: whereĤ|p = E p |p is assumed. This equation represents the Green's function by integral in time domain [37], and can be evaluated by the Fourier transform of the correlation function g|de −it(Ĥ−iǫ) S|e obtained by solving time-dependent Schrödinger equation. For our purpose, the vibrational (nuclear) Hamiltonian on the A or B ′′ electronic state is used for H.
As a test, we evaluated the right hand side and the left hand side of Eq. (A1) separately for the I 2 A ′ → A → X process using explicit summation of the vibrational states, and the wavepacket propagation on the A state, respectively. The same DVR basis parameters were used as in Sec. 3.3. The wave packet was propagated by expanding the time evolution operator in terms of Chebyshev polynomials [38], with ∆t = 0.1 fs. Since many bound vibrational eigenstates contribute to the time evolution of the wavepacket, we need relatively long total propagation time, 7.0 ps. In addition, we had to leave small but finite value of ǫ, 10 −4 a.u., in order to converge the Fourier transform in the right hand side of Eq. (A1). As shown in Fig. A1, numerical values of the left hand side and right hand side of Eq. (A1) agree perfectly well in case of the A electronic state.  For the wavepacket calculation on the B ′′ state, the same DVR basis parameters were used. Since the B ′′ state is repulsive, we need to put the absorbing potential from R = 5.0 20/24 to 7.0Å to prevent artificial reflection of the wavepacket at the boundary. This absorbing potential can be regarded as decay of the wavepacket into R = ∞. The wave packet on the B ′′ potential energy curve with the specific initial state S|A ′ (v) , where A ′ (v) represents the vibrational state on the A ′ state, is propagated by expanding the time evolution operator in terms of Chebyshev polynomials [38], with ∆t = 0.1 fs. The total propagation time is 400 fs, which is much shorter than the time needed in the wave packet propagation on the A electronic state. This short propagation time can be attributed to the fact that the wave packet is not reflected, but just absorbed in 400 fs at the right-hand side of the B ′′ potential energy curve. The diploe and spin transition moments, necessary to evaluate Eq. (A1) are shown in Fig. A2.  The amplitudes for the A ′ → X process with the A and B ′′ intermediate electronic states, obtained by the wavepacket calculations, are compared in Fig. A3. When the initial vibrational state is A ′ (v = 0) and the final vibrational state is X(v ′ = 20), the amplitude of the A state at ω = 0.4 eV is about 60 % larger than that of the B ′′ state. This means that the molecular factor including the A intermediate state is approximately 2.5 times larger than that including the B ′′ intermediate state. Although the contribution of the B ′′ state is not negligible, for the purpose of the present study, the molecular factor may be safely approximated by including only the A electronic state. 22/24