Modeling the Beam of Gravitational Radiation from a Cosmic String Loop

The space of the tangent vectors $\mathbf{\dot{X}_{\pm}}$ is the surface of a unit sphere and $\mathbf{\dot{X}_{\pm}}$ trace out two separate curves on this unit sphere. If two tangent curves intersect they give rise to a cusp. When that happens a part of the string loop momentarily doubles back on itself in spacetime. Suppose the two tangent curves intersect at

$$\mathbf{\dot{X}_{-}}(\sigma^{c}_{-})=\mathbf{\dot{X}_{+}}(\sigma^{c}_{+}).$$

(8)

Without loss of generality, define $\tau^{c}=\left(\sigma_{+}^{c}+\sigma_{-}^{c}\right)/2$ and $\sigma^{c}=\left(\sigma^{c}_{+}-\sigma^{c}_{-}\right)/2$. It follows from Eq. (3) that

$$\mathbf{\dot{X}}^{2}(\tau^{c},\sigma^{c})=1,$$

(9)

where $\mathbf{\dot{X}}=d\mathbf{X}/d\tau$. The above result implies that the string moves at the speed of light at the cusps in the tangent vector direction. The high Lorentz boost in the regions near the cusp leads to beamed gravitational radiation and quite possibly emission of other particles for non-minimally coupled strings [26, 27, 28, 29].

A cosmic string loop of length $l$ has fundamental frequency $f_{1}=1/{T_{1}}=2/l$ where $T_{1}$ is the oscillation period of the loop. The undamped loop motion is periodic in the center of mass frame and all relevant functions may be expanded in terms of sinusoidal modes with frequencies $f_{m}=mf_{1}$ where $m\in\mathbb{Z}^{+}$. Since we work to lowest non-vanishing order in the string tension the waveform is a linear sum over the sinusoidal modes.

One way to characterize the gravitational wave emission is in terms of the power emitted. The time-averaged power emitted is quadratic in the wave contributions mode-by-mode since all cross terms between different modes average to zero. We write the time-averaged differential power emitted at mode $m$ per solid angle as $dP_{m}/d\Omega$ which is a function of $\hat{k}$, the direction of emission. Figure 1 illustrates the celestial sphere for emission direction $\hat{k}$ overlaying the tangent sphere for tangent curves $\mathbf{\dot{X}_{\pm}}$. Usually we will use spherical polar coordinates and write $\hat{k}=\left(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta\right)$. The power emitted in a single mode, $P_{m}$, is the integral of this differential power over the celestial sphere

$$P_{m}=\int d\Omega\frac{dP_{m}}{d\Omega}.$$

(10)

The details of the calculation of power emitted have been laid out in [30, 31, 32, 33]. For a cosmic string loop of length $l$, define the one-dimensional integrals,

$$I^{\mu}_{\pm}\equiv\frac{1}{l}\int_{-l/2}^{l/2}d\sigma_{\pm}\dot{X}^{\mu}_{\pm}e^{-\frac{i}{2}\omega_{m}k.X_{\pm}},$$

(11)

where $k^{\mu}=\left(1,\hat{k}\right)$ and $\omega_{m}=4\pi m/l$ for a given mode number $m$. The vectors $\hat{u}$ and $\hat{v}$

	$\displaystyle\hat{u}$	$\displaystyle=\left(-\sin\phi,\cos\phi,0\right)$		(12)
	$\displaystyle\hat{v}$	$\displaystyle=\left(-\cos\theta\cos\phi,-\cos\theta\sin\phi,\sin\theta\right)$		(13)

form a right-handed orthonormal basis with $\hat{k}=\hat{u}\times\hat{v}$. We will label the 3 elements by $\left(n_{1},n_{2},n_{3}\right)=\left(\hat{k},\hat{u},\hat{v}\right)$ below.

The Fourier Transform of the energy-momentum tensor ³³3For the energy-momentum tensor of the string $T_{\mu\nu}$, the Fourier Transform is obtained as $\tau_{\mu\nu}(k^{\lambda})=\frac{1}{T_{1}}\int_{0}^{T_{1}}dt\;\int d^{3}\mathbf{x}e^{-ik.x}T_{\mu\nu}(x^{\lambda})$ where $x^{\lambda}$ denotes the spacetime point. of the string loop is

$$\tau_{ij}=\frac{\mu l}{2}\left[I_{-}^{(n_{i})}I_{+}^{(n_{j})}+I_{-}^{(n_{j})}I_{+}^{(n_{i})}\right]$$

(14)

where

$$I_{\pm}^{(n_{i})}=\mathbf{I}_{\pm}.n_{i}.$$

(15)

The power emitted per unit solid angle in mode $m$ along direction $\hat{k}$ is [31]

$\displaystyle\frac{dP_{m}}{d\Omega}=\frac{G\omega_{m}^{2}}{\pi}\left[\tau^{*}_{pq}\tau_{pq}-\frac{1}{2}\tau^{*}_{qq}\tau_{pp}\right]$

(16)

