  A variational study of two-nucleon systems with lattice quantum
chromodynamics is performed using a wide range of interpolating operators:
dibaryon operators built from products of momentum-projected nucleons,
hexaquark operators built from six spatially localized quarks, and quasi-local
operators inspired by two-nucleon bound-state wavefunctions in nuclear
effective field theories. Correlation-function matrices involving products of
these operators are constructed by computing timeslice-to-all quark propagators
with sparsening techniques. Comparisons between results obtained using the same
gauge-field ensemble but different interpolating-operator sets demonstrate that
interpolating-operator dependence can lead to significant effects on the
two-nucleon energy spectra obtained using both variational and non-variational
methods.
