Ab initio calculations of the isotopic dependence of nuclear clustering

TL;DR

This paper uses lattice Monte Carlo calculations based on chiral effective field theory to explore nuclear clustering, revealing the shape and correlations of alpha clusters in various isotopes and introducing a new computational method called the pinhole algorithm.

Contribution

It presents a novel computational approach, the pinhole algorithm, and applies it to analyze clustering and nucleon correlations in multiple isotopes, advancing understanding of nuclear structure.

Findings

Alpha clusters' shape and entanglement are characterized.

Density distributions and cluster geometries are determined for isotopes.

Structural similarities suggest analogous excitations in carbon isotopes.

Abstract

Nuclear clustering describes the appearance of structures resembling smaller nuclei such as alpha particles (4He nuclei) within the interior of a larger nucleus. While clustering is important for several well-known examples, much remains to be discovered about the general nature of clustering in nuclei. In this letter we present lattice Monte Carlo calculations based on chiral effective field theory for the ground states of helium, beryllium, carbon, and oxygen isotopes. By computing model-independent measures that probe three- and four-nucleon correlations at short distances, we determine the shape of the alpha clusters and the entanglement of nucleons comprising each alpha cluster with the outside medium. We also introduce a new computational approach called the pinhole algorithm, which solves a long-standing deficiency of auxiliary-field Monte Carlo simulations in computing density correlations relative to the center of mass. We use the pinhole algorithm to determine the proton and neutron density distributions and the geometry of cluster correlations in 12C, 14C, and 16C. The structural similarities among the carbon isotopes suggest that 14C and 16C have excitations analogous to the well-known Hoyle state resonance in 12C.

Figures

Tables

Table 1. Table 1: Observed charge radii from electron scattering and proton and neutron radii at leading order.

nucleus	observed charge radius	proton radius (LO)	neutron radius (LO)
12C	2.472(16) fm Schaller et al. (1982), 2.481(6) fm Sick (1982)	2.40(7) fm	2.39(5) fm
14C	2.497(17) fm Schaller et al. (1982)	2.43(7) fm	2.56(7) fm
16C	—	2.46(10) fm	2.65(8) fm

Equations

$$\exp \left[-\frac{{\sum }_i{{\mathbf{r}}_{\mathbf{i}}}^2}{2(2.6\mathrm{fm}{)}^2}\right]\prod _{j>k}\theta (|{\mathbf{r}}_{\mathbf{j}}-{\mathbf{r}}_{\mathbf{k}}|-1.7\mathrm{fm}).$$

$$\prod _{j>k}\exp \left[-\frac{{({\mathbf{r}}_{\mathbf{j}}-{\mathbf{r}}_{\mathbf{k}})}^2}{2(4.5\mathrm{fm}{)}^2}\right]\theta (|{\mathbf{r}}_{\mathbf{j}}-{\mathbf{r}}_{\mathbf{k}}|-1.7\mathrm{fm}).$$

$$a_{\mathrm{NL}}(\mathbf{n})$$

$$a_{\mathrm{NL}}^{†}(\mathbf{n})$$

$$\rho (\mathbf{n})$$

$${\rho }_S(\mathbf{n})$$

$${\rho }_I(\mathbf{n})$$

$${\rho }_{S,I}(\mathbf{n})$$

$${\rho }_{\mathrm{NL}}(\mathbf{n})$$

$${\rho }_{S,\mathrm{NL}}(\mathbf{n})$$

$${\rho }_{I,\mathrm{NL}}(\mathbf{n})$$

$${\rho }_{S,I,\mathrm{NL}}(\mathbf{n})$$

$$V_0=\frac{c_0}{2}\sum _{{\mathbf{n}}^{\prime },\mathbf{n},{\mathbf{n}}^{\prime \prime }}:{\rho }_{\mathrm{NL}}({\mathbf{n}}^{\prime })f_{s_{\mathrm{L}}}({\mathbf{n}}^{\prime }-\mathbf{n})f_{s_{\mathrm{L}}}(\mathbf{n}-{\mathbf{n}}^{\prime \prime }){\rho }_{\mathrm{NL}}({\mathbf{n}}^{\prime \prime }):$$

$$f_{s_{\mathrm{L}}}(\mathbf{n})$$

$$=s_L\mathrm{for}|\mathbf{n}|=1,$$

$$=0\mathrm{otherwise}.$$

$$V_{\mathrm{OPE}}=-\frac{g_A^2}{8f_{\pi }^2}\sum _{{\mathbf{n}}^{\prime },\mathbf{n},S^{\prime },S,I}:{\rho }_{S^{\prime },I}({\mathbf{n}}^{\prime })f_{S^{\prime }S}({\mathbf{n}}^{\prime }-\mathbf{n}){\rho }_{S,I}(\mathbf{n}):,$$

$$f_{S^{\prime }S}({\mathbf{n}}^{\prime }-\mathbf{n})=\frac{1}{L^3}\sum _{\mathbf{q}}\frac{\exp [-i\mathbf{q}\cdot ({\mathbf{n}}^{\prime }-\mathbf{n})-b_{\pi }{\mathbf{q}}^2]q_{S^{\prime }}{q}_S}{{\mathbf{q}}^2+m_{\pi }^2},$$

$$H_{\mathrm{free}}=$$

$$+\frac{3}{40m}\sum _{\mathbf{n},i}\sum _{⟨⟨{\mathbf{n}}^{\prime }\mathbf{n}⟩{⟩}_i}{a}^{†}({\mathbf{n}}^{\prime })a(\mathbf{n})-\frac{1}{180m}\sum _{\mathbf{n},i}\sum _{⟨⟨⟨{\mathbf{n}}^{\prime }\mathbf{n}⟩⟩{⟩}_i}{a}^{†}({\mathbf{n}}^{\prime })a(\mathbf{n}).$$

$$H_{\mathrm{B}}=H_{\mathrm{free}}+V_0+V_{\mathrm{OPE}},$$

