  We present results for the matrix elements of the additional $\Delta S=2$
operators that appear in models of physics beyond the Standard Model (BSM),
expressed in terms of four BSM $B$-parameters. Combined with experimental
results for $\Delta M_K$ and $\epsilon_K$, these constrain the parameters of
BSM models. We use improved staggered fermions, with valence HYP-smeared quarks
and $N_f=2+1$ flavors of "asqtad" sea quarks. The configurations have been
generated by the MILC collaboration. The matching between lattice and continuum
four-fermion operators and bilinears is done perturbatively at one-loop order.
We use three lattice spacings for the continuum extrapolation: $a\approx 0.09$,
$0.06$ and $0.045\;$fm. Valence light-quark masses range down to $\approx
m_s^{\rm phys}/13$ while the light sea-quark masses range down to $\approx
m_s^{\rm phys}/20$. Compared to our previous published work, we have added four
additional lattice ensembles, leading to better controlled extrapolations in
the lattice spacing and sea-quark masses. We report final results for two
renormalization scales, $\mu=2\;\text{GeV}$ and $3\;\text{GeV}$, and compare
them to those obtained by other collaborations. Agreement is found for two of
the four BSM $B$-parameters ($B_2$ and $B_3^\text{SUSY}$). The other two ($B_4$
and $B_5$) differ significantly from those obtained using RI-MOM
renormalization as an intermediate scheme, but are in agreement with recent
preliminary results obtained by the RBC-UKQCD collaboration using RI-SMOM
intermediate schemes.
