Variational calculation of the ground state of closed-shell nuclei up to $A=40$

TL;DR

This paper uses variational methods with realistic nuclear potentials to compute ground-state properties of light nuclei up to calcium-40, revealing insights into three-nucleon interactions and high-momentum nuclear behavior.

Contribution

It introduces a variational approach with cluster expansions and nonlinear optimization to accurately calculate nuclear properties, including effects of three-nucleon forces.

Findings

Three-nucleon interaction becomes repulsive for A≥16.

Including three-nucleon forces improves descriptions of 16O properties.

High-momentum distributions show universal behavior across nuclei.

Abstract

Variational calculations of ground-state properties of $^4$He, $^{16}$O, and $^{40}$Ca are carried out employing realistic phenomenological two- and three-nucleon potentials. The trial wave function includes two- and three-body correlations acting on a product of single-particle determinants. Expectation values are evaluated with a cluster expansion for the spin-isospin dependent correlations considering up to five-body cluster terms. The optimal wave function is obtained by minimizing the energy expectation value over a set of up to 20 parameters by means of a nonlinear optimization library. We present results for the binding energy, charge radius, one- and two-body densities, single-nucleon momentum distribution, charge form factor, and Coulomb sum rule. We find that the employed three-nucleon interaction becomes repulsive for $A\geq16$. In $^{16}$O the inclusion of such a force provides a better description of the properties of the nucleus. In $^{40}$Ca instead, the repulsive behavior of the three-body interaction fails to reproduce experimental data for the charge radius and the charge form factor. We find that the high-momentum region of the momentum distributions, determined by the short-range terms of nuclear correlations, exhibit a universal behavior independent of the particular nucleus. The comparison of the Coulomb sum rules for $^4$He, $^{16}$O, and $^{40}$Ca reported in this work will help elucidate in-medium modifications of the nucleon form factors.

Tables

Table 1. Table 1: Cluster contributions to the energy per nucleon and point radius in 4 He when using the AV18 potential. Energies are in MeV/ A 𝐴 A and radii in fm 2 . Averages are calculated using 𝒩 c = 2 × 10 6 subscript 𝒩 𝑐 2 superscript 10 6 \mathcal{N}_{c}=2\times 10^{6} configurations. In this and the following tables, the Monte Carlo statistical errors of the last digits are given in parentheses. A (0) indicates an error of less than 5 in the following digit.

observable	1b	2b	3b	4b	sum
$T$	$12.22\left(2\right)$	$12.51\left(1\right)$	$-1.12\left(1\right)$	$0.23\left(1\right)$	$23.84\left(3\right)$
$v_{ij}^{1-6}$		$-30.07\left(3\right)$	$0.52\left(1\right)$	$-0.22\left(1\right)$	$-29.77\left(3\right)$
$v_{ij}^{7-14}$		$-1.00\left(0\right)$	$1.07\left(0\right)$	$-0.08\left(0\right)$	$-0.01\left(0\right)$
$v_{ij}^{\gamma }$		$0.20\left(0\right)$	$0.01\left(0\right)$	$0.00\left(0\right)$	$0.21\left(0\right)$
$v_{ij}$		$-30.86\left(3\right)$	$1.58\left(1\right)$	$-0.29\left(1\right)$	$-29.57\left(3\right)$
$T+v_{ij}$	$12.22\left(2\right)$	$-18.35\left(2\right)$	$0.47\left(1\right)$	$-0.06\left(0\right)$	$-5.73\left(1\right)$
$r_{\mathrm{pt}}^2$	$2.813\left(2\right)$	$-0.568\left(1\right)$	$0.046\left(0\right)$	$-0.003\left(0\right)$	$2.289\left(1\right)$

Table 2. Table 2: Cluster contributions to the energy per nucleon and point radius in 4 He when using the AV18+UIX potential. Energies are in MeV/ A 𝐴 A and radii in fm 2 . Averages are calculated using 𝒩 c = 2 × 10 6 subscript 𝒩 𝑐 2 superscript 10 6 \mathcal{N}_{c}=2\times 10^{6} configurations.

observable	1b	2b	3b	4b	sum
$T$	$13.12\left(2\right)$	$15.55\left(2\right)$	$-2.30\left(1\right)$	$0.44\left(1\right)$	$26.80\left(3\right)$
$v_{ij}^{1-6}$		$-33.74\left(3\right)$	$1.45\left(1\right)$	$-0.30\left(1\right)$	$-32.60\left(3\right)$
$v_{ij}^{7-14}$		$-1.31\left(0\right)$	$1.66\left(0\right)$	$-0.19\left(0\right)$	$0.17\left(3\right)$
$v_{ij}^{\gamma }$		$0.21\left(0\right)$	$0.01\left(0\right)$	$0.00\left(0\right)$	$0.22\left(0\right)$
$v_{ij}$		$-34.83\left(3\right)$	$3.10\left(2\right)$	$-0.48\left(1\right)$	$-32.21\left(3\right)$
$V_{ijk}^{2\pi ,A}$			$-1.67\left(0\right)$	$0.04\left(0\right)$	$-1.64\left(0\right)$
$V_{ijk}^{2\pi ,C}$			$-1.02\left(0\right)$	$0.04\left(0\right)$	$-0.98\left(0\right)$
$V_{ijk}^R$			$1.36\left(0\right)$	$-0.03\left(0\right)$	$1.33\left(0\right)$
$V_{ijk}$			$-1.34\left(0\right)$	$0.05\left(0\right)$	$-1.29\left(0\right)$
$T+v_{ij}$	$13.12\left(2\right)$	$-19.29\left(2\right)$	$0.80\left(0\right)$	$-0.04\left(1\right)$	$-5.41\left(1\right)$
$T+v_{ij}+V_{ijk}$	$13.12\left(2\right)$	$-19.29\left(2\right)$	$-0.54\left(1\right)$	$0.01\left(1\right)$	$-6.70\left(1\right)$
$r_{\mathrm{pt}}^2$	$2.646\left(2\right)$	$-0.600\left(1\right)$	$0.070\left(0\right)$	$-0.005\left(0\right)$	$2.111\left(1\right)$

Table 3. Table 3: Cluster contributions to the energy per nucleon and point radius in 16 O when using the AV18 potential. Energies are in MeV/ A 𝐴 A and radii in fm 2 . Averages are calculated using 𝒩 c = 6 × 10 5 subscript 𝒩 𝑐 6 superscript 10 5 \mathcal{N}_{c}=6\times 10^{5} configurations.

observable	1b	2b	3b	4b	5b	6-16b	sum
$T$	$19.78\left(2\right)$	$16.67\left(2\right)$	$-6.13\left(3\right)$	$2.24\left(3\right)$	$-0.44\left(3\right)$	$0.07\left(1\right)$	$32.19\left(5\right)$
$v_{ij}^{1-6}$		$-48.30\left(4\right)$	$10.69\left(3\right)$	$-2.42\left(4\right)$	$-0.59\left(5\right)$	$-0.19\left(2\right)$	$-40.81\left(6\right)$
$v_{ij}^{7-14}$		$-0.46\left(0\right)$	$3.48\left(1\right)$	$-1.70\left(1\right)$	$0.69\left(1\right)$	$-0.20\left(2\right)$	$1.82\left(3\right)$
$v_{ij}^{\gamma }$		$0.91\left(0\right)$	$0.03\left(0\right)$	$-0.02\left(0\right)$	$0.00\left(0\right)$	$0.00\left(0\right)$	$0.91\left(0\right)$
$v_{ij}$		$-47.85\left(4\right)$	$14.20\left(3\right)$	$-4.15\left(4\right)$	$0.10\left(5\right)$	$0.00\left(0\right)$	$-37.70\left(3\right)$
$T+v_{ij}$	$19.78\left(2\right)$	$-31.18\left(2\right)$	$8.07\left(2\right)$	$-1.91\left(3\right)$	$-0.34\left(3\right)$	—	$-5.51\left(2\right)$
$r_{\mathrm{pt}}^2$	$6.235\left(2\right)$	$-0.672\left(2\right)$	$0.277\left(1\right)$	$-0.074\left(1\right)$	$-0.010\left(1\right)$	$-0.002\left(0\right)$	$5.754\left(3\right)$

Table 4. Table 4: Cluster contributions to the energy per nucleon and point radius in 16 O when using the AV18+UIX potential. Energies are in MeV/ A 𝐴 A and radii in fm 2 . Averages are calculated using 𝒩 c = 6 × 10 5 subscript 𝒩 𝑐 6 superscript 10 5 \mathcal{N}_{c}=6\times 10^{5} configurations.

observable	1b	2b	3b	4b	5b	6-16b	sum
$T$	$16.29\left(2\right)$	$16.34\left(2\right)$	$-3.79\left(2\right)$	$0.90\left(2\right)$	$-0.05\left(3\right)$	$0.00\left(0\right)$	$29.70\left(4\right)$
$v_{ij}^{1-6}$		$-41.90\left(3\right)$	$6.61\left(2\right)$	$-0.95\left(4\right)$	$-0.36\left(3\right)$	$-0.22\left(4\right)$	$-36.83\left(8\right)$
$v_{ij}^{7-14}$		$-0.82\left(0\right)$	$2.93\left(1\right)$	$-1.01\left(1\right)$	$0.27\left(1\right)$	$-0.06\left(1\right)$	$1.31\left(2\right)$
$v_{ij}^{\gamma }$		$0.83\left(0\right)$	$0.03\left(0\right)$	$-0.02\left(0\right)$	$0.00\left(0\right)$	$0.00\left(0\right)$	$0.84\left(0\right)$
$v_{ij}$		$-41.89\left(3\right)$	$9.57\left(2\right)$	$-1.98\left(4\right)$	$-0.10\left(4\right)$	$-0.01\left(0\right)$	$-34.41\left(2\right)$
$V_{ijk}^{2\pi ,A}$			$-2.49\left(1\right)$	$1.12\left(0\right)$	$-0.29\left(1\right)$	$0.06\left(1\right)$	$-1.59\left(2\right)$
$V_{ijk}^{2\pi ,C}$			$-1.69\left(0\right)$	$0.79\left(0\right)$	$-0.26\left(0\right)$	$0.06\left(0\right)$	$-1.10\left(0\right)$
$V_{ijk}^R$			$2.83\left(1\right)$	$-0.70\left(0\right)$	$0.12\left(1\right)$	$-0.02\left(0\right)$	$2.23\left(1\right)$
$V_{ijk}$			$-1.35\left(1\right)$	$1.21\left(0\right)$	$-0.43\left(1\right)$	$0.11\left(1\right)$	$-0.45\left(2\right)$
$T+v_{ij}$	$16.29\left(2\right)$	$-25.54\left(2\right)$	$5.78\left(2\right)$	$-1.08\left(2\right)$	$-0.14\left(2\right)$	—	$-4.70\left(2\right)$
$T+v_{ij}+V_{ijk}$	$16.29\left(2\right)$	$-25.54\left(2\right)$	$4.43\left(2\right)$	$0.13\left(2\right)$	$-0.57\left(2\right)$	—	$-5.15\left(2\right)$
$r_{\mathrm{pt}}^2$	$7.353\left(3\right)$	$-0.680\left(2\right)$	$0.217\left(1\right)$	$-0.035\left(2\right)$	$-0.007\left(1\right)$	$-0.002\left(0\right)$	$6.846\left(3\right)$

Table 5. Table 5: Cluster contributions to the energy per nucleon and point radius in 40 Ca when using the AV18 potential. Energies are in MeV/ A 𝐴 A and radii in fm 2 . Averages are calculated using 𝒩 c = 5 × 10 5 subscript 𝒩 𝑐 5 superscript 10 5 \mathcal{N}_{c}=5\times 10^{5} configurations.

