  We analyze neutrinoless double beta decay ($0\nu\beta\beta$) within the
framework of the Standard Model Effective Field Theory. Apart from the
dimension-five Weinberg operator, the first contributions appear at dimension
seven. We classify the operators and evolve them to the electroweak scale,
where we match them to effective dimension-six, -seven, and -nine operators. In
the next step, after renormalization group evolution to the QCD scale, we
construct the chiral Lagrangian arising from these operators. We develop a
power-counting scheme and derive the two-nucleon $0\nu\beta\beta$ currents up
to leading order in the power counting for each lepton-number-violating
operator. We argue that the leading-order contribution to the decay rate
depends on a relatively small number of nuclear matrix elements. We test our
power counting by comparing nuclear matrix elements obtained by various methods
and by different groups. We find that the power counting works well for nuclear
matrix elements calculated from a specific method, while, as in the case of
light Majorana neutrino exchange, the overall magnitude of the matrix elements
can differ by factors of two to three between methods. We calculate the
constraints that can be set on dimension-seven lepton-number-violating
operators from $0\nu\beta\beta$ experiments and study the interplay between
dimension-five and -seven operators, discussing how dimension-seven
contributions affect the interpretation of $0\nu\beta\beta$ in terms of the
effective Majorana mass $m_{\beta \beta}$.
