text:
# On the accuracy of the post-Newtonian approximation11footnote 1In “2001: a relativistic spacetime odyssey”, Proc. of the 25th Johns Hopkins workshop, edited by Ignacio Ciufolini, Daniele Dominici and Luca Lusanna, World Scientific, p. 411 (2001).

## 1 Introduction

The “standard” post-Newtonian approximation, or expansion when the speed of light $`c+\mathrm{}`$, is at the basis of an important body of research, which has provided us in the past with our best picture of the motion of compact objects in general relativity. We quote the pionneering works of Einstein, Droste, and DeSitter, the landmark analysis due to Einstein, Infeld and Hoffmann of the dynamics of $`N`$ separated bodies at the first post-Newtonian order (1PN, or $`1/c^2`$), and the seminal papers by Chandrasekhar and collaborators concerning the equations of motion of extended fluid systems, up to the 2.5PN level (at which order appears the first effect due to the reaction to the emission of gravitational radiation). In the case of two compact objects (neutron stars or black holes), we possess the 2.5PN equations of motion of the binary pulsar, and the 3PN equations of motion of inspiralling compact binaries. Regarding the problem of the gravitational radiation emitted by inspiralling compact binaries, we have under control most of the gravitational-wave flux up to the 3.5PN order, including the specific effects of wave tails.

In this paper we focus our attention on the question of the dynamics of black-hole binary systems (henceforth we assume that the compact objects are black holes) and its recent resolution up to the 3PN approximation, corresponding to the formal order $`1/c^6`$ beyond the Newtonian force law. On the one hand the ADM-Hamiltonian formalism of general relativity has been applied at the 3PN order by Jaranowski and Schäfer, and Damour, Jaranowski and Schäfer. On the other hand a direct 3PN iteration of the equations of motion — instead of a Hamiltonian — in harmonic coordinates (extending the method proposed in Ref.) has been implemented by Blanchet and Faye, and de Andrade, Blanchet and Faye. These two independent approaches have succeeded; it has been shown that there exists a unique transformation of the particle’s dynamical variables that changes the 3PN harmonic-coordinates Lagrangian into a Lagrangian whose Legendre transform is exactly identical to the 3PN ADM-coordinates Hamiltonian.

In the previous approaches the two black holes are modelled by point particles, solely described by their masses $`m_1`$ and $`m_2`$. This is consistent with the very spirit of the post-Newtonian method. The point-particle description has to be supplemented by a process of regularization of the self field of point particles. The standard regularization à la Hadamard is at the basis of the works (an extended version of this regularization has been defined in Refs.). Unfortunately it was shown that at the 3PN order the equations of motion of black holes contain a particular coefficient, i.e. the “static” ambiguity $`\omega _\mathrm{s}`$ in the ADM-Hamiltonian formalism or the parameter $`\lambda `$ in the harmonic-coordinates approac222See Eq. (8) below for the relation linking together $`\lambda `$ and $`\omega _\mathrm{s}`$., which cannot be fixed by the Hadamard regularization. The value of this coefficient has been obtained by means of a dimensional regularization instead of the Hadamard one within the ADM-Hamiltonian formalism. We shall discuss some implications of this result (i.e. $`\omega _\mathrm{s}=0`$) regarding the validity of the post-Newtonian expansion.

Let us look at the conserved energy of the black-hole binary in the center-of-mass frame at the 3PN order. Technically this energy follows from the 3PN harmonic-coordinates Lagrangian or equivalently from the 3PN ADM-coordinates Hamiltonian. The energy is conserved only when we neglect the radiation reaction damping at the 2.5PN order. The time derivative of the energy is equal to the radiation reaction effect, but for the present discussion we are interested only in the conservative part of the dynamics which is composed of the Newtonian, 1PN, 2PN and 3PN approximations. Specializing to the case of orbits which are circular (apart from the gradual radiation-reaction inspiral)333We know that most inspiralling compact binaries have circular orbits because of radiation-reaction effects. yields then the center-of-mass energy $`E`$ at the 3PN order for circular orbits as a function of the frequency $`\omega `$ of the orbital motion.

Here we want to assess the validity of the post-Newtonian approximation. More precisely we address, and to some extent we answer, the following questions. How accurate is the post-Newtonian expansion for describing the dynamics of binary black hole systems ? Is the innermost circular orbit (ICO) of binary black holes, defined by the minimum of the energy function $`E(\omega )`$, accurately determined at the highest currently known post-Newtonian order ? This question is pertinent because the ICO represents a point in the late stage of evolution of the binary which is very relativistic (orbital velocities of the order of 50% of the speed of light). How well does the 3PN approximation as compared with the prediction provided by numerical relativity ? What is the validity of the various post-Newtonian resummation techniques which aim at “boosting” the convergence of the standard post-Newtonian approximation ?

The previous questions are very interesting but difficult to settle down rigorously. Indeed the very essence of an approximation method is to cope with our ignorance of the higher-order terms in some expansion, but the higher-order terms are precisely the ones which would be needed for a satisfying answer to these problems. So we shall be able to give only some educated guesses and/or plausible answers, that we cannot justify rigorously, but which seem very likely from the standard point of view on the post-Newtonian theory, in particular that the successive orders of approximation get smaller and smaller as they should (in average), with only few accidents occuring at high orders where a particular approximation would be abnormally large with respect to the lower-order ones. Admittedly, in addition, our faith in the estimation we shall give regarding the accuracy of the 3PN order for instance, comes from the historical perspective, thanks to the many successes achieved in the past by the post-Newtonian approximation when confronting the theory and observations (e.g. of the Solar system dynamics and the binary pulsar). It is indeed beyond question, from our past experience, that the post-Newtonian method does work.

There are many other related issues that we shall not address. For instance : is the best presently known post-Newtonian approximation, i.e. 3PN, sufficient for the problem of the coalescence of two black holes in a foreseeable future444We have in mind the 3PN approximation considered as an input for the numerical calculation of the black-hole coalescence (assuming that a clear method would exist for implementing the post-Newtonian initial conditions into a numerical scheme). ? A related problem is the qualitative and quantitative influence of the high-order radiation reaction effects which have to be superposed to the conservative part of the dynamics. Finally we shall discuss only the binary’s equations of motion, and leave aside the problem of the radiation field. Certainly the accuracy of the 3.5PN gravitational-radiation flux of black-hole binaries which has been computed in Refs. should be discussed in a similar way.

Basically, the point we would like to emphasize555We are following the recent study in Ref.. is that the post-Newtonian approximation, in standard form (without using the resummation techniques advocated in Refs.), is able to located the ICO of two black holes with an excellent accuracy of the order of 1%, and perhaps much better than that. At first sight this statement is rather surprising, because the dynamics of two black holes at the point of the ICO is so relativistic. Indeed one often ears about the “bad convergence”, or the “fundamental breakdown”, of the post-Newtonian series in the regime of the ICO. However our estimates do show that the 3PN approximation is very good in this regime, and they are also confirmed by the remarkable agreement, we shall detail below, with a numerical calculation by Gourgoulhon et al. of the ICO in the case of two black holes moving on exactly circular orbits (“helical symmetry”). When comparing the post-Newtonian approximation with the numerical simulation we face an interesting problem : since the numerical work has been done for corotating black holes, which spin with the orbital velocity $`\omega `$, the effects of spins, appropriate to two Kerr black holes, are to be taken into account within our post-Newtonian framework.

Another point we shall develop is that most probably the general-relativistic dynamics of two objects with comparable masses is qualitatively different (in addition of being far more complicated), than the dynamics of a “test” particle with geodesic motion on a fixed background Schwarzschild space-time. Our argument has something to do with the value $`\omega _\mathrm{s}=0`$ obtained in Ref. for the 3PN regularization ambiguity parameter. Indeed this value, if we believe that it is really the prediction of general relativity, suggests that the two-body interaction is not “Schwarzschild-like”, in the sense that it does not seem to admit a light-ring singularity similar to the one of the Schwarzschild space-time. We then argue that this fact sheds a doubt on the validity, at high post-Newtonian orders, of some resummation techniques (like Padé approximants) that are based on the assumption that the field generated by two bodies of comparable sizes is a “deformation” of the Schwarzschild metric. A contrario we shall find that the value $`\omega _\mathrm{s}=0`$ gives some convincing evidence that the standard post-Newtonian approximation, based on Taylor rather than Padé approximants, is very accurate.

## 2 The binding energy of two black-holes