observable	1b	2b	3b	4b	5b	6-40b	sum
$T$	$20.80\left(1\right)$	$18.00\left(1\right)$	$-8.60\left(3\right)$	$4.04\left(3\right)$	$-1.22\left(10\right)$	$0.28\left(6\right)$	$32.29\left(15\right)$
$v_{ij}^{1-6}$		$-54.52\left(3\right)$	$15.80\left(2\right)$	$-5.09\left(6\right)$	$0.05\left(9\right)$	$0.00\left(0\right)$	$-43.76\left(8\right)$
$v_{ij}^{7-14}$		$-0.16\left(0\right)$	$4.45\left(1\right)$	$-2.78\left(1\right)$	$1.47\left(2\right)$	$-0.51\left(5\right)$	$2.47\left(7\right)$
$v_{ij}^{\gamma }$		$1.86\left(0\right)$	$0.07\left(0\right)$	$-0.08\left(0\right)$	$0.00\left(0\right)$	$0.00\left(0\right)$	$1.85\left(0\right)$
$v_{ij}$		$-52.81\left(3\right)$	$20.31\left(2\right)$	$-7.95\left(6\right)$	$1.52\left(10\right)$	$-0.24\left(3\right)$	$-39.17\left(11\right)$
$T+v_{ij}$	$20.80\left(2\right)$	$-34.81\left(2\right)$	$11.71\left(3\right)$	$-3.91\left(5\right)$	$0.30\left(9\right)$	—	$-5.88\left(10\right)$
$r_{\mathrm{pt}}^2$	$11.204\left(3\right)$	$-0.920\left(4\right)$	$0.506\left(3\right)$	$-0.188\left(5\right)$	$0.006\left(7\right)$	$0.000\left(0\right)$	$10.609\left(7\right)$

Table 6. Table 6: Cluster contributions to the energy per nucleon and point radius in 40 Ca when using the AV18+UIX potential. Energies are in MeV/ A 𝐴 A and radii in fm 2 . Averages are calculated using 𝒩 c = 5 × 10 5 subscript 𝒩 𝑐 5 superscript 10 5 \mathcal{N}_{c}=5\times 10^{5} configurations.

observable	1b	2b	3b	4b	5b	6-40b	sum
$T$	$17.35\left(1\right)$	$17.86\left(1\right)$	$-5.65\left(2\right)$	$1.83\left(3\right)$	$-0.54\left(9\right)$	$0.12\left(5\right)$	$30.97\left(14\right)$
$v_{ij}^{1-6}$		$-47.20\left(3\right)$	$10.20\left(2\right)$	$-2.33\left(5\right)$	$0.09\left(8\right)$	$0.00\left(0\right)$	$-39.24\left(7\right)$
$v_{ij}^{7-14}$		$-0.67\left(0\right)$	$3.68\left(1\right)$	$-1.60\left(1\right)$	$0.52\left(2\right)$	$-0.13\left(1\right)$	$1.80\left(3\right)$
$v_{ij}^{\gamma }$		$1.72\left(0\right)$	$0.06\left(0\right)$	$-0.05\left(0\right)$	$0.00\left(0\right)$	$0.00\left(0\right)$	$1.73\left(0\right)$
$v_{ij}$		$-46.15\left(2\right)$	$13.94\left(2\right)$	$-3.99\left(5\right)$	$0.62\left(8\right)$	$-0.08\left(2\right)$	$-35.66\left(9\right)$
$V_{ijk}^{2\pi ,A}$			$-2.98\left(0\right)$	$1.86\left(0\right)$	$-0.71\left(1\right)$	$0.20\left(1\right)$	$-1.63\left(2\right)$
$V_{ijk}^{2\pi ,C}$			$-2.10\left(0\right)$	$1.33\left(0\right)$	$-0.61\left(1\right)$	$0.19\left(1\right)$	$-1.19\left(2\right)$
$V_{ijk}^R$			$3.59\left(1\right)$	$-1.25\left(0\right)$	$0.28\left(1\right)$	$-0.05\left(0\right)$	$2.57\left(1\right)$
$V_{ijk}$			$-1.48\left(0\right)$	$1.93\left(1\right)$	$-1.05\left(1\right)$	$0.37\left(3\right)$	$-0.23\left(4\right)$
$T+v_{ij}$	$17.35\left(1\right)$	$-28.29\left(2\right)$	$8.30\left(2\right)$	$-2.16\left(3\right)$	$0.08\left(8\right)$	—	$-4.69\left(9\right)$
$T+v_{ij}+V_{ijk}$	$17.35\left(1\right)$	$-28.29\left(2\right)$	$6.81\left(2\right)$	$-0.23\left(3\right)$	$-0.97\left(8\right)$	—	$-4.92\left(10\right)$
$r_{\mathrm{pt}}^2$	$13.025\left(3\right)$	$-0.903\left(4\right)$	$0.393\left(2\right)$	$-0.091\left(5\right)$	$-0.022\left(7\right)$	$-0.007\left(4\right)$	$12.394\left(11\right)$

Table 7. Table 7: Total energies (in MeV/ A 𝐴 A ) and charge radii (in fm) in 4 He for different potentials. VMC and GFMC results are taken from Refs.  Pieper  et al. ( 2001 ); Pieper and Wiringa ( 2001 ) .

obs	potential	CVMC	VMC	GFMC	exp
$⟨E⟩$	AV18	$-5.73\left(1\right)$	$-5.93\left(1\right)$	$-6.02\left(1\right)$	$-7.07$
	AV18+UIX	$-6.70\left(1\right)$	$-6.95\left(1\right)$	$-7.08\left(1\right)$
$\sqrt{⟨r_{\mathrm{ch}}^2⟩}$	AV18	$1.725\left(3\right)$	$1.734\left(3\right)$	—	$1.676\left(3\right)$
	AV18+UIX	$1.673\left(3\right)$	$1.665\left(3\right)$	$1.661\left(3\right)$

Table 8. Table 8: Total energies (in MeV/ A 𝐴 A ) and charge radii (in fm) in 16 O for different potentials.

obs	potential	CVMC	exp
$⟨E⟩$	AV18	$-5.51\left(2\right)$	$-7.98$
	AV18+UIX	$-5.15\left(2\right)$
$\sqrt{⟨r_{\mathrm{ch}}^2⟩}$	AV18	$2.538\left(2\right)$	$2.699\left(5\right)$
	AV18+UIX	$2.745\left(2\right)$

Table 9. Table 10: Total energies (in MeV/ A 𝐴 A ) and charge radii (in fm) in 4 He and 16 O for the AV6 ′ potential. The electromagnetic term v i ​ j γ superscript subscript 𝑣 𝑖 𝑗 𝛾 v_{ij}^{\gamma} is not included. AFDMC energies are taken from Ref.  Gandolfi  et al. ( 2014a ) .

obs	nucleus	CVMC	AFDMC
$⟨E⟩$	4He	$-6.53\left(1\right)$	$-6.77\left(1\right)$
	16O	$-6.79\left(3\right)$	$-7.23\left(2\right)$
$\sqrt{⟨r_{\mathrm{ch}}^2⟩}$	4He	$1.678\left(3\right)$	$1.674\left(9\right)$
	16O	$2.580\left(2\right)$	$2.52\left(3\right)$

Table 10. Table 11: Variational parameters for 4 He, 16 O, and 40 Ca.

param	4He		16O		40Ca
	AV18	AV18+UIX	AV18	AV18+UIX	AV18	AV18+UIX
$V_S$	$-47.824\mathrm{MeV}$	$-51.459\mathrm{MeV}$	$-44.860\mathrm{MeV}$	$-43.032\mathrm{MeV}$	$-47.261\mathrm{MeV}$	$-45.606\mathrm{MeV}$
$R_s$	$2.174\mathrm{fm}$	$2.039\mathrm{fm}$	$3.325\mathrm{fm}$	$3.542\mathrm{fm}$	$4.592\mathrm{fm}$	$4.872\mathrm{fm}$
$a_s$	$0.371\mathrm{fm}$	$0.340\mathrm{fm}$	$0.439\mathrm{fm}$	$0.641\mathrm{fm}$	$0.690\mathrm{fm}$	$0.930\mathrm{fm}$
${\alpha }_s$	$0.126$	$0.285$	$-0.056$	$-0.091$	$-0.056$	$-0.091$
${\rho }_s$	$1.643\mathrm{fm}$	$1.131\mathrm{fm}$	$1.847\mathrm{fm}$	$1.009\mathrm{fm}$	$1.847\mathrm{fm}$	$1.009\mathrm{fm}$
$k_F$	$1.864{\mathrm{fm}}^{-1}$	$1.7036{\mathrm{fm}}^{-1}$	$1.604{\mathrm{fm}}^{-1}$	$1.309{\mathrm{fm}}^{-1}$	$1.604{\mathrm{fm}}^{-1}$	$1.309{\mathrm{fm}}^{-1}$
$\alpha$	$0.759$	$0.745$	$0.787$	$0.893$	$0.787$	$0.893$
${\beta }_c$	$0.976$	$0.994$	$1.078$	$1.172$	$1.078$	$1.172$
${\beta }_t$	$1.206$	$1.341$	$1.130$	$1.194$	$1.130$	$1.194$
$d_S$	$3.361\mathrm{fm}$	$3.539\mathrm{fm}$	$2.787\mathrm{fm}$	$2.502\mathrm{fm}$	$2.787\mathrm{fm}$	$2.502\mathrm{fm}$
$d_P$	$4.711\mathrm{fm}$	$4.039\mathrm{fm}$	$2.867\mathrm{fm}$	$3.212\mathrm{fm}$	$2.867\mathrm{fm}$	$3.212\mathrm{fm}$
$d_t$	$6.449\mathrm{fm}$	$6.716\mathrm{fm}$	$4.655\mathrm{fm}$	$4.312\mathrm{fm}$	$4.655\mathrm{fm}$	$4.312\mathrm{fm}$
$t_1$	$5.792$	$5.769$	$5.165$	$4.097$	$5.165$	$4.097$
$t_2$	$4$	$4$	$4$	$4$	$4$	$4$
$t_3$	$0.127{\mathrm{fm}}^{-1}$	$0.117{\mathrm{fm}}^{-1}$	$0.252{\mathrm{fm}}^{-1}$	$0.202{\mathrm{fm}}^{-1}$	$0.252{\mathrm{fm}}^{-1}$	$0.202{\mathrm{fm}}^{-1}$
${\epsilon }_{2\pi ,A}$		$-9.60\cdot {10}^{-4}$		$-8.77\cdot {10}^{-4}$		$-8.77\cdot {10}^{-4}$
${\epsilon }_R$		$-8.22\cdot {10}^{-4}$		$-7.87\cdot {10}^{-4}$		$-7.87\cdot {10}^{-4}$
$\eta$		$0.693$		$1.005$		$1.005$
$c_y$		$1.337{\mathrm{fm}}^{-2}$		$1.601{\mathrm{fm}}^{-2}$		$1.619{\mathrm{fm}}^{-2}$
$c_t$		$1.811{\mathrm{fm}}^{-2}$		$1.616{\mathrm{fm}}^{-2}$		$1.734{\mathrm{fm}}^{-2}$

Equations

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

$$v_{ij}=v_{ij}^{\gamma }+v_{ij}^{\pi }+v_{ij}^R.$$

$$v_{ij}=\sum _{p=1}^{18}{v}^p(r_{ij}){\mathcal{O}}_{ij}^p.$$

$$\displaystyle\mathcal{O}^{p=1,6}_{ij}=\big{[}\mathbbm{1},\bm{\sigma}_{i}\cdot\bm{\sigma}_{j},S_{ij}\big{]}\otimes\big{[}\mathbbm{1},\bm{\tau}_{i}\cdot\bm{\tau}_{j}\big{]},$$

$$S_{ij}=3({\mathbf{\sigma }}_i\cdot {\hat{\mathbf{r}}}_{ij})({\mathbf{\sigma }}_j\cdot {\hat{\mathbf{r}}}_{ij})-({\mathbf{\sigma }}_i\cdot {\mathbf{\sigma }}_j),$$

$${\mathcal{O}}_{ij}^{p=7,14}$$

$$\displaystyle\,\otimes\big{[}\mathbbm{1},\bm{\tau}_{i}\cdot\bm{\tau}_{j}\big{]},$$

