Head-on collisions and orbital mergers of Proca stars

Proca stars are self-gravitating Bose-Einstein condensates obtained as numerical stationary solutions of the Einstein-(complex)-Proca system. These solitonic can be both stable and form dynamically from generic initial data by the mechanism of gravitational cooling. In this paper we further explore the dynamical properties of these solitonic objects by performing both head-on collisions and orbital mergers of equal mass Proca stars, using fully non-linear numerical evolutions. For the head-on collisions, we show that the end point and the gravitational waveform from these collisions depends on the compactness of the Proca star. Proca stars with sufficiently small compactness collide leaving a stable Proca star remnant. But more compact Proca stars collide to form a transient ${\it hypermassive}$ Proca star, which ends up decaying into a black hole, albeit temporarily surrounded by Proca quasi-bound states. The unstable intermediate stage can leave an imprint in the waveform, making it distinct from that of a head-on collision of black holes. The final quasi-normal ringing matches that of Schwarzschild black hole, even though small deviations may occur, as a signature of sufficiently non-linear and long-lived Proca quasi-bound states. For the orbital mergers, the outcome also depends on the compactness of the stars. For most compact stars, the binary merger forms a Kerr black hole which retains part of the initial orbital angular momentum, being surrounded by a transient Proca field remnant; in cases with lower compactness, the binary merger forms a massive Proca star with angular momentum, but out of equilibrium. As in previous studies of (scalar) boson stars, the angular momentum of such objects appears to converge to zero as a final equilibrium state is approached.


