Ab initio short-range-correlation scaling factors from light to medium-mass nuclei

TL;DR

This paper uses ab initio nuclear theory to calculate short-range correlation scaling factors across light to medium-mass nuclei, linking high-energy scattering data with low-energy nuclear structure predictions.

Contribution

It introduces a method to predict SRC scaling factors from effective field theory, connecting high-energy scattering phenomena with ab initio nuclear calculations.

Findings

SRC scaling factors can be predicted from low-energy nuclear theory.

EMC effect and SRC factors have minimal isovector corrections.

Calculations extend to medium-mass nuclei.

Abstract

High-energy scattering processes, such as deep inelastic scattering (DIS) and quasielastic (QE) scattering provide a wealth of information about the structure of atomic nuclei. The remarkable discovery of the empirical linear relationship between the slope of the European Muon Collaboration (EMC) effect in DIS and the short-range-correlation (SRC) scaling factors $a_2$ in QE kinematics is naturally explained in terms of scale separation in effective field theory. This explanation has powerful consequences, allowing us to calculate and predict SRC scaling factors from ab initio low-energy nuclear theory. We present ab initio calculations of SRC scaling factors for a nucleus $A$ relative to the deuteron $a_2(A/d)$ and relative to $^3\rm He$ $a_2(A/^3\rm He)$ in light and medium-mass nuclei. Our framework further predicts that the EMC effect and SRC scaling factors have minimal or negligible isovector corrections.

Tables

Table 1. Table 1: Results for the SRC scaling factor a 2 ​ ( A / d ) subscript 𝑎 2 𝐴 𝑑 a_{2}(A/d) (upper table) and a 2 ​ ( A / \isotope ​ [ 3 ] ​ H ​ e ) subscript 𝑎 2 𝐴 \isotope delimited-[] 3 𝐻 𝑒 a_{2}(A/\isotope[3]{He}) (lower table) obtained via a linear fit to the Monte Carlo results for different nuclear interactions (see text for details). Both statistical (first) and systematic (second) uncertainties are reported in the parentheses. The latter include both the uncertainty coming from the fit of a 2 subscript 𝑎 2 a_{2} and the uncertainty associated with the truncation of the chiral expansion (for local chiral interactions). AV18 + UIX AV18 UIX \text{AV18}+\text{UIX} results are from GFMC calculations  [ 19 ] , while the other results are obtained using the AFDMC method. The last column shows the available experimental data from  [ 32 , 88 ] .

System	$\text{N}\text{<sup><span>2</span></sup>}\text{LO}E\tau$		$\text{N}\text{<sup><span>2</span></sup>}\text{LO}E\mathrm{𝟙}$		$\text{AV18}+\text{UIX}$	${\text{AV4}}^{\prime }+{\text{UIX}}_{\text{c}}$	Exp
	$R_0=1.0\text{ fm}$	$R_0=1.2\text{ fm}$	$R_0=1.0\text{ fm}$	$R_0=1.2\text{ fm}$
\isotope[3]H	$2.18\left(2\right)\left(27\right)$	$2.27\left(2\right)\left(8\right)$	$2.15\left(2\right)\left(28\right)$	$2.46\left(2\right)\left(8\right)$	$2.12\left(1\right)\left(22\right)$	$-$	$-$
\isotope[3]He	$2.13\left(2\right)\left(27\right)$	$2.23\left(1\right)\left(8\right)$	$2.10\left(2\right)\left(28\right)$	$2.38\left(2\right)\left(8\right)$	$2.07\left(1\right)\left(21\right)$	$2.09\left(2\right)\left(21\right)$	$2.13\left(4\right)$
\isotope[4]He	$3.76\left(2\right)\left(46\right)$	$4.16\left(1\right)\left(27\right)$	$3.77\left(2\right)\left(46\right)$	$5.31\left(2\right)\left(27\right)$	$4.00\left(1\right)\left(40\right)$	$3.83\left(2\right)\left(39\right)$	$3.60\left(10\right)$
\isotope[6]He	$3.23\left(2\right)\left(30\right)$	$3.51\left(1\right)\left(41\right)$	$3.14\left(1\right)\left(30\right)$	$4.04\left(2\right)\left(41\right)$	$-$	$-$	$-$
\isotope[6]Li	$3.46\left(1\right)\left(29\right)$	$3.78\left(1\right)\left(39\right)$	$3.33\left(1\right)\left(29\right)$	$4.18\left(2\right)\left(39\right)$	$-$	$-$	$-$
\isotope[12]C	$4.67\left(2\right)\left(82\right)$	$-$	$-$	$-$	$-$	$-$	$4.75\left(16\right)$
\isotope[16]O	$4.30\left(1\right)\left(1.75\right)$	$8.55\left(1\right)\left(4.40\right)$	$4.02\left(1\right)\left(1.75\right)$	$5.47\left(1\right)\left(4.40\right)$	$-$	$4.62\left(2\right)\left(47\right)$	$-$
\isotope[40]Ca	$-$	$-$	$-$	$-$	$-$	$5.15\left(2\right)\left(67\right)$	$-$
\isotope[63]Cu	$-$	$-$	$-$	$-$	$-$	$-$	$5.21\left(20\right)$
\isotope[197]Au	$-$	$-$	$-$	$-$	$-$	$-$	$5.16\left(22\right)$
\isotope[4]He	$1.76\left(2\right)\left(4\right)$	$1.87\left(1\right)\left(4\right)$	$1.80\left(2\right)\left(4\right)$	$2.23\left(2\right)\left(9\right)$	$1.83\left(2\right)\left(19\right)$	$1.94\left(1\right)\left(20\right)$	$1.93\left(2\right)\left(14\right)$
\isotope[6]He	$1.52\left(2\right)\left(5\right)$	$1.58\left(1\right)\left(6\right)$	$1.50\left(2\right)\left(4\right)$	$1.70\left(2\right)\left(7\right)$	$-$	$-$	$-$
\isotope[6]Li	$1.62\left(2\right)\left(7\right)$	$1.70\left(1\right)\left(6\right)$	$1.59\left(2\right)\left(6\right)$	$1.76\left(2\right)\left(6\right)$	$-$	$-$	$-$
\isotope[12]C	$2.19\left(2\right)\left(18\right)$	$-$	$-$	$-$	$-$	$-$	$2.41\left(2\right)\left(17\right)$
\isotope[16]O	$2.02\left(2\right)\left(33\right)$	$3.91\left(2\right)\left(1.01\right)$	$1.91\left(2\right)\left(33\right)$	$2.27\left(2\right)\left(1.01\right)$	$-$	$2.27\left(2\right)\left(23\right)$	$-$
\isotope[40]Ca	$-$	$-$	$-$	$-$	$-$	$2.55\left(2\right)\left(33\right)$	$-$
\isotope[56]Fe	$-$	$-$	$-$	$-$	$-$	$-$	$2.83\left(3\right)\left(18\right)$

Table 2. Table 2: Binding energies (in MeV) and charge radii (in fm) for A = 4 , 16 , 40 𝐴 4 16 40 A=4,16,40 with the AV4 ′ + UIX c superscript AV4 ′ subscript UIX c \text{AV4}^{\prime}+\text{UIX}_{\text{c}} potential. Energy results are from the AFDMC unconstrained evolution  [ 17 ] . Experimental results are shown for comparison.

