Ab initio description of monopole resonances in light- and medium-mass nuclei: III. Moments evaluation in ab initio PGCM calculations

The paper is the third of a series dedicated to the ab initio description of monopole giant resonances in mid-mass closed- and open-shell nuclei via the so-called projected generator coordinate method. The present focus is on the computation of the moments $m_k$ of the monopole strength distribution, which are used to quantify its centroid energy and dispersion. First, the capacity to compute low-order moments via two different methods is developed and benchmarked for the $m_1$ moment. Second, the impact of the angular momentum projection on the centroid energy and dispersion of the monopole strength is analysed before comparing the results to those obtained from consistent quasi-particle random phase approximation calculations. Next, the so-called energy weighted sum rule (EWSR) is investigated. First, the appropriate ESWR in the center-of-mass frame is derived analytically. Second, the exhaustion of the intrinsic EWSR is tested in order to quantify the (unwanted) local-gauge symmetry breaking of the presently employed chiral effective field theory ($\chi$EFT) interactions. Finally, the infinite nuclear matter incompressibility associated with the employed $\chi$EFT interactions is extracted by extrapolating the finite-nucleus incompressibility computed from the monopole centroid energy.


Introduction
The study of giant resonances (GRs) provides valuable insights into the structural and dynamical properties of atomic nuclei.In particular, the characteristics of the isoscalar giant monopole resonance (ISGMR or GMR for brevity here) and of the isovector giant dipole resonance (IVGDR) not only deepen our comprehension of nuclear structure but also have implications for the modelisation of several astrophysical systems.This is the case, for instance, of the description of core-collapse supernovae explosions and neutron stars mergers, both phenomena being associated to the nucleosynthesis of heavy elements and the behavior of nuclear matter under extreme conditions.This article is the third (Paper III) of a series of four addressing the properties of the GMR in closed-and open-shell nuclei from an ab initio standpoint using the so-called projected generator coordinate method (PGCM).While the first paper (Paper I) [1] detailed the uncertainty budget associated to PGCM calculations of monopole and quadrupole responses, the second paper (Paper II) [2] focused on the GMR properties of 16 O, 24 Mg, 28 Si and 46 Ti.Two-dimensional PGCM calculations were shown to account well for the fragmented monopole response of (rather) light doubly open-shell nuclei thanks to their capacity (i) to capture the impact of the intrinsic static quadrupole deformation and of its fluctuations on the position of the breathing mode (typically at play in spherical nuclei), (ii) to describe in a refined way the coupling between the GMR and the giant quadrupole resonance (GQR) mechanism responsible for the appearance of an additional component in the GMR of intrinsically-deformed nuclei and (iii) to seize anharmonic effects that were shown to be significant in light systems.
The present paper focuses on the computation of the moments m k of the monopole strength distribution in order to quantify its main characteristics such as its centroid energy and dispersion.Furthermore, the first mo-ment m 1 leads to the so-called energy-weighted sum rule (EWSR) that is used to extract experimental strength functions.Also, the inverse-energy weighted sum rule (IEWSR) associated with the moment m −1 delivers, when applied to the dipole response, the so-called dipole polarizability that is relevant to the computation of radiative capture cross sections.Finally, the centroid energy of the monopole strength distribution gives access to the nucleus-dependent nuclear compressibility K A that can eventually be linked to the nuclear matter incompressibility K ∞ .The latter quantity is a key characteristic of the nuclear equation of state (EOS) and, as such, has a clear interest for several astrophysical applications.
The moments of a strength function can be computed in two ways.The first one involves an explicit sum over excited states and matrix elements of the simple one-body excitation operator F .The second one does involve the expectation value of a complicated manybody operator, but in the sole ground state.The first approach is presently denoted as the sum over excited states (SOES) method whereas the second one is referred to as the ground-state expectation value (GSEV) method.For a given many-body method, the agreement between the two approaches constitutes an internal-consistency test to pass1 in order to correctly describe the excitation mode defined by the operator F .
In this context, the formal capacity to compute low-order moments via the GSEV approach is developed in Sec. 2 and Appendix D. Based on such an advancement, and after briefly introducing the numerical setting in Sec. 3, the SOES and GSEV approaches to m 1 are compared in Sec. 4 using the PGCM monopole responses of 16 O, 24 Mg, 28 Si and 46 Ti.Sec. 5 further discusses the impact of angular-momentum projection in the SOES approach.Next, PGCM moments are compared in Sec.6 to those obtained via the quasi-particle random phase approximation (QRPA), the goal being to complement the study of Paper II dedicated to the monopole strength function by focusing on its global characteristics.Section 7 focuses on the EWSR.First, it is demonstrated that the textbook expression of the EWSR must be corrected for the fact that nuclear excitations of interest are intrinsic excitations in the center-of-mass frame.Second, the exhaustion of the intrinsic EWSR is tested in order to quantify the (unwanted) local-gauge symmetry breaking of the presently employed chiral Hamiltonian.Eventually, Sec. 8 is dedicated to accessing K A in 16 O, 24 Mg, 28 Si and 46 Ti.The computed values are then employed to extract K ∞ and verify if the result thus obtained is consistent with empirical expectations.This constitutes an important test for the chiral Hamiltonian under present use.The main findings of this work are summarised in Sec. 9 whereas a set of technical appendices complement the main body of the paper.