The binding energy, in the center-of-mass frame, is defined as 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). Restricting our consideration to circular orbits, the energy is a function of a single variable, which can be chosen to be the distance $`r`$ between the two particles in a given coordinate system, or, better, the frequency $`\omega =\frac{2\pi }{P}`$ of the orbital motion ($`P`$ is the orbital period). We introduce for convenience the particular frequency-related parameter

$$x\left(\frac{GM\omega }{c^3}\right)^{2/3}.$$
(1)
The mass-energies of the black holes are $`m_1`$ and $`m_2`$ (they take notably into account the rotational energies); the total mass is $`M=m_1+m_2`$. It is important to express the energy in terms of the frequency-related parameter $`x`$, instead of some coordinate distance $`r`$, because the function $`E(x)`$ then takes an invariant expression (the same in different coordinate systems).

Having in hands the circular-orbit energy, we then define the innermost circular orbit (ICO) as the minimum, when it exists, of the energy function $`E(x)`$. The definition is motivated by our comparison with the numerical calculation, because this is precisely that minimum which is computed numerically. In particular, we do not define the ICO as a point of dynamical general-relativistic unstability. The circular-orbit energy, developed to the 3PN order, is of the form

$$E(x)=Mc^2\frac{\mu c^2x}{2}\left\{1+a_1(\nu )x+a_2(\nu )x^2+a_3(\nu )x^3+𝒪(x^4)\right\}.$$
(2)
The first term is the rest-mass, the next one, proportional to $`x`$, is the Newtonian term, and then we have many post-Newtonian corrections, the coefficients of which are known at present only up to the 3PN order, and given by

$`a_1(\nu )`$ $`=`$ $`{\displaystyle \frac{3}{4}}{\displaystyle \frac{\nu }{12}},`$ (3)
$`a_2(\nu )`$ $`=`$ $`{\displaystyle \frac{27}{8}}+{\displaystyle \frac{19}{8}}\nu {\displaystyle \frac{\nu ^2}{24}},`$ (4)
$`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.`$ (5)
We make use of the useful ratio between the reduced and total masses :

$$\nu =\frac{\mu }{M}\text{where}\mu =\frac{m_1m_2}{M}.$$
(6)
This ratio is interesting because of its range of variation,

$$0<\nu \frac{1}{4},$$
(7)
where $`\nu =\frac{1}{4}`$ in the equal-mass case and $`\nu 0`$ in the test-mass limit for one of the bodies. We shall investigate the value of the ICO as predicted by the 3PN energy (2)-(5) in Section 4.

The 3PN coefficient $`a_3(\nu )`$ involves the regularization-ambiguity parameter $`\lambda `$ introduced in Refs. and due to a physical incompleteness in the Hadamard method for regularizing the self-field of point particles. This parameter is equivalent to the parameter $`\omega _\mathrm{s}`$ in Refs. and related to it by :

$$\lambda =\frac{3}{11}\omega _\mathrm{s}\frac{1987}{3080}.$$
(8)
It has been argued in Ref. that the numerical value of $`\omega _\mathrm{s}`$ could be $`9`$, because for such a value some different resummation techniques, when they are implemented at the 3PN order, give approximately the same result for the ICO. Even more, it was suggested that $`\omega _\mathrm{s}`$ might be precisely equal to $`\omega _\mathrm{s}^{}`$, with

$$\omega _\mathrm{s}^{}=\frac{47}{3}+\frac{41}{64}\pi ^2=9.34\mathrm{}.$$
(9)
However, a more recent computation performed with the help of a dimensional regularization instead of the Hadamard regularization, within the ADM-Hamiltonian formalism, has yielded

$$\omega _\mathrm{s}=0\lambda =\frac{1987}{3080}.$$
(10)
Here we adopt $`\omega _\mathrm{s}=0`$ as the “correct” value predicted by general relativity for the ambiguity parameter. Note that the result (10) is quite different from $`\omega _\mathrm{s}^{}`$ given by Eq. (9) : this suggests, according to Ref., that different resummation techniques, viz Padé approximants and effective-one-body methods, which 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).

The appearance of one ambiguity parameter at the 3PN order is interesting in connection with the so-called effacing principle satisfied by general relativity, according to which the equations of motion and the radiation field of gravitationally condensed objects should depend only on their relativistic masses and not on the detailed features of their internal structure — we neglect the tidal effects between the objects. In general relativity this principle follows from the strong equivalence principle (which differs from the Einstein equivalence principle by the inclusion of bodies with self-gravitational interactions and of experiments involving gravitational forces). Indeed, in the freely falling frame of one of the (spherically symmetric) compact objects, we can apply the Birkhoff theorem and find that the external vacuum field depends only on the mass. In consequence we should expect that the parameter $`\lambda `$ is a pure number (e.g. a rational fraction), independent of the internal structure of the compact objects.

It would be interesting to confirm Eq. (10) by an independent calculation, hopefully without resort to any regularization scheme. We have in mind a calculation valid for extended (“fluid”) systems, taking a priori into account the internal structure of the fluids. Such a calculation would provide, in addition to the determination of $`\lambda `$, an explicit verification of the effacing principle of general relativity. Notice that from the point of view of a calculation valid for extended fluids, it is hard to believe that the ambiguity parameter could depend on $`\pi ^2`$ like in Eq. (9).

In the limiting case $`\nu 0`$, the energy (2)-(5) reduces to the 3PN approximation of the exact expression for a test particle in the Schwarzschild background,

$$E^{\mathrm{Sch}}(x)=\mu c^2\frac{12x}{\sqrt{13x}}.$$
(11)
The minimum of that function or ICO occurs at $`x_{\mathrm{ICO}}^{\mathrm{Sch}}=\frac{1}{6}`$, and we have $`E_{\mathrm{ICO}}^{\mathrm{Sch}}=\mu c^2\left(\sqrt{\frac{8}{9}}1\right)`$. We know that the ICO in the case of the Schwarzschild metric is also an innermost stable circular orbit (or ISCO), i.e. it corresponds to a point of dynamical unstability. Another important feature of Eq. (11) is the singularity at the value $`x_{\mathrm{light}\mathrm{ring}}^{\mathrm{Sch}}=\frac{1}{3}`$ which corresponds to the famous circular orbit of photons in the Schwarzschild metric (“light-ring” singularity). This orbit can also be viewed as the last unstable circular orbit. We can check that the post-Newtonian coefficients $`a_n^{\mathrm{Sch}}a_n(0)`$ corresponding to Eq. (11) are given by

$$a_n^{\mathrm{Sch}}=\frac{3^n(2n1)!!(2n1)}{2^n(n+1)!}.$$
(12)
They increase with $`n`$ by roughly a factor 3 at each order. This is simply the consequence of the fact that the radius of convergence of the post-Newtonian series is given by the Schwarzschild light-ring singularity at the value $`\frac{1}{3}`$. We may therefore recover the light-ring orbit by investigating the limit

$$\underset{n+\mathrm{}}{lim}\frac{a_{n1}^{\mathrm{Sch}}}{a_n^{\mathrm{Sch}}}=\frac{1}{3}=x_{\mathrm{light}\mathrm{ring}}^{\mathrm{Sch}}.$$
(13)
## 3 Accuracy of the post-Newtonian expansion