where each of $p$ and $q$ take the values 2 and 3 (repeated indices are summed over). In terms of $I_{\pm}$, the differential power is

	$\displaystyle\frac{dP_{m}}{d\Omega}$	$\displaystyle=8\pi G\mu^{2}m^{2}\left[\left(|I_{-}^{(u)}|^{2}+|I_{-}^{(v)}|^{2}\right)\left(|I_{+}^{(u)}|^{2}+|I_{+}^{(v)}|^{2}\right)\right.$
		$\displaystyle\left.+4\operatorname{Im}(I_{-}^{(u)}I_{-}^{(v)*})\operatorname{Im}(I_{+}^{(u)}I_{+}^{(v)*})\right].$		(17)

Generic oscillating string loops radiate not only energy but also momentum and angular momentum. The radiation of net momentum results in a gravitational rocket effect on the string loops which may play an important role in their dynamical motions [35, 36]. The differential rate of radiation of momentum from a cosmic string loop is

$\displaystyle\frac{d\mathbf{\dot{p}}_{m}}{d\Omega}$

$\displaystyle=\frac{dP_{m}}{d\Omega}\hat{k}.$

(18)

Note that the momentum density follows directly from the flux of energy.

The calculation of the rate of radiation of angular momentum from a cosmic string loop is more complicated and has been discussed in detail in [31]. Introducing the integrals

$$M_{\pm}^{\mu\nu}\equiv\frac{1}{l}\int_{-l/2}^{l/2}d\sigma_{\pm}\dot{X}^{\mu}_{\pm}X^{\nu}_{\pm}e^{-\frac{i}{2}\omega_{m}k.X_{\pm}}$$

(19)

and defining

$$M_{\pm}^{(n_{i},n_{j})}\equiv\frac{1}{l}\int_{-l/2}^{l/2}d\sigma_{\pm}\left(\mathbf{\dot{X}}_{\pm}.n_{i}\right)\Bigl{(}\mathbf{X}_{\pm}.n_{j}\Bigr{)}e^{-\frac{i}{2}\omega_{m}k.X_{\pm}},$$

(20)

the Fourier Transform of the first moment of the string energy-momentum tensor is

	$\displaystyle\tau_{ijk}$	$\displaystyle=\frac{\mu l}{4}\left[I_{-}^{(n_{i})}M_{+}^{(n_{j},n_{k})}+I_{-}^{(n_{j})}M_{+}^{(n_{i},n_{k})}\right.$
		$\displaystyle\left.+I_{+}^{(n_{i})}M_{-}^{(n_{j},n_{k})}+I_{+}^{(n_{j})}M_{-}^{(n_{i},n_{k})}\right].$		(21)

The radiated angular momentum is orthogonal to $\hat{k}$. The differential rate of radiated angular momentum per unit solid angle of mode number $m$ is

$$\frac{d\mathbf{\dot{L}}_{m}}{d\Omega}=\frac{d\dot{L}_{m,u}}{d\Omega}\hat{u}+\frac{d\dot{L}_{m,v}}{d\Omega}\hat{v},$$

(22)

where

	$\displaystyle\frac{d\dot{L}_{m,u}}{d\Omega}$	$\displaystyle=\frac{G}{2\pi}\left[-i\omega_{m}\left(3\tau^{*}_{13}\tau_{pp}+6\tau^{*}_{3p}\tau_{p1}\right)\right.$
		$\displaystyle\left.-\omega_{m}^{2}\left(2\tau^{*}_{3pq}\tau_{pq}-2\tau^{*}_{3p}\tau_{pqq}-\tau^{*}_{pq3}\tau_{pq}+\frac{1}{2}\tau^{*}_{qq3}\tau_{pp}\right)\right.$
		$\displaystyle\left.\hskip 156.49014pt+c.c.\right],$		(23)
	$\displaystyle\frac{d\dot{L}_{m,v}}{d\Omega}$	$\displaystyle=\frac{G}{2\pi}\left[i\omega_{m}\left(3\tau^{*}_{12}\tau_{pp}+6\tau^{*}_{2p}\tau_{p1}\right)\right.$
		$\displaystyle\left.+\omega_{m}^{2}\left(2\tau^{*}_{2pq}\tau_{pq}-2\tau^{*}_{2p}\tau_{pqq}-\tau^{*}_{pq2}\tau_{pq}+\frac{1}{2}\tau^{*}_{qq2}\tau_{pp}\right)\right.$
		$\displaystyle\left.\hskip 156.49014pt+c.c.\right].$		(24)

Finally, the total radiated energy, momentum and angular momentum are

	$\displaystyle P$	$\displaystyle=$	$\displaystyle\sum_{m}P_{m}=\sum_{m}\int d\Omega\frac{dP_{m}}{d\Omega},$		(25)
	$\displaystyle\mathbf{\dot{p}}$	$\displaystyle=$	$\displaystyle\sum_{m}\mathbf{\dot{p}}_{m}=\sum_{m}\int d\Omega\frac{d{\mathbf{\dot{p}}_{m}}}{d\Omega},$		(26)
	$\displaystyle\mathbf{\dot{L}}$	$\displaystyle=$	$\displaystyle\sum_{m}\mathbf{\dot{L}}_{m}=\sum_{m}\int d\Omega\frac{d{\mathbf{\dot{L}}_{m}}}{d\Omega}.$		(27)

