HAL QCD Method for Hadron Interactions

• HAL QCD method is a lattice-QCD framework that extracts hadron interactions through equal-time NBS wave functions and an energy-independent non-local potential.
• It uses a time-dependent formulation and derivative expansion to compute effective potentials without isolating a single ground-state plateau.
• The approach is validated by reproducing scattering phase shifts, finite-volume spectra, and bound-state energies across various baryon and meson systems.

The HAL QCD method is a lattice-QCD framework for extracting hadron–hadron interactions from equal-time Nambu–Bethe–Salpeter (NBS) wave functions, or from their time-dependent correlator analogues, by defining an energy-independent and generally non-local interaction kernel and then solving a Schrödinger-like equation for scattering phase shifts, bound states, and finite-volume spectra. Its distinctive feature is the use of spatial correlations and the information carried by all elastic states, so that the method does not require isolation of a single ground-state plateau in Euclidean time, a condition that is often numerically prohibitive in multi-hadron systems (Iritani et al., 2017, Aoki, 2013).

1. Formal basis

At the core of the method is the equal-time NBS wave function for a two-hadron state. In the two-body case, the asymptotic form of the NBS wave function outside the interaction range carries the physical phase shift, so the wave function can be used as the bridge between QCD and a Schrödinger description. In the original HAL QCD construction below inelastic threshold, one introduces an energy-independent non-local kernel $U(\mathbf{r},\mathbf{r}')$ through

$$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$

with $\mu$ the reduced mass. Solving the corresponding Schrödinger equation in infinite volume then reproduces the same scattering information encoded in the NBS asymptotics (Aoki, 2013).

This potential is non-local in general, but it is defined to be energy-independent in the elastic region. In practical lattice calculations, it is not treated as a unique microscopic object in an operator-independent sense. Rather, the potential depends on the choice of interpolating operators, while different choices are phase-equivalent in the sense that they are constructed to reproduce the same scattering observables. This operator dependence is explicit in HAL QCD studies of meson–baryon systems, where different interpolating choices define different schemes for the same underlying phase shifts (Murakami et al., 2020).

2. Time-dependent formulation and derivative expansion

The form used most widely in contemporary calculations is the time-dependent HAL QCD method. Instead of isolating a single eigenstate, one works with a normalized four-point correlator $R(\mathbf{r},t)$, which in the elastic region can be written as a superposition of NBS wave functions with Euclidean-time weights. For identical baryons, the correlator satisfies

$$\left[ \frac{1}{4m_B}\frac{\partial^2}{\partial t^2} - \frac{\partial}{\partial t} - H_0 \right] R(\mathbf{r},t) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,R(\mathbf{r}',t),$$

so only suppression of inelastic contributions is required; ground-state saturation is not (Iritani et al., 2017).

Because $U(\mathbf{r},\mathbf{r}')$ is non-local, the method employs a derivative expansion,

$$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$

or channel-specific variants such as

$$U(\mathbf{r},\mathbf{r}') \simeq \{V_{LO}(r)+V_{NLO}(r)\nabla^2\}\delta(\mathbf{r}-\mathbf{r}').$$

At leading order, one obtains an effective local potential directly from the correlator and its derivatives. In the $\Xi\Xi(^1S_0)$ case, for example,

$$V_{\mathrm{eff}}(r) = \frac{1}{4m_B}\frac{\partial_t^2R(\mathbf{r},t)}{R(\mathbf{r},t)} - \frac{\partial_tR(\mathbf{r},t)}{R(\mathbf{r},t)} - \frac{H_0R(\mathbf{r},t)}{R(\mathbf{r},t)}.$$

A quantum-mechanical study using a highly non-local separable potential sharpened the interpretation of this expansion: it concluded that the derivative expansion in HAL QCD is best regarded as a method to extract physical observables such as phase shifts and binding energies, rather than as a literal reconstruction of the microscopic kernel; in that setting, the observable-level convergence was substantially faster than in a formal derivative expansion of the exact non-local potential itself (Aoki et al., 2021).

3. Direct method, fake plateaux, and the “mirage problem”

The most persistent methodological contrast is with the direct method based on temporal plateaux. There one forms a two-hadron correlator $$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$0, divides by single-hadron correlators, and interprets a plateau in the effective energy shift

$$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$1

as the ground-state energy shift. In multi-baryon systems, however, the elastic level spacing scales as $$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$2, so low-lying elastic excitations decay very slowly; plateaux appearing at accessible Euclidean times can therefore be pseudo-plateaux generated by admixtures of nearby states rather than by ground-state dominance (Iritani et al., 2016).

HAL QCD analyses of $$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$3 made this point quantitatively. In that system, the direct-method effective energy shift exhibited strong source dependence between wall and smeared sources, so the apparent plateau could not be a physical eigenenergy. Using eigenfunctions and eigenvalues obtained from the HAL QCD potential, the correlator was decomposed into a few low-lying modes, and the reconstructed effective energy shift reproduced the plateau-like behavior seen in the direct method. For the smeared source, the first excited-state contamination was about $$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$4, and the time needed to reach the true ground-state plateau was estimated as $$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$5 fm, far beyond practical reach; the plateau near $$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$6 was therefore identified as a mirage produced by elastic excited states (Iritani et al., 2017).

Finite-volume eigenmode projection later generalized this diagnostic. For both $$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$7 and $$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$8 in the $$(E_W-H_0)\,\phi_W(\mathbf{r}) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,\phi_W(\mathbf{r}'), \qquad H_0=-\frac{\nabla^2}{2\mu},$$9 channel, eigenmode-projected correlators yielded source-independent plateaux for the ground and first excited states that matched the finite-volume Hamiltonian eigenvalues obtained from the HAL QCD potential. In particular, the $\mu$0 correlator was dominated by the first excited state rather than the ground state, yet the HAL QCD potential still reproduced the correct ground-state binding energy, providing explicit evidence that excited-state dominance in the correlator does not invalidate the potential method (Doi et al., 2021).

4. Systematic uncertainties and internal validation

A central issue for the method is the control of the derivative expansion. In the $\mu$1 analysis, the effective leading-order potential from the wall source was nearly time independent and closely matched the extracted $\mu$2, while the smeared-source effective potential showed larger time dependence that could be interpreted as higher-order terms in the derivative expansion. The resulting low-energy phase shifts from $\mu$3, $\mu$4, and $\mu$5 agreed within statistical errors at small $\mu$6, indicating good low-energy convergence (Iritani et al., 2017).

This low-energy convergence was tested more directly in finite volume for $\mu$7 and $\mu$8. There, leading-order HAL QCD potentials generated finite-volume eigenmodes whose projected correlators reproduced both the ground and first excited spectra. The fact that the same low-order potential correctly described distinct low-lying eigenstates in two systems with shallow bound states was taken as evidence that the derivative expansion is well converged in the relevant low-energy range (Doi et al., 2021).

On the lattice, cubic symmetry rather than $\mu$9 introduces another systematic: an $R(\mathbf{r},t)$0 projection contains not only $R(\mathbf{r},t)$1 but also higher partial waves such as $R(\mathbf{r},t)$2. Applying Misner’s approximate partial-wave decomposition to a $R(\mathbf{r},t)$3 NBS wave function showed that the $R(\mathbf{r},t)$4 component can be extracted successfully from $R(\mathbf{r},t)$5-projected data. The study further found that the comb-like structures often seen in conventional HAL QCD potentials arise from higher partial-wave contamination amplified by the usual second-order difference Laplacian; when the Laplacian was evaluated analytically on each partial wave after Misner decomposition, those structures disappeared, while the fitted potentials and resulting phase shifts remained almost identical to the conventional analysis. This indicated that higher partial-wave contamination in S-wave HAL QCD studies was already well controlled at the level of observables (Miyamoto et al., 2019).

Consistency with finite-volume quantization is another recurrent check. In the $R(\mathbf{r},t)$6 case, the HAL QCD potential produced finite-volume energy shifts proportional to $R(\mathbf{r},t)$7 and a $R(\mathbf{r},t)$8 curve consistent with Lüscher’s formula. In the moving-frame formulation for $R(\mathbf{r},t)$9$$\left[ \frac{1}{4m_B}\frac{\partial^2}{\partial t^2} - \frac{\partial}{\partial t} - H_0 \right] R(\mathbf{r},t) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,R(\mathbf{r}',t),$$0, effective leading-order potentials and phase shifts extracted in laboratory frames with non-zero total momentum agreed, within larger statistical errors, with both center-of-mass HAL QCD results and the finite-volume method (Aoki et al., 2021).

5. Generalizations and analytic refinements

The formalism was extended above inelastic thresholds and to genuine multi-particle systems by deriving the asymptotic behavior of multi-particle NBS wave functions in terms of the on-shell $$\left[ \frac{1}{4m_B}\frac{\partial^2}{\partial t^2} - \frac{\partial}{\partial t} - H_0 \right] R(\mathbf{r},t) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,R(\mathbf{r}',t),$$1-matrix and by constructing energy-independent coupled-channel potentials in the non-relativistic approximation. In that extension, the large-distance behavior of $$\left[ \frac{1}{4m_B}\frac{\partial^2}{\partial t^2} - \frac{\partial}{\partial t} - H_0 \right] R(\mathbf{r},t) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,R(\mathbf{r}',t),$$2-particle NBS wave functions carries generalized phase shifts and mixing angles, providing the necessary basis for coupled-channel and multi-particle HAL QCD analyses (Aoki, 2013).

The method has also been reformulated in moving frames. For non-zero total momentum, a laboratory-frame NBS wave function can be related to the energy-independent HAL QCD potential defined in the center-of-mass frame, and a corresponding time-dependent extraction formula can be written directly in terms of laboratory-frame correlators. In the %%%%76$U(\mathbf{r},\mathbf{r}')$76%%%%4 system, the potentials extracted for $$\left[ \frac{1}{4m_B}\frac{\partial^2}{\partial t^2} - \frac{\partial}{\partial t} - H_0 \right] R(\mathbf{r},t) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,R(\mathbf{r}',t),$$5 and $$\left[ \frac{1}{4m_B}\frac{\partial^2}{\partial t^2} - \frac{\partial}{\partial t} - H_0 \right] R(\mathbf{r},t) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,R(\mathbf{r}',t),$$6 were noisier than in the center-of-mass frame, but the corresponding phase shifts agreed with center-of-mass HAL QCD and with Lüscher’s method. This extension enlarges the applicability of HAL QCD to channels where moving frames are advantageous, including mesonic resonances such as $$\left[ \frac{1}{4m_B}\frac{\partial^2}{\partial t^2} - \frac{\partial}{\partial t} - H_0 \right] R(\mathbf{r},t) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,R(\mathbf{r}',t),$$7 and $$\left[ \frac{1}{4m_B}\frac{\partial^2}{\partial t^2} - \frac{\partial}{\partial t} - H_0 \right] R(\mathbf{r},t) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,R(\mathbf{r}',t),$$8 (Aoki et al., 2021).

A separate line of development applies the method to finite-temperature heavy quarkonium. Using NRQCD correlation functions of non-local mesonic S-wave operators on anisotropic FASTSUM ensembles, the HAL QCD method was used to extract channel-dependent bottomonium potentials and a spin-averaged central potential

$$\left[ \frac{1}{4m_B}\frac{\partial^2}{\partial t^2} - \frac{\partial}{\partial t} - H_0 \right] R(\mathbf{r},t) = \int d\mathbf{r}'\,U(\mathbf{r},\mathbf{r}')\,R(\mathbf{r}',t),$$9

Because NRQCD propagates as an initial-value problem, there are no backward movers, which simplifies the finite-temperature time dependence. In both the preliminary coordinate-space analysis and the later momentum-space/FFT implementation, the central potential was found to flatten as temperature increases, consistent with thermal weakening of the interaction and medium-induced screening (Spriggs et al., 2021, Spriggs et al., 2023).

Analytic questions have also been examined within the HAL QCD framework. A study of the left-hand cut problem considered Yukawa potentials with an infrared cutoff $U(\mathbf{r},\mathbf{r}')$0 and showed that the $U(\mathbf{r},\mathbf{r}')$1-matrix and phase shifts are well-defined at finite $U(\mathbf{r},\mathbf{r}')$2 for complex momentum $U(\mathbf{r},\mathbf{r}')$3. In the $U(\mathbf{r},\mathbf{r}')$4 limit, the phase shift approaches the analytically continued result for $U(\mathbf{r},\mathbf{r}')$5, while it differs for $U(\mathbf{r},\mathbf{r}')$6 except at the binding momentum $U(\mathbf{r},\mathbf{r}')$7. The analysis suggests that, when a bound state lies below the left-hand-cut branch point, the appropriate procedure in HAL QCD is to solve the Schrödinger equation with the extracted potential directly, rather than relying on naive analytic continuation of $U(\mathbf{r},\mathbf{r}')$8 (Aoki et al., 28 Jan 2025).

The conceptual scope of the method has even been expanded beyond color-singlet two-hadron systems. In an extended HAL QCD treatment of a static-quark–diquark baryonic system, a scalar-diquark mass was determined self-consistently by requiring that a $U(\mathbf{r},\mathbf{r}')$9-wave spectrum from two-point correlators be reproduced by a potential extracted from the $$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$0-wave system. The resulting quark–diquark potential had Cornell form, and the scalar diquark mass was found to be roughly $$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$1 (Kelvin-Lee et al., 15 Jan 2026).

6. Representative applications

The method has been applied across baryon–baryon, meson–baryon, meson–meson, and heavy-quark systems. Selected results are summarized below.

System	Setting	Representative outcome
$$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$2	Nearly physical quark masses, $$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$3 MeV	The interaction is attractive at all distances and produces a quasi-bound state with binding energy $$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$4 MeV; with Coulomb interaction, $$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$5 has binding energy $$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$6 MeV (Iritani et al., 2018)
S-wave $$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$7	$$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$8 MeV, all-to-all propagators with one-end trick	The $$U(\mathbf{r},\mathbf{r}') = \left[ V_0(r)+V_2(r)\nabla^2+V_4(r)\nabla^4+\cdots \right]\delta(\mathbf{r}-\mathbf{r}'),$$9 channel is more strongly repulsive than $$U(\mathbf{r},\mathbf{r}') \simeq \{V_{LO}(r)+V_{NLO}(r)\nabla^2\}\delta(\mathbf{r}-\mathbf{r}').$$0, and the phase shifts qualitatively reproduce the energy dependence of experimental phase shifts (Murakami et al., 2020)
$$U(\mathbf{r},\mathbf{r}') \simeq \{V_{LO}(r)+V_{NLO}(r)\nabla^2\}\delta(\mathbf{r}-\mathbf{r}').$$1, $$U(\mathbf{r},\mathbf{r}') \simeq \{V_{LO}(r)+V_{NLO}(r)\nabla^2\}\delta(\mathbf{r}-\mathbf{r}').$$2 in $$U(\mathbf{r},\mathbf{r}') \simeq \{V_{LO}(r)+V_{NLO}(r)\nabla^2\}\delta(\mathbf{r}-\mathbf{r}').$$3	Nearly physical light quark masses and physical charm	Finite-volume eigenmode analysis showed shallow bound states and demonstrated that the HAL QCD potential reproduces both ground and first excited energies even when the unprojected correlator is dominated by an excited scattering state (Doi et al., 2021)
%%%%78$\Xi\Xi(^1S_0)$78%%%%5 in moving frames	Non-zero total momentum	Effective leading-order potentials and phase shifts in laboratory frames agree with center-of-mass HAL QCD and with the finite-volume method, albeit with larger statistical errors (Aoki et al., 2021)
%%%%78$\mathcal{O}(1/L^2)$78%%%%7 and %%%%78$10\%$78%%%%9 P-wave	Decuplet-baryon study	The two meson–baryon potentials show quite similar behaviors; the phase shifts indicate $\Xi\Xi(^1S_0)$0 and $\Xi\Xi(^1S_0)$1 bound states in this lattice setup, with binding energies consistent with those from two-point functions (Murakami et al., 2021)
Bottomonium at non-zero temperature	NRQCD on FASTSUM ensembles	The extracted spin-averaged central potential becomes flatter as temperature increases, consistent with screening in the medium (Spriggs et al., 2023)

These applications show that the HAL QCD method is not confined to a single spectroscopy problem. It has been used to analyze weakly bound dibaryons, repulsive meson–baryon channels, shallow bound states in doubly heavy systems, moving-frame meson scattering, decuplet-baryon structure from meson–baryon interactions, and thermal interquark potentials. Across these settings, the recurring methodological pattern is the same: construct NBS-based correlators, extract an energy-independent interaction kernel within a derivative expansion, and validate the result through phase shifts, finite-volume spectra, or bound-state observables.