Clustering of Four-Component Unitary Fermions

TL;DR

This paper uses quantum Monte Carlo methods to study four-component fermions at unitarity, revealing clustering behaviors and near-degenerate ground states, with implications for nuclear physics and ultracold atomic experiments.

Contribution

It introduces novel trial wave functions for four-component fermions and demonstrates their effectiveness in analyzing clustering and energy properties at unitarity.

Findings

Ground state energy nearly equals two four-particle systems

Identification of clustering properties in the system

Extrapolation to the zero-range limit achieved

Abstract

Ab initio nuclear physics tackles the problem of strongly interacting four-component fermions. The same setting could foreseeably be probed experimentally in ultracold atomic systems, where two- and three-component experiments have led to major breakthroughs in recent years. Both due to the problem's inherent interest and as a pathway to nuclear physics, in this Letter we study four-component fermions at unitarity via the use of quantum Monte Carlo methods. We explore novel forms of the trial wave function and find one which leads to a ground state of the eight-particle system whose energy is almost equal to that of two four-particle systems. We investigate the clustering properties involved and also extrapolate to the zero-range limit. In addition to being experimentally testable, our results impact the prospects of developing nuclear physics as a perturbation around the unitary limit.

Equations

$$\hat{H}=-\frac{{\hslash }^2}{2m}\sum _i{\nabla }_i^2+\sum _{i<j}{V}_{i,j}+\sum _{i<j<k}{V}_{i,j,k}$$

$$V_{i,j}=-V_2{\mu }^2\frac{2{\hslash }^2}{m}\exp [-(\mu r_{ij}{)}^2/2]$$

$$\begin{split}V_{i,j,k}&=V_{3}\Big{(}\frac{\mu}{4}\Big{)}^{2}\frac{2\hbar^{2}}{m}\exp[-(\mu R_{ijk}/4)^{2}/2]\end{split}$$

$$E_V\approx \frac{1}{M}\sum _{i=1}^M{E}_L({\text{R}}_i)$$

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

$$\propto |{\psi }_0⟩\underset{\tau \to \infty }{\lim }{e}^{-(E_0-E_T)\tau }$$

$$\psi_{T}^{A}=f_{J}\big{[}\Phi^{{\rm I},{\rm II}}_{BCS}\Phi^{{\rm III},{\rm IV}}_{BCS}+\Phi^{{\rm I},{\rm III}}_{BCS}\Phi^{{\rm II},{\rm IV}}_{BCS}+\Phi^{{\rm I},{\rm IV}}_{BCS}\Phi^{{\rm II},{\rm III}}_{BCS}\big{]}$$

$$\begin{split}\psi_{T}^{B}=&\mathcal{A}\big{[}e^{-\beta_{0}\sum_{\begin{subarray}{c}i=1,3,5,7\end{subarray}}(\mathbf{r}_{i}-\mathbf{r}_{CM}^{1,3,5,7})^{2}}\times\\ &\qquad e^{-\beta_{0}\sum_{\begin{subarray}{c}j=2,4,6,8\end{subarray}}(\mathbf{r}_{j}-\mathbf{r}_{CM}^{2,4,6,8})^{2}}\times\\ &\qquad e^{-\beta_{1}(\mathbf{r}_{CM}^{1,3,5,7}-\mathbf{r}_{CM}^{2,4,6,8})^{2}}(\mathbf{r}_{CM}^{1,3,5,7}-\mathbf{r}_{CM}^{2,4,6,8})^{n}\big{]}\end{split}$$

$$\begin{array}{rl}{\psi }_T^C& =\mathcal{A}[F({\mathbf{r}}_{CM}^{1,3,5,7}-{\mathbf{r}}_{CM}^{2,4,6,8})\times f_J({\mathbf{r}}_1,{\mathbf{r}}_3,{\mathbf{r}}_5,{\mathbf{r}}_7)\times \\ & f_J({\mathbf{r}}_2,{\mathbf{r}}_4,{\mathbf{r}}_6,{\mathbf{r}}_8)\times \prod _{\begin{array}{c}n=1,3,5,7\\ m=2,4,6,8\end{array}}g(r_{nm})]\end{array}$$

