format: text

text:
# Valence quark distributions in nucleon at low 𝑄² in QCD.

## Abstract

Valence $`u`$\- and $`d`$-quarks distributions in proton are calculated in QCD at low $`Q^2`$ and intermediate $`x`$, basing on the operator product expansion (OPE). The imaginary part of the virtual photon scattering amplitude on quark current with proton quantum numbers is considered. The initial and final virtualities $`p_1^2`$ and $`p_2^2`$ of the currents are assumed to be large, negative and different, $`p_1^2p_2^2`$. The OPE in $`p_1^2`$, $`p_2^2`$ up to dimension 6 operators was performed. Double dispersion representations in $`p_1^2,p_2^2`$ of the amplitudes in terms of physical states contributions are used. Putting them to be equal to those calculated in QCD, the sum rules for quark distributions are found. The double Borel transformations are applied to the sum rules. Leading order perturbative corrections are accounted. Valence quark distributions are found: $`u(x)_v`$ at $`0.15<x<65`$$`d(x)_v`$ at $`0.25<x<0.55`$ with an accuracy $`30\%`$ in the middles and $`50\%`$ at the ends of these intervals. The quark distributions obtained are in agreement with those found from the analysis of hard processes data.

PACS: 11.15.Tk; 11.55.Hx; 12.38.Lg.

Keywords: Quantum chromodynamics; Quark distributions; Sum rules.

1.Introduction

Quark and gluon distributions in hadrons are not fully understood in QCD. QCD predicts the evolution of these distributions with $`Q^2`$ in accord with the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) - equations, but not the initial values from which this evolution starts. The standard way of determination of quark and gluon distributions in nucleon is the following \[4-8\] (for the recent review see ). At some $`Q^2=Q_0^2`$ (usually, at low or intermediate $`Q^225GeV^2`$) the form of quark (valence and sea) and gluon distributions is assumed and characterized by the number of free parameters. Then, by using DGLAP equations, quark and gluon distributions are calculated at all $`Q^2`$ and $`x`$ and compared with the whole set of the data on deep inelastic lepton-nucleon scattering (sometimes also with prompt photon production, jets at high $`p_{}`$ etc). The best fit for the parameters is found and, therefore, quark and gluon distributions are determined at all $`Q^2`$, including their initial values $`q(Q_0^2,x)`$, $`g(Q_0^2,x)`$. Evidently, such an approach is not completely satisfactory from theoretical point of view - it would be desirable to determine the initial distribution directly from QCD. Also, if the form of initial distribution is not divined properly, then the fitting procedure will give the relative minimum, but not the absolute ones and the results, especially at low $`Q_0^2`$ could be wrong. (This danger does not, probably, exists for nucleon, but such a situation happens, e.g., for transversally polarized $`\rho `$-meson, where unusual form of initial valence quark distribution was found ). Finally, extrapolation from high $`Q^2`$, to which the main part of the data belongs, to low $`Q^2`$ is not a quite stable procedure and may introduce some errors. For all these reasons it is desirable to find quark and gluon distribution in hadrons at low $`Q^225GeV^2`$ basing directly on QCD.

