text:
# Time-symmetric initial data for binary black holes in numerical relativity

## I Motivation and relation to other works

The numerical computation of the collision of two black holes is of paramount importance for the observation of gravitational waves by the network of laser-interferometric detectors. When investigating this problem the ten Einstein field equations are separated into: (i) four constraint equations that are to be satisfied by some initial data given on an initial three-dimensional Cauchy hypersurface; (ii) six hyperbolic-like equations describing the dynamical evolution of the gravitational field on neighbouring hypersurfaces. The Bianchi identities guarantee that the constraint equations are satisfied on neighbouring hypersurfaces if they are on the initial hypersurface. There are infinitely many ways that the initial data can be chosen to represent the starting state of the evolution of black holes. It is widely admitted that the problem of choosing physically realistic initial conditions for the collision of two black holes has not yet been solved. There has been a lot of concern in the litterature cook  for knowing what would really motivate physically a particular choice of initial data.

Let us consider the problem of time-symmetric initial data, which are physically appropriate to two momentarily static black holes, i.e. with zero initial velocities. The dynamical evolution of time-symmetric data describes the subsequent head-on collision of the two black holes. In this situation the constraint equations reduce to the Hamiltonian or scalar contraint equation $`R=0`$ (in vacuum — the case appropriate to black holes), with $`R`$ being the three-dimensional scalar curvature. Considering as usual lichne ; brill ; york  a conformal decomposition of the spatial metric (spatial indices $`i,j,\mathrm{}=1,2,3`$),

$$\gamma _{ij}=\mathrm{\Psi }^4\stackrel{~}{\gamma }_{ij},$$
(1)
where $`\gamma _{ij}`$ is the physical metric and $`\stackrel{~}{\gamma }_{ij}`$ denotes the conformal (unphysical) metric, we obtain the Lichnerowicz lichne  equation, which is an elliptic-type equation to be satisfied by the conformal factor $`\mathrm{\Psi }`$. In the time-symmetric case that equation becomes

$$\stackrel{~}{\mathrm{\Delta }}\mathrm{\Psi }=\frac{\mathrm{\Psi }}{8}\stackrel{~}{R}.$$
(2)
The scalar curvature $`\stackrel{~}{R}`$ and Laplacian $`\stackrel{~}{\mathrm{\Delta }}`$ are the ones associated with the conformal metric $`\stackrel{~}{\gamma }_{ij}`$. The interest of the conformal decomposition is that by solving Eq. (2) we generate a physical solution $`\gamma _{ij}`$ of the constraint equation starting from any choice for the conformal metric $`\stackrel{~}{\gamma }_{ij}`$ (a priori). Therefore the problem of initial conditions resides in the choice of a physically well-motivated conformal metric $`\stackrel{~}{\gamma }_{ij}`$.

The simplest choice of initial conditions (one motivated by simplicity rather than by physics) is the one for which $`\stackrel{~}{\gamma }_{ij}`$ equals the flat metric $`\delta _{ij}`$. In that case Eq. (2) reduces to the flat-space Laplace equation. Some exact solutions, appropriate to (momentarily static) black holes, have been obtained by Misner misner  and Lindquist lindquist , and Brill and Lindquist BL . The solution of Brill and Lindquist \[for which the conformal factor $`\mathrm{\Psi }^{\mathrm{BL}}`$ takes the form of Eq. (12b) below\] is particularly interesting: it describes the “geometrostatics” of two black holes, consisting of three asymptotically flat regions connected by two Einstein-Rosen bridges (actually, the solution is known for $`N`$ black holes). The one region containing the two throats is supposed to represent our universe, while the two sheets expanding behind it are associated with the two black holes. The beauty of the Brill-Lindquist solution is that not only the total ADM mass-energy $`M`$ of the binary, but also the two individual masses $`m_1`$ and $`m_2`$ of the black holes, are computed “at infinity”. One can use for instance standard surface integrals extending on topological two-spheres at infinity. The binary’s geometrostatic energy (i.e. the gravitational interacting or binding energy, in the center-of-mass frame) is therefore computed unambiguously as $`\frac{E}{c^2}=Mm_1m_2`$.

Yet the solution of Brill and Lindquist, despite its undeniable interest, is not “physically realistic” in the sense that it differs, in the limit $`c+\mathrm{}`$, from the three-metric found for a post-Newtonian solution. Indeed the post-Newtonian metric generated by two point-particles is known to deviate from conformal flatness at the 2PN order (see e.g. Ref. BFP98 ) 111In this paper $`c`$ and $`G`$ denote the speed of light and the gravitational constant. As usual the $`n`$PN order means the terms of order $`\frac{1}{c^{2n}}`$ when $`c+\mathrm{}`$.. In consequence both the Brill-Lindquist metric and the associated binding energy $`E`$ disagree with the post-Newtonian results from the 2PN order.

A compelling motivation for constructing physically realistic initial data is the agreement with the post-Newtonian approximation when $`c+\mathrm{}`$. Recall that the post-Newtonian theory provides some explicit expressions for the metric, equations of motion and energy of binary systems of point-particles. The post-Newtonian metric is valid in the binary’s “near-zone” (size of near-zone is much less than a gravitational wavelength), but has been proved to come from the re-expansion when $`c+\mathrm{}`$ of a “global” (post-Minkowskian-type) solution, defined everywhere in space-time including the wave zone livrev . Furthermore, in the post-Newtonian approach the modelling of black holes by point-like particles — i.e. technically by Dirac delta-functions in the stress-energy tensor — is rather well justified. We shall provide below some further evidence that the “post-Newtonian masses” of point-particles are indeed identical to the black-hole masses 222See also, in a similar context, the matching of a 1PN solution for the orbital motion of point-particles to two perturbed Schwarzschild black holes alvi ..

It has been suggested diener  (see also BCLT ) that, in order of taking into account the post-Newtonian physics into the initial data, one should adopt for the conformal metric $`\stackrel{~}{\gamma }_{ij}`$ directly a post-Newtonian solution. In such a proposal one expects that by correcting the post-Newtonian solution by means of a conformal factor $`\mathrm{\Psi }^4`$ (computed numerically), one will somehow be able to “compensate” in the physical metric $`\gamma _{ij}`$ (i.e., in fact, to cancel out) the systematic higher-order post-Newtonian error terms that are neglected in $`\stackrel{~}{\gamma }_{ij}`$.

There has been other proposals for realistic initial data, in particular built on the relaxation of the assumption of conformal flatness cook . One of these is to use for $`\stackrel{~}{\gamma }_{ij}`$ a linear combination of two (boosted versions of) Kerr-Schild metrics matzner ; BMNM . It is not known if the physical metric which is generated in this way from the numerically-computed conformal factor, is consistent with post-Newtonian (e.g. 2PN) calculations.

On the other hand, we should remark that to which extent the hypothesis of conformal flatness introduces some unphysical spurious (and numerically important) effects remains an issue. Indeed, relaxing this hypothesis may not always be a panacea. A quite different idea for settling the initial conditions of two black holes is to solve numerically a subset of the Einstein field equations (the four constraint equations plus one evolution equation), under the two premices of conformal flatness and the existence of a helical Killing symmetry GGB1 ; GGB2 . Despite the conformal flatness of the spatial metric, the latter calculation gives results in good agreement with post-Newtonian predictions GGB2 ; B02ico ; DGG02 . Nevertheless, we think that it is important, at some stage, to get rid of the hypothesis of conformal flatness.

We emphasize furthermore that the agreement between numerical relativity GGB1 ; GGB2  and post-Newtonian theory GGB2 ; B02ico ; DGG02  holds up to the very relativistic regime of the innermost circular orbit (ICO), where the orbital velocities are of the order of 50% of the speed of light 333The ICO is defined by the minimum of the binary’s energy function for circular orbits. It represents a useful reference point for the numerical GGB1 ; GGB2  and 3PN B02ico  calculations because for both of them the ICO is well-defined and can be meaningfully compared. However, the radiation reaction terms are neglected in its definition, so the ICO probably does not have a rigorous physical meaning in a context of exact radiative solutions.. In Refs. B02ico  it is suggested that the result for the ICO of two black holes with comparable masses, at the 3PN approximation 444We mean the standard Taylor-expanded form of the approximation — without using any post-Newtonian resummation techniques., is likely to be close to the “exact” solution, within 1% of fractional accuracy or better. This constitutes a motivation for advocating that the black-hole initial conditions, which are to be set (presumably) around the location of the ICO, should be in agreement with post-Newtonian theory.

In the present paper we propose an alternative way for incorporating the post-Newtonian information into the initial data of black-hole binaries (in the time-symmetric case). We find, in Section II, a simple expression for the conformal metric $`\stackrel{~}{\gamma }_{ij}`$, which is such that the corresponding physical metric $`\gamma _{ij}`$ is isometric (i.e. differs by a coordinate transformation) to the standard post-Newtonian spatial metric at the 2PN order in the limit $`c+\mathrm{}`$. At the same time, the solution is defined globally in space, with a global structure similar to the one of Brill and Lindquist. Our solution is not “exact”, but exists as a certain non-linear perturbation, investigated in Section III, of the Brill-Lindquist solution, playing here the role of a “background” metric. Most importantly, in Section IV we investigate the asymptotic structure of the solution, and compute “geometrically” the masses $`M`$ and $`m_1`$, $`m_2`$, i.e. by surface integrals performed in their respective domains at infinity. In Section V the binary’s interacting energy, deduced from the previous masses, is proved to be in agreement with the known post-Newtonian energy up to the 2PN order. (We shall find, however, that the definition we adopt for the two individual masses $`m_1`$ and $`m_2`$ makes sense only up to the 2PN order.)

## II Conformal decomposition of the spatial metric

### II.1 Definition of the conformal metric

The conformal metric we propose, appropriate to two black holes of masses $`m_1`$ and $`m_2`$ (to agree later with the “geometrical” masses), located at the singular points $`𝐲_1`$ and $`𝐲_2`$, and momentarily at rest ($`𝐯_1=𝐯_2=\mathrm{𝟎}`$), is

$$\stackrel{~}{\gamma }_{ij}=\delta _{ij}\frac{8G^2m_1m_2}{c^4}\frac{^2g}{y_1^{}}.$$
(3)
The function $`g`$ introduced here represents an elementary “kernel” playing an important role in post-Newtonian calculations S87 ; DI91 ; BDI95 ; WWi96 ; BFP98 ; AF97 . It depends on the “field” point $`𝐱`$ on the one hand, and on the pair of “source” points $`𝐲_1,𝐲_2`$ on the other hand; it is defined by

$$g(𝐱;𝐲_1,𝐲_2)=\mathrm{ln}\left(r_1+r_2+r_{12}\right),$$
(4)
where $`r_1=|𝐱𝐲_1|`$ and $`r_2=|𝐱𝐲_2|`$ are the distances to the black holes, and $`r_{12}=|𝐲_1𝐲_2|`$ is their separation. The function $`g`$ satisfies, in the sense of distributions (i.e. for all values of $`𝐱`$, including the singular points $`𝐲_1`$ and $`𝐲_2`$),

$$\mathrm{\Delta }g=\frac{1}{r_1r_2},$$
(5)
where $`\mathrm{\Delta }`$ denotes the usual flat-space Laplacian with respect to the field point $`𝐱`$. In Eq. (3) the derivatives are taken with respect to the two source points 555The following explicit formula holds:
$${}_{i}{}^{}g_{j}^{}\frac{^2g}{y_1^iy_2^j}=\frac{n_{12}^in_{12}^j\delta ^{ij}}{r_{12}(r_1+r_2+r_{12})}+\frac{(n_{12}^in_1^i)(n_{12}^j+n_2^j)}{(r_1+r_2+r_{12})^2}.$$
Here, $`n_1^i=(x^iy_1^i)/r_1`$, $`n_2^i=(x^iy_2^i)/r_2`$ and $`n_{12}^i=(y_1^iy_2^i)/r_{12}`$. See Ref. BFP98  for further discussion and formulas concerning the function $`g`$.. The carets around the indices refer to the symmetric and trace-free (STF) projection: $`T_{<ij>}\frac{1}{2}\left(T_{ij}+T_{ji}\right)\frac{1}{3}\delta _{ij}T_{kk}`$; so the trace of the metric (3) is normalized to be $`\stackrel{~}{\gamma }_{ii}=3`$. It is of course quite natural (and useful in practice) to impose that the deviation of the conformal metric from flat space be trace-free.