$$\begin{array}{rl}\mathbf{L}& =\frac{1}{2i}\left({\mathbf{r}}_i-{\mathbf{r}}_j\right)\times \left({\nabla }_i-{\nabla }_j\right),\\ \mathbf{S}& =\frac{1}{2}\left({\mathbf{\sigma }}_i+{\mathbf{\sigma }}_j\right),\end{array}$$

$$\displaystyle\begin{aligned} \mathcal{O}^{p=15,17}_{ij}&=\big{[}\mathbbm{1},\bm{\sigma}_{i}\cdot\bm{\sigma}_{j},S_{ij}\big{]}\otimes T_{ij},\\ \mathcal{O}^{p=18}_{ij}&=\tau_{z_{i}}+\tau_{z_{j}},\end{aligned}$$

$$V_{ijk}=V_{ijk}^{2\pi ,A}+V_{ijk}^{2\pi ,C}+V_{ijk}^R.$$

$$V_{ijk}^{2\pi ,A}$$

$$V_{ijk}^{2\pi ,C}$$

$$X_{ij}$$

$$Y_{\pi }(x)$$

$$T_{\pi }(x)$$

$$\xi (r)=1-e^{-c{r}^2}.$$

$$V_{ijk}^R=A_R\sum _{cyc}{T}_{\pi }^2({\mu }_{\pi }{r}_{ij})T_{\pi }^2({\mu }_{\pi }{r}_{ik}).$$

$$|{\Psi }_V⟩$$

$$\displaystyle\quad\;\times\Bigg{[}1+\sum_{i<j}U^{7-8}_{ij}\Bigg{]}|\Psi_{J}\rangle,$$

$$|{\Psi }_J⟩$$

$$|\Phi ⟩=\left\{D_{p\uparrow }{D}_{p\downarrow }{D}_{n\uparrow }{D}_{n\downarrow }\right\},$$

$${\varphi }_{nlm}(\mathbf{r})=R_{nl}(r)Y_{lm}(\theta ,\phi ),\mathbf{r}=(r,\theta ,\phi ),$$

$${\tilde{\mathbf{r}}}_i={\mathbf{r}}_i-{\mathbf{R}}_{c.m.},{\mathbf{R}}_{c.m.}=\frac{1}{A}\sum _i{\mathbf{r}}_i,$$

$$\displaystyle V(r)=V_{s}\Bigg{[}\frac{1}{1+e^{(r-R_{s})/a_{s}}}-\alpha_{s}e^{-(r/\rho_{s})^{2}}\Bigg{]},$$

$$U_{ij}=\sum _{p=2}^8{\beta }_p{u}_p(r_{ij}){\mathcal{O}}_{ij}^p,$$

$${\bar{v}}_{ij}=\sum _{p=1}^{14}{\alpha }_p{v}^p(r_{ij}){\mathcal{O}}_{ij}^p,$$

$$\begin{array}{rl}{f}_c(r\ge d_1)& =1,\\ u_p(r\ge d_p)& =0.\end{array}$$

$$\begin{array}{rlrl}{\beta }_{p=2-4,7-8}& ={\beta }_c,& {\alpha }_{p=1,9,13-14}=& 1,\\ {\beta }_{p=5-6}& ={\beta }_t,& {\alpha }_{p=2-8,10-12}=& \alpha ,\end{array}$$

$$\begin{array}{rl}{d}_{p=1-4,7-8}& =d_{\mathrm{S}},d_{\mathrm{P}},\\ d_{p=5-6}& =d_t,\end{array}$$

$$U_{ij}\to \prod _{k\ne i,j}{f}_3(r_{ij};r_{ik},r_{jk})U_{ij},$$

$$\begin{array}{rl}& f_3(r_{ij};r_{ik},r_{jk})=1-t_1{\left(\frac{r_{ij}}{R_{ijk}}\right)}^{t_2}{e}^{-t_3{R}_{ijk}},\\ & R_{ijk}=r_{ij}+r_{ik}+r_{jk}.\end{array}$$

$$U_{ijk}={\epsilon }_{2\pi ,A}{\tilde{V}}_{ijk}^{2\pi ,A}+{\epsilon }_R{\tilde{V}}_{ijk}^R,$$

$$E_V=\frac{⟨{\Psi }_V|H|{\Psi }_V⟩}{⟨{\Psi }_V|{\Psi }_V⟩}\ge E_0,$$

Full text

I Introduction

Atomic nuclei are self-bound systems of strongly interacting fermions. Understanding their structure, reactions, and electroweak properties in terms of the individual interactions among their constituents, protons and neutrons, has been a long-standing goal of theoretical nuclear physics. Ab initio approaches are aimed at solving the many-body Schrödinger equation associated with the nuclear Hamiltonian. This is made particularly difficult by the strong coupling of spin and spatial degrees of freedom which characterize nuclear forces. In addition, the nuclear many-body solution has to feature a self-emerging shell structure and should be able to encompass clusters of highly correlated nucleons.

One of the key advantages of ab initio approaches is that they allow the disentanglement of the theoretical uncertainty coming from modeling the nuclear potential and currents from that due to the approximations inherent in other many-body techniques. This is crucial for performing a comprehensive study of nuclear forces and properly assessing the theoretical uncertainty of the calculation.

Light nuclei, i.e., those with $A\leq 12$, where $A$ is the number of nucleons, have proven to be an effective laboratory to test a variety of nuclear interaction models. In this realm, quantum Monte Carlo (QMC) methods have been extensively used to compute binding energies for both the ground- and the low-lying excited states at $\simeq 1\%$ accuracy level (see Ref. Carlson et al. (2015) for a recent review).

The definition of the potential describing three-nucleon ($3N$) interactions is a central issue in nuclear theory. These forces are known to yield attractive contributions to the energy per particle of light nuclei. On the other hand, a repulsive contribution is needed for the stability of neutron stars against gravitational collapse and to reproduce the equilibrium properties of isospin-symmetric nuclear matter (SNM) Akmal et al. (1998); Carbone et al. (2014); Logoteta et al. (2016).

The most accurate phenomenological Hamiltonian for $A\leq 12$ nuclei comprises the Argonne $v_{18}$ (AV18) Wiringa et al. (1995) two-nucleon ($N\!N$) potential and the Illinois-7 (IL7) Pieper et al. (2001); Pieper (2008a) $3N$ potential. This provides a good description of the spectrum of nuclei up to 12C Pieper (2008b) but yields a pathological equation of state of pure neutron matter Maris et al. (2013). On the other hand, when constraints on the $3N$ interaction are inferred from saturation properties of symmetric nuclear matter, the resulting predictions for neutron stars are compatible with astrophysical observations Gandolfi et al. (2012); Steiner and Gandolfi (2012). However $p$-shell light nuclei turn out to be underbound compared to experiment by about $0.25-0.75\,{\rm MeV}/A$ Pieper et al. (2001).

Elucidating the role of $3N$ forces in the region of medium-mass nuclei, such as 16O and 40Ca, is of paramount importance. Studying these two nuclei will help us to understand the mass region where the $3N$ contribution might already become repulsive. This aspect is strongly connected to the long-standing problem of the oxygen and calcium drip lines, which will be a major experimental focus of the Facility for Rare Isotope Beams Facility for Rare Isotope Beams .

An accurate description of 16O, in particular its interaction with neutrinos, is also of immediate importance for the detection of supernova neutrinos Ankowski et al. (2016). The large water-Cherenkov detectors require precise determination of their backgrounds, especially the one involving neutron knockout through neutral-current scattering of atmospheric neutrinos on 16O Ankowski and Benhar (2013). The computation of the electromagnetic responses of 16O using realistic nuclear interactions is a first step in this direction. In addition, studying the Coulomb sum rules of both 16O and 40Ca allows the investigation of putative in-medium modifications of the nucleon electromagnetic form factors Cloët et al. (2016).

Highly advanced nuclear many-body techniques, such as the coupled cluster method Hagen et al. (2014), the no-core shell model Barrett et al. (2013), the similarity renormalization group Hergert et al. (2016), and the self-consistent Green’s function Dickhoff and Barbieri (2004), have been successfully employed to study oxygen and calcium isotopes. In this work we use nuclear quantum Monte Carlo methods, which are capable of dealing with a wider range of momentum and energy, and allow the use of nuclear interactions characterized by high-momentum components.

Standard quantum Monte Carlo techniques, namely variational Monte Carlo (VMC) and Green’s function Monte Carlo (GFMC), work in the complete spin-isospin space, which grows exponentially with $A$ Carlson et al. (2015). As a consequence, these methods are currently limited to $A\leq 12$ nuclei by available computational resources. Over the last two decades, the auxiliary field diffusion Monte Carlo (AFDMC) method Schmidt and Fantoni (1999); Gandolfi et al. (2014a); Carlson et al. (2015), which uses Monte Carlo to also sample the spin-isospin degrees of freedom, has emerged as a more efficient algorithm for dealing with larger nuclear systems, but so far only for somewhat simplified interactions. Within cluster variational Monte Carlo (CVMC) Pieper et al. (1990, 1992), expectation values are evaluated with a cluster expansion for the spin-isospin dependent correlations. The cluster expansion drastically reduces the computational effort necessary for the study of an $A$-body system, and it enables the study of medium-mass nuclei. Another approach based on a cluster expansion of nuclear correlations has been recently used to study the high-momentum components of nuclear wave functions (see Alvioli et al. (2016) and references therein). This work, not based on Monte Carlo techniques, has been carried out employing two-body nuclear interactions only and limiting the cluster expansion to the leading order.

In this work we employ CVMC to perform variational calculations of three closed-shell nuclei, 4He, 16O, and 40Ca. We use as input a realistic phenomenological Hamiltonian, capable of describing the nucleon-nucleon data, both in scattering and bound states, with remarkable accuracy. The binding energy of the $3N$ system and the saturation density of isospin-symmetric nuclear matter are also well reproduced. We present results for the binding energy, charge radius, point density, single-nucleon momentum distribution, charge form factor, and Coulomb sum rule, fully taking into account the high-momentum components of the nuclear interaction.

In Sec. II we briefly introduce the nuclear Hamiltonian and many-body wave functions used here. Section III is devoted to the description of the cluster variational Monte Carlo technique. In Sec. IV we present our results for 4He, 16O, and 40Ca. Finally, our conclusions are summarized in Sec. V.

II Nuclear Hamiltonian and wave functions

Over a substantial range of energy and momenta, atomic nuclei can be described as collections of point-like particles of mass $m$, whose dynamics is dictated by a nonrelativistic Hamiltonian

[TABLE]

Phenomenological $N\!N$ potentials include electromagnetic and one-pion-exchange terms at long range, and parametrize the intermediate- and short-distance region with phenomenological contributions that reproduce nucleon-nucleon elastic scattering data up to the pion-production threshold:

[TABLE]

A standard version in this class of potentials is the AV18 Wiringa et al. (1995) interaction. In AV18, the electromagnetic term $v_{ij}^{\gamma}$ includes one- and two-photon-exchange Coulomb interactions, vacuum polarization, Darwin-Foldy, and magnetic moment terms, with appropriate form factors that keep terms finite at $r_{ij}=0$, where $r_{ij}=|{\bm{r}}_{i}-{\bm{r}}_{j}|$ is the interparticle distance. The one-pion-exchange and phenomenological contributions can be written as a sum of 18 operators,

[TABLE]

The first six operators, corresponding to the static components of the $N\!N$ interaction, are

[TABLE]

where $\bm{\sigma}_{i}$ and $\bm{\tau}_{i}$ are Pauli matrices acting in spin and isospin space, respectively, and

[TABLE]

is the tensor operator. The operators $p=7,\ldots,14$ are associated with the non-static components of the $N\!N$ force. They have the form

[TABLE]

where

[TABLE]

are the relative angular momentum and the total spin of the pair $ij$, respectively. Overall, the first 14 operators of AV18 describe the charge-independent part of the $N\!N$ interaction. The last four operators account for small violations of isospin symmetry, and are grouped into charge-dependent $(p=15,17)$ and charge-symmetry breaking $(p=18)$ components,

