Parity violating magnetization at neutrino pair emission using trivalent lanthanoid ions

A new detection method using magnetization generated at triggered radiative emission of neutrino pairs (RENP), $ |e \rangle \rightarrow | g \rangle + \gamma + \sum_{ij}\nu_i \bar{\nu}_j $ (atomic de-transition from state $|e \rangle $ to state $|g \rangle$ emitting sum of neutrino pairs $\sum_{ij}\nu_i\bar{\nu}_j$ accompanied by a photon $\gamma$), is investigated in order to determine unknown neutrino properties; absolute neutrino masses of $\nu_i$ and Majorana/Dirac distinction. Magnetization associated with RENP events has parity violating component intrinsic to weak interaction enforced by crystal field effect in solids, and greatly helps background rejection of quantum electrodynamic (QED) origin even when these backgrounds are amplified. In proposed experiments we prepare a coherently excited body of trivalent lanthanoid ions, Er$^{3+}$ (a best candidate ion so far found), doped in a transparent dielectric crystal. The magnetic moment $\mu \langle \vec{S}\cdot\vec{k} \rangle/k $ arising from generated electron spin $\vec{S}$ parallel to trigger photon direction $\vec{k}/k$ is parity odd, and is absent in QED processes. The generated magnetic field of order nano gauss is stored in crystals long after pair emission event till spin relaxation time. An improved calculation method of coherent rate and angular distribution of magnetization is developed in order to incorporate finite size effect of crystal target beyond the infinite size limit in previous calculations.


I. INTRODUCTION
Neutrino oscillation experiments, despite their remarkable discovery and measurement of finite neutrino mass and mixing in interactions, cannot determine remaining important neutrino properties and physical parameters; absolute neutrino masses and Majorana/Dirac distinction. These are key ingredients to explain the matter anti-matter imbalance in cosmology and construct the unified theory of particle physics. We proposed more than ten years ago an idea to solve these issues by using atoms/ions of typical eV energy release much closer to expected neutrino masses instead of atomic nuclei of typical MeV energy release. A mechanism of rate enhancement is crucial for weak processes, and we initiated experiments that verify the principle of enhanced event rates in weak QED (Quantum Electro Dynamics) process [1]. The enhancement region extends to a whole body of excited target consisting of many atoms/ions within which macro-coherence prevails. The spatial region of macro-coherence for multiple particle emission is much larger than the coherence region of single photon emission in related Dicke's super-radiance scheme [2], [3], which is limited by photon wavelength 1/k = λ/2π. Our achieved rate enhancement reaches of order 10 18 close to expectation in QED two-photon emission from excited vibrational state of para-hydrogen [4], [5], [6].
The process we have considered for neutrino study is either radiative neutrino pair νν emission (RENP) |e → |g + γ + νν stimulated by trigger laser, or Raman stimulated radiative neutrino pair emission (RANP) [7] in atomic de-excitation between two states, |e → |g . One measures either dependence on signal (stimulated by trigger) photon energy or angular distribution to extract neutrino properties. A large number of target atoms/ions close to those in solid environment is required for realistic RENP experiments. We shall consider a special crystal in the present work; trivalent lanthanoid ion Er 3+ of low concentration doped in host crystals, known to have narrow radiative widths.
A potentially serious problem against RENP detection is QED backgrounds. Measurement of parity odd quantity, as first discussed in [8], in triggered neutrino pair emission is of great use to distinguish the process involving weak interaction from purely QED processes, even if they are macro-coherently amplified. Random QED backgrounds that mimic parity violating effects are separated, since parity violating quantities in QED average to zero in statistically meaningful data. An interesting quantity of this nature in RENP is magnetization arising from magnetic moments parallel to the trigger photon direction generated at neutrino pair emission: the expectation value in the final state |f in RENP, gµ B f | S·k|f , wherek is the unit vector along the signal photon momentum k and S is the electron spin operator to be multiplied by gµ B , the g-factor times Bohr magneton for magnetization. This quantity is parity odd and time reversal (T) even. The generated magnetization due to these magnetic moments persists in crystals till spin relaxation time, of order msec or even larger in appropriate trivalent lanthanoid ions doped in host crystals. Measurement of an accumulated bulk quantity such as magnetization gives more freedom in experiments than direct event detection that rarely occurs [9].
Parity odd quantity emerges from interference term of parity even and odd amplitudes in the weak hamiltonian H W of neutrino pair emission. Writing in terms of the neutrino mixing matrix elements, U ei , i = 1, 2, 3 , ν e = i U ei ν i , the weak hamiltonian density consists of two terms of different parities, where ν i denotes a neutrino of definite mass m i and θ w is the weak mixing angle. Electron axial vector current couples to neutrino-pair emission current N α b , while electron vector current couples to pair emission current N α c . Parameters b ij , c ij are known except CP phases that appear in U ei , one phase in the Dirac neutrino case and three phases in the Majorana case. Our strategy of neutrino mass spectroscopy is to assume unknown value m 1 of smallest neutrino mass and CP phases for fixing all other variables from neutrino oscillation measurements, and calculate observable quantities for experimental comparison. This way one does not have to commit to any particular model beyond the standard theory, except finite neutrino masses and CP violation phases. If CP conservation is assumed, all b ij , c ij 's are real and already known, leaving m 1 as a single unknown parameter. Decomposition of flavor states ν e into neutrino mass eigenstates ν i is achieved in our proposed experiments by precision of laser frequencies, usually less than µeV, and not by measurable accuracy of photon energy in detectors that may be larger than 1 meV.
The dominant interference term of parity violation arises from product of spatial part of axial vector current of electron,ē γγ 5 e which is the spin operator S in the nonrelativistic limit of atomic electrons, and spatial part of vector current, the velocity operatorē γe ∼ v = p /m e (nearly equivalent to position operator r times energy difference between two involved atomic states, often used in atomic physics). Matrix elements of spin S are usually of order unity, while those of velocity operator v are less than 10 −3 , hence interference terms are at least smaller by 10 −3 than rate ∝ S 2 of a parity even quantity. We shall quantify this ratio for lanthanoid ion we use as a target. Magnetization caused by (ij) (ij is abbreviation for ν iνj ) neutrino-pair emission is proportional to neutrino mixing parameters, ℜ(b ij c ij ) (assuming CP conservation for simplicity), while parity conserving quantity of RENP event rate is proportional to ℜ(b 2 ij ). A great advantage of 4f electrons in lanthanoid ions is that they are not sensitive to environmental effects of host crystals due to filled 5s and 5p electrons in outer shells working as a shield. This circumstance gives both a narrow transition width and a large spin relaxation time, the latter being very important to magnetization measurement. Lanthanoid ions of low concentration doped in host crystals exhibit paramagnetic property at room temperature. For an experimental purpose as explained below, we shall use Kramers degenerate ions in which two states of different time reversal quantum numbers ± are energetically degenerate.
We develop in the present work an improved calculation method of macro-coherent rates and magnetization by incorporating finite size effect of crystal target in integration over neutrino momenta. In the previous works we took the infinite volume limit giving rise to the momentum conservation, or more properly the phase matching condition. But it is difficult to justify this infinite size limit for experimentally used size of targets, of order ≈ 1 cm. This improved method gives quantitatively different results of rate and magnetization, although the phase matching condition may be still useful as a technical guide. Angular distribution of magnetization, the most important physical quantity to neutrino mass spectroscopy, is greatly influenced by change of this calculation method.
The measurement method we propose here is different from those in previous proposals: detection of accumulated magnetization remaining in target medium [9] rather than measuring direct individual events that rarely occur. If the accumulated magnetization is above the sensitivity level of detectors, for instance, a high quality SQUID, then the method is found to be sensitive to absolute neutrino mass determination and to Majorana/Dirac distinction. Proposed experiments can be conducted irrespective of the nature of neutrino masses, Majorana or Dirac, and experiments themselves can determine whether neutrinos are of either type by measuring the shape and the magnitude of angular distribution of generated magnetization. The principle of this distinction is due to existence of interference term intrinsic to identical fermions (particle = anti-particle) in the Majorana neutrino case [10]. In the Dirac neutrino case anti-neutrinos are distinguishable, hence there is no interference.
The present paper is organized as follows. In the next section we first present amplitude of individual atomic process, and in the next second section we spell out improved calculation method of phase factors over entire target atoms incorporating finite size effect of target. In Section IV we explain trivalent Er ion used as a target candidate and present results of magnetization calculation. Effect of crystal field is quantified and used in calculation of atomic matrix elements. A host crystal was chosen from the point of rich available optical data necessary for detailed computations. In Summary and Outlook of Section V we mention points to be studied for realistic experimental design. In Appendix we explain four-level optical Bloch equation and its solutions necessary to estimate populations in energy levels and coherence parameters.
We use the natural unit of = c = 1 throughout the present work unless otherwise stated.
Feynman diagram based on old-fashioned non-relativistic perturbation theory of RENP |e → |g + γs + νν. Relevant phase factors including imprint to state |e given by the wave vector, k1 + k2, at two-laser excitation are attached to ion and particle states.