Our proposal is to generate, by means of numerical techniques (i.e. with the help of elliptic solvers), a conformal factor $`\mathrm{\Psi }`$ solving the constraint equation (2) corresponding to the particular choice of conformal metric (3). The metric $`\gamma _{ij}=\mathrm{\Psi }^4\stackrel{~}{\gamma }_{ij}`$ we obtain in this way will incorporate the post-Newtonian (2PN) physics of the initial state of the head-on collision of two black holes 666See Refs. Smarr ; AHSSS  for numerical calculations of the head-on collision of black holes; see also Ref. SPW95  for a post-Newtonian calculation.. The latter assertion will now be proved, for the rest of the paper, with the help of analytic perturbation methods.

### II.2 Relation with the post-Newtonian metric

The first result shows that the “near-zone” behaviour of the solution (i.e. $`r/c0`$) is physically sound.

###### Theorem 1

The post-Newtonian expansion (when $`c+\mathrm{}`$) of the solution $`\gamma _{ij}=\mathrm{\Psi }^4\stackrel{~}{\gamma }_{ij}`$ of the constraint equation $`(\text{2})`$, where $`\stackrel{~}{\gamma }_{ij}`$ is defined by Eq. $`(\text{3})`$, differs from the standard post-Newtonian spatial metric, calculated by standard post-Newtonian methods, by a mere change of coordinates at the $`2\mathrm{P}\mathrm{N}`$ order, i.e.

$$\gamma _{ij}=\mathrm{g}_{ij}^{2\mathrm{P}\mathrm{N}}{}_{|_{𝐯_1=𝐯_2=\mathrm{𝟎}}}{}^{}+_i\xi _j+_j\xi _i+𝒪\left(\frac{1}{c^6}\right).$$
(6)
To be precise, by $`\mathrm{g}_{ij}^{2\mathrm{P}\mathrm{N}}`$ we mean the spatial metric in harmonic coordinates, when truncated at the 2PN order, that is given by Eq. (7.2c) in Ref. BFP98 . As indicated in Eq. (6), we must set the particles’ velocities $`𝐯_1`$ and $`𝐯_2`$ to zero in the post-Newtonian metric in order to conform with the assumption of time-symmetry. The remainder $`𝒪\left(\frac{1}{c^6}\right)`$ stands for the neglected 3PN and higher-order terms.

The proof of Theorem 1 is easily achieved on the basis of a post-Newtonian iteration of Eq. (2). At the 2PN order, this equation becomes

$$\mathrm{\Delta }\mathrm{\Psi }=\frac{G^2m_1m_2}{c^4}_{ij}\left({}_{}^{}\right)+𝒪\left(\frac{1}{c^6}\right),$$
(7)
where we denote $`{}_{i}{}^{}g_{j}^{}\frac{^2g}{y_1^iy_2^j}`$. The most general solution reads as

$$\mathrm{\Psi }=\psi \frac{G^2m_1m_2}{2c^4}\mathrm{D}\left(\frac{g}{3}+\frac{r_1+r_2}{2r_{12}}\right)+𝒪\left(\frac{1}{c^6}\right),$$
(8)
where $`\psi `$ represents a solution of the homogeneous equation (i.e. $`\mathrm{\Delta }\psi =0`$) 777We employ the notational shorthand $`\mathrm{D}\frac{^2}{y_1^iy_2^i}`$ (so that $`\mathrm{D}g={}_{i}{}^{}g_{i}^{}`$). Notice the useful relations    $`\mathrm{D}g`$ $`=`$ $`{\displaystyle \frac{1}{2r_1r_2}}{\displaystyle \frac{1}{2r_1r_{12}}}{\displaystyle \frac{1}{2r_2r_{12}}},`$    $`_{ij}\left({}_{i}{}^{}g_{j}^{}\right)`$ $`=`$ $`\mathrm{D}\left[{\displaystyle \frac{1}{2r_1r_2}}+{\displaystyle \frac{1}{2r_1r_{12}}}+{\displaystyle \frac{1}{2r_2r_{12}}}\right].`$  . By comparing with the post-Newtonian metric (see Eq. (7.2c) in Ref. BFP98  in which $`𝐯_1=𝐯_2=\mathrm{𝟎}`$), we readily find that the latter solution is uniquely specified, at the 2PN order, as being

$$\psi =1+\frac{Gm_1}{2c^2r_1}\left(1\frac{Gm_2}{2c^2r_{12}}\right)+\frac{Gm_2}{2c^2r_2}\left(1\frac{Gm_1}{2c^2r_{12}}\right)+𝒪\left(\frac{1}{c^6}\right).$$
(9)
Also uniquely determined is the expression of the vector $`\xi _i`$ in Eq. (6), which represents an infinitesimal gauge transformation. At the 2PN order it takes the expression

$$\xi _i=\frac{G^2m_1^2}{4c^4}\frac{n_1^i}{r_1}+\frac{G^2m_2^2}{4c^4}\frac{n_2^i}{r_2}+𝒪\left(\frac{1}{c^6}\right),$$
(10)
where $`𝐧_1(𝐱𝐲_1)/r_1`$ and $`𝐧_2(𝐱𝐲_2)/r_2`$. The system of spatial coordinates employed in the present paper, i.e. which corresponds to a conformal metric of the particular form displayed by Eq. (3), is thereby determined, and related up to the 2PN order to the harmonic coordinate system used in Ref. BFP98 , by the coordinate change $`\delta x^i=\xi ^i`$ (with $`\xi ^i=\delta ^{ij}\xi _j`$).

Summarizing this section, we have found a conformal decomposition of the metric $`\gamma _{ij}`$ that is in agreement with the post-Newtonian metric up to the 2PN order. Notice, however, that the conformal metric we propose is not unique, because we can always add higher-order terms (e.g. 3PN) to Eq. (3). A contrario, this means that with the metric (3) one should not expect Theorem 1 to hold at the 3PN order and beyond. At the 3PN order for instance the conformal metric will be more complicated than the simple expression (3).

## III Perturbation of the Brill-Lindquist solution

Having checked the near-zone structure of the metric, let us investigate some of its global properties, when it is viewed as a solution of the constraint equation (2) that is a priori valid everywhere on a space-like hypersurface. For this purpose, we impose that $`\gamma _{ij}`$ is a non-linear perturbation of the exact solution of Brill and Lindquist BL  — which is conformally flat ($`\stackrel{~}{\gamma }_{ij}^{\mathrm{BL}}=\delta _{ij}`$). This will mean that the topology of our solution is identical to the topology of the Brill-Lindquist solution, i.e. be “three-sheeted”, in contrast with the two-sheeted topology of the Misner-Lindquist solution misner ; lindquist .

### III.1 Hierarchy of perturbation equations

Let us write Eq. (3) in a more transparent form, $`\stackrel{~}{\gamma }_{ij}=\delta _{ij}+\beta s_{ij}`$, where

$`\beta `$ $``$ $`{\displaystyle \frac{8G^2m_1m_2}{c^4}},`$ (11a)
$`s_{ij}`$ $``$ $`{\displaystyle \frac{^2g}{y_1^{<i}y_2^{j>}}}\text{(such that }s_{ii}=0\text{)}.`$ (11b)
It is helpful to view the non-linear perturbation we want to consider, as being generated by the “seed” or “generating” function $`s_{ij}`$, and to interpret the parameter $`\beta `$ as the magnitude of that perturbation. In the following, $`\beta `$ will play the role of a “book-keeping” parameter allowing us to label the successive non-linear perturbation orders. Besides the non-conformally flat piece of the metric brought about by $`s_{ij}`$, it is evident that we must also introduce a perturbation in the conformal factor. We pose

$`\mathrm{\Psi }`$ $`=`$ $`\psi +\sigma ,`$ (12a)
$`\psi `$ $`=`$ $`1+{\displaystyle \frac{\alpha _1}{r_1}}+{\displaystyle \frac{\alpha _2}{r_2}},`$ (12b)
where $`\psi `$ denotes the Brill-Lindquist conformal factor BL , parametrized by two constants $`\alpha _1`$ and $`\alpha _2`$ (in the case of two particles). The quantity $`\sigma `$ denotes a certain perturbation of $`\psi `$.

With full generality — within the present perturbative framework —, we look for the expression of $`\sigma `$ in the form of an infinite power series in $`\beta `$, solving the equation (2):

$$\sigma =\underset{n=1}{\overset{+\mathrm{}}{}}\beta ^n\sigma _{(n)}.$$
(13)
We insist that in our approach the Brill-Lindquist solution plays the role of the background metric (it depends solely on $`\alpha _1`$ and $`\alpha _2`$), while the function $`s_{ij}`$ yields a non-linear deformation of this background, perturbatively ordered by the book-keeping parameter $`\beta `$. For a given choice of $`s_{ij}`$, we expect that the resulting solution is unique (at least in a sense of formal power series in $`\beta `$). Because $`\beta `$ involves two mass factors $`m_1`$ and $`m_2`$, so it is proportional to $`G^2`$, the perturbation series we look for is like a “double” post-Minkowskian expansion, going “twice as fast” as the usual post-Minkowskian expansion when $`G0`$. Considering this series on the point of view of a post-Newtonian re-expansion, we see that each non-linear order brings in a new factor $`1/c^4`$, so that our perturbation series can be said to go by — quite efficient indeed — steps of 2PN orders.

Let us compare our definitions (12)-(13) with the result (8)-(9) provided by the agreement with the 2PN metric. We find that the constants $`\alpha _1,\alpha _2`$ are determined with relative 1PN accuracy,

$$\alpha _1=\frac{Gm_1}{2c^2}\left[1\frac{Gm_2}{2r_{12}c^2}+𝒪\left(\frac{1}{c^4}\right)\right]\text{and }12,$$
(14)
and that with this accuracy they agree with the prediction of the Brill-Lindquist solution 888We recall that $`\alpha _1^{\mathrm{BL}}`$ and $`\alpha _2^{\mathrm{BL}}`$ in the case of the Brill-Lindquist solution are related to the masses by the exact relations    $`\alpha _1^{\mathrm{BL}}\left[1+{\displaystyle \frac{\alpha _2^{\mathrm{BL}}}{r_{12}}}\right]`$ $`=`$ $`{\displaystyle \frac{Gm_1}{2c^2}}\text{and }12,`$    $`\alpha _1^{\mathrm{BL}}+\alpha _2^{\mathrm{BL}}`$ $`=`$ $`{\displaystyle \frac{GM}{2c^2}}.`$  . In addition, we find that the perturbation $`\sigma `$ in the conformal factor is given by the second term in the right side of Eq. (8), which comes in only at the 2PN order — it is purely linear in $`\beta `$. Hence,

