text:
# Post-Newtonian Theory and its Application11footnote 1To appear in the Proc. of the 12th Workshop on General Relativity and Gravitation, M. Shibata (ed.)

## 1 Introduction

Recent years have shown a tremendous revival of interest in an Old Lady: the post-Newtonian approximation (or expansion when the speed of light $`c+\mathrm{}`$), which is surely the most important technique in the arsenal of General Relativity for drawing firm predictions for the outcome of experimental facts related to gravitation. The post-Newtonian approximation is ideally suited for describing the adiabatic phase of the very interesting astrophysical systems known as inspiralling compact binaries. These systems constitute our only known-to-exist (for sure) source to hunt for in the current network of laser-interferometric detectors of gravitational waves, today composed of the large-scale interferometers VIRGO and LIGO, and the medium-scale ones GEO and TAMA.

Two compact (i.e. gravitationally-condensed) objects — neutron stars or black holes — orbit an inward spiral, with decreasing orbital radius $`r`$, decreasing orbital period $`P`$, and increasing orbital frequency $`\omega =\frac{2\pi }{P}`$. The inspiral is driven by the loss of energy associated with the gravitational-wave emission, or, equivalently, by the action of radiation forces. There is a long phase of adiabatic inspiral, with associated dimensionless adiabatic parameter

$$\frac{\dot{\omega }}{\omega ^2}=𝒪\left(\frac{1}{c^5}\right).$$
(1)
The order of magnitude of the right-hand-side is that of the radiation-reaction force, namely $`1/c^5`$ or equivalently 2.5PN order beyond the standard Newtonian acceleration (following the usual post-Newtonian jargon). The binary’s dynamics being essentially “aspherical”, inspiralling compact binaries are strong emitters of gravitational radiation.

The main point about the theoretical description of inspiralling compact binary is that a model made of two structureless point-particles, characterized solely by two mass parameters $`m_1`$ and $`m_2`$ (and possibly two spins), is sufficient. Most of the non-gravitational effects usually plaguing the dynamics of binary star systems: a magnetic field, an interstellar medium, etc., are dominated by gravitational forces. However, the real justification for a model of point particles is that the effects due to the finite size of the compact bodies are small. In particular, the tidal interactions between the two compact objects are expected to play a little role during most of the inspiral phase; the mass transfer (in the case of neutron stars) does not occur until very late, near the final coalescence. Thus the inspiralling compact binaries are very clean systems, essentially dominated by gravitational forces during most of their life. This is the reason why these systems are so interesting for relativity theorists!

Our fascination for the old post-Newtonian lady is that it was recognized that improved waveform modelling is crucial for constructing efficient templates for searching and measuring the gravitational waves from inspiralling compact binaries in the LIGO/…/TAMA network. By improved we mean in fact much improved with respect to what is known from the “Newtonian” approximation, which corresponds in the case of the radiation field to the Einstein quadrupole formalism. And one needs much improved modelling because the orbital motion of the inspiralling compact binary — which is responsible for the gravitational-wave emission — becomes very relativistic when the binary reaches the so-called Innermost Circular Orbit or ICO (defined below in this article). At that point the orbital velocity is of the order of 50% of the speed of light. The price to be paid for applying the post-Newtonian approximation is that we must go to very high post-Newtonian order. The status of the field nowadays is the 3PN approximation (i.e. $`1/c^6`$) — or even better the 3.5PN one —, which is likely to be sufficient for practical purpose (at least in the case of neutron stars binaries).

After the two objects have passed the ICO they will plunge together and merge to form a single black hole, which will subsequently settle down into a stationary configuration, by emission of gravitational waves in the quasi-normal mode channel. Of course we should not expect the post-Newtonian approximation (and the description of the compact objects by point-particles) to be valid after the ICO, and we shall have to replace this model by a fully relativistic numerical computation of the phase of plunge and merger of two black holes . But, up to the point of the ICO, the post-Newtonian approximation, when carried out to 3PN order, is physically valid, and is probably very accurate. This last point may appear to be a little surprising (it contradicts some statements in the litterature about the “failure” of the post-Newtonian approximation in the regime of the ICO), but we shall present arguments supporting it. In particular the 3PN order is probably able to locate the ICO to within an accuracy of the order of 1% or better (for binary systems with comparable masses).

## 2 Post-Newtonian iteration for isolated systems

Before embarking into strong statements concerning the “accuracy” of the post-Newtonian description of binary systems, it is wise to examine in a somewhat general way the well-definiteness of the approximation itself. This is especially important in view of the fact that our old lady the post-Newtonian approximation has been unjustly accused of being plagued with some apparently inherent difficulties, and furthermore which crop up around the 3PN order we are interested in. The problems arise in the general case of regular (singularity-free) matter sources. Up to the 2.5PN order the approximation can be worked out without problems, and often at the 3PN order the problems can be solved specifically for each case at hands. However, it must be admitted that these difficulties, even if appearing at higher approximations, have cast doubt in the past on the actual soundness, on the theoretical point of view, of the post-Newtonian expansion. They pose the practical question of the reliability of the approximation when comparing the theory’s predictions with very precise experimental results. Here we discuss the nature of the problems – are they purely technical or linked with some fundamental drawback of the approximation? – and outline their resolution recently proposed in Ref. .

The first problem is that in higher approximations some divergent Poisson-type integrals appear. Recall that the post-Newtonian expansion replaces the resolution of an hyperbolic-like d’Alembertian equation by a perturbatively equivalent hierarchy of elliptic-like Poisson equations. Rapidly it is found during the post-Newtonian iteration that the right-hand-side of the Poisson equations acquires a non-compact support (it is distributed over all space), and that the standard Poisson integral diverges because of the bound of the integral at spatial infinity, i.e. $`r|𝐱|+\mathrm{}`$, with $`t=`$ const. For instance some of the potentials occuring at the 2PN order in Chandrasekhar’s work  are divergent, so the corresponding metric is formally infinite. In fact, Kerlick  showed that the post-Newtonian computation à la Chandrasekhar, following the iteration scheme of Anderson and DeCanio , can be made well-defined up to the 2.5PN order, but that this does not solve the problem at the next 3PN order, which has been found to involve some inexorably divergent Poisson integrals.

These divergencies come from the fact that the post-Newtonian expansion is actually a singular perturbation, in the sense that the coefficients of the successive powers of $`1/c`$ are not uniformly valid in space, since they typically blow up at spatial infinity like some positive powers of $`r`$. For instance, Rendall  has shown that the post-Newtonian expansion cannot be “asymptotically flat” starting at the 2PN or 3PN level, depending on the adopted coordinate system. The result is that the Poisson integrals are in general badly-behaving at infinity. Trying to solve the post-Newtonian equations by means of the standard Poisson integral does not a priori make sense. This does not mean that there are no solution to the problem, but simply that the Poisson integral does not constitute the correct solution of the Poisson equation in the context of post-Newtonian expansions. So the difficulty is purely of a technical nature, and will be solved once we succeed in finding the appropriate solution to the Poisson equation333The problem is somewhat similar to what happens in Newtonian cosmology. Here we have to solve the Poisson equation $`\mathrm{\Delta }U=4\pi G\rho `$, where the density $`\rho `$ of the cosmological fluid is constant all over space: $`\rho =\rho (t)`$. Clearly the Poisson integral of a constant density does not make sense, as it diverges at the bound at infinity like the integral $`r𝑑r`$. This nonsensical result has occasionally been referred to as the “paradox of Seeliger”. However the problem is solved once we realize that the Poisson integral does not constitute the appropriate solution of the Poisson equation in the context of Newtonian cosmology. A well-defined solution is simply given by $`U=\frac{2}{3}\pi G\rho r^2`$..

To cure the problem of divergencies we have introduced , at any post-Newtonian order, a generalized solution of the Poisson equation with non-compact support source, in the form of an appropriate finite part of the usual Poisson integral: namely we regularize the bound at infinity of the Poisson integral by means of a process of analytic continuation, analogous to the one already used to regularize the retarded integrals in Refs. . At any post-Newtonian order $`n`$ we have to solve a Poisson equation with non-compact-support “source-term” $`\overline{\tau }_n^{\mu \nu }`$, where the overbar indicates that we are considering a (formal) post-Newtonian expansion, and pick up the coefficient of $`1/c^n`$ in that expansion ($`\mu ,\nu `$ are space-time indices). We multiply the source-term by a regularization factor $`r^B`$, where $`r=|𝐱|`$ and $`B`$, and then apply the standard Poisson integral. The result,

$$\mathrm{\Delta }^1\left[r^B\overline{\tau }_n^{\mu \nu }\right](𝐱,t)=\frac{1}{4\pi }_^3\frac{d^3𝐲}{|𝐱𝐲|}|𝐲|^B\overline{\tau }_n^{\mu \nu }(𝐲,t),$$
(2)
is a well-defined integral (i.e. convergent at infinity, $`r+\mathrm{}`$) in some region of the $`B`$-complex plane, given by $`\mathrm{}(B)<a_{\mathrm{max}}2`$, where $`a_{\mathrm{max}}`$ denotes the maximal power of $`r`$ in the behaviour of the source-term $`\overline{\tau }_n^{\mu \nu }`$ when $`r+\mathrm{}`$, i.e. it corresponds to the maximal growth of the source-term at infinity (notice that $`a_{\mathrm{max}}`$ gets larger and larger when we increase the post-Newtonian order $`n`$). Next, we can prove that the latter function of $`B`$ generates a (unique) analytic continuation down to a neighbourhood of the origin $`B=0`$, except at $`B=0`$ itself, around which value it admits a Laurent expansion with multiple poles up to some finite order. Then, we consider the Laurent expansion of that function when $`B0`$ and compute the finite part ($`𝒫`$), or coefficient of the zero-th power of $`B`$, of that expansion. This defines our generalized Poisson integral:

$$\stackrel{~}{\mathrm{\Delta }^1}\left[\overline{\tau }_n^{\mu \nu }\right]=𝒫\mathrm{\Delta }^1\left[r^B\overline{\tau }_n^{\mu \nu }\right].$$
(3)
The resulting “generalized” Poisson integral constitutes an appropriate solution of the post-Newtonian equation. However this is only a particular solution of the Poisson equation, and the most general solution will be the sum of that particular solution and the most general solution of the corresponding homogeneous equation. At this stage, considering the post-Newtonian iteration scheme alone, we cannot do more and therefore we leave the homogeneous solution unspecified (we can see a posteriori that it is associated with radiation-reaction effects).

The second problem has to do with the near-zone limitation of the post-Newtonian approximation. Indeed the post-Newtonian expansion assumes that all retardations $`r/c`$ are small, so it can be viewed as a formal near-zone expansion when $`r0`$, which is valid only in the region surrounding the source that is of small extent with respect to the typical wavelength of the emitted radiation: $`r\lambda `$ (if we locate the origin of the coordinates $`r=0`$ inside the source). Therefore, the fact that the coefficients of the post-Newtonian expansion blow up at spatial infinity, when $`r+\mathrm{}`$, has nothing to do with the actual behaviour of the field at infinity. The serious consequence is that it is not possible, a priori, to implement within the post-Newtonian iteration the physical information that the matter system is isolated from the rest of the universe. Most importantly, the no-incoming radiation condition, imposed at past null infinity, cannot be taken into account, a priori, into the scheme. In a sense the post-Newtonian approximation is not “self-supporting”, because it necessitates some information taken from outside its own domain of validity.

The solution of the problem of the near-zone limitation of the post-Newtonian expansion resides in the matching of the near-zone field to the exterior field, a solution of the vacuum equations outside the source which has been developed in previous work  using some post-Minkowskian and multipolar expansions. In the case of post-Newtonian sources, the near zone, i.e. $`r\lambda `$, covers entirely the source, because the source’s radius itself is such that $`a\lambda `$. Thus the near zone overlaps with the exterior zone where the multipole expansion is valid. What we do is to impose a matching condition resulting from the numerical equality between the multipolar and post-Newtonian fields, verified in the external part of the near-zone (say $`a<r\lambda `$, where $`a`$ is the size of the source). The matching equation reads

$$\overline{(h^{\mu \nu })}(\overline{h}^{\mu \nu }),$$
(4)
where $``$ symbolizes the multipolar (actually multipolar-post-Minkowskian) series, and as before the overbar refers to the post-Newtonian or near-zone expansion. We emphasize that the matching equation is physically justified only for post-Newtonian sources, for which the exterior near-zone exists.

The requirement of matching to the post-Newtonian solution has been shown (in a previous work ) to entirely determine, up to any post-Newtonian order, the multipole moments parametrizing the exterior field. In our more recent paper  we have proved that the still undetermined homogeneous solutions in the post-Newtonian iteration alluded to before are also fully determined by the matching equation. These homogeneous solutions are associated with radiation-reaction effects — for instance they incorporate the dominant radiation-reaction force at the 2.5PN order, as well as the tail contribution in the radiation-reaction force which is known to arise at the 4PN order ). In conclusion, the post-Newtonian expansion of the field inside an isolated system is fully determined (it can be indefinitely reiterated, up to any post-Newtonian order) by the matching to the exterior solution satisfying the correct boundary condition at infinity — notably the absence of incoming radiation from past null infinity.

## 3 Dynamics and radiation of compact binaries

By equations of motion we mean the explicit expression of the accelerations of the compact bodies in terms of the positions and velocities. In Newtonian gravity, writing the equations of motion for a system of $`N`$ particles is trivial; in general relativity, even writing the equations in the case $`N=2`$ is difficult. The first relativistic term, at the 1PN order, was derived by Einstein, Infeld and Hoffmann  by means of their famous “surface-integral” method, by which the equations of motion are deduced from the vacuum field equations, and so which is applicable to any compact objects (they be neutron stars, black holes, or, maybe, naked singularities).

Concerning the 2PN and 2.5PN approximations the result for the equations of binary motion in harmonic coordinates was obtained by Damour and Deruelle and collaborators . The corresponding result for the ADM-Hamiltonian of two particles at the 2PN order was given in Ref. . By ADM-Hamiltonian we refer to the Fokker-type Hamiltonian, which is obtained from the matter-plus-field Arnowitt-Deser-Misner Hamiltonian by eliminating the field degrees of freedom. Kopeikin  derived the 2.5PN equations of motion for two extended compact objects. The 2.5PN-accurate harmonic-coordinate equations as well as the complete gravitational field (namely the metric $`g_{\mu \nu }`$) were computed by Blanchet, Faye and Ponsot , and by Pati and Will .

It is important to realize that the 2.5PN equations of motion are known to hold in the case of binary systems of strongly self-gravitating bodies. This is via an “effacing” principle for the internal structure of the bodies. As a result, the equations depend only on the “Schwarzschild” masses, $`m_1`$ and $`m_2`$, of the compact objects, as has been explicitly verified up to the 2.5PN order . The 2.5PN equations of motion have also been established by Itoh, Futamase and Asada , who use a variant of the surface-integral approach , valid for compact bodies.

The present state of the art is the 3PN order. The equations of motion at this order have been worked out independently by two groups, by means of different methods, and with equivalent results. On one hand, Jaranowski and Schäfer , and Damour, Jaranowski and Schäfer , employ the ADM-Hamiltonian formalism of general relativity; on the other hand, Blanchet and Faye , de Andrade, Blanchet and Faye , and Blanchet and Iyer , founding their approach on the post-Newtonian iteration initiated in Ref. , compute directly the equations of motion from which they infer the Lagrangian (instead of the Hamiltonian) in harmonic coordinates. The end results have been shown  to be physically equivalent in the sense that there exists a unique “contact” transformation of the dynamical variables, that changes the harmonic-coordinates Lagrangian obtained in Ref.  into a new Lagrangian, whose Legendre transform coincides exactly with the Hamiltonian given in Ref. . The 3PN equations of motion, however, depend on one unspecified numerical coefficient, $`\omega _{\mathrm{static}}`$ in the ADM-Hamiltonian formalism and $`\lambda `$ in the harmonic-coordinates approach, which is due to some incompleteness of the Hadamard self-field regularization method. This coefficient has been fixed by means of a dimensional regularization in Ref.  (see Section 5).

So far the status of the post-Newtonian equations of motion is quite satisfying. There is mutual agreement between all the results obtained by means of different approaches and techniques, whenever it is possible to compare them: point-particles described by Dirac delta-functions, extended post-Newtonian fluids, surface-integrals methods, mixed post-Minkowskian and post-Newtonian expansions, direct post-Newtonian iteration and matching, harmonic coordinates versus ADM-type coordinates, different processes or variants of the regularization of the self field of point-particles.

Let us remark that the 3PN equations of motion are merely “Newtonian” as regards the radiative aspects of the problem, because with that precision they contain the radiation reaction force at only the lowest 2.5PN order. From the conservative part of the 3PN dynamics (neglecting the 2.5PN radiation reaction) we shall obtain the binary’s center-of-mass energy $`E`$ at the 3PN order. We now want to compute the variation of $`E`$ because of the emission of gravitational radiation. For this purpose we replace the knowledge of the radiation reaction force in the local equations of motion by the computation of the total energy flux in gravitational waves, say $``$, and we apply the energy balance equation

$$\frac{dE}{dt}=.$$
(5)
Therefore in our approach the computation of the center-of-mass energy $`E`$ constitutes only one “half” of the solution of the problem. The second “half’, that of the computation of the energy flux $``$, is to be carried out by application of a wave-generation formalism.

Following earliest computations at the 1PN level  (at a time when the post-Newtonian corrections in $``$ had a purely academic interest), the energy flux of inspiralling compact binaries was completed to the 2PN order by Blanchet, Damour and Iyer , using the 2PN wave-generation formalism of Ref. , and, independently, by Will and Wiseman , using their own formalism (see Refs.  for joint reports of these calculations). The preceding approximation, 1.5PN, which represents in fact the dominant contribution of tails in the wave zone, had been obtained in Refs.  by application of a formula for tail integrals given in Ref. . Higher-order tail effects at the 2.5PN and 3.5PN orders, as well as a crucial contribution of tails generated by the tails themselves (the so-called “tails of tails”) at the 3PN order, were obtained by Blanchet . However, unlike the 1.5PN, 2.5PN and 3.5PN orders that are entirely composed of tail terms, the 3PN approximation involves also, besides the tails of tails, many non-tail contributions coming from the relativistic corrections in the (source) multipole moments of the binary. These have been completed by Blanchet, Iyer and Joguet , based on the expressions for the multipole moments given in Ref. , except for one single unknown numerical coefficient, due to the use of the Hadamard regularization, which is a combination of the parameter $`\lambda `$ in the equations of motion, and a new parameter $`\theta `$ coming from the computation of the 3PN quadrupole moment.

The post-Newtonian flux $``$ is in complete agreement, up to the 3.5PN order, with the result given by the very different technique of linear black-hole perturbations, valid in the “test-mass” limit where the mass of one of the bodies tends to zero (limit $`\nu 0`$, where $`\nu =\mu /m`$). Linear black-hole perturbations, triggered by the geodesic motion of a small mass around the black hole, have been applied to this problem . This technique has culminated with the beautiful analytical methods of Sasaki, Tagoshi and Tanaka , who solved the problem up to the extremely high 5.5PN order.

## 4 Post-Newtonian templates for compact binary inspiral

The orbital phase of the binary — i.e. the integral of the frequency $`\omega (t)`$:

$$\mathrm{\Phi }(t)=\omega (t)𝑑t,$$
(6)
constitutes the crucial quantity to be monitored (and therefore to be predicted) in the detectors. The templates built from the theoretical prediction for the orbital phase should be accurate enough over most of the inspiral phase, within the frequency bandwidth of the detectors, with reduced cumulative phase lags, so that the phasing errors are not significant when one attempts to extract the values of the binary’s parameters (essentially the masses and spins) from the data.

