Simulations of binary collisions involving compact objects require initial data that satisfy the constraint equations of general relativity. For binary boson star simulations it is common practice to use a superposition of two isolated star solutions to construct an approximate solution to the constraint equations. Such superposed data is simple to set up compared to solving these equations explicitly, but also introduces extra constraint violations in the time evolution. In this work we investigate how physical observables depend on the quality of initial data in the case of head-on boson star collisions. In particular we compare results obtained from data prepared using four different methods: the standard method to superpose isolated stars, a heuristic improvement to this superposition technique and two versions of this data where excess constraint violations were removed through a conformal thin-sandwich solver. We find that differences in the time evolutions are dominated by differences in the way the two superposition methods differ, whereas additionally constraint solving the superposed data has smaller impact. The numerical experiments are conducted using the pseudo-spectral code bamps. Our work demonstrates that bamps is a code suited for generating high accuracy numerical waveforms for boson star collisions due to the exponential convergence in the polynomial resolution of the numerical approximation.
