On the final fate of compact boson star mergers

TL;DR

This paper investigates the outcomes of binary mergers of compact boson stars using numerical simulations, revealing that head-on collisions produce a boson star, while boson-antiboson mergers lead to annihilation, and orbiting mergers form transient rotating bars.

Contribution

It introduces a detailed numerical study of boson star mergers with various parameters using a modified CCZ4 formalism, providing new insights into their final states.

Findings

Head-on collisions result in a boson star.

Boson-antiboson mergers lead to annihilation.

Orbiting mergers form a rotating bar that relaxes to a non-rotating boson star.

Abstract

Boson stars, self-gravitating objects made of a complex scalar field, have been proposed as simple models for very different scenarios, ranging from galaxy dark matter to black hole mimickers. Here we focus on a very compact type of boson stars to study binary mergers by varying different parameters, namely the phase shift, the direction of rotation and the angular momentum. Our aim is to investigate the properties of the object resulting from the merger in these different scenarios by means of numerical evolutions. These simulations, performed by using a modification of the covariant conformal Z4 (CCZ4) formalism of the Einstein Equations that does not require the algebraic enforcing of any constraint, indicate that the final state after a head-on collision of low mass boson stars is another boson star. However, almost complete annihilation of the stars occurs during the merger of a boson-antiboson pair. The merger of orbiting boson stars form a rotating bar that quickly relaxes to a non-rotating boson star.

Figures

Equations

$$R_{ab}$$

$$ds^{2}=-\alpha^{2}\,dt^{2}+\gamma_{ij}\bigl{(}dx^{i}+\beta^{i}dt\bigr{)}\bigl{(}dx^{j}+\beta^{j}dt\bigr{)},$$

$${\tilde{\gamma }}_{ij}$$

$${\tilde{A}}_{ij}$$

$$\tilde{\gamma }=1,tr\tilde{A}\equiv {\tilde{\gamma }}^{ij}{\tilde{A}}_{ij}=0,$$

$$tr\hat{K}$$

$${\hat{\Gamma }}^i$$

$${\partial }_t{\tilde{\gamma }}_{ij}$$

$${\partial }_t{\tilde{A}}_{ij}$$

$${\partial }_t\chi$$

$${\partial }_{t}tr\hat{K}$$

$${\partial }_t\Theta$$

$${\partial }_t{\hat{\Gamma }}^i$$

$${{}^{(3)}R}_{ij}$$

$${\hat{R}}_{ij}^{\chi }$$

$${{}^{(3)}\hat{R}}_{ij}$$

$${\nabla }_i{\nabla }_j\alpha$$

$$\tau$$

$$S_i$$

$${\tilde{S}}_i$$

$${\tilde{S}}_{ij}$$

$${\partial }_t\alpha$$

$${\partial }_t{\beta }^i$$

$${\partial }_t{B}^i$$

$$d{s}^2=-{\alpha }_0^{2}d{t}^2+{\tilde{\gamma }}_{ij}^0(d{x}^i+{\beta }_0^{i}dt)(d{x}^j+{\beta }_0^{j}dt),$$

$$g_{ab}-g_{ab}^0=e^{i{\omega }_k{x}^k}{\bar{g}}_{ab}(\omega ,t),$$

$$K_{ab}=i\omega e^{i{\omega }_k{x}^k}{\bar{K}}_{ab}(\omega ,t),$$

$${\partial }_t\bar{u}=-i\omega (\mathcal{A}-{\bar{\beta }}_0^s\mathrm{I})\bar{u},$$

$${\partial }_t\left(\begin{array}{c}\tilde{\gamma }\\ tr\tilde{A}\end{array}\right)\approx \alpha \left(\begin{array}{cc}0& 2({\lambda }_0-1)\\ 0& 0\end{array}\right)\left(\begin{array}{c}\tilde{\gamma }\\ tr\tilde{A}\end{array}\right)+\dots ,$$

$$g^{ab}{\nabla }^a{\nabla }_b\varphi =\frac{dV}{d|\varphi {|}^2}\varphi =V^{\prime }\varphi ,$$

$${\partial }_t\varphi$$

$${\partial }_t\Pi$$

$$T_{ab}$$

$$N\equiv {\int }_{{\Sigma }_t}(-n_a{J}^a)\sqrt{\gamma }{d}^{3}x,$$

$$J^a=i{g}^{ab}({\varphi }^{\ast }{\nabla }_b\varphi -\varphi {\nabla }_b{\varphi }^{\ast }).$$

$$d{s}^2=-H(x-t)d{t}^2+H(x-t)d{x}^2+d{y}^2+d{z}^2,$$

$$d{s}^2=-{\alpha }^2(r)d{t}^2+{\psi }^4(r)(d{r}^2+r^{2}d{\Omega }^2),$$

Full text

I Introduction