Let us now discuss the accuracy of the 3PN approximation when estimating the ICO of binary black holes. First of all we make a few order-of-magnitude estimates. At the location of the ICO we shall find (see the next Section) that the frequency-related parameter $`x`$ defined by Eq. (1) 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 that order. 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 which 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}x^{1/2}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 3PN will be even better, with the accuracy $`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 can reasonably expect that the neglected higher-order approximations, 4PN and so on, are in fact 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 estimate makes sense only if the numerical values of the post-Newtonian coefficients in Eqs. (3)-(5) stay 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. It has often been emphasized in the litterature (see e.g. Refs.) that in the test-mass limit $`\nu 0`$ the post-Newtonian series converges slowly, so the post-Newtonian approximation is not very good in the regime of the ICO. Indeed we have seen that when $`\nu =0`$ the radius of convergence of the series is $`\frac{1}{3}`$ (not so far from $`x_{\mathrm{ICO}}^{\mathrm{Sch}}=\frac{1}{6}`$), and that accordingly the post-Newtonian coefficients increase by a factor $`3`$ at each order. So it is perfectly correct that in the case of test particles in the Schwarzschild background the post-Newtonian approximation is to be carried out to a high order in order to locate the turning point of the ICO.

Let us immediately remark that this negative conclusion does not matter : indeed we shall never use the post-Newtonian approximation when $`\nu 0`$ simply because we know the exact result which is given by Eq. (8)666The exact result for the radiation field is also known, albeit only numerically.. Therefore we should not worry 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.

What happens when the two masses are comparable ($`\nu =\frac{1}{4}`$) ? It is clear that the accuracy of the post-Newtonian approximation depends crucially on how rapidly the post-Newtonian coefficients increase with $`n`$. We have seen that in the case of the Schwarzschild metric the latter increase is in turn related to the existence of a light-ring orbit. For continuing the discussion we shall say that the relativistic interaction between two bodies of comparable masses is “Schwarzschild-like” if the post-Newtonian coefficients $`a_n(\frac{1}{4})`$ increase when $`n+\mathrm{}`$. If this is the case this signals the existence of something like a light-ring singularity which could be interpreted as the deformation, when the mass ratio $`\nu `$ is “turned on”, of the Schwarzschild light-ring orbit. By analogy with Eq. (13) we can estimate the location of this “pseudo-light-ring” orbit by

$$\frac{a_{n1}(\nu )}{a_n(\nu )}x_{\mathrm{light}\mathrm{ring}}(\nu )\text{with }n=3.$$
(14)
Here $`n=3`$ is the highest known post-Newtonian order. If the two-body problem is “Schwarzschild-like” then the right-hand-side of Eq. (14) is small (say around $`\frac{1}{3}`$), the post-Newtonian coefficients typically increase with $`n`$, and most likely it should be difficult to get a reliable estimate by post-Newtonian methods of the location of the ICO. So we ask : is the gravitational interaction between two comparable masses Schwarzschild-like ?

In Table 1 we present the values of the coefficients $`a_n(\nu )`$ in the test-mass limit $`\nu =0`$ \[see Eq. (12) for the analytic expression\], and in the equal-mass case $`\nu =\frac{1}{4}`$ when the ambiguity parameter takes the “uncorrect” value $`\omega _\mathrm{s}^{}`$ \[defined by Eq. (9)\] and the correct one $`\omega _\mathrm{s}=0`$ predicted by general relativity. When $`\nu =0`$ we clearly see the expected increase of the coefficients by roughly a factor 3 at each step. Now when $`\nu =\frac{1}{4}`$ and $`\omega _\mathrm{s}=\omega _\mathrm{s}^{}`$ we notice that the coefficients increase approximately in the same manner as in the test-mass case $`\nu =0`$. This indicates that the gravitational interaction in the case of $`\omega _\mathrm{s}^{}`$ looks like that in a one-body problem. The associated pseudo-light-ring singularity is estimated using Eq. (14) as

$$x_{\mathrm{light}\mathrm{ring}}(\frac{1}{4},\omega _\mathrm{s}^{})0.32.$$
(15)
The pseudo-light-ring orbit seems to be a very small deformation of the Schwarzschild light-ring orbit given by Eq. (13). In this Schwarzschild-like situation, we should not expect the post-Newtonian series to be very accurate.

Now in the case $`\nu =\frac{1}{4}`$ but when the ambiguity parameter takes the correct value $`\omega _\mathrm{s}=0`$, we see that the 3PN coefficient $`a_3(\frac{1}{4})`$ is of the order of minus one instead of being $`10`$. This suggests (unless 3PN happens to be quite accidental) that the post-Newtonian coefficients in general relativity do not increase very much with $`n`$. This is an interesting finding because it indicates that the actual two-body interaction in general relativity is not Schwarzschild-like. There does not seem to exist something like a light-ring orbit which would be a deformation of the Schwarzschild one. Applying Eq. (14) we obtain as an estimate of the “light-ring” :

$$x_{\mathrm{light}\mathrm{ring}}(\frac{1}{4},\text{G.R.})2.88.$$
(16)
It is clear that 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 problem. Or, one might say (pictorially speaking) that the light-ring orbit gets hidden inside the horizon of the final black-hole formed by coalescence. Furthermore, if we apply Eq. (14) using the 2PN approximation $`n=2`$ instead of the 3PN one $`n=3`$, we get the value $`0.28`$ instead of Eq. (16). So at the 2PN order the metric seems to admit a light ring, while at the 3PN order it apparently does not admit any. This erratic behaviour 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 when the masses are equal.

It is impossible of course to be thoroughly confident about the validity of the previous statement because we know only the coefficients up to 3PN order. Any tentative conclusion based on 3PN 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 gravitational field generated by two bodies is more complicated than the Schwarzschild space-time. This appraisal should look cogent to relativists and is in accordance with the author’s respectfulness of the complexity of the Einstein field equations.

We want next to comment on a possible implication of our conclusion as regards 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. (2)-(5)\] would admit a singularity at some reasonable value of $`x`$ (e.g. $`0.5`$). In the Schwarzschild case, for which Eq. (13) holds, the Padé series converges rapidly : the Padé constructed only from the 2PN approximation of the energy — keeping only $`a_1^{\mathrm{Sch}}`$ and $`a_2^{\mathrm{Sch}}`$ — already coincide with the exact result given by Eq. (11). On the other hand, the EOB method maps the post-Newtonian two-body dynamics (at the 2PN or 3PN orders) 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 post-Newtonian resummation techniques is questionable because the value $`\omega _\mathrm{s}=0`$ suggests that the two-body interaction is not Schwarzschild-like. Our doubt is confirmed by the finding of Ref. (already alluded to above) that in the case of the “wrong” ambiguity parameter $`\omega _\mathrm{s}^{}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 light-ring singularity seems to exist with $`\omega _\mathrm{s}^{}`$. By contrast, in the case of general relativity ($`\omega _\mathrm{s}=0`$), the Padé and EOB methods give quite different results (cf. the figure 2 in Ref.). Such a disagreement, we argue, is due to the fact that the assumptions underlying the various resummation techniques are probably not fulfilled, so they may converge toward different solutions. 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}}(\frac{1}{4},\text{Padé})0.44.$$
(17)
This value is rather close to Eq. (15) but strongly disagrees with Eq. (16). Our explanation is that the Padé series converges toward a theory having $`\omega _\mathrm{s}\omega _\mathrm{s}^{}`$ and so which is different from general relativity.

Finally we come to the good news that, if really the (absolute value of the) post-Newtonian coefficients when $`\nu =\frac{1}{4}`$ stay of the order of one as it seems to, this means that the standard post-Newtonian approach, based on the standard Taylor approximants, is probably very accurate. The post-Newtonian series is likely 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 excellent 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.. Therefore the order-of-magnitude estimates we proposed 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.

Compare the situation with the case of Schwarzschild, for which we have seen that $`x_{\mathrm{ICO}}^{\mathrm{Sch}}=\frac{1}{6}`$ is rather close to the value of the convergence radius of the series given by $`x_{\mathrm{light}\mathrm{ring}}^{\mathrm{Sch}}=\frac{1}{3}`$, hence the poor convergence of the post-Newtonian expansion at the location of the ICO. By contrast, when the two masses have the same size, we have $`x_{\mathrm{ICO}}(\frac{1}{4})0.2`$ (see the Figure 1) which is quite far from the “convergence radius” of the series, that we argued is likely to be of the order of one. We thus expect the post-Newtonian expansion (in the case of two black holes) to be well appropriate near the ICO.

## 4 Comparison with a numerical simulation

We confront the prediction of the standard (Taylor-based) post-Newtonian approach with a recent result of numerical relativity by Gourgoulhon, Grandclément and Bonazzola. These authors computed numerically the energy of binary black holes 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. In the numerical approach there are no gravitational waves, the field is periodic in time, and the gravitational potentials tend to zero at spatial infinity within a restricted model equivalent to solving five out of the ten Einstein field equations. Considering an evolutionary sequence of equilibrium configurations Gourgoulhon et al. obtained numerically the circular-orbit energy $`E(\omega )`$ and looked for the innermost circular orbit or ICO.

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 our post-Newtonian formalism the effects of spins appropriate to two Kerr black holes rotating at the orbital rate $`\omega `$. The importance of the effect of spins in corotating systems of neutron stars, for which the ICO is usually determined by the hydrodynamical instability rather than by the effect of general relativity, is well known. We expect that these effects play some role in the case of black holes as well.

The total relativistic mass of the Kerr black hole is given by888In the formulas (18)-(25) we pose $`G=1=c`$.

