Nuclear matrix elements for the resonant neutrinoless double electron capture

The rate of the neutrinoless double electron capture ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\ensuremath 0\nu$\end{document} ECEC) decay with a resonance condition depends sensitively on the mass difference between the initial and final nuclei of decay. This is where the JYFLTRAP Penning-trap measurements at the JYFL become invaluable in estimation of the half-lives of these decays. In this work the resonant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\ensuremath 0\nu$\end{document} ECEC decay is discussed from the point of view of its theoretical aspects, in particular regarding the resonance condition and the involved nuclear matrix elements (NME). The associated decay amplitudes are derived and the calculations of the NMEs by the microscopic many-body approach of the multiple-commutator model are outlined. The resonant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\ensuremath 0\nu$\end{document} ECEC decays of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\ensuremath {\rm ^{74}Ge}$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\ensuremath {\rm ^{136}Ce}$\end{document} are discussed as applications of the theory framework.


Introduction
At present the neutrinoless double-beta (0νββ) decay is considered to be the most easily accessible means of extracting information on the possible Majorana mass of the neutrino. The reason for this is that the 0νββ decay occurs only if the neutrino is a massive Majorana particle. The search for the 0νββ decay is mostly concentrated on the 0νβ − β − decays due to their favorable Q values [1]. Assuming that the neutrino carries Majorana mass also the positron emitting modes of 0νββ decay, 0νβ + β + , 0νβ + EC and 0νECEC, are possible but they are hard to detect owing to their small decay Q values [2]. The neutrinoless double electron capture, 0νECEC, can only be realized as a resonant decay [3] or a radiative process with or without a resonance condition [4]. The resonant 0νECEC decay has attracted a lot of experimental attention recently [5][6][7][8][9].
The resonance condition -close degeneracy of the initial and final (excited) atomic states-can enhance the decay rate by a factor as large as 10 6 [3]. Possible candidates for the resonant decays are many [3,4]. Verification of the fulfillment of the resonance condition is of utmost importance before expensive experiments are conducted in the aim to detect the resonant 0νECEC. The half-life of the resonant 0νECEC depends sensitively on the mass difference between the initial and final nuclei of decay. The verification of the resonance enhancement thus calls for accurate measurements of these mass differences. At the JYFL the way to achieve the necessary accuracy -in the 100 eV range-is to engage the Penning ion trap, JYFLTRAP, a e-mail: jouni.suhonen@phys.jyu.fi for these measurements. In fact there is a campaign at the JYFL to run the Penning trap to sort out potential candidates for measurements of resonant 0νECEC decays by dedicated underground experiments.
In this work two examples -the decay of 74 Se to the second 2 + state in 74 Ge and the decay of 136 Ce to the fourth 0 + state in 136 Ba-are discussed from the point of view of the nuclear-structure aspects of resonant 0νECEC decays.
2 Theory framework of the resonant neutrinoless double electron capture

Decay half-life
In this work we study a particular type of neutrinoless double electron capture (0νECEC) decay, namely the resonant 0νECEC process of the form e − +e − +(A, Z) → (A, Z −2) * → (A, Z −2)+γ +2X, (1) where the capture of two atomic electrons leaves the final nucleus in an excited state that decays by one or more gamma-rays and the atomic vacancies are filled by outer electrons with emission of X-rays. The daughter state (A, Z − 2) * is a virtual state with energy including the nuclear excitation energy and the binding energies of the two captured electrons. For the half-life of the parent atom the resonance condition can be written [3,10] in the Breit-Wigner form where M contains the leptonic phase space and the nuclear matrix element, Γ denotes the combined nuclear and atomic radiative widths (few tens of electron volts [11]) and |Q − E| is the so-called degeneracy parameter containing the energy (2) of the virtual final state and the difference between the initial and final atomic masses (the Q value of resonant 0νECEC). The Q value needs to be measured very accurately in order to judge whether the resonant 0νECEC is detectable or not. For this, the JYFLTRAP Penning trap is suited in a perfect manner.
To be able to estimate the half-life we have to evaluate the quantity M in (3). In the mass mode of the resonant 0νECEC, it can be written as where m ν is the effective neutrino mass of the neutrinoless double-beta decay [2] and G ECEC 0ν is the leptonic phase-space factor with G F the Fermi coupling constant, θ C the Cabibbo angle, g A = 1.25 the axial-vector weak coupling constant, Z the charge of the mother nucleus, α the fine-structure constant, m e the electron rest mass and R A the radius of a nucleus with mass number A. The quantity η x in (5) is a suppression factor depending on the atomic orbitals where the two electrons are captured from. Solution of the Schrödinger equation for a point-charge nucleus produces the values This approximation is reasonable and compares well with the Dirac solution for a homogeneously charged spherical nucleus [12]. Expressions for the involved nuclear matrix element M ECEC 0ν are given in sect. 2.2 and 2.3.

