Converged ab initio calculations of heavy nuclei

Our starting Hamiltonian in second-quantized form is

where $t_{p^{\prime}p}$, $V^{\rm NN}_{p^{\prime}q^{\prime}pq}$, and $V^{\rm 3N}_{p^{\prime}q^{\prime}r^{\prime}pqr}$ are the one-, two-, and three-body matrix elements, respectively. The index $p$ labels the single-particle orbit with quantum numbers $\{n_{p},\ell_{p},j_{p},m_{p},t_{z_{p}}\}$ corresponding to the radial quantum number, orbital angular momentum, total angular momentum, total angular momentum projection, and isospin projection, respectively. Performing normal ordering with respect to a reference state characterized by a one-body density matrix $\rho_{p^{\prime}p}=\langle a^{\dagger}_{p^{\prime}}a_{p}\rangle$ and discarding the residual 3N part, we obtain the NO2B Hamiltonian:

where the braces $\{\ldots\}$ indicate that the enclosed string of creation and annihilation operators are normal-ordered with respect to the used reference state. The Hamiltonian is now expressed in terms of a zero-body part

a normal-ordered one-body part

and a normal-ordered two-body part

The accuracy of the NO2B approximation has been investigated for ground state energies [38, 42, 43], where it was found that by ¹⁶O the error is at the level of 1% of the binding energy. With increasing mass number, this error should decrease as a fraction of the total binding energy¹¹1 The approximation also breaks translational invariance [43], but this is only important for light nuclei (i.e. $A\lesssim 16$), where the NO2B truncation is not necessary and convergence in $E_{3\rm max}$ can be obtained by conventional methods..

If the one-body density matrix $\rho_{pp^{\prime}}$ is rotationally invariant and conserves parity and isospin projection, it must satisfy $(\ell_{p^{\prime}},j_{p^{\prime}},m_{p^{\prime}},t_{z_{p^{\prime}}})=(\ell_{p},j_{p},m_{p},t_{z_{p}})$, and the number of required 3N matrix elements is drastically smaller than that of the original full set. This condition is satisfied for single-reference calculations (e.g. coupled cluster, self-consistent Green’s function, IMSRG, HF-MBPT) with a closed-shell reference, as well as for particle-attached and particle-removed methods [6], and the ensemble normal ordering reference used in the valence-space IMSRG [44]. On the other hand, for broken-symmetry [45, 46] or multi-configurational [47, 48] references necessary to describe e.g. well-deformed nuclei, this would constitute an additional approximation.

Furthermore, in the practically used $JT$-coupled representation, we can sum up the 3N total angular momentum dependence. This can be seen in the $J$-coupled expression for the normal-ordered matrix element, the full expressions for which are provided in Appendix A. Here we show only the contributions from the three-nucleon interactions $V^{\rm 3N}$ (with the notation $[x]\equiv 2x+1$, and using un-normalized matrix elements):

We can see that in Eq. (6) all terms depend on $V^{\rm 3\textbf{}N}$ through the quantity

where the symbol $\tilde{\delta}_{r^{\prime}r}$ is shorthand for all the quantum numbers which are conserved by the one-body density matrix $\rho_{r^{\prime}r}$. If, instead of the full $V^{\rm 3N}$, we only store the quantity (7), we obtain the curve “NO2B” in Fig. 1, allowing us to access $E_{3\rm max}=26$ $(28)$ using single- (half-) precision floating point numbers.

Note that for a Hartree-Fock (HF) calculation, we only need the combination

The number of matrix elements (8) is sufficiently low that we can store the full $\overline{\mathcal{V}}_{p^{\prime}q^{\prime}r^{\prime}pqr}$ without any $E_{3\rm max}$ truncation. However, we find that the HF part of the calculation converges at lower $E_{3\rm max}$ than the beyond-mean-field corrections, which is why we store the quantity $\mathcal{V}^{J_{pq}}_{p^{\prime}q^{\prime}r^{\prime}pqr}$.