$$m^2=m_{\mathrm{irr}}^2+\frac{S^2}{4m_{\mathrm{irr}}^2},$$
(18)
where $`S`$ is the spin, related to the usual Kerr parameter by $`S=ma`$, and $`m_{\mathrm{irr}}`$ is the irreducible mass given by $`m_{\mathrm{irr}}=\frac{\sqrt{A}}{4\pi }`$ ($`A`$ is the hole’s surface area). The angular velocity of the corotating black hole is $`\omega =\frac{m}{S}`$ hence, from Eq. (18),

$$\omega =\frac{S}{2m^3\left[1+\sqrt{1\frac{S^2}{m^4}}\right]}.$$
(19)
Physically this angular velocity is the one of the outgoing photons that remain for ever at the location of the light-like horizon. Combining Eqs. (18) and (19) we obtain $`m`$ and $`S`$ as functions of $`m_{\mathrm{irr}}`$ and $`\omega `$ :

$`m`$ $`=`$ $`{\displaystyle \frac{m_{\mathrm{irr}}}{\sqrt{14m_{\mathrm{irr}}^2\omega ^2}}},`$ (20)
$`S`$ $`=`$ $`{\displaystyle \frac{4m_{\mathrm{irr}}^3\omega }{\sqrt{14m_{\mathrm{irr}}^2\omega ^2}}}.`$ (21)
This is the right thing to do since $`\omega `$ is the basic variable describing each equilibrium configuration calculated numerically, and because the irreducible masses are the ones which are held constant along the numerical evolutionary sequences. In the limit of small rotations we have

$$S=I\omega +𝒪\left(\omega ^3\right),$$
(22)
where $`I=4m_{\mathrm{irr}}^3`$ is the moment of inertia of the Kerr black hole. Next the total mass-energy is

$$m=m_{\mathrm{irr}}+\frac{1}{2}I\omega ^2+𝒪\left(\omega ^4\right).$$
(23)
It involves the standard kinetic energy of the spin.

To take into account all the spin effects our first task is to replace all the masses entering the energy function (2)-(5) by their equivalent expressions in terms of $`\omega `$ and the two irreducible masses. It is clear that the leading contribution is that of the spin kinetic energy given by Eq. (23) and comes from the replacement of the rest mass-energy $`Mc^2`$ (where $`M=m_1+m_2`$). From Eq. (23) this effect in the case of corotating binaries is of order $`\omega ^2`$, which means by comparison with Eqs. (1)-(2) that it is equivalent to an “orbital” effect at the 2PN order. Higher-order corrections in Eq. (23), which behave at least like $`\omega ^4`$, will correspond at least to the orbital 5PN order and are negligible for the present purpose. In addition there will be a subdominant contribution, of the order of $`\omega ^{8/3}`$ equivalent to 3PN order, which comes from the replacement of the masses into the “Newtonian” part, proportional to $`x\omega ^{2/3}`$, of the energy $`E`$ \[see Eq. (2)\]. With the 3PN accuracy we do not need to replace the masses that enter into the post-Newtonian corrections in $`E`$, so these masses can be considered to be the irreducible ones in these terms.

Our second task is to include the specific relativistic effects due to the spins, namely the spin-orbit (S.O.) interaction and the spin-spin (S.S.) one. In the case of spins $`S_1`$ and $`S_2`$ aligned parallel to the orbital angular momentum (and right-handed with respect to the sense of motion) the S.O. energy reads

$$E_{\mathrm{S}.\mathrm{O}.}=\mu (M\omega )^{5/3}\left[\left(\frac{4}{3}\frac{m_1^2}{M^2}+\nu \right)\frac{S_1}{m_1^2}+\left(\frac{4}{3}\frac{m_2^2}{M^2}+\nu \right)\frac{S_2}{m_2^2}\right].$$
(24)
Here we are employing the formula given by Kidder, Will and Wiseman who have computed the S.O. contribution as expressed by means of the orbital frequency $`\omega `$ (this is what we need in view of the comparison with the numerical work). The basis of their computation is the work of Barker and O’Connell who obtained the formula given in terms of the orbital separation $`r`$. The derivation of Eq. (24) in Ref. takes correctly into account the fact that the relation between the orbital separation $`r`$ (in a given coordinate system) and the frequency $`\omega `$ depends on the spins. We immediately infer from Eq. (22) that in the case of corotating black-holes the S.O. effect is equivalent to a 3PN orbital effect and thus must be retained with the present accuracy \[with this approximation, the masses in Eq. (24) are the irreducible ones\]. As for the S.S. interaction (still in the case of spins aligned with the orbital angular momentum) it is given by

$$E_{\mathrm{S}.\mathrm{S}.}=\mu \nu (M\omega )^2\frac{S_1S_2}{m_1^2m_2^2}.$$
(25)
The S.S. effect can be neglected here because it is of the orbital order 5PN for corotating systems.

Summaryzing, the contributions due to the corotating spins at the 3PN order are three in all : the main one is that of the spin kinetic energy and arises at the 2PN order; then we have two subdominant contributions at the 3PN order coming respectively from a mixing between the spin kinetic energy and the Newtonian orbital energy, and from the S.O. interaction (24). We can neglect all the other terms. Summing up these three contributions we find that the energy of the corotating spins is (coming back to the notation of Section 2)

$$\mathrm{\Delta }E^{\mathrm{corot}}(x)=Mc^2x\left\{(26\nu )x^2+\left(\frac{18}{3}\nu +13\nu ^2\right)x^3+𝒪(x^4)\right\}.$$
(26)
The complete 3PN energy of the corotating binary is the sum of Eqs. (2)-(5) and (26). Notice that we must now understand all the masses in (2)-(5) and (26) as being the irreducible masses — we no longer indicate the superscripts “irr” —, which for the comparison with the numerical work must be assumed to stay constant when the binary evolves.

The Figure 1 (issued from the work) presents our results for $`E_{\mathrm{ICO}}`$ in the case of irrotational and corotational binaries. Since $`\mathrm{\Delta }E^{\mathrm{corot}}`$, given by Eq. (26), is at least of order 2PN, the result for 1PNcorot is the same as for 1PN in the irrotational case; then, obviously, 2PNcorot takes into account only the leading 2PN corotation effect \[i.e. the spin kinetic energy given by Eq. (23)\], while 3PNcorot involves also, in particular, the corotational S.O. coupling at the 3PN order. In addition we present in Figure 1 the numerical point obtained by numerical relativity under the assumptions of conformal flatness and of helical symmetry. As we can see the 3PN points, and even the 2PN ones, are rather close to the numerical value. The fact that the 2PN and 3PN values are so close to each other is excellent, and confirms the good accuracy of the approximation we concluded in Section 3. In fact one might say that the role of the 3PN approximation is merely to “prove” the value already given by the 2PN one (but of course, had we not computed the 3PN term, we would not be able to trust very much the 2PN value). As expected, the best agreement we obtain is for the 3PN approximation and in the case of corotation : i.e. the point 3PNcorot. However, the 1PN approximation is clearly not precise enough, but this is not surprising in the highly relativistic regime of the ICO.

In conclusion, we find that the location of the ICO computed by numerical relativity, under the helical-symmetry and conformal-flatness approximations, is in good agreement with the post-Newtonian prediction. (See Ref. for the results calculated within the EOB method at the 3PN order, which are close to the ones reported in Figure 1.) 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 disagreed with each other, and do not match with the present 3PN results (see Ref. for further discussion).

Notice in Figure 1 the difference between the energies of the corotational and irrotational configurations which amounts to approximately $`3.5\times 10^3`$ units of $`Mc^2`$. We find that this amount is mainly made of the (positive) kinetic energy of the spins, which is about $`+5.0\times 10^3`$, but that the S.O interaction at the 3PN order, equal to $`1.3\times 10^3`$, reduces slightly the effect, while the coupling between spin kinetic and orbital terms is quite negligible, of the order of $`0.2\times 10^3`$. We feel that these numerical values are physically satisfying, and well in accordance with our order-of-magnitude estimates in Section 3 which showed that in the regime of the ICO the post-Newtonian series “converges well”, with the successive post-Newtonian approximations becoming numerically smaller and smaller. We have checked that the higher-order spin effects like the S.S. interaction \[see Eq. (25)\] which arises at the 5PN order for corotating systems are completely negligible numerically (the S.S. effect is about $`+0.02\times 10^3`$).

