  We present a new fully first order strongly hyperbolic representation of the
BSSN formulation of Einstein's equations with optional constraint damping
terms. We describe the characteristic fields of the system, discuss its
hyperbolicity properties, and present two numerical implementations and
simulations: one using finite differences, adaptive mesh refinement and in
particular binary black holes, and another one using the discontinuous Galerkin
method in spherical symmetry. The results of this paper constitute a first step
in an effort to combine the robustness of BSSN evolutions with very high
accuracy numerical techniques, such as spectral collocation multi-domain or
discontinuous Galerkin methods.