II. PROBABILITY AMPLITUDE OF RENP PROCESS AND PARITY VIOLATING MAGNETIZATION
A particular class of RENP process we consider in the present work consist of individual events occurring each at a space point x. Its diagram is shown in Fig(1). Individual atomic (or ionic) transition probability amplitude is given by [11] A ij e −iK0t+i K· x , Neutrino field operators ν i ,ν j in Introduction should be replaced here by their plane wave functions, but we shall keep the same notation for simplicity. We prepare excited state |e by two-photon laser excitation of their photon 4-momenta (ω i , k i ) , i = 1, 2 from the ground state |g . Emitted signal photon γ s is stimulated by trigger laser of photon 4-momentum (ω t , k t ) = (ω s , k s ).
, are 4-momenta of emitted neutrino pairs with their masses given by m i . The second order perturbation theory originally gives an energy denominator, 1/(−ǫ ep + E 1 + E 2 ), which is transformed as above by using the exact rule of energy conservation ǫ eg (= ω 1 + ω 2 ) = ω s + E 1 + E 2 , in order to eliminate neutrino energy dependence and make explicit dependence on the signal photon energy ω s . Two transition dipoles, magnetic ( µ pg ) and electric ( d pg ), are included in trigger field coupling ∝ E t , B t where laser intensity I t is given by the power Watt divided by an area, being related to field strength by | E t | = | B t | = √ I t . In Fig(2) we show momentum relation for target ion we shall consider in the present work. Three momenta, k 1 , k 2 , k s = k t , those for excitation and trigger, and we assume by experimental design that they are in a plane, while individual neutrino momenta, p i , i = 1, 2, may be out of this plane.
An implicit assumption for RENP amplitude is to be noted: emitted signal photons do not suffer from scattering or any reaction with surrounding atoms in solid after their emission. Even if there are interactions with surrounding atoms leaving modifications of amplitude formula, the subsequent formula for magnetization is related to RENP amplitude here, since magnetization is generated immediately after RENP and prior to interactions.
When the trigger laser frequency ω t = ω s is tuned to the resonance energy ǫ pg , the amplitude becomes A remarkable feature of this formula is that the RENP amplitude A ij is equal to neutrino-pair decay amplitude of |e → |p + ν iνj , simply multiplied by 2h r in which effect of stimulated photon emission is coded. Let us consider in more detail the factor h r multiplying the neutrino-pair emission amplitude, which may be related to various observable quantities. We consider unpolarized target and linear polarization of trigger laser consisting of an equal mixture of two circular polarizations. In this case In the case of trivalent Er ion scheme later proposed, this combination of parameters is numerically estimated as if one takes radiative decay rate, 57.3 sec −1 , for γ ep [12]. If an inhomogeneously broadened width of order 10 ∼ 100 MHz ×2π is taken for γ ep , this number decreases to 0.69 × (10 −3 ∼ 10 −4 ) I t /10mWcm −2 . We propose to use lowest Stark levels in each lanthanoid Jmanifold, since uncontrollable phonon relaxation is minimal for transitions between these levels. We may leave this uncertainty by inserting (57.3sec −1 /γ ep ) 2 in subsequent formulas of rate and magnetization. Neutrino variables, their helicities and momenta, are extremely difficult to measure, especially at these low energies, hence one integrates over these unobservables in atomic experiments. Neutrino helicity sum gives both for electron spin-current and for velocity current, [1] h1,h2 where dots · · · give irrelevant T-odd terms. If CP-odd phases of neutrino mixing matrix elements are taken into account, T-odd quantities such as M · k s × k 1 (up-down asymmetry of magnetization M relative to plane spanned by two wave vectors k s , k 1 ) become observable. We ne-glected T-odd terms for simplicity in the present work, assuming real mixing matrix. Electron matrix elements of (ij) neutrino pair emission are given by j e =ē( γ c ij − γγ 5 b ij )e with e,ē being replaced by relevant atomic wave function and conjugate of bound electrons states, |p , |e . δ M = 0 for Dirac neutrino and δ M = 1 for Majorana neutrino. Directional average over electron current j e may be taken for unpolarized targets, which gives the last term in eq.(10) of the form, p 1 · p 2 /(3E 1 E 2 ). Neutrino variable dependence in |A ij | 2 is then a function of p 1 · p 2 depending on the opening angle θ 12 and magnitudes p i , i = 1, 2.
Presence of crystal field is important to make it possible to measure parity violating effect of fundamental interaction. Both of RENP initial and final states, |e and |f , are eigenstates of QED interaction, including static Coulomb interaction of target ion with host crystal ions, which is a crystal field effect. Relevant parity violating quantity of our interest f |k · S e |f is transformed by parity operation P to − f |P −1k · S e P |f . Hence crystal field effect must produce parity mixture in the state |f ; |f = |f + + |f − such that f |k · S e |f = + f |k· S e |f − + − f |k· S e |f + in order to have non-vanishing value of f |k · S e |f = − f |P −1k · S e P |f .
The crystal field V c acting on a lanthanoid 4f n ion con-tains mixture of parity different states due to interaction with surrounding host crystal ions: by decomposing state vector into |4f = |4f 0 +δ|4f with |4f 0 defined as state without crystal field, one has as pointed out by [13]. This peculiar situation, which does not occur for isolated atoms in the free space, arises naturally when lanthanoid ions are placed in crystals. It is known that low-lying lanthanoid ions of 4f n system are ordinarily activated by magnetic dipole transitions, but can often simultaneously have forced electric dipole transitions of different parity, being caused by crystal field [13]: its calculation method was formulated in [14], [15].