I. INTRODUCTION
The recent spectacular detections of gravitational waves [1][2][3][4][5][6] opened up a new window into the strong-field regime of gravity [7]. Even though the data so far is well fitted by the expected physics-i.e., collision of Kerr black holes (BHs) or neutron stars-it is important to understand if alternative, nonconventional models of compact objects can also fit the data (i.e., the level of degeneracy) or how much these can be ruled out by current/future data-see [8] for such a discussion. In the best-case scenario, obtaining gravitational waveforms of such nonconventional objects could lead to their future discovery, and the exciting prospect of unveiling new surprising physics via the gravitational-wave window. Amongst such exotic models, self-gravitating solitons composed of complex (boson stars [9]) or real (oscillatons [10]) scalar fields are some of the most dynamically studied cases-see, e.g., [11]. Dynamical studies include the generation of waveforms for head-on collisions of boson stars [12,13] and oscillatons [14][15][16] as well as orbital mergers of boson stars [17][18][19][20][21]. Even if the orbital mergers describe the astrophysically more likely scenario [22,23], and constitute a possible mechanism to form spinning boson stars or Kerr BHs with scalar hair [24,25], the headon collisions represent a first step to compute the gravitational waveforms produced by these objects, allowing a simpler comparison to those produced in head-on collisions of BHs or neutron stars.
In has been recently found that a complex Proca field can form Proca stars (PSs), vector analogues of the scalar boson stars. PSs were constructed as static or stationary solutions of the Einstein-(complex)Proca system [26]-see also [27][28][29][30] for generalizations. Similarly to their scalar cousins, they can be regarded as single frequency, macroscopic, self-gravitating (vector) Bose-Einstein condensates. This frequency appears as a harmonic time dependence for the Proca potential and the domain of existence of PSs is very similar to that of scalar boson stars-see, e.g., Fig. 1 in [31]-except that the latter have a smaller maximal mass and a wider frequency range.
Several dynamical aspects of the gravitational interaction of Proca fields have been addressed in recent years, using numerical studies. Long-lived, quasibound states of Proca fields around Schwarzschild BHs have been considered in [32,33] (see also [34]). The superradiant instability of Kerr BHs was triggered by Proca fields in [35,36], leading, in particular, to the formation of Kerr BH with Proca hair [27,37]. In [38] fully nonlinear evolutions of PSs were performed, to assess their stability. The simulations showed the existence of a stable branch (connecting the vacuum with the solution with maximal ADM mass) and an unstable branch. Solutions belonging to the latter may have different fates, depending on the initial perturbation and the sign of their binding energy. As their scalar cousins, with and without a self-interacting term [39][40][41][42], unstable solutions with positive binding energy migrate to the stable branch, whereas unstable solutions with a negative binding energy (excess energy) undergo fission; i.e., they disperse entirely. Both cases can also collapse to form a Schwarzschild BH. More recently, additional numerical simulations [43] have shown that PS can also form dynamically from generic initial data describing a non-compact "cloud" of the Proca field, through the so-called gravitational cooling mechanism, first described in the 1990s in the context of scalar boson stars [44].
In this work we shall continue the exploration of the dynamics of PSs by performing numerical evolutions describing both head-on collisions and orbital mergers of PS binaries and computing the gravitational radiation emitted. In the case of head-on collisions, the results we obtain parallel qualitatively those for head-on collisions of scalar boson (or even oscillaton) stars. For sufficiently small compactness, or equivalently, low mass in units of the Proca field mass-see Fig. 1 below-the collisions form a more massive, but still stable, PS, not a BH. The final star is, however, perturbed and in the timescale of our simulations it only partially relaxes to equilibrium. Sufficiently compact PSs, on the other hand, form a horizon when colliding, but only after an intermediate phase that could be described as a hypermassive PS. This intermediate stage leaves an imprint in the waveform, making it distinct from that of a head-on collision of Schwarzschild BHs. After horizon formation the BH ringdown can be seen, which matches well that of a Schwarzschild BH. But in the cases where a larger Proca remnant remains outside the horizon, in the form of Proca quasibound states, we observe a difference with the BH ringdown. A similar observation was reported recently in the study of head-on collisions of oscillaton stars [16].
For the orbital mergers, similar results to the head-on case are also found. The mergers are eccentric and cover less than an orbit. They can lead to the formation of a Kerr BH wherein part of the initial orbital angular momentum of the configuration is deposited. In such cases, a Proca field remnant can be seen outside the horizon. This remnant is a quasi-bound state that decays exponentially. Thus, with the initial parameters chosen we do not see the formation of a Kerr BH with Proca hair. On the other hand, we obtain the formation of a massive PS with angular momentum for small compactness. The solitonic remnant is, however, out of equilibrium, and the angular momentum decreases significantly during and after the merger, approaching zero. In this process, the system emits gravitational waves continuously.
This paper is organized as follows. In Sec. II we present the equations of the Einstein-(complex)Proca model that will be used for the numerical evolutions. In Sec. III we present the initial data that will be used in our numerical evolutions. A brief description of the numerical techniques is given in Sec. IV and our results are presented in Sec. V. Final remarks are presented in Sec. VI. A brief assessment of the numerical code is given in Appendix.

II. BASIC EQUATIONS
We shall investigate the dynamics of a complex Proca field by solving numerically the Proca equations coupled to the Einstein equations. The system is described by the where the Lagrangian density depends on the Proca potential A, and field strength F ¼ dA; it is given by: where the bar denotes complex conjugation, R is the Ricci scalar, G is Newton's constant, and μ is the Proca field mass. The stress-energy tensor of the Proca field reads Using the standard 3 þ 1 split (see, e.g., [38] for more details) the Proca field is split into 3 þ 1 quantities: where X i is the vector potential and X ϕ is the scalar potential. The fully nonlinear Einstein-Proca system reads: where α is the lapse function, β is the shift vector, γ ij is the spatial metric, K ij is the extrinsic curvature (with K ¼ K i i ), D i is the covariant 3-derivative, and L β is the Lie derivative. Moreover, the three-dimensional "electric" E i and "magnetic" B i fields are also introduced in the previous equations. The system is solved using the time-evolution numerical code from [33] (see Sec. IV).