$$\begin{array}{r}F({\mathbf{r}}_{CM}^{1,3,5,7}-{\mathbf{r}}_{CM}^{2,4,6,8})={\left(1-\gamma e^{-({\mathbf{r}}_{CM}^{1,3,5,7}-{\mathbf{r}}_{CM}^{2,4,6,8}{)}^2/{\alpha }^2}\right)}^{-1}\end{array}$$

$$\begin{array}{r}{f}_J=\prod _i{e}^{-{\alpha }_J{r}_i^2}\prod _{i<j}\frac{K\tanh ({\mu }_J{r}_{ij})\cosh ({\gamma }_J{r}_{ij})}{r_{ij}}\\ \times \prod _{i<j<k}{e}^{u_J{e}^{-R_{ijk}^2/(2r_J^2)}}\end{array}$$

$$g(r_{nm})={\left(1-{\gamma }_g{e}^{-r_{nm}^2/{\alpha }_g^2}\right)}^{-1}$$

$$\frac{E_8}{E_4}=c_0+\frac{c_1}{\mu R_4(\bar{\mu })}+\frac{c_2}{[\mu R_4(\bar{\mu }){]}^2}.$$

Full text

The study of strongly interacting, ultracold fermionic atoms has witnessed a large number of exciting developments, including the BEC-BCS crossover, the physics of polarons, optical lattices, three- or many-component experiments, and lower-dimensional systems, whether at zero or finite temperature Braaten:2004rn ; Giorgini:2008 ; Levinsen:2014 ; Naidon:2017 . Experimentally the systems are trapped fluids, where pairing correlations typically have a considerable impact on ground-state properties. Some of the most accurate theoretical studies of these systems have involved the use of Quantum Monte Carlo (QMC) or related methods. While early work focused on the two-component problem at or near the unitary regime Carlson:2003 ; Astrakharchik:2004 ; Stecher:2008 ; McNeil:2011 ; Bertaina:2011 ; Shi:2015 ; Anderson:2015 ; Galea:2016 ; Shi:2016 ; Lacroix:2016 , progress has also been made on the study of three- or more-fermionic components Blume:2008 ; Dawkins:2017 ; Drut:2018 . A natural step is to tackle the four-component problem, where one could envision cold-atom experiments that directly probe the strongly interacting regime to address questions such as clustering in a many-particle system; crucially, a four-component cluster is bosonic in nature, thereby increasing the likelihood that heavier systems will be bound.

Experiments with cold fermionic atoms hold the promise (in part already borne out Horikoshi:2017 ) of constraining aspects of nuclear physics which are not amenable to terrestrial experiments. Most obviously, pairing in two-component Fermi gases is directly related to the physics of neutron-star crusts Gezerlis:2008 ; Gezerlis:2010 ; Carlson:2012 ; Gandolfi:2015 ; Lacroix:2017 ; Strinati:2018 . In low-density neutron matter, the components involved are neutrons with spin up ($n\uparrow$) and spin down ($n\downarrow$). Atomic nuclei involve two additional components of protons ($p\uparrow$, $p\downarrow$). Although the two-nucleon system is not exactly at unitarity, it has been argued that the properties of few-nucleon systems can be obtained in perturbation theory around this limit Konig:2017 ; Konig:2016iny ; Kolck:2017zzf ; Kievsky:2018 ; Rupak:2018gnc ; Gattobigio:2019 ; Konig:2019xxk . Four-component fermionic systems thus might shed light on issues like clustering in light nuclei and the convergence to the thermodynamic limit, i.e., nucleonic matter.

