Unequal-mass boson-star binaries: Initial data and merger dynamics

We present a generalization of the curative initial data construction derived for equal-mass compact binaries in Helfer {\it et al} (2019 Phys. Rev. D 99 044046; 2022 Class. Quantum Grav. 39 074001) to arbitrary mass ratios. We demonstrate how these improved initial data avoid substantial spurious artifacts in the collision dynamics of unequal-mass boson-star binaries in the same way as has previously been achieved with the simpler method restricted to the equal-mass case. We employ the improved initial data to explore in detail the impact of phase offsets in the coalescence of equal- and unequal-mass boson star binaries.


Introduction
Perhaps no concept is more central to modern physics than that of the field. Fields are building blocks of our universe and they play a key role in most paradigms of modern cosmology and theories that extend the Standard Model (SM) of particle physics. In recent years, inflationary [1][2][3] and dark matter (DM) models [4,5] have given an important role to scalar fields, which also naturally arise from string theory [6]. If given mass, scalar fields coupled to gravitational field can theoretically form astrophysical, compact, star-like objects. One of the examples of such stars include boson stars (BSs) described by a complex, massive scalar field; see [7,8] for a review. The constituents of a BS are bosonic particles, or bosons (hence, the name), whose mass in the range of 10 −22 − 10 −3 eV has been considered in cosmological and astrophysical settings [9,10]. While the existence of fermionic compact objects -such as neutron stars or white dwarfs -is supported by a plethora of observational evidence (e.g. [11][12][13]), the search for localized solitons made of bosons is still ongoing.
Self-gravitating bosonic fields were first studied in the form of Wheeler's gravitational electromagnetic entities or geons [14]. The concept of BSs, in the sense of equilibrium solutions to the Einstein equations, followed about a decade later with Kaup's pioneering work [15] on self-gravitating configurations of massive complex scalar field, dubbed as Klein-Gordon geons. Originally these configurations were devised from fundamental scalar (spin 0) fields [16,17] and later on extended to vector (spin 1) fields (aka Proca stars) [18][19][20][21][22][23][24][25][26][27][28][29] or high-spin fields [30,31]. The nature of the scalar field can be real -resulting in potentially long lived but not strictly stable compact objects commonly referred to as oscillatons (OSs) [32][33][34][35] -or complex for BSs; this latter complex case is the focus of our work. The first calculations of BSs employed free massive scalar fields, resulting in so-called mini boson stars. Extensive studies over the years, however, have uncovered a rich variety of other BS models, most notably through more elaborate scalar potential functions: self-interacting [36,37], solitonic [38,39] or axionic potentials [9,35,40]. The self-interaction terms lead to significantly more compact BSs, comfortably exceeding the compactness of neutron stars, and also increase the maximum mass BSs may acquire without forming a black hole (BH) [36,41]. Further BS models include charged stars [42,43], BSs comprised of multi-fields [44,45] and rotating models [46,47], where the nature of the spin is quantised. The stability of various bosonic configurations has been addressed in Refs. [48][49][50] and numerous numerical relativity (NR) simulations have demonstrated the robustness of the models [34,[51][52][53].
Due to their (potentially) very high compactness, BSs belong to the category of exotic compact objects and are even regarded as candidates for ultracompact BH mimickers in the sense of possessing a light ring [54,55]. More generally, thanks to their comparatively simple and mathematically regular nature but rich phenomenology, BSs are ideal proxies to study fundamental properties of compact objects using analytic and numerical methods. In this spirit, BSs are also intriguing probes in our search for evidence of extra degrees of freedom in theories of gravity extending general relativity (GR). From an observational viewpoint, BSs have been suggested as alternatives to primordial BHs [56] and supermassive BHs in the centres of galaxies [57]. Last but not least, BSs may contribute to the dark-matter sector of the universe and are an important target for gravitational-wave (GW) observations with the LIGO-Virgo-KAGRA (LVK) network [58][59][60][61][62], as well as future detectors like the Einstein Telescope and Laser Interferometer Space Antenna (LISA) [63,64].
Searches for BS signatures with these GW detectors require accurate waveform models [65,66] whose construction, in turn, relies on extensive high-precision numerical simulations of binary systems involving BSs. The numerical exploration of orbiting BS binary systems is still a relatively young field, but has already demonstrated the potentially rich phenomenology of the GW signals generated by these systems. To our knowledge, the first investigation dates back to Palenzuela et al [67] who find that the BSs' phase offset affects the merger phase more strongly than the inspiral. The GW signal generated by the merger remnant is furthermore mainly governed by the fundamental oscillation frequency of the remnant as it settles down into a non-spinning configuration [68]. Quite remarkably, the GW signal from BS binary mergers can be exceptionally long-lived, resulting in a GW afterglow that decays at a much slower rate than the remnant's angular momentum [53]. The inspiral of unequal-mass BS binaries has been studied in Ref. [52] and can result in large kicks of thousands of km/s which, however, is due to the asymmetric ejection of bosonic matter rather than that of GWs. So-called dark boson-star binaries with purely gravitational interaction have been found to generate GW signatures distinguishable from other astrophysical objects like black holes, neutron stars and even "normal" BSs [69]. Further simulations of compact binaries involving BSs include the piercing of bosonic clouds by a BH [70] and the inspiral of neutron stars with bosonic dark cores [71].
In spite of the tremendous progress made in these numerical explorations, our understanding of the GW emission across the BS binary parameter space remains very limited, both in terms of coverage and precision. One key ingredient indispensable for the systematic construction of GW waveforms forms the central focus of this paper: the generation of accurate initial data representing plausible physical configurations with negligible violations of the Einstein constraint equations. The importance of initial data for binary BS star evolutions in the equal-mass case has been previously addressed by Helfer et al. [72,73], who demonstrate how inaccurate initial superposition of boson stars can lead to substantial spurious features in the resulting gravitational waveforms; to overcome these issues they further propose a new binary superposition that we dub the equal-mass fix. This binary initial data is also applied to the case of equal-mass binary neutron star initial data in the FUKA code [74]. We see here an example how BS studies serve as a valuable proxy well beyond the immediate scope of BS physics. A key limitation of the above cure, however, is its restriction to equal-mass binaries. In this work, we develop a generalized version of this method that achieves the same benefits for binaries with arbitrary mass ratios and contains the equal-mass fix as a limiting case in the choice of two free parameters. The proposed methodology can be applied to head-on collisions as well as systems with angular momentum, as for example in Ref. [53]. Here we employ this improved initial data construction in the simulation of equal-and unequal-mass BS binary head-on collisions studying systematically the impact of the BSs' phase offset on the collision dynamics and resulting GW signals.
The outline of this work is as follows. We start by introducing the theoretical framework and the BS model of interest in Section 2. In Section 3 we summarise the 3+1 split of our equations of motion and the code infrastructures used for our simulations. Section 4 opens with a brief review of the improved construction of initial data in the equal-mass case and proceeds with the generalisation to unequal mass ratios. We explore the parameter space of this initial data construction in Section 5. The results from our exploration of the BSs' phase parameter are presented in Section 6 and we conclude in Section 7. Throughout this work, we set = c = 1 and use M to denote the total mass of the binary. We label spacetime indices by Greek letters running from 0 to 3 and spatial indices by Latin indices running from 1 to 3.

Theoretical framework
Mathematically, boson stars (BSs) are localized, soliton-like 1 solutions of the coupled system of the Einstein and general relativistic Klein-Gordon equations for a complex scalar field ϕ. The action is given by the Einstein-Hilbert term for four-dimensional gravity and a minimally coupled complex scalar field where V (ϕ) is the scalar potential. Varying this action, we recover the Einstein and matter evolution equations and the energy-momentum tensor reads In this work we will focus on the solitonic potential first proposed in [75] V sol = µ 2 |ϕ| 2 1 − 2 where µ is the mass of the scalar field and σ 0 quantifies the field's self-interaction. Note that the solitonic potential has multiple roots in |ϕ|: |ϕ| = 0, which corresponds to the true vacuum state and |ϕ| = σ 0 / √ 2, which represents a "false" or "degenerate" vacuum state. The potential (5) can result in highly compact stars and allows us to span a wider range of mass ratios. Furthermore, solitonic potentials produce some particularly interesting solutions; for example, in the case of ϕ ∼ σ 0 / √ 2 thin-wall configurations have been found, where the scalar field profile acquires a shape almost like a Heaviside function [76]. The resulting soliton profile is then split into three different regions: the interior solution where ϕ ∼ σ 0 / √ 2 (i.e. a false vacuum state), a transition region with a sharp drop from ϕ ∼ σ 0 / √ 2 to ϕ = 0 and the exterior true vacuum state ϕ = 0. In this work, we focus on time evolutions of head-on BS collisions. In general, the outcome of the collision is a non-spinning boson star or a black hole. However, a scenario where the two BSs "pass through" each other [77] is also possible. In our head-on collisions, the resulting remnant is always a BH. We model single BSs as stationary solutions in spherical symmetry, where our ansatz splits the complex solution into amplitude A(r), constant frequency ω ∈ R and phase-offset δφ ∈ [0, 2π) With this ansatz we construct single equilibrium BS solutions using a shooting algorithm. For details of this construction see Section 2.3 of Ref. [73].

The 3+1 decomposition and computational infrastructure
The simulations of BS collisions in this work have been performed with two independent numerical relativity codes, grchombo [78,79] and lean [80]. Both codes evolve the Einstein-Klein-Gordon equations using conformal variants of the 3+1 formalism of Arnowitt, Deser and Misner (ADM) [81] as reformulated by York [82]; cf. also [83].
Here, the spacetime metric is decomposed into the spatial metric γ ij , the shift vector β i and the lapse function α in adapted coordinates x α = (t, x i ) according to and the extrinsic curvature is given by the spatial projection of the covariant derivative of the time like unit normal n α of the foliation, K αβ = −(δ µ β + n µ n β )∇ µ n α . In analogy to the extrinsic curvature, we define a time derivative for the scalar field 2 so that the energy density, momentum density and stress-tensor can be written as These sources appear on the right-hand side of the Einstein equations as given in full detail in Eqs. (8)- (12) of Ref. [73], and the time evolution of the scalar field is given by the first-order-in-time system of equations (15) and (17) in Ref. [73]. Besides the evolution equations, Einstein's equations imply the Hamiltonian, H, and momentum, M i , constraints given by Both codes evolve lapse α and shift β i according to the moving puncture gauge [84,85], i.e. using 1+log slicing and the Γ driver condition as given in Eq. (18) of Ref. [73]. The equations are implemented in the form of finite differencing; more specifically, we use fourth-order spatial differencing with a fourth-order Runge-Kutta method of lines integration in time [86]. Other key ingredients of the codes differ in evolving the Einstein equations, so we summarise them in the following sections and list the grid setups employed for our runs.

Lean
The lean code is based on the cactus computational toolkit [87] and employs mesh refinement in the form of moving boxes as provided by carpet [88]. The code evolves the Einstein equations using the Baumgarte-Shapiro-Shibata-Nakamura-Oohara-Kojima (BSSNOK) formulation [89][90][91], i.e. describes the spacetime in terms of conformally rescaled and trace-split variables whereΓ i mn are the Christoffel symbols associated with the conformal metricγ ij . Apparent horizons are computed using Thornburg's ahfinderdirect [92,93].
The computational domain for all lean simulations consists of 7 nested refinement levels: 4 outer levels centered on the origin and 3 inner levels, each consisting of 2 boxes centered on the two BSs. The box (edge) size decreases from each outer level inwards by a factor of 2 except for level 4 to 5 where it decreases by a factor of 8. The grid spacing dx decreases by a factor 2 on each consecutive inner level. We can thus describe a grid in terms of two numbers, the edge L 1 and grid spacing dx 1 of the outermost level. For the lean simulations of this paper we use 3 L 1 ≈ 1024 M , dx 1 ≈ 2.67 M which implies boxes of size L 7 ≈ 4 M with spacing dx 7 ≈ M/24 on the innermost level.

GRChombo
The GRChombo code is built on the Chombo [94] adaptive mesh refinement (AMR) libraries and evolves the Einstein equations using the covariant and conformal Z4 (CCZ4) formulation [95]. The full Einstein equations in the CCZ4 formulation can be found in Section III F of [96], where we choose κ 1 → κ 1 /α, κ 1 = 0.1, κ 2 = 0 and κ 3 = 1. We set up a grid of length L 1 ≈ 512 M with 7 additional AMR levels such that on the finest level, the spatial grid spacing is dx 7 ≈ M/32. Finally, we use a tagging criterion based on second derivatives of the complex scalar field and the conformal factor (see Section 3.5 of [79] for more details).

Convergence testing
The full details of the convergence testing are provided in Appendix C. In summary, the total error budget, including finite radius extraction and discretisation errors, is 2 % for Lean and 3.7 % for GRChomobo. All of the results reported here use extraction radius R Lean ex ≈ 200 M for Lean and R GRChombo ex ≈ 120 M for GRChombo.

Boson-star binary initial data construction
The construction of binary initial data is a challenging task in GR, mainly due to the non-linearity of the Einstein equations and the gauge dependence of the variables describing the spacetime. For boson stars, we encounter two additional challenges not present for neutron stars or black holes: (i) exponentially growing modes of singlestar solutions and (ii) the lack of a consistent framework for binary initial data that are conformally flat; cf. for example the great simplification afforded by Bowen-York data [97,98]. In this section, we present an ansatz for computing BS-binary initial data that significantly reduce constraint violations relative to the superposition methods used in the literature, and that we also expect to be a valuable preconditioner reducing unphysical features in a full constraint solving process.

Revisiting the initial data construction for equal-mass binaries
Before describing our proposed methodology, we briefly outline the equal-mass initial data construction proposed in [72], which our method then generalises to the unequal mass case. In the remainder of the paper, the initial set-up of our BS binary configurations is as follows. We start off with two BS stars, star A and star B, initially located at x i A and x i B , and therefore separated by an initial distance d = ||x i A − x i B ||. We then boost the stars through Lorentz transformations with initial velocities v i A and v i B towards each other. The details of the 3 + 1 variables with the Lorentz boost can be found in Ref. [73]. We note that the initial positions x i A , x i B and initial boost velocities v i A , v i B are chosen such that the BSs are initially located in the centre of mass frame; these values will be given in our specification of the simulations in Table 2.
The most common procedure for constructing binary initial data is to superpose individual star solutions in a point-wise fashion. In terms of the 3 + 1 ADM variables this is written as This method and, in particular, equation (13) will henceforth be referred to as the method of plain superposition; here the value of δ ij is subtracted from the two superposed individual metrics to ensure the Minkowski metric is recovered in the farfield limit. Whilst the asymptotic flatness condition is thus satisfied, it has been shown in Refs. [72,73] that plain superposition can induce large deviations from the equilibrium values of the volume element, det(γ), at the centres of the stars. This effect arises from the fact that our binary system no longer contains isolated stars and simply superposing metric solutions induces a change in the volume element near the center of each BS due to the influence of its companion (see Appendix A1 of Ref. [72] for more details). To account for such a change in the volume element for the equal-mass binary stars, Ref. [72] proposes to modify the plain superposition by replacing (13) with This modification recovers the equilibrium volume element at the centres of stars A and B as in the case of isolated stars and from now on we will refer to equation (17) as the equal-mass fix. The equal-mass fix has been shown to significantly reduce constraint violations at the centres of the stars relative to the plain superposition procedure and also mitigate spurious physical features such as premature collapse of the stars to a BH and/or altered GW signals [73]. We stress, however, that (17) is only applicable for equal-mass binaries, since in that case In the unequal-mass case, this no longer holds true: the volume element change invoked by each star on its companion will no longer be the same for both stars.

Construction of the generalised unequal-mass initial data
Now, the question arises how to generalise the proposed modification (17) to unequalmass BS binaries. In the spirit of Eq. (17), we have to meet only two conditions: fix the volume element at the centres of each of the stars, x i A and x i B . One way to do so is to work with the 3-metric components γ ij directly and introduce spatially varying corrections that recover the required volume element at the centres of the stars. However, this leads to an under-determined problem, since each metric has 6 components, whilst we have only 2 conditions to satisfy. Instead, we choose to work with the conformal factor, which is a scalar density and thus makes it possible to satisfy both conditions at x i A and x i B using only 2 parameters.
We start with the plainly superposed metric (13) and conformally decompose it with a conformal factor λ defined bỹ where detγ ij = 1 by construction. We note that the conformal factor λ is related to the standard BSSN/CCZ4 variable χ by λ −1 = χ. Appendix A discusses more general choices for the conformal factor and illustrates why the exponent of 1/3 in Eq. (18) is a particularly convenient choice. In our procedure, we leave the conformally rescaled metricγ ij unchanged at the values it takes on in the plain superposition. So far this construction does not remedy the main error inherited from the plain superposition consisting in the change of the volume element at the centres of the stars and the resulting perturbation in their central energy densities. However, we can control the volume element by adjusting the superposed conformal factor λ with a correction δλ, such that we recover the equilibrium volume element at both stars' centres. At the centers of the stars this correction must satisfy where λ A and λ B are the unperturbed conformal factors for stars A and B. In contrast to the equal-mass case, this approach necessitates a spatially varying correction to the conformal factor. For this purpose, we construct weight functions w A (x i ) and w B (x i ) around the centres of the stars, which will be specified later in this Section. We then propose the following ansatz for the new conformal factor on the entire initial hypersurface Here h A and h B are determined by imposing our target conditions (19)- (20), which reduce to a (2 × 2) system of linear equations. The required values of h A and h B are then found to be .
(22) For these constants h A and h B the newly corrected conformal factor λ new allows us to recover the desired volume element at the centres of the stars, as if they were isolated. This can be seen by considering the newly re-scaled metric where γ −1/3 γ ij has unit determinant by construction and therefore ensures that  Table 1. Solitonic BS models with √ Gσ 0 = 0.2 considered in this work. A(0) denotes the central scalar field amplitude, M BS the mass of the BS, ω the frequency of the ground state solution, r 99 the areal radius containing 99% of the BS mass, and max(m(r)/r) our measure of compactness. Note that for this potential the maximum mass of a BS is µM BS = 0.7212 .
We are now left to choose what weight functions to use around the stars in Eq. (21). Focusing here on asymptotically flat spacetimes, we wish to obtain metric corrections that fall off ∝ 1/r. To guarantee such behaviour we construct the following weight functions where J ∈ {A, B}, r J . . = ||x i − x i J || and R J are freely specifiable constants that control the width of the functions. Our initial data method therefore consists of the following main steps: (i) Construct the plainly superposed metric γ ij according to (13).
(ii) Construct the conformal factor λ from a plainly superposed metric according to (18).
(iii) Choose a suitable parameter pair (R A , R B ).
(iv) Compute corrections at the stars' centres, δλ(x A ) and δλ(x B ), according to Eqs. (v) Correct the conformal factor λ to λ new in Eq. (21) to recover the equilibrium volume element at the stars' centres.
The choice of the two parameters, R A and R B , in step (ii) will be explored in more detail through our numerical simulations in the next Section.

Set-up and exploration of the parameter space
In this section, we start by listing the BS models we study in this work in Table 1. We focus on the solitonic BSs, with √ Gσ 0 = 0.2, which allows us to access a variety of mass ratios, including heavy and relatively compact BSs. To construct binary systems of various mass ratios q = M A /M B , where M A < M B , we superpose combinations of the BS models from Table 1. The resulting binary configurations are detailed in Table 2.

Run
Model for star A  Table 2. Binary initial data configurations considered in this work. Each run has a suffix dX-pX, where dX is a wildcard for the initial separation and pX is a wildcard for the off-phase parameter in degrees, 180 • π δφ ∈ [0 • , 360 • ). In our convention, the off-phase parameter, δφ, is added to star A, whilst for star B it remains zero. Here v J x for J ∈ {A, B} denote the initial boost velocities with associated Lorentz factors γ J , BS is the total mass of the binary.

Equal-mass binaries
For equal-mass binaries, the equal-mass fix (17) has been shown to remedy the spurious effects of plain superposition [73]. It is therefore important to verify that our method recovers these improvements in the equal-mass limit. First, we recall that our method requires the choice of free parameters (R A , R B ). As shown in more detail in Appendix B, we can recover both, plain superposition and the equal-mass fix, as limiting cases of this choice. Specifically, in the limit R A , R B → 0, our initial data reconstruction recovers plain superposition, whilst for the equal-mass case and in the limit R A , R B → ∞, we recover exactly the equal-mass fix (17). Our method thus provides a direct generalization for constructing BS binary initial data with arbitrary mass ratios that includes plain superposition and the equal-mass fix as limiting cases. We next verify this claim empirically by evolving in time the binary configuration q1-d11-p000 using four types of superposition: plain, the equal-mass fix and our generalized method using very small and very large (R A , R B ), namely R A = R B = 0.1 and 1000. Figure 1 shows the gravitational waveforms and the central scalar field value as functions of time obtained for these four cases. The figure demonstrates excellent agreement of our proposed method with plain superposition and the equal-mass fix, respectively, for R A = R B = 0.01 and R A = R B = 1000.

Unequal-mass binaries
In the unequal-mass case, the choice of radial parameters (R A , R B ) is more complex. Changing the profile shape of the radial functions affects the extent to which the corrections are applied to the conformal factor in Eq. (21) around the centres of the stars: very small (R A , R B ) will result in smaller corrections around the stars and vice versa. As such, there exists a region of suitable radial parameters, which we find numerically by calculating the L2-norm of the Hamiltonian constraint (10) in the parameter space of (R A , R B ). Figure 2 demonstrates the dependence of the L2-norm of the Hamiltonian  (20)-mode of the Newman-Penrose scalar Ψ 4 for binary sequence q1-d11-p000, obtained for plain superposition, the equal-mass fix and our improved method. The waveforms have been shifted by t peak , the time at which the maximum GW amplitude is reached. Right: Maximum of the scalar field amplitude for the same binary configuration. Notably, plain superposition forms a BH well before the equal-mass fix, therefore resulting in altered gravitational waveform.
Here we indicate the merger time by the vertical lines.
constraint violations on the choice of radial parameters for a binary configuration with q = 0.75. We show the corresponding "heat-maps" for other mass ratios in Appendix D. It is important to note that there are regions in the parameter space, where the initial constraint violations become very large -these regions are indicated by the grey shaded area for readability. The diverging behavior of the constraints for certain pairs of (R A , R B ) can be attributed to a zero crossing of λ new away from the stars' centers, which results in a singular 3-metric. Regions in the parameter space where this happens are evidently not suitable for time evolutions.
In the following, we base our specific choices for radial parameters (R A , R B ) on two criteria. First, the pair has to result in reduced constraint violations relative to plain superposition, and second, they should remedy spurious oscillations in the time evolutions of the scalar field profiles that result from plainly superposed initial data 4 . In this section we focus on constraint violations, whilst we discuss the scalar field profiles in more detail in Section 6.3. For unequal-mass BS binaries, we find that radii R A = 10 − 100 (light BS) and R B = 1 − 100 (heavy BS) work generally well. For smaller mass ratios q 0.5 we find that (R A , R B ) differing by at most one order of magnitude are optimal. We note that different pairs (R A , R B ) can result in small global time-shifts in the GW signals due to gauge effects; besides this time shift, however, the phase and amplitude of the waveforms remain unaffected. We summarise our choices of radial parameters in Table 3 for all unequal-mass binary configurations studied here.
For the radial parameters thus chosen, we observe a particularly pronounced reduction in the constraint violations at the centres of the stars relative to plain superposition. This is as desired by construction of our method, since the correction to the volume element is imposed exactly at the stars' centres. As for the equal-mass fix  studied in Ref. [73], we find this effect to be particularly pronounced in the Hamiltonian constraint. This is illustrated in Figure 3, which shows the constraint violations along the collision axis (x-axis) for the binary q05-d15-p000. Other binary configurations result in qualitatively similar constraint violation profiles.
6. Results: time evolutions of off-phase BS binaries

Equal-mass collisions
We begin the discussion of our results with the simplest case of equal-mass collisions simulated using plain superposition and the equal-mass fix. It has been shown in Ref. [73] that for q = 1 binaries with δφ = 0, plain superposition results in two crucial spurious features: (i) premature BH formation as indicated by a sudden drop in the scalar-field amplitude at the BS center and AH formation; i.e. the two BSs collapse to a BH prior to merger (cf. Figure 9 of Ref. [73]), which can be attributed to the spurious oscillations of the scalar field's central amplitude, (ii) energy dependence on the initial separation, which was found to be less pronounced for solitonic BSs.
Extending the argument of Ref. [73] to the case of arbitrary δφ, we span the dephasing parameter space over the range 5 δφ ∈ [0, π]. Similar to Ref. [73], we find that plain superposition results in premature BH formation regardless of the dephasing parameter δφ. Therefore, in case of plain superposition, the radiated energy becomes independent of the dephasing angle, as is demonstrated in the left panel of Figure 4. In contrast, the equal-mass fix avoids premature BH formation for all choices of δφ and the radiated energy takes on a non-trivial dependence on the dephasing parameter as displayed in the right panel of Fig. 4. This panel furthermore demonstrates that for large δφ the second observation in the above list no longer holds: for δφ 1, the radiated energy does vary significantly with initial separation. As we will discuss in more detail below, the choice of the dephasing parameter can significantly affect the dynamics and GW energy emission for equal as well as for unequal-mass BS binaries.
The role of the phase offset has been studied before in the context of head-on collisions of Proca stars and was found to significantly affect the GW emission and mode structure of the stars [99]. Similarly, in Ref. [100] the phase parameter was found to impact the merger dynamics of real-scalar-field solitons, aka oscillatons (OSs). In particular, over a considerable range of compactness values they find anti-phase (δφ = π) OS head-on collisions to bounce whereas equal-phase (δφ = 0) collisions result  in dispersal of the scalar field; cf. their Fig. 1. A similar repulsive effect has been found for δφ = π in BS head-on collisions in Ref. [101,102], where in particular Appendix B of Ref. [101] gives an explanation in terms of an effective interaction potential. This feature may be connected to the fact that anti-phase collisions produce destructive interference, as shown in the case of Newtonian gravity in Ref. [103]. Our collision sequences q1-dX-pX involve highly compact solitonic stars, and so a BH forms postcollision. We therefore do not observe bounces in the anti-phase collisions, but the scalar field's repulsive character still manifests itself in a weaker signal and reduced radiated energy. Equal-phase configurations form a BH most 'efficiently' and result in the largest energy burst. The energy dependence on the phase off-set is most naturally modelled as a single sinusoidal function 6 ; cf. Eq. (16) in Ref. [103]. This is confirmed by our numerical results in Figure 4, which are well fitted by where amplitude A 1 , frequency f 1 , phase p 1 and shift s are determined using a leastsquares algorithm. From the right panel of Figure 4, it is clear that in the evolutions starting from the equal-mass fix, some energy discrepancy occurs between various separations. Whereas for δφ 1 the radiated energy varies only mildly with d/M , it increases significantly with initial separation for larger dephasing parameters. In all cases, however, we observe a gradual convergence of the energy for large d/M , albeit at distinct rates for different δφ. These effects can be attributed to the increase in the collision velocity that results from larger separations and enhances the merger dynamics.
We can likewise attribute the decrease of E rad for larger δφ to the off-phase scalar fields' repulsion and a consequential weakening of the merger dynamics. However, a higher collision velocity (equivalent to larger initial separation) appears to 'break down' the repellent nature of the scalar field, thus narrowing the energy discrepancy as d → ∞. For the collisions starting from plain superposition, we also observe convergence of the radiated energy in the limit of large separations, but here the energy is a constant function of δφ for each given d.
The increase of energy with separation due to higher collision velocity is a plausible interpretation, supporting our results. But is it correct? We quantitatively test the hypothesis as follows. Using the Newtonian approximation, we estimate that in the evolution starting from initial separation d/ In simple terms, we should obtain a transition of the dark blue curve for d/M = 11.2 in Fig. 4 into the teal colored one for d/M = 22.3 by fixing d/M = 11.2 and increasing the initial velocity from v = 0.1 to v = 0.18. The left panel of Figure 5 illustrates this transition by displaying the energy obtained for equal-mass collisions with initial separation d/M = 11.2 and varying initial boost velocity in the range v ∈ [0.1, 0.18]. For δφ 1 we observe a gradual increase in the energy with rising initial boost velocity. This increase is most pronounced for the anti-phase sequence whereas the velocity change has almost no effect on the radiated energy for in-phase binaries. Over the entire range δφ ∈ [0, π], however, we recover with high accuracy the radiated energies of the d/M = 22.3 binaries by starting from d/M = 11.2 with larger velocity v = 0.18 as predicted by the Newtonian calculation. We are likewise able to recover comparable energies for separations of d/M = 33.5 and d/M = 44.7 by starting the sequence with d/M = 11.2 with yet higher initial boost velocities; see the right panel of Figure 5.

Unequal-mass collisions
In the remainder of this section, we focus on the unequal-mass binary simulations and present a direct comparison between the results obtained for plain superposition (13) and our improved method for q = 1 binaries developed in Section 4.2. There are two limits, in which we would expect the two methods to give comparable results: (i) With decreasing mass ratio, the metric of the lighter star will approach Minkowski, hence the volume factor change induced by it on the heavier star will be negligible. For plainly superposed data, overall constraint violations will be reduced, however, the volume factor change induced on the lighter star by its heavier companion would be inevitable, resulting in spurious star excitations. Again, starting a BS binary from d/M = 11.2 with a larger velocity, as obtained from a Newtonian calculation, the energy functions for v = 0.1 but larger initial separation are recovered. The initial separation d is given in units of M , even though this factor has been omitted in the legends for presentation purposes.
(ii) In the limit of infinite separation (d → ∞), both stars will be isolated and therefore even start in their 'equilibrium' state for plain superposition.
In practice, we are limited to finite mass ratios and initial separations, and as we will see later, the regimes, where plain superposition can give comparable results to our method require very large initial separations, often impractical due to the ensuing computational costs. For this comparison, it is important to realize that in the unequal-mass case, the initial dephasing parameter δφ no longer represents the dephasing at merger, δφ merger . This is a consequence of the two stars' different oscillation frequencies which introduce a "natural" time-dependent phase offset. Introducing a constant phase offset parameter δφ to one of the stars adds a constant phase difference to this time dependent dephasing in a controlled manner. Using multiple values of the initial dephasing parameter δφ for otherwise identical configurations allows us to cover a complete range of dephasing at merger, δφ merger ∈ [0, 2π).
As will be shown in the following analysis, the key effects of varying the off-set parameter δφ for a given BS binary configuration are as follows:  First, however, we test our improved initial data construction by exploring the time evolution of the BSs' central scalar field amplitude.

Scalar field profiles
The deficiencies of the plain-superposition procedure are diagnosed most directly in the time evolution of the scalar-field amplitude at the centres of the two BSs. For equal mass head-on collisions of BSs and OSs, respectively, this has been shown in Fig. 9 of Ref. [73] and Fig. 7 of Ref. [72]. A closer analysis of the scalar-field evolutions we obtain from our equal-mass BS collisions of Sec. 6.1 confirms this picture for all values of the dephasing parameter δφ. In this section, we demonstrate that for sufficiently compact BSs plain superposition results in the same spurious effects as in the unequal-mass BS collisions, namely spurious oscillations of the two stars' central scalar amplitudes around their equilibrium values that subsequently trigger premature collapse to a BH.
In this analysis of unequal-mass collisions, we encounter one minor difficulty, the apriori unknown dephasing at merger. To ensure that our comparison of binaries starting from different initial separations is not adversely affected by possible variations in the dephasing at merger, we choose for each separation the initial phase-parameter that maximizes the radiated GW energy (cf. Fig. 9 below). For the mass ratio q = 0.5 and our two superposition types the specific values of the dephasing parameter can be found in Table 4.
The resulting time evolutions |ϕ c (t)| are shown for all four initial separations and both superposition methods in Fig. 6 and exhibit the same features as mentioned above: plain superposition (left panel) results in significant unphysical oscillations of |ϕ c | which for d/M = 14.8 and d/M = 29.7 also cause a collapse of the heavier star B into a BH well before a common horizon forms as marked by the vertical dashed lines. For the larger separations d/M = 44.5 and d/M = 59.3, the premature BH formation is avoided, but plain superposition still results in significant pulsations of the BSs. In contrast, the scalar-field amplitude of both BSs remains very close to its equilibrium value for our improved superposition method throughout the infall as demonstrated in the right panel of Fig. 6. As expected, significant dynamics in the scalar field, including the eventual collapse to a single BH, are only encountered around merger: the local maximum in the scalar field coincides with common horizon formation for small and large initial separations alike. We have repeated this analysis for different mass ratios and different choices of the dephasing parameter. As it turns out, the dephasing parameter has no significant effect on the results shown in Fig. 6 and our above concern about choosing δφ appropriately has been unnecessary. The mass ratio, however, does affect the results to some extent. The spurious effects of plain superposition are even more pronounced for q = 0.75 (where premature BH formation occurs for all but the largest d/M ) and less pronounced for q = 0.38 (where only the smallest initial separation results in premature BH formation). This q dependence is fully consistent with the above observation (i) in Sec. 6.2 that plain superposition becomes viable for q → 0.

Gravitational waveforms
As we have seen in Sec. 6.1 and, in particular, in Fig. 4, the initial binary separation d can have a significant effect on the magnitude of the GW signal due to the corresponding differences in the collision velocity around merger. This variation arises additionally to the impact due to the choice of the initial dephasing parameter δφ. We can still check the consistency of our evolutions, however, in the limit of large separation d and simultaneously selecting δφ such that it maximizes the GW energy. Specifically, we show the resulting GW signals in Fig. 7 for the q075-dX-pX binaries with initial separation from d/M = 12.7 to d/M = 50.9 and the q038-dX-pX binaries with d/M = 16.2 to d/M = 64.8. The dephasing parameters δφ maximising the GW energy for these configurations are shown in Table 5.
The Figure clearly demonstrates that for mass ratio q = 0.75 (q = 0.38), our improved superposition method results in comparable GW amplitudes for initial   ). For small initial separations, plain superposition systematically results in weaker GW amplitude when compared to our improved superposition method. In fact, this is a distinct feature of premature BH formation that occurs precisely for these plainly superposed configurations; cf. Sec. 6.3. As expected, the agreement in the GW amplitude between the two superposition methods improves as we increase the separation and/or decrease the mass ratio, i.e. the very limits described in items (i) and (ii) of Sec. 6.2.
Quite remarkably, we obtain the largest GW amplitudes for the smallest mass ratio q = 0.38. This mass ratio also exhibits the most interesting waveforms. For a wide range of the dephasing angle δφ, we find the GW radiation to be quadrupole dominated with a merger pulse reminiscent of BH head-on collisions; see e.g. Fig.8 of Ref. [80]. However, certain phase off-set parameters result in fainter and aberrant GW signals with signatures of tidal deformation of the binary constituents. We attribute this more complex shape of the fainter signals to the lighter BS in the binary becoming more prone to tidal effects from its more compact and heavier companion. For a given self-interation constant σ 0 and with decreasing mass, it has been shown that tidal deformability of the star increases [60], and as a result the gravitational waveform can significantly depart from the BH-like form [104]. We illustrate this phenomenon in more detail in Figure 8, which shows the q = 0.38 waveforms for phase off-set parameters δφ = 60 • , 90 • , 150 • , 210 • . The former two give the smallest GW amplitude and exhibit significant deformation, whilst the latter two show a clear, 'black hole' like signal. In the case of fainter GW signals, the importance of higher modes becomes more prominent: as indicated in the same Figure, = 3 modes become almost comparable in amplitude to the = 2 modes for binaries with δφ = 60 • , 90 • .