III. IMPROVED CALCULATION METHOD OF MACRO-COHERENT RENP RATE
We re-examine and improve the calculation method of macro-coherent RENP rate and persistent magnetization, directly summing over plane-wave factors of relevant absorbed and emitted particles: eq.(4) summed over ion positions x. Suppose that one wants to determine energy and momentum to an accuracy level, (100 ∼ 1)µeV, whose inverse corresponds to 6.6 × (10 −12 ∼ 10 −10 ) sec in time and 2 ∼ 200 cm in length. Time scale given here is much smaller than the usual measurement time of atomic experiments, which implies that the Fermi golden rule of the infinite time limit should be valid and the energy conservation holds. The length scale, on the other hand, is larger than a typical spatial region of actual experiments, < 1 cm. Hence one cannot take the infinite size limit giving the momentum conservation of macrocoherent rates. What is achieved in the present improved method is incorporation of finite size effect by which it becomes possible to deal with intricate double resonance region of dynamical variables.
The squared amplitude relevant for rate and observables is decomposed into a sum of parity-even and parityodd terms. Differential and total rates arises from parityeven part, while magnetization discussed in the present work arises from parity-odd part. In the present section we discuss parity-even rate and in later sections we proceed to parity-odd contributions which shall be denoted by notations such as I P V . We investigate total event rate emerging from N = nV atoms or ions in a target volume of cylinder, V = πR 2 L, irradiated by excitation and trigger lasers. Event of an individual target atom/ion at a position x contributes to the total amplitude a piece proportional to product of plane waves of absorbed and emitted particles, e i K· x A ij ( K = sum of an imprint wave vector + relevant momenta of emitted particles. In solids momentum vectors should be multiplied by refractive indexes of host crystals denoted by n r .), with the atomic amplitude A ij independent of x position. We need to calculate an elementary integral of squared amplitude over neutrino momenta, for the triggered process |e → |g + γ s + νν. Here n is the target ion number density assumed constant. The averaged coherence |χ| may be estimated by solving an optical Bloch equation, as discussed in Appendix.
In the previous simplified treatment [1], [7], [9], we took the infinite volume limit which gives the space integral of the form, How this result is changed for a large, but finite volume is a subject of this section. For a target region of cylinder, radius R, length L, and its volume V = πR 2 L (having in mind R ≪ L = O(1) cm), the spatial integration of phase factors in eq.(12) gives where K ⊥ (K z ) is the magnitude of transverse (longitudinal) component to cylinder axis of excitation (parallel to the first excitation laser photon k 1 ). J 1 (z) (∼ z/2 as z → 0) is the Bessel function of the first order. Both functions U (x) , W (y) decrease as arguments |x| or |y| increases, showing damped oscillatory behaviors. It is difficult to accurately incorporate these damped oscillations in numerical integration over neutrino momenta, d 3 p 1 d 3 p 2 . We thus replace these by a smoothed out function, either Lorentzian function or Gaussian function. It was found that Lorentzian functions are easier to handle from the point of numerical computations. Difficulty of numerical simulations based on Gaussian approximation arises from its rapid decrease in large argument region. Small, but non-negligible contributions in signal light angular distribution come from a combination of one peak region and the other tail region of Gaussian functions, which greatly enhances computer tension.
The neutrino-pair phase space integral necessary for rate calculation of RENP process, |e → |g + γ s + νν, depends on how spin matrix element S ep (and smaller v ep , too) in the amplitude, eq.(7), is oriented. In order to simplify unnecessary complications, we shall assume for the sake of discussion here the case in which contributions proportional to S ep · p 1 S ep · p 2 averages out to give vanishing contribution. More realistic spin and velocity orientations shall be discussed in Section 4. Assuming the azimuthal symmetry, RENP rate is given by with n r the index of refraction of host crystal (∼ 1.45 in host crystal of Section IV). R ij is a function of (E i , p i ) , i = 1, 2 up to bi-linear orders, and shall be explicitly given below in eq.(56) for each neutrino (ij) pair. We introduced here notations for light directions: θ e refers to direction of the second excitation laser relative to the first excitation laser, while θ s (equal to the trigger one θ t ) does to signal photon direction relative to the first excitation laser. Lorentzian parameters x l , y t were determined by requiring that Lorentzian functions coincide with exact plane-wave integrals of U (x) = (sin x/x) 2 , W (y) = (2J 1 (y)/y) 2 in values at x = 0 , y = 0 and integrated values. All angles are measured from zaxis of the first excitation laser. For large size values of L , R, two Loretzian functions, U 0 (x) , W 0 (y), have sharp peaks in neutrino variables at x = 0 , y = 0, which restricts the region of largest contribution in the neutrino phase space integral given by three independent variables, p 1 , θ 1 , θ 2 . The largest region is given by K = 0 where K = k 1 + k 2 − p 1 − p 2 − k s , and momentum vector relation is depicted in Fig(2). Thresholds of neutrino-pair production are given by taking vanishing neutrino momenta p i = 0. The threshold condition gives w t = 0 , w l = 0, or This determines both θ s , θ e : in resonant trivalent Er scheme, energies ω i , ω s are all fixed, and two angles are determined as cos θ e = −0.7838 , cos θ s = 0.5182 (142 • and 59 • ). Thus, introduction of non-aligned second excitation laser of θ e = 0, π is required to experimentally approach neutrino-pair thresholds which are most sensitive to Majorana/Dirac distinction.
An important technical question concerns how the infinite volume limit leading to the momentum conservation is approached, ultimately giving a method of how to determine the size appropriate for realistic crystal growth. This and related problems are postponed and discussed after we solve how to compute magnetization in Section IV.