The relevant model for describing the inspiral phase consists of two point-masses moving under their mutual gravitational attraction. As a simplification for the theoretical analysis, the orbit of inspiralling binaries can be considered to be circular, apart from the gradual inspiral, with a good approximation. The templates are based on the energy-balance equation (5), from which one deduces the orbital phase as

$$\mathrm{\Phi }_c\mathrm{\Phi }=_{\omega _c}^\omega \frac{\omega dE}{},$$
(7)
where $`\mathrm{\Phi }_c`$ and $`\omega _c`$ denote the values at the instant of coalescence. The number of gravitational-wave cycles left from the current time till the coalescence instant (we consider only the dominant harmonics at twice the orbital frequency) is

$$𝒩=\frac{\mathrm{\Phi }_c\mathrm{\Phi }}{\pi }.$$
(8)
It is clear — because of Eq. (1) — that $`𝒩`$ is of the order of the inverse of radiation-reaction effects, hence the formal “post-Newtonian” order is $`𝒩=𝒪\left(c^{+5}\right)`$.

As a matter of fact, $`𝒩`$ will be a large number, approximately equal to $`1.6\times 10^4`$ in the case of two neutron stars between 10 and 1000 Hz (roughly the frequency bandwidth of the detector VIRGO). Data analysts have estimated that, in order not to suffer a too severe reduction of signal-to-noise ratio, one should monitor the phase evolution with an accuracy comparable to one gravitational-wave cycle (i.e. $`\delta 𝒩1`$), over the whole detector’s bandwidth. From a strict post-Newtonian point of view, we see that the consequence, since the “Newtonian” number of cycles is formally $`𝒪\left(c^{+5}\right)`$, is that any post-Newtonian correction therein that is larger than the order $`c^5`$ is expected to contribute to the phase evolution more than what is allowed by the previous estimate. Therefore, one expects that in order to construct accurate enough templates it will be necessary to compute the orbital phase up to at least the 2.5PN order at a minimum. This back-on-the-envelope estimate has been confirmed by measurement-accuracy analyses  which showed that in advanced generations of detectors the 3PN approximation (or, even better, the 3.5PN one) is required in the case of inspiralling neutron star binaries.

The first ingredient in the theoretical calculation is the binding energy $`E`$ (in the center-of-mass frame), defined as being the invariant energy associated with the conservative part of the binary’s 3PN dynamics (we ignore the radiation reaction effect at the 2.5PN order). The center-of-mass energy $`E`$ is deduced from the Hamiltonian in ADM-type coordinates , or equivalently from the Lagrangian in harmonic-coordinates . Restricting our consideration to circular orbits, the energy is a function or a single variable, the radial distance $`r`$ between the two particles in a given coordinate system. In fact it is better to express the energy in terms of the frequency $`\omega =\frac{2\pi }{P}`$ of the orbital motion, or, rather, in terms of the particular frequency-related parameter

$$x\left(\frac{GM\omega }{c^3}\right)^{2/3}.$$
(9)
The individual masses of the black holes are denoted by $`m_1`$ and $`m_2`$, and the total mass is $`M=m_1+m_2`$. The interest of using the parameter $`x`$, instead of some coordinate distance $`r`$, is that the energy function $`E(x)`$ then takes an invariant invariant (the same in different coordinate systems). The result we get consists of the Newtonian contribution, proportional to $`x`$, followed by post-Newtonian corrections up to the 3PN order:

$`E`$ $`=`$ $`{\displaystyle \frac{\mu c^2x}{2}}\{1+({\displaystyle \frac{3}{4}}{\displaystyle \frac{1}{12}}\nu )x+({\displaystyle \frac{27}{8}}+{\displaystyle \frac{19}{8}}\nu {\displaystyle \frac{1}{24}}\nu ^2)x^2`$ (10)
$`+`$ $`({\displaystyle \frac{675}{64}}+[{\displaystyle \frac{209323}{4032}}{\displaystyle \frac{205}{96}}\pi ^2{\displaystyle \frac{110}{9}}\lambda ]\nu {\displaystyle \frac{155}{96}}\nu ^2{\displaystyle \frac{35}{5184}}\nu ^3)x^3\}.`$
This expression involves the useful ratio between reduced and total masses:

$$\nu =\frac{\mu }{M}\text{where}\mu =\frac{m_1m_2}{M}.$$
(11)
This ratio is interesting because of its range of variation: $`0<\nu \frac{1}{4}`$, where $`\nu =\frac{1}{4}`$ in the equal-mass case and $`\nu 0`$ in the test-mass limit for one of the bodies. The parameter $`\lambda `$ denotes a point-mass regularization-constant and will be discussed in Section 5.

The second ingredient in this analysis concerns the gravitational-wave luminosity $``$, calculated, in the post-Newtonian approximation, from a wave-generation formalism valid for extended “fluid” systems , and then specialized to binary systems of point-particles . The calculation takes properly into account the relativistic corrections linked with the description of the source (multipole moments), as well as the non-linear effects in the propagation of the waves from the source to the far zone. The crucial input of any post-Newtonian computation of the flux is the mass quadrupole moment (because the post-Newtonian precision required for the higher multipole moments is smaller). The 3PN quadrupole moment for circular binary orbits in a harmonic coordinate system is of the form444We neglect a 2.5PN term and denote e.g. $`\widehat{x}_{ij}=x_ix_j\frac{1}{3}\delta _{ij}r^2`$, where $`x_i`$ is the harmonic-coordinate separation between the two bodies and $`r=|𝐱|`$.

$$I_{ij}=\mu \left(A\widehat{x}_{ij}+B\frac{r^2}{c^2}\widehat{v}_{ij}\right),$$
(12)
and has been obtained recently by Blanchet, Iyer and Joguet , with result

$`A`$ $`=`$ $`1+\left({\displaystyle \frac{1}{42}}{\displaystyle \frac{13}{14}}\nu \right)\gamma +\left({\displaystyle \frac{461}{1512}}{\displaystyle \frac{18395}{1512}}\nu {\displaystyle \frac{241}{1512}}\nu ^2\right)\gamma ^2`$
$`+`$ $`\left({\displaystyle \frac{395899}{13200}}{\displaystyle \frac{428}{105}}\mathrm{ln}\left({\displaystyle \frac{r}{r_0}}\right)+\left[{\displaystyle \frac{139675}{33264}}{\displaystyle \frac{44}{3}}\mathrm{ln}\left({\displaystyle \frac{r}{r_0^{}}}\right){\displaystyle \frac{44}{3}}\xi {\displaystyle \frac{88}{3}}\kappa \right]\nu +{\displaystyle \frac{162539}{16632}}\nu ^2+{\displaystyle \frac{2351}{33264}}\nu ^3\right)\gamma ^3,`$
$`B`$ $`=`$ $`{\displaystyle \frac{11}{21}}{\displaystyle \frac{11}{7}}\nu +\left({\displaystyle \frac{1607}{378}}{\displaystyle \frac{1681}{378}}\nu +{\displaystyle \frac{229}{378}}\nu ^2\right)\gamma `$ (14)
$`+`$ $`\left({\displaystyle \frac{357761}{19800}}+{\displaystyle \frac{428}{105}}\mathrm{ln}\left({\displaystyle \frac{r}{r_0}}\right)+\left[{\displaystyle \frac{75091}{5544}}+{\displaystyle \frac{44}{3}}\zeta \right]\nu +{\displaystyle \frac{35759}{924}}\nu ^2+{\displaystyle \frac{457}{5544}}\nu ^3\right)\gamma ^2,`$
where the post-Newtonian ordering parameter (defined in harmonic coordinates) reads

$$\gamma =\frac{GM}{rc^2}.$$
(15)
Notice the logarithms entering these formulas at the 3PN order. One type involves a constant length scale $`r_0`$ coming from the general wave-generation formalism of Refs. , and which corresponds to some “infrared” cut-off in the bound at infinity of the integrals defining the multipole moments. The constant $`r_0`$ should and will cancel out when considering the complete multipole expansion of the field exterior to the source. The other type of logarithm contains a different length scale $`r_0^{}`$, which should rather be viewed as an “ultra-violet” cut-off associated with the singular behaviour of the metric at the location of the point-particles; $`r_0^{}`$ is related to the two constants $`r_1^{}`$ and $`r_2^{}`$ — one for each particles — appearing in the 3PN equations of motion of Refs.  by

$$\mathrm{ln}r_0^{}=\frac{m_1}{M}\mathrm{ln}r_1^{}+\frac{m_2}{M}\mathrm{ln}r_2^{}.$$
(16)
As we know that $`r_1^{}`$ and $`r_2^{}`$ are gauge constants, i.e. they can be eliminated by a gauge transformation at the 3PN order, $`r_0^{}`$ will necessarily disappear from our physical, gauge-invariant, result at the end. Besides the harmless constants $`r_0`$ and $`r_0^{}`$, there are three point-mass regularization-constants in Eqs. (LABEL:13)-(14): $`\xi `$, $`\kappa `$ and $`\zeta `$ (see Section 5). However, we shall see that the energy flux for circular orbits depends on one combination only of these constants: $`\theta =\xi +2\kappa +\zeta `$.

Through 3.5PN order, the result is decomposed into an “instantaneous” part, i.e. generated solely by the multipole moments of the source, and a tail part. What we call here the tail part is in fact a complicated sum of “tails”, “tail squares”, and “tails of tails”, as computed in Ref. . We find, respectively,