\isotope[A]Z	$E_{\text{AFDMC}}$	$E_{\text{Exp}}$	$r_{\text{ch}}^{\text{AFDMC}}$	$r_{\text{ch}}^{\text{Exp}}$
\isotope[4]He	$-26.00\left(2\right)$	$-28.296$	$1.74\left(1\right)$	$1.680\left(4\right)$ [89]
\isotope[16]O	$-113\left(2\right)$	$-127.619$	$2.61\left(6\right)$	$2.699\left(5\right)$ [90]
\isotope[40]Ca	$-321\left(3\right)$	$-342.052$	$3.25\left(8\right)$	$3.478\left(2\right)$ [90]

Equations

$${\rho }_{2,\text{mathbbm}1}(A,r)\propto ⟨\Psi |\sum _{i<j}^A\delta (r-r_{ij})|\Psi ⟩,$$

$$\frac{{\mathrm{d}}^2\sigma }{\mathrm{d}x\mathrm{d}{Q}^2}=\frac{2\pi {\alpha }^2}{x{Q}^4}{F}_2(x,Q^2)\left[1+(1-y{)}^2\right].$$

$$R_{\text{EMC}}(A,x)\equiv \frac{2F_2^A(x,Q^2)}{A{F}_2^d(x,Q^2)}$$

$$\sigma (x,Q^2)=\sum _{j=2}^A\frac{1}{j}{a}_j(A){\sigma }_j(x,Q^2),$$

$$a_2(A/d)\equiv {\frac{2{\sigma }_A}{A{\sigma }_d}|}_{1.5<x<2},$$

$$F_2^A(x,Q^2)/A\simeq F_2^N(x,Q^2)+g_2(A,\Lambda )f_2(x,Q^2,\Lambda ),$$

$$R_{\mathrm{EMC}}(A,x)-1\simeq C(x)\left[a_2(A)-1\right],$$

$$C(x)$$

$$a_2(A)$$

$$R_{\text{EMC}}(A,1<x<2)\simeq a_2(A),$$

$$g_2(A,\Lambda )\equiv \frac{1}{2A}⟨A|:(N^{†}N{)}^2:|A{⟩}_{\Lambda },$$

$$\begin{array}{rl}{F}_2^A(x,Q^2)& \simeq Z{F}_2^p(x,Q^2)+N{F}_2^n(x,Q^2)+A{g}_2(A,\Lambda )f_2(x,Q^2,\Lambda )+\cdots ,\end{array}$$

$$⟨x^n{⟩}_A(Q)={\int }_{-A}^A\mathrm{d}x{x}^n{q}_A(x,Q),$$

$$⟨A;p|{\mathcal{O}}^{{\mu }_0\cdots {\mu }_n}|A;p⟩=2⟨x^n{⟩}_A(Q)p^{({\mu }_0}\cdots p^{{\mu }_n)},$$

$${\mathcal{O}}^{{\mu }_0\cdots {\mu }_n}=\overline{q}{\gamma }^{({\mu }_0}\mathrm{i}{D}^{{\mu }_1}\cdots \mathrm{i}{D}^{{\mu }_n)}q,$$

$$\begin{array}{rl}{\mathcal{O}}^{{\mu }_0\cdots {\mu }_n}\to & :2⟨x^n{⟩}_N{M}_N^{n+1}{v}^{({\mu }_0}\cdots v^{{\mu }_n)}{N}^{†}N\left[1+{\alpha }_n{N}^{†}N\right]\\ & +⟨x^n{⟩}_{\pi }{\pi }^{\alpha }\mathrm{i}{\partial }^{({\mu }_0}\cdots \mathrm{i}{\partial }^{{\mu }_n)}{\pi }^{\alpha }+\cdots :,\end{array}$$

$$\displaystyle\begin{split}\langle x^{n}\rangle_{A}(Q)&=\langle x^{n}\rangle_{N}(Q)\Bigl{[}A+\alpha_{n}(\Lambda,Q)\langle A|:\mathrel{\mkern 2.0mu(N^{\dagger}N)^{2}\mkern 2.0mu}:|A\rangle_{\Lambda}\Bigr{]}\,\\ &+\delta_{n=1}\text{ term}\,,\end{split}$$

$$q_A(x,Q)/A\simeq q_N(x,Q)+g_2(A,\Lambda ){\tilde{q}}_2(x,Q,\Lambda ),$$

$${\mathcal{O}}_3^{{\mu }_0\cdots {\mu }_n}=\overline{q}{\tau }_3{\gamma }^{({\mu }_0}\mathrm{i}{D}^{{\mu }_1}\cdots \mathrm{i}{D}^{{\mu }_n)}q,$$

$$\begin{array}{rl}{\mathcal{O}}_3^{{\mu }_0\cdots {\mu }_n}\to & :2⟨x^n{⟩}_{N,3}{M}_N^{n+1}{v}^{({\mu }_0}\cdots v^{{\mu }_n)}{N}^{†}{\tau }_{3}N\left[1+{\gamma }_n{N}^{†}N\right]\\ & +2{\delta }_n{M}_N^{n+1}{N}^{†}{S}^{({\mu }_0}{v}^{{\mu }_1}\cdots v^{{\mu }_n)}{\pi }^{\alpha }[{\tau }^{\alpha },{\tau }_3]N\\ & +⟨x^n{⟩}_{\pi ,3}\mathrm{i}{ϵ}^{3\alpha \beta }{\pi }^{\alpha }\mathrm{i}{\partial }^{({\mu }_0}\cdots \mathrm{i}{\partial }^{{\mu }_n)}{\pi }^{\beta }+\cdots :.\end{array}$$

$$H=-\frac{{\hslash }^2}{2M_N}\sum _i^A{\nabla }_i^2+\sum _{i<j}^A{v}_{ij}+\sum _{i<j<k}^A{V}_{ijk},$$

$$v_{ij}=\sum _{p=1}^7{v}_p(r_{ij}){\mathcal{O}}_{ij}^p,$$

