## Subspace-projected multireference covariant density functional theory

X. Zhang, 1 C. C. Wang, 1 C. R. Ding, 1 and J. M. Yao 1, ∗

1 School of Physics and Astronomy, Sun Yat-sen University, Zhuhai 519082, P.R. China (Dated: February 11, 2025)

Multireference density functional theory (MR-DFT) has been a pivotal method for studying nuclear lowlying states and neutrinoless double-beta (0 νββ ) decay. However, quantifying their theoretical uncertainties has been a significant challenge due to the computational demands. This study introduces a subspace-projected covariant density functional theory (SP-CDFT), which e ffi ciently emulates MR-CDFT calculations for nuclear low-lying states. This approach leverages the eigenvector continuation method combined with the quantumnumber projected generator coordinate method, based on a relativistic energy density functional (EDF). We apply SP-CDFT to investigate the correlations among the physical quantities of nuclear matter, nuclear lowlying spectroscopy, and the nuclear matrix elements (NMEs) of 0 νββ decay in the two heaviest candidate nuclei. Our findings reveal generally strong correlations between the NMEs of 0 νββ decay and the excitation energy of the 2 + 1 state, as well as the E 2 transition strength, although these correlations vary significantly among nuclei. This work also paves the way for refining nuclear EDF parameters using spectroscopic data.

Introduction. Understanding nuclear spectroscopy is essential not only for probing the intricacies of strongly correlated quantum many-body systems [1] but also for constraining fundamental symmetries and interactions at the nuclearenergy scale [2, 3]. Theoretical calculations of nuclear matrix elements (NMEs) are pivotal for interpreting and designing experiments aimed at exploring new physics at the highintensity frontier [4-7]. Among the various nuclear models, nuclear density functional theory (DFT) is uniquely versatile, successfully applied to study nuclear properties across a wide range of mass regions [8-12]. This includes not only groundstate bulk properties [13-16] but also nuclear low-lying spectroscopy [17, 18], nuclear Schi ff moments [19, 20], and 0 νββ decay [21-27]. Despite its success, nuclear DFT faces significant challenges due to discrepancies in predictions made by di ff erent energy density functionals (EDFs). These discrepancies contribute to considerable uncertainties, especially for neutron-rich nuclei with sparse experimental data [14] and for the equation of state of nuclear matter at densities far from the saturation point [28]. As facilities for radioactive ion beams and precision measurements continue to advance [29], there is an increasing need to rigorously quantify the theoretical uncertainties of nuclear DFT. This e ff ort is crucial not only for making meaningful comparisons with other models and available data, but also for conducting systematic studies of nuclear low-lying states approaching the driplines [30] and determining NMEs relevant to rare processes [31, 32].

NMEsof 0 νββ decay [21, 23, 38]. To date, the uncertainties in MR-DFT calculations have not been assessed due to the prohibitive computational cost of performing numerous repeated calculations with varying EDF parameters. This challenge also impedes the optimization of EDF parameters using nuclear spectroscopic data.

In the past decade, significant progress has been made in quantifying the uncertainties of DFT predictions for nuclear ground-state bulk properties [33-35], and in identifying possible correlations between nuclear-matter quantities [36] and neutron-star observables [37]. Extending DFT to nuclear spectroscopic properties is more challenging as it usually requires going beyond the mean-field approximation [38]. This is typically achieved through the implementation of quantum-number projections and generator coordinate method (PGCM). This framework is known as multireference (MR) DFT [8, 11, 12, 38, 39], which has been successfully applied to study nuclear low-lying states [11, 12, 40-43] and

In this work, we develop a subspace-projected covariant density functional theory (SP-CDFT), which merges multireference (MR)-CDFT [44] with the eigenvector continuation (EC) method [45] for nuclear low-lying states. The EC method, a specialized version of reduced basis methods [46, 47], has recently been frequently employed to emulate highfidelity calculations in nuclear physics [48-50]. However, it is challenging to emulate nuclear excited states [51, 52], especially in heavy nuclei where excitation energies are several orders of magnitude smaller than ground-state binding energy. We demonstrate that our SP-CDFT can e ff ectively reproduce MR-CDFT results for nuclear low-lying states with greatly reduced computational e ff ort. This approach allows to explore the correlations among the NMEs of 0 νββ decay, nuclear-matter properties, and nuclear low-lying states starting from a universal nuclear EDF, linking nuclear weak-decay properties to those of infinite nuclear matter and finite nuclei. Our findings indicate that nuclear-matter properties at saturation density are weakly correlated with nuclear spectroscopy and 0 νββ decay, whereas the latter two exhibit much stronger correlations, albeit varying with the specific nucleus. Using nuclear matter and spectroscopic data, we quantify the statistical uncertainties of the NMEs, which are found to be much smaller than the systematic discrepancies among di ff erent nuclear models.

Emulating MR-CDFT with SP-CDFT. Westart from a relativistic EDF composed of the standard kinetic energy τ ( r ), electromagnetic energy E em ( r ), as well as the nucleon-nucleon ( NN ) interaction energy [53],

$$\text{nce} \\ \text{all} \quad \ E [ \tau, \rho, \nabla \rho ; \mathbf C ] = \int d ^ { 3 } r \left [ \tau ( r ) + \mathcal { E } ^ { \text{em} } ( r ) + \sum _ { \ell = 1 } ^ { 9 } c _ { \ell } \mathcal { E } ^ { N N } _ { \ell } ( r ) \right ], \ \ ( 1 )$$

