Constraints on the finite volume two-nucleon spectrum at $m_π\approx 806$ MeV

The Big Picture

Imagine trying to understand a knot by looking through a foggy window, knowing only that what you see is the loosest possible tangle. The real knot could be tighter, but never looser. That is roughly where physicists stand when they extract energies from lattice quantum chromodynamics (LQCD), the main computational method for studying nuclear physics from first principles.

The strong nuclear force holds together every atomic nucleus and confines quarks inside protons and neutrons. Yet no one has managed to calculate from first principles whether two neutrons can form a stable pair. That gap has real consequences: it affects how we interpret dark matter experiments, how we understand heavy-element formation in stars, and how we design searches for rare radioactive decays that might reveal new physics.

A team from the NPLQCD collaboration, with members at MIT, Fermilab, and the University of Barcelona, has now taken on one of the persistent sources of confusion in this field. Does the answer you extract from LQCD depend on how you ask the question? Their study tests how the choice of mathematical “operators” (quantum-mechanical probes used to interrogate two-nucleon systems) affects the extracted energy levels.

Key Insight: The variational bounds you extract from lattice QCD depend on which operators you use to probe the system. Some choices give much weaker bounds and can be misleading. This work sorts out which choices are trustworthy and which are not.

How It Works

Lattice QCD discretizes spacetime into a grid and numerically computes how quarks and gluons behave. To study two-nucleon systems, physicists extract energy eigenvalues: the discrete energy levels a proton-neutron or neutron-neutron pair can occupy when confined in a finite computational box. From those levels, they infer scattering properties using the Lüscher method, which connects the finite-volume spectrum to the infinite-volume physics we actually care about.

Energy levels cannot be measured directly, though. Physicists construct interpolating operators, mathematical objects that mimic the two-nucleon state of interest, and then watch how signals built from these operators decay over a computational analog of time. The decay rate encodes the energy.

The variational method goes further. You deploy a whole basis of operators at once and set up a generalized eigenvalue problem (GEVP). The lowest energy found this way is guaranteed to be an upper bound on the true ground-state energy. A tighter bound means you are closer to the truth.

The researchers tested several distinct operator families:

• Plane-wave dibaryon operators: products of two nucleon operators, each carrying definite momentum. The standard workhorse of the field.
• Local hexaquark operators: all six quarks placed at the same spatial point, with no assumption that they form two separate nucleons.
• Negative-parity nucleon operators: built from a specific component of the relativistic quantum description of each nucleon. This type had not been studied in prior work on this system.

The calculation used a lattice with spatial side-length L ≈ 3.4 fm and an artificially heavy pion mass of mπ ≈ 806 MeV, about six times the physical value. Heavier quarks reduce computational cost and sharpen the signal, making this a controlled testbed for method development.

Both the isosinglet channel (I = 0, proton-neutron with spin 1, analogous to the deuteron) and the isotriplet channel (I = 1, two identical nucleons, analogous to the dineutron) were studied. For each channel, the team tested many operator subsets, including bases with as many as 46 operators at once.

Why It Matters

The central result is clean. Operator bases dominated by hexaquark-only operators, or built exclusively from negative-parity nucleon operators, produced variational bounds that were noticeably weaker. Those sets did not capture the low-energy physics efficiently.

All other combinations (plane-wave dibaryon operators, quasi-local operators, and mixtures) produced bounds that agreed within statistical uncertainties. That consistency points to a stable picture of the spectrum, not artifacts of a bad operator choice.

Do bound two-nucleon states exist at this heavy quark mass? The variational method can only deliver upper bounds, and the results show no evidence for a bound state in either channel. They cannot rule one out, either. But what the study accomplishes may matter more in the long run: it sorts out which operator strategies are trustworthy and which are not.

Getting the two-body problem right comes first. The long-term goal of computing nuclear matrix elements from QCD, needed as inputs for dark matter searches and neutrino oscillation experiments, depends on it.

The paper carries a practical warning, too. Not all operator bases are equal. Hexaquark-only results being weaker is not a curiosity; it signals that operator sets with poor overlap onto the physical states of interest can lead you astray. Future LQCD calculations of two-nucleon systems should include diverse, physically motivated operator families and verify consistency across choices.

Open questions remain. The calculation used an unphysical quark mass, and pushing to the physical value mπ ≈ 140 MeV is computationally demanding but necessary. Whether bound dineutron or deuteron states appear at large quark masses, and what that would imply about how nuclear binding changes with quark mass, is still unresolved.

Bottom Line: By stress-testing a wide variety of operator choices in lattice QCD calculations of two-nucleon systems, this work identifies which computational strategies produce reliable energy bounds and which do not. That is a necessary step before anyone can claim trustworthy first-principles nuclear physics from the lattice.