[TABLE]

where $T_{ij}=3\,\tau_{z_{i}}\tau_{z_{j}}-\bm{\tau}_{i}\cdot\bm{\tau}_{j}$ is the isotensor operator.

AV18 fits the 1993 Nijmegen database Stoks et al. (1993), which includes $4301$ $N\!N$ scattering data up to $E_{\rm lab}=350\,\rm MeV$, with a $\chi^{2}/N_{data}\simeq 1.1$, as well as the deuteron binding energy and $nn$ scattering length. It is also found to be qualitatively good to much higher energies (up to $600\,\rm MeV$) Gandolfi et al. (2014b).

The inclusion of the $3N$ interaction $V_{ijk}$ is needed to explain the binding energies of the $3N$ systems and the saturation properties of SNM. The derivation of $V_{ijk}$ was first discussed in the pioneering work of Fujita and Miyazawa Fujita and Miyazawa (1957), who argued that its main contribution originates from the two-pion-exchange process in which the $N\!N$ interaction leads to the excitation of one of the participating nucleons to a (virtual) $\Delta$ resonance, which then decays by interacting with a third nucleon.

In this work we use a phenomenological model of the $3N$ force, namely the Urbana IX (UIX) potential Pudliner et al. (1995), which is written as a sum of three contributions:

[TABLE]

The Fujita-Miyazawa anticommutator and commutator terms are

[TABLE]

where $cyc$ denotes a cyclic sum over the three particle indexes and

[TABLE]

with $\mu_{\pi}=m_{\pi}/\hbar c$ the pion mass, and $Y_{\pi}(x)$ and $T_{\pi}(x)$ the Yukawa and tensor Yukawa functions respectively, with cutoffs

[TABLE]

The purely phenomenological repulsive term is given by

[TABLE]

The parameters $A_{2\pi}$ and $A_{R}$ are adjusted to reproduce the ground-state energy of the $3N$ systems and the SNM saturation density when used in conjunction with the AV18 $N\!N$ interaction. The IL7 $3N$ potential also includes multi-pion-exchange components. The resulting AV18+IL7 Hamiltonian leads to predictions of $\simeq 100$ ground- and excited-state energies up to $A=12$ nuclei in very good agreement with the corresponding empirical values Carlson et al. (2015). However, when used to compute the neutron star matter equation of state, IL7 does not provide sufficient repulsion to guarantee the stability of observed stars against gravitational collapse Maris et al. (2013). We have therefore used the simpler UIX interaction in this study.

We note that the local $N\!N$ potentials recently derived within chiral perturbation theory Gezerlis et al. (2013, 2014); Lynn et al. (2014); Piarulli et al. (2015, 2016) are written in the same fashion as in Eq. 3. Because local versions of the chiral $3N$ potentials Tews et al. (2016); Lynn et al. (2016); Logoteta et al. (2016) have spin-isospin structure analogous to that of UIX, the formalism developed in this paper can be readily applied to this class of interactions.

Variational Monte Carlo exploits the stochastic Metropolis algorithm Metropolis et al. (1953) to evaluate the expectation value of a given many-body operator using a suitably parametrized trial wave function $\Psi_{V}$. The nuclear potential introduces spin-isospin correlations into the nuclear wave function so the variational wave function should, to the extent possible, contain operator correlations of $v_{ij}$ and $V_{ijk}$. In the same spirit of Ref. Pieper et al. (1992), in this work we assume that a good variational wave function for the ground state of a closed-shell nucleus can be expressed as the product of two- and three-body correlation operators acting on a Jastrow wave function $\Psi_{J}$:

[TABLE]

In the above equations, $U_{ij}$ and $U_{ijk}$ are correlations depending upon the spin and isospin of particles $ij$ and $ijk$, respectively. The $U^{2-6}_{ij}$ are static correlations (they contain no derivatives) while $U^{7-8}_{ij}$ are $\bm{L}\cdot\bm{S}$ correlations, fully defined following Eq. 23. The first term in the parentheses comes from the approximation of the independent triplet product of $(1+U_{ijk})$ to the linear term only. The symmetrization operator $\mathcal{S}$ is needed for the wave function to be fully antisymmetric, because $[U^{2-6}_{ij},U^{2-6}_{jk}]\neq 0$. To avoid multiple-order derivatives, the spin-orbit correlations $U^{7-8}_{ij}$ are done as a sum and act first on just the Jastrow wave function. In the Jastrow wave function, $f_{c}(r_{ij})$ denotes a central pair correlation function, $\mathcal{A}$ is the antisymmetrization operator, and $\Phi$ is an independent-particle wave function.

For doubly closed-shell nuclei, we can use a single product of four determinants $D_{\tau\sigma}$, one each for protons and neutrons, spin up and spin down, for $\Phi$:

[TABLE]

where each determinant contains $A/4$ nucleons. It follows that $\mathcal{A}|\Phi\rangle$ of Eq. 18 is a sum over all the possible partitions of the $A$ nucleons into four groups of $A/4$ nucleons.

Each determinant is constructed from single-particle radial wave functions

[TABLE]

calculated on the relative coordinates $\tilde{\bm{r}}_{i}$,

[TABLE]

in order to make $\Phi$ translationally invariant. $Y_{lm}(\theta,\varphi)$ is the spherical harmonic. The radial wave functions $R_{nl}(r)$ are obtained from the bound-state solutions of the Woods-Saxon wine-bottle potential,

[TABLE]

where the five parameters $V_{S}$, $R_{s}$, $a_{s}$, $\alpha_{s}$, and $\rho_{s}$ are determined variationally.

As stated above, the two-body correlation operator $U_{ij}$ should reflect the spin-isospin structure of the underlying $N\!N$ potential. In this work we consider only the first eight spin-isospin operators, which capture the dominant features in the $N\!N$ phase shifts,

[TABLE]

with $u_{p}(r_{ij})=f_{p}(r_{ij})/f_{c}(r_{ij})$. The radial correlation functions $f_{c,p}(r_{ij})$ are obtained by minimizing the two-body cluster contribution to the energy per particle of SNM at the Fermi momentum $k_{F}$. Euler-Lagrange (EL) equations are solved in a partial-wave $(S,T)$ basis for a quenched potential,

[TABLE]

by imposing the boundary conditions Lagaris and Pandharipande (1981):

[TABLE]

In the present calculations we assume

[TABLE]

and we consider three independent healing distances,

[TABLE]

where $d_{\rm S}\neq d_{\rm P}$ are used in order to differentiate $s$-wave (${}^{1}S$ and ${}^{3}S-{{}^{3}D}$) from $p$-wave (${}^{1}P$ and ${}^{3}P-{{}^{3}F}$) channels, and the general relation $d_{\rm S}<d_{\rm P}<d_{t}$ should hold. The functions $f_{c,p}(r_{ij})$ are projected from the solutions of the ($S,T$) partial-wave EL equations. The pair correlation functions are thus fully specified by a total of seven variational parameters: $k_{F}$, $\alpha$, $\beta_{c}$, $\beta_{t}$, $d_{\rm S}$, $d_{\rm P}$, and $d_{t}$.

In a many-body system it has been found advantageous to screen the spin- and isospin-dependent pair correlation functions when other particles are nearby Lomnitz-Adler et al. (1981); Pudliner et al. (1997). This can be achieved by multiplying $U_{ij}$ by three-body correlation factors,

[TABLE]

where

[TABLE]

The three parameters $t_{1}$, $t_{2}$, and $t_{3}$ are found variationally.

Explicit triplet correlations significantly improve the variational energy for Hamiltonians including a $3N$ interaction. In this work we employed the form

[TABLE]

where $\tilde{V}_{ijk}$ have the structures of Eqs. 10 and 16 but the two-particle distances are rescaled by a factor $\eta$, and two different constants $c_{y}$ and $c_{t}$ are used for the cutoff function $\xi(r)$ of Eq. 15 used in Eqs. 13 and 14. The triplet correlation functions are then given in terms of five variational parameters: $\varepsilon_{2\pi,A}$, $\varepsilon_{R}$, $\eta$, $c_{y}$, and $c_{t}$.

We did not include correlations arising from the commutator of Eq. 11 because it is significantly more computationally expensive to evaluate than the anticommutator of Eq. 10. However, it has been shown that most of the correlations induced by the commutator can be effectively obtained by an appropriate choice of the coefficient $\varepsilon_{2\pi,A}$ Pudliner et al. (1997).

III Cluster variational Monte Carlo

In VMC, once the form for the trial wave function is assumed, one optimizes the variational parameters, typically by minimizing the expectation value and/or the variance of the total energy with respect to the variations of the parameters. The energy expectation value is given by

[TABLE]

and it is always greater than or equal to the ground-state energy with the same quantum numbers as $\Psi_{V}$. By minimizing $E_{V}$ the optimal $\Psi_{V}$ is obtained, and it is used to evaluate other quantities of interest.

In general, for spin-isospin dependent interactions, the wave function is a sum of complex amplitudes for each spin-isospin state. The number of these components grows exponentially with the number of particles. This scaling can be mitigated by considering charge conservation and by assuming that the nucleus has good isospin $T$. However, for $A>12$ nuclei, quantum Monte Carlo calculations employing the complete many-body wave function currently represent a computational challenge Carlson et al. (2015).

One way to overcome the scaling problem and perform calculations for larger systems is to employ a cluster expansion scheme. The expectation value $\langle\Psi_{V}|H|\Psi_{V}\rangle$ as well as $\langle\Psi_{V}|\Psi_{V}\rangle$ can be expanded according to the number of nucleons connected by the spin-isospin correlations $U_{ij}$ and $U_{ijk}$. The resulting cluster expansion for the expectation value $E_{V}$, which is constructed according to Ref. Pandharipande and Wiringa (1979), has been used up to four-body cluster for the VMC study of 16O Pieper et al. (1992) and ${}_{\leavevmode\nobreak\ \Lambda}^{17}$O Usmani et al. (1995) with earlier versions of the phenomenological $N\!N$+$3N$ potentials. In this work the calculations have been performed including up to five-body cluster contributions and considering closed-shell nuclei as large as 40Ca. The modern AV18 $N\!N$ potential plus the UIX $3N$ force has been employed.

III.1 Cluster expansion

The trial wave function of Eq. 17 contains a large number of terms because there are many ways of partitioning $A$ nucleons into four groups of $A/4$ nucleons that preserve the antisymmetrization of $\Psi_{V}$. However, since $H$ is a symmetric operator, we can reduce the problem by considering a trial wave function $\Psi_{R}$ not fully antisymmetric,

[TABLE]

and by re-defining the energy expectation value as

[TABLE]

The cluster expansion adopted in this work is the one associated with expectation values of the form (34). In the reference work Pieper et al. (1992) this cluster expansion is referred to as “CEA.”

Let us consider the expectation value of a symmetric one-body operator $\mathcal{O}_{i}$:

[TABLE]

The numerator $N$ and denominator $D$ can be expanded as a sum of $n$-body contributions,

[TABLE]

Obviously extending the sums to $A$-body contributions gives the exact expectation value. We define the generic expectation value $\langle X\rangle$, to be used for both $N$ and $D$ terms in Eq. 35, as

[TABLE]

The contributions $n_{ij\ldots}$ and $d_{ij\ldots}$ then take the following form:

[TABLE]

The expansions (36) and (37) for $N$ and $D$ are divergent. On the other hand, a convergent expansion is achieved by considering the linked cluster expansion

[TABLE]

whose coefficients can be obtained from the equation $C\cdot D=N$ by equating terms containing the same number of particles,

[TABLE]

The cluster expansion for the expectation value of two-body operators $\mathcal{O}_{ij}$ and three-body operators $\mathcal{O}_{ijk}$, such as $v_{ij}$ and $V_{ijk}$, resembles the one for the one-body operator $\mathcal{O}_{i}$. However, in the case of $\mathcal{O}_{ij}$, there are no one-body terms $n_{i}$, nor terms such as $n_{i,jk}$ in the numerator (36). Therefore the cluster expansion (41) only contains terms of the kind $c_{ij},\,c_{ijk},\,c_{ij,kl},\,c_{ijkl},\,\ldots$. In a similar fashion, the cluster expansion for $\sum_{ijk}\mathcal{O}_{ijk}$ only comprises terms like $c_{ijk},\,c_{ijkl},\,c_{ijk,lm},\,c_{ijklm},\,\ldots$.