In the unitary limit where the two-body scattering length diverges and the binding energy is vanishingly small, the two-body system is scale invariant and parameter-free. At or close to the unitarity limit, weakly bound systems can be described with an effective field theory (EFT), contact (or pionless) EFT Hammer:2019poc , where all interactions are of zero range. For bosons and three- or more-component fermions, the three-body system collapses  Thomas:1935zz unless a three-body interaction prevents it Bedaque:1999a ; Bedaque:1999b ; as a result, heavier systems will also collapse. The three-body force contains a single dimensionful parameter and gives rise to the Efimov effect Efimov:1970zz thanks to a remaining discrete scale invariance. Systems with more bosons also display the consequences of discrete scale symmetry Platter:2004he ; Hammer:2006ct ; Schmidt:2009kq ; Deltuva:2010xd ; Stecher:2010 ; Nicholson:2012 ; Kievsky:2014 ; Yan:2015 ; Bazak:2016wxm ; Bazak:2018qnu and saturate at finite density Carlson:2017 .

While the four-nucleon system is well described in pionless EFT Platter:2004zs ; Stetcu:2006ey ; Kirscher:2009aj ; Lensky:2016djr ; Konig:2017 ; Contessi:2017 ; Bansal:2018 , results for heavier nuclear systems (6Li, 16O, 40Ca) have so far been somewhat disappointing: although binding energies are within the error of the leading-order calculation, systems are unstable with respect to breakup into 4He and 2H clusters Stetcu:2006ey ; Contessi:2017 ; Bansal:2018 . As one pushes pionless EFT to larger systems, the issue arises of whether the details of the interactions matter. In recent decades the strongest contender to account for these details has been chiral EFT Hammer:2019poc , where contact interactions are supplemented by exchanges of the lightest mesons (pions). Nuclear many-body calculations now routinely employ interactions from chiral EFT with much success (for example Refs. Hebeler:2010 ; Gezerlis:2013 ; Coraggio:2013 ; Hagen:2014 ; Gezerlis:2014 ; Carbone:2014 ; Roggero:2014 ; Wlazlowski:2014 ; Soma:2014 ; Elhatisari:2015 ; Tews:2016 ; Piarulli:2018 ; Lonardoni:2018 ; Lynn:2019 ; Stroberg:2019 ), but they do not saturate properly at leading order Ekstrom:2017koy ; Sammarruca:2018bqh ; Machleidt:2009bh . Properties of nuclear binding and clustering have recently been explored in the context of a lattice approach using chiral or minimal interactions Elhatisari:2016 ; Elhatisari:2017 ; Lu:2019 .

Here we consider for the first time four-component unitary systems with more than four fermions. Like the corresponding bosonic systems Carlson:2017 , these systems are expected to be universal in the sense that all energies are given by dimensionless numbers times the energy of the three-body ground state, $E_{3}$, or alternatively the four-body ground-state energy $E_{4}=4.611E_{3}$ Deltuva:2010xd . We report on novel quantum Monte Carlo calculations where the antisymmetrization required to respect the Pauli exclusion principle is carried out explicitly. We use a variety of trial wave functions as we have had to explicitly check and minimize the effect of the fermion-sign problem on our results. We find that the form of the trial wave function plays a considerable role. Thankfully, our approach is variational, meaning that if a trial wave function leads to a low energy, we can discard other guesses for the trial wave function which gave higher values for the energy.

Unlike bosonic systems, four is the maximum number of four-component fermions that can be found in relative $S$ waves. We might then expect clustering when an integer number of four particles are considered. We focus on the simplest such system, containing eight bodies. This is the analog of the 8Be nucleus, whose ground state is observed to be a narrow resonance very close to the two-4He threshold Benn:1967owz ; Wustenbecker:1992 . As we describe in the remainder of this Letter, we find that eight four-component fermions at unitarity do cluster into two four-particle subsystems. This suggests that clustering is a universal feature of weakly bound, multicomponent fermion systems, with tantalizing implications to the wider program of producing nuclear observables as small corrections to the corresponding cold-atom ones. Since the binding energies of nuclei up to 52Fe are close to the energies of the corresponding number of independent alpha particles, a description of the eight-particle system holds the promise of being extensible to heavier systems.

Our Hamiltonian is:

[TABLE]