A similar idea was employed in Ref. [20]. However, in that work an iterative procedure was adopted in which a HF calculation was performed at a manageable $E_{3\rm max}=14$, and then the lab-frame matrix elements necessary for the NO2B approximation for larger $E_{3\rm max}$ were computed from the relative-basis matrix elements, and the procedure was iterated until self consistency was attained. In our approach, the transformation to the lab frame is performed once, and the resulting matrix elements $\mathcal{V}$ are written to disk and can be used for future calculations of any desired nucleus, without the need for iteration.

Although we can compress the file size by calculating only the NO2B relevant matrix elements via Eq. (7), we still need an efficient way to perform the transformation from the three-body Jacobi basis to single-particle basis. Originally, this transformation was derived for the three-nucleon single-particle $m$-scheme basis [12, 49]. Memory requirements for the $m$-scheme storage limited calculations to $E_{3\rm max}{=}9$. Later, a $j$-coupled storage scheme was introduced [50, 51] that allowed calculations with $E_{3\rm max}\lesssim 18$ as discussed in the introduction with file size requirements illustrated in Fig. 1. Here, we specify the angular momentum coupling in detail. The antisymmetrized three-body states in the lab frame are defined as

with the antisymmetrizer

defined in terms of the permutation operator $P_{ij}$. In Eq.(9), the symbol $\mathcal{C}$ indicate a Clebsch-Gordan coefficient. A state in the antisymmetrized Jacobi basis is denoted $|NiJ_{\rm rel}\rangle$, with total oscillator quanta $N$, total angular momentum $J_{\rm rel}$, and an additional quantum number $i$ to distinguish the states. The transformation from the Jacobi basis to the lab frame may be expressed as

The quantity $\langle N_{\rm cm}L_{\rm cm}NiJ_{\rm rel}:JT|pqr:J_{pq}T_{pq}JT\rangle$ denotes the transformation coefficient. The quantum numbers $N_{\rm cm}$ and $L_{\rm cm}$ are the radial nodes and orbital angular momentum of the center-of-mass (c.m.) motion. The summations over $N,i,N^{\prime},i^{\prime}$ can be performed efficiently by matrix-matrix multiplication, and the remaining summations over $N_{\rm cm}$, $L_{\rm cm}$ and $J_{\rm rel}$ can be computed manually.

The transformation coefficient can be calculated through the non-antisymmetrized Jacobi state:

The index $\alpha$ labels the set of Jacobi quantum numbers $\alpha=\{n_{12},l_{12},s_{12},j_{12},t_{12},n_{3},l_{3},j_{3}\}$. The quantum numbers $\{n_{12},l_{12},s_{12},j_{12},t_{12}\}$ are used for the relative motion of nucleons 1 and 2, i.e., the nodal, orbital angular momentum, spin, total angular momentum, and total isospin quantum numbers, respectively. Similarly, the quantum numbers $\{n_{3},l_{3},j_{3}\}$ correspond to the motion of nucleon 3 with respect to the c.m. of nucleons 1 and 2. Since the antisymmetrized state $|NiJ_{\rm rel}\rangle$ is an eigenstate of the antisymmetrizer $\mathcal{A}$, the coefficient $\langle NiJ_{\rm rel}|N\alpha J_{\rm rel}\rangle$ is also known as the coefficient of fractional parentage [52, 53]. The coefficient $\langle N_{\rm cm}L_{\rm cm}N\alpha J_{\rm rel}:JT|pqr:J_{pq}T_{pq}JT\rangle$ is known as the $T$-coefficient [49, 51], and is the bottleneck of the calculation. It turns out one can sum up three of the angular momentum sums in Eq. (B11) in Ref. [49] ($S_{3},L_{3},{\cal L}$) and obtain a significantly more efficient expression for the $T$-coefficient:

