  We investigate neutrinoless double beta decay ($0\nu\beta\beta$) in the
presence of sterile neutrinos with Majorana mass terms. These gauge-singlet
fields are allowed to interact with Standard-Model (SM) fields via
renormalizable Yukawa couplings as well as higher-dimensional gauge-invariant
operators up to dimension seven in the Standard Model Effective Field Theory
extended with sterile neutrinos. At the GeV scale, we use Chiral effective
field theory involving sterile neutrinos to connect the operators at the level
of quarks and gluons to hadronic interactions involving pions and nucleons.
This allows us to derive an expression for $0\nu\beta\beta$ rates for various
isotopes in terms of phase-space factors, hadronic low-energy constants,
nuclear matrix elements, the neutrino masses, and the Wilson coefficients of
higher-dimensional operators. The needed hadronic low-energy constants and
nuclear matrix elements depend on the neutrino masses, for which we obtain
interpolation formulae grounded in QCD and chiral perturbation theory that
improve existing formulae that are only valid in a small regime of neutrino
masses. The resulting framework can be used directly to assess the impact of
$0\nu\beta\beta$ experiments on scenarios with light sterile neutrinos and
should prove useful in global analyses of sterile-neutrino searches. We perform
several phenomenological studies of $0\nu\beta\beta$ in the presence of sterile
neutrinos with and without higher-dimensional operators. We find that
non-standard interactions involving sterile neutrinos have a dramatic impact on
$0\nu\beta\beta$ phenomenology, and next-generation experiments can probe such
interactions up to scales of $\mathcal O(100)$ TeV.
