Radiation from a D-dimensional collision of shock waves: first order perturbation theory

We study the spacetime obtained by superimposing two equal Aichelburg-Sexl shock waves in D dimensions traveling, head-on, in opposite directions. Considering the collision in a boosted frame, one shock becomes stronger than the other, and a perturbative framework to compute the metric in the future of the collision is setup. The geometry is given, in first order perturbation theory, as an integral solution, in terms of initial data on the null surface where the strong shock has support. We then extract the radiation emitted in the collision by using a D-dimensional generalisation of the Landau-Lifschitz pseudo-tensor and compute the percentage of the initial centre of mass energy epsilon emitted as gravitational waves. In D=4 we find epsilon=25.0%, in agreement with the result of D'Eath and Payne. As D increases, this percentage increases monotonically, reaching 40.0% in D=10. Our result is always within the bound obtained from apparent horizons by Penrose, in D=4, yielding 29.3%, and Eardley and Giddings, in D>4, which also increases monotonically with dimension, reaching 41.2% in D=10. We also present the wave forms and provide a physical interpretation for the observed peaks, in terms of the null generators of the shocks.


Introduction
The ongoing physics runs at the Large Hadron Collider (LHC) with 7 TeV centre of mass energy, are starting to set bounds on physics beyond the standard model. For the particular case of TeV gravity models [1], microscopic black hole formation and evaporation is predicted for scattering experiments with partonic centre of mass energy beyond the TeV [2,3]. The first set of bounds for these models was recently released by the CMS [4] and ATLAS [6,7] collaborations. The current analysis for the 7 TeV data, is however extremely dependent on regions of parameter space where (at best) black holes with masses close to the unknown Planckian regime would be produced (see [8] for a discussion). Otherwise if strong cuts are imposed, such that the objects produced are in the semi-classical regime (which is the calculable one), the cross-sections become negligible at 7 TeV, and only after the upgrade of the beam energy to 14 TeV, planned to take place in 2013, will the scenario be properly tested. Any improvement in the phenomenology of these models is therefore quite timely.
Two event generators, charybdis2 [9] and blackmax [10] are being used at the LHC to look for signatures of black hole production and evaporation. The event rates they produce depend sensitively on various input parameters, amongst which two fundamental quantities for the modeling of the black hole production phase are: the critical impact parameter for black hole formation in parton-parton scattering; and the energy lost into gravitational radiation in this process. The latter, if large, could be an important signature for discovery or exclusion, if dominating the distribution of missing energy and if calculated with enough precision in the semi-classical regime. Furthermore, if the experimental searches continue to rely on events with little (or under-estimated) missing energy as in [4][5][6], then this information will be crucial to determine if there is enough phase space left with little missing energy.
As first argued by t'Hooft [11], classical gravity, described by general relativity, should dominate transplanckian scattering. Following this rationale, the best estimates, so far, for the two mentioned quantities, come from apparent horizon computations, initiated by Penrose in four spacetime dimensions. He found an apparent horizon on the past light cone of the collision between two shock waves, representing the gravitational field of two particles boosted to the speed of light. Thus, Penrose concluded that no more than 29.3% of the centre of mass energy could be emitted in gravitational radiation in a head-on-collision.
This analysis was generalised to D dimensions by Eardley and Giddings [16]. They obtained the bound on the percentage of the centre of mass energy, ǫ radiated , emitted in a head on collision, which increases with dimension and is given by where Ω n is the volume of the unit n-sphere. This bound increases monotonically with dimension approaching 50% in the limit of infinite D. Using the results in [16] for nonhead on collisions, the critical impact parameters were numerically evaluated by Nambu and Yoshino [17]. According to their estimate, a black hole will form if at least half of the effective gravitational radius of each of the colliding particles overlap. The state of the art results were obtained by Yoshino and Rychkov [20], who found an apparent horizon on the future light cone of the collision. Their analysis coincides with the one in [16] for the energy emitted in head on collisions; but the critical impact parameters for black hole formation become larger, yielding a larger cross section for black hole formation. Events generated by charybdis2 rely on a selection of random points from the configuration space allowed by the bounds determined in [20]. The event horizon, however, is located outside the apparent horizon, and the latter will evolve and settle to the former in the future of the collision. Thus, both the critical impact parameter and ǫ radiated should be exactly determined by knowledge of the metric in the future of the collision 1 . There are in principle two methods to compute it, both of which have been developed in four spacetime dimensions.
Modeling the gravitational field of the colliding partons as shock waves, one may use the Aichelburg-Sexl metric [15] to describe the process. Using a method first developed by D'Eath [18] and D'Eath and Payne [12][13][14], one may compute the metric to the future of the collision analytically, as a perturbative expansion. In four dimensions, second (first) order perturbation theory yields 16.3% (25.0%) for the energy emitted in gravitational waves in a head on collision.
Modeling the gravitational field of the colliding partons as highly boosted black holes, one may use numerical relativity techniques to perform fully non-linear numerical simulations of the collisions. In four dimensions the high energy collision of two black holes was performed by Sperhake et. al. [21], obtaining that 14±3% of the energy was converted into gravitational radiation. This overlaps with the result obtained by D'Eath and Payne in second order perturbation theory [14], which indicates that such an approximation provides a good estimate. Moreover, introducing an impact parameter [24], one observes that the critical impact parameter increases from about 0.8 in [20] to about 1.2, in units of the Schwarzschild radius of the centre of mass energy. This is an increase of about 50% which, if verified in the higher dimensional case as well, leads to significantly larger cross-sections, thus becoming easier to exclude regions in the model parameter space.
Phenomenologically interesting TeV gravity models occur in dimension D ≥ 6. It would therefore be useful to develop each of the two above methods in higher dimensions. Numerical relativity for higher dimensional spacetimes is being developed (see eg. [29,32]), and the first ever results for black hole collisions in higher dimensions (albeit low energy collisions) have been produced [30,31]. In this paper we shall perform a higher dimensional collision of shock waves and determine, for even D, the energy radiated by computing the geometry to the future of the collision to first order in perturbation theory, building up on the technique of [12]. The results we obtain for the radiated energy are summarised in the following table, where the apparent horizon bounds obtained in [16] are also given for comparison: Like the apparent horizon bound, our results indicate that the energy radiated increases monotonically with D. Moreover, our result is always below the bound, as expected. If this reduction trend is verified by higher order perturbation theory (or other methods, such as numerical black hole collisions) this result has phenomenological implications. More energy lost into gravitational radiation means a less massive final black hole is produced. Thus, this result suggests that the final black hole will be more massive than previously estimated, making it more consistent with the semi-classical analysis used for estimating the potentially observable Hawking radiation. This paper is organised as follows. In Section 2 we discuss the spacetime geometry of one and two shock waves (before the collision) in both Brinkmann and Rosen coordinates. In Section 3 the perturbative computation is setup and the integral solution for the metric in the future of the collision is provided. Some details concerning the derivation of this solution are provided in Appendix A. In Section 4 the extraction of gravitational radiation is studied. A formula in D dimensions is first derived in terms of the metric perturbation in the future of the collision starting from the Landau-Lifschitz pseudo tensor (in D dimen-sions). Then, using results from Appendix B, a workable expression is obtained to extract the energy carried out by the gravitational radiation (formulas (4.7) and (4.8)). These expressions are evaluated numerically to yield the results presented in the table above. The wave forms are exhibited and an interpretation is given for the domain where they have support and peak. We also discuss the approximations used in our derivation and further comment on our results in Section 5.