$$M=:\exp [-H{\alpha }_t]:.$$

$$Z_{f,i}(L_t)=⟨{\Psi }_f|M^{L_t}|{\Psi }_i⟩.$$

$$Z_{f,i}(L_t)=⟨{\Psi }_f|M_{\ast }^{L_t^{\prime }}{M}^{L_t}{M}_{\ast }^{L_t^{\prime }}|{\Psi }_i⟩.$$

$$:\exp \left(-\frac{c{\alpha }_t}{2}{\rho }^2\right):=\sqrt{\frac{1}{2\pi }}{\int }_{-\infty }^{\infty }ds:\exp \left(-\frac{1}{2}{s}^2+\sqrt{-c{\alpha }_t}s\rho \right):.$$

$$V_s^{(n_t)}=\sqrt{-c_0}\sum _{\mathbf{n},{\mathbf{n}}^{\prime }}{\rho }_{\mathrm{NL}}(\mathbf{n})f_{s_{\mathrm{L}}}(\mathbf{n}-{\mathbf{n}}^{\prime })s({\mathbf{n}}^{\prime },n_t),$$

$$V_{ss}^{(n_t)}=\frac{1}{2}\sum _{\mathbf{n}}{s}^2(\mathbf{n},n_t).$$

$$V_{\pi }^{(n_t)}=\frac{g_A}{2f_{\pi }}\sum _{\mathbf{n},{\mathbf{n}}^{\prime },S,I}{\rho }_{S,I}({\mathbf{n}}^{\prime })f_S^{\pi }({\mathbf{n}}^{\prime }-\mathbf{n}){\pi }_I(\mathbf{n},n_t),$$

$$V_{\pi \pi }^{(n_t)}=\frac{1}{2}\sum _{\mathbf{n},{\mathbf{n}}^{\prime },I}{\pi }_I({\mathbf{n}}^{\prime },n_t)f^{\pi \pi }({\mathbf{n}}^{\prime }-\mathbf{n}){\pi }_I(\mathbf{n},n_t),$$

$$f_S^{\pi }({\mathbf{n}}^{\prime }-\mathbf{n})=\frac{1}{L^3}\sum _{\mathbf{q}}\exp [-i\mathbf{q}\cdot ({\mathbf{n}}^{\prime }-\mathbf{n})]q_S,$$

$$f^{\pi \pi }({\mathbf{n}}^{\prime }-\mathbf{n})=\frac{1}{L^3}\sum _{\mathbf{q}}\exp [-i\mathbf{q}\cdot ({\mathbf{n}}^{\prime }-\mathbf{n})+b_{\pi }{\mathbf{q}}^2]({\mathbf{q}}^2+m_{\pi }^2).$$

$$M=\int D{s}^{(n_t)}D{\pi }^{(n_t)}{M}^{(n_t)},$$

$$M^{(n_t)}=:\exp \left(-H_{\mathrm{free}}{\alpha }_t-V_s^{(n_t)}\sqrt{{\alpha }_t}-V_{ss}^{(n_t)}-V_{\pi }^{(n_t)}{\alpha }_t-V_{\pi \pi }^{(n_t)}{\alpha }_t\right):.$$

$$E(t)=E_0+c\exp [-\Delta Et],$$

$$:\exp \left[\sum _{\mathbf{n}}ϵ(\mathbf{n})\rho (\mathbf{n})\right]:$$

Full text

There have been many exciting recent advances in ab initio nuclear structure theory Romero-Redondo et al. (2014); Maris et al. (2014); Dytrych et al. (2016); Duguet et al. (2016); Stroberg et al. (2017); Garcia Ruiz et al. (2016); Hagen et al. (2016); Elhatisari et al. (2016) which link nuclear forces to nuclear structure in impressive agreement with experimental data. However, we still know very little about the quantum correlations among nucleons that give rise to nuclear clustering and collective behavior. The main difficulty in studying alpha clusters in nuclei is that the calculation must include four-nucleon correlations. Unfortunately in many cases this dramatically increases the amount of computer memory and computing time needed in calculations of heavier nuclei. Nevertheless there is promising work in progress using the symmetry-adapted no-core shell model Launey et al. (2016), antisymmetrized molecular dynamics Yoshida and Kanada-En’yo (2016), fermionic molecular dynamics Feldmeier and Neff (2016), the alpha-container model Schuck et al. (2016), Monte Carlo shell model Yoshida et al. (2014), and Green’s function Monte Carlo Lovato et al. (2016).

Lattice calculations using chiral effective field theory and auxiliary-field Monte Carlo methods have probed alpha clustering in the 12C and 16O systems Epelbaum et al. (2011, 2012, 2013, 2014). However these lattice simulations have encountered severe Monte Carlo sign oscillations in cases where the number of protons $Z$ and number of neutrons $N$ are different. In this work we solve this problem by using a new leading-order lattice action that retains a greater amount of symmetry, thereby removing nearly all of the Monte Carlo sign oscillations. The relevant symmetry is Wigner’s SU(4) spin-isospin symmetry Wigner (1937), where the four nucleon degrees of freedom can be rotated as four components of a complex vector. Previous attempts using SU(4) symmetry had failed due to the tendency of nuclei to overbind in larger nuclei. However recent progress has uncovered important connections between local interactions and nuclear binding, as well as the significance of the alpha-alpha interaction Elhatisari et al. (2015, 2016); Rokash et al. (2016). Following this approach, we have constructed a leading-order lattice action with highly-suppressed sign oscillations and which reproduces the ground-state binding energies of the hydrogen, helium, beryllium, carbon, and oxygen isotopes to an accuracy of 0.7 MeV per nucleon or better. The lattice results are shown in panel a of Fig. 1 in comparison with the observed ground state energies. The astonishingly good agreement at leading order in chiral effective field theory with only three free parameters is quite remarkable and bodes well for future calculations at higher orders. We use auxiliary-field Monte Carlo simulations with a spatial lattice spacing of 1.97 fm and lattice time spacing 1.97 fm$/c$. We comment that the results for these ground state energies are equally good when including Coulomb repulsion and a slightly more attractive nucleon-nucleon short-range interaction. The full details of the lattice interaction, nucleon-nucleon phase shifts, simulation methods, and results are given in the Supplemental Materials.

