Gravitational Wave Signatures of Highly Compact Boson Star Binaries

Solitonic boson stars are stable objects made of a complex scalar field with a compactness that can reach values comparable to that of neutron stars. A recent study of the collision of identical boson stars produced only non-rotating boson stars or black holes, suggesting that rotating boson stars may not form from binary mergers. Here we extend this study to include an analysis of the gravitational waves radiated during the coalescence of such a binary, which is crucial to distinguish these events from other binaries with LIGO and Virgo observations. Our studies reveal that the remnant's gravitational wave signature is mainly governed by its fundamental frequency as it settles down to a non-rotating boson star, emitting significant gravitational radiation during this post-merger state. We calculate how the waveforms and their post-merger frequencies depend on the compactness of the initial boson stars and estimate analytically the amount of energy radiated after the merger.


I. INTRODUCTION
The direct detection of gravitational waves (GWs) by LIGO has begun a new era for strong-field gravity. Four events so far [1][2][3][4] appear to represent the inspiral, merger, and ring-down of binaries composed of two black holes (BHs), and an additional one is consistent with a binary neutron star system [5] (or a low mass BHneutron star system [5,6]). More detections of compact binary mergers will surely increase our understanding of compact objects generally. While BHs and neutron stars now represent the standard model of compact objects, exploring the extent to which alternatives differ in their GW signatures remains an important test.
Among the many exotic alternative compact objects that have been proposed (e.g., fuzzballs [7], gravastars [8], wormholes [9], etc), boson stars (BSs) [10] are arguably among the better motivated and very likely the cleanest to model (see Ref. [11] for a recent review on exotic compact objects). Indeed, even if perhaps not realized in nature, they are excellent test-beds to explore possible phenomenology associated with highly compact objects.
BSs are composed of a complex scalar field with harmonic time dependence [10,12] (see [13] for a review). These systems possess a U(1) symmetry that gives rise to a conserved Noether charge. The stability properties of a BS resemble those of neutron stars, and, in particular, they are stable below a critical mass. The discovery of the Higgs boson in 2012 [14,15] shows that at least one scalar field exists in nature; if other (stable) bosonic particles exist in the universe they might clump together to form self-gravitating objects.
Motivated by the existent and future observations of compact object mergers, we study here the inspiral and merger of two BSs initially in a tight, quasi-circular or-bit. The remnant of this merger can generally be either a BS or a BH. For the merger to result in a long-lived BS, however, it is required that (i) the remnant star be stable, which in turn implies that its mass is less than the maximum mass allowed for the model, and (ii) the angular momentum left over after the merger satisfy the quantization condition for BSs [13] and hence either vanish or be larger than the minimum angular momentum for a rotating BS.
We are interested in these BS signals within the context of binary BHs and binary neutron stars, and we will therefore consider that range of compactnesses for the initial binary components. For increasing compactness, the GW signal of a BS binary is expected to more closely resemble that of a BH binary. For small compactness, the structure and tidal deformability [16,17] of the star will play a significant role, similar to the signals from binary neutron stars.
This work extends the study begun in [18] that focused on the dynamics of solitonic BS binaries. Solitonic BSs [19] are a specific family of BS with a potential (see Eq. (4) below) that yields compactness C ≡ M/R (M and R being the mass and radius of the BS, to be defined below) comparable and even higher than that of neutron stars. In particular, the maximum compactness of the stable configurations can reach values C max ≈ 0.33. As shown in [20], a Schwarzschild photon sphere -which is crucial in order to get a behavior similar to that of a BH under linear perturbations -appears outside the star at such compactness.
In particular, in [18] the constituents of these binaries were solitonic BSs with a fixed compactness C = 0.12 but with different phases and signs in their time oscillations. For example, the head-on collision of two BSs differing only by a relative phase difference generically produced arXiv:1710.09432v1 [gr-qc] 25 Oct 2017 a massive and significantly perturbed BS. However, a BS that collides head-on with another BS with the opposite frequency but otherwise identical (i.e., an anti-BS) annihilates the Noether charge of the system such that all the scalar field disperses to infinity.
Surprisingly, the addition of angular momentum does not seem to change this qualitative picture [18]. Although the merger of a pair of orbiting BSs produces a rotating remnant, the final object eventually settles down into a stationary non-rotating BS by radiating all its angular momentum via GW radiation and scalar field dispersion. Notice that similar scenarios were studied previously [21,22] in the context of mini BSs, which are much less compact than solitonic ones. However, in those cases, the long dynamical timescales of the stars prevented definite conclusions about the end state of the remnants.
Here we focus on the dynamics and GW signatures of the inspiral of binaries consisting of solitonic BSs of increasing compactness. As expected, the frequency of the GWs after the merger increases with compactness, an aspect they share with neutron stars and other compact objects. We also compare these signals with those of BH and neutron-star binaries and discuss the prospects for differentiating these objects with LIGO/Virgo observations.
This work is organized as follows. In Sec. II the structure of the initial data and evolution system are described in some detail. In particular, the Einstein-Klein-Gordon equations are presented, considering a particular potential leading to solitonic BSs. The procedure to construct initial data consisting of an isolated solitonic BS is described, followed by a summary of our initial configurations of binary BS systems. This section culminates with a short summary of the evolution formalism adopted in this work. In Sec. III we study the coalescence of binary BSs for different compactness, focusing on the dynamics and the corresponding GW signal. Finally, we present our conclusion in Sec. IV. Throughout this paper, Roman letters from the beginning of the alphabet (a, b, c, . . .) denote space-time indices ranging from 0 to 3, while letters near the middle (i, j, k, . . .) range from 1 to 3, denoting spatial indices. We also use geometric units in which G = 1 and c = 1, unless otherwise stated.