where ρ represents di ff erent types of densities and currents. The NN interaction energy is parameterized using di ff erent powers of ρ with nine contact terms. The low-energy constants c ℓ s are collectively denoted as C = { α , β S S , γ S , δ S , α V , γ V , δ V , α TV , δ TV } . The subscripts ( S , V ) indicate the scalar and vector types of NN interaction vertices in Minkowski space, respectively, and T for the vector in isospin space. For open-shell nuclei, the EDF (1) has an additional term from pairing correlation between nucleons. For the sake of simplicity, the strength parameters of pairing energy terms are fixed throughout this work. These model parameters are usually optimized to the ground-state bulk properties of finite nuclei, and infinite nuclear-matter properties around saturation density at the mean-field level [53, 54].

The wave function of nuclear low-lying state in the MRCDFT is constructed as a linear combination of mean-field configurations with projection onto good quantum numbers

$$| \Psi _ { \nu } ^ { J N Z } ( \mathbf C ) \rangle = \sum _ { \mathbf q } ^ { N _ { \mathbf q } } f _ { \nu } ^ { J N Z } ( \mathbf q, \mathbf C ) \, | J N Z ; \mathbf q, \mathbf C \rangle \,, \quad ( 2 ) \, \quad \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \ \, \$$

where the basis of many-body wave functions are given by

$$| J N Z ; \mathbf q, \mathbf C \rangle = \hat { P } ^ { J } _ { M K } \hat { P } ^ { N } \hat { P } ^ { Z } \left | \Phi ( \mathbf q, \mathbf C ) \right \rangle. \quad \ \ ( 3 ) \quad \ e p \varepsilon _ { \mathbf C 1 }$$

Here, ˆ P J MK and ˆ P N Z , are projection operators that select the component with angular momentum J and its z -projection K , neutron number N , and proton number Z , respectively [1]. The mean-field wave functions | Φ ( q C , ) ⟩ are generated from a self-consistent CDFT calculation based on the EDF (1), with constraints on the collective coordinates collectively labeled as q for a given set of parameters C . Following our previous studies [23, 24], we only include axially-deformed and parityconserving mean-field configurations, which are reasonable for the even-even nuclei of all concerned. In this case, we have K = 0 for all configurations and thus omit this symbol. Pairing correlation in each configuration state is treated using the Bardeen-Cooper-Schrie ff er (BCS) theory. The number of basis functions is denoted by N q . The mixing weight f JNZ ν ( q C , ) is determined with the variational principle, leading to the Hill-Wheeler-Gri ffi n equation [55, 56]. For each parameter set C , there are about N 2 q PGCM kernels,

$$\mathcal { N } ^ { C } ( \mathbf q, \mathbf q ^ { \prime } ) \ = \ \langle J N Z ; \mathbf q, C | J N Z ; \mathbf q ^ { \prime }, \mathbf C \rangle, \quad \ \ ( 4 a ) \quad \text{ios re}$$

$$\mathcal { H } ^ { \text{C} } ( \mathbf q, \mathbf q ^ { \prime } ) \ = \ \langle J N Z ; \mathbf q, \text{C} | \hat { H } ( \text{C} ) | J N Z ; \mathbf q ^ { \prime }, \text{C} \rangle \,, \ \ ( 4 \text{b} ) \quad \text{spee} \,.$$

to be determined numerically. The Hamiltonian kernels H C ( q q , ′ ) are evaluated with the generalized Wick's theorem [57]. In particular, the energy overlap is determined with the mixed-density prescription [58, 59]. The MR-DFT framework usually meets the problems of spurious divergences [60, 61] and finite steps [62, 63]. Although this issue has yet to be fully resolved [64], its impact is monitored by checking the convergent behaviors of the energies of the projected mean-field states that are ultimately used in the configuration mixing calculation for the nuclei of interest. We

FIG. 1. (Color online) The comparison of the excitation energies Ex (2 + 1 ) of the first 2 + state in 150 Nd from both MR-CDFT and SPCDFT(2 , k max) calculations with di ff erent values of k max, as a function of the scaling factors f ( c ℓ ) of the scalar and vector coupling constants α S and α V . The green shaded area indicates the ranges of these values in the samples.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00691v4-figure-68f66626-9dd6-4f1f-b310-781120bf0b86.jpg)

find that its impact on the low-lying states in MR-CDFT calculations is minor, consistent with the findings in light nuclei [42, 65].

To emulate the high-fidelity MR-CDFT calculation for nuclear low-lying states, we propose a SP-CDFT( Nt , k max) based on the EC method [45]. In this method, the wave function | Ψ JNZ k ( C ⊙ ) ⟩ of the k -th state for a target EDF E [ ρ, ∇ ρ ; C ⊙ ] is expanded in terms of the wave functions | Ψ JNZ ν ( C t ) ⟩ (called EC basis) of the first k max states by the Nt sampling EDFs,