Decays to 0 + final states
The nuclear matrix element (NME) involved in (4) for decays to the 0 + final states is defined as where R A = 1.2A 1/3 fm is the nuclear radius, g V (g A ) is the vector (axial-vector) weak coupling constant, M GT is the Gamow-Teller and M F is the Fermi NME. The two matrix elements are given by Here the summation runs over all the states J π of the intermediate nucleus and r mn = |r m − r n | is the relative coordinate between the nucleons m and n. The isospin raising operators τ + n convert the n-th proton to a neutron and σ are the usual Pauli spin matrices. Here 0 + i (0 + f ) is the ground state of the even-even mother (daughter) nucleus and the neutrino potential h K (r mn , E a ), K = F, GT, is defined as where j 0 is the spherical Bessel function. The term h K (q 2 ) in (10) includes the contributions arising from the shortrange correlations, nucleon form factors and higher-order terms of the nucleonic weak current [13][14][15].
Next, we write the nuclear matrix elements explicitly. They are given by where k 1 and k 2 label the different nuclear-model solutions for a given multipole J π , the set k 1 stemming from the calculation based on the final nucleus and the set k 2 stemming from the calculation based on the initial nucleus. The operators O K inside the two-particle matrix element derive from (8) and (9) and they can be written as where E k is the average of the k-th eigenvalues of the nuclear-model calculations based on the initial and final nuclei of the decay. Here the one-body transition densities , and they are given separately for the different nuclear-structure formalisms in sect. 3. The quantity J π k1 |J π k2 is the overlap factor and its purpose is to match the the two sets of J π states that stem from the two different nuclear-model calculations -one starting from the mother nucleus (the k 2 states) and the other starting from the daughter nucleus (the k 1 states). Its explicit form for the pnQRPA nuclear model is given in sect. 3.1.
The two-particle matrix element of (11) can be written as M λ (n 1 lN L; n n l n n n l n )M λ (n 2 lN L; n p l p n p l p ) where j = √ 2j + 1, the neutrino potentials are those of (10) and we have defined The quantities M λ are the Moshinsky brackets that mediate the transformation from the laboratory coordinates r 1 and r 2 to the center-of-mass coordinate R = 1 √ 2 (r 1 + r 2 ) and the relative coordinate r = 1 √ 2 (r 1 − r 2 ). In this way, the short-range correlations of the two decaying protons are easily incorporated in the theory. The wave functions φ nl (r) are taken to be the eigenfunctions of the isotropic harmonic oscillator.

Decays to 2 + final states
The 0νECEC transition to the 2 + f final states is mediated by a spherical thensor of rank two. Denoting this tensor by T 2 one can write the transition matrix element as The most straightforward way to build T 2 is to use a double Gamow-Teller type of operator where the neutrino potential is given in (10). Then the two-nucleon operator O 2 in (15) becomes where E k is again the average of the k-th eigenvalues of the two nuclear-model calculations. The resulting two-particle matrix element is given by where the neutrino potential is the usual one (10) and the other quantities have been defined in sect. 2.2.

Nuclear-structure models
The starting point in the present calculations is the theoretical framework of the quasiparticle random-phase approximation (QRPA). This framework is based on the quasiparticles that are created in a BCS calculation via the Bogoliubov-Valatin transformation [16]. We use two kinds of QRPA approaches. The proton-neutron QRPA (pnQRPA) is used to produce the J π states of the intermediate odd-odd nucleus and the charge-conserving QRPA (ccQRPA) is used to produce the excited states of the even-even daughter nucleus. The states of the pnQRPA and the ccQRPA are connected by decay amplitudes that are calculated in the higher-QRPA framework called the multiple-commutator model (MCM). All these approaches are discussed below.

