  In order to obtain an evolution system which is robust against the violation
of constraints, we present a new set of evolution systems based on the
so-called Baumgarte-Shapiro-Shibata-Nakamura (BSSN) equations.The idea is to
add functional derivatives of the norm of constraints, $C^2$, to the evolution
equations, which was proposed by Fiske (2004) and was applied to the ADM
formulation in our previous study. We derive the constraint propagation
equations, discuss the behavior of constraint damping, and present the results
of numerical tests using the gauge-wave and polarized Gowdy wave spacetimes.
The construction of the $C^2$-adjusted system is straightforward. However, in
BSSN, there are two kinetic constraints and three algebraic constraints; thus,
the definition of $C^2$ is a matter of concern. By analyzing constraint
propagation equations, we conclude that $C^2$ should include all the
constraints, which is also confirmed numerically. By tuning the parameters, the
lifetime of the simulations can be increased as 2-10 times as longer than those
of the standard BSSN evolutions.
