  This is the first paper of our systematic efforts on lepton number violating
(LNV) hadronic decays in the effective field theory approach. These decays
provide information complementary to popular nuclear neutrinoless
double-$\beta$ ($0\nu\beta\beta$) decay in that they can probe LNV interactions
involving heavier quarks and charged leptons. We may call them hadronic
$0\nu\beta\beta$ decays in short, though $\beta$ refers to all charged leptons.
In this work we investigate the decays $K^\pm\rightarrow\pi^\mp l^{\pm}l^{\pm}$
that arise from short-distance or contact interactions involving four quark
fields and two charged lepton fields, which have canonical dimension nine
(dim-9) at leading order in low energy effective field theory (LEFT). We make a
complete analysis on the basis of all dim-9 operators that violate lepton
number by two units, and compute their one-loop QCD renormalization effects. We
match these effective interactions in LEFT to those in chiral perturbation
theory ($\chi$PT) for pseudoscalar mesons, and determine the resulting hadronic
low energy constants (LECs) by chiral symmetry and lattice results in the
literature. The obtained decay rate is general in that all physics at and above
the electroweak scale is completely parameterized by the relevant Wilson
coefficients in LEFT and hadronic LECs in $\chi$PT. Assuming the standard model
effective field theory (SMEFT) is the appropriate effective field theory
between some new physics scale and the electroweak scale, we match our LEFT
results to SMEFT whose leading effective interactions arise from LNV dim-7
operators. This connection to SMEFT simplifies significantly the interaction
structures entering in the kaon decays, and we employ the current experimental
bounds to set constraints on the relevant Wilson coefficients in SMEFT.
