  We investigate the dynamics of bow shock nebulae created by pulsars moving
supersonically through a partially ionized interstellar medium. A fraction of
interstellar neutral hydrogen atoms penetrating into the tail region of a
pulsar wind will undergo photo-ionization due to the UV light emitted by the
nebula, with the resulting mass loading dramatically changing the flow dynamics
of the light leptonic pulsar wind. Using a quasi 1-D hydrodynamic model of both
non-relativistic and relativistic flow, and focusing on scales much larger than
the stand-off distance, we find that if a relatively small density of neutral
hydrogen, as low as $10^{-4}\,\text{cm}^{-3}$, penetrate inside the pulsar
wind, this is sufficient to strongly affect the tail flow. Mass loading leads
to the fast expansion of the pulsar wind tail, making the tail flow
intrinsically non-stationary. The shapes predicted for the bow shock nebulae
compare well with observations, both in H$\alpha$ and X-rays.