All three conserved quantities depends upon one dimensional integrals Eq. (11) and Eq. (19). Our focus is on the calculation of these basic quantities. It is important to recognize that the one-dimensional integrals are gauge-dependent but that the combinations giving the physical quantities above are gauge-independent.

To sketch the larger context, we examine the power emitted $P_{m}$, the rate of emission of momentum $\mathbf{\dot{p}}_{m}$ and the rate of emission of angular momentum $\mathbf{\dot{L}}_{m}$ as a function of the mode number $m$ for an example cosmic string loop.

For the calculations we select a particular family of loops (which we will refer to as “Vachaspati-Vilenkin loops” to distinguish from other, more symmetric classes of loops introduced in the later sections) with three frequencies and characterized by two parameters $\alpha$ and $\Phi$, introduced in [36],

	$\displaystyle X_{-}(\sigma_{-})$	$\displaystyle=\frac{l}{2\pi}\left\{\frac{2\pi}{l}\sigma_{-},\right.$
		$\displaystyle\left.-\left(1-\alpha\right)\sin\left(\frac{2\pi\sigma_{-}}{l}\right)+\frac{\alpha}{3}\sin\left(\frac{6\pi\sigma_{-}}{l}\right),\right.$
		$\displaystyle\left.-\left(1-\alpha\right)\cos\left(\frac{2\pi\sigma_{-}}{l}\right)-\frac{\alpha}{3}\cos\left(\frac{6\pi\sigma_{-}}{l}\right),\right.$
		$\displaystyle\left.-\sqrt{\alpha\left(1-\alpha\right)}\sin\left(\frac{4\pi\sigma_{-}}{l}\right)\right\},$		(28a)
	$\displaystyle X_{+}(\sigma_{+})$	$\displaystyle=\frac{l}{2\pi}\left\{\frac{2\pi}{l}\sigma_{+},\sin\left(\frac{2\pi\sigma_{+}}{l}\right),\right.$
		$\displaystyle\left.-\cos\Phi\cos\left(\frac{2\pi\sigma_{+}}{l}\right),-\sin\Phi\cos\left(\frac{2\pi\sigma_{+}}{l}\right)\right\}.$		(28b)

For the purposes of demonstration, let us select $\alpha=1/2$ and $\Phi=\pi/4$. The integrals Eq. (11) and Eq. (19) are performed numerically and combined to calculate $\frac{dP_{m}}{d\Omega},\frac{d|\mathbf{\dot{p}}_{m}|}{d\Omega}$ and $\frac{d|\mathbf{\dot{L}}_{m}|}{d\Omega}$. The differential quantities are integrated over the celestial sphere numerically to give $P_{m}$, $|\mathbf{\dot{p}}_{m}|$ and $|\mathbf{\dot{L}}_{m}|$ according to Eqs. (25), (26) and (27). The left hand panel of Fig. 2 shows the radiated power and the rate of radiated momentum and angular momentum per mode number for this loop up to $m=50$. The calculated power, rate of radiation of momentum and rate of radiation of angular momentum from the first $50$ modes are $53.8$ $G\mu^{2}$, $5.67G\mu^{2}$ and $4.18G\mu^{2}l$, respectively, in general agreement with [31]. The right hand panel shows the cumulative quantities up to and including mode $m$. Evidently, most of the total radiated power, momentum and angular momentum originates from the lower modes. These will be responsible for the dominant secular change in the loop’s invariant length, center of mass acceleration and spin changes. This particular loop has zero x-components of both $\mathbf{\dot{p}}_{m}$ and $\mathbf{\dot{L}}_{m}$.

We used numerical results to estimate the total power, momentum and angular momentum emitted by the loop from all the modes. We performed quadratures using a regular grid on the celestial sphere for all modes with $m\leq 200$ and using Monte Carlo sampling on the celestial sphere for modes $m=100-200$ and $207$, $248$, $298$, $358$, $429$, $515$, $619$, $743$ and $891$ (successive values increase by $\sim 1.2$). We used linear interpolation of the Monte Carlo results to infer the quadratures for $200<m\leq 891$. We extrapolated the results for $m>891$ based on the asymptotic scaling of the integrals for a dominant cusp (large $m$ quadratures $\propto m^{-4/3}$). We observed and verified the power law scaling for the numerically calculated quadratures with $300\lesssim m\leq 891$. Summing these three contributions (explicit, interpolated and extrapolated quadratures) allowed us to estimate the loop’s total power, total rate of radiation of momentum and total rate of radiation of angular momentum: $P=(74.0-74.2)G\mu^{2}$, $|\mathbf{\dot{p}}|=(9.6-9.7)G\mu^{2}$ and $|\mathbf{\dot{L}}|=(5.32-5.33)G\mu^{2}l$, respectively. The quoted range is derived by comparing regular versus Monte Carlo quadrature estimates in the overlap region $101\leq m\leq 200$. Approximately 80% of the total power, momentum magnitude and angular momentum magnitude radiated is sourced by modes with $m<121-122$, $527-549$ and $60$, respectively, for this particular loop (the range is as above). Although the upper limits for $m$ may appear variable they can be understood in terms of 3 qualitative factors: the intrinsic magnitude of radiated quantities, the possibility of sign cancellations in radiated quantities and the onset of asymptotic variation. Details are quantitatively discussed in Appendix A.

