text:
# Reconstruction of Black Hole Metric Perturbations from Weyl Curvature

## I Introduction

The spherically symmetry of a Schwarzschild black hole background allows for a multipole decomposition of metric perturbations, even in the time domain. These were studied originally by Regge and Wheeler Regge57  for odd-parity perturbations and by Zerilli Zerilli70a  for the even-parity case. Moncrief Moncrief74  has given a gauge-invariant formulation of the problem, in terms of the three-metric perturbations. The two degrees of freedom of the gravitational field are described in terms of two waveforms, $`\varphi ^\pm `$ satisfying a simple wave equation

$$\frac{^2\varphi _{(\mathrm{}m)}^\pm }{t^2}+\frac{^2\varphi _{(\mathrm{}m)}^\pm }{r^2}V_\mathrm{}^\pm (r)\varphi _{(\mathrm{}m)}^\pm =0.$$
(1)
Here $`r^{}r+2M\mathrm{ln}(r/2M1)`$, and $`V_\mathrm{}^\pm (r)`$ are the Zerilli and Regge–Wheeler potentials respectively.

Given the solution to the wave equation (1) one can reconstruct explicitly both the even and odd parity metric perturbations of a Schwarzschild background in the Regge–Wheeler gauge Lousto99b ; Nicasio:2000ge . This permits a complete description of the perturbative spacetime. But the Regge–Wheeler gauge is unfortunately not asymptotically flat, and in order to extract radiation information for a second order perturbative expansion one has to perform a new gauge transformation Gleiser:1998rw . Moreover, the desirable properties of the the Regge–Wheeler gauge being unique and invertible are lost in the case the background is a rotating black hole, i.e. represented by the Kerr metric, where no effective multipole decomposition is yet known.

There is an independent formulation of the perturbation problem derived from the Newman-Penrose formalism Newman62a  that is valid for perturbations of rotating black holes.Teukolsky73  This formulation fully exploits the null structure of black holes to decouple the curvature perturbation equations into a single wave equation that, in Boyer-Lindquist coordinates $`(t,r,\theta ,\phi )`$, can be written as:

$`\left\{\right[a^2\mathrm{sin}^2\theta {\displaystyle \frac{(r^2+a^2)^2}{\mathrm{\Delta }}}]_{tt}{\displaystyle \frac{4Mar}{\mathrm{\Delta }}}_{t\phi }`$
$`2s\left[(r+ia\mathrm{cos}\theta ){\displaystyle \frac{M(r^2a^2)}{\mathrm{\Delta }}}\right]_t`$
$`+\mathrm{\Delta }^s_r\left(\mathrm{\Delta }^{s+1}_r\right)+{\displaystyle \frac{1}{\mathrm{sin}\theta }}_\theta \left(\mathrm{sin}\theta _\theta \right)`$
$`+\left({\displaystyle \frac{1}{\mathrm{sin}^2\theta }}{\displaystyle \frac{a^2}{\mathrm{\Delta }}}\right)_{\phi \phi }+\mathrm{\hspace{0.17em}2}s\left[{\displaystyle \frac{a(rM)}{\mathrm{\Delta }}}+{\displaystyle \frac{i\mathrm{cos}\theta }{\mathrm{sin}^2\theta }}\right]_\phi `$
$`(s^2\mathrm{cot}^2\theta s)\}\psi =4\pi \mathrm{\Sigma }T,`$ (2)
where $`M`$ is the mass of the black hole, $`a`$ its angular momentum per unit mass, $`\mathrm{\Sigma }r^2+a^2\mathrm{cos}^2\theta `$, and $`\mathrm{\Delta }r^22Mr+a^2`$. The source term $`T`$ is built up from the energy-momentum tensor Teukolsky73 . Gravitational perturbations, corresponding to $`s=\pm 2`$, are compactly described in terms of contractions of the Weyl tensor with a null tetrad. The components of the tetrad (also given in Ref. Teukolsky73 ) are appropriately chosen along the repeated principal null directions of the background spacetime \[see Eq. (4) below\]. The resulting gauge and (infinitesimally) tetrad invariant components of the Weyl curvature are given by

