The Quadrupole Formula

1. Introduction — Why Quadrupole and Not Dipole?

In electromagnetism, the dominant radiation comes from the electric dipole moment. The dipole radiation power scales as $P_{\text{EM}} \propto |\ddot{\mathbf{d}}|^2$, where$\mathbf{d} = \sum q_i \mathbf{x}_i$ is the electric dipole moment. One might expect gravity to work similarly, with the "gravitational dipole" $\mathbf{d}_{\text{grav}} = \sum m_i \mathbf{x}_i$playing the analogous role. But there is a fundamental difference.

The gravitational "charge" is mass-energy, and there is only one sign of it (no negative masses). The mass dipole moment is related to the center of mass:

Its second time derivative is $\ddot{\mathbf{d}}_{\text{grav}} = M \ddot{\mathbf{x}}_{\text{cm}}$. For an isolated system, the total momentum $\mathbf{P} = M \dot{\mathbf{x}}_{\text{cm}}$ is conserved, so $\ddot{\mathbf{d}}_{\text{grav}} = \dot{\mathbf{P}} = 0$. The mass dipole radiation vanishes identically due to conservation of momentum.

Similarly, there is a "current dipole" moment associated with angular momentum$\mathbf{J} = \sum m_i \mathbf{x}_i \times \dot{\mathbf{x}}_i$. Its time derivative vanishes for an isolated system: $\dot{\mathbf{J}} = 0$. Therefore current dipole radiation also vanishes.

This has profound consequences: gravitational radiation is intrinsically weaker than electromagnetic radiation. While EM dipole radiation scales as $\omega^4$, GW quadrupole radiation scales as $\omega^6$, making it much harder to produce and detect. Only the most extreme astrophysical events — merging black holes, colliding neutron stars — produce detectable gravitational waves.

2. Multipole Expansion of the Retarded Solution

We start from the linearized Einstein equations in the Lorenz gauge. The trace-reversed perturbation $\bar{h}_{\mu\nu} = h_{\mu\nu} - \frac{1}{2}\eta_{\mu\nu}h$ satisfies the wave equation:

The retarded solution (outgoing waves only) is given by the Green's function integral:

where the retarded time is $t_{\text{ret}} = t - |\mathbf{x} - \mathbf{x}'|/c$. We now perform a multipole expansion for the far zone, where the observer distance $r = |\mathbf{x}|$ is much larger than the source size $r' = |\mathbf{x}'|$.

Step 1: Expand the Denominator

For $r \gg r'$, we expand:

where $\mathbf{n} = \mathbf{x}/r$ is the unit vector pointing from the source to the observer. Similarly, the retarded time becomes:

We can Taylor-expand the stress-energy tensor around the leading retarded time $t_r = t - r/c$:

Step 2: Identify the Multipole Terms

Substituting the expansion back into the integral and collecting terms by powers of $1/r$, the leading (monopole) contribution at order $1/r$ is:

For the $00$-component, $\int T^{00} d^3x = Mc^2$ is the total mass-energy, which is conserved. Thus the monopole term is static — it does not radiate. For the $0i$-component, $\int T^{0i} d^3x = c P^i$is the total momentum, also conserved in the center-of-mass frame.

Step 3: The Dipole Terms Vanish

The next-order (dipole) contribution involves integrals of the form$\int T^{\mu\nu} x'^k \, d^3x'$. Using the conservation law$\partial_\mu T^{\mu\nu} = 0$, one can show by integration by parts:

The second derivative gives $\ddot{D}^i = c^2 \dot{P}^i = 0$ by momentum conservation. The mass dipole moment $D^i = \int T^{00} x'^i d^3x' / c^2$ has a vanishing second derivative, producing no radiation. Similarly, the angular momentum (current dipole)$J^{ij} = \int (x'^i T^{0j} - x'^j T^{0i}) d^3x'$ is conserved, so $\dot{J}^{ij} = 0$.

This remarkable identity connects the spatial stress integral to the second time derivative of the second moment of the energy distribution. The proof follows from applying$\partial_\mu T^{\mu\nu} = 0$ twice. We define the mass quadrupole moment tensor and proceed to the radiation formula.

3. The Quadrupole Radiation Formula

Step 1: Define the Quadrupole Moments

The mass quadrupole moment tensor is defined as:

For a system of point masses: $I_{ij} = \sum_A m_A\, x_A^i\, x_A^j$. The reduced (trace-free) quadrupole moment is:

where $I_{kk} = \text{tr}(I) = \sum_k I_{kk}$. The trace-free part is what matters for radiation because the trace part, being proportional to $\delta_{ij}$, is removed by the transverse-traceless (TT) projection.

