  We calculate for the first time all four-quark hadronic matrix elements of
local operators possibly contributing to $K\to\pi\pi$ decays and in particular
to the ratio $\epsilon'/\epsilon$ beyond the Standard Model (BSM). To this end
we use the Dual QCD (DQCD) approach. In addition to 7 new mirror operators
obtained from the SM ones by flipping the chirality, we count 13 BSM four-quark
operators of a given chirality linearly independent of each other and of the
aforesaid 14 operators for which hadronic matrix elements are already known. We
present results in two bases for all these operators, one termed DQCD basis
useful for the calculation of the hadronic matrix elements in the DQCD approach
and the other called SD basis suited to the short distance renormalization
group evolution above the 1~GeV scale. We demonstrate that the pattern of long
distance evolution (meson evolution) matches the one of short distance
evolution (quark-gluon evolution), a property which to our knowledge cannot be
presently achieved in any other analytical framework. The highlights of our
paper are chirally enhanced matrix elements of tensor-tensor and scalar-scalar
BSM operators. They could thereby explain the emerging $\epsilon'/\epsilon$
anomaly which is strongly indicated within DQCD with some support from lattice
QCD. On the other hand we do not expect the BSM operators to be relevant for
the $\Delta I=1/2$ rule.
