Comparative Analysis of Nuclear Matrix Elements of 0νβ+β+ Decay and Muon Capture in 106Cd

Comparative analyses of the nuclear matrix elements (NMEs) related to the 0νβ+β+ decay of 106Cd to the ground state of 106Pd and the ordinary muon capture (OMC) in 106Cd are performed. This is the first time the OMC NMEs are studied for a nucleus decaying via positron-emitting/electron-capture modes of double beta decay. All the present calculations are based on the proton-neutron quasiparticle random-phase approximation with large no-core single-particle bases and realistic two-nucleon interactions. The effect of the particle-particle interaction parameter gpp of pnQRPA on the NMEs is discussed. In the case of the OMC, the effect of different bound-muon wave functions is studied.


INTRODUCTION
Neutrinoless double beta (0νββ) decay is a process in which a nucleus (A, Z), with mass number A and proton number Z, decays to a daughter nucleus with two more or two less protons. In the 0νβ − β − mode the final nucleus is (A, Z + 2), and two electrons are emitted. In the case of the 0νβ + β + mode the final nucleus is (A, Z − 2), and two positrons are emitted. In the latter case also the electron capture (EC) is possible through the mode 0νβ + EC. In this article we denote also this mode by the generic symbol 0νβ + β + . In the case of 106 Cd also the 0νECEC (neutrinoless double electron capture) [1] is possible, but it goes to an excited state, and thus is beyond the scope of the present work, as we analyze here only the ground-state-to-ground-state transition. In addition, it should be noted that the same nuclear matrix elements (NMEs) are involved in the 0νβ + β + and 0νβ + EC modes. In any case, the neutrinoless double beta decay would immediately provide striking new-physics vistas beyond the standard model, since it not only violates lepton-number conservation, but also requires the neutrino to be of Majorana character. After the discovery of neutrino oscillations [2][3][4], this process has become even the more of vital interest because its discovery could potentially provide us with information on the yet-unknown absolute mass scale of neutrinos.
While neutrinoless double beta decay remains hypothetical, the two-neutrino decay mode (2νββ), which is allowed by the standard model, has been observed in several isotopes. Most of the observed decays are of two-neutrino double-beta minus (2νβ − β − ) type, and there are only six isotopes capable of 2νβ + β + decaying: 78 Kr, 96 Ru, 106 Cd, 124 Xe, 130 Ba, and 136 Ce. Naturally, these isotopes can also decay via the 2νβ + EC and 2νECEC (two-neutrino double electron capture) modes. Of these isotopes, 106 Cd is a particularly promising candidate for the 2νβ + β + -decay searches since it has the biggest decay energy, Q ββ = 2775.39 (10) keV, as well as other experimentally favorable features. At present, there are three running experiments searching for the β + β + decay of 106 Cd, namely COBRA [5,6], TGV-2 [7], and 106 CdWO 4 crystal scintillator [8].
Ordinary muon capture (OMC) on nuclei is a weakinteraction nuclear process, in which a negative muon µ − is captured by a nucleus (A, Z) resulting in an atomic-number reduction by one and emission of a muon neutrino. It can significantly extend the kinematic region of ordinary beta decay, owing to the high energy release and large momentum transfer associated with the process. The energy release in the nuclear capture process is about 100 MeV, of which the largest fraction is donated to the neutrino, being the lightest particle in the process. Large mass of the captured muon allows highly-forbidden transitions and high excitation energies of the final states. These features make the OMC a particularly promising probe for the 0νββ decay. In fact, there are several completed, ongoing and planned experiments aiming to study OMC in double-beta-decay triplets. In [9], partial OMC rates to numerous excited states of intermediate nuclei of ββ-decay triplets, including the A = 106 triplet we are studying here, were extracted from γ-ray spectra. In [10], on the other hand, OMC strength function and the associated giant resonances in 100 Nb were studied for the first time. There is an ongoing joint program pursued at RCNP, J-PARC, and the Paul Scherrer Institute (PSI) aiming to extend these studies to a wide range of nuclei from sd-shell nuclei such as 24 Mg up to as heavy nuclei as 240 Pu [11].
In the work of Kortelainen et al. [12] the OMC rates were compared against the 2νβ − β − -decay NMEs for light nuclei using the nuclear shell model. It was found that there was a clear correlation between the energy-distributed OMC rates to 1 + states and the energy-based decomposition of the NMEs for the 2νβ − β − decays of the sd-shell nuclei 36 Ar, 46 Ca, and 48 Ca. In [13], we extended these studies to 0νβ − β − decays of medium-heavy and heavy nuclei by computing the average matrix elements corresponding to the OMC transitions to the intermediate nuclei of 0νβ − β − decays up to some 50 MeV using the pnQRPA formalism. We then compared these matrix elements with the energy-multipole decompositions of the NMEs of 0νβ − β − decays computed using the same formalism and model spaces. We found that there are clear correspondencies between the 0νβ − β − -decay NMEs and the average OMC matrix elements, especially for the J π = 3 + , 3 − , 4 + , and 4 − states.
In [14], double beta decays of 106 Cd were studied in the pnQRPA framework using 40 Ca as the inert core. Here we extend those studies, for the ground-state-to-ground-state transition, by making a comparative analysis of the 0νβ + β + -decay and OMC NMEs of 106 Cd in the pnQRPA formalism with large no-core single-particle bases, in a manner pursued in [13]. The OMC on 106 Cd leads to excited states of 106 Ag which, on the other hand, act as virtual intermediate states of the β + β + decay of 106 Cd (see Figure 1). Hence, we hope the comparison between OMC and 0νβ + β + -decay matrix elements will help improve the accuracy of the 0νβ + β + -decay NMEs by using the data of future muon-capture experiments. Particularly, in the case of 106 Cd, a measured OMC-strength spectrum would help pin down the value of the particle-particle parameter g pp of pnQRPA, which in this case cannot be adjusted to 2νββ-decay data. We also study the effect of different bound-muon wave function on the OMC matrix elements. This is the first time when such calculations are being done on the positron-decay side of the nuclear chart.

