Gravitational Radiation Reaction

Our SPH code has been used extensively to study the properties of colliding neutron stars. In such systems, the general relativistic phenomenon of gravitational radiation losses can be very important. We include in our code a Newtonian-style treatment of radiation reaction, derived from the formalism of Blanchet et al. (1990), hereafter referred to as BDS, used by us (Faber & Rasio, 2000; Faber et al., 2001; Faber & Rasio, 2002) and other groups (Ruffert et al., 1997a; Ruffert et al., 1997b; Ruffert et al., 1996; Ayal et al., 2001). It is consistent to lowest order, and more accurate than many slow-motion approaches, which are completely incapable of modeling the spatial dependence of the radiation reaction force.

To include radiation reaction, we make several changes to the code. The BDS form for the radiation reaction form is given by

$$\vec{F}_{reac}=m_i\vec{a}_{reac}=m_i \frac{1}{c^5}\nabla U_5,$$ (27)

where $c$ is the speed of light in code units, and the relativistic potential $U_5$ is given by

$$U_5=\frac{2}{5}G\left(R-Q^{[3]}_{ab}x^a \frac{\partial\Phi}{\partial x^b}\right),$$ (28)

where $G$ is the gravitational constant, set equal to $1.0$ in our code by definition. Taking the derivative of the radiation reaction potential requires a second derivative of the standard gravitational potential $\Phi$, which we compute by finite differencing during a modified version of the gravity solver. Here, $a$ and $b$ are spatial indices to be summed over in a standard tensor product, i.e. they run from 1 to 3. The third time derivative of the traceless quadrupole moment, $Q^{[3]}_{ab}$, is given by taking the time derivative of the second derivative, which itself can be expressed, as shown in Rasio & Shapiro (1992), in SPH terms as

$$Q^{[3]}_{ab}=\frac{d}{dt}Q^{[2]}_{ab}=\frac{d}{dt} {\rm STF}\left [\sum_j 2m_j(v_j^a v_j^b +x_j^a a_{grav}^b)\right] .$$ (29)

Note that this quantity is a scalar. ``STF'' indicates to symmetrize and subtract away the trace of the tensor. The gravitational acceleration on a particle is denoted $a_{grav}$. The quantity $R$ which appears in the radiation reaction potential is a quantity whose value is given by solving a Poisson equation of the form

$$\nabla^2 R=-4\pi Q^{[3]}_{ab}x^a\frac{\partial \rho}{\partial x^b}.$$ (30)

To solve this equation numerically, we finite difference the density field on the grid directly to get the gradient of the density, rather than import directly to the grid the gradient which can be calculated using SPH techniques. We have found that this approach leads to a smoother and more realistic source term distribution

Besides the computation of the radiation reaction force, we also have to adjust very slightly our velocity equation. The specific momentum $\vec{v}$, representing the quantity whose derivative with respect to time is given by the hydrodynamic, AV, gravitational, and radiation reaction forces, is not the same quantity in the BDS formalism as the velocity $\vec{u}\equiv d\vec{x}/dt$, as it is in Newtonian physics. Instead, we find

$$u^a=v^a+\frac{4}{5c^5} Q^{[3]}_{ab}v^b.$$ (31)

Unfortunately, in order to keep our code simple, the notation we use does not match up with that used in the BDS paper. What we call ``$\vec{u}$'', they call ``$\vec{v}$'', and what we call ``$\vec{v}$'', they call ``$\vec{w}$''. Please be careful if you feel the need to edit these routines.

To include radiation reaction in a run, you must set the input parameter NGRAVRAD=1. Compute the speed of light in the code units, since it is often not unity, and use this to set the parameter ``SOL''.