  Binary boson stars can be used to model the nonlinear dynamics and
gravitational wave signals of merging ultracompact, but horizonless, objects.
However, doing so requires initial data satisfying the Hamiltonian and momentum
constraints of the Einstein equations, something that has not yet been
addressed. In this work, we construct constraint-satisfying initial data for a
variety of binary boson star configurations. We do this using the conformal
thin-sandwich formulation of the constraint equations, together with a specific
choice for the matter terms appropriate for scalar fields. The free data is
chosen based upon a superposition of isolated boson star solutions, but with
several modifications designed to suppress the spurious oscillations in the
stars that such an approach can lead to. We show that the standard approach to
reducing orbital eccentricity can be applied to construct quasi-circular binary
boson star initial data, reducing the eccentricity of selected binaries to the
$\sim 10^{-3}$ level. Using these methods, we construct initial data for
quasi-circular binaries with different mass-ratios and spins, including a
configuration where the spin is misaligned with the orbital angular momentum,
and where the dimensionless spins of the boson stars exceeds the Kerr bound. We
evolve these to produce the first such inspiral-merger-ringdown gravitational
waveforms for constraint-satisfying binary boson stars. Finally, we comment on
how equilibrium equations for the scalar matter could be used to improve the
construction of binary initial data, analogous to the approach used for
quasi-equilibrium binary neutron stars.