We find at least a factor of 2 variation among different Turok and Vachaspati-Vilenkin loops in terms of the total rates of emission but the qualitative conclusion that lower modes contribute the most remains true.

In general, the integrals $I^{\mu}_{\pm}$ cannot be reduced to an explicit closed form (except for a few special cases [36, 31, 30]). At low $m$ it is straightforward to use a direct numerical evaluation to derive an “exact” numerical result. The simple trapezoidal rule is exponentially convergent [37] for smooth periodic functions but there are two practical difficulties as $m$ grows large. First, the integrand oscillates more and more rapidly and the number of function evaluations needed scales $\propto m$. Second, the cancellation of signed quantities of order unity in the sum accumulates errors that may exceed the exact, small final result. Eventually, higher precision arithmetic must be used in the direct evaluation because the residual answer is so small. Numerical schemes exist that can integrate rapidly oscillating functions [38, 39] but in the next section we will take advantage of the special features of the integrand and develop a customized approach.

Direct numerical schemes related to the Fast Fourier Transform (FFT) were explored in [24, 22]. The methodology may be succinctly summarized as follows: Assume that the integrand, a periodic function in $\sigma_{\pm}$, is represented by $N$ discrete samples. A transformation of the original independent variable to a new form absorbs part of the periodic function leaving a pure harmonic (dependent on $m$) times a weight at the expense of introducing a new grid with uneven spacing. Adopt an interpolation scheme for the weight. Let the new uniform grid have dimension $M=cN$ with $c>>1$ (values $c=8-16$ were typical in [24]). Finally, use an FFT of length $M$ to perform the quadrature sum. FFT-related methods have the great advantage of computational speed for evaluating a range of mode numbers but suffer from aliasing effects because power at high frequencies is redistributed to low frequencies. As a practical matter, the contamination at low $m$ is small because the true power at high frequencies is small compared to the true power at low frequencies, i.e. the induced relative error is small at low frequencies; on the other hand, large $m$ results are significantly impacted in terms of relative error. We will see by direct calculation that the accuracy degrades as $m$ increases and, in general, the largest mode $m$ at which good accuracy can be achieved will satisfy $m<M/c$.

Let us focus on the complex analytic properties of the integration of $I^{\mu}_{-}$. We simplify the notation and consider one spacetime component $\mu$ and one harmonic $m\in\mathbb{Z}^{+}$ of the fundamental frequency. Let $h=I^{\mu}_{-}$, $z=\sigma_{-}/l$, $g=\dot{X}^{\mu}_{-}$ and $f=\omega_{1}k\cdot X_{-}/2$. The integral has the form

$$h=\int_{-1/2}^{1/2}dz\ g(z)e^{-imf(z)}$$

(29)

where $f$ and $g$ are periodic functions of $z$. For convenience, let $F=if$ so that the exponential above is $e^{-mF}$.

If we were interested in finding an asymptotic approximation to the integral at large $m$ we might follow the method of steepest descent [40]. The integrand is a smooth complex function without poles and with various critical points $z_{c}$ in the complex plane where $F^{\prime}(z_{c})=0$. The method of steepest descent proceeds by deforming the integration contour to pass through critical point(s) and also requiring that the imaginary part of $F(z)$ be constant over the deformed contour. The latter condition suppresses the oscillations of the integrand. Laplace’s method may be used to provide an approximation for the integral along a path for which the integrand is a maximum at the critical point. The technique is commonly used to generate asymptotic approximations (e.g. [40, 41]).

We apply the same ideas here to deform the integration path in the complex plane. A numerical calculation of $h$ along the new path will exactly match the direct real calculation along the old path. The new path also provides an asymptotic approximation to $h$ by way of the contributions at each of the critical points visited.