One of the most prominent opportunities in the rising era of gravitational wave (GW) astronomy is to study the strong-gravity regime through the signals produced during the coalescence of compact objects. Very recently, LIGO detectors observed the first GW signals consistent with binary black hole mergers GW150914 ; GW151226 , although other sources could produce comparable waveforms 2016PhRvL.116q1101C ; 2017arXiv170101116C ; car . Consequently, these observations are already setting bounds both on the nature of these objects and on alternative theories of gravity Yunes:2016jcc ; TheLIGOScientific:2016src .The next most likely candidate to be detected by GWs is the coalescence of binary neutron stars. Neutron stars are very compact objects, usually possessing strong magnetic field, which existence has been confirmed with the observations of several binary pulsar systems on the electromagnetic (EM) band for several decades already (see for instance doi:10.1146/annurev-nucl-102711-095018 ).

In addition to these known sources, gravitational waves can allow us to find unexpected astrophysical compact objects with low brightness, known generically as Exotic Compact Objects (ECOs). One of the most plausible ECO candidates are the boson stars (BSs), self-gravitating objects made of complex scalar field liebpa ; mace . Even if their existence is still under debate due to the lack of any observational evidence, boson stars provide a simple and useful model to study compact bodies in very different scenarios, ranging from dark matter candidates to black hole mimickers. The maximum compactness $GM/c^{2}R$ of the boson star depends on the self-interaction terms of its associated potential, and ranges from $\mathcal{O}(10^{-3})$ for mini-boson stars to $\mathcal{O}(10^{-1})$ for non-topological solitonic BS frie ; 1992PhR…221..251L .

Despite the simplicity of these smooth solutions, there are only few studies on binary boson star collisions within General Relativity. Preliminary head-on collisions of mini-boson stars were first studied in 1999PhDT……..44B within a 3D code. The dynamics of the merger, which showed an interesting interference pattern, was further analyzed in 2004PhDT…….230L with and axisymmetric code. Ultra-relativistic collisions were considered in 2010PhRvL.104k1101C , and head-on and orbital mergers of non-identical boson stars in pale1 ; pale2 . Other related works include the study of the orbital case within the conformally flat approximation instead of full GR 2010PhDT…….390M , and head-on collision of oscillatons brito , a solution analogous to BS but using just a real scalar field. Much more recently, collisions of solitonic boson stars has been numerically performed car , leading to dynamics qualitatively similar to the observed for mini-boson stars.

This work aims to extend these recent numerical studies of solitonic boson stars by exploring a wider parameter space on the initial head-on configurations. We also consider binaries of identical solitonic boson stars with angular momentum to study whether the formation of a rotating boson star as a final state is possible. Our results indicate that, at least for the low mass boson stars considered here111The merger of high mass boson stars leads to the formation of a black hole car , the head-on collision generically produce another boson star if the phase shift is not too close to $\pi$. However, the orbiting binary will not lead to a rotating boson star, but to a non-rotating one. This is probably due to the quantization of the angular momentum, and implies that the massive boson star formed by the merger must shed its excess of angular momentum by emitting GW and, in some cases, very rapid blobs of scalar field.

Our simulations will be performed by using a novel modification of the CCZ4 formalism alic that treats all the constraints in the same manner. A common feature of the current conformal formulations, like the different flavors of BSSN1995PhRvD..52.5428S ; BSSN and CCZ4 bernu ; alic , is that a subset of the constraints of the system must be enforced after each time step of the simulation in order to obtain a stable evolution. Although this feature does not present a problem when using explicit time integrators, it might be not so straightforward for more sophisticated numerical methods or for automatically generated codes Arbona20132321 . We introduce new terms in some of the equations to ensure that the full system is strongly hyperbolic (and well posed) and that the constraints are dynamically enforced during the evolution. Hence, our modified CCZ4 formalism does not require the algebraic enforcing of any constraint.

The paper is organized as follows. In section II the Einstein-Klein-Gordon (EKG) evolution system is described in some detail. In particular, after a short summary on the CCZ4 formulation, we introduce novel modifications to avoid the algebraic enforcing of the constraints. In section III the implementation of the evolution equations is briefly discussed, together with several numerical spacetimes –robust stability test, gauge waves and single solitonic boson star– to test our evolution system. We study the dynamics of binary boson stars in Section IV by analyzing our numerical simulations of head-on and orbiting cases. Finally, we present our conclusion in section V. 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=c=1$, unless otherwise stated.

II Evolution system

The interaction between scalar-field matter and gravity, required to study the dynamics of boson stars, is given by the Einstein-Klein-Gordon equations. We adopt the CCZ4 formalism of the Einstein equations, which is briefly summarized next. We stress the modifications with respect to previous works and perform a characteristic analysis of the resulting system. The evolution equations for the complex scalar field are also described.

II.1 CCZ4 formalism