$$| \bar { \Psi } _ { k } ^ { J N Z } ( \mathbf C _ { \odot } ) \rangle = \sum _ { \nu = 1 } ^ { k _ { \max } } \sum _ { t = 1 } ^ { N _ { t } } \bar { f } _ { k, \mathbf C _ { \odot } } ^ { J N Z } ( \nu, \mathbf C _ { t } ) \, | \Psi _ { \nu } ^ { J N Z } ( \mathbf C _ { t } ) \rangle \,, \quad ( 5 )$$

where k ∈ [1 , 2 , · · · , k max]. The dimension of the EC basis is thus N EC = Ntk max. The expansion coe ffi cient ¯ f JNZ k , C ⊙ ( ν, C t ) is determined by the following generalized eigenvalue equation,

$$\sum _ { v ^ { \prime } = 1 } ^ { k _ { \max } } \sum _ { t ^ { \prime } = 1 } ^ { N _ { t } } \left [ \mathcal { H } _ { t t ^ { \prime } } ^ { \nu \nu ^ { \prime } } ( \mathbf C _ { \odot } ) - \bar { E } _ { k } ^ { \mathbf C _ { \odot } } \mathcal { N } _ { t t ^ { \prime } } ^ { \nu, \nu ^ { \prime } } \right ] \bar { f } _ { k, \mathbf C _ { \odot } } ^ { J N Z } ( \nu ^ { \prime }, \mathbf C _ { t ^ { \prime } } ) = 0, \quad ( 6 )$$

where the norm and Hamiltonian kernels of EC for a target EDF E [ ρ, ∇ ρ ; C ⊙ ] are defined as [66]

$$\mathcal { N } _ { t t ^ { \prime } } ^ { \nu \nu ^ { \prime } } \ = \ \langle \Psi _ { \nu } ^ { J N Z } ( \mathbf C _ { t } ) | \, \Psi _ { \nu ^ { \prime } } ^ { J N Z } ( \mathbf C _ { t ^ { \prime } } \rangle, \quad \quad ( 7 \mathbf a )$$

$$\mathcal { H } _ { t t ^ { \prime } } ^ { \nu \nu ^ { \prime } } ( \mathbf C _ { \odot } ) \ = \ \langle \Psi _ { \nu } ^ { J N Z } ( \mathbf C _ { t } ) | \, \hat { H } ( \mathbf C _ { \odot } ) \, | \Psi _ { \nu ^ { \prime } } ^ { J N Z } ( \mathbf C _ { t ^ { \prime } } ) \rangle \,. \quad ( 7 b )$$

The e ffi ciency of the SP-CDFT is demonstrated in scenarios requiring numerous repeated MR-CDFT calculations. The speed-up factor increases almost linearly up to 10 4 when the number of sampling EDFs reaches 10 6 , a typical value used in the assessment of statistical uncertainty [49]. In other words, the SP-CDFT enables us to predict nuclear low-lying states for millions of EDFs within half an hour using a PC, a task that would otherwise take years with the high-fidelity MR-CDFT. See Ref. [66] for more details.

Figure 1 displays the excitation energy Ex (2 + 1 ) of the first 2 + state in Nd from both the MR-CDFT and SP150 CDFT( Nt , k max) calculations based on the covariant EDF (1) with di ff erent values of α S and α V . To facilitate the sampling of parameter sets, we introduce a scaling factor f ( c ℓ ) = c ℓ / c 0

ℓ

FIG. 2. (Color online) The correlation relation among di ff erent quantities, including the quantities Θ sat = { ρ , 0 E A E / , sym , L K , } of infinite nuclear matter at the saturation density ρ 0 , the low-lying states and NME M 0 ν of (a) Nd and (b) 150 Xe from the SP-CDFT(14 3) cal136 , culation. The red distributions are refined by the empirical values of Θ sat [66]. The ground-state energies are in GeV, excitation energies in MeV, and B E ( 2) in e b 2 2 . In the bottom row, the Pearson coe ffi cients are derived from the red distributions, while those in parentheses are based on the blue distributions. See the main text for details.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00691v4-figure-8cbf1184-eaab-47e0-9b5a-2008ba159382.jpg)

for each parameter in C , where c 0 ℓ is the value of the coupling constant c ℓ in the PC-PK1 set [54], which has shown remarkable success in describing the nuclear masses of a total of 4829 isotopes [16]. Since nuclear properties are mainly sensitive to these two coupling constants [66], we vary their values within 0.2% during the sampling procedure. We find that the SP-CDFT( Nt , k max ≥ 3) reasonably reproduces the results of MR-CDFT calculations and generally performs much better than the SP-CDFT( Nt , 1). This finding confirms the advantage of using the extended EC scheme in Ref. [67], which helps reduce the number Nt of required training samples. Of particular interest is the observation that the excitation energy in Nd 150 increases smoothly with α V and decreases with α S . Notably, these trends vary across di ff erent nuclei [66]. We also checked the global performance of the SP-CDFT for other spectroscopic quantities, including ground-state energies and proton radii, excitation energies of 2 + 1 states, and B E ( 2 : 0 + 1 → 2 ) + 1 for Nd, 150 Xe, and their daughter nuclei in 0 136 νββ decay. To this end, we sampled 14 training parameter sets and 64 testing sets in the vicinity of the PC-PK1 set using the Latin hypercube sampling method [68]. We checked the results by sampling the parameters in the vicinity of the PC-F1 [53], and obtained similar conclusions. Therefore, we mainly show the results by the samples around the PC-PK1 set in this paper. For the 64 testing sets, the relative emulator errors are generally less than 6% [66].