$$\psi =\{\begin{array}{cc}\rho ^4\psi _4\rho ^4C_{n\overline{m}n\overline{m}}\hfill & \mathrm{for}s=2\hfill \\ \psi _0C_{lmlm}\hfill & \mathrm{for}s=+2\hfill \end{array},$$
(3)
where an overbar means complex conjugation and $`\rho `$ is given in Eq. (I) below. Asymptotically, the leading behavior of the field $`\psi `$ represents either the outgoing radiative part of the perturbed Weyl tensor, ($`s=2`$), or the ingoing radiative part, ($`s=+2`$).

The components of the Boyer–Lindquist null tetrad for the Kerr background are given by

$`(l^\alpha )`$ $`=`$ $`({\displaystyle \frac{r^2+a^2}{\mathrm{\Delta }}},1,0,{\displaystyle \frac{a}{\mathrm{\Delta }}}),`$ (4a)
$`(n^\alpha )`$ $`=`$ $`{\displaystyle \frac{1}{2(r^2+a^2\mathrm{cos}^2\theta )}}(r^2+a^2,\mathrm{\Delta },0,a),`$ (4b)
$`(m^\alpha )`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}(r+ia\mathrm{cos}\theta )}}(ia\mathrm{sin}\theta ,0,1,i/\mathrm{sin}\theta ).`$ (4c)
With this choice of the tetrad the non-vanishing spin coefficients are

$`\rho `$ $`=`$ $`{\displaystyle \frac{1}{(ria\mathrm{cos}\theta )}},\beta =\overline{\rho }{\displaystyle \frac{\mathrm{cot}\theta }{2\sqrt{2}}},`$
$`\pi `$ $`=`$ $`ia\rho ^2{\displaystyle \frac{\mathrm{sin}\theta }{\sqrt{2}}},\tau =ia\rho \overline{\rho }{\displaystyle \frac{\mathrm{sin}\theta }{\sqrt{2}}},`$
$`\mu `$ $`=`$ $`\rho ^2\overline{\rho }{\displaystyle \frac{\mathrm{\Delta }}{2}},\alpha =\pi \overline{\beta },`$
$`\gamma `$ $`=`$ $`\mu +\rho \overline{\rho }{\displaystyle \frac{(rM)}{2}},`$ (5)
and the only non-vanishing Weyl scalar in the background is

$$\psi _2=M\rho ^3.$$
(6)
Analogously to the Zerilli-Regge-Wheeler waveforms, $`\psi _4`$ can be directly used to compute energy and momentum radiated at infinity, but it remains to relate it to metric perturbations. Chandrasekhar Chandrasekhar83  studied a way to to obtain metric perturbations of the Kerr metric, but it was proved by Price and Ipser Price91  that this choice is not a proper gauge, namely is an incomplete constraint on the coordinates.

Chrzanowski Chrzanowski:1975wv  generalized work of Cohen and Kegeles Cohen:1974  on Hertz potentials to the gravitational perturbations of the Kerr metric. In Ref. Chrzanowski:1975wv  explicit expressions are given for homogeneous metric perturbations in the frequency domain. Wald Wald:1978vm  subsequently showed that the expressions given in Ref. Chrzanowski:1975wv  do not lead to real metric perturbations. Cohen and Kegeles Kegeles:1979an  then corrected their expressions and gave explicit equations (see Sec. II) relating metric perturbations to a gravitational Hertz potential, $`\mathrm{\Psi }`$ that fulfills Eq. (2), but that is different from $`\psi _4`$ or $`\psi _0`$. In those works no explicit method was given for determining $`\mathrm{\Psi }`$ in any specific astrophysical problem.

In Sec. III we provide the explicit formulae relating a gravitational Hertz potential to $`\psi _4`$ or $`\psi _0`$ in the time domain, hence to the given initial data defining the astrophysical model one wants to evolve (See Ref. Campanelli98c  for the 3+1 decomposition of the Weyl scalars). These allow one to define the outgoing and ingoing radiation gauges that are asymptotically flat at future infinity and regular on the event horizon respectively. Such gauges have been found Campanelli99  especially well suited for computing second order perturbations of a Kerr hole and, once generalized in presence of matter, can provide a first step toward computing radiation reaction (self-force) corrections Lousto99b  to the trajectory of a particle orbiting a rotating black hole.

Vacuum second order perturbations have a direct application to the close limit expansion of the final merger stage of binary black holes with comparable masses. Perturbation theory in conjunction with limited full numerical simulations has proved to be an extremely powerful tool to compute waveforms from binary black holes from near the innermost stable circular orbit Baker:2001nu ; Baker:2002qf . The combination of radiation reaction correction plus second order perturbations provides a formidable tool for computing gravitational radiation from binary black hole–neutron star systems, and is particularly relevant to the computation of template banks for ground based interferometers such as LIGO/VIRGO/GEO about to enter on-line, as well as space missions such as LISA, sensitive to supermassive black hole binary systems. Thus in order to incorporate these improvements to our theoretical predictions it is imperative to know how to build up explicit metric perturbations around a Kerr background.

## II Formulation of the problem

We use the notation

$$g_{\mu \nu }=g_{\mu \nu }^{Kerr}+h_{\mu \nu }$$
(7)
to describe metric perturbations.

### II.1 ingoing and outgoing radiation gauges

Chrzanowski Chrzanowski:1975wv , and Cohen and Kegeles Kegeles:1979an  found two convenient gauges that allow one to invert the metric perturbations in terms of potentials $`\mathrm{\Psi }_{IRG}`$ or $`\mathrm{\Psi }_{ORG}`$ satisfying the same wave equations as the Weyl scalars $`\rho ^4\psi _4`$ or $`\psi _0`$ respectively.

In the ingoing radiation gauge (IRG) we have

$`h_{ll}`$ $`=`$ $`h_{\mu \nu }l^\mu l^\nu =0;h_{ln}=h_{\mu \nu }l^\mu n^\nu =0,`$
$`h_{lm}`$ $`=`$ $`h_{\mu \nu }l^\mu m^\nu =0;h_{l\overline{m}}=h_{\mu \nu }l^\mu \overline{m}^\nu =0,`$
$`h_{m\overline{m}}`$ $`=`$ $`h_{\mu \nu }m^\mu \overline{m}^\nu =0,`$ (8)
and the homogeneous (for vacuum) metric components can be written, in the time domain, in terms of solutions to the wave equation for $`\rho ^4\psi _4`$ only, as follows:111We have made the choice that $`ϵ=0`$.

$`h_{nn}^{IRG}`$ $`=`$ $`\left\{(\delta +\overline{\alpha }+3\beta \tau )(\delta +4\beta +3\tau )\right\}`$ (9a)
$`(\mathrm{\Psi }_{IRG})+c.c.`$
$`h_{\overline{m}\overline{m}}^{IRG}`$ $`=`$ $`\left\{(D\rho )(D+3\rho )\right\}(\mathrm{\Psi }_{IRG}),`$ (9b)
$`h_{(n\overline{m})}^{IRG}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\{(\delta \overline{\alpha }+3\beta \overline{\pi }\tau )(D+3\rho )+`$ (9c)
$`(D+\overline{\rho }\rho )(\delta +4\beta +3\tau )\}(\mathrm{\Psi }_{IRG}),`$
where $`c.c.`$ stands for the complex conjugate part of the whole object to ensure that the metric be real Wald:1978vm ; Kegeles:1979an . Note that in this gauge the metric potential has the properties of being transverse $`(h_{\mu \nu }l^\mu =0)`$ and traceless $`(h_\mu ^\mu =0)`$ at the future horizon and past null infinity. It is thus a suitable gauge to study gravitational radiation effects near the event horizon.

The complementary (adjoint) gauge to the ingoing radiation gauge is the outgoing radiation gauge (ORG), which can be obtained from the IRG upon exchange of the tetrad vectors $`ln`$, $`\overline{m}m`$ and the appropriate renormalization. It satisfies:

$$h_{nn}=0=h_{ln}=0=h_{nm}=0=h_{n\overline{m}}=0=h_{m\overline{m}}.$$
(10)
The metric potential has now the property of being transverse $`(h_{\mu \nu }n^\mu =0)`$ and traceless $`(h_\mu ^\mu =0)`$ at the past horizon and future null infinity. It is thus an example of a suitable asymptotically flat gauge in which to directly compute radiated energy and momenta at infinity.

In this gauge, the homogeneous metric components can be written in terms of solutions to the wave equation for $`\psi _0`$, as

$`h_{ll}^{ORG}`$ $`=`$ $`\rho ^4\left\{(\overline{\delta }3\alpha \overline{\beta }+5\pi )(\overline{\delta }4\alpha +\pi )\right\}`$ (11a)
$`(\mathrm{\Psi }_{ORG})+c.c.`$
$`h_{mm}^{ORG}`$ $`=`$ $`\rho ^4\left\{(\widehat{\mathrm{\Delta }}+5\mu 3\gamma +\overline{\gamma })(\widehat{\mathrm{\Delta }}+\mu 4\gamma )\right\}`$ (11b)
$`(\mathrm{\Psi }_{ORG}),`$
$`h_{(lm)}^{ORG}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\rho ^4\{(\overline{\delta }3\alpha +\overline{\beta }+5\pi +\overline{\tau })(\widehat{\mathrm{\Delta }}+\mu 4\gamma )`$ (11c)
$`+(\widehat{\mathrm{\Delta }}+5\mu \overline{\mu }3\gamma \overline{\gamma })(\overline{\delta }4\alpha +\pi )\}`$
$`(\mathrm{\Psi }_{ORG}),`$
where the directional derivatives are $`D=l_{}^\mu _\mu ,\widehat{\mathrm{\Delta }}=n_{}^\mu _\mu ,\delta =m_{}^\mu _\mu `$, and the rest of the Greek letters represent the usual notation for spin coefficients.Newman62a

### II.2 4th order equations for the potential

From the evolution of the Teukolsky equation we can obtain $`\psi _4`$ and $`\psi _0`$, the two gauge and tetrad invariant objects representing outgoing and ingoing radiation respectively. To relate the unknown potential $`\mathrm{\Psi }`$ to them we use their definitions (3) to obtain

$`\psi _0`$ $`=`$ $`DDDD\left[\overline{\mathrm{\Psi }}_{IRG}\right],\mathrm{and}`$ (12)
$`\rho ^4\psi _4`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left[\overline{}\overline{}\overline{}\overline{}\overline{\mathrm{\Psi }}_{IRG}12\rho ^3\psi _2_t\mathrm{\Psi }_{IRG}\right],`$ (13)
where $`\overline{}=\overline{ð}+ia\mathrm{sin}\theta _t`$ and $`\overline{ð}=[_\theta +s\mathrm{cot}\theta i\mathrm{csc}\theta _\phi ]`$, while for the outgoing radiation gauge we have

