Bound-state double-$\beta$ decay

We consider new modes of two-neutrino and neutrinoless double-$\beta$ decays in which one $\beta$ electron goes over to a continuous spectrum and the other occupies a vacant bound level of the daughter ion. We calculate the corresponding phase-space factors of the final states, estimate the partial decay rates, and derive the one- and two-electron energy spectra using relativistic many-electron wave functions of atoms provided by the multiconfiguration Dirac-Hartree-Fock package GRASP2K. While the bound-state neutrinoless double-$\beta$ decays are strongly suppressed, their two-neutrino counterparts can be observed in the next-generation double-$\beta$-decay experiments, most notably SuperNEMO.


I. INTRODUCTION
Among the most challenging problems of the modern neutrino physics are the mechanism of neutrino mixing and the nature of neutrino masses (Dirac or Majorana). If the diagonal neutrinos ν i (i = 1, 2, 3) are Majorana fermions, then the flavor neutrinos ν α (α = e, µ, τ ) are identical to their charge-conjugated states, as a result of which the total lepton number is not conserved (see, e.g., [1]). Observation of the neutrinoless double-β (0νββ) decay can provide evidence for the Majorana nature of massive neutrinos, which would be of great value for extensions of the Standard Model [2]. Measurement of the half-life of a 0νββ decay could provide a key to the absolute scale of neutrino masses and also shed light on the leptonic CP violation mechanism required to explain the observed baryon asymmetry of the Universe [2,3]. Given the opportunity to get answers to so many fundamental questions, the 0νββ decay has attracted much attention of theorists and experimentalists in the recent decades.
The neutrinoless (two-neutrino) double-β decay of a parent nucleus A Z X into a daughter nucleus A Z+2 Y, denoted 0ν(2ν)ββ, involves an emission of two electrons e − (and a pair of electron antineutrinos ν e ) from the atom: The 2νββ decay occurs in the 2nd order of weak interaction and as such it conserves the total lepton number: ∆L = 0. It forms the dominant decay channel of beta radioactivity of the even-even isotopes for which the single-β decay into the odd-odd intermediate nucleus is either energetically forbidden or suppressed by spin selection rules. The double-β decay has so far been observed in 11 out of 35 candidate isotopes, with half-lives T 2νββ 1/2 ∼ 10 19 -10 21 y, making it the rarest known spon-taneous decay in nuclear physics. In contrast, the 0νββ decay violates the total lepton number by two units: ∆L = +2, and requires a Majorana mass term. Such process could be observed as a monoenergetic peak at the 2νββ spectrum endpoint in calorimetric measurements of the sum of electron energies. The current limits on the half-lives set T 0νββ 1/2 > (0.18-1.07) × 10 26 y at 90% C.L. for the 136 Xe and 76 Ge isotopes [4][5][6].
In 1961, Bahcall [7] developed a formalism for description of bound-state β decays in which the β-electron is produced in an atomic K or L shell, while the monochromatic antineutrino ν e carries away the entire energy of the decay. The bound-state β decay was observed on bare 163 66 Dy 66+ ions collected in the heavy-ion storage ring ESR at GSI, Darmstadt, with a half-life of 47 d for the otherwise stable nuclide [8].
The neutrinoless double-β decay with two bound electrons in the final state, denoted 0νEPEP (where EP stands for the electron production): was discussed in Ref. [9] as an inverse process to the neutrinoless double-electron capture. A resonant enhancement of the 0νEPEP decay probability can occur in the case of quasi-degeneracy of the initial-and final-state atomic energies. The ground-state 0 + −→ 0 + nuclear transition of 148 Nd to an 1.921 MeV excited state of 148 Sm * fulfills the resonance condition with the experimental accuracy of ∼ 10 keV. The estimated half-life, however, was found to be beyond the reach of experiments at the present stage.
In this paper, we develop a formalism for description of the bound-state two-neutrino and neutrinoless double-β  Figure 1. A schematic view of the 0ν(2ν)EPβ decays. The final state involves a daughter nucleus A Z+2 Y, a bound electron e − b produced above the subshells occupied by Z atomic electrons, and a single free electron e − (and a pair of electron antineutrinos νe) emitted from the atom. Upon the deexcitation, the bound electron e − b radiates photons of energy 10 eV.
decays, denoted 0ν(2ν)EPβ: The process is shown schematically in Fig. 1. The appearance of the first β-electron in the continuous energy spectrum is accompanied by a production of the second β-electron in a vacant discrete ns 1/2 or np 1/2 level above the valence shell of the daughter ion A Z+2 Y 2+ . The inclusion of atomic levels with higher angular momenta is not required because their wave functions exhibit only a negligible overlap with the nucleus. Since the 0νEPβ, 0νββ, 2νEPβ and 2νββ decay modes constitute 1-, 2-, 3-and 4-body decays, respectively, they could be distinguished by their one-and two-electron energy distributions.
The outline of the paper is as follows. In Sec. II, the relativistic electron wave functions as one-particle solutions to the Dirac equation are described and expressions for the relativistic Fermi function and its bound-state analog are derived. In Sec. III, the 0ν(2ν)EPβ decay rates are derived within the V − A weak interaction theory including the mixing of Majorana neutrinos. We restrict ourselves to the ground-state 0 + −→ 0 + nuclear transitions and obtain the phase-space factors entering the decay rates. Section IV describes the evaluation of relativistic bound-electron wave functions at short distances via the multiconfiguration Dirac-Hartree-Fock package Grasp2K [10]. Numerical estimates of the half-lives and the 0ν(2ν)EPβ to 0ν(2ν)ββ decay-rate ratios are given in Sec. V, in addition to the one-and two-electron energy spectra. In Sec. VI, we finally draw conclusions regarding a possible experimental observation of the boundstate double-β decays and provide a motivation for further studies.