In this paper we calculate $`u`$ and $`d`$ valence quark distributions in proton. The idea of the method was suggested in  and developed in \[12-14\]. Recently, the method had been improved and valence quark distributions in pion  and transversally and longitudinally polarized $`\rho `$-meson  had been calculated, what was impossible in the initial version of the method. The idea of the approach (in the improved version) is to consider the imaginary part (in $`s`$-channel) of a four-point correlator $`\mathrm{\Pi }(p_1,p_2,q,q^{})`$ corresponding to the non-forward scattering of two quark currents, one of which has the quantum numbers of hadron of interest (in our case – of proton) and the other is electromagnetic (or weak). It is supposed that virtualities of the photon $`q^2,q^2`$ and hadron currents $`p_1^2,p_2^2`$ are large and negative $`|q^2|=|q^2||p_1^2|,|p_2^2|R_c^2`$, where $`R_c`$ is the confinement radius. It was shown in  that in this case the imaginary part in $`s`$-channel $`[s=(p_1+q)^2]`$ of $`\mathrm{\Pi }(p_1,p_2;q_1,q^{})`$ is dominated by a small distance contribution at intermediate $`x`$. (The standard notation is used: $`x`$ is the Bjorken scaling variable, $`x=q^2/2\nu `$ , $`\nu =p_1q)`$. The proof of this statement is given in ref.. So, in the mentioned above domain of $`q^2,q^2`$$`p_1^2,p_2^2`$ and intermediate $`x`$ $`Im\mathrm{\Pi }(p,p_2;q,q^{})`$ can be calculated using the perturbation theory and the operator product expansion in both sets of variables $`q^2=q^2`$ and $`p_1^2,p_2^2`$. The approach is inapplicable at small $`x`$ and $`x`$ close to 1. This can be easily understood for physical reasons. In deep inelastic scattering at large $`|q^2|`$ the main interaction region in space-time is the light-cone domain and longitudinal distances along the light-cone are proportional to $`1/x`$ and become large at small $`x`$ . For OPE validity it is necessary for these longitudinal distances along light-cone to be also small, that is not the case at small $`x`$. At $`1x1`$ another condition of applicability of the method is violated. The total energy square $`s=Q^2(1/x1)+p_1^2`$ $`Q^2=q^2`$ is not large at $`1x1`$. Numerically, the typical values to be used below are $`Q^25GeV^2`$, $`p_1^21GeV^2`$. Then, even at $`x0.7`$, $`s1GeV^2`$, i.e., at such $`x`$ we are in the resonance, but not in the scaling region. So, one may expect beforehand, that our method could work only up to $`x0.7`$. The inapplicability of the method at small and large $`x`$ manifests itself in the blow-up of higher order terms of OPE. More precise limits on the applicability domain in $`x`$ will be found from the magnitude of these terms.

The further procedure is common for QCD sum rules. On one hand the four-point correlator $`\mathrm{\Pi }(p_1,p_2;q,q^{})`$ is calculated by perturbation theory and OPE.On the other hand, the double dispersion representation in $`p_1^2,p_2^2`$ in terms of physical states contributions is written for the same correlator and the contribution of the lowest state is extracted using the Borel transformaion. By equalling these two expression the desired quark distribution is found.

In Sec.2 the general outline of the method is presented. It is explained here, why it is necessary to consider at the beginning the nonequal $`p_1^2`$ and $`p_2^2`$ and only at the end of the calculation to go to the forward scattering amplitude. Sec.3 presents the results of the bare loop calculation ($`d=0`$ term in OPE) and the LO perturbative corrections to it. Sec.4 is devoted to the calculation of power corrections to the QCD side of the sum rules: the contribution of gluonic condensate (d=4), $`\alpha _s\overline{\psi }\psi ^2`$ term (d=6). The valence $`u`$ and $`d`$-quark distributions are numerically calculated in Sec.5, the errors are estimated and the results are compared with those obtained from the deep inelastic scattering data using the evolution equations. Sec.6 contains our conslusion.

2. The outline of the method

Consider the 4-current correlator which corresponds to the virtual photon scattering on the quark current with quantum number of proton:

$$T^{\mu \nu }(p_1,p_2,q,q^{})=id^4xd^4yd^4ze^{i(p_1x+qyp_2z)}$$

$$0|T\{\eta (x),j_\mu ^{u,d}(y),j_\nu ^{u,d}(0),\overline{\eta }(z)\}|0,$$
(1)
where $`\eta (x)`$ is the three-quark current (Ioffe current ). Choose the currents in the form $`j_\mu ^u=\overline{u}\gamma _\mu u`$$`j_\mu ^d=\overline{d}\gamma _\mu d`$, i.e. as an electromagnetic current which interacts only with $`u(d)`$ quark (with unit charges). Such a choice allows us to get sum rules separately for distribuion functions of $`u`$ and $`d`$ quarks. The general method of calculation of distribution functions from consideration of the 4-point correlator in the QCD sum rules for light quarks was used in . Unfortunately, this approach being in the form used in  did not give completely satisfactory sum rules for quark distribution functions (especially for $`d`$-quark in nucleon) and did not allow one to calculate quark distributions in mesons. Recently, in ref.  was suggested a generlized method of calculation of quark distribution functions, which made it possible to get much more reliable sum rules, what, in particular, made it possible to find valence quark distributions in $`\pi `$ and $`\rho `$\- mesons. The main difference of this method is that the hadronic current momenta are put to be unequal, $`p_1^2p_2^2`$ and then independent borelization over $`p_1^2`$ and $`p_2^2`$ is performed and only at the very end the Borel the parameters $`M_1^2`$ and $`M_2^2`$ are put to be equal. (Unlike the approach of ref. where $`p_1`$ was put to be equal to $`p_2`$ and single borelization is performed from the very beginning). The described procedure allows one to kill nondiagonal transitions of the type

$$0|j^h|h^{}h^{}|j_\mu ^{el}(y)j_\nu ^{el}(0)|hh|j^h|0$$
(2)
and thus makes it possible to separate the diagonal transition of interest

$$0|j^h|hh|j_\mu ^{el}(y)j_\nu ^{el}(0)|hh|j^h|0,$$
(3)
without using additional tricks like differentiation over Borel parameter $`M^2`$, what strongly worsens the accuracy of the sum rules.

Let us briefly remind the main points of derivation of the sum rules in the method under investigaion (for details see . The most general form of the double dispersion relation (in $`p_1^2`$ and $`p_2^2`$) for imaginary part of the correlator (1) (we omit for some time all indeces assuming the necessary invariant amplitude to be chosen) has the form

$$ImT(p_1^2,p_2^2,q^2,s)=a(q^2,s)+\underset{0}{\overset{\mathrm{}}{}}\frac{\psi (q^2,s,u)}{up_1^2}𝑑u+\frac{\psi (q^2,s,u)}{up_2^2}𝑑u$$

$$+\underset{0}{\overset{\mathrm{}}{}}𝑑u_1\underset{0}{\overset{\mathrm{}}{}}𝑑u_2\frac{\rho (q^2,s,u_1,u_2)}{(u_1p_1^2)(u_2p_2^2)}$$
(4)
(without loss of generality one may put $`q^2=q^2`$, $`t=(p_1p_2)^2=0`$). The double Borel transformation in $`p_1^2`$ and $`p_2^2`$ eliminates three first terms and we have

$$B_{M_1^2}B_{M_2^2}ImT(p_1^2,p_2^2,q^2,s)=\underset{0}{\overset{\mathrm{}}{}}𝑑u_1\underset{0}{\overset{\mathrm{}}{}}𝑑u_2\rho (q^2,s,u_1,u_2)exp\left[\frac{u_1}{M_1^2}\frac{u_2}{M_2^2}\right]$$
(5)
where $`M_1^2`$ and $`M_2^2`$ are the squared Borel mass. The integration region with respect to $`u_1,u_2`$ may be divided into four areas

I. $`u_1<s_0`$$`u_2<s_0`$;

II. $`u_1<s_0`$, $`u_2>s_0`$;

III. $`u_1>s_0`$, $`u_2<s_0`$;

IV. $`u_1`$, $`u_2>s_0`$.

Here $`s_0`$ is the continuum threshold in the standard QCD sum rule model of the hadronic spectrum with one lowest resonance plus continuum. Area I obviously corresponds to the resonance contribution and spectral density in this area can be written as

$$\rho (u_1,u_2,x,Q^2)=g_h^22\pi F_2(x,Q^2)\delta (u_1m_h^2)\delta (u_2m_h^2),$$
(6)
where $`g_h`$ are coupling constants of the corresponding hadronic current. As for the areas II-IV, they are exponentially suppressed, and, using the standard hypothesis of quark-hadron duality, we may estimate them as a bare loop contribution in the same integration region.

Thus, the general form of the sum rule is

$$ImT_{QCD}^0+\text{Power correction}=2\pi F_2(x,Q^2)g_h^2e^{m_h^2(\frac{1}{M_1^2}+\frac{1}{M_2^2})}$$

$$ImT_{QCD}^0=\underset{0}{\overset{s_0}{}}\underset{0}{\overset{s_0}{}}\rho ^0(u_1,u_2,x)𝑑u_1𝑑u_2e^{(\frac{u_1}{M_1^2}+\frac{u_2}{M_2^2})}$$
(7)
where one should go to the limit $`M_1^2=M_2^2`$.

We choose $`M_1^2=M_2^2=2M^2`$ following  where it was shown that the value of the Borel mass square in the double sum rules is approximately twice as large than in single ones.

Choose now invariant amplitude for the case of the proton current of interest in the correlator (1) (cf.). Its imaginary (in $`s`$) part can be written as

$$ImT_{\mu \nu }^{(p)}=\lambda _n^2\frac{1}{p_1^2m^2}\underset{r,r^{}}{}v^r(p_1)\times $$

$$\times Im\{id^4xe^{iqy}p_1,r|T\{j_\mu (y),j_\nu (0)\}|p_2,r^{}\}\overline{v}^r^{}(p_2)\frac{1}{p_2^2m^2}$$
(8)
where $`v^r(p)`$ is proton spinor with polarization $`r`$ and momentum $`p`$, $`\lambda _N`$ is the coupling constant of proton with the current $`0|\eta |p,r=\lambda _Nv^r(p)`$.

In order to choose the most suitable invariant amplitude, rewrite eq.(8) in the form

$$ImT_{\mu \nu }=\lambda _N^2\frac{1}{(p_1^2m^2)(p_2^2m^2)}\underset{r,r^{}}{}v^r(p_1)$$

$$\left(\overline{v}^r(p_1)Im\stackrel{~}{T}_{\mu \nu }^{(p)}v^r^{}(p_2)\right)\overline{v}^r(p_2)$$
(9)
where $`Im\stackrel{~}{T}_{\mu \nu }^{(p)}`$ is the amplitude before averaging in proton spin, $`m`$ is the proton mass.

The general form of this amplitude is

$$Im\stackrel{~}{T}_{\mu \nu }^{(p)}(p_1,p_2)=C_1P_\mu P_\nu +C_2(P_\mu \gamma _\nu +P_\mu \gamma _\nu )+C_3\widehat{P}P_\mu P_\nu +\mathrm{}$$

$$+(\text{terms with}r)$$
(10)
where $`P=\frac{p_1+p_2}{2}`$; $`r=p_1p_2`$.

Let us now take into account that we are interested in such a combination of invariant amplitudes which in the spin-averaged matrix element

$$\underset{r}{}\overline{v}^r(p_1)Im\stackrel{~}{T}_{\mu \nu }^{(p)}v^r(p_2)$$
(11)
appears at the kinematical structure $`P^\mu P^\nu `$ (since this structure in the limit $`p_1p_2p`$ transforms into $`p^\mu p^\nu `$, the coefficient at which is just $`F_2(x)`$).

Using the equation of motion

$$\widehat{p}_{1,2}v^r(p_{1,2})=\frac{1}{2}(\widehat{p}_{1,2}+m)v^r(p_{1,2})$$
(12)
and

$$\underset{r}{}v_\alpha ^r(p_{1,2})\overline{v}_\beta ^r(p_{1,2})=(\widehat{p}_{1,2}+m)_{\alpha \beta }$$
(13)
it can be seen that the combination of invariant amplitudes in eq.(10), which appears in eq.(9) at the kinematical structure $`\widehat{P}P^\mu P^\nu `$ coincides (up to numerical factor) with the combination of invariant functions at the structure $`P_\mu P_\nu `$ in (11).

Thus, we come to a conclusion that the sum rules should be written for invariant amplitude at the kinematical structure $`\widehat{P}P^\mu P^\nu `$ (in what follows we will denote it $`ImT/\nu `$).

So, the sum rules for nucleon have the form

$$\frac{2\pi }{4M^4}\frac{\overline{\lambda }_N^2}{32\pi ^4}xq^{u,d}(x)e^{m^2/M^2}=ImT_{u,d}^0+\text{Power  corrections}$$
(14)
where $`\overline{\lambda }_N^2=32\pi ^4\lambda _N^2`$; $`q^{u,d}(x)`$ are distribution functions of $`u(d)`$ quark in nucleon, $`Im\stackrel{~}{T}^0`$ is perturbative contribution, i.e. of a bare loop with perturbative corrections. (The continuum contribution, i.e., of II, III, IV regions should be subtracted from $`Im\stackrel{~}{T}`$. Note, that really, the contribution of regions II and III to the bare loop is zero, since $`\rho ^0\delta (u_1u_2)`$).

3. Bare loop contribution and leading order perturbative corrections.

Bare loop contribution to the sum rules is represented in Fig.1. In the calculation we use the technique described in ref.15, Appendix. The following formulae are exploited:

$$\frac{d^4k}{(p_1k)^2(p_2k)^2}\delta [(p_1+qk)^2]\theta (k^2)=\frac{\pi }{4\nu x}(1x)𝑑u\frac{u}{(p_1^2u)(p_2^2u)}$$

$$\frac{d^4kk^2}{(p_1k)^2(p_2k)^2}\delta [(p_1+qk)^2]\theta (k^2)=\frac{\pi }{8\nu x}(1x)^2𝑑u\frac{u^2}{(p_1^2u)(p_2^2u)}$$

$$\frac{d^4k(p_1k)p_2}{(p_1k)^2(p_2k)^2}\delta [(p_1+qk)^2]\theta (k^2)=\frac{\pi }{16\nu x}(1x^2)𝑑u\frac{u^2}{(p_1^2u)(p_2^2u)}$$
(15)
(The terms, vanishing at the double Borel transformation are omitted). The results after the double Borel transformation are the same as in the case for equal $`p_1=p_2`$ :

$$ImT_{u(d)}^0=\phi _0^{u(d)}(x)\frac{M^2}{32\pi ^3}E_2(s_0/M^2)$$
(16)
where

$$\phi _0^u(x)=x(1x)^2(1+8x),\phi _0^d(x)=x(1x)^2(12x),$$
(17)
$`s_0`$ is the continuum threshold

$$E_2(z)=1(1+z+z^2/2)e^z$$
(18)
The substitution of eq.(16) into the sum rules (14) results in

$$xq(x)_0^{u(d)}=\frac{2M^6e^{m^2/M^2}}{\overline{\lambda }_N^2}\phi _0^{u(d)}(x)E_2(\frac{s_0}{M^2})$$
(19)
In the bare loop approximation the moments of the quark structure function are equal to

$$\underset{0}{\overset{1}{}}q_0^d(x)𝑑x=\frac{M^6e^{m^2/M^2}}{\overline{\lambda }_N^2}E_2$$

$$\underset{0}{\overset{1}{}}q_0^u(x)𝑑x=2\frac{M^6e^{m^2/M^2}}{\overline{\lambda }_N^2}E_2$$
(20)
Making use of relation $`\overline{\lambda }_N^2e^{m^2/M^2}=M^6E_2`$ which follows from the sum rule for the nucleon mass (see ) in the same approximation, we get

$$\underset{0}{\overset{1}{}}q_0^d(x)𝑑x=1$$

$$\underset{0}{\overset{1}{}}q_0^u(x)𝑑x=2$$
(21)
In the bare loop approximation there also appears the sum rule for the second moment:

$$\underset{0}{\overset{1}{}}x(q_0^u(x)+q_0^d(x))𝑑x=1$$
(22)
Analogously to  one can show that relations (21),(22) hold also when taking into account power corrections proportional to the quark condensate square in the sum rules for the 4-point correlator (Fig.2) and in the sum rules for the nucleon mass. Relations (21) reflect the fact that proton has two $`u`$-quarks and one $`d`$-quark. Relation (22) expresses the momentum conservation law – in the bare loop approximation all momentum is carried by valence quarks. Therefore, the sum rules (21),(22) demonstrate that the zero order approximation is reasonable. In the real physical theory the regions $`x1`$ and $`1x1`$ are off the frames of our consideration. However, in the noninteracting quark model which corresponds to the bare loop approximation, the whole region $`0x1`$ should be considered and relations (21),(22) should take place.

Let us calculate the perturbative corrections to bare loop and restrict ourselves by the leading order (LO) corrections proportional to $`lnQ_0^2/\mu ^2`$, where $`Q_0^2`$ is the point, where the quark distributions $`q(x,Q_0^2)`$ is calculated and $`\mu ^2`$ is the normalization point. In our case it is reasonable to choose $`\mu ^2`$ to be equal to the Borel parameter $`\mu ^2=M^2`$. The results take the form:

$$d^{LO}(x)=d_0(x)\{1+\frac{4}{3}ln(Q_0^2/M^2)\frac{\alpha _s(Q_0^2)}{2\pi }$$

$$[1/2+x+ln((1x)^2/x)+\frac{517x+16x^2+12x^3}{6(1x)(1+2x)}\frac{(32x)x^2ln(1/x)}{(1x)^2(1+2x)}]\}$$
(23)
$$u^{LO}(x)=u_0(x)\{1+\frac{4}{3}\frac{\alpha _s(Q_0^2)}{2\pi }ln(Q_0^2/M^2)[1/2+x+ln(1x)^2/x+$$

$$\frac{759x+46x^2+48x^3}{6(1x)(1+8x)}\frac{(158x)x^2ln(1/x)}{(1x)^2(1+8x)}]\}$$
(24)
where $`u_0(x)`$ and $`d_0(x)`$ are bare loop contributions, given by (19).

4.Power corrections

In this Section the power corrections to the sum rules are calculated.

The lowest dimension (d=4) power corrections are the corrections due to gluon condensate $`0|(\alpha _s/\pi )G_{\mu \nu }^nG_{\mu \nu }^n|0`$. Some examples of the corresponding diagrams are given in Fig.3. The calculations are performed in the fixed point (Fock-Schwinger) gauge. The corresponding expression for quark propagator is given in ref.12, eq.(24). (The error in the coefficient in front of the last term was corrected: it should be 1/288 instead of 1/96). In order to be sure in the final results the fixed point was chosen in two ways: at the upper and lower left-hand vertices of Fig.1 diagrams. The results for the sum of diagrams coincide in these two cases, as it should be. However, the contributions of separate diagrams are different. Particularly, if the fixed point is chosen at the upper left-hand vertex, then all diagrams where the soft gluon is emitted from the upper horizontal line (i.e., by active quark) in Fig.1, are zero in the lowest twist approximation. However, if the fixed point is at the lower left-hand vertex, the diagrams with soft gluon, emitted by active quark are nonzero. 111May be, this observation is related to the Brodsky et al.  statement, that structure functions are influent by active quark interactions with spectators. Therefore, it is generally untrue the folklore statement, that in the lowest twist approximation one may neglect in the forward scattering amplitude the active quark interaction between absorption and emission of virtual photon – this statement is gauge dependent.

All other technique is the same as in ref.15. The contributions of gluon condensate to $`u`$ and $`d`$-quarks distribution were found to be (the ratios to bare loop contributions are presented):

$$\frac{u(x)_{G^2}}{u_0(x)}=\frac{(\alpha _s/\pi )G^2}{M^4}\frac{\pi ^2}{12}\frac{(11+4x31x^2)}{x(1x)^2(1+8x)}(1e^{s_0/M^2})/E_2(\frac{s_0}{M^2})$$
(25)
$$\frac{d(x)_{G^2}}{d_0(x)}=\frac{(\alpha _s/\pi )G^2}{M^4}\frac{\pi ^2}{6}\frac{(12x^2)}{x^2(1x)^2(1+2x)}(1e^{s_0/M^2})/E_2(\frac{s_0}{M^2})$$
(26)
(The factor $`(1e^{s_0/M^2})`$ appears since gluonic condensate also contributes to continuum and this contribution should be subtracted).

Consider now the contributions of d=6 operators.

We start from studying the diagrams of Fig.4 and Fig.5, the contribution of which is proportional to the vacuum mean values of the type $`\overline{\psi }_\alpha \psi _\beta D_\rho G_{\mu \nu }^n`$, $`\overline{\psi }_\alpha _D\psi _\beta G_{\mu \nu }^n`$, $`\overline{\psi }\psi `$ ets. All of them are expressed through $`\overline{\psi }\psi ^2`$ (see ). Here it should be necessary to make the following remark: since the approach is inapplicable at $`x1`$, then the diagrams of the type of Fig.2, the imaginary part of which is proportional to $`\delta (1x)`$, should not to be taken into account. But the diagrams of Fig.4 arising by evolution of the diagram Fig.2, should not be taken into account too. The diagrams which contribute to intermediate region $`x`$ (and not being evolution of the corresponding non-loop ones) are given in Fig.5 (cf. Fig.5b in ). These diagrams are one-loop ones and for this reason it may be expected that their contribution will be dominant among contributions of d=6 operators. The contribution of Fig.5 diagram into sum rules for quark distribution in proton appeared to be equal to:

$$\frac{u(x)_{\alpha _s\overline{\psi }\psi ^2}}{u_0(x)}=\frac{\alpha _sa^2(215867x+172x^2+288(1x)ln2)}{M^681\pi 8x(1x)^3(1+8x)}\frac{1}{E_2(s_0/M^2)}$$
(27)
$$\frac{d(x)_{\alpha _s\overline{\psi }\psi ^2}}{d_0(x)}=\frac{\alpha _sa^2(1943x+36x^2)}{M^681\pi x(1x)^3(1+2x)}\frac{1}{E_2(s_0/M^2)},$$
(28)
where

$$a=(2\pi )^20|\overline{\psi }\psi |0$$
(29)
Besides the considered above contribution of the quark condensate square $`\alpha _s0|\overline{\psi }\psi |0^2`$, there exist contributions of d=6 operators described by two-loop diagrams. Examples of these diagrams are shown in Fig.6. The contributions of these diagrams are expressed via $`\alpha _s^20|\overline{\psi }\psi |0^2`$ and $`0|g^3f^{abc}G_{\mu \nu }^aG_{\nu \lambda }^bG_{\lambda \mu }^c|0`$. The diagrams of the first type are suppressed in comparison with (27), (28) by $`\alpha _s/\pi `$ and by numerical factor $`1/2\pi `$, because they are two loop diagrams. The diagrams of the second type, proportional to and $`0|g^3f^{abc}G_{\mu \nu }^aG_{\nu \lambda }^bG_{\lambda \mu }^c|0`$ are also suppressed by two-loop factor $`1/2\pi `$, but vacuum expectation value of $`G^3`$ operator is poorly known: there is only instanton estimate

$$g^3f^{abc}G^uG^bG^c=\frac{48\pi ^2}{5}(1/\rho _c^2)0|(\alpha _s/\pi )G^2|0$$
(30)
where $`\rho _c`$ is the effective radius of instanton. The magnitude of $`\rho _c`$ is not well known: various estimations result in essenially different values, varying from 1/3  to 1 fm , Also, one may have doubts  if instantons quantitatively describe physical processes at the scale $`1GeV^2`$, where we are working. Moreover, some of these diagrams are infrared divergent (for examples see Fig.6b).

For these reasons we restrict ourselves to the estimate of $`G^3`$ contribution only. Such estimate shows, that it may be essential in the domain of small $`x\begin{array}{c}<\hfill \\ \hfill \end{array}0.2`$ and is positive in both cases – for $`u`$\- and $`d`$-quarks.

The final result for valence quark distribution in proton are of the form

$$xu(x)=\frac{M^6e^{m^2/M^2}}{\overline{\lambda }_N^2}2x(1x)^2(1+8x)E_2(\frac{s_0}{M^2})\{[1+\frac{u^{LO}(x,Q_0^2)}{u_0(x)}]+$$

$$+\frac{1}{u_0(x)}[u(x)_{G^2}+u(x)_{\alpha _s\overline{\psi }\psi ^2}]\}$$
(31)
$$xd(x)=\frac{M^6e^{m^2/M^2}}{\overline{\lambda }_N^2}2x(1x)^2(1+2x)E_2(\frac{s_0}{M^2})\{[1+\frac{d^{LO}(x,Q_0^2)}{d_0(x)}]+$$

$$+\frac{1}{d_0(x)}[d(x)_{G^2}+d(x)_{\alpha _s\overline{\psi }\psi ^2}]\}$$
(32)
5. Valence $`u`$\- and $`d`$-quark distributions at intermediate $`x`$.

We are now in a position to calculate numerically valence $`u`$\- and $`d`$-quark distributions using eqs.(31),(32). The following values of the parameters were taken: $`\stackrel{~}{\lambda }_N^2=2.1GeV^6`$, $`s_0=2.3GeV^2`$ , $`\lambda _{QCD}=250MeV`$. The latter is the effective one-loop value, which gives the same $`\alpha _s(Q^2)`$ as many-loop calculations at low $`Q^2=15GeV^2`$ and particularly, $`\alpha _s(m_\tau ^2)=0.355`$ found in . The parameter $`\alpha _s0|\overline{\psi }\psi |o^2`$ was varied from $`8.10^5GeV^6`$ to $`\mathrm{2.2.10}^4GeV^6`$, i.e, $`\alpha _sa^2`$ varied from $`0.13GeV^6`$ to $`0.34GeV^6`$. (The lower limit is the old value , the upper limit corresponds to the recent determination of this parameter from $`\tau `$-decay data ,). The important parameter is the value of gluon condensate $`b=0|(\alpha _s/\pi )G_{\mu \nu }^2|0`$. We allow it to vary from $`0`$ up to $`0.012GeV^4`$. (The recent determination of $`0|(\alpha _s/\pi )G_{\mu \nu }^2|0`$ from the analysis of the $`\tau `$-decay data results in: $`0|(\alpha _s/\pi )G_{\mu \nu }^2|0=0.009\pm 0.007GeV^4`$ , the value $`0.012GeV^4`$ is the old result ). All results are presented at $`Q_0^2=5GeV^2`$.

Let us start with studying of the Borel mass dependence of quark distributions at various $`x`$. (Fig.7; the chosen values of the parameters are: $`a^2=0.34GeV^6`$, $`b=0.006GeV^4`$ – our favourite values – see below.) As is seen from Fig.7, for $`u`$-quark the stability in $`M^2`$ is good or satisfactory at $`0.1\begin{array}{c}<\hfill \\ \hfill \end{array}x\begin{array}{c}<\hfill \\ \hfill \end{array}0.5`$, at $`x=0.6`$ the $`M^2`$ dependence is rather strong. For $`d`$-quark distribution the stability interval is narrower; $`0.25\begin{array}{c}<\hfill \\ \hfill \end{array}x\begin{array}{c}<\hfill \\ \hfill \end{array}0.55`$.

Check now how much are the contributions of power corrections in comparison with $`d=0`$ OPE term (bare loop). For $`u`$-quark the applicability domain at low $`x`$ is limited by gluon condensate contribution – it comprises about $`30\%`$ at $`x=0.15`$ (and $`M^2=1.1GeV^2`$). At the same $`x=0.15`$ the contribution of $`G^3`$ term becomes large, although uncertain. At large $`x`$ the limits for $`u`$-quark come from perturbative corrections, they comprise about 40% at $`x=0.65`$. For $`d`$-quark the corresponding limits are $`0.25\begin{array}{c}<\hfill \\ \hfill \end{array}x\begin{array}{c}<\hfill \\ \hfill \end{array}0.60`$. (All the mentioned above values are given for $`\alpha _sa^2=0.34GeV^6`$, $`b=0.006GeV^4`$, $`M^2=1.1GeV^2`$).

Fig.8a,b shows the dependence of $`u(x)`$, $`d(x)`$ on the magnitude of gluon condensate. One can see from Fig.8, that the data following from the analysis of deep inelastic scattering are best described at $`b=0.006GeV^4`$. The zero value of gluon condensate cannot be excluded with certainty, but higher values, $`b\begin{array}{c}>\hfill \\ \hfill \end{array}0.012GeV^4`$ give much worse description of the data. The curves in Fig.8a,b were calculated at $`\alpha _sa^2=0.34GeV^6`$ (and $`M^2=1.1.GeV^2`$).The dependence of quark distributions on the magnitude of quark condensate square times $`\alpha _s`$$`\alpha _sa^2=(2\pi )^4\alpha _s|\overline{\psi }\psi |0^2`$ is plotted in Fig.9a,b (at $`b=0.006GeV^4`$). The best fit is obtained at $`\alpha _sa^2=0.34GeV^6`$. Therefore, as follows from the analysis, our favourite values of gluon condensate and $`\alpha _s`$ times quark condensate square are

$$0|\frac{\alpha _s}{\pi }G_{\mu \nu }^2|0=0.006GeV^4$$

$$\alpha _s0|\overline{\psi }\psi |0^2=2.2.10^4GeV^6$$
(33)
which agree with the recent result from the $`\tau `$-decay analysis , , .

Besides, the uncertainties arising from the spread of possible values of $`b`$ and $`a^2`$, there are additional sourses of errors – the uncertainty in $`\overline{\lambda }_N^2`$ and those coming from continuum contribution. The value of $`\overline{\lambda }_N^2=2.1GeV^6`$ was found in  by the best fit of the proton QCD sum rules. However, in  the old value of $`\alpha _s0|\psi \psi |0^2=8.10^5GeV^6`$ was used and $`\alpha _s`$ corrections were not accounted. This may result in $`1020\%`$ deviation of $`\overline{\lambda }_N^2`$ from the above accepted value. Continuum contribution comprises about $`30\%`$ of the total in $`u`$-quark case and about $`60\%`$ in $`d`$-quark case. So, one may estimate possible errors from uncertainty of continuum model as 10% and 20%, correspondingly, for $`u`$\- and $`d`$\- quarks. Therefore, our estimates of total errors are: for $`u`$-quarks distributions, in the middle of $`x`$-interval – $`0.25<x<0.45`$, about 25%, at the ends of intervals – $`x=0.15`$ and $`0.65`$ – about 50%; for $`d`$ quarks distributions, in the middle of $`x`$-interval – $`0.3<x<0.45`$ about 30%, at the ends of intervals – $`x=0.25`$ and $`0.55`$ by a factor of 2. Taking into account these errors, the agreement with $`u`$\- and $`d`$-quark distributions, found in  from experimental data (solid curve in Fig.8) is satisfactory. It should be mentioned, that the account of $`G^3`$ terms would improve the agreement in the domain of low $`x`$.

A few additional remarks are in order. In Figs.8,9 we performed the comparison with valence quark distributions found from hard processes, we used only one group results . We limited our comparison to these data not because we consider these results as more confident or reliable, than those of other groups \[5-9\]. Quite opposite, we believe, that the precision of new analysis (see especially ,) is better than the older ones. The reason is that theoretical errors in determination of $`u`$\- and $`d`$\- quark distribution exceed the differences in various treatment of the data – all of them are in the limit of our theoretical uncertainties. We chosed for comparison the LO results of , because our perturbative calculations are done also in LO. Our results can be compared directly with experimental data by considering the difference

$$F_2^p(x)F_2^n(x)=\frac{1}{3}[u(x)d(x)]_{val}+\frac{1}{3}[\overline{u}(x)\overline{d}(x)]$$
(34)
$`F_2^p(x)F_2^n(x)`$ was measured by NMC . At $`x>0.30`$ one may expect, that sea quark contribution – the last term in (36) – is negligible and the data can be compared with valence quark distributions. The result is that at $`0.3<x<0.55`$ the theoretical curve (at $`b=0.006`$ GeV6) is about 30-40% higher than the data points. Such disagreement is in the limit of estimated errors. (At $`b=0`$ the agreement is better.)

6. Conclusion

The distributions of valence $`u`$\- and $`d`$-quarks at low $`Q^2`$ and intermediate $`x`$ were theoretically calculated basing on the first principles of QCD: perturbation theory and operator product expansion (OPE). No experimental data and no fitting parameters were used. New technique, suggested in ref.15, was exploited: the double Borel transformation in virtualities of two hadronic currents in 4-point function eq.(1), which allows one to kill the background non-diagonal transition amplitudes. In OPE the operators of dimension 4 (gluon condensate) and dimension 6 (quark condensate square times $`\alpha _s`$) were accounted. The theoretical analysis of the obtained valence quark distributions at $`Q^2=5GeV^2`$ showed, that $`u`$-quark distribution is reliable at $`0.15<x<0.65`$, its accuracy is about 25% in the middle of this interval and decreases to 50% at the ends of interval; $`d`$-quark distribution is reliable at $`0.25<x<0.55`$ with an accuracy of about 30% in the middle and is given by a factor of 2 at the ends of interval. In the limit of these accuracies the theoretically calculated valence quark distributions are in agreement with those found from deep inelastic scattering and other hard processes data. The account of $`G^3`$ contribution, probably, improves the agreement, especially at small $`x`$.

Acknowledgement

The authors are thankful to Wu-Ki Tung for his kind presentation of the CTEQ analysis data. The research described in this publication was made possible in part by Award No.RP2-2247 of the US Civilian and Development Foundation for the Independent States of Former Soviet Union (CRDF), by INTAS Grant 2000, Project 587 and by the Russian Found of Basic Research, Grant No. 00-02-17808.


text:
# I Introduction

## I Introduction

A compact object of solar-mass size orbiting a supermassive black hole is one of promising candidates for the source of gravitational waves. Since the internal structure of such a compact object may be neglected in this situation, we may adopt the black hole perturbation approach with the compact object being regarded as a point particle. In the black hole perturbation approach, we consider the metric perturbation induced by a point particle of mass $`\mu `$ orbiting a black hole of mass $`M`$, where $`\mu M`$. At the lowest-order in the mass ratio ($`\mu /M`$), the motion of the particle follows a geodesic of the background spacetime. In the next order, however, the particle moves no longer along a geodesic of the background because of its interaction with the self-field.

Although this deviation from a background geodesic is small for $`\mu /M1`$ at each instant of time, after a large lapse of time, it accumulates to become non-negligible. For example, a circular orbit will not remain circular but becomes a spiral-in orbit and the orbit eventually plunges into the black hole.

If the time scale of the orbital evolution due to the self-force is sufficiently long compared to the characteristic orbital time, we may adopt the so-called adiabatic approximation in which the orbit is assumed to be instantaneously geodesic with the constants of motion changing very slowly with time. In the Schwarzschild background case, we may assume the orbit to lie on the equatorial plane and the geodesic motion is determined by the energy and (the $`z`$-component of) angular momentum of the particle. In this case, the time variation of the energy and angular momentum can be determined from the energy and angular momentum emitted to infinity and absorbed into the black hole horizon by using the conservation law.

However, there are cases when the adiabatic approximation breaks down. For example, in the case of an extremely eccentric orbit or an orbit close to the inner-most stable circular orbit, the orbital evolution will not be adiabatic because the stability of the orbit is strongly affected by an infinitesimally small reaction force. Furthermore, in the Kerr background, there is an additional constant of motion, known as the Carter constant. Intuitively, it describes the total orbital angular momentum, but unlike the case of spherical symmetry, it has nothing to do with the Killing vector field of the Kerr geometry. The lack of its relation to the Killing vector makes us impossible to evaluate the time change of the Carter constant from the gravitational waves emitted to infinity and to event horizon, even in the case when the adiabatic approximation is valid. Thus it is in any case necessary to derive the self-force of a particle explicitly.

The gravitational self-force $`F^\mu `$ is formally given as

$$\frac{d^2z^\alpha }{d\tau ^2}+\mathrm{\Gamma }_{\mu \nu }^\alpha \frac{dz^\mu }{d\tau }\frac{dz^\nu }{d\tau }=F^\alpha ,$$

where $`\{z^\alpha (\tau )\}`$ represents the orbit with $`\tau `$ being the proper time measured in the background geometry and $`\mathrm{\Gamma }_{\mu \nu }^\alpha `$ is the connection of the background. The self-force arizes from the metric perturbation $`h_{\mu \nu }`$ induced by the particle:

$$\stackrel{~}{g}_{\mu \nu }=g_{\mu \nu }+h_{\mu \nu },$$

and it is expressed as

$$F^\alpha [h]=\mu P_\beta ^\alpha (\overline{h}_{\beta \gamma ;\delta }\frac{1}{2}g_{\beta \gamma }\overline{h}^ϵ{}_{ϵ;\delta }{}^{}\frac{1}{2}\overline{h}_{\gamma \delta ;\beta }+\frac{1}{4}g_{\gamma \delta }\overline{h}^ϵ{}_{ϵ;\beta }{}^{})u^\gamma u^\delta ,$$

where $`P_{\alpha }^{}{}_{}{}^{\beta }=\delta _{\alpha }^{}{}_{}{}^{\beta }+u_\alpha u^\beta `$, $`\overline{h}_{\alpha \beta }=h_{\alpha \beta }\frac{1}{2}g_{\alpha \beta }h`$ and $`u^\alpha =dz^\alpha /d\tau `$.

The metric perturbation diverges at the location of the particle and so does the self-force. Thus the above formal expression is in fact meaningless. Fortunately, however, it is known that the metric perturbation in the vicinity of the orbit can be divided into two parts under the harmonic gauge condition; the direct part which has support only on the past light-cone emanating from the field point $`x^\mu `$ and the tail part which has support inside the past light-cone, and the physical self-force is given by the tail part of the metric perturbation which is regular as we let the field point coincide with a point on the orbit; $`x^\mu z^\mu (\tau )`$ Mino:1996nk ; Quinn:1999kj . It must be noted that the direct part can be evaluated by local analysis, i.e., only with the knowledge of local geometrical quantities. Therefore the physical self-force can be calculated as

$$\underset{xz(\tau )}{lim}F_\alpha [h^{\mathrm{tail}}(x)]=\underset{xz(\tau )}{lim}\left(F_\alpha [h(x)]F_\alpha [h^{\mathrm{dir}}(x)]\right).$$

Furthermore, it has been revealed recently by Detweiler and Whiting Detweiler:2002mi  that the above devision of the metric can be slightly modified so that the new direct part, called the S part, satisfies the same Einstein equations as the full metric perturbation does, and the new tail part, called the R part, satisfies the source-free Einstein equations, and that the R part gives the identical, regular self-force as the tail part does. The important point is that the S part can be still evaluated locally near the orbit without knowing the global solution.

When we perform this subtraction, we must evaluate the full self-force and the direct part under the same gauge condition. But the direct part is, by definition, defined only in the harmonic gauge. On the other hand, the full metric perturbation is directly obtainable only by the Regge-Wheeler-Zerilli or Teukolsky formalism Regge:1957td ; Zerilli:wd ; Teukolsky:1973ha ; Chrzanowski:wv . Therefore one must find a gauge transformation that brings both the full metric perturbation and the direct part of the metric perturbation to those in the same gauge. This is called the gauge problem 111Furthermore, the gravitational self-force is, because of the equivalence principle, a gauge-variant notion. To give a genuinely physical meaning to it, one must solve the second order metric perturbation completely. This is however beyond the scope of the present paper..

For the direct part, methods to obtain it under the harmonic gauge condition were proposed Mode-sum ; Mino:2001mq ; Barack:2001gx ; Barack:2002mh . However, it seems extremely difficult to solve the metric perturbation under the harmonic gauge because the metric components couple to each other in a complicated way. This is one of the reasons why the gauge problem is difficult to solve.

Recently, Barack and Ori Barack:2001ph  gave a useful insight into the gauge problem. They proposed an intermediate gauge approach in which only the direct part of the metric in the harmonic gauge is subtracted from the full metric perturbation in the RW gauge. They then argued that the gauge-dependence of the self-force is unimportant when averaged over a sufficiently long lapse of time. Using this approach, the gravitational self-force for an orbit plunging into a Schwarzschild black hole was calculated by Barack and Lousto Barack:2002ku . But they also pointed out that the RW gauge is singular in the sense that the resulting self-force will still have a direction-dependent limit for general orbits. The situation becomes worse in the Kerr background where the only known gauge in which the metric perturbation can be evaluated is the radiation gauge Chrzanowski:wv , but the metric perturbation becomes ill-defined in the neighborhood of the particle, i.e., the Einstein equations are not satisfied there Barack:2001ph .

In this paper, as a direct approach to the gauge problem, we consider a formalism to calculate the metric perturbation in the harmonic gauge. We focus on the Schwarzschild background. Instead of directly solving the metric perturbation in the harmonic gauge, we consider the metric perturbation in the RW gauge first, and then transform it to the one in the harmonic gauge. Namely, we derive a set of equations for gauge functions that transform the metric perturbation in the RW gauge to the one in the harmonic gauge.

The paper is organized as follows. In Sec. II, we formulate the gauge transformation from the RW gauge to the harmonic gauge. First we decompose the gauge transformation generators into the Fourier-harmonics components. Then the generators are devided into three parts; the odd parity part and the even parity part which is further devided into scalar and divergence free parts. By the above procedure, we find a set of decoupled equations for the gauge functions. In Sec. III, we summarize our formulation and discuss remaining issues. In Appendix A, we recapitulate the equations for the Regge-Wheeler-Zerilli formalism by correcting typos in Zerilli’s paper Zerilli:wd .

## II Formulation

We consider a metric perturbation $`h_{\mu \nu }`$ in the Schwarzschild background,

$`g_{\mu \nu }dx^\mu dx^\nu =f(r)dt^2+f(r)^1dr^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2);f(r)=1{\displaystyle \frac{2M}{r}}.`$ (1)
We express the gauge transformation from the RW gauge to the harmonic gauge as

$`x_{\mathrm{RW}}^\mu `$ $``$ $`x_\mathrm{H}^\mu =x_{\mathrm{RW}}^\mu +\xi ^\mu ,`$ (2)
$`h_{\mu \nu }^{\mathrm{RW}}`$ $``$ $`h_{\mu \nu }^\mathrm{H}=h_{\mu \nu }^{\mathrm{RW}}\xi _{\mu ;\nu }\xi _{\nu ;\mu },`$ (3)
where the suffix RW stands for the RW gauge and H for the harmonic gauge.

Substituting Eq. (3) into the harmonic gauge condition $`\overline{h}_{\mu \nu }^\mathrm{H}{}_{}{}^{;\nu }=0`$, we obtain the equations for $`\xi ^\mu `$:

$`\xi _\mu {}_{}{}^{;\nu }{}_{;\nu }{}^{}=\overline{h}_\mu ^{\mathrm{RW}}{}_{}{}^{\nu }{}_{;\nu }{}^{}=h_\mu ^{\mathrm{RW}}{}_{}{}^{\nu }{}_{;\nu }{}^{}{\displaystyle \frac{1}{2}}h_{;\mu }^{\mathrm{RW}}.`$ (4)
Using the static and spherical symmetry of the background, we perform the Fourier-harmonic expansion of the above equation and consider the equations for the expansion coefficients for $`\xi ^\mu `$ and $`h_{\mu \nu }^{\mathrm{RW}}`$. We use the tensor harmonics introduced by Zerilli Zerilli:wd  which are recapitulated in Appendix A. Then according to the property under the parity transformation $`(\theta ,\varphi )(\pi \theta ,\varphi +\pi )`$, we devide $`\xi ^\mu `$ into the odd and even parity parts,

$`\xi ^\mu =\xi _{(\mathrm{odd})}^\mu +\xi _{(\mathrm{even})}^\mu ,`$ (5)
and the even part is further decomposed into the scalar and divergence-free part,

$`\xi _{(\mathrm{even})}^\mu =\xi ^{;\mu }+\xi _{(\mathrm{v})}^\mu ,`$ (6)
where $`\xi _{(\mathrm{v});\mu }^\mu =0`$.

### II.1 Odd part

First, we consider the odd part which has the odd parity $`(1)^{\mathrm{}+1}`$ under the parity transformation. The gauge transformation generators and the metric perturbation are given in the Fourier-harmonic expanded form as

$`\xi _\mu ^{(\mathrm{odd})}`$ $`=`$ $`{\displaystyle 𝑑\omega \underset{\mathrm{}m}{}e^{i\omega t}\mathrm{\Lambda }_{\mathrm{}m\omega }(r)\{0,0,\frac{1}{\mathrm{sin}\theta }_\varphi Y_\mathrm{}m(\theta ,\varphi ),\mathrm{sin}\theta _\theta Y_\mathrm{}m(\theta ,\varphi )\}},`$ (7)
$`h_{\mu \nu }^{(\mathrm{odd})}`$ $`=`$ $`{\displaystyle 𝑑\omega \underset{\mathrm{}m}{}\frac{\sqrt{2\mathrm{}(\mathrm{}+1)}}{r}e^{i\omega t}}`$ (8)
$`\times \left[h_{0\mathrm{}m\omega }(r)𝒄_{\mathrm{}m\mu \nu }^{(0)}+ih_{1\mathrm{}m\omega }(r)𝒄_{\mathrm{}m\mu \nu }+{\displaystyle \frac{\sqrt{(\mathrm{}1)(\mathrm{}+2)}}{2r}}h_{2\mathrm{}m\omega }(r)𝒅_{\mathrm{}m\mu \nu }\right],`$
where $`𝒄_{\mathrm{}m\mu \nu }^{(0)}`$,$`𝒄_{\mathrm{}m\mu \nu }`$, $`𝒅_{\mathrm{}m\mu \nu }`$ are tensor harmonics with odd parity Zerilli:wd . By substituting Eqs. (7) and (8) into Eq. (4), we obtain the equation for $`\mathrm{\Lambda }_{\mathrm{}m\omega }(r)`$:

$`^{(\mathrm{odd})}\mathrm{\Lambda }_{\mathrm{}m\omega }(r)=2R_{\mathrm{}m\omega }^{(\mathrm{odd})}(r){\displaystyle \frac{16i\pi r^2}{\sqrt{2\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)}}}D_{\mathrm{}m\omega }(r),`$ (9)
where $`R_{\mathrm{}m\omega }^{(\mathrm{odd})}(r)`$ is the Regge-Wheeler gauge invariant variable, $`D_{\mathrm{}m\omega }(r)`$ is the Fourier-harmonic coefficient of the stress-energy tensor, and the differential operator $``$ is defined as

$`^{(\mathrm{odd})}`$ $``$ $`f(r){\displaystyle \frac{d^2}{dr^2}}+f^{}(r){\displaystyle \frac{d}{dr}}+\left({\displaystyle \frac{\omega ^2}{f(r)}}{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r^2}}\right)`$ (10)
$`=`$ $`{\displaystyle \frac{1}{f(r)}}{\displaystyle \frac{d^2}{dr^{\mathrm{\hspace{0.17em}2}}}}+\left({\displaystyle \frac{\omega ^2}{f(r)}}{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r^2}}\right),`$
where $`r^{}=r+2M\mathrm{log}(r/2M1)`$ and $`{}_{}{}^{}=d/dr`$. The form of the differential operator $``$ is slightly different from the radial part of the scalar d’Alembertian, but Eq. (9) may be solved by the standard Green function method. The homogeneous solutions to construct the Green function can be obtained by applying the method developed by Mano et al. Mano1 . Here we note that the retarded causal boundary condition should be imposed for the Green function.

### II.2 Even scalar part

Next we consider the even scalar part of $`\xi ^\mu `$ which has the even parity $`(1)^{\mathrm{}}`$ and is expressed by a gradient of a scalar function $`\xi `$. It is expressed as

$`\xi _\mu ^{(\mathrm{s})}=\xi _{;\mu };\xi ={\displaystyle 𝑑\omega \underset{\mathrm{}m}{}\frac{1}{r}\stackrel{~}{\xi }_{\mathrm{}m\omega }(r)e^{i\omega t}Y_\mathrm{}m(\theta ,\varphi )}.`$ (11)
The gauge equation for $`\xi _\mu ^{(\mathrm{s})}`$ is derived from $`\xi _{\mu ;\nu }^{(\mathrm{s})}{}_{}{}^{;\nu }=J_{;\mu }^{(\mathrm{s})}`$, which, because of the vanishing of the background Ricci tensor, gives

$`\xi ^{;\nu }{}_{;\nu }{}^{}=J^{(\mathrm{s})},`$ (12)
where the source term $`J^{(\mathrm{s})}`$ is determined from the equation,

$`J_{}^{(\mathrm{s});\mu }{}_{;\mu }{}^{}`$ $`=`$ $`\overline{h}_{\mathrm{RW};\mu \nu }^{\mu \nu }.`$ (13)
The Fourier-harmonic expanded form of the above equation is derived as follows. First, we take the divergence of the metric perturbation under the RW gauge;

$`\overline{h}_{\mu \nu }^{\mathrm{RW};\nu }`$ $`=`$ $`{\displaystyle 𝑑\omega \underset{\mathrm{}m}{}e^{i\omega t}}`$ (14)
$`\times \{[{\displaystyle \frac{2f(r)}{r}}H_{1\mathrm{}m\omega }^{\mathrm{RW}}(r){\displaystyle \frac{8\pi irB_{\mathrm{}m\omega }^{(0)}(r)}{\sqrt{\mathrm{}(\mathrm{}+1)/2}}}{\displaystyle \frac{8\pi i\omega r^2F_{\mathrm{}m\omega }(r)}{\sqrt{\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)/2}}}]Y_\mathrm{}m(\theta ,\varphi )(𝒆_t)_\mu `$
$`+\left[{\displaystyle \frac{2}{r}}(H_{2\mathrm{}m\omega }^{\mathrm{RW}}(r)K_{\mathrm{}m\omega }^{\mathrm{RW}}(r))+{\displaystyle \frac{8\pi rB_{\mathrm{}m\omega }(r)}{\sqrt{\mathrm{}(\mathrm{}+1)/2}}}{\displaystyle \frac{8\pi r^2\frac{d}{dr}F_{\mathrm{}m\omega }(r)}{\sqrt{\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)/2}}}\right]Y_\mathrm{}m(\theta ,\varphi )(𝒆_r)_\mu `$
$`+{\displaystyle \frac{8\pi r^2}{\sqrt{\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)/2}}}F_{\mathrm{}m\omega }(r)(𝒆_3)_\mu \},`$
where $`(𝒆_t)_\mu =(1,0,0,0)`$, $`(𝒆_r)_\mu =(0,1,0,0)`$ and $`(𝒆_3)_\mu =(0,0,_\theta Y_\mathrm{}m,_\varphi Y_\mathrm{}m)`$. Then we take the divergence of the above once more to obtain

$`\overline{h}_{\mu \nu }^{\mathrm{RW};\mu \nu }`$ $`=`$ $`{\displaystyle 𝑑\omega \underset{\mathrm{}m}{}e^{i\omega t}Y_\mathrm{}m(\theta ,\varphi )}`$ (15)
$`\times \left[^{(\mathrm{s})}\stackrel{~}{H}_{\mathrm{}m\omega }(r)+8\pi \left({\displaystyle \frac{A_{\mathrm{}m\omega }^{(0)}(r)}{f(r)}}f(r)A_{\mathrm{}m\omega }(r)\sqrt{2}G_{\mathrm{}m\omega }^{(\mathrm{s})}\right)\right]`$
where $`\stackrel{~}{H}_{\mathrm{}m\omega }(r)`$ is the trace of the metric perturbation given by

$`\stackrel{~}{H}_{\mathrm{}m\omega }(r)={\displaystyle \frac{H_{0\mathrm{}m\omega }^{\mathrm{RW}}(r)H_{2\mathrm{}m\omega }^{\mathrm{RW}}(r)}{2}}K_{\mathrm{}m\omega }^{\mathrm{RW}}(r).`$ (16)
The d’Alembertian of $`J^{(\mathrm{s})}`$ is expanded as

$`J_{}^{(\mathrm{s});\mu }{}_{;\mu }{}^{}={\displaystyle 𝑑\omega \underset{\mathrm{}m}{}e^{i\omega t}Y_\mathrm{}m(\theta ,\varphi )^{(\mathrm{s})}\stackrel{~}{J}_{\mathrm{}m\omega }^{(\mathrm{s})}(r)},`$ (17)
where $`\stackrel{~}{J}_{\mathrm{}m\omega }^{(\mathrm{s})}`$ is the Fourier-harmonic coefficient of $`J^{(\mathrm{s})}`$ and $`^{(\mathrm{s})}`$ is the radial part of the scalar d’Alembertian,

$`^{(\mathrm{s})}`$ $``$ $`f(r){\displaystyle \frac{d^2}{dr^2}}+f^{}(r){\displaystyle \frac{d}{dr}}+\left({\displaystyle \frac{\omega ^2}{f(r)}}{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r^2}}{\displaystyle \frac{f^{}(r)}{r}}\right)`$ (18)
$`=`$ $`{\displaystyle \frac{1}{f(r)}}{\displaystyle \frac{d^2}{dr^{\mathrm{\hspace{0.17em}2}}}}+\left({\displaystyle \frac{\omega ^2}{f(r)}}{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r^2}}{\displaystyle \frac{2M}{r^3}}\right).`$
Hence Eq. (13) becomes

$`^{(\mathrm{s})}\stackrel{~}{J}_{\mathrm{}m\omega }^{(\mathrm{s})}(r)=^{(\mathrm{s})}\stackrel{~}{H}_{\mathrm{}m\omega }(r)+8\pi \left({\displaystyle \frac{A_{\mathrm{}m\omega }^{(0)}(r)}{f(r)}}f(r)A_{\mathrm{}m\omega }(r)\sqrt{2}G_{\mathrm{}m\omega }^{(\mathrm{s})}\right).`$ (19)
Once the source term is obtained, we obtain from Eq. (12) the equation for $`\stackrel{~}{\xi }_{\mathrm{}m\omega }(r)`$ as

$`^{(\mathrm{s})}\stackrel{~}{\xi }_{\mathrm{}m\omega }(r)`$ $`=`$ $`r\stackrel{~}{J}_{\mathrm{}m\omega }^{(\mathrm{s})}(r).`$ (20)
Note that this equation as well as Eq. (19) are just the radial part of the scalar d’Alembertian, or equivalently the spin $`s=0`$ Teukolsky equation, hence can be also solved by the Green function method. The detail analysis of the homogeneous solutions which satisfy the $`s=0`$ Teukolsky equation is discussed in Mano1 .

### II.3 Even vector part

Since the even vector part $`\xi _{(\mathrm{v})}^\mu `$ satisfies the divergence free condition, $`\xi _\mu ^{(\mathrm{v})}{}_{}{}^{;\mu }=0`$, this part has two degrees of freedom. Therefore $`\xi _{(\mathrm{v})}^\mu `$ is expressed in terms of two independent radial functions,

$`\xi _\mu ^{(\mathrm{v})}={\displaystyle 𝑑\omega \underset{\mathrm{}m}{}e^{i\omega t}}`$ $`[`$ $`{\displaystyle \frac{1}{r}}M_{0\mathrm{}m\omega }(r)Y_\mathrm{}m(\theta ,\varphi )(𝒆_t)_\mu +{\displaystyle \frac{1}{r}}M_{1\mathrm{}m\omega }(r)Y_\mathrm{}m(\theta ,\varphi )(𝒆_r)_\mu `$ (21)
$`+{\displaystyle \frac{1}{\mathrm{}(\mathrm{}+1)}}\{{\displaystyle \frac{i\omega r}{f(r)}}M_{0\mathrm{}m\omega }(r)+_r\left(rf(r)M_{1\mathrm{}m\omega }(r)\right)\}(𝒆_3)_\mu ].`$
To solve the gauge transformation equation (4) for this part, we introduce an auxiliary field $`F_{\mu \nu }:=\xi _{\nu ;\mu }^{(\mathrm{v})}\xi _{\mu ;\nu }^{(\mathrm{v})}`$. Then we have

$`F_{\mu \nu }{}_{}{}^{;\nu }=J_\mu ^{(\mathrm{v})};J_\mu ^{(\mathrm{v})}\overline{h}_{\mu \nu }^{\mathrm{RW};\nu }J_{;\mu }^{(\mathrm{s})},`$ (22)
where we have again used the fact that the background Ricci tensor vanishes. Note that the Fourier-harmonic expansion of $`\overline{h}_{\mu \nu }^{\mathrm{RW};\nu }`$ is given by Eq. (14). The above equation is the Maxwell equation, so it can be solved by using the Teukolsky formalism Teukolsky:1973ha  for the electromagnetic field ($`s=\pm 1`$). Namely, we introduce the following variables:

$`\varphi _0=F_{\mu \nu }l^\mu m^\nu ,\varphi _2=F_{\mu \nu }m^\mu n^\nu ,`$ (23)
where $``$ denotes the complex conjugate, and $`l^\mu `$, $`n^\mu `$, $`m^\mu `$ are the Kinnersley null tetrad defined by

$`l^\mu =\{{\displaystyle \frac{1}{f(r)}},1,0,0\},n^\mu =\{{\displaystyle \frac{1}{2}},{\displaystyle \frac{f(r)}{2}},0,0\},m^\mu ={\displaystyle \frac{1}{\sqrt{2}r}}\{0,0,1,{\displaystyle \frac{i}{\mathrm{sin}\theta }}\}.`$ (24)
The variables $`\varphi _0`$ and $`\varphi _2`$ satisfy the $`s=\pm 1`$ Teukolsky equation,

$`_s^{(\mathrm{Teuk})}\mathrm{\Psi }_s=4\pi r^2T_s,`$ (25)
where $`\mathrm{\Psi }_1=\varphi _0`$ and $`\mathrm{\Psi }_1=\varphi _2/r^2`$, and $`_s^{(\mathrm{Teuk})}`$ is the Teukolsky differential operator defined as

$`_s^{(\mathrm{Teuk})}`$ $`:=`$ $`{\displaystyle \frac{r^2}{f(r)}}_t^2+2s\left({\displaystyle \frac{M}{f(r)}}r\right)_t+\left(r^2f(r)\right)^s_r\left[\left(r^2f(r)\right)^{s+1}_r\right]`$ (26)
$`+{\displaystyle \frac{1}{\mathrm{sin}\theta }}_\theta (\mathrm{sin}\theta _\theta )+{\displaystyle \frac{1}{\mathrm{sin}^2\theta }}_\varphi ^2+{\displaystyle \frac{2is\mathrm{cos}\theta }{\mathrm{sin}^2\theta }}_\varphi s(s\mathrm{cot}^2\theta 1).`$
The source terms $`T_s`$, which are calculated from $`J_\mu ^{(\mathrm{v})}`$, are

$`T_1`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}r}}\left(_\theta +{\displaystyle \frac{i}{\mathrm{sin}\theta }}_\varphi \right)\left[J_\mu ^{(\mathrm{v})}l^\mu \right]\left(_r{\displaystyle \frac{i\omega }{f(r)}}+{\displaystyle \frac{3}{r}}\right)\left[J_\mu ^{(\mathrm{v})}m^\mu \right],`$ (27)
$`T_1`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}r}}\left(_\theta {\displaystyle \frac{i}{\mathrm{sin}\theta }}_\varphi \right)\left[J_\mu ^{(\mathrm{v})}n^\mu \right]{\displaystyle \frac{f(r)}{2}}\left(_r+{\displaystyle \frac{i\omega }{f(r)}}+{\displaystyle \frac{3}{r}}\right)\left[J_\mu ^{(\mathrm{v})}m^\mu \right].`$ (28)
The solution of this Teukolsky equation is also analyzed in Mano1 .

