  In this paper a new double-domain spectral method to compute binary black
hole excision initial data is presented. The method solves a system of elliptic
partial differential equations in the exterior of two excised spheres. At the
surface of these spheres, boundary conditions need to be imposed. As such, the
method can be used to construct arbitrary initial data corresponding to binary
black holes with specific boundary conditions at their apparent horizons. We
give representative examples corresponding to initial data that fulfill the
requirements of the quasi-stationary framework, which combines the thin
sandwich formulation of the constraint equations with the isolated horizon
conditions for black holes in quasi-equilibrium. For all examples considered,
numerical solutions with extremely high accuracy were obtained with moderate
computational effort. Moreover, the method proves to be applicable even when
tending toward limiting cases such as large radius ratios for the black holes.
