  We report on the first application of the stochastic Laplacian Heaviside
method for computing multi-particle interactions with lattice QCD to the
two-nucleon system. Like the Laplacian Heaviside method, this method allows for
the construction of interpolating operators which can be used to construct a
positive definite set of two-nucleon correlation functions, unlike nearly all
other applications of lattice QCD to two nucleons in the literature. It also
allows for a variational analysis in which optimal linear combinations of the
interpolating operators are formed that couple predominantly to the eigenstates
of the system. Utilizing such methods has become of paramount importance in
order to help resolve the discrepancy in the literature on whether two nucleons
in either isospin channel form a bound state at pion masses heavier than
physical, with the discrepancy persisting even in the $SU(3)$-flavor symmetric
point with all quark masses near the physical strange quark mass. This is the
first in a series of papers aimed at resolving this discrepancy. In the present
work, we employ the stochastic Laplacian Heaviside method without a hexaquark
operator in the basis at a lattice spacing of $a\sim0.086$~fm, lattice volume
of $L=48a\simeq4.1$~fm and pion mass $m_\pi\simeq714$ MeV. With this setup, the
observed spectrum of two-nucleon energy levels strongly disfavors the presence
of a bound state in either the deuteron or dineutron channel.