The Z4 formalism was first proposed as a covariant extension of Einstein equations to achieve an hyperbolic evolution system free of elliptic constraints Z44 ; Z45 . The equations of motion, which might also be derived from a Palatini-type variation Z41 , are

[TABLE]

where $R_{ab}$ is the Ricci tensor associated to the spacetime metric $g_{ab}$, $T_{ab}$ is the stress-energy tensor (with trace $trT\equiv g^{ab}T_{ab}$) and $Z_{a}$ is a new four-vector which measures the deviation from Einstein’s solutions. Although the original formulation, corresponding to the choice $\kappa_{z}=0$, was completely covariant, additional damping terms were included to enforce a dynamical decay of the constraint violations associated to $Z_{a}$. As it is shown in gundlach , all the physical constraint modes are exponentially damped if $\kappa_{z}>0$. However, since the damping terms are proportional to the unit normal of the time slicing $n_{a}$, the full covariance of the system is lost due to the presence of a privileged time vector.

A conformal and covariant version of the $Z4$ (CCZ4) can be obtained from the 3+1 decomposition of the evolution equations by using conformal variables alic (i.e., see also bernu for other conformal but non-covariant Z4 formulations). Since this formulation is the starting point to our modifications, we first briefly summarize the derivation of the equations.

The first step involves writing the line element by using the 3+1 decomposition, namely

[TABLE]

where $\alpha$ is the lapse function, $\beta^{i}$ is the shift vector and $\gamma_{ij}$ is the induced metric in each spatial foliation. In this foliation geometry we can define the normal to the hypersurfaces $n_{a}=(-\alpha,0)$ and the extrinsic curvature $K_{ij}\equiv-\frac{1}{2}\mathcal{L}_{n}\gamma_{ij}$, where $\mathcal{L}_{n}$ is the Lie derivative along $n^{a}$. Therefore, the Z4 formalism given by eq. (1), together with the metric decomposition eq. (2) and these definitions, lead to evolution equations for the evolved fields $\{\gamma_{ij},K_{ij},Z_{i},\Theta\}$, where we have defined $\Theta\equiv-n_{a}Z^{a}$.

In the second step a conformal decomposition is applied to the evolved fields. A conformal metric ${\tilde{\gamma}}_{ij}$ with unit determinant and a conformal trace-less extrinsic curvature ${\tilde{A}}_{ij}$ can be defined as

[TABLE]

where $trK\equiv\gamma^{ij}K_{ij}$. These definitions lead to the following new constraints

[TABLE]

which will be denoted as conformal constraints from now on to distinguish them from the physical constraints associated to $Z_{a}$. Notice that now the evolved fields are $\{\chi,{\tilde{\gamma}}_{ij},trK,{\tilde{A}}_{ij},Z_{i},\Theta\}$.

Instead of using $trK$ and $Z_{i}$, it is more convenient to use the following quantities

[TABLE]

so that the evolution equations are closer to those in the BSSN formulation 1995PhRvD..52.5428S ; BSSN , where the quantity ${\tilde{\Gamma}}^{i}={\tilde{\gamma}}^{jk}\,{\tilde{\Gamma}}^{i}{}_{jk}=-\partial_{j}{\tilde{\gamma}}^{ij}$ is directly evolved. Therefore, the final list of evolved fields become $\{\chi,{\tilde{\gamma}}_{ij},{tr\hat{K}},{\tilde{A}}_{ij},{\hat{\Gamma}}^{i},\Theta\}$, following the evolution equations

[TABLE]

where the expression $[\ldots]^{\rm TF}$ indicates the trace-free part with respect to the metric $\tilde{\gamma}_{ij}$. The non-trivial terms inside this expression can be written as

[TABLE]

The matter terms are computed by contracting the stress-energy tensor, namely

[TABLE]

The evolution equations (8-13) are equivalent to those obtained in alic , by defining the conformal factor $\chi=\gamma^{-1/3}$ instead of $\phi=\gamma^{-1/6}$, except by two significant differences. First, there is a new term proportional to $tr\tilde{A}$. This term, which was already suggested in alic , is crucial to obtain a well-posed evolution system if the algebraic conformal constraints $\ln{\tilde{\gamma}}=tr{\tilde{A}}=0$ are not enforced during the evolution. Second, damping terms proportional to $\kappa_{c}$ have been included in order to dynamically control the conformal constraints, exactly in the same way as it is done with the physical ones.

In order to close the system of equations, coordinate (or gauge) conditions for the evolution of the lapse and shift must be supplied. We use the Bona-Massó family of slicing conditions BM and the Gamma-driver shift condition alcub , namely

[TABLE]

being $f$ and $g$ arbitrary functions depending on the lapse and the metric, and $\eta$ a constant parameter.

II.2 Characteristic structure