Terms such as $c_{i,jk}$ are referred to as semifactorizable. They are typically small because of the large cancellation between $n_{i,jk}$ and $c_{i}d_{jk}$, but they are finite. It is not necessary to treat them separately from the others. For example it is possible to define cluster contributions $\tilde{c}_{ijk}$ as the sum of all those that contain particles $ijk$ so that

[TABLE]

The corresponding $\tilde{n}_{ijk}$ can also be directly computed without separating their semifactorizable contributions. The total $n$-body cluster contribution $C_{n}$ is then obtained from the sum

[TABLE]

and Eq. 41 can be simply rewritten as

[TABLE]

In the present work the cluster expansion is carried out up to five-body cluster, $n=5$. Since the operators in the expectation value $n_{ij\ldots l}$ or $d_{ij\ldots l}$ only contain the spin and isospin of particles $ij\ldots l$, the spin and isospin of the other particles are unchanged and can be ignored. If $ij\ldots l$ are in a single determinant $D_{\tau\sigma}$ in $|\Phi\rangle$, then only the term $\langle\Phi|$ in $\langle\Phi|\mathcal{A}$ contributes, and the rest can be ignored. If $i$ is in $D_{\tau^{\prime}\sigma^{\prime}}$ and $j\ldots l$ are in $D_{\tau\sigma}$ in $|\Phi\rangle$, then only the direct term $\langle\Phi|$ and those obtained by exchanging $i$ with $j\ldots l$ in $\langle\Phi|\mathcal{A}$ need to be considered. This implies a large reduction of the number of contributions to be calculated at each order, allowing for a full evaluation up to five-body cluster.

All the expectation values $n_{i},\,n_{ij},\,\ldots$ and $d_{ij},\,d_{ijk},\,\ldots$ are calculated up to four-body cluster. Five-body cluster contributions are instead sampled according to the probability

[TABLE]

where $x=\sum_{i<j}r_{ij}$, and typical values are $P_{min}=0.02$, $b=35\,\rm fm$, and $a=3.2\,\rm fm$. If $P(x)$ is larger than $\xi$, where $\xi$ is a random number in the interval $[0,1]$, then the five-body contribution is calculated. For 16O it has been verified that sampling five-body cluster terms yields an energy expectation value that is compatible to the one obtained with the full five-body cluster calculation $(P_{\min}=1)$. In 16O the sampling procedure speeds up the evaluation of $E_{V}$ by a factor of 1.7 when using the $N\!N$ potential only, and by a factor of 2.2 when also $3N$ interactions are included. This is crucial for the calculation of 40Ca, in particular when using the full AV18+UIX potential. In 16O there are 4368 quintuplets, while in 40Ca there are 658008 quintuplets, making the full five-body cluster calculation extremely time demanding.

Further simplifications can be made by looking at the structure of the employed trial wave function. $|\Phi\rangle$ is a product of four determinants in which particle $(1:d)$, $(d+1:2d)$, $(2d+1:3d)$, and $(3d+1:4d)$, with $d=A/4$ are, respectively, $p\uparrow$, $p\downarrow$, $n\uparrow$ and $n\downarrow$. $\langle\Phi|\mathcal{A}$ is instead fully antisymmetric, so that when particle $i$ and $i^{\prime}$ belong to the same determinant, the following equivalences among expectation values apply:

[TABLE]

By neglecting the effects of the Coulomb potential on the wave function, for the isospin-symmetric nuclei considered in this work it follows that, for instance, there are only four nonequivalent classes of $n_{ij}$ contributions:

[TABLE]

We note that the employed cluster expansion treats exactly all the exchanges and central correlations among the $A$ nucleons. Every term in the cluster expansion (III.1) and (40) contains the complete product of central correlations. In the conventional cluster expansions Pandharipande and Wiringa (1979), one also expands in powers of $f_{c}^{2}(r)-1$ and this does not necessarily keep all the exchange terms.

The current work includes the $\bm{L}\cdot\bm{S}$ correlations and $\bm{L}\cdot\bm{S}$, $\bm{L}^{2}$, and $(\bm{L}\cdot\bm{S})^{2}$ potentials in all cluster expansion orders. Reference Pieper et al. (1992) included these in only the two-body clusters, arguing that their total contribution is small. However we find a large, repulsive, three-body contribution from these potential terms.

Note that in the process of expanding the numerator and the denominator of the Hamiltonian’s expectation value of Eq. 34, the variational principle is not guaranteed to hold. However, since summing up to the $A$-body contribution gives the exact expectation value, the convergence of the cluster expansion itself will restore the validity of the variational principle. For this reason, during the optimization of the variational parameters, the convergence of the cluster expansion has been carefully checked for each of the analyzed cases.

III.2 VMC sampling

The spatial integrals in Eq. 38 are evaluated using Metropolis Monte Carlo techniques Metropolis et al. (1953). The Metropolis method allows one to sample points in large-dimensional spaces according to a probability distribution $W(\bm{R})$, where $\bm{R}=\{\bm{r}_{1},\ldots,\bm{r}_{A}\}$. The algorithm generates a sequence of points (random walk) in the $3A$-dimensional space. This is achieved by a sequence of moves that can either be accepted or rejected depending upon the ratio of the function $W$ computed at the original and proposed points. According to the central limit theorem, the generic expectation value $\langle I\rangle$ can be written as

[TABLE]

where $\mathcal{N}_{c}$ is the number of configurations $\bm{R}_{i}$ sampled with probability proportional to $W(\bm{R})$. The Monte Carlo statistical error associated to $\langle I\rangle$ can be estimated with $\epsilon_{I}=\sqrt{\sigma_{I}/\mathcal{N}_{c}}$, where $\sigma_{I}$ is the variance of $I$.

The weight function $W(\bm{R})$ must be positive definite and normalizable. The choice adopted in this work is to use the Jastrow part of the trial wave function $\Psi_{R}$

[TABLE]

The expectation value $\langle X\rangle$ is

[TABLE]

and the function to evaluate at a sampled $\bm{R}_{i}$ is $\Phi^{*}(\bm{R})\mathcal{A}\,X\,\Phi(\bm{R})$ [with the normalization factor $|\Phi(\bm{R})|^{2}F(\bm{R})$], where the spin-isospin summations are implicit. In the present case $\Phi(\bm{R})$ is real, so that $\Phi^{*}(\bm{R})=\Phi(\bm{R})$.

The factor $F(\bm{R})$ is introduced in the weight function $W(\bm{R})$ in order to prevent the quantity $\Phi^{*}(\bm{R})\mathcal{A}\,X\,\Phi(\bm{R})/|\Phi(\bm{R})|^{2}$ from becoming very large. It is chosen so that $\Phi^{*}(\bm{R})\mathcal{A}\,X\,\Phi(\bm{R})/\!\left[|\Phi(\bm{R})|^{2}F(\rm R)\right]$ is finite at all $\bm{R}$. All the exchanges that contribute to $\Phi^{*}(\bm{R})\mathcal{A}\,X\,\Phi(\bm{R})$ are included in $|\Phi(\bm{R})|^{2}F(\bm{R})$ so that

[TABLE]

where $\mathcal{P}_{ij}$ is the exchange operator acting on particles $i$ and $j$. In the present work contributions up to four-body exchanges are considered in Eq. 52. The function $\omega(r)$ is chosen to be proportional to the sum of the squares of $u_{p}(r)$, since the exchange of particles $i$ and $j$ with different spin-isospin states in $\Phi(\bm{R})$ must be accompanied by a $v_{ij}$, $U_{ij}$, $V_{ijk}$, or $U_{ijk}$. The use of the importance function $F(\bm{R})$ drastically reduces the variance on the expectation values. For instance, in 16O the same statistical error for the energy expectation value can be achieved by using just half of the configurations when $F(\bm{R})\neq 1$.

The $\Phi(\bm{R})\mathcal{A}\,X\,\Phi(\bm{R})$ for a given cluster is calculated with methods developed for few-body systems Wiringa (1991); Lomnitz-Adler et al. (1981). The terms in $\Phi(\bm{R})\mathcal{A}$ that can contribute are summed, and $\Phi(\bm{R})\mathcal{A}$ is represented as a vector whose components give the amplitudes of the spin-isospin states of the nucleons in the cluster. The corresponding vector representing $\Phi(\bm{R})$ has only one nonzero component since all particles have definite values of $\tau_{z}$ and $\sigma_{z}$ in $\Phi(\bm{R})$. The $v_{ij}$, $V_{ijk}$, $U_{ij}$, and $U_{ijk}$ operate on these vectors as discussed in Refs. Wiringa (1991); Lomnitz-Adler et al. (1981). The expectation values of the kinetic energy operators are obtained by computing $\Psi_{R}$ at slightly shifted positions and using finite differences to evaluate terms in $\nabla^{2}\Psi_{R}$.

Due to the tremendous increase in computer power of the last decades, many of the approximations implemented in the reference work Pieper et al. (1992) are no longer necessary. For instance, in the current calculations the three-body correlation operators $U_{ijk}$ act last in $\Psi_{R}$ and $\Psi_{V}$ of Eqs. 17 and 32, as in the original formulation of the trial wave function. In Ref. Pieper et al. (1992), because of the computational limitations of the time, $\Psi_{R}$ and $\Psi_{V}$ were approximated acting first with the $U_{ijk}$ on the sparse vectors representing $\Phi(\bm{R})$ and $\mathcal{A}\Phi(\bm{R})$, and then operating with the two-body correlations $U_{ij}$. The latter are now implemented in all orders, including the spin-orbit correlations that were previously calculated at the two-body level only.

Moreover, the calculation of the contribution of the kinetic energy, $N\!N$, and $3N$ potential operators is fully carried out at each order of the cluster expansion. For the five-body cluster, all the one-, two-, and three-body operators are evaluated, although their contributions to $C_{5}$ are sampled as previously discussed.

III.3 Optimization

The trial wave functions of Eqs. 32 and 17 contain a total of 15 variational parameters when only two-body correlations are considered, and up to 20 parameters if three-body correlations are also included. In order to perform the minimization of the energy expectation value with respect to these sets of parameters, we used the NLopt optimization tool, as recently done in other standard VMC calculations Piarulli et al. (2016).

NLopt is a free/open-source library for nonlinear optimization developed at the Massachusetts Institute of Technology Johnson . It provides a common interface for a number of different free optimization routines available online as well as original implementations of various other algorithms, including both global and local optimization algorithms, both derivative-free and user-supplied gradients algorithms, and algorithms for unconstrained optimization, bound-constrained optimization, and general nonlinear inequality/equality constraints.

In this work we implemented different local derivative-free algorithms, and in particular we made extensive use of the COBYLA (constrained optimization bY linear approximations) Powell ; *powell:1998 and Nelder-Mead simplex Nelder and Mead (1965); *box:1965; *richardson:1972 algorithms. It has been observed that both algorithms perform well in the case of 4He Piarulli et al. (2016). For heavier systems, Nelder-Mead simplex seems instead to be the optimal algorithm, providing better convergence and reliability of the minimization search. This is probably related to the fact that the minimization was done using correlated energy differences Wiringa (1991) for $A=4$ but not for the larger nuclei.

For both 4He and 16O the energy minimization was carried out in the full parameter space, with the energy expectation value calculated up to the highest cluster contribution for the system under study. In order to reduce the computational cost of the optimization process, the spin-orbit correlations are turned off during the variational search. However, once the optimal set of parameters is found, the full two-body correlations of Eq. 23 are employed in the calculation of the expectation values. In the case of 40Ca, the computation of five-body cluster contributions to the total energy is quite demanding, even with the sampling procedure. Each CVMC run for $A=40$ requires approximately two hours on 18 32-core Intel Haswell 2.3 GHz nodes to obtain a statistical error of $\simeq 0.5\,{\rm MeV}/A$ for the energy of the full $N\!N$ plus $3N$ Hamiltonian. The variational search over the entire 20-dimensional parameter space would have required at least $\sim 140$ h on the same hardware configuration, i.e., more than $80\cdot 10^{3}$ CPU hours.