III. INITIAL DATA
PSs were obtained in [26] as stationary solutions to the model described by the action (1). Five illustrative examples of spherically symmetric PSs will be taken as the initial data for our time evolutions. Their basic physical properties, frequency, w, ADM mass, M ADM , Noether charge Q and the Proca "electric" potential at the origin, Φ c , all in units of the vector field mass, can be found in Table I and their distribution in an ADM mass vs. Proca field frequency diagram is shown in Fig. 1.
When computed as static solutions [26], spherically symmetric PSs are given by the line element where F 0 , F 1 are radial functions and r; θ; φ correspond to isotropic coordinates. The Proca field ansatz is given in terms of another two real functions ðV; H 1 Þ which depend also on r where w > 0 is the frequency of the field. The translation between the four radial functions above, F 0 ; F 1 ; V; H 1 , and the initial value for the metric and the 3 þ 1 Proca field variables described is given as follows: We follow [12] to construct appropriate initial data to study the head-on collision of these compact objects. We take a superposition of two PS solutions: where superindex ðiÞ labels the stars and AEx 0 indicates their initial positions. The stars are initially separated by Δx ¼ 39, which corresponds to x 0 ¼ AE19.5 (in G ¼ c ¼ 1 units). The solutions are not boosted. These initial data introduce constraint violations. However, they are small and do not grow during the evolution, as we discuss in Appendix.
For the construction of initial data describing orbiting PS binary we extended the same method of the superposition of two isolated PSs solution, with the only difference that in this case the two stars are boosted along the y-axis which is perpendicular to the line segment linking them, with velocity AEv y , following [23,45]. In this case, 5he stars are initially separated by Δx ¼ 30, which corresponds to x 0 ¼ AE15.0. We call B ðiÞ b ¼v ðiÞ y =c and Γ the Lorentz factor, and the matrix associated with the transformation has the following form with the features that Λ T ¼ Λ and Λ −1 can be obtained using opposite velocity v y . Note that Γ The line element of each star is in Cartesian coordinates given by We perform a Lorentz transformation t ¼ Γ b ðt 0 þ v y y 0 Þ and y ¼ Γ b ðy 0 þ v y t 0 Þ and obtain from (20) The extrinsic curvature is computed from: We have to transform the Proca fields. The Lorentz transformation resembles the one for the common electromagnetic fields. We consider that in the rest frame the magnetic field is zero, due to the spherical symmetry. The boosted fields are obtained as follows: The initial data for the PS binary is a superposition of the two boosted solution as in the head-on case.

IV. NUMERICS
To perform the numerical evolutions we use the freely available EINSTEIN TOOLKIT [46,47], which uses the CACTUS framework and mesh refinement. The methodof-lines is employed to integrate the time-dependent differential equations. In particular, we use a fourth-order Runge-Kutta scheme for this task. The left-hand-side of the Einstein equations is solved using the MACLACHLAN code [48,49], which is based on the 3 þ 1 Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation.
The Proca evolution equations, Eqs. (6)- (11), are solved using the code developed independently by M. Zilhão and H. Witek [33]. We have extended this code to take into account a complex field. All technical details, assessment of the code and convergence tests can be found in [33].
We use a numerical grid with 6 refinement levels for the head-on collisions, its structure being fð512; 64; 64; 32; 32; 8Þ; ð4; 2; 1; 0.5; 0.25; 0.125Þg, where the first set of numbers indicates the spatial domain of each level and the second set indicates the resolution. Figure 2 shows the grid structure at t ¼ 0. The outer boundary is located at Δx ¼ 392 from the gravitational-wave extraction radii. This ensures that the total time of the simulations is shorter than twice the light-crossing time, which prevents numerical reflections at the boundary from affecting the extraction. Due to the geometry of head-on collisions, we consider equatorial-plane symmetry and reflection symmetry with respect to the x-z plane.
Correspondingly, Fig. 3 shows the grid structure we employ for the orbital mergers, which contains 7 refinement levels with spatial extent and resolution given by fð412;64;64;32;32;8;4Þ;ð4;2;1;0.5;0.25;0.125;0.0625Þg. In this case there is a reflection symmetry along the equatorial-plane but no reflection symmetry with respect to the x-z plane.

