Progress in Multibaryon Spectroscopy

TL;DR

This paper reviews recent progress in multi-baryon spectroscopy within lattice QCD, highlighting CalLat's strategy and latest results in understanding nuclear interactions from fundamental theory.

Contribution

It introduces CalLat's innovative approach to multi-nucleon spectroscopy and reports new results advancing the understanding of nuclear forces from QCD.

Findings

Progress in multi-nucleon spectroscopy techniques

CalLat's latest results on multi-baryon systems

Advancement towards understanding nuclear interactions from QCD

Abstract

Anchoring the nuclear interaction in QCD is a long-outstanding problem in nuclear physics. While the lattice community has made enormous progress in mesonic physics and single nucleon physics, continuum-limit physical-point multi-nucleon physics has remained out of reach. I will review CalLat's strategy for multi-nucleon spectroscopy and our latest results.

Figures

Equations

$$\det \left[{\left(\mathcal{M}\right)}^{-1}+\delta {\mathcal{G}}^V\right]=0,$$

$$C^{ij}(t)=⟨\Omega |{\mathcal{O}}_{{\Lambda }^{\prime }{\mu }^{\prime },I^{\prime }{m}_I^{\prime }}^{i[J^{\prime }{\ell }^{\prime }{S}^{\prime }]}(t){\mathcal{O}}_{\Lambda \mu I{m}_I}^{j[J\ell S]†}(0)|\Omega ⟩$$

$$C^{ij}(t)=\sum _n{z}_n^i{z}_n^{j†}{e}^{-E_{n}t}$$

$$m(t)$$

$$y(t+\tau )=\hat{T}(\tau )y(t).$$

$$M$$

$$⟨\Omega |\mathcal{O}(\vec{P},\vec{p};t){\mathcal{O}}^{†}(\Delta \vec{x};0)|\Omega ⟩$$

Full text

1 Introduction