IV. TRIVALENT ER ION SCHEME
In the present work we shall restrict to cases of neutrino pair emission, in which two-photon excitation from the ground state to excited state |e is followed by triggered RENP: |g + γ 1 + γ 2 → |e ; |e → |g + γ s + νν. Experimental layout is shown in Fig(3). Excitation laser and trigger laser frequencies are matched to level spacings of trivalent Er ion as shown in Fig(4).
Typical candidate targets we may consider are trivalent Er and Ho ions which have odd (11) and even (10) number of 4f valence electrons, respectively. Trivalent 4f 11 Er ion has Kramers degeneracy, degeneracy due to time reversal invariance, while Ho ion does not have this degeneracy automatically. Symmetry of host crystal may however induce accidental Kramers degeneracy in nonsinglets. This occurs in the ground state of Ho 3+ in a special host crystal of higher symmetry [16].
Lantanoid ions in the free space, or vacuum, are classified using the LS coupling scheme, giving multiplets 2S+1 L J in which J is the total angular momentum combining the total spin S and the total orbital L angular momentum. In the ground state of trivalent Er ion, S = 3/2 and L = 6 , J = 15/2, and is denoted by a term notation 4 I 15/2 . These multiplets, called J-manifolds, are split into Stark states under crystal field of host crystals, as illustrated in Fig(5). Typical J-manifold splitting is of order 0.5 ∼ 1 eV, while Stark splitting is of order several tens of meV [17]. The number of separated Stark states in J-manifold is J + 1/2 for half-odd J, each forming Kramers doublets, whose degeneracy may be lifted, with amount of order 60µeVB/Tesla, by application of magnetic field. Relaxation due to phonons in solids, which might become serious obstacles against RENP, is the smallest at the lowest Stark state in each J-manifold, hence we shall use the lowest Stark states of J-manifold in optical transitions of our scheme.
The choice of host crystals is important for success of RENP project. We choose LiYF 4 (usually called YLF) crystal as host. This crystal has a point group symmetry S 4 at the site of F, and target ion Er 3+ replaces a fraction of F sites. Energy levels of Stark states in doped crystal Er 3+ :YLF are well known, but we further need magnetic and electric dipole transition rates denoted by MD and ED. For these we have to use available theoretical calculation of different, but similar host crystal Er 3+ :Cs 2 NaYF 6 . According to [18], calculated radiative decay rates of 10% doped Er 3+ in host Cs 2 NaYF 6 and other data are given in Table 1 Coexistence of magnetic and electric dipole transitions, a necessary ingredient for magnetization measurement in RENP, is a manifestation of forced electric dipole moment in lanthanoid crystals [13], [14], [15]. The quantity v/S in the last column contributes to ratio of magneti-zation to rate which shall be precisely defined and calculated later.
Our RENP scheme of trivalent Er ion uses the following J-manifolds as in Fig(4 The expected position of neutrino pair thresholds calculated from eq.(23) appears at signal (=trigger) direction given by corresponding to ±1.026(±59 • ) radians. The second excitation laser should be irradiated with aligned angle at corresponding to ±2.47(±140 • ) radians.
Optical and neutrino-pair emission processes in solids are governed by the selection rule of time reversal (T) quantum number: its change must be T-odd, since two major single-photon operators, electric dipole (E1) p/m e and magnetic dipole (M1) µ, and neutrino-pair emission operators, S and v, are all T-odd. (Use of the length gauge e r for E1 obscures this nature of T-oddness.) On the other hand, two-photon excitation is governed by Teven operator, hence |e and |g are both T-even or both T-odd [19]. The intermediate state |p that participates in neutrino-pair emission has different T quantum number from these.
States constructed in the LS coupling scheme are transformed under time reversal to In the leading order of our scheme, (−1) J+MJ = 1. The second terms ∝ b J ′ are sub-leading compared to ∝ a J terms. States of T = ± are constructed as Matrix elements of T-odd operator V are then given by to leading orders. The first term dominates over the second term due to smaller ∆M J . Crystal field acting on trivalent Er ion (sum of Coulomb interaction caused by surrounding host ions) splits J-manifold states into J + 1/2 degenerate Kramers doublets, consisting of T-even and -odd linear combinations of MJ c MJ >0 (|J, M J ±(−1) J+MJ (|J, −M J ) with quantization axis taken parallel to c-axis of tetragonal host crystal. Magnetic quantum numbers of lowest Stark levels for relevant RENP J-manifolds are mixtures of different M J values, and their weights are calculated using tabulated crystal field parameters [20]. Our calculation shows that dominant components have probabilities, For this choice of M J = 1/2 , M ′ J = 3/2, other spin components, matrix elements of S 0 , S − vanish. This means that vector matrix elements are parallel to crystal c-axis.
For estimate of v ep = ±i d ep ǫ ep /e we need to consider the crystal field effect of host crystals, since parity forbids dipole elements between the same parity J-states. The leading effect of host crystal surrounding target ion is to introduce Coulomb-induced potential derivative of the form, with a the lattice constant. Thus, a direct product of two vector operators for instance, is non-vanishing, where ∆E is a typical energy difference of 5d-4f ionic levels (roughly of order 10 eV). We defined components of vector operator V parallel to quantization axis as V 0 , and perpendicular components as V ± , depending on their angular momentum. In Judd-Ofelt theory higher order terms O(r 2n+1 ) in crystal field expansion and higher levels than 5d are effectively taken into account. RENP amplitude is finally given in terms of reduced matrix elements which are related to observables. The following general formula for calculation of reduced matrix element of a vector operator is useful: For application to the spin operator, we first note a relation in LS coupling scheme, These give S 2 ep ≈ 0.0157 for lowest Stark states in our scheme.
In order to work out parity violating interference term, we start from squared amplitude when neutrino helicities are specified (not summed over), although neutrino helicities cannot be measured. This calculation gives more insight than providing helicity summed squared amplitude from the outset. Neutrino pair emission of two planewave modes specified by their momenta and helicities, ( p 1 h 1 ) , ( p 2 h 2 ), 1 and 2 referring to anti-neutrino (distinguishable from neutrino only for Dirac case) and neutrino variables, respectively, gives interference term of squared amplitude at |e → |p [22] ( where θ 12 is the opening angle of emitted neutrino-pair, andp i is neutrino propagation vector of unit length. Two important limiting cases of quantity (I P V ) 12 are listed in the following tables; the case of relativistic (R) neutrinos and the case of non-relativistic (NR) neutrinos. The relativistic case gives dominant helicity combination of left-handed neutrino and right-handed anti-neutrino (opposite helicity combination for Majorana neutrino pair), all other combinations giving vanishing contribution. The non-relativistic case gives a different picture of nearly all non-vanishing helicity combination, as shown in the table.  Table 2 Contribution to (I P V ) 12 of neutrino helicity combinations: relativistic case. The following definition was used; Dirac 3 4 (Majorana 1) 0(0) Table 3 Contribution to (I P V ) 12 of neutrino helicity combinations: non-relativistic case. The factor δ M = 0 for Dirac neutrino and = 1 for Majorana neutrino, assuming an equal mass pair of m i = m , i = 1, 2.
The helicity-summed interference term relevant to actual experiments is given by h1,h2 This result is equivalent to eq.(10) in the preceding section. Neutrino pair phase space integration over momenta of (I P V ) 12 further gives PV interference contribution, This gives contribution from a single neutrino-pair, and one should further add all six pair contributions with appropriate weights. On the other hand, RENP rate is given by parity conserving quantity taking dominant spin current contribution,F Details of RENP rate and magnetization depend on how crystals are placed relative to directions of excitation and trigger laser irradiation. Crystal c-axis defines quantization axis and matrix element S ep is oriented along this axis due to the selection rule, J ′ , M + 1|S i |J, M = 0 unless S i = S + . For simplicity we take direction of the first laser excitation k 1 , either parallel or orthogonal to c-axis. There are two important cases of interest: (1) c-axis is parallel to k 1 , (2) c-axis is not parallel to k 1 , but in the plane spanned by two vectors, k 1 , k 2 , (3) c-axis orthogonal to the plane spanned by these two vectors. We assume the signal wave-vector k s in k 1 , k 2 plane. In the neutrino phase space integral ofJ ij 2 for case (3) p i · S ep and p i · v ep give odd functions to integrands, leading to vanishingJ ij 2 contribution. HenceF P V =J 1 v ep · S ep and F P C =J 1 S 2 ep . On the other hand, in case (1) where both v ep , S ep are parallel to the first excitation axis ∝ k 1 . The quantity R ij of eq.(17) is thus RENP rates are numerically We turn to parity violating magnetization. Amount of generated magnetization by neutrino pair emission is parity violating amplitude multiplied by magnetic moment in the state |f , hence its magnitude is We use the relevant component of g-factor, g ⊥ of 5.9 for Er 3+ in the first excited J-manifold 4 I 13/2 Γ 56 of trivalent Er ion [23]. The number density of Er ions in 0.1 % Er doped crystal (YLF) is 1.4×10 19 cm −3 , which we shall use as a reference value. We further introduce a conversion factor to magnetic field from magnetic moment, ξ, since magnetization is measured outside the target, which gives µ eff = ξgµ B with ξ is related to how near the target one measures magnetic field. To convert the generated magnetic moment µ to the magnetic field caused by these moments, we use the formula of magnetic field due to magnetic moment µ, to be integrated over the target cylinder region (n( r ′ ) is the number density of RENP affected ions). When the measurement site is chosen to be close to the target region, the magnetic field generated by RENP may be effectively expressed as ξµn, ξ being a fraction of RENP affected ions ≈ 0.01 for reference parameters. Actual value of ξ to be taken in magnetization depends on how experiment is done. Magnetization generation rate caused by neutrino pair emission is calculated using eq.(58) with k · S = cos θ s (signal photon directional factor) and replacementS 2 ep in rate by v ·S ep for magnetization. Weight of individual pairs is also changed to b ij c ij for parity-odd magnetization. These lead to The angular cosine cos θ s arises from k · S ep /k ∝ cos θ s . Using | v ep || S ep | = 0.0157 × 2.49 × 10 −6 /3.47 = 1.1 × 10 −8 , calculated magnetization at neutrino pair emission is numerically The number in front is for 1000/57 msec of 1/γ ep , and if the value γ ep = 2π× 10 MHz is taken, it is changed to 1 × 10 −3 . A typical value for computed ij G ij is of order 10 −8 eV 5 , implying generated magnetization of order, 10 5 |χ| 2 nG sec −1 L = 2R = 2cm. See below on |χ| 2 .
We work out results of magnetization angular distribution for case (1); the case of c-axis parallel to excitation laser wave vector k 1 . Rate is proportional to F ij ∝ b 2 ij , while magnetization is to cos θ s G ij ∝ b ij c ij , assuming CP conservation of all vanishing phase δ = α = β = 0. Different weights in rate and magnetization distributions in six neutrino-pair thresholds are calculable from neutrino oscillation data [24], and they are listed in the following table, Table 4, assuming CP conservation case described by real number b ij , c ij .  Table 4 Relative weights of neutrino-pair thresholds as determined by neutrino oscillation data, assuming CP conservation.
We now address the postponed question of how the size or volume of crystal targets affects the angular distribution, and wish to estimate an optimized target crystal size for experiments. For definiteness we take the trivalent Er scheme identified by eq.(25) ∼ (29), and focus on excited region of either cylinder or disc shape whose boundary lies near the aspect ratio of R/L = 1/2. It is always easy to experimentally study excited region of smaller aspect ratio, 2R < L by laser focusing, hence we concentrate to the aspect ratio 1/2, V = πL 3 /4, for subsequent study. We may classify small, medium and large crystals with varying volumes, taken here in a range L = 0.5 ∼ 5cm, since much smaller and much larger crystal sizes have problems either of bad angular resolution or of crystal growth in actual experiments. Macro-coherent rate (magnetization) is proportional to angular function, given by a unit of eV 5 cm 6 .
We first plot magnetization strength given by ij G ij (θ s ) in Fig(6), in two cases of small and intermediate crystal sizes, L = 0.5 and 2 cm's. The unit of angular resolution in this figure is 0.36 • covering a region of 21.6 • centered at 59 • . The resolution 0.36 • seems smaller than experimentally attainable angular resolution, hence it is better for comparison with experimentally obtained data to theoretically average over discretized angles by a smoothing procedure. We adopt here a sum over four adjacent angular points (1.44 • ), taking two overlapping angular 0.72 • points with neighboring regions. This gives result of Fig(6b) corresponding to result of Fig(6a). As expected, smoothed angular distribution does not show complicated wiggling structures. We confirmed that central peaking increases as the crystal size increases, to 5 cm ∼ 10 cm, approaching the infinite volume limit.
This kind of smoothing should be incorporated in detailed design study of experiments, but there are other factors such as spatial laser profile defining the excited target region and angular resolution of measuring device of magnetization, and consideration of these is outside the present scope of work. We shall therefore present numerical results without smoothing in subsequent results. For simplicity we give results for L = 2R = 2 cm size of crystal in subsequent figures.
The first physics issue we would like to discuss is neutrino-pair opening angle (θ 2ν = θ 1 − θ 2 ) distribution given byF P C , integrated over θ i , i = 1, 2 with the constraint of fixed θ 2ν around π, which is illustrated in Fig(7). It shows a remarkable back-to-back configuration, consistent with the momentum conservation. We assumed a single neutrino-pair of mass as small as 10 meV.
Decomposition into individual (ij) neutrino-pair contribution may be illuminating, although the task is not rewarded in experimental analysis. Result is illustrated in Fig(8) showing different angular structures for different pair contributions. Note weight factors given in b ij c ij row of Table 4, showing two largest (12) and (33) pair contributions.
Magnetization angular distributions given by cos θ s ij G ij (θ s ) are illustrated for Dirac NH cases of three smallest neutrino masses, 10, 30, and 50 meV's in Fig(9a) and Fig(9b). These mass differences are expected to be distinguishable in experimental analysis.
Majorana-Dirac distinction in terms of angular distribution is not easy as shown in Fig(10), which gives Dirac and Majorana functions of ij G ij (θ s ), taking a (eV 5 ) ij G ij n θ [θ s = 1.026 + n θ π/500 (rad)] n θ [θ s = 1.026 + n θ π/500 (rad)] NH smallest mass as large as 600 meV. Ratio of angular point sums is 1.06 for m 1 =500 meV and 1.11 for 600 meV (Dirac values larger than Majorana ones). In order to detect Majorana/Dirac distinction for smaller neutrino masses, it is better to use hyperfine levels of odd Er isotopes along with magnetic field, which also requires a joint microwave spectroscopy. This subject is beyond the present scope of work. Let us discuss the remaining important factor |χ| 2 . We assume CW (continuous wave) laser operation, and CW trigger laser is irradiated with a time lag after two-photon excitation. RENP events start to occur at trigger irradiation, and generated magnetization remains after events till spin relaxation time comes. Time evolution of popu- n θ [θ s = 1.026 + n θ π/60 (rad)] n θ [θ s = 1.026 + n θ π/300 (rad)] lation ρ ee in the upper excited state |e and coherence ρ ep between levels of neutrino-pair emission is described by four-level optical Bloch equation, as given in Appendix.
The coherence ρ eg and population ρ ee develop prior to trigger irradiation, and the coherence is partially transfered to ρ ep after trigger irradiation. Time dependent |χ| 2 is identified to be |ρ ee ρ ep | 2 . Generation rates of magnetization are proportional to n|χ| 2 ∝ ρ 3 ee |ρ ep | 2 , which has a time structure, as illustrated in Fig(11). Note that event rate has a different power dependence ρ 2 ee |ρ ep | 2 . In our example of numerical computation, generated value ρ 3 ee |ρ ep | 2 becomes of order, 10 −2 ∼ 10 −3 . It is important to note that this large value is realized for weak intensity I 1 , I 2 , I t of three CW lasers, which however requires a small relaxation rate γ ep . The largest RENP event occurs right after the trigger irradiation, accompanied by side events, which is somewhat similar to ringing phenomena of super-radiance events [3]. As stressed, signal of RENP events persists as magnetization till spin relaxation time, hence time-accumulated events show increasing profile of magnetization with time as in Fig(11). We confirmed that laser turn off leads to gradual decrease of ρ 3 ee |ρ ep | 2 . In actual experiments time-integrated magnetization till spin relaxation time is expected to be a readily observable quantity.
Our calculated magnetization rate per second is applicable at cycle end of three laser operation, and cycles are repeated with a cycle repetition time to be taken ideally of order spin relaxation time. Thus, the quoted rate per second so far given should be multiplied by the number of repetition times per second. If the repetition cycle is 100 per second, the effective rate per second would be increased by 100. It would be nice if one can make dedicated simulations including both spin relaxation and coherence development at the same time. For this simulation data of spin relaxation is needed. This is left to future works.

V. QED BACKGROUND
Major obstacle of QED background is the problem of extinction of initial ions in state |e for RENP. We shall discuss amplified processes first and at the end of this section we discuss effect of spontaneous electric dipole transition to the ground state.
We consider potentially dominant amplified QED event of two photon emission, |e → |g + γ s + γ b via |p , with extra background photon γ b of energy ǫ ep . Three vectors, k 1 , k 2 , − k s , are arranged to be in the same plane, hence four momentum vectors, k 1 , k 2 , − k s , − k b , form closed quadrangle to leading order of momentum conservation. One can directly prove that there is no extra photon vector k b satisfying the quadrangle relation, using parameter sets of proposed scheme. Still, the improved method of plane-wave summation in this work gives a finite contribution. Our estimate of finite size effect gives extinction rate per ion of order 3 × 10 −4 , tak-Generated and accumulated magnetization FIG. 11. Time structure of magnetization generation rates ∝ ρ 3 ee |ρep| 2 (not to be confused with RENP event rate ∝ ρ 2 ee |ρep| 2 ) in solid black, when trigger is irradiated at 0.2 × lifetime of state 2 H 11/2 with time lag after irradiation of two-photon excitation lasers. Also shown is time profile of accumulated magnetization (ten times values of directly integrated result) in dashed red, which remain in crystals within this time scale. The same set of parameters as in Fig(12) is used here, with the green curve indicating trigger laser time profile.
ing inhomogeneous width γ e to be 1 MHz (effect ∝ γ −2 e ). This rate is sufficiently small to ignore the process.
The next sub-dominant contribution is amplified threephoton emission, |e → |g +γ s +γ b +γ ′ b with γ b +γ b ′ = ǫ ep . In order to suppress this process, we resort to symmetry and apply a weak magnetic field of order 100 G to distinguish by irradiated laser frequency two Zeeman states of separated doublet components. Zeeman split levels are given in terms of states of definite time reversal quantum numbers ±, Optical transitions between two states of |e , |p given by |J, ±J z and |J ′ , ∓J ′ z , proceed via either single photon T-odd operator V o or two-photon T-even operator V e via paths of different energy, |J, ±J z → |J ′ , ∓J ′ z for T-odd case, |J, ±J z → |J ′ , ±J ′ z for T-even case. This situation occurs due to cancellation of two terms, for instance between |J, J z , ± → |J ′ , J ′ z , ± for T-odd case. One can thus single out T-odd transition by laser frequency, rejecting two photon emission. On the other hand, when Kramers degeneracy exists in the zero field, both T-even and T-odd transitions occur since lasers cannot distinguish two modes of transitions.
This leaves amplified four-photon emission, |e → |g + γ s + γ 1 + γ 2 + γ 3 , of QED origin as the main extinction process. Estimate of total rate of this process is as follows. Except a common factor related to signal photon part, the major difference arises from coupling factors and phase space integration either over two-body or three-body states. Couping factors of QED background is of order, (dipole 2 /3) 3 ≈ 10 −39 eV −6 , while RENP coupling is of order G 2 F ≈ 10 −46 eV −4 . Three-body phase space of background QED photons have where energy denominators ∆ i , i = 1, 2 that appear in perturbation theory are linear combinations of level energy differences ǫ ab and one of emitted photon energies ω i , i = 1, 2, 3, with∆ i its appropriate average value, while two-neutrino phase space has Thus, the ratio of three-to two-body phase space is of order, 3 × 10 −7 (ǫ 6 ep /∆ 2 1∆ 2 2 ), which gives roughly comparable magnitude of amplified QED background to RENP events. Even if the four-photon process depletes to some extent the excited state |e , parity violating magnetization can single out RENP process.
Finally, we discuss how magnetization may be measured under a spontaneous electric dipole transition from state |e to the ground J-manifold of lifetime 1/518 sec (from Table 1). Despite of small branching ratio of order 10 −12 (proportional to laser power, taken here 10 mW cm −2 , and other factors, hence controllable to some extent) RENP events occur with reasonably large total event number (assuming reference target ion density and target size) during lifetime of |e . Accumulated magnetization is measurable by taking sufficient measurement time with a fast repetition cycle of laser excitation. Although leakage due to spontaneous decay occurs to different time reversal state, this is not a problem since CW two-excitation to |e followed by another spontaneous emission takes the excited state back the original path.

VI. SUMMARY AND OUTLOOK
Even if the experimentally confirmed macro-coherence amplification [4], [5], [6], works, the major problem for a successful RENP measurement of process |e → |g + γ s + νν from an excited state |e is how to reject similarly amplified QED backgrounds. We proposed in this paper a method of overcoming this difficulty by measuring parity violating observable present in weak interaction and absent in QED. The proposed magnetization along the signal photon direction is parity violating, and emerges with reasonably large magnitudes, if one adopts a clever way of choosing excitation and trigger directions and their photon energies. The best candidate target we have found so far is trivalent lanthanoid ion of odd number of 4f electrons: Er 3+ doped in host crystals such as YLF. Host crystals should provide to the ion a sharp radiative width and a large spin relaxation time. Moreover, the Kramers degeneracy of odd number 4f systems such as Er 3+ brings about an experimental bonus. Excitation at |g → |e is done by T(time reversal)-even twophoton process to imprint a phase to target ions, while macro-coherence generation between states of neutrino pair emission is achieved via single photon T-odd trigger laser irradiation with a time lag. Numerical simulations using four-level optical Bloch equation supports a sizable generation of coherence and population among relevant states.
Calculation of magnetization suffers from uncertainty, most notably due to a lack of optical experimental data of separated magnetic and electric dipole transitions between two relevant states of neutrino pair emission.
Coexistence of magnetic and electric dipole transitions is a key element for parity-violating magnetization measurement. Crystal field of surrounding host ions introduces an environmental parity violation, a necessary ingredient for a successful measurement of RENP originated magnetization. Generated magnetization may become of order 10 2 nG sec −1 for 0.1% concentration of trivalent Er ion, depending on coherence factor generated at excitation. Numerical analysis of RENP rate and magnetization have been carried out by incorporating finite size effect of target crystal, which is important due to realistic crystal sizes of order 1 cm hopefully used in experiments.
If the angular resolution of magnetization measurement is reasonably good, of order a few degrees, magnetization has a good sensitivity to the smallest neutrino mass down to of order 20 meV. Majorana/Dirac distinction is difficult from measurement of angular distribution of magnetization in the proposed Er 3+ scheme. One possibility for improvement is to use another crystal such as Ho 3+ :YLF.
We did not study sensitivity to CP violation parameters, but this may open a new possibility of Majorana/Dirac determination assuming presence of CPV phases, because finite CPV phases may enhance Majorana/Dirac difference term proportional to m i m j ℜ(b ij c ij ). There are two ways to determine CP violation parameters in the mixing matrix: one is to incorporate CP phases in time reversal even quantities such as rate and magnetization rate, which is a straightforward extension of our formalism here. This gives even functions of CP phases. The other is to measure time reversal violation directly and, for this purpose, to calculate time reversal odd measurable which is given by odd function of CP phases. These studies are left to future work.
A rich variety of angular structures resulting from incorporation of finite size effect as seen in our figures is not really what we can precisely measure in actual experiments, since detector angular resolutions, in particular for magnetization measurement using SQUID, should be taken into account. We may expect to observe smoothed out distributions, but details of the smooth-out should be worked out taking into account realistic angular resolution of adopted detector. It is rather better to first set up experimental goals for neutrino parameters, and to design and choose detector system matching this objective by detailed simulations using the methods developed in the present work. In the meantime it is desirable to measure electric and magnetic dipole transition rates between two relevant ion states in crystals such as Er 3+ :YLF.

VII. APPENDIX: OPTICAL BLOCH EQUATION FOR ESTIMATE OF POPULATION AND COHERENCE
Before we present four-level optical Bloch equation and its numerical simulations, we give an intuitive picture of what may be realized in ideal situations. Suppose that trigger laser irradiation is followed by two-laser excitation after a time lag. All lasers are assumed to be operated by CW (continuous wave) mode, which effectively means that lasers are irradiated during the entire lifetime scale. One may take Rabi frequencies much larger than radiative decay rates, Ω ab = d 2 ab + µ 2 ab E ≫ γ ab , using a reasonable range of laser intensity, I t = E 2 = B 2 .
CW steady state solution derived by analysis of laddertype three-level optical Bloch equation after various lifetimes has the following form of population and coherence, in the zero detuning limit, [25], [26], It is thus possible to generate a maximal coherence, |ρ eg | = 1/2 between |e and |g , with population ratio of 1 : 1 : 0 between three levels |e , |g , |q when Ω qg = Ω eq . A nearly complete inverted population to |e is also possible when Ω qg ≫ Ω eq , the inequality sometimes felt counter-intuitive. Suppose that trigger laser matched to energy difference of |g and |p is irradiated with a time lag when the steady state caused by two-photon excitation is being formed. One can then effectively consider three level Vscheme in which two lasers irradiated between |g , |e are replaced by an effective single laser, forming V-type laser irradiation along with trigger |g → |p . The steady-state solution of V-type is given by [26] One can thus generate a large coherence between states, |p , |e , provided that Rabi frequencies, an effective Ω eg and Ω pg , are of comparable magnitudes. For instance, ρ ep = 1/4 for Ω pg = Ω eg . This is important, because neutrino pairs are emitted between these states and large coherence is required for ρ ep . The picture here is a kind of coherence transfer from ρ eg at excitation to ρ ep at trigger. We would like to check whether this picture is realized in four-level problem under three laser excitation advocated in the text.
We are bound to use for relaxation radiative decay rates due to lack of inhomogeneous widths in actual host crystals. This might be a proper choice when we use lowest Stark levels, but experimental data of relaxation rates of lowest Stark states in each J-manifold are needed to settle this question. Time unit used in simulations is the lifetime 1.7 msec of |e =Er 3+ 2 H 11/2 . Calculated Rabi frequencies Ω ab = d ab E = 2πν ab under CW operation of powers ∼ 10mWcm −2 = 1.26 × 10 −3 eV 2 are of orders, using calculated A-coefficients in the text. Divided by respective radiative decay rates, these Ω/γ are corresponds to use of CW laser power of order µW cm −2 or less during a few msec irradiation. Solution of large coherence ρ ep generation is illustrated in Fig(12), which shows a kind of coherence transfer at the trigger irradiation, as suggested by heuristic argument using three-level optical Bloch equation. Thus, there are ranges of laser parameter that give rise to sufficiently large coherence and population.
RENP magnetization emerges with time profile as illustrated in Fig(11) of the text, since rates are proportional to ρ 3 ee |ρ ep | 2 .