$`_{\mathrm{inst}}`$ $`=`$ $`{\displaystyle \frac{32c^5}{5G}}\nu ^2\gamma ^5\{1+({\displaystyle \frac{2927}{336}}{\displaystyle \frac{5}{4}}\nu )\gamma +({\displaystyle \frac{293383}{9072}}+{\displaystyle \frac{380}{9}}\nu )\gamma ^2`$ (17)
$`+`$ $`({\displaystyle \frac{53712289}{1108800}}{\displaystyle \frac{1712}{105}}\mathrm{ln}\left({\displaystyle \frac{r}{r_0}}\right)`$
$`+`$ $`[{\displaystyle \frac{332051}{720}}+{\displaystyle \frac{123}{64}}\pi ^2+{\displaystyle \frac{110}{3}}\mathrm{ln}\left({\displaystyle \frac{r}{r_0^{}}}\right)+44\lambda {\displaystyle \frac{88}{3}}\theta ]\nu {\displaystyle \frac{383}{9}}\nu ^2)\gamma ^3\},`$
$`_{\mathrm{tail}}`$ $`=`$ $`{\displaystyle \frac{32c^5}{5G}}\nu ^2\gamma ^5\{4\pi \gamma ^{3/2}+({\displaystyle \frac{25663}{672}}{\displaystyle \frac{109}{8}}\nu )\pi \gamma ^{5/2}`$ (18)
$`+`$ $`\left({\displaystyle \frac{116761}{3675}}+{\displaystyle \frac{16}{3}}\pi ^2{\displaystyle \frac{1712}{105}}C{\displaystyle \frac{856}{105}}\mathrm{ln}(16\gamma )+{\displaystyle \frac{1712}{105}}\mathrm{ln}\left({\displaystyle \frac{r}{r_0}}\right)\right)\gamma ^3`$
$`+`$ $`({\displaystyle \frac{90205}{576}}+{\displaystyle \frac{3772673}{12096}}\nu +{\displaystyle \frac{32147}{3024}}\nu ^2)\pi \gamma ^{7/2}\}.`$
The Newtonian result has been factorized out in front. Here $`C=0.577\mathrm{}`$ denotes the Euler constant. As we can see, the constant $`r_0`$ drops out from the sum of Eqs. (17) and (18). However, the gauge constant $`r_0^{}`$ does not seem to disappear at this stage, but that is simply due to our use of the post-Newtonian parameter $`\gamma `$ defined by (15), and which depends via the equation of motion on the choice of harmonic coordinates. After substituting the frequency-related parameter $`x`$ in place of $`\gamma `$ (with the help of the 3PN equations of motion), we find that $`r_0^{}`$ does cancel as well — this nicely illustrates the consistency between our two computations, in harmonic-coordinates, of the equation of motion on one hand and the multipole moments on the other hand. Finally we obtain

$``$ $`=`$ $`{\displaystyle \frac{32c^5}{5G}}\nu ^2x^5\{1+({\displaystyle \frac{1247}{336}}{\displaystyle \frac{35}{12}}\nu )x+4\pi x^{3/2}`$ (19)
$`+`$ $`\left({\displaystyle \frac{44711}{9072}}+{\displaystyle \frac{9271}{504}}\nu +{\displaystyle \frac{65}{18}}\nu ^2\right)x^2+\left({\displaystyle \frac{8191}{672}}{\displaystyle \frac{535}{24}}\nu \right)\pi x^{5/2}`$
$`+`$ $`({\displaystyle \frac{6643739519}{69854400}}+{\displaystyle \frac{16}{3}}\pi ^2{\displaystyle \frac{1712}{105}}C{\displaystyle \frac{856}{105}}\mathrm{ln}(16x)`$
$`+`$ $`[{\displaystyle \frac{11497453}{272160}}+{\displaystyle \frac{41}{48}}\pi ^2+{\displaystyle \frac{176}{9}}\lambda {\displaystyle \frac{88}{3}}\theta ]\nu {\displaystyle \frac{94403}{3024}}\nu ^2{\displaystyle \frac{775}{324}}\nu ^3)x^3`$
$`+`$ $`({\displaystyle \frac{16285}{504}}+{\displaystyle \frac{176419}{1512}}\nu +{\displaystyle \frac{19897}{378}}\nu ^2)\pi x^{7/2}\}.`$
The last test (but not the least) is that Eq. (19) is in perfect agreement, in the test-mass limit $`\nu 0`$, with the result of linear black-hole perturbations obtained by Tagoshi and Sasaki .

## 5 Problem of regularization ambiguities

A model of structureless point masses is expected to be sufficient to describe the inspiral phase of compact binaries. Thus we want to compute the metric (and its gradient needed in the equations of motion) at the 3PN order for a system of two point-like particles. Applying general expressions for the metric valid in the case of continuous (smooth) matter sources, we find that most of the integrals become divergent at the location of the particles. Consequently we must supplement the calculation by a prescription for how to remove (i.e. to regularize) the infinite part of these integrals.

The “standard” Hadamard regularization yields some ambiguous results for the computation of certain integrals at the 3PN order, as Jaranowski and Schäfer  noticed in their computation of the equations of motion within the ADM-Hamiltonian formulation of general relativity. They showed that there are two and only two types of ambiguous terms in the 3PN Hamiltonian, which were then parametrized by two unknown numerical coefficients $`\omega _{\mathrm{static}}`$ and $`\omega _{\mathrm{kinetic}}`$.

Blanchet and Faye , motivated by the previous result, introduced their “improved” Hadamard regularization, based on a theory of pseudo-functions and generalized distributional derivatives. This new regularization is mathematically well-defined and free of ambiguities; in particular it yields unique results for the computation of any of the integrals occuring in the 3PN equations of motion. Unfortunately, this regularization turned out to be in a sense incomplete, because it was found  that the 3PN equations of motion involve one and only one unknown numerical constant, called $`\lambda `$, which cannot be determined within the method. The comparison of this result with the work of Jaranowski and Schäfer , on the basis of the computation of the invariant energy of binaries moving on circular orbits, showed  that

$`\omega _{\mathrm{kinetic}}`$ $`=`$ $`{\displaystyle \frac{41}{24}},`$ (20)
$`\omega _{\mathrm{static}}`$ $`=`$ $`{\displaystyle \frac{11}{3}}\lambda {\displaystyle \frac{1987}{840}}.`$ (21)
Therefore, the ambiguity $`\omega _{\mathrm{kinetic}}`$ is fixed, while $`\lambda `$ is equivalent to the other ambiguity $`\omega _{\mathrm{static}}`$. The value of $`\omega _{\mathrm{kinetic}}`$ given by Eq. (20) was recovered by Damour, Jaranowski and Schäfer , who proved that this value is the unique one for which the global Poincaré invariance of their formalism is verified. By contrast, the harmonic-coordinate conditions preserve the Poincaré invariance, and therefore the associated equations of motion should be Lorentz-invariant, as was indeed found to be the case by Blanchet and Faye , thanks in particular to their use of a Lorentz-invariant regularization  (hence their determination of $`\omega _{\mathrm{kinetic}}`$).

The appearance of one and only one physical undeterminacy $`\lambda \omega _{\mathrm{static}}`$ in the equations of motion constitutes a quite striking fact, specifically related to the use of an Hadamard-type regularization. Mathematically speaking, the presence of $`\lambda `$ is (probably) related to the fact that it is impossible to construct a distributional derivative operator satisfying the Leibniz rule for the derivation of the product. The Einstein field equations can be written into many different forms, by operating some terms by parts with the help of the Leibniz rule. All these forms are equivalent in the case of regular sources, but they become inequivalent for point particles if the derivative operator violates the Leibniz rule.

In Ref.  it has been argued that the numerical value of the parameter $`\omega _{\mathrm{static}}`$ could be $`9`$, because for such a value some different resummation techniques, viz Padé approximants  and effective-one-body (EOB) method , when they are implemented at the 3PN order, give approximately the same numerical result for the location of the last stable circular orbit. Even more, it was suggested  that $`\omega _{\mathrm{static}}`$ might be precisely equal to $`\omega _{\mathrm{static}}^{}`$, with

$$\omega _{\mathrm{static}}^{}=\frac{47}{3}+\frac{41}{64}\pi ^2.$$
(22)
(We have $`\omega _{\mathrm{static}}^{}=9.34\mathrm{}`$.) However, the value of $`\omega _{\mathrm{static}}`$ in general relativity has been computed by Damour, Jaranowski and Schäfer  by means of a dimensional regularization, instead of an Hadamard-type one, within the ADM-Hamiltonian formalism. Their result is

$$\omega _{\mathrm{static}}=0\lambda =\frac{1987}{3080}.$$
(23)
As Damour et al  argue, clearing up the ambiguity is made possible by the fact that the dimensional regularization, contrary to the Hadamard regularization, respects all the basic properties of the algebraic and differential calculus of ordinary functions555Note that the result (23) is very different from $`\omega _{\mathrm{static}}^{}`$ given by Eq. (22): this suggests that the resummation techniques (Padé approximants and EOB method), though they are designed to “accelerate” the convergence of the post-Newtonian series, do not in fact converge toward the same exact solution (or, at least, not as fast as expected). See Section 7 for discussion on this point..

Let us comment that the use of a self-field regularization in this problem, it be dimensional or based on the Hadamard partie finie, signals a somewhat unsatisfactory situation on the physical point of view, because we would like to perform, ideally, a complete calculation valid for extended bodies, taking into account the details of the internal structure of the bodies (energy density, pressure, internal velocities, etc.). By considering the limit where the radii of the objects tend to zero, one should recover the same result as obtained by means of the point-mass regularization, and determine the value of the regularization parameter $`\lambda `$. Because of considerable difficulties arising at the 3PN order this program has not yet been achieved.

Concerning the 3PN radiation field of two point masses — the second half of the problem, besides the 3PN equations of motion —, Blanchet, Iyer and Joguet  used the (standard) Hadamard regularization and found necessary to introduce three additional regularization constants $`\xi `$, $`\kappa `$ and $`\zeta `$. However the total gravitational-wave flux, in the case of circular orbits, depends only on the linear combination

$$\theta =\xi +2\kappa +\zeta .$$
(24)
Furthermore, this $`\theta `$ comes in at the same level as $`\lambda `$, so there is in fact only one unknown constant in the flux \[see Eqs. (17)-(19)\]. Notice that the improved version of the Hadamard regularization proposed in Refs.  should be able, in principle, to fix the value of the constant $`\zeta `$.

