Assume initial data is given at some space-like hypersurface which satisfies the constraints there. We use the new evolution system and get a solution in the domain of dependence of the system, (which, if gauges propagate at speeds greater than light, might be smaller than the domain given by the metric). But using the harmonic gauge I know that there is a solution to the Einstein equation on a maximally extended domain of dependence. If one can diffeomorphically transform metric corresponding to that solution into one satisfying the gauge used for the evolution with the new system, then, since it satisfies all the equations, including the constraints, it follows that it will also satisfy the equations of the new system. Uniqueness of solutions of the new system implies it must be the one found initially and so it also satisfies the constraints. Thus we see that no particular consideration for the constraint equations is needed.
© Max Planck Society and the author(s)