Figure 3 illustrates part of the complex $z$ plane for an example ⁴⁴4The particulars are not crucial at this point. This example is discussed in greater detail in a following section. It is Case 1 of the $X_{-}(\sigma_{-})$ mode of a Turok loop with $\alpha=1/10$ and with $\hat{k}$ direction implied by $\theta=7/10$ and $\phi=6/5$.. The original integration path over real $z$ is the light gray line that extends from one gray dot at $z=-1/2$ to the other gray dot at $z=1/2$. The dashed black, blue and green lines are contours for different fixed $\operatorname{Im}[F(z)]$. The red dots are critical points with $F^{\prime}(z_{c})=0$, a subset of all the critical points in the plane for this particular example. The original path is deformed to a connected set of 4 particular segments. The first segment starts at $z=-1/2$, proceeds to large imaginary $z$ (near vertical), then links to colored paths that begin and end at $\operatorname{Im}[z]=+\infty$ and the last segment returns to the real axis at $z=1/2$. The first segment (dashed) has $\operatorname{Im}[F(z)]=\operatorname{Im}[F(z=-1/2)]$, ascending in the direction $\operatorname{Im}[z]\to+\infty$ with large positive $\operatorname{Re}[F]$ (the integrand is zero). The last segment (also dashed) descends from infinite $\operatorname{Im}[z]$ and large positive $\operatorname{Re}[F]$ to $z=1/2$ along the contour with $\operatorname{Im}[F]=\operatorname{Im}[F(z=1/2)]$. These two segments are identical except for a shift in $z$ by $1$. Since $f$ and $g$ are periodic in $z$ the two integrations, one up and one down, exactly cancel each other. The remaining two segments form two looping jumps, from and to $\operatorname{Im}[z]=+\infty$ with large positive $\operatorname{Re}[F]$ at each end; each segment passes through a critical point at finite $z=z_{c}$ with $F^{\prime}(z_{c})=0$ and $F^{\prime\prime}(z_{c})>0$.

Carrying out an exact evaluation in the complex plane involves three tasks: locating the critical points, finding the complete closed path with piecewise-constant $\operatorname{Im}[F]$ and integrating along the path. The deformed contour need not be found exactly (in the sense that $\operatorname{Im}[F]$ is precisely constant) because any closed contour gives the same integrated result. The advantage of locating a contour with constant $\operatorname{Im}[F]$ is most significant when $m$ is large. We have compared numerical results for complex integration along the deformed path to the real integration along the original contour for $1\leq m\leq 10^{4}$ and achieved agreement at the level of machine precision ⁵⁵5We used the simple trapezoidal rule and also NIntegrate (the generic Mathematica integration recipes) for the real direct method. We used NIntegrate along a numerically determined path for the complex direct integration..

It is worth pointing out a few practical aspects of the numerical integration over the deformed path. The path is independent of $m$ so once the path is determined then large batches of $m$ can be done efficiently. The integration itself is relatively easy to carry out because, by design, the phase of $e^{-mF}$ doesn’t oscillate (the prefactor $g$ varies slowly). The maximum of the integrand for each segment occurs near the associated critical point $z_{c}$. The concentration of the integrand about the peak is related to $\operatorname{Re}[F^{\prime\prime}(z_{c})]$ – large, positive values imply sharper peaks near the critical point. The peak size $e^{-mF(z_{c})}$ sets the scale for additive arithmetic for any integration scheme. If the relative error is $\epsilon$ (stemming from the truncation error of the integration scheme) then the inferred absolute error $\sim\epsilon e^{-mF}$. Other than its dependence on $\epsilon$ the minimum absolute error that can be achieved numerically is constrained by the smallest number that can be represented ⁶⁶6For IEEE-754 double precision normalized floating point we estimate the absolute error $\max\{\epsilon e^{-mF},10^{-308}\}$.. This is quite different than the residual errors seen in the direct real calculations which involves cancellation of signed quantities of order unity.

The direct numerical schemes for real integration are suitable at small $m$. Here we describe a new method, which turns out to be appropriate for intermediate mode numbers $m$ where intermediate means larger than those needing direct methods and smaller than those most accurately evaluated by asymptotic means.

All previous approximations [36, 21, 22] to high mode number integrals exploit the existence of small neighborhoods where the integral builds up because the complex phase doesn’t oscillate too rapidly. These methods are ideal for describing the emission from cusps where the derivatives of both the left and right-moving oscillatory phases $k.X_{\pm}$ vanish. Write the mode conformal coordinates at the cusp as $\sigma_{\pm}=\sigma_{\pm}^{(c)}$. At the cusp the two tangent vectors $\mathbf{\dot{X}}_{\pm}(\sigma_{\pm}^{(c)})$ are identical. These pick out a special direction which becomes the center of expansion for the approximation. At large $m$ the emission is strongest along $\mathbf{\dot{X}}_{\pm}(\sigma_{\pm}^{(c)})$ and falls rapidly away from this special direction on the celestial sphere.

This approach can be modified to yield an improved description in two senses: it will treat emission even when the special direction associated with a cusp does not exist or is irrelevant and it will improve the fidelity of the beam shape. As we will see the improvements also will enable us to extend the description to somewhat lower mode frequencies than previous methods.

For each mode the integral Eq. (11) gets its maximum contribution if $\hat{k}$ aligns with some $\mathbf{\dot{X}}^{*}_{\pm}=\mathbf{\dot{X}}_{\pm}(\sigma_{\pm}^{*})$. But this situation does not hold for either mode for most emission directions. In general, $\hat{k}$ will point away from $\mathbf{\dot{X}}^{*}_{\pm}$,