$$\sigma _{(1)}=\frac{1}{16}\mathrm{D}\left(\frac{g}{3}+\frac{r_1+r_2}{2r_{12}}\right)+𝒪\left(\frac{1}{c^2}\right),$$
(15)
with all higher-order $`\sigma _{(n)}`$’s being negligible with this approximation. Actually, Eqs. (14)-(15) correspond to a particular “sharing” of the terms between $`\alpha _1,\alpha _2`$ on the one hand, and $`\sigma _{(1)}`$ on the other hand, since we can always add to $`\sigma _{(1)}`$ some “homogeneous” terms $`1/r_1`$ and $`1/r_2`$ without changing the equation for $`\sigma _{(1)}`$. It is obvious that such a sharing is physically irrelevant. However we shall forbid a different sharing of terms by adopting below the prescription (which represents simply a convenient choice) that $`\sigma _{(1)}`$ and all subsequent iterations solve the equation for the conformal factor in the sense of distributions. Some results more complete than (14)-(15) will be obtained below.

By inserting the perturbation ansatz (13) into the constraint equation (2), and by identifying each of the coefficients of the successive powers of $`\beta `$ in both sides of the equation, we obtain a hierarchy (indexed by $`n`$) of Poisson-type equations:

$$\mathrm{\Delta }\sigma _{(n)}=\mathrm{\Sigma }_{(n)}[\sigma _{(1)},\mathrm{},\sigma _{(n1)}],$$
(16)
where we recall that $`\mathrm{\Delta }_i_i`$. The $`n`$th-order source term $`\mathrm{\Sigma }_{(n)}`$ depends on the solutions of the preceding iterations $`\sigma _{(1)},\mathrm{},\sigma _{(n1)}`$, on the “generating” function $`s_{ij}`$ and on the “background” conformal factor $`\psi `$. Notice that $`\psi `$ satisfies the equation $`\mathrm{\Delta }\psi =4\pi (\alpha _1\delta _1+\alpha _2\delta _2)`$, where e.g. $`\delta _1\delta (𝐱𝐲_1)`$ denotes the Dirac function at the point $`𝐲_1`$. Obviously we may include the equation for the background conformal factor into our hierarchy of equations by posing $`\mathrm{\Psi }=_{n=0}^+\mathrm{}\beta ^n\sigma _{(n)}`$, with $`\sigma _{(0)}=\psi `$ and $`\mathrm{\Sigma }_{(0)}=4\pi (\alpha _1\delta _1+\alpha _2\delta _2)`$.

### III.2 Analytic closed form of the linearized solution

To the linearized order the perturbation equation reads explicitly

$$\mathrm{\Delta }\sigma _{(1)}=\frac{1}{8}\psi _{ij}s_{ij}+_js_{ij}_i\psi +s_{ij}_{ij}\psi .$$
(17)
Interestingly, this equation turns out to be solvable in analytic closed form. We find that the unique solution of Eq. (17), that is valid in the sense of distributions and tends to zero at spatial infinity (i.e. when $`r|𝐱|+\mathrm{}`$), is

$`\sigma _{(1)}`$ $`=`$ $`{\displaystyle \frac{1}{16}}\mathrm{D}\left({\displaystyle \frac{g}{3}}+{\displaystyle \frac{r_1+r_2}{2r_{12}}}\right)`$ (18)
$`+`$ $`\alpha _1\{{\displaystyle \frac{H_1}{2}}+{\displaystyle \frac{K_1}{32}}{\displaystyle \frac{1}{32}}\mathrm{D}\left({\displaystyle \frac{1}{r_2}}\mathrm{ln}\left[{\displaystyle \frac{r_1}{r_{12}}}\right]\right)+{\displaystyle \frac{9}{32}}\mathrm{D}\left({\displaystyle \frac{\mathrm{ln}r_1}{r_{12}}}\right)+{\displaystyle \frac{\mathrm{D}g}{12r_1}}`$
$`+{\displaystyle \frac{1}{4}}\mathrm{\Delta }_1\left({\displaystyle \frac{g}{r_{12}}}\right){\displaystyle \frac{1}{32}}\mathrm{\Delta }_2\left({\displaystyle \frac{g}{r_{12}}}\right){\displaystyle \frac{1}{24r_1r_{12}^2}}+{\displaystyle \frac{1}{32r_2r_{12}^2}}\}`$
$`+`$ $`\alpha _2\{12\}.`$
Besides the already met shorthand $`\mathrm{D}\frac{^2}{y_1^iy_2^i}`$, we denote the Laplacians with respect to the source points by $`\mathrm{\Delta }_1\frac{^2}{y_1^iy_1^i}`$ and $`\mathrm{\Delta }_2\frac{^2}{y_2^iy_2^i}`$. In Eq. (18), the two terms proportional to $`\alpha _1`$ and $`\alpha _2`$ are deduced from each other by label exchange $`12`$. The solution involves the special functions $`H_1`$ and $`K_1`$ (and $`12`$) which were introduced in Refs. BDI95 ; BFP98  for solving some elementary equations (in the sense of distributions) in the problems of equations of motion and wave generation at the 2PN order. We have,

$`\mathrm{\Delta }H_1`$ $`=`$ $`2{}_{i}{}^{}g_{j}^{}_{ij}\left({\displaystyle \frac{1}{r_1}}\right),`$ (19a)
$`\mathrm{\Delta }K_1`$ $`=`$ $`2\mathrm{D}^2\left({\displaystyle \frac{\mathrm{ln}r_1}{r_2}}\right),`$ (19b)
in which $`{}_{i}{}^{}g_{j}^{}\frac{^2g}{y_1^iy_2^j}`$. The functions $`H_1`$ and $`K_1`$ admit the closed-form expressions given by Eqs. (3.48)-(3.51) in Ref. BDI95 , or equivalently by Eqs. (6.3)-(6.5) in Ref. BFP98 ; for the present purpose we adopt the formulas 999We notice here (since it was not noticed in Refs. BDI95 ; BFP98 ) that $`K_1`$ can also be given an interesting form in connection with a simple elementary kernel $`k_1`$, viz    $`\mathrm{\Delta }k_1`$ $`=`$ $`{\displaystyle \frac{1}{r_1^2r_2}},`$    $`K_1`$ $`=`$ $`\mathrm{D}\left({\displaystyle \frac{1}{r_2}}\mathrm{ln}\left[{\displaystyle \frac{r_1}{r_{12}}}\right]k_1\right).`$

$`H_1`$ $`=`$ $`\mathrm{\Delta }_1\left[{\displaystyle \frac{g}{2r_{12}}}+\mathrm{D}\left({\displaystyle \frac{r_1+r_{12}}{2}}g\right)\right]\mathrm{D}\left({\displaystyle \frac{\mathrm{ln}r_{12}}{r_1}}\right){\displaystyle \frac{3}{2}}\mathrm{D}\left({\displaystyle \frac{\mathrm{ln}r_1}{r_{12}}}\right)`$ (20a)
$``$ $`{\displaystyle \frac{r_2}{2r_1^2r_{12}^2}}+{\displaystyle \frac{1}{2r_1^2r_{12}}}{\displaystyle \frac{1}{2r_1r_{12}^2}},`$
$`K_1`$ $`=`$ $`\mathrm{D}\left({\displaystyle \frac{1}{r_2}}\mathrm{ln}\left[{\displaystyle \frac{r_1}{r_{12}}}\right]\right){\displaystyle \frac{1}{2r_1^2r_2}}+{\displaystyle \frac{1}{2r_2r_{12}^2}}+{\displaystyle \frac{r_2}{2r_1^2r_{12}^2}}.`$ (20b)
Anyway, we find that the exact solution of the linearized perturbation equation is

$`\sigma _{(1)}`$ $`=`$ $`{\displaystyle \frac{1}{16}}\mathrm{D}\left({\displaystyle \frac{g}{3}}+{\displaystyle \frac{r_1+r_2}{2r_{12}}}\right)`$ (21)
$`+`$ $`\alpha _1\{{\displaystyle \frac{1}{2}}\mathrm{\Delta }_1[{\displaystyle \frac{g}{r_{12}}}+\mathrm{D}\left({\displaystyle \frac{r_1+r_{12}}{2}}g\right)]{\displaystyle \frac{1}{2}}\mathrm{D}\left({\displaystyle \frac{\mathrm{ln}r_{12}}{r_1}}\right){\displaystyle \frac{15}{32}}\mathrm{D}\left({\displaystyle \frac{\mathrm{ln}r_1}{r_{12}}}\right)+{\displaystyle \frac{\mathrm{D}g}{12r_1}}`$
$`{\displaystyle \frac{1}{32}}\mathrm{\Delta }_2\left({\displaystyle \frac{g}{r_{12}}}\right){\displaystyle \frac{15r_2}{64r_1^2r_{12}^2}}+{\displaystyle \frac{1}{4r_1^2r_{12}}}{\displaystyle \frac{1}{64r_1^2r_2}}{\displaystyle \frac{7}{24r_1r_{12}^2}}+{\displaystyle \frac{3}{64r_2r_{12}^2}}\}`$
$`+`$ $`\alpha _2\{12\}.`$
Let us quote also, for completeness, the fully explicit form obtained by expanding all the derivatives in this result:

$`\sigma _{(1)}`$ $`=`$ $`{\displaystyle \frac{1}{96r_1r_2}}+{\displaystyle \frac{r_1+r_2}{64r_{12}^3}}+{\displaystyle \frac{1}{192r_{12}}}\left({\displaystyle \frac{1}{r_1}}+{\displaystyle \frac{1}{r_2}}\right){\displaystyle \frac{1}{64r_{12}^3}}\left({\displaystyle \frac{r_1^2}{r_2}}+{\displaystyle \frac{r_2^2}{r_1}}\right)`$ (22)
$`+`$ $`\alpha _1\{{\displaystyle \frac{1}{4r_1^3}}{\displaystyle \frac{13}{64r_{12}^3}}{\displaystyle \frac{1}{24r_1r_{12}^2}}{\displaystyle \frac{5}{192r_1^2r_{12}}}+{\displaystyle \frac{5}{192r_1^2r_2}}{\displaystyle \frac{r_1}{32r_2r_{12}^3}}+{\displaystyle \frac{3}{64r_2r_{12}^2}}`$
$`{\displaystyle \frac{1}{24r_1r_2r_{12}}}+{\displaystyle \frac{r_2}{4r_1r_{12}^3}}{\displaystyle \frac{15r_2}{64r_1^2r_{12}^2}}+{\displaystyle \frac{r_2}{4r_1^3r_{12}}}+{\displaystyle \frac{15r_2^2}{64r_1^2r_{12}^3}}+{\displaystyle \frac{r_2^2}{4r_1^3r_{12}^2}}{\displaystyle \frac{r_2^3}{4r_1^3r_{12}^3}}\}`$
$`+`$ $`\alpha _2\{12\}.`$
The terms in the first line of Eqs. (21) or (22) contribute at the 2PN order in the conformal factor and they are in agreement with Eq. (15). The other terms, proportional to $`\alpha _1`$ or $`\alpha _2`$, will not contribute before the 3PN order \[because $`\alpha _1`$ and $`\alpha _2`$ carry a factor $`1/c^2`$ in front\]. Therefore, only a small part of the linearized approximation, the one given by the first term in (21), is necessary in the proof of Theorem 1. As a matter of fact, the non-linear perturbation we consider contains much more information than a mere post-Newtonian (2PN) expansion. Theorem 1 does not constitute a very stringent requirement on the non-linear solution of the perturbation equations.

A point we make by writing the linear solution $`\sigma _{(1)}`$ into the primary form (18), involving the intermediate functions $`H_1`$ and $`K_1`$, is that the latter functions do enter in the post-Newtonian metric at the 3PN order — with the same numerical coefficients as predicted by Eq. (18). This can be inferred from the expression of the 3PN spatial metric given by Eq. (111) in Ref. livrev , which contains the particular non-linear potential called $`\widehat{X}`$, together with the way that the potential $`\widehat{X}`$ contains the functions $`H_1`$ and $`K_1`$ as calculated by Eq. (6.11) in BFP98 . Thus, although our solution agrees with post-Newtonian calculations up to the 2PN order only, it does contain some correct 3PN features.