Moreover, the SP-CDFT allows us to perform a global sensitivity analysis of the low-lying states and the NME M 0 ν to the EDF parameters. The M 0 ν is computed using the transition operators based on the standard mechanism, see Refs. [23, 24] for details. It is important to note that the contributions from the contact transition operator [69] and two-body currents are not included in this analysis. While their omission may introduce some errors, these contributions are expected to be sub-leading for the following reasons. Recent studies have demonstrated that the renormalizability of the transition amplitude is automatically ensured at leading order within the relativistic framework [70], classifying the contact transition operator as a next-to-next-to-leading order (N2LO) e ff ect in this context. Additionally, two-body currents, which appear at next-to-next-to-next-to-leading order (N3LO), are known to reduce the matrix elements by approximately 10% [71]. Moreover, we assume that all parameters in C are independent of each other and uniformly distributed. We find that approximately 90% of the variance in the quantities of interest can be attributed to the parameters α S and α V . Secondary contributions come from β s and γ s , which account for most of the remaining variance. Interestingly, the sensitivity indices are consistent across di ff erent nuclei [66].

Application of the SP-CDFT to statistical analysis. We sampled totally about 1 3 . × 10 6 parameter sets of the EDF by varying the nine parameters C around their optimal values using quasi Monte-Carlo sampling with a uniform distribution. To explore correlations among di ff erent quantities, we compute the properties Θ sat = { ρ , 0 E A E / , sym , L K , } of infinite nuclear matter at the saturation density ρ 0 , the low-lying states and NMEs of two heaviest candidate nuclei for 0 νββ decay, where the symmetry energy and its slope are defined as E sym ≡ 1 2 ∂ 2 ( E A / ) ∂η 2 GLYPH&lt;12&gt; GLYPH&lt;12&gt; GLYPH&lt;12&gt; GLYPH&lt;12&gt; η = 0 and L ≡ 3 ρ ∂ E sym ∂ρ , respectively, with η = ( ρ ( n ) -ρ ( p ) ) /ρ . The total nucleon-number density is a summation of those for neutrons and protons ρ = ρ ( n ) + ρ ( p ) . Figure 2 shows the correlation relations among these quantities. In addition, we select out 457,380 samples that are able to reproduce the empirical data of the nuclear matter simultaneously [66], and display these results in Fig. 2 as well. The correlations of nuclear spectroscopic quantities with E sym and L are generally weak in both cases, especially for the nearspherical nucleus Xe. 136 Specifically, the excitation energy Ex (2 + 1 ) of Nd is somewhat correlated to the saturation den150 sity ρ 0 , average nucleon energy E A / and incompressibility K . However, these correlations disappear in Xe. 136 Moreover, one observes a strong anti-correlation between the proton radius Rp and ρ 0 , and a linear correlation between the groundstate energy E (0 + 1 ) and E A / , consistent with the finding in a previous study based on Skyrme EDFs [72]. This finding justifies the use of E A / and ρ 0 of nuclear matter to infer the prediction of EDFs for the binding energies and radii of finite nuclei, respectively [53, 54].

In particular, we observe that the NME M 0 ν of 0 νββ decay is generally strongly correlated with the Ex (2 + 1 ) and B E ( 2 : 0 + 1 → 2 ), + 1 even though the correlation relations are di ff er-

ent for the prolate deformed nucleus Nd and spherical nu150 cleus Xe. 136 With the refinement from the empirical values of nuclear matter at saturation density, the correlations are slightly increased, and the Pearson coe ffi cients for 150 Nd are 0 93 and . -0 93, respectively. . For 136 Xe, the Pearson coe ffi -cients become -0 91 and 0 67, respectively. . . These numbers are much larger than the values 0 65 and . -0 18, respectively, . found in the interacting shell model (ISM) [73]. In the ISM for Ca, these numbers are 0 23 and 0 43, respectively [74]. 48 . . It is worth noting that the Pearson coe ffi cients between the M 0 ν and Ex (2 + 1 ) in Ge and 76 Se from the 76 ab initio valence-space in-medium similarity renormalization group (VS-IMSRG) are found to be 0 69 and 0 79, respectively [75]. . . In short, there are strong correlations between the NME M 0 ν and the properties of nuclear low-lying states, but these correlations may vary with nuclei. The observed correlations provide a basis for employing the Bayesian model averaging method for 0 νββ studies [76].

Figure 3 shows the probability distribution of excitation energies of 2 + 1 states and transition strengths B E ( 2 : 0 + 1 → 2 ) + 1 in Nd, 150 Sm, 150 Xe, and 136 Ba. In the Bayesian analysis, 136 the priors of the parameters are taken in the Gaussian form of χ 2 which reflects the deviation from the optimal value [36]. Besides, we also use the results (including energy per particle, pressure, speed of sound, nuclear symmetry energy and its slope) of ab initio many-body perturbation theory calculations for pure-neutron and symmetric nuclear matter with density ρ &lt; 0.16 fm -3 based on the chiral NN and 3 N potentials up to N 3 LO [82] to weight the 457,380 samples through Gaussian form of χ 2 [66]. With the constraints from nuclear matter and the inclusion of emulator errors of SP-CDFT, independently, we obtain the posterior probability distributions for Ex (2 + 1 ) and B E ( 2). It is seen from Fig. 3 that the excitation energies and E 2 transition strengths of Nd and 150 Sm 150 are reasonably reproduced, while the quadrupole collectivity of Xe is overestimated and this is expected to be improved 136 by including non-collective particle-hole excitation configurations in the model space. Moreover, we find that the statistical uncertainties (92% C.L.) from the nine parameters are within 21% for the excitation energies and 12% for the E 2 transition strengths.