Now that the system is complete, it is possible to calculate its characteristic structure. Here we use the concept of pseudo-hyperbolicity Z43 ; CCC , that relies on a plane-wave analysis applied to the linearized equations around a background metric 222Notice also the work in Ref. nagy that extend these ideas using pseudo-differential operators.. Hence, we consider the line element

[TABLE]

and study the dynamics of perturbations over this background spacetime which propagates along a given normalized direction $s_{i}$ (i.e., such that $\gamma_{ij}^{0}s^{i}s^{j}=1$). The perturbation for the metric fields $\{\alpha,\beta^{k},\tilde{\gamma}_{ij},\chi\}$ has a plane-wave form,

[TABLE]

where $\omega^{k}$ is the wavenumber and $\omega\equiv\omega_{k}\,s^{k}$. An additional factor $i\omega$ appears in the perturbations of the fields $\{\tilde{A}_{ij},\hat{K},\Theta,\hat{\Gamma}^{i},B^{i}\}$, which are first derivatives of the metric, namely

[TABLE]

Replacing the above mentioned definitions in (8 - 13) one can obtain the following system:

[TABLE]

where $\bar{u}$ is a vector containing the perturbation of the fields, $\mathcal{A}$ is the characteristic matrix and $\mathrm{I}$ the identity one. The index $s$ means a contraction with the propagation direction $s_{i}$ (i.e., $\bar{\beta}^{s}_{0}=s_{i}\bar{\beta}_{0}^{i}$). The projection orthogonal to $s_{i}$ will be denoted by the index $\perp$.

The system (24) is pseudo-hyperbolic if and only if the characteristic matrix $\mathcal{A}$ has real eigenvalues and a complete set of eigenvectors. Therefore, with this definition, hyperbolicity of the system translates into a set of algebraic conditions reula . The analysis of the characteristic structure can be simplified by splitting the perturbations in different sectors which do not interact (i.e., or at least, not strongly) with the others.

It is instructive to analyze first the effect of the term proportional to $\lambda_{0}$. There is a sector, involving only the perturbations of $\tilde{\gamma}$ and $tr\tilde{A}$, given by:

[TABLE]

where $\approx$ means that only the principal part is considered. Obviously, for the original choice $\lambda_{0}=0$, there is not a complete set of eigenvectors. Then, the system is only weakly pseudo-hyperbolic system and, consequently, the problem is ill-posed. The same problem appears for any other value except for $\lambda_{0}=1$. Only for this choice both $\{\tilde{\gamma},\tilde{A}\}$ are standing modes, implying that this sector has a complete set of eigenvectors. As it is shown next, the other sectors are also complete, meaning that the full system is strongly pseudo-hyperbolic. Notice that this lack of strong hyperbolicity (together with the unbound growth of the conformal constraints) prevents to evolve directly the unconstrained CCZ4, unless the conformal constraints are algebraically enforced during the evolution hilditch .

We can now study the characteristic structure of the other modes. The lapse sector, constituted by $\{\bar{\alpha},tr\bar{K}\}$, has a complete set of eigenvectors with eigenvalues $-\bar{\beta}^{s}_{0}\pm\alpha_{0}\sqrt{f}$. The longitudinal shift and energy modes form another closed sector with a complete set of eigenvectors including $\{\bar{\chi},\bar{\Theta},\bar{\Gamma}^{s},\bar{\beta}^{s},\bar{B}^{s}\}$ with characteristic speeds given by $\{-\bar{\beta}^{s}_{0},-\bar{\beta}^{s}_{0}\pm\sqrt{4\,g/3\chi_{0}},-\bar{\beta}^{s}_{0}\pm\alpha_{0}\}.$ The transverse shift sector, including $\{\bar{\beta}^{\perp},\bar{\Gamma}^{\perp},\bar{B}^{\perp}\}$, is also complete with characteristic speeds $\{-\bar{\beta}^{s}_{0},-\bar{\beta}^{s}_{0}\pm\sqrt{g/\chi_{0}}\}$. The light sector also has a complete set of eigenvectors including the projections $\{\overline{\gamma}_{\perp\perp},\bar{A}_{\perp\perp},\overline{\gamma}_{s\perp},\bar{A}_{s\perp}\}$, with characteristic speeds $\{-\bar{\beta}^{s}_{0}\pm\alpha_{0}\}$.

Finally, notice that the choice $g=3/4$ is especially delicate because the characteristic velocities of the longitudinal shift modes collapse to light speed: no complete set of eigenvectors can be found, and the strong pseudo-hyperbolicity of the system is spoiled. This might be a problem, at least around Minkowski spacetimes, as it has been reported previously by several authors Z43 ; beyer ; gundlach2 .

II.3 Klein-Gordon equation

The evolution of a complex scalar field $\phi$ is described by the Klein-Gordon equation

[TABLE]

where $V(|\phi|^{2})$ is the potential depending only on the scalar field magnitude. By using the 3+1 decomposition, and introducing $\Pi\equiv-\mathcal{L}_{n}\phi$ as a new evolved field, the Klein-Gordon equation (26) can be written as