Here, as above, we use the notation $[x]\equiv 2x+1$ and the usual 6-$j$ and 9-$j$ symbols are used. In addition, we use a 12-$j$ symbol of the first kind [54], and $\langle\ldots|\ldots\rangle_{d}$ is Talmi-Moshinsky bracket with mass ratio $d$ [55]. For efficiency, 12-$j$ symbols are calculated on the fly from cached 6-$j$ symbols [54]. We have also used $s_{p}=s_{q}=s_{r}=1/2$. While (13) is a complicated expression, it involves four nested summations (including the expansion of the 12-$j$ symbol; the sum over $N_{12}$ is trivial by energy conservation), rather than the six needed for the expression in Refs. [49, 51].

Our implementation of the above expressions allows us to generate all the lab-frame three-body matrix elements (with half-precision floating point numbers) needed for a calculation employing the NO2B approximation with a spherical reference state up to $E_{3\rm max}=28$ using $\sim 10^{5}$ CPU hours with 187 GB RAM per node. Importantly, this step only needs to be done once for a given interaction. Subsequent many-body calculations for different nuclei, or using different methods can be performed using the same file.

Here, we investigate large-$E_{3\rm max}$ calculations around ¹³²Sn using the well-established NN+3N interaction 1.8/2.0 (EM) [39], which accurately reproduces binding energies to $A\sim 100$ [66, 17, 67]. We employ an oscillator basis with frequency $\hbar\omega=16$ MeV, which is near the optimal value giving the most rapid $e_{\rm max}$ convergence for the ground-state energies and radii of the medium-mass nuclei (converged results are independent of $\hbar\omega$) [66].

One important feature of the 1.8/2.0 (EM) interaction for our purposes is that, while the NN force is softened by a free-space similarity renormalization group (SRG) evolution to a scale $\lambda_{\rm SRG}=1.8~{}{\rm fm}^{-1}$, the corresponding 3N interactions are not SRG evolved. Instead, the cutoff is chosen to be $\Lambda_{\rm 3N}=2.0~{}{\rm fm}^{-1}$ and the short-range low-energy constants $c_{D}$ and $c_{E}$ are refit to the triton binding energy and ⁴He radius. This means that we can avoid SRG evolution of the 3N interaction, which introduces additional challenges due to basis truncations (we address these in section IV.2).

In Fig. 2 we show the ground-state energy of ¹³²Sn calculated with HF-MBPT and IMSRG as a function of $E_{3\rm max}$. The non-variational nature mentioned above is evident, and is present even at the mean-field level. We see that truncations at $E_{3\rm max}=22$ or $24$ are sufficient to obtain convergence within a few MeV. For all points in Fig. 2, the 3N matrix elements are stored and read in using half-precision floating point numbers to reduce the memory footprint. Up to $E_{3\rm max}=24$, we can use single-precision numbers to check the impact of this choice. At $e_{\rm max}=14$, $E_{3\rm max}=24$, the half-precision calculation yields HF energies shifted by $-2.14$ MeV, while the second- and third-order MBPT corrections are changed by 0.68 MeV and 0.11 MeV, respectively, yielding a total difference up to third order of $-1.35$ MeV. This is completely negligible compared with uncertainties arising from many-body truncations (which we expect to be on the order of 20 MeV here³³3This estimate is based on the difference between the MBPT(2), MBPT(3) and IMSRG(2) energies, and is consistent with Ref. [60] where the error at MBPT(3) for similarly soft interactions was found to be 0.1-0.2 MeV per particle. We have further corroborated this estimate with MBPT(4) calculations in a smaller $e_{\rm max}$ space.) and the interaction itself. We also show in Fig. 2 the convergence with respect to $e_{\rm max}$. At $E_{3\rm max}=28$, the third-order energies for $e_{\rm max}=14,16,18$, are $-1115.85$ MeV, $-1117.61$ MeV, and $-1118.16$ MeV, respectively, demonstrating convergence at the $1$ MeV level.