TWO-NEUTRINO DOUBLE-BETA DECAY
The half-life of a ground-state-to-ground-state two-neutrino double-beta decay can be written in the form where g eff A is the effective value of the weak axial-vector coupling strength. The factor G 2ν is a leptonic phase-space factor (in units of inverse years) defined in [15,16]. The ground states of the initial and final nuclei are denoted by 0 + i and 0 + f , correspondingly. The Gamow-Teller NME involved in Equation (1) can be written as with the energy denominator where is the nuclear mass difference between the initial and final 0 + ground states, M i the mass of the initial nucleus, and m e the electron rest mass.Ẽ(1 + m ) and E(1 + m ) are the (absolute) energies of the mth 1 + state in a pnQRPA calculation based on the left-and right-side ground states.
In principle, the expression in Equation (2) could also contain a Fermi NME. However, the ground states of the mother and daughter nuclei belong to different isospin multiplets, and due to the isospin symmetry, the Fermi contribution to the 2νββdecay NME should vanish, leaving the Gamow-Teller NME in Equation (2) as the sole contributor to the 2νββ-decay rate.
In [16], the phase-space factors for the 2νβ + β + decay, as well as for the competing modes 2νβ + EC and 2νECEC were computed. In [8], the experimentally extracted lower limits for half-lives of the different modes for 106 Cd were given. These values, together with the corresponding phase-space factors and resulting experimental matrix elements, are listed in Table 1.

NEUTRINOLESS DOUBLE-BETA DECAY
We assume that the 0νββ decay is dominated by the light-Majorana-neutrino-exchange mechanism, and exploit the formalism presented in [17]. Here we are only interested in the ground-state-to-ground-state transitions. The half-life for such a 0νββ transition can be written as where G 0ν is a phase-space factor for the final-state leptons in units of inverse years (see [15,16]), defined here without including the axial-vector coupling g A . The effective lightneutrino mass, m ν , of Equation (4) is defined as with m j being the mass eigenstates of light neutrinos. The amplitudes U ej are the components of the electron row of the light-neutrino-mass mixing matrix. The 0νββ-decay NME M (0ν) in Equation (4) is defined as where we adopt the CVC value g V = 1.0 for the weak vector coupling strength. The definitions for the double Fermi, Gamow-Teller, and tensor NMEs M (0ν) F , M (0ν) GT , and M (0ν) T can be found e.g., in [17].