One D-dimensional shock wave in Rosen form
The D-dimensional Tangherlini solution with mass M is [25] and Ω D−2 are the line element and volume of the unit D − 2 sphere. The Aichelburg-Sexl solution [15] is found by boosting this black hole and then taking simultaneously the limit of infinite boost and vanishing mass, keeping the total energy µ fixed. We use a coordinate system (u, v, x i ), where the retarded and advanced times (u, v) are (t − z, t + z) in terms of Minkowski coordinates, and x i are the remaining Cartesian coordinates on the plane of the shock, i = 1 . . . D − 2, such that the transverse radius is ρ = x i x i . The resulting geometry for a particle moving in the +z direction is The function Φ depends only on ρ and takes the form [16] (2. 3) The coordinates used in (2.2) are of Brinkmann type [26]. Presenting the shock wave in this chart, geodesics and their tangent vectors appear discontinuous across the shock. One can, however, introduce a new coordinate system defined by [16] where θ is the Heaviside step function and Φ and its derivative Φ ′ are evaluated atρ. In this new chart both geodesics and their tangents are continuous across the shock atū = u = 0; φ a are the angles on the (D − 3)-sphere and a = 1...D − 3. Using (2.4) the metric becomes [18][19][20]  In the centre of mass frame, the radiation is symmetric under z → −z. Only the radiation that propagates along a small solid angle around θ = π in the centre of mass frame (red area, delimited by the polar angle tan θ = β −2 − 1) is still propagating in the −z direction in the boosted frame. All the radiation in this solid angle is redshifted in the boosted frame.
which we dub the Rosen form [27] of the shock wave. The geometry for an identical shock wave traveling in the −z direction is obtained by changing z → −z in (2.2) or equivalently exchangev ↔ū in (2.5). Following [12] it is simple to see that in a boosted frame (moving with respect to the (u, v) chart with velocity β in the −z direction), the oppositely directed shock waves keep their form, but with new energy parameters, respectively, where e α = (1 + β)/(1 − β). By causality, the spacetime where two shock waves travel in opposite directions along the z coordinate, is described everywhere by superimposing the two shock wave metrics, except in the future light cone of the collision. Thus, in the boosted frame the spacetime geometry reads which is valid everywhere except in the future light cone ofū =v = 0. As a consequence of the collision a strong burst of gravitational radiation (followed by a tail) is expected to be emitted. An illustration of this signal in the centre of mass/energy and in the boosted frame is presented in Fig. 1.