$$\displaystyle\mathcal{O}_{ij}^{p=1,\ldots,7}=\Big{\{}$$

$$V_{ijk}=V_{2\pi }^P+V_{2\pi }^S+V_D+V_E.$$

$$V_{D2}=\frac{g_A{c}_D{m}_{\pi }^2}{96\pi {\Lambda }_{\chi }{F}_{\pi }^4}\sum _{i<j<k}\sum _{\text{cyc}}{\mathbf{\tau }}_i\cdot {\mathbf{\tau }}_k\left[X_{ik}({\mathbf{r}}_{ik})-\frac{4\pi }{m_{\pi }^2}{\mathbf{\sigma }}_i\cdot {\mathbf{\sigma }}_k{\delta }_{R_{3N}}(r_{ik})\right][{\delta }_{R_{3N}}(r_{ij})+{\delta }_{R_{3N}}(r_{kj})],$$

$$V_{E\tau }$$

$$V_{E\text{mathbbm}1}$$

$$|\Psi (\tau )⟩$$

$$\underset{\tau \to \infty }{\lim }|\Psi (\tau )⟩$$

$$\displaystyle|\Psi_{T}\rangle=\big{[}F_{C}+F_{2}+F_{3}\big{]}|\Phi\rangle_{J^{\pi}T}\,,$$

$$⟨\mathcal{O}⟩=\frac{1}{\mathcal{N}}\sum _{i=1}^{\mathcal{N}}\frac{⟨R_i{S}_i|\mathcal{O}|{\Psi }_T⟩}{⟨R_i{S}_i|{\Psi }_T⟩},$$

$$⟨\mathcal{O}⟩\approx 2\frac{⟨{\Psi }_T|\mathcal{O}|\Psi (\tau )⟩}{⟨{\Psi }_T|\Psi (\tau )⟩}-\frac{⟨{\Psi }_T|\mathcal{O}|{\Psi }_T⟩}{⟨{\Psi }_T|{\Psi }_T⟩}.$$

Full text

Keywords: short-range correlations, short-range-correlation scaling factors, EMC effect, quantum Monte Carlo

1 Introduction

The accurate description and prediction of the structure and behavior of atomic nuclei remains an important problem in physics. In spite of decades of experimental, theoretical, and computational research and the fact that quantum chromodynamics (QCD) is widely understood to provide the underlying field theoretic description, the strong interaction between protons and neutrons still surprises us with its subtlety. Because of the nature of the nonabelian gauge interactions, QCD resists perturbative treatments at low energies [1, 2]. Explicit solutions at these energies are possible via the computational framework of lattice QCD, wherein observables are calculated directly in QCD but on a finite Euclidean space-time lattice [3, 4, 5].

While lattice QCD promises a fundamental explanation of nuclear physics phenomena, the computational difficulties it faces grow rapidly with the system size, so that current simulations are limited to few-nucleon systems [3, 4, 6, 5, 7, 8, 9]. This means that for many interesting nuclear systems, other methods are needed at present. Low-energy ab initio nuclear theory, working with protons and neutrons as degrees of freedom and fixing the parameters of the theory with results from either experimental data or lattice QCD, naturally fills this role, and the field has made significant progress in recent years in terms of working with systematically improvable Hamiltonians derived from chiral effective field theory (EFT) and in terms of the size of the nuclear systems that can be accurately handled [10, 11, 12, 13, 14, 15, 16, 17].

In particular, in recent years advances made in accurate quantum Monte Carlo (QMC) methods and their combination with interactions derived from chiral EFT has provided many new insights in low-energy nuclear theory (see [18] for a review). One such insight to arise from the use of QMC methods with EFT techniques is that, while commonly calculated two-body central densities

[TABLE]

with $r_{ij}$ the internucleon separation for a nucleus with $A$ nucleons, are scheme and scale dependent, their ratios are largely scheme and scale independent [19]. Moreover, these ratios at small internucleon separation correspond to short-range-correlation (SRC) observables in quasielastic (QE) lepton-nucleus scattering: in short, $\lim_{r\to 0}\frac{\rho_{2,\mathbbm{1}}(A,r)}{\rho_{2,\mathbbm{1}}(d,r)}\propto a_{2}(A/d)$, where $a_{2}$ is the so-called SRC scaling factor, and $d$ stands for the deuteron [19].

In this paper, we exploit this unique convergence of advances in QMC methods and EFT to confirm this relationship in light nuclei up to \isotope[12]C by comparing with existing experimental data. We then make predictions for several light systems (\isotope[6]He, \isotope[6]Li, and \isotope[16]O) and for the medium-mass nucleus \isotope[40]Ca, which could be tested in existing and near-term future experimental facilities.

The structure of this article is as follows. In what remains of Section 1 we present some background (Section 1.1), the main EFT arguments (Section 1.2), and details on the EFT power counting (Section 1.3). In Section 2 we briefly discuss the nuclear Hamiltonian and our QMC methods. In Section 3 we present our main results, discussing how best to extract the SRC scaling factors from our QMC results. Finally, in Section 4 we summarize our results and provide an outlook for this novel framework.

1.1 Background

Deep-inelastic scattering (DIS) of leptons on nuclear targets has been one of the most valuable experimental tools for learning about the structure of nucleons and nuclei. In DIS, a highly energetic ($Q^{2}\sim 5\,\text{GeV}^{2}$) leptonic probe with four momentum $p$ is scattered from a hadronic target with four momentum $P$, transferring four momentum $q$ to the struck quark, see Fig. 1. The cross section can be written in terms of the dimensionless Bjorken $x\equiv\frac{Q^{2}}{2P\cdot q}$, with $Q^{2}=-q^{2}$, because $q$ is spacelike, the dimensionless variable $y\equiv\frac{P\cdot q}{P\cdot p}$, and the structure function $F_{2}(x,Q^{2})$:

[TABLE]

In 1983, the European Muon Collaboration (EMC) announced their unexpected results for the measurement of the structure functions $F_{2}^{A}(x,Q^{2})$ in leptonic DIS for iron (\isotope[56]Fe) and deuterium [20]. The surprise came because, given that the typical binding energy per nucleon is so small (i.e. $\lesssim 1\%$) compared to the nucleon mass and the energy transfer in the DIS process, the expectation was that the cross section would have only trivial dependence on the nuclear target. Instead, in the region $0.2\lesssim x\lesssim 0.7$, the ratio $2F_{2}^{\isotope{Fe}}/AF_{2}^{d}$ was observed to fall off linearly to a significant reduction of $\sim 10\%$ at $x\sim 0.7$. This reduction in the ratio

[TABLE]

has come to be known as the EMC effect. Since then, significant experimental and theoretical effort has been invested to understand this effect (see [21, 22, 23, 24] for reviews).

As part of this effort to further understand the implications of the EMC effect, more experiments were carried out for smaller values of $x$, at different $Q^{2}$ [25], and for various nuclei [26], and more recently in QE scattering at higher $x$, $1\lesssim x\lesssim 2$ [27]. The picture that emerges is that the ratio of nuclear structure functions $R_{\text{EMC}}(A,x)$ has very little $Q^{2}$ dependence, and for isoscalar nuclei, the $A$ and $x$ dependence of $R_{\text{EMC}}-1$ factorizes. That is, the shape of the deviation from unity of the ratio $R_{\text{EMC}}(A,x)$ is independent of $A$, while the maximum magnitude only depends on $A$. Figure 2 shows an example of the universal $x$ dependence of the data. The different regions are labeled with the favored explanation for the behavior of the ratio in that region (see [22] for a more detailed explanation of the history of attempts at explaining the EMC effect).

In this work, we are interested in the EMC region $(0.35<x<0.7)$ and the SRC region $(1<x<2)$. The strength of the effect in the former region is usually characterized by the slope $|{\rm d}R_{\text{EMC}}/{\rm d}x|$ (see again Fig. 2), which ranges from $\sim 0.07$ in \isotope[3]He up to $\sim 0.5$ in \isotope[108]Ag, showing a trend towards saturation as the mass number $A$ increases. In the latter region, based on an impulse-approximation argument, Frankfurt et al [30] cast the inclusive cross section as

[TABLE]

where the $a_{j}(A)$ are proportional to the probabilities to find a nucleon in a $j$-nucleon SRC, and $\sigma_{j}(x,Q^{2})=0$ for $x>j$. This framework correctly predicted the scaling behavior ($x$ and $Q^{2}$ independence) in the ratio of cross sections:

[TABLE]

where $a_{2}$ is the SRC scaling factor introduced earlier. (Note that Fermi motion pushes the onset of the plateau from $x\sim 1$ to $x\sim 1.5$). In some of the more recent experiments at Thomas Jefferson National Accelerator Facility (Jefferson Lab) [27], these plateaus have been observed for nuclear targets from \isotope[3]He to \isotope[197]Au (see also Fig. 2).

Recently, a fascinating empirical discovery was made: the slope of the EMC effect in the EMC region is linearly correlated with the SRC scaling factor [31, 32], see Fig. 3. This remarkable result has motivated a series of experiments attempting to further understand this phenomenon, as well as many theoretical proposals. As discussed in [19], the physics behind this correlation is naturally explained in the EFT approach used here.

1.2 Effective field theory

Effective field theory is a model-independent approach that relies on the symmetries and the separation of scales in a given system. Effective field theory has been successfully applied to many aspects of meson [33], single- [34], and multi-nucleon systems [35, 36, 37, 38, 39, 40]. In particular, chiral EFT has been applied to parton distribution functions (PDFs) in the meson, single-nucleon [41, 42, 43, 44, 45, 46, 47, 48, 49], and multi-nucleon sectors [50, 51], as well as to other light-cone dominated observables [52, 53, 54, 55, 56, 57].

In 2005, using EFT, Chen and Detmold [50] found that, up to higher order corrections, the $F_{2}$ structure function of an isoscalar nucleus has the form

[TABLE]

where $F_{2}^{N}$ is the isoscalar combination of the nucleon structure function, which receives the nuclear modification from the second term in which the $x$ and $A$ dependence factorizes. The $A$ dependence comes from momenta smaller than the ultraviolet momentum cutoff of the EFT $\Lambda\sim 0.5\,\rm GeV$, while the $x$ dependence comes from scales larger than $\Lambda$.

An immediate consequence of Eq. 6 is that

[TABLE]

with the $x$ and $A$ dependence factorized, and

[TABLE]

The deviation of $R_{\text{EMC}}$ from unity in Eq. 7 means that the nuclear modification to the structure functions has a universal shape ($x$ dependence), while its maximum magnitude depends only on $A$ [50]. This feature describes experimental data with $x<1$ for many nuclei, ranging from \isotopeHe to \isotopePb very well [58, 59].

Because $F_{2}^{A}(x)$ has support for $0<x<A$, if DIS experiments were carried out at $1<x<2$, where $F_{2}^{N}(x)=0$ but $F_{2}^{d}(x)\neq 0$, then Eqs. 7 and 8 yields

[TABLE]

which is an $x$-independent plateau. Experimentally, the measurements at $x>1$ are performed not in the DIS region, but in the QE region at lower $Q^{2}$ because of the larger associated rate. Generalizing the analysis to the QE region by including all the higher twist effects does not change the plateau value of Eq. 10 [19]. The plateau is observed experimentally at $1.5<x<2$, possibly because Fermi motion, which is a higher-order effect in the EFT, extends the contribution of the single-nucleon PDF to $x$ slightly above 1, so that the onset of the plateau is also pushed to larger $x$.

From Eqs. 7 and 10, the observed linear relation between $-{\rm d}R_{\text{EMC}}/{\rm d}x$ and the SRC scaling factor $a_{2}(A)$ is easily obtained. Equation 10 demands that the scaling factor, which comes from the ratio of two cross sections, be independent of the cutoff $\Lambda$. Therefore, the $\Lambda$ dependence on the right-hand side of Eq. 9 should cancel. This provides a nontrivial test of EFT, because it implies that, although $g_{2}(A,\Lambda)$ depends on the renormalization scheme and scale ($\Lambda$) of the EFT, $a_{2}$ is scheme and scale independent. This occurs if the $\Lambda$ and $A$ dependence factorize in $g_{2}$, which is defined as

[TABLE]

where $N$ is the nucleon field and $:\mathrel{\mkern 2.0mu\ \cdots\ \mkern 2.0mu}:$ indicates normal ordering of the enclosed operators with respect to the vacuum state.

The above analysis is for isoscalar operators. Including isovector corrections, one has

[TABLE]

with $N$ $(Z)$ the number of neutrons (protons) in the nucleus. The isovector counterpart of the $g_{2}$ term is neglected because it is $\mathcal{O}((N-Z)/AN_{c})$ smaller than $g_{2}$, with the number of colors $N_{c}=3$. This implies that, even with isovector corrections, the SRC plateaus still exist, and the plateau values of $a_{2}$ remain unchanged. Also, for the EMC effect, recent experimental results including nonisoscalar nuclei are well described by Eq. 12 [60].

1.3 EFT power counting

In DIS, the structure functions $F_{2}^{A}(x,Q^{2})$ can be expressed in terms of nuclear PDFs $q_{i}^{A}(x,Q)$ as $F_{2}^{A}(x,Q^{2})=\sum_{i}Q_{i}^{2}xq_{i}^{A}(x,Q)$, where the sum runs over quarks and antiquarks of flavor $i$ with charge $\pm Q_{i}$ in a nucleus $A$. In what follows, we first focus on the power counting for isoscalar PDFs, $q_{A}=q_{A,0}=q_{A,u}+q_{A,d}$, then we discuss the isovector correction from $q_{A,3}=q_{A,u}-q_{A,d}$. The dominant (leading-twist) PDFs are determined by target matrix elements of bilocal light-cone operators. Applying the operator product expansion, the Mellin moments of the PDFs,

[TABLE]

are determined by matrix elements of local operators,

[TABLE]

with

[TABLE]

where $(\cdots)$ indicates that the enclosed indices have been symmetrized and made traceless, D^{\mu}\equiv(\vec{D}^{\mu}-\mathchoice{\mkern 3.0mu\reflectbox{\displaystyle\vec{\reflectbox{$\mkern-3.0mu\displaystyle D\mkern 3.0mu$}}}\mkern-3.0mu}{\mkern 3.0mu\reflectbox{\textstyle\vec{\reflectbox{$\mkern-3.0mu\textstyle D\mkern 3.0mu$}}}\mkern-3.0mu}{\mkern 2.0mu\reflectbox{\scriptstyle\vec{\reflectbox{$\mkern-2.0mu\scriptstyle D\mkern 2.0mu$}}}\mkern-2.0mu}{\mkern 2.0mu\reflectbox{\scriptscriptstyle\vec{\reflectbox{$\mkern-2.0mu\scriptscriptstyle D\mkern 2.0mu$}}}\mkern-2.0mu}^{\mu})/2 is the covariant derivative, and a sum over flavors $q=u,d$ is implied. The negative $x$ distribution is the antiquark distribution: $q_{A}(-x)=-\bar{q}_{A}(x)$.

In nuclear matrix elements of these operators, there are other relevant momentum scales below the hard scattering scale $Q$: $\Lambda\sim 0.5\,\rm GeV$ is the range of validity of the EFT, and $P\sim m_{\pi}$ is a typical momentum inside the nucleus ($m_{\pi}$ is the pion mass). These scales satisfy $Q\gg\Lambda\gg P$, and the ratio $\Lambda/Q$ is the small expansion parameter in the twist expansion, while the ratio $\epsilon\sim P/\Lambda\sim 0.2\text{--}0.3$ is the small expansion parameter for the chiral expansion.

In EFT, each of the QCD operators is matched to a sum of all possible hadronic operators of the same symmetries at the scale $\Lambda$ [50]

[TABLE]

where $\pi$ ($N$) is the pion (nucleon) field, $v$ is the nucleon four-velocity, and $\langle x^{n}\rangle_{N\,(\pi)}$ is the $n$th moment of the isoscalar quark PDF in a free nucleon (pion). There are an infinite number of terms on the right-hand side of Eq. 16, whose importance will be estimated by power counting. The $\langle x^{n}\rangle_{N\,(\pi)}$ terms are one-body operators acting on a single hadron, whose prefactors can be determined by taking the nucleon (pion) matrix element of Eq. 16. The $\alpha_{n}$ terms are two-body operators. Here we have only kept the SU(4) (spin and isospin) singlet two-body operator $\propto\left(N^{\dagger}N\right)^{2}$ and neglected the SU(4) nonsinglet operator $\propto(N^{\dagger}\bm{\sigma}N)^{2}-(N^{\dagger}\bm{\tau}N)^{2}$, which changes sign when interchanging the spin ($\bm{\sigma}$) and isospin ($\bm{\tau}$) matrices [61]. The latter operator has an additional $\mathcal{O}(1/N_{c}^{2})\sim 0.1$ suppression in its prefactor [62]. We also replace the nucleon velocity by the nucleus velocity and include the correction $\mathrm{i}\partial_{0}/M_{N}$ at higher orders.

In Weinberg’s power counting scheme, the typical nucleon momenta $|{\bf q}|$ are counted as $\mathcal{O}(\epsilon)$, while their energies $q^{0}$ are $\mathcal{O}(\epsilon^{2})$. Two-nucleon contact operators $(N^{\dagger}N)^{2}$ are counted as $\mathcal{O}(\epsilon^{0})$, while the three-body contact operator $(N^{\dagger}N)^{3}$ is counted as $\mathcal{O}(\epsilon^{3})$, both according to their mass dimension. We will focus on the twist-2 operators with all $\mu_{i}=0$ in Eq. 16. Because $v^{0}=1$, the $v^{(\mu_{0}}\cdots v^{\mu_{n})}(N^{\dagger}N)$ operator is $\mathcal{O}(\epsilon^{-3})$ and $v^{(\mu_{0}}\cdots v^{\mu_{n})}(N^{\dagger}N)^{2}$ is $\mathcal{O}(\epsilon^{0})$. The one-derivative operator $N^{\dagger}\partial^{(\mu_{0}}v^{\mu_{1}}\cdots v^{\mu_{n})}N$ is $\mathcal{O}(\epsilon^{-1})$, but its net effect is to shift the value of $p^{0}$ on the right-hand side of Eq. 14 from $AM_{N}$ to $M_{A}$. This can be seen from the special case of $n=0$. The vector current operator $\mathcal{O}^{\mu_{0}}$ is matched to the operator $2M_{N}N^{\dagger}(v^{\mu_{0}}+\mathrm{i}\partial^{\mu_{0}}/M_{A})N$. The nuclear matrix element of the first term yields $2AM_{N}$. The relative coefficient between the two terms are fixed by reparametrization invariance [63], and the nuclear matrix element of the sum yields $2M_{A}$.

The two-derivative operator given by $N^{\dagger}\partial^{(\mu_{0}}\partial^{\mu_{1}}v^{\mu_{2}}\cdots v^{\mu_{n})}N$, (again with $\mu_{i}=0$ for all $i$) is $\mathcal{O}(\epsilon)$, and it can cause $q_{N}(x)$ or $F^{N}_{2}(x)$ to “spill” into $x>1$. This is related to Fermi motion. Although it is higher order than the two-body operator, if $f_{2}(x)$ of Eq. 6 is very small when $x$ is just above one, then the Fermi-motion effect could be dominant and explain why the $a_{2}$ plateau only sets in at $x\gtrsim 1.5$. It is important to note that, in the EFT approach, off-shell effects that enter through Fermi motion can be absorbed into many-body operators through a field redefinition [64, 65]. Therefore the separation between “Fermi motion” and “two-body effects” is meaningful only after the theory is clearly specified.

The pion one-body operator $\pi^{a}\mathrm{i}\partial^{(\mu_{0}}\cdots\mathrm{i}\partial^{\mu_{n})}\pi^{a}$ inserted in the one-pion-exchange diagram contributes at $\mathcal{O}(\epsilon^{n-1})$. Because $\langle x^{n}\rangle_{\pi}=0$ for even $n$ due to charge conjugation symmetry, the $n=1$ pion operator enters at $\mathcal{O}(\epsilon^{0})$, but for higher $n$ the contributions either vanish or are higher order compared with the other operators in Eq. 16. This means that, at $\mathcal{O}(\epsilon^{0})$, the pion contribution to Eq. 6 is proportional to $\delta(x)/x$ and breaks the factorization of the $x$ and $A$ dependence of $F^{2}_{A}$ of Eq. 6, but only at $x=0$.

All the other operators in the matching are found to be higher order than $\epsilon^{0}$ in this power counting. Using nucleon number conservation, $\langle A|:\mathrel{\mkern 2.0muN^{\dagger}N\mkern 2.0mu}:|A\rangle=A$, the isoscalar nuclear matrix element of Eq. 16 is

[TABLE]

where $\alpha_{n}$ is $A$ independent but $\Lambda$ dependent, and is completely determined by the two-nucleon system. After an inverse Mellin transform, except at $x=0$ as explained above, the isoscalar PDFs satisfy

[TABLE]

where $\tilde{q}_{2}(x,Q,\Lambda)$ is an unknown function independent of $A$ whose Mellin moments are determined by the low-energy constants $\alpha_{n}$. This result also holds at the level of the structure function, which leads to Eq. 6.

The isovector operator

[TABLE]

is matched to hadronic operators as

[TABLE]

The $\langle x^{n}\rangle_{N,3}$ term is $\mathcal{O}(\epsilon^{-3})$. The $\gamma_{n}$ term is $\mathcal{O}(\epsilon^{0})$, like the $\alpha_{n}$ operator of Eq. 16, but it has an additional $1/N_{c}$ suppression in its prefactor [62] and an $(N-Z)/A$ suppression in its nuclear matrix element compared with the $\alpha_{n}$ term and, hence can be neglected. $S^{\mu}$ is the nucleon spin vector. Using $\pi^{\alpha}[\tau^{\alpha},\tau_{3}]\propto(\pi^{+}\tau_{+}-\pi^{-}\tau_{-})$, the $\delta_{n}$ term involves a charged pion exchange, which can only happen between $np$ states in two-nucleon systems. However, $\tau_{3}$ for $np$ states (which have isospin zero) vanishes, therefore, there is no net two-nucleon contribution from this term. The $\langle x^{n}\rangle_{\pi,3}$ term contributes at $\mathcal{O}(\epsilon^{n-1})$. However, $\langle x^{n}\rangle_{\pi,3}$ vanishes for odd $n$ by charge conjugation. The $\langle x^{0}\rangle_{\pi,3}$ term is the isospin charge, which is protected from nuclear modifications. The other terms $\langle x^{n\geq 2}\rangle_{\pi,3}$ are $\mathcal{O}(\epsilon)$ and higher and can be neglected. The leading three-body operator $v^{(\mu_{0}}\cdots v^{\mu_{n})}(N^{\dagger}N)^{2}N^{\dagger}\tau_{3}N$ is $\mathcal{O}(\epsilon^{3})$, which can also be neglected.

We remark that in the large $N_{c}$ limit, the nucleon and delta resonances are degenerate, hence one should explicitly include the deltas in the $1/N_{c}$ expansion. In the real world with $N_{c}=3$, the mass difference between delta and nucleon $\Delta m$ is much larger than the typical Fermi energy $E_{F}$ in a nucleus. Therefore, one can choose to integrate out the delta degrees of freedom, as done in this work. The effect is that the nucleon operators studied here will receive $\mathcal{O}(E_{F}/\Delta m)$ corrections, but their $N_{c}$ scalings remain the same.

In summary, up to $\mathcal{O}(\epsilon^{0})$, only the one-body operator $\langle x^{n}\rangle_{N,3}$ contributes to isovector corrections. Therefore, the nuclear effects are dominated by the isoscalar PDF contributions, while the isovector PDFs are relatively unaltered by the nuclear environment, leading to Eq. 12.

2 Hamiltonian and Quantum Monte Carlo Methods

In ab initio methods, nuclei are treated as a collection of point-like particles of mass $M_{N}$ interacting via two- and three-body potentials according to the nonrelativistic Hamiltonian

[TABLE]

where the two-body interaction $v_{ij}$ also includes the Coulomb force.

In this work, we adopt the local chiral nucleon-nucleon ($N\!N$) interactions at next-to-next-to-leading order (N2LO) in Weinberg counting of [66, 67], with coordinate-space cutoffs $R_{0}=1.0\,\rm fm$ and $R_{0}=1.2\,\rm fm$. Such interactions include long-range pion-exchange contributions, determined by pion-nucleon couplings, and shorter-range contributions, defined by low-energy couplings (LECs) that are fit to reproduce $N\!N$ scattering data. The local chiral $N\!N$ potentials are written in coordinate space as a sum of spin/isospin operators

[TABLE]

with

[TABLE]

where $r_{ij}=|{\bf r}_{i}-{\bf r}_{j}|$ is the $N\!N$ relative distance, $S_{ij}=3\,\bm{\sigma}_{i}\cdot\hat{{\bf r}}_{ij}\,\bm{\sigma}_{j}\cdot\hat{{\bf r}}_{ij}-\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}$ is the tensor operator, and ${\bf L}=({\bf r}_{i}-{\bf r}_{j})\times(\bm{\nabla}_{i}-\bm{\nabla}_{j})/2\mathrm{i}$ and ${\bf S}=(\bm{\sigma}_{i}+\bm{\sigma}_{j})/2$ are the relative angular momentum and the total spin of the pair $ij$, respectively.