Let $\rho(\bf n)$ be the total nucleon density operator on lattice site $\bf n$. We will use short-distance three- and four-nucleon operators as probes of the nuclear clusters. To construct a probe for alpha clusters, we define $\rho_{4}$ as the expectation value of $:\rho^{4}({\bf n})/4!:$ summed over $\bf n$. The $::$ symbols indicate normal-ordering where all annihilation operators are moved to the right and all creation operators are moved to the left. For nuclei with even $Z$ and even $N$, there are likely no well-defined 3H or 3He clusters since their formation is not energetically favorable. Therefore we can use short-distance three-nucleon operators as a second probe of alpha clusters. We define $\rho_{3}$ as the expectation value of $:\rho^{3}({\bf n})/3!:$ summed over $\bf n$. A 3H or 3He cluster may form in nuclei with odd $Z$ or odd $N$. In these cases we can use spin- and isospin-dependent three-nucleon operators to probe the 3H and 3He clusters. As we consider only nuclei with even $Z$ and even $N$ here, we focus on $\rho_{3}$ and $\rho_{4}$ for the remainder of the discussion. We note that another measure of clustering in nuclei by measuring short-distance correlations has been introduced in Ref. Zhang et al. (2016).

Due to divergences at short distances, $\rho_{3}$ and $\rho_{4}$ will depend on the short-distance regularization scale, which in our case is the lattice spacing. However the regularization-scale dependence of $\rho_{3}$ and $\rho_{4}$ does not depend on the nucleus being considered. Therefore if we let $\rho_{3,\alpha}$ and $\rho_{4,\alpha}$ be the corresponding values for the alpha particle, then the ratios $\rho_{3}/\rho_{3,\alpha}$ and $\rho_{4}/\rho_{4,\alpha}$ are free from short-distance divergences and are model-independent quantities up to contributions from higher-dimensional operators in an operator product expansion. The derivations of these statements are given in the supplemental materials. We have computed $\rho_{3}$ and $\rho_{4}$ for the helium, beryllium, carbon, and oxygen isotopes. As our leading-order interactions are invariant under an isospin mirror flip that interchanges protons and neutrons, we focus here on neutron-rich nuclei. The results for $\rho_{3}/\rho_{3,\alpha}$ and $\rho_{4}/\rho_{4,\alpha}$ are presented in panel b of Fig. 1. As we might expect, the values for $\rho_{3}/\rho_{3,\alpha}$ and $\rho_{4}/\rho_{4,\alpha}$ are roughly the same for the different neutron-rich isotopes of each element.

Since $\rho_{4}$ involves four nucleons, it couples to the center of the alpha cluster while $\rho_{3}$ gets a contribution from a wider portion of the alpha-cluster wave function. Therefore, a value larger than 1 for the ratio of $\rho_{4}/\rho_{4,\alpha}$ to $\rho_{3}/\rho_{3,\alpha}$ corresponds to a more compact alpha-cluster shape than in vacuum, and a value less than 1 corresponds to a more diffuse alpha-cluster shape. In panel b of Fig. 1 we observe that the ratio of $\rho_{4}/\rho_{4,\alpha}$ to $\rho_{3}/\rho_{3,\alpha}$ starts at 1 or slightly above 1 when $N$ is comparable to $Z$, and the ratio gradually decreases as the number of neutrons is increased. This is evidence for the swelling of the alpha clusters as the system becomes saturated with excess neutrons. The effect has also been seen in 6He and 8He in Green’s Function Monte Carlo calculations Wiringa et al. (2000).

We comment here that if one wants to study the swelling of alpha clusters in detail, then there are other local operators that provide more direct geometrical information such as the operators $:\!\rho^{3}({\bf n})\rho({\bf n}^{\prime})\!:$ and $:\!\rho^{2}({\bf n})\rho^{2}({\bf n}^{\prime})\!:$, where ${\bf n}^{\prime}$ is a nearest-neighbor site to ${\bf n}$. These local operators have the advantage of measuring four-nucleon correlations directly rather than inferring them from the ratio of four-body and three-body correlations, which may not work well for cases with very large isospin imbalance.

The traditional approach to nuclear clustering usually involves a variational ansatz where the wave function is expanded in terms of some chosen set of alpha-cluster wave functions. However the answer obtained this way may depend strongly on the details of the interactions and the choice of alpha-cluster wave functions. This problem of model dependence is solved by calculating short-range multi-nucleon quantities. Even though we use only short-range operators, the quantities $\rho_{3}/\rho_{3,\alpha}$ and $\rho_{4}/\rho_{4,\alpha}$ act as high-fidelity alpha-cluster detectors. Their values are strongly enhanced if the nuclear wave function has a well-defined alpha-cluster substructure. As shown in the supplemental materials, the enhancement factor for $\rho_{3}/\rho_{3,\alpha}$ is $(R_{A}/R_{\alpha})^{6}$, where $R_{A}$ is the nuclear radius and $R_{\alpha}$ is the alpha-particle radius. The enhancement factor for $\rho_{4}/\rho_{4,\alpha}$ is an even larger factor of $(R_{A}/R_{\alpha})^{9}$.