MUON CAPTURE
Ordinary muon capture (OMC) is a semileptonic weak interaction process, quite like electron capture (EC). The OMC process we are interested in here can be written as where the negatively charged muon (µ − ) is captured by the 0 + ground state of the even-even nucleus X of mass number A and atomic number Z. The process leads to the J π multipole states of Y, the odd-odd isobar of the mother nucleus, of atomic number Z − 1; here J is the angular momentum and π the parity of the final state. At the same time a muon neutrino ν µ is emitted.

Bound-Muon Wave Functions
The wave function of the muon bound on an atomic orbit of the initial nucleus can be written as an expansion in terms of the normalized spherical spinors χ κµ where G κ and F κ are the large and small radial components of the wave functions of the bound state [21]. Here κ denotes the atomic orbit in the following manner l = κ and j = l − 1 2 , for κ > 0 l = −κ − 1 and j = l + 1 2 , for κ < 0.
After being stopped in the outermost shell of an atom, the negatively charged muon undergoes a cascade of transitions to lower atomic orbitals, leaving it finally on the lowest, K atomic orbit. Hence, the captured muon can be assumed to be initially bound in the lowest state, 1s 1/2 , corresponding to κ = −1 and µ = ± 1 2 . Making this assumption, we can estimate the boundmuon wave function by the Bethe-Salpeter point-like-nucleus approximation formulae [22] where γ is defined as where α is the fine-structure constant and Z the atomic number of the nucleus. The Bohr radius of the µ-mesonic atom is where we have adopted the valuesh = c = 1. The is the reduced muon mass in the µ-mesonic atom. If we assume that αZ is very small, γ ≈ 1, and therefore Alternatively, we can reconstruct a realistic bound-muon wave function by solving from the Dirac wave equations the large, G −1 , and small, F −1 , parts of the wave function (9) in the Coulomb field created by the nucleus. If we assume that the muon is in the lowest state 1s 1/2 (κ = −1), the components satisfy the coupled differential equations (see, e.g., [23], but note that they use different notations for the large and small parts) Furthermore, "pl" refers to a point-like nucleus and "fs" to a finite-size nucleus.
Assuming finite nuclear size and uniform distribution of the nuclear charge within a charge radius R c = r 0 A 1/3 with r 0 = 1.2 fm, the potential energy V(r) in Equation (14) can be written in the form similarly as in [15,16,24] in the case of bound-electron wave functions in the context of double beta decay. Equation (14) can then be solved by means of the package RADIAL [25] by using a piecewise-exact power-series expansion of the radial functions, which then are summed up to a prescribed accuracy.
In Figure 2, we compare the large component, G −1 , of the bound-muon wave function, computed using this method (blue line), with the approximate wave function (black line) of Equation (13). For the sake of comparison, we have also plotted the exact solution of the Dirac equation corresponding to pointlike nucleus (red line). The exact solution for the point-like nucleus is notably close to the Bethe-Salpeter approximation, whereas the finite-size-nucleus solution differs significantly from the point-like-nucleus solution, especially at r ≤ 7 fm.