At N2LO, in addition to the $N\!N$ interactions specified above, three-nucleon ($3N$) interactions enter [68, 69, 16, 17], see also [70, 71] for earlier formulations, often used as non-local interactions in momentum space. The employed $3N$ forces include two-pion-exchange (TPE) contributions in $P$ and $S$ waves, plus shorter-range components parametrized by two contact terms, usually referred to as $V_{D}$ and $V_{E}$:

[TABLE]

The TPE components are characterized by the LECs $c_{1}$, $c_{3}$, and $c_{4}$ from the pion-nucleon sector. The LECs of the contact terms, $c_{D}$ and $c_{E}$, have been fit to the $\alpha$ particle binding energy and to the spin-orbit splitting in the neutron-$\alpha$ $P$-wave phase shifts [69, 17]. We employ the form

[TABLE]

for $V_{D}$, and we consider two choices for $V_{E}$, namely $E\tau$ and $E\mathbbm{1}$:

[TABLE]

where $g_{A}$ is the axial vector coupling constant, $m_{\pi}$ is the pion mass, $\Lambda_{\chi}=700\,\text{MeV}$, $F_{\pi}$ is the pion decay constant, $X_{ij}({\bf r}_{ij})=[S_{ij}({\bf r}_{ij})T(r_{ij})+\bm{\sigma}_{i}\cdot\bm{\sigma}_{j}]Y(r_{ij})$ is the coordinate-space pion propagator, with the tensor and Yukawa functions defined as $T(r)=1+3/(m_{\pi}r)+3/(m_{\pi}r)^{2}$ and $Y(r)=\mathrm{e}^{-m_{\pi}r}/r$, respectively, and $\delta_{R_{3N}}(r)=\frac{\mathrm{e}^{-(r/R_{3N})^{4}}}{\pi\Gamma(3/4)R_{3N}^{3}}$ is a smeared-out delta function with $3N$ coordinate-space cutoff $R_{3N}$. We take this $3N$ cutoff equal to the $N\!N$ cutoff $R_{3N}=R_{0}$. The notation $\sum_{\text{cyc}}$ indicates a cyclic summation over the indices $\{ijk\}$. See [69, 16, 17] for more details including values for $c_{D}$ and $c_{E}$.