Once we obtain $`\varphi _0`$ and $`\varphi _2`$, the two radial functions for the gauge transformation are calculated from Eq. (23) which reads

$`\varphi _0`$ $`=`$ $`{\displaystyle }d\omega {\displaystyle \underset{\mathrm{}m}{}}e^{i\omega t}{\displaystyle \frac{{}_{1}{}^{}Y_{\mathrm{}m}^{}(\theta ,\varphi )}{\sqrt{2\mathrm{}(\mathrm{}+1)}rf(r)}}[{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r}}(M_{0\mathrm{}m\omega }(r)f(r)M_{1\mathrm{}m\omega }(r))`$ (29)
$`i\omega \{{\displaystyle \frac{i\omega r}{f(r)}}M_{0\mathrm{}m\omega }(r)+{\displaystyle \frac{d}{dr}}(rf(r)M_{1\mathrm{}m\omega }(r))\}+f(r){\displaystyle \frac{d}{dr}}\{{\displaystyle \frac{i\omega r}{f(r)}}M_{0\mathrm{}m\omega }(r)+{\displaystyle \frac{d}{dr}}(rf(r)M_{1\mathrm{}m\omega }(r))\}],`$
$`\varphi _2`$ $`=`$ $`{\displaystyle }d\omega {\displaystyle \underset{\mathrm{}m}{}}e^{i\omega t}{\displaystyle \frac{{}_{1}{}^{}Y_{\mathrm{}m}^{}(\theta ,\varphi )}{2\sqrt{2\mathrm{}(\mathrm{}+1)}r}}[{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r}}(M_{0\mathrm{}m\omega }(r)+f(r)M_{1\mathrm{}m\omega }(r))`$ (30)
$`+i\omega \{{\displaystyle \frac{i\omega r}{f(r)}}M_{0\mathrm{}m\omega }(r)+{\displaystyle \frac{d}{dr}}(rf(r)M_{1\mathrm{}m\omega }(r))\}+f(r){\displaystyle \frac{d}{dr}}\{{\displaystyle \frac{i\omega r}{f(r)}}M_{0\mathrm{}m\omega }(r)+{\displaystyle \frac{d}{dr}}(rf(r)M_{1\mathrm{}m\omega }(r))\}],`$
where $`{}_{s}{}^{}Y_{\mathrm{}m}^{}(\theta ,\varphi )`$ are the spin-weighted spherical harmonics. Performing the Fourier and spin-weighted spherical harmonic expansion for $`\varphi _0`$ and $`\varphi _2`$,