Muon-Capture Matrix Elements
We compute the OMC matrix elements using the formalism that was originally developed by Morita and Fujii [21]. This formalism takes into account both the genuine and induced vector and axialvector weak nucleon currents. The OMC rate from a J i initial state to a J f final state can be written as where (17) and the Q-value of the OMC process can be computed from [21]. Here M f (M i ) is the nuclear mass of the final (initial) nucleus, m e the rest mass of an electron and E X the excitation energy of the final J π state. The expectation value in Equation (16) can be written as where j ′ is the angular momentum of the bound muon. The definitions of the matrix elements M (i) νu and the corresponding coefficients C (i) can be found in Table 5.1 of [26] (or in Table  1 of [21]). However, note that we use different notation for the coupling constants.
The factors C (i) contain the usual weak vector and axial-vector couplings g V ≡ g V (q) and g A ≡ g A (q) at finite momentum transfer q > 0. The conserved vector current (CVC) and partially conserved axial-vector current (PCAC) hypotheses give the values g V (0) = 1.00 and g A (0) = 1.27 for a free nucleon at zero momentum transfer, and for finite momentum transfer we can use the dipole approximation [27]. For the induced pseudoscalar coupling g P the Goldberger-Treiman PCAC relation [28] gives g P /g A = 7.0. However, at zero momentum transfer deviations from the CVC and PCAC values have been obtained in several analyses [29][30][31][32].
The matrix elements M (i) νu in Equation (19) consist of radial integrals of different integrands containing spherical harmonics, geometric factors and components of neutrino and muon wave functions. The different terms are listed in Table 5.1 of [26]. We assume that the muon is bound on the κ = −1 orbit and that the small component of the bound muon wave function is negligible, which simplifies the expressions of the matrix elements notably [see Equation (5.26) and Table 5.2 of [26] for the Bethe-Salpeter approximation and Equation (5.37) and Table 5.3 for a general muon wave function].
Here we define an average OMC matrix element as and compare this quantity, instead of OMC rate, with the 0νβ + β + -decay NME in order to reduce the phase-space effects.
In this work we choose the slightly quenched values of g A (0) = g eff A = 1.0 and g P (0) = 7.0 and keep the CVC value g V (0) = 1.0 for all the studied cases. In general, OMC serves as a probe of the effective values of these parameters at the momentumexchange region g ≈ 100 MeV, which is particularly relevant for 0νββ decay. These parameter values could be constrained by the measured capture rates to individual excited nuclear states, especially in light nuclei which are well-described by shell-model and different ab initio methods. However, for the lack of the OMC data on the OMC in 106 Cd, this is not possible in the present study.

SPHERICAL PROTON-NEUTRON QRPA AND ITS HAMILTONIAN PARAMETERS
The results reported in the present study are based on a spherical proton-neutron quasiparticle random-phase approximation (pnQRPA), which describes nuclear excitations in odd-odd nuclei (such as 106 Ag) as proton-neutron quasiparticle pairs. In order to reach wide region of excitationenergies up to 50 MeV, we use large no-core single-particle bases consisting of proton and neutron states from the oscillator shells. As a starting point, the single-particle energies were generated by a spherical Coulomb-corrected Woods-Saxon (WS) potential with the parametrization of [33]. This basis is denoted as "WS" in this study. Furthermore, we modified the WS energies in the same way as in [34] in order to better reproduce the spectra of the neighboring odd-A nuclei. This basis, in turn, is denoted by "Adj.".
For the muon-capture calculations, we generate the pnQRPA excitations in 106 Ag starting from the even-even mother nucleus, 106 Cd. As for the 0νβ + β + decay, we generate two sets of pnQRPA excitations for each J π multipole-one based on the mother nucleus 106 Cd, and one based on the daughter nucleus 106 Pd. We call these sets the right-and left-hand pnQRPA sets, correspondingly. In the 0νβ + β + -decay calculations, we then use the average of the left-and right-hand-side excitation energies as the excitation energy of a given intermediate state. We also take into account the overlap of these two sets in the definition of the matrix element.
The quasiparticle spectra for protons and neutrons, needed in the pnQRPA diagonalization, are obtained by solving the BCS equations for protons and neutrons in the even-even reference nuclei. The calculated BCS pairing gaps are adjusted to the phenomenological proton and neutron pairing gaps in a way described in detail in [35].
The X and Y amplitudes in the pnQRPA equations are calculated by diagonalizing the pnQRPA matrix separately for each multipole J π . We adopt as the two-body interaction the one derived from the Bonn-A one-boson-exchange potential, introduced in [36]. The particle-hole part was scaled by a common factor g ph fixed by fitting the centroid of the Gamow-Teller giant resonance (GTGR) in the 1 + channel of the calculations.
As for the particle-particle parameter g pp , we follow the partial isospin-restoration scheme introduced in [37], and multiply the isoscalar (T = 0) and isovector (T = 1) parts of the particle-particle G-matrix elements by factors g (T=0) pp and g (T=1) pp , respectively. The isovector parameter g (T=1) pp is adjusted such that the Fermi part of the corresponding two-neutrino 2 | Adopted values of the pairing parameters of the BCS and the particle-hole and particle-particle parameters of the pnQRPA. The selection procedures of these values are explained in the text.
double beta (2νβ + β + ) NME vanishes, leading to partial isospinsymmetry restoration. Usually, for the 2νβ − β − decays, the value of the isoscalar parameter g (T=0) pp is determined by fitting the corresponding experimental half-life. Here we do not have a measured half-life and the corresponding experimental NME available (see Table 1) so that the value of g (T=0) pp is a matter of choice. In this work we choose to adopt the rather wide range of values g (T=0) pp = 0.6 − 0.8 for this parameter, where the upper limit is at a safe distance from the collapse point of the pnQRPA for both single-particle bases.
All the parameter values resulting from the above-described procedures are listed in Table 2.

