We present a new numerical scheme to solve the initial value problem for
black hole-neutron star binaries. This method takes advantage of the
flexibility and fast convergence of a multidomain spectral representation of
the initial data to construct high-accuracy solutions at a relatively low
computational cost. We provide convergence tests of the method for both
isolated neutron stars and irrotational binaries. In the second case, we show
that we can resolve the small inconsistencies that are part of the
quasi-equilibrium formulation, and that these inconsistencies are significantly
smaller than observed in previous works. The possibility of generating a wide
variety of initial data is also demonstrated through two new configurations
inspired by results from binary black holes. First, we show that choosing a
modified Kerr-Schild conformal metric instead of a flat conformal metric allows
for the construction of quasi-equilibrium binaries with a spinning black hole.
Second, we construct binaries in low-eccentricity orbits, which are a better
approximation to astrophysical binaries than quasi-equilibrium systems.