II. RELATIVISTIC ELECTRON WAVE FUNCTIONS IN CENTRAL FIELD
The electronic structure of atoms is described by the shell-model relativistic wave functions obtained as solutions to the Dirac equation in a self-consistent centrally symmetric potential, which represents a superposition of the nuclear Coulomb potential and the screening potential of the electron shell. The corresponding bispinors with separated radial (r = |r|) and angular (n = r/|r|) variables take the form (see, e.g., [11]): where κ = (l − j)(2j + 1) = ±1, ±2, . . . labels combinations of the orbital l = 0, 1, . . . and spin s = 1/2 angular momenta (κ = −1, +1 for ns 1/2 and np 1/2 states, respectively), while µ = −j, . . . , +j denotes the projection of the total angular momentum j = l + s onto the z-axis.
The radial functions f κ (r) and g κ (r) in the continuum further depend on the electron energy E. In the double-β decays, the leading s 1/2 term of the partial-wave expansion which enters the nuclear matrix elements reads [12]: wherep is a unit vector in the direction of the electron momentum p. Note that the continuum radial functions are normalized to the δ function in p = |p| [11], in contrast to the bound states which obey [13]: dr r 2 (f 2 + g 2 ) = 1. The Fermi function F (Z, E), introduced to correct the short-distance behavior of the β-electron plane waves due to the Coulomb potential, is defined in terms of the radial wave functions f −1 (E, r) and g +1 (E, r) evaluated at the nuclear surface at r = R 1.2 fm A 1/3 : where γ = κ 2 − (αZ) 2 , ν = αZE/p, p = E 2 − m 2 e , m e is the electron mass, and α 1/137 is the finestructure constant. We remark that F (Z, E) → 1 for Z → 0. For αZ 1 and l = 0, the Fermi function F (Z, E) coincides with the Gamow-Sommerfeld factor [14][15][16].
The Fermi function in Eq. (7) is given by a standard approximation [12] in which the relativistic electron wave function for a uniform charge distribution in the nucleus is considered and only the lowest-order terms in the power expansion in r are taken into account. The exact Dirac electron wave function accounting for a finite nuclear size and electron-shell screening effects [17] modifies the 0νββ-decay phase-space factor for 150 Nd by 30% (see Ref. [18] and Table 1 therein), which results in an increase of the 0νββ-decay half-life. The 0ν(2ν)EPβ decay rate with one electron in the continuous spectrum is thus less sensitive to the details of the Dirac electron wave function, since only one Fermi function enters the corresponding phase-space factor. We therefore restrict ourselves to the continuous-spectrum solutions of the Coulomb problem for V (r) = −α(Z + 2)/r, where Z + 2 is the atomic number of the daughter nucleus A Z+2 Y. In the discrete spectrum, the radial wave functions f nκ (r) and g nκ (r) in the Coulomb potential correspond to the energy eigenvalues (see, e.g., [11]): where n = 1, 2, . . . is the principal quantum number and n r = n − |κ| is the radial quantum number which counts the number of radial nodes. At small distances nκ , the leading term from the series expansion of the radial wave functions f nκ (r) and g nκ (r) for a point-like source can be found in Ref. [9].
The radial wave functions enter the bound-state βdecay probabilities in the combination: which formally coincides with the relativistic Fermi function (7). Note, however, that the first and second terms in the right-hand side of Eq. (9) originate from the production of β-electrons in the ns 1/2 and np 1/2 orbits, respectively. For αZ 1 and l = 0, we have f nκ (r) ≈ R nl (r) and g nκ (r) ≈ 0, where R nl (r) is the nonrelativistic radial wave function obtained by solving the Schrödinger wave equation for a hydrogen-like atom. The screening of the Coulomb potential modifies the short-distance behavior of the bound-state wave functions. This effect is taken into account in Sec. IV via the relativistic atomic structure package Grasp2K.