Since the second-order correction of $\sim-300$ MeV is much larger than third-order correction of $\sim-20$ MeV, the correlation energy is dominated by second-order correction. This supports the claim that the extrapolation formula Eq. (21) based on the second-order energy correction is applicable in the case of the HF-MBPT(3) and IMSRG, which includes correlations beyond second order. In panel (a) of Fig. 3, we show $n=2,4,6$ curves of Eq. (21) fitted with the HF-MBPT(2) and IMSRG energy results at $e_{\max}=14$, indicated by the solid symbols in the panel. We see that Eq. (21) works for IMSRG energies as well. Panels (b) and (c) show the extrapolated energies to $E_{3\rm max}=28$, which is the largest value we can calculate. Since the extrapolated point is finite, the uncertainty of all the fitting parameters can propagate to uncertainty of the extrapolated energies. The uncertainty of the energy is estimated as the standard deviation of the 10000 samples generated with the covariance matrix from the fit. Comparing the extrapolated and calculated energies, we see that $n=2$ (Gaussian) reproduces the energies for both HF-MBPT(2) and IMSRG cases, and $n=2$ is the most likely to reproduce the convergence behavior in this case. With $n=2$ formula, we observed that the extrapolated energy to $E_{3\max}=42$ is $-1110.57(2)$ [$-1097.13(2)$] MeV using the IMSRG [HF-MBPT(2)] data $18\leq E_{3\max}\leq 23$.

As already mentioned in Sec. I, we have observed a lack of convergence with respect to $E_{3\rm max}$ in some calculations of heavier systems. One particular example is ¹²⁷Cd as discussed in Ref. [21]. We revisit the calculations in that work, obtained with the VS-IMSRG, and extend them to larger $E_{3\rm max}$. Here, our single-particle basis truncation is $e_{\max}=14$, and we take the valence space as $\{1p_{3/2},1p_{1/2},0f_{5/2},0g_{9/2}\}$ for protons and $\{1d_{5/2},2s_{1/2},1d_{3/2},0g_{7/2},0h_{11/2}\}$ for neutrons above ⁷⁸Ni core. As seen in Fig 4, by $E_{3\rm max}=28$ we obtain convergence in excitation energies at the level of 5 keV. With the previous limit of $E_{3\rm max}=18$ there is no sign of convergence. This behavior can be understood by noting that the $h_{11/2}$ orbit with $e\geq 5$, is impacted by the $E_{3\rm max}$ cut differently than the other neutron valence orbits which have $e\geq 4$, and that the parity of a state is driven by the occupation of the $h_{11/2}$. An analogous argument applies to the proton orbits.

In contrast, when two states have the same number of oscillator quanta in their naive configurations, we expect that their convergence with respect to $E_{3\rm max}$ will be similar and so the energy difference will be less sensitive to the $E_{3\rm max}$ truncation. To illustrate this, we present in the top panel of Fig 5 the first $2^{+}$ excitation energies of even-mass tin isotopes, obtained with the VS-IMSRG. The valence space is $\{1p_{3/2},1p_{1/2},0f_{5/2},0g_{9/2}\}$ for protons and $\{2s_{1/2},1d_{3/2},0h_{11/2},1f_{7/2}\}$ for neutrons above a ⁹²Ni core, indicated by the open symbols in the figure. During the IMSRG evolution, the center-of-mass (c.m.) motions are separated with the Glöckner-Lawson prescription [69] with the coefficient $\beta=3$, and we observe the stability with respect to $\beta$ (see [70] for a detailed discussion). The ground-state energies are converged within approximately 2 MeV. It is clear from the figure that the $2^{+}$ energies show convergence as $E_{3\rm max}$ is increased and $E_{3\rm max}=18$ is sufficient to see the systematics of the $2^{+}$ energies. We note that the $2^{+}$ energy of ¹³²Sn at $E_{3\rm max}=18$ and $24$ differ by 200 keV. We successfully reproduce the $A$-independent excitation energies of the open-shell nuclei, consistent with the seniority picture. In fact, the analysis of the calculated wave function of ^126-130Sn reveals that our valence-space wave functions of ground and first $2^{+}$ states are dominated more than $70$% by the seniority $v=0$ and $v=2$ states, respectively⁴⁴4While such quantitative details about the wave function will in general depend on the details of the IMSRG transformation, this is a relatively simple way to understand the convergence behavior.. The relatively fast convergence of the excitation energies with respect to $E_{3\rm max}$ reflects the fact that both the ground and excited states are dominated by configurations with the same occupancies in the oscillator basis. On the other hand, the excitation at $N=82$ is dominated by a single neutron excitation $0h_{11/2}\to 1f_{7/2}$. As these orbits have the same naive number of oscillator quanta, dependence on the $E_{3\rm max}$ is still mild. The predicted excitation energy is about 1.5 MeV above the experimental value, which is attributed to the IMSRG(2) approximation, as seen in earlier works [66, 71]. Efforts to go beyond the IMSRG(2) approximation are underway [72].