In addition to nuclear matter properties, we also incorporate data on both Ex (2 + 1 ) and B E ( 2) to determine the posterior distribution of M 0 ν using the Bayesian method [36], where Pearson coe ffi cients are applied to weight their relative contributions. The median, 4th percentile and 96th percentile of M 0 ν , derived from the probability distribution with and without the refinement from nuclear matter, as well as from the posterior distribution, are compared with those of other nuclear model calculations in Fig. 3(e). The results indicate that incorporating nuclear matter and nuclear low-lying state information into the Bayesian method does not significantly alter the NMEs. The median value and uncertainty (92% C.L.) for the NMEs are 4 34 . + 0 09 . -0 10 . for Xe and 5 51 136 . + 0 24 . -0 26 . for Nd, which align closely with the values of 4.32 and 150

FIG. 3. (Color online) Probability distributions for the excitation energy of the 2 + 1 state, B E ( 2 : 0 + 1 → 2 ), and the NME + 1 M 0 ν of 0 νββ decay in Nd, 150 Sm, shown without (blue) and with (red) refine150 ment of Θ sat, as well as the posterior distribution (grey) considering Θ low and emulator errors ( ϵ em). The median values and uncertainties (92% C.L.) of these distributions are displayed alongside the experimental data for comparison. The results of M 0 ν are compared to the values by the MR-CDFT based on the relativistic EDF PCPK1 [24] (REDF) and those of other model calculations, including the interacting-boson model (IBM2) [77], MR-DFT based on a nonrelativistic Gogny force (NREDF) [21], quasiparticle random-phase approximation (QRPA) based on an e ff ective interaction of G matrix [78] and a Skyrme EDF [22], quasiparticle vacua shell model (QVSM) [79], triaxial projected shell model (TPSM) [80], and ab initio calculation (VS-IMSRG) [81].

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00691v4-figure-6c2c6923-a697-45c6-a880-3a9b4babb5e4.jpg)

5.60 obtained using the PC-PK1 functional [24]. We also perform a statistical analysis that accounts for the correlations among di ff erent observables, obtaining an NME of 4 33 . + 0 09 . -0 11 . for 136 Xe and 5 52 . + 0 22 . -0 26 . for 150 Nd. For the samples around the PC-F1, these values become 4 21 . + 0 07 . -0 10 . for 136 Xe and 6 07 . + 0 30 . -0 32 . for 150 Nd. These two sets of results are consistent with each other when the uncertainties are taken into account. By using the Bayesian model averaging (BMA) method for these two popular EDFs based on their predictions on the B E ( 2) values of candidate nuclei, we obtain the NMEs of 4 34 . + 0 09 . -0 11 . for 136 Xe and 5 52 . + 0 23 . -0 26 . for 150 Nd. These results demonstrate that the statistical uncertainties in the NMEs calculated with the relativistic EDF are much smaller than the discrepancies observed among di ff erent nuclear models. It is worth noting that the statistical uncertainty shown in Fig. 3 does not include the

contribution from the pairing part [25, 83], which may further increase the overall uncertainty.

Conclusions. In summary, we have developed a SPCDFT to emulate MR-CDFT for nuclear low-lying states by integrating an extended eigenvector continuation method into the PGCM. Using SP-CDFT, we have demonstrated that the NMEs of 0 νββ decay in the two heaviest candidate nuclei are correlated with the properties of nuclear low-lying states, although these correlations vary among di ff erent nuclei. We derived the posterior distribution of NMEs by utilizing empirical values and ab initio calculations of nuclear-matter properties at and below saturation density, alongside data from nuclear low-lying spectroscopy. By using the Bayesian model averaging (BMA) method in the analysis of samples in the vicinity of two popular EDFs, we obtain the NMEs of 4 34 . + 0 09 . -0 11 . for 136 Xe and 5 52 . + 0 23 . -0 26 . for 150 Nd. These statistical uncertainties are notably smaller than the discrepancies observed among di ff erent nuclear models. It is important to acknowledge, however, that several sources of uncertainty remain, with systematic uncertainties being particularly challenging to quantify within the current framework. Future work should include additional generator coordinates and explore variations in both coupling constants and energy density functional forms. The correlations observed across di ff erent nuclear models provide a foundation for applying the BMA to determine the NMEs of 0 νββ decay [6]. This approach can also be extended to PGCM calculations based on non-relativistic EDFs and chiral Hamiltonians, o ff ering a pathway to refine nuclear EDF parameters beyond the mean-field level using nuclear spectroscopy in the future. Moreover, our method provides a tool of choice to examine the correlations between the NMEs with the observables in high-energy collisions of candidate nuclei in the future [84, 85].