Energy radiated by unequal-mass binaries
Similar to the equal-mass binaries, in the unequal-mass case we also observe some discrepancy in the GW energy with varying initial separation and dephasing parameter. We recall that these effects (see Fig. 4) have been attributed to the differences in the collision velocity and the degree of the repellent nature of the scalar field. However, in the unequal-mass case the dependency of the energy on the dephasing parameter and For each mass ratio, initial separation and superposition method, we choose the phase off-set parameter that maximises the GW energy; cf. Table 5. The straight line shows the waveform obtained for our improved method, whilst the dashed line shows that for plain superposition. As we increase the separation, both superpositions give comparable results. This is most notable in the case of smaller mass ratio q = 0.38, where good agreement is already reached at smaller distances d/M 32.4.
separation becomes even more complex. This is due to the fact that the phase difference at merger for unequal-mass binaries is no longer the initial dephasing δφ we apply to one of the BSs and as we vary the initial separation d we further change the phase difference at merger. In Figure 9, we illustrate this dependence of the radiated energy on the dephasing parameter, E rad (δφ), for varying initial separations and all the mass ratios considered here. As shown in the left column, plain superposition results in flat energy profiles for smaller separations. We see here once again a manifestation of premature BH formation which largely eliminates the effect of the scalar field's dephasing on the merger dynamics. This is unlike the energy profiles of our improved superposition displayed in the right column of Fig. 9. Here the energy dependence on the dephasing parameter takes on an approximately sinusoidal shape for all separations d. The distinct horizontal shift between these profiles can be attributed to the infall-time dependent contribution to the dephasing δφ. Except for the smallest separation and modulo the horizontal shift, the energy profiles exhibit comparable maxima and minima as we change the initial separation. However, similar to the equal-mass case, the results for the smallest separation differ significantly, presumably due to differing collision velocity as illustrated in Fig. 5.
Since the dependence of the radiated energy on δφ is time-dependent in the unequalmass case, E rad (δφ) is no longer described by a single sinusoidal fit (25). In fact, we find that a two-mode sinusoidal fit well approximates the data. This two-mode fit applies to all configurations using the improved initial data construction and the data of plain superposition at larger initial separations. Only in the case of small initial separations, where plain superposition results in premature BH formation, the energy is well fitted with the one-mode fit (25); here the BH formation eliminates the overall effect of the dephasing as well as any complications arising from its time dependence during the merger stage. In summary, our results demonstrate that plain superposition not only results in quantitative changes in the emitted GW signals, but also leads to a significant over-simplification of the merger dynamics.