### III.3 Quadratic and higher-order approximations

At the level of the second-order perturbation ($`\beta ^2`$) the equation reads

$`\mathrm{\Delta }\sigma _{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{8}}\sigma _{(1)}_{ij}s_{ij}+_js_{ij}_i\sigma _{(1)}+s_{ij}_{ij}\sigma _{(1)}`$ (23)
$`+`$ $`{\displaystyle \frac{1}{8}}\left(s_{ij}\mathrm{\Delta }s_{ij}2s_{ij}_i_ks_{jk}_js_{ij}_ks_{ik}+{\displaystyle \frac{3}{4}}_ks_{ij}_ks_{ij}{\displaystyle \frac{1}{2}}_ks_{ij}_is_{jk}\right)\psi `$
$`+`$ $`\left(s_{ij}_ks_{jk}s_{jk}_ks_{ij}+{\displaystyle \frac{1}{2}}s_{jk}_is_{jk}\right)_i\psi s_{ik}s_{jk}_{ij}\psi .`$
At the next level, cubic-order, the equation will be made of the same terms as in Eq. (23) but with the replacements $`\sigma _{(1)}\sigma _{(2)}`$ and $`\psi \sigma _{(1)}`$, together with many other terms that are purely cubic in $`s_{ij}`$. And so on for the higher-order equations.

We shall see in Section IV that in order to prove the agreement with the 2PN binary’s energy (in contrast with the 2PN metric), we need the full information content about the linearized solution given by Eqs. (21)-(22), and also a crucial piece coming from the second-order perturbation, solution of Eq. (23). In the case of non-linear perturbations ($`n2`$), it is in general impossible to find a solution in analytic closed form. Fortunately, what we shall need is only to control the expansion of $`\sigma _{(2)}`$ at spatial infinity (i.e. in the far zone, $`r+\mathrm{}`$). And this can be achieved, without disposing of the exact expression of $`\sigma _{(2)}`$, from the knowledge of the far-zone expansion of the corresponding source-term $`\mathrm{\Sigma }_{(2)}`$. More generally, the far-zone expansion of the solution $`\sigma _{(n)}`$ for any $`n`$ can be obtained on condition that the far-zone expansion of its source $`\mathrm{\Sigma }_{(n)}`$ has been determined beforehand (say, by induction on $`n`$).

The method is issued from the investigation, in Refs. B98mult ; PB02 , of the multipole expansion of the solution of a Poisson equation with non-compact-support source (i.e. whose support is $`^3`$). In the present context, the multipole expansion is completely equivalent to the far-zone expansion, when $`r+\mathrm{}`$. Following Eq. (C.9) in Ref. PB02 , we obtain the (formal) multipole expansion of $`\sigma _{(n)}`$ in the form 101010Technical notations in Eq. (24) are $`\mathrm{\Delta }^1`$ for the standard Poisson integral; $`Li_1\mathrm{}i_{\mathrm{}}`$ for a multi-index composed of $`\mathrm{}`$ indices; $`_L_{i_1}\mathrm{}_i_{\mathrm{}}`$ for the product of $`\mathrm{}`$ partial derivatives; $`y_Ly_{i_1}\mathrm{}y_i_{\mathrm{}}`$ for the product of $`\mathrm{}`$ spatial vectors. We do not write the $`\mathrm{}`$ summation symbols over the $`\mathrm{}`$ indices composing $`L`$.

$`(\sigma _{(n)})`$ $`=`$ $`FP\left\{\mathrm{\Delta }^1\left[(\mathrm{\Sigma }_{(n)})\right]{\displaystyle \frac{1}{4\pi }}{\displaystyle \underset{\mathrm{}=0}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{()^{\mathrm{}}}{\mathrm{}!}}_L\left({\displaystyle \frac{1}{r}}\right){\displaystyle _^3}d^3𝐲y_L\mathrm{\Sigma }_{(n)}(𝐲)\right\},`$ (24)
where the calligraphic letter $``$ refers to the multipole or equivalently the far-zone expansion. The first term in the right side of (24) represents the effect of integrating “term by term” the multipole expansion of the source $`\mathrm{\Sigma }_{(n)}`$ (this term exists only in the case of non-compact-support sources, because the multipole expansion of a compact-support function is zero). The second term in Eq. (24) is parametrized by the “multipole moments” associated with the solution, which are given by some definite integrals extending over the non-compact support $`^3`$ of $`\mathrm{\Sigma }_{(n)}`$. The symbol FP acts on both terms of (24), and stands for a certain operation of taking the Finite Part, defined in Refs. B98mult ; PB02  by means of a process of complex analytic continuation (with parameter $`B`$). The finite part has proved to play a crucial role, in the case of non-compact-support sources like $`\mathrm{\Sigma }_{(n)}`$, in order to ensure the well-definiteness of the integrals giving the multipole moments in (24). Notice that Eq. (24) has been proved, in Refs. B98mult ; PB02 , to hold in the case of a regular source \[i.e. $`C^{\mathrm{}}(^3)`$\]. Nevertheless, in the presence of the singular points $`𝐲_1`$ and $`𝐲_2`$, the formula (24) can still be applied, but only in the case of that solution $`\sigma _{(n)}`$ for which $`\mathrm{\Delta }\sigma _{(n)}=\mathrm{\Sigma }_{(n)}`$ is satisfied in the sense of distributions.

## IV Total mass and individual black-hole masses

### IV.1 Asymptotic structure of the solution

A central result of the present paper concerns the “asymptotics” of our solution — both at spatial infinity and in the vicinity of the two particles.

###### Theorem 2

(i) The metric $`\gamma _{ij}`$ is asymptotically flat at infinity (when $`r|𝐱|+\mathrm{}`$), i.e.

$$\gamma _{ij}=\delta _{ij}+𝒪\left(\frac{1}{r}\right).$$
(25)
(ii) The two singular points $`𝐲_1`$ and $`𝐲_2`$ are the images, via some appropriate inversions of the radial coordinates: $`\rho _11/r_1`$ and $`\rho _21/r_2`$, of two “asymptotically finite” regions (when $`\rho _1+\mathrm{}`$ and $`\rho _2+\mathrm{}`$), in the sense that the metric, in coordinates $`(\rho _1,𝐧_1)`$, behaves like

$$\mathrm{\Gamma }_{ij}(\rho _1,𝐧_1)=\delta _{ij}+\pi _{ij}(𝐧_1)+𝒪\left(\frac{1}{\rho _1}\right),$$
(26)
where $`\pi _{ij}`$ denotes a certain function of the angles.

Theorem 2 says that the solution is composed of an asymptotically flat universe connected by continuity, via some Einstein-Rosen-like bridges, to two other regions which are asymptotically finite in the sense of Eq. (26). Note that, though the metric (26), unlike the corresponding Brill-Lindquist metric, is not asymptotically flat in the vicinity of the two particles (stricto-sensu), the violation of asymptotic flatness is not very severe, because the metric tends toward a constant with respect to $`\rho _1`$, and does not involve any divergency “at infinity”: it remains asymptotically finite — hence the name. On the other hand, what is very important is that the real universe, as depicted by this solution, is asymptotically flat at spatial infinity in the usual sense of Eq. (25). The global structure described by Theorem 2 represents an attractive feature, we argue, for considering the solution as a physically well-motivated initial state of binary black holes.

The asymptotic flatness when $`r+\mathrm{}`$ follows from the far-zone expansion of the generating function (11b), which is easily checked to start at $`s_{ij}=𝒪\left(\frac{1}{r}\right)`$. Using this fact we find that the source term for the linearized perturbation — i.e. the right side of Eq. (17) — behaves like $`\mathrm{\Sigma }_{(1)}=𝒪\left(\frac{1}{r^3}\right)`$. With the help of Eq. (24) we readily obtain $`\sigma _{(1)}=𝒪\left(\frac{1}{r}\right)`$, as can also be checked directly with our exact results (21)-(22). Then, similarly, we deduce, by induction on the non-linear order $`n`$ \[using Eq. (24) at each step\], that $`\sigma _{(n)}=𝒪\left(\frac{1}{r}\right)`$ for any $`n`$. So $`\sigma =𝒪\left(\frac{1}{r}\right)`$, and the result follows.

The main point about Theorem 2 is the behaviour of the solution in the vicinity of the two particles. When $`r_10`$, we find that $`s_{ij}`$ admits a bounded expansion \[i.e. $`s_{ij}`$ does not blow up when $`r_10`$\], of the type

$$\beta s_{ij}=\epsilon _{ij}(𝐧_1)+𝒪(r_1),$$
(27)
where $`\epsilon _{ij}`$ is a function of $`𝐧_1=\frac{𝐱𝐲_1}{r_1}`$, the direction of approach to the singularity \[$`\epsilon _{ij}`$ is given by Eq. (55a)\]. From Eq. (27), the linearized source-term $`\mathrm{\Sigma }_{(1)}`$ when $`r_10`$ is like some $`\mathrm{\Sigma }_{(1)}=𝒪\left(\frac{1}{r_1^3}\right)`$. Now by integrating “term by term” the expansion of $`\mathrm{\Sigma }_{(1)}`$, we get a similar expansion, but which starts at the order $`𝒪\left(\frac{1}{r_1}\right)`$. Clearly, the solution $`\sigma _{(1)}`$, when $`r_10`$, will be composed of the latter expansion, which represents a particular solution, and augmented by a possible homogeneous solution, solving the (source-free) Laplace equation. However, because our original Poisson equation $`\mathrm{\Delta }\sigma _{(1)}=\mathrm{\Sigma }_{(1)}`$ is satisfied in the sense of distributions, the only possible homogeneous solution we can add is one of the type “regular at the origin”, taking the form of a sum of STF products $`\widehat{x}_L\mathrm{STF}(x_{i_1}\mathrm{}x_i_{\mathrm{}})`$. Because it is regular when $`r_10`$, this homogeneous solution does not modify the leading singular behaviour of $`\sigma _{(1)}`$, that we have therefore proved to be given by $`\sigma _{(1)}=𝒪\left(\frac{1}{r_1}\right)`$. And, again, this can be checked with Eqs. (21)-(22). The argument is easily generalized to any non-linear order $`n`$ (by induction on $`n`$), thus we conclude that the non-linear perturbation $`\sigma `$ diverges when $`r_10`$, but not faster than $`\sigma =𝒪(\frac{1}{r_1})`$. As a result, the expansion of the conformal factor $`\mathrm{\Psi }=\psi +\sigma `$ is of the type

$$\mathrm{\Psi }=\frac{\zeta (𝐧_1)}{r_1}+𝒪\left(r_1^0\right),$$
(28)
where $`\zeta `$ depends on the unit direction $`𝐧_1`$ but not on $`r_1`$. From this we deduce that $`\gamma _{ij}=\mathrm{\Psi }^4\stackrel{~}{\gamma }_{ij}`$ behaves dominantly when $`r_10`$ like

$$\gamma _{ij}(𝐱)=\left(\frac{\zeta (𝐧_1)}{r_1}\right)^4\left[\delta _{ij}+\epsilon _{ij}(𝐧_1)+𝒪\left(r_1\right)\right].$$
(29)
See Eq. (56a) for the expression of $`\zeta (𝐧_1)`$ at linearized order.

Let us now perform some inversion of the radial coordinate $`r_1`$. We consider the coordinate change, valid in a neighbourhood of the particle 1, that is defined by $`𝐱(\rho _1,𝐧_1)`$, where

$$\rho _1=\frac{\zeta ^2(𝐧_1)}{r_1}.$$
(30)
First we have $`dx^i=dr_1^i`$, where $`r_1^i(r_1,𝐧_1)`$, because the particle is at rest: $`v_1^i=0`$. Then, the coordinate transformation $`r_1^i\rho _1^i(\rho _1,𝐧_1)`$ involves simply the change of the radial variable given by Eq. (30), with the angular part $`𝐧_1`$ being unchanged 111111More general coordinate transformations, involving a change of the unit direction $`𝐧_1`$, are also possible.. We compute