Proton-neutron quasiparticle random-phase approximation
We use the pnQRPA to form the wave functions that represent states of an odd-odd nucleus. The pnQRPA state is obtained by acting by a pnQRPA creation operator on the QRPA ground state |QRPA . This action can be expressed in the form Here J π is the multipolarity of the nuclear state and and k enumerates the phonons. The operator a † p (a † n ) creates a proton (neutron) quasiparticle in the orbital p (n). The sum runs over all proton-neutron configurations in the chosen valence space. Making use of the ansatz (19) one can derive (see, e.g., [16]) the pnQRPA equations of motion in the quasiboson approximation. The equations can be cast in the matrix form where J is the angular momentum and the matrix elements of the A(J) and B(J) matrices read and with v and u corresponding to the occupation and unoccupation amplitudes stemming from the BCS calculation. Above V is the normalized, J-coupled and antisymmetrized two-body interaction matrix element (see eq. (8.16) of [16]) and the particle-hole matrix element V RES of the residual interaction is obtained from V by the Pandya transformation The quantities E p and E n in (21) are the quasiparticle energies for the proton orbital p and neutron orbital n, respectively. The j p,n in (23) are the total angular momenta of the proton or neutron single-particle orbitals.
In (21) and (22) the particle-hole and particle-particle parts of the pnQRPA matrices are separately scaled by the particle-hole parameter g ph and particle-particle parameter g pp [17][18][19]. The particle-hole parameter affects the position of the Gamow-Teller giant resonance and its value was fixed by the available systematics [16] on the location of the giant state. The parameter g pp affects only weakly the neutrinoless double electron capture rates and it value was set to the default g pp = 1.0.
The pnQRPA transition densities are given by where v (v) and u (ū) correspond to the BCS occupation and unoccupation amplitudes of the initial (final) eveneven nucleus. The amplitudes X and Y (X andȲ ) come from the pnQRPA calculation starting from the initial (final) nucleus of the 0νECEC decay. The overlap factor in (11) and (15) is given by and it takes care of the matching of the corresponding states in the two sets of states based on the initial and final even-even reference nuclei.

The charge-conserving QRPA and the multiple-commutator model
The multiple-commutator model (MCM) [19,20] is designed to connect excited states of an even-even reference nucleus to states of the neighboring odd-odd nucleus. Earlier the MCM has been used extensively in the calculations of double-beta-decay rates [21,22]. The states of the oddodd nucleus are given by the pnQRPA in the form (19). The excited states of the even-even nucleus are generated by the (charge conserving) quasiparticle randomphase approximation (ccQRPA) described in detail in [16].
Here the symmetrized form of the phonon amplitudes is adopted contrary to ref. [16] so that the k-th J π state can be written as a QRPA phonon in the form The symmetrized amplitudes Z and W are obtained from the usual ccQRPA amplitudes X and Y [16] through the following transformation: and similarly for W in terms of Y . The amplitudes X and Y are obtained by solving the QRPA equations of motion that formally look like the matrix equation (20). In the case of the ccQRPA the matrix elements of the A(J) and B(J) matrices read and where N ab (J) is the quasiparticle-pair normalization constant and the amplitudes v and u correspond to the occupation and unoccupation amplitudes of the BCS. Definitions of the two-body interaction matrix elements appearing in (29) and (30) were given in sect. 3.1.
The ccQRPA phonon defines a state in the final nucleus of the resonant 0νECEC decay. In particular, if this final state is the k-th J + state the related transition density, to be inserted in (11) or (15), becomes instead of the expression (24) for the ground-state transition. Again v (v) and u (ū) correspond to the BCS occupation and unoccupation amplitudes of the initial (final) even-even nucleus. The amplitudes X and Y (X andȲ ) come from the pnQRPA calculation starting from the initial (final) nucleus of the 0νECEC decay. The amplitudes Z andW are the amplitudes of the k-th J + state in the 0νECEC daughter nucleus. One can take a J π k = 2 + 1 phonon of (27) and build an ideal two-phonon J + state of the form An ideal two-phonon state consists of partner states J π = 0 + , 2 + , 4 + that are degenerate in energy, and exactly at an energy twice the excitation energy of the 2 + 1 state. In practice, this degeneracy is always lifted by the residual interaction between the one-and two-phonon states [23]. The related transition density of the MCM, that can be inserted in (11) or (15), attains the form where, as usual, the barred quantities denote amplitudes obtained for the 0νECEC daughter nucleus.

Nuclear-structure calculations
In this section we discuss as examples the nuclearstructure calculations related to the 0νECEC decays of two nuclei, 74 Se and 136 Ce, that are located in quite different mass regions. Since 0νECEC decays involve chargechanging decay transitions it is advantageous to perform also an auxiliary analysis by engaging theoretical description of the lateral beta-decay feeding of the low-lying states of the final nuclei of the resonant 0νECEC decays.