where $m$ is the particle mass, $V_{i,j}$ is a two-body attractive potential which acts between particles belonging to distinct components, and $V_{i,j,k}$ is a three-body repulsive potential where, again, each particle within a triplet must belong to a distinct component. At unitarity observables should be insensitive to the form of the potentials, as long as their ranges are small compared to interparticle distances. We take Gaussian forms with a common range $\mu^{-1}$:

[TABLE]

[TABLE]

where $r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$ and $R_{ijk}=(r_{ij}^{2}+r_{ik}^{2}+r_{jk}^{2})^{-1/2}$. The strengths $V_{2}$ and $V_{3}$ are adjusted to ensure, respectively, two-body unitarity and a nonzero four-body energy $E_{4}$; as a preliminary check, we ensured that our results for $E_{4}$ matched the four-boson values from Ref. Carlson:2017 . To produce dimensionless quantities, we employ $E_{4}$ or, alternatively, the corresponding length $R_{4}=(-2mE_{4}/\hbar^{2})^{-1/2}$.

To solve for the ground-state energy of this Hamiltonian we use a combination of the variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC) methods. The VMC method is based on using trial wave functions that attempt to capture the physics of the system being studied. The VMC method evaluates the expectation value of the Hamiltonian using the trial wave functions by computing:

[TABLE]

where $E_{L}(\textbf{R})=\psi^{-1}_{T}(\textbf{R})\hat{H}\psi_{T}(\textbf{R})$ is the local energy, $M$ is the number of sample points, and $\mathbf{R}$ (or $\mathbf{R}_{i}$) encapsulates all the particle positions; this is the result of attacking a multidimensional integral via the Monte Carlo approach. The trial wave functions $\psi_{T}$ contain variational parameters that are adjusted to find a lower VMC energy, carried out via automated optimization techniques Sorella:2001 . Once the VMC energy cannot be lowered further by adjusting these parameters we move onto the Diffusion Monte Carlo method. In DMC we are still making random walks in coordinate space via so-called walkers. However, rather than evaluating a variational estimate of the ground-state energy, we are propagating through imaginary time to project to the ground state $|\psi_{0}\rangle$:

[TABLE]

where we have used imaginary time $\tau=it$ and have decomposed the trial wave function into the energy eigenstates using the completeness relation. The trial energy $E_{T}$ is included as an offset to the Hamiltonian. Of course, when propagating in imaginary time one has to also address the fermion-sign problem; we employ the fixed-node approximation, which implies that our final DMC answers are upper bounds to the true ground-state energy.

In this work specifically, we study the ground-state energy of an eight-particle system. Such a system is made up of eight fermions, with two fermions each belonging to four components (I, II, III, IV). We investigated three possibilities of increasing physical content, so we note ahead of time that Eq. (8) is the main new idea in this Letter. First, we extended the two-component BCS wave function (successfully used in DMC calculations of spin-$1/2$ Fermi gases Carlson:2003 ) to the problem of four components:

[TABLE]

where $\Phi^{m,n}_{BCS}$ is the BCS function that pairs components $m$ and $n$ and $f_{J}$ is a nodeless Jastrow function. $\Phi^{m,n}_{BCS}$ is simply the form used when a system contains only two components; Eq. (6) applies it to all permutations of component pairing. The motivation behind this choice is the success of $\Phi^{m,n}_{BCS}$ in describing the two-component system: via the use of 10 variational parameters, the two-component unitary Fermi gas is described extremely well Carlson:2003 ; McNeil:2011 .

Second, we employed a “cluster wave function,” motivated by the forms used for 8Be in the mid-20th century Wildermuth:1958 :

[TABLE]

