1 Introduction↩︎

The nuclear radius, which characterizes the size of the nucleus, is one of the most fundamental properties of atomic nuclei. With the advancement of experimental techniques, there has been a rapid increase in available data on nuclear charge and matter radii for both stable and unstable isotopes [1]–[4]. On the one hand, these data have significantly enhanced our understanding of nuclear structure, revealing novel phenomena such as nuclear halos [5], [6], neutron skins [7], [8], and nuclear shape phase transitions [9], [10]. On the other hand, they pose considerable challenges to nuclear theory, particularly to ab initio approaches.

Despite the remarkable progress in nuclear theory, achieving accurate descriptions of nuclear charge radii remains a significant challenge [11]–[17]. A notable example is the oxygen isotope chain. Most ab initio calculations underestimate the charge radius of \(^{16}\)O by roughly 7–15% [11], [18]. Furthermore, while \(^{16}\)O and \(^{17}\)O exhibit similar charge radii, a pronounced increase in the charge radius is observed experimentally for \(^{18}\)O in comparison to \(^{16}\)O and \(^{17}\)O [2]. To our knowledge, this intriguing feature is not captured by any existing theoretical approach. Although the chiral interaction \({\rm NNLO}_{\rm sat}\) successfully reproduces the charge radius of \(^{16}\)O by including it among the fitted observables [18], it still fails to give the substantial increase seen in \(^{18}\)O [13], [16].

Nuclear lattice effective field theory (NLEFT) [19], [20] offers a different perspective on this problem. Recently, the wave function matching method [21] was proposed to mitigate the notorious Monte Carlo sign problem in the imaginary-time evolution of quantum many-body simulations, which introduces strong cancellations between positive and negative amplitudes. In combination with high-fidelity N\(^3\)LO chiral interactions, the binding energies of both light and medium-mass nuclei up to \(A=58\) are well reproduced. While charge radii are generally well described, including the pronounced increase in \(^{18}\)O compared to \(^{16,17}\)O, they are often accompanied by significantly larger statistical uncertainties. This is largely due to the commonly used pinhole algorithm, which inserts an \(A\)-body density operator during the imaginary-time evolution [22]. This insertion leads to unpaired nucleons which results in additional strong sign-oscillations, particularly for nuclei with larger mass numbers and/or pronounced proton–neutron asymmetry, thereby complicating the investigation of subtle isotopic trends in nuclear radii.

In this work, we propose an extension of the pinhole algorithm, referred to as the partial pinhole algorithm, that reduces the rank of the many-body density operator in order to mitigate the sign problem in nuclear radius calculations within the framework of NLEFT. This reduced operator is implemented via the recently developed rank-operator method [23]. We combine the partial pinhole algorithm with the wave function matching method using high-fidelity N\(^3\)LO chiral interactions. As a first application, we investigate both charge and matter radii of a number of oxygen isotopes.

2 Formalism↩︎

The full details of our lattice calculations are shown in the Supplemental Material [24]. To illustrate the idea of the partial pinhole algorithm and how it is combined with NLEFT in calculating nuclear radii, we adopt a simplified Hamiltonian with an attractive two-body contact short range interaction (\(C<0\)), \[H=K+\frac{1}{2}C\,\sum_{\boldsymbol{n}}:\rho^2(\boldsymbol{n}):,\] where \(K\) is the kinetic energy term with nucleon mass \(m = 938.9\) MeV, the colons indicate normal ordering, and \(\rho(\boldsymbol{n})=\sum_{\sigma\tau}\rho_{\sigma\tau}(\boldsymbol{n})=\sum_{\sigma\tau}a^\dagger_{\sigma\tau}(\boldsymbol{n})a_{\sigma\tau}(\boldsymbol{n})\). Here, \(\rho_{\sigma\tau}(\boldsymbol{n})\) is the one-body density operator at lattice site \(\boldsymbol{n}\) for spin \(\sigma\) and isospin \(\tau\).

The ground state wave functions of \(H\) can be obtained by applying imaginary time projectors to a trial wave function \(\Psi_0\), \(|\Psi\rangle=\lim_{\tau\to\infty}e^{-H\tau/2}|\Psi_0\rangle\). In Monte Carlo simulations, \(\tau\) is divided into \(L_t\) slices with temporal spacing \(a_t\) such that \(\tau= L_ta_t\). For each time slice, the two-body interaction is defined as nucleons propagated in a fluctuating background auxiliary field using a Hubbard-Stratonovich transformation, \[\exp\left(-\frac{a_tC}{2}\rho^2\right) = \sqrt{\frac{1}{2\pi}} \int ds~\exp\left(-\frac{s^2}{2}+\sqrt{-a_tC}s\rho\right),\] where \(s\) is the auxiliary field at a lattice site. Therefore, the wave function \(|\Psi\rangle\) can be written as an auxiliary field path-integral expression, \[\label{Eq:wfs} |\Psi\rangle = \int \mathcal{D}{s_1}\cdots \mathcal{D}{s_{L_t/2}} M(s_{Lt/2})\cdots M(s_1) |\Psi_0\rangle~,\tag{1}\] where the initial wave function \(|\Psi_0\rangle\) is a \(A\)-nucleon Slater determinant, such as alpha clusters states or shell-model wave functions, and the normal-ordered transfer matrix \(M\) is defined as, \[M(s_{n_t}) = :\exp\left[-a_tK+\sqrt{-a_tC}\sum_{\boldsymbol{n}}s_{n_t}(\boldsymbol{n})\rho(\boldsymbol{n})\right]:.\] Therefore, \(|\Psi\rangle\) is expressed as a linear combination of \(A\)-nucleon Slater determinants, each of which is associated with a specific configuration of auxiliary fields \(\vec{s}\).