Determination of the model parameters
In the case of the decay of 74 Se the nuclear-structure calculation was performed in the 1p-0f-2s-1d-0g (9 orbitals) single-particle valence space for both protons and neutrons. For the decay of 136 Ce the proton single-particle space was chosen as 1p-0f-2s-1d-0g-0h 11/2 (10 orbitals) and that of the neutrons consisted of 2s-1d-0g-2p-1f-0h (11 orbitals). The single-particle energies of all these orbitals were generated by the use of a spherical Coulombcorrected Woods-Saxon (WS) potential with a standard parametrization [24], optimized for nuclei near the line of beta stability. As the two-body interaction we have used the Bonn-A G-matrix and we have renormalized it in the standard way [17][18][19]: The quasiparticles are treated in the BCS formalism and the pairing matrix elements are scaled by a common factor, separately for protons and neutrons. In practice these factors are fitted such that the lowest quasiparticle energies obtained from the BCS match the experimental pairing gaps for protons and neutrons, respectively. The average behavior of both the proton and neutron pairing gaps is given by the three-point formulae (see, e.g., [16]) where S p and S n are the experimental separation energies obtained from the measured nuclear masses (see, e.g., [25]) for protons and neutrons, respectively, and Z and N are the proton and neutron numbers. The intermediate J π states were generated both from the mother and daughter ground states by the use of the pnQRPA, discussed in sect. 3.1. The two obtained sets of J π states were matched by using the overlap (26). The particle-hole strength parameter g ph of the pnQRPA, shown in (21) and (22), determines the energy of the Gamow-Teller giant resonance (GTGR) and its value was fixed by fitting the empirical location of the GTGR [16]. The particle-particle parameter g pp has little effect on the β + /EC type of transitions and we have used the default value g pp = 1.0.
The value of the NME in (4) is affected by both the value of the axial-vector coupling constant g A and the nucleon-nucleon short-range correlations (SRC). In this work the quenched value g A = 1.0 is assumed and we use the unitary correlation operator method (UCOM) to account for the SRC. The UCOM has recently been applied successfully in calculations of the neutrinoless double-beta decay [13][14][15]26]. The effect of the UCOM SRC is to slightly reduce the magnitude of the NME by gently shifting the wave function of the relative motion away from the touching point of the two decaying nucleons. Also the dipole form factors of the nucleons and higher-order nucleonic currents have been included in the present calculations in the way described in [13][14][15].
The final states of the presently discussed decays are the second 2 + state, 2 + 2 , in 74 Ge and the fourth 0 + state, 0 + 4 , in 136 Ba (the 0 + ground state is counted as the first 0 + state). The 2 + 2 state in 74 Ge is assumed to be a member of a two-phonon triplet with a wave function (33). The experimental energy of the involved first 2 + state was reproduced in the ccQRPA calculation by varying the value of the particle-hole strength g ph that scales the particlehole matrix elements in the way presented in eqs. (29) and (30). The value of the particle-particle strength can safely be taken to be g pp = 1.0 since the properties of the ccQRPA states depend only very weakly on the value of this parameter. The 0 + states are a special case in a ccQRPA calculation: The first ccQRPA root is spurious and has to be removed by setting its value to zero. It has to be noted that in the ccQRPA the other 0 + states are not contaminated by this spuriosity [27]. To bring the first ccQRPA root to zero one has to change the values of both g ph and g pp . The values of these parameters can be fixed uniquely by reproducing at the same time the low-energy spectrum of the 0 + states. The low-energy spectra of 74 Ge and 136 Ba are discussed in sect. 4.2.
Finally, it should be stressed that the present calculations assume spherical shape of the involved nuclei. In fact many of the nuclei under study here possess a nonzero static oblate or prolate deformation in their ground state. The corresponding deformation parameter is, however, not very large, of the order of |β| ≤ 0.25. Thus one should expect to encounter effects arising from the nuclear deformation itself or the possible deformation differences between different nuclei participating the weak decay processes. These effects should show up already at the level of single beta decays and that is why it pays to study the lateral beta-decay feeding of the nuclei involved in the double-beta transitions. This aspect, among others, is addressed in sect. 4.2.