Relying on the observation that short-range correlations for medium-heavy systems should be independent of $A$, the energy minimization for 40Ca has been carried out in a subset of the parameter space. When using AV18 only, the optimal parameters for the two-body correlations found in 16O for the same potential were employed, as were the wine-bottle coefficients and the induced three-body correlations of Eq. 28. The variational search was performed for only the three parameters defining the Wood-Saxon potential of Eq. 22. When the $3N$ potential is also included, from the best set of parameters for 16O with the same interaction, we minimized over the three parameters of the Wood-Saxon potential and over two of the five parameters of the three-body correlations. For the latter, $c_{y}$ and $c_{t}$ appear to be the most effective to produce appreciable changes in the total energy. All the parameters for the systems under study for both AV18 and AV18+UIX are listed the Appendix.

IV Results

The expectation values of all observables are calculated for each nucleus by summing all cluster contributions up to five-body cluster (four-body in the case of 4He). For 16O and 40Ca, the full expansion should consider contributions up to 16- and 40-body clusters, respectively. Under the observation that the ratio between the last successive cluster contributions is small and approximately constant, we can estimate the cluster contribution $C_{6-A}$ by assuming uniform convergence, i.e., by using the relation

[TABLE]

Given the cluster expansion at order $k$, the total extrapolated result is obtained by summing over all the cluster contributions, including the extrapolated ones

[TABLE]

Under the assumption of uniform convergence, the $C_{n}^{\rm ext}$ form a geometric progression, and we can then recast the total $C_{\rm ext}$ using the sum of the geometric series

[TABLE]

Note that in the employed cluster expansion, successive cluster contributions $C_{k}$ and $C_{k-1}$ have decreasing magnitude and opposite sign, so that $|x|<1$.

Equations 54 and 55 give consistent results for all the observables under study. In the following, unless otherwise specified, we will report results for 16O and 40Ca using the extrapolation of Eq. 54 for contributions above the five-body cluster. Errors on $\sum_{n=k+1}^{A}C_{n}^{\rm ext}$ are estimated by propagating the CVMC statistical errors from the previous cluster contributions.

IV.1 Energies, radii, and densities

The contributions of the cluster expansion to the kinetic energy $T$, to the $N\!N$ and $3N$ potentials, and to the point radius are listed in Tables 1 and 2 for 4He, in Tables 3 and 4 for 16O, and in Tables 5 and 6 for 40Ca. The expectation value of $v_{ij}^{15-18}$ is zero for all the systems under study, as we are assuming pure $T=0$ ground states. Since one-, two-, and three-body operators exhibit different convergence patterns in the cluster expansion, for $A>4$ the total energy is estimated as the sum of the extrapolated results for $T$, $v_{ij}$, and $V_{ijk}$, and it is italicized in the tables.

Let us first consider the 4He nucleus. By comparing the energy per nucleon obtained with AV18 and AV18+UIX Hamiltonians, both reported in Table 7, it is apparent that the $3N$ force gives overall $\simeq 1\,{\rm MeV}/A$ more binding. This result is consistent with VMC and GFMC calculations Pieper et al. (2001); Pieper and Wiringa (2001) for the same interactions. It is interesting to note that the best wave function for the full Hamiltonian including UIX sacrifices $\simeq 0.3\,{\rm MeV}/A$ from the $T+v_{ij}$ contribution which is made up by increasing the attraction from UIX.

The more sophisticated two-body correlations employed in the VMC wave function for $s$-shell nuclei yield $\simeq 0.2\,{\rm MeV}/A$ additional binding compared to the CVMC results. Nevertheless, the CVMC energies are within $5\%$ from the GFMC values, and charge radii are remarkably close for all the three quantum Monte Carlo methods. This corroborates the accuracy of the wave functions employed in this work to describe the ground state of closed-shell nuclei, which, combined with the cluster expansion technique, allows reliable variational calculations for nuclei as heavy as 40Ca.

The total energies of 16O and 40Ca for both AV18 and AV18+UIX are reported in Tables 9 and 9, respectively. Our variational calculations show that the AV18+UIX Hamiltonian underbinds both 16O and 40Ca, by $2.83(3)\,{\rm MeV}/A$ and $3.63(10)\,{\rm MeV}/A$, respectively. The results obtained for 40Ca are consistent with variational calculations for SNM performed with the same interaction Akmal et al. (1998), which yield $-11.85\,{\rm MeV}/A$, to be compared to the empirical value of $\simeq-16\,{\rm MeV}/A$. This underbinding can be only partly ascribed to deficiencies of the variational wave function, which has proven to be accurate for describing infinite matter properties Lovato et al. (2011). To gauge the accuracy of the CVMC wave function in describing closed-shell nuclei, we performed a benchmark calculation with AFDMC using the AV6′ potential. This is a reprojection of the full AV18 onto the first six operators that preserves the deuteron binding energy and many of the properties of elastic $N\!N$ scattering Wiringa and Pieper (2002). To obtain 16O energies that are bound against $\alpha$-particle break up, the Coulomb interaction was omitted. The results listed in Table 10 show a $\simeq 0.25\,\rm MeV/A$ and a $\simeq 0.45\,\rm MeV/A$ energy difference in 4He and 16O respectively between CVMC and AFDMC results. This is expected for a variational versus a diffusion Monte Carlo calculation. Charge radii are instead compatible between the two methods, confirming the quality of the employed wave function. Therefore, a large fraction of the missing binding in 16O and 40Ca is due to limitations of the AV18+UIX Hamiltonian. Note that these results show significantly less binding per nucleon for both 16O and 40Ca than for 4He, i.e., they predict that 16O and 40Ca would break apart into 4He nuclei. However, Table 10 indicates that 16O is stable against breakup with the AV6′ interaction if the Coulomb interaction is omitted.

For both 16O and 40Ca the expectation value of $V_{ijk}^{2\pi}$ is negative, and that of $V_{ijk}^{R}$ is positive, leading to an overall attractive contribution of the $3N$ force, as for 4He. However, by comparing the total energies for AV18 and AV18+UIX, it turns out that 16O and 40Ca are less bound when the $3N$ force is included. This is particularly evident in 40Ca, where the UIX potential reduces the binding energy of $\simeq 1\,{\rm MeV}/A$. This is somewhat consistent with the fact that the UIX force is repulsive in SNM. Finally, it is interesting to notice how, within a variational approach, the change in the behavior of the employed $3N$ force—from attractive to repulsive—is already manifest in relatively small nuclear systems, like 16O.

Three-nucleon forces significantly affect quantities other than the energy, such as point radii and point densities. The latter are related to the charge density, which can be extracted from electron-nucleus scattering data, but they are not observables themselves, as many-body currents and single-nucleon electromagnetic form factors need to be accounted for.

Neglecting small effects, the charge radius $\left\langle r_{\rm ch}^{2}\right\rangle$ can be expressed in terms of the point proton radius $\left\langle r_{\rm pt}^{2}\right\rangle$ Friar et al. (1997)

[TABLE]

where $\left\langle R_{p}^{2}\right\rangle=0.770(9)\,\rm fm^{2}$ is the proton radius Beringer et al. (2012), $\left\langle R_{n}^{2}\right\rangle=-0.116(2)\,\rm fm^{2}$ is the neutron radius Beringer et al. (2012), and $3\hbar^{2}/(4M_{p}^{2}c^{2})\simeq 0.033\,\rm fm^{2}$ is the Darwin-Foldy term. Charge radii in 4He, 16O, and 40Ca for both AV18 and AV18+UIX are reported in Tables 7, 9 and 9, respectively. In 4He AV18 produces a charge radius larger than the experimental value. However, the $3N$ force shrinks the nucleus, improving the agreement with experiment. Both CVMC values are reasonably consistent with the VMC ones. In 16O and 40Ca instead, the $N\!N$ interaction alone results in too small radii, while the UIX potential increases them towards and above their experimental values. This is consistent with the observation that, as opposed to 4He, for $A=16$ and $A=40$ the net effect of the UIX potential is to make the systems more loosely bound.

The single-nucleon, two-nucleon, and two-nucleon operator point densities are defined as

[TABLE]

where $N=p,n$, $\mathcal{P}_{N_{i}}=(1\pm\tau_{z_{i}})/2$ are isospin projection operators, and the operators $\mathcal{O}_{ij}^{p}$ are given in Eq. 4. With these definitions, $\rho_{N}(r)$ is normalized to the number of protons or neutrons, and $\rho_{\rm{NN}}$ to the number of $pp$, $np$ or $nn$ pairs. Note that for some of the alternative expansion schemes normalization is ensured order-by-order by construction (either defining a “number conserving” expansion Alvioli et al. (2005, 2008, 2013, 2016), or requiring a normalization factor Ryckebusch et al. (2015)). In our expansion scheme, the central one- and two-body densities are properly normalized order-by-order. The first term of the corresponding cluster expansion carries the full normalization, and higher order contributions integrate to zero within Monte Carlo statistical errors. This reflects the fact that at every order a $3A$ dimensional integral is performed. The normalization of the two-body operator densities is instead recovered only at convergence, and each order of the cluster expansion contributes to it.

Figures 1, 2 and 3 show the point proton densities of 4He, 16O, and 40Ca, respectively, obtained with the AV18 and AV18+UIX interactions. They are compared to the values obtained from the “Sum-of-Gaussians” parametrization of the charge densities given in Ref. de Vries et al. (1987) by unfolding the nucleon form factors and subtracting the small contribution of the neutrons. As discussed at length in Sec. IV.3, neglecting two-body meson exchange currents (MECs) is likely to have little effect in 16O and 40Ca. On the other hand, MECs are important in the description of the 4He elastic form factor, from which the charge densities are extracted. Hence, the discrepancy between theory and experiment of Fig. 1 does not have to be ascribed to deficiencies of the CVMC wave function. In fact, the 4He point proton density obtained within CVMC for AV18+UIX agrees very well with the GFMC result for the same interaction.

For the lightest system the effect of the $3N$ force on the density is not dramatic, as expected by looking at the small difference in the charge radii of Table 7. In oxygen and calcium, instead, the addition of the UIX potential pushes the nucleons far away from the center of mass. For both systems the density at small distances is substantially depleted, with a $\simeq 25\%$ reduction of both the peak in 16O at $1.4\,\rm fm$ and the plateau in 40Ca around $2\,\rm fm$. Remarkably, this effect results in a better description of the structure of 16O, for which the AV18+UIX prediction of the charge radius is less than $2\%$ different from the experimental value, as shown in Table 9. However, the situation is different in the case of 40Ca, for which the employed $3N$ force is too repulsive and pushes the nucleons towards the surface of the nucleus yielding an excessively large charge radius, as in Table 9. The point proton density of 40Ca turns out to be $0.8\,\rm fm^{-3}$ at $1\,\rm fm$, and $0.7\,\rm fm^{-3}$ in the plateau after $1.6\,\rm fm$. These values are consistent with the saturation density of SNM obtained employing the same Hamiltonian, while AV18 alone significantly overpredicts the saturation density Akmal et al. (1998).