With the wave function \(|\Psi\rangle\) in Eq. 1 , the root-mean-square (RMS) point-proton radius of a nucleus with proton number \(Z\), neutron number \(N\), and mass number \(A\) can be evaluated by \[\langle r_p^2\rangle = \frac{1}{Z} \sum_{\boldsymbol{n},\sigma,\tau=p} \langle \rho_{\sigma\tau}(\boldsymbol{n}) (\boldsymbol{n} - \hat{\boldsymbol{r}}_{\rm c.m.})^2 \rangle,\] where \(\hat{\boldsymbol{r}}_{\rm c.m.} = \sum_{\boldsymbol{n},\sigma,\tau} \boldsymbol{n} \rho_{\sigma\tau}(\boldsymbol{n}) / A\) is the center-of-mass operator. The expression can be rewritten in terms of a set of two-body correlation functions \(G_{\alpha\beta}\) as \[\langle r_p^2\rangle = \frac{N}{ZA^2} G_{pp} + \frac{N-Z}{ZA^2} G_{pn} - \frac{1}{A^2} G_{nn}~,\] with \[G_{\alpha\beta} = \sum_{\substack{\boldsymbol{n}_1,\sigma_1,\tau_1 = \alpha \\ \boldsymbol{n}_2,\sigma_2,\tau_2 = \beta}} \langle :\rho_{\sigma_1\tau_1}(\boldsymbol{n}_1) \rho_{\sigma_2\tau_2}(\boldsymbol{n}_2): \rangle (\boldsymbol{n}_1 - \boldsymbol{n}_2)^2.\] Similarly, the matter radius is given by \[\begin{align} \langle r_m^2\rangle &= \frac{1}{A} \sum_{\boldsymbol{n},\sigma,\tau} \langle \rho_{\sigma\tau}(\boldsymbol{n}) (\boldsymbol{n} - \hat{\boldsymbol{r}}_{\rm c.m.})^2 \rangle \nonumber \\ &= \frac{1}{2A^2} G_{pp} + \frac{1}{A^2} G_{pn} + \frac{1}{2A^2} G_{nn}. \end{align}\] Thus, nuclear radii can be expressed in terms of just three two-body correlation functions. It is important to note that since the matter radius corresponds to the second radial moment of the intrinsic nuclear density distribution (or equivalently, the squared nuclear wave function), its numerical value is sensitive to both the theoretical framework used and the specific methods adopted when extracting it from experimental observables.

Similar to the previous studies on computing the RMS charge radii [18], [21], [25], we use the standard relation [26] \[r_{\rm ch}^2 = \langle r_p^2 \rangle + R_p^2 + \frac{N}{Z} R_n^2 + \frac{3}{4m_p^2},\] where \(R_p^2 = 0.7056~{\rm fm}^2\) [27], [28], \(R_n^2 = -0.105~{\rm fm}^2\) [29], and \(m_p = 938.27\) MeV. Relativistic spin-orbit corrections and two-nucleon current contributions are not included in this work.

A direct calculation of the correlation functions \(G_{\alpha\beta}\) is computationally expensive, especially when the corrections from high-fidelity chiral interactions are included perturbatively. To evaluate them more efficiently, we introduce the partial pinhole algorithm. We define a normal-ordered \(M\)-body density operator, \[\label{eq:rhoM} \rho_M(\vec{c}\,) = :\rho_{\sigma_1\tau_1}(\boldsymbol{n}_1) \cdots \rho_{\sigma_M\tau_M}(\boldsymbol{n}_M):,\tag{2}\] using the notation \(\vec{c} = (c_1, \cdots, c_M)\). Here, \(c_i = (\boldsymbol{n}_i, \sigma_i, \tau_i)\) specifies the quantum numbers of the \(i\)-th density operator, with \(\boldsymbol{n}_i\) denoting the lattice coordinate, and \(\sigma_i\) and \(\tau_i\) representing spin and isospin. The summation over \(\vec{c}\) satisfies \[\begin{align} \sum_{\vec{c}} \rho_M(\vec{c}\,) &= \hat{N}(\hat{N} - 1) \cdots (\hat{N} - M + 1) \nonumber \\ &= A(A - 1) \cdots (A - M + 1), \end{align}\] where \(\hat{N}\) is the particle-number operator. The equality of the second line stems from the fact that we are working in \(A\)-nucleon subspace. Therefore, this summation is proportional to the identity operator.