We denote the number of alpha clusters as $N_{\alpha}$. A simple counting of protons gives $N_{\alpha}=1$ for neutron-rich helium, $N_{\alpha}=2$ for neutron-rich beryllium, $N_{\alpha}=3$ for neutron-rich carbon, and $N_{\alpha}=4$ for neutron-rich oxygen. However the alpha clusters are immersed in a complex many-body system, and it is useful to quantify the entanglement of the nucleons comprising each alpha cluster with the outside medium. The observables $\rho_{3}/\rho_{3,\alpha}$ and $\rho_{4}/\rho_{4,\alpha}$ are useful for this purpose. Let us define $\delta^{\rho_{3}}_{\alpha}$ as the difference $\rho_{3}/\rho_{3,\alpha}-N_{\alpha}$ divided by $N_{\alpha}$. Since $\delta^{\rho_{3}}_{\alpha}$ measures the deviation of the nuclear wave function from a pure product state of alpha clusters and excess nucleons, we call it the $\rho_{3}$-entanglement of the alpha clusters. In an analogous manner, we can also define the $\rho_{4}$-entanglement $\delta^{\rho_{4}}_{\alpha}$ as the difference $\rho_{4}/\rho_{4,\alpha}-N_{\alpha}$ divided by $N_{\alpha}$. $\delta^{\rho_{4}}_{\alpha}$ turns out to be quantitatively similar to $\delta^{\rho_{3}}_{\alpha}$, though with more sensitivity to the shape of the alpha clusters.

In panel b of Fig. 1, we show $N_{\alpha}$ along with the ratios $\rho_{3}/\rho_{3,\alpha}$ and $\rho_{4}/\rho_{4,\alpha}$. The relative excess of $\rho_{3}/\rho_{3,\alpha}$ compared to $N_{\alpha}$ gives $\delta^{\rho_{3}}_{\alpha}$, and the relative excess of $\rho_{4}/\rho_{4,\alpha}$ compared to $N_{\alpha}$ gives $\delta^{\rho_{4}}_{\alpha}$. We see that $\delta^{\rho_{3}}_{\alpha}$ is negligible for 6He and 8He, indicating an almost pure product state of alpha clusters and excess neutrons. For the beryllium isotopes, $\delta^{\rho_{3}}_{\alpha}$ is about $0.18$111In this leading-order calculation the 8Be ground state is about 1 MeV below the two-$\alpha$ threshold. The addition of the Coulomb interaction and other corrections should push this energy closer to threshold, and one expects $\delta^{\rho_{3}}_{\alpha}$ to decrease as a result. for 8Be and rises to about $0.34$ for 14Be. For the carbon isotopes, it is about $0.28$ for 12C and rises to a maximum of about $0.50$ near the drip line. For the oxygen isotopes, $\delta^{\rho_{3}}_{\alpha}$ is about $0.50$ for 16O and increases with neutron number up to $0.73$. For such high values of the $\rho_{3}$-entanglement, we expect a simple picture in terms of alpha clusters and excess neutrons will break down. $\delta^{\rho_{3}}_{\alpha}$ should be much lower for excited cluster-like states of the oxygen isotopes. With $\rho_{3}$-entanglement, we have a model-independent quantitative measure of nuclear cluster formation in terms of entanglement of the wave function. Our results show that the transition from cluster-like states in light systems to nuclear liquid-like states in heavier systems should not be viewed as a simple suppression of multi-nucleon short-distance correlations, but rather an increasing entanglement of the nucleons involved in the multi-nucleon correlations.

Despite the many computational advantages of auxiliary-field Monte Carlo methods, one fundamental deficiency is that the simulations involve quantum states that are superpositions of many different center-of-mass positions. Therefore density distributions of the nucleons cannot be computed directly. To solve this problem we have developed a new method called the pinhole algorithm. In this algorithm an opaque screen is placed at the middle time step with pinholes bearing spin and isospin labels that allow nucleons with the corresponding spin and isospin to pass. We use $A$ pinholes for a simulation of $A$ nucleons, and the locations as well as the spin and isospin labels of the pinholes are updated by Monte Carlo importance sampling. From the simulations, we obtain the expectation value of the normal-ordered $A$-body density operator $:\rho_{i_{1},j_{1}}({\bf n}_{1})\cdots\rho_{i_{A},j_{A}}({\bf n}_{A}):$, where $\rho_{i,j}$ is the density operator for a nucleon with spin $i$ and isospin $j$.

Using the pinhole algorithm, we have computed the proton and neutron densities for the ground states of 12C, 14C, and 16C. In order to account for the nonzero size of the nucleons, we have convolved the point-nucleon distributions with a Gaussian distribution with root-mean-square radius 0.84 fm, the charge radius of the proton Belushkin et al. (2007); Pohl et al. (2010). The results are shown in Fig. 2 along with the experimentally observed proton densities for 12C and 14C Kline et al. (1973), which we define as the charge density divided by the electric charge $e$. From Fig. 2 we see that the agreement between the calculated proton densities and experimental data for 12C and 14C is rather good. We show data for $L_{t}=7,9,11,13,15$ time steps. The fact that the results have little dependence on $L_{t}$ means that we are seeing ground state properties. As we increase the number of neutrons and go from 12C to 16C, the shape of the proton density profile remains roughly the same. However there is a gradual decrease in the central density and a broadening of the proton density distribution. We see also that the excess neutrons in 14C and 16C are distributed fairly evenly, appearing in both the central region as well as the tail.

We now study the alpha-cluster structures of 12C, 14C, and 16C in more detail. In order to probe the alpha cluster geometry, we use the fact that there is only one spin-up proton per alpha cluster. Using the pinhole algorithm, we consider the triangular shapes formed by the three spin-up protons in the carbon isotopes. This correlation function is free of short-distance divergences, and so, up to the contribution of higher-dimensional operators, it provides a model-independent measure that serves as a proxy for the geometry of the alpha-cluster configurations.

The three spin-up protons form the vertices of a triangle. When collecting the lattice simulation data, we rotate the triangle so that the longest side lies on the $x$-axis. We also rescale the triangle so the longest side has length one, and flip the triangle, if needed, so that the third spin-up proton is in the upper half of the $xy$-plane. Histograms of the third spin-up proton probability distributions for 12C, 14C, and 16C are plotted in panel a, b, c of Fig. 3 using the data at $L_{t}=15$ time steps. The data for other values of $L_{t}$ are almost identical. There is some jaggedness due to the discreteness of the lattice, but we see quite clearly that the histograms for 12C, 14C, and 16C are very similar. While there is some increase in the overall radius of the nucleus, the rescaled cluster geometry of the three carbon isotopes remain largely the same. In each case we see that there is a strong preference for triangles where the largest angle is less than or equal to 90 degrees. We should note that idea that the ground state of 12C has an acute triangular alpha-cluster structure has a long history dating back to Ref. Hafstad and Teller (1938).