Conclusion
In this work, we have extended previous studies of boson-star (BS) binaries in two principle directions: (i) we have generalised the initial data construction for equal-mass BS binaries of Ref. [73] to the unequal-mass case, (ii) we have systematically explored the effect of the scalar field's dephasing δφ on the merger dynamics and GW emission. For all mass ratios and dephasing parameters, we have shown how plainly superposed initial data can result in spurious physical effects, such as increased constraint violations, oscillations of the scalar fields' central amplitudes and premature BH formation. The key drawback of plain superposition that leads to these spurious effects is its failure to recover equilibrium values of the volume element at the stars' centres. By appropriately correcting the conformal factor of the spatial metric, our method exactly recovers the equilibrium volume element at the centres of the stars and thus circumvents the spurious features of plainly superposed data for all mass ratios and δφ. Notably, our improved method explicitly incorporates the equal-mass fix and plain superposition as limiting cases. Similar to previous studies (e.g. [68,100,102]), we find that the choice of the dephasing parameter significantly affects the merger dynamics of BS binaries, most notably through the repellent character of merging scalar-field solitons with large phase differences. Crucially, for unequal-mass binaries, this phase difference is not equal to the initial dephasing parameter δφ but also acquires a time-dependent contribution due to the individual BSs' different oscillation frequencies.
To leading order, we find that the radiated energy and, hence, GW amplitude depends approximately sinusoidally on the dephasing parameter.
More specifically, we find equal-mass collisions to result in the strongest GW signals (and maximal energy) for δφ = 0, i.e. configurations with equal phase at merger. As the dephasing δφ is increased towards π, the energy and GW amplitude decrease monotonically. For unequal-mass binaries, we observe the same behaviour, accompanied, however, by a constant offset in δφ due to the differing oscillation frequencies and the subsequent additional phase offset at merger; cf. Figs. 4 and 9.
In general, the BSs' phase offset introduces considerable complexity to the merger dynamics and GW emission. In the equal-mass case, this manifests itself most prominently in a significant variation of the GW energy from off-phase BS binaries as we change the initial separation; cf. Fig. 4. We attribute this variation to differences in the binaries' binding energies and the corresponding differences in the collision velocities as shown in Fig. 5. We observe the same phenomenon for unequal-mass binaries but with maximal and minimal radiation occurring for shifted values of δφ which we ascribe to the time dependent nature of the dephasing. The extent of the energy discrepancy depends on the dephasing parameter. This may be connected to complex interactions between the collision velocity and the scalar fields' repulsion, which we leave for future study.
Additionally our unequal-mass collisions exhibit several distinct features which we summarize as follows.
(1) As shown in Fig. 9, smaller mass ratios produce larger GW energy than the equalmass case. This is in contrast to the BH case, where E rad ∼ η 2 . . = q 2 /(1 + q) 4 [106], where η is a monotonically decreasing function of q.
(2) As shown in Fig. 8, certain dephasing parameters result in weaker GW signals with signatures of tidal deformation. For these configurations, higher modes, such as = 3, exhibit almost comparable magnitude as their quadrupolar counterparts, especially in the case of binary with q = 0.38. This is unlike the black hole case, where = 3 mode is 5 times smaller than the = 2 mode [105].
(3) The numerical data of the radiated energy E rad (δφ) in Fig. 9 displays some deviation from the pure sinusoidal fit of (25) and is better described by the two-mode sinusoidal fit (26). We believe this is a feature of the time-dependent phase difference during the merger.
Clearly the enhanced parameter space leads to a rich structure in the merger dynamics of unequal-mass BS binaries which necessitates further systematic exploration, especially in the case of inspiralling binaries. Furthermore, the proposed initial data superposition in this work is still an approximation to the ultimate goal of fully solving the constraint equations. We leave these explorations for future efforts.