$`dr_1^i`$ $`=`$ $`\left({\displaystyle \frac{\zeta }{\rho _1}}\right)^2\left(\delta ^{ij}2n_1^in_1^j+2n_1^i\chi ^j\right)d\rho _1^j,`$ (31a)
$`\chi ^j`$ $``$ $`\left(\delta ^{jk}n_1^jn_1^k\right){\displaystyle \frac{\mathrm{ln}\zeta }{n_1^k}}.`$ (31b)
(Notice that $`n_1^j\chi ^j=0`$.) We then find that, in the new coordinate system $`\rho _1^i=(\rho _1,𝐧_1)`$, the metric, say $`\mathrm{\Gamma }_{ij}`$, admits when $`\rho _1+\mathrm{}`$ the bounded expansion announced in Eq. (26), in which the quantity $`\pi _{ij}(𝐧_1)`$ reads

$$\pi _{ij}4n_1^{(i}\chi ^{j)}+4\chi ^i\chi ^j+(\delta ^{ki}2n_1^kn_1^i+2n_1^k\chi ^i)(\delta ^{lj}2n_1^ln_1^j+2n_1^l\chi ^j)\epsilon _{kl}.$$
(32)
It is such that $`n_1^in_1^j\epsilon _{ij}=n_1^in_1^j\pi _{ij}`$. \[See also Eqs. (57b) and (58).\] This completes the proof of Theorem 2.

### IV.2 The ADM mass

Because the metric is asymptotically flat at spatial infinity (Theorem 2), the binary’s total ADM mass is given by the usual surface integral on a topological 2-sphere at infinity:

$$M=\frac{c^2}{16\pi G}\underset{r+\mathrm{}}{lim}𝑑S^i\left(_j\gamma _{ij}_i\gamma _{jj}\right),$$
(33)
where $`dS^i`$ is the outward surface element on the surface at infinity ($`dS^i=d\mathrm{\Omega }r^2n^i`$ in the case of a coordinate sphere). To be more explicit about $`M`$, we recall from Section IV.1 that some functions $`A_{ij}`$ and $`X`$ of the unit direction $`𝐧=𝐱/r`$ exist so that, when $`r+\mathrm{}`$,

$`\beta s_{ij}`$ $`=`$ $`{\displaystyle \frac{A_{ij}(𝐧)}{r}}+𝒪\left({\displaystyle \frac{1}{r^2}}\right),`$ (34a)
$`\sigma `$ $`=`$ $`{\displaystyle \frac{X(𝐧)}{r}}+𝒪\left({\displaystyle \frac{1}{r^2}}\right).`$ (34b)
Here $`A_{ij}`$ is simply linear in $`\beta `$, while $`X`$ is in the form of a full non-linear series in $`\beta `$: $`X=\beta X_{(1)}+\beta ^2X_{(2)}+\mathrm{}`$ (see also the Appendix). In terms of these functions $`M`$ is given by

$$M=\frac{c^2}{G}\left\{2(\alpha _1+\alpha _2)+\frac{d\mathrm{\Omega }}{4\pi }\left[2X+\frac{1}{2}n^in^jA_{ij}\right]\right\}.$$
(35)
The first term is identical to the result in the Brill-Lindquist solution.

The computation of Eq. (35) is quite straightforward at the level of the linearized approximation, thanks to the closed-form expression (21). However, our aim will be the computation of the 2PN energy in Section V, and we can ascertain beforehand (by counting the required powers of $`G`$), that among all the terms contributing to the 2PN energy there must be one coming from the second-order perturbation, solution of Eq. (23). We find that this particular non-linear term is to be computed only in $`M`$ (not in the individual masses $`m_1,m_2`$). It comes from that part of $`\sigma _{(2)}`$ — or, more precisely, of its leading order coefficient $`X_{(2)}`$ (when $`r+\mathrm{}`$) — which does not involve the constants $`\alpha _1`$ and $`\alpha _2`$, i.e. the part which would be the analogue of the first term in $`\sigma _{(1)}`$ as given by Eq. (21). The other parts, proportional to $`\alpha _1`$ or $`\alpha _2`$, appear at the 3PN order at least.

We succeeded in obtaining this non-linear term thanks to the method described by Eq. (24). First, one can check that the source-term behaves like $`\mathrm{\Sigma }_{(2)}=𝒪\left(\frac{1}{r^4}\right)`$, so there is no direct contribution, computed “term-by-term” from the expansion of the source-term, to the multipole expansion at the order $`𝒪\left(1/r\right)`$: i.e. the first term in the right side of (24) is at least $`𝒪\left(1/r^2\right)`$. The looked-for non-linear contribution is therefore given directly by the value of the “mass monopole” of $`\sigma _{(2)}`$, i.e. the integral appearing in the second term of Eq. (24) for $`\mathrm{}=0`$. So,

$$\sigma _{(2)}=\frac{1}{4\pi r}_^3d^3𝐲\mathrm{\Sigma }_{(2)}(𝐲)+𝒪\left(\frac{1}{r^2}\right).$$
(36)
The integral (36) can be computed analytically 121212This integral is convergent, thus it is unnecessary to include a finite-part operation FP (for notational simplicity we skip the multipole-expansion symbol $``$)., and yields the following contribution to the ADM mass:

$$2X_{(2)}=\frac{1}{2\pi }_^3d^3𝐲\mathrm{\Sigma }_{(2)}(𝐲)=\frac{1}{128r_{12}^3}+𝒪\left(\alpha \right),$$
(37)
where the remainder $`𝒪\left(\alpha \right)`$ indicates that the terms proportional to $`\alpha _1`$ or $`\alpha _2`$ (in this second-order perturbation $`\beta ^2`$), are not to be considered for the present calculation. We thereby obtain $`M`$ with sufficient accuracy for controlling the energy at the 2PN order 131313This means, by the way, that the relative accuracy on $`M`$ itself is actually 3PN — because of the rest-mass contribution.. The result reads

$$M=\frac{c^2}{G}\left\{2(\alpha _1+\alpha _2)+\frac{\beta (\alpha _1+\alpha _2)}{24r_{12}^2}+\frac{\beta ^2}{r_{12}^3}\left[\frac{1}{128}+𝒪\left(\alpha \right)\right]+𝒪\left(\beta ^3\right)\right\}.$$
(38)
The cubic and higher-order perturbations $`𝒪\left(\beta ^3\right)`$ are neglected. See the Appendix for the explicit expansion coefficients needed in this computation.

### IV.3 The black-hole masses

The situation as concerns the two individual masses $`m_1`$ and $`m_2`$ is less easy than with the ADM mass because we dispose only of the notion of “asymptotic finiteness” in the vicinity of the particles, described by the fall-off property (26). Nevertheless, we wish to find an appropriate concept for the black-hole masses. What we shall do is to define $`m_1`$ by the same formula as for the ADM mass $`M`$, but using the coordinate system $`\rho _1^i=(\rho _1,𝐧_1)`$ in the limit where $`\rho _1+\mathrm{}`$. Accordingly we pose

$$m_1=\frac{c^2}{16\pi G}\underset{\rho _1+\mathrm{}}{lim}𝑑S_1^i\left(\frac{\mathrm{\Gamma }_{ij}}{\rho _1^j}\frac{\mathrm{\Gamma }_{jj}}{\rho _1^i}\right)\text{and }12,$$
(39)
where $`\mathrm{\Gamma }_{ij}(\rho _1,𝐧_1)`$ denotes the metric in the coordinates $`\rho _1^i`$, with $`\rho _1`$ being shown in Eq. (30), and where the outward surface element reads $`dS_1^i=d\mathrm{\Omega }_1\rho _1^2n_1^i`$ in the case of a coordinate sphere.

Because $`\mathrm{\Gamma }_{ij}`$ is not asymptotically flat in the usual sense, due to the term with $`\pi _{ij}`$ in Eq. (26), the mass defined by the previous integral will be typically unbounded (i.e. $`m_1`$ will in general tend toward infinity like $`\rho _1`$). Therefore, the definition (39) does not a priori make sense. A possibility would be to discard the divergent part ($`\rho _1`$) of the mass, and thereby to consider only the finite part of the integral in Eq. (39), i.e. the coefficient of the zeroth power of $`\rho _1`$ in the expansion at infinity. This could represent an appropriate postulate for the mass in a general situation \[the finite part prescription would be similar to the FP present in Eq. (24)\]. However such a “finite part” process imposed into the definition of the mass appears to represent (untill more convincing justification is proposed) a somewhat artificial and ad-hoc recipe 141414The ADM mass, given by the surface term that arises in the Hamiltonian, has been obtained even for space-times that are not asymptotically flat HH95 . However, it seems difficult to apply this result in the present case, notably because of our lack of knowledge of the lapse function N in the vicinity of the particles (after radial inversion $`r_1\rho _1`$ of the coordinates).. Gladly enough, we shall not need to invoke any ad-hoc finite part because we shall prove that if we restrict ourselves to the computation of the binary’s energy at the 2PN order, which represents anyway the maximal order at which our solution is physically relevant, the masses $`m_1`$ and $`m_2`$ needed in this computation are perfectly well-defined.

To compute $`m_1`$ and $`m_2`$, we must control the expansion of the metric when $`r_10`$ to one order beyond Eqs. (27)-(28). Let us pose

$`\beta s_{ij}`$ $`=`$ $`\epsilon _{ij}(𝐧_1)+r_1\mu _{ij}(𝐧_1)+𝒪(r_1^2),`$ (40a)
$`\mathrm{\Psi }`$ $`=`$ $`{\displaystyle \frac{\zeta (𝐧_1)}{r_1}}+\eta (𝐧_1)+𝒪\left(r_1\right).`$ (40b)
\[The expressions of the coefficients are relagated to Eqs. (55)-(56) in the Appendix.\] It is straightforward to derive from this the expression of the metric $`\mathrm{\Gamma }_{ij}`$ in the coordinate system $`(\rho _1,𝐧_1)`$, taking into account the term $`1/\rho _1`$ in Eq. (26):

$$\mathrm{\Gamma }_{ij}=\left(1+4\frac{\zeta \eta }{\rho _1}\right)\left[\delta _{ij}+\pi _{ij}\right]+\frac{\zeta ^2}{\rho _1}\kappa _{ij}+𝒪\left(\frac{1}{\rho _1^2}\right),$$
(41)
where $`\kappa _{ij}=(\delta ^{ki}2n_1^kn_1^i+2n_1^k\chi ^i)(\delta ^{lj}2n_1^ln_1^j+2n_1^l\chi ^j)\mu _{kl}`$. Then, by inserting it into Eq. (39), we obtain

$$m_1=\frac{c^2}{G}\underset{\rho _1+\mathrm{}}{lim}\frac{d\mathrm{\Omega }_1}{4\pi }\left[\frac{1}{4}\rho _1\pi _{jj}+2\eta \zeta \left(1+n_1^in_1^j\epsilon _{ij}\right)+\frac{1}{2}\zeta ^2n_1^in_1^j\mu _{ij}\right].$$
(42)
Here, $`\pi _{jj}`$ denotes the trace of the quantity (32), and we have used the facts that $`n_1^in_1^j\pi _{ij}=n_1^in_1^j\epsilon _{ij}`$ and $`n_1^in_1^j\kappa _{ij}=n_1^in_1^j\mu _{ij}`$, together with the useful cancellation of the angular integral

$$\frac{d\mathrm{\Omega }_1}{4\pi }n_1^in_1^j\epsilon _{ij}(𝐧_1)=0,$$
(43)
which can be checked from Eq. (56a) in the Appendix.