Finally, we consider the convergence behavior of point-proton and point-neutron radii through the Hartree-Fock, second-order HF-MBPT, and IMSRG(2). The diagrams taken into account in the second-order HF-MBPT are listed in Appendix B. The charge radii of several isotopes including ¹³²Sn were recently computed with the self-consistent Green’s function method using chiral forces up to $e_{\rm max}=13$ and $E_{3\rm max}=16$ [18]. We compute point-proton and point-neutron root-mean-squared radii and the neutron skin of ¹³²Sn as a function of $E_{3\rm max}$, and plot the result in Fig. 6. We see that convergence is achieved by $E_{3\rm max}\sim 22$. The corresponding converged charge radii are $r_{ch}=4.43$ fm and $4.42$ fm with IMSRG and second-order HF-MBPT, respectively, demonstrating that the effect of the many-body truncation is controllable for radii. Eq. (21) with $n=2$ reasonably captures the asymptotic convergence behavior of radii. Also we note that the often-used N²LO_sat interaction is harder than the interaction employed here, and thus we would expect the calculations with N²LO_sat will show slower convergence with respect to $E_{\rm 3max}$. Our converged neutron skin with IMSRG(2) is 0.2202(4)—where the quoted uncertainty only accounts for the $E_{3\rm max}$ truncation—consistent with the (model-dependent) extraction [73] of 0.24(4).

The NN and 3N contributions to the 1.8/2.0 (EM) interaction, used for the calculations discussed in the previous section, are defined at different cutoff and resolution scales. For a more systematic convergence study of the calculations it would be desirable explore the resolution-scale and cutoff dependence of observables. Using interactions with a higher cutoff, observables in heavy nuclei will be impossible to converge in the largest feasible model spaces, even with the advances discussed in this work. Therefore, these interactions first need to be softened via a free-space SRG evolution [74] (or some other procedure [75, 76, 77]). For the following calculations we evolve NN and 3N sectors consistently in the harmonic-oscillator basis space. For the NN sector, the evolution is done within the space spanned by the principle quantum number of the relative motion up to 200. Assuming our typical basis frequency of a few tens MeV, the UV scale of this space is a few GeV$/c$—sufficiently larger than the typical momentum scale of $\sim 500$ MeV$/c$ of the bare NN interaction from the chiral EFT, and we can safely evolve the NN Hamiltonian.

For the 3N sector, we evolve the 3N Hamiltonian within the space defined by the three-body principle quantum number $N_{3\max}$, the sum of the principle quantum numbers of the motions for corresponding Jacobi variables. Since the 3N evolution is computationally demanding compared to the NN evolution, we cannot handle a value of $N_{3\max}$ well beyond the typical nuclear interaction scale. We therefore need to investigate the $N_{3\max}$ dependence as we move to heavier systems, as done in Ref. [20]. In the following, we use the chiral N³LO NN interaction from Entem and Machleidt [78] and the N²LO 3N interaction with both local and non-local regulators developed in Ref. [40] denoted as NN+3N(lnl).