Moments of the strength function
In many respects, the present section follows Ref. [3].By convention, all operators at play are redefined in such a way that their expectation value in the ground state is subtracted, i.e. for a given operator Q its rescaled companion is introduced as
The k -th moment of the strength distribution 3 associated with the operator F is defined as

Mean energy(ies) and dispersion
Two sets of quantities having the dimension of an energy are introduced according to 2 Notice that Eq. ( 1) is no longer relevant as soon as one deals with commutators, since it is immediate to show that 3 From a mathematical standpoint, the moment m k constitutes the k-th moment of a discretized probability distribution associated with the transition generated by F .Moments of a physical strength function are not guaranteed to be finite.The fact that it is indeed the case or not depends on mathematical characteristics of the Hamiltonian, e.g. of inter-nucleon interactions [3]. 4 Except for m 0 , and as it will become evident below, there is in fact no difference in using F or F in Eq. (8).
They coincide for all k's if the strength distribution is concentrated in a single peak.The degree to which they differ reflects the fragmentation of the distribution.By definition, the average value of the energy distribution is given by In this work the following energy averages are also employed Compared to the centroid energy Ē1 , the scaled (constrained) energy Ẽ3 ( Ẽ1 ) is more sensitive to the high (low) energy part of the strength.
As shown in Appendix A, the moments entertain the set of inequalities providing a practical tool to set boundaries on a specific moment in case it cannot be easily computed 5 .Thanks to these inequalities, the variance of the strength distribution is shown to satisfy

SOES formulation
Inserting Eq. ( 2) into Eq.( 3) delivers the expression requiring the knowledge of excited states of the system.Equation ( 8) constitutes the SOES approach to the moments computation. 5From a practical standpoint, Eq. ( 6) holds if the involved moments are all computed within the same approximation scheme.

GSEV formulation
By means of the identity resolution on the A-body Hilbert space Eq. ( 8) can be rewritten as a ground-state expectation value Computing moments via Eq.( 10) constitutes the GSEV method based on the expectation value of a complicated operator in the sole ground state.
Clearly, the complexity of the operator at play in Eq. ( 10) increases with |k|.For k ≥ 0 the many-body rank increases with k whereas for k < 0 it further involves a non-trivial inversion.

Moment operators
Positive moments can be re-expressed in more convenient forms by invoking the appropriate definition of moment operators.As shown in Appendix B moments with k ≥ 0 can be further rewritten as with and where i and j are any pair of integers fulfilling i + j = k.By definition C 0 ≡ F .
For odd moments, Eq. ( 11) can be further expressed in terms of a commutator The last step provides a useful simplification to the structure of the operator whose ground-state expectation value is to be computed.Indeed, taking F to be a one-body operator, while the product C i C j contains up to [(n−1)k+2]-body operators, n being the highest-rank component of H, the commutator contains only up to [(n − 1)k + 1]-body operators.Because even moments can only be written in terms of anti-commutators that have the same many-body rank as the product C i C j , this simplification does not occur in this case.
Eventually, two sets of moment operators are introduced according to whose expectation value in |Ψ σ0 0 ⟩ delivers m k .Based on a Hamiltonian H containing up to three-body operators, the algebraic expressions of the tensors defining M 1 (1, 0) are explicitly derived in Appendix D .The result is used to numerically compute the PGCM m 1 moment associated with the monopole operator F = r 2 via the GSEV approach in Sec.4.2.

Alternative formulation
It is possible to access the operator M k (j + 1, j) associated with the odd positive moment m k in an alternative way.To do so, the similarity-transformed Hamiltonian is introduced, where the expansion in powers of the parameter η results from the application of Baker-Campbell-Hausdorff's identity.To match the expression given in Eq. (14b) one takes i = j + 1 and such that Based on a Hamiltonian H containing up to three-body operators, the algebraic expressions of the tensors defining M 1 (1, 0) are also derived in Appendix D as a way to validate the correctness of the expressions obtained via the more direct commutator approach laid down in Sec.2.5.