[TABLE]

or, by using the conformal fields of the CCZ4 formalism,

[TABLE]

In order to study self-gravitating boson stars we need to define the stress-energy tensor produced by this complex scalar field333Notice that there is a difference of a factor $1/2$ with respect to other definitions. mace

[TABLE]

where $\phi^{*}$ is the complex conjugate of $\phi$. One can also take advantage of the $U(1)$ symmetry to define the Noether charge, given by:

[TABLE]

where

[TABLE]

The Noether charge can be interpreted as the number of bosonic particles liebpa .

III Numerical tests

Our modification of the CCZ4 formalisms is examined by evolving known solutions –robust stability, gauge waves and single boson star–, some of which are included as standard testbeds apples . Here we describe first the setup of our simulations and present the tests results.

III.1 Setup

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 ander . To ensure sufficient resolution within the compact objects (i.e., boson stars in this work) in an efficient manner, we employ adaptive mesh refinement (AMR) via the had computational infrastructure that provides distributed, Berger-Oliger style AMR had ; lieb with full sub-cycling in time, together with an improved treatment of artificial boundaries lehner . We adopt a Courant parameter of $\lambda_{c}\approx 0.25$ such that $\Delta t_{l}=\lambda_{c}\,\Delta x_{l}$ on each refinement level $l$ to guarantee that the Courant-Friedrichs-Levy condition is satisfied. This code has been used extensively for a number of other projects and it has already been rigorously tested.

III.2 Robust stability test

We first carry out the robust stability test-bed alcub2 ; CCC to confirm the characteristic structure of our system (8-13). The test consists on a Minkowski background metric plus a small random perturbation in each of the evolution fields, such that only the principal part and the linear terms are significant (i.e., the scalar field potential $V$ and all the damping coefficients are set to zero in this test). A linear growth on any field indicates a weakly hyperbolic system. We set a 2D domain $[-0.5,0.5]^{2}$ with periodic boundaries conditions and $N=100$ grid points in each direction. No artificial Kreiss-Oliguer dissipation is included for this test. As it was shown in subsection II.2, the hyperbolicity of the system depends on the parameter $\lambda_{0}$, so we analyze the effect of this parameter. Besides, we use $f=2/\alpha$, $g=3/4$ and $\eta=2$, that is, the $1+\log$ slice with common values for the Gamma-driver shift condition.

The $L_{2}$-norms of some constraints are displayed in Fig. 1 for three cases. In the first one (CCZ4e) the conformal constraints are algebraically enforced after each timestep, as it is currently done in all the flavors of BSSN and conformal Z4. The other two cases (CCZ4) corresponds to $\lambda_{0}=0$ and $\lambda_{0}=1$ without any algebraic enforcing. Our simulations show a linear growth on $||\tilde{\gamma}-1||$ that propagates to $||Z_{x}||$ and, eventually, to all the other fields. This linear growth indicates a lack of strong hyperbolicity of the system for $\lambda_{0}=0$. All the norms are constant both for CCZ4e and for CCZ4 with $\lambda_{0}=1$, as it is expected for a well-posed system. Henceforth, we are going to use the choice $\lambda_{0}=1$ for all the forthcoming simulations with CCZ4.

III.3 Gauge waves

A family of non-trivial exact solutions can be constructed by performing a coordinate transformation on $\{x,t\}$ to Minkowski spacetime. The resulting line element can be written as alcub2 ; CCC ; alic ; apples

[TABLE]

where $H(x-t)=1-A\sin\left[k(x-t)\right]$. We set an amplitude $A=0.1$ and a wave number $k=2\pi/L$, being $L$ the size of the domain. This solution is exact with harmonic slicing and zero shift, corresponding to the choice $f=1$ and $g=0$. The domain for this one-dimensional test is $[0,1]$ with periodic boundary conditions and $100$ grid points.

$L_{\infty}$-norm for some constraints are shown in Fig. 2 for four different cases, which can be compared with figure 1 in alic . In the first two cases (i.e., BSSN and CCZ4e), the conformal constraints are algebraically enforced after each time step. The other two cases correspond to CCZ4 with either $\kappa_{c}=1/L$ or $\kappa_{c}=0$ (in both cases $\kappa_{z}=1/L$). It is clear that the BSSN formulation fails this test, as both the conformal and the Hamiltonian constraint suffer exponential growth. In contrast, all the constraints remain under control when using the CCZ4e. The most important outcome of this test is the fact that the CCZ4 formulation with $\kappa_{c}=1$ is also stable, meaning that it is not required to enforce algebraically the conformal constraints to keep them under control. The last case, CCZ4 with $\kappa_{c}=0$, presents a linear growth in the conformal constraint $|\tilde{\gamma}-1|$, which unavoidably will lead to a failure due to the propagation to other fields. The same behavior is observed in simulations on generic spacetimes, indicating that the choice of the damping coefficients $\{\kappa_{z},\kappa_{c}\}$ is crucial to achieve accurate and stable solutions.