Recently Damour et al. computed the effects of spin in corotating black hole binaries within the effective-one-body (EOB) approach. They find a result which differs from our computation presented above. Basically they obtain that the kinetic energy of the spins is about $`+5.0\times 10^3`$, in agreement with ours, but that the S.O. effect is much larger than in our computation, and nearly compensates in their approach the spin kinetic energy. The net result given in Ref. for the total effect of spins is then about $`+1\times 10^3`$, much smaller than our own result. The surprising fact (in the author’s opinion) is that in the EOB method a relativistic effect like the S.O. coupling, which according to the a priori expectation based on the estimates of Section 3 should be small, can nearly cancel out a “Newtonian” term, which one would physically consider to be the dominant one. Although Refs. and are in very good agreement for the ICO of corotational binaries, and in rather good agreement for irrotational binaries, their discrepancy concerning the numerical value of the S.O. effect remains unexplained. But let us try a guess. As we argued in Section 3 the resummation techniques are probably justified only when the post-Newtonian coefficients in the energy function get bigger and bigger as the order of approximation increases. Our guess would be that perhaps the EOB method is somehow forced to attribute in front of the S.O. term a big post-Newtonian coefficient in order to maintain its internal consistency. With an abnormally big coefficient one might explain the discrepancy.

## Acknowledgements

The author is grateful to an anonymous referee of his paper for pointing out the interest of looking at the numerical values of the post-Newtonian coefficients.


format: text

text:
# Large anisotropy in the optical conductivity of YNi2B2C

## Abstract

The optical properties of YNi2B2C are studied by using the first-principles full-potential linearized augmented plane wave (FLAPW) method within the local density approximation. Anisotropic behavior is obtained in the optical conductivity, even though the electronic structure shows 3D character. A large peak in $`\sigma _z`$ is obtained at 2.4 eV. The anisotropic optical properties are analyzed in terms of interband transitions between energy levels and found that the Ni site plays an important role. The electronic energy loss spectroscopy (EELS) spectra are also calculated to help elucidate the anisotropic properties in this system.

Quaternary transition metal borocarbide superconductors with formula RM2B2C (R= Y, rare-earths; M=transition metal) have attracted much attention due to their relatively high superconducting transition temperatures in the intermetallic compounds including exchange enhanced transition metal elements such as Ni, Pd, and Ru . Of special interest is the coexistence and interplay between superconductivity and magnetism in the systems containing rare-earth elements. YNi2B2C with $`T_c`$=15.6 K serves as a reference system for understanding the superconducting mechanism in this family, because there are no complications introduced by the presence of magnetic rare-earth elements.

YNi2B2C crystallizes in a body-centered tetragonal structure with a ThCr2Si2-type (space group $`I4/mmm`$, $`D_{4h}^{17}`$). The crystal structure is layer-like in the $`c`$ axis, reminiscent of the high-$`T_c`$ cuprate superconductors, consisting of two alternating layers of Ni2B2 and YC. Electronic structure studies show three dimensional (3D) character and so it is thought to be a conventional BCS type superconductor with a relatively high density of states (DOS) at the Fermi level . Despite the layered anisotropic crystal structure, some physical properties of YNi2B2C are often considered to be isotropic. Note that the 3D character of the electronic structure does not necessarily mean the isotropy of physical properties. In fact, it has not yet been settled whether the physical properties of YNi2B2C are isotropic or not. Civale et al.  observed an anisotropic effective mass ratio $`\gamma =(m_{max}/m_{min})^{1/2}1.1`$, while Johnston-Halperin et al.  observed an isotropic result, $`\gamma 1.005`$. Fisher et al.  reported isotropic resistivity for non-magnetic YNi2B2C and anisotropic resistivity for magnetic borocarbide superconductors like ErNi2B2C and HoNi2B2C. Rathnayaka et al.  observed an isotropic upper critical field, $`H_{c2}`$, for YNi2B2C, but a small anisotropic $`H_{c2}`$ for LaNi2B2C. In contrast, anisotropic properties were reported for YNi2B2C thin-films in the upper critical field, thermal conductivity, critical field for vortex lattice transition , and paramagnetic susceptibility from nuclear magnetic resonance (NMR) experiments. von Lips et al.  used X-ray absorption spectroscopy (XAS) and the linear combination of atomic-like orbitals (LCAO) band method for YNi2B2C to determine the anisotropic nature of the unoccupied DOS.

The optical properties of YNi2B2C have been studied several times since the first measurement by Widder et al., who obtained anisotropic results by linking together optical reflectance and electron energy loss spectroscopy (EELS) measurements. On the other hand, Bommeli et al.  reported that the structural anisotropy would not affect the optical data based on the observation that both the polycrystalline and single crystal samples produced equivalent experimental results. Therefore, the results were conflicting, and until recently, the anisotropy problem in the optical properties of Ni borocarbides has not been resolved. It is important to understand correct optical properties since the superconducting energy gap is obtained from optical measurements.

Interestingly, all the above optical groups utilize the 3D electronic structures of Ni borocarbides for their experiments. It thus indicates that the electronic structure information only may not be enough to determine definitely the optical anisotropy. In this Letter, to resolve this controversial issue, the optical conductivity and the electron energy loss spectroscopy (EELS) have been studied theoretically for the first time by using first principles calculations. We have found that the optical properties of YNi2B2C are really anisotropic.

Calculations are performed by employing the highly precise first principles local density full-potential linearized augmented plane wave (FLAPW) method, which makes no shape approximations to the charge density or the potential. The Hedin-Lundqvist form is employed for the LDA exchange-correlation energy. We used lattice parameters determined by experiment, 240 special $`𝐤`$-points for the self-consistency cycle and 8619 $`𝐤`$-points for the calculation of the optical properties inside the irreducible Brillouin zone. The optical conductivity arises from the optical transitions of the occupied valence electrons to unoccupied states above EF. Each interband transition is broadened in the calculation by introducing a finite phenomenological inverse life time, $`\mathrm{\Delta }`$=0.05 eV. The Drude term due to free carriers is taken into account by using parameters from experiment.

Figure 1 shows the crystal structure and the first Brillouin zone of YNi2B2C. In Fig. 1(b), points a,b,c, and d are not symmetry points since the symmetry is not a local maximum at those points. Points a and d (and b and c likewise) are the same points by the translational symmetry in $`𝐤`$-space. X has the highest symmetry as $`\mathrm{\Gamma }`$ with 16 symmetry elements.

Figure 2 shows the calculated energy band structures along the symmetry lines. Band dispersions agree well with existing calculations. However, in the present full-potential band calculation, the 19th band at $`\mathrm{\Gamma }`$, which gives a small calabash shaped Fermi surface in the LMTO-ASA calculation, is lifted up to higher energy.

Figure 3 compares the calculated optical conductivities for YNi2B2C with experiments, in which $`\sigma _x`$ and $`\sigma _z`$ represent theoretical conductivities along the $`x`$ and $`z`$ directions, respectively. One can see a clear difference between $`\sigma _x`$ and $`\sigma _z`$, implying the anisotropic nature of the optical conductivity. Anisotropic optical conductivities in the calculation are consistent with the anisotropy observed in the experiment by Widder et al. , represented by $`\sigma _x`$\_w and $`\sigma _z`$\_w in Fig. 3. Peaks are located at 2.48, 7.6, 8.9, and 9.6 eV for $`\sigma _x`$, and at 2.4, 5.0, 8.16, and 9.5 eV for $`\sigma _z`$ in the calculation. The most prominent and anisotropic peak in the calculation is the one at 2.4 eV of $`\sigma _z`$, while $`\sigma _x`$ has a peak at 2.48 eV with a smaller magnitude than $`\sigma _z`$. However, the intensity and position of the peak of $`\sigma _z`$ are somewhat different from those in experiment: the highest peak is located at 1.6 eV, which is lower by 0.8 eV than the theoretical one. It is unusual that the experimental peak is located at lower frequency considering that, in most LDA calculations, theoretical peaks are found at lower frequencies than experimental ones. The disagreement is to be clarified further theoretically or experimentally. Other optical properties such as the reflectivity $`R`$ and the dielectric constants $`ϵ`$ show similar anisotropies which will be discussed elsewhere.

