For highly compact neutron stars, constructing numerical initial data for black hole–neutron star binary evolutions is very difficult. We describe improvements to an earlier method that enable it to handle these more challenging cases. These improvements were found by invoking a general relaxation principle that may be helpful in improving robustness in other initial data solvers. We examine the case of a 6:1 mass ratio system in inspiral close to merger, where the star is governed by a polytropic Γ =2, an SLy, or an LS220 equation of state (EOS). In particular, we are able to obtain a solution with a realistic LS220 EOS for a star with compactness 0.26 and mass 1.98 M_⊙, which is representative of the highest reliably determined neutron star masses. For the SLy EOS, we can obtain solutions with a comparable compactness of 0.25, while for a family of polytropic equations of state, we obtain solutions with compactness up to 0.21, the largest compactness that is stable in this family. These compactness values are significantly higher than any previously published results
