Gravitational wave recoils in non-axisymmetric Robinson–Trautman spacetimes

We examine the gravitational wave recoil and the associated kick velocities in non-axisymmetric Robinson–Trautman (RT) spacetimes. Characteristic initial data used for the dynamics correspond to non-head-on collisions of two black holes. We make a parameter study of the kick distributions for an extended range of the incidence angle ρ0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _0$$\end{document} in the initial data. In the range examined 3∘≤ρ0≤125∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3^{\circ }\le \rho _0 \le 125^{\circ }$$\end{document} the kick distribution Vk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_k$$\end{document} as a function of the symmetric mass η\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} satisfies an empirically modified Fitchett law, with a parameter C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C$$\end{document} that accounts for the nonzero net gravitational wave momentum flux in the equal-mass case. The law fits accurately the numerical data with a normalized rms error ≤0.3%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\le } 0.3~\%$$\end{document}. The maximum kick velocity is ≃190km/s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\simeq }190~\hbox {km}/\hbox {s}$$\end{document} for ρ0≃55∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _0 \simeq \! 55^{\circ }$$\end{document}. A recent modification of the Fitchett law motivated by the effective-one-body formalism (Nagar in Phys Rev D 88:121501R, 2013) is also examined and, with the needed introduction of C\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C$$\end{document} to account for non-head-on collisions, fits accurately the RT data with a normalized rms error ≤3×10-5%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\le }3\times 10^{-5}~\%$$\end{document}. We construct the surface Vk(η,ρ0)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_k(\eta ,\rho _0)$$\end{document} in the parameter space of RT initial data, which gives an overall view of the behavior of Vk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_k$$\end{document} as the parameters change. The angular patterns of the gravitational waves emitted are analyzed and include the two polarization modes present in the radiative field of a non-head-on collision.