$`\varphi _0`$ $`=`$ $`{\displaystyle 𝑑\omega \underset{\mathrm{}m}{}\stackrel{~}{\varphi }_{0\mathrm{}m\omega }(r)e^{i\omega t}{}_{1}{}^{}Y_{\mathrm{}m}^{}(\theta ,\varphi )},`$ (31)
$`\varphi _2`$ $`=`$ $`{\displaystyle 𝑑\omega \underset{\mathrm{}m}{}\stackrel{~}{\varphi }_{2\mathrm{}m\omega }(r)e^{i\omega t}{}_{1}{}^{}Y_{\mathrm{}m}^{}(\theta ,\varphi )},`$ (32)
we obtain the equations,

$`\stackrel{~}{\varphi }_{0\mathrm{}m\omega }(r)+{\displaystyle \frac{2}{f(r)}}\stackrel{~}{\varphi }_{2\mathrm{}m\omega }(r)`$ $`=`$ $`{\displaystyle \frac{\sqrt{2/\mathrm{}(\mathrm{}+1)}}{rf(r)}}\left[\left({\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r}}{\displaystyle \frac{\omega ^2r}{f(r)}}\right)M_{0\mathrm{}m\omega }(r)i\omega {\displaystyle \frac{d}{dr}}(rf(r)M_{1\mathrm{}m\omega }(r))\right],`$ (33)
$`\stackrel{~}{\varphi }_{0\mathrm{}m\omega }(r){\displaystyle \frac{2}{f(r)}}\stackrel{~}{\varphi }_{2\mathrm{}m\omega }(r)`$ $`=`$ $`{\displaystyle \frac{\sqrt{2/\mathrm{}(\mathrm{}+1)}}{rf(r)}}[{\displaystyle \frac{\mathrm{}(\mathrm{}+1)f(r)}{r}}M_{1\mathrm{}m\omega }(r)`$ (34)
$`f(r){\displaystyle \frac{d^2}{dr^2}}(rf(r)M_{1\mathrm{}m\omega }(r))+i\omega f(r){\displaystyle \frac{d}{dr}}\left({\displaystyle \frac{r}{f(r)}}M_{0\mathrm{}m\omega }(r)\right)].`$
We can eliminate $`M_{1\mathrm{}m\omega }(r)`$ from Eqs. (33) and (34) to obtain a decoupled equation for $`M_{0\mathrm{}m\omega }(r)`$,