$`\rho ^4\psi _4`$ $`=`$ $`\mathrm{\Delta }^2\widehat{\mathrm{\Delta }}\widehat{\mathrm{\Delta }}\widehat{\mathrm{\Delta }}\widehat{\mathrm{\Delta }}\left[\mathrm{\Delta }^2\overline{\mathrm{\Psi }}_{ORG}\right],\mathrm{and}`$ (14)
$`\psi _0`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left[\overline{\mathrm{\Psi }}_{ORG}+12\rho ^3\psi _2_t\mathrm{\Psi }_{ORG}\right].`$ (15)
In order to obtain an expression for the potentials in terms of the known quantities $`\psi _4`$ or $`\psi _0`$, one has to invert a fourth order differential equation where $`\psi _4`$ or $`\psi _0`$ act as source terms. This will be the central task of our paper.

## III Explicit solution for nonrotating black holes

### III.1 Master equation

The key observation here is that for $`a=0`$ the differential operator $``$ acting on the potentials in Eqs. (13) and Eqs. (15) contains only angular derivatives. Since for the spherically symmetric background we can decompose the potentials into spherical harmonics of spin weight $`s`$ and, from Goldberg:1967 , we have

$`{}_{s}{}^{}\overline{Y}_{\mathrm{}m}^{}`$ $`=`$ $`()^{m+s}{}_{s}{}^{}Y_{\mathrm{},m}^{},`$ (16a)
$`\overline{ð}\overline{ð}\overline{ð}\overline{ð}\left[{}_{2}{}^{}\overline{Y}_{\mathrm{}m}^{}\right]`$ $`=`$ $`()^{m2}{\displaystyle \frac{(\mathrm{}+2)!}{(\mathrm{}2)!}}{}_{2}{}^{}Y_{\mathrm{},m}^{},`$ (16b)
$`ðððð\left[{}_{2}{}^{}\overline{Y}_{\mathrm{}m}^{}\right]`$ $`=`$ $`()^{m+2}{\displaystyle \frac{(\mathrm{}+2)!}{(\mathrm{}2)!}}{}_{2}{}^{}Y_{\mathrm{},m}^{},`$ (16c)
we obtain a first order relationship among $`\psi _4`$ or $`\psi _0`$ and the IRG or ORG potentials decomposed into multipoles

$`\left[\rho ^4\psi _4\right]^\pm `$ $`=`$ $`\pm {\displaystyle \frac{(\mathrm{}+2)!}{4(\mathrm{}2)!}}\mathrm{\Psi }_{IRG}^\pm 3M_t\mathrm{\Psi }_{IRG}^\pm ,`$ (17)
$`\psi _0^\pm `$ $`=`$ $`\pm {\displaystyle \frac{(\mathrm{}+2)!}{4(\mathrm{}2)!}}\mathrm{\Psi }_{ORG}^\pm +3M_t\mathrm{\Psi }_{ORG}^\pm ,`$ (18)
where we have used the notation

$`\psi ^\pm ={\displaystyle \frac{1}{2}}[\psi ^{\mathrm{},m}\pm ()^m\overline{\psi ^{\mathrm{},m}}]`$ (19)
for all fields decomposed into multipoles.

Since $`\mathrm{\Psi }_{IRG}`$ and $`\mathrm{\Psi }_{ORG}`$ satisfy the master equation (2) for spin $`s=2`$ respectively we can eliminate from this equation all time derivatives by replacing Eqs. (17) or (18) and its time derivatives into the Teukolsky equation. This leads to the following equation for the IRG / ORG potentials, both represented here by $`\mathrm{\Psi }`$

$`\mathrm{\Delta }^s_r[\mathrm{\Delta }^{s+1}_r\mathrm{\Psi }^\pm ]{\displaystyle \frac{r^4}{\mathrm{\Delta }}}(\mathrm{\Omega }_{AS})^2\mathrm{\Psi }^\pm `$
$`\pm 2sr(\mathrm{\Omega }_{AS})({\displaystyle \frac{Mr}{\mathrm{\Delta }}}1)\mathrm{\Psi }^\pm `$
$`(\mathrm{}s)(\mathrm{}+s+1)\mathrm{\Psi }^\pm =F^\pm `$ (20)
where

$$\mathrm{\Psi }=\{\begin{array}{cc}\mathrm{\Psi }_{IRG}\hfill & \mathrm{for}s=2\hfill \\ \mathrm{\Psi }_{ORG}\hfill & \mathrm{for}s=+2\hfill \end{array},$$
(21)
$$\mathrm{\Omega }_{AS}=\frac{1}{12M}\frac{(\mathrm{}+2)!}{(\mathrm{}2)!},$$
(22)
and the source term is

$`F^\pm `$ $`:={\displaystyle \frac{r^4}{3M\mathrm{\Delta }}}(_t\psi ^\pm ){\displaystyle \frac{r^4}{3M\mathrm{\Delta }}}(\mathrm{\Omega }_{AS})(\psi ^\pm )`$ (23)
$`+2s{\displaystyle \frac{r}{3M}}(Mr/\mathrm{\Delta }1)(\psi ^\pm ).`$
The equation (III.1) is our fundamental equation to solve for the potential in terms of the known fields $`\psi _4`$ or $`\psi _0`$ that appear in the source terms. The key observations here are to note that, separately for the “plus” and “minus” parts: i) the left hand side of this equation contains the Teukolsky operator in the frequency domain for the algebraically special frequencies, $`\omega =\pm i\mathrm{\Omega }_{AS}`$, and ii) the source terms in (23) are the Cauchy data for the Teukolsky operator in terms of $`\psi _4`$ or $`\psi _0`$, precisely as they would appear in a Laplace transform approach to Eq. (2) \[See Eqs. (A2-A3) of Ref. Campanelli98a \]. Notably, we have arrived to this equation working in the time domain, without any frequency decomposition.

### III.2 solution

The single frequency appearing on the right hand side of Eq. (2) is precisely an algebraic special frequency. Algebraically special perturbations of black holes excite gravitational waves which are either purely ingoing or outgoing. Hence only one of $`\psi _4`$ or $`\psi _0`$ is non-zero while the other vanishes. Chandrasekhar Chandrasekhar:1984  has obtained the explicit form of the algebraic special perturbations of the Kerr black hole. This is a remarkable fact, because we can then use these analytic solutions to the Teukolsky equation for algebraic special perturbations to construct explicit solutions to (III.1).

For Schwarzschild, Chandrasekhar Chandrasekhar:1984  gives two algebraic special solutions

$`y_1(r)`$ $`=`$ $`\left[1+{\displaystyle \frac{\lambda r}{M}}+{\displaystyle \frac{\lambda ^2r^2}{3M^2}}+{\displaystyle \frac{\lambda ^2(\lambda +1)r^3}{9M^3}}\right]e^{\mathrm{\Omega }_{AS}r^{}},`$ (24)
$`y_2(r)`$ $`=`$ $`\left[1{\displaystyle \frac{(\lambda +1)r}{3M}}\right]r^2e^{+\mathrm{\Omega }_{AS}r^{}},`$ (25)
where $`\lambda =(\mathrm{}1)(\mathrm{}+2)/2`$. Note that $`y_2`$ gives rise to an algebraically special solution for $`\psi _4`$, ($`s=2`$), at frequency $`\omega _{AS}^{}=i\mathrm{\Omega }_{AS}`$ and for $`\mathrm{\Delta }^2\psi _0`$, ($`s=+2`$), at frequency $`\omega _{AS}^+=+i\mathrm{\Omega }_{AS}`$, while $`y_1`$ gives rise to an algebraically special solution for $`\psi _4`$ at frequency $`\omega _{AS}^+`$ and for $`\mathrm{\Delta }^2\psi _0`$ at frequency $`\omega _{AS}^{}`$.

These two solutions satisfy the “plus” and “minus” parts of Eq. (III.1) for $`\mathrm{\Psi }_{IRG}^\pm `$ ($`s=+2`$) and vice-versa for $`\mathrm{\Psi }_{ORG}^{}`$ ($`s=2`$). A second set of independent solutions can be found:

$`z_1(r)`$ $`=`$ $`y_1(r){\displaystyle _{2M}^r}{\displaystyle \frac{W(r^{})}{y_1(r^{})^2}}𝑑r^{},`$ (26)
$`z_2(r)`$ $`=`$ $`y_2(r){\displaystyle _{\mathrm{}}^r}{\displaystyle \frac{W(r^{})}{y_2(r^{})^2}}𝑑r^{},`$ (27)
where the Wronskian of the solutions is

$`W(r):`$ $`=`$ $`y_{1,2}(r)_rz_{1,2}(r)z_{1,2}(r)_ry_{1,2}(r)`$ (28)
$`=`$ $`\mathrm{\Delta }(r)=r(r2M).`$
Note: $`z_1`$ and $`z_2`$ are not algebraically special solutions, although they are each a homogeneous solution to their respective Teukolsky equation at an algebraically special frequency.

Making use of these solutions, the explicitly expression for the potential can be written as follows:

$`\mathrm{\Psi }_{IRG}^+(r)`$ $`=`$ $`y_1(r){\displaystyle _{2M}^r}{\displaystyle \frac{z_1(r^{})F^+(r^{})}{\mathrm{\Delta }(r^{})^2}}𝑑r^{}`$ (29)
$``$ $`z_1(r){\displaystyle _r^{\mathrm{}}}{\displaystyle \frac{y_1(r^{})F^+(r^{})}{\mathrm{\Delta }(r^{})^2}}𝑑r^{},`$
$`\mathrm{\Psi }_{IRG}^{}(r)`$ $`=`$ $`z_2(r){\displaystyle _{2M}^r}{\displaystyle \frac{y_2(r^{})F^{}(r^{})}{\mathrm{\Delta }(r^{})^2}}𝑑r^{}`$ (30)
$`+`$ $`y_2(r){\displaystyle _r^{\mathrm{}}}{\displaystyle \frac{z_2(r^{})F^{}(r^{})}{\mathrm{\Delta }(r^{})^2}}𝑑r^{},`$
valid for each hypersurface where $`\psi _0`$ or $`\psi _4`$, entering in $`F`$ given by Eq. (III.1), are evaluated.

The equations for the ORG are obtained by exchanging “plus” and “minus” parts above (cf. Eq. (43) and (44)), and by adopting the corresponding dependence of $`F`$ on the components of $`\psi `$ \[See Eq. (3)\].

## IV Applications

### IV.1 Metric perturbations (IRG)

Using the Regge-Wheeler notation Regge57  for metric perturbations, conditions (II.1) read

$`h_0^{\mathrm{}}(r,t)^{(\mathrm{even},\mathrm{odd})}=(12M/r)h_1^{\mathrm{}}(r,t)^{(\mathrm{even},\mathrm{odd})},`$
$`G^{\mathrm{}}(r,t)={\displaystyle \frac{2}{\mathrm{}(\mathrm{}+1)}}K^{\mathrm{}}(r,t),`$
$`H_0^{\mathrm{}}(r,t)=H_2^{\mathrm{}}(r,t)=H_1^{\mathrm{}}(r,t).`$ (31)
We can now explicitly compute the metric perturbations (9) or (11) in terms of the computed $`\psi _0`$ or $`\psi _4`$:222For Newman-Penrose scalars, as in Eqs. (IV.1) – (36), even and odd really refer to contributions to the real and imaginary parts of the whole scalar, not specifically to the individual coefficients in a mode sum decomposition.

$`[h_{(n\overline{m})}]^{\mathrm{}}`$
$`=\left\{h_0^{\mathrm{}}(r,t)^{(\mathrm{even})}ih_0^{\mathrm{}}(r,t)^{(\mathrm{odd})}\right\}{\displaystyle \frac{\sqrt{\mathrm{}(\mathrm{}+1)}}{\sqrt{2}r}}{}_{1}{}^{}Y_{\mathrm{}m}^{}`$
$`=\{({\displaystyle \frac{y_1^{}}{ry_1}}+{\displaystyle \frac{C_{\mathrm{}}}{12M}}{\displaystyle \frac{1}{(12M/r)}}{\displaystyle \frac{2}{r^2}})[\mathrm{\Psi }_{IRG}^+]`$
$`{\displaystyle \frac{\mathrm{\Delta }}{ry_1}}{\displaystyle _r^{\mathrm{}}}{\displaystyle \frac{y_1}{\mathrm{\Delta }^2}}[F^+]𝑑r^{}{\displaystyle \frac{[\rho ^4\psi _4^+]}{3M(r2M)}}`$
$`+\left({\displaystyle \frac{y_2^{}}{ry_2}}{\displaystyle \frac{C_{\mathrm{}}}{12M}}{\displaystyle \frac{1}{(12M/r)}}{\displaystyle \frac{2}{r^2}}\right)[\mathrm{\Psi }_{IRG}^{}]`$
$`+{\displaystyle \frac{\mathrm{\Delta }}{ry_2}}{\displaystyle _{2M}^r}{\displaystyle \frac{y_2}{\mathrm{\Delta }^2}}[F^{}]dr^{}{\displaystyle \frac{[\rho ^4\psi _4^{}]}{3M(r2M)}}\}`$
$`\times \sqrt{2(\mathrm{}+2)(\mathrm{}1)}{}_{1}{}^{}Y_{\mathrm{}m}^{},`$ (32)
where $`C_{\mathrm{}}=(l+2)!/(l2)!`$, and

$`[h_{nn}]^{\mathrm{}}`$ $`=`$ $`(12M/r)H_0^{\mathrm{}}(r,t){}_{0}{}^{}Y_{\mathrm{}m}^{}`$ (33)
$`=`$ $`{\displaystyle \frac{\rho ^2}{2}}[\mathrm{\Psi }_{IRG}^+]\sqrt{{\displaystyle \frac{(\mathrm{}+2)!}{(\mathrm{}2)!}}}{}_{0}{}^{}Y_{\mathrm{}m}^{},`$
$`[h_{\overline{m}\overline{m}}]^{\mathrm{}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left\{G(r,t)^{\mathrm{}}+i{\displaystyle \frac{h_2(r,t)^{\mathrm{}}}{r^2}}\right\}\sqrt{{\displaystyle \frac{(\mathrm{}+2)!}{(\mathrm{}2)!}}}{}_{2}{}^{}Y_{\mathrm{}m}^{}`$ (34)
$`=`$ $`[h_{\overline{m}\overline{m}}]^{(\mathrm{even})}+[h_{\overline{m}\overline{m}}]^{(\mathrm{odd})}`$
where

$`[h_{\overline{m}\overline{m}}]^{(\mathrm{even})}=\{1/6{\displaystyle \frac{[\mathrm{\Psi }_{IRG}^+]\left(12M^2+r^2C_l\right)y_1^{}}{y_1r\left(r2M\right)M}}`$
$`1/6{\displaystyle \frac{\left(\mathrm{\Delta }\right)\left(12M^2+r^2C_l\right)}{y_1r\left(r2M\right)M}}{\displaystyle _r^{\mathrm{}}}{\displaystyle \frac{y_1}{\mathrm{\Delta }^2}}[F^+]𝑑r^{}`$
$`2/3{\displaystyle \frac{r^2[\rho ^4_t\psi _4^+]}{\left(r2M\right)^2M}}2/3{\displaystyle \frac{r\frac{}{r}[\rho ^4\psi _4^+]}{\left(r2M\right)M}}`$
$`1/18{\displaystyle \frac{[\rho ^4\psi _4^+]\left(36r+84M+r^2C_l\right)}{\left(r2M\right)^2M}}`$
$`+{\displaystyle \frac{1}{72}}(r^3C_{l}^{}{}_{}{}^{2}+12r(6l^2+6l+7C_l12)M^2`$
$`144(l+2)(l1)M^336r^2C_lM)`$
$`\times {\displaystyle \frac{[\mathrm{\Psi }_{\mathrm{𝐼𝑅𝐺}}^+]}{\left(r2M\right)^2M^2r}}\}_2Y_\mathrm{}m,`$ (35)
and

$`[h_{\overline{m}\overline{m}}]^{(\mathrm{odd})}=\{1/6{\displaystyle \frac{[\mathrm{\Psi }_{\mathrm{𝐼𝑅𝐺}}^{}]\left(12M^2r^2C_l\right)y_2^{}}{y_2r\left(r2M\right)M}}`$ (36)
$`+1/6{\displaystyle \frac{\left(\mathrm{\Delta }\right)\left(12M^2r^2C_l\right)}{y_2r\left(r2M\right)M}}{\displaystyle _{2M}^r}{\displaystyle \frac{y_2}{\mathrm{\Delta }^2}}[F^{}]𝑑r^{}`$
$`2/3{\displaystyle \frac{r^2[\rho ^4_t\psi _4^{}]}{\left(r2M\right)^2M}}2/3{\displaystyle \frac{r\frac{}{r}[\rho ^4\psi _4^{}]}{\left(r2M\right)M}}`$
$`1/18{\displaystyle \frac{[\rho ^4\psi _4^{}]\left(36r+84Mr^2C_l\right)}{\left(r2M\right)^2M}}`$
$`+{\displaystyle \frac{1}{72}}(r^3C_{l}^{}{}_{}{}^{2}+12r(6l^2+6l7C_l12)M^2`$
$`144(l+2)(l1)M^3+36r^2C_lM)`$
$`\times {\displaystyle \frac{[\mathrm{\Psi }_{IRG}^{}]}{\left(r2M\right)^2M^2r}}\}_2Y_\mathrm{}m.`$
In writing these we have explicitly lowered the order of the derivatives of the potential by making use of the following identities

$`_r\left([\mathrm{\Psi }^+]\right)=`$
$`\left({\displaystyle \frac{_ry_1}{y_1}}\right)[\mathrm{\Psi }^+]{\displaystyle \frac{\mathrm{\Delta }}{y_1}}{\displaystyle _r^{\mathrm{}}}\left({\displaystyle \frac{y_1}{\mathrm{\Delta }^2}}\right)[F^+]𝑑r^{},`$ (37)
$`_r\left([\mathrm{\Psi }^{}]\right)=`$
$`\left({\displaystyle \frac{_ry_2}{y_2}}\right)[\mathrm{\Psi }^{}]+{\displaystyle \frac{\mathrm{\Delta }}{y_2}}{\displaystyle _{2M}^r}\left({\displaystyle \frac{y_2}{\mathrm{\Delta }^2}}\right)[F^{}]𝑑r^{},`$ (38)
and directly from Eq. (17)

$`_t\mathrm{\Psi }_{IRG}^\pm =\pm {\displaystyle \frac{(\mathrm{}+2)!}{12M(\mathrm{}2)!}}\mathrm{\Psi }_{IRG}^\pm {\displaystyle \frac{1}{3M\rho ^4}}\psi _4^\pm .`$ (39)
### IV.2 Metric perturbations (ORG)

Using the Regge-Wheeler notation Regge57  for metric perturbations conditions (10) read

$`h_0^{\mathrm{}}(r,t)^{(\mathrm{even},\mathrm{odd})}=(12M/r)h_1^{\mathrm{}}(r,t)^{(\mathrm{even},\mathrm{odd})},`$
$`G^{\mathrm{}}(r,t)={\displaystyle \frac{2}{\mathrm{}(\mathrm{}+1)}}K^{\mathrm{}}(r,t),`$
$`H_0^{\mathrm{}}(r,t)=H_2^{\mathrm{}}(r,t)=H_1^{\mathrm{}}(r,t).`$ (40)
The explicit metric perturbations can be found directly from the previous subsection, Eqs. (IV.1)-(36) upon exchanges of the tetrad contractions $`ln`$ and $`m\overline{m}`$, and the consequent change in normalizations, $`s`$, and potentials, as described throughout the paper.

## V Discussion

We have explicitly computed the metric perturbations of a nonrotating black hole in terms of the Weyl scalars $`\psi _4`$ and $`\psi _0`$ which can be computed directly by solving the Teukolsky equation for any appropriate astrophysical scenario, given the corresponding initial data. In doing so we had to invert Eqs. (13) or (15). This was performed making explicit use of the multipole decomposition of the potential, Weyl scalars and metric. The extension of this procedure to the rotating background is not straightforward, but we have learned some key features: The algebraic special solutions of the Teukolsky equation will play a crucial role in finding the solution for the Hertz potential in terms of $`\psi _4`$ or $`\psi _0`$. One can see this as follows. In order to invert Eq. (13) for the potential we first seek out solutions of the homogeneous equation. Hence $`\psi _4`$ should vanish, but for the solution to be nontrivial, $`\psi _0`$, given by Eq. (12), must not vanish Wald:1973 . These two conditions are precisely the conditions that define the algebraically special solutions for the potential satisfying the vacuum Teukolsky equation (2). These two solutions (in the time domain) should allow one to build up the Kernel that inverts the fourth order equation (13). An identical argument applies for the ORG potential when working with Eqs. (14) and (15).

One main application of this formalism is to go beyond first order perturbation theory and compute second order perturbations of rotating black holes. In Ref. Campanelli99  there was developed a formalism for the second order Teukolsky equation that takes the form of the first order wave operator acting on the second order piece of the Weyl scalar $`\psi _4`$, and a source term build up as a quadratic combination of first order perturbations. It is precisely to compute this source term that one needs the explicit form of the metric perturbations. In Ref. Campanelli99  it was found that to describe the emitted gravitational radiation the ORG gauge is specially well suited. Hence one has to solve the the first order Teukolsky equation for $`\psi _4^{(1)}`$ and $`\psi _0^{(1)}`$, then later to build up the source of the second order piece of the Weyl scalar $`\psi _4^{(2)}`$ (See Eq. (9) in Ref. Campanelli99 ).

A second important application of the reconstruction of metric perturbations around Kerr background is to compute the self force of a particle orbiting a massive black hole Mino:1997nk ; Quinn:1997am  and to compute corrected trajectories Lousto99b  depending on the perturbed metric and connection coefficients along the particle world line. This task is left for a future paper. While we know the form of the Teukolsky equation in the presence of of perturbative matter around a Kerr hole (see Eq. (2)), we need to generalize the equation satisfied by the potential and the relationship between this potential and the metric perturbations, i.e. the generalization of Eqs. (9) and (11). In particular we know that not all of conditions (II.1) or (10) can hold in the presence of matter since they are then incompatible with the Einstein equations Barack:2001ph .

###### Acknowledgements.

We wish to thank H. Beyer, M. Campanelli and S. Detweiler for helpful discussions. This research has been supported in part by NSF Grant No. PHY-9800977 (B.F.W.) with the University of Florida. C.O.L developed part of this work in the Albert Einstein Institut, Germany, which B.F.W. also thanks for hospitality.

## Appendix A Relation between IRG/ORG potentials

We give here a relation expressing the potential $`\mathrm{\Psi }_{ORG}`$ in terms of the result for $`\mathrm{\Psi }_{IRG}`$. This was not given in Chrzanowski:1975wv ; Cohen:1974 ; Wald:1978vm ; Kegeles:1979an . We begin by defining a field $`\chi `$, with spin weight $`2`$, through:

$$\mathrm{\Psi }_{IRG}=\frac{1}{4}\left[\overline{}\overline{}\overline{}\overline{}\overline{\chi }+12\rho ^3\psi _2_t\chi \right],$$
(41)
analogous to Eq. (15). The solution for $`\mathrm{\Psi }_{ORG}`$ is:

$$\mathrm{\Psi }_{ORG}=DDDD\left[\overline{\chi }\right],$$
(42)
also the relation (12) between $`\psi _0`$ and $`\mathrm{\Psi }_{IRG}`$. A potential degeneracy for algebraically special modes was discussed briefly in Whiting:1989vc . The explicit solution for the modes $`\chi ^\mathrm{}m`$ can be written, in terms of $`\chi ^\pm `$ given by Eq. (19), as:

$`\chi ^+(r)`$ $`=`$ $`z_2(r){\displaystyle _{2M}^r}{\displaystyle \frac{y_2(r^{})V^+(r^{})}{\mathrm{\Delta }(r^{})^2}}𝑑r^{}`$ (43)
$`+`$ $`y_2(r){\displaystyle _r^{\mathrm{}}}{\displaystyle \frac{z_2(r^{})V^+(r^{})}{\mathrm{\Delta }(r^{})^2}}𝑑r^{},`$
$`\chi ^{}(r)`$ $`=`$ $`y_1(r){\displaystyle _{2M}^r}{\displaystyle \frac{z_1(r^{})V^{}(r^{})}{\mathrm{\Delta }(r^{})^2}}𝑑r^{}`$ (44)
$``$ $`z_1(r){\displaystyle _r^{\mathrm{}}}{\displaystyle \frac{y_1(r^{})V^{}(r^{})}{\mathrm{\Delta }(r^{})^2}}𝑑r^{},`$
cf. Eqs. (29) and (30) above. Here $`V^\pm `$ are given by:

$`V^\pm `$ $`:={\displaystyle \frac{r^4}{3M\mathrm{\Delta }}}(_t\mathrm{\Psi }_{IRG}^\pm ){\displaystyle \frac{r^4}{3M\mathrm{\Delta }}}(\mathrm{\Omega }_{AS})(\mathrm{\Psi }_{IRG}^\pm )`$ (45)
$`+4{\displaystyle \frac{r}{3M}}(Mr/\mathrm{\Delta }1)(\mathrm{\Psi }_{IRG}^\pm ),`$
similar to Eq. (23), but note the sign differences.


format: text

text:
# Parallel J-W Monte Carlo Simulations of Thermal Phase Changes in Finite-size Systems

## 1 Introduction

The interest in phase transitions in finite systems is two-fold: firstly, bulk properties of the material can be simulated if the system is studied under periodic boundary conditions  and proper account for rounding and shifting of the measurable quantities is taken in the analysis. Secondly, free finite systems, such as molecular clusters, are of interest due to their peculiar properties, which are not observed in bulk systems of the same substance . These are important for the synthesis of nanomaterials and nanodevices .

We use a canonical Monte Carlo (MC) method to investigate temperature-induced critical behavior in clusters containing 59 $`TeF_6`$ molecules.

Our model system contains 59 octahedral rigid $`TeF_6`$ molecules which are allowed to rotate and translate. Their interaction is reliably described by a Lennard-Jones and a Coulomb atom-atom potential:

$`V(𝐪)`$ $`=`$ $`{\displaystyle \underset{i<j}{}}U(r_{ij})`$ (1)
$`V(𝐪)`$ $`=`$ $`{\displaystyle \underset{i<j}{}}{\displaystyle \underset{\alpha ,\beta }{}}4ϵ_{\alpha \beta }\left[\left({\displaystyle \frac{\sigma _{\alpha \beta }}{r_{ij}}}\right)^{12}\left({\displaystyle \frac{\sigma _{\alpha \beta }}{r_{ij}}}\right)^6\right]+`$
$`+`$ $`{\displaystyle \underset{i<j}{}}{\displaystyle \underset{\alpha ,\beta }{}}\left({\displaystyle \frac{q_{i\alpha }q_{i\beta }}{r_{ij}}}\right)`$
where $`𝐪`$ is a generalized coordinate; $`r_{ij}`$ is the distance between the $`ith`$ and $`jth`$ atom. The indices $`\alpha ,\beta `$ denote either a fluorine or a tellurium atom. The parameter values are taken from .

This orientation-dependent potential generates rugged potential energy surfaces. How to tackle the behavior of a system having a rugged potential energy surface in MD simulations is discussed in .

Here we tackle the problem with the help of the jump-walking Monte Carlo method, , which is a smart development of the Metropolis algorithm. We have parallelized the code to speed up the computations .

## 2 J-walking algorithm

In its general form, the method generates trial moves from a higher-temperature $`(T_J)`$ equilibrium distribution with a probability specified by the variable $`P_J`$:

$$P=min\{1,q[x^{}|x]\}$$
(2)
where the Metropolis probability function is:

$$q(x^{}|x)=e^{(\frac{1}{T}\frac{1}{T_J})\left\{V(x^{})V(x)\right\}}$$
(3)
$`V(x)`$ is the potential energy of the system at a temperature $`T`$ and $`V(x^{})`$ is the potential energy of a new configuration at the temperature $`T_J`$.

The remaining trial moves are of conventional Metropolis character.

In our study we use multi-temperature generalizations  of the basic approach.

We develop a parallel J-walking code that enables us to carry out a Monte Carlo simulation efficiently in a multiprocessor computing environment. We apply the code in the study of the thermodynamic properties of $`(TeF_6)_n`$, $`n=59`$ clusters.

The J-walking technique can be implemented in two ways. The first approach is to write the configurations from the simulation at the J-walking temperature to an external file and access these configurations randomly, while carrying out a simulation at the lower temperature. It is necessary to access the external files randomly to avoid correlation errors . The large storage requirements limit the application of the method to small systems. The second approach uses tandem walkers, one at a high temperature where Metropolis sampling is ergodic and multiple walkers at lower temperatures.

The best features of these two approaches can be combined into a single J-walking algorithm  with the use of multiple processors and the Message Passing Interface (MPI) library. We incorporate MPI functions into the MC code to send and receive configuration geometries and potential energies of the clusters. Instead of generating external distributions and storing them before the actual simulation, we generate the required distributions during simulation and pass them to the lower-temperature walkers.

Parallel J-walking algorithm:

* Step 1. For each t make Metropolis MC steps:

+ Rotate each molecule.
    + Translate each molecule.
    + Reject or accept step.
    + Go to step 1.
* Step 2. After $`S_1`$ steps - collect statistics:

+ Potential energy histogram.
    + Energy average and deviation.
    + Heat capacity $`C_V`$.
    + Save current configuration.
* Step 3. After $`S_2`$ steps - make jump-walking step by exchanging the configurations using MPI.
* Step 4. Go to step 1.

The diagram of model J-walking is shown in fig.1. Each square is a Metropolis MC simulation at a particular temperature. The set of boxes on the right-hand side represents the array of previous configurations of the system, which are stored in the memory to avoid correlations between the lower and the higher temperature. At each trial jump we randomly choose one of the 4 systems. When a configuration is transmitted to a lower-temperature process, it is a configuration randomly chosen from the array of higher-temperature walkers. The current configuration of the walker then replaces the configuration just passed from the array to another temperature. In fig.2 we show the parallel decomposition of computation. Each process computes part of one of 4 multistage J-walking chains, and exchanges configurations and energy with others.

We use array sizes of 2500 configurations. The number of configurations is limited by the processor RAM used in simulation. In the parallelization implemented in our code, the arrays are small and do not inhibit applications of the method to large systems. For the computations in the present work the number of MC passes for each temperature is $`6.10^5`$ for each cluster containing 59 $`TeF_6`$ molecules (413 atoms). We make J-walk jump attempts at every 50 Metropolis MC steps during the thermalisation and at every 150 steps during the main computation.

The computer code has been ported, tested and optimized on SUN 3500 and CRAY-T3E machines. In our program dynamic memory management has been implemented for optimal usage of the memory. We find that memory and performance requirement make CRAY-T3E more suitable for such computations. In our runs we use 64 processors each with 64 MB RAM. Each run of $`6.10^5`$ steps takes approximately 11h per CPU.

## 3 Results and Conclusions

Using the parallel code described in the previous section we make sets of production runs for 59 molecule clusters. The previous MD analysis of the temperature behavior of $`TeF_6`$ clusters pointed out a two-step structural transformation process  from an orientationally disordered bcc structure below cluster solidification to an orientationally oriented ordered monoclinic structure detected at low $`T`$ (below $`20K`$). The first step involves lattice reconstruction (bcc to monoclinic) and a partial order of one of the molecular axes, when the cluster is cooled down to its freezing point. This transition has been proved to be a first-order phase change  by detecting coexisting phases. A further temperature decrease causes complete orientational order of the three molecular axes. This transformation is continuous. The diagnostic method developed in  has been implemented to animate the solid-solid transformations .

To distinguish between the different phases, which the clusters adopt at different temperatures, we have computed the Lindemann index :

$$\delta _{lin}=\frac{2}{N(N1)}\underset{i,j(>i)=1}{\overset{N}{}}\frac{\sqrt{<r_{ij}^2><r_{ij}>^2}}{<r_{ij}>}$$
(4)
where $`r_{ij}`$ is the distance between the centers of the molecules ($`Te`$ atoms). For the values of $`\delta _{lin}>0.08`$ the system is in liquid phase and for $`\delta _{lin}<0.08`$ is in the solid phase. Fig.3 presents the Lindemann criterion for the 4 systems. The solidification occurs in interval $`100K115K`$. In the interval $`75K80K`$ the angle of the curve has changed.

We calculate the radial distribution function to reconstruct the crystal lattice. Fig.4 and fig.5 show the normalized radial distribution function for distances between centers of the $`TeF_6`$ molecules ($`Te`$ atoms). Fig.4 is at $`157`$. The peaks of the curve correspond to disordered liquid structure. Fig.5 is at $`76K`$ where the peaks correspond to an ordered structure (bcc lattice).

In this study a J-walking Monte Carlo algorithm has been implemented to study $`(TeF_6)`$ clusters. A parallel J-walking code has been developed. Thermodynamic properties have been observed in the range from $`60K`$ to $`150K`$. We can compare the results to the results obtained from MD simulations.

Acknowledgments

The work was supported under the EC (contract No. HPRI-CT-1999-00026) TRACS-EPCC.