We investigate the appearance of di-neutron bound states in pure neutron matter within the Brueckner-Hartree-Fock approach at zero temperature. We consider Argonne $v_{18}$ and Paris bare interactions as well as chiral two- and three-nucleon forces. Self-consistent single-particle potentials are calculated controlling explicitly singularities in the $g$ matrix associated with bound states. Di-neutrons are loosely bound, with binding energies below $1$ MeV, but are unambiguously present for Fermi momenta below $1$ fm$^{-1}$ for all interactions. Within the same framework we are able to calculate and characterize di-neutron bound states, obtaining mean radii as high as $sim 110$ fm. The resulting equations of state and mass-radius relations for pure neutron stars are analyzed including di-neutron contributions.