$`^{(\mathrm{s})}M_{0\mathrm{}m\omega }(r)`$ $`=`$ $`{\displaystyle \frac{r^2}{\sqrt{2\mathrm{}(\mathrm{}+1)}}}[f(r)^2{\displaystyle \frac{d^2}{dr^2}}\stackrel{~}{\varphi }_{0\mathrm{}m\omega }(r)+f(r)({\displaystyle \frac{2+2f(r)}{r}}+f^{}(r)+i\omega ){\displaystyle \frac{d}{dr}}\stackrel{~}{\varphi }_{0\mathrm{}m\omega }(r)`$ (36)
$`+{\displaystyle \frac{rf^{}(r)(\mathrm{}1)(\mathrm{}+2)f(r)+i\omega r(1+2f(r))}{r^2}}\stackrel{~}{\varphi }_{0\mathrm{}m\omega }(r)+2f(r){\displaystyle \frac{d^2}{dr^2}}\stackrel{~}{\varphi }_{2\mathrm{}m\omega }(r)`$
$`+2({\displaystyle \frac{2+2f(r)}{r}}f^{}(r)i\omega ){\displaystyle \frac{d}{dr}}\stackrel{~}{\varphi }_{2\mathrm{}m\omega }(r)2{\displaystyle \frac{rf^{}(r)+(\mathrm{}1)(\mathrm{}+2)+3i\omega r}{r^2}}\stackrel{~}{\varphi }_{2\mathrm{}m\omega }(r)],`$
$`=`$ $`{\displaystyle \frac{2r^2f(r)}{\sqrt{2\mathrm{}(\mathrm{}+1)}}}[4\pi (2r^2\stackrel{~}{T}_{1\mathrm{}m\omega }(r)+f(r)\stackrel{~}{T}_{1\mathrm{}m\omega }(r))`$
$`f(r)(i\omega +f^{}(r)){\displaystyle \frac{d}{dr}}\stackrel{~}{\varphi }_{0\mathrm{}m\omega }(r){\displaystyle \frac{(1+i\omega r)(rf^{}(r)+i\omega r)}{r^2}}\stackrel{~}{\varphi }_{0\mathrm{}m\omega }(r)`$
$`+{\displaystyle \frac{2}{r}}(ir\omega 17f(r)){\displaystyle \frac{d}{dr}}\stackrel{~}{\varphi }_{2\mathrm{}m\omega }(r)+2({\displaystyle \frac{\omega ^2+i\omega /r}{f(r)}}+{\displaystyle \frac{5f(r)1}{r^2}})\stackrel{~}{\varphi }_{2\mathrm{}m\omega }(r)],`$
where $`\stackrel{~}{T}_{s\mathrm{}m\omega }(r)`$ is the Fourier-harmonic coefficient of the source term in the Teukolsky equation (25). Interestingly, Eq. (36) has the same form as the $`s=0`$ Teukolsky equation. Once we obtain $`M_{0\mathrm{}m\omega }(r)`$ by solving the above equation, $`M_{1\mathrm{}m\omega }(r)`$ is derived from