The correlation functions \(G_{\alpha\beta}\) can then be evaluated as \[\begin{align} \label{Eq:corr95G} G_{\alpha\beta} &= \frac{A(A - 1)}{M(M - 1)} \nonumber \\ &\quad \times \frac{\displaystyle\sum_{\vec{c}} \left[ \langle \Psi | \rho_M(\vec{c}\,) | \Psi \rangle \sum_{i<j} \delta_{\tau_i\alpha} \delta_{\tau_j\beta} (\boldsymbol{n}_i - \boldsymbol{n}_j)^2 \right]}{ \displaystyle\sum_{\vec{c}} \langle \Psi | \rho_M(\vec{c}\,) | \Psi \rangle }. \end{align}\tag{3}\] The details about computing \((\boldsymbol{n}_i - \boldsymbol{n}_j)^2\) are presented in the Supplemental Material [24]. Unlike the standard pinhole algorithm [22], the value of \(M\) is much smaller than the total number of nucleons \(A\). As a result, \(\rho_M\) is no longer a projection operator, and acting with \(\rho_M\) on a single \(A\)-nucleon Slater determinant yields a linear combination of multiple \(A\)-nucleon Slater determinants.

Since \(\rho_{\sigma\tau}(\boldsymbol{n})\) is a rank-one operator, that is, \(:\rho_{\sigma\tau}^2(\boldsymbol{n}): = 0\), the rank-one operator method [23] allows us to express \(\rho_M(\vec{c}\,)\) as \[\rho_M(\vec{c}\,) = \frac{1}{\varepsilon^M} :\exp[ \varepsilon \rho_{\sigma_1\tau_1}(\boldsymbol{n}_1) + \cdots + \varepsilon \rho_{\sigma_M\tau_M}(\boldsymbol{n}_M) ]:\] with \(\varepsilon \to \infty\). It consists of a string of one-body operators which can be applied directly to each single-particle wave function in the Slater determinant.