RESULTS
In order to investigate the possibility of using the OMC as a probe of 0νββ decay, we have studied the 0νβ + β + -decay matrix elements and the average OMC matrix elements of 106 Cd in the pnQRPA framework in detail. The results are presented in the following subsections. In order to make the comparison between the two processes meaningful, we need to adjust the excitation energies in 106 Ag (being the intermediate nucleus of 0νβ + β + decay and the final nucleus of OMC) in a consistent manner. For muon capture, the excitation energy of the lowest pnQRPA excited state for each J π multipole is adjusted to the measured excitation energy, when available. For the 0νβ + β + decay we adjust the right-hand pnQRPA set of states in a similar manner. Hence, we can compare the two processes as a function of the excitation energy in 106 Ag in a consistent way.

Multipole Decompositions of the Total 0νββ-Decay and OMC Matrix Elements
In contrast to the 2νββ decay which has only the J π = 1 + states active in the process, in the case of the 0νββ decay all the multipole states J π of the intermediate nucleus are active. On the other hand, in the OMC the large mass of the captured muon allows highly forbidden transitions to all possible J π final states up to highest excitation energy. Hence, by studying the OMC on relevant nuclei, one can access the intermediate states of 0νββ decay by complementary means. In this section we investigate the multipole decompositions of the 0νβ + β + -decay and OMC matrix elements.
In Figures 3, 4, we show the multipole decompositions of the total 0νβ + β + -decay NME and the average OMC matrix element of 106 Cd, respectively. The compositions correspond to the parameter value g (T=0) pp = 0.7, which is at a safe distance from the pnQRPA breaking point. In the case of the 0νβ + β + decay, the 1 + multipole plays a dominant role both for the bare Woods-Saxon basis and for the adjusted basis. This ensues from the dominating role of the Gamow-Teller type of transitions. The 1 − contribution is the second largest, whereas the 0 + and 0 − contributions are negligible. The contributions coming from the higher multipoles decrease rather smoothly as a function of J.
As for the muon capture (see Figure 4), the major part of the average matrix element consists of transitions to the states with 1 ≤ J ≤ 4, the leading multipoles being those with J π = 2 − , 2 + , 1 − , and 3 + . In contrast to 0νβ + β + decay, the strength is more evenly distributed among the few leading multipoles. On the other hand, the multipoles with J ≥ 6 play a minor role compared with the 0νβ + β + decay.