## 6 Accuracy of the post-Newtonian approximation

In this section we discuss (some aspects of) the accuracy of the post-Newtonian approximation as regards the determination of the binary’s innermost circular orbit (ICO). The ICO will be defined by the minimum, when it exists, of the energy function $`E(x)`$ given by Eq. (10). Notice that we do not define the ICO as an ISCO, i.e. as a point of dynamical general-relativistic unstability. See Ref.  for a discussion of the dynamical stability of circular binary orbits at the 3PN order.

Let us first confront the prediction of the standard (Taylor-based) post-Newtonian approximation at the 3PN order for the ICO — as given by the minimum of Eq. (10) — with a recent numerical calculation by Gourgoulhon, Grandclément and Bonazzola . These authors obtained numerically the energy $`E(\omega )`$ of binary black holes along evolutionary sequences of equilibrium configurations under the assumptions of conformal flatness for the spatial metric and of exactly circular orbits. The latter restriction is implemented by requiring the existence of an “helical” Killing vector, time-like inside the light cylinder associated with the circular motion and space-like outside. The numerical calculation  has been performed in the case of corotating black holes, which are spinning with the orbital angular velocity $`\omega `$. For the comparison we must therefore include within the post-Newtonian formalism yielding Eq. (10) the effects of spins appropriate to two Kerr black holes rotating at the orbital rate $`\omega `$.

The Figure 1 (issued from Ref. ) presents the post-Newtonian results for $`E_{\mathrm{ICO}}`$ in the case of irrotational and corotational binaries. The points indicated by 1PN, 2PN and 3PN are defined from the obvious truncation of Eq. (10). The points 1PNcorot, 2PNcorot and 3PNcorot take into account the spin effects of corotational binaries and are computed in Ref. . Notice that the irrotational and corotational configurations differ only from the 2PN order. As we can see the 3PN points, and even the 2PN ones, are rather close to the numerical value (marked by an asterisk in Figure 1). As expected, the best agreement is for the 3PN approximation and in the case of corotation: i.e. the point 3PNcorot. So our first conclusion is that the location of the ICO as computed by numerical relativity, under the helical-symmetry approximation , is in good agreement with post-Newtonian predictions. This constitutes an appreciable improvement of the previous situation, because we recall that the earlier estimates of the ICO in post-Newtonian theory  and numerical relativity  strongly disagree with each other, and do not match with the present 3PN results (see Ref.  for further discussion).

Our second conclusion comes from the fact that the 2PN and 3PN values in Figure 1 — either for irrotational or corotational binaries — are so close to each other. Indeed it seems that the 3PN points are useful only to confirm the result given by the 2PN ones. This is a quite satisfying state of affairs because it indicates that the post-Newtonian approximation converges very well (in the sense of Cauchy’s criterion). Figure 1 shows that the post-Newtonian approximation is well qualified to accurately locate the ICO (of course for this purpose one must go to the 3PN order — the 1PN approximation is clearly not accurate enough).

Let us elaborate more on this point. First we make a few order-of-magnitude estimates. At the location of the ICO we find (see Figure 1) that the frequency-related parameter $`x`$ defined by Eq. (9) is approximately of the order of 20%. Therefore, we might a priori expect that the contribution of the 1PN approximation to the energy at the point of the ICO should be of the order of 20%. For the present discussion we take the pessimistic view that the order of magnitude of an approximation represents also the order of magnitude of the higher-order terms that are neglected. We see that the 1PN approximation should yield a rather poor estimate of the “exact” result, but this is quite normal at this very relativistic point where the orbital velocity is $`\frac{v}{c}\sqrt{x}50\%`$. By the same argument we infer that the 2PN approximation should do much better, with fractional errors of the order of $`x^25\%`$, while the 3PN approximation will be even better, with the precision $`x^31\%`$.

The simple order-of-magnitude estimate suggests therefore that the 3PN order should be close to the “exact” solution for the ICO to within 1% of fractional accuracy. We think that this is very good, and we should even remember that this estimate is pessimistic, because we could reasonably expect that the neglected higher-order approximations, 4PN and so on, are in fact much smaller numerically (e.g. of the order of $`x^40.2\%`$). But let us keep for the present discussion the 1% guess for the accuracy of the 3PN approximation.

Now the previous estimates make sense only if the numerical values of the post-Newtonian coefficients in Eq. (10) are roughly of the order of one. If this is not the case, and if the coefficients increase dangerously with the post-Newtonian order $`n`$, one sees that the post-Newtonian approximation might in fact be very bad. So let us look at the values of the coefficients of the 1PN, 2PN and 3PN approximations in Eq. (10), say

$`a_1(\nu )`$ $`=`$ $`{\displaystyle \frac{3}{4}}{\displaystyle \frac{\nu }{12}},`$ (25)
$`a_2(\nu )`$ $`=`$ $`{\displaystyle \frac{27}{8}}+{\displaystyle \frac{19}{8}}\nu {\displaystyle \frac{\nu ^2}{24}},`$ (26)
$`a_3(\nu )`$ $`=`$ $`{\displaystyle \frac{675}{64}}+\left[{\displaystyle \frac{209323}{4032}}{\displaystyle \frac{205}{96}}\pi ^2{\displaystyle \frac{110}{9}}\lambda \right]\nu {\displaystyle \frac{155}{96}}\nu ^2{\displaystyle \frac{35}{5184}}\nu ^3.`$ (27)
We present in Table 1 the values of these coefficients in the test-mass limit $`\nu =0`$, and in the equal-mass case $`\nu =\frac{1}{4}`$ when the ambiguity parameter takes the “uncorrect” value $`\omega _{\mathrm{static}}^{}`$ \[see Eq. (22)\] and the correct one $`\omega _{\mathrm{static}}=0`$ predicted by general relativity .

Our first comment is that when $`\nu =0`$ there is an increase of the coefficients by roughly a factor 3 at each step. This behaviour is fairly easy to understand. It comes from the existence in the Schwarzshild metric of the famous light-ring orbit: a geodesics of photon which is a circular orbit located at $`R=3M`$ in Schwarzshild coordinates. As a result, the energy of the test particle in the Schwarzshild metric exhibits a singularity at $`x_{\mathrm{light}\mathrm{ring}}=\frac{1}{3}`$,

$$E^{\mathrm{test}}(x)=\mu c^2\left(\frac{12x}{\sqrt{13x}}1\right).$$
(28)
The consequence is that the radius of convergence of the post-Newtonian series is $`\frac{1}{3}`$, so the post-Newtonian coefficients increase by a factor $`3`$. So the post-Newtonian series is not very accurate in the case $`\nu =0`$. This fact has motivated several statements in the literature (see e.g. ), according to which the post-Newtonian approximation would be “poorly convergent”, or that there should be a “fundamental breakdown” of its validity in the regime of the ICO. This is indeed true in the test-mass limit ($`\nu 0`$), where the post-Newtonian series converges slowly666Let us remark that this negative conclusion does not matter: indeed we shall never use the post-Newtonian approximation in the case $`\nu 0`$ simply because we know the exact result which is given by Eq. (10). The exact result for the radiation field is known as well, albeit numerically only . Therefore we should not worry too much about the poor convergence of the post-Newtonian series in the test-mass limit. The post-Newtonian method is useless and even one might say irrelevant when considering the motion of a test particle around a Schwarzschild black hole..

On the other hand, what happens in the equal-mass case ($`\nu =\frac{1}{4}`$)? When $`\nu =\frac{1}{4}`$ and we have the value $`\omega _{\mathrm{static}}^{}`$ we notice that the coefficients increase approximately like in the test-mass case $`\nu =0`$. This indicates that the gravitational interaction in the case of $`\omega _{\mathrm{static}}^{}`$ looks like that in a one-body problem. We shall say that when the post-Newtonian coefficients rapidly increase with the order of approximation, the interaction is “Schwarzschild-like”.

Now when $`\nu =\frac{1}{4}`$ and the ambiguity parameter takes the correct value $`\omega _{\mathrm{static}}=0`$, we see that the 3PN coefficient $`a_3(\frac{1}{4})`$ is of the order of minus one instead of $`10`$. We think that this strongly suggests, unless 3PN happens to be quite accidental, that the post-Newtonian coefficients in general relativity do not increase very much with $`n`$, and stay rather of the order of one. This is interesting as it indicates that the actual general-relativistic two-body interaction is not Schwarzschild-like.

It is impossible of course to be very confident about the validity of the previous statement because we know only the coefficients up to the 3PN order. Any tentative conclusion based on the 3PN order can be “falsified” when we obtain the next 4PN order. Nevertheless, we feel that the mere fact that $`a_3(\frac{1}{4})=0.97`$ in Table 1 is sufficient to motivate our (tentative) conclusion that the field of two bodies is more complicated than the Schwarzschild space-time. This conclusion is in accordance with the present author’s respectfulness of the complexity of the Einstein field equations.

The nice consequence is that because the post-Newtonian coefficients when $`\nu =\frac{1}{4}`$ stay of the order of one, the standard post-Newtonian approach, based on the standard Taylor approximants, is probably very accurate. The post-Newtonian series seems to “converge well”, with a “convergence radius” of the order of one777Actually, the post-Newtonian series could be only asymptotic (hence divergent), but nevertheless it should give good results provided that the series is truncated near some optimal order of approximation. In this discussion we assume that the 3PN order is not too far from that optimum.. A convincing support of this view is provided by Figure 1 itself. The order-of-magnitude estimates we did at the beginning of this Section are probably correct. In particular the 3PN order should be close to the “exact” solution even in the regime of the ICO.

It is also interesting to look at the numerical values of the post-Newtonian coefficients in the total flux $``$ obtained in Eq. (19). For this case we do not have (like for the ICO) a clear point at which we can compare with a result of numerical relativity. Furthermore things concerning the convergence of the post-Newtonian series are less clear than with the energy function. A possible reason for that could be that the flux, in contrast to the energy, includes the contribution of tails at the 1.5PN, 2.5PN and 3.5PN orders, and even of tails-of-tails at the 3PN order.

