  We investigate numerically the propagation of density waves excited by a
low-mass planet in a protoplanetary disk in the nonlinear regime, using 2D
local shearing box simulations with the grid-based code Athena at high spatial
resolution (256 grid points per scale height h). The nonlinear evolution
results in the wave steepening into a shock, causing damping and angular
momentum transfer to the disk. On long timescales this leads to spatial
redistribution of the disk density, causing migration feedback and potentially
resulting in gap opening. Previous numerical studies concentrated on exploring
these secondary phenomena as probes of the nonlinear wave evolution. Here we
focus on exploring the evolution of the basic wave properties, such as its
density profile evolution, shock formation, post-shock wave behavior, and
provide comparison with analytical theory. The generation of potential
vorticity at the shock is computed analytically and is subsequently verified by
simulations and used to pinpoint the shock location. We confirm the theoretical
relation between the shocking length and the planet mass (including the effect
of the equation of state), and the post-shock decay of the angular momentum
flux carried by the wave. The post-shock evolution of the wave profile is
explored, and we quantitatively confirm its convergence to the theoretically
expected N-wave shape. The accuracy of various numerical algorithms used to
compute the nonlinear wave evolution is also investigated: we find that higher
order spatial reconstruction and high resolution are crucial for capturing the
shock formation correctly.