The summations over \(\vec{c}\) in Eq.@eq:Eq:corr95G and the path integral over \(\vec{s}\) in Eq.@eq:Eq:wfs are evaluated via Monte Carlo importance sampling. We generate an ensemble of \({\vec{c},\vec{s}}\) configurations according to the relative probability weight, \[\begin{align} P(\vec{c},\vec{s}\,) = |\langle \Psi_0 |&M(s_{L_t})\cdots M(s_{L_t/2+1})\rho_M(\vec{c}\,)\nonumber\\ &\times M(s_{L_t/2})\cdots M(s_1)| \Psi_0 \rangle|. \end{align}\] The strings of transfer matrices are generated by updating the auxiliary fields \(\vec{s}\) using the shuttle algorithm [30]. At the middle time slice, the auxiliary field updates are paused, and the pinhole configuration \(\vec{c}\) is updated by performing standard Metropolis accept/reject steps. Specifically, to update \(\vec{c}\), we randomly select an index \(i\) in \(\rho_M\) and either move it to a neighboring site, \[c_i = \{\boldsymbol{n}_i, \sigma_i, \tau_i\} \rightarrow c_i' = \{\boldsymbol{n}_i', \sigma_i, \tau_i\},\] or reassign its spin and isospin, \[c_i = \{\boldsymbol{n}_i, \sigma_i, \tau_i\} \rightarrow c_i' = \{\boldsymbol{n}_i, \sigma_i', \tau_i'\}.\] The new configuration \(\vec{c}\,'\) is accepted if \[{ P(\vec{c}\,',\vec{s}\,)}/{ P(\vec{c},\vec{s}\,)} > r,\] where the random number \(r\) is uniformly distributed between 0 and 1.

In this work, we combine the partial pinhole algorithm with the wave function matching method. The original chiral Hamiltonian is unitarily transformed into a new high-fidelity Hamiltonian \(H\), whose wave function matches that of a computationally simple Hamiltonian, \(H_S\), at a given radius. This ensures rapid convergence of the perturbative expansion in \(H - H_S\). The details can be found in Ref. [21]. The formalism for applying the partial pinhole algorithm in perturbative calculations is provided in [24].

3 Results and discussion↩︎

In the following calculations, we adopt the N\(^3\)LO chiral interactions from Ref. [21]. We also use a minimal sets of three-body force terms, where, in addition to nuclear binding energies, the charge radius of \(^4\)He is used to determine the low-energy constants (LECs) [31]. Although the number of involved free LECs is reduced from eight to six, the overall description of nuclear binding energies is similar to Ref. [21], as detailed in [24]. Our simulations are performed with a spatial lattice spacing of \(a = 1.32~{\rm fm}\), corresponding to a momentum cutoff \(\Lambda = \pi/a \approx 471~{\rm MeV}\). Additionally, we use a temporal lattice spacing of \(a_t = 0.001~{\rm MeV}^{-1}\) and a cubic box of length \(L = 13.2~{\rm fm}\). Since the radii of the oxygen isotopes under investigation are all below 3 fm, this box size is clearly sufficient to suppress finite-volume effects.

In the partial pinhole algorithm, we take \(M = 4\) for \(\rho_M\) in Eq.@eq:eq:rhoM . We have verified that the results are largely insensitive to the specific choice of \(M\) [24]. To enhance computational efficiency, we fix half of the isospin indices in \(\rho_M\) as protons and the other half as neutrons. This setup ensures that each configuration \(\vec{c}\) contributes simultaneously to all three correlation functions \(G_{pp}\), \(G_{pn}\), and \(G_{nn}\).

In [24], we benchmark the charge radii of \(^{16}\)O and \(^{17}\)O using both the standard pinhole and the partial pinhole algorithms with the simple Hamiltonian \(H_S\). Although we find that both methods yield consistent charge radii, the partial pinhole algorithm maintains a significantly larger average phase factor and substantially reduces statistical uncertainties, particularly in the case of \(^{17}\)O. This improvement becomes even more substantial for heavier nuclei such as \(^{40}\)Ca, and is further amplified when applying the full N\(^3\)LO chiral interactions.

Figure 1 shows the charge radii of \(^{16}\)O, \(^{17}\)O, \(^{18}\)O, and \(^{20}\)O calculated using the partial pinhole algorithm together with the chiral interactions at N\(^3\)LO. Due to the substantial suppression of the sign oscillations, we are now able to reach significantly larger imaginary times as compared to earlier works. For \(^{16}\)O, \(^{17}\)O, and \(^{18}\)O, the calculated radii exhibit a clear exponential decay behavior. A double-exponential fit is used to account for residual excited-state contamination and extrapolate the results to \(\tau \to \infty\); see the Supplemental Material [24] for details. The uncertainties of the extrapolation are quantified using the covariance matrix obtained from the nonlinear fit. The resulting 1\(\sigma\) errors are found to be small, consistent with the near-convergence of the data at \(\tau = 0.48~{\rm MeV}^{-1}\). For \(^{20}\)O, however, the stronger sign problem arising from its larger proton–neutron asymmetry prevents the simulation from reaching very large projection times, and the extrapolation is terminated at \(\tau = 0.4~{\rm MeV}^{-1}\). As a result, the extrapolated charge radius exhibits a larger uncertainty.

The extrapolated charge radii are summarized in Table 1, where results using nine additional sets of three-nucleon forces are included to estimate theoretical uncertainties [31]. The results show that \(^{16}\)O and \(^{17}\)O have very similar charge radii, while \(^{18}\)O exhibits a noticeable increase, consistent with experimental data [2]. We note that, although \(^{16}\)O is commonly described as a spherical doubly magic nucleus within shell-model-based approaches, NLEFT simulations predict a distinct tetrahedral arrangement of alpha clusters [32]. Supporting this picture, previous NLEFT studies [33] demonstrated that the ground state of \(^{12}\)C robustly converges to an equilateral triangular arrangement of alpha clusters independent of the choice of the initial state. By analogy, assuming that \(^{16}\)O similarly forms a stable tetrahedral cluster configuration as indicated by NLEFT calculations, this intrinsic clustering could explain the success of NLEFT in reproducing oxygen isotope charge radii.

We also compute the matter radii of oxygen isotopes using the partial pinhole algorithm. The results for \(^{16}\)O, \(^{17}\)O, \(^{18}\)O, and \(^{20}\)O are shown in Fig. 2. Similar to the charge radii shown in Fig. 1, we perform a double-exponential fit to extrapolate the data to \(\tau \to\infty\); see [24] for details. In contrast to the charge radii, the matter radii increase monotonically with mass number. The extrapolated matter radii are listed in Table 1, along with values obtained using nine additional sets of three-nucleon forces [31].

In Table 1, we also compare our results with experimental matter radii extracted from electron and proton scattering [\((e,e)\) and \((p,p)\)] [13], interaction cross sections (\(\sigma_I\)) [1], and charge-exchange cross sections (\(\sigma_{cc}\)) [16]. Unlike charge radii, the extraction of matter radii often involves some degree of model dependence, leading to discrepancies among different experimental methods. For example, matter radii derived from low-energy proton elastic scattering tend to be systematically larger than those extracted from \(\sigma_I\) and \(\sigma_{cc}\) measurements. The most notable case is \(^{20}\)O, where differences up to 0.2 fm are observed among the three different measurements. As discussed in Ref. [13], the \((e,e)\) and \((p,p)\) data provide a more reliable descriptions, since the correlations in the target are usually not included in the other methods. From Table 1, we can find our calculated matter radii agree well with the values extracted from \((e,e)\) and \((p,p)\) data, and thus support the statement given in Ref. [13].

4 Summary and perspectives↩︎

We have performed an ab initio calculation of the charge and matter radii of oxygen isotopes from \(^{16}\)O to \(^{20}\)O using nuclear lattice effective field theory (NLEFT) with the high-fidelity N\(^3\)LO chiral interactions. To mitigate the Monte Carlo sign problem in nuclear radius calculations, we introduce the partial pinhole algorithm, which reduces the rank of the many-body density operator in the standard pinhole approach and significantly lowers the statistical uncertainties. This method enables access to much larger projection times and extends the applicability to more neutron-rich and proton-rich nuclei. Our results accurately reproduce the experimental charge radii of \(^{16}\)O, \(^{17}\)O, and \(^{18}\)O, demonstrating the capability of NLEFT in describing nuclear structure properties. The charge radius of \(^{20}\)O is predicted to be \(2.810(32)\) fm. The calculated matter radii show excellent agreement with values extracted from low-energy electron and proton scattering, while deviating from those inferred from interaction and charge-exchange cross sections. These results underscore the need for further experimental clarification, particularly for neutron-rich isotopes, with improved precision and reduced model dependence.

As a general and versatile method, the partial pinhole algorithm holds broad potential future applications. It can be integrated with other Monte Carlo techniques, such as the perturbative quantum Monte-Carlo method [34], and has already been successfully applied to the calculation of additional observables, including electromagnetic transitions [35]. Its ability to suppress sign oscillations while maintaining high accuracy suggests wide applicability not only across nuclear systems, but also to broader classes of quantum many-body problems.

5 Supplemental Materials↩︎

5.1 Lattice Hamiltonian and wave function matching with the N\(^3\)LO interaction↩︎

In our lattice calculations, we combine the partial pinhole algorithm with the wave function matching method introduced in Ref. Elhatisari:2022zrb_sm?. In this approach, the original chiral Hamiltonian is unitarily transformed into a new high-fidelity Hamiltonian \(H\), whose wave function matches those of a computationally simpler Hamiltonian \(H_S\) at short distances, up to a specified matching radius. This construction ensures rapid convergence of the perturbative expansion in \(H - H_S\).

The full formalism of the wave function matching method is given in Ref. Elhatisari:2022zrb_sm?, and here we summarize the key elements relevant to our study.

The simplified Hamiltonian \(H_S\) is constructed from a leading-order chiral effective field theory (\(\chi\)EFT) interaction, \[\label{eq:H95S} H_S = K + \frac{1}{2}c_{\rm SU(4)}\sum_{\boldsymbol{n}}:\tilde{\rho}^2(\boldsymbol{n}): + V_{\rm OPE}^{\Lambda_\pi=180},\tag{4}\] where \(K\) denotes the kinetic energy term with nucleon mass \(m = 938.92\) MeV, and the colons \(::\) indicate normal ordering. The smeared density operator \(\tilde{\rho}(\boldsymbol{n})\) includes both local and nonlocal components and is defined as \[\tilde{\rho}(\boldsymbol{n}) = \sum_{\sigma\tau}\tilde{a}^\dagger_{\sigma\tau}(\boldsymbol{n})\tilde{a}^{}_{\sigma\tau}(\boldsymbol{n}) + s_L\sum_{|\boldsymbol{n}'-\boldsymbol{n}|=1}\tilde{a}^\dagger_{\sigma\tau}(\boldsymbol{n}')\tilde{a}^{}_{\sigma\tau}(\boldsymbol{n}'),\] with smeared creation and annihilation operators given by \[\tilde{a}^{(\dagger)}_{\sigma\tau}(\boldsymbol{n}) = a^{(\dagger)}_{\sigma\tau}(\boldsymbol{n}) + s_{NL}\sum_{|\boldsymbol{n}'-\boldsymbol{n}|=1}a^{(\dagger)}_{\sigma\tau}(\boldsymbol{n}').\] Throughout our calculations, we use a local smearing parameter \(s_L = 0.07\) and a nonlocal smearing parameter \(s_{NL} = 0.5\). The term \(V_{\rm OPE}^{\Lambda_\pi=180}\) denotes the one-pion-exchange (OPE) potential at leading order, regulated with a cutoff \(\Lambda_\pi = 180\) MeV.

The high-fidelity \(\chi\)EFT Hamiltonian \(H\) at N\(^3\)LO is given by \[\begin{align} \label{eq:H95N3LO} H = K + V_{\rm OPE}^{\Lambda_\pi=300} + V_{\rm Coulomb} + V_{\rm 3N}^{Q^3} + V_{\rm 2N}^{Q^4} + W_{\rm 2N}^{Q^4} + V_{\rm 2N,WFM}^{Q^4} + W_{\rm 2N,WFM}^{Q^4}. \end{align}\tag{5}\] Here, the OPE potentials \(V_{\rm OPE}^{\Lambda_\pi=300}\) has a larger cutoff \(\Lambda_\pi=300\) MeV compared to \(H_S\). The full Hamiltonian also includes the Coulomb interaction (\(V_{\rm Coulomb}\)), the three-nucleon (3N) interaction at order \(Q^3\) (\(V_{\rm 3N}^{Q^3}\)), two-nucleon (2N) short-range interactions at order \(Q^4\) (\(V_{\rm 2N}^{Q^4}\)), Galilean-invariance-restoring (GIR) terms (\(W_{\rm 2N}^{Q^4}\)), and wave function-matching interaction (\(V_{\rm 2N,WFM}^{Q^4}\) as well as the GIR correction of the wave function matching interaction \(W^{Q^4}_{\rm 2N,WFM}\). The contribution of \(V_{\rm 2N,WFM}^{Q^4}\) to nuclear radii is found to be negligible, while its inclusion significantly increases computational cost. Therefore, this term is omitted in our simulations.

All simulations are performed in periodic cubic lattices with a spatial lattice spacing of \(a = 1.32\) fm. The low-energy constants (LECs) of the two-nucleon interaction are taken from Ref. Elhatisari:2022zrb_sm?. For the 3N potential, we adopt a minimal set of three-body force terms (3NFs) that reduce the number of free LECs from eight to six while preserving the overall description of nuclear binding energies Elhatisari:pre_SM?. See Fig. 3 for details.

5.2 Perturbative formalism for the partial pinhole algorithm↩︎

To explain the perturbative formalism, we begin by expressing the two-body correlation functions \(G_{\alpha\beta}\) in terms of the transfer matrix, \[\begin{align} \label{eq:Gab95sm} G_{\alpha\beta} = \frac{A(A - 1)}{M(M - 1)} \frac{ \displaystyle\sum_{\vec{c}} \left[ \langle \Psi_0 |\hat{M}^{L_t/2} \rho_M(\vec{c}\,) \hat{M}^{L_t/2}|\Psi_0\rangle \sum_{i<j} \delta_{\tau_i\alpha} \delta_{\tau_j\beta} (\boldsymbol{n}_i - \boldsymbol{n}_j)^2 \right] }{ \displaystyle\sum_{\vec{c}} \langle \Psi_0 | \hat{M}^{L_t/2}\rho_M(\vec{c}\,) \hat{M}^{L_t/2}| \Psi_0 \rangle }, \end{align}\tag{6}\] with the definition, \[\hat{M} = :\exp(-a_t H):.\] To formulate the perturbative expansion, we decompose the high-fidelity Hamiltonian \(H\) in Eq. 5 as \[H = H_S + (H - H_S),\] where the leading-order Hamiltonian \(H_S\) is used as the nonperturbative part, while the difference \((H - H_S)\) is treated as a perturbation. Accordingly, the transfer matrix can be expanded to first order as \[\begin{align} \hat{M} = :\exp(-a_t H): = :\exp(-a_t H_S): - a_t :\exp(-a_t H_S)(H - H_S): \equiv \hat{M}^{(0)} + \hat{M}^{(1)}. \end{align}\]

Keeping terms up to \(\mathcal{O}(\hat{M}^{(1)})\), the correlation function becomes \[G_{\alpha\beta} = \frac{A(A - 1)}{M(M - 1)}\frac{\mathcal{M}_{\alpha\beta}^{(0)}+\mathcal{M}_{\alpha\beta}^{(1)}}{\mathcal{M}^{(0)}+\mathcal{M}^{(1)}}=\frac{A(A - 1)}{M(M - 1)}\left(\frac{\mathcal{M}_{\alpha\beta}^{(0)}}{\mathcal{M}^{(0)}} +\frac{\mathcal{M}_{\alpha\beta}^{(1)}}{\mathcal{M}^{(0)}} -\frac{\mathcal{M}_{\alpha\beta}^{(0)}\mathcal{M}^{(1)}}{(\mathcal{M}^{(0)})^2} \right),\] The nonperturbative amplitudes are defined as \[\begin{align} &\mathcal{M}^{(0)} = \sum_{\vec{c}} \langle \Psi_0 | (\hat{M}^{(0)})^{L_t/2}\rho_M(\vec{c}\,) (\hat{M}^{(0)})^{L_t/2}| \Psi_0 \rangle,\\ &\mathcal{M}^{(0)}_{\alpha\beta} = \sum_{\vec{c}}{\bigg [} \langle \Psi_0 | (\hat{M}^{(0)})^{L_t/2}\rho_M(\vec{c}\,) (\hat{M}^{(0)})^{L_t/2}| \Psi_0 \rangle\sum_{i<j} \delta_{\tau_i\alpha} \delta_{\tau_j\beta} (\boldsymbol{n}_i - \boldsymbol{n}_j)^2{\bigg ]}~. \end{align}\] The perturbative amplitudes \(\mathcal{M}^{(1)}\) and \(\mathcal{M}_{\alpha\beta}^{(1)}\) are obtained by replacing one of the \(\hat{M}^{(0)}\) factors with \(\hat{M}^{(1)}\) at each time slice and summing over all possible insertion points, \[\begin{align} &\mathcal{M}^{(1)} = \sum_{\vec{c}}\sum_{n_t=0}^{L_t/2-1} \langle \Psi_0 | (\hat{M}^{(0)})^{L_t/2-n_t-1}M^{(1)}(\hat{M}^{(0)})^{n_t}\rho_M(\vec{c}\,) (\hat{M}^{(0)})^{L_t/2}| \Psi_0 \rangle+{\rm c.c.},\\ &\mathcal{M}^{(1)}_{\alpha\beta} = \sum_{\vec{c}}\sum_{n_t=0}^{L_t/2-1} {\bigg [}\langle \Psi | (\hat{M}^{(0)})^{L_t/2-n_t-1}M^{(1)}(\hat{M}^{(0)})^{n_t}\rho_M(\vec{c}\,) (\hat{M}^{(0)})^{L_t/2}| \Psi \rangle\sum_{i<j} \delta_{\tau_i\alpha} \delta_{\tau_j\beta} (\boldsymbol{n}_i - \boldsymbol{n}_j)^2{\bigg ]} +{\rm c.c.}. \end{align}\] In the calculation of all amplitudes \(\mathcal{M}^{(0)}\), \(\mathcal{M}^{(0)}_{\alpha\beta}\), \(\mathcal{M}^{(1)}\), and \(\mathcal{M}^{(1)}_{\alpha\beta}\), \(:\exp(-a_tH_S):\) will be used in the auxiliary-field formalism, and the other terms involving sophisticated \(\chi\)EFT interactions are evaluated via Jacobi formulas.

5.3 The distance between two lattice sites↩︎

To evaluate the correlation function \(G_{\alpha\beta}\), one needs to calculate the squared distance between two lattice sites, \((\boldsymbol{n}_i-\boldsymbol{n}_j)^2\). Due to the periodic boundary conditions used in our lattice calculations, the default manner is to minimize \((\boldsymbol{n}_i-\boldsymbol{n}_j)^2\) by considering all periodic images of the lattice sites. We denote the minimal-distance on the lattice as \((\widetilde{\boldsymbol{n}_i-\boldsymbol{n}_j})^2\). However, this conventional method introduces additional finite volume effects. This can be seen in Fig. 4, the physical distance between \(\boldsymbol{n}_i\) and \(\boldsymbol{n}_j\) will be replaced by the distance between \(\boldsymbol{n}_i\) and \(\boldsymbol{n}_j'\) due to the periodic boundary conditions.

Inspired by the fact that the nuclear radii are defined with respect to the center of mass, we calculate \((\boldsymbol{n}_i-\boldsymbol{n}_j)^2\) by introducing a reference point chosen to be close to the nuclear center of mass, denoted as \(\boldsymbol{r}_{\rm ref}\). From Fig. 4, we can find that the relative vector \((\widetilde{\boldsymbol{n}_i-\boldsymbol{r}_{\rm ref}})\) is always physically correct. We further require that \(\boldsymbol{r}_{\rm ref}\) should coincide exactly with the center of the periodic box so that it is symmetric with respect to all lattice directions. Consequently, the coordinates of \(\boldsymbol{r}_{\rm ref}\) depend on whether the box size \(L\) is odd or even, e.g. along the \(x\)-direction, \[x_{\rm ref} = \begin{cases} n_x, &L={\rm odd},\\ n_x +1/2, &L={\rm even}, \end{cases}\] with \(n_x\) an integer number. The reference points \(\boldsymbol{r}_{\rm ref}\) are then determined by minimizing \[\sum_{\boldsymbol{n}}\langle\rho(\boldsymbol{n})\rangle(\widetilde{\boldsymbol{n}-\boldsymbol{r}_{\rm ref}})^2,\] where the expectation of one-body density \(\langle\rho(\boldsymbol{n})\rangle\) is defined as, \[\langle\rho(\boldsymbol{n})\rangle = \frac{\sum_{\sigma\tau}\langle\Psi|\rho_{\sigma\tau}(\boldsymbol{n}):\rho_{\sigma_i\tau_i}(\boldsymbol{n}_i)\rho_{\sigma_j\tau_j}(\boldsymbol{n}_j):|\Psi\rangle}{\langle\Psi|:\rho_{\sigma_i\tau_i}(\boldsymbol{n}_i)\rho_{\sigma_j\tau_j}(\boldsymbol{n}_j):|\Psi\rangle}.\] With the determined reference point, we first determine the relative vector of \(\boldsymbol{n}_{i,j}\) with considering the periodicity property, and \((\boldsymbol{n}_i-\boldsymbol{n}_j)^2\) is then evaluated as, \[[(\widetilde{\boldsymbol{n}_i-\boldsymbol{r}_{\rm ref}})-(\widetilde{\boldsymbol{n}_j-\boldsymbol{r}_{\rm ref}})]^2= (\widetilde{\boldsymbol{n}_i-\boldsymbol{r}_{\rm ref}})^2-2(\widetilde{\boldsymbol{n}_i-\boldsymbol{r}_{\rm ref}})\cdot(\widetilde{\boldsymbol{n}_j-\boldsymbol{r}_{\rm ref}}) + (\widetilde{\boldsymbol{n}_j-\boldsymbol{r}_{\rm ref}})^2.\]

5.4 Benchmark with the simple Hamiltonian \(H_{S}\)↩︎

In this section, we benchmark the nuclear radii calculated using the partial pinhole algorithm against those obtained from the standard pinhole algorithm Elhatisari:2017eno_SM?. In the standard pinhole algorithm, the positions of \(A\)-nucleons, denoted as \(\boldsymbol{n}_i\), are sampled according to the \(A\)-body density amplitudes. With these \(A\) pinhole coordinates, one can compute density correlations relative to the center of mass, and thus the nuclear radii can be obtained Lu:2018bat_SM?.

For benchmarking purposes, we use the simplified Hamiltonian given in Eq. 4 . To facilitate perturbation theory, we split \(H_S\) into a nonperturbative part \(H^{(0)}\) and a perturbative correction \(H^{(1)}\), given explicitly by, \[\begin{align} &H^{(0)} = K + (1-\lambda)\cdot\frac{1}{2}c_{\rm SU(4)}\sum_{\boldsymbol{n}}:\tilde{\rho}^2(\boldsymbol{n}): + V_{\rm OPE}^{\Lambda_\pi=180},\\ &H^{(1)} = \lambda\cdot\frac{1}{2}c_{\rm SU(4)}\sum_{\boldsymbol{n}}:\tilde{\rho}^2(\boldsymbol{n}):, \end{align}\] with parameter \(\lambda=0.1\). Such a value is chosen to ensures the validity of first-order perturbation theory.

As specific benchmarks, we compare the charge radii of \(^{16}\)O and \(^{17}\)O calculated using the partial pinhole and the standard pinhole algorithms in Fig. 5. We can find that both algorithms give us almost identical charge radii. However, the partial pinhole algorithm maintains a significantly large phase factor, especially noticeable for \(^{17}\)O at projection time longer than \(0.4~{\rm MeV}^{-1}\). As a result, the charge radii from partial pinhole algorithm are always shown with much smaller statistical uncertainties. For heavier nuclei, such as \(^{40}\)Ca, the difference becomes even more pronounced. As presented in Fig. 6, the average phase factor in the standard pinhole algorithm drops dramatically with the increase of \(\tau\), a phase factor of less than \(0.01\) at a projection time of \(0.4~{\rm MeV}^{-1}\), while the partial pinhole algorithm achieves a phase factor as high as \(0.95\). As a results, the statistical uncertainties of the charge radius are greatly reduced by the partial pinhole algorithm. The advantage of the partial pinhole approach becomes even greater when employing the full N\(^3\)LO chiral interactions.

To further illustrate the effectiveness of the partial pinhole algorithm, Fig. 7 displays how charge radii depend on the number \(M\) of nucleons included in the sampled density operator \(\rho_M\), taking \(^{16}\)O as example. We can find that calculated charge radii are independent on the choice of \(M\). The average phase factors decrease gradually as \(M\) increases, approaching the smaller values characteristic of the standard pinhole method in the limit \(M = A\). This also explains how the partial pinhole algorithm mitigates the Monte Carlo sign problem encountered in nuclear radius calculation.

5.5 Imaginary time extrapolation↩︎

In Figs. 1 and 2 of the main text, we perform extrapolations to the \(\tau \to \infty\) in order to extract the ground-state charge and matter radii. For each observable \(O^{(i)}\), we fit the calculated values at finite projection time \(\tau\) using a double-exponential ansatz, \[\label{eq:fit95eq} O^{(i)} = O^{(i)}({\infty}) + O_1^{(i)}e^{-d\tau/2} + O_2^{(i)}e^{-d\tau},\tag{7}\] where \(O^{(i)}({\infty})\), \(O_1^{(i)}\), \(O_2^{(i)}\), and the decay constant \(d\) are fit parameters. The observables considered in the fit include both charge and matter radii, computed at leading order and at N\(^3\)LO. For each nucleus, all observables share the same decay parameter \(d\) to ensure consistency in the treatment of excited-state contaminations.

The function form in Eq. 7 is motivated by the assumption that at sufficiently large imaginary time, and the wave function is dominated by the ground state and first excited state. However, there is no strict a priori criterion for how large \(\tau\) must be for this assumption to hold. We therefore adopt a practical strategy in which we vary the lower bound \(\tau_c\) of the fit imaginary-time window to identify a stable region. Specifically, data points with \(\tau < \tau_c\) are the suppressed by assigning them an artificially large uncertainty (e.g., 0.05 fm), thereby reducing their influence on the fit. By scanning a range of \(\tau_c\) and monitoring the stability of the extrapolated quantities, we identify a plateau region where the results become insensitive to \(\tau_c\). This behavior signals a reliable window for extracting the ground-state properties.

The dependence of the extrapolated charge and matter radii on \(\tau_c\) is presented in Tables 2–5 for \(^{16}\)O, \(^{17}\)O, \(^{18}\)O, and \(^{20}\)O, respectively. In all cases, we observe stable plateaus in the extrapolated charge and matter radii as \(\tau_c\) increases. Once the plateau is reached, the fitted parameter \(O^{(i)}({\infty})\) in Eq. 7 exhibits a negligible variation. Consequently, the fitted curves shown in Figs. 1 and 2 of the main text remain virtually unchanged. This confirms the robustness of our extrapolated nuclear radii. In the discussion of the main text, we adopt \(\tau_c = 0.18~\mathrm{MeV}^{-1}\) for all extrapolations, as it lies within the plateau region for the investigated oxygen isotopes and ensures stable and consistent extrapolation results.

References↩︎