Numerical tests of evolution systems, gauge conditions, and boundary conditions for 1D colliding gravitational plane waves

We investigate how the accuracy and stability of numerical relativity simulations of 1D colliding plane waves depends on choices of equation formulations, gauge conditions, boundary conditions, and numerical methods, all in the context of a first-order $3+1$ approach to the Einstein equations, with basic variables some combination of first derivatives of the spatial metric and components of the extrinsic curvature tensor. Hyperbolic schemes, specifically variations on schemes proposed by Bona and Massó and Anderson and York, are compared with variations of the Arnowitt-Deser-Misner formulation. Modifications of the three basic schemes include raising one index in the metric derivative and extrinsic curvature variables and adding a multiple of the energy constraint to the extrinsic curvature evolution equations. Redundant variables in the Bona-Massó formulation may be reset frequently or allowed to evolve freely. Gauge conditions which simplify the dynamical structure of the system are imposed during each time step, but the lapse and shift are reset periodically to control the evolution of the spacetime slicing and the longitudinal part of the metric. We show that physically correct boundary conditions, satisfying the energy and momentum constraint equations, generically require the presence of some ingoing eigenmodes of the characteristic matrix. Numerical methods are developed for the hyperbolic systems based on decomposing flux differences into linear combinations of eigenvectors of the characteristic matrix. These methods are shown to be second-order accurate, and in practice second-order convergent, for smooth solutions, even when the eigenvectors and eigenvalues of the characteristic matrix are spatially varying.