## 7 On Schwarzschild-like templates for binary inspiral

Let us finally comment about a possible implication of our conclusion as regards the validity of the so-called post-Newtonian resummation techniques, i.e. Padé approximants , which aim at “accelerating” the convergence of the post-Newtonian series in the pre-coalescence stage, and effective-one-body (EOB) methods , which attempt at describing the late stage of the coalescence of two black holes. These techniques are based on the idea that the gravitational two-body interaction is a “deformation” — with $`\nu \frac{1}{4}`$ being the deformation parameter — of the Schwarzschild space-time.

The Padé approximants are valuable tools for giving accurate representations of functions having some singularities. In the problem at hands they would be justified if the “exact” expression of the energy \[whose 3PN expansion is given by Eqs. (25)-(27)\] would admit some Schwarzschild-like features, and in particular a light-ring singularity. In the Schwarzschild case the Padé series converges rapidly toward the solution : the Padé constructed only from the 2PN approximation of the energy — keeping only $`a_1(0)`$ and $`a_2(0)`$ in Eqs. (2) — already coincide with the exact result given by Eq. (8). On the other hand, the EOB method maps the post-Newtonian two-body dynamics (at the 2PN or 3PN order) on the geodesic motion on some effective metric which happens to be a $`\nu `$-deformation of the Schwarzschild space-time. In the EOB method the effective metric looks like Schwarzschild by definition, and we might expect the two-body interaction to own some Schwarzschild-like features.

Our comment is that the validity of these resummation techniques (and the templates for binary inspiral built on them) is questionable, because as we have seen in Section 6 the value $`\omega _{\mathrm{static}}=0`$ suggests that most probably the two-body interaction is not Schwarzschild-like. In particular, there does not seem to exist something like a light-ring orbit which would be a deformation of the Schwarzschild one.

To test the previous idea, let us come back to the Schwarzschild limit, where we observed that the radius of convergence of the post-Newtonian series is given by the light-ring singularity at the value $`\frac{1}{3}`$. Now the radius of convergence is given by d’Alembert’s criterion as the limit when $`n+\mathrm{}`$ of the ratio of successive coefficients, i.e. $`a_{n1}/a_n`$. We might therefore compute the light-ring orbit by investigating the limit

$$x_{\mathrm{light}\mathrm{ring}}^{\mathrm{test}}=\underset{n+\mathrm{}}{lim}\frac{a_{n1}^{\mathrm{test}}}{a_n^{\mathrm{test}}}=\frac{1}{3}.$$
(29)
To test for the possible existence of a light-ring singularity in the case of comparable masses, we consider the ratio between the two highest known post-Newtonian coefficients, that are $`a_2(\nu )`$ and $`a_3(\nu )`$. This ratio will give an estimate of the “light-ring” singularity for non-infinitesimal mass ratios,

$$x_{\mathrm{light}\mathrm{ring}}(\nu )\frac{a_2(\nu )}{a_3(\nu )}.$$
(30)
Using the values given in Table 1 we obtain for equal masses in the case of the “wrong” ambiguity parameter $`\omega _{\mathrm{static}}^{}9.34`$,

$$x_{\mathrm{light}\mathrm{ring}}(\frac{1}{4},\omega _{\mathrm{static}}^{})0.32.$$
(31)
As we see there seems to be in this case a (pseudo-)light-ring orbit which is a small deformation of the Schwarzschild light-ring orbit given by Eq. (29). In this Schwarzschild-like situation, we expect that Padé approximants and EOB-type methods to be appropriate.

But, in the case of general relativity ($`\omega _{\mathrm{static}}=0`$), we obtain a drastically different result:

$$x_{\mathrm{light}\mathrm{ring}}(\frac{1}{4},G.R.)2.86.$$
(32)
If we believe the correctness of this estimate we must conclude that there is in fact no notion of a light-ring orbit in the real two-body case. Or, one might say (pictorially speaking) that the light-ring orbit gets hidden inside the horizon of the final black-hole formed by coalescence. If we consider the ratio between the 1PN and 2PN coefficients (instead of between 2PN and 3PN), we get the value $`0.28`$ instead of Eq. (32). So at the 2PN order the field seems to admit a light ring, while at the 3PN order it does not. This reinforces our idea that it is meaningless (with our present 3PN-based knowledge, and untill fuller information is available), to assume the existence of a light-ring singularity in the equal-mass case. Our expectation, therefore, is that the conditions under which the Padé and EOB methods should be legitimate might in fact not be fulfilled.

This doubt is confirmed by the finding of Ref.  (already alluded to above) that in the case of the “wrong” ambiguity parameter $`\omega _{\mathrm{static}}^{}9.34`$ the Padé approximants and the EOB method at the 3PN order give the same result for the ICO. From the previous discussion we see that this agreement is to be expected because a deformed (pseudo-)light-ring singularity seems to exist with $`\omega _{\mathrm{static}}^{}`$. By contrast, in the case of general relativity, where $`\omega _{\mathrm{static}}=0`$, the Padé and EOB methods give quite different results (cf the figure 2 in Ref. ). Another confirmation comes from the light-ring singularity which is determined from the Padé approximants at the 2PN order (see Eq. (3.22) in Ref. ) as $`x_{\mathrm{light}\mathrm{ring}}^{\mathrm{P}2\mathrm{P}\mathrm{N}}(\frac{1}{4})0.44`$. This value is rather close to Eq. (31) but strongly disagrees with Eq. (32). Our explanation is that the Padé series might converge toward a theory having $`\omega _{\mathrm{static}}=\omega _{\mathrm{static}}^{}`$ and therefore which is different from general relativity.


format: text

text:
# 1 Introduction

## 1 Introduction

In the last five years there have been a number of neutron stars discovered to have magnetic fields of magnitude only imagined by a few theorists. Previous to these recent discoveries, the highest measured neutron star fields were around $`10^{13}`$ G. The Parkes Multibeam Survey (Manchester et al. 2001, Lorimer 2003) has now detected a number of new radio pulsars having surface magnetic fields, deduced from their dipole spin-down rates, of near $`10^{14}`$ G. Even more remarkable are the Soft Gamma-Ray Repeaters (SGRs) (Kouveliotou 2003) and Anomalous X-Ray Pulsars (AXPs) (Gavril 2003), both now believed to be the predicted magnetars (Thompson & Duncan 1992) with surface fields possibly as high as $`10^{15}`$ G. The fields of these new sources are well above the quantum critical magnetic field of $`B_{\mathrm{cr}}=m^2c^3/e\mathrm{}`$($`4.413\times 10^{13}`$ G for electrons), where physical processes are profoundly changed and new processes unique to strong fields dominate. Observations of these sources can not only provide direct evidence for the high fields but will begin to study in detail the exotic processes which can take place only at these field strengths.

This review will discuss the physics of some of the radiative processes taking place in ultrastrong magnetic fields and their possible observational signatures in high-energy wavebands. Other recent and useful reviews of physics in strong fields have been given by Duncan (1999) and of atomic processes and matter in strong magnetic fields by Lai (2001).

## 2 Basics of Physics in Strong Fields

Magnetic fields can affect radiative processes in a number of ways. The critical magnetic field strength, $`B_{\mathrm{cr}}=m^2c^3/e\mathrm{}`$ where particle’s cyclotron energy equals its rest mass, is a useful measure of when relativistic effects become important. In neutron star fields, relativistic effects are always important for electrons, whose critical field is $`4.413\times 10^{13}`$ G, but never for protons, whose critical field is $`1.43\times 10^{20}`$ G. Charged particles occupy Landau states, which arise from the solutions to the Dirac Equation in a homogeneous magnetic field, with energy (in rest mass units)

$$E_n=(1+p^2+2nB^{})^{1/2},n=l+\frac{1}{2}(s+1)=0,1,2,\mathrm{}$$
(1)
where $`n`$, $`l`$, and $`s`$ are the principal, orbital and spin quantum numbers, respectively and $`B^{}=B/B_{\mathrm{cr}}`$. Particles may have either spin up ($`s=1`$) or down ($`s=1`$), except in the ground state $`n=0`$, where only the spin-down state is allowed. The two spin states in levels $`n>0`$ are thus degenerate. This degeneracy is however broken by higher order interaction of the particles with the radiation field (Herold et al. 1982, Pavlov et al. 1991) which causes a fine-structure splitting of all the $`n>0`$ levels. The momentum component parallel to the field, $`p`$, is continuous but momentum perpendicular to the field is quantized. In magnetic fields $`B^{}>1`$, the transverse motion becomes relativistic. In the non-relativistic limit ($`p1,B^{}1`$), the transverse energy becomes $`nmcB^{}=neB/mc=nϵ_{\mathrm{cr}}`$, multiples of the cyclotron energy.

Since the magnetic field can absorb or supply momentum, momentum perpendicular to B is not conserved, although parallel momentum and total energy are strictly conserved. This allows a number of processes to take place that would be forbidden in free-space.

A magnetic field also strongly affects the photon polarization states and propagation modes. In a magnetized plasma the two normal modes of propagation are termed the extraordinary (X) mode, which is mostly polarized perpendicular to the plane containing the field and photon wavevector $`𝐤`$, and the ordinary (O) mode, which is mostly polarized parallel to the $`𝐤𝐁`$ plane (e.g., Meszaros 1992). The character of these modes is a strong function of the angle between $`𝐤`$ and $`𝐁`$, the plasma density and field strength. In very strong fields, the vacuum electron-positron pairs are polarized (see Section 4.2) and this effect dominates the photon propagation modes, which become almost purely linear (in contrast to the plasma-dominated modes that are elliptical).

In neutron star magnetic fields, the cyclotron decay rate is high enough so that nearly all particles occupy the ground state. Radiative (rather than collisional) processes thus control the Landau state populations. This means that resonant scattering (where a photon excites an electron to a higher Landau state from which it then spontaneously decays) dominates over absorption (where a photons excites an electron to a higher Landau state from which it collisionally de-excites).