Given the rich cluster structure of the excited states of 12C, this raises the interesting possibility of similar cluster states appearing in 14C and 16C. In particular, the bound $0^{+}_{2}$ state at 6.59 MeV above the ground state of 14C may be a bound-state analog to the Hoyle state resonance in 12C at 7.65 MeV. It may also have a clean experimental signature since low-lying neutron excitations are suppressed by the shell closure at eight neutrons. There is also a bound $0^{+}_{2}$ in 16C, however in this case one expects low-lying two-neutron excitations to be important, thereby making the analysis more complicated. We note that there is ample experimental evidence for the cluster properties of the neutron-rich beryllium and carbon isotopes Bohlen et al. (2007, 2003); Freer (2008); Marin-Lambarri et al. (2014).

In order to analyze what we are seeing in the lattice data, we can make a simple Gaussian lattice model of the distribution of the spin-up protons. We consider a probability distribution $P({\bf{r}_{1}},{\bf{r}_{2}},{\bf{r}_{3}})$ on our lattice grid for the positions of the protons ${\bf{r}_{1}}$, ${\bf{r}_{2}}$, and ${\bf{r}_{3}}$. We take the probability distribution to be a product of Gaussians with root-mean-square radius 2.6 fm (charge radius of 14C) and unit step functions which vanish if the magnitude of ${\bf{r}_{1}}-{\bf{r}_{2}}$, ${\bf{r}_{2}}-{\bf{r}_{3}}$, or ${\bf{r}_{3}}-{\bf{r}_{1}}$ is smaller than 1.7 fm (charge radius of 4He),

[TABLE]

We can factor out the center-of-mass distribution of the three spin-up protons and recast the Gaussian factors as a product of Gaussians for the separation vectors ${\bf{r}_{1}}-{\bf{r}_{2}}$, ${\bf{r}_{2}}-{\bf{r}_{3}}$, or ${\bf{r}_{3}}-{\bf{r}_{1}}$ with root-mean-square radius 4.5 fm,

[TABLE]

In panel d of Fig. 3 we show the third spin-up proton probability distribution corresponding to this model. Despite the simplicity of this model with no free parameters, we note the good agreement with the lattice data for 12C, 14C, and 16C. The only discrepancy is that the model overpredicts the probability of producing obtuse triangular configurations. This indicates that there are some additional correlations between the clusters that go beyond this simple Gaussian lattice model.

In this letter we have presented a number of novel approaches to computing and quantifying clustering and entanglement in nuclei. We hope that this work may help to accelerate progress in theoretical and experimental efforts to understand the correlations that produce nuclear clustering and collective behavior.

Acknowledgement

We are grateful for the hospitality of the Kavli Institute for Theoretical Physics at UC Santa Barbara for hosting E.E., H.K., and D.L. We are indebted to Ingo Sick for providing the experimental data tables on the electric form factor for 12C. We acknowledge partial financial support from the CRC110: Deutsche Forschungsgemeinschaft (SGB/TR 110, “Symmetries and the Emergence of Structure in QCD”), the BMBF (Verbundprojekt 05P2015 - NUSTAR R&D), the U.S. Department of Energy (DE-FG02-03ER41260), and U.S. National Science Foundation grant No. PHY-1307453. Further support was provided by the Magnus Ehrnrooth Foundation of the Finnish Society of Sciences and Letters and the Chinese Academy of Sciences (CAS) President’s International Fellowship Initiative (PIFI) grant no. 2017VMA0025. The computational resources were provided by the Jülich Supercomputing Centre at Forschungszentrum Jülich, RWTH Aachen, and North Carolina State University.

Supplemental Materials

Lattice interactions

In our lattice simulations the spatial lattice spacing is taken to be $a=1.97$ fm, and the time lattice spacing is $a_{t}=1.97$ fm$/c$. The axial-vector coupling constant is $g_{A}=1.29$, pion decay constant is$\ f_{\pi}=92.2$ MeV, pion mass is $m_{\pi}=m_{\pi^{0}}=134.98$ MeV, and nucleon mass is $m=938.92$ MeV. We write $\sigma_{S}$ with $S=1,2,3$ for the spin Pauli matrices, and $\tau_{I}$ with $I=1,2,3$ for the isospin Pauli matrices. We use dimensionless lattice units, where the physical quantities are multiplied by powers of the spatial lattice spacing $a$ to make dimensionless combinations. We write $\alpha_{t}$ for the ratio $a_{t}/a$.

The notation $\sum_{\langle{\bf n^{\prime}\,n}\rangle}$ represents the summation over nearest-neighbor lattice sites of n. We use $\sum_{\langle{\bf n^{\prime}\,n}\rangle_{i}}$ to indicate the sum over nearest-neighbor lattice sites of n along the $i^{\rm th}$ spatial axis. Similarly, $\sum_{\langle\langle{\bf n^{\prime}\,n}\rangle\rangle_{i}}$ is the sum over next-to-nearest-neighbor lattice sites of n along the $i^{\rm th}$ axis, and $\sum_{\langle\langle\langle{\bf n^{\prime}\,n}\rangle\rangle\rangle_{i}}$ is the sum over next-to-next-to-nearest-neighbor lattice sites of n along the $i^{\rm th}$ axis. Our lattice system is defined on an $L\times L\times L$ periodic cube, and so the summations over ${\bf n^{\prime}}$ are defined with periodic boundary conditions.

