  We present non-perturbative results for beyond the standard model kaon mixing
matrix elements in the isospin symmetric limit ($m_u=m_d$) of QCD, including a
complete estimate of all dominant sources of systematic error. Our results are
obtained from numerical simulations of lattice QCD with $N_f = 2+1$ flavours of
dynamical domain wall fermions. For the first time, these quantities are
simulated directly at the physical pion mass $m_\pi$~$\sim$~$139\,\mathrm{MeV}$
for two different lattice spacings. We include data at three lattice spacings
in the range $a = 0.11 $ - $ 0.07\,\mathrm{fm}$ and with pion masses ranging
from the physical value up to 450$\,\mathrm{MeV}$. Compared to our earlier
work, we have added both direct calculations at physical quark masses and a
third lattice spacing making the removal of discretisation effects
significantly more precise and eliminating the need for any significant mass
extrapolation beyond the range of simulated data. We renormalise the lattice
operators non-perturbatively using RI-SMOM off-shell schemes. These schemes
eliminate the need to model and subtract non-perturbative pion poles that
arises in the RI-MOM scheme and, since the calculations are performed with
domain wall fermions, the unphysical mixing between chirality sectors is
suppressed. Our results for the bag parameters in the $\overline{\mathrm{MS}}$
scheme at $3\,\mathrm{GeV}$ are $B_K~\equiv~\mathcal{B}_1 = 0.5240(17)(54)$,
$\mathcal{B}_2 = 0.4794(25)(35)$, $\mathcal{B}_3 = 0.746(13)(17)$,
$\mathcal{B}_4 = 0.897(02)(10)$ and $\mathcal{B}_5 = 0.6882(78)(94)$, where the
first error is from lattice uncertainties and the second is the uncertainty due
to the perturbative matching to $\overline{\mathrm{MS}}$.