III. PHASE-SPACE FACTORS
The double-β decay is a 2nd-order process governed by the effective beta-decay Hamiltonian: Here, G β = G F cos θ C includes the Fermi coupling constant G F 1.166 × 10 −5 GeV −2 together with the Cabibbo angle θ C 13 • due to the quark mixing [19], e and ν e are the electron and electron-neutrino fields, respectively, and the baryon charged current j µ = p γ µ g V − g A γ 5 n couples the proton and neutron fields via the vector g V = 1 and (unquenched) axial-vector g A 1.27 coupling constants. The V−A structure of H β ensures that only the left-handed leptons enter the weak interaction. The flavor-and diagonal-neutrino fields are related by the unitary 3×3 Pontecorvo-Maki-Nakagawa-Sakata (PMNS) matrix U : The neutrinoless double-β decay is assumed to be related to a light Majorana-neutrino exchange between nucleons in the parent nucleus. The inverse 0νββ and 2νββ half-lives (see, e.g., [17]): factorize in terms of the kinematical phase-space factors G 0ν(2ν)ββ (Z, Q), the nuclear matrix elements (NMEs) M 0ν(2ν)ββ , and the effective Majorana neutrino mass: where m i are the masses of diagonal neutrinos. Since the absolute scale of neutrino masses and the Majorana phases are unknown, the value of |m ββ | is treated as a parameter. The experimental lower bounds on T 0νββ 1/2 set an upper limit on |m ββ |. The most stringent limit has so far been obtained in the KamLAND-Zen experiment [6]: |m ββ | < 61-165 meV at 90% C.L., where the range of values accounts for the uncertainties inherent to the nuclear-structure models. In case of the inverted hierarchy of neutrino masses, the effective mass is constrained by cosmology: |m ββ | = 20-50 meV. We estimate the 0νEPβ and 0νββ half-lives assuming |m ββ | = 50 meV. Since the 2νββ half-life is unambiguously defined within the Standard Model, the measured values of T 2νββ 1/2 can be used to fix the phenomenological parameters, improve the predictions of the nuclear-structure models for M 0νββ and probe the possible quenching of g A .
The energy conservation in the 0ν(2ν)ββ decays implies: where M i and M f are the masses of the parent and daughter nuclei and E 1,2 and ω 1,2 are the total energies of the emitted electrons and antineutrinos, respectively. The total released kinetic energy in both scenarios equals: Due to indistinguishability of the final-state leptons, the NMEs contain a superposition of two (four) energy denominators [20]: where E n denotes the nth energy level of the intermediate nucleus and q = (q 0 , q) is the four-momentum of the exchanged Majorana neutrino. Since q 0 = q 2 + m 2 i ≈ |q| ∼ 200 MeV, the difference between the lepton energies can be safely neglected: In case of the 0ν(2ν)EPβ decay modes, a similar approximation ensures that the corresponding NMEs remain essentially unchanged: M 0ν(2ν)EPβ ≈ M 0ν(2ν)ββ and the distinction between the 0ν(2ν)EPβ and 0ν(2ν)ββ decay modes is fully reflected in the phase-space factors G 0ν(2ν)EPβ (Z, Q).
The phase-space factors of the 0ν(2ν)EPβ decays can be found to be: where n min is the principal quantum number of the lowest available electron shell (this can in principle be different for the s 1/2 and p 1/2 states). Equations (15)-(16) can be derived from G 0νββ and G 2νββ using the substitution: and taking into account the identity of the electrons: the integrated phase space of the 0ν(2ν)ββ decays contains a statistical factor of 1/2!, which is not present in case of the 0ν(2ν)EPβ decay modes since the bound and free electrons occupy complementary regions of the phase space. The corresponding rule for the integrated phase space reads: In the bound-state double-β decays, the binding energy of the produced electron 10 eV can be neglected. Such an approximation does not affect the required accuracy, but greatly simplifies the computation since the infinite sum of integrals in Eq. (18) is factorized into the Fermi sum ∞ n=nmin B n (Z + 2) and just one integral independent of n. The energy conservation in the 0νEPβ decay implies that the free electron carries away the entire energy released in the decay: E = m e + Q, whereas in the 2νEPβ decay the energy is distributed between the electron and the two antineutrinos:

IV. BOUND-STATE WAVE FUNCTIONS OF ELECTRONS IN DIRAC-HARTREE-FOCK METHOD
The multiconfiguration Dirac-Hartree-Fock package Grasp2K solves a stationary N -body Dirac equation with the separable central atomic Hamiltonian [10]: where α = γ 0 γ and β = γ 0 . The first two kineticenergy terms are followed by the potential-energy terms which account for electron-nucleus Coulomb attraction and electron-electron Coulomb repulsion, respectively, where the latter is approximated by a mean field V (r i ) generated by the surrounding electron cloud. The separability ensures that the energy eigenvalues are additive: while the many-electron wave functions are expressed in terms of the Slater determinants: where ψ i (r j ) = ψ niκiµi (r j ) and P is a permutation of the quantum numbers with parity (−1) P . The nuclear part of the total wave function is disregarded by virtue of the Born-Oppenheimer approximation. The self-consistent field procedure then varies the radial orbitals f nκ (r) and g nκ (r) in iterative cycles until a convergence is achieved.
The radial orbitals f n,−1 (R) and g n,+1 (R) are computed in the nuclear Coulomb potential of the daughter nucleus A Z+2 Y for the ground-state electron configuration of the parent atom A Z X with the additional β-decay electron occupying an empty orbit. Since the convergence cannot be always guaranteed and the program only provides the electron-shell wave functions up to n = 9, we employ a combined approach: 1. The radial orbitals f n,−1 (R) and g n,+1 (R) are calculated based on initial estimates provided by the Thomas-Fermi model.
2. If the convergence cannot be achieved within the specified number of iterations, the radial orbitals f n,−1 (R) and g n,+1 (R) are calculated based on initial estimates provided by the non-relativistic Hartree-Fock approximation.
3. If both methods fail for the charge Z, we are looking for the values of Z = Z for which the calculation can be completed. The squares of the radial orbitals are then determined by fitting the available values for a fixed orbit using a power-law function: Finally, the squares of the radial orbitals with the principal quantum numbers above n = 9 are estimated for a given isotope from a fit of the available values for n ≤ 9 using a power-law function:  In the atomic spectroscopy, power functions are often used to fit the dependence of observables on the atomic number Z (see, e.g., [21]). On the other hand, the power law of the principal quantum number n is motivated by the fact that, in the absence of shielding, the squares of nonrelativistic radial orbitals ns 1/2 decrease at the origin as: R 2 n0 (0) ∝ n −3 . The simple power law enables us to explicitly perform the summation in ∞ n=nmin B n (Z + 2) over the vacancies in the electron shell. The sum is expressed in terms of the Riemann zeta function ζ(z) = ∞ n=1 1/n z . In average, the radial orbitals with n > 9 contribute to the decay rates at the level of only ∼ 4% of the total value. Figure 2 shows the results for a power-law fitting of the squared radial orbitals f 2 n,−1 (R) and g 2 n,+1 (R) at the nuclear radius r = R as functions of the initial nuclear charge Z. The example of 8s 1/2 and 8p 1/2 subshells is considered. In general, the convergence cannot be achieved for all nuclei. The power-law dependence is in excellent agreement with the observed behavior of the computed radial wave functions. The radial orbitals at the nuclear radius r = R as functions of the principal quantum number n are shown in Fig. 3 for the isotope A simple qualitative explanation for the dependence of the bound-electron radial wave functions on Z and n at r = R follows from the following considerations. The nodes of the radial part of a nonrelativistic wave function with l = 0 are localized partially outside the atom at r 1 (in a.u.) and partially inside the atom at r 1. The number of nodes inside the atom can be estimated for highly excited states using a semiclassical approximation, which is justified for Z 1 and r 1. At the boundary of the atom, the phase of the radial wave function is estimated to be:  the Coulomb potential, the radial wave function for small r behaves like ∼ 1/n 3 . The atomic radius ∼ 1 is small compared to the average radius ∼ n 2 of the bound βelectron. The ratio R n0 (1)/R n0 (0) is independent of n for large n and tends to 0.283 at the infinity. Since n a nodes moved inside the atom, the square of wave function at the atomic boundary becomes: R 2 n0 (1) ∼ 1/(n − n a ) 3 . The matching at r ∼ 1 of the outer part of the wave function with the semiclassical wave function at r 1 leads to the appearance at r ∼ 0 of an additional factor Z (see, e.g., [22]), so finally: The same result follows from the requirement of orthogonality of the wave function of the bound β-electron to the electron wave functions in the atom. The number of electrons occupying the atomic levels up to the principal quantum number n s with completely filled outer shell is expressed as follows: In agreement with the semiclassical arguments given above, n s ∼ (3Z/2) 1/3 . To ensure orthogonality, the bound β-electron should have one more node inside the atom compared to n s −1. One can verify that for n a ∼ n s Eq. (21) reproduces the qualitative behavior of the upper radial function for r = R. The dependence on Z for n = 6, shown in Fig. 2, appears reasonable for Z 20. In the case of 82 34 Se, shown in Fig. 3, the approximation (21) works reasonably well for n 7. We remark that Eq.  atoms, so the results are applicable directly to gaseous substances such as krypton or xenon. We expect that the presented calculations yield reasonable estimates also for solids.