Dependence of the Matrix Elements on the Single-Particle Bases and g pp
As mentioned in section 5, the particle-particle strength parameter g pp strongly affects the ββ-decay rates. In the muoncapture studies of e.g., [38,39] it was found that g pp affects also the muon capture rates. Hence, in this section we study the effect of g pp on both the 0νβ + β + -decay and OMC matrix elements of 106 Cd in detail. The exploration also paves the way for possible future adjustments of g pp using, e.g., the shape of the OMC strength function.
The total 0νβ + β + -decay NME of 106 Cd is plotted as a function of the isoscalar part of the particle-particle parameter pp in different single-particle bases in Figure 5. For the isovector part g (T=1) pp , we adopt the value that was adjusted so that the Fermi part of the 2νβ + β + -decay NME vanishes, as explained in section 5. It is evident from the plots that the value of M (0ν) is sensitive to the value of g (T=0) pp : increasing g (T=0) pp decreases the value of the matrix element. In both bases, varying g pp from 0.6 to 0.8 reduces the value of M (0ν) by some 25%. Hence, constraining the value of g pp is of utmost importance in the 0νβ + β + -decay studies. In the absence of a measured 2νβ + β + -decay half-life, adjusting g pp to OMC data, once measured, would help reduce the large uncertainty related to g pp . It is seen that in the adjusted ("Adj.") basis, the matrix element is consistently about 20% larger than in the bare Woods-Saxon ("WS") basis. It is to be noted that soon after the g (T=0) pp values shown in the x axis, between g (T=0) pp 0.82 − 0.84, depending on the single-particle basis, the ground state of pnQRPA collapses and the value of M (0ν) blows up: the value rapidly increases by some 10%. Since the values after the pnQRPA breaking point are not physically meaningful, they are not shown in the figure.
In Figure 6, we plot the average OMC matrix element |M (µ) | av as a function of g (T=0) pp . Contrary to the 0νβ + β + -decay NME, the OMC matrix element is not sensitive to the small adjustments of the single-particle bases or to the value of the   The different parts of the 0νβ + β + -decay NMEs in different single-particle bases in the adopted ranges of the parameter g (T=0) pp (see Table 2) are listed in Table 3. In Table 4, we list the corresponding values obtained in [14] in a smaller single-particle basis, which corresponds to the "Adj." basis in the present study.
The total 0νβ + β + -decay NMEs computed in the present study ( Table 3) are consistently larger than those computed in the smaller, core-based single-particle bases, used in [14]. The values computed with the UCOM short-range correlations are closer to the values computed in the present study, which is natural since the presently adopted short-range correlations are FIGURE 5 | Dependence of the 0νβ + β + -decay NME on the isoscalar part of the particle-particle parameter g pp in the two different single-particle bases.
FIGURE 6 | Dependence of the average OMC matrix element on the isoscalar part of the particle-particle parameter g pp in the two different single-particle bases, and for the exact Dirac finite-nucleus and approximate Bethe-Salpeter point-like-nucleus bound-muon wave functions. a parametrization of the UCOM correlator. The dependence of the NMEs on the size of the single-particle bases is in accordance with the findings of our previous work [40], where we noticed that the size of the single-particle bases affects the 0νβ − β − -decay NMEs much more than the different adjustment procedures of the particle-hole parameter g ph of the pnQRPA. This is most likely due to the fact that in the smaller bases we cannot reach the highly-excited intermediate states which play a non-negligible role in the 0νββ-decay process.
According to Figures 3, 4 the multipoles J π = 1 + and 2 + are among the leading ones for both the 0νβ + β + -decay and OMC The matrix elements are computed in different single-particle bases with different values of gpp.
TABLE 4 | Nuclear matrix elements for 0νβ + β + decay of 106 Cd corresponding to g A (0) = g eff A = 1.0 in the single-particle bases used in [14], corresponding to the "Adj." basis of the present study. The results are adopted from [14].
matrix elements of 106 Cd. Hence, it is illuminating to study the effect of g pp on these multipoles in more detail. In Figures 7, 8 we plot the total 0νβ + β + -decay NME on the positive y axis and the average OMC matrix element on the negative y axis in the cases of J π = 1 + and 2 + , respectively. We decompose the average OMC NME for the J π = 1 + , 2 + multipole states within MeV energy bins while for the 0νβ + β + decay the energy-multipole decomposition entails division of the NMEs into multipoles J π = 1 + , 2 + and their energy distributions binned by MeV-energy intervals. We have chosen to plot only the absolute values of the matrix elements since they carry the essential information needed in the present comparison of the basic features of the OMC and 0νβ + β + decay. Figures 7, 8 show similar behavior for both the multipoles J π = 1 + and 2 + : decreasing the value of g (T=0) pp from 0.8 to 0.6 shifts the spectrum to lower energies for both the 0νβ + β + decay and OMC. Note that the present comparison does not reflect the results of Figure 5, since here we are considering the absolute values of the 0νβ + β + -decay NME for each bin. However, these figures show that even though the average OMC matrix element is quite independent of the value of g pp , the shape of the OMC spectrum depends on g pp . This, in turn, raises interest of studying the possibility of adjusting g pp to the locations of OMC giant resonances, once measured.

