  We present three-dimensional, nonrelativistic, hydrodynamic simulations of
bow shocks in pulsar wind nebulae. The simulations are performed for a range of
initial and boundary conditions to quantify the degree of asymmetry produced by
latitudinal variations in the momentum flux of the pulsar wind, radiative
cooling in the postshock flow, and density gradients in the interstellar medium
(ISM). We find that the bow shock is stable even when travelling through a
strong ISM gradient. We demonstrate how the shape of the bow shock changes when
the pulsar encounters density variations in the ISM. We show that a density
wall can account for the peculiar bow shock shapes of the nebulae around PSR
J2124-3358 and PSR B0740-28. A wall produces kinks in the shock, whereas a
smooth ISM density gradient tilts the shock. We conclude that the anisotropy of
the wind momentum flux alone cannot explain the observed bow shock morphologies
but it is instead necessary to take into account external effects. We show that
the analytic (single layer, thin shell) solution is a good approximation when
the momentum flux is anisotropic, fails for a steep ISM density gradient, and
ap- proaches the numerical solution for efficient cooling. We provide analytic
expressions for the latitudinal dependence of a vacuum-dipole wind and the
associated shock shape, and compare the results to a split-monopole wind. We
find that we are unable to distinguish between these two wind models purely
from the bow shock morphology.