The operator structure of the employed local chiral interactions is suited for QMC calculations. QMC methods are a family of ab initio many-body techniques that allow one to solve the many-body Schrödinger equation in a nonperturbative fashion with high accuracy. In particular, imaginary-time projection algorithms, also known as diffusion Monte Carlo (DMC) algorithms, have proven to be remarkably successful in the description of nuclei and their global properties, e.g. binding energies, radii, transitions, and reactions, and in the prediction of properties of neutron star matter (for a review of QMC methods see [14]).

In this work, we employ two different DMC techniques, namely the Green’s function Monte Carlo (GFMC) method [72] and the auxiliary field diffusion Monte Carlo (AFDMC) method [73]. Both approaches rely on the application of an imaginary-time propagator to an initial trial wave function in order to project out the true many-body ground state of the system:

[TABLE]

The trial wave function is given in terms of a variational state of the form

[TABLE]

where $F_{C}$ accounts for all of the spin/isospin-independent correlations, and $F_{2}$ and $F_{3}$ are spin/isospin-dependent two- and three-body correlations, respectively. The term $|\Phi\rangle$ is taken to be a shell-model-like state with total angular momentum $J$, parity $\pi$, and total isospin $T$. Its wave function is constructed using single-particle orbitals that depend on the nucleon spatial coordinates, spin, and isospin. An initial optimization procedure is applied to the trial state of Eq. 28 in order to find the optimal parameters providing the best, i.e. lowest, variational energy. The optimized wave function is then repetitively evolved in small imaginary-time steps until the ground state of the system is reached (more details can be found in [14, 17]).

