Bow shocks generated by pulsars moving through weakly ionized interstellar medium (ISM) produce emission dominated by non-equilibrium atomic transitions. These bow shocks are primarily observed as H
$\alpha$ nebulae. We developed a package, named Shu, that calculates non-LTE intensity maps in more than 150 spectral lines, taking into account geometrical properties of the pulsars’ motion and lines of sight. We argue here that atomic (C i, N i, O i) and ionic (S ii, N ii, O iii, Ne iv) transitions can be used as complementary and sensitive probes of ISM. We perform self-consistent 2D relativistic hydrodynamic calculations of the bow shock structure and generate non-LTE emissivity maps, combining global dynamics of relativistic flows, and detailed calculations of the non-equilibrium ionization states. We find that though typically
$\text{H}_\alpha$ emission is dominant, spectral fluxes in [O iii], [S ii] and [N ii] may become comparable for relatively slowly moving pulsars. Overall, morphology of non-LTE emission, especially of the ionic species, is a sensitive probe of the density structures of the ISM.