where $\beta_{0}$, $\beta_{1}$ and $n$ are variational parameters. In this notation we have arbitrarily labeled the positions of the particles that belong to the component I as $\mathbf{r}_{1}$ and $\mathbf{r}_{2}$, component II has $\mathbf{r}_{3}$ and $\mathbf{r}_{4}$ and so on and we have written the center of mass for particles $a,b,c,d$ as $\mathbf{r}_{CM}^{a,b,c,d}$. In Eq. (7), $\mathcal{A}$ is the antisymmetrization operator, needed to make the form of this wave function appropriately antisymmetric under exchange of identical particles, belonging to the same component. Here, clustering is captured in the first two exponential terms where a positive $\beta_{0}$ causes a decay in the wave function as particles move away from the center of mass of their cluster. Intercluster interactions are dictated by the remaining exponential and polynomial terms, and can be tuned with $\beta_{1}$ and $n$ to favor close mixing of the two clusters, or for the two to remain separate.

Third, we improved the trial wave function of Eq. (7) by extending it to also allow for more complicated correlations:

[TABLE]

where:

[TABLE]

and:

[TABLE]

and, finally:

[TABLE]

The $K$ and $\gamma_{J}$ are chosen such that the function $K\tanh(\mu_{J}r_{ij})\cosh(\gamma_{J}r_{ij})/r_{ij}$ goes to 1 and its derivative goes to 0 at $r_{ij}=d$, where $d$ is referred to as the “healing distance.” In this process, $\alpha$, $\gamma$, $\alpha_{J}$, $\mu_{J}$, $d$, $u_{J}$, $r_{J}$, $\gamma_{g}$ and $\alpha_{g}$ are variational parameters which are adjusted during VMC simulations to find an upper bound on the energy.

The physical motivation behind $\psi_{T}^{C}$ was to capture the behavior seen in the four-particle boson case, which is equivalent to the four-particle fermion case when all fermions belong to distinct components. The $f_{J}$ functions are nodeless and have been used previously in simulations of bosonic clustering. In this context we are using this function to account for the formation of four-particle clusters (where each particle in the cluster belongs to each one of the distinct components), that we suspect will occur in the eight-particle case. Equations (9) and (11) attempt to account for cross-cluster interactions, Eq. (9) is a function of the distance between the centers of mass of the two clusters, while Eq. (11) is simply a function of the separation distances between individual particles belonging to different clusters; both go to unity at large distances. The intercluster pair correlations allow for a deformation of the individual clusters to enhance the attractive interaction between different components between clusters, and to reduce the impact of the Pauli repulsion between like particles. It is this additional correlation that might allow the system to bind. In a Born-Oppenheimer picture, it also suggests that if two $N$-body clusters are bound (here with $N=4$), problems with more than four components will also bind two clusters since the ratio of attractive unlike interactions to Pauli repulsion is increasing with the number of components.

We can qualitatively interpret the physical content of $\psi_{T}^{C}$ by plotting out the paths of the particles over some number of VMC steps. Figure 1 shows the system from the simulation using Eq. (8) over 5000 VMC steps, taken after equilibration. We can see that the system forms two clusters, rather than a larger single cluster containing all of the particles; this is similar to what was seen in 8Be using the nuclear Green’s function Monte Carlo method Wiringa:2000 . While the two clusters here appear to be distinct, it is worth keeping in mind that the length scales involved in $\psi_{T}^{C}$ are such that the two regions communicate with each other; in the DMC method, the two clusters tend to drift apart. Note that in our case the two-body interactions have been tuned at unitarity, namely our calculations correspond to a deuteron with vanishing binding energy. In the figure, we have divided all positions with $R_{4}(\bar{\mu})$, the length corresponding to the range at which $\bar{\mu}R_{4}(\bar{\mu})=13.64$.

A simple benchmark for these calculations is to compare the energy calculated for the eight-particle system to that of a four-particle bosonic system with the same potential parameters. If our simulations have successfully reached the ground state of the eight-particle system we expect the energy found to be a factor of at least two times that of the energy of the corresponding four-particle system. We expect this upper limit because it should be possible for the system to produce two four-particle clusters (that are identical to those found in the four-particle simulation) which do not interact with one another, due to the tiny range of our potentials. If this is the case, the energy of the eight-particle system should be exactly two times the four-particle system, due to the two independent energy contributions of each cluster.