Bommeli et al.  obtained isotropic optical conductivity (denoted by $`\sigma _x`$\_b in Fig. 3). The overall shape of $`\sigma _x`$\_b is seemingly close to $`\sigma _x`$\_w, except that some fluctuations occur when the relative magnitude of $`\sigma _z`$\_w with respect to $`\sigma _x`$\_w becomes large. This feature can be understood if $`\sigma _x`$\_b corresponds to the average of the anisotropic optical conductivity. Although it is not clear in $`\sigma _x`$\_b, there is a peak at around the 2.4 eV in $`\sigma _x`$\_w. The 2.4 eV peak in $`\sigma _x`$ is also observed in a recent experiment by Lee et al.. Hence the experimental peak at 2.4 eV can be assigned to the calculated $`\sigma _x`$ peak at 2.48 eV,

The optical conductivity is obtained directly from the optical transitions between energy levels in the solid. Hence the peaks in the optical conductivity can be assigned to specific interband transitions. However, peaks in derived quantities from the optical conductivity, such as the reflectivity $`R`$ and the refractive index $`n`$, may be shifted from the peaks of the optical conductivity. In order to identify the peaks manifested in the optical conductivity, we plot in Fig. 4 the momentum resolved optical conductivity $`\sigma (𝐤,\omega )`$ along the symmetry lines in the Brillouin zone. Integrating $`\sigma (𝐤,\omega )`$ over $`𝐤`$-space gives rise to the optical conductivity $`\sigma (\omega )`$.

Figures 4 (a) and (b) provide the optical transitions along the $`x`$ and $`z`$ directions, respectively. The darker intensity in Fig. 4 represents the stronger interband transitions. It is evident from Figs. 4(a) and (b) that the optical transitions are anisotropic between the $`x`$ and $`z`$ directions. The anisotropy is due to different matrix element between the $`x`$ and $`z`$ directions. Strong transitions in the $`z`$-direction occur around P and N, while they occur at points near $`\mathrm{\Gamma }`$ in the $`x`$-direction. The anisotropy is most prominent at the 2.4 eV transition near P; there is a strong transition in the $`z`$-direction, but a negligible transition in the $`x`$-direction. The intensity of the 2.4 eV transition is further enhanced by the transition near N. In Fig. 2, corresponding transitions for 2.4 eV are marked as small solid arrows at N and P. Other strong transitions identified in Fig. 4 are also marked as arrows in Fig. 2. It is noteworthy that the unoccupied bands (19th band) for the 2.4 eV transition near P have similar shapes to the occupied bands (18th band), reminiscent of the nesting of bands.

In order to see the detailed nature of the transition, we have calculated contributions from each atomic species for the bands related to the 2.4 eV transition at P. The 18th band consists of Ni-$`d_{3z^2r^2}`$(36%) and Ni-$`d_{xy}`$(35%) states, while the unoccupied band (19th band) consists of Ni-$`d_{xy}`$(16%), Ni-$`s`$(12%), Ni-$`p_z`$(4.4%), B-$`p`$(6.9%), and C-$`p`$(16%) states. Since the occupied band originates mainly from Ni-$`d`$ states, the transition matrix element is determined by Ni site. By the selection rule about the symmetry of the wave functions, one can assign the 2.4 eV peak in $`\sigma _z`$ to the transition between Ni-$`d_{3z^2r^2}`$ and Ni-$`p_z`$. The transition matrix element for $`\sigma _x`$ is zero in this case, which explains the negligible intensity at 2.4 eV near P in $`\sigma _x(𝐤,\omega )`$. Thus Ni states play a vital role in yielding a large anisotropic peak in $`\sigma _z`$. This is indeed consistent with the dominant Ni-$`d`$ DOS near the Fermi level .

Figure 5 depicts the loss function, Im$`(1/ϵ(\omega ))`$, of the EELS measurement with zero momentum transfer. There is no EELS data available for single crystals. Hence, in the figure, the calculated result is compared with Widder et al.’s experiment for polycrystals . There is no plasma peak due to free carriers which is expected to occur between 3.75 eV and 4.25 eV. It is damped due to interband transitions (Landau damping), in agreement with the experimental EELS result. The calculated loss function shows a two peak structure in both the $`x`$ and $`z`$-direction, and the spectrum in the $`z`$-direction is narrower and higher than that of $`x`$ direction. Two peaks are located at 22.7 eV and 26 eV in the z-direction. In YNi2B2C, there are 33 valence electrons per formula unit. The simple plasma frequency relation, $`\omega _p^2=4\pi ne^2/m`$, gives 26.3 eV which agrees well with the above plasma frequency of 26 eV. Therefore, one can explain the plasma frequency at 26 eV as originating from the response of all valence electrons to the external electric field. Then the peak at 22.7 eV can be ascribed to the valence electrons near EF: Ni-$`d`$, Y-$`d`$, B-$`p`$, and C-$`p`$ states. In the experiment, two peaks are observed at $``$ 20 eV and $``$ 35 eV. Although the two peak structure in the calculation agrees qualitatively with experiment, the peak positions are different. The EELS experiment with single crystal is required to clarify the discrepancy between theory and experiment.

In conclusion, we have investigated the optical properties of YNi2B2C by using the first-principles FLAPW method. The optical conductivity is found to be anisotropic between the $`x`$ and the $`z`$ direction. A strong peak at 2.4 eV is obtained in the $`z`$ direction, which is assigned to the transition near P between Ni-$`d_{3z^2r^2}`$ and Ni-$`p_z`$ states, indicating that the Ni-site plays an important role in the anisotropy of YNi2B2C. The calculated EELS spectra are also found to be anisotropic with a two peak structure.

###### Acknowledgements.

The authors are grateful to KISTEP/CRI of Korean ministry of Science and Technology for financial support. BIM is also supported by the KOSEF through the eSSC at POSTECH. SJY thanks to Prof. J. H. Kim for helpful discussions. Work at Northwestern University supported by the U.S. Department of Energy (Grant No. DE-F602-88ER45372).


text:
# Phase-fluctuation induced reduction of the kinetic energy at the superconducting transition

## I Introduction

The key idea of the phase-fluctuation scenario in the high-$`T_c`$ superconductors is the notion that the pseudogap observed in a wide variety of experiments arises from phase fluctuations of the superconducting gap em.ki.95 ; ra.tr.92 ; fr.mi.98 ; kw.do.99 ; herb.02 ; fr.te.01 ; ec.sc.02 . In this scenario, below a mean-field temperature scale $`T_c^{MF}`$, a $`d_{x^2y^2}`$-wave gap amplitude is assumed to develop. However, the superconducting transition is suppressed to a considerably lower transition temperature $`T_c`$ by phase fluctuations em.ki.95 ; ec.sc.02 . In the intermediate temperature regime between $`T_c^{MF}`$ and $`T_c`$, phase fluctuations of the superconducting order parameter give rise to the pseudogap phenomena. Recently, we have shown that indeed a two-dimensional BCS-like Hamiltonian with a $`d_{x^2y^2}`$-wave gap and phase fluctuations, which were treated by a Monte-Carlo simulation of an $`XY`$ model, yields results which compare very well with scanning tunneling measurements over a wide temperature range ec.sc.02 ; ku.fi.01 . Thus, they support the phase-fluctuation scenario for the pseudogap.

There is also increasing evidence from a number of recent experiments for the relevance of phase fluctuations such as the Corson et al. measurements co.ma.99  of the high-frequency conductivity that track the phase correlation time. $`XY`$ vortices are probably responsible for the large Nernst effect wa.xu.00 ; wa.on.02 . The evolution of $`T_c`$ with electron irradiation found very recently ru.al.03u  also emphasizes the importance of phase fluctuations. In this paper, we argue that phase fluctuations should have yet another rather unexpected consequence: they induce a reduction of the kinetic energy at the superconducting transition. This reduction is due to a transition from a “disordered”, i. e. phase-incoherent Cooper-pair motion in the pseudogap regime above $`T_c`$ to an “ordered”, i. e. phase-coherent motion at $`T_c`$. Comparison of our results, based on the BCS phase-fluctuation model, with optical experiments mo.pr.02 ; sa.lo.01u  support this idea.

In ordinary BCS superconductors the optical conductivity is suppressed at frequencies within a range of about twice the superconducting gap. The corresponding low-frequency spectral weight $`W_{low}`$ is transfered to the zero-frequency delta peak $`W_D`$ tinkham , associated with the dissipationless transport (and the superfluid weight $`D`$) in the superconducting state. This is the Glover-Ferrell-Tinkham (GFT) sum rule. On the other hand, the total frequency integral of the optical conductivity is conserved, when decreasing the temperature across the superconducting transition, due to the f-sum rule tinkham , i. e. $`W_{tot}^{sc}=W_{tot}^n`$.