A. Head-on collisions
We evolve the initial data described in Sec. III for four equal-mass PS binaries, using the four types of PSs described in Table I and highlighted in Fig. 1. We label these four models as PS0-PS0, PS1-PS1, PS2-PS2, and PS3-PS3. We also evolve the latter including an initial perturbation which consists in multiplying by a number slightly larger than one the initial values of the Proca variables, as we did in [38], a model that we dub PS3b-PS3b. By including a perturbation the stars are forced to collapse before the collision, since the PS3 model is unstable against radial perturbations. Figure 4 shows equatorial-plane snapshots of the evolutions of the PS0-PS0, PS1-PS1, PS2-PS2, and PS3-PS3 binaries, from top to bottom. Our computational grid only extends from y ¼ 0 to 512, therefore the data for negatives values of y are mirrored by the corresponding positive y values, due to axisymmetry. Animations illustrating these collisions can be found in [50].

Visualisation of the collisions
The collision of the two Proca stars happens at t ∼ 160 for models PS0-PS0 and PS3-PS3, and at t ∼ 192 for PS1-PS1 and PS2-PS2. The merged object oscillates in the x and y directions. For PS0-PS0, part of the field is ejected from the polar caps and approaches spherical symmetry.
Model PS2-PS2 collapses quickly after the collision at t collapse ∼ 245, while models PS1-PS1 and PS3-PS3 oscillate during a short period of time before an apparent horizon (AH) appears at t collapse ∼ 375 and t collapse ∼ 425, respectively. Moreover, the PS1-PS1 collision forms a perturbed massive PS that oscillates and produces the first part of the gravitational-wave signal, as we show below. In the case of PS3-PS3, the stars migrate and expand before the collision. For all models, the morphology and amplitude of the gravitational waves are clearly influenced by the dynamics during the collision (cf. Sec. VA 3).

Proca energy and horizon formation
In Fig. 5 we plot the time evolution of the Proca field energy "scalar" potential) for model PS0-PS0. In this case we see no horizon formation. The stars have sufficiently low mass so that the final state is still a star. The amplitude of the scalar potential stabilizes within the timescale of our simulation. It has been reported in [14,16] that some collapsed oscillatons from head-on collisions can lose part of the field and form a stable, less compact oscillaton, not collapsing to a BH even if the total mass is larger than the maximum mass of the equilibrium configuration. This seems to be precisely what is occurring for model PS0-PS0. The mass of the PS0 Proca star is μM ADM ¼ 0.693, so that twice that mass is larger than the maximum mass for spherical PSs, which is μM max ADM ¼ 1.058 [26]. Still, a PS forms, as a result of the inelasticity of the collision. Figure 6 displays the time evolution of the Proca field energy and of the irreducible mass of the AH for models PS1-PS1, PS2-PS2 and PS3-PS3. For these three models, the collapse of the solutions is triggered after the collision and an AH forms. After the collapse there is still a Proca field remnant outside the horizon, as Fig. 7 shows. For PS1-PS1 and PS3-PS3 this remnant has a visible dynamics, reminiscent of a beating pattern, a signature of the presence of more than one quasibound state outside the AH. In the case of the PS3-PS3 model, the Proca remnant seems particularly long lived. We shall observe a possible impact of this feature in the gravitational waveform below.  , for the modes l ¼ 2, m ¼ f0; þ2g, extracted at radii r ext ¼ f100; 120g. The waves, conveniently shifted and rescaled, overlap, as expected in the wave zone. Nonaxisymmetric modes l ¼ 2, m ¼ AE1 are consistent with zero. In agreement with the results for head-on collisions of spherical boson stars in [12], we find that the coefficients C 2;m of the different modes are related in the following way: For our sample of models, the amplitude of the gravitational waves increases monotonically with decreasing vector-field frequency. The largest amplitude is achieved for model PS3-PS3 which is more than an order of magnitude larger than for model PS0-PS0. All waveforms are of the burst-type, of which PS2-PS2 and PS3-PS3 are clear examples. Morphological differences are apparent in model PS0-PS0 and, in particular in model PS1-PS1. For the latter, plotted in the middle panels of Fig. 8, there is a time delay Δt ∼ 100 between the first negative peak of the waveform (at around t ∼ 275) and the collapse (at around t ∼ 375). The waveform can be regarded as composed by   two contributions: the first part, from t ∼ 250 to t ∼ 375, would correspond to the collision of the stars, forming a "hypermassive" PS, a remnant that oscillates and eventually collapses to a BH at t ∼ 375, triggering the second part of the wave. The collapse nearly coincides with a peak of the waveform and the corresponding time is highlighted by a vertical dashed line in the figure. For PS2-PS2 and PS3-PS3, these dynamics are absent, the AHs form promptly and their times of formation are very close to the first positive peak in the waveforms (vertical dashed lines).
The bottom panels of Fig. 9 show the gravitational waveform for the PS3b-PS3b model together with those of models PS2-PS2 and PS3-PS3. Recall that PS2 and PS3 have almost the same mass (less than 1% difference), therefore the BH formed will have a similar mass. Thus, we can compare the quasinormal modes (QNMs) of the three cases. For PS2-PS2, the QNMs are in good agreement in at least four of the peaks, while for PS3-PS3 there are some differences in the frequency. This comparison becomes more clear in the QNM ringdown plots shown in Fig. 10, where we plot the waveforms in logarithmic scale and we fit them to the QNMs of a Schwarzschild BH, following the results of [51]. PS2-PS2, PS3-PS3, and PS3b-PS3b have almost the same total mass, therefore the fit is the same. The agreement for PS2-PS2 and PS3b-PS3b is very good. The masses that we obtain from the fits have an error of 7% for PS1-PS1 and 4% for PS2-PS2, PS3-PS3, and PS3b-PS3b, with respect to the masses computed with the AHFINDER algorithm of the EINSTEIN TOOLKIT. However, for the PS3-PS3 waveform the frequency does not match the PS3b-PS3b nor the QNM fit. A possible interpretation is that this is due to the fact that there is still a rich Proca field environment around the newly form BH, cf.