It is interesting to compare the densities of these nuclei in which $\alpha$-clustering can potentially occur. In Fig. 4 we collect the CVMC results for the point proton densities of 4He, 16O, and 40Ca obtained with AV18+UIX together with those for 8Be and 12C coming from VMC calculations using the same interaction R. B. Wiringa . The 4He density shows a very large point density at small distance. When integrated over the volume, about half the nucleons reside inside $1.25\,\rm fm$, where the density is above $0.08\,\rm fm^{-3}$. The 8Be density has a low, broad peak with half the nucleons residing inside $2.25\,\rm fm$, consistent with a two-$\alpha$ cluster structure as observed in Fig. 15 of Ref. Wiringa et al. (2000). The 12C density peaks at a slightly smaller distance and noticeably higher value, with a larger dip at the center. This is consistent with a more tightly bound three-$\alpha$ cluster—either in a triangular configuration with a low-density region at the center of mass, or alternatively with one $\alpha$ in the $s$-shell and two $\alpha$’s in the $p$-shell. Similarly, 16O can be viewed as a tetrahedral four-$\alpha$ cluster with the $\alpha$’s at somewhat greater distance from the center of mass, or as one $s$-shell and three $p$-shell $\alpha$’s with a larger dip-peak difference than in 12C. The 40Ca density is more complicated, but might be thought of as two $s$-shell $\alpha$’s giving a larger central peak, while three $p$-shell and five $d$-shell $\alpha$’s give a broad shoulder at $1-3\,\rm fm$.

The two-nucleon point densities of 4He, 16O, and 40Ca are reported in Figs. 5, 6 and 7, respectively, for both AV18 and AV18+UIX. Upper and lower curves refer to $np$ and $pp$ pairs, respectively. The fact that $\rho_{\rm{NN}}$ is very small for $r\simeq 0$ is a consequence of the repulsive core of the $N\!N$ potential. As observed for the point-proton densities, the effect of the $3N$ force on the two-nucleon densities is appreciably different in light- and medium-heavy systems. In 4He the $pp$ density is almost unchanged, while the $np$ density is enhanced around the peak at $1.1\,\rm fm$. In heavier systems there is a severe depletion of both $pp$ and $np$ densities, again due to the peculiar repulsive effect of the UIX potential that tends to push nucleons apart.

Both figures and tables for the CVMC single-nucleon and two-nucleon point densities for 16O and 40Ca, together with the VMC results for $A\leq 12$, are available online R. B. Wiringa ; R. B. Wiringa .

The two-nucleon operator point densities are shown in Figs. 8, 9 and 10 for 4He, 16O, and 40Ca, respectively. It can be observed that the larger the system, the wider the range of central two-body density. In fact, the central two-body operator density is just the sum of $pp$, $np$, and $nn$ densities in Figs. 5, 6 and 7. On the other hand, spin-isospin densities are appreciably nonvanishing only for $r\lesssim 3.5\,\rm fm$, and are largely independent of the nucleus, with the position of the peaks situated around $1\,\rm{fm}$. This extends the results of Ref. Feldmeier et al. (2011), where the two-body densities normalized at short distances in $A=3$ and $A=4$ systems exhibit a universal behavior up to about $1\,\rm{fm}$ in all nuclei. Among the spin-isospin densities, $\rho_{2,\sigma\tau}$ and $\rho_{2,t\tau}$ are characterized by longer ranges and larger amplitudes, as they arise from the one-pion-exchange part of the $N\!N$ interaction. These results are qualitatively consistent with the findings of Ref. Alvioli et al. (2005), although no $3N$ forces were employed in that work. The peak values of the $\rho_{2,\sigma\tau}$ and $\rho_{2,t\tau}$ scale as $1:4:10$ for ${{}^{4}\rm He}:{{}^{16}\rm O}:{{}^{40}\rm Ca}$, or just as the number of $\alpha$-particle clusters.

IV.2 Momentum distributions

The probability of finding a proton or neutron with momentum $\bm{k}$ is proportional to the momentum distribution,

[TABLE]

which is normalized as

[TABLE]

$\mathcal{N}_{N}$ being the number of protons or neutrons ($\mathcal{N}_{p}\equiv Z$). In this work we present results for symmetric nuclei implying $n_{p}(\bm{k})=n_{n}(\bm{k})$. Equation (60) can be rewritten as

[TABLE]

The Fourier transform can be computed by Monte Carlo integration. Spatial configurations are sampled as explained in Sec. III.2. The average over all particles $i$ in each configuration is then performed, and for each particle, a grid of Gauss-Legendre points $x_{i}$ is used to compute the Fourier transform. The polar angle $d\Omega_{x}$ is also sampled by Monte Carlo integration, with a randomly chosen direction for each particle in each configuration. For all the nuclei under study we calculated $n(\bm{k})$ up to $k=10\,\rm fm^{-1}$, integrating to $x_{\max}=20\,\rm fm$ using 200 Gauss-Legendre points.

As reported in Ref. Pieper et al. (1992), with the employed expansion the three-body clusters give small contribution to the momentum distribution. In this work, $n_{N}(\bm{k})$ is evaluated up to three-body cluster and then extrapolated using Eq. 55. In order to save computing time, spin-orbit correlations are turned off in the calculation of the momentum distribution. This approximation, also used in standard VMC calculations Wiringa et al. (2014), is justified by the small effect of spin-orbit correlations on $n_{N}(\bm{k})$ compared to the first six operators of the two-body correlations. The results for 4He, 16O, and 40Ca are shown in Figs. 11, 12 and 13, respectively. For $A=4$ the VMC result for AV18+UIX Wiringa is also displayed for comparison. The proton momentum distributions are reported for both AV18 and AV18+UIX. The $3N$ force makes only small changes to $n_{N}(\bm{k})$. Near $k=2\,\rm fm^{-1}$ the momentum distribution manifests a sharp change in slope, as previously observed in both light- Wiringa et al. (2014) and medium-mass Alvioli et al. (2008) nuclei. This is attributed to the strong tensor correlations induced by the one-pion-exchange part of the $N\!N$ potential, further enhanced by the two-pion-exchange part of the $3N$ potential, when included. At higher momentum, the tail of $n_{N}(\bm{k})$ manifests the expected universal behavior determined by the short-range correlations, i.e., by the short-range structure of the employed Hamiltonian, as shown in Fig. 14 and discussed at length in a number of other works Feldmeier et al. (2011); Alvioli et al. (2012, 2013, 2016). Such universality refers to the independence of the tail with respect to the specific nucleus. On the other hand, the high-momentum tail strongly depends on the nuclear interaction model. The recently developed local chiral interactions, which are significantly softer than the phenomenological interactions employed in this work, yield a momentum distribution characterized by weaker high-momentum components Gandolfi et al. than those of Figs. 11, 12 and 13.

Compared to the other local or nearly local operators, like the kinetic energy, the potential energy, and the densities, the momentum distribution is strictly a nonlocal operator. In order to check the convergence of the cluster expansion for such operator we computed the kinetic energy by integrating the momentum distribution

[TABLE]

for each order of the expansion. The contributions for 16O with AV18+UIX up to $k=10\,\rm fm^{-1}$ are $16.3(2)\,\rm MeV/A$ for one-body cluster, $16.0(5)\,\rm MeV/A$ for two-body cluster, and $-4.4(4)\,\rm MeV/A$ for three-body cluster. The integration of the extrapolated $n_{N}(\bm{k})$ leads to $28.9(6)\,\rm MeV/A$. This is compatible with the cluster contributions reported in the first line of Table 4. The missing 4- to 16-body cluster contributions to the integrated kinetic energy, that account for $\simeq 1\,\rm MeV/A$, are fully recovered by the extrapolation of $n_{N}(\bm{k})$. This validates the convergence of the expansion and confirms the negligible effect of spin-orbit correlations on the momentum distribution. Similar outcomes are found for the other nuclei considered in this work. The errors on the integrated kinetic energies are larger than those of the direct calculation because of the propagation of uncertainties in the integration of $n_{N}(\bm{k})$, which above $5\,\rm fm^{-1}$ has large statistical errors due to the cancellation of positive and negative small cluster contributions. However, as discussed in the next paragraph, the integrated strength of the momentum distribution saturates before $5\,\rm fm^{-1}$. Simulations for $n_{N}(\bm{k})$ have been thus carried out with good statistics up to that momentum value, using up to $10^{7}$ Monte Carlo configurations.

The momentum distribution integrated strength as a function of $k$ is reported in Fig. 15 for AV18+UIX. At low momentum it decreases as $A$ increases, because the nuclei become more tightly bound, and the fraction of nucleons at low momentum decreases. At $k=2\,\rm fm^{-1}$ for all the systems analyzed, the integrated strength is already $\simeq 93\%$ of the total, and it becomes $\simeq 99\%$ at $k=4\,\rm fm^{-1}$. Less than $1\%$ of the total strength is given by the tail of the momentum distribution above $4\,\rm fm^{-1}$.

The figures and the tables for the CVMC momentum distributions for 16O and 40Ca, together with the VMC results for $A\leq 12$, are available online Wiringa .

IV.3 Charge form factors and Coulomb sum rules

The double differential cross section of the inclusive electron-nucleus scattering process in which an electron of initial four-momentum $k_{e}=(\bm{k}_{e},E_{e})$ scatters off a nuclear target to a state of four-momentum $k_{e}^{\prime}=(\bm{k}_{e}^{\prime},E_{e}^{\prime})$, the hadronic final state being undetected, can be written in the one-photon-exchange approximation as

[TABLE]

where

[TABLE]

and

[TABLE]

is the Mott cross section. In the above expressions, $\alpha\simeq 1/137$ is the fine structure constant, $d\Omega_{e^{\prime}}$ is the differential solid angle in the direction of $\bm{k}_{e^{\prime}}$, $q=k_{e}-k_{e^{\prime}}=(\bm{q},\omega)$ is the four-momentum transfer, and $Q^{2}=-q^{2}=\bm{q}^{2}-\omega^{2}$. The longitudinal and transverse response functions are defined as

[TABLE]

where $|0\rangle$ and $|f\rangle$ represent the nuclear initial and final states of energies $E_{0}$ and $E_{f}$, and $j_{L}({\bm{q}},\omega)$ and $j_{T}({\bm{q}},\omega)$ are the electromagnetic charge and current operators, respectively.

Recently, the quasielastic electromagnetic response functions of 4He and 12C have been computed within GFMC using realistic nuclear two- and three-body forces and consistent one- and two-body electroweak currents Lovato et al. (2015, 2016). Besides the transverse enhancement brought about by two-body current contributions, the authors of Ref. Lovato et al. (2016) have found no evidence of in-medium modification of the nucleon form factor in the analysis of the longitudinal response function of 12C. This is at variance with the findings of Ref. Cloët et al. (2016), where changes to the proton Dirac form factor induced by the nuclear medium leads to a dramatic quenching of the Coulomb sum rule,

[TABLE]

where $\omega_{\rm th}$ is the energy transfer corresponding to the inelastic threshold, and $G_{E}^{p}(Q^{2})$ is the proton electric form factor evaluated at four-momentum transfer $Q^{2}$.

The one-body charge operator employed in the GFMC calculations has the standard expressions obtained from a relativistic reduction of the time component of the covariant single-nucleon current,

[TABLE]

with

[TABLE]

In this work we adopted Kelly’s parametrization Kelly (2004) for the nucleon electric and magnetic form factors $G^{p,n}_{E,M}$.

In $R_{L}(\bm{q},\omega)$, the $\omega$ dependence enters via the energy-conserving $\delta$ function and the four-momentum transfer $Q^{2}$ of the electroweak form factors of the nucleon. The latter can be removed by evaluating these form factors at $Q^{2}_{\rm qe}=\bm{q}^{2}-\omega^{2}_{\rm qe}$, where $\omega_{\rm qe}$ is the energy transfer corresponding to the quasielastic peak, and by dividing the response by the factor $G_{E}^{p\,2}(Q_{\rm qe}^{2})$. Therefore, the Coulomb sum rule can be very well approximated by the following ground-state expectation value

[TABLE]

where $\mathcal{O}_{L}({\bm{q}})=j_{L}({\bm{q}},\omega_{\rm qe})/G_{E}^{p}(Q_{\rm qe}^{2})$ Lovato et al. (2013). The Coulomb sum rule defined in Eq. 68 only includes the inelastic contribution to $R_{\alpha}(\bm{q},\omega)$, i.e., the elastic contribution represented by the second term on the right-hand side of Eq. 71, where $|0;{\bm{q}}\rangle$ denotes the ground state of the nucleus recoiling with total momentum ${\bm{q}}$, has been removed. This term is proportional to the longitudinal elastic form factor, which is given by

