Bow shocks generated by pulsars moving through weakly ionised interstellar medium (ISM) produce emission dominated by non-equilibrium atomic transitions. These bow shocks are primarily observed as Hα 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 ionisation states. We find that though typically Hα 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.