  We determine the neutral kaon mixing matrix element $B_K$ in the continuum
limit with 2+1 flavors of domain wall fermions, using the Iwasaki gauge action
at two different lattice spacings. These lattice fermions have near exact
chiral symmetry and therefore avoid artificial lattice operator mixing.
  We introduce a significant improvement to the conventional NPR method in
which the bare matrix elements are renormalized non-perturbatively in the
RI-MOM scheme and are then converted into the MSbar scheme using continuum
perturbation theory. In addition to RI-MOM, we introduce and implement four
non-exceptional intermediate momentum schemes that suppress infrared
non-perturbative uncertainties in the renormalization procedure. We compute the
conversion factors relating the matrix elements in this family of RI-SMOM
schemes and MSbar at one-loop order. Comparison of the results obtained using
these different intermediate schemes allows for a more reliable estimate of the
unknown higher-order contributions and hence for a correspondingly more robust
estimate of the systematic error. We also apply a recently proposed approach in
which twisted boundary conditions are used to control the Symanzik expansion
for off-shell vertex functions leading to a better control of the
renormalization in the continuum limit.
  We control chiral extrapolation errors by considering both the NLO SU(2)
chiral effective theory, and an analytic mass expansion. We obtain
$B_K^{\msbar}(3 GeV) = 0.529(5)_{stat}(15)_\chi(2)_{FV}(11)_{NPR}$. This
corresponds to $\hat{B}_K = 0.749(7)_{stat}(21)_\chi(3)_{FV}(15)_{NPR}$. Adding
all sources of error in quadrature we obtain $\hat{B}_K =
0.749(27)_{combined}$, with an overall combined error of 3.6%.