The perturbative metric in the future of the collision
We shall now compute the geometry in the future of the collision to first order in perturbation theory. The main idea is that in the boosted frame one shock wave can carry much more energy than the other, since Thus, one may face the wave traveling in the −z direction as a small perturbation of the wave traveling in the +z direction. Since the geometry of the latter is flat forū > 0, we make a perturbative expansion of the Einstein equations around flat space, in order to compute the metric in the future of the collision (v,ū > 0). To first order this amounts to solving the linearised Einstein equations around flat Minkowski spacetime, subject to a boundary condition given by the metric (2.7) in the limit u = 0 + . The computation will be performed in the more intuitive Brinkmann type coordinates. To obtain the boundary condition we transform back to these coordinates on u = 0 + . Using (2.4) with κ replaced by ν we find whereḡ αβ is the metric for the weak shock by takingū = 0 in (2.7). Observe that the dimension of the parameters ν and λ is [Length] D−3 , so we choose to perform the following rescaling to dimensionless coordinates (u, v, Performing the matrix multiplication and expressing everything in the unbarred dimensionless coordinates, we write the geometry, on u = 0 + , as where the flat metric is The perturbations 2 are Observe that the step function jumps at v = Φ/ √ 2. To understand why, recall that, from (2.7), the collision occurs atū = 0 =v, in Rosen coordinates. From (2.4), taking into account the rescaling to dimensionless coordinates performed above, the collision takes place at in Brinkmann coordinates. The future light cone of the collision has two branches:ū = 0,v > 0 andū > 0,v = 0. In terms of the Brinkmann coordinates these conditions read: In order to understand the interaction, let us follow the null generators of the shock with support atv = 0 (which we call the weak shock). Before the collision these can be parametrised as To go beyond the collision we look at Brinkmann coordinates. From (2.4): In terms of these coordinates the focusing effect on the weak shock null generators as they cross u = 0 is clear - Fig. 2, 3 and 4. Observe that there is a qualitative difference between the behaviour in D = 4 and D > 4: the discontinuous jump in v of the weak shock null generators is always positive in D > 4 but it becomes negative for large ρ in D = 4. This will give rise to a different domain in the time integration when computing the overall energy emitted. In all cases, nevertheless, the generators will focus, generating a caustic, at √ 2Λ = ξ D−2 , corresponding to After crossing at the caustic the generators enter the curved region of the spacetime. Thus, until they cross, the light cone structure is that of the Minkowski spacetime (3.3). It is simple to see that the deflection angle α undergone by a null shock generator at distance ρ from the axis obeys, in Brinkmann coordinates, Thus the focusing increases (decreases) with D for short (long) distances, as expected from the behaviour of the gravitational force. u v x weak shock null generators u = 0 strong shock Figure 3. Evolution of the weak shock null generators (blue arrows) from the viewpoint of Brinkmann coordinates in the boosted frame in D > 4. For u < 0 they are at v = 0; then the generators undergo a discontinuity in v at u = 0, which is ρ dependent but always positive. They jump to the collision surface (green lines, at , gain shear and focus along the caustic (red line).

Future development of the metric
To the future of u = 0 the metric will have a perturbative expansion of the type  After this event, which is at the caustic, spacetime becomes curved therein, and therefore the continuation of the rays beyond the axis is merely illustrative and is represented as a dotted line. For points outside the axis, such as point P, this diagram suggests an initial radiation signal, associated to ray 1, followed by a later burst, associated to ray 2. We shall see this interpretation matches the wave forms computed below. and we should be able to find h (i) µν by successively solving the Einstein equations at each order in λ/ν, with the boundary condition expressed before. To first order we therefore have to solve the linearised D-dimensional Einstein equations around Minkowski space, so in the remainder we consider the linear perturbation h µν in (3.10) and for notational simplicity drop the "(1)" label. We denote the trace reversed perturbation by barring ith µν = h µν − η µν h/2, and perform a coordinate transformation of the form This can be chosen such that the De Donder gauge condition is obeyed in u > 0 so that the linearised Einstein equations become simply a set of decoupled wave equations in Minkowski space for each component: In [12], and integral solution for the wave equation 2F = 0 in u ≥ 0 with initial data on u = 0 was used. In appendix A, we show that in D dimensions, the corresponding solution is where, for each x ′ , v ′ defines points at the u = 0 hypersurface, which are on the past light cone of the event (u, v, x i ) (cf. Fig. 5): 15) and the derivative operator is suitably defined in Fourier space (with respect to v) for odd dimensions (see appendix A). Similarly to [12], the gauge condition (3.12) is obeyed in To see this we first contract the wave equation with another derivative to get 2h N ,ν µν = 0. It then follows from the integral solution (3.14), that if then the De Donder condition holds everywhere. In particular (3.16) is sufficient, since it implies (3.17).

Gauge fixing at u = 0
Condition (3.16) can be written, in terms of the perturbations in the original coordinate system, using (3.11) (on u = 0). Then, using (3.13) and the fact that the trace of the first order metric perturbation h vanishes on u = 0, we find, that (3.16) implies, on u = 0 where the right hand side of (3.18) and (3.19) was computed using (3.4). We look for a solution of (3.18) and (3.19) which has a power series expansion around u = 0: One such solution is ξ Applying this gauge transformation, only the uu first order perturbation changes in this new (De Donder) gauge, on u = 0, which becomes Since the metric perturbation is traceless on u = 0, it is traceless everywhere, due to (3.14) The general solution of (3.23) in u > 0 is again obtained using (3.14) with the replacement F → h N µν . Observe that the components of the metric perturbation that vanish at u = 0, will vanish when u > 0. We can boost back to the centre of mass frame simply by using the replacements {u, v} → {e −α u, e α v}. At linear order, we can use x N µ = x µ for the coordinates in the new gauge. Finally, we can change from units of ν and obtain the linearised, centre of mass frame metric in u > 0

Extracting the gravitational radiation
To extract the gravitational radiation produced in the collision we shall use the Landau-Lifshitz pseudo-tensor [28], which was generalised to higher dimensions in [32]. In terms of our perturbations h N µν (u, v, x i ), taking into account that they are traceless, it reads (we omit from now on the superscript N for notational simplicity) [32]: Despite not being unique or gauge-invariant it is well known that the integral computed on a 'distant' surface with area element dS outward unit normal n i , is a gaugeinvariant well defined energy [33].
Following [12], we shall compute the power emitted along the vicinity of the θ = π direction (in the CM frame), corresponding to the negative z direction. This maps to the region where perturbation theory should be valid in the boosted frame (cf. Fig. 1). So we shall need the flux along the z direction which is given by ( 3) The first (second) term corresponds to the flux across a v = constant (u = constant) surface. Thus, in order to determine the flux in the θ = π direction we need only the second term. Then, from (4.1) (observe that the first two terms are zero by the De Donder gauge condition): The total radiation emitted will be computed, to first order, under the assumption that dE/d cos θ is isotropic (we shall further discuss this approximation below). Thus we integrate the power emitted inside the narrow cone around theθ ≡ π − θ = 0 axis. Using dS = r D−2 dΩ D−2 , taking the limit close to the axis and multiplying by the area of the sphere of radius r, this energy is As expected, only the ij components (which did not change with our gauge choice) determine the energy emitted in gravitational radiation. In appendix B we show that, in fact, there is only one independent quantity that determines the integrand in the last equation, denoted E = E(u, v, ρ). In terms of E, we can express the radiated energy as We identify, inside the parenthesis, the analogous of Bondi's news function used in [12]. The coordinates used in the last formula are dimensionful. The fraction of the energy radiated in gravitational waves ǫ radiated defined in (4.6) is given by, and, from appendix B where All coordinates in (4.7) and (4.8) are now the dimensionless coordinates used in Section 3.
Observe that at ρ = 0, the argument of the delta function has no x dependence and the x integral can be immediately performed to yield zero. This can be interpreted as the result of a destructive interference phenomenon. To see this consider for instance the h 11 component of the initial data (3.4), which can be written as The angular part of the perturbation is a scalar harmonic with ℓ = 2 (or ℓ = −2 in D = 4) on the (D − 3)-sphere. Integrating such scalar harmonic on the (D − 3)-sphere gives zero. This is the integral appearing in (4.8) at the axis and may be seen as a cancellation between the different phases coming from different points on the integration circle of radius ρ ′ . The vanishing of the power radiated at the axis in a head on collision of two equal objects (like particles or black holes) can also be physically interpreted as follows. Gravitational radiation emission is determined by the variation of the gravitational quadrupole. Taking the observation point at the collision axis, and behind one of the objects, no quadrupole variation is observed. 3

Integration limits
We now discuss the integration domain for the time integration in (4.7). Equation (4.8) suggests that the gravitational radiation seen at a spacetime point P to the future of the collision, originates from the points on the collision surface wherein the delta function has support; these are the points at the intersection of the past light cone of P with the collision surface, as expected. 4 Instead of using the spacetime coordinates (u, v, ρ) it is convenient to introduce a retarded time coordinate τ = t − r, where the radial coordinate is r = z 2 + ρ 2 and an angular coordinate θ by ρ = r sin θ. Thus we wish to know the first instant τ 1 for which the point specified by (r, θ) receives gravitational radiation. This event occurs when the intersection of its past light cone with u = 0 becomes tangent to the collision surface and is determined in a similar way to the determination of the caustic in Section 3. One finds that τ 1 = τ 1 (r, θ), together with the auxiliary variableρ, is obtained by solving the system of equations for specified ρ and θ and s = +1. The intersection of the past light cone of points (τ 1 , r, θ) (points 1 to 3 in Fig. 6) with u = 0 is plotted in Fig. 7, for θ = π and Fig. 8, outside the axis. As claimed these curves are tangent to the collision surface. 3 We thank U. Sperhake for this observation. 4 Observe, however, that an integration in v ′ was already performed to get (4.8). For even D this integration is performed using a delta function which enforces the support of the integrand on the light cone. But for odd D the initial integral has support also inside the light cone, cf. Appendix A. . This is tangent to the collision surface u = 0 and v = Φ/ √ 2, represented by the solid (green) lines. The generators of the shock traveling along u that emerge from the intersection points will focus and converge at 3. For an observer at fixed z, following the worldline represented by the τ axis, no gravitational radiation will be observed before τ = τ 1 . In Fig. 7  Thus the initial time for observing gravitational radiation, τ 1 , is positive, for D > 4 and goes to zero as r goes to infinity (cf. Fig. 3). By contrast, in D = 4, τ 1 becomes negative and approaches negative infinity as r goes to infinity (cf. Fig. 2). Thus, the time The previous computation of τ 1 amounts to computing the retarded time at which ray 1 intersects point P in Fig. 4. The retarded time, τ 2 , for which ray 2 intersects P is computed by solving (4.11) with s = −1. This coincides with a second peak in the wave forms exhibited in the next section.

Numerical evaluation
In this section we finally obtain the wave forms for several dimensions and integrate the radiated power to estimate the amount of gravitational radiation emitted in the collision. First we perform the angular integral in (4.8) (see appendix B where the polynomials of degree M + 2, P (M +2) and Q (M +2) are defined), to obtain in odd dimensions and in even dimensions. D is defined such that where U = τ + 2r sin 2 (θ/2), T = τ + 2r cos 2 (θ/2) − ρ 2 /U and we are now expressing the result in coordinates {τ, r, θ} as in Eq. (4.11).
In the remaining numerical analysis, we present only results for even D. Concerning odd D, we have attempted to evaluate (4.12); the wave forms obtained (numerically),  Figure 9. Wave forms: These plots show the even D wave forms (4.13) for fixed θ, close to the axis of symmetry and for various r. Note the initial step which is due to the impulsive nature of the initial conditions associated with the Aichelburg-Sexl shocks. This step is exactly at τ = τ 1 , whereas the sharp peak is at τ = τ 2 , for each curve. however, contained non-integrable tails when squared to get the radiated power. Our best hint, at the moment, is that these divergences are similar to those present in the second order computation [13]. Nevertheless, our aim is to investigate the variation of the radiated energy with D. Thus the even case already shows the general trend. It seems plausible that the result for D odd will interpolate between the even D cases we shall exhibit.
To evaluate (4.13), we wrote firstly a test code in mathematica7 and then a code in the c++ language, using the numerical integration routines of the Gnu Standard Library (GSL). The purpose was to check the two codes in some cases and use the (lower level) c++ code for generating all the data in a practical amount of time. The domain of integration D was determined by looking at the roots of two polynomials constructed from x ⋆ and using the root finding routines in mathematica or the GSL library. In general, depending on the spacetime point, we found either an empty integration domain (giving a value of exactly zero to the wave form), or a non-empty domain which can be simply connected or the union of two disconnected domains. The integrable singularities at the points x ⋆ = ±1 where removed explicitly through a change of variable around each singularity. In general, for the wave forms, we have demanded a relative error of 10 −8 , and 10 −5 for the final τ integral of the radiated power.
In Fig. 9, we present some wave forms (D even) as a function of τ . Since we are interested in the limit near the negative z axis (where our approximations are justified), we use examples with θ away from the axis by 0.01 radians. A first observation is that we have shifted τ by τ 1 for all wave forms such that zero corresponds to the beginning of the radiation burst. A consistency check is that, indeed, the numerically determined integration domain agrees exactly (within very small numerical errors) with this expectation. Thus all the wave forms have a sudden step at zero which arises naturally from the numerical code. Another feature which is verified numerically, is that the (integrable) singular peak of radiation observed in all plots, occurs exactly at τ = τ 2 . This is explicitly exhibited in Fig. 10, where besides the shift of τ by τ 1 we have rescaled the axis by ∆τ = τ 2 − τ 1 . An advantage of our numerical method compared to [12] is that no approximations were made in the integral. In fact in [12], the limit r → +∞ was taken inside the integral before actually performing the integration. This removes the ρ ′2 term in (4.14), which is responsible for excluding regions of the domain of integration with large ρ ′ . In four dimensions, due to the logarithmically growing Φ profile, it turns out that such regions of large ρ ′ are indeed cut off from the integral by the ln ρ ′ term. Thus in the large r limit, the wave form does approach the one provided in Fig. 4 of [12]. This can be seen in the first plot as r becomes large. However the approximation in [12] makes it less clear when the radiation burst starts and peaks. In our approach we simply use the full numerical integral without further approximations, for finite r, and find the limit numerically. Finally, for D > 4, the procedure in [12] is not valid, because it would amount to neglecting the only term growing with ρ ′ in the definition of x ⋆ , so there would be no regularisation of the domain and the integral would not converge. For D > 4, there is a notable difference as we take the large r limit, namely that the radiation pulse becomes more concentrated around τ = τ 2 . In fact, the same happens in four dimensions, however more slowly (as r increases) due to the logarithmic nature of the profile of the shock wave 5 . The other difference is that the wave form acquires one more oscillation in each jump to the next D even. Finally we note (without plotting explicitly) that the effect of decreasing θ (with r fixed), is the same as when we decrease r with θ fixed 6 , that is, the radiation burst becomes sharper and occurs in a smaller period ∼ ∆τ as we move closer to the axis.
After we obtain the wave forms for a certain r, θ, we can compute the estimate for the fraction of radiated energy, by computing the integral in (4.7) and taking the limit r → +∞,θ → 0. In Fig. 11 we show plots of the radiated fraction of energy as a function of r for several θ. We have performed a fit of each numerical curve to extract the limit, by  Figure 11. Limiting fractions: These plots show ǫ(r), whose limit for r large and small π − θ gives the estimate ǫ radiated . The best estimate for the limit can be read from the constant term in the asymptotic form used to fit the numerical points of the red curves (for which the angle is the smallest).
using an asymptotic expansion of the form ǫ(r) ≃ ǫ radiated + a/r (indicated next to each curve). Since we have computed this function for very large values of r we obtain a fit which is basically indistinguishable from the numerical tails, so ǫ radiated is extracted with a relative error of less than 0.001. Note that the difference in the shape of the curves for different θ, is again smaller for the D = 4 case due to the logarithmic dependence of the shock wave profile. Since we have used very small angles, all the asymptotic values in the different curves are consistent with the limiting values quoted in the introduction (within the relative precision of 0.001 and the error from the a/r correction).

Discussion of the method
The computation performed in this paper relies on perturbation theory. Why should perturbation theory be valid in describing a non-linear effect (like black hole formation) and an associated strong field effect, such as the generation of a strong burst of gravitational radiation? In [12] this is justified as follows. Mathematically, the boost provides a small parameter for expanding the geometry (the ratio of energies of weak to strong shock waves). Analysing (3.4), one concludes that in the boosted frame the initial conditions are indeed perturbative except in a small vicinity of ρ = 0 and sufficiently large v > 0. This has a physical picture. The null generators of the weak shock wave (in the boosted frame) are not only bent towards the axis of the collision, as they cross the strong shock, but they also suffer a redshift which becomes larger as the null generator approaches the axis and infinitely large for generators at the axis. This is manifest in Fig. 2 and 3. The radiation that hits a spacetime point (not at the axis), therefore, comes first from the far field region of the collision (ray 1 in Fig. 4), which is given a large head start as compared to its near field counterpart (ray 2 in Fig. 4). Since gravity is weak in the far field region, perturbation theory should therefore be accurate in describing this part of the signal. Moreover, due to focusing, the overall amount of radiation measured near the axis, at very large distances, coming from the far field region need not be small. Therefore, it is concluded in [12] that this problem provides one example where perturbation theory can tackle a strong field effect.
The interpretation of the wave forms we have provided in this paper shows, however, that the above argument must be considered with care. The wave forms exhibited in Fig.  9 and 10 are dominated by the signal near τ 2 , which is the signal that crosses the near field region. This indicates therefore, that there is a large contribution to the signal from outside the region where linear theory is clearly applicable and a correction from higher order terms in perturbation theory is expected. But the near field region experienced by ray 2 becomes smaller at very large r and smallθ. Furthermore, as we have learned in D = 4, the first order perturbation provides a better estimate than the trapped surface argument and, more importantly, the second order perturbation provides a result consistent with the numerical simulations. We find this an important motivation to pursue this computation. It will be interesting to see if the agreement between numerical results and perturbation theory still holds in higher dimensions.
In order to compute the total radiation emitted, we have determined the radiation in the vicinity of θ = π and extrapolated it under an isotropy assumption for dE/d cos θ, following [12]. This is partly justified by the Zero Frequency Limit approximation [34], which predicts an isotropic angular distribution of gravitational radiation (with two blind spots along the symmetry axis) for the collision of two equal mass black holes as the speed of light is approached. Nevertheless, the second order computation of D'Eath and Payne introduces an angular dependence [14], which seems to be crucial for a more accurate estimate.
One final remark concerning odd D. It is well known that the Green's functions of the D'Alembertian operator have support on the light cone in even D and both on and inside the light cone for odd D, cf. (A.12) and (A.13) (see, e.g. [35][36][37]). This leads to the unfamiliar property that odd D dimensional Minkowski space behaves as a dispersive medium. It may seem that the integral solution (3.14) is missing this property for odd D, since it has support solely at the intersection of the past light cone with the initial data surface (as illustrated in Fig. 5). But it is not so. The equivalence between the integral solution and the Green's function method is shown in appendix A.

Conclusions
In this paper we have provided an estimate for the energy radiated into gravitational waves in a collision of shock waves in a D dimensional spacetime. These shock waves provide an intuitive description, due to the infinite Lorentz boost, of the gravitational field of ultrarelativistic particles. Moreover, they do capture their gravitational interactions. Indeed, the scattering amplitude for two (scalar) particles in the eikonal regime of perturbative quantum gravity [38], can be equivalently computed by analysing the scattering of a (test) plane wave in a shock wave background [11] (see also [39]). Thus our computation applies to particle collisions in a regime wherein their interaction is dominated by gravity, and well described by classical general relativity. This is the transplanckian scattering regime [11].
If the fundamental Planck scale is of the order of the TeV scale, a key issue for phenomenology is how good the classical and semi-classical approximations are to describe the black holes that could be produced in particle collisions. In other words how much of the initial centre of mass energy of the process stays in the final black hole. The answer essentially amounts to understand how much energy is lost in gravitational radiation. The best estimates for this radiated energy, to date, come from apparent horizon computations [16]. In four dimensions the apparent horizon estimates are off by a factor of two as compared to improved values obtained either from perturbative computations in colliding shock waves backgrounds [12][13][14] or numerical relativity results [21]. The final black hole is therefore more massive (and hence more classical) than what has been anticipated by apparent horizon estimates. A similar conclusion in higher dimensions would have phenomenological relevance.
In this paper we have computed the metric in the future of the collision of two shock waves in first order perturbation theory and extracted the corresponding gravitational radiation. Although our result should only be faced as an estimate, what we know from D = 4 indicates that it is an improved estimate as compared to that obtained from the analysis of trapped surfaces [16].
The apparent horizon estimate gives the same trend for the D dependence as our result. This trend is also observed in numerical relativity results for a low energy collision of black holes in D = 4 and D = 5: these show that the energy radiated increases, from D = 4 to D = 5 [30]. Curiously, another bound that applies to these low energy collisions exhibits an opposite trend. From the area law, that bound may be obtained, as first shown by Hawking [40], for the energy emitted into gravitational radiation in a collision of nonspinning, equal mass black holes starting from rest at infinite distance. This bound starts at precisely the same value, 29.3%, as the one discussed herein for high energy collisions, but it decreases with dimension. Another context where a similar decrease with D occurs, though in a different regime -that of extreme mass ratio, was studied in [41][42][43], where the gravitational radiation emitted by point particles falling into a black hole was computed.
This paper also clarifies two aspects of the first order computation by D'Eath and Payne [12]. The first is the approximation procedure in determining the emitted power, namely neglecting ρ ′2 terms in the computation of (4.8). Since we have evaluated (4.8) numerically, no such approximation was used, and our result is consistent with that of D'Eath and Payne within less than one percent. This shows that the approximation used in [12] is robust in D = 4. The second is the time integration domain. It seems strange, at first sight, that there is a non-vanishing gravitational signal for all τ ∈ R in D = 4, in the limit r → ∞ for the observation point. This is a consequence of the nature of the coordinate transformation from Brinkmann to Rosen in D = 4 and the jump of the null generators of the weak shock to negative values of v as they cross the strong shock (cf. Fig. 2), for large ρ. In D > 4 the coordinate transformation is different and so is the time integration, which has support only for τ ∈ R + . Moreover, the observation we have made concerning the matching of the peaks in the obtained wave forms with a simple ray analysis, is useful in understanding the shape of the wave forms, which was not discussed in [12].
Finally, we hope that this paper will help setting the stage for further use of this technique, since it can be applied to various generalisations of the problem we are considering. In particular, one immediate goal is to consider the second order in perturbation theory. We observe that the initial conditions are exact to second order which might be at the basis of the good agreement between the perturbative computation and the numerical relativity result in D = 4.

A General integral solution
In this section we generalise the integral solution found in [12] to higher dimensions. The most standard method is to use Green's functions. In our case, however, there is a simple way of finding the solution by using an ansatz inspired by the four dimensional case. First we provide such a derivation and after we check it using the well know Green's functions method in D-dimensional Minkowski spacetime.

A.1 Integral operator ansatz
Before attempting to generalise the integral solution of the wave equation F = 0 in u ≥ 0 with initial data at u = 0 provided in [12], let us start by investigating how it works in four dimensions. In that case the solution is where, for each x ′ , v ′ defines points, at u = 0, on the past light cone of the event (u, v, x i ): We can compute partial derivatives of F by acting directly on the explicit u dependence, and the implicit dependence on u, v, x i through v ′ . Then using 3) we can compute F and check it is indeed zero. Furthermore, we need to check the initial condition. If we try to take the limit u → 0 directly in (A.1) we realise that there is no easy way to deal with the v ′ derivative. To remove it we can work in Fourier space by taking the Fourier transform with respect to v (we denote itF and replace the corresponding argument by q) where we have identified the parameter ǫ = −iu/q and a function δ ǫ which becomes the 2-dimensional Dirac delta distribution in the ǫ → 0 limit. So, in Fourier space, the initial condition is obeyed and the same holds after Fourier inversion. In higher dimensions we try a similar form where for the moment we leave the overall normalisation N arbitrary, v ′ is as before, and α is to be determined as well as the O v ′ operator acting on F . We can again compute F = 0 to obtain Thus this integral ansatz must have α = (D − 2)/2. Incidentally, this is exactly the correct number of u powers needed to obtain a D − 2 dimensional delta function by repeating the reasoning for the initial condition. Again, we Fourier transform with respect to v to obtaiñ So in the u → 0 limit we obtaiñ In higher dimensions, for the right hand side to reduce to the left hand side, we need a generalisation of the Fourier transform of the operator . This is defined through its Fourier transform when acting on a function f (v) (so it is a distributional operator) Using (A.9) and choosing N = (2π) (D−2)/2 , then (A.8) is solved. In even dimensions O v is just a partial derivative. However in odd dimensions we have a fractional partial derivative, thus formally, the general integral solution is