$$k^{\mu}(\theta,\phi)=\dot{X}^{*\mu}_{\pm}+\delta_{\pm}^{\mu}$$

(30)

where $\delta_{\pm}^{\mu}$ corresponds to the deviation for each mode. With our definitions $k^{\mu}=(1,{\hat{k}})$ and $\dot{X}^{*\mu}=(1,{\hat{n}})$ for unit vectors ${\hat{k}}$ and ${\hat{n}}$ so that $\delta^{\mu}=(0,{\hat{k}}-{\hat{n}})$. For any $\sigma_{\pm}^{*}$, there is a corresponding $\delta_{\pm}^{\mu}$, which is fixed by the choice of $\theta$, $\phi$ and $\sigma_{\pm}^{*}$. Previous authors have generally worked in the limit of small $\delta_{\pm}^{\mu}$, but we will retain $\delta_{\pm}^{\mu}$ exactly for the moment. The four-vectors $X_{\pm}$ and $\dot{X}_{\pm}$ can be Taylor-expanded in a small region around $\sigma_{\pm}^{*}$. The largest contribution to the variation of the phase $k.X_{\pm}$ comes from the first three powers of $\sigma_{\pm}$. At the cusp, the leading contribution comes from the cubic term $\left(\mathcal{O}\left(\sigma_{\pm}^{3}\right)\right)$ while the linear $\left(\mathcal{O}\left(\sigma_{\pm}\right)\right)$ and quadratic terms $\left(\mathcal{O}\left(\sigma_{\pm}^{2}\right)\right)$ are zero. So, for the description of the phase, we shall retain terms up to cubic order in $\sigma_{\pm}$ whereas for the prefactor $\dot{X}_{\pm}$, we will retain the linear term only. With these assumptions, the Taylor expansion of $X_{\pm},k.X_{\pm}$ and $\dot{X}_{\pm}$ are

	$\displaystyle X_{\pm}(\sigma_{\pm})$	$\displaystyle=X_{\pm}^{*}+\dot{X}_{\pm}^{*}(\sigma_{\pm}-\sigma^{*}_{\pm})$
		$\displaystyle+\frac{1}{2}\ddot{X}_{\pm}^{*}(\sigma_{\pm}-\sigma^{*}_{\pm})^{2}$
		$\displaystyle+\frac{1}{6}X^{(3)*}_{\pm}(\sigma_{\pm}-\sigma^{*}_{\pm})^{3}$		(31)
	$\displaystyle k_{\mu}X^{\mu}_{\pm}(\sigma_{\pm})$	$\displaystyle=\left(\dot{X}^{*}_{\pm,\mu}+\delta_{\pm,\mu}\right)X^{\mu}_{\pm}(\sigma_{\pm})$
		$\displaystyle=\dot{X}_{\pm,\mu}^{*}X^{*\mu}_{\pm}-\frac{1}{6}|\ddot{X}_{\pm}^{*}|^{2}(\sigma_{\pm}-\sigma_{\pm}^{*})^{3}$
		$\displaystyle+\delta_{\pm,\mu}X^{*\mu}_{\pm}+\delta_{\pm,\mu}\dot{X}^{*\mu}_{\pm}(\sigma_{\pm}-\sigma_{\pm}^{*})$
		$\displaystyle+\frac{1}{2}\delta_{\pm,\mu}\ddot{X}^{*\mu}_{\pm}(\sigma_{\pm}-\sigma_{\pm}^{*})^{2}$
		$\displaystyle+\frac{1}{6}\delta_{\pm,\mu}X^{(3)*\mu}_{\pm}(\sigma_{\pm}-\sigma_{\pm}^{*})^{3},\quad\text{to cubic order}$		(32)
	$\displaystyle\dot{X}_{\pm}(\sigma_{\pm})$	$\displaystyle=\dot{X}_{\pm}^{*}+\ddot{X}_{\pm}^{*}(\sigma_{\pm}-\sigma^{*}_{\pm}),\quad\text{to linear order}$		(33)

where the terms in the expansion of the phase, Eq. (32), have been rewritten using Eq. (7). We define the coefficients,

	$\displaystyle A_{m}$	$\displaystyle\equiv-\frac{2\pi m}{l}\left(-\frac{1}{6}|\ddot{X}^{*}_{\pm}|^{2}+\frac{1}{6}\delta_{\pm,\mu}X^{(3)*\mu}_{\pm}\right),$		(34)
	$\displaystyle B_{m}$	$\displaystyle\equiv-\frac{2\pi m}{l}\frac{1}{2}\delta_{\pm,\mu}\ddot{X}^{*\mu}_{\pm},$		(35)
	$\displaystyle C_{m}$	$\displaystyle\equiv-\frac{2\pi m}{l}\delta_{\pm,\mu}\dot{X}^{*\mu}_{\pm},$		(36)
	$\displaystyle D_{m}$	$\displaystyle\equiv-\frac{2\pi m}{l}\left(\dot{X}^{*}_{\pm,\mu}X^{*\mu}_{\pm}+\delta_{\pm,\mu}X^{*\mu}_{\pm}\right),$		(37)

