$\gamma\gamma$ decay as a probe of neutrinoless $\beta\beta$ decay nuclear matrix elements

We study double gamma ($\gamma\gamma$) decay nuclear matrix elements (NMEs) for a wide range of nuclei from titanium to xenon, and explore their relation to neutrinoless double-beta ($0\nu\beta\beta$) NMEs. To favor the comparison, we focus on double-magnetic dipole transitions in the final $\beta\beta$ nuclei, in particular the $\gamma\gamma$ decay of the double isobaric analog of the initial $\beta\beta$ state into the ground state. For the decay with equal-energy photons, our large-scale nuclear shell model results show a good linear correlation between the $\gamma\gamma$ and $0\nu\beta\beta$ NMEs. Our analysis reveals that the correlation holds for $\gamma\gamma$ transitions driven by the spin or orbital angular momentum due to the dominance of zero-coupled nucleon pairs, a feature common to $0\nu\beta\beta$ decay. Our shell-model findings point out the potential of future $\gamma\gamma$ decay measurements to constrain $0\nu\beta\beta$ NMEs, which are key to answer fundamental physics questions based on $0\nu\beta\beta$ experiments.


Introduction and main result
The observation of the decay of an atomic nucleus emitting only two electrons, neutrinoless double-beta (0νββ) decay, is the process experimentally most feasible to demonstrate that neutrinos are their own antiparticles [1]. Moreover, 0νββ decay is one of the most promising probes of physics beyond the standard model (BSM) of particle physics [2]. For instance, the observation of change in lepton number in 0νββ decay could help to explain the prevalence of matter in the universe [3,4]. Because of this unique potential, a very active program aims to detect 0νββ decay. Currently the most stringent constraints reach half-lives longer than 10 26 years [5][6][7][8][9][10][11][12], and next generation ton-scale experiments are being proposed, among others, for 76 Ge, 100 Mo, 130 Te and 136 Xe nuclei.
Since 0νββ decay changes lepton number-no antineutrinos are emitted to balance the two electrons-its decay rate depends on some unknown BSM parameter(s). In the standard scenario that 0νββ is triggered by the exchange of known neutrinos, this role is played by a combination of absolute neutrino masses and mixing matrix elements, m ββ . The decay rate also depends on a calculable phase-space factor [13,14], and quadratically on the nuclear matrix element (NME) that involves the initial and final nuclear states [15]. Thus, 0νββ NMEs are needed to interpret current experimental half-life limits and to anticipate the reach of future searches. However, typical NME calculations disagree up to a factor 3 [16][17][18][19][20][21][22][23][24][25][26][27][28], about an order of magnitude on the decay rate. Furthermore, first attempts to obtain more controlled 0νββ NMEs using ab initio techniques suggest smaller NME values than most previous studies [29][30][31].
A widely explored approach to reduce the uncertainty in 0νββ analyses is to study related nuclear observables. Nuclear structure [32] or muon capture [33] data are very useful to test nuclear models used to calculate NMEs, but they are not directly related to 0νββ decay. Likewise, ββ decay with neutrino emission shows no apparent correlation with 0νββ, in spite of both being second-order weak processes sharing initial and final states [34]. Useful insights could be gained from nuclear reactions, in the same spirit of the β decay information obtained in charge-exchange experiments [35,36]. Recent efforts include the measurement of nucleon pair transfers [37,38] and double charge-exchange reactions [39]. The good correlation found between 0νββ and double Gamow-Teller transitions [40] could in principle be exploited in double charge-exchange reactions, but the analyses are challenged by tiny cross sections [41,42] and involved reaction mechanisms [43,44].
In this Letter we study the correlation between the NMEs of 0νββ and second-order electromagnetic (EM) transitions emitting two photons (γγ). In fact, the latter were first studied in atoms by Goeppert-Mayer [45,46], and it was the extension to the weak interaction which led her to propose the ββ decay [47]. To ensure that nuclear-structure aspects are as similar as possible in the γγ and ββ sectors, we focus on EM double-magnetic dipole decays-which depend, like the 0νββ operator, on the nuclear spin. In addition, isospin symmetry assures a good correspondence between the γγ and ββ nuclear states if we consider the decay of the double isobaric analogue state (DIAS) of the initial ββ state-an excited state of the final ββ nucleus-into the final ββ state-the ground state (GS) of that nucleus. Our proposal expands the connections between first-order weak and EM transitions involving isobaric analogue states exploited in the past [48][49][50]. Figure 1 summarizes the main result of this work. We find a good linear correlation between γγ and 0νββ NMEs obtained with the nuclear shell model, valid across the nuclear chart. The upper panel presents results for decays in nineteen nuclei comprising titanium, chromium and iron isotopes with nucleon number 46 ≤ A ≤ 60. The lower panel covers twenty five nuclei comprising zinc, germanium, selenium, krypton, tellurium, xenon and barium isotopes with 72 ≤ A ≤ 136. The correlation is independent on the nuclear interaction used. Therefore, our findings call for γγ calculations with other many-body methods to test to what extent the shell-model correlation in Fig. 1 is universal or depends on the theoretical approach. Indeed the 0νββ correlation with double Gamow-Teller matrix elements is common for approaches as different as the nuclear shell model and energy-density functional theory, but it is apparently not fulfilled by the quasiparticle randomphase approximation method.
Second-order EM decays are naturally suppressed with respect to first-order ones. Nevertheless, γγ transitions have been measured between 0 + first-excited states and GSs, where single-γ decay is forbidden [51][52][53], and, recently, among general nuclear states in competition with γ transitions [54,55]. Future DIAS to GS γγ decay measurements, combined with the good linear correlation between NMEs presented in this work, show as a promising tool to give insights on 0νββ NMEs. A linear regression analysis supports this potential. Assuming a ±15% uncertainty in the branching ratio measurement as in Refs. [54,55] would lead to a relatively moderate error around ±(30% − 40%), dominated by the correlation. This would imply a clear improvement over the large spread in current M 0νββ calculations [15] if the same correlation is found to be valid for other many-body methods used to study 0νββ decay.

