scieee AI-readable full text Open interactive document viewer

Impact of ultrasoft neutrinos and realistic radial wavefunctions on neutrinoless double-beta decay matrix elements

Soriano Fajardo, Pablo

Full text

Impact of ultrasoft neutrinos and realistic radial wavefunctions on neutrinoless double-beta decay matrix elements Master’s degree in Nuclear Physics Master’s thesis by Pablo Soriano Fajardo Advisor: Dr. Javier Men´endez S´anchez (UB, Barcelona) University of Seville (US) University of Barcelona (UB) November 2020 1 Contents 1 Introduction 2 1.1 Neutrinoless double-beta decay and neutrino physics . . . . . . . . . . . . 2 1.2 Thenuclearshellmodel ............................ 5 1.2.1 From the independent particle model to the interacting shell model 5 1.2.2 TheLanczosmethod.......................... 7 1.2.3 The choice of the basis . . . . . . . . . . . . . . . . . . . . . . . . . 8 2 Contribution of ultrasoft neutrinos to 0νββ decay 10 2.1 Conventional formalism of the 0νββ decay amplitude . . . . . . . . . . . . 10 2.2 Contribution of ultrasoft neutrinos (|k|  kF) ................ 12 2.3 Numerical calculations of M0ν usoft in ββ decaying nuclei . . . . . . . . . . . . 14 2.3.1 Calculating transition matrix elements . . . . . . . . . . . . . . . . 14 2.3.2 Calculating M0ν usoft ............................ 15 2.3.3 Sensitivity of M0ν usoft to the nuclear interaction . . . . . . . . . . . . 19 2.3.4 Computationally complex nuclei: 0νββ decay of 128Te and 130Te . . 19 2.4 Energy dependence: comparison to 2νββ decay ............... 21 3 Realistic radial wavefunctions in NME calculations 23 3.1 Harmonic oscillator wavefunctions . . . . . . . . . . . . . . . . . . . . . . . 23 3.2 Woods-Saxon wavefunctions . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.2.1 Separationenergies........................... 32 3.3 Transitiondensities............................... 33 3.4 Radial distribution of the NME components . . . . . . . . . . . . . . . . . 36 4 Summary and conclusions 41 A Appendix A: Correction to the electron energies 43 B Appendix B: Convergence of S(nmax)and M0ν usoft 44 C Appendix C: Associated Laguerre polynomials 47 D Appendix D: Oscillator amplitudes A(nlj) ν48 E Appendix E: Talmi-Moshinsky transformation for Woods-Saxon wavefunctions in HO basis 52 2 1 Introduction Precise calculations of the neutrinoless double-beta 0νββ decay nuclear matrix elements (NME) used to calculate the 0νββ decay rate are relevant to the study of neutrino physics, since this decay rate also depends on a combination of the neutrino masses and mixing matrix elements. Our goal will be to improve the calculations of these NME by including previously neglected precision components, such as the contribution of very low-momentum (ultrasoft) neutrinos (Chapter 2) and the use of realistic radial wavefunctions (Chapter 3). In this starting chapter, we will introduce and discuss the basic concepts that will make up the foundations of this work, as well as its objectives. In the first section (1.1) we will deal with the concept of neutrinoless double-beta decay and its significance and interest towards neutrino physics. In the second section (1.2), we will provide an insight of the nuclear shell model, which will be the framework of our calculations. 1.1 Neutrinoless double-beta decay and neutrino physics Double-beta decay is a transition between isobaric nuclei in which two neutrons simultaneously decay into protons, meaning that the parent nucleus decays into a daughter nucleus with two fewer neutrons and two more protons. Being a second-order weak-interaction process, it is strongly suppressed and only observable for isotopes in which single beta decay is forbidden. This results in typical half-lives ranging from 1018 to 1024 yr, as observed in 11 different nuclei undergoing ββ decay [1] and 3 undergoing double electron capture εε [2]. So far, all measured double-beta decay processes have been two-neutrino double-beta (2νββ) decay (corresponding to Figure 1a), which can be expressed as A ZX−−→ A Z+2Y + e− 1+ e− 2+ ¯νe,1+ ¯νe,2,(1) where A ZX and A Z+2Y are the parent and daughter nuclei; A,Z and N are the number of nucleons, protons and neutrons (such that N=A−Z), e− 1,2are the emitted electrons and ¯νe;1,2are the emitted antineutrinos. One particular detail about neutrinos is that they are the only neutral fermions we know of. As such, neutrinos are the only known particles that may be Majorana fermions, which means that they would be their own antiparticle. This would imply the existence of a yet-hypothetical alternate decay mode, the neutrinoless double-beta (0νββ) decay A ZX−−→ A Z+2Y + e− 1+ e− 2(2) As their own antiparticle, Majorana neutrinos could virtually annihilate each other in the decay, as depicted in Figure 1b. The inverse half-life of a 0νββ decay between JP= 0+states of the parent and daughter nuclei can be written as T0ν 1/2(0+ i→0+ f)−1=G0ν(Qββ, Z)g4 A|M0ν|2m2 ββ ,(3) 3 (a) Two-neutrino double-beta decay (b) Neutrinoless double-beta decay Figure 1: Feynman diagrams for (a) 2νββ and (b) 0νββ decay. Two neutrons (n) decay into two protons (p), emitting two electrons (e) and (a) two antineutrinos (¯ν) or (b) no neutrinos, implying that they are Majorana particles (νM) in the second case. From [3]. where G0ν(Qββ, Z) is a phase space factor that can be calculated with great precision [4], Qββ =Ei−Ef−2me(4) is the Q value of the reaction, M0νis the ”nuclear matrix element” (NME) M0ν=M0ν GT −g2 V g2 A M0ν F+M0ν T,(5) with each of the terms corresponding to the Gamow-Teller, Fermi and tensor NME contributions, respectively. gV= 1 and gA≃1.27 are the vector and axial coupling, and mββ a combination of the neutrino masses mjand mixing matrix elements Uej, defined as mββ ≡X j mjU2 ej =m1|Ue1|2+m2|Ue2|2ei(α2−α1)+m3|Ue3|2ei(−α1−2δ).(6) where mjare the neutrino mass eigenstates (not the same as the neutrino flavor eigenstates: m1leaning heavily towards electron flavor, m2being a more even blend of the three lepton flavors and m3mostly muon and tau flavor), δis the so-called Dirac phase and α1,2are Majorana phases that vanish if neutrinos are Dirac particles [3]. Likewise, we can express the inverse half-life for 2νββ decay as T2ν 1/2−1=G2ν(Qββ, Z)g4 A M2ν GT −g2 V g2 A M2ν F 2 ≃G2ν(Qββ, Z)g4 AM2ν GT  2.(7) where the Fermi contribution can be safely neglected when the isospin Tof the initial and final nuclei is different. 4 Obtaining reliable values of M0νis key to obtain the value of mββ, since experimentally observing a 0νββ decay would provide its half-life, T0ν 1/2. This would provide vital information about the absolute scale of the neutrino masses and, potentially, about its ordering. As illustrated in Figure 2, knowing the value of mββ along with the mass of the lightest neutrino (a parameter βdecay experiments and cosmological observations like the KATRIN [5] and Planck [6] collaborations are sensible to) could confirm which of the neutrino-mass hierarchies is correct: the normal hierarchy (m1< m2< m3) or the inverted hierarchy (m3< m1< m2). This would have fundamental implications for neutrino physics and cosmology [7], since it could shed some light on the mechanism of neutrino mass generation. 0νββ decay would violate the conservation of the total lepton number, as well as the conservation of the ”B−L” quantity (baryon-lepton number), the latter being a fundamental symmetry of the Standard Model [1]. Accordingly, the detection of this process would imply Physics beyond the Standard Model. Additionally, the existence of the aforementioned violation of the lepton number could be very useful to explain the generation of the matter-antimatter asymmetry in our universe via a process known as ”leptogenesis” [8,9]. Figure 2: Left panel: bands for the value of the parameter mββ as a function of the mass of the lightest neutrino, for the case of normal (NH, red band) and inverted (IH, green band) neutrino-mass hierarchies. The present best experimental upper limits on mββ are shown in the blue band. Right panel: present best upper limits, with uncertainty bars, on mββ from experiments performed on each ββ emitter, as a function of mass number A. The uncertainty bands and bars include experimental uncertainties and ranges of calculated nuclear matrix elements. From [10]. 5 1.2 The nuclear shell model 1.2.1 From the independent particle model to the interacting shell model The nuclear shell model was originally introduced by M. Goeppert-Mayer [11] and H. Jensen et al. [12] to explain the regularities of the nuclear properties associated with magic numbers. They proposed an independent particle model, assuming that the main effect of the two-body nucleon-nucleon (NN) and three-body (3N) interactions was to generate a mean field [13]. This nuclear mean field was constructed as a surface corrected harmonic oscillator whose main novelty was the very strong spin-orbit splitting needed to explain these magic numbers [14]. It is described as follows: U(r) = 1 2mω2r2+Dl2−Cl·s,(8) where mis the mass of the particle, ωits the angular frequency, rits radial component, land sare the angular momentum and spin operators (respectively) and Cand Dare parameters to be fit for best results. Figure 3shows the single particle levels of the nuclear mean field. From left to right, we first find the shell structure of the harmonic oscillator, then the splitting due to the l2 term in Equation (7) and, in the middle of the figure, the actual single particle levels, that take into account the spin-orbit splitting [14]. On the right hand side we find three rows of numbers corresponding to, respectively, the maximum occupancy of the level, the accumulated occupancy and the predicted magic numbers for protons and neutrons (2, 8, 20, 28, 50, 82 and 126). Alternately to the notation in Figure 3, the principal quantum number nof the lowest level is sometimes defined as n= 0 instead of n= 1, such that N= 2n+l. This is the convention we will be following. The success of the independent particle model strongly suggests that the very singular free NN interaction can be regularized in the nuclear medium. Starting with this bare interaction, the exact solution of the many-body problem in the infinite Hilbert space built on the mean field orbits is approximated in the large scale shell model calculations by the solution of the Schr¨odinger equation in the valence space, using an effective interaction [13,14] such that observables are preserved: H|Ψi=E|Ψi → Heff|Ψeff i=E|Ψeffi.(9) Note that Heff includes interactions between nucleons, therefore describing physics beyond the mean field of Equation (8). In general, effective operators have to be introduced to account for the restrictions of the Hilbert space, such as 6 Figure 3: Structure of the spherical mean field, see text for details. From [14]. hΨ|O|Ψi=hΨeff|Oeff |Ψeffi.(10) Once we adopt a regularized interaction that is compatible with the experimental mean field (in this case, magic numbers) we can proceed using the spherical mean field orbits as the basis for the occupation number space (Fock space). We will have states i, j, k, ... with energies i, j, k, ... that will bunch in shells, giving rise to magic numbers when the energy difference between them is large enough. This motivates the separation of the full space intro three different regions: •Inert core: encompasses the orbits that are always full, therefore not undergoing any change. If this core is made up of Zcprotons and Ncneutrons, there will remain zv=Z−Zcvalence protons and zv=N−Ncvalence neutrons. •Valence space: formed by the orbits that are available to the aforementioned valence particles, which will partially occupy them according to the dictates of the effective interaction Heff. •External space: refers to all the remaining orbits, which always remain empty. In the valence space, a formal solution to the A-body problem can be obtained through the following procedure: •First, we choose a single particle basis: a+ i|0i. 7 •We will then proceed to build the A-particle wavefunctions as Slater determinants: |Φαi= A Y i a+ i|0i.(11) •The physical states are then expressed as a linear combination of these Slater determinants: |Ψeffi=X α Cα|Φαi.(12) •Lastly, the solution of the many-body problem Heff |Ψeffi=E|Ψeffiis given by the eigenvalues and eigenvectors of the many-body matrix hΦα|Heff|Φ0 αi. 1.2.2 The Lanczos method If one aims to carry out large scale shell model calculations, standard diagonalization methods where CPU times increase as N3 dim [15,16] (Ndim being the dimension of the many-body matrix) are not suitable. Taking into account that, in general, only a small number of eigenvalues and eigenvectors are needed, we can try to use another approach to this problem. This is where the Lanczos algorithm comes in handy. The Lanczos method is based in the building of an orthogonal basis in which the manybody matrix has a tridiagonal structure. We initialize the algorithm with a normalized vector Φ1(that we call the ”pivot” state) and apply the effective hamiltonian Heff operator. We will then get a parallel and an orthogonal component of the pivot Φ1: Heff|Φ1i=E11|Φ1i+E12|Φ2i,(13) where E11 =hΦ1|Heff|Φ1iand E12|Φ2i=Heff |Φ1i−E11|Φ1i. Applying Heff on Φ2, we will generate a third vector Φ3orthogonal to Φ1and Φ2: Heff|Φ2i=E21|Φ1i+E22|Φ2i+E23|Φ3i,(14) where the hermicity of Heff implies E21 =E12. Analogously to the previous step, E22 = hΦ2|Heff|Φ2i, and E23 is obtained through normalization: E23|Φ3i= (Heff −E22)|Φ2i−E21|Φ1i.(15) Continuing this process, at iteration n, we obtain the diagonal energy of the vector |Φni, a new vector |Φn+1iand the non-diagonal energy En,n+1: Enn =hΦn|Heff|Φni,(16) Heff|Φni=En,n−1|Φn−1i+Enn|Φni+En,n+1|Φn+1i,(17) 8 En,n+1|Φn+1i= (Heff −Enn)|Φni−En,n−1|Φn−1i,(18) always remembering that En,n−1=En−1,n. Due to Heff being hermitic, the construction of the Lanczos matrix ensures that Eij = 0 if |i−j|>1, thus obtaining a tridiagonal matrix:          E11 E12 0 0 ··· E21 E22 E23 0··· 0E32 E33 E34 ··· 0 0 E43 E44 ··· . . .. . .. . .. . ....          (19) This matrix is diagonalized every certain number of steps. When the eigenvalues obtained through the diagonalization meet an established convergence criterion, the calculation will stop. 1.2.3 The choice of the basis For a given valence space, the choice of the basis is simply driven by convenience. Depending on the properties we want to describe, one or another basis may be more appropriate. There are two possible choices: the m-scheme and the J-coupled scheme. In the m-scheme, the basis is composed of all the Slater determinants made from all the possible partitions of the valence particles among every valence orbit |nljmτi[14], such that |Φαi=Y i=nljmτ a+ i|0i.(20) where n,l,j,mand τare the principal, orbital, total angular momentum, magnetic and isospin quantum numbers, respectively. The principal advantage of this representation is the simplicity of the calculation of the many particle matrix elements, since they reduce to the two-body matrix elements of H in m-scheme with a phase [16]. Its major drawback, however, lies in the fact that only Jz and Tzare good quantum numbers, meaning that the basis takes into account all possible (J,T) states. This causes the dimensions of the matrix to be maximal, being proportional to Ndim ∝dπ zv·dν nv,(21) where dπand dνare the total degeneracies of the proton and neutron valence spaces and zvand nvare the valence protons and neutrons (as explained on section 1.2.1). 15 •Last, we will calculate the overlaps of the nmax |nistates with |GT, fi, obtaining hn|GT, fi. With this, we have everything we need to calculate the (nmax) overlap matrix elements as hf|στ+|nihn|στ+|ii=NiNfhf, GT|nihn|GT, ii.(44) It is important to note that Equation (36) takes into account the complete set of infinite |niintermediate states. As we obviously cannot carry out infinite calculations, we will need to find an appropriate value of nmax that guarantees the convergence of the result. As an additional check, S(nmax)≡ nmax X n=1 hf|στ+|nihn|στ+|ii(45) should converge towards Sfi ≡ hf|στ+στ+|ii=NiNfhf, GT|GT, ii,(46) since the closure relation Pn|nihn|=1implies nmax X n=1 |nihn|nmax→∞ −−−−−→ 1.(47) 2.3.2 Calculating M0ν usoft With the transition matrix elements we just calculated, we only need the corresponding energies from Equation (36) in order to obtain values for M0ν usoft. We will calculate the energy of the initial state (Ei) using the binding energies per nucleon (B/A) that can be found in [20,21]. This way, Ei=Zmp+Nmn−(B/A)A ZW·A , (48) where A ZW is the parent nucleus. To calculate the ground state energy of the intermediate states (that is, the energy of the lowest-lying intermediate state, E1+ 1), we will also need an EGS→1+ 1excitation energy, which we will be taking from [22]. So, it follows that E1+ 1= (Z+ 1)mp+ (N−1)mn−(B/A)A Z+1X·A+EGS→1+ 1 A Z+1X,(49) 16 where A Z+1X is the intermediate nucleus. For experimentally observed double-beta decays, we will take Qββ from the literature [23,24]. Otherwise, we will simply take Qββ =Ei−Ef−2me= 2(mn−mp) + [(B/A)A Z+2Y−(B/A)A ZW]−2me,(50) where Ef= (Z+ 2)mp+ (N−2)mn−(B/A)A Z+2Y·A; (51) A Z+2Y being the daughter nucleus. Now that we have everything we need to calculate M0ν usoft, we will begin by representing S(nmax) for a certain range of values of nmax in order to find an appropriate nmax that guarantees the convergence S(nmax)→Sfi, as we discussed in the previous subsection (2.3.1). In Figures 4and 5we have represented S(nmax) with 10 ≤nmax ≤60 for the double-beta decay of 48Ca and 136Xe using the KB3 and GCN5082 interaction, respectively. As we can see, S(nmax) is well converged to Sfi for nmax = 60. We can also appreciate this convergence if we represent M0ν usoft as a function of nmax, as shown in Figures 6and 7. As a note, we have performed this convergence check for all nuclei and interactions considered in this chapter. For the sake of brevity, we will present its corresponding figures in Appendix B(Figures 25-30). Figure 4: S(nmax) for the double-beta decay of 48Ca, using the KB3 interaction. The value of Sfi is depicted with a red dashed line. 17 Figure 5: S(nmax) for the double-beta decay of 136Xe, using the GCN5082 interaction. Figure 6: Matrix element M0ν usoft with respect to the number 1+of intermediate states nmax for the double-beta decay of 48Ca, using the KB3 interaction. 18 Figure 7: Matrix element M0ν usoft with respect to the number 1+of intermediate states nmax for the double-beta decay of 136Xe, using the GCN5082 interaction. Thus, we will calculate M0ν usoft for two different experimentally observed double-beta decaying nuclei: 48Ca, where we will work in the valence space consisting on 0f7/2, 1p3/2, 0f5/2 and 1p1/2(also known as the pf-shell), using the KB3 interaction; and 136Xe, working in the 0g7/2, 1d5/2, 1d3/2, 2s1/2and 0h11/2valence space and using the GCN5082 interaction. The results are presented in Table 1below. Now, in order to calculate the M0ν usoft/M0νratio, we will use the shell model code NATHAN. Since the intermediate 1+states do not play any role in the calculation, only J= 0 states will be involved, an optimal condition for J-coupled scheme code like NATHAN, as we discussed in subsection 1.2.3. Parent Interaction Ei(MeV) E1+ 1(MeV) Ef(MeV) Qββ (MeV) M0ν usoft 48Ca KB3 44657.259 44658.669 44651.970 4.263 [23] 1.12 ·10−2 136Xe GCN5082 126569.136 126569.306 126565.657 2.458 [24] 9.86 ·10−2 Table 1: Energies, Qββ value and M0ν usoft (calculated for nmax = 60) for the ββ emitters 48Ca and 136Xe. The energies Ei,En,GS and Efare taken from [20,21,22]. Parent Interaction M0ν usoft M0νM0ν usoft/M0ν 48Ca KB3 1.12 ·10−20.93 1.2 ·10−2 136Xe GCN5082 9.86 ·10−22.30 4.3 ·10−2 Table 2: M0ν usoft,M0νand its ratio for the 0νββ decay of 48Ca and 136Xe. 19 As we can see, the ratio M0ν usoft/M0ν∼10−2holds up for both of the calculated cases. In particular, the ultrasoft contribution for the 0νββ decay of 48Ca and 136Xe represents a correction to the conventional NME of around 1.2% and 4.3% respectively. This, in turn, implies a correction to the 0νββ decay rate (Equation (3)) of ∼2.4% and ∼8.6%, which should not be neglected if one wants precise values of the NME. 2.3.3 Sensitivity of M0ν usoft to the nuclear interaction In the previous subsection, for the calculations involving the double-beta decay of 136Xe, we used the GCN5082 interaction, which is an effective hamiltonian based in a renormalized G-matrix obtained from the Bonn-C potential and constructed through a fit to about 300 energy levels from ∼90 nuclei [25], making it a more realistic and sophisticated interaction. However, it might be interesting to repeat the calculations using three different G-matrices based on the different Bonn potentials to see how much they differ from the GCN5082 calculations. In Table 3below, we present the calculated values of M0ν usoft,M0νand its ratio using these Bonn interactions (A through C), including the previously calculated case for the GCN5082 interaction. Parent Interaction M0ν usoft M0νM0ν usoft/M0ν 136Xe Bonn-A 2.88 ·10−21.11 2.6 ·10−2 136Xe Bonn-B 3.52 ·10−21.18 3.0 ·10−2 136Xe Bonn-C 4.05 ·10−21.25 3.2 ·10−2 136Xe GCN5082 9.86 ·10−22.30 4.3 ·10−2 Table 3: M0ν usoft,M0νand its ratio for the 0νββ decay of 136Xe using different nuclear interactions. Table 3shows that the results are pretty similar for the three variants of the Bonn interaction. Even though both M0ν usoft and M0νincrease by around a factor 2 for the GCN5082 interaction with respect to Bonn-C, the ratio between the two values roughly stays the same for all interactions, leading to a very robust M0ν usoft correction. 2.3.4 Computationally complex nuclei: 0νββ decay of 128Te and 130Te In addition to 48Ca and 136Xe, we would have also wanted to study other experimentally observed double-beta decaying nuclei, like 128Te and 130Te. However, as we mentioned in subsection 1.2.3, the dimension of the many-body matrix increases proportionally to the product of two binomial coefficients that depend on the total degeneracy of the valence spaces and the number of valence nucleons. 20 The maximum dimension that can be handled by the version of ANTOINE we used is Ndim,max = 2 ·108, while the dimensions of the matrices needed to calculate the daughter nuclei 128Xe and 130Xe are Ndim(128Xe) = 32 4·32 24∼4·1011 ,(52) Ndim(130Xe) = 32 4·32 26∼3·1010 .(53) We could, however, do a preliminary study and carry out calculations for similar nuclei that are computationally feasible (like 132Te and 134Te) in order to check how the relative contribution of M0ν usoft varies along an isotopic chain. To do this, we will first study Ca isotopes around 48Ca like 46Ca and 50Ca using the same nuclear interaction. Parent Interaction Qββ (MeV) M0ν usoft M0νM0ν usoft/M0ν 46Ca KB3 0.988 [20,21] 1.04 ·10−21.20 8.7 ·10−3 48Ca KB3 4.263 [23] 1.12 ·10−20.929 1.2 ·10−2 50Ca KB3 5.921 [20,21] 9.97 ·10−41.063 9.4 ·10−4 Table 4: M0ν usoft,M0νand its ratio for different Ca isotopes, alongside it’s Qββ value. Table 4shows that the obtained values of M0ν usoft are reasonably similar for 46Ca and 48Ca, but decrease an order of magnitude for 50Ca. However, as can be seen in Figures 4,6, 27 and 28, the individual contributions of each term in the summation are comparable for the three cases, with the exception that, for 50Ca, they cancel each other more often, resulting in an overall lower total value. Therefore, our results suggest that calculations for certain nuclei can qualitatively give us an idea of the results for similar isotopes. Parent Interaction Qββ (MeV) M0ν usoft M0νM0ν usoft/M0ν 132Te GCN5082 4.090 [20,21] 6.34 ·10−22.55 2.5 ·10−2 134Te GCN5082 5.592 [20,21] 8.58 ·10−21.94 4.4 ·10−2 Table 5: M0ν usoft,M0νand its ratio for different Te isotopes, alongside it’s Qββ value. Table 5shows the results for 132Te and 134Te, which are comparable to those obtained for 136Xe. As such, we can expect a similar contribution of M0ν usoft for the 0νββ decay of 128Te and 130Te. 21 2.4 Energy dependence: comparison to 2νββ decay The transition matrix elements hf|Jµ|nihn|Jµ|iifrom Equation (36) that we calculated in the previous subsections also control the 2νββ decay amplitude, although through a different En-dependent weighted sum [9]. In particular, the Gamow-Teller matrix element M2ν GT [26] can be written as M2ν GT =meX n hf|στ+|nihn|στ+|ii En−(Ei+Ef)/2.(54) We will analyze the dependence with respect to Enof M2ν GT and M0ν usoft. In particular, we will represent the running sum with respect to the energy difference between the corresponding intermediate state and the initial state, En−Ei. This way, we can visualize the contribution of each state to M2ν GT and M0ν usoft. We define S2ν GT (n) = mehf|στ+|nihn|στ+|ii En−(Ei+Ef)/2,(55) S0ν usoft (n) = R 2πhf|στ+|nihn|στ+|ii×2Qββ 2+En−Ei ln mπ 2Qββ 2+En−Ei+ 1 , (56) such that M2ν GT =PnS2ν GT (n) and M0ν usoft =PnS0ν usoft (n). Figure 8: Running sum of M2ν GT and M0ν usoft with respect to En−Eifor the double-beta decay of 48Ca, using the KB3 interaction. 22 Figure 8shows that, as expected due to its dependence on En,M2ν GT is more sensitive to low-energy intermediate states (En−Ei<5 MeV) than to the ones with a higher energy difference like those in the 5 MeV < En−Ei<10 MeV region (depicted between red dashed lines for clarity). This implies higher values of S2ν GT (n) for the low-energy region and thus a more abrupt change in the value of the running sum. After that, the values of S2ν GT (n) practically vanish for En−Ei>10 MeV. On the other hand, and following the same reasoning, M0ν usoft is much more sensitive to states in the 5 MeV < En−Ei<10 MeV region than M2ν GT , after which S0ν usoft (n) also rapidly decays for high energies. This is even clearer for the double-beta decay of 136Xe. In Figure 9we can distinguish how M0ν usoft is more sensitive to states between 5 MeV < En−Ei<10 MeV than M2ν GT . This is consistent with the estimation of subsection 2.1, where we said that the energy differences En−Eibetween bound nuclear states have a typical size of O(5 −10) MeV, while the transition matrix elements hf|Jµ|nihn|Jµ|ii(and, consequently, S2ν GT (n) and S0ν usoft (n)) rapidly die out for En−Ei>10 MeV. Figure 9: Running sum of M2ν GT and M0ν usoft with respect to En−Eifor the double-beta decay of 136Xe, using the GCN5082 interaction. 23 3 Realistic radial wavefunctions in NME calculations In order to solve Heff|Ψeffi=E|Ψeff i, the radial part of the wavefunctions is not explicitly assumed, since only matrix elements of Heff are needed and Ψeff is given as a combination of Slater determinants. However, in order to evaluate NME of r-dependent operators, the explicit form of the radial wavefunction is needed. While the harmonic oscillator is the conventional choice, it might be too simple if we want to obtain accurate results. In the following chapter we will study the impact of realistic radial wavefunctions on the value of nuclear matrix elements, as well as the choice of the valence space. In particular, we will present the differences between harmonic oscillator and Woods-Saxon radial wavefunctions and compare them to the results obtained through the more accurate variational Monte Carlo (VMC) calculations. Following an article by X.B. Wang et al. [27], we will study the effects of these radial wavefunctions towards the calculation of the NME of the double-beta decay of light nuclei (A= 10 and A= 12), as well as the differences that arise from using a simple valence space (like the p-shell) or an extended valence space (psd-shell). We extend this study taking into account all three contributions to M0ν(M0ν F,M0ν GT and M0ν T, the last of which is neglected) in order to produce a more accurate result. In addition, for the first time we extend this study towards the double-beta decay of heavier nuclei (A= 48). 3.1 Harmonic oscillator wavefunctions Typically, instead of going through the steps that lead to self-consistency, one simply selects a particular type of mean-field potential. The use of such phenomenological potentials is a practical shortcut taken at the expense of theoretical preciseness [28]. The simplest frequently used potential is the three-dimensional harmonic oscillator (HO) potential: vHO(r) = −V1+kr2=−V1+1 2mNω2r2,(57) where V1and kare parameters to be fitted for best results. Given this potential, the harmonic oscillator wavefunctions gnl(r) are solutions of the radial Schr¨odinger equation −~2 2mN∇2 r−l(l+ 1) r2gnl(r) + −V1+1 2mNω2r2gnl(r) = εnlgnl(r),(58) where 2n+l=N= 0,1,2,3, ... (59) are the radial quantum numbers and 24 εnl =−V1+N+3 2~ω=−V1+2n+l+3 2~ω(60) are the energy eigenvalues. As we previously mentioned in section 1.2, we will be using the convention where the principal quantum number n= 0,1,2,3, ... indicates the number of nodes of the wavefunction. The HO radial wavefunction can explicitly be written as [29] gnl(r) = s2n! b3Γ(n+l+3 2)r ble−r2/2b2L(l+1 2) nr2 b2,(61) where bis the oscillator length and L(l+1 2) n(x) is the associated Laguerre polynomial [30]. The oscillator length, as the name suggests, characterizes the width of the oscillator potential, and it can be expressed as b≡r~ mNω=~c p(mNc2)(~ω)≈197.33 √940 ·~ωfm,(62) where we will take the value of ~ωfrom the Blomqvist-Molinari formula [31], based on a fit to charge radii across the nuclear chart: ~ω=45A−1/3−25A−2/3MeV.(63) Instead of directly using the associated Laguerre polynomials, numerical values of gnl(r) can be obtained via auxiliary functions, as we explain in Appendix C. This way, we can generate the HO radial wavefunctions gnl(r), which we can use, for example, to obtain the nucleon density profile (either for protons or neutrons) ρ(r), defined as ρ(r) = |Ψ(r)|2=X nl a2 nl |gnl(r)|2,(64) where a2 nl is the occupancy of orbitals with quantum numbers nl and the sum runs over all orbitals of the inert core and the chosen valence space. We will use the shell model code NATHAN to calculate these occupancies. Following Wang’s paper [27], we will use 2 different valence spaces for the study of light nuclei: the p-shell (0p3/2and 0p1/2) and the psd-shell (0p3/2, 0p1/2, 0d5/2, 1s1/2and 0d3/2), where the PSDMWK shell model hamiltonian [32,33] will be used, adding a correction that takes 31 We have included the Woods-Saxon energies and oscillator amplitudes of all the other cases we will be treating in this chapter in Appendix D(Tables 15-24). Note that this only refers to neutrons of the parent nuclei and protons of the daughter nuclei, as we are interested in the study of double-beta decay. Figure 14: Proton density profile of 48Ca. The dashed red (blue) line corresponds to shell model calculations in the pf-shell configuration space using HO (WS) radial wavefunctions. The asymptotic behaviour of the density profile has been included. Figure 15: Neutron density profile of 48Ca. The dashed red (blue) line corresponds to shell model calculations in the pf-shell configuration space using HO (WS) radial wavefunctions. The asymptotic behaviour of the density profile has been included. 32 3.2.1 Separation energies Up until now, we have followed the usual parametrization for the Woods-Saxon potential (Eqs. (66)-(68)). However, Wang et al. [27] introduced a different parametrization: they took the p-shell proton and neutron asymptotic wavefunction behaviour to be determined by the corresponding separation energies, while the unbound sd-shell particles are assumed to be bound by ∼0.05 MeV. HO wavefunctions decay rapidly at great distances, with an asymptotic form e−(r/b)2. WS wavefunctions, on the other hand, decay less dramatically, like e−kr [36], with k=√2µBs ~,(76) where µis the reduced mass of the last nucleon and the resulting nucleus and Bs=B(AX) −B(A−1Y) (77) is the separation energy of the last nucleon. We calculate Bsfrom the binding energies of the original (AX) and the resulting nucleus (A–1Y), taken from [20,21]. Thus, using Equation (76), we can find the kparameter for each case and fit the asymptotic behaviour of our p-shell WS wavefunctions by modifying the V0parameter of the WS potential. From now on, the usual parametrization will be referred as Suhonen (S) parametrization, while the new one based on the separation energies will be called Wang (W) parametrization. AX N Bs(MeV) k (fm−1)V0,S (MeV) kS(fm−1)V0,W (MeV) kW(fm−1) 10Be n 6.812 0.543 44.4 0.822 33.7 0.639 10C p 4.007 0.417 44.4 0.782 42.3 0.748 12Be n 3.171 0.374 40.0 0.647 29.7 0.574 12C p 15.957 0.839 51.0 0.909 46.2 0.839 Table 11: Suhonen and Wang WS parametrizations for the nuclei involved in the doublebeta decays 10Be →10C and 12Be →12C. As we can see in Table 11 above, it is not possible to exactly fit kWto kfor most of the cases. This is because, as V0decreases, the single-particle states become more unbound. As such, we fitted kWas close as we could to kwhile keeping all the single-particle states of the p-shell bound by ∼0.05 MeV. However, for the case of the 12C protons, we have been able to fit k=kWsince none of the p-shell states became unbound. We extended this parametrization to the double-beta decay of heavier nuclei like 48Ca. In this case, we fitted the asymptotic behaviour of the wavefunction up to the pf-shell, while making sure that all single-particle states are bound by ∼0.05 MeV. Here, each orbital of the pf-shell has its own V0(which is what we did for the sd-shell in the previous case), 33 which is the reason why we did not present any single V0values in Table 12 as we did in Table 11. AX N Bs(MeV) k (fm−1)kS(fm−1)kW(fm−1) 48Ca n 9.952 0.685 0.900 0.685 48Ti p 11.878 0.749 0.993 0.832 Table 12: Suhonen and Wang WS parametrizations for the nuclei involved in the doublebeta decay 48Ca →48Ti. The Woods-Saxon energies and oscillator amplitudes obtained through Wang’s parametrization have been included in Appendix D. 3.3 Transition densities Prior to begin studying the nuclear matrix elements, we will first compare the wavefunction transition densities, defined as C(r) = hf|X a<b δ(r−rab)τ+ aτ+ b|ii.(78) Shell model calculations generally use HO radial wavefunctions, since the HO basis provides the advantage of a simple well-defined method of separating into center-of-mass and relative coordinates via a Talmi-Moshinsky transformation [37,38]. For that reason, we built a new option which, introducing the oscillator amplitudes we obtained in section 3.2, finally allows us to obtain results with WS wavefunctions. To do this, we also generalized the Talmi-Moshinsky transformation from HO wavefunctions to WS wavefunctions in HO basis. Following Equation (72), we express the two-body WS wavefunction in HO basis as |n1l1j1t1, n2l2j2t2;JiWS =X ν1ν2 A(n1l1j1) ν1A(n2l2j2) ν2|ν1l1j1t1, ν2l2j2t2;JiHO (79) so we can expand the corresponding Talmi-Moshinsky transformation as WShn0 1l0 1j0 1p, n0 2l0 2j0 2p;J|Oα|n1l1j1n, n2l2j2n;JiWS = (80) =X ν1ν2ν0 1ν0 2 A0 ν0 1 (n0 1l0 1j0 1)A0 ν0 2 (n0 2l0 2j0 2)A(n1l1j1) ν1A(n2l2j2) ν2HOhν0 1l0 1j0 1p, ν0 2l0 2j0 2p;J|Oα|ν1l1j1n, ν2l2j2n;JiHO where Oα=τ+ 1τ+ 2SαHα(r)˜ Hα(R), with r=|r1−r2|the relative radial coordinate (distance between decaying neutrons) and R=|r1+r2|the center-of-mass radial coordinate. The full expansion of this transformation can be found in Appendix E. 34 Figure 16: The normalization densities C(r), for (a) A= 10, and (b) A= 12. The N= 1 model space (p-shell only) and extended model space of psd shell model calculation results are shown. The different choices of radial wavefunctions, HO, and WS, are also shown. VMC results with shell-model-like wavefunctions are labeled as “VMC-1”, and those with cluster-like wave functions are labeled as “VMC-2”. Taken from [27]. Figure 17: Our rendition of the normalization densities C(r) for A= 10 and A= 12. The green dotted and dashed lines would correspond to our version of the WS data in Figure 16, using the same parametrization as [27]. The blue lines would correspond to the usual WS parametrization. 35 Figure 16 shows the transition density of 10,12Be →10,12C calculated with the shell model in the pand psd-shell configurations spaces, compared to VMC calculations as obtained in reference [27]. On the other hand, Figure 17 aims to reproduce the shell model results of Figure 16 by using the Suhonen and Wang Woods-Saxon parametrization sets described in subsection 3.2.1. Figure 16 shows that the WS wavefunctions provide a better description at large distances, as we previously checked in section 3.2. The normalization of C(r) for the A= 10 transition is R∞ 0C(r)dr = 1 due to the parent and daughter nuclei being mirror nuclei (10Be and 10C), while being R∞ 0C(r)dr = 0 for the A= 12 transition. However, the normalizations we obtained in our first WS results ranged around 0.90-0.98 for A= 10 and 10−3-10−2for A= 12, making us realize our calculation wasn’t correct. This made us acknowledge the isospin symmetry breaking introduced by the oscillator amplitudes A(nlj) ν. Accordingly, we upgraded our code in order to distinguish between protons and neutrons, and modified the effective hamiltonian files to substitute the general nucleon-nucleon interactions for proton-proton, neutron-neutron and proton-neutron interactions. This way, the normalizations correctly converged to 1 for A= 10 and 0 for A= 12. An important thing to note is that, for A= 12, there is a node around r∼5 fm that appears for p-shell calculations and does not for psd and VMC, highlighting the impact that different model spaces can have for the calculation. Figure 18: Normalization density C(r) for A= 48. The color code is the same as Figures 16 and 17, but for the pf-shell instead. 36 Even though our results are qualitatively similar to Wang et al.’s, there exist a couple differences regarding our Wang WS parametrization: the second peak for A= 10 appears to be smoothed out, and the positive peak of the p-shell calculation for A= 12 seems to be slightly lower. We have tried to pinpoint the source of this difference without success. We also calculated the normalization density for A= 48 using the WS parametrizations described in subsection 3.2.1. Figure 18 indicates that, although the differences are not as obvious as for lighter nuclei, there still exist non-negligible discrepancies between HO and WS wavefunctions, especially for Wang’s parametrization. 3.4 Radial distribution of the NME components Just as we calculated the normalization densities in the previous section, we can also represent the radial NME distributions CGT (r), CF(r) and CT(r), which fulfill the property M0ν GT =Z∞ 0 CGT (r)dr , M0ν F=Z∞ 0 CF(r)dr , M0ν T=Z∞ 0 CT(r)dr , (81) where M0ν GT ,M0ν Fand M0ν Tare explained in section 1.1 and defined in Equations (29)-(31). As we have seen in the previous section, the extended psd-shell configuration space provides a better description than p-shell for light (A= 10,12) nuclei. Thus, we have represented the radial distributions from Equation (81) in Figures 19-24 and the integrated matrix elements in Tables 13 and 14 working in the psd-shell for A= 10,12 and the pf-shell for A= 48, using HO and WS radial wavefunctions. Figures 19 through 22 show a lower first peak in the radial distribution of the NMEs for the WS results in light nuclei (except for the tensor component in the A= 10 case, which presents a higher peak for WS calculations). As we can see in Table 13, this results in lower values of the total NME M0ν, which implies corrections ranging between ∼13%- 16% for A= 10 and ∼7%-15% for A= 12, depending on the chosen WS parametrization. On the other hand, for heavier nuclei (A= 48), Figures 23 and 24 show that the Suhonen WS parametrization gives off radial distributions very similar to the HO case, while the Wang parametrization presents slightly more prominent differences to both of them. Table 14 shows us this true for M0ν Fand M0ν T, presenting differences below 1% between Suhonen WS and HO results. However, even while having a visually similar profile, this difference goes up to more than 4% for M0ν GT. This is due to different cancellations happening beyond the first peak: as we can see in Figure 23, in the second peak HO results present more negative values than Suhonen WS’s. If we zoom in at the tail of the profile (r > 4fm), we can distinguish how CGT (r) still presents a negative value until r∼7 fm for the HO results, while Suhonen WS’s CGT (r) stays positive beyond r∼4.5 fm. Both of this differences add up, resulting in a more important cancellation for the HO case, and finally resulting in a ∼4.5% difference for the total NME M0ν. Wang’s parametrization, on the other hand, presents roughly half that difference, at ∼2.2%. 37 Figure 19: Radial distribution CGT (r) (Equation (77)) for A= 10, calculated in the psd-shell configuration space and using different radial wavefunctions. Figure 20: Radial distributions CF(r) and CT(r) (Equation (77)) for A= 10, calculated in the psd-shell configuration space and using different radial wavefunctions. 38 Figure 21: Radial distribution CGT (r) (Equation (77)) for A= 12, calculated in the psd-shell configuration space and using different radial wavefunctions. Figure 22: Radial distributions CF(r) and CT(r) (Equation (77)) for A= 10, calculated in the psd-shell configuration space and using different radial wavefunctions. 39 Figure 23: Radial distribution CGT (r) (Equation (77)) for A= 48, calculated in the pfshell configuration space and using different radial wavefunctions. The tail is shown in detail to observe the different cancellations. Figure 24: Radial distributions CF(r) and CT(r) (Equation (77)) for A= 10, calculated in the pf-shell configuration space and using different radial wavefunctions. 40 Regarding the introduced corrections, Table 13 shows that WS radial wavefunctions reduce the total NME M0νfor A= 10,12, while Table 14 shows that they may either reduce or enhance M0νfor A= 48, depending on the chosen parametrization. In the end, the more realistic WS radial wavefunctions introduce a correction to the 0νββ decay rate of up to ∼30% for light nuclei (like 10Be and 12Be) and ∼9% for the decay of a heavier nuclei like 48Ca when compared to the decay rate obtained with HO radial wavefunctions. 10Be →10C12Be →12C M0ν GT M0ν FM0ν TM0νM0ν GT M0ν FM0ν TM0ν psd : HO 4.90 -2.11 -0.038 6.16 1.03 -0.390 -0.042 1.22 psd : WS (S) 4.26 -1.90 -0.053 5.39 0.867 -0.315 -0.024 1.04 psd : WS (W) 4.11 -1.83 -0.059 5.18 0.951 -0.321 -0.024 1.13 Table 13: M0νand its Gamow-Teller, Fermi and tensor components for the double-beta decay 10Be →10C (A= 10) and 12Be →12C (A= 12), calculated in the psd-shell configuration space and using different radial wavefunctions. 48Ca →48Ti M0ν GT M0ν FM0ν TM0ν pf : HO 0.845 -0.228 -0.058 0.929 pf : WS (S) 0.882 -0.230 -0.054 0.970 pf : WS (W) 0.825 -0.217 -0.050 0.909 Table 14: M0νand its Gamow-Teller, Fermi and tensor components for the double-beta decays 48Ca →48Ti (A= 48), calculated in the pf-shell configuration space and using different radial wavefunctions. 47 C Appendix C: Associated Laguerre polynomials Back in Equation (61), we defined the harmonic oscillator radial wavefunction gnl(r) using the associated Laguerre polynomials L(l+1 2) n(x) [30]. The first three associated Laguerre polynomials are L(l+1 2) 0= 1 ,(84) L(l+1 2) 1=l−x+3 2,(85) L(l+1 2) 2=1 2l+3 2l+5 2−2l+5 2x+x2,(86) while further polynomials can be obtained through the recursion relation L(l+1 2) n=L(l+3 2) n−L(l+3 2) n−1.(87) Instead of directly using these polynomials, numerical values for gnl(r) can be obtained through the auxiliary functions vnl(r) [39], defined by gnl(r) = s2−n+l+2(2n+ 2l+ 1)!! b3√πn![(2l+ 1)!!]2r ble−r2/2b2vnl r2 b2,(88) and then using the recursion relations vn,l−1(x) = vn−1,l−1(x)−2xvn−1,l(x) 2l+ 1 ,(89) vnl(x) = (2l+ 1)vn,l−1(x)+2nvn−1,l(x) 2n+ 2l+ 1 .(90) 48 D Appendix D: Oscillator amplitudes A(nlj) ν nlj εnlj ν= 0 ν= 1 ν= 2 ν= 3 ν= 4 ν= 5 0s1/2-17.661 0.997 -0.063 0.037 -0.022 0.004 -0.004 0p3/2-5.682 0.980 -0.154 0.110 -0.062 0.026 -0.015 0p1/2-0.806 0.940 -0.252 0.189 -0.110 0.061 -0.030 0d5/24.596 0.870 -0.336 0.279 -0.192 0.115 -0.059 1s1/23.985 0.075 0.771 -0.458 0.358 -0.228 0.106 0d3/28.151 0.616 -0.487 0.445 -0.346 0.231 -0.115 Table 15: Woods-Saxon energies εnlj and oscillator amplitudes A(nlj) νfor proton singleparticle states in 10C, calculated using the WS parameters from Eqs. (66)-(68). nlj εnlj ν= 0 ν= 1 ν= 2 ν= 3 ν= 4 ν= 5 0s1/2-19.868 0.997 -0.068 0.030 -0.021 0.003 -0.003 0p3/2-8.229 0.983 -0.146 0.094 -0.055 0.021 -0.012 0p1/2-3.977 0.960 -0.212 0.154 -0.086 0.044 -0.022 0d5/22.071 0.898 -0.309 0.244 -0.164 0.094 -0.048 1s1/21.775 0.076 0.799 -0.439 0.331 -0.208 0.095 0d3/25.818 0.687 -0.465 0.409 -0.307 0.200 -0.099 Table 16: Woods-Saxon energies εnlj and oscillator amplitudes A(nlj) νfor neutron singleparticle states in 12Be, calculated using the WS parameters from Eqs. (66)-(68). nlj εnlj ν= 0 ν= 1 ν= 2 ν= 3 ν= 4 ν= 5 0s1/2-24.512 0.999 -0.004 0.019 -0.015 -0.001 -0.003 0p3/2-11.768 0.996 -0.046 0.062 -0.033 0.007 -0.008 0p1/2-6.005 0.987 -0.104 0.111 -0.050 0.023 -0.014 0d5/20.186 0.965 -0.168 0.165 -0.096 0.049 -0.028 1s1/21.224 0.014 0.884 -0.340 0.270 -0.159 0.069 0d3/26.429 0.766 -0.410 0.370 -0.266 0.173 -0.086 Table 17: Woods-Saxon energies εnlj and oscillator amplitudes A(nlj) νfor proton singleparticle states in 12C, calculated using the WS parameters from Eqs. (66)-(68). 49 nlj εnlj ν= 0 ν= 1 ν= 2 ν= 3 ν= 4 ν= 5 0s1/2-13.934 0.988 -0.137 0.059 -0.034 0.011 -0.006 0p3/2-3.307 0.942 -0.267 0.166 -0.103 0.052 -0.026 0p1/2-0.055 0.810 -0.427 0.311 -0.212 0.128 -0.061 0d5/2-0.064 0.953 -0.197 0.191 -0.110 0.061 -0.033 1s1/2-0.064 -0.003 0.882 -0.340 0.278 -0.158 0.070 0d3/2-0.053 0.964 -0.115 0.213 -0.082 0.064 -0.027 Table 18: Woods-Saxon energies εnlj and oscillator amplitudes A(nlj) νfor neutron singleparticle states in 10Be, calculated using Wang’s parametrization (see subsection 3.2.1). nlj εnlj ν= 0 ν= 1 ν= 2 ν= 3 ν= 4 ν= 5 0s1/2-16.159 0.995 -0.082 0.041 -0.024 0.005 -0.004 0p3/2-4.551 0.972 -0.183 0.122 -0.071 0.032 -0.017 0p1/2-0.045 0.925 -0.281 0.207 -0.125 0.070 -0.034 0d5/2-0.085 0.971 -0.137 0.166 -0.087 0.046 -0.027 1s1/2-0.066 -0.040 0.922 -0.265 0.244 -0.128 0.056 0d3/2-0.044 0.474 -0.518 0.495 -0.406 0.278 -0.138 Table 19: Woods-Saxon energies εnlj and oscillator amplitudes A(nlj) νfor proton singleparticle states in 10C, calculated using Wang’s parametrization (see subsection 3.2.1). nlj εnlj ν= 0 ν= 1 ν= 2 ν= 3 ν= 4 ν= 5 0s1/2-12.429 0.982 -0.172 0.064 -0.039 0.013 -0.007 0p3/2-2.733 0.929 -0.298 0.176 -0.112 0.058 -0.028 0p1/2-0.056 0.877 -0.364 0.249 -0.163 0.093 -0.045 0d5/2-0.052 0.948 -0.217 0.188 -0.115 0.061 -0.033 1s1/2-0.051 0.022 0.877 -0.352 0.274 -0.161 0.070 0d3/2-0.066 0.498 -0.522 0.487 -0.393 0.266 -0.131 Table 20: Woods-Saxon energies εnlj and oscillator amplitudes A(nlj) νfor neutron singleparticle states in 12Be, calculated using Wang’s parametrization (see subsection 3.2.1). 50 nlj εnlj ν= 0 ν= 1 ν= 2 ν= 3 ν= 4 ν= 5 0s1/2-20.843 0.999 -0.040 0.024 -0.019 0.001 -0.003 0p3/2-8.776 0.991 -0.100 0.078 -0.045 0.014 -0.010 0p1/2-3.715 0.973 -0.167 0.136 -0.071 0.035 -0.019 0d5/2-0.043 0.968 -0.160 0.161 -0.093 0.047 -0.027 1s1/2-0.065 0.049 0.834 -0.403 0.309 -0.191 0.086 0d3/2-0.058 0.706 -0.448 0.402 -0.300 0.197 -0.098 Table 21: Woods-Saxon energies εnlj and oscillator amplitudes A(nlj) νfor proton singleparticle states in 12C, calculated using Wang’s parametrization (see subsection 3.2.1). nlj εnlj ν= 0 ν= 1 ν= 2 ν= 3 ν= 4 ν= 5 0s1/2-24.718 0.986 -0.165 -0.013 -0.010 0.006 0.002 0p3/2-17.203 0.988 -0.150 -0.005 -0.025 0.006 0.002 0p1/2-15.566 0.993 -0.117 0.004 -0.028 0.004 0.001 0d5/2-9.180 0.988 -0.148 0.027 -0.046 0.010 -0.002 1s1/2-6.641 0.157 0.952 -0.235 0.080 -0.076 0.022 0d3/2-5.906 0.988 -0.130 0.064 -0.057 0.015 -0.009 0f7/2-1.039 0.972 -0.190 0.099 -0.090 0.033 -0.017 1p3/2-0.051 0.127 0.882 -0.360 0.215 -0.160 0.064 0f5/2-0.044 0.982 -0.097 0.133 -0.085 0.033 -0.024 1p1/2-0.053 0.069 0.904 -0.319 0.220 -0.155 0.058 Table 22: Woods-Saxon energies εnlj and oscillator amplitudes A(nlj) νfor neutron singleparticle states in 48Ca, calculated using Wang’s parametrization (see subsection 3.2.1). 51 nlj εnlj ν= 0 ν= 1 ν= 2 ν= 3 ν= 4 ν= 5 0s1/2-33.426 0.995 -0.095 -0.042 -0.012 0.003 0.004 0p3/2-25.211 0.998 -0.038 -0.040 -0.025 -0.001 0.004 0p1/2-22.864 0.999 0.009 -0.025 -0.026 -0.005 0.001 0d5/2-16.228 0.999 0.013 -0.021 -0.035 -0.007 0.001 1s1/2-12.060 0.094 0.994 -0.013 -0.001 -0.046 -0.006 0d3/2-11.324 0.997 0.064 0.018 -0.031 -0.008 -0.005 0f7/2-6.745 0.998 0.042 0.016 -0.044 -0.008 -0.005 1p3/2-2.431 0.037 0.991 -0.073 0.072 -0.077 0.004 0f5/20.943 0.993 0.035 0.097 -0.050 0.010 -0.018 1p1/2-0.018 -0.009 0.984 -0.094 0.120 -0.087 0.013 Table 23: Woods-Saxon energies εnlj and oscillator amplitudes A(nlj) νfor proton singleparticle states in 48Ti, calculated using the WS parameters from Eqs. (66)-(68). nlj εnlj ν= 0 ν= 1 ν= 2 ν= 3 ν= 4 ν= 5 0s1/2-24.595 0.989 -0.147 -0.028 -0.010 0.005 0.003 0p3/2-16.980 0.993 -0.111 -0.025 -0.024 0.003 0.004 0p1/2-14.929 0.997 -0.070 -0.013 -0.027 -0.001 0.001 0d5/2-8.690 0.996 -0.084 -0.001 -0.040 0.002 0.001 1s1/2-5.238 0.142 0.975 -0.154 0.038 -0.062 0.011 0d3/2-4.510 0.997 -0.051 0.037 -0.045 0.004 -0.007 0f7/2-0.054 0.992 -0.091 0.051 -0.066 0.012 -0.009 1p3/2-0.060 0.063 0.975 -0.158 0.103 -0.098 0.019 0f5/2-0.026 0.993 0.061 0.092 -0.043 0.006 -0.017 1p1/2-0.044 -0.009 0.985 -0.093 0.120 -0.087 0.012 Table 24: Woods-Saxon energies εnlj and oscillator amplitudes A(nlj) νfor proton singleparticle states in 48Ti, calculated using Wang’s parametrization (see subsection 3.2.1). 52 E Appendix E: Talmi-Moshinsky transformation for Woods-Saxon wavefunctions in HO basis Defining the two-body WS wavefunction as we did in Equation (75), we can expand the Talmi-Moshinsky transformation for double-beta decay as follows WShn0 1l0 1j0 1p, n0 2l0 2j0 2p;J|Oα|n1l1j1n, n2l2j2n;JiWS = =X ν1ν2ν0 1ν0 2 A0 ν0 1 (n0 1l0 1j0 1)A0 ν0 2 (n0 2l0 2j0 2)A(n1l1j1) ν1A(n2l2j2) ν2HOhν0 1l0 1j0 1p, ν0 2l0 2j0 2p;J|Oα|ν1l1j1n, ν2l2j2n;JiHO = =X S,S0,λ,λ0qˆ j0 1ˆ j0 2ˆ S0ˆ λ0qˆ j1ˆ j2ˆ Sˆ λ       l0 11 2j1 l0 21 2j2 λ0S0J               l11 2j1 l21 2j2 λ S J        hl0 1l0 2λ0, S0, J|Sα|l1l2λ, S, Ji X ν1ν2ν0 1ν0 2 A0 ν0 1 (n0 1l0 1j0 1)A0 ν0 2 (n0 2l0 2j0 2)A(n1l1j1) ν1A(n2l2j2) ν2hν0 1l0 1, ν0 2l0 2|Hα(r)˜ Hα(R)|ν1l1, ν2l2i= =X S,S0,λ,λ0qˆ j0 1ˆ j0 2ˆ S0ˆ λ0qˆ j1ˆ j2ˆ Sˆ λ       l0 11 2j1 l0 21 2j2 λ0S0J               l11 2j1 l21 2j2 λ S J        X ν,ν0,l,l0,N,N0,L,L0hl0L0λ0, S0, J|Sα|lLλ, S, JiX ν1ν2ν0 1ν0 2 A0 ν0 1 (n0 1l0 1j0 1)A0 ν0 2 (n0 2l0 2j0 2)A(n1l1j1) ν1A(n2l2j2) ν2 hν0 1l0 1, ν0 2l0 2|ν0l0, N0L0iλ0hν1l1, ν2l2|νl, NLiλhν0l0, N0L0|Hα(r)˜ Hα(R)|νl, NLi =X S,S0,λ,λ0qˆ j0 1ˆ j0 2ˆ S0ˆ λ0qˆ j1ˆ j2ˆ Sˆ λ       l0 11 2j1 l0 21 2j2 λ0S0J               l11 2j1 l21 2j2 λ S J        X ν,ν0,l,l0,N,N0,L,L0hl0L0λ0, S0, J|Sα|lLλ, S, JiX ν1ν2ν0 1ν0 2 A0 ν0 1 (n0 1l0 1j0 1)A0 ν0 2 (n0 2l0 2j0 2)A(n1l1j1) ν1A(n2l2j2) ν2 hν0 1l0 1, ν0 2l0 2|ν0l0, N0L0iλ0hν1l1, ν2l2|νl, NLiλhν0l0|Hα(r)|νlihN0L0|˜ Hα(R)|NLi =X S,S0,λ,λ0qˆ j0 1ˆ j0 2ˆ S0ˆ λ0qˆ j1ˆ j2ˆ Sˆ λ       l0 11 2j1 l0 21 2j2 λ0S0J               l11 2j1 l21 2j2 λ S J        X ν,ν0,l,l0,N,N0,L,L0hl0L0λ0, S0, J|Sα|lLλ, S, JiX ν1ν2ν0 1ν0 2 A0 ν0 1 (n0 1l0 1j0 1)A0 ν0 2 (n0 2l0 2j0 2)A(n1l1j1) ν1A(n2l2j2) ν2 hν0 1l0 1, ν0 2l0 2|ν0l0, N0L0iλ0hν1l1, ν2l2|νl, NLiλhν0l0|Hα(r)|νliδN,N0δL,L0. Since most of the times ˜ Hα(R) = 1, in the last step of the development we have considered hN0L0|˜ Hα(R)|NLi=δN,N0δL,L0. 53 References [1] L. Cardani. Neutrinoless double beta decay overview. SciPost Physics Proceedings, 02 2019. [2] C. Patrignani et al. Review of Particle Physics. Chin. Phys. C, 40(10):100001, 2016. [3] Jonathan Engel and Javier Men´endez. Status and future of nuclear matrix elements for neutrinoless double-beta decay: a review. Reports on Progress in Physics, 80(4):046301, Mar 2017. [4] J. Kotila and F. Iachello. Phase-space factors for double-βdecay. Phys. Rev. C, 85:034316, Mar 2012. [5] KATRIN Collaboration and KATRIN Collaboration. Katrin design report 2004. Technical report, Forschungszentrum, Karlsruhe, 2005. 51.54.01; LK 01. [6] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. B. Barreiro, N. Bartolo, and et al. Planck 2018 results. Astronomy Astrophysics, 641:A6, Sep 2020. [7] M. Biassoni and O. Cremonesi. Search for neutrino-less double beta decay with thermal detectors. Progress in Particle and Nuclear Physics, 114:103803, 2020. [8] Sacha Davidson, Enrico Nardi, and Yosef Nir. Leptogenesis. Physics Reports, 466(4):105 – 177, 2008. [9] Vincenzo Cirigliano, Wouter Dekens, Emanuele Mereghetti, and Andr´e Walker-Loud. Neutrinoless double-βdecay in effective field theory: The light-majorana neutrinoexchange mechanism. Phys. Rev. C, 97:065501, Jun 2018. [10] A. Gando, Y. Gando, T. Hachiya, A. Hayashi, S. Hayashida, H. Ikeda, K. Inoue, K. Ishidoshiro, Y. Karino, M. Koga, S. Matsuda, T. Mitsui, K. Nakamura, S. Obara, T. Oura, H. Ozaki, I. Shimizu, Y. Shirahata, J. Shirai, A. Suzuki, T. Takai, K. Tamae, Y. Teraoka, K. Ueshima, H. Watanabe, A. Kozlov, Y. Takemoto, S. Yoshida, K. Fushimi, T. I. Banks, B. E. Berger, B. K. Fujikawa, T. O’Donnell, L. A. Winslow, Y. Efremenko, H. J. Karwowski, D. M. Markoff, W. Tornow, J. A. Detwiler, S. Enomoto, and M. P. Decowski. Search for majorana neutrinos near the inverted mass hierarchy region with kamland-zen. Phys. Rev. Lett., 117:082503, Aug 2016. [11] Maria Goeppert Mayer. On closed shells in nuclei. ii. Phys. Rev., 75:1969–1970, Jun 1949. [12] Otto Haxel, J. Hans D. Jensen, and Hans E. Suess. On the ”magic numbers” in nuclear structure. Phys. Rev., 75:1766–1766, Jun 1949. [13] Etienne Caurier. Shell model and nuclear structure. Progress in Particle and Nuclear Physics, 59(1):226 – 242, 2007. International Workshop on Nuclear Physics 28th Course. 54 [14] Alfredo Poves and Frederic Nowacki. The nuclear shell model, pages 70–101. Springer Berlin Heidelberg, Berlin, Heidelberg, 2001. [15] E. Caurier, G. Mart´ınez-Pinedo, F. Nowacki, A. Poves, and A. P. Zuker. The shell model as a unified view of nuclear structure. Rev. Mod. Phys., 77:427–488, Jun 2005. [16] Present status of shell model techniques. [17] Main website of the shell model code ANTOINE. http://www.iphc.cnrs.fr/ nutheo/code_antoine/menu.html. [18] Masaru Doi, Tsuneyuki Kotani, and Eiichi Takasugi. Double Beta Decay and Majorana Neutrino. Progress of Theoretical Physics Supplement, 83:1–175, 03 1985. [19] Mattias Blennow, Enrique Fernandez-Martinez, Jacobo Lopez-Pavon, and Javier Men´endez. Neutrinoless double beta decay in seesaw models. Journal of High Energy Physics, 2010(7), Jul 2010. [20] Meng Wang, G. Audi, F. G. Kondev, W.J. Huang, S. Naimi, and Xing Xu. The AME2016 atomic mass evaluation (II). tables, graphs and references. Chinese Physics C, 41(3):030003, mar 2017. [21] Website of the Atomic Mass Data Center. http://amdc.impcas.ac.cn/. [22] Evaluated Nuclear Structure Data File of the National Nuclear Data Center. http: //amdc.impcas.ac.cn/. [23] Matthew Redshaw, Georg Bollen, Maxime Brodeur, Scott Bustabad, David L. Lincoln, Samuel J. Novario, Ryan Ringle, and Stefan Schwarz. Atomic mass and doublebeta-decay Q value of Ca-48. Phys. Rev. C, 86:041306, 2012. [24] M. Redshaw, E. Wingfield, J. McDaniel, and E.G. Myers. Mass and double-betadecay Q value of Xe-136. Phys. Rev. Lett., 98:053003, 2007. [25] E. Caurier, F. Nowacki, and A. Poves. Shell Model description of the ββ decay of 136Xe. Physics Letters B, 711(1):62–64, May 2012. [26] Fedor ˇ Simkovic, Rastislav Dvornick´y, Du ˇsan ˇ Stef´anik, and Amand Faessler. Improved description of the 2νββ-decay and a possibility to determine the effective axial-vector coupling constant. Phys. Rev. C, 97:034315, Mar 2018. [27] X.B. Wang, A.C. Hayes, J. Carlson, G.X. Dong, E. Mereghetti, S. Pastore, and R.B. Wiringa. Comparison between variational monte carlo and shell model calculations of neutrinoless double beta decay matrix elements in light nuclei. Physics Letters B, 798:134974, Nov 2019. [28] Jouni Suhonen. From Nucleons to Nucleus: Concepts of Microscopic Nuclear Theory. Theoretical and Mathematical Physics. Springer, Berlin, Germany, 2007. [29] P. J Brussaard and (joint author.) Glaudemans, P. W. M. Shell-model applications in nuclear spectroscopy. Amsterdam ; New York : North-Holland Pub. Co. ; New York : sole distributors for the USA and Canada, Elsevier/North-Holland, 1977. 55 [30] F. Oberhettinger W. Magnus and R. P. Soni. Formulas and theorems for the special functions of mathematical physics. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift f¨ur Angewandte Mathematik und Mechanik, 47(8):554–554, 1967. [31] J. Blomqvist and A. Molinari. Collective 0−vibrations in even spherical nuclei with tensor forces. Nuclear Physics A, 106(3):545 – 569, 1968. [32] E. K. Warburton and B. A. Brown. Effective interactions for the 0p1s0d nuclear shell-model space. Phys. Rev. C, 46:923–944, Sep 1992. [33] E.K. Warburton, B.A. Brown, and D.J. Millener. Large-basis shell-model treatment of a = 16. Physics Letters B, 293(1):7 – 12, 1992. [34] D.H. Gloeckner and R.D. Lawson. Spurious center-of-mass motion. Physics Letters B, 53(4):313 – 318, 1974. [35] Alessandro Lovato. Private communication. [36] Kenneth M. Nollett and R. B. Wiringa. Asymptotic normalization coefficients from ”ab initio” calculations. Physical Review C, 83(4), Apr 2011. [37] I. Talmi. Helv. Phys, Acta 25:185, 1952. [38] Marcos Moshinsky. Transformation brackets for harmonic oscillator functions. Nuclear Physics, 13(1):104 – 116, 1959. [39] Amos de Shalit, Igal Talmi, and H S W Massey. Nuclear shell theory. Pure Appl. Phys. Academic Press, New York, NY, 1963.