Dependence of the OMC Matrix Elements on the Bound-Muon Wave Function
As mentioned in section 4.1, the OMC matrix elements have usually been computed by approximating the bound-muon wave function by a point-like-nucleus approximation. In our previous works [13,41,42] we used the Bethe-Salpeter pointlike-nucleus approximation formula for the muon wave function. Here we study the effect of the exact muon wave function, solved from the Dirac equation by taking into account the finite size of the nucleus with uniform charge distribution, on the OMC matrix elements. The effects of the two muon wave functions (see Figure 2) on the OMC matrix element are clearly seen in Figure 6, where the g pp dependence of the OMC matrix element is displayed in the two different single-particle bases and for the exact Dirac finite-nucleus and approximate Bethe-Salpeter point-like-nucleus bound-muon wave functions. The difference between the matrix elements is considerable. However, one has to keep in mind that in the OMC-rate calculations this difference is to a major part compensated by the use of a phenomenological effective charge Z eff (the so-called Primakoff approximation [43]) in the calculations using the Bethe-Salpeter approximation.
Next we study the effect of the different bound-muon wave functions on the multipole decomposition of the average OMC matrix elements in the cases of the multipoles J π = 1 + , 2 + , which are among the leading ones for both the 0νβ + β + decay and the OMC. In Figures 9, 10 we plot the energy-decompositions of the average OMC matrix element for the transitions to J π = 1 + and 2 + states, respectively. In both figures, the Bethe-Salpeter point-like-nucleus approximation (blue bars) results in notably larger values of the average OMC matrix element than the Dirac wave function (black and white bars). All in all, the use of the Dirac wave function results in about 50-60% reduction of the matrix elements in all the energy bins. This makes sense, since looking at Figure 2, especially at r ≤ 7 fm, the behavior of the Dirac wave function, taking into account the finite size of the nucleus, differs significantly from the Bethe-Salpeter approximation. The finding is also in keeping with results for the total OMC matrix element, depicted in Figure 6.

Comparison of the 0νβ + β + and OMC Matrix Elements
Here we finally compare the absolute values of the 0νβ + β +decay and average OMC matrix elements in the same manner as in [13]. We analyze the summed absolute values of the matrix elements in the same way as we did in Figures 7, 8. We plot the summed absolute values of the 0νβ + β + -decay NMEs and the average OMC matrix elements for J π = 0 + , 1 + , 2 + , 3 + , 4 + , 1 − , 2 − , 3 − , and 4 − in Figure 11. The matrix elements are computed in the adjusted Woods-Saxon basis ("Adj.") with the parameter value g (T=0) pp = 0.7. The OMC matrix elements are computed with the exact Dirac muon wave function. FIGURE 8 | Dependence of the energy distributions of the 0νβ + β + -decay and OMC matrix elements for the J π = 2 + multipole on the particle-particle parameter g Looking at Figure 11, one can see clear correspondences between the 0νβ + β + -decay and OMC matrix elements, especially in the cases of J π = 3 + (Figure 11D), J π = 4 + ( Figure 11E) and J π = 4 − ( Figure 11I). This observation is in accordance with our earlier study in the 0νβ − β − side of double beta decays [13]. There are also notable similarities in the distributions corresponding to multipoles J π = 2 − ( Figure 11G) and J π = 3 − (Figure 11H). For the rest of the multipoles, the correspondencies are not so well visible. Especially, in the case of J π = 1 + (Figure 11B), the major part of the 0νβ + β + -decay NME is coming from the first energy bin E ≤ 1 MeV, while the OMC distribution is clearly more spread to higher energies. This is also the most notable difference between the present results and those of our earlier study [13], where the 1 + contributions to the 0νβ − β − -decay matrix elements were more evenly distributed to higher excitation energies.