We observe that Eq. (42) contains a term proportional to $`\rho _1`$, with $`\pi _{jj}`$ as a factor, which can and does make the mass to become infinite in the limit $`\rho _1+\mathrm{}`$ \[indeed, see Eq. (59)\]. As suggested above, such an infinite term could be removed by some procedure of taking the finite part (whose meaning would a priori be quite unclear), but we find that this infinite term is in fact “negligible” for the present purpose, because $`\pi _{jj}`$ starts only at non-linear order in $`\beta `$,

$$\pi _{jj}=4\chi ^j\chi ^j+4\left(n_1^k\chi ^l+n_1^kn_1^l\chi ^j\chi ^j\right)\epsilon _{kl}=𝒪\left(\beta ^2\right).$$
(44)
\[See also the expression (58).\] Therefore we get, at the linearized level, a finite expression for the mass (formally):

$$m_1=\frac{c^2}{G}\frac{d\mathrm{\Omega }_1}{4\pi }\left[2\eta \zeta \left(1+n_1^in_1^j\epsilon _{ij}\right)+\frac{1}{2}\zeta ^2n_1^in_1^j\mu _{ij}\right]+𝒪\left(\beta ^2\right).$$
(45)
This expression is sufficient to control the 2PN energy. The non-linear corrections $`𝒪\left(\beta ^2\right)`$ are infinite but will not be needed in this paper. In a sense, they do not belong to the “realm” of the present solution, which is limited to the physics at the 2PN approximation: i.e. agreement with the 2PN metric in the near-zone, and internal consistency of the asymptotics up to the 2PN order as regards the energy content of the solution (see Section V) 151515Recall that our proposal for binary black-hole initial data is to adopt the conformal metric (3) and to deduce the conformal factor $`\mathrm{\Psi }`$ from it by numerical methods. The present analytic investigation \[including the theoretical definition of the mass (39)\] is aimed at verifying the physical soundness of this proposal.. As for the linearized terms in Eq. (45), they are dealt with thanks to the explicit result (21), and we get

$$m_1=\frac{2c^2\alpha _1}{G}\left[1+\frac{\alpha _2}{r_{12}}+\frac{\beta }{48r_{12}^2}\left(1\frac{25\alpha _2}{r_{12}}\right)+𝒪\left(\beta ^2\right)\right]\text{and }12.$$
(46)
By setting $`\beta =0`$ into Eq. (46), we reproduce the result valid in the case of the Brill-Lindquist solution.

## V The binary’s geometrostatic energy

With the masses being defined and computed “geometrically” in Section IV, we get the opportunity of an interesting consistency check of our solution, concerning the “geometrostatic” energy that is associated with those masses. Indeed, in Theorem 1 we recovered the 2PN metric in the near zone, therefore we know the coordinate transformation $`\delta x^i=\xi ^i`$ between our presently used coordinate system and the harmonic one. The gauge transformation vector has been obtained in Eq. (10). Therefore it is possible to control, without ambiguity — because we determined the coordinate transformation —, the binary’s interacting energy $`E`$ up to the same 2PN order. If our solution has anything cogent (physically speaking), the latter energy should be in complete agreement with the post-Newtonian result at the 2PN order known from Refs. DD81 ; D82 ; B96 ; BFeom .

###### Theorem 3

The binary’s interacting energy, deduced from those masses $`(\text{33})`$ and $`(\text{39})`$ computed by surface integrals at infinity, i.e.

$$\frac{E}{c^2}=Mm_1m_2,$$
(47)
gives back the energy calculated by standard post-Newtonian methods at the $`2\mathrm{P}\mathrm{N}`$ order (after invoking the same coordinate transformation as in Theorem 1).

The ADM mass has been previously obtained in the form (38), while the masses $`m_1`$ and $`m_2`$ follow from the result (46). The corresponding energy $`E`$ reads then

$$E=\frac{c^4}{G}\left\{\frac{4\alpha _1\alpha _2}{r_{12}}\left(1\frac{25\beta }{48r_{12}^2}\right)+\frac{\beta ^2}{r_{12}^3}\left[\frac{1}{128}+𝒪\left(\alpha \right)\right]+𝒪\left(\beta ^3\right)\right\}.$$
(48)
As we said above, this formula is accurate enough to get full control of the 2PN approximation. We also remarked that at this order the energy is well-defined (i.e. finite).

Our task is to expand that energy when $`c+\mathrm{}`$, using the facts that $`\alpha _1,\alpha _2=𝒪\left(\frac{1}{c^2}\right)`$ and $`\beta =𝒪\left(\frac{1}{c^4}\right)`$ \[actually, we already applied the post-Newtonian approximation when arguing that the remainder $`𝒪\left(\alpha \right)`$ is negligible\]. Also, we know that $`\beta `$ is given by Eq. (11a). We first invert Eq. (46), in an iterative post-Newtonian way, so as to obtain $`\alpha _1,\alpha _2`$ with relative 2PN accuracy. The result is

$$\alpha _1=\frac{Gm_1}{2c^2}\left[1\frac{Gm_2}{2r_{12}c^2}+\frac{G^2m_2(5m_1+3m_2)}{12r_{12}^2c^4}+𝒪\left(\frac{1}{c^6}\right)\right]\text{and }12.$$
(49)
At the 1PN order, we are consistent with the Brill-Lindquist prediction and with what is given by Eq. (14). But, quite naturally, we find that Eq. (49) differs from the Brill-Lindquist result at the 2PN order 161616In the case of the Brill-Lindquist solution we have
$$\alpha _1^{\mathrm{BL}}=\frac{Gm_1}{2c^2}\left[1\frac{Gm_2}{2r_{12}c^2}+\frac{G^2m_2(m_1+m_2)}{4r_{12}^2c^4}+𝒪\left(\frac{1}{c^6}\right)\right].$$
.

Substituting (49) back into (48), we then arrive, after suitable post-Newtonian expansion, at the expression of the 2PN energy in terms of the two “physical” individual masses:

$$E=\frac{Gm_1m_2}{r_{12}}+\frac{G^2m_1m_2(m_1+m_2)}{2r_{12}^2c^2}\frac{G^3m_1m_2(m_1^2+19m_1m_2+m_2^2)}{4r_{12}^3c^4}+𝒪\left(\frac{1}{c^6}\right).$$
(50)
It is composed of the standard Newtonian potential energy (for static black holes), augmented by 1PN and 2PN corrections. The result, however, is not yet the one we want to prove, because the particles’ separation $`r_{12}`$ corresponds to our particular coordinate system. If we want to compare $`E`$ with the post-Newtonian prediction, we must take into account the coordinate transformation (in the “near zone”) that was determined in Theorem 1.

The required link between $`r_{12}`$ and the particles’ separation $`R_{12}`$ in harmonic coordinates is readily obtained with the help of Eq. (10), which permits to compute the shift of the world-lines that is induced by the coordinate transformation 171717We discard an infinite “self-term” when considering the coordinate-transformation vector $`\xi _i(𝐱)`$ at the singular location of the two particles. Thus, from Eq. (10), which is valid a priori only outside the singularities, we obtain the shift vector
$$\xi _i(𝐲_1)=\frac{G^2m_2^2}{4c^4}\frac{n_{12}^i}{r_{12}}+𝒪\left(\frac{1}{c^6}\right).$$
In other words we consider only the finite part of $`\xi _i(𝐱)`$ when $`𝐱𝐲_1`$ (in Hadamard’s sense).. We get

$$r_{12}=R_{12}\frac{G^2(m_1^2+m_2^2)}{4R_{12}c^4}+𝒪\left(\frac{1}{c^6}\right).$$
(51)
The result we find after replacing Eq. (51) into Eq. (50), and effecting the post-Newtonian re-expansion \[the replacement is to be made only into the Newtonian term of (50)\], is

$$E=\frac{Gm_1m_2}{R_{12}}+\frac{G^2m_1m_2(m_1+m_2)}{2R_{12}^2c^2}\frac{G^3m_1m_2(2m_1^2+19m_1m_2+2m_2^2)}{4R_{12}^3c^4}+𝒪\left(\frac{1}{c^6}\right).$$
(52)
Most satisfactorily, we discover that (52) is in complete agreement with the prediction for the 2PN energy in harmonic coordinates, calculated in Refs. DD81 ; D82 ; B96 ; BFeom . The expression we end up with is the same as given by Eq. (B6) in Ref. B96  — in which of course one must set the two particles’ velocities to zero.

The latter agreement is interesting, but we said that it is quite mandatory if our solution makes sense. Technically speaking, it necessitated the control of the metric up to second-order perturbation theory on the Brill-Lindquist background \[cf. the crucial contribution of Eq. (37) to the ADM mass\]. Furthermore, the result (52) can be said to check the global character of the solution (notably the asymptotics therein) — because the physical masses have been computed by surface integrals at infinity. Alternatively, it shows the relevance at 2PN order of the definition (39) we adopted for the black-hole individual masses.

Finally, let us comment that Theorem 3 tells us something about the physical tenets at the basis of the usual post-Newtonian approximation when it is applied to the description of black holes. Indeed, the “post-Newtonian” masses $`m_1`$ and $`m_2`$, which parametrize the post-Newtonian iteration, are introduced as being the coefficients of Dirac delta-functions in the Newtonian density of point-like particles livrev . Now, we have just seen that in fact these post-Newtonian masses agree (within the present approach at least, i.e. in the time-symmetric situation, and up to the 2PN order) with the “geometrostatic” masses which are associated with some Einstein-Rosen-like bridges. We view this as a confirmation that the post-Newtonian calculations, which treat formally the compact objects by means of delta-function singularities 181818For instance, the calculation of the 3PN equations of motion of compact binaries reported in Ref. BFeom ., are appropriate to the description of systems of black holes (as long as the orbital motion of the black holes can be considered to be “slow” in the post-Newtonian sense, i.e. in the so-called inspiralling phase of black-hole binaries).

###### Acknowledgements.

It is a pleasure to thank Gilles Esposito-Farèse and Eric Gourgoulhon for discussions and remarks. This work was partially supported by the International Academic Grant of the Hirosaki University. *

## Appendix A Relevant expansion formulas

In this Appendix we provide the explicit expressions of the various coefficients in the expansion formulas introduced in Section IV for the computation of the masses.

For the ADM mass we need the dominant term, when $`r+\mathrm{}`$, in the generating function $`s_{ij}`$ and the solution of the linear-order perturbation for the conformal factor \[see Eqs. (34)\]. They are given by

$`A_{ij}`$ $`=`$ $`{\displaystyle \frac{\beta }{2r_{12}}}n_{12}^{<i}n_{12}^{j>},`$ (53a)
$`\beta X_{(1)}`$ $`=`$ $`{\displaystyle \frac{\beta }{r_{12}}}\{{\displaystyle \frac{1}{32}}[(𝐧.𝐧_{12})^2{\displaystyle \frac{1}{3}}]+{\displaystyle \frac{1}{48}}{\displaystyle \frac{\alpha _1+\alpha _2}{r_{12}}}\},`$ (53b)
where $`𝐧=𝐱/r`$ and $`𝐧_{12}=(𝐲_1𝐲_2)/r_{12}`$ (with $`𝐧.𝐧_{12}`$ denoting the ordinary scalar product). The second-order perturbation coefficient was already computed in Eq. (37):

$$\beta ^2X_{(2)}=\frac{\beta ^2}{r_{12}^3}\left[\frac{1}{256}+𝒪\left(\alpha \right)\right].$$
(54)
For the masses $`m_1`$ and $`m_2`$ we require the expansions when $`r_10`$ of $`s_{ij}`$ and the conformal factor to one order beyond the dominant term \[see Eqs. (40)\]. For $`s_{ij}`$ we have