## 3 Processes Modified by Strong Fields

### 3.1 Cyclotron/Synchrotron radiation

At low field strengths, cyclotron radiation (and synchrotron radiation at relativistic energies) is a process in which a charged particle radiates as it spirals around the magnetic field. At high fields, the particle motion perpendicular to the field is quantized and the emission must be described as transitions between Landau states. When the cyclotron energy, and thus the energies of the emitted photon, is a significant fraction of the particle rest mass (in fields $`B\stackrel{>}{}0.1B\mathrm{cr}`$), the particle recoil is important. The classical synchrotron emission formula will violate conservation of energy when the critical frequency, $`ϵ_{\mathrm{cr}}=(3/2)\gamma ^2B^{}\mathrm{sin}\psi `$, exceeds the particle’s kinetic energy, or when

$$\gamma ^2B^{}\mathrm{sin}\psi >\left(\gamma 1\right)$$
(2)
where $`\psi `$ is the pitch angle. Modified formulae (i.e. Sokolov & Ternov 1968, Harding & Preece 1987) for synchrotron radiation, based on QED, should be used in this regime. The quantum synchrotron spectra have a cutoff at the particle kinetic energy (see Figure 1), which limits the emission at high energies. When $`\mathrm{\Gamma }\gamma B^{}\mathrm{sin}\psi >0.1`$ there is a decrease in the energy loss rate to $`\dot{\gamma }\gamma ^{2/3}`$ (Erber et al. 1968), compared to the classical formula where $`\dot{\gamma }\gamma ^2`$. A more extreme spectral effect, the dominance of ground-state transitions, occurs for high values of $`\mathrm{\Gamma }\stackrel{>}{}1.0`$. In this regime, it is most likely for the particle in an excited Landau state to make a transition directly to the ground state, rather than the classical mode of decaying in a series of many single harmonic number transitions. This produces a “tip” at the upper end of the spectrum just before the cutoff, as is evident in the $`B^{}=10`$ case in Figure 1. Such a feature would be a signature of synchrotron radiation in very high magnetic fields, but would most likely be smeared out by variations in field strength and particle energy in a real source. For emission from low Landau states in strong fields, quantum cyclotron transition rates should be used and the spin of the particle is important in determining the spectrum.

### 3.2 Cyclotron absorption and scattering

Cyclotron absorption, the inverse of cyclotron emission, is a first-order process in which a photon excites a particle to a higher Landau state. For particles in the ground state, the required energy for excitation to state $`n`$ is (in rest mass units)

$$ϵ_n=[(1+2nB^{}\mathrm{sin}^2\theta )^{1/2}1]/\mathrm{sin}^2\theta $$
(3)
for a photon propagating at angle $`\theta `$ to the field (Daugherty & Ventura 1978). Because of the recoil of the particle, the cyclotron harmonics are actually anharmonic in high magnetic fields, so that the energy difference between successive harmonics decreases. As noted in Section 2, radiative transitions dominate over collisions in astrophysical sources having strong magnetic fields. A particle that absorbs a cyclotron photon will almost always de-excite by emitting another photon, rather than collisionally de-exciting, so that the end result of the process is a scattering of the photon rather than a true absorption. Thus in strong fields cyclotron absorption is most accurately treated as a second-order process, as resonances in the Compton scattering cross section, where the excited state of the particle is a virtual state for which energy and momentum are not strictly conserved. The cyclotron lines are thus broadened by the intrinsic width which is equal to the inverse of the decay rate from that state (Harding & Daugherty 1991), although Doppler broadening usually dominates at most angles. The cross section for Compton scattering in a magnetic field was first studied in the non-relativistic limit by Canuto, Lodenquai & Ruderman (1971) and the full QED cross section has been computed by Herold (1979), Daugherty & Harding (1986) and Bussard, Meszaros & Alexander (1986). Since the non-relativistic treatment is limited to dipole radiation, only scattering at the cyclotron fundamental is allowed. In the relativistic (QED) treatment scattering at higher harmonics is allowed, including Raman scattering, in which the state of the particle after scattering is higher than the initial state.

#### 3.2.1 Electron cyclotron lines

When continuum radiation is transmitted through a thermal plasma before reaching an observer, the radiation at the cyclotron resonant frequencies is scattered most strongly (and out of the line-of-sight), forming a distinctive series of line features. In most neutron star fields, the electron cyclotron features occur in the hard X-ray region and can be uniquely identified because of their (near) harmonic structure and also because this spectral region is absent of atomic and nuclear lines. Cyclotron lines are thus an excellent diagnostic tool for providing both a direct in-situ measurement of magnetic field strength and physical conditions of the radiating plasma. An electron cyclotron line at 35 keV was first discovered in the accreting pulsar Her X-1 (Trumper et al. 1978) and such lines have since been detected in about a dozen accreting X-ray pulsars (see Orlandini & dal Fiume 2001, for review). Magnetic fields have thus been measured for these sources in the range of $`10^{12}10^{13}`$ G. The most recent and sensitive detections of cyclotron lines have been made by BeppoSax, which detected four features (presumably the fundamental and three harmonics) in the spectrum of X0115+63 (Santangelo et al. 1999). Figure 2 shows this detected spectrum along with recent theoretical models of cyclotron resonance lines, illustrating the diagnostic power of such spectral features. By comparing the observed spectra to even the simplest homogeneous slab and cylinder models of the scattering plasma (where the magnetic field is either perpendicular or parallel to the direction of highest optical depth, respectively) it is possible to identify both the geometry of the plasma (probably slab) and the angle of the observer to the field (about $`50^070^0`$).

As was noted in Section 2, the interaction with the radiation field breaks the spin degeneracy of the excited Landau states, causing a splitting of the levels. The level splitting is quite small (Pavlov et al. 1991),

$$\mathrm{\Delta }\epsilon _{sf}\frac{\alpha }{2\pi }B^{}10^3\epsilon _{cyc},$$
(4)
$`\alpha `$ being the fine-structure constant, and produces a fine structure of the cyclotron harmonic features that would be potentially observable with very high resolution. There are also electron spin-flip transitions possible in the same Landau state which would produce lines at the above energies, in the 1-10 keV range for magnetar fields. However, the rate of such transitions (Parle 1996, Geprags 1994), $`\mathrm{\Gamma }_{sf}\alpha ^6B^3`$, is extremely small, being a third-order process, so that these lines would not be detectable.

#### 3.2.2 Proton cyclotron lines

In the very high fields of SGRs and AXPs, the electron cyclotron fundamental would fall well above 1 MeV and therefore above the observed continuum. However the proton cyclotron fundamental, $`ϵ_{\mathrm{cyc}}=6.3(B/10^{15}`$G) keV would fall right in the range of the observed spectra. Proton cyclotron features have been predicted in a number of spectral models (Ho & Lai 2001, Zane et al. 2001, Ozel 2002, Lai & Ho 2002). Very recently, absorption-like features have been detected around 5 keV in the burst spectra of SGR 1806-20 (Ibrabim et al. 2002) which could be proton cyclotron lines if the surface field is really as high as $`10^{15}`$ G. A feature in the spectrum of the AXP 1RXS J170849-400910 at around 8 keV has also been reported (Rea et al. 2003), again implying a magnetic field around $`10^{15}`$ G if the feature is due to proton cyclotron scattering. One interesting result of the theoretical models of radiation transfer in very high magnetic fields is that vacuum polarization (see Section 4.2) effectively suppresses the strength of the proton cylotron features which would otherwise be much deeper than the observed spectral features.

### 3.3 Atomic processes

When the Lorentz force on electrons in an atom is comparable to or greater than the Coulomb force, which occurs for fields

$$\frac{B^{}mc^2}{2Ry}>1B>2\times 10^9G,$$
(5)
where $`Ry=`$ is the Rydberg constant, the atom is distorted into a cylindrical shape and the electron energy levels are closer to Landau states than to atomic states. Atomic transitions thus become more cyclotron-like. Among the many interesting consequences (see Wunner & Herold 1986, for review) are an increase in the binding energy of the atom. For example, the ionization potential of Hydrogen,

$$V_o\frac{e^2}{\text{}\lambda }\left(\frac{B}{B_{cr}}\right)^{1/2}=3.8keV\left(\frac{B}{B_{cr}}\right)^{1/2}$$
(6)
moves from the UV into the X-ray range in magnetar fields. Atoms are thus more tightly bound so that the ionization fraction is lowered by the magnetic field. Even in the relatively low fields ($`10^8`$ G) of some White Dwarfs, atomic emission and absorption lines can provide an accurate measure of the magnetic field of these objects (Angel et al. 1985).

### 3.4 Two-photon pair production and annihilation

The second-order process in which two photons of energy $`ϵ_1`$ and $`ϵ_2`$ interact to produce an electron-positron pair, and its inverse process, pair annihilation into two photons, are modified in magnetic fields $`B\stackrel{>}{}10^{13}`$ G. Due to the relaxation of momentum conservation perpendicular to the field, the threshold condition (Daughterty & Bussard 1980) for two-photon pair production becomes

$$(ϵ_1\mathrm{sin}\theta _1+ϵ_2\mathrm{sin}\theta _2)^2+2ϵ_1ϵ_2[1\mathrm{cos}(\theta _1\theta _2)]4,$$
(7)
where $`\theta _1`$ and $`\theta _2`$ are the angles the photons are propagating with respect to the field. The free-space threshold condition is obtained by dropping the first term on the left-hand side. It is thus possible for two photons traveling parallel to each other, and also at an angle to the magnetic field, to produce a pair. The cross section is modified by the requirement that each member of the pair must occupy Landau states, producing a series of spikes at the thresholds for each new state combination (Kozlenkov & Mitrofanov 1987).