III.4 Single solitonic boson star

The initial data for complex scalar field configurations in spherical symmetry can be solved numerically for the static metric liebpa ; mace

[TABLE]

by adopting the following harmonic ansatz for the scalar field

[TABLE]

Within these assumptions the EKG system reduces to a set of ordinary differential equations (ODEs) that can be solved by imposing appropriate boundary conditions (i.e., regularity at the origin and asymptotically flat at large distances). The equations are further simplified by using polar-areal coordinates. Therefore, the standard procedure is to solve the equations in these coordinates and then perform a (numerical) coordinate transformation into isotropic coordinates, which can be transformed easily to Cartesian ones 2004PhDT…….230L .

Different interaction potentials $V(|\phi|^{2})$ lead to boson stars with different compactness liebpa . We are interested in a particular family of very compact boson stars, commonly known as nontopological solitonic boson stars frie ; 1992PhR…221..251L , where the potential is given by

[TABLE]

Here $\sigma_{0}$ is a constant that determines the compactness of the star and $\mu$ is related to the scalar field mass. By setting $\sigma_{0}=0.05$ and the scaling factor $\mu\sigma_{0}\sqrt{8\pi}=1$, a suitable stable equilibrium configuration can be obtained for $\phi_{0}(r=0)=0.0364$ and $\omega=1.0666$. The resulting star has mass $M=0.36$ and radius $R=3.08$, so its compactness is $C=M/R=0.118$. This configuration is well inside the stable branch, since the most massive stable star has a mass $M_{max}\approx 1.84$. The profiles of $\alpha(r)$, $\psi(r)$ and $\phi_{0}(r)$ for this particular solution are plotted in Fig. 3.

This configuration is evolved in a domain $[-16,16]^{3}$ with radiative boundary conditions. There are $60$ grid points in each direction and three refinement levels, such that the highest resolution is $\Delta x=0.1$. The simulations are performed by using the CCZ4 formulation with $\kappa_{z}=0.1$ and either $\kappa_{c}=1$ or $\kappa_{c}=0$. We also include the solutions obtained with the BSSN formulation for comparison purposes.

The evolution of the scalar field real part $\phi_{R}(t,r=0)$ is displayed in Fig. 4, together with the expected analytical behavior $\phi_{0}(r_{0})\cos(\omega t)$. The solutions for all the cases considered, either with BSSN or CCZ4, show a very good agreement with the analytical expectation. Differences arise however in the $L_{2}$-norm of some constraints, plotted in Fig. 5. Both the conformal and the physical constraint remain under control by using either BSSN or CCZ4 with $\kappa_{z}=0.1$ and $\kappa_{c}=1$. However, notice that the errors of the physical constraints obtained with the CCZ4 and this parameter choice are several orders of magnitude smaller than the ones obtained by using BSSN.

IV Dynamics of solitonic boson stars

Here the dynamics of binary boson stars is studied, focusing on the final state of the system after the merger. We first describe how to construct initial data for a binary boson star system with generic angular momentum. We consider several head-on cases, and finish with orbiting binary systems.

IV.1 Initial data and setup

We extend the procedure describe in shi to construct accurate boosted initial data from an spherically symmetric solution. Our starting point is the line element given by

[TABLE]

where $\alpha_{0}=\alpha_{0}(r_{0})$ and $\psi_{0}=\psi_{0}(r_{0})$, being $r_{0}=\sqrt{x_{0}^{2}+y_{0}^{2}+z_{0}^{2}}.$ By performing a Lorentz transformation $t=\Gamma(t_{0}+vx_{0}),\,x=\Gamma(x_{0}+vt_{0})$ (being $\Gamma=1/\sqrt{1-v^{2}}$), one can obtain

[TABLE]

The gauge choice is given by

[TABLE]

with $B_{0}=\Gamma\,\sqrt{\left(1-\frac{v^{2}\alpha_{0}^{2}}{\psi_{0}^{4}}\right)}$. Notice that $r_{0}$ can be written in terms of the new coordinates, namely $r_{0}=\sqrt{\Gamma^{2}(x-vt)^{2}+y^{2}+z^{2}}$.

Now we only have to perform the Lorentz transformation to the scalar field quantities. First, the harmonic ansatz given by eq. (33) can be generalized, to allow for non-identical boson stars, by including a phase shift $\theta$ and the direction of rotation $\epsilon=\pm 1$liebpa , namely:

[TABLE]

We can compute the field $\Pi(t,r_{0})$ from its definition, calculated in the boosted frame. The final expressions, evaluated at $t=0$, are

[TABLE]