Practical merits and limitations
The great practical advantage of the GSEV approach is to access strength function's moments based on the sole knowledge of the nuclear ground state.This indeed is a tremendous simplification given that accessing a complete-enough set of excited states constitutes a challenge within any state-of-the-art ab initio many-body method 6Such a benefit however comes at the price of evaluating the ground-state expectation value of operators (Eq.( 14)) whose many-body complexity increases with the moment order.The set of moment operators indeed involve the hierarchy of operators whose many-body rank increases with j due to the new commutator involved at each step.With F a one-body operator and H containing up to three-body operators, C 1 contains up to three-body operators, C 2 up to fivebody operators, C 3 up to seven-body operators, i.e.C j contains up to (2j + 1))-body operators.As a result, Mk (i, j) contains up to (2k + 2)-body operators and M k (i, j) up to (2k + 1) operators.For example, M0 (0, 0) contains up to two-body operators and M 1 (1, 0) contains up to three-body operators.Knowing that dealing with three-body operators constitutes the current computational limit, it makes possible to compute both m 0 and m 1 exactly via the GSEV approach in PGCM calculations 7 .Moving further, M2 (1, 1) contains up to six-body operators and M 3 (2, 1) contains up to seven-body operators, which makes them beyond reach 8 .While it is in principle possible to design approximations to M2 (1, 1) and M 3 (2, 1) based on rank-reduction techniques [4], this avenue is not pursued in the present work and PGCM moments such as m 2 and m 3 are accessed via the SOES approach.
2.8 Pseudo-GSEV approach to m −1 The GSEV delivers an alternative strategy to compute moments based on the introduction of moment operators.However, this strategy does not apply to m k with k ≤ 0.
The way to evaluate m −1 , which in the SOES approach reads as via a pseudo-GSEV approach relies on time-independent perturbation theory.Perturbing the system by the external field F , the Hamiltonian becomes and the associated Schrödinger equation for the ground state is9 In the small-λ limit, perturbation theory allows one to expand |Ψ 0 (λ)⟩ in powers of λ according to [5] |Ψ 0 (λ The variation of the ground-state expectation values of a generic operator Q and of the Hamiltonian H read as where the term linear in λ disappears in Eq. (22b) due to the fact that {|Ψ σ ν ⟩, ν = 0, . . ., ν max } constitutes an orthonormal eigenbasis of H.It is easy to see that both expressions provide a direct link to Notice that the first contribution to the variation of the ground-state energy is of order λ 2 , which makes in general Eq.(23a) a more reliable numerical option to compute m −1 .While this pseudo-GSEV approach can be rather straightforwardly implemented within PGCM calculations, it is postponed to a later work such that numerical values of m −1 presented below actually rely on the SOES approach (Eq.( 18)).

Numerical setting
The PGCM formalism and the characteristics of the numerical applications were detailed in Paper I. In particular, the definition of the three mean-square-radiuslike operators r 2 , r 2 lab and r 2 int under use all throughout the present paper can be found in Paper I. All calculations presented here use the same setting as in Paper II.
A one-body spherical harmonic oscillator basis characterised by the optimal frequency ℏω = 12 MeV is employed.All states up to e max ≡ max(2n + l) = 10 are included, with n the principal quantum number and l the orbital angular momentum.The representation of three-body operators is further restricted by only employing three-body states up to e 3max = 14.
A Hamiltonian based on chiral effective field theory (χEFT) and built at next-to-next-to-next-to-leadingorder (N 3 LO) [6] is employed.It contains consistent two-(2N) and three-nucleon (3N) interactions and is further evolved via similarity renormalization group (SRG) transformations [7] to the low-momentum scale λ = 1.88 fm −1 (i.e.flow parameter α=0.08 fm 4 ) and truncated at the three-body operator level.The resulting three-body force is approximated via the rank-reduction method developed in Ref. [4].Two-dimensional (2D) PGCM calculations mix a set of constrained HFB states with axial symmetry using the root-mean-square radius r ≡ ⟨r 2 lab ⟩ and the axial mass quadrupole deformation parameter β 2 as generator coordinates.The QRPA calculations are performed at the HFB minimum via the quasi-particle finite amplitude method (QFAM) [8].The QFAM monopole moments are computed via the contour integration of the response function in the complex energy plane [9].
Eventually, the present analysis is based on the (P)GCM and QFAM monopoles responses of 16 O, 24 Mg, 28 Si and 46 Ti.Given that the present paper only focuses on spectral moments, the reader is referred to Paper II for a detailed analysis of the corresponding strength functions.4 SOES and GSEV approaches to m 1

Rationale
From a formal standpoint, the equivalence between the SOES and GSEV approaches relies on the completeness assumption from Eq. ( 9) allowing one to use the identity resolution.While the GSEV value of a given moment can be considered to be the formal value of reference, the SOES value is the one corresponding to the strength function actually computed in practice on the basis of a necessarily incomplete set of excited states 10 .In this context, the agreement between the two values constitutes an internal-consistency test for the employed many-body method relative to the excitation operator F of interest.The agreement tests whether the vector F|Ψ σ 0 ⟩ belongs to the subspace S spanned by the set of computed eigenstates {|Ψ σ ν ⟩, ν = 0, . . ., ν max } explicitly at play in the SOES approach.
Presently employed (P)GCM eigenstates are linear combinations of the non-orthogonal (projected) Bogoliubov vacua defining the set {(P σ )|Φ(r 2 , β 2 )⟩} with r ∈ [r min , r max ] and β 2 ∈ [β min , β max ] (σ ∈ IRREPs).Consequently, the (P)GCM subspace S (σ) is nothing but the span of that set.In Appendix C of Ref. [10], the monopole and quadrupole operators were shown to be indeed exhausted for a GCM calculation based on Slater determinants built out of the lowest eigenstates of axially deformed harmonic oscillators, the two generator coordinates being the corresponding axial and perpendicular oscillator frequencies.While realistic (P)GCM calculations rely on more general Bogoliubov vacua (and include particlenumber and angular-momentum projections), such a proof gives some confidence that the monopole operator might be well exhausted in present 2D (P)GCM calculations using r 2 and β 2 as generator coordinates.It is the goal of the present section to test quantitatively to which extent this is indeed the case for m 1 .

Results
The (P)GCM m 1 values obtained from both evaluation methods are reported in Tab. 1. Furthermore, their difference (rescaled according to their expected A 5/3 scaling; see Eq. ( 25)) is displayed in Fig. 1 along with the difference in percentage.
Results obtained via the SOES approach are about 6 − 7% smaller than their GSEV counterpart across the five cases under consideration.The underestimation of the SOES approach is stable from A = 16 to A = 46 once the A 5/3 scaling has been removed.The small but systematic improvement of the PGCM over the GCM is attributed to the benefit of the symmetry restoration, i.e. symmetry contaminants are removed by the angular momentum projection (AMP) on J = 0 such that the operator r 2 is better exhausted by the corresponding subspace S P .
Eventually, the operator r 2 is exhausted, within a few percents, by the (P)GCM subspace S (P ) .This translates

18976050
Table 2 Monopole moments computed using the SOES approach for GCM and PGCM calculations of 16 O, 24 Mg and 46 Ti.Numbers in between GCM and PGCM results indicate the variation between the former and the latter in percentage.into the fact that the SOES approach to m 1 can be safely used within a few percent uncertainty 11 .

Angular-momentum projection
The effect of AMP on the monopole moments m k , k = −1, 0, 1, 2, 3, evaluated via the SOES approach is presently quantified by comparing results from GCM and PGCM calculations.As seen in Tab. 2, the AMP systematically enlarges m k in a way that increases with k.In fact, while the increase with the moment order is rather marked in 16 O, it is limited in 24 Mg and has entirely disappeared in 46 Ti.Thus, and even though the range of nuclei presently tested is too limited to draw general conclusions, the impact of the AMP seems to decrease with A.
While the behavior of specific moments is interesting, it is more pertinent to investigate how this translates into the modification of physically-relevant quantities, e.g. the mean value and the dispersion of the monopole strength function.As visible from Tab. 3, the impact of the AMP on the centroid energy Ē1 decreases from 3.4% in 16 O to 1.5% in 24 Mg, and eventually down to 1.1% 11 The resulting uncertainty for a moment m k can be conjectured to increase with k.Indeed, the energy weight E k entering m k accentuates the importance of higher-energy states as k increases while the truncation of the completeness relation in the SOES approach probably affects more this higher-energy domain.Given that m 1 is the highest moment that can be computed exactly within the GSEV approach, this conjecture cannot be presently tested.
in 46 Ti.Except in 16 O, where it amounts to 750 keV, the GMR energy shift due to AMP is thus essentially negligible.The situation is similar for Ẽ1 that is used as an alternative to evaluate the GMR energy.
The impact on the dispersion is typically more significant.Again, the set of nuclei is too limited to draw general conclusions.Still, the dispersion varies by as much as 19.3% in 16 O and 7.1% in 46 Ti.In 24 Mg, the impact of AMP on σ is small but the strongly fragmented monopole strength is in fact significantly modified as can be seen in Paper II, which reflects the fact that Ē1 and σ are anyway insufficient to characterize the behavior of the strength in such a case.13].This result demonstrates the internal consistency of the QRPA as far as strength functions are concerned.While the GCM does not strictly share this property as discussed above, the GCM ground-state is necessarily a better approximation of the exact ground state than the HFB state, such that GCM moments based on the GSEV approach are necessarily better than QRPA ones 12 .This is testified by the larger values of the GCM m 1 moment reflecting the beneficial impact of (static) correlations associated with fluctuations of r 2 and β 2 leading to slightly larger GCM mean-square radii compared to HFB ones 13 .
The trend with A of the difference between GCM and QRPA m 1 values is better inferred from Fig. 2. Given the hypothesis at the heart of the QRPA, such a difference is expected to increase with the degree of anharmonicity of the system.As expected, and as discussed in Sec.6 of Paper II, larger systems are more harmonic than lighter ones.This is indeed consistent with the fact that the difference with GCM values decreases with A. This interpretation is further supported by Fig. 3 where the difference is shown to grow with the size of the cubic coefficient 14 a 3 extracted in Sec.6 of Paper II.
Figure 2 also displays the deviation of QFAM m 1 values from GCM results based on the SOES approach.In this case, QFAM values are systematically a few percent above GCM ones; i.e. they are located in between the two sets of GCM values.Eventually, the disagreement between QRPA and GCM is smaller than the uncertainty in the evaluation of the GCM values.Contrary to values based on the GSEV approach, GCM values obtained from the SOES approach do not converge towards QRPA as A increases and thus do not scale as expected with the harmonic character of the system.
Eventually, the centroid and the dispersion of the QRPA monopole strength function are compared to GCM val- 12 This recalls that the capacity of a method to fully exhaust the strength of a given excitation operator is not a sufficient condition to deliver a better approximation of the exact moments than a method that does not fully exhaust it. 13The argument qualitatively relates to the EWSR expressing m 1 in terms of the ground-state mean-square matter radius. 14The cubic coefficient is rescaled by A −3/2 to remove its trivial A dependence due to the use of the rms radius as the variable in the fitted function; see Paper II for details.ues based on the SOES approach in Tab. 5.The GCM centroid energy is typically 4 − 6% below the QRPA one 15 , which amounts to less than 1 MeV difference in the studied nuclei.The QRPA and GCM dispersions are also very consistent, especially in view of the remaining many-body uncertainty.
7 Energy weighted sum rule