The local chiral interactions considered in this work can be efficiently implemented in both the GFMC and AFDMC methods. The GFMC method, which includes a sum over all possible spin/isospin states at each step in the diffusion, scales exponentially with the number of nucleons $A$. This limits current calculations to around $A=12$. The AFDMC method, on the other hand, samples the sum over all spin/isospin states, and therefore exhibits a much gentler, polynomial scaling with $A$. The two algorithms are thus complementary, and they allow one to vastly extend the region of applicability of QMC calculations. Results employing local chiral forces are now available for several quantities (binding energies, charge radii, charge form factors, single- and two-nucleon radial distributions, and single- and two-nucleon momentum distributions) in light and medium-mass nuclei [69, 16, 74, 17, 75], and for properties of pure neutron systems [76, 77, 78], including pure neutron matter [66, 67, 68, 69].

In QMC methods, the expectation value of an observable $\mathcal{O}$ is calculated as

[TABLE]

where $\{R_{i},S_{i}\}$ are spatial and spin/isospin configurations typically sampled using the Metropolis algorithm [79], and $\mathcal{N}$ is the (large) number of configurations in the simulation. In the AFDMC method, both spatial and spin/isospin degrees of freedom are sampled during the imaginary-time propagation, the latter through the so-called Hubbard-Stratonovich transformation. In the GFMC approach, all possible spin/isospin configuration are included in the trial many-body wave function, and only configurations in coordinate space are sampled. The above expression is valid only for observables that commute with the Hamiltonian. For other observables, such as radii and densities, expectation values are extracted from so-called mixed estimates

[TABLE]

In the above expression, the first term is the mixed estimate (propagated wave function on one side, trial wave function on the other side), and the second term is the variational estimate. This relationship can be derived under the assumption that the variational trial wave function is a good starting point, i.e. that $|\Psi(\tau\to\infty)\rangle=|\Psi_{T}\rangle+|\delta\Psi_{T}\rangle$, with $|\delta\Psi_{T}\rangle$ small. Then, if we calculate the expectation value of an operator between two propagated wave functions and discard terms of $\mathcal{O}(\delta\Psi_{T}^{2})$, we arrive at Eq. 30. Additional details, including the sampling procedure and the calculation of statistical errors, can be found, e.g. in [80].

The SRC scaling factors can be expressed in terms of the central two-nucleon distribution (two-body point-nucleon density) [19]:

[TABLE]

where the central two-nucleon distribution is defined as

[TABLE]

The normalization is such that $\rho_{2,\mathbbm{1}}(A,r)$ integrates to the number of nucleon pairs. Equation 32 involves a mixed estimate and is evaluated according to Eq. 30. In this work, we ensure that the difference between the mixed and variational estimates of the distributions is $\lesssim 10\%$.