$`M_{1\mathrm{}m\omega }(r)`$ $`=`$ $`{\displaystyle \frac{ir^2}{\omega \sqrt{2\mathrm{}(\mathrm{}+1)}}}\left[f(r){\displaystyle \frac{d}{dr}}\stackrel{~}{\varphi }_{0\mathrm{}m\omega }(r)+{\displaystyle \frac{1+i\omega r}{r}}\stackrel{~}{\varphi }_{0\mathrm{}m\omega }(r)+2{\displaystyle \frac{d}{dr}}\stackrel{~}{\varphi }_{2\mathrm{}m\omega }(r)+2{\displaystyle \frac{f(r)i\omega r}{rf(r)}}\stackrel{~}{\varphi }_{2\mathrm{}m\omega }(r)\right]`$ (37)
$`{\displaystyle \frac{i}{\omega }}\left({\displaystyle \frac{d}{dr}}M_{0\mathrm{}m\omega }(r){\displaystyle \frac{M_{0\mathrm{}m\omega }(r)}{r}}\right),`$
which also follows from Eqs. (33) and (34).

## III Summary and Discussion

In this paper, to solve the gauge problem of the gravitational self-force, we have considered the gauge transformation from the Regge-Wheeler gauge to the harmonic gauge and have presented a formalism to obtain the infinitesimal displacement vector of this transformation, $`\xi ^\mu `$. First, we have performed the Fourier-harmonic expansion of $`\xi ^\mu `$ and divided it into the odd and even parity parts. The odd part has only one degree of freedom and it turns out that the gauge transformation can be found by solving a single second-order differential equation for the radial function. As for the even parity part, we have further divided it into scalar and vector parts where the scalar part is given by the gradient of a scalar function and the vector part is divergence-free. The scalar part has by definition only one degree of freedom, and we have found that it can be obtained by solving two second-order differential equations consecutively. These two equations are found to be identical to the $`s=0`$ Teukolsky equation. The vector part has two degrees of freedom, and the gauge transformation equations give equations that are coupled in a complicated way. However, by introducing two auxiliary variables which satisfy the $`s=\pm 1`$ Teukolsky equations, we have succeeded in deriving a decoupled second-order equation for one of the gauge functions with the source term given by the auxiliary variables. Interestingly, this second-order equation has the same form as the $`s=0`$ Teukolsky equation. The other gauge function is then simply given by applying a differential operator to the first.

Since all the equations to be solved have the form analogous to or equal to the Regge-Wheeler equation, we can derive analytic expressions for their homogeneous solutions by using the Mano-Suzuki-Takasugi method Mano1 , and construct the Green function from these homogeneous solutions. So we conclude that the gauge transformation can be solved by using the Green function method, and we can construct the metric perturbation in the harmonic gauge. In practice, however, it may not be easy to solve for the gauge transformation since it involves products of Green functions with double integrals. Derivation of the gauge transformation functions in a closed, practically tractable form is left for future study.

Another approach to the gauge problem is to consider the self-force in a gauge different from the harmonic gauge, similar to (but very different in principle from) the intermediate gauge approach proposed by Barack and Ori Barack:2001ph . Here the recent result by Detweiler and Whiting Detweiler:2002mi  becomes crucial. Their observation that the S part and the R part play the identical roles as the direct part and the tail part, respectively, and that the S part satisfies the same inhomogeneous Einstein equations as the full metric perturbation enables us to define the S part and the R part of the metric perturbation unambiguously in an arbitrary gauge as long as the gauge condition is consistent with the Einstein equations. For example, given the S part of the metric perturbation in the harmonic gauge, one can perform the gauge transformation of it to the RW gauge and the resulting metric perturbation which satisfies the Einstein equations can be identified as the S part of the metric perturbation in the RW gauge. Then, after solving the Regge-Wheeler-Zerilli equations to obtain the full metric perturbation, it is straightforward to derive the R part of the metric perturbation in the RW gauge Mino . The calculation of the self-force in the RW gauge in this manner is in progress SNS2 .

Finally, we comment on the self-force in the case of the Kerr background. In the Schwarzschild case, it was possible to use the Regge-Wheeler-Zerilli formalism to obtain the metric perturbation in the RW gauge. However, in the Kerr case, there is no known gauge in which the full metric perturbation can be calculated. The Chrzanowski method Chrzanowski:wv  based on the Teukolsky formalism can give the metric perturbation in the (ingoing or outgoing) radiation gauge, but only outside the range of radial coordinates the orbit resides in. One possible way to circumvent this difficulty is to consider first the regularization of the Weyl scalar $`\mathrm{\Psi }_4`$. Given an orbit, $`\mathrm{\Psi }_4`$ can be calculated by the Teukolsky formalism, and the S part of it, $`\mathrm{\Psi }_4^\mathrm{S}`$, can be calculated from the S part of the metric perturbation in the harmonic gauge, $`h_{\mu \nu }^{\mathrm{S},\mathrm{H}}`$,

$`\mathrm{\Psi }_4^\mathrm{S}=\widehat{\mathrm{\Psi }}_4[h_{\mu \nu }^{\mathrm{S},\mathrm{H}}],`$ (38)
where $`\widehat{\mathrm{\Psi }}_4`$ is the operator to derive the Wely scalar from a given metric perturbation. Then the R part of $`\mathrm{\Psi }_4`$ can be derived by subtracting the S part from the Weyl scalar,

$`\mathrm{\Psi }_4^\mathrm{R}=\mathrm{\Psi }_4\mathrm{\Psi }_4^\mathrm{S}.`$ (39)
Now $`\mathrm{\Psi }_4^\mathrm{R}`$ satisfies the homogeneous Teukolsky equation. Hence using the Chrzanowski method, we may construct the R part of the metric perturbation in the radiation gauge and derive the self-force. Since this procedure involves many derivative operations, the metric perturbation $`h_{\mu \nu }^{\mathrm{S},\mathrm{H}}`$ has to be evaluated to with a sufficiently high accuracy which may be practically a difficult task, if not impossible. Feasibility of this method should surely be investigated.

###### Acknowledgements.

We would like to thank T. Tanaka and Y. Mino for useful discussions. HN and MS would like to thank all the participants at the Radiation Reaction Focus Session and the 5th Capra Ranch Meeting, PSU, State College in U.S.A., for invaluable discussions. This work was supported in part by Monbukagaku-sho Grant-in-Aid for Scientific Research Nos. 1047214 and 12640269, and by the Center for Gravitational Wave Physics, which is funded by the National Science Foundation under Cooperative Agreement PHY 0114375.

## Appendix A The field equation on the Regge-Wheeler gauge

In this appendix, we recapitulate the equations for the Regge-Wheeler-Zerilli formalism. In doing so, we correct some minor errors in Zerilli’s paper Zerilli:wd . Here an equation number given as (Z:1) denotes the equation (1) in Zerilli’s paper for comparison, and a label $`[CRTD]`$ to an equation means it is corrected.

We consider the linearized Einstein equations for the perturbed metric,

$$\stackrel{~}{g}_{\mu \nu }=g_{\mu \nu }+h_{\mu \nu },$$

where $`g_{\mu \nu }`$ is the background metric. Then the Einstein tensor and the stress-energy tensor up to the linear order can be expanded as