Electromagnetic DIAS to GS transitions
The γγ decay of a nuclear excited state is an EM process where two photons are emitted simultaneously: where N i , N f are the initial and final nuclear states with four-momenta p i and p f , respectively, and photons have four-momenta k, k and helicities λ, λ . The theoretical framework of nuclear two-photon decay is presented in detail in Refs. [52,56,57]. The nonrelativistic interaction Hamiltonian is given bŷ  where A µ (x) denotes the EM field,Ĵ µ (x) the nuclear current, andB µν (x, y) is a contact (seagull) operator which represents intermediate nuclear-state excitations not captured by the nuclear model, such as nucleon-antinucleon pairs. Perturbation theory up to second order in the photon field leads to the transition amplitudes where ε µλ (k) is the photon polarization vector. The initial (|i ), intermediate (|n ) and final (|f ) nuclear states have energies E i , E n and E f , respectively. The amplitude M (2) can be neglected for DIAS to GS transitions, in the absence of subleading two-nucleon currents, because it involves a one-nucleon operator in isospin space [52]. It is very useful to perform a multipole decomposition of the γγ amplitude, because nuclear states have good angular momentum. The expansion involves electric (E) and magnetic (M ) multipole operators with angular momentum L, denoted as X. The transition amplitude sums over multipoles, which factorize into a geometrical (phase space) factor and the generalized nuclear polarizability, P J , containing all the information on the nuclear structure and dynamics [52]: where the 6j-symbols depend on the total angular momenta of the initial, intermediate, and The reduced matrix elements of the EM multipole operators involve the photon energy: O(X) ∝ k L 0 , as well as the nucleon radial r, angular spherical harmonics Y L , orbital angular momentum l, and spin s operators.
Double EM and weak decays involve different nuclei: Z the neutron and proton number. In order to study the correlation between 0νββ and γγ NMEs, we focus on the γγ decay of the DIAS of the initial ββ state. This is an excited state, with isospin T = T z + 2, of the ββ daughter nucleus with isospin third component T z = (N −Z)/2. The DIAS γγ decay to the GS-the final ββ state-with T = T z , thus connects states with the same isospin structure as ββ decay: an initial state with isospin T i = T f + 2 with a final one with T f = T z . Since isospin symmetry holds very well in nuclei, we expect the nuclear structure aspects of DIAS to GS γγ and 0νββ transitions to be very similar. Altogether, the γγ decay involves the following positive-parity J i = J f = 0 nuclear states: with K a normalization constant and T − = A i t − i the nucleus isospin lowering operator, which only changes T z .
Angular momentum and parity conservation impose that transitions between 0 + states just involve the zeromultipole polarizability P 0 , with two EL or M L operators. In the long wave approximation k 0 |x| 1, satisfied when Q = E i − E f ∼ 1 − 10 MeV, dipole (L = 1) decays dominate. Since the nuclear spin is key for 0νββ decay, we focus on double-magnetic dipole (M 1M 1) processes, governed by the operator with µ N the nuclear magneton, and the neutron (n) and proton (p) spin and orbital g-factors: g s n = −3.826, g s p = 5.586, g l n = 0, g l p = 1. For M1M1 transitions Eq. (5) factorizes and the expression can be written in terms of a single NME M γγ (M 1M 1, ∆ε): which depends on the energy difference between the two photons, ∆ε = k 0 −k 0 , and on ε n = E n −(E i +E f )/2. The M 1 operator demands 1 + intermediate states.
In order to avoid the dependence on the photon energies, which depends on the nucleus, we require that the two photons share the transition energy: k 0 = k 0 = Q/2 and thus ∆ε = 0. In this case we can define the following nuclear matrix element M γγ (M 1M 1) = M γγ (M 1M 1, 0): In the following, we calculate the NMEs in Eq. (11). Note that the γγ measurement to constrain these NMEs requires an experimental setup for k 0 k 0 . This is needed because in transitions between the DIAS and the ground state Q and hence ∆ε can easily be of the order of 10 MeV, exceeding the value of ε n .
First, we calculate the final γγ state and the initial ββ one, which we rotate in isospin to obtain its DIAS as in Eq. (6). Next, we build a finite set of intermediate states {1 + n } with the Lanczos strength function method, taking as doorway state the isospin T n = T f + 1 projection of the isovector M 1 operator applied to the final state:

This guarantees intermediate states with correct angular momentum and isospin.
We evaluate the energy denominator ε n using experimental energies when possible [68,69]. For 48 Ti, E f , E i (DIAS) and also the energy of a T = T z + 1 state (6 + ) are known. We use the latter, together with the calculated energy difference between the 6 + and 1 + states with T = T z + 1, to fix the energy of the intermediate states E n . With this experimental input, M γγ (M 1M 1) only varies the result obtained with calculated energies by 0.2%. Isospin-breaking effects cancel in ε n to a very good approximation [70]. Therefore, in nuclei with unknown energy of the DIAS or T = T z +1 states, we use experimental data on states of the same isospin multiplet in neighboring nuclei: the ββ parent to fix E i , and the ββ intermediate nucleus-when available-for E 1 . Using these experimental energies modifies M γγ (M 1M 1) results by less than 5%.

Results
With these ingredients we evaluate Eq. (11). Figure 2 shows M γγ (M 1M 1) as a function of the excitation energy of the intermediate states, for nuclei covering the three configuration spaces: 48 Ti, 82 Se and 128 Te. The Lanczos strength function gives converged results to ∼ 1% after 50 − 100 iterations. Figure 2 illustrates that, in general, intermediate states up to ∼ 15 MeV can contribute to the double-magnetic dipole NME, and that only a few states dominate each transition. The comparison between weak and EM decays needs to take into account that while 0νββ changes N and Z by two units, they are conserved in γγ decay. This is achieved by comparing isospin-reduced NMEs or, alternatively, by including the ratio of Clebsch-Gordan coefficients dictated by the Wigner-Eckart theorem [71]: Figure 1 shows the good linear correlation between 0νββ NMEs and double-magnetic dipole NMEs obtained with bare spin and orbital g-factors. We observe essentially the same correlation when using effective g-factors that give slightly better agreement with experimental magnetic dipole moments and transitions: g s i (eff) = 0.9g s i , g l p (eff) = g l p +0.1, g l n (eff) = g l n −0.1 in the pf shell [72]; and g s i (eff) = 0.7g s i for pfg nuclei [73]. We have performed a linear regression analysis to the data leading to the correlation in Fig. 1. A fit to the function M 0νββ = a + bM γγ gives best-fit parameters a = 0.872, b = 0.459 for 46 ≤ A ≤ 60 (top panel), and a = 1.29, b = 1.11 for 72 ≤ A ≤ 136 (bottom). The best fit is shown with a solid line, while prediction bands at 90% confidence level (CL) are given in dashed lines. The correlation coefficients for the top and bottom panels are ρ = 0.83 and ρ = 0.84, respectively. The 90% CL bands could be combined with a hypothetical measurement of the γγ M1M1 decay to obtain a M 0νββ NME. A branching ratio measurement with a ±15% uncertainty [54,55] combined with the linear correlation would lead to a relatively moderate error around ±(30% − 40%).
The slope of the linear correlation between γγ and 0νββ NMEs in Fig. 1 only depends mildly on the mass number, being slightly larger in the pf shell than for pfg and sdgh nuclei. This distinct behaviour is due to the energy denominator in M γγ (M 1M 1): when only the numerator in Eq. (10) is considered,M γγ , the same linear correlation is common to all nuclei. Ultimately, this mild dependence on the energy denominator is key for the good correlation between M γγ and M 0νββ .
The small dependence on the energy denominator is illustrated by Fig. 2: the intermediate states that contribute more to M γγ (M 1M 1) lie systematically at lower energies in pf -shell nuclei, compared to A ≥ 72 systems. In fact, the ratio of average energy of the dominant states contributing to M γγ (M 1M 1) in the pf shell over the pfg−sdgh spaces matches very well the ratio of the slopes in the top and bottom panels of Fig. 1. Also, in the bottom panel of Fig. 1 heavier nuclei (Te, Xe, Ba) calculated with the GCN5082 interaction appear in the upper part of the correlation band. This is partly due to a mildly smaller energy denominator and also because of a slightly larger contribution of the orbital angular momentum component of the M 1M 1 operator. We can gain additional insights on the γγ − 0νββ correlation by decomposing the double-magnetic dipole NME into spin, orbital and interference parts. Since the energy denominator plays a relatively minor role, we focus on the changes in the numerator matrix element: M γγ =M γγ ss +M γγ ll +M γγ ls . Figure 3 shows the decomposition for the γγ decay of several nuclei. In some cases like 72 Zn, the spin part dominates. Here, sinceM γγ ss is proportional to the double Gamow-Teller operator, a very good correlation with 0νββ is expected [40]. In contrast, the or-bitalM γγ ll part dominates in 134 Xe or 136 Ba, sdgh nuclei with an l = 5 orbital. Remarkably, these nuclei follow the common trend in Fig. 1, which means that the correlation with 0νββ decay is not limited to operators driven by the nuclear spin. The interferenceM γγ ls is generally smaller, and can be of different sign to the dominant terms. In fact, Fig. 3 also shows that the spin and orbital contributions to γγ decay always have the same sign, preventing a cancellation that would blur the correlation with 0νββ decay. Figure 4 investigates further the relation between spin and orbital γγ contributions, decomposing the NMEs in terms of the two-body angular momenta J of the two nucleons involved in the transition. Analogously to 0νββ NMEs [17,18],M γγ is dominated by the contribution of J = 0 pairs, partially canceled by that of J > 0 ones. This behaviour is common toM γγ ss andM γγ ll , with a more marked cancellation in the spin part, as expected due to the spin-isospin SU(4) symmetry of the isovector spin operator [22,74]. The J = 0 dominance suggests that spin and orbital S = L = 0 pairs are the most relevant in γγ DIAS to GS transitions, implying that s 1 s 2 = (S 2 − 3/2)/2 < 0, and likewise l 1 l 2 < 0. Since the spin and orbital isovector g-factors also share sign, the hierarchy in Fig. 4 explains the absence of cancellations that leads to the γγ correlation with 0νββ decay. The shell-model 0νββ nuclear matrix elements in Fig. 1 have been obtained with axial coupling g A = 1.27. While the nuclear shell model is known to overestimate β [75] and two-neutrino ββ [76][77][78] matrix elements-a feature usually known as "g A quenching"-the need and amount of "quenching" required by shell-model 0νββ NMEs is uncertain. For instance, the larger 0νββ-decay momentum transfer may imply different sensitivity to missing nuclear correlations and two-body currents, the main aspects that cause the deficiencies in β and two-neutrino ββ calculations [15,79]-two-body currents are expected to be less relevant for 0νββ NMEs [80,81]. In contrast, ab initio calculations in light [79,82,83] and middle-mass nuclei [79,84] describe well β-decay matrix elements without additional adjustments. A comparison to the first ab initio 0νββ NMEs [29][30][31] and also to a recently proposed hybrid approach that combines ab initio short-range correlations with the nuclear shell model [85] suggests that 0νββ NMEs in the shell model may be moderately overestimated-by several tens of percent-in a relatively similar way for all ββ emitters. This would imply a similar correlation to the one presented in Fig. 1 but with a and b parameters modified according to the possible overestimation of the shell-model 0νββ NMEs.

Summary
We have observed a good linear correlation between 0νββ NMEs and γγ ones when the two photons share the energy of the decay. For our shell model calculations, the correlation holds across the nuclear chart, independently on the nuclear interaction used. While the correlation should also be tested with other leading manybody methods used to study 0νββ decay, this suggests a new avenue to reduce 0νββ NME uncertainties if doublemagnetic dipole DIAS to GS γγ transitions can be mea-sured, especially on the most relevant 0νββ nuclei. In fact, first steps in this direction are underway: Valiente-Dobón et al. [95] recently proposed a flagship experiment to determine the conditions of a future program to measure the γγ decay of the 48 Ca DIAS in 48 Ti. Even though these experiments are challenging due to the competition with single-γ, E1E1 γγ and nucleon-emission channels, their potential should not be underestimated. Next generation 0νββ experiments imply a significant investment with the promise to fully cover the inverted neutrino-mass hierarchy region [96], but current NME uncertainties may limit the reach of the proposals under discussion.