In addition to Monte Carlo statistical errors, the use of chiral interactions allows one to estimate the theoretical uncertainties coming from the truncation of the chiral expansion. In this work, we consider results for $\rho_{2,\mathbbm{1}}(A,r)$ at leading-order (LO), next-to-leading-order (NLO), and N2LO, and we estimate the truncation errors on the ratio $X=\frac{2\rho_{2,\mathbbm{1}}(A,r)}{A\rho_{2,\mathbbm{1}}(2,r)}$ entering the definition of the SRC scaling factor of Eqs. 31a and 31b following [81]:

[TABLE]

where we take $Q=m_{\pi}/\Lambda_{b}$ with $m_{\pi}\approx 140\,\rm MeV$ and $\Lambda_{b}=600\,\rm MeV$, as in [74, 17].

Auxiliary field diffusion Monte Carlo calculations for nuclei employing local chiral interactions have been carried out up to $A=16$ [74, 17, 75]. Preliminary results for heavier systems suggest that improved wave functions are necessary to obtain ground-state properties with the same accuracy as for lighter systems. However, such a prescription will increase the computational cost by a factor proportional to $A^{2}$ (see Ref. [17] for details), making calculations for $A\gtrsim 20$ no longer feasible. One way to move beyond oxygen is to use a simplified interaction, capable of capturing most of the ground-state physics of nuclei, for which the employed wave function still gives an accurate description of larger nuclei, thus maintaining the good computational scaling of the current implementation of the AFDMC algorithm. We consider the phenomenological two-body Argonne $v_{4}^{\prime}$ potential ($\text{AV4}^{\prime}$) [82], a simplified version of the more sophisticated Argonne $v_{18}$ (AV18) potential [83], obtained by reprojecting the full potential onto the first four operator channels in order to preserve the phase shifts of lower partial waves and the deuteron binding energy. We note that this potential is very simple and excludes, for example, tensor forces. The Coulomb interaction is, however, still included. Such a potential typically overbinds light nuclei [82]. The inclusion of a repulsive three-body force can be used to compensate for the excessive attraction. As done in other works [84, 85, 86], we consider the central component of the Urbana IX (UIX) interaction [87] as a source of repulsion. In the following, this simplified potential will be referred to as $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$.

3 Results

3.1 Fitting $a_{2}$

In [19], the SRC scaling factors were obtained by taking the limit $r\to 0$ of the ratio of two-body distributions as in Eqs. 31a and 31b. However, this is precisely the region where the Monte Carlo statistical uncertainties become large, see, e.g. Fig. 5. Nuclear potentials are generally repulsive at short distances, and therefore the likelihood of finding two nucleons at small separations is small, giving rise to large statistical uncertainties as $r\to 0$.) In this work, we exploit the fact that, as pointed out in [19], in EFT, “locality” means a shorter distance than the resolution scale. Thus, we expect that we can replace $r\to 0$ in Eqs. 31a and 31b by simply smearing in some region $r<R$, where $R$ is set by the cutoff scale $R_{0}$ (but $R$ is not necessarily equal to $R_{0}$), and still obtain the same $a_{2}$ value. Therefore, we fit a horizontal line to the ratio of two-body distributions $2\rho_{2,\mathbbm{1}}(A,r)/A\rho_{2,\mathbbm{1}}(d,r)$ and $3\rho_{2,\mathbbm{1}}(A,r)/A\rho_{2,\mathbbm{1}}(\isotope[3]{He},r)$ in the region with $0\leq r\leq R$, and we take $R=0.7\,\text{fm}$. This region is chosen as the empirical region where the expected plateau sets in. (Note that for the systems \isotope[3]H, \isotope[3]He, and \isotope[4]He, the results from this linear fit agree with our previous results from [19] using the limit $r\to 0$). We have further checked that varying $R$ from $0.4$ to $1.0\,\text{fm}$ makes a 1%–3% difference in our extracted values of $a_{2}$ for the local chiral interactions. For results with the simplified $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$ and the $\text{AV18}+\text{UIX}$ potentials, varying $R$ in this range makes a $\sim 10\%$ difference (up to 13% for \isotope[40]Ca). For the phenomenological potentials, we use this variation as an estimate for the systematic uncertainty coming from the fit of $a_{2}$.

If each of the discrete values $\{\rho_{2,\mathbbm{1}}(A,r_{i})\}$ obtained from Monte Carlo calculations were equally likely, this procedure would be entirely equivalent to taking the average in the region $0\leq r\leq R$. However, as discussed above, the statistical uncertainties in $\rho_{2,\mathbbm{1}}(A,r)$ grow rapidly as $r\to 0$. In short, our fitting problem is heteroskedastic, and therefore, we use a weighted linear least squares fitting procedure,

[TABLE]

where we take the weight matrix diagonal and equal to the inverse of the Monte Carlo variances for each point:

[TABLE]

In our case $X_{ij}$ reduces to a vector of $1$’s, and the $\{y_{i}\}$ are the set of values $\{a_{2}(r_{i})\}$ from the Monte Carlo simulations. Then our procedure amounts to

[TABLE]

where $\mathcal{M}$ is taken such that $r_{i}\in[0.0,0.7]\,\rm fm$.

3.2 Results for selected nuclei up to \isotope[40]Ca

We first present results for the two-body distributions Eq. 32 for selected nuclei up to $\isotope[40]{Ca}$ in Fig. 4. The left panel shows results for the deuteron, \isotope[6]He, and \isotope[16]O using local chiral interactions at N2LO with the $E\mathbbm{1}$ parameterization of the $3N$ interaction and both cutoffs $R_{0}=1.0,\,1.2\,\rm fm$, whereas the right panel shows results for the deuteron, \isotope[4]He, \isotope[16]O, and \isotope[40]Ca for the simplified nuclear potential $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$. The figure shows the definite scheme and scale dependence of these distributions. This is especially clear in the left panel where the distributions are calculated in an EFT framework at two cutoff scales. The softer cutoff $R_{0}=1.2\,\rm fm$ resembles more a mean-field calculation at short distances where what are typically referred to as SRCs are reduced, leading to a higher probability to find two nucleons separated by very short distances $r\lesssim 1.0\,\rm fm$. The right panel, which utilizes phenomenological potentials where the effective cutoff is much harder (though a particular value is not identified) shows a significantly lower probability to find a pair of nucleons separated by $r\lesssim 1.0\,\rm fm$.

In contrast to the two-body distributions shown in Fig. 4, their ratios to the deuteron and \isotope[3]He two-body distributions, i.e. Eqs. 31a and 31b are largely scheme and scale independent. In Fig. 5, we show the ratio $a_{2}(A/d)$ at short internucleon distances $0\leq r\leq 1.0$ fm for \isotope[6]He and \isotope[12]C using chiral EFT interactions at N2LO with the $E\tau$ parameterization of the $3N$ interaction (left and middle panels) and for \isotope[40]Ca using the simplified $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$ potential (right panel). For \isotope[6]He and \isotope[12]C, we show the Monte Carlo results with statistical uncertainties using $R_{0}=1.0\,\rm fm$ as the blue squares with error bars, while we use green squares for \isotope[40]Ca using the simplified $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$ potentials. The blue (red) band represents the combined statistical and systematic uncertainties coming from the truncation of the chiral expansion for the $R_{0}=1.0\,\rm fm$ ($1.2\,\rm fm$) cutoff. The light blue band in the right panel shows the uncertainty in the fit for \isotope[40]Ca. The light blue horizontal lines indicate the weighted linear fits to the Monte Carlo results as described above. Of the three cases shown here, there is currently only an experimental result with which to compare for \isotope[12]C. This is shown in the middle panel as the black dashed line with the gray band representing the experimental uncertainty.