In our notation $a_{\rm NL}$ is a four-component spin-isospin column vector while $a^{\dagger}_{\rm NL}$ is a four-component spin-isospin row vector. For real parameter $s_{\rm NL}$, we define the nonlocal annihilation and creation operators for each spin and isospin component of the nucleon,

[TABLE]

For spin indices $S=1,2,3,$ and isospin indices $I=1,2,3$, we define point-like densities,

[TABLE]

and also the smeared nonlocal densities,

[TABLE]

For the leading-order short-range interactions we use

[TABLE]

where $f_{s_{\rm L}}$ is defined for real parameter $s_{L}$ as

[TABLE]

The :: symbol indicates normal ordering, where the annihilation operators are on the right-hand side and the creation operators are on the left-hand side.

The one-pion exchange interaction is given by

[TABLE]

where $f_{S^{\prime}S}$ is defined as

[TABLE]

and each $q_{S}$ is an integer multiplied by $2\pi/L$. The parameter $b_{\pi}$ removes short-distance lattice artifacts in the one-pion exchange interaction, and in this work we use the value $b_{\pi}=0.700$. We take the free lattice Hamiltonian to have the form Epelbaum et al. (2010)

[TABLE]

The full leading-order (LO) lattice Hamiltonian can be written as

[TABLE]

with $s_{\rm NL}=0.0800$, $s_{\rm L}=0.0800$, and $c_{0}=-0.1850$. In tuning our interactions here, we fit the parameters $s_{\rm NL}$, $s_{\rm L}$, and $c_{0}$ to the average inverse scattering length and effective range of the two $s$-wave channels, as well as the finite-volume energies of 8Be. The finite-volume energies for 8Be give a measure of the alpha-alpha scattering length, which was emphasized in Ref. Elhatisari et al. (2016) as a sensitive indicator correlated with the binding energies of medium-mass nuclei.

Nucleon-nucleon scattering

The details of the nucleon-nucleon scattering calculations are given in Ref. Elhatisari et al. (2016). In Fig. S1 we show the LO lattice phase shifts for proton-neutron scattering versus the center-of-mass relative momentum. For comparison we also present phase shifts from the Nijmegen partial wave analysis Stoks et al. (1993). In the first row, the data in panels a, b, c, d correspond to ${{}^{1}s_{0}},{{}^{3}s_{1}},{{}^{1}p_{1}},{{}^{3}p_{0}}$ respectively. In the second row, panels e, f, g, h correspond to ${{}^{3}p_{1}},{{}^{3}p_{2}},{{}^{1}d_{2}},{{}^{3}d_{1}}$ respectively. In the third row, panels i, j, k, l correspond to ${{}^{3}d_{2}},{{}^{3}d_{3}},{\varepsilon_{1}},{\varepsilon_{2}}$ respectively. As can been seen from Fig. S1, the ${}^{1}s_{0}$ phase shift requires significant higher-order corrections. These leading-order results are just the first step in the chiral effective field theory expansion, and the phase shifts would be systematically improved at each higher order, NLO, NNLO, etc. While the behavior of the ${{}^{1}s_{0}}$ phase shift near threshold seems rather poor, it requires only a rather small higher-order correction to reproduce the proper ${{}^{1}s_{0}}$ phase shift. We have checked this explicitly and it is also one of the central themes in a recent paper on nuclear physics expanded around the unitarity limit König et al. (2017). The key point is that the ${{}^{1}s_{0}}$ phase shift depends strongly on small changes in the ${{}^{1}s_{0}}$ coupling strength because it sits very close to the unitarity limit where the scattering length diverges.

Euclidean time projection and auxiliary-field Monte Carlo

The Euclidean time transfer matrix $M$ is defined as the normal-ordered exponential of the lattice Hamiltonian $H$ over one time lattice step,

[TABLE]

We use an initial state $|\Psi_{i}\rangle$ and final state $|\Psi_{f}\rangle$ that have nonzero overlap with the ground state nucleus of interest. By multiplying by powers of $M$ upon $|\Psi_{i}\rangle$, we can project out the ground state. We compute projection amplitudes of the form

[TABLE]

By calculating the ratio $Z_{f,i}(L_{t})/Z_{f,i}(L_{t}-1)$ for large $L_{t}$ we can determine the ground state energy.

It is useful however to first prepare the initial state using a simpler transfer matrix $M_{*}$ that is a good approximation to $M$. We choose $M_{*}$ to be invariant under Wigner’s SU$(4)$ symmetry Wigner (1937). The SU(4) symmetry eliminates sign oscillations from auxiliary-field Monte Carlo simulations of $M_{*}$ Chen et al. (2004); Lee (2007). $M_{*}$ has the same form as $M,$ but the operator coefficients that violate SU(4) symmetry are turned off. We use $M_{*}$ as an approximate low-energy filter by multiplying the initial and final states by $M_{*}$ some fixed number of times, $L_{t}^{\prime}$,

[TABLE]

We use auxiliary fields to generate the lattice interactions.  The auxiliary field method can be viewed as a Gaussian integral formula which relates the exponential of the two-particle density, $\rho^{2}$, to the integral of the exponential of the one-particle density, $\rho$,

[TABLE]

The normal ordering symbol :: ensures that the operator products of the creation and annihilation operators behave as classical anticommuting Grassmann variables Lee (2009). We use this integral identity to introduce auxiliary fields at every lattice site Hubbard (1959); Stratonovich (1958); Koonin (1986). The pion fields are treated in a manner similar to the auxiliary fields.

We couple the auxiliary field $s$ at time step $n_{t}$ to $\rho_{\rm NL}$ through a convolution with the smearing function $f_{s_{\rm L}}$. The linear term in the auxiliary field is

[TABLE]

and the quadratic term in the auxiliary field is

[TABLE]

For the one-pion exchange interaction, the gradient of the pion field $\pi_{I}$ is coupled to the point-like density $\rho_{S,I}$,

[TABLE]

[TABLE]

where $f^{\pi}_{S}$ ($S=1,2,3$) and $f^{\pi\pi}$ are defined as

