  The low-energy spectrum and scattering of two-nucleon systems are studied
with lattice quantum chromodynamics using a variational approach. A wide range
of interpolating operators are used: dibaryon operators built from products of
plane-wave nucleons, hexaquark operators built from six localized quarks, and
quasi-local operators inspired by two-nucleon bound-state wavefunctions in
low-energy effective theories. Sparsening techniques are used to compute the
timeslice-to-all quark propagators required to form correlation-function
matrices using products of these operators. Projection of these matrices onto
irreducible representations of the cubic group, including spin-orbit coupling,
is detailed. Variational methods are applied to constrain the low-energy
spectra of two-nucleon systems in a single finite volume with quark masses
corresponding to a pion mass of 806 MeV. Results for S- and D-wave phase shifts
in the isospin singlet and triplet channels are obtained under the assumption
that partial-wave mixing is negligible. Tests of interpolating-operator
dependence are used to investigate the reliability of the energy spectra
obtained and highlight both the strengths and weaknesses of variational
methods. These studies and comparisons to previous studies using the same
gauge-field ensemble demonstrate that interpolating-operator dependence can
lead to significant effects on the two-nucleon energy spectra obtained using
both variational and non-variational methods, including missing energy levels
and other discrepancies. While this study is inconclusive regarding the
presence of two-nucleon bound states at this quark mass, it provides robust
upper bounds on two-nucleon energy levels that can be improved in future
calculations using additional interpolating operators and is therefore a step
toward reliable nuclear spectroscopy from the underlying Standard Model of
particle physics.