Introduction
The collision and merger of two black holes is presently considered to be an important astrophysical configuration where processes of generation and emission of gravitational waves take place (cf. [1,2] and references therein). The radiative transfer involved in these processes, evaluated in the full nonlinear regime of general relativity, shows that gravitaa e-mail: rfaranha@cbpf.br b e-mail: ivano@cbpf.br c e-mail: tonini@cefetes.br tional waves extract mass, momentum, and angular momentum of the source, and may turn out to be fundamental for the astrophysics of the collapse of stars and the formation of black holes. The process of momentum extraction and the associated recoils in the system can have important consequences for astrophysical scenarios, as the evolution and the population of massive black holes in galaxies or in the intergalactic medium [3][4][5]. Observational evidence of black hole recoils have been reported in [6,7] and references therein.
Gravitational wave recoils and the associated kick processes of two black holes have been investigated within several approaches, most of them connected to binary black hole in-spirals. Post-Newtonian approximations (cf. [8] and references therein) estimated the kick velocity accumulated during the adiabatical in-spiral of the system plus the kick velocity accumulated during the plunge phase. Sopuerta et al. [9] computed the recoil velocity based on the close limit approximation (CLA) supplemented with post-Newtonian (PN) calculations. The first full numerical relativity (NR) evaluation of the recoil in nonspinning black hole binaries was reported by Baker et al. [3] for a mass ratio 0.667, while González et al. [10] and Campanelli et al. [11] simultaneously obtained much larger recoils for black hole binaries with antialigned spins. González et al. [12] undertook a more complete NR treatment of kicks in the merger of black hole binaries by contemplating a larger parameter domain. For the case of small mass ratios in the interval 0.01 ≤ α ≤ 0.1 full numerical relativity evaluations bridged with perturbative techniques were implemented in Refs. [13][14][15][16]. Le Tiec et al. [17], combining PN+CLA methods, recently evaluated the gravitational wave recoil in black hole binaries and showed that the ringdown phase produces a significant antikick. In the same vein Choi et al. [18] examined recoils in head-on collisions of black holes, considering the head-on case as a model problem which can be seen as an approximation to the final plunge to merger and allow one to iso-late kick effects from the orbital in-spiral motion. Finally Rezzolla et al. [19,20] obtained an important injective relation between the kick velocities and the effective curvature parameter of the global apparent horizon in head-on collisions, using initial data derived in [21,22]. In spite of the enormous progress achieved until now using approximation methods and numerical techniques, the information on wave form patterns and radiative transfer processes in the dynamics of gravitational wave emission is far from being complete [23].
In this paper we examine the distribution of kicks for a large domain of the collision angle ρ 0 , a parameter of the initial data for non-head-on collisions of two black holes used in our simulations. The work completes Ref. [24] where we examined the energy and momentum extraction in the post-merger phase of non-head-on collisions, in the realm of Robinson-Trautman (RT) spacetimes [25,26]. Our treatment is based on the Bondi-Sachs (BS) conservation laws [27][28][29][30][31] that regulate the radiative transfer processes involved in the emission of gravitational waves. The characteristic initial data constructed for the dynamics already present a global apparent horizon so that the dynamics covers the post-merger phase of the system up to the final configuration of the remnant black hole. The novel feature of this topic is connected to the presence of a nonzero net gravitational wave momentum flux in the non-head-on case, absent in the cases of nonspinning black hole binary in-spirals or head-on collisions examined in the literature of numerical relativity simulations. The resulting kick distributions, for relatively large angles of collision ρ 0 , have a monotonous behavior with the symmetric mass-ratio parameter, also not seen yet in the literature.
We also examine the angular patterns of the gravitational waves. Our analysis includes the two polarization modes present in the radiative field of a non-head-on collision. Previous treatments of the angular patterns in the emission of gravitational waves-connected to the luminosity and the angular distribution of the emitted energy (antenna pattern)-were done in Refs. [32][33][34][35][36][37] where the authors used a post-Newtonian approach and showed the universality of gravitational bremsstrahlung in the low deflection encounters of two massive bodies. However, the use of black hole initial conditions and distinct incidence angles were not contemplated due to the approximations used. In Ref. [21] angular patterns of gravitational waves emitted in small mass-ratio head-on collisions of two black holes were examined, presenting short bursts of gravitational bremsstrahlung. Also Ref. [38] addressed the angular pattern of the same system examined here; however, the analysis there is incomplete, since its authors did not realize the presence of the additional polarization mode B in the wave zone curvature of non-axisymmetric Robinson-Trautman spacetimes.
Due to the presence of a global apparent horizon the initial data effectively represents an initial single distorted black hole which is evolved via the RT dynamics. Similarly to the case of the CLA-where the perturbation equations of a black hole [39][40][41] are feeded either with numerically generated, or with Misner, or Bowen-York-type initial datawe feed the (nonlinear) RT equation with the above mentioned characteristic data. It is in this sense that we denote the dynamics thus generated as "the post-merger phase of two colliding black holes". The interpretation of the outcomes of the RT dynamics should be considered with the above caveats.
In the next section we give a brief review of the basic properties of RT spacetimes which will be necessary for the discussions in the paper.

Robinson-Trautman spacetimes
RT spacetimes [25,26] are asymptotically flat solutions of Einstein's vacuum equations that describe the exterior gravitational field of a bounded system radiating gravitational waves. In a suitable coordinate system the metric can be expressed as where The Einstein vacuum equations for (1) result in Subscripts u, θ and φ, preceded by a comma, denote derivatives with respect to u, θ and φ, respectively. m 0 > 0 is the only dimensional parameter of the geometry, which fixes the mass and length scales of the spacetime. Equation (3), the RT equation, governs the dynamics of the system and allows one to evolve the initial data K (u 0 , θ, φ), given in the characteristic surface u = u 0 , for times u > u 0 . For sufficiently regular initial data RT spacetimes exist globally for all positive u and converge asymptotically to the Schwarzschild metric as u → ∞ [42,43]. Once the initial data K (u 0 , θ, φ) is specified, a unique apparent horizon (AH) solution is fixed for that u 0 [44,45]. Since the AH is the outer past marginally trapped surface, the closest of a white hole definition (the remnant black hole will form as u → ∞), only the exterior and its future development via RT dynamics with outgoing gravitational waves is of interest. We note that all the BS quantities, measured at the future null infinity J + , are constructed and well defined under the outgoing radiation condition [19,20,31]. The field equations have a stationary solution that will play an important role in our discussions, wherex = (sin θ cos φ, sin θ sin φ, cos θ) is the unit vector along an arbitrary direction x, and n = (n 1 , n 2 , n 3 ) is a constant unit vector (satisfying n 2 1 + n 2 2 + n 2 3 = 1). Also K 0 and γ are constants. We note that (4) yields λ = 1/K 2 0 , showing its stationary character. This solution can be interpreted [27] as a boosted black hole along the axis determined by the unit vector n with boost parameter γ , or equivalently, with velocity parameter v = tanh γ . The Bondi mass function associated with (4) is m(θ, φ) = m 0 K 3 (θ, φ) and the total mass-energy of this gravitational configuration is given by the Bondi mass The interpretation of (4) as a boosted black hole is relative to the asymptotic Lorentz frame which is the rest frame of the black hole when γ = 0. In the paper we use units such that 8π G = c = 1; c is however restored in the definition of the kick velocity. All the numerical results were done for the boost parameter γ = 0.5. In our computational work we used m 0 = 10 but the results are given in terms of u/m 0 . We should note that we can always set m 0 = 1 in the RT equation (3) by the transformation u →ũ = u/m 0 .