Definition
The EWSR is a standard quantity in GRs studies and provides a good indicator of the degree of collectivity of nuclear excitations.Furthermore, the EWSR is used to extract strength functions from experimental data as briefly recalled in Appendix E. 15 While it is true for the centroid of the actually computed GCM strength function (SOES value), the formal (GSEV) value not computed here is probably higher in view of the behavior of m 1 studied above.The ESWR relies on an analytical evaluation of m 1 via Eq.( 14), i.e. using the GSEV approach.Targeting the first moment of the isoscalar monopole strength function, the similarity transformation of the Hamiltonian in Eq. ( 15) computed with the isoscalar local operator G 0 = F = r 2 is nothing but a local-gauge transformation.In case inter-nucleon interactions are local-gauge invariant, they do not contribute to the quadratic term in η in Eq. ( 15) which is nothing but the operator M 1 (1, 0) delivering m 1 .Under such an assumption only the (laboratory-frame) kinetic-energy operator contributes to M 1 (1, 0) such that m 1 is obtained analytically under the form [3,14] EWSR lab (r which constitutes the textbook EWSR formula for the isocalar monopole mode.Interestingly, EWSR lab (r 2 ) is proportional to the (laboratory-frame) ground-state mean-square matter radius.Thus, accessing it only requires the computation of that mean-square radius from the many-body method of interest, e.g. using the Bogoliubov state at the HFB minimum in QFAM calculations or the (P)GCM ground state in (P)GCM calculations.
However, it happens that nuclear excitations of interest are intrinsic excitations.Consequently, present manybody calculations employ the intrinsic Hamiltonian H containing the intrinsic kinetic-energy operator T int ≡ T lab − T cm , with the subtracted center-of-mass kineticenergy operator reading as In this context, the monopole m 1 moment reads as with ⃗ R cm the center-of-mass position vector.The derivation of the correction term from T cm is provided in Appendix F. Eventually, the two terms can be combined such that the appropriate, i.e. intrinsic, ESWR is given by and thus amounts to using the intrinsic mean-square radius rather than the laboratory-frame one.
The EWSR int (r 2 ) must in principle be fulfilled in ab initio calculations given that χEFT-based 2N and 3N interactions are meant to be local-gauge invariant, which is a necessary condition to achieve a consistent coupling to the electromagnetic field [15].However, enforcing the local-gauge invariance is not straightforward in practice.First, it cannot be exactly fulfilled if the same EFT truncation level is applied to both nuclear interactions and currents, even in the case of dimensional regularization.Second, the use of (nonlocal) cutoff regulators makes its fulfillment even more challenging [16].Eventually, existing χEFT-based 2N and 3N interactions are not strictly local-gauge invariant and it is our goal to quantify such a feature by testing the exhaustion of EWSR int (r 2 ) by the computed m 1 .
The potential breaking of the local gauge invariance can be straightforwardly formulated by schematically expressing the intrinsic Hamiltonian as where V lgi ≡ V − V lgi formally defines the departure of the nuclear interactions from their local-gauge invariant formulation.Given Eq. ( 29), the monopole m 1 moment effectively reads in practice as where δm lgi 1 quantifies the effective breaking of EWSR int (r 2 ).

EWSR lab versus EWSR int
Values of EWSR int (r 2 ) (EWSR lab (r   While PGCM values are systematically larger, the difference is eventually very small.These features reflect the behavior of the point-matter radii reported in Tab. 8 16 .
The relative difference between EWSR int and EWSR lab is plotted as a function of A in Fig. 4. Results are identical for QFAM and PGCM calculations.The A dependence of the difference relates to the −A −1 scaling driven by the center-of-mass-correction entering EWSR int , as is analytically demonstrated in Appendix F for its onebody component.Again, the trend reflects directly how 16 The difference between HFB and PGCM radii relates to the impact of so-called static correlations beyond the mean-field included into the PGCM ansatz.In general, static correlations have little impact on radii, the exceptions being light spherical nuclei in which they can non-negligibly increase radii and transitional nuclei in which they can strongly reduce them.
In nuclei displaying a sharp total energy surface around the deformed HFB minimum, as is the case here, the impact of static correlations on the mean-square radius is typically very small [17].The presently employed χEFT Hamiltonian typically delivers good radii such that the further addition of missing dynamical correlations to the PGCM ansatz [4] is expected to enlarge radii.
the difference of the mean-square point matter radii computed in the intrinsic frame and in the laboratory frame decreases with A.

Exhaustion of Sum Rules
The actual exhaustion of the values provided in Tab.6 is now tested.Introducing the deviation in percentage the exhaustion (100 + ε) of EWSR int (r 2 ) is reported in Tab. 9 and Fig. 5  Overall, the violation of EWSR int (r 2 ) due to the breaking δm lgi 1 of local-gauge invariance by the presently employed χEFT interactions is small and remains below 3% in the present calculations.Still, it manifests slightly differently depending on the (approximate) many-body method, the nucleus or the eigenstate under consideration.
To illustrate this point more transparently, ε is plotted in Fig. 6 as a function of A. The difference to EWSR int (r 2 ) evolves with A for the ground states under consideration and is systematically larger for PGCM than for QRPA.One observes that the trend with A is flatter for the PGCM and that QFAM results seem to approach PGCM ones as the mass increases.One may conjecture that this is a sign of better convergence of the PGCM groundstate.
Eventually, a thorough investigation of the violation of EWSR int (r 2 ) requires a larger set of nuclei and excited states as well as to employ an expansion many-body method at various truncation orders.Furthermore, δm lgi 1 must be studied as a function of the chiral order and for various regularizations of the employed χEFT interactions.Such a systematic study is left to a future work.

Nuclear incompressibility
The monopole breathing mode probes the compressibility of nuclear matter.Consequently, the infinite matter incompressibility modulus K ∞ has been extracted based on microscopic calculations of E GMR , typically within the frame of the nuclear energy density functional (EDF)    method [18].As a matter of fact, the extraction procedure is not unambiguous in itself.Furthermore, while originally applying it to a couple of doubly closed-shell nuclei ( 208 Pb and 92 Zr) led to consistent values of K ∞ , the more recent use of open-shell nuclei produced conflicting results.
The goal is to extract the value of K ∞ associated with χEFT-based interactions via ab initio calculations.In EDF calculations, it has become customary to extract K ∞ by computing directly the symmetric nuclear matter EOS, while checking that E GMR is well reproduced in a selected set of finite nuclei on the basis of the same EDF parameterization.Another approach, presently in use, consists of extracting K ∞ from the leptodermous expansion of the finite-nucleus compressibility modulus computed microscopically [19].While the former approach typically carries smaller uncertainties, the latter bypasses the need to compute the infinite matter EOS.
The second approach was recently employed to extract K ∞ for NNLO sat [20] and NNLO opt [21] χEFT-based Hamiltonians via symmetry-adapted no core shell model (SA-NCSM) calculations of 4 He, 16 O, 20 Ne and 40 Ca [22].The extracted result for NNLO sat (K ∞ = 297) was shown to be consistent, within the rather large extrapolation uncertainties, with the value (K ∞ = 253) based on the computation of the EOS with the same Hamiltonian.
Following the same protocol but only relying on a set of intrinsically-deformed nuclei, i.e. 24 Mg, 28 Si and 46 Ti, the compressibility modulus K ∞ associated with the N 3 LO Hamiltonian under use [6] is presently estimated based on PGCM and QRPA calculations.

Finite-nucleus compression modulus
The first step consists of accessing the finite-system compression modulus given by [19] which thus requires the ground-state mean-square matter radius and the GMR energy as inputs.In finite, especially light and deformed, nuclei the GMR strength is not concentrated into a single peak.Consequently, the choice of E GMR to be used in Eq. ( 32) is neither unique nor obvious.Specific derivations support the use of Ẽ1 or Ẽ3 whereas general arguments also motivate the use of the centroid energy Ē1 [19].In the following, all three cases are tested 17 .
Based on the GMR energies provided in Tab. 10, the set of K A values are given in Tab.11 and displayed in Fig. 7 as a function of A. The higher values of K A in QRPA than in PGCM reflects the characteristics of the GMR energies pointed out earlier on whenever computing PGCM moments via the SOES approach as presently done.The spread of K A values depending on the definition of E GMR is the manifestation that Ẽ1 ( Ẽ3 ) is more sensitive to the part of the strength located at lower (higher) energies than Ē1 .Eventually, K A can typically vary by as much as 30% in 24 Mg depending on that choice.However, this variation quickly decreases with A to reach 14% in 46 Ti.Such a trend is encouraging in view of extracting K ∞ .

Extraction of K ∞
The method to extract K ∞ is based on the leptodermous expansion of K A given by [19]
where K vol , K surf , K Coul and K sym are the volume, surface, Coulomb and symmetry contributions to the compression modulus, respectively.The parameter β characterizes the isospin asymmetry where N (Z) denotes the neutron (proton) number.Equation ( 33) is fitted based on the values of K A given in Tab.11 and K vol is interpreted as the infinite nuclear matter incompressibility K ∞ .Given that the Coulomb and symmetry terms do not significantly impact the asymptotic behaviour of K A for very large A [22], K ∞ can be obtained via a simple linear fit in the variable While the linear fits are displayed in Fig. 8, the corresponding values of K ∞ and K surf are reported in Tab. 12 along with the uncertainties associated with the fit.The extracted incompressibility is K ∞ ≈ 290.While QRPA central values are a few MeV higher than PGCM ones, they only differ by about 3.3% and 4.2% when using E GMR ≡ Ē1 and E GMR ≡ Ẽ3 , respectively.Eventually, QRPA and PGCM values are consistent within extrapolation uncertainties, which are significantly larger for QRPA than for PGCM results 18 .
Interestingly, while the hierarchy K A ( Ẽ1 ) < K A ( Ē1 ) < K A ( Ẽ3 ) is systematically valid for all computed nuclei with A ≤ 46, the trends are such that the extrapolation to very large A values leads to K ∞ being the smallest for E GMR ≡ Ẽ3 .Eventually, the nuclear matter incompressibility varies by 6.6% (7.5%) for QRPA (PGCM) between the two extreme values obtained for Ẽ3 and Ē1 .This confirms the trend observed above for K A as a function of A.
of K ∞ fall, within extrapolation uncertainties, into this region 19 .
Uncertainties of the present theoretical predictions (partially) evaluated in Paper I are not presently propagated to K ∞ .While they are not negligible, they are typically subleading compared to the uncertainties associated with the choice of E GMR and with the extrapolation based on the leptodermous expansion.While the range of masses presently used in the fit allows one to make quantitative statement, the use of (much) heavier systems in the future will help reducing the extrapolation uncertainty and ensure the stability of the fit.In any case, and as already stipulated in Ref. [22], the present work demonstrates that extrapolating the finite-nuclei compressibility modulus for a large enough set of nuclei can be complementary to the computation of the EOS in order to extract the nuclear matter compressibility.

Conclusions
The present paper focused on the ab initio computation of the monopole strength's moments m k .As a first step, the formal capacity to compute low-order moments in PGCM calculations via the ground-state expectation value of moment operators was achieved.This development was then exploited to validate the use, within a few percent uncertainty, of the approach based on the explicit sum over excited states for the first moment m 1 in 16 O, 24 Mg, 28 Si and 46 Ti.
With this at hand, the momentum projection was shown to have little impact on the centroid but to affect significantly the dispersion of the monopole strength distribution.Next, the centroid energy obtained in GCM calculations was demonstrated to be typically 4 − 6% below QRPA results, which amounts to less than 1 MeV difference in the nuclei under study.The QRPA and GCM dispersions were seen to be also very consistent. 19Values of K surf are also in qualitative agreement with systematic studies [24].
The next part of the study focused on the EWSR and first demonstrated that its textbook expression must be corrected for the fact that nuclear excitations of interest are intrinsic excitations in the center-of-mass frame.Having derived the appropriate intrinsic analytical EWSR, its exhaustion was shown to be violated on the level of 3% as a result of the (unwanted) localgauge symmetry breaking of the employed χEFT-based Hamiltonian [6].
Eventually, the finite-nucleus compressibility K A was computed in 24 Mg, 28 Si and 46 Ti in order to extract the infinite matter nuclear incompressibility K ∞ = 290(15) that happens to be consistent, within uncertainties, with empirical expectations.
is positive definite.For f (E) = E and g(E) = 1 Eq.(A.1) reads which provides the sequence of inequalities in Eq. ( 6).

Appendix B: Commutator approach
The introduction of moment operators first relies on expressing the moments in terms of the commutators C l .This step relies on rewriting Eq. ( 10) as with i + j = k and where the property H |Ψ σ0 0 ⟩ = 0 has been used.Since [H n , F] = [H n , F ] for n ∈ N, the bold notation can in fact be omitted.The needed commutators can be rewritten [25] as ) by virtue of H |Ψ σ0 0 ⟩ = 0 (even though the bold notation could be omitted in the meantime).This finishes to prove Eq. (11).

Appendix C: Second-quantized operators
Given an arbitrary orthonormal basis of the one-body Hilbert space H 1 represented by the particle annihilation and creation operators {c p , c † p }, a generic (particlenumber conserving) operator O containing up to threebody operators reads as O ≡ O [0] + O [2] + O [4] + O [6]  ≡ O 00 where where σ(P ) refers to the signature of the permutation P .The notation P (. . .| . ..) denotes a separation into the k particle-creation operators and the k particleannihilation operators such that permutations are only considered between members of the same group.
The algebraic expressions of the matrix elements defining the operator M 1 (1, 0) allowing to m 1 via the GSEV approach are presently derived.All notations are consistent with Appendix C for operators expressed in normal order with respect to the particle vacuum.
There are two equivalent ways to obtain the odd-moment operators, namely given by Eqs.(14b) and ( 17), respectively.They are explored separately below.
Appendix D.1: Similarity-transformed H Using Eq. ( 17) for k = 1, the operator is given by As shown below, the similarity transformation, F being a one-body operator, does not change the rank of the operator, such that M 1 (1, 0) has the same rank as H.
Introducing the identity operator in-between each pair of creation and/or annihilation operators under the form 1 = e ηF e −ηF , (D.12) the similarity transformation is separately performed on each creation (annihilation) operator.The elementary commutator together with Baker-Campbell-Hausdorff's formula allows one to obtain  The matrix elements of M 1 can also be derived from Eq. ( 14) for k = 1 (e.g.i = 0 and j = 1).This is achieved by applying Wick's theorem with respect to the particle vacuum |0⟩.In this case the only non-vanishing contraction at play is The commutator C 1 = [H, F ] is computed separately for the various components of H.The operator F being a one-body operator, the commutator preserves the n-body nature of the component H [H [1] , F ] ≡ ab c 11 1,ab c † a c b , (D.20b) [H [2] , F ] ≡ The derivation of the matrix elements from Eqs. (D.20) relies on the tool developed in Ref. [26].This tool allows one to compute the antisymmetrized matrix elements of the normal-ordered operator obtained via the commutator of any two normal-ordered operators.While the development was originally done with respect to a Bogoliubov vacuum |Φ HFB ⟩ and expressing normalordered operators in the associated quai-particle basis, it can be readily exploited here by simply substituting quasi-particle operators β † (β) with particle operators c † (c) and by using the particle vacuum |0⟩ instead of the Bogoliubov one.Naturally the particle formalism only needs to retain particle-number-conserving components.
Eventually, the matrix elements of the elementary commutator from Eqs. (D.20) can be expressed as The above result is exploited to readily compute the nested commutator needed to obtain the m 1 operator The extended writing of Eqs.(D.24) is provided in Sec.4.3.2 of Ref. [13] and is found to be identical to the similarity-evolved derivation from Appendix D.1.

Appendix E: Strength function extraction
The actual relation of the strength function to scattering observables is hereby briefly discussed.At first order in perturbation theory, the transition rate w 0→ν from the ground state |Ψ σ0 0 ⟩ to an excited state |Ψ σ ν ⟩ mediated by the time-independent operator F is provided by Fermi's golden rule The corresponding cross section σ 0→ν is obtained normalising the transition rate by the flux of incident particles and the number of scattering centers The total cross section is computed by summing over all possible final states ν so that it can be expressed as In practice, double-differential cross sections are experimentally measured to perform a multipole-decomposition analysis (MDA), allowing the extraction of the multipole strength distributions [18].In the MDA process, the experimental cross-sections at each angle are binned into small (typically, ≤ 1 MeV) excitation energy intervals.The laboratory angular distributions for each excitation-energy bin are then converted into the centreof-mass frame using standard Jacobian and relativistic kinematics.For each excitation energy bin, the experimental angular distributions are fitted by means of the least-square method with the linear combination of the calculated double-differential cross sections associated to different multipoles: where a L (E x ) is the m 1 sum rule fraction for the L-th component.The cross sections used for the fit procedure correspond to the 100% of m 1 for the L-th multipole at excitation energy E x calculated using the distorted-wave Born approximation (DWBA).In such calculations an optical potential is used as the scattering potential.The fractions of m 1 , a L (E x ), for various multipole components are determined by minimising χ 2 error.Eventually, the strength distributions for different multipolarities are obtained by multiplying the extracted a L (E x )'s by the strength corresponding to 100% m 1 at the given energy Traditionally, the energy-weighted sum rules m L,1 employed in the above procedure for different L's are always the textbook EWSR lab rather than the appropriate intrinsic one discussed in Sec. 7.

Appendix F: Intrinsic EWSR
The monopole EWSR from Eq. ( 25) is evaluated under the assumption that only the kinetic energy T lab from Eq. ( 24) contributes to Eq. (14b), such that (F.30) The correction to EWSR lab (r 2 ) due to the subtraction of the center-of-mass kinetic energy T cm (Eq.( 26)) from the Hamiltonian is given by δm cm 1 (r The above result demonstrates that the subtraction of T cm in H leads to replacing the laboratory-frame mean-square radius ⟨r 2 lab ⟩ by the intrinsic one ⟨r 2 int ⟩.The last line splits δm cm 1 (r 2 ) into its one-and a two-body contributions to demonstrate that the one-body part of T cm leads to a simple A-dependent renormalization of EWSR lab (r 2 ) [27].

Fig. 1
Fig. 1 Difference between monopole m 1 values obtained via the GSEV and SOES approaches in (P)GCM calculations as a function of A. Upper panel: difference in percentage.Lower panel: absolute difference multiplied by A −5/3 to remove the expected trivial A dependence (see Eq. (25)).

Fig. 2
Fig.2Difference between the GCM and the QRPA monopole m 1 values as a function of A. The GCM moment is evaluated both through the SOES and GSEV approaches.Upper panel: absolute difference in percentage.Lower panel: difference multiplied by A −5/3 to remove the expected trivial A dependence (see Eq. (25)).

Fig. 3
Fig.3Relative difference between the QFAM and the GSEV GCM monopole m 1 values as a function of the cubic coefficient a 3 (see Sec. 6 of Paper II).The factor A −3/2 is included to remove the trivial A dependence of a 3 .See text for details.

Fig. 4
Fig.4 Relative difference between the monopole EWSR int and EWSR lab as a function of A for QFAM and PGCM calculations.

Fig. 6
Fig. 6 Percent variation of the computed monopole m 1 compared to the corresponding EWSR int as a function of A for QFAM and PGCM calculations.The PGCM values are based on the GSEV approach.

Table 1
GCM and PGCM m 1 monopole moments computed via the SOES and GSEV approaches.All quantities are in fm4MeV.

Table 3
Average energies and dispersion computed using the SOES for GCM and PGCM calculations of16O,24Mg and 46 Ti.All results are expressed in MeV units.Numbers in between GCM and PGCM results indicate the variation between the former and the latter in percentage.

6
Comparison to QRPAIn Paper II, the QFAM et GCM monopole strengths of16O,24Mg, 28 Si and 46 Ti were discussed at length.
analytically that, within the quasi-boson approximation, odd-k QRPA moments computed with the GSEV and the SOES approaches are strictly identical, the state at play in the GSEV being the HFB ground state[11,12,

Table 5
Centroid energy and dispersion from QFAM and GCM(SOES) calculations of16O,24Mg and 46 Ti.All results are expressed in MeV units.Numbers in between QFAM and GCM results indicate the variation between the former and the latter in percentage.

Table 8
Ground-state expectation value of point-matter nuclear radii in the laboratory and intrinsic frames for HFB and PGCM calculations.All results are expressed in fm.

Table 9
Exhaustion of the monopole EWSR int in PGCM and QFAM calculations.

Table 10
Finite-nuclei compression modulus K A as a function ofA for PGCM and QFAM calculations.Different definitions of the average GMR energy E GMR entering Eq. (32) are used, see Eqs. (5) for the notation.Average GMR energies in MeV computed from QFAM and PGCM calculations according to Eqs. (4).

Table 11
Finite-nucleus compression modulus K A computed from QFAM and PGCM calculations.Values are categorised according to the definition of the GMR energy (see Eqs. (4)) employed to compute K A via Eq.(32).
O 00 is a number.Given that O is presently taken to be particle-number conserving, the k-body class O[2k]contains a single operator O kk characterized by the equal number k of particle-creation and annihilation operators.Such an operator is obviously in normal order with respect to the particle vacuum.