Comparison with BH collisions
We also perform head-on collisions of Schwarzschild BHs with the same mass as the PSs in each model. This  allows to compare the resulting gravitational waveforms with those from the PS binaries. The real part of rΨ l¼2;m¼0 4 for both cases is shown in Fig. 11. Observe that for model PS3b-PS3b, we obtain essentially the same gravitational waveform. This is the expected result; in this case, the PSs collapse to BHs before the actual collision due to the initial perturbation we have imposed, thus leading to a head-on collision of BHs. In the case of PS2-PS2, the ringdown phase overlaps almost perfectly that same phase for the BH collision in agreement with the QNMs analysis in Fig. 10. The first part of the wave, however, is different: the frequency is higher but the amplitude is slightly smaller at the peak of the emission.
On the other hand, as already shown in the previous subsection, for PS1-PS1 and PS3-PS3, the comparison indicates that the waveforms depart from those of the Schwarzschild BH collisions. The amplitudes are different for both models. Some of the QNM oscillations can be fitted with the Schwarzschild BH-BH waveforms but not all, as the frequency changes, see Fig. 12. The deviation in the frequency could be a signature of the presence of quasibound states around the BH and might be in particular related to the compactness of the remaining Proca field around the BH. The gravitational radiation induced by accreting shells of matter evolving in fixed BH backgrounds was first studied in [52] by numerically solving the linearized curvature perturbation equations (see also [53] for a dynamical spacetime study). It was found that the excitation of the BH QNM ringing strongly depends on the shell thickness, becoming increasingly clear with progressively more compact shells (see Fig. 10 in [52]). In the infinitesimally thin limit, the gravitational energy asymptotes to a finite value, about a third of the point particle upper limit. Those findings confirmed earlier ideas about the QNM excitation mechanism made by [54], namely that the strong excitation is induced by curvature profiles that have spatial wavelengths comparable to the width of the BH potential.

B. Orbital binary simulations
We turn now to describe the orbital binary merger simulations of models PS00-PS00, PS1-PS1, PS2-PS2, and PS3-PS3. The stars are boosted in the y direction. In this case, the objects are initially separated by Δx ¼ 30.