Auxiliary analysis through lateral beta decays
The available data on beta-decay rates allow for studies of the lateral beta-decay feeding of the lowenergy states in the 0νECEC daughter nuclei 74 Ge and 136 Ba. The transitions are first-forbidden for 74 As(2 − 1 ) → 74 Ge(0 + f , 2 + f , 4 + f ) and allowed Gamow-Teller for 136 La(1 + 1 ) → 136 Ba(0 + f , 2 + f ). The allowed Gamow-Teller beta-decay transitions of interest in this work are of the type 1 + → 0 + , 2 + . For them the log ft value is defined as [16] log ft = log(f 0 t 1/2 ) = log 6147 B GT , for the initial 1 + and final J + = 0 + , 2 + states. Here is the the leptonic phase-space factor for the allowed β + /EC decays defined in [16]. The first-forbidden decay transitions can be divided into two categories: those with an angular-momentum change of two units (first-forbidden unique, FFU) and those with an angularmomentum change of at most one unit (first-forbidden non-unique, FFNU). For the first-forbidden unique transitions 2 − → 0 + , 4 + we can define [16] log ft = log(f 1u t 1/2 ) = log 6147 where for the initial 2 − and final J + = 0 + , 4 + states.
Here f 1u = f + 1u + f EC 1u is the the leptonic phase-space factor for the FFU β + /EC decays as given in [16]. For the first-forbidden non-unique transitions 2 − → 2 + one defines [16,28,29] log ft = log(f 0 t 1/2 ) = log f 0 6147 S 1 , where S 1 is the integrated shape factor that is a complex combination of the various nuclear matrix elements and phase-space factors, see ref. [29].
Results of the beta-decay calculations for 74 Ge are shown in table 1. There the experimental and computed energies of the 2 + 1 state and the two-phonon triplet in In particular, the theoretical decay rate to the two-phonon 2 + state, 2 + 2 , is too high pointing to slight overestimation of the magnitude of the corresponding 0νECEC NME.
The 0 + 4 state in 136 Ba is presumed to be a ccQRPA phonon of the form (27). The experimental and computed energies of this state and a number of other low-energy 0 + and 2 + states in 136 Ba are shown in table 1 in columns two and four. As can be seen the computed energies agree very well with the experimental ones. In the same table we show also the log ft values of the β + /EC-decay feeding of these states via allowed Gamow-Teller transitions from the 1 + ground state of 136 La. As can be seen the computed log ft values are usually a bit too small indicating that the asso-ciated transitions are predicted to be too fast by the MCM formalism. In particular, the theoretical decay rate to the 0 + 4 state is too high pointing to a possible overestimation of the magnitude of the corresponding 0νECEC NME.
One possible source of the noticed differences between the computed and measured log ft values could be the omission of the deformation degree of freedom in the present calculations. In particular the deformation differences between the mother and daughter nuclei could drive the suppression of beta-decay transition rates. Such tendencies are visible in tables 1 and 2, especially for the decay of 136 La in table 2. However, it seems that the deformation effects are not so strong as to ruin the overall compatibility of the trends visible in the computed and measured log ft values.

Results for the NMEs of the resonant 0νECEC decays
We can use relations (3) and (4) to write the half-life of the resonant 0νECEC in the form where the effective neutrino mass has to be inserted in units of eV. The 0νECEC-decay half-life of 74 Se can now be obtained directly by the use of eqs. (41), (4), (15), (25) and (34), where the second leg of the decay is evaluated by the use of the MCM-calculated transition density (34). The corresponding value for the matrix element in (15) is M ECEC 0ν (0 + → 2 + 2 ) = 0.624 MeV. However, additional considerations, discussed extensively in [32], reduce this value by some two orders of magnitude. Such a small value of M ECEC 0ν results in an extremely long 0νECEC half-life [32] irrespective of the fulfillment of the resonance condition (3). Also the uncertainty in the value of the final NME should be accounted for. In this case, the uncertainty is considerable due to the recoil structure of the NME. The decay rate is further suppressed by the JYFLTRAP-measured Q value that yields a rather large value |Q − E| min = 2.4 keV for the degeneracy parameter. An additional suppression is imposed by the necessary capture from the L 2 (2p 1/2 ) and L 3 (2p 3/2 ) atomic orbitals (see (6)). Summing up all these effects produces a half-life value (41) with C ECEC ≈ 2 × 10 42 − 1 × 10 45 being indicative of the large uncertainty in the recoil structure of the NME, not so much of the uncertainties in the nuclear wave functions.

Conclusions
A general theoretical framework of the resonant neutrinoless double electron capture is presented from the point of view of the involved nuclear matrix elements. Use of this theory framework together with the Q value measurements by the JYFLTRAP produces information about the nuclei whose 0νECEC decays could be detected in dedicated underground experiments. This, in turn, is important in pinning down the value of the possible Majorana mass of the neutrino.
In this work explicit expressions for the amplitudes of the resonant 0νECEC decays have been derived. The nuclear matrix elements have been evaluated within the higher-QRPA framework of the multiple commutator model. The associated transition densities have been derived and presented in the article. An auxiliary analysis has been performed by the use of the lateral betadecay feeding of the final nuclei of the double electron capture decays. The related allowed and first-forbidden transitions have been defined and computed by the use of the multiple-commutator model. As an example, the developed formalism is applied to the decays of 74 Se and 136 Ce. The values of the computed matrix elements and the Q values measured by the JYFLTRAP suggest that the resonant double electron captures in 74 Se and 136 Ce are impossible to detect in foreseeable future.