II. SETUP
BSs are self-gravitating compact objects made of a complex scalar field satisfying the Einstein-Klein-Gordon equations (for a review, see [13]). In this section we present these equations, describe the initial data, detail the evolution formalism, and finally summarize the numerical techniques used to evolve a BS binary.
A. Einstein-Klein-Gordon equations BSs are described by the Einstein-Klein-Gordon theory (1) where R is the Ricci scalar associated with the metric g ab , Φ is a minimally coupled, complex scalar field, and V |Φ| 2 is the potential for the scalar field. Variations with respect to the metric and the scalar field yield the following evolution equations where T ab is the scalar stress-energy tensor Different BS models are classified according to their scalar self-potential V |Φ| 2 . Here we focus on the socalled solitonic potential [19] where m b and σ 0 are two free parameters. The bare mass of the scalar field reads as m b whereas, in our units, m b has the dimensions of an inverse mass. After restoring physical units, the maximum mass of nonspinning BS solutions in this model reads as M max ≈ 5 10 −12 /σ 0 2 500 GeV m b M , when σ 0 1. Therefore, depending on the model parameters, solitonic BSs can be as massive as (super)massive BHs and, when σ 0 1, they can be slightly more compact than ordinary neutron stars [19].
Owing to the U(1) invariance of the action (1), BSs admit a Noether charge current satisfying the conservation law The spatial integral of the time component of this current defines the conserved Noether charge, N , which can be interpreted as the number of bosonic particles in the star [13].