V. RESULTS AND DISCUSSION
In Table I, the double-β-decaying isotopes A Z X are listed together with (a) the Q values obtained from the recent evaluation of atomic masses [23], (b) the Fermi sums ∞ n=nmin B n (Z + 2) computed using the Grasp2K package, (c) the phase-space factors G 0ν(2ν)EPβ and G 0ν(2ν)ββ associated with the ground-state 0 + −→ 0 + nuclear transitions, and (d) the decay-rate ratios: which are independent of the NMEs and m ββ , and hence are free of uncertainties inherent to the nuclear-structure models and the neutrino masses. The Fermi sum ∞ n=nmin B n (Z + 2), shown in Fig. 4, increases with the initial atomic number Z and drops whenever the valence shell becomes fully occupied; the very large value of 2.199 × 10 3 for the isotope 198 78 Pt with n min = 6 is out of bounds. The decay-rate ratios Γ 0ν(2ν)EPβ /Γ 0ν(2ν)ββ , shown in Figs. 5-6, achieve their maximum for the isotopes with very low Q values: 98 Mo, 80 Se and 146 Nd, and decrease with increasing both Z and Q. The two-neutrino channels exhibit decay-rate ratios by one order of magnitude higher than the neutrinoless channels. The overall suppression is mainly attributed to the presence of other electrons in the atom: the lowlying electron states (which would otherwise provide a dominant contribution) are already occupied, while the shielding effect of nuclear charge substantially reduces the bound-electron wave functions on the surface of the nucleus.
In Table II, the double-β-decaying isotopes A Z X with available NMEs are listed together with their half-lives T   Figure 6. The decay-rate ratio Γ 2νEPβ /Γ 2νββ as a function of the atomic number Z of the parent nucleus and the Q value. The two-neutrino mode exhibits a behavior similar to the neutrinoless mode (see Fig. 5), but the absolute values are by one order of magnitude higher.
approach including the CD-Bonn nucleon-nucleon potential with short-range correlations and the partial isospinsymmetry restoration [24], except for the isotope 150 Nd which was treated separately within the deformed pn-QRPA model [25]. We estimate the neutrinoless half-lives assuming the unquenched value of the axial-vector coupling constant g A = 1.27 and the effective Majorana neutrino mass at the top of the allowed inverted-hierarchy region: |m ββ | = 50 meV. The half-lives T 2νEPβ 1/2 are derived  [23], the Fermi sums ∞ n=n min Bn(Z +2) over the vacant electron shells of the daughter ion (in a.u.), the bound-state decay phase-space factors G 0ν(2ν)EPβ from Eqs. (15)- (16), the standard phase-space factors G 0ν(2ν)ββ , and the relative frequencies of bound-state to continuum-state decays Γ 0ν(2ν)EPβ /Γ 0ν(2ν)ββ .  [26]; these are further used to extract the NMEs listed for g A = 1.27. While the 0νEPβ decay mode is strongly suppressed and can hardly be experimentally observed in a near future, the half-lives of its 2νEPβ counterpart are already comparable to the present sensitivity to the 0νββ decay. Figures 7 and 8 show the neutrinoless and twoneutrino double-β-decay half-lives for the isotopes listed in Table II. The 0ν(2ν)EPβ and 0ν(2ν)ββ one-electron spectra are described by the differential decay rates 1/Γ dΓ/dε, conventionally normalized to unity and expressed as functions of the dimensionless portion of the electron kinetic energy ε = (E − m e )/Q:  Shown in Figs. 9-10 are the one-electron spectra for the neutrinoless and two-neutrino double-β decays of 82 34 Se. The 0νEPβ peak consists of a large number of discrete contributions, each shifted above the Q value by the electron binding energy ( 10 eV); however, these are indistinguishable under any realistic energy resolution. The 2νEPβ spectrum covers the entire energy range, which could lead to a slight deformation of the measured 2νββ data.
The one-electron spectra are studied with unprecedented accuracy in the tracking-and-calorimetry doubleβ decay experiments based on the external-source technique at the Modane Underground Laboratory (LSM). The NEMO-3 detector [27], which operated during 2003-2011, exploited a cylindrical geometry and observed more than 7 × 10 5 positive 2νββ events with a high signal-tobackground ratio for 7 kg of its primary source isotope 100 Mo during 3.5 y of data taking (the low-radon phase) [28]. The next-generation detector SuperNEMO [29] [26] for the isotopes with double-β decays observed experimentally. The 2νEPβ and 0νββ decay rates are comparable in magnitude.
which is currently under construction, will deploy the source modules comprising 20 thin foils totalling in 100 kg of enriched and purified 82 Se, with possible addition of 48 Ca or 150 Nd isotopes. The tracking chamber will consist of nine planar high-granularity drift cells operating in Geiger regime in a magnetic field of 2.5 mT, and thus enable charge-sign particle identification and vertex reconstruction, secure enhanced background rejection, and provide means to study angular correlations in addition to the one-electron spectra. Calorimeter walls will be composed of segmented low-Z organic-scintillator blocks connected to photomultiplier tubes, striving to achieve the energy resolution: FWHM/Q = 7%/ Q/MeV in the region of interest (ROI) 2.8-3.2 MeV around the endpoint Q 2.998 MeV. The first planar SuperNEMO module "Demonstrator" with 7 kg of the source isotope 82 Se is currently in its final stages of the development.
For the data analysis, it is often desirable to specify the ratio: between the integrated one-electron (two-electron) spectra as a function of the interval [ε (12)min , ε (12)max ] in order to identify the ROIs in which the 2νEPβ decay is best visible relative to its 2νββ counterpart. While the oneelectron ratios are maximal in a small ROI at the spectrum endpoint Q, the two-electron ratios reveal the highest 2νEPβ sensitivity near the opposite end of the energy domain. In these ROIs, the 2νEPβ decay mode could for the given isotopes account for as much as ∼ 100 ppm of the registered events. The ratios from Eq. (27) for the one-and two-electron spectra associated with the decays of 82 Se and 76 Ge, respectively, are shown in Fig. 12.
At temperatures T α 2 Z 2 m e ∼ 10 8 (Z/34) 2 K, atoms become fully ionized and the β-electrons can occupy all discrete levels, provided that the Debye screening length λ D is sufficiently large. In such case, the Fermi sum ∞ n=1 B n (Z + 2) is enhanced by 3-5 orders of magnitude and some of the decay-rate ratios Γ 0ν(2ν)EPβ /Γ 0ν(2ν)ββ exceed unity. The effect can be interpreted as follows: the sum ∞ n=nmin R 2 n0 (0) ∼ Z/(n min −n a ) 2 from Eq. (21) is replaced due to the full ionization by its hydrogenlike analog ∞ n=1 R 2 n0 (0) ∼ Z 3 . For the parent isotope 82 34 Se with n min = 5 (for the ns 1/2 states) and n a = (3Z/2) 1/3 , the enhancement factor can be estimated to give ∼ 2 × 10 3 , and it increases with Z. The 0νEPβ decay channel becomes the only possible one for The one-electron ratio is given for 82 Se (left panel) and the two-electron ratio refers to 76 Ge (right panel). The ROIs with the highest 2νEPβ sensitivity belong to the opposite sides of the ε (12) -interval at ε = 1 and ε12 = 0, respectively. the fully ionized atoms of 98 Mo and 146 Nd, in addition to 80 Se, 114 Cd, 122 Sn, 134 Xe and the rest of doubleβ-decaying isotopes starting from 170 Er in case of the 2νEPβ decay.
In plasma conditions, there is a shift and broadening of the atomic levels which affects the bound-state decay rates [31]. In an extreme case when the Debye screening length λ D decreases below the Bohr radius a 0 , the discrete levels of atoms are pushed to the continuum and, as a result, the bound states cease to exist. This phenomenon is known as the Mott transition [32]. In the cores of the Sun and Sun-like stars where λ D a 0 , the discrete levels of hydrogen are nonexistent. A similar situation occurs in the inner layers of white dwarfs. In the radiative zone of the Sun, e.g., where λ D = (0.7-4) a 0 , the lowest discrete levels of hydrogen become a discrete part of the spectrum but remain vacant because of the ionization. The bound-state double-β decays can thus occur in the outer layers of stars where the screening length is sufficiently large.