Acknowledgments We thank K. Hagino, H. Hergert, W.G. Jiang, C.F. Jiao, D. Lee, G. Li, and X.L. Zhang for fruitful discussions at di ff erent stages, Y.K. Wang for carefully reading the manuscript, and C. Drischler for providing nuclearmatter properties from ab initio calculations. This work is partly supported by the National Natural Science Foundation of China (Grant Nos. 12405143, 12375119 and 12141501), the Guangdong Basic and Applied Basic Research Foundation (2023A1515010936), and the Fundamental Research Funds for the Central Universities, Sun Yat-sen University. We also acknowledge the Beijing Super Cloud Computing Center (BSCC) for providing HPC resources.

∗ Corresponding author: yaojm8@sysu.edu.cn

- [1] P. Ring and P. Schuck, The nuclear many-body problem (Springer-Verlag, New York, 1980).
- [2] J. Engel, M. J. Ramsey-Musolf, and U. van Kolck, Prog. Part. Nucl. Phys. 71 , 21 (2013), arXiv:1303.2371 [nucl-th].
- [3] X. F. Yang, S. J. Wang, S. G. Wilkins, and R. F. Garcia Ruiz, Prog. Part. Nucl. Phys. 129 , 104005 (2023), arXiv:2209.15228 [nucl-ex].
- [4] T. Chupp, P. Fierlinger, M. Ramsey-Musolf, and J. Singh, Rev. Mod. Phys. 91 , 015001 (2019), arXiv:1710.02504 [physics.atom-ph].
- [5] R. Alarcon et al. , in Snowmass 2021 (2022) arXiv:2203.08103 [hep-ph].
- [6] V. Cirigliano et al. , (2022), arXiv:2203.12169 [hep-ph].
- [7] M. Agostini, G. Benato, J. A. Detwiler, J. Menéndez, and F. Vissani, Rev. Mod. Phys. 95 , 025002 (2023), arXiv:2202.01787 [hep-ex].
- [8] M. Bender, P.-H. Heenen, and P.-G. Reinhard, Rev. Mod. Phys. 75 , 121 (2003).
- [9] D. Vretenar, A. Afanasjev, G. Lalazissis, and P. Ring, Physics Reports 409 , 101 (2005).
- [10] J. Meng, H. Toki, S. G. Zhou, S. Q. Zhang, W. H. Long, and L. S. Geng, Prog. Part. Nucl. Phys. 57 , 470 (2006), arXiv:nuclth 0508020. /
- [11] T. Nikši´ c, D. Vretenar, and P. Ring, Prog. Part. Nucl. Phys. 66 , 519 (2011).
- [12] J. A. Sheikh, J. Dobaczewski, P. Ring, L. M. Robledo, and C. Yannouleas, Journal of Physics G: Nuclear and Particle Physics 48 , 123001 (2021).
- [13] S. Goriely, S. Hilaire, M. Girod, and S. Peru, Phys. Rev. Lett. 102 , 242501 (2009).
- [14] J. Erler, N. Birge, M. Kortelainen, W. Nazarewicz, E. Olsen, A. M. Perhac, , and M. Stoitsov, Nature 486 , 509 (2012).
- [15] A. V. Afanasjev, S. E. Abgemava, D. Ray, and P. Ring, Phys. Lett. B 726 , 680 (2013), arXiv:1309.3289 [nucl-th].
- [16] P. Guo et al. (DRHBcMassTable), Atom. Data Nucl. Data Tabl. 158 , 101661 (2024), arXiv:2402.02935 [nucl-th].
- [17] M. Bender, G. F. Bertsch, and P. H. Heenen, Phys. Rev. Lett. 94 , 102503 (2005), arXiv:nucl-th 0410023. /
- [18] T. R. Rodríguez, A. Arzhanov, and G. Martínez-Pinedo, Phys. Rev. C 91 , 044315 (2015), arXiv:1407.7699 [nucl-th].
- [19] J. Engel, M. Bender, J. Dobaczewski, J. H. De Jesus, and P. Olbratowski, Phys. Rev. C 68 , 025501 (2003), arXiv:nuclth 0304075. /
- [20] J. Dobaczewski, J. Engel, M. Kortelainen, and P. Becker, Phys. Rev. Lett. 121 , 232501 (2018), arXiv:1807.09581 [nucl-th].
- [21] T. R. Rodríguez and G. Martínez-Pinedo, Phys. Rev. Lett. 105 , 252503 (2010).
- [22] M. T. Mustonen and J. Engel, Phys. Rev. C 87 , 064302 (2013).
- [23] L. S. Song, J. M. Yao, P. Ring, and J. Meng, Phys. Rev. C 90 , 054309 (2014).
- [24] J. M. Yao, L. S. Song, K. Hagino, P. Ring, and J. Meng, Phys. Rev. C 91 , 024316 (2015).
- [25] C. R. Ding, X. Zhang, J. M. Yao, P. Ring, and J. Meng, Phys. Rev. C 108 , 054304 (2023).
- [26] Y. K. Wang, P. W. Zhao, and J. Meng, Phys. Lett. B 855 , 138796 (2024), arXiv:2403.06455 [nucl-th].
- [27] Y. Wang, P. Zhao, and J. Meng, Sci. Bull. 69 , 2017 (2024), arXiv:2304.12009 [nucl-th].
- [28] B.-A. Li, L.-W. Chen, and C. M. Ko, Phys. Rept. 464 , 113 (2008), arXiv:0804.3580 [nucl-th].
- [29] G. Arrowsmith-Kron et al. , Rept. Prog. Phys. 87 , 084301 (2024), arXiv:2302.02165 [nucl-ex].
- [30] T. R. Rodriguez and J. L. Egido, Phys. Rev. Lett. 99 , 062501 (2007).
- [31] J. Engel and J. Menéndez, Rep. Prog. Phys. 80 , 046301 (2017).
- [32] A. Belley et al. , Phys. Rev. Lett. 132 , 182502 (2024), arXiv:2308.15634 [nucl-th].
- [33] J. Dobaczewski, W. Nazarewicz, and P. G. Reinhard, J. Phys. G 41 , 074001 (2014), arXiv:1402.4657 [nucl-th].
- [34] J. D. McDonnell, N. Schunck, D. Higdon, J. Sarich, S. M. Wild, and W. Nazarewicz, Phys. Rev. Lett. 114 , 122501 (2015),