Appendix A. Choice of conformal factor
In Eq. (18), rewritten here as we have expressed the conformal factor in terms of the variable λ = γ 1/3 . The exponent 1/3 is, of course, a free choice in this equation and we will now explore the implications of choosing a different exponent and why 1/3 turns out to work particularly well in our practical applications. Let us start with the conformal factor motivated by the Schwarzschild metric in (Cartesian) isotropic coordinates, Note that the spatial metric is now written as which is Eq. (18) -for the special case of conformal flatness -written in terms of the alternative variable ψ and an exponent 4. We generalize this freedom of writing the conformal factor by introducing the variable Λ = ψ n = γ n/12 with γ . . = det γ ij , (A.4) where 7 n ∈ Z. For n = −4, for example, we recover the customary BSSN/CCZ4 conformal function Λ = χ, whilst for the choice of n = 4 we recover the conformal factor Λ = λ = γ 1/3 of Eq. (18). In the general case, we conformally decompose the spatial metric according to In our construction of BS binary initial data, we start with the spatial metric γ ij obtained from plain superposition and then conformally rescale this metric in order to correct the volume element at the individual BSs' centers. This correction gives us a new spatial metric By construction, this correction will recover the correct volume element at the centers of both BSs for any choice of n. The metric corrections thus introduced in the neighbourhood of the BS centers, however, will differ for different choices of n. And Table A1. A list of the conformal factor functions used in our discussion of binary initial-data construction. Note that in the main text we fix n = 4 in the newly defined conformal factor, so that γ ij = λ newγij .