and make a change of variables,

	$\displaystyle\sigma_{\pm}$	$\displaystyle\rightarrow t-\frac{B_{m}}{3A_{m}},$		(38)
	$\displaystyle p$	$\displaystyle\rightarrow\frac{3A_{m}C_{m}-B_{m}^{2}}{3A_{m}^{2}},$		(39)
	$\displaystyle q$	$\displaystyle\rightarrow\frac{2B_{m}^{3}-9A_{m}B_{m}C_{m}+27A_{m}^{2}D_{m}}{27A_{m}^{3}}.$		(40)

Substituting for these quantities in the integral Eq. (11), shifting the origin $(\sigma_{\pm}-\sigma^{*}_{\pm})\rightarrow\sigma_{\pm}$ and extending the integration limits to $\pm\infty$ yields

	$\displaystyle I^{\mu}_{\pm}$	$\displaystyle=\left(\dot{X}^{*\mu}_{\pm}-\frac{B_{m}}{3A_{m}}\ddot{X}^{*\mu}_{\pm}\right)\int_{-\infty}^{\infty}dt\;e^{iA_{m}(t^{3}+pt+q)}$
		$\displaystyle+\ddot{X}^{*\mu}_{\pm}\int_{-\infty}^{\infty}dt\;t\;e^{iA_{m}(t^{3}+pt+q)}.$		(41)

The integrals

	$\displaystyle I_{1}$	$\displaystyle\equiv\int_{-\infty}^{\infty}dt\;e^{iA_{m}(t^{3}+pt+q)},$		(42a)
	$\displaystyle I_{2}$	$\displaystyle\equiv\int_{-\infty}^{\infty}dt\;t\;e^{iA_{m}(t^{3}+pt+q)}$		(42b)

can be evaluated analytically. See Appendix B for the explicit expressions.

If ${\hat{k}}$ aligns exactly with one of the tangent vectors $\mathbf{\dot{X}}_{\pm}(\sigma_{\pm})$ then $\delta_{\pm}^{\mu}=0$ and the results for that mode simplify as follows:

	$\displaystyle A_{m}$	$\displaystyle\rightarrow-\frac{2\pi m}{l}\left(-\frac{1}{6}|\ddot{X}^{*}_{\pm}|^{2}\right),$		(43)
	$\displaystyle B_{m}$	$\displaystyle\rightarrow 0$		(44)
	$\displaystyle C_{m}$	$\displaystyle\rightarrow 0$		(45)
	$\displaystyle D_{m}$	$\displaystyle\rightarrow-\frac{2\pi m}{l}\left(\dot{X}^{*}_{\pm,\mu}X^{*\mu}_{\pm}\right),$		(46)
	$\displaystyle\sigma_{\pm}$	$\displaystyle\rightarrow t$		(47)
	$\displaystyle p$	$\displaystyle\rightarrow 0$		(48)
	$\displaystyle q$	$\displaystyle\rightarrow\frac{D_{m}}{A_{m}}$		(49)
	$\displaystyle I^{\mu}_{\pm}$	$\displaystyle\rightarrow\dot{X}^{*\mu}_{\pm}I_{1}+\ddot{X}^{*\mu}_{\pm}I_{2}$		(50)
	$\displaystyle I_{1}$	$\displaystyle\rightarrow\frac{-2\pi e^{iD_{m}}}{A_{m}^{1/3}\Gamma(-1/3)}$		(51)
	$\displaystyle I_{2}$	$\displaystyle\rightarrow\frac{-i\pi e^{iD_{m}}}{A_{m}^{2/3}\Gamma(-2/3)}.$		(52)

In this case $I_{1}\propto m^{-1/3}$ and $I_{2}\propto m^{-2/3}$.

In the case that the two tangent vectors coincide forming a cusp write $l^{\mu}\equiv\dot{X}_{+}^{\mu}=\dot{X}_{-}^{\mu}$. Exact alignment of viewing direction is $k^{\mu}=l^{\mu}$. Now dropping constants for clarity $I^{\mu}_{\pm}\sim l^{\mu}I_{1}+\ddot{X}^{*\mu}_{\pm}I_{2}$. From above $I^{\mu}_{\pm}$ contain terms $\propto m^{-1/3}$ and $m^{-2/3}$. The stress energy tensor involves symmetrized products like $I_{+}^{(\mu}I_{-}^{\nu)}$ and such expressions involve powers $m^{-2/3}$, $m^{-1}$ and $m^{-4/3}$. Now [21] showed that all but the last were gauge dependent and could be removed from the stress tensor by a suitable coordinate transformation. In our treatment we do not make any explicit gauge transformations but rely on the fact that we are calculating physical observables even though they are written in terms of gauge dependent quantities. Such observables depend on the appropriate symmetrized combination of the one-dimensional integrals, e.g. the differential power in Eq. (17). Using the Virasoro conditions in Section II we find all the gauge-dependent terms identified by [21] explicitly vanish when we substitute $I^{\mu}_{\pm}\propto l^{\mu}I_{1}+\ddot{X}^{*\mu}_{\pm}I_{2}$. The final differential cusp power is $dP_{m}/d\Omega\propto m^{2}(m^{-4/3})^{2}\propto m^{-2/3}$.