However, recent measurements of the in-plane optical conductivity mo.pr.02 ; sa.lo.01u  have indicated a *violation* of the GFT optical sum rule for frequencies up to $`2eV`$ in underdoped $`Bi_2Sr_2CaCu_2O_{8+\delta }`$ (Bi2212). By entering the superconducting state, not only spectral weight $`W_{low}`$ from the microwave and far-infrared, but also from the visible optical spectrum i. e. high-frequency spectral weight $`W_{high}`$ contributes to the superfluid condensate $`W_D`$. That is, in contrast to ordinary BCS-superconductors, a “*color* change” is introduced at the superconducting transition. The interpretation of this unusual result may require the inclusion of local-field effects and other (such as excitonic) many-body effects. They are known to play a crucial role already in weakly-correlated systems (such as semiconductors), and introduce a shift of order of the Coulomb correlation energy between single-particle and two-particle, i. e. optical excitations ha.sh.80 . Therefore, they may partly account for the “high-energy” features observed in $`\sigma (\omega )`$. On the other hand, within a tight-binding one-band model, the anomalously large energy scale, which contributes to the superfluid weight, and the corresponding *color* change can be attributed to a *reduction of kinetic energy* hirs.02  at the superconducting transition. This is rather surprising, since one would expect that in a conventional (BCS) pairing process, it is the potential energy which is reduced at the expense of the kinetic energy, with the latter being increased due to particle-hole mixing.

The full optical integral, when integrated over all frequencies and *energy bands*, is proportional to the carrier density ($`n`$) over the bare mass ($`m`$)

$$W_{tot}W_{low}+W_D+W_{high}=_0^{\mathrm{}}Re\sigma _{xx}(\omega )𝑑\omega =\frac{ne^2}{2m},$$
(1)
and, thus, is conserved. When the optical integral is restricted over a finite (low) range of frequencies $`\mathrm{\Omega }`$, in the HTSC typically of the order of eV, one may consider the weight $`W_{low}+W_D`$ as being essentially due to a single band around the Fermi energy, i. e.

$$W_{low}+W_D=_0^{\mathrm{}}Re\stackrel{~}{\sigma }_{xx}(\omega )𝑑\omega =(\pi e^2a^2/2\mathrm{}^2V)E_K,$$
(2)
where $`\stackrel{~}{\sigma }`$ is the single-band conductivity, $`a`$ the lattice constant and $`V`$ the unit cell volume. With this single-band assumption, the frequency integral of the optical conductivity is proportional to the the inverse mass tensor ($`\frac{^2ϵ_k}{k_x^2}`$, $`\widehat{x}`$ being the direction in which the conductivity is measured) weighted with the momentum distribution $`n_k`$ no.pe.02u ; kubo.57 :

$$E_K=(2/a^2N)\underset{k}{}\frac{^2ϵ_k}{k_x^2}n_k,$$
(3)
with $`N`$ the number of k points. This quantity depends upon the bare single-particle band structure $`ϵ_k`$, being proportional to minus the kinetic energy $`E_K=E_{kin}`$ for a (nearest-neighbor) tight-binding (TB) model, while for free electrons it is a constant given by the electron density divided by the effective mass.

## II Phase Fluctuation Scenario for Kinetic Energy Reduction

In this paper, we propose as a mechanism for a kinetic-energy reduction phase-fluctuations. That is, in order to have condensation into the superconducting state, one needs, in addition to the binding of charge carriers into Cooper pairs, long-range phase coherence among the pairs. Since superconductors with low superconducting carrier density (such as the organic and underdoped high-$`T_c`$ superconductors) are characterized by a relatively small phase *stiffness*, this implies a significantly larger role for phase fluctuations, than in conventional superconductors em.ki.95 ; ca.ki.99 ; klei.00 . As a consequence, in these materials the transition to the superconducting state does not display a typical mean-field (BCS) behavior, and phase fluctuations, both classical and quantum, may have a significant influence on low temperature properties. When coherence is lost due to thermal fluctuations of the phase at and above the transition temperature $`T_c`$, pairing remains, together with short-range phase correlations. These phase fluctuations can cause the pseudogap phenomena observed e. g. in tunneling experiments re.re.p.98 ; mi.za.99 ; ku.fi.01 ; ec.sc.02  in the underdoped HTSC.

We show here that indeed phase fluctuations contribute to a significant reduction of the in-plane kinetic energy upon entering into the superconducting phase below $`T_c`$, with a magnitude comparable to recent experimental results. The physical reason for this kinetic energy lowering is that, due to phase fluctuations and to the associated incoherent motion of Cooper pairs (cf. Fig. 1), the pseudogap region has a higher kinetic energy than the simple BCS mean-field state. When long-ranged phase coherence finally develops at Tc , the Cooper-pair motion becomes phase coherent and the kinetic energy decreases. The onset of the coherent motion can be seen, for example, from the development of coherence peaks in the tunneling spectrum of BiSrCaCuO compounds (see e. g. ku.fi.01 ; ec.sc.02  and Fig. 3). The initial *cost* of kinetic energy, which is needed for pairing, is payed at a mean-field temperature $`T_c^{MF}`$ considerably higher than $`T_c`$. Therefore, the reduction of kinetic energy observed experimentally mo.pr.02 ; sa.lo.01u  can be attributed to a transition from a phase-disordered pseudogap to a phase-ordered superconducting state. We stress that this effect is independent of the particular mechanism leading to pair formation, as long as the superconductor considered is characterized by a small phase stiffness kivelson .

Our starting Hamiltonian is of a simple BCS form given by

$$H=K\frac{1}{4}\underset{i\delta }{}(\mathrm{\Delta }_{i\delta }\mathrm{\Delta }_{i\delta }^{}+\mathrm{\Delta }_{i\delta }^{}\mathrm{\Delta }_{i\delta }),$$
(4)
with the nearest-neighbor hopping term

$$K=t\underset{ij,\sigma }{}(c_{i\sigma }^{}c_{j\sigma }+c_{j\sigma }^{}c_{i\sigma }),$$
(5)
where $`c_{i\sigma }^{}`$ creates an electron of spin $`\sigma `$ on the $`i^{th}`$ site and $`t`$ denotes an effective nearest-neighbor hopping. The $`ij`$ sum is over nearest-neighbor sites of a 2D square lattice, and, in the pairing term, $`\delta `$ connects $`i`$ to its nearest-neighbor sites. The local $`d`$-wave gap,

$$\mathrm{\Delta }_{i\delta }^{}=\frac{1}{\sqrt{2}}c_i^{}c_{i+\delta }^{}c_i^{}c_{i+\delta }^{}=\mathrm{\Delta }e^{i\mathrm{\Phi }_{i\delta }},$$
(6)
is characterized by the *fluctuating* phases