B. Initial data of binary solitonic BSs
Our initial data is constructed as a superposition of two boosted, isolated, solitonic BSs. Let us first describe the solution of a spherically symmetric BS in isolation. The ansatz for the metric and the scalar field are given by g ab = diag −e Γ(r) , e Λ(r) , r 2 , r 2 sin 2 θ , where ω is a real frequency. Despite the phase oscillation of the scalar field in time, the Einstein-Klein-Gordon equations yield solutions that are otherwise static and that are described by the following set of ordinary differential equations where a prime denotes derivative with respect to r. The equilibrium configurations are found by integrating numerically Eqs. (9)-(11) along with suitable boundary conditions. Namely, we impose regularity at the center, Γ(0) = Γ c , Λ(0) = 0, φ(0) = φ c , φ (0) = 0, and at infinity we require the metric to be Minkowski and the scalar field to vanish. For a given value φ c of the scalar at the center of the star, the problem is then reduced to an eigenvalue problem for the frequency ω, which we solve through a standard shooting method. The value of Γ c can be arbitrarily chosen so that Γ(r → ∞) = 0. The ADM mass of the BS is M = m(r → ∞), where m(r) is defined by e −Λ(r) ≡ 1 − 2m(r)/r. In the solitonic model, the scalar field has a very steep profile which makes the numerical integration of the equilibrium equations very challenging, requiring very finetuned shooting parameters [23]. On the other hand, the steepness of the scalar field also alleviates a customary problem in defining the BS radius; strictly speaking, the scalar field has non-compact support, and thus BSs do not have a hard surface. Following previous work, we define the effective radius R M as the radius within which 99% of the total mass is contained, i.e. m(R M ) = 0.99M . In the model under consideration, since the scalar field is very steep, choosing a higher threshold does not affect the mass significantly. Accordingly, we define the compactness as C ≡ M/R M , so that C = 0.5 for a Schwarzschild BH and C ≈ 0.1 − 0.2 for a neutron star. In numerical simulations, it is more practical to estimate the radius of the final remnant through the radius that contains 99% of the Noether charge, R N , so we will also adopt this definition when needed. For the initial configurations considered in this work, the two definitions of the radius differ by 5−10%, more compact configurations displaying the smaller difference.
With the above shooting procedure, we compute a sequence of isolated BS solutions characterized by the central value of the scalar field φ c . In particular, here we restrict ourselves to σ 0 = 0.05 which is chosen because it allows for very compact, stable configurations. In Fig. 1, a number of relevant quantities of this family are shown. Among these, we show the compactness C as a function of φ c which achieves a maximum of roughly 0.33. The markers shown in Fig. 1 denote the four representative initial configurations investigated in this study. The relevant parameters of these solutions are given in Table I. Following [23], we can rewrite Eqs. (9)- (11) in terms of the following dimensionless quantities (12) where λ = σ 0 √ 8π. Doing so produces equations independent of m b , and hence m b serves to set the units of the physical solution. As such, we choose m b so that the BS mass is M = 0.5, and the total mass of binary systems constructed with these solutions (described below) is roughly unity for any compactness of the binary components. Note that while the scaling with m b shown in Eq. (12) is exact, the scaling with respect to λ is not. Only in the σ 0 1 limit does the scaling hold. For the value σ 0 = 0.05 considered here, the scaling is only approximate.
With the isolated BS configurations in hand, we can now discuss construction of initial data describing an orbiting BS binary. The method, described in detail in [18,24,25] is based on a superposition of two boosted isolated BS solutions, and can be summarized as follows: 1. construct the initial data for two identical, spherically symmetric BSs as discussed in the previous section. We denote the metric and scalar field as {g 2. extend the solution to Cartesian coordinates, with the center of the star located at a given position x j c , that is, {g

perform a Lorentz transformation to the solution
for each BSs along the x-direction with velocity v x , namely {g . superpose the two boosted solutions for each of the two stars: where η ab is the Minkowski metric.
Notice that this superposition in only an approximate solution that does not satisfy exactly the constraints at the initial time. However, our evolution scheme enforces the exponential decay of these constraint violation dynamically (e.g., see Figure 10 in [18]).
The positions and velocities of each binary system considered in this work, together with other parameters of our simulations, are presented in Table II.  Table I] whereas squared markers in the right panel refer to the final remnant produced by the merger of stars in an initial configuration indicated by the same color [cf. Table II]. Two squares corresponding to two configurations are not shown; the remnant of the black configuration is not well-enough resolved and did not reach a quasi-stationary state, and the green configuration produces a BH (with C > 0.5) instead of a BS. In the left panel, the horizontal line denotes M = Mmax/2, i.e. the merger of two identical BSs laying above this line is expected to produce a BH (neglecting energy dissipation during the coalescence and non-linear effects). The radius RM is defined as that containing 99% of the mass of the star, except for the radius of the remnant which is instead defined as that containing 99% of the Noether charge, RN . or RN , respectively) and angular frequency of the phase of φ in the complex plane, in dimensionless units on the left and in units such that M = 0.5 on the right. Note that high-compactness configurations require a very fine tuning in ω. Here we show only the first nine decimal figures. In the last two columns, we give the normalized Newtonian moment of inertia (where I = dmL 2 , L being the distance from the axis of rotation) and dimensionless tidal Love number (k tidal ) of the corresponding configuration as computed in [16,17]. As a reference, k tidal ≈ 200 for a neutron star with an ordinary equation of state, and k tidal = 0 for a BH.
C y x , the initial total ADM mass M0, the initial total orbital angular momentum J0 of the system, the time of contact of the two stars tc, the final remnant, the final total ADM mass Mr, the averaged final radius of the remnant star R N r (i.e., containing 99% of the total Noether charge), the frequency fr of the fundamental mode of the remnant, its dimensionless value Mrωr (where ωr = 2πfr), the total radiated energy in gravitational waves for each simulation E rad (i.e., integrated from the beginning and extrapolated to large times after the merger) and the one estimated analytically E rad as described in appendix A. The final angular momentum of the BS remnant tends to zero quite rapidly. The final (dimensionless) angular momentum of the BH obtained in the C = 0.22 case is Jr/M 2 r ≈ 0.64.

C. Evolution formalism
We adopt the covariant conformal Z4 formulation (CCZ4) (for full details see [18,26]). The Z4 formulation [27,28] extends the Einstein equations by introducing a new four-vector Z a which measures the deviation from Einstein's solutions, namely Notice the presence of damping terms, proportional to the parameter κ z , which enforce the dynamical decay of the constraint violations associated with Z a [29] when κ z > 0. To write Eq. (15) as an evolution system, one needs to split the spacetime tensors and equations into their space and time components by means of the 3 + 1 decomposition. The line element can be decomposed as (16) where α is the lapse function, β i is the shift vector, and γ ij is the induced metric on each spatial foliation, denoted by Σ t . In this foliation we can define the normal to the hypersurfaces Σ t as n a = (−α, 0) and the extrinsic curvature Via a conformal transformation, one transforms to a conformal metricγ ij with unit determinant. In terms of a real, positive conformal factor χ, one then obtains a conformal trace-less extrinsic curvatureÃ ij where trK = γ ij K ij . These definitions leads to new constraints, called conformal constraints, which can also be enforced dynamically by including additional damping terms to the evolution equations proportional to κ c > 0. For further convenience we can redefine some of the evolved quantities as where Θ ≡ −n a Z a . The final set of evolution fields is {χ,γ ij , trK,Ã ij ,Γ i , Θ}, and the equations can be found in [18]. These equations must be supplemented with gauge conditions for the lapse and shift. We use the 1+log slicing condition [30] and the Gamma-driver shift condition [31]. Our version of the CCZ4 formalism, together with the dynamical gauge conditions, is a strongly hyperbolic system, as demonstrated in [18].

D. Numerical setup and analysis
We adopt finite difference schemes based on the Method of Lines on a regular Cartesian grid. A fourth order accurate spatial discretization satisfying the summation by parts rule, together with a third order accurate (Runge-Kutta) time integrator, are used to achieve stability of the numerical implementation [32][33][34]. To ensure sufficient resolution within the BSs in an efficient manner, we employ adaptive mesh refinement (AMR) via the had computational infrastructure that provides distributed, Berger-Oliger style AMR [35,36] with full subcycling in time, together with an improved treatment of artificial boundaries [37]. We adopt a Courant parameter of λ c ≈ 0.25 such that ∆t l = λ c ∆x l on each refinement level l to guarantee that the Courant-Friedrichs-Levy condition is satisfied. Previous work with this code [18,21,22] established that the code is convergent and consistent for the evolution of BSs.
Our simulations are performed in a domain [−264, 264] 3 with a coarse resolution of ∆x 1 = 4 and either 7 or 8 levels of refinement, the last one only covering each star, such that the finest resolution is ∆x 8 = 0.03125.
We analyze several relevant physical quantities from our simulations, such as the ADM mass, the angular momentum, and the Noether charge, computed as described in [18]. To analyze the gravitational radiation emitted during the coalescence, we have considered the Newman-Penrose scalar Ψ 4 . As is customary, we expand Ψ 4 in terms of spin-weighted spherical harmonics (with spin weight s = −2) The components h l,m (t) of the strain (per unit distance) can be computed by integrating twice in time the Ψ 4 . However, this procedure introduces integration constants which might cloud its purely harmonic behavior. These constants do not appear in the frequency domain, is the Fourier transform of the l, m mode of Ψ 4 (t). This relation is valid for the physical frequencies f ≥ f 0 , where f 0 is the initial orbital frequency. A way to enforce this condition is by setting f = f 0 for f < f 0 in Eq. (23) [38], which acts as a high-pass filter and attenuates signals with frequencies lower than the cutoff frequency f 0 . The strain in the time domain can be obtained by performing the inverse The GWs, the ADM mass, and the angular momentum are computed as spherical surface integrals at different extraction radii to check the consistency of the results.

III. COALESCENCE OF BINARY BSS
In this section we analyze the dynamics of binary BSs as a function of their individual compactness. For concreteness, we consider four binary systems composed of BSs on the stable branch with compactness ranging from C ≈ 0.06 to C ≈ 0.22 (i.e., see Table I). Notice that all the stars have been rescaled such that their individual masses in isolation are M = 0.5, so that the binary has approximately a total initial mass M 0 ≈ 1.

A. Dynamics
Some snapshots of the Noether charge density and the norm of the scalar field during the binary coalescence are displayed in Figs. 2 and 3, respectively. The early part of the inspiral is qualitatively similar for all cases, with the stars completing at least an orbit before they make contact. Notice, however, that the initial orbital velocities differ among the binaries, mainly because of two reasons. The first reason is that the binaries are constructed of BSs with different masses, although the scaling is such that the binary mass is approximately unity. The second reason is that for the smaller compactness cases, the initial separation is chosen larger than that for the larger compactness to accommodate the larger size of the stars. Consequently, the GWs produced during the earlier stages of the inspiral will be weaker than the corresponding stage of the larger compactness cases.
Once the stars make contact, scalar field interactions play a significant role in the dynamics of the remnant and nonlinear effects due to the scalar potential becoming significant. The combination of gravity forces and matter interactions produces a compact, rotating remnant immediately after the merger. The final fate of this remnant will depend both on the potential and on the location of the solution within the stable branch, that is, on how far from the unstable branch this specific configuration is. This distance can be parametrized in different ways, such as the fraction M/M max (where M max is the maximum allowed mass on the stable branch) or, equivalently, on the initial stellar compactness.
Our results indicate a transition between binaries with small compactnesses with those of larger compactness. In particular, the critical compactness, C T , appears to be somewhere in the range of 0.18 − 0.22. This transition is roughly estimated by the fact that twice the individual BS mass above the critical value exceeds the maximum stable BS mass, suggesting that any remnant (less any radiated scalar field) would be unstable. We can distinguish two different behaviors depending on the initial compactness, separated by this transition value: • C C T : the remnant is a significantly perturbed BS, the angular momentum of which decreases through dispersion of scalar field and gravitational radiation, settling down into a non-rotating BS.
Furthermore, for the case with C = 0.12, as already observed in Ref. [18], the angular momentum of the remnant is further reduced through the ejection of two "rather cohesive" scalar field blobs soon after the merger (see third panel, second row in Fig. 2).
• C C T : the remnant mass exceeds the maximum mass and promptly collapses to a BH with approximately the mass and angular momentum of the system at the merger time. There is some scalar field surrounding the BH that carries angular momentum and is being either slowly dispersed to larger distances or falling into the BH.
The evolution of the ADM mass, angular momentum, and Noether charge are illustrated in Fig. 4. The binaries show only a significant loss of ADM mass near the merger due to scalar field dispersion/ejection and energy carried away by GWs. Similar behavior is reflected in the total Noether charge -when the remnant does not collapse to a BH -. Since the Noether charge is a conserved quantity, the fact that it remains mostly constant further supports the reliability of our simulations (see also Section II D for further discussion of numerical tests).
Notice however that the case C = 0.18 shows peculiar behavior in its mass and angular momentum after merger when compared to the other cases. A close inspection of the dynamics of this case shows that large gradients develop that are not accurately tracked by the finest resolution we have allowed our adaptive grid to achieve. Indeed, it appears the system explores a nearthreshold regime which is not correctly captured by our simulations and we consider the post-merger period of this case to be unreliable. Nevertheless, for reference and comparison purposes we include it in the overall discussions since, in any case, its pre-merger behavior is informative.
Angular momentum during the early inspiral is radiated mainly through GWs. Near the merger stage there is dispersion (and in some cases, also ejection) of scalar field, which also carries away a significant fraction of the angular momentum (notice the sharp decrease in the middle panel of Fig. 4). After the merger, the final object is much more compact than the initial stars and rotates rapidly, emitting GWs more copiously than the late inspiral. When the remnant is not a BH, the system radiates angular momentum until settling down to a non-rotating BS.
As discussed in [18], that the system approaches a nonrotating BS might be a consequence of two combined effects: (i) the quantization of the angular momentum J z = kN (where k an integer) of rotating BSs that prevents stationary solutions with an arbitrary angular momentum, and (ii) the rigid structure of the rotation boson star may present difficulties for the scalar field to organize itself into a stable, rotating configuration. In particular, the rotating boson star is harmonic both in time and azimuth with the level sets of its magnitude being toroidal. This structure contrasts with a rotating neutron star which can have a range of rotational profiles either rigid or differential.
Also, another possibility is that the first rotating configuration (k = 1) is unstable in the high-compactness regime explored here (see Ref. [39] for a discussion of the stability of less compact rotating BSs). As shown in Fig. 1, the mass and the compactness are very steep functions of the central scalar field, and it is therefore possible that spinning solutions exceed the maximum mass when their non-spinning counterparts are close to such a maximum. If this is the case, the unstable spinning solution would not be formed dynamically and we expect the outcome to be a non-spinning BS or a BH, in agreement with the results of our simulations.
Of particular interest is the case C = 0.12 due to the presence of two blobs of scalar field that are ejected from the remnant at speed v ≈ 0.6c. The formation of two peaks in the Noether charge density during the merger is common for all studied cases, but only for this compactness are they able to detach from the star while maintaining -at least temporarily -their character (i.e., its Noether charge). As can be surmised from Fig. 4, these blobs carry away little mass (M blob = m b N blob ≈ 0.025) but a significant fraction of angular momentum. A Newtonian estimate, assuming the distance from the blobs to the plane of symmetry is L ≈ 7.5, yields J z ≈ 2M blob vL ≈ 0.2, consistent with the additional decrease of angular momentum displayed in the C = 0.12 case with respect to the C = 0.06 configuration, as shown in the middle panel of Fig. 4.
We can gain some insight into the ejection of these blobs by examining some characteristic speeds in the problem. A simple calculation shows that the Newtonian angular orbital frequency Ω c when the two identical stars first make contact is where C is the compactness of the individual stars and M 0 is the total initial mass of the system. Notice that the blob velocity (0.6c) is considerably larger than the maximum velocity v c = Ω c (2R) ≈ 0.35c predicted by this rotational rate for solid body rotation. Similarly one can compute the velocity associated with rotation at the angular frequency of the remnant 1 . By using the values of table II for the case C = 0.12 one can obtain that the orbital frequency of the remnant is Ω r = πf r = 0.098 and its radius is R r = 4.6. The velocity of the remnant, v r = Ω r R r ≈ 0.45c, also does not reach the level of the blob velocity. However, these blobs are ejected during the time when the binary is transitioning from first contact to quasinormal ringing. As such, the characteristic angular frequency and radius change Ω c → Ω r and 2R → R r . At some point during this transition the frequency and the radius might be large enough so that some of the scalar field might move with a speed larger than the escape velocity v esc = √ 2C of the star. If this occurs, then it is conceivable that some amount of scalar field may be ejected from the remnant at such a speed.

B. GW signal
The merging binaries produce GWs measured by the Newman-Penrose Ψ 4 scalar, as displayed in Fig. 5. In Fig. 6 we also show the corresponding strain. Notice that the amplitude and the time scale of the strain has been rescaled with the total initial mass, and the time has been shifted such that the contact time occurs at t = 0.
With these waveforms, we can look for the effect of compactness on the gravitational wave signal. Starting with the least compact case (C = 0.06), it radiates the least in the inspiral. The weakness of its inspiral signal results because its stellar constituents have the largest radii and thus they make contact at the smallest frequency of these four cases. As the compactness increases, the late inspiral occurs at higher frequencies and the signal becomes stronger.
Once the stars make contact, both scalar field interactions and gravitational forces determine the dynamics and tend to homogenize the scalar field profile. This period is very dynamical producing a rapidly rotating compact remnant that radiates strongly in gravitational radiation with an amplitude and frequency much larger than the inspiral. This contrast between the pre-and post-merger signals is particularly marked for the two low compactness binaries, but becomes less so with increasing compactness. This trend indicates that this contrast likely results from the disparity between the initial compactness of the boson stars and the compactness of the remnant (see Table II). Notice also that the strain amplitude (Fig. 6) does not show such disparate scales as the Newman-Penrose Ψ 4 scalar (Fig. 5) due to the additional frequency dependence (see Eq. 23). A simple estimate of the maximum amount of total energy the system can radiate follows from a model of energy balance presented in [41], and discussed in detail in Appendix A. Within some approximations, when the final object is a non-rotating BS, and the total radiated energy in GWs is estimated to be where M and C are the mass of the system and compactness of the initial stars respectively. This estimate is largely consistent (i.e., within a factor of two) with the results of our simulations, given in Table II, obtained by integrating the gravitational wave luminosity displayed in the bottom panel of Fig. 5. Notice the energy emitted in gravitational waves for the case with C = 0.12 exceeds the 5% of the total mass M 0 emitted during the analogous coalescence of a binary BH system; thus, boson star binaries in a suitable range can be considered super-emitters in the terminology of [41].
The most compact cases considered here are also interesting in the context of the recent observations of GWs by the LIGO detectors. A simple Newtonian calculation shows that the GW frequency at the contact of the two The time has been rescaled by the initial total ADM mass M0 and shifted such that t = 0 corresponds to the maximum of the norm of the mode. The amplitude has been also rescaled with the mass of the system. We have chosen the same range in the axes to make clear the increase in frequency as the stellar compactness also increases. The different cases are compared with a recent version of the EOB approximation of a quasi-circular binary BH coalescence [40] by matching the waveforms at the early inspiral. For the highest compactness C = 0.22 we have also matched to the EOB waveform at the merger time (dotted red curve).
stars is f c Ω c /π. This relation can be contrasted to the frequency at which the analogous case of binary BHs would make a transition from inspiral to plunge. This frequency is well approximated by the ISCO frequency of the resulting BH produced through the merger [42]. For non-spinning binary BHs, a handy expression for the ISCO frequency is provided in [41], which indicates that non-rotating BSs would have a contact frequency higher than the corresponding ISCO frequency for binary BHs provided C 0.27 (i.e., a compactness higher than any of the ones considered in this work). The gravitational waveforms for the binaries considered here display the expected "point-particle" behavior at the low frequencies but evidence a sharp transition around contact.
For cases not collapsing to a BH the gravitational waveforms have a rather simple structure with principal modes which can be tied to quasi-normal modes of boson stars (see appendices C and B). For cases collapsing to a BH, the post merger gravitational wave signatures are captured well by the familiar ring-down behavior of a BH. We however note an argument that the angular momentum of the remnant BH is slightly less than that of the analogous binary BH merger. Because the BH pair merges well within the system's ISCO, much of the angular momentum is trapped within the remnant. In contrast, the BS binary, being less compact and merging at a lower frequency than the BH binary, allows for the radiation of more angular momentum during the merger. Recall also that tidal effects introduce modulations in the (late) inspiral waveforms (e.g. [17,43]) but such modulations become smaller for higher compactness (see the tidal Love numbers in Table I).
To examine in more detail the after-merger behavior, we analyze the strain in the frequency domain. The Fourier transform of h 22 (t > t merger ) is shown in Fig. 7, where t merger is defined as the time where the strength of the GW is maximum. For C < C T the final remnant rotates and oscillates while settling down to a (nonrotating) stationary configuration, producing GWs at certain frequencies. Clearly, the frequency of the main peak increases with the compactness of initial objects in the binary. These post-merger frequencies, reported in Table II, can be fit in terms of the contact frequencies, namely M r ω r = 0.064 + 1.72 M 0 ω c (26) which is valid even when the remnant is not a BH. On the other hand, we can compare with the fit obtained for neutron stars [44,45], that in the same units reads These relations (i.e., with M = 2.7M in Eq. 27), together with the observed frequencies, are displayed in Fig. 8. The best fit lines have quite different slopes and intercepts, but for high compactness stars they produce similar frequencies. We expect that the difference in these frequencies implies that remnant BSs and NSs are potentially distinguishable with GW spectroscopy if either (i) a large enough SNR is achieved by increasingly sensitive detectors or (ii) a sufficient number of events can be combined (i.e., stacked) [46,47]. For comparison, we include the neutron star cases studied in [44]. The case C = 0.18 is included for reference as an unfilled square.
Quite interesting is the comparison of the main gravitational wave mode (i.e., l = m = 2) with the quasinormal modes of single isolated stars, displayed in Fig. 9 and discussed in more detail in Appendix C. Clearly, the frequencies of the remnant agree very well with the frequencies of the fundamental quasi-normal mode of single non-rotating boson stars, providing further evidence that these cases produce non-rotating, remnant BSs.
This agreement contrasts with the mergers of neutron stars for which the strongest radiating mode is produced by a rotating quadrupole distribution. Here, the BS binary remnants are similar to those of binary BHs in that the remnant rings down in accordance with QNM dominated by the m = 2 mode.
For C > C T the final remnant is a rotating BH, and its post-merger signal shows the characteristic ring-down signal. For such a case, a significant amount of energy corresponding to the rotational energy of the BH is retained after merger, which contrasts with the cases producing a remnant non-rotating boson star. For the BH case we can calculate both the frequency and the decay rate of the gravitational wave signature, which can be obtained by fitting the post-merger strain signal to We find that σ = 0.049±0.003 and f = 0.056±0.002. The final mass and angular momentum of the BH, calculated asymptotically at a spherical surface of radii R = 50, are M r ≈ 1.42 and J r ≈ 1.3 respectively. Therefore, one can also calculate from linear theory the quasi-normalmode frequencies for a BH with final dimensionless spin a r = J r /M 2 r ≈ 0.64. For this spin, a perturbative calculation yields M r ω QNM = 0.50819 − 0.082748i [48,49] Comparison between the frequencies of the fundamental quasi-normal mode of single BSs in isolation (circles) and the gravitational frequencies of the merger remnant (squares), as a function of the compactness C ≡ M/RN (the case C = 0.18 is included for reference as an unfilled square) where RN is the radius containing 99% of the Noether charge. The good agreement between these frequencies suggests that the remnant is indeed a perturbed non-rotating BS ringing down to a quiescent one.
very good agreement with the value obtained from the fit M r (2πf − iσ) ≈ 0.5 − 0.07i. The small deviations might be due either to inaccuracies in the extraction of the final BH mass or to a sub-leading effect of the scalar field during the post-merger.
Notice that at the end of our simulation there is still some rotating and dispersing scalar field around the BH, but the system does not seem to relax to the hairy BH solution found in [50], at least in the timescales considered here. This is rather natural, if one recalls the quantized tidal locking relation ω = mΩ H , which is necessary for the existence of such solutions (with m = ±1, ±2, ...). This condition relates the internal frequency of the BS, ω, with the angular velocity of the BH horizon, Ω H ≡ J r /(2M 2 r R H ) (where R H is the radius of the horizon). In our simulations Ω H ≈ 0.087 and ω ∈ [0.4 − 1.8], and thus m should be between 4 − 20. Since the equal mass binary systems drive primarily the m = 2 mode, it is difficult to fulfill the tidal locking condition unless higher modes are nonlinearly generated. Though, these would be significantly sub-leading in strength.

IV. CONCLUSIONS
In this work we have studied in detail the dynamics of binary BSs for different initial compactness in the range [0.06, 0.22]. Our analysis reveals several interesting features. We confirm and extend the results obtained in [18], namely that the merger of two solitonic BSs leads either to a BH or to a non-rotating BS. The latter case is due to the inherent difficulty to achieving, simultaneously, the quantized angular momentum condition, the specific profile for the scalar field required in rotating boson stars, and stability of highly-compact rotating BS.
Gravitational waves emitted during the coalescence of these objects show that for low-compactness BSs (C < C T ), the maximum strength achievable in the inspiral phase is rather weak but it rises rapidly during merger with a significant amount of energy during that phase. For high-compactness BSs (C > C T ), a more monotonic transition -as judged in the rate of upward frequency sweep-is observed in between inspiral and merger phases. The final object promptly collapses to a BH and the post-merger gravitational wave is dominated by the typical ring-down of a spinning BH.
For the less compact cases with C < C T , the main mode in the post-merger gravitational radiation is given by the fundamental quasi-normal frequency of the nonrotating boson star. This is in contrast with the behavior manifested in binary neutron stars (that do not collapse promptly to a BH), where the main mode is linked to the rotation of the newly formed (hyper) massive neutron star. Nevertheless, in both cases, the main mode can be linked to the contact frequency by a rather simple linear relation. Importantly for efforts to try and distinguish binary boson stars from binary neutron stars with gravitational waves, the relation is sufficiently distinct to be probed by 3rd-generation detectors and/or the combination of multiple events in aLIGO/VIRGO. Moreover, it is clear that the overall gravitational wave signals from solitonic boson stars are unlikely degenerate with those from either binary neutron stars or binary BHs except possibly in the highest compactness cases (see related discussion in [51]). That is, if the primary and secondary objects are not compact enough, the late inspiral phase will differ significantly from that of two BHs and the remnant will have peak gravitational modes at significantly higher frequencies in general than those expected in binary neutron star mergers (which do not collapse to a BH promptly). Thus, sensitive searches could in principle distinguish rather easily these two systems. For compactnesses of the order of C T 0.2 a BH is formed but one with less angular momentum than a corresponding BH binary would have formed. As such, with a sufficiently sensitive detection of the post-merger spin, one may be able to differentiate a compact BS binary from a BH binary by comparing to the inspiral angular momentum. Notice however that there could still be a degeneracy with a sufficiently compact binary neutron star merger, the presence of an electromagnetic counterpart would favor the merger involving at least one neutron star as opposed to binary boson stars [52].
We also notice in passing that in all cases considered in this work (and likely in all cases allowed by the solitonic model) the final BS remnant is not compact enough to produce echoes in the post-merger ring-down phase at late time [53], since this effect requires the exotic compact object to have an external, unstable, photon sphere, which is not the case for the BS models discussed here.
We expect our results to be qualitatively similar in other scenarios including self-gravitating, compact scalar configurations. For instance, if the scalar field is real, action (1) admits compact, self-gravitating, oscillating solutions known as oscillatons [54]. These solutions are metastable, but their decay time scale can significant exceed the age of the universe and we expect their properties to be very similar to those of BSs.
Although not likely to have astrophysical consequences, it would nevertheless be interesting to investigate the transition at C = C T which separates markedly different behaviors in the post-merger phase.
ior observed in the merger of binary BSs. Namely, we observe that, in the case where collapse to a BH is avoided, the final result is a BS remnant with no angular momentum and with radius R r having mass M r that is roughly the total initial mass M 0 (i.e., M r ≈ M 0 = M 1 + M 2 ) . The merger takes place at "contact," that is, when the stars are separated by R 1 + R 2 . Thus, the energy of the system at such an instant is roughly where we have included the binding energy of each star and have assumed that the stars have constant density. The moment of inertia I c with respect to the center of mass, assuming the stars are irrotational, can be written at contact time as Notice that at contact time the orbital frequency can be estimated as Ω 2 c = M1+M2 (R1+R2) 3 . Following the discussion in [41], for the equal mass case (M 1 = M 2 = M and R 1 = R 2 = R ) this energy can be expressed as a function of the compactness C = C 1 = C 2 as As the collapse takes place, the system ultimately settles into a non-rotating BS. The energy left in the system (beyond the rest mass) is given by the potential energy of the BS which, assuming a spherical object of uniform density, can be estimated as where we have considered an upper bound M r ≈ 2M . Assuming no scalar radiation, we can now estimate the radiated energy in gravitational waves during different states of the system. In particular, the total amount of energy radiated E rad and radiated after contact E ac rad are, and and we have estimated the ratio R/R r ≈ 0.9 from our simulations (alternatively, assuming the effective density of the individual BSs is similar to that of the merged BS, one has R/R r = 1/2 1/3 ≈ 0.8). Notice that for C 0.1, if no collapse to a black hole takes place, this implies that highly compact binary boson star systems are examples of "super-emitters," in the terminology introduced in [41] (as the analogous BH binary system emits ≈ 5% of their total mass).
Appendix B: Estimate of after-merger frequency of gravitational waves We can also attempt to estimate the frequency of gravitational waves soon after merger has taken place. For this, we begin assuming angular momentum is nearly conserved around the moment where the collision (contact) takes place. For the case of irrotational binaries, the angular momentum at contact can be approximated by the expression with I c and Ω c as defined in Appendix A. Soon after contact, we have instead where I r is the moment of inertia of the newly formed object (i.e, the remnant) assuming no prompt collapse to a BH takes place. Let us now assume a relation of these angular momentum given by L r = κL c with κ ∈ [0, 1] a factor introduced to account for loss of angular momentum during the merger. Now, since the angular frequency of gravitational waves ω = 2Ω, we have (specializing to the equal mass case) with C denoting the compactness of the individual stars, M their individual masses, and M r and I r the mass of the remnant and its moment of inertia respectively soon after contact. Rearranging, the gravitational wave frequency soon after contact implies (where we have definedĪ r = I r M −3 r ). In the case of binary neutron star mergers, extensive studies already indicate a small amount of angular momentum is radiated during the merger, so we can adopt κ 1. For realistic equations of state, we can make use of thorough investigations of the values ofĪ r for isolated stars for a wide range of compactness and mass (e.g. [55]) to evaluate Eq. (B4). This gives ω r ≈ 2.8 ω c , in excellent agreement with results from numerical simulations (e.g. [45]). This is not surprising since the angular momentum in the system right before contact is primarily transferred to the object formed after merger. This is not the case for binary boson stars. As opposed to the case for neutron stars, general values ofĪ r for boson stars are not available but we can estimate them either from our isolated BS solutions or by considering that they behave as constant density spheres, namely I = (2/5)M R 2 . Since the result will not change significantly, we use the latter approximation which yields where again we estimated R/R f ≈ 0.9 from our simulations. The corresponding estimate would give an upper bound factor (taking κ = 1) in [2 − 2.5] which is higher than the fit in Eq. (26). However, as we have seen, a large amount of angular momentum is lost through the merger. From Fig. 4 one would expect κ ≈ 1/4 making the "expected" frequencies from this naive estimate much too low when compared with the measured peak frequency in gravitational waves. This is in strong contrast to what is observed in the case of binary neutron star mergers.
As the next section illustrates, the after-merger radiation of boson stars is determined by the quasi-normal modes of the produced boson star or by the BH in the case of collapse. This is a consequence of our observation that the merger of boson stars does not produce a rotating boson star and that the speed of propagation of perturbations in boson stars is faster than that in neutron stars.
Appendix C: Quasi-normal modes of isolated solitonic boson star Quasi-normal modes for isolated solitonic boson stars can be computed by evolving numerically a perturbed star and analyzing the gravitational wave radiation. The formalism and numerical schemes are the same as the ones used in this work for binaries, such that only the initial data differs. We have chosen stable boson stars with total mass M = 1 for different compactnesses ranging from C = 0.06 − 0.22. These equilibrium configurations are deformed by adding a small perturbation on the conformal factor which introduces constraint violations below the truncation error of the unperturbed configuration, and so we need not re-solve the constraints. In order to ensure the excitation of gravitational modes, the perturbation has a toroidal shape with a m = 2 dependence in the axial direction.
The top panel of Fig. 10 shows the Fourier transform of the main gravitational wave mode (i.e., l = m = 2) as a function of frequency. Although there are several peaks in the spectra, we focus only on the two strongest modes at the lowest frequencies. Clearly, the frequencies of the fundamental and the secondary quasi-normal modes increase with compactness. The bottom panel displays the adimensional frequency of these two modes as a function 0.02 0.04 0.06 0.08 0.10 0.12 0.14 The circles and the diamonds correspond,respectively, to the frequencies ω 0 QNM and ω 1 QNM of the lowest quasi-normal modes (i.e., fundamental and secondary peaks), as a function of the compactness CN ≡ M/RN . Squares represent the gravitational frequencies of the remnant resulting from a binary merger with that initial compactness. Notice that we have included the case with C = 0.18, which we do not trust completely, and the case C = 0.22 that ends up in Kerr BH. For comparison purposes, we have included also the QNM of a Schwarzschild BH (triangle). of the compactness. We have also included the frequencies of the remnant after the merger of the case studied here, given in Table II. There is a good agreement between the QNM of the single stars and the fundamental mode of the remnant of the binary.