Name
Variable Relation toγ ij Standard notation for isotropic Schwarzschild Conformal variable for n = 4 (our preferred choice) λ γ ij = λγ ij we see in Eq. (A.6) that the conformal rescaling is a linear function of our conformal variable only for n = 4. If we Taylor expand the rescaling factor Λ new /Λ around either BS center (where, we recall, it gives us the exact correction), then any choice n = 4 will result in complicated additional terms in the Taylor expansion of the factor (Λ new /Λ) 4/n . We cannot rigorously prove that this non-linearity in Eq. (A.6) inevitably leads to a significant deterioration of our initial data construction, but this is exactly what we observe in all our numerical experiments; choices n = 4 systematically result in significantly larger constraint violations compared to those of Fig. 3 and generally more so the further n deviates from 4. As a summary, we list in Table A1 the different variables for the conformal factor discussed in our initial data construction.
Appendix B. Improved superposition in the limits R A , R B → 0 or R A , R B → ∞ In Section 5.2 we have demonstrated that in the equal mass case our improved superposition (23) with R A , R B → 0 results in the same gravitational waveform as obtained with plain superposition, whilst with R A , R B → ∞ we recover the waveform from the equal mass fix of Refs. [72,73]. Here we derive analytically that our proposed method indeed reduces to plain superposition and the equal-mass fix in the respective limits.
Appendix B.1. The limit R A , R B → 0 In the limit R A , R B → 0, the weight functions of Eq. (24) become where r A = ||x i − x i A || and r B = ||x i − x i B ||. This implies, in particular, that in Eq. (22) the terms w A (x i A ) → ∞ and w B (x i B ) → ∞ diverge at the respective stars' centres, whereas w A (x i B ) and w B (x i A ) remain finite, so that . (B. 2) The correction δλ(x i ) applied to the conformal factor in Eq. (21) then becomes We thus recover δλ = 0, i.e. plain superposition, everywhere except at the isolated points x i A and x i B . The Dirac δ function like correction at x i A and x i B is a consequence of our condition (20) but is lost in numerical evolutions due to finite resolution, so that for R A , R B → 0 we expect to obtain the same results as for plain superposition.
Appendix B.2. The limit R A , R B → ∞ For R A , R B → ∞, we can write the weight functions (24) as where J = A, B and we have Taylor expanded to first order in is simply the separation d of the two BSs, we obtain, again to leading order, (B.5) We likewise expand, to leading order in (B.6) In the equal-mass case, the corrections at the BSs' centers are the same, δλ( which, after multiplication withγ ij , is the equal mass fix (16).  Table 2. The upper panel shows the differences in the energy results obtained for different resolutions, measured here in terms of the grid spacing on the innermost refinement level. The high-resolution difference has been rescaled by Q 1 = 0.8 and Q 4 = 1.57 corresponding to first and fourth-order convergence. While the small contribution due to noisy junk radiation converges only at first order, the total radiated energy exhibits convergence close to fourth order. The bottom panel displays the radiated energy as a function of time for the three resolutions as well as a fourth-order Richardson extrapolation.
We have analyzed in the same way the quadrupole of the Newman-Penrose scalar Ψ 4 and observe the same order of convergence but find a total error about twice as large, 4 % in the GW amplitude.  Figure C2. Convergence analysis of the GW energy computed with at R ex = 120 for the GRChombo simulation of BS binary q075 d12 p000 of Table 2. The difference between high and medium resolutions has been rescaled by factors Q 1 to Q 4 corresponding to convergence of first to fourth order. At early times, the order of convergence fluctuates between second and third, whilst at later times it settles to just below second order of convergence. uncertainty of 2.6%. This thus gives us a total error budget for GRChombo of 3.7%.  Figure D1. The logarithm of the L2-norm of the Hamiltonian constraint violations is shown for the binary configurations q05-d14-p000 (upper) and q038-d16-p000 (lower panels). For each case, the right panel shows a zoom-in on the region (R A , R B ) = [(0, 100) × (0, 100)]. The grey region indicates the parameter space of (R A , R B ), where constraint violations are too large; note that this region increases in size for smaller mass ratios.