$$\mathrm{\Phi }_{i\delta }=\{\begin{array}{cc}(\phi _i+\phi _{i+\delta })/2\mathrm{for}\hfill & \delta \text{ in x-direction}\hfill \\ (\phi _i+\phi _{i+\delta })/2+\pi \mathrm{for}\hfill & \delta \text{ in y-direction,}\hfill \end{array}$$
(7)
and by a spatially constant amplitude $`\mathrm{\Delta }`$. We neglect the relative bond phase fluctuations between $`\delta =\widehat{x}`$ and $`\widehat{y}`$ as well as amplitude fluctuations. Thus, we consider only center of mass pair phase fluctuations, which are the relevant low-energy degrees of freedom, in a situation in which the superfluid density is small, like in the underdoped cuprates. In our two-dimensional (2D) model $`T_c`$ corresponds to the Kosterlitz-Thouless transition temperature $`T_{KT}`$, where the phase correlation length $`\xi `$ diverges.

With this Hamiltonian it is straightforward to show, that the optical sum rule yields sc.wh.93 ; kubo.57 ; delta ,

$$_0^{\mathrm{}}Re\stackrel{~}{\sigma }_{xx}(\omega )𝑑\omega =e^2\pi k_x/2,$$
(8)
in units where $`\mathrm{}=c=1`$, with $`k_x`$ being the expectation value of the nearest neighbor hopping kin  in the $`x`$-direction, i. e.

$$k_x=t\underset{\sigma }{}c_{i\sigma }^{}c_{i+x\sigma }+c_{i+x\sigma }^{}c_{i\sigma }.$$
(9)
Since we are only interested in the temperature region $`TT_c`$, one can safely assume that the fluctuations of the phase $`\phi _i`$ are predominantly determined by a classical $`XY`$ free energy em.ki.95 ; ca.ki.99 ,

$$F\left[\phi _i\right]=J\underset{ij}{}\mathrm{cos}\left(\phi _i\phi _j\right).$$
(10)
Our physical picture here is that the XY-action arises from integrating out the shorter wavelength fermion degrees of freedom, including those responsible for the formation of the local pair amplitude and the internal $`d_{x^2y^2}`$ structure of the pair. Thus, the *scale* of the XY-lattice spacing is actually set by the pair coherence length $`\xi _0`$. In our work, we have chosen $`\mathrm{\Delta }`$ so that $`\xi _0\frac{v_F}{\pi \mathrm{\Delta }}1`$. In this case, the phase configurations $`\phi _i`$ calculation can be carried out on the same $`L\times L`$ $`(L=32)`$ lattice that is used for the diagonalization of the Hamiltonian. This allows the Kosterlitz-Thouless phase correlation length $`\xi `$ to grow over a sufficient range as $`T`$ approaches $`T_{KT}`$ and minimizes finite-size effects. Thus, we are always in the limit where the phase correlation length $`\xi `$ is larger than the Cooper pair size $`\xi _0`$, when the temperature T is below the mean-field critical temperature $`T_c^{MF}`$.

In principle, also the coupling energy $`J`$ can be considered as arising from integrating out the high-energy degrees of freedom of the underlying microscopic system. Here, we will proceed phenomenologically, neglecting the temperature dependence of $`J`$ and simply use it to set the Kosterlitz-Thouless transition temperature $`T_{KT}`$ equal to some fraction of $`T_c^{MF}`$. Specifically, for the present calculations we will set $`T_{KT}\frac{1}{4}T_c^{MF}`$. This choice is motivated by the recent scanning-tunneling results in $`Bi_2Sr_2CuO_{6+\delta }`$, where $`T_c10K`$ and the pseudogap regime extends to about 50K, which we take as $`T_c^{MF}`$.

In a previous paper ec.sc.02  we have presented a detailed numerical solution of the 2D-BCS like Hamiltonian of Eq. 4 with a $`d`$-wave gap and phase fluctuations. This is a minimal model but, nevertheless, contains the key ideas of the cuprate phase fluctuation scenario: that is, a $`d`$-wave BCS gap amplitude forms below a mean-field temperature $`T_c^{MF}`$, but phase fluctuation suppress the actual transition to a considerably lower temperature $`T_c`$. In the intermediate temperature regime between $`T_c^{MF}`$ and $`T_c`$, the phase fluctuations of the gap give rise to pseudogap phenomena. Comparison of these results with recent scanning tunneling spectra of Bi-based high-$`T_c`$ cuprates supports the idea that the pseudogap behavior observed in these experiments can be understood as arising from phase fluctuations ec.sc.02 .

In the present calculations, where we assume a BCS temperature dependence of the pairing gap $`\mathrm{\Delta }(T)`$, we have therefore set $`\mathrm{\Delta }(T=0)=1.0t`$ corresponding to $`T_c^{MF}0.42t`$ and selected $`J`$ so that $`T_{KT}=0.1t`$ fe.fe.86n . The condition that $`\xi >\xi _0`$ is thus always fulfilled, if we are not too close to $`T_c^{MF}`$. The calculation of the kinetic energy for an $`L\times L`$ $`(L=32)`$ periodic lattice now proceeds as follows ref1 ; da.yu.98n : a set of phases $`\{\phi _i\}`$ is generated by a Monte Carlo importance sampling procedure, in which the probability of a given configuration is proportional to $`\mathrm{exp}(F[\phi _i]/T)`$ with $`F`$ given by Eq. (10). With $`\{\phi _i\}`$ given, the Hamiltonian of Eq. (4) is diagonalized and the kinetic energy $`E_{kin}(T,\{\phi _i\})=k_x_{\{\phi _i\}}`$ is extracted. Further Monte-Carlo $`\{\phi _i\}`$ configurations are generated and an average kinetic energy $`E_{kin}(T)=k_x`$, at a given temperature, is determined.

Fig. 1 displays the kinetic energy $`k_x`$ as a function of temperature for non-interacting tight-binding electrons, for BCS electrons, and for our phase-fluctuation model, respectively.

We can clearly see that pairing, as expected, produces an overall increase of kinetic energy (indicated as vertical arrows) with respect to the free-electron case. We observe that in the phase-fluctuation model the kinetic energy is further increased (small vertical arrows) due to the incoherent motion of the paired electrons. The kinetic energy is a smoothly decreasing function of temperature for $`T0`$. This is expected from the fact that, at high temperature, more electrons are transferred to higher kinetic energies, and is in agreement with the experimental results mo.pr.02 ; sa.lo.01u . What we are especially interested in, is the rather pronounced change (magnified in Fig. 2 by using a different scale for the kinetic energy) near $`T_cT_{KT}`$, where the kinetic energy of our phase-fluctuation model rather suddenly reduces to the BCS value. This sudden deviation from the $`TT_c`$ behavior is also obtained in experiments, which show a kink in the temperature dependence of the low-frequency spectral weight $`W_{low}+W_D`$ at $`T_c`$ mo.pr.02 .

This pronounced change of in-plane kinetic energy can be better observed in Fig. 2, where we plot the difference between the BCS kinetic energy and the kinetic anergy of our phase-fluctuation model $`\delta k_x=k_x_{PP}k_x_{BCS}`$.

As discussed above, this reduction is due to the onset of phase coherence of the Cooper pairs below the superconducting transition temperature $`T_cT_{KT}`$. This is signaled by the appearance of sharp coherence peaks in the single-particle spectral function upon developing long-range phase coherence ec.sc.02 . The corresponding result for the density of states $`N(\omega )`$ displaying these coherence peaks is shown in Fig. 3.

Notice that this argument for the reduction of kinetic energy at Tc due to a phase ordering transition is quite robust. For example, we expect it to be valid (and actually to be stronger) in a true three-dimensional system. As a matter of fact, it has been argued em.ki.95 ; ca.ki.99  that even small interplane couplings play an important role due to the infinite-order nature of the the KT transition.

In order to get a rough estimate of the kinetic condensation energy, we calculate the reduction in kinetic energy near $`T_c`$, i. e.

$$\mathrm{\Delta }E_k=\frac{2}{e^2\pi }_0^{\mathrm{}}(Re\stackrel{~}{\sigma }_{xx}^{sc}(\omega )Re\stackrel{~}{\sigma }_{xx}^n(\omega ))𝑑\omega ,$$
(11)
as indicated by the energy change $`\mathrm{\Delta }E_k`$ in figure 2. Assuming that $`t250meV`$, we get a condensation energy estimate of $`1.5meV`$ per Copper site, which is in order of magnitude agreement with the experimental results (again assuming a one-band TB analysis).

Up to now, to refrain from further approximations, we have set the chemical potential $`\mu `$ equal to zero and have only considered nearest-neighbor hopping. We have checked to some extent, how robust these results are with respect to finite doping ($`n0.9`$) and the inclusion of a next-nearest neighbor-hopping term $`t^{}`$ in our Hamiltonian Eq. 4. Notice that in this case $`E_k`$ is no longer proportional to the kinetic energy. For $`t^{}0.3t`$, our results for the sum rule *violation* are reduced only by about $`20\%30\%`$.

## III Summary

In conclusion, we have shown that the recently observed *violation* of the low-frequency optical sum rule in the superconducting state, associated with a reduction of kinetic energy, can be related to the role of phase fluctuations. The decrease in kinetic energy is due to the sharpening of the quasiparticle peaks close to the superconducting transition at $`T_cT_{KT}`$, where the phase correlation length $`\xi `$ diverges. We suggest that this sum rule violation should also appear in other superconductors with low charge carrier density (phase stiffness) such as the organic superconductors.

## Acknowledgments

We would like to acknowledge useful discussions and comments by S. A. Kivelson and D. J. Scalapino. This work was supported by the DFG under Grant No. Ha 1537/16-2 and by a Heisenberg fellowship (AR 324/3-1), by the Bavaria California Technology Center (BaCaTeC), the KONWHIR projects OOPCV and CUHE. The calculations were carried out at the high-performance computing centers HLRS (Stuttgart) and LRZ (München).