$`G_{\mu \nu }[\stackrel{~}{g}_{\mu \nu }]`$ $`=`$ $`G_{\mu \nu }[g_{\mu \nu }]+\delta G_{\mu \nu }[h_{\mu \nu }]+O(h^2),`$ (40)
$`\stackrel{~}{T}_{\mu \nu }`$ $`=`$ $`T_{\mu \nu }+\delta T_{\mu \nu },`$ (41)
where $`f_\mu =h_{\mu \alpha }^{}{}_{}{}^{;\alpha }`$ and

$`\delta G_{\mu \nu }[h_{\mu \nu }]`$ $`=`$ $`{\displaystyle \frac{1}{2}}h_{\mu \nu ;\alpha }{}_{}{}^{;\alpha }+f_{(\mu ;\nu )}R_{\alpha \mu \beta \nu }h^{\alpha \beta }{\displaystyle \frac{1}{2}}h_{;\mu ;\nu }+R^\alpha {}_{(\mu }{}^{}h_{\nu )\alpha }^{}`$ (42)
$`{\displaystyle \frac{1}{2}}g_{\mu \nu }(f_\lambda {}_{}{}^{;\lambda }h_{;\lambda }{}_{}{}^{;\lambda }){\displaystyle \frac{1}{2}}h_{\mu \nu }R+{\displaystyle \frac{1}{2}}g_{\mu \nu }h_{\alpha \beta }R^{\alpha \beta }.`$
When the background is Ricci flat, $`R_{\mu \nu }^{(\mathrm{b})}=0`$, the above equation is rewritten as

$`{\displaystyle \frac{1}{2}}h_{\mu \nu ;\alpha }{}_{}{}^{;\alpha }+f_{(\mu ;\nu )}R_{\alpha \mu \beta \nu }h^{\alpha \beta }{\displaystyle \frac{1}{2}}h_{;\mu ;\nu }{\displaystyle \frac{1}{2}}g_{\mu \nu }(f_\lambda {}_{}{}^{;\lambda }h_{;\lambda }{}_{}{}^{;\lambda })=8\pi \delta T_{\mu \nu }.`$ (43)
We apply the above to the case of the Schwarzschild background, and expand $`h_{\mu \nu }`$ ((Z:D2a) and (Z:D2b)) and $`\delta T_{\mu \nu }`$ in tensor harmonics,

$`𝒉`$ $`=`$ $`{\displaystyle \underset{\mathrm{}m}{}}[f(r)H_{0\mathrm{}m}(t,r)𝒂_\mathrm{}m^{(0)}i\sqrt{2}H_{1\mathrm{}m}(t,r)𝒂_\mathrm{}m^{(1)}+{\displaystyle \frac{1}{f(r)}}H_{2\mathrm{}m}(t,r)𝒂_\mathrm{}m`$ (44)
$`{\displaystyle \frac{i}{r}}\sqrt{2\mathrm{}(\mathrm{}+1)}h_{0\mathrm{}m}^{(e)}(t,r)𝒃_\mathrm{}m^{(0)}+{\displaystyle \frac{1}{r}}\sqrt{2\mathrm{}(\mathrm{}+1)}h_{1\mathrm{}m}^{(e)}(t,r)𝒃_\mathrm{}m`$
$`+\sqrt{{\displaystyle \frac{1}{2}}\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)}G_\mathrm{}m(t,r)𝒇_\mathrm{}m+\left(\sqrt{2}K_\mathrm{}m(t,r){\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{\sqrt{2}}}G_\mathrm{}m(t,r)\right)𝒈_\mathrm{}m`$
$`{\displaystyle \frac{\sqrt{2\mathrm{}(\mathrm{}+1)}}{r}}h_{0\mathrm{}m}(t,r)𝒄_\mathrm{}m^{(0)}+{\displaystyle \frac{i\sqrt{2\mathrm{}(\mathrm{}+1)}}{r}}h_{1\mathrm{}m}(t,r)𝒄_\mathrm{}m`$
$`+{\displaystyle \frac{\sqrt{2\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)}}{2r}}h_{2\mathrm{}m}(t,r)𝒅_\mathrm{}m][CRTD],`$
$`\delta 𝑻`$ $`=`$ $`{\displaystyle \underset{\mathrm{}m}{}}[A_\mathrm{}m^{(0)}𝒂_\mathrm{}m^{(0)}+A_\mathrm{}m^{(1)}𝒂_\mathrm{}m^{(1)}+A_\mathrm{}m𝒂_\mathrm{}m+B_\mathrm{}m^{(0)}𝒃_\mathrm{}m^{(0)}+B_\mathrm{}m𝒃_\mathrm{}m+Q_\mathrm{}m^{(0)}𝒄_\mathrm{}m^{(0)}+Q_\mathrm{}m𝒄_\mathrm{}m`$ (45)
$`+D_\mathrm{}m𝒅_\mathrm{}m+G_\mathrm{}m^{(s)}𝒈_\mathrm{}m+F_\mathrm{}m𝒇_\mathrm{}m],`$
where we use $`h_{0\mathrm{}m}^{(e)}`$ and $`h_{1\mathrm{}m}^{(e)}`$ for the even part coefficients instead of $`h_{0\mathrm{}m}^{(m)}`$ and $`h_{0\mathrm{}m}^{(m)}`$, respectively, in Zerilli’s paper, and the coefficient $`G_\mathrm{}m^{(s)}`$ instead of the Zerilli’s notation $`G_\mathrm{}m`$ for the energy-momentum tensor, and $`𝒂_\mathrm{}m^{(0)}`$, $`𝒂_\mathrm{}m,\mathrm{}`$ are the ten tensor harmonics (Z:A2a-j) defined as

$`𝐚_\mathrm{}m^{(0)}`$ $`=`$ $`\left(\begin{array}{cccc}Y_\mathrm{}m& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),`$ (50)
$`𝐚_\mathrm{}m^{(1)}`$ $`=`$ $`(i/\sqrt{2})\left(\begin{array}{cccc}0& Y_\mathrm{}m& 0& 0\\ Sym& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),`$ (55)
$`𝐚_\mathrm{}m`$ $`=`$ $`\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& Y_\mathrm{}m& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),`$ (60)
$`𝐛_\mathrm{}m^{(0)}`$ $`=`$ $`ir[2\mathrm{}(\mathrm{}+1)]^{1/2}\left(\begin{array}{cccc}0& 0& (/\theta )Y_\mathrm{}m& (/\varphi )Y_\mathrm{}m\\ 0& 0& 0& 0\\ Sym& 0& 0& 0\\ Sym& 0& 0& 0\end{array}\right),`$ (65)
$`𝐛_\mathrm{}m`$ $`=`$ $`r[2\mathrm{}(\mathrm{}+1)]^{1/2}\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& (/\theta )Y_\mathrm{}m& (/\varphi )Y_\mathrm{}m\\ 0& Sym& 0& 0\\ 0& Sym& 0& 0\end{array}\right),`$ (70)
$`𝐜_\mathrm{}m^{(0)}`$ $`=`$ $`r[2\mathrm{}(\mathrm{}+1)]^{1/2}\left(\begin{array}{cccc}0& 0& (1/\mathrm{sin}\theta )(/\varphi )Y_\mathrm{}m& \mathrm{sin}\theta (/\theta )Y_\mathrm{}m\\ 0& 0& 0& 0\\ Sym& 0& 0& 0\\ Sym& 0& 0& 0\end{array}\right),`$ (75)
$`𝐜_\mathrm{}m`$ $`=`$ $`ir[2\mathrm{}(\mathrm{}+1)]^{1/2}\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& (1/\mathrm{sin}\theta )(/\varphi )Y_\mathrm{}m& \mathrm{sin}\theta (/\theta )Y_\mathrm{}m\\ 0& Sym& 0& 0\\ 0& Sym& 0& 0\end{array}\right),`$ (80)
$`𝐝_\mathrm{}m`$ $`=`$ $`ir^2[2\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)]^{1/2}\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& (1/\mathrm{sin}\theta )X_\mathrm{}m& \mathrm{sin}\theta W_\mathrm{}m\\ 0& 0& Sym& \mathrm{sin}\theta X_\mathrm{}m\end{array}\right),`$ (85)
$`𝐠_\mathrm{}m`$ $`=`$ $`(r^2/\sqrt{2})\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& Y_\mathrm{}m& 0\\ 0& 0& 0& \mathrm{sin}^2\theta Y_\mathrm{}m\end{array}\right),`$ (90)
$`𝐟_\mathrm{}m`$ $`=`$ $`r^2[2\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)]^{1/2}\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& W_\mathrm{}m& X_\mathrm{}m\\ 0& 0& Sym& \mathrm{sin}^2\theta W_\mathrm{}m\end{array}\right)[CRTD].`$ (95)
Here the angular functions $`X_\mathrm{}m`$ and $`W_\mathrm{}m`$ are given by

$`X_\mathrm{}m`$ $`=`$ $`2{\displaystyle \frac{}{\varphi }}\left({\displaystyle \frac{}{\theta }}\mathrm{cot}\theta \right)Y_\mathrm{}m,`$ (96)
$`W_\mathrm{}m`$ $`=`$ $`\left({\displaystyle \frac{^2}{\theta ^2}}\mathrm{cot}\theta {\displaystyle \frac{}{\theta }}{\displaystyle \frac{1}{\mathrm{sin}^2\theta }}{\displaystyle \frac{^2}{\varphi ^2}}Y_\mathrm{}m\right).`$ (97)
For a point particle moving along a geodesic, the stress-energy tensor takes the form,

$`T^{\mu \nu }`$ $`=`$ $`\mu {\displaystyle _{\mathrm{}}^+\mathrm{}}\delta ^{(4)}(xz(\tau )){\displaystyle \frac{dz^\mu }{d\tau }}{\displaystyle \frac{dz^\nu }{d\tau }}𝑑\tau `$ (98)
$`=`$ $`\mu \gamma {\displaystyle \frac{dz^\mu }{dt}}{\displaystyle \frac{dz^\nu }{dt}}{\displaystyle \frac{\delta (rR(t))}{r^2}}\delta ^{(2)}(\mathrm{\Omega }\mathrm{\Omega }(t)),`$
where the following notation for the particle orbit is used.

$`z^\mu `$ $`=`$ $`z^\mu (\tau )=\{T(\tau ),R(\tau ),\mathrm{\Theta }(\tau ),\mathrm{\Phi }(\tau )\},`$ (99)
$`\gamma `$ $`=`$ $`{\displaystyle \frac{dT(\tau )}{d\tau }}.`$ (100)
This stress-energy tensor is expressed in terms of the tensor harmonics as given in Table I (correspond to (Z:Table III)).

Substituting Eqs. (44) and (45) into Eq. (43), we obtain the field equations for each harmonic mode. For the odd part which has the odd parity $`(1)^{\mathrm{}+1}`$, in the RW gauge in which $`h_2=0`$, the following three equations (Z:C6a-c) are derived.

$`{\displaystyle \frac{^2h_0}{r^2}}{\displaystyle \frac{^2h_1}{tr}}{\displaystyle \frac{2}{r}}{\displaystyle \frac{h_1}{t}}+\left[{\displaystyle \frac{4M}{r^2}}{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r}}\right]{\displaystyle \frac{h_0}{r2M}}`$ $`=`$ $`{\displaystyle \frac{8\pi }{\sqrt{\mathrm{}(\mathrm{}+1)/2}}}{\displaystyle \frac{r^2}{r2M}}Q_\mathrm{}m^{(0)}[CRTD],`$ (101)
$`{\displaystyle \frac{^2h_1}{t^2}}{\displaystyle \frac{^2h_0}{tr}}+{\displaystyle \frac{2}{r}}{\displaystyle \frac{h_0}{t}}+{\displaystyle \frac{(\mathrm{}1)(\mathrm{}+2)(r2M)}{r^3}}h_1`$ $`=`$ $`{\displaystyle \frac{8\pi i(r2M)}{\sqrt{\mathrm{}(\mathrm{}+1)/2}}}Q_\mathrm{}m[CRTD],`$ (102)
$`{\displaystyle \frac{}{r}}\left[\left(1{\displaystyle \frac{2M}{r}}\right)h_1\right]{\displaystyle \frac{r}{r2M}}{\displaystyle \frac{h_0}{t}}`$ $`=`$ $`{\displaystyle \frac{8\pi ir^2}{\sqrt{\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)/2}}}D_\mathrm{}m[CRTD].`$ (103)
Fo the even part which has the even parity $`(1)^{\mathrm{}}`$, we have seven equations (Z:C7a-g) in the RW gauge in which $`h_0^{(e)}=h_1^{(e)}=G=0`$.