Visualization of the binary mergers
The results are similar to the head-on collision case. The stars do not complete a full orbit since the mergers are eccentric; nonetheless, the final object has nonzero angular momentum. Figure 13 shows snapshots of the evolution of the energy density of models PS1-PS1, PS2-PS2, and PS3-PS3 in the equatorial plane for an initial velocity v y ¼ 0.050. The collision happens at t ∼ 128 for models  PS1-PS1 and PS3-PS3, and at t ∼ 144 for PS2-PS2. The merged object has still some angular momentum left. An AH appears at t ∼ 350 for model PS1-PS1, at t ∼ 200 for model PS2-PS2 and at t ∼ 300 for PS3-PS3.
In Fig. 14 we exhibit the evolution of the energy density for the model PS00-PS00, with three different initial velocities, v y ¼ 0.0125, v y ¼ 0.050 and v y ¼ 0.10. For PS00-PS00, the result of the merger does not lead to BH formation; instead the final object is a Proca star with angular momentum. The larger the initial velocity, the larger the fraction of the initial Proca mass and angular momentum that is ejected during and after the merger. With the largest initial velocity, v y ¼ 0.10, almost all the Proca field is dispersed away at the end of the simulation.

Proca energy and horizon formation
In Fig. 15 we plot the time evolution of the Proca field energy, Eq. (31), and the Proca field angular momentum together with the BH mass and the BH spin J for models PS1-PS1, PS2-PS2, and PS3-PS3 (bottom panels). In all these cases there is AH formation. These models collapse after the merger and form a Kerr BH with angular momentum. There is a remnant Proca field outside the rotating horizon, forming a quasibound state. We do not see evidence of infinitely long-lived Proca remnants (hair) forming around the BH. The binary takes longer to collapse for a larger initial velocity v y , but for the models PS2-PS2 and PS3-PS3 the difference is small for the velocities we have chosen. A large part of the total angular momentum is lost before the collapse for these models. On the other hand, for model PS1-PS1 the difference in collapse time is more noticeable: it takes about 50 times longer to collapse for the largest velocity and the final BH stores almost all of the initial angular momentum. Figure 16 shows that a quasibound state appears outside the horizon after the collapse and formation of a BH. As in the head-on collision case, the final amount of Proca field remaining is smaller for model PS2-PS2 than for models PS1-PS1 and PS3-PS3. For PS00-PS00 there is no BH formation. In Fig. 17, we show the time evolution of the energy and angular momentum of this model that we evolve with five different initial boost velocities, therefore increasing the initial angular momentum, namely v y ¼ f0.0125; 0.025; 0.050; 0.075; 0.10g. The PS formed after the collision has nonzero residual angular momentum. The total Proca energy is larger than the maximum mass of a rotating PS with m ¼ 1 (M ∼ 1.124, see Fig. 6 in [37]). To prevent the collapse, the star can lose energy through two mechanisms: gravitational cooling (ejecting Proca particles) and gravitational-wave emission (which also carries angular momentum). One may ask if there is formation of a rotating Proca star after the merger; we see that if we increase the initial velocity the angular momentum is rapidly lost and goes below the energy of the Proca field. This result could be due to the constraint-violating initial data we use; thus, it is interesting to revisit this problem once constraint-satisfying initial data becomes available, as it would provide a more reliable answer. For very large initial velocities, the final star can be completely dispersed away. For the largest velocity, we see that the amplitude of the scalar potential is decreasing with time.