The results of this test are displayed in Fig. 2. In this plot we see the running averages over imaginary time in DMC simulations for all three trial wave functions for a system where $V_{3}=3.0$ and $\bar{\mu}R_{4}(\bar{\mu})=13.64$. We divide our $E_{8}$’s with the large-imaginary-time (central) value of $E_{4}$. As can be seen in the figure, the running averages equilibrate over imaginary time and eventually converge to an answer with an acceptable statistical uncertainty. In this test, the trial wave function from Eq. (8) converged just below the $2E_{4}$ value. This test was carried out for multiple variations of adjusting the potential parameters and in all of these the same outcome was observed. More specifically, we varied the magnitude of $V_{3}$ (see the Supplemental Material supple ): this impacts the value of $E_{4}$, but each time we produced the same ratio $E_{8}/E_{4}$ within statistical error; this suggests that our findings are independent of the details of the short-range interaction.

It is worth highlighting that the computations employing $\psi_{T}^{C}$ were the only ones that led to an eight-particle system whose energy is one or two standard deviations away from breakup into two four-particle clusters. As a matter of fact, the runs with $\psi_{T}^{A}$ did not even come close to the $2E_{4}$ value: despite having ten variational parameters at its disposal, the BCS determinant is designed for a gas and therefore does not do a good job of capturing clustering physics. The $\psi_{T}^{B}$ wave function, on the other hand, that is cluster aware does bring the energy to within 5% of $2E_{4}$; of course, this is still noticeably different than the answer for two four-particle clusters, similarly to what was found in Ref. Contessi:2017 . While physically we know that the system could dissolve into two disconnected clusters, it is worth reiterating that DMC computations impose a nodal Ansatz: if this is sufficiently constricting, then the intuitively expected state of matter may not materialize. This is another way of saying that $\psi_{T}^{C}$ effectively captures both the clustering physics and more involved correlations.

Once the testing of the trial wave function was completed, calculations were performed over a set of varying potentials. We are interested in the zero-range interaction limit, therefore in our simulations we proceeded to also vary $\mu$; we make this dimensionless by forming the product $\mu R_{4}(\bar{\mu})$. In Fig. 3 we plot the ratio of the energies of the eight-particle system and the four-particle system. We can see that as the range of the interactions become smaller, the absolute value of the energy of the eight-particle system becomes larger; this is consistent with the phase-shift expansion, which tells us that the effective range tends to reduce the overall attraction. In order to extrapolate, we fit our DMC results to the form:

[TABLE]

In the limit of $\mu R_{4}(\bar{\mu})$ going to infinity, namely a zero-range interaction, the ratio of the eight-particle to four-particle system goes to $2.04\pm 0.05$ where we carried out standard error propagation. We also checked that employing a higher-degree (or lower-degree) model does not qualitatively change the main result: $E_{8}$ is always within one standard deviation of $2E_{4}$.

In summary, we have extended the DMC approach, which has in the past been applied to the two-component unitary Fermi gas, to four-component unitary fermions. Employing microscopic interactions containing solely two- and three-body central potentials, we have investigated three different types of trial wave function, two of which allow for clustering to emerge. For our most general form of the wave function, we have produced an eight-particle state which is very close to decaying into four-particle clusters. These results constitute an example of pushing the applicability of pionless EFT to heavier systems. One could thus add in a small perturbation which would bring 8Be to the physical point. More generally, our findings could be experimentally tested in the future, when it becomes possible to manipulate four components of strongly interacting fermionic atoms in the lab.

The authors wish to acknowledge several conversations with S. König and D. Lee. This work was supported by the Natural Sciences and Engineering Research Council (NSERC) of Canada, the Canada Foundation for Innovation (CFI), the Early Researcher Award (ERA) program of the Ontario Ministry of Research, Innovation and Science, the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Contracts No. DE-AC52-06NA25396 and No. DE-FG02-04ER41338, and by the European Union Research and Innovation program Horizon 2020 under Grant No. 654002. Computational resources were provided by SHARCNET and NERSC.