For the cases where $\hat{k}$ aligns with only one of or neither of the tangent vectors the cross terms $\dot{X}_{-}.\dot{X}_{+}$ do not cancel out. Then $dP_{m}/d\Omega$ includes contributions from the lower powers of $m$. The gauge argument in [21] is valid only for the case of a cusp and in the direction of the cusp. It cannot be applied to the more general case. It is important to recognize that the asymptotic behavior as $m\to\infty$ will involve both exponentials and powers of $m$, i.e. one should not conclude that the drop off with $m$ is necessarily slower than $m^{-2/3}$. We will see, however, that at intermediate $m$ the emission in the exact cusp direction need not be dominant.

A similar procedure can be employed to estimate $M_{\pm}^{\mu\nu}$, with the prefactor $\dot{X}_{\pm}^{\mu}X_{\pm}^{\nu}$ expanded to linear order and the phase expanded to cubic order.

The cubic expansion Eq. (32) obviously ignores higher order terms. Consider the expansion of the phase up to the 5${}^{\text{th}}$ order, i.e.,

	$\displaystyle-\frac{2\pi m}{l}k_{\mu}X_{\pm}^{\mu}\left(\sigma_{\pm}\right)$	$\displaystyle=D_{m}+C_{m}(\sigma_{\pm}-\sigma^{*}_{\pm})+B_{m}(\sigma_{\pm}-\sigma^{*}_{\pm})^{2}$
		$\displaystyle+A_{m}(\sigma_{\pm}-\sigma^{*}_{\pm})^{3}+G_{m}(\sigma_{\pm}-\sigma^{*}_{\pm})^{4}$
		$\displaystyle+H_{m}(\sigma_{\pm}-\sigma^{*}_{\pm})^{5}$		(53)

where $A_{m},B_{m},C_{m},D_{m}$ are defined in Eq. (34) - Eq. (37) and

	$\displaystyle G_{m}$	$\displaystyle=-\frac{2\pi m}{l}\frac{1}{24}k_{\mu}X^{(4)*\mu}_{\pm},$		(54)
	$\displaystyle H_{m}$	$\displaystyle=-\frac{2\pi m}{l}\frac{1}{120}k_{\mu}X^{(5)*\mu}_{\pm}.$		(55)

Assume that the main contribution for the cubic phase occurs over width $|\sigma_{\pm}-\sigma^{*}_{\pm}|\sim 1/|A_{m}|^{1/3}$. Then the cubic term is dominant if both

$$\left(\frac{|A_{m}|}{|G_{m}|^{3/4}},\frac{|A_{m}|}{|H_{m}|^{3/5}}\right)\gtrsim\mathcal{O}(1).$$

(56)

The above condition sets an approximate lower bound on $m$ given by,

	$\displaystyle m$	$\displaystyle\gtrsim\frac{3l}{64\pi}\text{max}\left(\frac{|k.X^{(4)*}|^{3}}{|k.X^{(3)*}|^{4}},\frac{8}{5^{3/2}}\frac{|k.X^{(5)*}|^{3/2}}{|k.X^{(3)*}|^{5/2}}\right)$		(57)
		$\displaystyle\equiv m_{low,4+5}$		(58)

Whenever any of the aforementioned conditions is violated there is no good reason to assume the cubic expansion will provide reliable results. More generally any Taylor series expansion may have a limited radius of convergence so that even if one wanted to extend the method beyond cubic order it might not yield more accurate results.

There are other considerations as well. The Taylor expansion of the prefactor, $\mathbf{\dot{X}_{\pm}}$ at linear order must be a good approximation over the contributing width $|\sigma_{\pm}-\sigma^{*}_{\pm}|\sim|A_{m}|^{-1/3}$. A necessary condition is

$$\left|\frac{1}{2}\mathbf{\dddot{X}}^{*}_{\pm}\left(\sigma_{\pm}-\sigma^{*}_{\pm}\right)^{2}\right|<\left|\mathbf{\dot{X}}^{*}_{\pm}+\mathbf{\ddot{X}}^{*}_{\pm}\left(\sigma_{\pm}-\sigma^{*}_{\pm}\right)\right|.$$

(59)

In addition, the width of the peak should be small compared to the size of the fundamental domain ($|A_{m}|>>1$) and small compared to the distance to nearby peaks.

Collectively, these inequalities imply the existence of a lower mode $m$ for which Eq. (41) is a good approximation. We summarize this by $m>m_{low}\equiv\max\left(m_{low,4+5},...\right)$.

We will work exclusively at cubic order for the phase and at linear order for the prefactor and quantify the errors of that approach with respect to exactly known answers.