Waveforms
In Figs. 18-21 we exhibit the gravitational waveforms produced in the orbital mergers, showing again the l ¼ 2, m ¼ f0; þ2g modes, extracted at two radii, r ext ¼ f100; 120g. The waveforms, conveniently shifted and rescaled, overlap in the wave zone. For the PS1-PS1 merger one observes that, as one increases the initial velocity, the non-axisymmetric m ¼ 2 mode grows visibly more than the axisymmetric m ¼ 0 mode. In fact the same trend occurs for the PS2-PS2 and PS3-PS3 mergers, but it is less pronounced. Therefore, the rescaling presented in Eq. (33) is no longer true. Moreover, the two parts of the waveforms shown in Figs. [18][19][20] seem to indicate that for models PS1-PS1 and PS3-PS3 a transient hypermassive PS forms before it collapses to a BH. Finally, for model PS00-PS00 (cf. 21) the gravitationalwave emission does not decay as the result of this merger is a highly perturbed PS. The waveform is filled with high frequency noise, probably coming from reflections with the outer boundary. The gravitational waveforms for this case are markedly similar to previous results in the scalar case [17,18,21].

VI. CONCLUSIONS
In this paper we have used numerical-relativity techniques to study both head-on collisions and orbital binary mergers of Proca stars of equal ADM mass and vector field frequency and we have extracted the gravitational waves produced in those collisions. This work continues our numerical exploration of the dynamics of PSs initiated in [38,43]. PSs are macroscopic, self-gravitating, Bose-Einstein condensates built out of a massive, complex, vector field [26]. Since they can achieve a compactness comparable to that of BHs, they are a type of BH mimicker with appealing dynamical properties, namely they are stable against perturbations [26,38] and they can form dynamically through a gravitational cooling mechanism [43].
Our investigation shows that the head-on collision of these spherically symmetric solutions may lead either to the formation of a more massive PS, which we dub "hypermassive" PS or, if the initial PSs are sufficiently massive/compact, to the formation of a Schwarzschild BH. Horizon formation, however, only occurs after an intermediate phase, which leaves an imprint in the waveform, making it distinct from that of a head-on collision of Schwarzschild BHs. After horizon formation the BH QNMs match those of a Schwarzschild BH. However, we have found that for those cases where the final BH is surrounded by a sufficiently extended and long-lived   quasistationary Proca cloud, differences with the BH ringdown are noticeable.
In the orbital binary case, we have also observed two fates: (i) the formation of a Proca star remnant, initially out of equilibrium and with angular momentum, but loosing such angular momentum as it approaches equilibrium; (ii) the formation of a Kerr BH, also initially surrounded by a cloud of quasi-bound states of the Proca field, which also tends to be absorbed/scattered during the evolution. We have not found evidence for the formation of either rotating Proca stars or infinitely long-lived Proca hair around a rotating horizon. In the case of the rotating Proca stars, since these are thought to be perturbatively stable in some region of parameter space (see [25] for a discussion in the analogue case of rotating boson stars), it may be that the island of initial conditions leading to their formation has not yet been scanned by our simulations. In the near future we plan to scan the space of initial data, in particular initial velocities. For larger velocities, obtaining reliable results requires constraint-preserving initial data. Obtaining such data and using it for performing further numerical evolutions is work underway.

APPENDIX: CODE ASSESSMENT
We briefly comment here on the standard analyses we carried out to assess the quality of our simulations with the EINSTEIN TOOLKIT. In Fig. 22 we show the behavior of the L 2 -norm of the Hamiltonian constraint for the head-on case. The floor of the violation of this constraint is at ∼4 × 10 −4 during most of the simulation. A significant peak in the violations of the constraint appears when the BH forms, but the errors quickly decrease to the precollision values and remain stable during the rest of the simulation.
Finally, in Figs. 24 and 25 we plot the Hamiltonian constraint of the head-on collision of the PS2-PS2 model at t ¼ 0, t ¼ 125, and t ¼ 500 and of the orbital merger at t ¼ 0, t ¼ 128, and t ¼ 304. The analysis yields a similar result for both types of collisions. At t ¼ 0 the constraint violations do not converge with resolution (our initial data do not satisfy the constraints). At intermediate times we obtain first-order convergence and, when the BH forms, a convergence of around second order is achieved. Therefore, we conclude that during the simulation we obtain between first-and second-order convergence. We note that the convergence order is greatly influenced by the linear second-order interpolation from the spherical grid of the Proca star solutions to the Cartesian grid we use to perform the evolutions with the EINSTEIN TOOLKIT.