[TABLE]

[TABLE]

Then the transfer matrix at leading order can be written as an path integral,

[TABLE]

where $Ds^{(n_{t})}$ is the path integral measure for $s$ at time step $n_{t}$, $D\pi^{(n_{t})}$ is the path integral measure for $\pi_{I}$ ($I=1,2,3$) at time step $n_{t}$, and

[TABLE]

In the projection Monte Carlo calculations we use the same procedure for the initial states as discussed in Ref. Elhatisari et al. (2016). Four nucleons are inserted at each time step. For neutron-rich nuclei we also insert pairs of spin-up and spin-down neutrons, and for proton-rich nuclei we insert pairs of spin-up and spin-down protons. For the calculations of 3H and 3He we use an $L\simeq 16\,{\rm fm}$ periodic box, and for the rest of the nuclei we use an $L\simeq 12\,{\rm fm}$ periodic box.

Results for the ground state energies

In Fig. S2 we show the energy versus projection time for 3H and 3He. The error bars indicate one standard deviation errors due to the stochastic noise of the Monte Carlo simulations. The lines are extrapolations to infinite projection time using the functional form

[TABLE]

where $E_{0}$ is the ground state energy that we wish to determine. The results for the helium isotopes are shown in Fig. S3, the beryllium isotopes in Fig. S4, the carbon isotopes in Fig. S5, and the oxygen isotopes in Fig. S6.

Results for $\rho_{3}$ and $\rho_{4}$

We compute $\rho_{3}$ by inserting the operator

[TABLE]

at the middle time step and taking three numerical derivatives with respect to $\epsilon({\bf n})$ for infinitesmally small $\epsilon({\bf n})$. We then divide by $3!$ and sum over n. For $\rho_{4}$ we compute four numerical derivatives with respect to $\epsilon({\bf n})$, divide by $4!$, and sum over n.

In Fig. S7 we show $\rho_{3}$ versus projection time for the neutron-rich helium, beryllium, and carbon isotopes. The error bars indicate one standard deviation errors due to the stochastic noise of the Monte Carlo simulations. The lines are extrapolations to infinite projection time using the functional forms

[TABLE]

where $\Delta E$ is determined from the ground state energy fit in Eq. (31). The factor of $t/2$ rather than $t$ comes from the fact that we are computing expectation values of $:\rho^{3}({\bf n})/3!:$ and $:\rho^{4}({\bf n})/4!:$ inserted at the middle time step. This leads to exponential corrections from matrix elements connecting the ground state to the first excited state, each of which are propagated for time duration $t/2$. In Fig. S8 we show $\rho_{4}$ versus projection time for the neutron-rich helium, beryllium, and carbon isotopes.

Local cluster operators and operator product expansion

Let us consider any short-distance three-nucleon operator of the form

[TABLE]

where $f$ is a spatially-localized function with width parameter $\Delta r$. This operator can be expanded as a sum of local operator products Wilson (1969); Zimmermann (1973),

[TABLE]

where $d_{n}$ is the momentum dimension of the operator $O_{n}(\vec{r})$, $\Lambda$ is the renormalization momentum scale, and $c_{n}$ gives the dependence on $\Lambda$ through quantum loop effects. The lowest possible value for $d_{n}$ is $9$ and is associated with the operator product $:\rho^{3}(\vec{r}):$. Other operators with the same quantum numbers as $X_{3}(\vec{r})$ have higher dimension.

Let us now consider the expectation value of $X_{3}(\vec{r})$ for the alpha particle and for an arbitrary nucleus which we label $A$. The ratio of these expectation values is then

[TABLE]

where the omitted terms are contributions from higher operators in Eq. (36) and therefore suppressed by powers of $\Delta r$. Similarly, we find that for any short-distance four-nucleon operator

[TABLE]

the ratio of expectation values is

[TABLE]

We can now turn Eq. (37) and Eq. (39) around and conclude that the ratios ${\rho_{3}}/{\rho_{3,\alpha}}$ and ${\rho_{4}}/{\rho_{4,\alpha}}$ are independent of renormalization scale up to higher dimension corrections,

[TABLE]

It would be interesting to check this statement of model independence in the future using a variety of different lattice and continuum ab initio methods.

Local cluster operators as a measure of clustering

Let us consider the short-distance operators $X_{3}(\vec{r})$ and $X_{4}(\vec{r})$ as defined in Eq. (35) and Eq. (38) respectively. As the width of the spatial distributions $\Delta r$ becomes small, the expectation values of these short-distance operators will depend very strongly on the amount of clustering present in the nucleus. If the nucleus is a homogeneous liquid of uncorrelated nucleons then

[TABLE]

where $R_{A}$ is the radius of the nucleus $A$. If on the other hand, the nucleus is comprised of non-overlapping alpha clusters, then

[TABLE]

where $R_{A}$ is the radius of an alpha particle $A$. Therefore if we measure ${\rho_{3}}/{\rho_{3,\alpha}}$ there is an enhancement by a factor of $(R_{A}/R_{\alpha})^{6}$ if the nucleus is comprised of alpha clusters. For ${\rho_{4}}/{\rho_{4,\alpha}}$ the enhancement factor is $(R_{A}/R_{\alpha})^{9}$.

Pinhole algorithm

Auxiliary-field Monte Carlo simulations are efficient for computing the quantum properties of systems with attractive pairing interactions. By the calculating the exact quantum amplitude for each configuration of auxiliary fields, we obtain the full set of correlations induced by the interactions. However, the exact quantum amplitude for each auxiliary field configuration involves quantum states which are superpositions of many different center-of-mass positions. Therefore information about density correlations relative to the center of mass is lost. The pinhole algorithm is a new computational approach that allows for the calculation of arbitrary density correlations with respect to the center of mass. As this was not possible in all previous auxiliary-field Monte Carlo simulations, adaptations of this technique should have wide applications to hadronic, nuclear, condensed matter, and ultracold atomic simulations.