We are mainly interested on binary systems. Along the lines described in pale1 , initial data for head-on binaries can be constructed as the superposition of two single solitonic boson star solutions, located at positions $r_{i}$, in the following way

[TABLE]

where superindex $(i)$ indicate each star. The initial data for orbiting binaries can be constructed in an analogous way, but using the boosted boson star –with velocity $v$ along the $x$-direction– as the basic solution and performing the superposition of all the non-trivial evolution fields.

For the head-on cases each solitonic BS is centered at $(\pm x_{c},0,0)$, while that for the orbiting case they are centered at $(0,\pm y_{c},0)$, being $x_{c}=y_{c}=4$. In both cases the cubical domain is given by $[-60,60]^{3}$ and $120$ grid points in each axis, leading to a coarsest resolution $\Delta x_{0}=1.0$. We set five refinement levels such that the last one, covering both stars, has a resolution $\Delta x_{4}=0.0625$ (i.e., there are approximately 96 points covering each star).

Besides the constraints, additional analysis quantities are evaluated in our binary boson stars simulations: the Noether charge or boson number, given by the volume integral eq. (30), and the ADM mass and angular momentum as described in the Appendix. The surface integrals, eq. (50,A), are computed in a sphere located at $R_{ext}=40$.

IV.2 Head-on cases

We consider initial data constructed as described in the previous subsection with the ansatz (38) and no boost. Even for a particular solution (i.e., the one described in subsection III.4), located at fixed initial positions, there is an infinite family of configurations depending on the parameters $\{\theta,\epsilon\}$. We differentiate two essentially different scenarios: (i) two identical boson stars, that is, with the same phase direction of rotation (i.e., $\epsilon=+1$), that will be denoted as B-B($\theta$), and (ii) a binary formed by a boson and an anti-boson star, denoted as B-aB($\theta$), such that their phases have opposite rotation direction (i.e., $\epsilon=-1$) and their Noether charges have opposite sign. We consider four different cases for B-B($\theta$) with phase shifts $\theta=\{0,\pi/2,\pi,3\pi/2\}$, and two for B-aB($\theta$) with $\theta=\{0,\pi\}$.

The evolution of the ADM mass and the Noether charge is presented in Fig. 6, and some snapshots of the evolution for all cases are displayed in Figure 7. Let us start by describing the dynamics of the binaries formed by two boson stars. The collision of the B-B([math]), B-B($\pi/2$) and B-B($3\pi/2$) cases leads to a merger, resulting into a single solitonic boson star with roughly the same total initial mass and Noether charge (i.e., except a small fraction that is emitted by gravitational waves and scalar field radiation, see Fig. 6). However, the scalar field interaction in the B-B($\pi$) produces a repulsive force that overcomes the gravitational attraction. Therefore, the binary suffers several inelastic collisions before relaxing to a system with two touching stars –which do not merge into a single one. The binaries formed by a boson and an antiboson star annihilates each other during the merger for the two extreme phase shifts considered (i.e., $\theta=0,\pi$). Most of the scalar field is radiated away to infinity, and only a small fraction remains near the region of the collision.

IV.3 Orbiting cases

The initial data for the orbiting cases is constructed as described in subsection IV for identical stars (i.e., with $\epsilon=1,\theta=0$). We have considered different Lorentz boost velocities $v=\{0,0.05,0.10,0.15\}$ along the $x$-direction. The last case corresponds to a binary almost in quasi-circular orbit.

The ADM mass, the Noether charge and the angular momentum are shown in Fig. 8, whereas few snapshots of the evolutions are displayed in Figure 9. All the cases –with angular momentum– merge and form a rotating bar that quickly losses angular momentum and settles down to a non-rotating boson star.

Although the system only losses a small fraction of the mass and the boson number during the coalescence, all the angular momentum is emitted by gravitational waves and scalar field radiation soon after the merger. This prompt radiation of angular momentum is most extreme for the system in quasi-circular orbit, corresponding to $v=0.15$: after the merger, two blobs of scalar field are ejected from the system at almost light speed, carrying a small fraction of the mass and boson number but a large amount of angular momentum (i.e., see $t=80$ of the $v=0.15$ in Figure 9, and the sudden drop of angular momentum in Fig. 8).

This behavior can possibly be explained by the quantization of (adimensional) angular momentum $J^{rot}_{z}=kN$ (being $k=$ an integer) for rotating boson star configurations, which might be difficult to achieve exactly by a dynamical merger. The resulting star from the merger should have at least $J^{rot}_{z}(k=1)\approx 0.62$ in order to correspond to a stable rotating boson star (see the bottom panel of Fig. 8 at $t\approx 50$). Although the (adimensional) angular momentum of the system $J_{z}/M^{2}\approx 0.78$ is larger than $J^{rot}_{z}(k=1)$, the system does not relax to that state but to the one with the lowest angular momentum $k=0$. Stated in a different way, the system does not dynamically evolve towards a rotating boson star configuration, but towards a non-rotating one.