Nuclear physics remains quantitatively disconnected from the input parameters of the Standard Model. A long-outstanding dream is to determine the nuclear interaction directly from the fundamental strong interaction of quarks and gluons[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. The difficulty, of course, is that QCD, the theory of those particles’ interactions is analytically intractable at the low energies of nuclear physics and nonperturbative methods are required. The only known method without uncontrollable uncertainties is lattice QCD.

While the experimental coverage of two-nucleon scattering processes is extensive and nuclei with more than two baryons and nuclear reactions are well-measured, a quantitative understanding of nuclear physics directly from QCD remains desirable. Three obvious reasons are, first, to extend the reach of the explanatory power of the Standard Model; second, to enable precision tests of the Standard Model; third, to understand matter in extreme conditions.

Without an understanding of what the Standard Model predicts, it is impossible to quantify any potentially-observed new physics in low-energy nuclear experiments. For example, a complete understanding of any experimental signal of neutrinoless double beta decay requires quantifying how the nucleon axial coupling is modified in a nuclear environment and quantifying how short-distance beyond-the-Standard-Model operators infect the nuclear interaction.

The lattice has long held out the promise of precision non-perturbative QCD calculations. In recent years we have seen the development of continuum-limit, physical-point calculations. In particular, the spectrum of single hadrons is well-determined and dynamical electromagnetism accounted for.

Multi-particle phenomena are much less constrained. Those cases with the greatest precision are entirely mesonic. At first glance, the mesonic case is simpler because the number of required Wick contractions is smaller. On the other hand, many channels have quark-line disconnected diagrams that require stochastic methods. More seriously, the baryon signal-to-noise problem currently encumbers attempts to reproduce the precision of mesonic scattering, and only worsens with baryon number.

Nevertheless, a variety of calculations studying the few-baryon sector exist[12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 9, 32, 33]. In the remainder of these proceedings, I will review the methods and calculations of the CalLat collaboration.

2 Applying the Lüscher Finite-Volume Formalism

The Lüscher finite-volume formalism is the method by which finite-volume energy levels can be transformed into information about infinite-volume scattering[34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 11, 53, 54]. It is quite general and is applied with enormous success in the meson sector. To extract infinite-volume scattering data one must solve the determinant equation

[TABLE]

the quantization condition that the interacting states exactly fit into the box, under the assumption that the box is substantially larger than the range of the interaction. The matrix $\delta\mathcal{G}^{V}$ contains information about the finite-volume spectrum and the finite volume itself, while $\mathcal{M}$ encodes the infinite-volume scattering amplitudes of interest at the respective energies. By changing the size of the box or boosting the center of mass we can access different kinematic points and build a map of the scattering data as a function of momentum. The formal development of the Lüscher method is quite advanced, with an understanding of how to handle inelasticities and coupled channels[43, 55], twisted boundary conditions[56, 57], and other complications and progress towards an understanding of the three-body sector advancing rapidly[58, 59, 60, 61, 62].

The nuclear force is not a central interaction, and orbital and spin angular momenta $L$ and $S$ are not conserved separately. Instead, only their total $J$ is conserved. However, finite-volume energy eigenfunctions must satisfy the quantization condition, which arises from the boundary conditions of the finite volume. The states therefore have the symmetry of the finite volume, rather than the SO(3) symmetry expected from the conservation of angular momentum in an infinite volume. In a cubic box, the states therefore fall into irreducible representations of the octahedral group $O_{H}$. Once the center of mass is boosted, relativistic length contraction further breaks this symmetry and the irreps split into their descendants in the remaining symmetry group. The reduced symmetry implies that finite-volume eigenstates will be superpositions of infinite-volume channels that would not otherwise mix. We say that the infinite-volume states with good angular momentum quantum numbers are subduced onto states labeled by $\mu$ in irreps labeled by $\Lambda$.

The required energy levels may be extracted by any reliable spectroscopic method. In lattice QCD we typically extract the spectrum by fitting Euclidean-time correlation functions. Correlation functions,

[TABLE]

where $i$ and $j$ run over some set of interpolating operators with the quantum numbers of interest and the primed indices label quantities at the sink and unprimed indices quantitites at the source, admit a spectral decomposition

[TABLE]

where the sum is over eigenstates with the quantum numbers of interest, $E_{n}$ the energy of the $n^{\text{th}}$ state in the spectrum, $t$ the Euclidean time, and $z_{n}$ are overlap factors that express how well the interpolating operator gives the $n^{\text{th}}$ state.

In the long-time limit the correlator is dominated by the lowest-lying state, and one can extract the energy by calculating the effective mass,

[TABLE]

where the derivative is typically computed as a finite difference between timeslices on the lattice. Other more sophisticated spectroscopic techniques are available.

Alternatives to the Lüscher formalism include the the potential method and the unitary isobar formalism. An enormous amount of progress both formal and numerical has been achieved, and new results and updates appear elsewhere in these proceedings[63, 64, 65, 66, 67, 68]. Another approach is to put an EFT of interest into a finite volume and tune its low-energy constants to match the spectrum directly, rather than converting the energy levels to infinite-volume observables[69]. The approach taken in this work is the application of Lüscher’s formalism.

In baryonic channels the signal-to-noise problem prevents a clean examination of the late-time limit. While new methods, such as phase reweighting[70], have appeared, their application so far has remained limited. One hope is to find interpolating operators which have very little excited state contamination, so that they plateau early before the noise overwhelms the signal. In the next two sections we discuss the construction of interpolating operators and what may be done to extend plateaus earlier in Euclidean time.

3 Single-Nucleon Improvements

Rather than fight the noise, one may try to remove excited state contamination to produce an early, reliable plateau, by taking linear combinations of interpolating operators. One worrisome source of excited state contamination comes from the inelastic relaxation of the nucleon to its ground state. By improving the individual nucleons and using those nucleons in a multi-nucleon calculation, we can get a better handle on the elastic processes of interest.

A large basis of operators for nucleons has long been known[71] which each give different overlap with the ground and tower of excited states. Since we are only interested in the scattering of nucleons and not any inelastic excited states, we can take the simpler approach of using a point-like operator and a similar, smeared operator. By changing the smearing we can get an operator with markedly different overlap factors from the point operator. By finding a good linear combination of the point and smeared operator we improve the ground state, throwing much of the tower of excited states into the garbage pail of the orthogonal linear combination.

In a truly variational method, the same improvement would be done at the source and the sink in a positive-definite way. However, the high cost of multinucleon contractions has led us to use fixed source interpolators and only improve the sink. To get a procedural improvement of the single nucleon, rather than a so-called ‘$\chi$-by-eye’ improvement, we apply the Matrix Prony method[72, 73], which we summarize here. We hope to extract the ground state of some vector of correlators $y(t)$ which should all have the same ground state. These $y$ might be generated, for example, by using a fixed smeared source and both point and smeared sinks, as we do here. We assume there exists a transfer operator $\hat{T}(\tau)$ such that

[TABLE]

By multiplying both sides on the right by $y{}^{\top}$ and assuming $\hat{T}=M^{-1}V$ we can solve via

[TABLE]

under the assumption that only as many excited states as entries in $y$ contribute from times $t_{i}$ to $t_{f}$. Knowing $M$ and $V$ then allow us to construct $\hat{T}$. Diagonalizing $\hat{T}$ gives orthogonal linear combinations of $y$ with an isolated exponential decay. We observe that the lower-energy state is typically reliable while the orthogonal higher-energy state is typically much more noisy and tends to capture violations of the few-excited-states assumption.

Taking the linear combination that corresponds to the lowest energy state allows us to construct what we call ‘calm baryons’, as they are less excited[74]. With a reliable linear combination the calm baryon correlator plateaus substantially earlier (see, for example, Fig. 1 of Ref. [74]). These linear combinations can be naturally incorporated into baryon block methods[75, 76, 77, 78]. For a single nucleon this doesn’t yield any substantive computational advantage, but in the multinucleon case it reduces substantially the number of contractions needed to eliminate these single-nucleon inelastic excited states. More importantly, incorporating these calm baryons into multinucleon calculations dramatically reduces the excited state contamination, as can be see in Fig. 3 of Ref. [74].

4 Spatially-Displaced Two-Nucleon Operators

The analysis of correlation functions $C^{ij}(t)$ can pose a nontrivial difficulty.

Because they are substantially cheaper in terms of contractions, mesonic calculations have had enormous success applying a variational method, where good combinations of interpolating operators yield very flat effective masses $m(t)$. To be positive-definite, these variational methods require the same set of operators at both source and sink.

The high cost of constructing a complete cross-correlator from a large set of baryonic interpolators has held back calculations with baryons, which have only recently begun applying any truly positive-definite variational method[79].

The variational method relies on interpolating operators projected to definite center-of-mass and definite irrep. In other words, the noninteracting states needed are superpositions of states with baryons of definite momenta. With point-to-all propagators we retain the ability to construct baryons with definite momenta at the sink, but the source is more problematic. While it is easy to project quark sources to a definite momentum, without a combinatorial explosion it’s not simple to project a baryon source to a definite momentum. Without an all-to-all method we’re stuck considering correlators where the source consists of baryons with definite spatial location and the sink consists of baryons with definite momenta and thus are not positive-definite. This indefinite correlation can cause problems of false plateaus[80].

Schematically, such a correlator looks like

[TABLE]

where $\vec{P}$ is the total center of mass momentum, $\vec{p}$ the relative momentum, $t$ the Euclidean time, and $\Delta\vec{x}$ a triplet of integers giving the relative spatial displacement of two baryons. The operators $\mathcal{O}$ themselves are built from two single-baryon operators.

By combining different relative momenta $\vec{p}$ we can project the sink to a noninteracting eigenstate. After selecting an irrep $\Lambda$ one then can select an $n^{2}$, where $\vec{n}$ is a triplet of integers that corresponds to the momentum $\vec{p}$. Increasing $n^{2}$ increases the energy of the target noninteracting state. A comprehensive set of states of good cubic symmetry was tabulated by Luu and Savage[81]111The arXiv:1101.3347 version contains substantially more tables and useful information, which the referee requested removed in the published version, according to private correspondence..

On the source side we can construct operators with good overlap with particular irreps. In the simplest case $\Delta\vec{x}=0$ we get overlap with $A_{1}^{+}$ of all $n^{2}$ and no other irrep. As $(\Delta x)^{2}$ increases one can construct more and more complicated operators and, therefore, irreps.

Spatial displacements exhibit different types of symmetry under the octahedral group. The zero displacement, of course, is completely symmetric; to access parity-odd partial waves one must avoid the accidental parity symmetry that occurs when two baryons are created on the same lattice site. Under the action of the octahedral group a displacement along a lattice direction $\Delta\vec{x}\sim(0,0,1)$ generates a total of six displacements, $(0,1,0)$, $(1,0,0)$, and all their additive inverses. Displacements proportional to $(0,1,1)$ yield a family of twelve displacements (three rearrangements, each with four possible sign assignments for the nonzero entries). The different shapes and their respective number of vertices, which corresponds to the cost of constructing a source of that shape are detailed in Table 3. The names we give to different displacements correspond to locations on a cube surrounding the origin, as shown in Figure 3.

By scaling up that cube we spread the nucleons apart. By using different sources we achieve different overlaps with different cubic irreps. We have previously studied the agreement between different displacements and face, edge, and corner sources[82, 32]. By studying the overlap of faces, edges, and corners onto noninteracting states as a function of cube size we can understand how to build sources that have good overlap with many different channels of interest.

For example, restricting a calculation to a local source gives very good overlap with $A_{1}^{+}$ but no overlap with other irreps (it is particularly clear that access to the parity-odd irreps is exactly zero by symmetry). If we picked the maximal possible displacement, $\Delta\vec{x}=(1,1,1)L/2$, where $L$ is the linear size of the cubic volume, we need only perform one additional solve, as that antipodal displacement goes to itself under the action of the octahedral group, once accounting for periodic boundary conditions. This source too gives overlap with $A_{1}^{+}$ but no overlap with parity-odd or other parity-even irreps.

Note that a cube around the origin is also a cube of a different size around the the antipode. This suggests a cost-saving measure: if we do ten solves—at the origin, the antipode, and on a set of corners—we can build a large variety of sources. Obviously, we can build 10 point sources, an antipodal source, a corner source around the origin on a cube of size $\Delta x$, a corner source around the antipode on a cube of size $L/2-\Delta x$. By combining two solves from a corner source we can construct the displacements for 2 faces, 1 edge, and $\frac{1}{2}$ a corner source on a cube of size $2\Delta x$, albeit not around a common point. However, as the contraction costs are not trivial, we have not yet taken advantage of these other sources.

By trying to maximize the overlap between our position-space displaced sources and the lowest noninteracting states in a variety of cubic irreps for both the cubic source around the origin (OC) and the cubic source around the antipode (AC) one concludes that $\Delta x=L/8$ (or equivalently $3L/8$) is best. Picking the symmetric $\Delta x=L/4$ often yields a zero overlap with parity-odd cubic states.

Using those displacements for a $24^{3}\times 64$ lattice described in Section 5 with a calm baryon interpolating operator as previously discussed we can study four different correlation functions in the isotriplet $A_{1}^{+}$ $n^{2}=0$ which corresponds to the ${}^{1}S_{0}$ dineutron channel, and the isosinglet $A_{1}^{+}$ $n^{2}=0$ (${}^{3}S_{1}$ deuteron) as shown in Figure 4. We see that even with calm baryons, the entirely local operator doesn’t plateau before hitting the noise. This may not be surprising, as putting two baryons on top of one another should cause large inelastic distortion, while the tuning of calm baryons kills the excited states of a single nucleon in isolation.

In contrast, the spatially displaced operators have visible, compatible plateaus. Moreover, the maximally-displaced and OC sources have similar correlation functions while AC has some additional excited state contamination, in comparision. Again, this is compatible with the idea that calming the baryons is more effective when they are far apart, as AC is on a cube with $\Delta x=L/8\sim 0.36$fm while the nucleons in the OC source are more widely separated, $\Delta x=3L/8\sim 1.08$fm.

5 Preliminary Results at $m_{\pi}\sim 350$MeV

We generated an ensemble of 10,000 $24^{3}\times 64$ HISQ gauge configurations[83, 84, 85, 86] with a pion mass of about 350 MeV and a lattice spacing of about 0.12 fm ($m_{\pi}L=5.1$). For the valence quark action we use the Möbius Domain Wall Fermion action described in Ref. [87], with which we previously determined neutrinoless double beta decay matrix elements[88, 89] and the nucleon axial coupling[90, 91] and have been studying the convergence of heavy-baryon $\chi$PT[92], hadronic CP violation[93], methodological improvements[94, 95], and an independent scale setting[96].

Mixed action calculations have a long history(see, eg. [97]); our mixed action has a well-developed EFT (see eg. [98, 99, 100] and references in Ref. [87]) and can take advantage of an excellent GPU solver in QUDA[101]. Moreover, the chiral properties of the MDWF valence action helps maintain chiral symmetry in observables.

Figure 5 shows effective masses with corner sources and exactly projected to the three lowest-energy $A_{1}^{+}$ states for both the isotriplet and isosinglet channels and how they compare to the noninteracting box eigenenergies. The correlation functions were generated from 1 complete corner source (meaning 9 solves—a center an 8 corners) per configuration and were fit with two exponentials. The plateau values are reported in Table 7.

Table 7 shows binding energies in the isosinglet and isovector channels at two nearby pion masses. We emphasize that the negative energy splitting we find should not be directly compared to those binding energies—without additional volumes or negative-energy points one cannot determine an infinite-volume binding energy; a crude phase-shift analysis of our results is consistent with no bound state but further study is warranted. The characteristic lengths $1/\sqrt{\mu B}$ where $\mu$ is the reduced mass and $B$ is the the $n^{2}=0$ energies are about 2.89 fm, while the box size is about 2.88 fm; such a close match in size may indicate a large finite-volume distortion of the wavefunction.

6 Outlook

In addition to the two-nucleon sector, there have been characterizations of the three-baryon and four-baryon sectors. In the calculations thus far performed the baryons are always created on a single site. Figure 4 suggests that these calculations may benefit dramatically from a combination of single-nucleon improvement and spatially displaced sources. As more baryons are added, however, the displaced sources become more expensive.

However, elsewhere in these proceedings the first proof-of-principle calculation of the three-neutron (isospin-$\frac{3}{2}$) system is described[104]. This channel differs from other many-baryon calculations in that the Pauli exclusion principle implies that spatially displaced sources are mandatory, rather than simply beneficial. To construct spin-$\frac{3}{2}$ interpolators, moreover, each baryon must be on its own site[105]. It is natural, then, to try to remove the inelastic single-nucleon excited states in these calculations.

As that calculation matures, the interpretation of the three-body spectrum will present an outstanding theoretical challenge. The formal development of the three-body formalism is rapidly advancing and connects finite-volume spectra with infinite-volume scattering observables[58, 59, 60, 61, 62, 63]. Rather than extracting the infinite-volume physics in a model-independent way and then fit an EFT of choice to those observables, one may instead put the EFT in a matching finite volume and tune its LECs to match the finite-volume spectrum[69]—once the LECs are fixed, one may then use the EFT in an infinite volume. Progress along these lines for analyzing the NPLQCD results at $m_{\pi}\sim 800$MeV with pionless EFT has been reported at the recent Few Body conference[106], while progress using the Harmonic-Oscillator Basis Effective Theory[107] formulated with periodic boundary conditions was described at the April APS meeting in 2017[108].

Whatever method used to go from finite to infinite volume, the extraction of a reliable spectrum is paramount. Ref.[109] details a truly variational calculation in the H-dibaryon channel using distillation at $m_{\pi}\sim$960, and demonstrates clean, reliable plateaus that begin almost immediately, where a simple point-to-all calculation on the same configurations suffers from bad excited state contamination. The impressive success of this calculation heralds exciting promise for future progress in the few-baryon sector.

Acknowledgements

An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program to CalLat (2016) and by the OLCF Director’s Discretionary Time (2017). This research used the NVIDIA GPU-accelerated Titan supercomputer at the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725, and the GPU-enabled Surface and RZHasGPU and BG/Q Vulcan clusters at LLNL. We thank the LLNL Multiprogrammatic and Institutional Computing program for Grand Challenge allocations on the LLNL supercomputers.

This work was supported by the NVIDIA Corporation (MAC), the DFG and the NSFC Sino-German CRC110 (EB), an LBNL LDRD (AWL), the RIKEN Special Postdoctoral Researcher Program (ER), the U.S. Department of Energy, Office of Science: Office of Nuclear Physics (EB, CCC, TK, HMC, AN, ER, BJ, PV, AWL); Office of Advanced Scientific Computing (EB, BJ, TK, AWL); and the DOE Early Career Award Program (DAB, CCC, HMC, AWL). This work was performed under the auspices of the U.S. Department of Energy by LLNL under Contract No. DE-AC52-07NA27344 (EB, ER, PV).

Calculations were performed with the chroma software suite and quda solvers built atop the USQCD stack[110, 111, 112] with the HDF5 library[113, 114] and managed with METAQ and mpi_jm[115, 116]. The a12m350 ensemble was generated with the MILC collaboration’s public lattice gauge theory code (see http://physics.utah.edu/~detar/milc.html).