  We investigate the application of Hybrid Effective Field Theory (HEFT) --
which combines a Lagrangian bias expansion with subsequent particle dynamics
from $N$-body simulations -- to the modeling of $k$-Nearest Neighbor Cumulative
Distribution Functions ($k{\rm NN}$-${\rm CDF}$s) of biased tracers of the
cosmological matter field. The $k{\rm NN}$-${\rm CDF}$s are sensitive to all
higher order connected $N$-point functions in the data, but are computationally
cheap to compute. We develop the formalism to predict the $k{\rm NN}$-${\rm
CDF}$s of discrete tracers of a continuous field from the statistics of the
continuous field itself. Using this formalism, we demonstrate how $k{\rm
NN}$-${\rm CDF}$ statistics of a set of biased tracers, such as halos or
galaxies, of the cosmological matter field can be modeled given a set of
low-redshift HEFT component fields and bias parameter values. These are the
same ingredients needed to predict the two-point clustering. For a specific
sample of halos, we show that both the two-point clustering \textit{and} the
$k{\rm NN}$-${\rm CDF}$s can be well-fit on quasi-linear scales ($\gtrsim 20
h^{-1}{\rm Mpc}$) by the second-order HEFT formalism with the \textit{same
values} of the bias parameters, implying that joint modeling of the two is
possible. Finally, using a Fisher matrix analysis, we show that including
$k{\rm NN}$-${\rm CDF}$ measurements over the range of allowed scales in the
HEFT framework can improve the constraints on $\sigma_8$ by roughly a factor of
$3$, compared to the case where only two-point measurements are considered.
Combining the statistical power of $k{\rm NN}$ measurements with the modeling
power of HEFT, therefore, represents an exciting prospect for extracting
greater information from small-scale cosmological clustering.