The total cross section for the inverse process, two-photon annihilation, begins to decrease around $`10^{13}`$ G (Bussard & Daugherty 1980). As the cross section decreases, the energy of a pair at rest, which is split equally between the two photons at low and zero field strength, are divided unequally so that one of the photons gets almost all of the initial pair energy. Also at $`B>10^{13}`$ G, the cross section for the first order process of one-photon annhilation (see Section 4.1) is increasing exponentially with increasing field strength. Thus two photon annihilation seems to transform into one-photon annhilation in magnetic fields above $`10^{13}`$ G. The spectral lines resulting from two-photon pair annihilation are therefore quite asymmetric in high magnetic fields (Harding & Baring 1992). The annihilation line begins to broaden in fields $`B^{}\stackrel{>}{}0.1`$ and at fields $`B^{}>0.3`$ the lines have become nearly a continuum below the kinematic cutoff at 1 MeV. Thus the characteristic signature of pair annihilation at low fields, a line at 511 keV, is much more difficult to detect at high field strengths.

## 4 Processes Unique to Strong Fields

In addition to the many changes and transformations that occur to well-known processes, strong magnetic fields allow new and exotic processes that occur only in strong field environments.

### 4.1 One-photon pair production and annihilation

The first-order processes of one-photon pair production and its inverse, one-photon pair annihilation, do not take place in field-free environments because energy and momentum cannot be simultaneously conserved. However, magnetic fields that are strong enough (i.e. near or above $`B_{\mathrm{cr}}`$) can absorb/supply the extra momentum that a photon must have to create or be created by a pair. The threshold for one-photon pair production is $`ϵ_{1\gamma }=2/\mathrm{sin}\theta `$, where $`\theta `$ is the angle between the photon momentum and the magnetic field. Although the threshold is purely kinematic, the attenuation coefficient is exponentially dependent on field strength (Erber 1966) so that the processes effectively does not take place unless $`B^{}ϵ\mathrm{sin}\theta \stackrel{>}{}0.2`$. In neutron star fields $`B^{}>0.1`$, the attenuation coeffient is so high that pair production occurs at threshold (Harding & Daugherty 1983), so that the quantum effects resulting from discreteness of the pair states causes sharp spikes in the attentuation coefficient at thresholds of pair states.

One-photon pair production has become a critically important process in models of isolated pulsars, where acceleration of particles to energies of 1-10 TeV results in the radiation of photons well above pair threshold in very strong fields. However, these photons initially have such small angles to the field ($`\theta 1/\gamma `$, where $`\gamma `$ is the Lorentz factor of the particle) that they are below threshold and must travel a distance through the curved magnetic field to attain large enough angles. Since the strength of a dipole (or multipole) field falls quickly as the photon travels away from the neutron star surface, many escape without pair producing. Observed spectra that originate at radius $`r`$ near the poles of a neutron star dipole field will thus have high-energy cutoffs at the pair escape energy (Harding et al. 1997, Harding 2002),

$$E_c2\mathrm{GeV}P^{1/2}\left(\frac{r}{R}\right)^{1/2}\mathrm{max}\{0.1,B_{0,12}^1\left(\frac{r}{R}\right)^3\}$$
(8)
where $`B_{0,12}B_0/10^{12}`$ G, and $`P`$ and $`R`$ are the period and radius of the neutron star. The field strength dependence of the energy of these cutoffs, which may have been detected in the spectra of $`\gamma `$-ray pulsars (Harding 2001, Thompson 2001), could provide a measure of neutron star magnetic fields and test high-energy pulsar emission models.

Production of electron-positron pairs has been an important element in pulsar models, since it is believed that pairs may be necessary for coherent radio emission (e.g. Melrose 2000). Much theoretical work has been devoted to understanding the absence of radio pulsars detected with long periods, defining a “death line” in $`\dot{P}`$-$`P`$ space, in terms of the inability of very old pulsars to create pairs (e.g. Harding et al. 2002, and references therein).

### 4.2 Vacuum polarization

Very strong fields are capable of polarizing the vacuum electron-positron pairs. This affects the photon progagation modes at high frequencies. The vacuum dominates the propagation modes at photon energies

$$E_V\stackrel{>}{}1\mathrm{k}\mathrm{e}\mathrm{V}\left(\frac{Y_e\rho }{1\mathrm{g}\mathrm{c}\mathrm{m}^3}\right)^{1/2}\left(\frac{B}{10^{14}G}\right)^1$$
(9)
(Ventura, Nagel, & Meszaros 1979) where $`Y_e=Z/A`$ is the electron fraction and $`\rho `$ is the plasma density. Vacuum polarization causes a number of interesting radiative transfer effects with possibly observable signatures.

When the condition of Eqn (9) above is satisfied, the photon progation modes are determined by the vacuum polarization, where the modes are predominantly linear, rather than by the plasma, where the modes are elliptical. The photon modes in the magnetized vacuum are thus linear combinations of the plasma modes. Photons propagating through either density or magnetic field gradients of a neutron star magnetosphere may undergo adiabatic evolution of their normal modes (Heyl & Shaviv 2000, Lai & Ho 2002, Ozel 2001, 2002). For example a photon that is in purely O-mode at relatively high density, where plasma modes dominate, may slowly evolve into the X-mode if it crosses the vacuum resonance and moves into a region where the vacuum modes dominate. Since the scattering opacities of the O and X modes are very different, this vacuum-induced mode conversion has significant transfer effects for X-ray photons in neutron star, and especially magnetar, magnetospheres.

Another transfer effect caused by vacuum polarization is enhanced mode coupling during scattering (Bulik & Miller 1997, Ozel 2002). When photons undergo scattering in a magnetized plasma, they can change their polarization mode, as well as their energy and direction. When scattering occurs near the vacuum resonance, the partial cross section for scattering into another mode (e.g. Meszaros & Ventura 1979) is greatly enhanced. This effect can significantly change the emergent spectrum and angular distribution of radiation from a neutron star by converting photons from modes with high/low opacity into modes with low/high opacity. Both the effects of adiabatic mode conversion (without scattering) and mode switching (with scattering) have been found to greatly suppress the strength of ion-cyclotron features near the vacuum resonance, by effectively reducing the opacity of the resonant photons.

### 4.3 Photon splitting

Photon splitting is a third-order process in which a photon splits into two photons of lower energy (Adler 1971). It cannot occur in field-free environments but the cross section is significant in high fields, where it is able to compete with one-photon pair production. Because photon splitting has no threshold, high-energy photons propagating at very small angles to the magnetic field in neutron star magnetospheres may split before reaching the threshold for pair production (Harding et al. 1997). Generally, splitting dominates over pair production in fields $`B^{}\stackrel{>}{}1.0`$. Thus it may be important in emission models of high-field radio pulsars (Harding et al. 1997) and magnetars (Thompson & Duncan 1995, Harding & Baring 1996, Zhang & Harding 2001).

Polarization signatures due to photon splitting may be observable in the spectra of some sources since kinematic selection rules in the low dispersion (low B) limit (Adler 1971) only allows splitting of photons with electric vectors perpendicular to the field (X mode) into two photons with parallel (O) mode polarization. As a result, parallel-mode photons are attenuated only by pair production while perpendicular-mode photons are attenuated mostly by splitting. Observed spectra will have polizarization-dependent cutoffs, since the escape energy for photon splitting is lower than the escape energy for pair production, causing the high end of the spectrum to be 100% polarized in the parallel mode (see Figure 4).

Since photon splitting can dominate, at least in the X-mode, it has the ability to suppress the production of pairs in neutron star magnetospheres. Baring & Harding (1998) suggested that the radio-quiescence of SGRs and AXPs could be explained by pair suppression in fields $`B^{}\stackrel{>}{}1`$, if pairs are required for detectable radio emission. However, complete suppression of pair creation requires that photon splitting be allowed for both photon polarization modes (Baring & Harding 2001), and thus a break-down of the kinematic selection rules. But recently, Usov (2002) determined through analysis of the linear vacuum polarization modes that Adler’s kinematic selection rules are valid to arbitrarily high field strengths. Thus, unless quadratic vacuum polarization contributions allow a violation of these selection rules, it seems that the parallel mode photons may produce pairs even in magnetar fields and some other reason for their radio quiecsence must emerge. It is possible that the parallel-mode photons produce bound pairs (Shabad & Usov 1982) rather than free pairs in high magnetic fields. If these bound pairs are stable to dissociation by either the radiation field or a strong electric field (Usov & Melrose 1995), then complete pair suppression would occur, at least in the strong-field region near the surface.

While the theoretical jury is still out on whether photon splitting (and bound-pair creation) can cause radio quiescence, the Parkes Multibeam radio pulsar survey (Manchester et al. 2001) is detecting a number of pulsars with very high magnetic fields (as inferred from dipole spin down) near and above the critical field. Several of these pulsars lie above the photon splitting death line in $`\dot{P}`$-$`P`$ for surface emission (Baring & Hardng 1998). However, Zhang & Harding (2002) have noted that in space-charge limited flow acceleration models, particles can keep accelerating to higher altitudes where fields are low enough for pair production to dominate over photon splitting. This will move the death line to higher $`\dot{P}`$.

## 5 Neutron star magnetic fields: Is there a limit?

The strength of an electric field is limited by the “Klein Catastrophe”, which occurs when the field is strong enough to cause a breakdown of the vacuum into free electron-positron pairs. A magnetic field is stable to vacuum breakdown because it cannot accelerate the particles in a vacuum pair (magnetic fields do no work). But there are several other physical limits that apply. Neutron stars cannot have magnetic fields with pressure high enough to blow them apart. This gives a limit

$$\frac{B^2R^3}{8\pi }<\frac{GM^2}{R}B<10^{18}\mathrm{G}.$$
(10)
There are also proposed particle physics limits to the strength of neutron star fields. Various interactions with photons and matter deplete energy and momentum from the field, limiting its strength to $`B<10^{16}10^{18}`$ G (Lerche & Shramm 1977). An ultimate limit of $`B<10^{20}`$ G may come from spontaneous production of bound pairs (Zaumen 1976).