  Recent observational and theoretical results have suggested that some of
ultraluminous X-ray (ULX) sources may contain neutron star (NS) accretors.
However, the formation channel and properties of donor stars of NS ULXs remain
uncertain. By adopting the non-conservative and rotation-dependent mass
transfer model in the primordial binary evolution, we investigate the way to
form pulsar ULXs like observed pulsar ULXs in a systematic way. Our simulation
results indicate that pulsar ULXs with Be stars and intermediate or/and high
mass donors match observed apparent luminosities, orbital periods and
observationally indicated donor masses of known pulsar ULXs. ULXs with Be and
intermediate donors are main contributors. The route of accretion-induced
collapse of WDs has 4.5% contribution to the NS ULXs, 4.0% of NSs in ULXs are
formed through electron-capture supernovae (SNe), and 91.5% of NSs in ULXs are
born with core-collapse SNe. We also studied the evolution of pulsar ULXs to
double compact star systems. We do not find NS-black hole systems (merging in a
Hubble time) that evolved from pulsar ULXs. Pulsar-white dwarf (WD) cases that
evolve through pulsar ULXs have significant contributions to the whole NS-WD
gravitational wave sources. Contributions of pulsar-WD and pulsar-NS cases that
experienced pulsar ULXs are $\sim$40% and 11% among all LISA NS-WD and NS-NS
sources, respectively. Monte Carlo simulation noise with different models give
a non-negligible uncertainty.