Step 2: The Far-Zone Field

Combining the retarded solution with the identity$\int T^{ij} d^3x' = \frac{1}{2}\ddot{I}_{ij}$, the spatial part of the trace-reversed perturbation in the far zone is:

where $t_r = t - r/c$ is the retarded time. All quantities on the right-hand side are evaluated at $t_r$.

Step 3: TT Projection

The physical gravitational wave is the transverse-traceless (TT) part of $\bar{h}_{ij}$. To extract it, we apply the TT projection operator$\Lambda_{ij,kl}$ (also written $\Pi_{ijkl}$). Given the propagation direction $\mathbf{n}$, define the transverse projector:

This projects any vector onto the plane perpendicular to $\mathbf{n}$. The full TT projector is:

This operator is symmetric in $(ij)$ and $(kl)$, transverse ($n^i \Lambda_{ij,kl} = 0$), traceless ($\delta^{ij}\Lambda_{ij,kl} = 0$), and idempotent ($\Lambda_{ij,mn}\Lambda_{mn,kl} = \Lambda_{ij,kl}$). The quadrupole radiation formula in TT gauge is therefore:

Step 4: Polarizations

For a wave propagating along $\mathbf{n} = \hat{z}$, only the $xy$-plane components survive the TT projection. The two polarizations are:

For a general direction $\mathbf{n} = (\sin\theta\cos\phi,\, \sin\theta\sin\phi,\, \cos\theta)$, one constructs an orthonormal basis $(\hat{e}_\theta, \hat{e}_\phi)$ in the transverse plane and projects accordingly.

4. Radiated Power — The GW Luminosity

The energy flux (power per unit area) carried by a gravitational wave is given by the Isaacson stress-energy tensor, averaged over several wavelengths:

The radiated power (luminosity) is obtained by integrating the energy flux over a sphere at distance $r$:

Substituting $h^{TT}_{ij} = \frac{2G}{c^4 r}\Lambda_{ij,kl}\ddot{Q}_{kl}$ and performing the angular integration over $d\Omega$ (which involves contracting the TT projectors over the sphere), the key integral identity is:

When contracted with the trace-free tensor $\dddot{Q}_{kl}\dddot{Q}_{mn}$(the third time derivative, since we took one more time derivative in going from$\ddot{Q}$ to $\dot{h}$), the trace terms drop out and we get:

Comparison with the Electromagnetic Larmor Formula

The electromagnetic analogue — the Larmor formula for dipole radiation — is:

Note the structural similarity: both are proportional to the square of time derivatives of multipole moments, divided by powers of $c$. The GW formula has $c^5$in the denominator (vs. $c^3$ for EM) because it involves one higher multipole order (quadrupole vs. dipole). The prefactor in the GW formula is:

Equivalently, in units where we factor out the mass and length scales, the characteristic GW power scale is:

This is the "Planck luminosity" — an astronomical number. The factor $G/c^5 \approx 2.76 \times 10^{-53}\;\text{W}^{-1}$in front of $\dddot{Q}^2$ explains why gravitational waves are so extraordinarily weak: you need enormous, rapidly accelerating masses to produce detectable radiation. Only compact binary systems with $v/c \sim \mathcal{O}(1)$ approach the Planck luminosity in their final merger.

5. Circular Binary System

Consider two point masses $m_1$ and $m_2$ in a circular orbit of separation $a$ and orbital angular frequency $\omega$. Define the total mass $M = m_1 + m_2$, the reduced mass $\mu = m_1 m_2/M$, and the symmetric mass ratio $\eta = \mu/M$.

Step 1: Compute the Quadrupole Moment

In the center-of-mass frame, the two bodies orbit at distances $r_1 = (m_2/M)a$and $r_2 = (m_1/M)a$ from the center. Their positions are:

The mass quadrupole moment is $I_{ij} = m_1 x_1^i x_1^j + m_2 x_2^i x_2^j$. Substituting and simplifying:

Using double-angle identities $\cos^2\omega t = \frac{1}{2}(1 + \cos 2\omega t)$and $\cos\omega t\,\sin\omega t = \frac{1}{2}\sin 2\omega t$, the trace-free quadrupole is:

Step 2: Time Derivatives

The second time derivative (needed for the strain):

The third time derivative (needed for the power):

Step 3: Strain Amplitudes

Note that the GW frequency is $f_{\text{GW}} = 2f_{\text{orb}} = \omega/\pi$(twice the orbital frequency, because the quadrupole has a $\cos 2\omega t$ dependence). For an observer along the $z$-axis (inclination $\iota = 0$, face-on):