A.2 Green's function solution
A more standard derivation, consists of reducing our Cauchy problem with initial conditions on the null hypersurface u = 0, to a wave equation problem for a distribution with a source term. Then the general solution in D dimensions can be found, for example, from Theorem 6.3.1 of [44] F (u, v, The Green's function for even D is where the superscript denotes denotes the distributional derivative of order (D − 2)/2, of the Dirac delta distribution. For odd D where this distribution is to be understood as giving the finite part of the integral when acting on functions. In even dimensions, we can perform the u ′ , v ′ integration immediately and use the action of the distributional derivative of the delta function to get exactly (A.10).
For odd D, we first note that by working in Fourier space, we can define the negative order of the delta distribution, (p < 0) as This definition can be used to extend (recursively) to positive orders, by acting with derivatives on δ (p) (x). For example, the fractional derivative of the delta function of order 1/2, is simply and any higher order fractional derivative can be defined similarly. We note that up to singular delta function terms which we must discard (due to the finite part prescription for its integral with a function), we can write the D odd Green's function exactly as (A.12) with the fractional order delta distribution. Thus the same integral solution, with a fractional partial derivative, holds in odd dimensions in agreement with the first derivation.

B Solution for the transverse components
In Fourier space (with respect to v), the solution of (3.23) for the spatial components of the metric perturbations reads For convenience we write the initial perturbation as to separate the v dependent term: Then the Fourier transform of the initial perturbation is obtained as 3) The integration measure and ρ ′ are invariant under rotations. So we define a rotation, such that y i = R ij x j , is aligned with the y 1 direction and has magnitude ρ. Similarly where θ ′ projects the vector y ′ onto the y ′ 1 direction. In D = 4, θ ′ ∈ [0, 2π] whereas for D > 4 θ ′ ∈ [0, π]. The parity of the integrand with respect to y ′ 2 , . . . , y ′ D−2 is always well defined (there is only no definite parity with respect to y ′ 1 due to the exponential). Thus only the ℓ = m components are non-zero, either due to the δ ℓm or because if ℓ = m, then y ′ ℓ y ′ m is odd when integrated over y ′ 2 , . . . , y ′ D−2 . Changing to hyperspherical coordinates where we have defined Ω 0 = 2 to account for the fact that θ ′ is a polar angle in D = 4. Note that the first angular integral (corresponding to y ′ 1 y ′ 1 ) is just an area factor and the other components can be obtained using the vanishing of the trace. This last property is satisfied by the final result (B.7) which works as an independent check of the angular integrals. Fourier inversion then yields The exponent of the delta denotes the order of the derivative, which can be fractional or negative as in (A. 15). Note that due to the vanishing of the trace, actually only one component of the ij is independent. If we factor out the rotations, we can define the rotated 11 perturbation as (B.9) Now if we define x ≡ cos θ ′ , take a v derivative and use the scaling properties of the derivative of the delta distribution where U ≡ √ 2u and T ≡ √ 2v − ρ 2 /U . To perform the integral, we can always integrate by parts M = [(D − 4)/2] times (the brackets indicate the integer part) without obtaining any boundary terms 7 , so we end up with where q = 0 for D even and q = 1 for D odd. Then in even dimensions (note D is defined such that −1 ≤ x ⋆ ≤ 1) (B.13) and in odd dimensions, using (A.15), (B.14) The x integral can be performed by expanding the polynomials and integrating by parts. Then if we define the following polynomials we obtain (4.12) in odd dimensions and (4.13) in even dimensions.