In Fig. 7 we show the ground-state energy of ¹³²Sn as a function of $N_{\rm 3max}$. Because the Hamiltonian is block diagonal in the relative angular momentum $J_{\rm rel}$, we can apply a different $N_{\rm 3max}$ cut to each $J_{\rm rel}$ block. We include all channels up to $J_{\rm rel}\leq E_{3\rm max}+3/2=47/2$, which is the highest value that can contribute. To simplify the analysis, we divide the $J_{\rm rel}$ blocks into low-$J$ ($J_{\rm rel}\leq 13/2$), middle-$J$ ($15/2\leq J_{\rm rel}\leq 21/2$), and high-$J$ ($J_{\rm rel}\geq 23/2$) partitions, and vary $N_{\rm 3max}$ for each partition. The SRG evolution is run to a scale of $\lambda_{\rm SRG}=2.0$ fm^-1, working with a basis frequency $\hbar\omega=30$ MeV. After the evolution, the frequency is converted to $\hbar\omega=15$ MeV for the many-body calculations (see Ref. [51] for details). The many-body calculations are done with third-order HF-MBPT with lab-frame truncations $e_{\rm max}=16$ and $E_{3\rm max}=22$. In Fig. 7 we see that the low-$J$ and high-$J$ partitions are converged within the level of a few MeV by $N_{\rm 3max}=48$ and $N_{\rm 3max}=42$, respectively, while the middle-$J$ region converges more slowly. To our knowledge, these are the largest $N_{\rm 3max}$ spaces explored in the literature. It appears that SRG evolution of the $J_{\rm rel}\gtrsim 15/2$ blocks such that heavy nuclei are converged with respect to $N_{\rm 3max}$ will not be possible in the near term without further technical developments.

An alternative truncation scheme we can explore is the maximum value of $J_{\rm rel}$. Since the nuclear interaction is short range, we naively expect that the high $J_{\rm rel}$ components are suppressed by the angular momentum barrier. In the top panel of Fig. 8, the ground-state energy of ¹³²Sn is shown as a function of the maximum $J_{\rm rel}$ included in the transformation Eq. (11). The points labeled “Full” use a uniform $N_{\rm 3max}$ for all blocks up to $J_{\rm rel}=47/2$. Again, energies are computed with HF-MBPT(3) at $\hbar\omega=15$ MeV, $e_{\max}=16$, and $E_{\rm 3max}=22$. We observe that as we increase $N_{\rm 3max}$, the contribution of channels with $J_{\rm rel}>9/2$ becomes essentially negligible.

In order to extrapolate to $N_{\rm 3max}\to\infty$, we fit the calculated energies with an exponential function $E=a\exp(-bN_{\rm 3max})+E_{\infty}$ as shown in the middle panel. In the fitting procedure, we used the energies at $N_{\rm 3max}=36,38,40,42,44$, and used the $N_{\rm 3max}=48$ points to validate the assumed functional form.⁵⁵5We also tried fitting with the Gaussian function $E=a\exp(-bN_{\rm 3max}^{2})+E_{\infty}$, and found this does not provide consistent results with the computed $N_{\rm 3max}=48$ energies. The $N_{\rm 3max}$ extrapolated energies are shown in the bottom panel. The agreement of calculated and extrapolated energies at $N_{\rm 3max}=48$ validates the fitting formula employed here. The final energy obtained by extrapolating the ”full” results to $N_{\rm 3max}\to\infty$ is $-1105.6(15)$ MeV, which agrees within the error bars with the extrapolated energies from $\max(J_{\rm rel})=9/2,11/2,13/2$ results. This reinforces the observation that the contributions from channels with $J_{\rm rel}>9/2$ are negligible in the $N_{\rm 3max}\to\infty$ limit.