For general inclination angle $\iota$ between the orbital angular momentum and the line of sight:

Using Kepler's third law $\omega^2 a^3 = GM$, we can re-express the amplitude as:

Step 4: The Chirp Mass

In the last expression, we introduced the chirp mass:

Step 5: Radiated Power and Orbital Decay

Computing $\dddot{Q}_{ij}\dddot{Q}^{ij}$ from the third derivative above:

(this is constant in time — no need for time-averaging). The radiated power is:

where in the last step we used Kepler's law to eliminate $\omega$. This energy loss comes at the expense of the orbital energy $E = -G m_1 m_2 / (2a)$. Setting $dE/dt = -L_{\text{GW}}$:

or equivalently $\dot{a}/a = -(64/5)(G^3 m_1 m_2 M)/(c^5 a^4)$. The orbit shrinks, the frequency increases, and the amplitude grows — the classic chirp signal.

Step 6: The Chirp Equation

Converting $\dot{a}$ to frequency evolution using Kepler's law$a = (GM/\omega^2)^{1/3}$ and $f_{\text{GW}} = \omega/\pi$:

Integrating the chirp equation gives the time to coalescence from frequency $f_0$:

6. Peters Formula for Eccentric Orbits

Real binary systems often have eccentric orbits, especially early in their evolution. Peters & Mathews (1963) generalized the quadrupole radiation results to Keplerian orbits with semi-major axis $a$ and eccentricity $e$.

Step 1: Orbit Parametrization

A Keplerian orbit has the separation $r(\phi) = a(1-e^2)/(1+e\cos\phi)$ where$\phi$ is the true anomaly. The orbital frequency is no longer constant; instead, Kepler's second law gives $\dot{\phi} = L/(\mu r^2)$ with angular momentum$L = \mu\sqrt{GMa(1-e^2)}$. To compute the time-averaged power, one must evaluate the quadrupole moment as a function of $\phi$ and average over one orbit.

Step 2: Orbit-Averaged Energy Loss

After extensive calculation (evaluating the quadrupole moment components in the orbital plane, taking three time derivatives, contracting, and averaging over one orbital period using the eccentric anomaly), Peters & Mathews obtained:

Step 3: Semi-Major Axis and Eccentricity Evolution

Using $E = -Gm_1 m_2/(2a)$ and $L^2 = G m_1 m_2 \mu a(1-e^2)$, Peters (1964) derived the orbit-averaged evolution equations:

Step 4: Circularization

A crucial feature of the Peters equations is that GW emission circularizes the orbit. To see this, note that $\langle de/dt \rangle < 0$for all $e > 0$ — eccentricity always decreases. More quantitatively, dividing the two Peters equations:

This can be integrated analytically (Peters 1964):

where $c_0$ is set by initial conditions. As $a \to 0$ (merger),$e \to 0$. Physical explanation: at pericenter (closest approach), velocities are highest and GW emission is strongest. Pericenter emission preferentially removes energy from the orbit while keeping the angular momentum relatively unchanged, making the orbit more circular.

7. Applications

The Hulse-Taylor Binary Pulsar PSR B1913+16

Discovered in 1974 by Russell Hulse and Joseph Taylor, PSR B1913+16 is a binary system consisting of two neutron stars with $m_1 \approx 1.4411\,M_\odot$ and$m_2 \approx 1.3874\,M_\odot$. The orbital parameters are:

The observed orbital decay matches the Peters formula prediction to better than 0.3%. This was the first indirect evidence for gravitational waves and earned Hulse & Taylor the 1993 Nobel Prize in Physics.

LIGO Waveforms

The quadrupole formula (and its post-Newtonian extensions) forms the basis of matched filtering — the primary detection technique used by LIGO/Virgo. Template waveforms are computed from the chirp equation and compared against detector output. The first direct detection, GW150914, was a binary black hole merger with $m_1 \approx 36\,M_\odot$, $m_2 \approx 29\,M_\odot$, and chirp mass $\mathcal{M}_c \approx 28.3\,M_\odot$. The signal swept from about 35 Hz to 250 Hz in the final ~0.2 seconds, exactly matching the chirp predicted by the quadrupole formula (with post-Newtonian corrections for the strong-field regime).

8. Historical Context

9. Python Simulation

The simulation below computes and plots three things: (1) the gravitational-wave chirp waveform $h_+(t)$ and $h_\times(t)$ from an inspiralling binary, showing the characteristic increase in both frequency and amplitude; (2) the time-frequency evolution of the chirp; and (3) the orbital decay of the Hulse-Taylor binary pulsar PSR B1913+16 compared to the GR prediction.