We study the $^1S_0$ superfluidity of $\Lambda$ hyperons in neutron star matter and neutron stars. We use the relativistic mean field (RMF) theory to calculate the properties of neutron star matter. In the RMF approach, the meson-hyperon couplings are constrained by reasonable hyperon potentials that include the updated information from recent developments in hypernuclear physics. To examine the $^1S_0$ pairing gap of $\Lambda$ hyperons, we employ several $\Lambda\Lambda$ interactions based on the Nijmegen models and used in double-$\Lambda$ hypernuclei studies. It is found that the maximal pairing gap obtained is a few tenths of a MeV. The magnitude and the density region of the pairing gap are dependent on the $\Lambda\Lambda$ interaction and the treatment of neutron star matter. We calculate neutron star properties and find that whether the $^1S_0 $ superfluidity of $\Lambda$ hyperons exists in the core of neutron stars mainly depends on the $\Lambda\Lambda$ interaction used.