$`\left(1{\displaystyle \frac{2M}{r}}\right)^2{\displaystyle \frac{^2K}{r^2}}+{\displaystyle \frac{1}{r}}\left(1{\displaystyle \frac{2M}{r}}\right)\left(3{\displaystyle \frac{5M}{r}}\right){\displaystyle \frac{K}{r}}`$
$`{\displaystyle \frac{1}{r}}\left(1{\displaystyle \frac{2M}{r}}\right)^2{\displaystyle \frac{H_2}{r}}{\displaystyle \frac{1}{r^2}}\left(1{\displaystyle \frac{2M}{r}}\right)(H_2K)`$
$`{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{2r^2}}\left(1{\displaystyle \frac{2M}{r}}\right)(H_2+K)`$ $`=`$ $`8\pi A_\mathrm{}m^{(0)}[CRTD],`$ (104)
$`{\displaystyle \frac{}{t}}\left[{\displaystyle \frac{K}{r}}+{\displaystyle \frac{1}{r}}(KH_2){\displaystyle \frac{M}{r(r2M)}}K\right]`$
$`{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{2r^2}}H_1`$ $`=`$ $`4\sqrt{2}\pi iA_\mathrm{}m^{(1)}[CRTD],`$ (105)
$`\left({\displaystyle \frac{r}{r2M}}\right)^2{\displaystyle \frac{^2K}{t^2}}{\displaystyle \frac{rM}{r(r2M)}}{\displaystyle \frac{K}{r}}{\displaystyle \frac{2}{r2M}}{\displaystyle \frac{H_1}{t}}`$
$`+{\displaystyle \frac{1}{r}}{\displaystyle \frac{H_0}{r}}+{\displaystyle \frac{1}{r(r2M)}}(H_2K)`$
$`+{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{2r(r2M)}}(KH_0)`$ $`=`$ $`8\pi A_\mathrm{}m[CRTD],`$ (106)
$`{\displaystyle \frac{}{r}}\left[\left(1{\displaystyle \frac{2M}{r}}\right)H_1\right]{\displaystyle \frac{}{t}}(H_2+K)`$ $`=`$ $`{\displaystyle \frac{8\pi ir}{\sqrt{\mathrm{}(\mathrm{}+1)/2}}}B_\mathrm{}m^{(0)}[CRTD],`$ (107)
$`{\displaystyle \frac{H_1}{t}}+\left(1{\displaystyle \frac{2M}{r}}\right){\displaystyle \frac{}{r}}(H_0K)+{\displaystyle \frac{2M}{r^2}}H_0`$
$`+{\displaystyle \frac{1}{r}}\left(1{\displaystyle \frac{M}{r}}\right)(H_2H_0)`$ $`=`$ $`{\displaystyle \frac{8\pi (r2M)}{\sqrt{\mathrm{}(\mathrm{}+1)/2}}}B_\mathrm{}m[CRTD],`$ (108)
$`{\displaystyle \frac{r}{r2M}}{\displaystyle \frac{^2K}{t^2}}+\left(1{\displaystyle \frac{2M}{r}}\right){\displaystyle \frac{^2K}{r^2}}+{\displaystyle \frac{2}{r}}\left(1{\displaystyle \frac{M}{r}}\right){\displaystyle \frac{K}{r}}`$
$`{\displaystyle \frac{r}{r2M}}{\displaystyle \frac{^2H_2}{t^2}}+2{\displaystyle \frac{^2H_1}{tr}}\left(1{\displaystyle \frac{2M}{r}}\right){\displaystyle \frac{^2H_0}{r^2}}`$
$`+{\displaystyle \frac{2(rM)}{r(r2M)}}{\displaystyle \frac{H_1}{t}}{\displaystyle \frac{1}{r}}\left(1{\displaystyle \frac{M}{r}}\right){\displaystyle \frac{H_2}{r}}{\displaystyle \frac{r+M}{r^2}}{\displaystyle \frac{H_0}{r}}`$
$`+{\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{2r^2}}(H_0H_2)`$ $`=`$ $`8\sqrt{2}\pi G_\mathrm{}m^{(s)}[CRTD],`$ (109)
$`{\displaystyle \frac{H_0H_2}{2}}`$ $`=`$ $`{\displaystyle \frac{8\pi r^2F_\mathrm{}m}{\sqrt{\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)/2}}}[CRTD].`$ (110)
We now consider the Fourier transform of the above field equations. The Fourier coefficients are defined, for example, as

$`h_{0\mathrm{}m\omega }(r)={\displaystyle _{\mathrm{}}^+\mathrm{}}𝑑th_{0\mathrm{}m}(t,r)e^{i\omega t}.`$ (111)
Then we derive the Regge-Wheeler-Zerilli equations and construct the metric perturbation under the RW gauge condition.

For the odd part, a new radial function $`R_{\mathrm{}m\omega }^{(\mathrm{odd})}(r)`$ is introduced, in terms of which the two radial functions $`h_{0\mathrm{}m\omega }`$ and $`h_{1\mathrm{}m\omega }`$ for the metric perturbation are expressed as

$`h_{1\mathrm{}m\omega }`$ $`=`$ $`{\displaystyle \frac{r^2}{r2M}}R_{\mathrm{}m\omega }^{(\mathrm{odd})},`$ (112)
$`h_{0\mathrm{}m\omega }`$ $`=`$ $`{\displaystyle \frac{i}{\omega }}{\displaystyle \frac{d}{dr^{}}}(rR_{\mathrm{}m\omega }^{(\mathrm{odd})}){\displaystyle \frac{8\pi r(r2M)}{\omega [\frac{1}{2}\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)]^{1/2}}}D_{\mathrm{}m\omega }.`$ (113)
The new radial function satisfies the Regge-Wheeler equation (Z:11),

$`{\displaystyle \frac{d^2R_{\mathrm{}m\omega }^{(\mathrm{odd})}}{dr^2}}`$ $`+`$ $`[\omega ^2V_{\mathrm{}}^{(\mathrm{odd})}(r)]R_{\mathrm{}m\omega }^{(\mathrm{odd})}`$ (114)
$`=`$ $`{\displaystyle \frac{8\pi i}{[\frac{1}{2}l(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)]^{1/2}}}{\displaystyle \frac{r2M}{r^2}}`$
$`\times \left(r^2{\displaystyle \frac{d}{dr}}[(1{\displaystyle \frac{2M}{r}})D_{\mathrm{}m\omega }]+(r2M)[(\mathrm{}1)(\mathrm{}+2)]^{1/2}Q_{\mathrm{}m\omega }\right)[CRTD],`$
where $`r^{}=r+2M\mathrm{log}(r/2M1)`$ and

$`V_{\mathrm{}}^{(\mathrm{odd})}(r)=\left(1{\displaystyle \frac{2M}{r}}\right)\left({\displaystyle \frac{\mathrm{}(\mathrm{}+1)}{r^2}}{\displaystyle \frac{6M}{r^3}}\right).`$ (115)
For the even part, a new radial function $`R_{\mathrm{}m\omega }^{(\mathrm{even})}(r)`$ is introduced, in terms of which the four radial functions (Z:13-16) are expressed as

$`K_{\mathrm{}m\omega }`$ $`=`$ $`{\displaystyle \frac{\lambda (\lambda +1)r^2+3\lambda Mr+6M^2}{r^2(\lambda r+3M)}}R_{\mathrm{}m\omega }^{(\mathrm{even})}+{\displaystyle \frac{r2M}{r}}{\displaystyle \frac{dR_{\mathrm{}m\omega }^{(\mathrm{even})}}{dr}}`$ (116)
$`{\displaystyle \frac{r(r2M)}{\lambda r+3M}}\stackrel{~}{C}_{1\mathrm{}m\omega }+{\displaystyle \frac{i(r2M)^2}{r(\lambda r+3M)}}\stackrel{~}{C}_{2\mathrm{}m\omega }[CRTD],`$
$`H_{1\mathrm{}m\omega }`$ $`=`$ $`i\omega {\displaystyle \frac{\lambda r^23\lambda Mr3M^2}{(r2M)(\lambda r+3M)}}R_{\mathrm{}m\omega }^{(\mathrm{even})}i\omega r{\displaystyle \frac{dR_{\mathrm{}m\omega }^{(\mathrm{even})}}{dr}}`$ (117)
$`+{\displaystyle \frac{i\omega r^3}{\lambda r+3M}}\stackrel{~}{C}_{1\mathrm{}m\omega }+{\displaystyle \frac{\omega r(r2M)}{\lambda r+3M}}\stackrel{~}{C}_{2\mathrm{}m\omega }[CRTD],`$
$`H_{0\mathrm{}m\omega }`$ $`=`$ $`{\displaystyle \frac{\lambda r(r2M)\omega ^2r^4+M(r3M)}{(r2M)(\lambda r+3M)}}K_{\mathrm{}m\omega }+{\displaystyle \frac{M(\lambda +1)\omega ^2r^3}{i\omega r(\lambda r+3M)}}H_{1\mathrm{}m\omega }+\stackrel{~}{B}_{\mathrm{}m\omega }[CRTD],`$ (118)
$`H_{2\mathrm{}m\omega }`$ $`=`$ $`H_{0\mathrm{}m\omega }16\pi r^2[{\displaystyle \frac{1}{2}}\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)]^{1/2}F_{\mathrm{}m\omega }[CRTD],`$ (119)
where we have introduced the symbol $`\lambda `$ for

$`\lambda ={\displaystyle \frac{1}{2}}(\mathrm{}1)(\mathrm{}+2),`$ (120)
and the local source terms ((Z:17), (Z:20) and (Z:21)) by

$`\stackrel{~}{B}_{\mathrm{}m\omega }`$ $`=`$ $`{\displaystyle \frac{8\pi r^2(r2M)}{\lambda r+3M}}\{A_{\mathrm{}m\omega }+[{\displaystyle \frac{1}{2}}\mathrm{}(\mathrm{}+1)]^{1/2}B_{\mathrm{}m\omega }\}{\displaystyle \frac{4\pi \sqrt{2}}{\lambda r+3M}}{\displaystyle \frac{Mr}{\omega }}A_{\mathrm{}m\omega }^{(1)},`$ (121)
$`\stackrel{~}{C}_{1\mathrm{}m\omega }`$ $`=`$ $`{\displaystyle \frac{8\pi }{\sqrt{2}\omega }}A_{\mathrm{}m\omega }^{(1)}+{\displaystyle \frac{1}{r}}\stackrel{~}{B}_{\mathrm{}m\omega }16\pi r[{\displaystyle \frac{1}{2}}\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)]^{1/2}F_{\mathrm{}m\omega }[CRTD],`$ (122)
$`\stackrel{~}{C}_{2\mathrm{}m\omega }`$ $`=`$ $`{\displaystyle \frac{8\pi r^2}{i\omega }}{\displaystyle \frac{[\frac{1}{2}\mathrm{}(\mathrm{}+1)]^{1/2}}{r2M}}B_{\mathrm{}m\omega }{}_{}{}^{(0)}{\displaystyle \frac{ir}{r2M}}\stackrel{~}{B}_{\mathrm{}m\omega }`$ (123)
$`+{\displaystyle \frac{16\pi ir^3}{r2M}}[{\displaystyle \frac{1}{2}}\mathrm{}(\mathrm{}+1)(\mathrm{}1)(\mathrm{}+2)]^{1/2}F_{\mathrm{}m\omega }[CRTD].`$
We note that the above radial functions for the metric perturbation have the local source terms which have the $`\delta `$-function behavior at the particle location. The new radial function obeys the wave equation,

$`{\displaystyle \frac{d^2R_{\mathrm{}m\omega }^{(\mathrm{even})}}{dr^2}}+\left[\omega ^2V_{\mathrm{}}^{(\mathrm{even})}(r)\right]R_{\mathrm{}m\omega }^{(\mathrm{even})}=S_{\mathrm{}m\omega },`$ (124)
where

$`V_{\mathrm{}}^{(\mathrm{even})}(r)=\left(1{\displaystyle \frac{2M}{r}}\right){\displaystyle \frac{2\lambda ^2(\lambda +1)r^3+6\lambda ^2Mr^2+18\lambda M^2r+18M^3}{r^3(\lambda r+3M)^2}},`$ (125)
and the source term is

$`S_{\mathrm{}m\omega }`$ $`=`$ $`i{\displaystyle \frac{r2M}{r}}{\displaystyle \frac{d}{dr}}\left[{\displaystyle \frac{(r2M)^2}{r(\lambda r+3M)}}\left({\displaystyle \frac{ir^2}{r2M}}\stackrel{~}{C}_{1\mathrm{}m\omega }+\stackrel{~}{C}_{2\mathrm{}m\omega }\right)\right]`$ (126)
$`+i{\displaystyle \frac{(r2M)^2}{r(\lambda r+3M)^2}}\left[{\displaystyle \frac{\lambda (\lambda +1)r^2+3\lambda Mr+6M^2}{r^2}}\stackrel{~}{C}_{2\mathrm{}m\omega }+i{\displaystyle \frac{\lambda r^23\lambda Mr3M^2}{r2M}}\stackrel{~}{C}_{1\mathrm{}m\omega }\right].`$
The above equation (124) is called the Zerilli equation. The Zerilli equation can be transformed to the Regge-Wheeler equation by the Chandrasekhar transformation chandra . So we may focus only on the Regge-Wheeler equation if desired. The Regge-Wheeler homogeneous solutions are discussed in detail by Mano et al. Mano1 . Using their method, one can construct the retarded Green function to solve the inhomogeneous Regge-Wheeler equation. Then the metric perturbation in the RW gauge is obtained from Eqs. (112) and (113) for the odd part and from Eqs. (116-119) for the even part.