- arXiv:1501.03572 [nucl-th].
- [35] S. E. Agbemava, A. V. Afanasjev, and A. Taninah, Phys. Rev. C 99 , 014318 (2019), arXiv:1811.11473 [nucl-th].
- [36] P. Giuliani, K. Godbey, E. Bonilla, F. Viens, and J. Piekarewicz, (2022), arXiv:2209.13039 [nucl-th].
- [37] M. Salinas and J. Piekarewicz, Phys. Rev. C 107 , 045802 (2023), arXiv:2301.09692 [nucl-th].
- [38] J. M. Yao, J. Meng, Y. F. Niu, and P. Ring, Prog. Part. Nucl. Phys. 126 , 103965 (2022), arXiv:2111.15543 [nucl-th].
- [39] L. M. Robledo, T. R. Rodríguez, and R. R. Rodríguez-Guzmán, J. Phys. G 46 , 013001 (2019), arXiv:1807.02518 [nucl-th].
- [40] M. Bender and P.-H. Heenen, Phys. Rev. C 78 , 024309 (2008).
- [41] T. R. Rodríguez and J. L. Egido, Phys. Rev. C 81 , 064323 (2010).
- [42] J. M. Yao, J. Meng, P. Ring, and D. Vretenar, Phys. Rev. C 81 , 044311 (2010).
- [43] B. Bally, B. Avez, M. Bender, and P. H. Heenen, Phys. Rev. Lett. 113 , 162501 (2014), arXiv:1406.5984 [nucl-th].
- [44] J. M. Yao, K. Hagino, Z. P. Li, J. Meng, and P. Ring, Phys. Rev. C 89 , 054306 (2014), arXiv:1403.4812 [nucl-th].
- [45] D. Frame, R. He, I. Ipsen, D. Lee, D. Lee, and E. Rrapaj, Phys. Rev. Lett. 121 , 032501 (2018).
- [46] A. Quarteroni and G. Rozza, Reduced Order Methods for Modeling and Computational Reduction (Springer Cham, 2013).
- [47] E. Bonilla, P. Giuliani, K. Godbey, and D. Lee, Phys. Rev. C 106 , 054322 (2022), arXiv:2203.05284 [nucl-th].
- [48] S. Yoshida, N. Shimizu, T. Togashi, and T. Otsuka, Phys. Rev. C 98 , 061301 (2018), arXiv:1810.03263 [nucl-th].
- [49] A. Ekström and G. Hagen, Phys. Rev. Lett. 123 , 252501 (2019), arXiv:1910.02922 [nucl-th].
- [50] W. G. Jiang, C. Forssén, T. Djärv, and G. Hagen, Phys. Rev. C 109 , 064314 (2024), arXiv:2212.13216 [nucl-th].
- [51] S. Yoshida and N. Shimizu, PTEP 2022 , 053D02 (2022), arXiv:2105.08256 [nucl-th].
- [52] Z. H. Sun, A. Ekström, C. Forssén, G. Hagen, G. R. Jansen, and T. Papenbrock (2024) arXiv:2404.00058 [nucl-th].
- [53] T. Burvenich, D. G. Madland, J. A. Maruhn, and P. G. Reinhard, Phys. Rev. C 65 , 044308 (2002), arXiv:nucl-th 0111012. /
- [54] P. W. Zhao, Z. P. Li, J. M. Yao, and J. Meng, Phys. Rev. C 82 , 054319 (2010).
- [55] D. L. Hill and J. A. Wheeler, Phys. Rev. 89 , 1102 (1953).
- [56] J. J. Gri ffi n and J. A. Wheeler, Phys. Rev. 108 , 311 (1957).
- [57] R. Balian and E. Brezin, Nuovo Cimento B Serie 64 , 37 (1969).
- [58] P. Bonche, J. Dobaczewski, H. Flocard, P. H. Heenen, and J. Meyer, Nucl. Phys. A 510 , 466 (1990).
- [59] A. Valor, P. H. Heenen, and P. Bonche, Nucl. Phys. A 671 , 145 (2000), arXiv:nucl-th 9910032. /
- [60] M. Anguiano, J. L. Egido, and L. M. Robledo, Nucl. Phys. A 696 , 467 (2001), arXiv:nucl-th 0105003. /
- [61] J. Dobaczewski, M. V. Stoitsov, W. Nazarewicz, and P. G. Reinhard, Phys. Rev. C 76 , 054315 (2007), arXiv:0708.0441 [nuclth].
- [62] M. Bender, T. Duguet, and D. Lacroix, Phys. Rev. C 79 , 044319 (2009), arXiv:0809.2045 [nucl-th].
- [63] T. Duguet, M. Bender, K. Bennaceur, D. Lacroix, and T. Lesinski, Phys. Rev. C 79 , 044320 (2009), arXiv:0809.2049 [nuclth].
- [64] J. A. Sheikh, J. Dobaczewski, P. Ring, L. M. Robledo, and C. Yannouleas, J. Phys. G 48 , 123001 (2021), arXiv:1901.06992 [nucl-th].
- [65] E. F. Zhou, X. Y. Wu, and J. M. Yao, Phys. Rev. C 109 , 034305 (2024), arXiv:2311.15305 [nucl-th].
- [66] X. Zhang, C. C. Wang, C. R. Ding, and J. M. Yao, (2024), arXiv:2408.13209 [nucl-th].
- [67] Q. Y. Luo, X. Zhang, L. H. Chen, and J. M. Yao, Phys. Rev. C 110 , 014309 (2024).
- [68] S. Dutta and A. H. Gandomi, in Handbook of Probabilistic Models , edited by P. Samui, D. Tien Bui, S. Chakraborty, and R. C. Deo (Butterworth-Heinemann, 2020) pp. 369-381.
- [69] V. Cirigliano, W. Dekens, J. de Vries, M. L. Graesser, E. Mereghetti, S. Pastore, and U. van Kolck, Phys. Rev. Lett. 120 , 202001 (2018).
- [70] Y. Yang and P. Zhao, Physics Letters B 855 , 138782 (2024).
- [71] J.-Q. Wang, Z.-C. Gao, Y.-J. Ma, and Y. S. Chen, Phys. Rev. C 98 , 021301 (2018).
- [72] P.-G. Reinhard and W. Nazarewicz, Phys. Rev. C 93 , 051303 (2016).
- [73] M. Horoi, A. Neacsu, and S. Stoica, Phys. Rev. C 107 , 045501 (2023), arXiv:2302.03664 [nucl-th].
- [74] M. Horoi, A. Neacsu, and S. Stoica, Phys. Rev. C 106 , 054302 (2022), arXiv:2203.10577 [nucl-th].
- [75] A. Belley, T. Miyagi, S. R. Stroberg, and J. D. Holt, in Matrix Elements for the Double-beta-decay EXperiments (2022) arXiv:2210.05809 [nucl-th].
- [76] V. Cirigliano et al. , J. Phys. G 49 , 120502 (2022), arXiv:2207.01085 [nucl-th].
- [77] J. Barea, J. Kotila, and F. Iachello, Phys. Rev. C 91 , 034304 (2015).
- [78] D.-L. Fang and A. Faessler, Phys. Rev. C 103 , 045501 (2021).
- [79] Y. Tsunoda, N. Shimizu, and T. Otsuka, Phys. Rev. C 108 , L021302 (2023), arXiv:2304.11780 [nucl-th].
- [80] Y. K. Wang, P. W. Zhao, and J. Meng, Phys. Rev. C 104 , 014320 (2021).
- [81] A. Belley, T. Miyagi, S. R. Stroberg, and J. D. Holt, (2023), arXiv:2307.15156 [nucl-th].
- [82] C. Drischler, R. J. Furnstahl, J. A. Melendez, and D. R. Phillips, Phys. Rev. Lett. 125 , 202702 (2020).
- [83] W. L. Lv, Y. F. Niu, D. L. Fang, J. M. Yao, C. L. Bai, and J. Meng, Phys. Rev. C 108 , L051304 (2023), arXiv:2302.04423 [nucl-th].
- [84] M. I. Abdulhamid et al. (STAR), Nature 635 , 67 (2024), arXiv:2401.06625 [nucl-ex].
- [85] J. Jia et al. , Nucl. Sci. Tech. 35 , 220 (2024), arXiv:2209.11042 [nucl-ex].
X. Zhang C. C. Wang C. R. Ding J. M. Yao
Subspace-projected multireference covariant density functional theory
Multireference density functional theory (MR-DFT) has been a pivotal method
for studying nuclear low-lying states and neutrinoless double-beta
($0\nu\beta\beta$) decay. However, quantifying their theoretical uncertainties
has been a significant challenge due to the computational demands. This study
introduces a subspace-projected covariant density functional theory (SP-CDFT),
which efficiently emulates MR-CDFT calculations for nuclear low-lying states.
This approach leverages the eigenvector continuation method combined with the
quantum-number projected generator coordinate method, based on a relativistic
energy density functional (EDF). We apply SP-CDFT to investigate the
correlations among the physical quantities of nuclear matter, nuclear low-lying
spectroscopy, and the nuclear matrix elements (NMEs) of $0\nu\beta\beta$ decay
in the two heaviest candidate nuclei. Our findings reveal generally strong
correlations between the NMEs of $0\nu\beta\beta$ decay and the excitation
energy of the $2_1^+$ state, as well as the $E2$ transition strength, although
these correlations vary significantly among nuclei. This work also paves the
way for refining nuclear EDF parameters using spectroscopic data.