Although this seems in contradiction with the results in pale2 , where a rotating boson star seemed to be produced, there are two important differences. First, these stars are much more compact, so the dynamics might be more dominated by non-linear gravitational effects. Second, due also to the high compactness, the dynamics is faster (i.e., the crossing time is shorter), so we can follow the evolution until a stationary state is achieved, which might not had been possible with the mini-boson stars.

Finally, the evolution of some constraints is displayed in Fig. 10, showing that they are small and kept under control during all the simulation.

V Conclusions

In this work we have studied the merger of binary compact boson stars, focusing on two different cases: head-on collisions of boson-boson/boson-antiboson for a wide range of shift phases, and orbital mergers of identical stars with different velocities.

These simulations have been performed by solving the Einstein-Klein-Gordon system with a novel modification of the CCZ4 formalism which does not require the algebraic enforcement of the conformal constraints after each step of the numerical integration. Our formulation treats these conformal constraints in the same way than the physical ones, and they are also kept under control by including damping terms. We have studied the pseudo-hyperbolicity of the evolution system using a linear plane-wave analysis to show that, only for the preferred choice $\lambda_{0}=1$, the system has a complete set of eigenvectors with real eigenvalues. Therefore, the linearized system is well posed, provided that $\lambda_{0}=1$. We have performed some standard numerical tests to check the robustness of our implementation, including the evolution of an isolated solitonic boson star. Our numerical results –stable evolutions with small and bounded constraints– confirm the analytical analysis.

Within this formulation we have studied binary solitonic boson stars systems with and without angular momentum. All the binary systems considered were constituted by two equal-mass solitonic boson stars. In the head-on cases we allowed for non-identical boson stars with different rotation and shift phase. Our simulations show that the merger of a boson-boson binary leads, in general, to another solitonic boson stars. However, when the phase shift approaches $\pi$, the scalar field interaction is repulsive and stronger than gravity, preventing the merger of the two stars. The merger of a boson star and an anti-boson star completely annihilates each other for any of the phase shift considered. This behavior, combined with the results described in pale1 ; pale2 , allow us to hypothesize that the generic behavior during the collision of a boson and an antiboson star is the annihilation of both, independently on the interaction potential and the phase shift, producing large amounts of unbound scalar field that is radiated to infinity.

In the scenario of orbital binaries, our studies with identical boson stars revealed that the merger always lead to the formation of a rotating bar which sheds quickly all its angular momentum by emitting scalar field and gravitational waves, to finally relax into a non-rotating boson star. This inability to form a rotating boson star from a merger might be due to the angular momentum quantization of the rotating solutions. Of particular interest is the case with the highest angular momentum considered, leading roughly to a system in quasi-circular orbits. In this case, soon after the merger, two blobs of scalar field, carrying away small amount of Noether charge but large fraction of angular momentum, were expelled from the remnant almost at light speed.

Future studies will further study orbital binary systems by considering different masses and analyzing the gravitational waves produced during the coalescence.

Acknowledgements

It is a pleasure to thank Vitor Cardoso, Luis Lehner, Steve Liebling and Paolo Pani for interesting discussions on the subject and helpful comments on the manuscript, as well as to Hiro Okawa for useful suggestions on the initial data. We acknowledge support from the Spanish Ministry of Economy and Competitiveness grants FPA2013-41042-P and AYA2016-80289-P (AEI/FEDER, UE). CP acknowledges support from the Spanish Ministry of Education and Science through a Ramon y Cajal grant. MB would like to thank CONICYT Becas Chile (Concurso Becas de Doctorado en el Extranjero) for financial support.

VI Appendix

Appendix A Analysis quantities

Hence, the three dimensional Christoffel symbols tensor can be written in function of conformal metric as:

[TABLE]

Besides, the BSSN conformal BSSN connections is defined by

[TABLE]

There are several analysis quantities in the CCZ4. The first term of the momentum constraint can be computed directly from this relation

[TABLE]

Assuming that $\theta=Z_{i}=0$, we have that Hamiltonian and momentum constraint are given by

[TABLE]

Other interesting quantities are ADM mass and angular momentum, they can be computed by performing a surface integral at spatial infinity

[TABLE]

in terms of the conformal quantities, where ${\hat{\phi}}_{j}=(-y,x,0)$ is the cartesian coordinate basis axial vector. The surface element is given by

[TABLE]

Other interesting analysis integral quantities, which are exact for stationary and stationary spacetimes, are the Komar mass and angular momentum gour :

[TABLE]

where $\zeta_{t}^{a}=(1,0,0,0)$ and ${\zeta}_{\varphi}^{a}=(0,{\hat{\phi}}^{i}/||{\hat{\phi}}^{i}||)$ are the time-like and the axial Killing vectors.