While Fig. 5 illustrates the method by which we extract the SRC scaling factors, in Fig. 6 we show our main predictions. The left panel shows results for $a_{2}(A/d)$ for selected nuclei from \isotope[3]H up to \isotope[40]Ca. The blue squares (red circles) show the results for chiral interactions up to N2LO with the $E\tau$ parameterization of the $3N$ interaction and the cutoff $R_{0}=1.0\,\rm fm$ ($R_{0}=1.2\,\rm fm$). The green upward-pointing (downward-pointing) triangles show the results using the $\text{AV18}+\text{UIX}$ ($\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$) potentials. The black stars show the experimental results from [32], where available. The light colored bands show the systematic uncertainties stemming from the truncation of the chiral expansion at N2LO for the $R_{0}=1.0\,\rm fm$ ($R_{0}=1.2\,\rm fm$) cutoff (where available), and coming from the fit of $a_{2}$ in the case of the phenomenological potentials. The gray region appearing at large $A$ represents an expected saturation region taken as the difference, including uncertainties, between $a_{2}(\isotope[197]{Au}/d)=5.16(22)$ and $a_{2}(\isotope[63]{Cu}/d)=5.21(20)$ [32], i.e. we estimate $\lim_{A\to\infty}a_{2}(A/d)\sim 4.94\text{--}5.41$. Similarly, the right panel of Fig. 6 shows results for $a_{2}(A/\isotope[3]{He})$ for selected nuclei from \isotope[4]He up to \isotope[40]Ca using the same color and symbol scheme as in the left panel. Note that the gray saturation region is provided by a single experimental value at large $A$: namely, $a_{2}(\isotope[56]{Fe}/\isotope[3]{He})=2.83(3)(18)$ [88]. (Table 1 collect these and more results using both parametrizations $E\tau$ and $E\mathbbm{1}$ of the $3N$ interaction and both cutoffs $R_{0}=1.0$, 1.2 fm, as well as results for phenomenological potentials and experimental results).

The results in Fig. 6 compare very well with experimental values, where available. In particular, we find 0.0%, 4.4%, and 1.7% relative agreement between our results for $a_{2}(A/d)$ using chiral interactions at N2LO with the cutoff $R_{0}=1.0\,\rm fm$ and experiment for \isotope[3]He, \isotope[4]He, and \isotope[12]C, respectively. Results using the softer cutoff $R_{0}=1.2\,\rm fm$ are typically higher than for the lower cutoff by $\sim 5\text{--}10\%$ (an exception occurs for \isotope[16]O, where the softer interaction with the $E\tau$ parametrization has already been found to exhibit significant overbinding [17]), but are always within the estimated systematic uncertainties. It is also interesting to note that the predicted values for $a_{2}$ for the $A=6$ systems fall below the values for \isotope[4]He, placing them between \isotope[3]He and \isotope[4]He along the fitted line in Fig. 3. As suggested in [28], what appears to dictate the strength of the EMC effect (and therefore the height of the SRC scaling plateaus through the EMC-SRC linear relation) is the local nuclear density. Given that \isotope[4]He is such a compact nucleus, and that both \isotope[6]Li and \isotope[6]He can be thought of as $\alpha$ particles with additional nucleons “orbiting,” one might expect that the strong attraction of the $\alpha$ core to the orbiting nucleons would tend to lower the local central two-nucleon density. These predictions for \isotope[6]Li could be tested already using existing experimental setups for $(e,e^{\prime})$ inclusive scattering in QE kinematics at Jefferson Lab. For \isotope[6]He, these predictions could be tested at future rare isotope facilities such as the Facility for Antiproton and Ion Research with experiments in inverse kinematics using a \isotope[6]He beam on a proton target inducing $(p,2p)$ reactions.

[FIGURE:]

We also make predictions for \isotope[16]O and \isotope[40]Ca in Figs. 6 and 1. While the latter is only calculated using the simplified phenomenological potential $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$, our expectation based on calculations for light systems with $3\leq A\leq 16$ is that this rather simplified Hamiltonian is capturing most of the important SRC physics: this can be seen by comparing the results using the realistic chiral EFT interactions at N2LO with $R_{0}=1.0\,\rm fm$ (blue squares in Fig. 6) with the results using $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$ (green downward-pointing triangles). The relative agreement between the results is 1.9%, 1.9%, and 7.4% for \isotope[3]He, \isotope[4]He, and \isotope[16]O, respectively. We also refer the reader to Table 2: Both the binding energies and radii for \isotope[4]He, \isotope[16]O, and \isotope[40]Ca are reasonably well reproduced using $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$. Nevertheless, given the relative agreement between our chiral interactions at N2LO with the cutoff $R_{0}=1.0\,\rm fm$ and the simplified potential $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$, and the slight systematic underbinding of the latter, we assign a conservative uncertainty to our $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$ calculations, e.g. $a_{2}(\isotope[40]{Ca}/d)=5.15(67)$ and $a_{2}(\isotope[40]{Ca}/\isotope[3]{He})=2.55(33)$. This 13% can be justified from our study of the sensitivity of the extracted $a_{2}$ to the chosen region $0\leq r\leq R$.

4 Summary

In this work, we have used DMC algorithms, namely the GFMC and AFDMC methods, to calculate the SRC scaling factors $a_{2}(A/d)$ and $a_{2}(A/\isotope[3]{He})$ for nuclei from $A=3$ to $A=40$. We have reviewed in detail the derivation of $a_{2}$ from EFT, arguing that isovector corrections are very small. We have then shown that fitting a constant to the ratio of two-body central densities in some empirical region $0\leq r\leq R$ reproduces the values from our previous work [19] and provides a reliable method to extract SRC scaling factors. Where experimental values exist, our calculations agree very well using both chiral EFT interactions at N2LO and phenomenological potentials, including the simplified $\text{AV4}^{\prime}+\text{UIX}_{\text{c}}$ potential, providing further evidence of the value of the novel framework first proposed in [19]. We also show the first ab initio predictions for SRC scaling factors for \isotope[6]Li, \isotope[6]He, \isotope[16]O, and \isotope[40]Ca. These predictions could be tested in future experiments, offering intriguing insights into the evolution of SRC scaling factors with the nuclear mass $A$. Our framework may also shed light on the proposed, but so far elusive, $3N$ SRC scaling. This topic is currently being investigated and we leave it for future work.

We thank N Fomin for providing us with the data from [27] and T Aumann, H-W Hammer, K Hebeler, O Hen, and I Tews for valuable discussions. The work of JL and AS was supported by the ERC Grant No. 307986 STRONGINT and the BMBF under Contract No. 05P18RDFN1. The work of DL was supported by the US Department of Energy, Office of Science, Office of Nuclear Physics, under Contract No. DE-SC0013617, and by the NUCLEI SciDAC program. The work of JC and SG was supported by the NUCLEI SciDAC program, by the US Department of Energy, Office of Science, Office of Nuclear Physics, under contract No. DE-AC52-06NA25396, and by the LDRD program at LANL. SG was also supported by the DOE Early Career research Program. J-W C is partly supported by the Ministry of Science and Technology, Taiwan, under Grant No. 105-2112-M-002-017-MY3 and the Kenda Foundation. WD is supported by the US Department of Energy under grant DE-SC0011090. WDUS is also supported within the framework of the TMD Topical Collaboration of the US Department of Energy, Office of Science, Office of Nuclear Physics, and by the SciDAC4 award DE-SC0018121. Computational resources have been provided by the Lichtenberg high performance computer of the TU Darmstadt, by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the US Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001 and by the National Energy Research Scientific Computing Center (NERSC), which is supported by the US Department of Energy, Office of Science, under contract No. DE-AC02-05CH11231.