VI. CONCLUSION
In this paper, we studied the bound-state two-neutrino and neutrinoless double-β decays. The corresponding phase-space factors were calculated in the framework of the V − A weak-interaction theory including the mixing of Majorana neutrinos. The continuum wave functions of the β-electrons were approximated by the solutions to the Dirac equation in the Coulomb potential of the daughter nucleus, while the relativistic bound-electron wave functions, which are sensitive to the electron-shell screening effects, were computed via the multiconfiguration Dirac-Hartree-Fock package Grasp2K. The ratios between the decay rates of the bound-state and continuumstate double-β decays, which are independent of the nuclear matrix elements and the effective Majorana neutrino mass, are maximal for the isotopes with very low Q values.
The bound-state double-β decays were found to be several orders of magnitude less probable than the continuum-state double-β decays. The bound-state neutrinoless channel is therefore not very suitable in the searches for lepton number violation. In contrast, the sensitivity of the modern 0νββ-decay experiments is already sufficient to observe the 2νEPβ decay mode. We propose to set experimental limits on the 0νEPβ peak and study the 2νEPβ one-electron spectra in the tracking-and-calorimetry double-β-decay experiment NEMO-3 and its next-generation successor SuperNEMO, and examine the two-electron spectra in the calorimetric experiments CUORE [30], EXO-200 [4] and GERDA [5], as well as their upcoming tonne-scale upgrades.
Since under the standard conditions for pressure and temperature most of the double-β-decaying isotopes are solids, it would be desirable to generalize the proposed formalism to the scenario in which the electron shells belong to atoms embedded in a crystal lattice.