$`\epsilon _{ij}`$ $`=`$ $`{\displaystyle \frac{\beta }{r_{12}^2}}\left(n_{12}^{<i}n_{12}^{j>}{\displaystyle \frac{1}{2}}n_1^{<i}n_{12}^{j>}\right),`$ (55a)
$`\mu _{ij}`$ $`=`$ $`{\displaystyle \frac{\beta }{r_{12}^3}}[n_{12}^{}(𝐧_1.𝐧_{12}{\displaystyle \frac{3}{4}})+n_1^{}({\displaystyle \frac{3}{4}}𝐧_1.𝐧_{12}+{\displaystyle \frac{3}{4}}){\displaystyle \frac{1}{4}}n_1^{}].`$ (55b)
For the conformal factor:

$`\zeta `$ $`=`$ $`\alpha _1[1+{\displaystyle \frac{\beta }{r_{12}^2}}({\displaystyle \frac{1}{2}}(𝐧_1.𝐧_{12})^2+{\displaystyle \frac{5}{24}}𝐧_1.𝐧_{12}+{\displaystyle \frac{1}{6}})]+𝒪\left(\beta ^2\right),`$ (56a)
$`\eta `$ $`=`$ $`1+{\displaystyle \frac{\alpha _2}{r_{12}}}+{\displaystyle \frac{\beta }{r_{12}^2}}[{\displaystyle \frac{1}{24}}𝐧_1.𝐧_{12}+{\displaystyle \frac{1}{48}}`$ (56b)
$`+{\displaystyle \frac{\alpha _1}{r_{12}}}({\displaystyle \frac{1}{4}}(𝐧_1.𝐧_{12})^3+{\displaystyle \frac{5}{32}}(𝐧_1.𝐧_{12})^2{\displaystyle \frac{5}{24}}𝐧_1.𝐧_{12}{\displaystyle \frac{5}{96}})`$
$`+{\displaystyle \frac{\alpha _2}{r_{12}}}({\displaystyle \frac{1}{24}}𝐧_1.𝐧_{12}{\displaystyle \frac{25}{48}})]+𝒪(\beta ^2).`$
The quantities $`\chi ^i`$ and $`\pi _{ij}`$ defined by Eqs. (31b) and (32) read (at linear order in $`\beta `$)

$`\chi ^i`$ $`=`$ $`{\displaystyle \frac{\beta }{r_{12}^2}}(𝐧_1.𝐧_{12}{\displaystyle \frac{5}{24}})[n_{12}^i+(𝐧_1.𝐧_{12})n_1^i]+𝒪\left(\beta ^2\right),`$ (57a)
$`\pi _{ij}`$ $`=`$ $`{\displaystyle \frac{\beta }{r_{12}^2}}(n_{12}^{<i}n_{12}^{j>}{\displaystyle \frac{1}{3}}n_1^{<i}n_{12}^{j>}{\displaystyle \frac{1}{6}}(𝐧_1.𝐧_{12})n_1^{<i}n_1^{j>})+𝒪\left(\beta ^2\right).`$ (57b)
Note that at linear order $`\pi _{ij}`$ is trace-free. Its trace arises only at second-order in $`\beta `$,

$$\pi _{jj}=\frac{\beta ^2}{r_{12}^4}(\frac{1}{6}𝐧_1.𝐧_{12}+\frac{5}{144})[(𝐧_1.𝐧_{12})^21]+𝒪\left(\beta ^3\right).$$
(58)
At this order the angular average of $`\pi _{ij}`$ is non-zero:

$$\frac{d\mathrm{\Omega }_1}{4\pi }\pi _{jj}=\frac{5}{216}\frac{\beta ^2}{r_{12}^4}+𝒪\left(\beta ^3\right).$$
(59)


format: text

text:
# Superconductivity near the vibrational mode instability in MgCNi3.

## Abstract

To understand the role of electron-phonon interaction in superconducting MgCNi3 we have performed density functional based linear response calculations of its lattice dynamical properties. A large coupling constant $`\lambda `$= 1.51 is predicted and contributing phonons are identified as displacements of Ni atoms towards octahedral interstitials of the perovskite lattice. Instabilities found for some vibrational modes emphasize the role of anharmonic effects in resolving experimental controversies.

The discovery of superconductivity in MgCNi3 Nature  has generated a new puzzle in the recent series of found superconductors Recent . Despite its relatively low Tc~8K, the presence of Ni signals the possible importance of correlation effects which makes the physics of the pairing mechanism relevant to the famous high Tc cuprates and brings the discussion of unconventional non-electron-phonon mechanism. The experimental information characterizes MgCNi3 as moderate Nature ; Lin  or strong Mao  coupling conventional superconductor through the analysis of specific–heat data, supports the s–wave pairing by CMR experiments CMR , and at the same time shows a zero-bias anomaly in tunnelling data Mao . A clear need for detailed information about the phonon spectra emerges from these controversial data in order to clarify the role of electron–phonon interaction (EPI) and understand the mechanism of superconductivity.

In this work, we perform theoretical studies of the strength of the electron–phonon coupling in MgCNi3 by using fully self–consistent density functional based linear response calculations PHNPRL  of the lattice dynamical properties as a function of phonon wavevector $`𝐪`$. This method was proven to provide reliable estimates for the phonon spectra and electron-phonon interactions in a large variety of systemsPHNPRL ; HTCPRL . We extract a large coupling constant $`\lambda =1.51`$ and identify phonons contributing to it as displacements of Ni atoms towards octahedral interstitials of the perovskite–like cubic structure. We find some of the lattice vibrations to be unstable in linear order which emphasizes the role of strong anharmonic effects in the possible resolution of the recent experimental puzzles.

The basic element of the perovskite MgCNi3 structure is given by a carbon atom placed at the center of the cube and octahedrally coordinated by 6 Ni atoms. Our electronic structure calculation using full potential linear muffin-tin orbital (LMTO) method FPLMTO  reveals Ni-d C-p hybridized valence bands $`ϵ_{𝐤j}`$ in accord with the previous studiesMazin ; Jarlborg ; Shim . The Fermi surface consists of several sheets such as rounded cube sections centered at $`\mathrm{\Gamma }`$, thin jungle gym area spanning from $`R`$ $`[\frac{1}{2}\frac{1}{2}\frac{1}{2}]\frac{2\pi }{a}`$ to $`M`$ $`[\frac{1}{2}\frac{1}{2}0]\frac{2\pi }{a}`$points, dimpled square shaped hole pockets centered around $`X`$ point $`[\frac{1}{2}00]\frac{2\pi }{a}`$ as well as little ovoids along $`\mathrm{\Gamma }R.`$The tight–binding picture discussed before Mazin  consists of carbon $`p_x`$ $`p_y`$ $`p_z`$ orbitals hybridized with $`d`$ states of three Ni atoms numerated accordingly as Nix Niy and Ni$`_z.`$For example (see Fig. 1) carbon $`p_z`$ state hybridizes with Níz $`d_{z^21^2}`$ and also with Níx $`d_{xz}`$y $`d_{yz}`$. Similar picture holds for carbon $`p_x,p_y`$ orbitals which are not shown. As it has been noted Mazin , this in particular results in two antibonding states crossing the Fermi level, which in the nearest neighbor approximation have no dispersion along some directions in the Brillouin zone (BZ).

The appearance of nearly flat areas of $`ϵ_{𝐤j}`$ gives rise to 2D van Hove singularity (vHS) placed 40 meV below the Fermi energy $`E_F`$ which is responsible for a strong narrow peak near $`ϵ_F`$. This has generated speculation about closeness of MgCNi3to ferromagnetic instability upon doping Mazin ; Jarlborg ; Shim ; Pickett . The narrowness of the vHS band is controlled by the second–nearest neighbor hopping integrals, which for the states shown in Fig.1 correspond to the overlap between Nix and Niy $`d_{xy}`$ orbitals. As we discuss in this work, nearly unstable phonon modes exists when each of the two Ni atoms moves toward octahedral interstitial sites. (This is shown in Fig. 1 by arrows for a phonon wave vector $`𝐪=[\frac{1}{2}\frac{1}{2}0]\frac{2\pi }{a}`$). We shall see that these distortions wipe out the narrow vHS peak and give rise to a large electron–phonon coupling.

To calculate lattice dynamics of MgCNi3 as a function of wave vector $`𝐪`$ we utilize the linear response methodPHNPRL ; FPLMTO , We use $`2\kappa `$ LMTO basis set, generalized gradient approximation for exchange–correlationGGA , experimental lattice constant a=7.206 a.u., as well as effective (40,40,40) grid in $`𝐤`$-space (total 1771 irreducible $`𝐤`$ points) to generate the phonon dispersions $`\omega _{𝐪\nu }`$ and electron–phonon matrix elements $`g_{𝐤+𝐪j^{}𝐤j}`$on a (10,10,10) grid of the $`𝐪`$ vectors (total 56 irreducible $`𝐪`$ points).

Our calculated phonon spectrum along major high symmetry lines of the cubic Brillouin zone is given on Fig. 2. The frequencies are seen to be span up to 900 K, with some of the modes showing significant dispersion. In general, we distinguish three panels where the top three branches around 900 K are carbon based, the middle three branches around 600 K are Mg based and 9 lower branches are all Ni based. For the $`\mathrm{\Gamma }`$ point z–polarized modes, in particular, consist of (i) Niz-C against Mg-Nix-Niy vibrations (186 K), (ii) pure Nix-Niy vibrations (262 K), (iii) Mg against Nix-Niy (561 K) (iv) C against Niz vibrations(857 K).

A striking feature of this phonon spectrum is the presence of a low-frequency acoustic mode which is very soft and is even seen to be unstable along $`(\xi \xi 0`$) direction in the BZ. This mode is essentially Ni based and corresponds to perpendicular movements of two Ni atoms towards octahedral interstitials of the perovskite structure. The latter is made of each of the four Ni atoms and two Mg atoms. For example, considering the xy plane (see Fig. 1) for the $`𝐪`$ point M such movements can be seen as a 2D breathing around this vacant interstitial. We find a similar situation for other wave vectors and in other directions of the BZ, where each pair of Ni atoms prefers such in-phase displacements perpendicular to each other. The softness and instability here can be understood as the octahedral interstitials are only places to escape for each Ni atom stressed between two surrounding carbons.

The discussed displacements affect the overlap integrals between nearest Ni t2g orbitals which control the width of the vHS band. For the xy plane these are the hoppings between $`d_{xy}`$ orbitals (see Fig. 1). It is therefore clear that these modes should have large EPI. As they bend the Ni–C bonds, this effect has already been previously emphasizedMazin  by considering rotations of the Ni based octahedra at the R point. Here, we point out movements for wave vectors away from R point, where, for example, Nix–C and Niy–C bonds can be seen as changing their relative angle.

To illustrate the crucial change in the electronic structure due to such distortions, Fig. 3 shows two one–electron spectra in the vicinity of the Fermi level corresponding to the M point frozen phonon involving Níx displacement along the y axis in phase with the Niy displacement along the x axis as illustrated on Fig. 1. The top panel of Fig. 3 corresponds to the undistorted energy bands drawn in the original cubic BZ for easier comparison with the published data Mazin ; Jarlborg ; Shim . As the point $`[\frac{1}{2}\frac{1}{2}0]\frac{2\pi }{a}`$ is now a reciprocal vector of the new doubled lattice, the bands are seen to be simply folded, and the narrow vHS band is readily recognized. When we introduce a distortion by 0.2 Å, the only essential difference is the width of the vHS band which now disperses as much as 1.5 eV. Given the smallness of the assumed displacement this emphasizes the large electron–phonon coupling. Its appearance cannot be understood from a simple geometric overlap between the two $`t_{2g}`$ orbitals ($`d_{xy}`$ states between Nix and Niy shown on Fig. 1). We therefore look for an electronic enhancement due to nesting–like features of the Fermi surface. Indeed, such nesting can be found for the two dimpled square shaped hole pockets centered at X separated exactly by the wave vector $`[\frac{1}{2}\frac{1}{2}0]\frac{2\pi }{a}.`$ We have confirmed that feature by corresponding calculation of the integral $`_{𝐤jj^{}}\delta (ϵ_{𝐤j}ϵ_F)\delta (ϵ_{𝐤+𝐪j}ϵ_F)`$ which provides the total phase space available for the electrons to scatter at given wave vector $`𝐪`$ with no energy change.