[TABLE]

where $Q^{2}_{\rm el}=\bm{q}^{2}-\omega_{\rm el}^{2}$, and $\omega_{\rm el}$ is the energy transfer corresponding to elastic scattering, $\omega_{\rm el}=\sqrt{q^{2}+m_{A}^{2}}-m_{A}$ ($m_{A}$ is the mass of the target nucleus).

Neglecting the small spin-orbit contribution of Eq. 69, the Coulomb sum rule and the elastic form factor can be expressed as

[TABLE]

where $\tilde{\rho}_{N}(q)$ and $\tilde{\rho}_{\rm{NN}}(q)$ are the Fourier transform of the densities defined in Eqs. 57 and 58).

Here we compute the Coulomb sum rules and the elastic form factors of 4He, 16O, and 40Ca, to provide a useful benchmark for current and future analysis of electron-nucleus scattering data. In particular, our results for 16O and 40Ca, when compared to experiment, should further elucidate the role of in-medium modification of the nucleon form factors.

The 4He longitudinal elastic form factor is compared to experimental data in Fig. 16. Our theoretical results for the AV18 interaction significantly overpredict the diffraction minimum and maximum positions. Inclusion of the $3N$ force brings theory closer to experiment, but it is known that MECs are needed to further shift the peaks of the longitudinal elastic form factor to lower values of the momentum transfer and achieve agreement with experiment Marcucci et al. (2016); Carlson et al. (2015). This is shown in Fig. 17 where the GFMC longitudinal elastic form factor with and without MEC contributions is displayed. Note that up to $\simeq 6\,\rm fm^{-1}$ the CVMC form factor perfectly matches the GFMC result obtained without MECs.

The longitudinal form factor of 16O is shown in Fig. 18. The experimental data are well reproduced by our calculations once the $3N$ force is included. In analogy to 12C Lovato et al. (2013), it is plausible that two-body current contributions are negligible at low $q$, and become appreciable only for $q>3\,\rm fm^{-1}$. In fact, in the high-momentum region MECs interfere destructively with the one-body contributions, bringing theoretical prediction of 12C into closer agreement with experiment. This is consistent with the findings of Ref. Mihaila and Heisenberg (2000), where MECs improve the description of 16O experimental data above $2.5\,\rm fm^{-1}$.

As for the 40Ca nucleus, a better agreement with experiments is achieved when AV18 only is present in the Hamiltonian (see Fig. 19). Assuming that, as for 12C and 16O, two-body current contributions have little effect for $q\leq 3\,\rm fm^{-1}$, we can infer that the UIX potential moves the diffraction peaks to excessively low values of $q$. This failure of the UIX interaction is directly related to the behavior of the point proton density displayed in Fig. 3, where nucleons are pushed too far away from the center of mass when UIX is employed.

The longitudinal sum rules of 4He, 16O, and 40Ca for AV18+UIX obtained from Section IV.3 are displayed in Fig. 20. The best GFMC estimates for $S_{L}(\bm{q})$ in 4He and 12C Lovato et al. (2013) are also shown for comparison (solid symbols). GFMC calculations have been carried out employing the AV18+IL7 potential and considering the full contribution of one- and two-body electromagnetic currents. The latter have only a relatively small effect on the longitudinal sum rule, mainly affecting the magnitude of the peak for 12C and the region above $3\,\rm fm^{-1}$. In this region, in addition to MECs, the discrepancies between CVMC and GFMC are due to the spin-orbit contribution in the charge operator, neglected in CVMC calculations but included in Ref. Lovato et al. (2013). In the large $q$ limit, the CVMC sum rules differ from unity because of relativistic corrections in the charge current, which gives the factor $1/(1+Q_{\rm qe}^{2}/(4m^{2}))$ of Section IV.3.

Extracting the Coulomb sum rules from the experimental response functions involves nontrivial difficulties. The experimental determination of $S_{L}(\bm{q})$ requires measuring the associated $R_{L}(\bm{q},\omega)$ from the inelastic threshold to infinity. However, inclusive electron scattering experiments can only explore the space-like region of the four-momentum transfer $\omega<q$. Therefore, a meaningful comparison between theory and experiment requires estimating the strength outside the region covered by electron-scattering experiments. Furthermore, the authors of Ref. Lovato et al. (2016) have shown that the transitions to the low-lying states of 12C give significant contributions to $\langle 0|\mathcal{O}^{\dagger}_{L}(\bm{q})\mathcal{O}_{L}(\bm{q})|0\rangle$ that are not present in the longitudinal response functions extracted from inclusive $(e,e^{\prime})$ cross sections. Therefore, before comparing experiment with the present theory, which computes the sum rule of the total inelastic response rather than just the quasielastic one, these contributions have to be explicitly removed from the theoretical sum rule. In the 12C case, the transition form factors to $J^{\pi}=2^{+}$, $0_{2}^{+}$ (Hoyle), and $4^{+}$ states were taken from experiments. However, this approach is not suitable to the present work because of the large numbers of low-lying transitions of 16O and 40Ca. For this reason we refrain from reporting experimental data in Fig. 20.

V Conclusions

A variational Monte Carlo analysis of the properties of three closed-shell nuclei, 4He, 16O, and 40Ca, has been performed. We employed the accurate phenomenological nuclear Hamiltonian AV18+UIX, which is capable of simultaneously describing two-nucleon bound and scattering states, the binding energy of 4He, and the saturation density of isospin-symmetric nuclear matter. The CVMC algorithm has been improved by including five-body terms in the cluster expansion of all the spin-isospin dependent correlations. Therefore, this work represents significant progress with respect to Ref. Pieper et al. (1990, 1992), in which the older AV14+UVII Hamiltonian was employed, the cluster expansion was limited to four-body terms only, spin-orbit correlations were treated only at two-body cluster level, and other approximations were made in the construction of the wave function and in estimating the variational expectation values.

In order to perform extensive searches for the optimal variational parameters in the multidimensional parameter space defined by the employed wave functions, we implemented in the CVMC program the open-source library for nonlinear optimization NLopt Johnson . The accuracy of the optimized wave function has been tested against standard VMC and GFMC calculations for 4He using both AV18 and AV18+UIX, and against AFDMC results for 4He and 16O employing the AV6′ potential.

We present results for the binding energy, charge radius, one- and two-body densities, single-nucleon momentum distribution, charge form factor, and Coulomb sum rule, fully accounting for the high-momentum components of the nuclear interaction. We find that the UIX three-body potential, known to be attractive for $A\leq 12$, becomes repulsive for $A\geq 16$. At variance with the 4He case, the addition of the UIX potential makes 16O and 40Ca less bound. This repulsive effect is not limited to the binding energies. In 16O and 40Ca nucleons are pushed far away from the center of mass when the $3N$ force is included, resulting in larger radii, broader densities, and a shift of the charge form factor diffraction peaks towards smaller momenta. Although relying on different interaction schemes, a similar behavior of three-body interactions is found in CC and IM-SRG calculations for medium-heavy nuclei (see Hagen et al. (2014); Hebeler et al. (2015); Hergert et al. (2016) and references therein). We note that within CVMC there is no need to soften the $N\!N$ potential and to employ either the normal ordering procedure or a two-body density dependent approximation for the three-body force.

Although the UIX three-nucleon interaction manifests a change in behavior—from attractive to repulsive—for $A\geq 16$, it appears to provide a better description of radii, densities and charge form factors, of nuclei at least up to $A=16$. For instance, the charge radius and the position of the first peak in the longitudinal elastic form factor of 16O are better reproduced by the full AV18+UIX interaction than by the AV18 potential alone. This is no longer true in 40Ca, where the inclusion of the $3N$ potential yields a too large charge radius and shifts the diffraction peaks of the charge form factor towards too small momenta. The experimental data for 40Ca lie in between the CVMC theoretical predictions for AV18 and AV18+UIX. The fact that the AV18+UIX Hamiltonian is not adequate to describe medium-mass nuclei is consistent with the deficiencies in the theoretical prediction of isospin-symmetric nuclear matter employing the same interaction. Although the correct saturation density is obtained, the binding energy per nucleon is too small Akmal et al. (1998). In this regard, as a follow up of this work, we will consider local $N\!N$ potentials recently derived in coordinate space within chiral perturbation theory Gezerlis et al. (2013, 2014); Piarulli et al. (2015, 2016); Tews et al. (2016); Lynn et al. (2016); Logoteta et al. (2016). The latter are characterized by a spin-isospin structure analogous to the one of AV18+UIX so the CVMC can be straightforwardly extended to this class of interactions. It will be interesting to see whether local chiral effective field theory Hamiltonians provide a satisfactory description of 16O, 40Ca, and light nuclei.

We also computed the single-nucleon momentum distributions of 16O and 40Ca. These extend the VMC collection of Refs. Wiringa et al. (2014); Wiringa obtained using realistic phenomenological Hamiltonians, which include both $N\!N$ and $3N$ interactions. Together with the inclusion of three-body and higher-order terms in the cluster expansion, this makes the calculations of $n(\bm{k})$ accurate in both the high- and low-momentum regions. The universality of the tail of the momentum distribution, i.e., the independence of the high-momentum component upon the specific nucleus, has been confirmed for the selected interaction scheme. The momentum distributions are of immediate use for the studies of the high-momentum structure of nuclei, which includes the EMC effect and the analysis of short-range correlations in nuclei Arrington (2016). For the latter, the analysis of two-nucleon momentum distributions derived employing realistic two- and three-body nuclear interactions will be of great interest. A future project will focus on the CVMC computation of two-nucleon momentum distributions in medium-heavy nuclei, extending the VMC collection of Refs. Wiringa et al. (2014); R. B. Wiringa and providing a comparison with the findings of Ref. Alvioli et al. (2016).

We plan to employ the momentum distributions, and the average separation energies, computed in this work to evaluate the electroweak response functions of 16O and 40Ca in the impulse approximation, with particular emphasis on the role of three-nucleon forces, extending the study of Ref. Bacca et al. (2009) to heavier nuclei. This will be relevant for neutrino-oscillation experiments, such as the Deep Underground Neutrino Experiment (DUNE) Deep Underground Neutrino Experiment , and to elucidate quark and gluon effects in nuclei, which have long been actively sought, but never unambiguously identified.

We computed the Coulomb sum rules for closed-shell nuclei ranging from $A=4$ to $A=40$. Our calculations show very little $A$ dependence of the sum rules for $A\geq 12$ for momentum transfers as low as $1\,\rm fm^{-1}$. These results are also consistent with the recent GFMC calculation for 12C Lovato et al. (2016).

Another future project will be to examine closed-shell nuclei $+/-$ one nucleon, e.g., 15N, 15O, 17O, 17F, to study various properties such as spin-orbit splitting, which was previously evaluated in 15N using CVMC in Ref. Pieper and Pandharipande (1993), charge-symmetry breaking Wiringa et al. (2013), and $\beta$ decay.

Acknowledgements.

We thank I. Sick for providing us with the compilation of the experimental longitudinal elastic form factors and for the interpretation of the data. We are also thankful to J. Carlson and S. Gandolfi for insightful discussions. This work was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under the FRIB Theory Alliance Grant Contract No. DE-SC0013617 titled “FRIB Theory Center—A path for the science at FRIB” (D.L.), under the NUCLEI SciDAC grant (D.L. and A.L.), and under Contract No. DE-AC02-06CH11357 (A.L., S.C.P, and R.B.W). Computing time was provided by Los Alamos Open Supercomputing via the Institutional Computing (IC) program, by the National Energy Research Scientific Computing Center (NERSC), which is supported by the U.S. Department of Energy, Office of Science, under Contract No. DE-AC02-05CH11231, and by the Laboratory Computing Resource Center (LCRC) at Argonne National Laboratory.

Appendix: Wave function details

Figures 21, 22, 23, 24 and 25 and Table 11 provide the radial correlation functions and all the variational parameters for the systems under study for both AV18 and AV18+UIX.