The Bondi-Sachs four-momentum for RT spacetimes and the initial data
Since RT spacetimes describe asymptotically flat radiating spacetimes and the initial data of its dynamics are prescribed on null characteristic surfaces, they are in the realm of the 2+2 Bondi-Sachs formulation of gravitational waves in general relativity [27][28][29][30]. Consequently we must use suitable physical quantities of this formulation appearing in the description of gravitational wave emission processes, as the BS fourmomentum and its conservation laws. A detailed derivation of the BS four-momentum conservation laws in RT spacetimes was given in [31]. We can show that, from the supplementary vacuum Einstein equations in the B-S integration scheme together with the outgoing radiation condition, the B-S four-momentum conservation laws for RT spacetimes are where is the Bondi-Sachs four-momentum and where m(u, θ, φ) is the Bondi mass function. In the above the four vector l μ = (1, − sin θ cos φ, − sin θ sin φ, − cos θ), defined in an asymptotic Lorentz frame, characterizes the generators l μ ∂/∂U of the BMS translations in the temporal and Cartesian x, y, z directions of the asymptotic Lorentz frame [30], and is the net flux of energy-momentum carried out by the gravitational waves. In (8), the quantities c (1) ,u and c (2) ,u are the news functions for RT spacetimes expressed as c (1) ,u (u, θ, φ) = where we have introduced the variable P ≡ 1/K , for notation convenience. We remark that c (1) ,u = 0 = c (2) ,u for the boosted Schwarzschild solution (4), as should be expected. The mass-energy conservation law (Eq. (6) for μ = 0) is the Bondi mass formula. Our main interest here is the analysis of the momentum conservation, (Eq. (6), for μ = x, y, z. Due to the planar nature of a general collision, namely, the motion of the two initial colliding black holes and the motion of the remnant are restricted to a plane, without loss of generality we will fix this plane as the (x, z)-plane so that the momentum conservation equations relevant to our discussion reduce to with wheren = (sin θ cos φ, 0, cos θ). Obviously P y is conserved, a consequence of P y W (u) = 0 for all u. The initial data to be used was derived in Ref. [24] and is interpreted as representing two instantaneously colliding Schwarzschild black holes in the (x, z) plane, at u = u 0 , In the derivation of (12) it turns out that α 2 /α 1 is the mass ratio of the Schwarzschild mass of the initial data, as seen by an asymptotic observer. It is also worth mentioning the following properties of (12): (i) for α 2 = 0 (α 1 = 0) or α 1 = 0 (α 2 = 0) the initial data (12) corresponds to a boosted Schwarzschild black hole along, respectively, the positive z-axis and along the direction defined by unit vector n = (− cos ρ 0 , 0, sin ρ 0 ) with respect to an asymptotic Lorentz frame (cf. Eq. (4)). (ii) The specific combination (12) of two black hole solutions (4) is not arbitrary but arises as the conformal factor of an asymptotically flat 3-geometry which is a solution of the constraint (3) R = 0. (iii) Additionally (12) results in a planar dynamics, namely, for all u the net gravitational wave momentum flux P W (u) is restricted to the plane determined by the unit vectors n z and n defining the direction of motion of the two black holes entering in (12); the momentum of the remnant black hole is also contained in this plane, as should be expected. The above properties reinforces the interpretation of (12) as related to the postmerger phase of a black hole collision. The parameter ρ 0 of the initial data, which we denote the incidence angle, defines the direction of the second initial black hole with respect to the z-axis of the asymptotic Lorentz observer. γ is the boost parameter of the black hole solutions (4) which enter in (12). As in the 1+3 numerical relativity approach the interpretation of the initial data parameters involves an approximation, namely, that the initial gravitational interaction in the data is neglected.
As mentioned already this data has a single apparent horizon so that the evolution covers the post-merger regime up to the final configuration, when the gravitational wave emission ceases. It is worth remarking here that, in the full Bondi-Sachs problem, further data (the news functions) are needed to determine the evolution of the system. However for the RT dynamics the news are specified once K (u, θ, φ) is given, cf. (9).

Numerical evolution
The initial data (12) is evolved numerically via the RT equation (3), which is integrated using a Galerkin method with a spherical harmonics projection basis space [46] adapted to the non-axisymmetric dynamics of RT spacetimes. The implementation of the Galerkin method, as well as its accuracy and stability for long time runs, is described in detail in Section V of Ref. [24]. The autonomous dynamical system derived with the Galerkin basis projection is integrated using a fourth-order Runge-Kutta recursive method (adapted to our constraints) together with a C++ integrator [24] for a truncation N = 7. Exhaustive numerical experiments show that after a sufficiently long time u ∼ u f all the modal coefficients of the Galerkin expansion become constant up to 12 significant digits, corresponding to the final time of computation u f . At u f the gravitational wave emission is considered to effectively cease. By reconstructing numerically K (u, θ, φ) for all u > u 0 we can obtain the time behavior of important physical quantities, as for instance the net gravitational wave flux and the associated total impulse imparted to the merged system by the emission of gravitational waves. From the final constant modal coefficients we obtain K (u f , θ, φ) that, in all cases, can be approximated as With the final parameters The values of the parameters of the remnant black hole are one of the basic results to be extracted from our numerical experiments, and are included in the tables of the next sections.

Gravitational wave net momentum fluxes and kicks in a non-head-on collision
We can now examine the processes of momentum extraction and the associated impulses imparted to the merged system by the emission of gravitational waves. Our starting point is the construction, of the curves of the net momentum fluxes carried out by gravitational waves, via the numerically integrated function K (u, θ, φ). Our numerical work in the present paper contemplates the parameter intervals where is the impulse imparted to the merged system due to the momentum carried out by the gravitational waves emitted up to the time u, withn = (sin θ cos φ, 0, cos θ). In Fig. 1 (left) we show the curves of the net momentum fluxes P x W (u) and P z W (u) for the mass ratio α = 0.25 and the incidence angle ρ 0 = 60 • , for u > u 0 . For this value of ρ 0 the net momentum flux is negative for all u, corresponding to a strong deceleration regime of the system by the emission of gravitational waves up to the final configuration of the remnant black hole, when the gravitational wave emission ceases. We can also see that these fluxes correspond to short pulses of gravitational bremsstrahlung in an interval Δu/m 0 ∼ 10. For incidence angles smaller than 55 • (and γ = 0.5) the net momentum flux P z W (u) is always positive for a short initial period.
The behavior of the associated impulses is illustrated in Fig. 1 (right). As expected I z W (u) and I x W (u) are negative for all u and tend to a constant negative value (a plateau) for large u ∼ u f corresponding to the final configuration of the system, that of the remnant black hole. The plateau is considered to be reached when |I W (u) − I W (u + h)| 10 −10 , where h is the stepsize of the integration used for the evaluation of I W (u). At this stage the remnant black hole has a momentum P = (n 1 f , 0, n 3 f ) P f , with The numerical tables include values of the parameters characterizing the remnant black hole for several ρ 0 . Typically the net total impulse imparted to the system has a dominant contribution from the deceleration regimes (where P z W (u) < 0 and P x W (u) < 0) and will correspond to a net kick on the merged system. As we will discuss later this net total impulse corresponds to the momentum of the remnant in a zeroinitial-Bondi-momentum frame.
From Eq. (14) we derive where the right-hand side of (17) are the nonzero components of the net total impulse I W (u f ) generated by the gravitational waves emitted. The values of I W (u f ) correspond to the final plateau which are present the impulse curves for any value of the initial data parameters, as illustrated for instance in Fig. 1 (right). We define the net kick velocity V k as proportional to the net momentum imparted to the system by the total impulse of the gravitational waves. This definition is based on the impulse function I W (u) evaluated at u = u f (cf. Eqs. (17)) and are in accordance with [12]. We obtain (restoring universal constants) with modulus where m 0 K 3 f is the rest mass of the remnant black hole. Taking into account the momentum conservation equations evaluated at u = u f we interpret (18) as the balance between the Bondi momentum of the system and the impulse of the gravitational waves in a zero-initial-Bondi-momentum frame, which can then be compared with the results of the literature. We remark that the zero-initial-Bondi-momentum frame is the inertial frame related to the asymptotic Lorentz Table 1 Summary of our numerical results corresponding to an incidence angle ρ 0 = 21 • and boost parameter γ = 0.5  Table 2 Summary of our numerical results corresponding an incidence angle ρ 0 = 60 • and boost parameter γ = 0.5 frame used in our computations by a velocity transformation with velocity parameter v B = P(u 0 )/m 0 K 3 f ; in the parameter domain of our numerical experiments the relativistic corrections in this transformation may be neglected. We note that the velocity (18) is directed along an axis making the angle Θ f = arctan(I x W (u f )/I z W (u f )) with the negative z-axis of the zero-initial-Bondi-momentum frame.
For our initial data (12) we have numerically evaluated V k contemplating an extended range of the parameters α and ρ 0 , with fixed γ = 0.5. The numerical results, illustrated in Tables 1 and 2 Fig. 2, for two separate domains of the incidence angle ρ 0 , the first corresponding to a domain of ρ 0 for which V k for α = 1 increases with ρ 0 , and the latter for which V k for α = 1 decreases with the increase of ρ 0 . The threshold between the two behaviors is ρ 0 55 • . The cases ρ 0 = 110 • , 115 • , 125 • were not included in the figures to avoid overcluttering, presenting a similar monotonous increase in η. The continuous curves are the best fit of the points to the empirical analytical formula with best-fit parameters given in Table 3. Equation (20) is an empirical modification of the Fitchett law [8,47], where  (20), with best-fit parameters given in Table 3   Table 3 Data for gravitational wave recoil of the equal-mass case (α = 1), and the best-fit parameters for the empirical law V = Aη 2 (1 − 4Cη) 1/2 (1 + Bη) × 10 3 km/s. For the case α = 1 we have the exact relation ρ f = (180 the additional parameter C was empirically introduced to account for the nonzero net gravitational wave momentum flux in the equal-mass case α = 1. (20) reduces to the Fitchett law for C = 1. The Fitchett law was derived from post-Newtonian analysis and used by a number of authors [12,17,47] to adjust the distribution of kick velocities in numerical relativity evaluations of the gravitational wave recoil in merging binary in-spirals of black holes and consistently yields a zero result for the equal-mass case. Therefore the results for kick distributions in the non-head-on case have no connection with black hole binary in-spirals, but rather possibly with two colliding black holes in premerger unbounded trajectories, or in hyperbolic encounters of two nonspinning black holes followed by a merger, the lat-ter configuration recently discussed by Gold and Brügmann [48]. We mention that the modification introduced in (20) is the only one that works to produce an accurate fit of our results, with normalized rms error of the order of, or smaller that 0.3 %, for the whole range of ρ 0 considered, this error decreasing as ρ 0 increases. Recently Nagar [49], motivated by the effective-one-body formalism, suggested a modification of the Fitchet law for binary black holes of the form where f (η) = (1 + B 1 η + B 2 η 2 + B 3 η 3 + B 4 η 4 ) and C = 1.
We verified that such a model, with the introduction of the factor C = 1 needed to account for non-head-on collisions, is also able to fit the RT data. In fact (21) yields a very accurate fit for RT data with a normalized rms error (between the curve and numerical points) of the order of, or smaller than 3 × 10 −5 % for the range of ρ 0 considered, this error decreasing as ρ 0 increases. We remark that the normalized rms error is four orders of magnitude smaller than the fit with the empirical law (20). The curves for each law, if plotted together, are indistinguishable on the scale of Fig. 2. We illustrate the case ρ 0 = 21 • for which the best-fit parameters for the RT data (cf. Table 1) In fact we verified that the upper bound of Δ(η) decreases monotonically with ρ 0 , from Δ ≤ 0.46 km/s for ρ 0 = 3 • to Δ ≤ 0.005 km/s for ρ 0 = 90 • . Consequently, with the introduction of the parameter C needed to account for non-head-on collisions, the Nagar distribution law (21) yields a very accurate function compatible to fit the RT data for kick velocities, as it does for numerical relativity simulations of binary black hole inspirals. We have a final comment that may be relevant for the present issue. As discussed in [49], the introduction of multipole terms in the wave form up to l = 8 is essential for obtaining (21). Now in our physical description with the RT dynamics the multipole expansion is not necessary since we have the exact expressions for the wave forms r Ψ 4 ∼ (D + i B) and equivalently for the total impulses I W (u f ), Eqs. (23) and (15) respectively. If we extend f (η) in (21) up to the sixth power in η we verify that the error of the fit of RT numerical points relative to the former fourth power fit is negligible, indicating a rapid convergence of f (η). In Fig. 3 we plot the parameters A, B and C (associated with the best fit of the law (20) to the numerical RT data) versus ρ 0 , as given in Table 3. The continuous curves are the best fit of the points through an eighth order polynomial least-squares method. A similar pattern is verified for the parameters A, B and C in the distribution (21). By using the best-fit curves A(ρ 0 ), B(ρ 0 ) and C(ρ 0 ) (cf. Fig. 3) and the kick velocity distributions we are able to construct the surface V k (η, ρ 0 ) in the parameter space of RT initial data for non-head-on collisions of black holes, as shown in Fig. 4. This gives us a global view of the behavior of V k as the massratio and incidence angle parameters vary, as for instance the absolute maximum of V k for (η = 0.25, ρ 0 55 • ).
The nonzero kick velocity for the non-head-on data with α = 1 deserves a further discussion that, without loss of generality, will be restricted to the case ρ 0 = 21 • . In this instance  (20), in terms of ρ 0 (cf. Table 3). The continuous curves are the best fit of the points through an eighth order polynomial least-squares method Fig. 4 Plot of the kick velocity as a surface V k (η, ρ 0 ) in the parameter space of RT initial data we can evaluate the components of the initial BS momentum to be, P x (0)/m 0 4.489093, P z (0)/m 0 0.832004 and P y (0)/m 0 0, with respect to an asymptotic Lorentz observer. This momentum vector, which lies in the right quadrant of the upper hemisphere z > 0 of the plane x − z, makes an angle Θ B = arctan |P x (0)/P z (0)| 1.387537 radians (or Θ B 79.5 • ) with the positive z-axis. This is also the direction of the nonzero momentum of the remnant with respect to the same asymptotic Lorentz frame, determined by the angle ρ f , which satisfies ρ f = (180 • − ρ 0 )/2 for α = 1 and any ρ 0 (cf. Table 3 and Ref. [24]). The axis determined by ρ f ≡ Θ B actually plays an important role in the dynamics. If we take a new frame with its z-axis coincid-ing with this axis the net gravitational wave momentum flux vector P W (u) lies along the new z-axis for all u. This was verified numerically by evaluating the ratios of the computed fluxes, sampled in the interval 0 < u/m 0 ≤ 890, yielding in all cases arctan |I x W (u)/I z W (u)| 1.387541 with a relative error of the order of 10 −6 . Still in this new frame the data will not be symmetric under θ → π − θ , leading to a nonzero net gravitational wave momentum flux, contrary to the case of merging binary in-spirals and head-on collisions. As expected the z-axis of the new frame is the direction of the kick velocity since arctan |I x W (u f )/I z W (u f )| 1.387547 rad or 79.5 • (within the precision of data in Table 1).
A remark is in order now concerning the balance between the total rest mass of the remnant and the total net impulse of the gravitational waves in the distributions of the net kick velocity, as observed from the numerical results displayed in the tables. In the domain 0 < ρ 0 < 55 • , as η increases from 0 to 0.25, both the parameter K f and the total net 2 increase; however, the increase of the rescaled Bondi rest mass of the remnant, K 3 f , is smaller than the increase of I W (u f ) up to η 0.225, implying that in this range the net kick velocity increases in accordance with (19). Beyond this point the increase of K 3 f is larger than the increase of the total net impulse leading to a decrease in the values of V k up to η = 0.25. On the other hand, in the domain 55 • < ρ 0 < 125 • , the above behavior is reversed leading to the monotonous increase in η shown in Fig. 2.
In general, in a zero-initial-Bondi-momentum frame, the Bondi momentum of the merged system satisfies P(u) = I W (u) so that an integral curve x(u) of the wave impulse vector field I W (u), defined as dx/du = P(u), can give a schematic picture of the motion of the merged system in this frame. In Fig. 5 we display this integral curve for α = 0.2 and ρ 0 = 21 • , generated with initial conditions x(u 0 ) = 0 = z(u 0 ) in the zero-initial-Bondi-momentum frame. The initial phase of positive momentum flux along z is responsible for the curved form of the trajectory in the semiplane z > 0.
For u → u f the curve approaches the asymptote with angle 15.11 • with respect to the negative z-axis of the zero-initial-Bondi-momentum frame, which is actually the direction of the kick velocity in this frame. In the case α = 1 the integral curve is a straight line with angle Θ f 79.5 • ≡ (180 • − 21 • )/2 with respect to the z-axis of the zero-initial-Bondi-momentum frame (cf. Table 3).

The angular wave pattern and the bremsstrahlung regime of the gravitational waves
The radiative character of RT spacetimes is given by the expression of its curvature tensor that in a suitable semi-null tetrad basis [24] assumes the form where the quantities N ABCD , III ABCD and II ABCD are of the algebraic type N , III and II, respectively, in the Petrov classification of the curvature tensor [50,51], and r is the distance parameter along the principal null direction ∂/∂r . Equation (22) displays the peeling property [52,53] of the curvature tensor, showing that indeed RT is the exterior gravitational field of a bounded source emitting gravitational waves. For large r we have so that at large r the gravitational field looks like a gravitational wave with propagation vector ∂/∂r . The nonvanishing of the N ABCD is therefore an invariant criterion for the presence of gravitational waves, and the asymptotic region where O(1/r )-terms are dominant defined as the wave zone. The curvature tensor components in the above basis that con- ,u (u, θ, φ) P , with the news c (1) ,u (u, θ, φ) and c (2) ,u (u, θ, φ) given in (9). From (23) we can see that the functions D and B contain all the information of the angular, and time dependence of the gravitational wave amplitudes in the wave zone, once K (u, θ, φ) is given. D and B actually correspond to the two polarization modes of the gravitational wave, transverse to the direction of propagation in the wave zone. We note that in the axisymmetric case B = 0, which is the case of head-on collisions. We will consider the particular combination where Ψ 4 is the Weyl spinor associated with N ABCD /r in a suitable Newman-Penrose null tetrad basis [53][54][55].
The quantity (25) is specified once we have the function K (u, θ, φ), which is numerically obtained via the numerical integration of the dynamics, as we have discussed. In Fig. 6 (left) we display the polar plots of √ D 2 + B 2 at early times u = 0.01 (dotted), u = 0.05 (dash-dotted) and u = 0.1 (continuous), and initial data parameters α = 0.2, ρ 0 = 55 • and γ = 0.5, with section by the plane φ = 0 • (corresponding to the plane of collision (x, z)). The plots show, for each time, a pattern with two dominant lobes in the forward direction of motion of the merged system. The direction of the Bondi momentum vector at u = 0.01 makes an angle Θ B 8, 43 • with the z axis. The pattern is typical of a bremsstrahlung process due to the deceleration of the merged system, analogous to the electromagnetic case of a charge decelerated along its direction of motion. As time increases we observe that the cone enveloping the dominant lobes opens up and the amplitudes decrease. For later times the pattern evolves to the expected quadrupole structure with a much smaller amplitude, as shown in Fig. 6 (right) for u = 5.0. We mention that the increase of the initial boost parameter γ would sharpen the forward cone enveloping of the two dominant lobes in the early regime, as expected in a In Fig. 7 we show the polar plot of √ D 2 + B 2 , with section by the plane φ = 90 • (corresponding to the plane (y, z), orthogonal to the plane of collision), at u = 0.1. The same initial data parameters of the previous Figures were used. As expected the pattern is symmetric about the z axis in accordance with the conservation of P y W (u) = 0. We then see that, although gravitational waves are emitted outside the plane of the collision, this radiation component has a zero net momentum flux; therefore it does not extract momentum of the system, consistent with the planar nature of the collision.

Conclusions and final comments
In the present paper we have examined the gravitational wave recoils and the associated net kick velocities in the postmerger phase of two black holes in non-head-on collision, by contemplating an extended domain of the incidence angle parameter ρ 0 present in the initial data. Our treatment is made in the realm of Robinson-Trautman dynamics and based on the Bondi-Sachs characteristic formulation of gravitational waves, and completes Ref. [24]. The kick velocity distributions are evaluated for an extended domain of the incident angle parameter, 3 • ≤ ρ 0 ≤ 125 • , of the initial data. We also obtain that in general the total net impulse is negative for all u which, from the Bondi-Sachs momentum conservation laws, corresponds to a dominant deceleration regime of the system due to the emission of gravitational waves. However, for relatively small values of the mass-ratio parameter α and of ρ 0 , an initial positive impulse in the z direction may be present, which will be responsible for an initial in-spiral branch in the motion of the system. By using the Bondi-Sachs conservation laws we evaluate the net kick velocity V k imparted to the system as proportional to the total net gravitational wave impulse in a zeroinitial-Bondi-momentum frame. For each value of the incidence angle ρ 0 considered the distributions of V k as a function of the symmetric mass parameter η are obtained.
The distributions of V k as a function of η were shown to be fitted, for the whole domain of ρ 0 considered, by the empirical law (20) obtained by a modification of the Fitchett ηscaling law, this modification corresponding to the introduction of an additional parameter C to account for the nonzero gravitational wave momentum fluxes in the equal-mass case (η = 0.25) of non-head-on collisions. The best fit of the points (V k , η) with the modified law is sufficiently accurate with a normalized rms error of the order of, or smaller than 0.3 % for all ρ 0 considered. For ρ 0 = 0 • (the case of headon collisions) we have C = 1 and this distribution reduces to the Fitchett law, as expected. For large incidence angles (e.g. ρ 0 > 55 • in the case of γ = 0.5) the distributions are monotonous in η. We also verified that the best-fit values of the parameter C decrease monotonically from C = 1 as the incidence angle ρ 0 increases. Finally we examined a kick distribution law recently proposed in [49] using the effective-one-body formalism (with the needed introduction of C to account for non-head-on collisions). We showed that this modified distribution yields a very accurate fit for RT numerical data, with a normalized rms error ≤3 × 10 −5 %, this error decreasing as ρ 0 increases. By using the best-fit curves A(ρ 0 ), B(ρ 0 ) and C(ρ 0 ) for the results in Table 3 and the kick velocity distributions we constructed the surface V k (η, ρ 0 ) in the parameter space of RT initial data for non-head-on collisions of black holes, giving an overall view of the behavior of V k as the parameters are changed.
A novel feature of non-head-on collisions (ρ 0 = 0) is a nonzero net gravitational wave flux for the equal-mass case, contrary to the cases of head-on collisions and merging of black hole in-spiral binaries. This implies that the net kick velocity for non-head-on collisions are nonzero for the equal-mass case. As a consequence the results for kick distributions in the non-head-on case have no connection with black hole binary in-spirals, but rather possibly with two colliding black holes in pre-merger unbounded trajectories, or in hyperbolic encounters of two black holes followed by a merger as recently discussed by Gold and Brügmann [48].
We examined the behavior of the integral curves of the gravitational wave impulse I W (u) that, in accordance with the Bondi-Sachs momentum conservation law, can describe the motion of the merged system in a zero-initial-Bondimomentum frame. This integral curve exhibits an initial inspiral branch in the positive z semiplane whenever an initial phase of P z W (u) > 0 is present. Finally we have examined the angular patterns of the radiation both in the initial regime, when the emission is typically bremsstrahlung, and in later times when the quadrupole pattern is already set up.