We now turn our discussion to the detailed dependence of the electron–phonon coupling across the entire BZ. This at least can be done for all stable phonons. Fig.2 shows the calculated phonon linewidths $`\gamma _{𝐪\nu }`$by widening some representative dispersion curves $`\omega _{𝐪\nu }`$ proportionally to $`\gamma _{𝐪\nu }`$. Each phonon linewidth is proportional to Allen72  $`_{𝐤jj^{}}|g_{𝐤+𝐪j^{}𝐤j}^{𝐪\nu }|^2\delta (ϵ_{𝐤j}ϵ_F)\delta (ϵ_{𝐤+𝐪j}ϵ_F)`$ where the electron–phonon matrix element is found self–consistently from the linear response theory PHNPRL . In particular, we see that some phonons have rather large linewidhts. This, for example holds, for all carbon based higher lying vibrational modes. The strength of the coupling, $`\lambda _{𝐪\nu },`$ for each mode can be obtained by dividing $`\gamma _{𝐪\nu }`$ by $`\pi N(ϵ_F)\omega _{𝐪\nu }^2,`$where $`N(ϵ_F)`$ is the density of states at the Fermi level equal to 5.3 st./\[eV\*cell\] in our calculation. Due to large $`\omega _{𝐪\nu }^2,`$ this unfortunately results in strongly suppressed coupling for all carbon modes which would favor high critical temperatures. The coupling, however, is relatively strong for the Ni based modes. For example, we can find $`\lambda `$’s of the order of 1–3 for the Ni based optical phonons around 250 K. Again, the analysis of the polarization vectors shows that these vibrations involve Ni movements towards octahedral interstitials. For example, (see Fig. 1) the movement of Nix and Niy atoms along the z direction either in-phase or out-of phase result in larger overlap between Niz Niy $`d_{yz}`$ orbitals and between Niz Nix $`d_{xz}`$ orbitals. Similar to what we find for the M point using the frozen phonon method (Fig.3), this again enlarges the vHS bandwidth resulting in large EPI. A extremely large coupling ($`\lambda \mathrm{~}70)`$ occurs for the soft acoustic mode involving the interstitial breathing (M–point). Here we point out a triple effect: (i) breathing of four Nix Niy atoms into the interstitial, which results in larger $`d_{xy}`$ overlap, (ii) nesting enhancement which helps wiping out the vHS peak, and, (iii) the smallness of $`\omega _{𝐪\nu }^2.`$As the linewidth of this particular phonon is so large, the concept of phonon itself has to be questioned, but due to a smallness of the phase space associated with this vibration, this has a little effect on integral characteristics such as $`\lambda .`$

Unfortunately, finding the integral value of $`\lambda `$ is another challenging problem due to the appearance of the imaginary frequencies. Neglecting the unstable mode completely results in the average coupling constant equal to 0.95 mainly due to the discussed Ni vibrations around 250 K (see Fig. 2).. It is however clear that the low–frequency mode has a large contribution to $`\lambda `$ and cannot be omitted. We have asked a question for which wave vectors $`𝐪`$ this mode is unstable? Since we know its dispersion across the entire BZ, we can determine a surface in $`𝐪`$ space which separates the real and imaginary frequencies. Fig. 4 shows the result of such an analysis. We see the area around the $`\mathrm{\Gamma }`$ point which continues along the lines towards the $`R,M,`$ and $`X`$ points. Here we find the stability of the mode. The area around the M point, where we find enormously large EPI, is seen to be very small . This is also clear from the dispersion relations shown on Fig. 2 for the ($`\xi \xi 0`$) direction where the frequency becomes real just near the M point itself.

This instability carefully avoids symmetry points which can be understood by keeping in mind the symmetry of the discussed distortions. Namely, anharmonicity is expected to be large for all in–plane movements involving two, three, or four Ni atoms towards octahedral interstitials. This is forbidden near the $`\mathrm{\Gamma }`$ point and for all wavevectors $`(\xi 00)`$ along $`\mathrm{\Gamma }X`$ as well as for $`(\xi \xi \xi )`$ along $`\mathrm{\Gamma }R`$ directions if, for example, xy plane is considered. It is however allowed for $`(\xi \xi 0)`$ along $`\mathrm{\Gamma }M`$and for $`(\frac{1}{2}\xi 0)`$ along $`XM.`$ Moreover, for the wavevectors with $`q_z0,`$ the displacements away from xy plane are allowed and the instability is quickly suppressed. This, e.g., is the case of $`(\frac{1}{2}\frac{1}{2}\xi )`$ along $`MR.`$

The persistence of the instability which does not occur for any of the high–symmetry point needs a non–trivial frozen–phonon analysis. As our polarization vectors promt that the largest anharmonicity is expected for the $`\mathrm{\Gamma }XM`$ plane with $`q_z=0`$, we have performed three such calculations for the points $`a=`$($`\frac{1}{4}\frac{1}{4}0)\frac{2\pi }{a}`$, $`b=`$($`\frac{1}{2}\frac{1}{4}0)\frac{2\pi }{a}`$, and $`M`$ shown on Fig. 4. This corresponds to the discussed displacements of two (point $`a`$), three (point $`b`$), or four (point $`M`$) Ni atoms within xy–plane. The results of these calculations reveal esentially anharmonic interatomic Ni potentials. A shallow double well with a depth of the order of 20–40 K and the curvature at the equillbrium of the order of 40–50i K exists for the $`a`$ and $`b`$ points which becomes vanishing at the M point. We have found such a behavior by both total energy and force calculations for the supercells of 20 atoms ($`a`$ and $`b`$ points) and 10 atoms ($`M`$ point). Clearly, such a small depth on the temperature scale indicates that the distortions are dynamical and zero point motions would prevent the appearance of the static long–range order. This picture is confirmed by our recent extended x–ray absorbtion fine spectroscopy experiments (EXAFS) which reveal dynamical distortions Ignatov  not seen in the diffraction experiments. Remarkably that the parameters of the double well (depth ~20K and the equillbrium curvature ~50i K) extracted from EXAFS are in accord with our frozen phonon data.

It is now clear that our instabile mode is not related to a statically distorted structure of MgCNi3 but should be resolved using an anharmonic theory of the EPI. Using a combination of our frozen phonon and linear response data we are able to map the behavior of the instable mode inside the entire Brillouin zone using a (4,4,4) grid of $`𝐪`$–vectors. This gives us a possibility to determine a contribution to the Eliashberg spectral function using a generalized formula Allen ; BKBO

$`\alpha ^2F(\omega )={\displaystyle \frac{1}{N(ϵ_F)}}{\displaystyle \underset{𝐪}{}}{\displaystyle \underset{𝐤j𝐤^{}j^{}}{}}{\displaystyle \underset{n}{}}\delta [\omega \omega _n(𝐪)]\times `$
$`{\displaystyle \frac{[f_{𝐤j}f_{𝐤^{}j^{}}]\delta [ϵ_{𝐤j}ϵ_{𝐤^{}j}+\omega _n(𝐪)]|G_{𝐤^{}j^{}𝐤j}^{(n)}(𝐪)|^2}{[\omega _n(𝐪)]^2}}`$ (1)
where the summation over the eigenstates $`n`$ of the anharmonic oscillator appears here with $`\omega _n(𝐪)=ϵ_n(𝐪)ϵ_0(𝐪)`$ being the excitation frequencies around the ground state level for the anharmonic well at a given wave vector $`𝐪`$. The generalized electron–phonon matrix element $`G_{𝐤^{}j^{}𝐤j}^{(n)}(𝐪)`$ involves transitions from the ground to $`n`$th excited state and also includes changes in the effective one–electron potential to all orders with respect to the displacements. The most important is the effect of the frequency renormalization which appears in denominator of Eq. (1). As we find, our anharmonic mode has the first excitation frequency which, depending on the $`𝐪`$–point, varies from 120 to 150 K. We also find that the summation over $`n`$ is fastly convergent in accord with the previous workAllen : the oscillator strength for $`n`$=2 is always less than 2% of that for $`n`$=1. The linear electron–phonon scattering matrix elements are known to us from the linear response calculations, and numerically small corrections due to second–order couplings have been extracted from our frozen phonon data using the band splittings technique BKBO .

Our resulting value of $`\lambda `$ for this anharmonic mode appears to be 0.56. Adding the result for $`\lambda =0.95`$ from all other modes our total $`\lambda `$ is now 1.51. This is consistent with the values of $`\lambda `$ extracted from specific heat measurements Nature ; Lin ; Walte  and would cause MgCNi3 to be a strongly coupled electron–phonon superconductor. To estimate the T$`_c,`$ we solve the Eliashberg equation with our total (harmonic+anharmonic) $`\alpha ^2F(\omega ).`$ The cutoff frequency $`\omega _{cut}`$ is taken to be 10 times the maximum phonon frequency. The Coulomb pseudopotential $`\mu ^{}(\omega _{cut})`$ is varied and the Tc within the range 7–20 K is obtained. To make connections with the values of $`\mu ^{}`$ which are standardly discussed with the Allen–Dynes modified McMillan Tc expression McMillan  we need to use a rescaling formula V  $`[\mu ^{}]^1=[\mu ^{}(\omega _{cut})]^1+\mathrm{ln}(\omega _{cut}/\omega _{\mathrm{log}})`$ with $`\omega _{\mathrm{log}}120K.`$ We found that that both the experimental $`T_c`$ and the superconducting energy gap $`\mathrm{\Delta }(0)`$ can be reproduced with the value of $`\mu ^{}=0.33`$. While the effect of spin fluctuations on superconductinvity and mass enchancement is better to discuss in terms of its own coupling constant $`\lambda _{spin}`$, at the absence of detailed calculation of the latter, we can only conclude that our enchanced $`\mu ^{}`$ can in principle appear due to localized nature of Ni $`d`$ orbitals. It is, for example, known that spin fluctuatlions are important in V and $`\mu ^{}\mathrm{~}0.3`$ is neededV . In fact, at the absence of spin fluctuations, $`\mu ^{}`$ is usually 0.1–0.15 and Tc in MgCNi3would be larger. We conclude that spin fluctuations partially suppress superconductivity, the result expected from the conventional theory. Note that the same conclusion has been reached based on the recent analysis of the specific heat data Walte . We finally calculated the phonon contribution to the specific heat and fitted it to the form $`C(T)=\beta T^3`$ at low temperatures. Our value of $`\beta =0.35`$ mJ/\[mol\*K4\] is close to the values 0.40–0.42 mJ/\[mol\*K4\] determined experimentally Lin ; Mao ; Walte .

In conclusion, by performing linear response calculations of the electron–phonon interaction in MgCNi3 we reported the value of $`\lambda =1.51`$ consistent with the strong–coupling limit of electron–phonon mechanism of superconductivity. The unusually large anharmonic correction to $`\lambda `$ for the lattice near instability is emphasized.

We thank O. V. Dolgov, I. I. Mazin, W. Pickett for useful comments, and S. Shulga for availability of his programs to solve Eliashberg’s equations. The work was supported by the grants NSF DMR No.9733862, 0209243, 0238188, 0342290, US DOE No. DE–FG02–99ER45761, NJSGC No. 02–42.