DISCUSSION
Double beta decay is one of the most intensively studied topics in neutrino, nuclear and particle physics. While the ordinary two-neutrino double beta decay mode has been observed in several isotopes, the neutrinoless decay mode remains hypothetical. Most of the observed decays are of β − β − type, and there are only six isotopes known to be capable of β + β + decaying. Here we have studied a particularly promising candidate: 106 Cd, for which currently only the lower limit of the 2νβ + β + -decay half-life has been extracted experimentally.
In the present work, we have made a comparative analysis of the 0νβ + β + -decay and average OMC matrix elements of 106 Cd in the pnQRPA framework using large no-core singleparticle bases. This comparison is the first ever done on the positron-emission side of the nuclear chart, and could potentially help improve the accuracy of the 0νβ + β + -decay matrix elements once access to the data of future muoncapture experiments is gained. In particular, adjusting the g pp parameter to future data on OMC giant resonances could help reducing the sizeable uncertainty related to the unknown value of g pp .
Analysis of the multipole decompositions of the total 0νβ + β + -decay matrix element and the average OMC matrix element shows that the J π = 1 + multipole has a dominating role in the 0νβ + β + -decay process, while the total OMC matrix FIGURE 9 | Average OMC matrix elements for the captures to J π = 1 + states computed with different bound-muon wave functions. "B-S(pl)" denotes the Bethe-Salpeter point-like-nucleus approximation and "Dirac(fs)" the exact wave function solved from the Dirac equation taking into account the finite size of the nucleus. The value g (T=0) pp = 0.7 is adopted for particle-particle interaction parameter. (A) Basis: "Adj.". (B) Basis: "WS". FIGURE 10 | Average OMC matrix elements for the captures to J π = 2 + states computed with different bound-muon wave functions. "B-S(pl)" denotes the Bethe-Salpeter point-like-nucleus approximation and "Dirac(fs)" the exact wave function solved from the Dirac equation taking into account the finite size of the nucleus. The value g (T=0) pp = 0.7 is adopted for particle-particle interaction parameter. (A) Basis: "Adj.". (B) Basis: "WS". element is more evenly distributed to a few leading multipoles. The multipoles J π = 1 + and 2 + play a major role in both processes, hence we have studied the transitions involving those multipoles in more detail: we have studied the effect of different particle-particle parameter values on both the 0νβ + β + -decay and the OMC matrix elements, and the effect of different bound-muon wave functions on the OMC matrix elements in these cases.
Our studies indicate that the 0νβ + β + -decay matrix element of 106 Cd strongly depends on the value of the isovector FIGURE 11 | Multipole decompositions in terms of relative 0νβ + β + -decay matrix elements (positive y axes) and average matrix elements of the OMC (negative y axes) for 106 Cd as functions of the excitation energy E in the intermediate nucleus ( 106 Ag) of the 0νβ + β + decay of 106 Cd. These matrix elements correspond the adjusted Woods-Saxon basis and g (T=0) pp = 0.7. For the bound-muon wave function we have used the realistic exact wave function solved from the Dirac equation.
Here J π refer to the angular momenta and parities of the virtual states in 106 Ag and all quantities have been summed within 1 MeV energy bins. The subfigures represent different J π values. The scale values of the y axes have been omitted, since they are not relevant for the current analysis. For more information see the text. part g (T=0) pp of the particle-particle interaction parameter of pnQRPA. Contrary to this, the average value of the OMC matrix element is less dependent on the g (T=0) pp . However, near the pnQRPA breaking point the average OMC matrix element becomes unstable and grows fast in magnitude. Furthermore, when comparing the 0νβ + β +decay matrix elements of the present work with those computed in (much) smaller single-particle bases in [14], we noticed that the matrix elements are sensitive to the size of the single-particle basis. This observation is in accordance with our earlier work on the β − β − type of decays [13].
Finally, we compared the energy distributions of the multipole-decomposed 0νβ + β + -decay matrix elements and the average OMC matrix elements, computed in the adjusted Woods-Saxon single-particle basis. We identified a clear correspondence between the absolute values of the 0νβ + β + -decay and OMC multipole contributions, especially in the cases of the J π = 3 + , 4 + , and 4 − multipoles. This finding is in accordance with our previous work [13], where we compared the energy distributions of the multipole-decomposed 0νβ − β − -decay and OMC matrix elements for several 0νβ − β − -decay triplets in a similar manner.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.