We let $\rho_{i,j}({\bf n})$ be the density operator for nucleons with spin $i$ and isospin $j$ at lattice site n,

[TABLE]

We construct the normal-ordered $A$-body density operator

[TABLE]

In the $A$-nucleon subspace, we note the completeness identity

[TABLE]

Using the transfer matrices $M$ and $M_{*}$ defined in Eq. (19) and Eq. (21), in the pinhole algorithm we work with the expectation value

[TABLE]

Due to the completeness identity Eq. (48), the sum of the expectation value in Eq. (49) over ${\bf n}_{1},\cdots{\bf n}_{A}$ and $i_{1},j_{1},\cdots i_{A},j_{A}$ gives $A!$ times the amplitude

[TABLE]

The quantities $Z_{f,i}(i_{1},j_{1},\cdots i_{A},j_{A};{\bf n}_{1},\cdots{\bf n}_{A})$ and $Z_{f,i}$ are computed using Monte Carlo simulations with auxiliary fields. Within the auxiliary-field framework, the pinhole locations ${\bf n}_{1},\cdots{\bf n}_{A}$ and spin-isospin indices $i_{1},j_{1},\cdots i_{A},j_{A}$ are sampled by Metropolis updates Hastings (1970), while the auxiliary fields are sampled by the hybrid Monte Carlo algorithm Duane et al. (1987); Gottlieb et al. (1987). In Fig. S9 we show a sketch of the pinhole locations and spin-isospin indices for the operator $\rho_{i_{1},j_{1},\cdots i_{A},j_{A}}({\bf n}_{1},\cdots{\bf n}_{A})$ inserted at time $t=L_{t}a_{t}/2$. We obtain the ground state expectation value by extrapolating to the limit of infinite projection time. We compute the path integrals

[TABLE]

where $DsD\pi$ is the path integral measure for all time steps of the auxiliary field $s$ and pion field $\pi$, and

[TABLE]

We perform importance sampling of the path integral in Eq. (51) according to the absolute value of the integrand,

[TABLE]

The complex phase of the integrand is treated as an observable that is accumulated to give a total sum over all selected configurations. In the pinhole algorithm we alternate the auxiliary field and pion field updates with updates of the spin-isospin indices and pinhole locations. For fixed indices $i_{1},j_{1},\cdots i_{A},j_{A}$ and pinhole locations ${\bf n}_{1},\cdots{\bf n}_{A}$ we use the hybrid Monte Carlo algorithm Gottlieb et al. (1987); Duane et al. (1987) to update the auxiliary field and pion field. This is the same method used in previous nuclear lattice simulations, and the details of the implementation can be found in Ref. Lee (2009, 2017).

The spin-isospin indices $i_{1},j_{1},\cdots i_{A},j_{A}$ and are updated using the Metropolis algorithm Metropolis et al. (1953). We propose a new set of indices $i^{\prime}_{1},j^{\prime}_{1},\cdots i^{\prime}_{A},j^{\prime}_{A}$ by randomly reassigning the spin and isospin for some of the nucleons. We select a random number $r$ uniformly distributed between [math] and $1$ and accept the new indices if

[TABLE]

We also choose new pinhole locations ${\bf n}^{\prime}_{1},\cdots{\bf n}^{\prime}_{A}$ by randomly displacing one of the pinhole locations by one lattice unit. We select another random number $r$ uniformly distributed between [math] and $1$ and accept the new pinhole locations if

[TABLE]

In this manner we update the auxiliary and pion fields, spin-isospin indices, and pinhole locations.

Density correlations

For spatial lattice spacing $a$, the coordinates ${\bf r}_{i}$ of each nucleon on the lattice is an integer vector ${\bf n}_{i}$ times $a$. We do not consider mass differences between protons and neutrons in these calculations. Since the center of mass is a mass-weighted average of $A$ nucleons with the same mass, the center-of-mass position ${\bf r}_{\rm CM}$ is an integer vector ${\bf n}_{\rm CM}$ times $a/A$. Therefore the density distribution has a resolution scale that is $A$ times smaller than the lattice spacing. In order to determine the center-of-mass position ${\bf r}_{\rm CM}$, we minimize the squared radius

[TABLE]

where each term $\left|{\bf r}_{\rm CM}-{\bf r}_{i}\right|$ is minimized with respect to all periodic copies of the separation distance on the lattice. We comment that the tails of the proton and neutron density distributions are determined from the asymptotic properties of the $A$-body wave function, which have been derived in a recent paper König and Lee (2017) for interactions with finite range.

As discussed in the main text, from the $A$-body density information we can view the triangular shapes formed by the three spin-up protons in the carbon isotopes. The positions of the three spin-up protons serve as a measure of the alpha cluster geometry. In Fig. S10 we sketch a typical configuration of the protons (red) and neutrons (blue) with the arrows indicating up and down spins in 12C. The three spin-up protons form the vertices of a triangle, and this is indicated by the orange triangle in Fig. S10. When collecting the lattice simulation data, we rotate the triangle so that the longest side lies on the $x$-axis. We also rescale the triangle so the longest side has length one, and flip the triangle, if needed, so that the third spin-up proton is in the upper half of the $xy$-plane.

.1 Form factors and radii

From the density distribution of the protons relative to the center of mass, we compute the Fourier transform to determine the electric form factor, $F(q)$, where $q$ is the momentum transfer. In order to reduce systematic errors due to the lattice spacing, we perform a least squares fit of the density distribution using a two-parameter Fermi model,

[TABLE]

and then Fourier transform to momentum space. The results are shown in Fig. S11.

From the density distribution of protons and neutrons, we also compute the root-mean-square (rms) radius for the proton and neutron distributions at leading order. The results are shown in Table 1. The shown error bars include Monte Carlo errors as well as errors due to extrapolation to infinite projection time. For comparison we show the rms charge radius observed in electron scattering experiments. We find reasonable agreement between the 12C and 14C proton radii at leading order and the corresponding observed charge radii.