This result is somewhat surprising as it suggests that we can obtain a more accurate result by neglecting the high-$J_{\rm rel}$ sector altogether than we can by evolving it in the largest space we can manage. Evidently, the main impact of the high-$J_{\rm rel}$ matrix elements is to introduce an artifact due to the $N_{\rm 3max}$ truncation, which is removed in the limit $N_{\rm 3max}\to\infty$. To investigate the origin of this artifact, in Fig. 9 we hold fixed $N_{\rm 3max}=48$ for the $J_{\rm rel}\leq 13/2$ partition, and plot the expectation values $\langle T+V^{\rm NN}\rangle$, $\langle V^{\rm 3N}_{\rm ind}\rangle$, $\langle V^{\rm 3N}_{\rm gen}\rangle$, and $\langle H\rangle$ as a function of the $N_{3\max}$ cut applied to the $J_{\rm rel}\geq 15/2$ partition. Here $T$ is the relative kinetic energy, $V^{\rm NN}$ is the evolved NN potential, $V^{\rm 3N}_{\rm ind}$ is the induced 3N potential, $V^{\rm 3N}_{\rm gen}$ is the evolved “genuine” 3N potential, and $H$ is the transformed Hamiltonian obtained by summing all of the kinetic and potential terms. The expectation values are taken in a naive harmonic oscillator ground state of ¹³²Sn.

At $N_{\rm 3max}=0$, corresponding to the $J=13/2$ point in Fig. 8, we obtain a bound energy. With increasing $N_{\rm 3max}$ the energy shoots up to 15 GeV, driven by the $\langle V^{\rm 3N}_{\rm ind}\rangle$ component, before converging back towards the $N_{\rm 3max}=0$ value. It appears that the impact the high-$J_{\rm rel}$ matrix elements are negligible. Similar behavior is found in ⁷⁸Ni, where the fully-converged and $N_{\rm 3max}=0$ HF energies differ by 0.3 MeV.

To further investigate the enormous contributions from induced 3N interactions, we decompose $\langle V^{\rm 3N}_{\rm ind}\rangle$ into terms induced by transforming the one-pion exchange, two-pion exchange, and contact parts of $V_{\rm NN}$, as well as the kinetic energy. We find that at $N_{\rm 3max}=16$, all four of these induced terms contribute several GeV to the energy in Fig. 9, indicating that this behavior is generic and not tied to the detailed structure of the NN interaction. Further understanding of the mechanism of this large induced component should be pursued, as it may point the way to a more efficient treatment.

Through this analysis, we conclude that one can perform a more accurate 3N SRG evolution with a truncation in $J_{\rm rel}$, rather than using all possible $J_{\rm rel}$ channels without fully achieving the convergence with respect to $N_{\rm 3max}$. We leave for future work the question of whether this holds for other operators.

Finally, we demonstrate that the asymptotic convergence in $E_{3\rm max}$ discussed in section III is also observed for a consistently SRG-evolved NN+3N interaction. In Fig. 10, we show the 3rd order HF-MBPT ground-state energy of ¹³²Sn as a function of $E_{3\rm max}$, at multiple values of $N_{\rm 3max}$ for the case $J_{\rm rel}\leq 13/2$ (similar behavior is observed for $J_{\rm rel}\leq 9/2$). In contrast to the unevolved case, we observe an increase in the energy for large $E_{3\rm max}$. This bump diminishes with increasing $N_{\rm 3max}$, indicating that the truncation artifact shows up most significantly in the large $E_{3\rm max}$ matrix elements, as would be expected. For each $E_{3\rm max}$, we extrapolate to $N_{\rm 3max}\to\infty$ using an exponential form, and we obtain the gray squares in Fig. 10 (the extrapolation uncertainties are smaller than the markers).

The extrapolated points still display a minimum as a function of $E_{3\rm max}$ before converging to the final answer from below. The decreasing trend below $E_{3\rm max}=20$ is driven by the convergence of the HF energy, while the increase above $E_{3\rm max}=20$ is driven by the second order MBPT correction. The fact that the energy converges from below in this case supports the assumption in section III that the asymptotic convergence in $E_{3\rm max}$ is driven by the $V_{NN}$-$V_{3N}$ cross-term , which can be either positive or negative. The asymptotic behavior is fit well with a Gaussian with similar parameters (aside from the overal sign) to those in the unevolved case.