text:
# The violent past of Cygnus X–2

## 1 Introduction

The $`9.84`$ d period neutron–star binary Cygnus X–2 has long been regarded as an archetypal long–period low–mass X–ray binary (Cowley et al. 1979, Webbink et al. 1983) where nuclear expansion of a Hayashi line giant drives mass transfer. Yet recent optical photometry and high–resolution spectroscopy, while confirming that the donor has a low mass (Casares et al. 1998; Orosz & Kuulkers 1999), unambiguously showed that the spectral type of the optical counterpart is A$`5\pm 2`$ (Casares et al. 1998), significantly too hot for a Hayashi line donor.

King & Ritter (1999; hereafter KR) and Podsiadlowski & Rappaport (2000; PR) argued that the system must be the descendent of an intermediate–mass X–ray binary (IMXB), with $`34\mathrm{M}_{}`$ as the likely initial donor mass. Such systems undergo a rapid mass transfer phase, previously regarded as fatal, with transfer rates exceeding the Eddington value by several orders of magnitude. For a narrow range of initial separations the systems never reach the Hayashi line, or evolve away from it, during the subsequent phase with slower mass transfer.

KR suggested that Cygnus X–2 is the product of a “case B” mass transfer sequence, where the donor star was already expanding towards the giant branch when mass transfer began (Kippenhahn & Weigert ). By contrast, PR preferred an evolution where mass transfer started while the donor was still on the main sequence, with core hydrogen burning terminating during the mass transfer phase. KR’s considerations are semi–analytical and entirely based on generalised main–sequences (Giannone et al. 1968), while PR performed detailed binary sequences with a full stellar code. In Sect. 2 and 3 we use our evolutionary code to reexamine critically the case B solution rejected by PR, and resolve the discrepancy between KR and PR. To see how common a Cygnus X–2–like evolution is we consider IMXBs with still higher initial donor masses in Sect. 4.

A further peculiarity of Cygnus X–2, addressed in Sect. 5, is its dynamical state. The system has a measured line-of-sight velocity of -208.6 km/s (Casares et al. 1998) and a Galactic latitude of -11.32 and longitude of 87.33. At a distance of 11.6 kpc from the sun (Smale 1998), this places Cygnus X–2 at a Galactocentric distance of 14.2 kpc, and a distance from the Galactic plane of 2.28 kpc. Integrating the equations of motion in the Galactic potential we use Monte Carlo techniques to investigate possible trajectories of Cygnus X–2.

## 2 Constraining the prehistory of Cygnus X–2

The observed location of the donor star of Cygnus X–2 in the middle of the Hertzsprung gap in the HR diagram, and its small mass ($`0.40.7\mathrm{M}_{}`$), suggest that the mass of the hydrogen–rich envelope remaining above the donor’s helium core is very small ($`M_\mathrm{H}/M_{\mathrm{He}}0.05`$). Stars with an initial mass $`34.5\mathrm{M}_{}`$ leave the core hydrogen burning phase with a helium core mass $`M_{\mathrm{He}}`$ in this range (e.g. Bressan et al. 1993). Mass transfer from such an intermediate–mass star on to a less massive neutron star is thermally unstable and involves an initial phase of rapid mass transfer with a highly super–Eddington transfer rate. This rapid phase lasts roughly until the mass ratio is reversed and the donor’s Roche lobe expands upon further mass transfer. The subsequent slower transfer phase proceeds on the donor’s thermal timescale, which is essentially given by the Kelvin–Helmholtz time of the donor when it left the main sequence.

PR identified two possible routes from the post–supernova binary, i.e. after the formation of the neutron star, to the present system configuration. One route is via a genuine case B evolution, the other via an evolution they labelled ‘case AB’. In the former the donor star has already left the main sequence when it fills its Roche lobe for the first time. The mass transfer phase is short–lived, with most of the mass transferred at a highly super–Eddington rate. Hence the neutron star mass increase is negligible, even if it accretes at the Eddington rate during the whole evolution. In the case AB solution discussed by PR mass transfer already starts during the donor’s main–sequence phase. The system briefly detaches when core hydrogen burning terminates, but mass transfer resumes when shell–burning is well established. As the rapid mass transfer phase terminates before the donor leaves the main sequence the subsequent slow phase proceeds on a much longer timescale than in the genuine case B solution — the thermal time of the now less massive donor at the terminal main–sequence. The transfer rate is sub–Eddington for some time and the neutron star grows in mass.

Given this we explore three different prescriptions to project the evolution of Cygnus X–2 backwards in time:

(Model 1) To represent a case B evolution we assume that $`M_1=`$const., and that any material lost from the system carries the specific orbital angular momentum of the neutron star. Then the initial orbital separation $`a_i`$ is

$$a_i=a\frac{M_1+M_2}{M_1+M_{2i}}\left(\frac{M_2}{M_{2i}}\right)^2\mathrm{exp}\left(\frac{2(M_{2i}M_2)}{M_1}\right),$$
(1)
where $`M_{2i}`$, $`M_2`$ denote the initial and present donor mass, and $`a`$ the present separation (e.g. KR). This solution is valid only if $`a_i>a_{\mathrm{TMS}}`$, the orbital separation of a binary with a Roche–lobe filling terminal main–sequence (TMS) star of mass $`M_{2,i}`$.

(Model 2) The same applies to a case AB evolution, except that we require $`a_i<a_{\mathrm{TMS}}`$.

(Model 3) To account for a possible prolonged phase with sub–Eddington mass transfer in a case AB sequence we allow mass transfer to be conservative (total binary mass and orbital angular momentum is constant) for the last part of the evolution. We assume that the neutron star mass at birth was $`1.4\mathrm{M}_{}`$, and that the present neutron star mass is $`>1.4\mathrm{M}_{}`$. Hence the evolution backwards in time consists of two branches. During the conservative phase, characterised by $`aM_1^2M_2^2=`$constant, $`M_1`$ reduces from its present value to $`1.4\mathrm{M}_{}`$. The second phase is calculated with constant neutron star mass ($`1.4\mathrm{M}_{}`$), as in models 1 and 2 above.

The present system parameters of Cygnus X–2 are: orbital period $`P_0=9.844`$ d, mass ratio $`q=0.34\pm 0.04`$ (Casares et al. 1998), and neutron star mass $`M_1=1.41.8\mathrm{M}_{}`$. This mass range accommodates the canonical neutron star mass at birth, as well as the estimate by Orosz & Kuulkers (1999) based on modelling the observed ellipsoidal variations. If we adopt a specific value $`M_1`$ for the present neutron star mass, then the present donor mass is constrained to the range $`0.30M_1M_20.38M_1`$ by the observed mass ratio. This translates into a range of allowed initial separations, as shown by the hatched regions in Fig. 1, for case B (model 1, wide spacing) and case AB (model 2, narrow spacing) evolution. The parameter space for case AB solutions with model 3 assumptions is tiny and hence not shown in the figure. For conservative mass transfer (model 3) the orbital separation of Cygnus X–2 today increases less steeply with decreasing donor mass than for the isotropic wind case (model 2). Hence the separation for model 3 sequences describing the past evolution of Cygnus X–2 is generally larger than for the corresponding model 2 sequences. In particular, model 3 solutions require larger initial separations. Most of them are in conflict with the limit $`a_i<a_{\mathrm{TMS}}`$, i.e. inconsistent with the assumption that the donor was on the main sequence when mass transfer started.

More generally, if the mass lost from the system in sequences with $`M_1=`$ const. carries more (less) specific angular momentum than that of the neutron star, the required initial separation is larger (smaller) than estimated by (1). A somewhat higher loss seems more likely than a smaller loss, hence this reduces the parameter space available for case AB solutions.

Even if Fig. 1 indicates a viable solution it is not clear if the corresponding evolutionary sequence reproduces Cygnus X–2. Although the donor will have the observed mass (and hence radius) at the observed orbital period $`P_0`$, the effective temperature is of course not constrained by the above considerations. To check this we need a detailed binary sequence with full stellar models, as presented in the next section. In particular, it appears that the case AB sequences described by PR require the donor star to be already fairly close to the end of core hydrogen burning at turn–on of mass transfer. This is likely to narrow the parameter space for case AB sequences even further.

We note that the limits shown in Fig. 1 depend via $`a_{\mathrm{TMS}}`$ somewhat on the stellar input physics, in particular on the amount of convective core overshooting (cf. the discussion in the next section).

## 3 Model calculations

We calculated several detailed early massive case B binary evolution sequences, using Mazzitelli’s stellar code in its 1989 version (see Mazzitelli 1989, and references therein, for a summary of the input physics). Mass transfer was treated as in Kolb (1998). Table 1 summarises initial and final system parameters. The neutron star mass was $`1.4\mathrm{M}_{}=`$ const. in each case, with angular momentum loss treated as in model 1 above. Specifically, we present in detail the evolution of a system with initial donor mass $`M_2=3.5\mathrm{M}_{}`$ and orbital separation $`a_i=13.3R_{}`$ (Sequence I) and $`11.1R_{}`$ (sequence II).

The figures 2-4 confirm the general behaviour of case B evolutionary sequences as described above. Sequence I gets very close to the observed state of Cygnus X–2, although it fails to fit all observed parameters simultaneously with high accuracy. The evolutionary track in the HR diagram passes through the Cygnus X–2 error box, but at a period slightly longer than the observed value $`P_0=9.844`$ d. This can be seen in Fig. 4, which shows that at $`P_0`$ the sequence I donor is somewhat too cool, while in sequence II it is somewhat too hot. Hence the initial separation for a simultaneous fit of $`P`$, $`T_{\mathrm{eff}}`$ (and $`L`$) lies between $`11.1R_{}`$ and $`13.3R_{}`$. The transfer rate at $`P_0`$ for sequence I is $`5\times 10^8\mathrm{M}_{}\mathrm{yr}^1`$, a factor $`23`$ higher than the Eddington rate.

Observational estimates place the actual accretion rate in Cygnus X–2 close to the Eddington limit. With no evidence of significant mass loss from the system we expect that the transfer rate is also of the order of the Eddington rate. As the transfer rate in the model sequence is set by the thermal time of the progenitor star, a somewhat smaller initial secondary mass (e.g. $`3.0\mathrm{M}_{}`$) would give a lower rate; this has already been noted by PR.

We note the following main differences between our sequence I and the case B sequence calculated by PR:
(1) PR chose $`21.3R_{}`$ as the initial separation. This leads to the combination $`1.4+0.66\mathrm{M}_{}`$ for the component masses at $`P=P_0`$, i.e. to a mass ratio $`0.47`$, significantly larger than observed. Hence the PR case B sequence represents a poor fit for Cygnus X–2.
(2) More importantly, the mass transfer rate in the slow phase of sequence I decreases from $`10^7\mathrm{M}_{}\mathrm{yr}^1`$ to $`10^8\mathrm{M}_{}\mathrm{yr}^1`$ over a period of 2.5 Myr, while in PR’s model the transfer rate is always larger than $`10^7\mathrm{M}_{}\mathrm{yr}^1`$ (which is inconsistent with observations). A similar calculation by Tauris et al. (2000) with essentially the same stellar code as the one used by PR gives a similarly high rate in the slow phase.

The initial separation for our sequence II is essentially the same as the one used by PR for their case AB evolution ($`11.5R_{}`$, corresponding to a stellar radius $`5.6R_{}`$ of the donor at turn–on of mass transfer). The most important difference between the two calculations is the degree of convective overshooting assumed during the main–sequence phase: none in our models, a very strong one in PR’s models. Hence PR’s main–sequence band is much wider, with $`6.6R_{}`$ as the maximum radius $`R_{2,\mathrm{max}}`$ that a $`3.5\mathrm{M}_{}`$ star reaches during core hydrogen burning, compared to $`3.6R_{}`$ in our case. A moderate extent of core overshooting is favoured in the literature, giving $`R_{2,\mathrm{max}}=5.0R_{}`$ (Schaller et al. 1992), $`5.7R_{}`$ (Bressan et al. 1993) and $`4.8R_{}`$ (Dominguez et al. 1999) for a $`3.5\mathrm{M}_{}`$ star with solar composition. Schröder et al. (1997) — who use the same code and input physics as PR — prefer a rather efficient overshooting leading to values up to $`R_{2,\mathrm{max}}=7.1R_{}`$, while in a subsequent paper (Pols et al. 1997) the same authors concluded that $`R_{2,\mathrm{max}}=5.9R_{}`$ (which corresponds to their overshooting parameter $`\delta =0.12`$) provides the best overall fit to observations. This is only barely larger than the donor’s radius at turn–on of mass transfer in PR’s case AB solution. Hence standard assumptions about overshooting favour a case B solution for Cygnus X–2 over a case AB solution.

We conclude that a case B solution is a viable fit for the evolutionary history of Cygnus X–2. Given the parameter space limitations for a case AB evolution it does appear as the more likely solution for Cygnus X–2, even though this implies that the presently observed state of Cygnus X–2 is rather short–lived, of order several million years.

PR pointed out that the surface composition of the donor in Cygnus X–2 would be significantly hydrogen–depleted ($`X0.1`$) and show signs of CNO–processing if it had undergone a case AB evolution. Unfortunately, this does not distinguish unambiguously between a case AB and case B evolution, as the same is in principle true for a donor that had undergone a case B mass transfer. The models in our sequence I close to the position of Cygnus X–2 have a surface hydogen mass fraction of $`X=0.29`$, while C/N and O/N are close to the equilibrium values for CNO burning at $`2\times 10^7`$ K.

## 4 Sequences with initially more massive donor stars

The calculations discussed so far represent early case B mass transfer solutions for systems with a moderate initial mass ratio, $`q=M_2/M_12.5`$. Systems with smaller mass ratio ($`q1`$) have no initial rapid mass transfer phase and can be understood semi–analytically (see e.g. Kolb 1998, Ritter 1999), while the fate of systems with yet larger initial mass ratio $`3`$ is unclear. Hjellming (1989) pointed out that sustained thermal–timescale mass transfer can lead to a delayed transition to dynamical–timescale mass transfer. The reason for this is that the adiabatic mass–radius index of initially radiative stars decreases significantly when the radiative envelope is stripped rapidly. This causes the Roche lobe to shrink faster than the star. Estimates for the critical initial mass ratio which just avoids the delayed dynamical instability give $`q34`$ (Hjellming 1989; Kalogera & Webbink 1996), although none of these are based on self–consistent mass transfer calculations.

This critical value is relevant for identifying possible descendants of neutron–star systems undergoing early massive case B mass transfer. These should appear as binary millisecond pulsars lying significantly below the orbital period–white dwarf mass relation found by Rappaport et al. (1995; see also Tauris & Savonije 1999 for calculations with updated input physics) for systems descending from Hayashi line low–mass X–ray binaries. Obvious candidates are systems with a fairly high–mass white dwarf but short orbital period. The most discrepant system is B0655+64 ($`P=1.03`$ d, $`M_{\mathrm{WD}}0.7`$; cf. KR). The fairly massive white dwarf implies a donor mass $`5\mathrm{M}_{}`$ in the progenitor binary.

We performed test calculations similar to our sequences I and II, but with initial donor mass $`3.75\mathrm{M}_{}`$, $`4.0\mathrm{M}_{}`$ and $`5.0\mathrm{M}_{}`$. While the $`3.75\mathrm{M}_{}`$ sequence was stable throughout, the $`4.0\mathrm{M}_{}`$ and $`5.0\mathrm{M}_{}`$ sequences encountered runaway mass transfer (where we stopped the calculations) rather early in the rapid mass transfer phase (see Tab.1). An additional test sequence with a $`3\mathrm{M}_{}`$ donor near the end of core hydrogen burning and initial mass ratio 4 (i.e. primary mass $`0.75\mathrm{M}_{}`$, assumed constant, as above) encountered runaway mass transfer at a donor mass $`2.6\mathrm{M}_{}`$. This is in perfect agreement with Hjellming’s prediction, based on his Fig. IV.1 and a Roche lobe curve corresponding to $`q=4`$.

Unless a neutron star binary can survive even such a dynamical–timescale mass transfer, the phase space available for a Cygnus X–2–like evolution is severely limited by this upper limit on the initial donor mass. Our calculations suggest that if most neutron stars form with $`1.4\mathrm{M}_{}`$ the maximum white dwarf mass in an endproduct of early massive case B evolution is $`0.55\mathrm{M}_{}`$. If this is true, neither B0655+64 nor the recently discovered system J1453-58 ($`P=12.422`$ d, $`M_{\mathrm{WD}}0.88\mathrm{M}_{}`$; cf. Manchester et al. 1999) can have formed in this way.

The occurrence of the delayed dynamical instability is intimately linked to the internal structure of the donor star. Therefore it is not surprising that the maximum initial donor mass for early case B mass transfer, just avoiding this instability, depends on the stellar input physics. Using an updated version of Eggleton’s stellar code (see e.g. Tauris & Savonije 1999), Tauris et al. (2000) found $`5\mathrm{M}_{}`$ for this limit, with a corresponding maximum white dwarf mass $`0.8\mathrm{M}_{}`$ in millisecond pulsar binaries formed in this way.

## 5 The trajectory of Cygnus X–2

As described in the introduction, Cygnus X–2 is in an unusual dynamical state. In this section we investigate possible trajectories of Cygnus X–2, given the observed line–of–sight velocity and the constraints on the evolutionary state as derived in Secs. 2 and 3.

In Fig. 5 we plot the cumulative velocity distribution for systems with initial separations of $`d=15R_{}`$, i.e. immediately after circularisation of the binary orbit following the SN explosion producing the neutron star. The secondary is taken to have a mass $`M_2=3.5`$ M while the mass of the helium star prior to the supernova is $`M_1=5.0`$ M. We consider two different distributions for the kick velocity imparted to the neutron star; namely those by Hansen & Phinney (1997) and Fryer (1999). Throughout the following section, when we refer to kick velocities we mean the kick imparted on the binary allowing for the effects of mass–loss and an asymmetric supernova, and not simply the kick imparted on the neutron star from the latter.

Given its present position and velocity, we integrate the trajectory of Cygnus X–2 backwards to the birth of the neutron star using a model for the Galactic potential suggested by Paczyński (1990) (reviewed in the Appendix). From sections 2 and 3 it is clear that the age $`t`$ of the neutron star is essentially the main–sequence lifetime of the progenitor star of the present donor. To account for the allowed range of the initial donor mass and for uncertainties from the width of the main–sequence we choose $`150<t/\mathrm{Myr}<275`$. As only the line-of-sight velocity, $`𝐯_{los}`$, is known, we considered a set of trajectories where the current velocity was given by

$$𝐯=𝐯_{}+𝐯_{los}+\alpha 𝐯_1+\beta 𝐯_2$$
(2)
where $`𝐯_{los}`$, $`𝐯_1`$, and $`𝐯_2`$ are mutually orthogonal, and $`𝐯_1`$ has a zero component in the $`z`$ direction. When computing trajectories, it is most convenient to work in units related to the mass and size of the galaxy. As a natural unit of velocity we use 207 km/s, the Kepler velocity at 1 kpc distance from a point mass $`10^{10}`$ M. In the above equation $`𝐯_1`$, and $`𝐯_2`$ are given unit lengths in our code units.

We integrated the equations (A9) for a range of values of $`\alpha `$ and $`\beta `$ and noted the radius, $`r`$, and time, $`t`$, whenever the trajectory cut the Galactic plane. We also computed the kick velocity, $`V_\mathrm{k}`$ the system must have had if the neutron star has formed at that time, assuming the system had previously been in a circular orbit in the Galactic plane.

In Fig. 6 we plot the values of $`r`$ and $`V_\mathrm{k}`$ for all values of $`\alpha `$ and $`\beta `$. We note that for $`r<10`$ kpc, 140 km/s $`<V_\mathrm{k}<`$ 600 km/s. This range of required values of $`V_\mathrm{k}`$ should be compared with that expected given the neutron star kick distributions of Hansen & Phinney (1997) and Fryer (1999) in Fig. 5. It is clear from Fig. 5 that the expected value of $`V_\mathrm{k}`$ is far lower than that required for many of the trajectories shown in Fig. 6. In other words, in a large fraction of systems formed at a radius $`r<10`$ kpc, the trajectory of the binary will confine the system to the inner regions of the Galaxy; Cygnus X–2 must have been relatively unusual in reaching out beyond the solar circle.

We plot the time of formation, $`t`$, as a function of $`r`$ for all values of $`\alpha `$ and $`\beta `$ in Fig. 7. As in Fig. 6, the data points where $`150`$ Myr $`<t<275`$ Myr and $`V_k200`$ km/s are plotted as larger dots. From the work of Sections 2 and 3 these are the values most likely to be applicable to a progenitor of Cygnus X–2. We note from Fig. 7 that the binary is unlikely to have originated within 5 kpc of the Galactic centre as trajectories cutting the Galactic plane at such small radii do so at the wrong time. Fig. 6 also shows they would require unreasonably large kicks.

In Fig. 8 we plot the values of $`\alpha `$ and $`\beta `$ for those trajectories which cut the Galactic plane $`t`$ Myr ago, where $`150<t<275`$, with system kick velocities required to take the binary from a circular orbit restricted by $`V_k200`$ km/s. Values of $`\alpha `$ and $`\beta `$ besides those shown in Fig. 8 were considered, but none satisfied the above conditions. We see that only a small number of possible trajectories are credible with these restrictions applied. Here we note the effect of the Galactic rotation. For $`\alpha <0.4`$ the trajectories are in the opposite direction to the Galactic rotation. In other words the kick the system receives on formation will oppose its initial velocity in a circular orbit, and will therefore have to be larger than for a similar trajectory in the same direction as the Galactic rotation. Here the effect of the kick is boosted by the initial velocity of the circular orbit.

For illustration, in Fig. 9 we plot the trajectory ($`R`$ and $`z`$) for $`\alpha =0.75`$ and $`\beta =0.2`$. In this case, the trajectory cuts the Galactic plane on two occasions within the last 250 Myr, the first occasion being 95 Myr ago, at a radius, $`R_{\mathrm{cut}}=14.1`$ kpc, and the second occasion being 172 Myr ago, at a radius , $`R_{\mathrm{cut}}=5.9`$ kpc. If the system was formed in the disc at this latter time, and was initially travelling on a circular orbit, it would have had to receive a kick, $`V_\mathrm{k}130`$ km/s, which sits comfortably within the range given in Fig. 5.

Thus far we have demonstrated that there are viable trajectories for Cygnus X–2 compatible with the observed line-of-sight velocity. We now ask the reverse question: what fraction of binaries originating at some radius will produce systems resembling Cygnus X–2? Given that we expect systems to have originated within the solar circle, as most massive stars are located within 10 kpc of the Galactic centre, we might suspect that Cygnus X–2 is relatively unusual in being at a larger distance from the Galactic centre and significantly away from the Galactic plane. By Monte Carlo simulation we were able to produce a large number of systems at various initial radii, assumed to be initially on circular orbits, with kick velocities drawn from the distributions given in Fig. 5, and follow their trajectories in the Galactic potential. Investigation demonstrated that in determining whether a particular trajectory would produce a Cygnus X–2 like system, the maximum radius reached by the system, $`R_{\mathrm{max}}`$ was a good diagnostic. This is illustrated in Fig. 10 where we plot the radii of 100 binaries as a function of time where the initial trajectories are drawn randomly, from the kick distribution of Hansen and Phinney (see Fig. 5), and choosing the initial radius to be between 2 kpc and 10 kpc. The trajectories fall into two categories: those where the binary remains at a radius similar to that at which it was located initially, and those where the binary is ejected significantly into the Galactic halo, reaching maximum radii of $`>10`$ kpc. Cygnus X–2 clearly belongs to the latter category. In order for a system to be at a radius today similar to that of Cygnus X–2, we require $`15<R_{\mathrm{max}}<20`$ kpc, providing the system originated in the Galactic disc somewhere within 10 kpc of the Galactic centre. A number of systems will be ejected to even larger radii. In such cases a binary of age similar to Cygnus X–2 would still be on the outward bound portion of its orbit.

In Fig. 11 we plot $`R_{\mathrm{max}}`$ as a function of initial radius $`R`$. In this case the initial radius was chosen randomly from 2 kpc to 10 kpc, and the kick velocity drawn from the distribution given by Hansen and Phinney (see Fig. 5). From this figure we note that Cygnus X–2–like objects may be produced when $`r>5`$ kpc, and that the relative frequency is relatively independent of formation radius. Systems having larger values of $`R_{\mathrm{max}}`$ will also be produced, the frequency increasing with radius of formation. The relative frequency for Cygnus X–2–like systems and those on longer period orbits is plotted as a function of $`R`$ in Fig. 12.

Assuming that the formation rate of binaries scales as the surface density of stars, which is given by $`\mathrm{\Sigma }\mathrm{\Sigma }_0e^{r/3\mathrm{k}\mathrm{p}\mathrm{c}}`$ (see Binney & Merrifield 1999), we computed the relative number of Cygnus X–2–like systems which would be produced by integrating over the entire disc within the solar circle. The fraction of systems resembling Cygnus X–2 ($`f_{\mathrm{cyg}}`$) and the total fraction of systems which will travel significantly further out than their formation radius ($`f_{\mathrm{out}}`$) are listed in Table 2 as a function of the initial separation within the post-supernova binary, $`d`$, the primary helium star mass prior to the supernova, $`M_1`$, and the mass of the secondary, $`M_2`$. We found that $`7`$% of all binaries will produce Cygnus X–2–like binaries on trajectories that would place them at a radius similar to Cygnus X–2 today. A further $`515`$% will be further out than Cygnus X–2, whilst the remainder will be located closer to the Galactic centre. These results apply equally to both velocity distributions plotted in Fig. 5.

## 6 Summary

Our binary evolution calculations with full stellar models have verified that Cygnus X–2 can be understood as the descendant of an intermediate–mass X–ray binary (IMXB). The most likely evolutionary solution is an early massive case B sequence, starting from a donor with mass $`3.5\mathrm{M}_{}`$. This implies that the presently observed state is rather short–lived, of order 3 Myr. This in turn points to a large IMXB formation rate of order $`10^710^6`$ yr<sup>-1</sup> if Cygnus X–2 is the only such system in the Galaxy. It is likely that there are more as yet unrecognized IMXBs, so the Galactic IMXB formation rate could be even higher. Hercules X–1 seems to be in the early stages of a similar case A or case AB evolution (e.g. van den Heuvel 1981). The prehistory of Cygnus X–2 is sensitive to the width of the main–sequence band in the HR diagram, i.e. to convective overshooting in that phase. The alternative evolutionary solution for Cygnus X–2 suggested by PR, a case A mass transfer followed by a case B phase, is viable only if overshooting is very effective, i.e. more effective than hitherto assumed in the literature.

Using Mazzitelli’s stellar code (Mazzitelli 1989) we found that neutron star IMXBs that start case B mass transfer with initial donor mass $`4\mathrm{M}_{}`$ will encounter a delayed dynamical instability. (Note that the maximum donor mass that just avoids this instability depends on stellar input physics). The components are likely to merge and perhaps form a low–mass black hole. Given the high formation rate of Cygnus X–2–like objects the formation rate of such black holes could be rather substantial.

We have shown that the large Galactocentric distance of Cygnus X–2 and its high negative radial velocity do not require unusual circumstances at birth of the neutron star. There are viable trajectories for Cygnus X–2 that firstly cut the Galactic plane inside the solar circle at a time consistent with the evolutionary age of Cygnus X–2, and secondly which require only a moderate neutron star kick velocity ($`<200`$km/s). We estimate the fraction of systems resembling Cygnus X–2, i.e. of systems on orbits that reach Galactocentric distances $`>15`$ kpc, as $`1020\%`$ of the entire IMXB population of systems.

### Acknowledgements

MBD gratefully acknowledges the support of a URF from the Royal Society. ARK thanks the UK Particle Physics & Astronomy Research Council for a Senior Fellowship. This work was partially supported by a PPARC short–term visitors grant. We thank the anonymous referee for a careful reading of the manuscript and for comments that helped to improve the paper.

APPENDIX: THE GALACTIC POTENTIAL

The Galactic potential can be modelled as the sum of three potentials. The spheroid and disc components are given by

$$\mathrm{\Phi }_s(R,z)=\frac{GM_s}{\left(R^2+[a_s+(z^2+b_s^2)^{1/2}]^2\right)^{1/2}}$$
(3)
$$\mathrm{\Phi }_d(R,z)=\frac{GM_d}{\left(R^2+[a_d+(z^2+b_d^2)^{1/2}]^2\right)^{1/2}}$$
(4)
where $`R^2=x^2+y^2`$. The component from the Galactic halo can be derived assuming a halo density distribution, $`\rho _\mathrm{h}`$, given by

$$\rho _h=\frac{\rho _c}{1+(r/r_c)^2}$$
(5)
where $`r^2=x^2+y^2+z^2`$. The above density distribution yields the potential

$$\mathrm{\Phi }_h=\frac{GM_c}{r_c}\left[\frac{1}{2}\mathrm{ln}\left(1+\frac{r^2}{r_c^2}\right)+\frac{r_c}{r}\mathrm{atan}\left(\frac{r}{r_c}\right)\right]$$
(6)
where $`M_c=4\pi \rho _cr_c^3`$. The total Galactic potential is the sum

$$\mathrm{\Phi }=\mathrm{\Phi }_s+\mathrm{\Phi }_d+\mathrm{\Phi }_h$$
(7)
Following Paczyński (1990), we use the following choice of parameters:

$$a_s=0,b_s=0.277\mathrm{kpc},M_s=1.12\times 10^{10}M_{},$$
(8)
$$a_d=3.7\mathrm{kpc},b_d=0.20\mathrm{kpc},M_d=8.07\times 10^{10}M_{},$$
(9)
$$r_c=6.0\mathrm{kpc},M_c=5.0\times 10^{10}M_{},$$
(10)
Because of the cylindrical symmetry of the potential, the integration of the trajectories can be simplified to consider the evolution of the $`z`$ and $`R`$ components only, as given below

$`{\displaystyle \frac{dR}{dt}}`$ $`=`$ $`v_R,{\displaystyle \frac{dz}{dt}}=v_z,`$
$`{\displaystyle \frac{dv_R}{dt}}`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Phi }}{R}}\right)_z+{\displaystyle \frac{j_z^2}{R^3}},{\displaystyle \frac{dv_z}{dt}}=\left({\displaystyle \frac{\mathrm{\Phi }}{z}}\right)_R`$ (11)
where the $`z`$ component of the angular momentum, $`j_z=Rv_\varphi `$.


format: text

text:
# Moments of the ARPES spectral function of an undoped Mott insulator.

## I The Calculation

### A Half-filled Hubbard model at large U

The large U Hubbard model at half filling is the simplest model that produces the large charge gap and the antiferromagnetic tendencies observed experimentally in CuO planes. The Hamiltonian,

$$\widehat{H}=\underset{\sigma ,i,j}{}t_{ij}C_{\stackrel{}{i},\sigma }^{}C_{\stackrel{}{j},\sigma }+U\underset{\stackrel{}{j}}{}n_\stackrel{}{j}^\sigma n_\stackrel{}{j}^\sigma $$
(1)
contains a strong onsite repulsion term $`U_\stackrel{}{j}n_\stackrel{}{j}^\sigma n_\stackrel{}{j}^\sigma `$ that freezes (together with the lattice potential) charge motion at half-filling provided that the kinetic energy term is small enough to be treated as a perturbation, $`|t_{ij}|U`$.

When $`t_{ij}/U=0`$ the ground state consists of singly occupied sites, has energy zero, is $`2^N`$-fold degenerate where N is the size of the system, and is completely characterized by the spin configuration, $`\{S_i\}`$.

The perturbative effects of the kinetic energy modify all of the above statements except the last one: each of the perturbed states, though it contains an admixture of doubly occupied sites, can still be labeled by the spin configuration of the unperturbed state from which it has evolved. Thus the expectation value of any operator in the ground state manifold can always be expressed in terms of spin variables. Formally this is accomplished by computing perturbatively in powers of $`t_{ij}/U`$ the unitary transformation $`\mathrm{exp}[i\widehat{X}]`$ that expresses the evolution of the low energy unpeturbed states $`|\{S_i\}>`$ as a function of increasing kinetic energy:

$`\stackrel{~}{|\{S_i\}>}`$ $`=`$ $`e^{i\widehat{X}}|\{S_i\}>`$ (2)
$`<\widehat{𝒪}>`$ $`=`$ $`<\{S_i\}|e^{i\widehat{X}}\widehat{𝒪}e^{i\widehat{X}}|\{S_i\}>,`$ (3)
where $`𝒪`$ is any observable.

In the familiar fashion, this transformation maps the low energy physics of the Hubbard Hamiltonian into an effective Heisenberg antiferromagnet whose leading term is:

$`\widehat{H}_{eff}`$ $`=`$ $`e^{i\widehat{X}}\widehat{H}e^{i\widehat{X}}`$ (4)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}{\displaystyle \frac{4t_{ij}^2}{U}}({\displaystyle \frac{1}{4}}\stackrel{}{S_i}\stackrel{}{S_j})+𝒪(t_{ij}^4/U^3).`$ (5)
### B The spectral function

The emission spectral function is

$`A(\stackrel{}{k},w)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma ,m,n}{}}e^{\beta E_n}|<m|C_{\stackrel{}{k},\sigma }^{}|n>|^2\delta (\omega +E_mE_n)`$ (6)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma }{}}{\displaystyle 𝑑te^{iwt}<C_{\stackrel{}{k},\sigma }^{}(t)C_{\stackrel{}{k},\sigma }(0)>}.`$ (7)
Frequency moments of $`A(\stackrel{}{k},w)`$ correspond to ground state (or thermodynamic) averages of the following operators:

$`n_\stackrel{}{k}`$ $``$ $`{\displaystyle \frac{dw}{2\pi }A(k,w)}={\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma }{}}<C_{\stackrel{}{k},\sigma }^{}C_{\stackrel{}{k},\sigma }>`$ (8)
$`A_1(\stackrel{}{k})`$ $``$ $`{\displaystyle \frac{dw}{2\pi }\omega A(\stackrel{}{k},w)}`$ (9)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma }{}}<C_{\stackrel{}{k},\sigma }^{}[\widehat{H},C_{\stackrel{}{k},\sigma }]>`$ (10)
$`A_2(\stackrel{}{k})`$ $``$ $`{\displaystyle \frac{dw}{2\pi }\omega ^2A(\stackrel{}{k},w)}`$ (11)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma }{}}<C_{\stackrel{}{k},\sigma }^{}[\widehat{H},[\widehat{H},C_{\stackrel{}{k},\sigma }]]>`$ (12)
A systematic evaluation of these averages, using the perturbative expression for $`\widehat{X}`$ in powers of $`t_{ij}/U`$, can be found in the Appendix.

## II The spectral weight, $`n_\stackrel{}{k}`$

We begin with a discussion of the spectral weight because it has an immediate physical interpretation as the occupation probability. As shown in the Appendix,

$$n_\stackrel{}{k}=\frac{1}{2}[1\frac{4\stackrel{~}{ϵ_k}}{U}+O(t^3/U^3)]$$
(13)
where

$$\stackrel{~}{ϵ_k}=\underset{j}{}S_{0i}t_{0j}e^{i\stackrel{}{k}\stackrel{}{R}_j}$$
(14)
and

$$S_{ij}=<\frac{1}{4}\stackrel{}{S_i}\stackrel{}{S_j}>$$
(15)
is the equilibrium spin correlation between spins i and j. $`\stackrel{~}{ϵ_k}`$ is a sort of renormalized band energy in which each hopping matrix element $`t_{ij}`$ is renormalized by a factor of $`S_{ij}`$. However, this renormalized energy does not correspond in any simple way to the energy of any elementary excitation of the system. Note that the non-interacting free electron band is given by $`ϵ_k=_jt_{0j}e^{i\stackrel{}{k}\stackrel{}{R}_j}`$.

This expression for $`n_\stackrel{}{k}`$ can be derived in a different, simpler manner, which is readily generalizible to more complicated situations, such as the three band Cu-O or Emery model. From the Hellman-Feynman theorem, it follows that

$$\underset{\sigma }{}<[c_{i\sigma }^{}c_{j\sigma }+\mathrm{H}.\mathrm{C}.]>=E/t_{ij}=[J_{ij}/t_{ij}]S_{ij}$$
(16)
where $`E`$ is the internal energy, which can be computed using the effective Hamiltonian $`\widehat{H}_{eff}`$. From the expression for $`J_{ij}`$ in terms of $`t_{ij}`$, the result in Eq.13 follows immediately.

### A The remnant fermi surface

At low temperatures, the short-range spin correlation functions are essentially temperature independent, and equal to their value in the ground state. To be concrete, let us consider the Hubbard model with nearest, second, and third neighbor hopping, $`t`$, $`t^{}`$ and $`t^{\prime \prime }`$, respectively; this sort of model was used in the numerical studies to fit the dispersions seen in ARPES. The zero temperature spin correlations of the corresponding spin 1/2 Heisenberg model have been computed fairly accurately in numerical studies. For the Heisenberg model with only nearest-neighbor exchange coupling, the spin correlation functions are $`S_{01}7/12`$ and $`S_{03}S_{02}\frac{1}{20}`$ between nearest, next-nearest, and third nearest neighbor sites, respectively. These correlations are, moreover, found to be relatively insenstitive to the inclusion of a modest amount of further neighbor exchange couplings, which anyway are expected to be quite small since $`J^{}/J=[t^{}/t]^2`$.

In computing $`\stackrel{~}{ϵ}_k`$ the antiferromagnetic correlations between neighboring spins imply a factor of 1/2 renormalization of $`t`$, compared to a factor of 1/20 renormalization of $`t^{}`$ and $`t^{\prime \prime }`$. Since in most cases of physical interest, $`|t^{}|,|t^{\prime \prime }|10|t|`$, even when $`t^{}`$ and $`t^{\prime \prime }`$ are large enough to make signicant shifts in the original Fermi surface defined by $`ϵ_k`$, the occupation probability is well approximated (Figure 1) as

$$n_\stackrel{}{k}\frac{1}{2}[1+\frac{7t(cos(k_x)+cos(k_y))}{3U}]$$
(17)
It is both $`\frac{1}{2}`$ and has the steepest slope along the Fermi surface of the non-interacting electrons with only the nearest neighbor hopping. In fact, ARPES experiments that observe such momentum dependence of $`n_\stackrel{}{k}`$ have been interpreted as an indication that there is a “remnant Fermi surface” in the undoped Mott insulator. By contrast, our result suggests that the observed $`n_\stackrel{}{k}`$ is reflective of the the spin physics of the strongly correlated Neel state rather than a vestige of the original Fermi surface.

### B Specific heat of the Neel transition

The natural connection we find between the spectral weight and the Neel state can be exploited further. Since, as we already pointed out, the Heisenberg Hamiltonian is dominated by the nearest neighbor term, and assuming that $`J_{n.n.}`$ is only weakly (if at all) temperature dependent, the specific heat is

$$C(T)=J\frac{}{T}<\frac{1}{4}\stackrel{}{S_1}\stackrel{}{S_0}>=J\frac{S_1}{T}.$$
(18)
Since

$$n_\stackrel{}{k}\frac{1}{2}[1+\frac{8t(cos(k_x)+cos(k_y))}{U}S_1],$$
(19)
by measuring the temperature dependence of $`S_1`$ (as extracted from $`n_\stackrel{}{k}`$) and differentiating it one gets the magnetic contribution to the specific heat. The specific heat extracted in this manner contains only the contribution from the spin fluctuations.

## III Higher Moments

Higher spectral moments correspond to the derivatives $`<C_{\stackrel{}{k},}^{}(t)C_{\stackrel{}{k},}(0)>`$ at $`t=0`$ and thus provide further insight into the problem of a hole in a Mott insulator. Our results for first and second moments, $`A_1(\stackrel{}{k})`$ and $`A_2(\stackrel{}{k})`$, are presented in a table below. The moments are written as sums over momentum space Fourier harmonics, each corresponding to a summation over sites a given near neighbor distance away on the square lattice.

$`A_n(\stackrel{}{k})`$ $`=`$ $`{\displaystyle \underset{j}{}}A_n^j\gamma _j(\stackrel{}{k}),`$ (20)
$`\gamma _j(\stackrel{}{k})`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{R}=j^{th}n.n.}{}}e^{i\stackrel{}{k}\stackrel{}{R_j}}`$ (21)
| i | $`A_1^j`$ | $`A_2^j`$ |
| --- | --- | --- |
| 1 | $`t(\frac{1}{2}S_1)`$ | $`(2tt^{}+tt^{\prime \prime })(S_11)`$ |
|  | $`\frac{2t^{}t^{\prime \prime }}{U}(S_2+2S_1)`$ |  |
|  | $`\frac{tt^{\prime \prime }}{U}(S_3+2S_1)`$ |  |
| 2 | $`t^{}(\frac{1}{2}S_2)`$ | $`(2t^{}t^{\prime \prime }+t^2)(S_21)`$ |
|  | $`\frac{2t^{}t^{\prime \prime }}{U}(S_3+2S_2)`$ |  |
|  | $`\frac{t^2}{U}(S_2+2S_1)`$ |  |
| 3 | $`t^{\prime \prime }(\frac{1}{2}S_1)`$ | $`((t^{})^2+\frac{t^2}{2})(S_31)`$ |
|  | $`\frac{(t^{})^2}{U}(S_3+2S_2)`$ |  |
|  | $`\frac{t^2}{2U}(S_3+2S_1)`$ |  |
| 4 | $`\frac{tt^{}}{U}(S_4+S_2+S_1)`$ | $`(t^{}t^{\prime \prime }+tt^{\prime \prime })(S_41)`$ |
|  | $`\frac{tt^{\prime \prime }}{U}(S_4+S_3+S_1)`$ |  |
| 5 | $`\frac{(t^{})^2}{U}(S_5+2S_2)`$ | $`(\frac{(t^{})^2}{2}+(t^{\prime \prime })^2)(S_51)`$ |
|  | $`\frac{(t^{\prime \prime })^2}{U}(S_5+2S_3)`$ |  |
| 6 | $`\frac{tt^{\prime \prime }}{U}(S_6+S_3+S_1)`$ | $`tt^{\prime \prime }(S_61)`$ |
| 7 | $`\frac{t^{}t^{\prime \prime }}{U}(S_7+S_3+S_2)`$ | $`t^{}t^{\prime \prime }(S_71)`$ |
| 10 | $`\frac{(t^{\prime \prime })^2}{2U}(S_{10}+2S_3)`$ | $`\frac{(t^{\prime \prime })^2}{2}(S_{10}1)`$ |

One notices that all computed moments are finite. Although the existance of the moment expansion is a general requirement of any physical system, it is more rule than exception that approximations lead to divergencies past some finite order. Though we haven’t constructed an explicit proof, there are indications that the expansion is well behaved for the Hubbard model.

In principle, the moment expansion can be used to study the QP dispersion ($`E_\stackrel{}{k}`$) directly: a coherent oscillation results in $`E_\stackrel{}{k}^n`$ contribution to $`A_n(\stackrel{}{k})`$. The crudest (single mode) approximation of this sort identifies $`\overline{E}_\stackrel{}{k}=A_1(\stackrel{}{k})/n_\stackrel{}{k}`$. We found it to be in a surprising agreement with previously obtained results for the momentum dependence of the QP energy in t-t’-t”-J as well as t-J models (one needs to assume that term proportional to t does not contribute to QP dispersion). Since $`\overline{E}_\stackrel{}{k}`$ rather seriously overestimates the overall bandwidth, it isn’t clear if one is justified claiming to have obtained even an approximate QP energy yet. Another likely use of our results (or rather their extension to higher moments and orders in perturbation theory) can be in comparing with spectral functions obtained by other means (either numerics, self-consistent Born approximation or other).

In conclusion, we have outlined a well controlled method for analysing the spectral moments of a hole in a Mott insulator. We find the occupation probability, $`n_\stackrel{}{k}`$, is in agreement with the well established experimental result, which as our calculation suggests is strongly constrained by the presence of AF order. We further propose that the temperature dependence of $`n_\stackrel{}{k}`$ can be used to study the specific heat of the Neel transition. The implications of our results for higher moments are yet to be understood properly.

I would like to thank Steven Kivelson for innumerable illuminating conversations and suggesting this problem in the first place. I would also like to greatfully acknowledge M.Z.Hasan, A. Lanzara, F.Ronning and Z.X.Shen for useful discussions and the hospitality of the Physics and Applied Physics Departments at Stanford University where this work was carried out. This work was supported by NSF grant number DMR98-08685.

## A Evaluation of moments

To compute the moments we first express them in terms of real space electronic correlations

$`n_k`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma }{}}<C_{\stackrel{}{k},\sigma }^{}C_{\stackrel{}{k},\sigma }>`$ (A1)
$`=`$ $`{\displaystyle \frac{1}{2}}\{1+{\displaystyle \underset{i,\sigma }{}}<C_{\stackrel{}{i},\sigma }^{}C_{\stackrel{}{0},\sigma }>\gamma _i(k)\}`$ (A2)
$`A_1(\stackrel{}{k})`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma }{}}<C_{\stackrel{}{k},\sigma }^{}[\widehat{H},C_{\stackrel{}{k},\sigma }]>`$ (A3)
$`=`$ $`ϵ_kn_k+{\displaystyle \frac{U}{2}}{\displaystyle \underset{i,\sigma }{}}\gamma _i(k)<C_{\stackrel{}{i},\sigma }^{}C_{\stackrel{}{0},\sigma }n_{0,\sigma }>`$ (A4)
$`A_2(\stackrel{}{k})`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma }{}}<C_{\stackrel{}{k},\sigma }^{}[\widehat{H},[\widehat{H},C_{\stackrel{}{k},\sigma }]]>`$ (A5)
$`=`$ $`ϵ_kn_\stackrel{}{k}\overline{E}_\stackrel{}{k}+{\displaystyle \frac{U^2}{2}}{\displaystyle \underset{i,\sigma }{}}\gamma _i(k)<C_{\stackrel{}{i},\sigma }^{}C_{\stackrel{}{0},\sigma }n_{0,\sigma }>`$ (A6)
$`+`$ $`{\displaystyle \frac{U}{2}}{\displaystyle \underset{i,j,\sigma }{}}\gamma _i(k)t_{j,0}<C_{\stackrel{}{i},\sigma }^{}(C_{\stackrel{}{j},\sigma }n_{0,\sigma }`$ (A7)
$`+`$ $`C_{\stackrel{}{j},\sigma }C_{\stackrel{}{0},\sigma }^{}C_{\stackrel{}{0},\sigma }+C_{\stackrel{}{j},\sigma }^{}C_{\stackrel{}{0},\sigma }C_{\stackrel{}{0},\sigma }>`$ (A8)
As before $`\gamma _i(\stackrel{}{k})=_{R_i}e^{i\stackrel{}{k}\stackrel{}{R_i}}`$,where the sum is over i’th nearest neighbors.

Next, these electronic correlations are evaluated perturbatively ($`S_{ij}=<\frac{1}{4}\stackrel{}{S}_i\stackrel{}{S}_j>`$):

$`{\displaystyle \underset{\sigma }{}}<C_{\stackrel{}{i},\sigma }^{}C_{\stackrel{}{0},\sigma }>`$ $`=`$ $`{\displaystyle \frac{4t_{i0}}{U}}S_{i0}+O(t^3/U^3)`$ (A9)
$`{\displaystyle \underset{\sigma }{}}<C_{\stackrel{}{i},\sigma }^{}C_{\stackrel{}{0},\sigma }n_{0,\sigma }>`$ $`=`$ $`{\displaystyle \frac{2t_{i0}}{U}}S_{i0}`$ (A10)
$`+{\displaystyle \underset{j,R_j}{}}{\displaystyle \frac{t_{iR_j}t_{R_j0}}{U^2}}\{3S_{iR_j}S_{R_j0}`$ $``$ $`S_{i0}+3i<\stackrel{}{S}_i(\stackrel{}{S}_{R_j}\times \stackrel{}{S}_0)>\}`$ (A11)
On a square lattice and in a state that doesn’t break time reflection invariance ($`<\stackrel{}{S}_i(\stackrel{}{S}_{R_j}\times \stackrel{}{S}_0)>=0`$) the second sum is (for different i):

| 1 | $`4tt^{}S_{20}+2tt^{\prime \prime }S_{30}`$ |
| --- | --- |
| 2 | $`2(2S_1S_2)t^2+4S_3t^{}t^{\prime \prime }`$ |
| 3 | $`(2S_1S_3)t^2+2(2S_2S_3)t^2`$ |
| 4 | 2$`(S_1+S_2S_4)tt^{}+2(S_1+S_3S_4)tt^{\prime \prime }`$ |
| 5 | $`(2S_2S_5)t^2+2(2S_3S_5)t^{\prime \prime 2}`$ |
| 6 | 2$`(S_1+S_3S_6)tt^{\prime \prime }`$ |
| 7 | 2$`(S_2+S_3S_7)t^{}t^{\prime \prime }`$ |
| 10 | $`(2S_3S_{10})t^{\prime \prime 2}`$ |

Substituting these terms into A1,A2,A3 yields results of Section III.


text:
# 1 Introduction

## 1 Introduction

Gravitational collapse of a stellar core to a black hole has been studied since many years. Efforts have been also devoted to the calculation of the gravitational radiation emitted in this process. Most studies have been based on a perturbative approach,?,?,?,?,? others on the numerical solution of the full Einstein equations ?,?. All these papers consider the collapse of a “naked” stellar core, described by a “dust” of particles or, at most, by a politropic equation of state, without taking into account the presence of the outer layers of the star, which are involved in the process. In particular, depending on the ratio between the energy released in the final explosion of a massive star and the binding energy of the ejected material, two different kinds of collapse have been outlined ?,?: the prompt collapse, in which a large fraction of the star collapses on a dynamical time-scale forming a massive black hole, and the delayed collapse, in which a low mass black hole or, alternatively, a neutron star forms at the beginning and later accretes matter, due to fall-back. If a neutron star is the initial outcome of the collapse, fall-back pushes its mass above the critical mass and the formation of a black hole takes place. In both cases, this light black hole continues to slowly accrete matter until the final mass is reached.

In the gravitational wave community the prompt collapse has been considered for a long time as representative of realistic collapse processes. The estimated mass (several solar masses) of the first black hole candidates (like $`CygnusX1`$) has led to the idea that black holes should often be born with a ”typical” mass of $`10M_{}`$. This assumption appears no more justified now: we know that the evolution of massive stars in binary systems (to which all observed black hole candidates belong) is different from that of single stars; in addition, more refined observational techniques have allowed to find black hole candidates of few solar masses. In the light of these results and of the recent numerical simulations which we will describe, the prompt birth of such massive black holes ($`10M_{}`$ or more) should be considered as a very rare event.

In this paper, we will discuss the delayed collapse model, the expected mass distribution of isolated black holes and the detectability of the emitted signal, by forthcoming interferometric detectors. Our main aim is to understand how our perspectives of detection change with respect to the ”naive” prompt collapse. The plan of the paper is as follows. In Sec.2 we will shortly describe the delayed collapse scenario. In Sec.3 we will estimate the contribution of fall-back to the total gravitational emission in the delayed collapse discussing the consequences. In Sec.4 we will discuss the detectability of the emitted signal by forthcoming interferometric detectors, and compare with the prompt collapse which could happen for progenitor stars above about $`40M_{}`$, if stellar winds were much less important in the evolution of massive stars than it is currently believed. In Sec.5 we will derive the theoretical final mass distribution function for isolated black holes, using the results of recent simulations of the collapse of massive single stars. Finally, in Sec.6 the results and their implications will be discussed.

## 2 Star Collapse to a Black Hole

The fate of a massive star ($`m_{prog}>9M_{}`$) is the core collapse with the formation of a compact object: a neutron star or a black hole. Numerical simulations show that the actual final product, and also the way in which it is formed, depend on many factors, among which the mass and the angular momentum of the progenitor, the explosion energy, the high density matter equation of state ($`EOS`$) and also the way in which the physics is implemented (regarding, for instance, neutrino physics or angular momentum transport). Moreover, results cannot be considered conclusive until fully relativistic 3D simulations will be performed. It has been shown that, for a wide range of progenitor masses and explosion energies, the shock cannot expell all the matter outside the collapsing core, so that part of the helium mantle and heavy elements may slow down below the escape velocity and be accreted by the just formed neutron star, with a timescale of minutes to hours. If the neutron star mass grows above a critical value, a black hole forms and we have a delayed collapse. On the other hand, if the explosion completely fails, or if it is too weak, a black hole immediately forms: this is the prompt collapse. According to recent simulations (starting from non-rotating progenitors) by Woosley $`\&`$ Weaver ? and by Fryer $`\&`$ Kalogera ?, typical collapses are always delayed, the prompt ones occurring only for very massive progenitor stars ($`M>40M_{}`$) if stellar winds are negligible, an assumption which appears rather unlikely. The $`EOS`$ of high density matter plays a basic role in the determination of neutron stars limiting mass, i.e. the critical mass for black hole formation, $`m_{min}`$. For conventional EOS the neutron star maximum mass ranges between $`1.7M_{}`$ and $`2.2M_{}`$ (e.g. ?, ?,?,?), with a $`10÷20\%`$ increase if rotation is taken into account. On the other hand, if $`\pi `$ or $`K`$ condensation, or formation of quark matter, takes place at very high densities (e.g. ?,?,?,?), the $`EOS`$ is softened and this reduces the maximum mass that the pressure of degenerate matter can sustain. In this case, the critical neutron star mass is $`1.5M_{}`$. The masses of 26 neutron stars, measured in pulsars, are compatible with a gaussian distribution with $`\overline{m}=1.35\pm 0.04`$, and then are consistent with that previous limit ?. Moreover, the lack of evidence for neutron stars with mass near the maximum derived from conventional $`EOS`$ is reinforced by the lack of evidence for a pulsar as a remnant of the supernova $`SN1987A`$ whose progenitor had a mass of $`18M_{}`$, which should have left behind a remnant mass of about $`1.5M_{}`$ ?,?. However, no definite conclusion can be still drawn?,?. Given the uncertainties in the neutron stars equation of state, in the following we will always refer to two different values of the black hole minimum mass: $`m_{min}=1.5M_{}`$, representative of soft equations of state (soft EOS), and $`m_{min}=2M_{}`$ for standard, with no phase transition, equations of state (conventional EOS).

## 3 Gravitational Radiation from Fall-Back

In this section, we want to estimate the contribution of fall-back to the emission of gravitational waves in the delayed collapse of a massive star to a black hole.

There are many studies on the capture by a black hole of finite-size shells of matter ?,?,?. A general result is that the amount of gravitational radiation emitted is always smaller than that emitted by the capture of a pointlike particle with the same mass of the shell. This is a consequence of destructive interference of the radiation emitted by different parts of the infalling extended matter. These results are confirmed also by the recent, more realistic, calculations by Papadopoulos $`\&`$ Font ?. In the case of an axisymmetric irrotational shell of matter with mass $`\mu `$, much smaller than the mass $`m`$ of the black hole, they find that the efficiency in the emission of gravitational waves decreases as a function of the radial width $`L`$ of the shell as follows:

$$ϵ_s=\frac{\mathrm{\Delta }E}{mc^2}=810^3\left(\frac{\mu }{m}\right)^2\left(\frac{m}{M_{}}\right)^{2.4}\left(\frac{L}{1km}\right)^{2.4}$$
(1)
We shall now extrapolate this equation for $`\mu >m`$ and compare to the efficiency we expect from core collapse which, for a maximally rotating core, is $`ϵ_c710^4`$ in the axisymmeric case ?. For instance, assuming that the initial mass of the newly formed black hole is $`m=2M_{}`$, and that the mass of the shell is $`\mu =10M_{}`$, we find that the condition $`ϵ_s<0.1ϵ_c`$ requires $`L>100km`$. This condition is largely verified for massive stars. Then, we can conclude that the main burst of gravitational radiation is emitted at the formation of the black hole while subsequent accretion of matter gives no important contributions. To this respect, two possibilities can occur. First, the $`Fe`$ core of the star has a mass greater than $`m_{min}`$. In such a case the black hole initial mass is equal to the core mass. Second, the $`Fe`$ core is lighter than $`m_{min}`$. In this case, a neutron star is produced at the beginning. It then accretes matter until the critical mass is reached, so that the black hole initial mass is equal to the minimum one. From the core masses given by Woosley $`\&`$ Weaver ? we see that for soft EOS ($`m_{min}=1.5M_{}`$) all stars with initial mass greater than $`18M_{}`$ produce cores that immediately collapse to a black hole of mass greater than $`m_{min}`$, but lower than $`1.8M_{}`$. The same happens according to Fryer ?. For conventional EOS ($`m_{min}=2M_{}`$) black holes are formed, after fall-back on a neutron star, from the collapse of progenitors of mass $`m_{prog}>26M_{}`$ (following Fryer $`\&`$ Kalogera ? this value is reduced to $`m_{prog}20M_{}`$). Then, all black holes form with an initial mass equal, or nearly equal, to the minimum one, presumably in the range $`1.5÷2.6M_{}`$$`^\text{e}`$footnotetext: $`^\text{e}`$The upper limit takes into account also the possible stabilizing effect of rotation., depending on the nuclear density matter equation of state and on the rotation rate. In the following we will assume, for simplicity of calculation, that all black holes form with the same mass $`m_{min}`$, so that their initial mass distribution function, for the delayed collapse scenario, can be written as a $`\delta `$-function: $`f(m)=\delta (mm_{min})`$. We will consider $`m_{min}=1.5M_{}`$ and $`m_{min}=2M_{}`$.

### 3.1 The rate of black hole formation

Let us now estimate the expected event rate for collapses. We assume a Galactic rate for Supernovae of type $`II`$ given by $`R_{SNII}=0.02yr^1`$. Such rate is the sum of the rate of collapses to neutron stars and, in the delayed scenario, of the rate of collapses to black hole (a supernova explosion is produced by all delayed collapses, indipendently of $`m_{min}`$). For instance, for soft EOS, we have seen that a black hole is produced by progenitor stars with mass $`18M_{}m_{prog}40M_{}`$. Their formation rate, $`R_{bh}`$, is a fraction $`\lambda `$ of the rate of formation of neutron stars $`R_{ns}`$, that are generated by progenitors with mass $`9M_{}m_{prog}18M_{}`$. Of course, this ratio depends on the initial mass function of the progenitor stars: $`\lambda =0.43`$ for a Salpeter law with exponent $`\alpha =2.35`$, while $`\lambda =0.32`$ for a Scalo law with $`\alpha =2.7`$. Then, we can write

$$R_{SNII}=R_{ns}+R_{bh}=R_{ns}(1+\lambda ),$$
(2)
from which we can get $`R_{ns}`$ and $`R_{bh}`$. In Tab.(1) we tabulate the galactic rates of black hole formation in the different cases, considering both delayed ($`R_{bh}`$) and prompt ($`R_{bh,p}`$) collapses. Moreover, the total expected rate within the Virgo cluster, $`R_{\mathrm{@}20Mpc}`$, is given in the last column.

From Tab.(1) we see that the fraction of delayed collapses producing a black hole is in the range $`10\%÷43\%`$ of those leading to a neutron star. These percentages increase to $`18\%÷60\%`$ if prompt collapse, for stars more massive than $`40M_{}`$, is also taken into account.

## 4 The Detectability

Gravitational signals emitted in star collapse to a black hole are characterized by an initial part, emitted during the in-fall phase and bounce, followed by an oscillating tail, which can be described as a superposition of damped sinusoids, corresponding to the black hole quasi-normal modes. The frequency and the damping time of the quasi-normal modes are a function of the black hole parameters, mass and angular momentum. For axisymmetric collapses, the main contribution (typically $`90\%`$ of the whole energy emitted) is given by the $`l=2`$ mode, for which the frequencies and the damping times are given, as a function of the rotation parameter $`a=J/\left(\frac{Gm^2}{c}\right)`$, in Tab.(2).

We use the waveforms calculated by Stark $`\&`$ Piran ?, who have computed, in the framework of the full non-linear theory, the gravitational signals emitted in the axisymmetric collapse to a black hole of a rotating core, for different values of the angular momentum. The maximum efficiency reached in their simulations is $`ϵ_{max}710^4`$. From the waveforms we have computed the corresponding one-sided power spectrum $`f(\nu )`$, i.e. the flux of energy per unit frequency which is given by

$$f(\nu )=\frac{\pi c^3\nu ^2}{2G}<h^2(\nu )>$$
(3)
where $`<h^2>`$ denotes the average of the squared gravitational signal with respect to its angular dependence. These computed energy spectra allow to evaluate the signal to noise ratio ($`SNR`$), according to the well-known formula given by Eq.(4) below, which holds if the matched filter is applied to the data. This optimum filtering procedures implies that we know the exact waveform emitted and this is a rather optimistic hypothesis even in the case of star collapses to a black hole. We know that in the simplest cases (low rotation, no hydrodynamical effects) most of the energy is emitted in the phase of quasi-normal ringing, which can be described as a superposition of damped sinusoids. On the other hand, if the collapse has a high degree of rotation, or if strong hydrodynamical effects take place, the emitted waveform may be no longer dominated by the quasi-normal ringing and would be less predictable ?. In this case, other, less optimum, data analysis procedures should be used ?. Then, our results have to be considered as upper limits on the $`SNR`$. The average squared $`SNR`$ can be expressed as

$$\overline{SNR^2}=\frac{8G}{5\pi c^3}_0^{\mathrm{}}\frac{f(\nu )}{\nu ^2S_h(\nu )}𝑑\nu .$$
(4)
where $`S_h(\nu )`$ is the detector noise power spectrum. In Eq.(4) an average on the source-detector relative position and on the polarization of the gravitational waves emitted is also performed. As discussed in Sec.3, the delayed collapse produces black holes that, immediately after birth, have mass very close to the minimum mass. Thus, the emitted signal depends essentially on the black hole angular momentum and on the distance at which the collapse takes place. In prompt collapses, on the contrary, the $`SNR`$ clearly depends on the mass of the black hole and, due to the typical sensitivity curve of ground-based interferometers, the radiation emitted by more massive black holes is more easily detectable, so that the delayed collapse process is definitely less favourable than the “naive” prompt one. In Tab.(3) we give the $`SNR=\sqrt{\overline{SNR^2}}`$ for delayed collapses taking place in the Galaxy (we fix $`r=10kpc`$), for different values of the minimum black hole mass and of the rotation parameter.

The values of the $`SNR`$ increase with the angular momentum of the black hole. This is the consequence of two effects: first, the efficiency of emission increases as $`a^4`$; second, for very high angular momentum (say, $`a>0.8`$) the bounce of the collapsing star produces a lower frequency component in the signal energy spectrum which fits better to the sensitivity curve of interferometers. On the other hand, the expected rate of detectable events is low: collapses to a black hole are detectable, essentially, only within the Local Group, with a total rate, strongly dominated by the Milky Way, which is less than $`1`$ event per century, see Tab.(1). We stress, however, that higher $`SNR`$ could be obtained if collapses had a higher degree of asymmetry and there are some observative indications supporting this hypothesis ?,?,?,?. Moreover, the initial core collapse to a neutron star which takes place, in the delayed case, if $`m_{min}=2M_{}`$, could be a promising process. If it is highly asymmetric, as observations seem to indicate, a large amount of gravitational radiation could be emitted in the range of frequencies where interferometric detectors reach their best sensitivity$`^\text{f}`$footnotetext: $`^\text{f}`$It should be noted, however, that according to the newtonian simulations by Rampp et al. ?, the amount of gravitational radiation emitted in non-axisymmetric collapses is comparable to that of the axisymmetric case..

We have repeated the calculation of the $`SNR`$ for the advanced $`LIGO`$ detector, using an approximation to its sensitivity curve given by Flanagan $`\&`$ Huges ?.

Results are given, for delayed collapses, in Tab.(4) where a distance of $`r=20Mpc`$ has been assumed. Within this distance, corresponding approximately to the Virgo Cluster, the expected black hole formation rate is $`1÷5yr^1`$, but we see that the detection perspectives are not much better because the $`SNR`$ is much smaller than one. The situation is different for prompt collapses. In Fig.(4) we have plotted the signal-to-noise ratio for prompt collapses as a function of the newborn black hole mass and it appears that the emitted gravitational radiation could be detectable if the rotation rate is high enough.

## 5 Black hole final mass distribution function

In this section we calculate the black hole final mass distribution function. We do not consider the chemical evolution of the galaxy, which affects the mass distribution of compact remnants (see, e.g., Timmes et al. ?). Such a theoretical distribution could be useful in statistical studies of recently formed black holes.

The black hole final mass distribution depends on the mass distribution of progenitor stars and on the relation between the mass of progenitors and that of the final remnants produced after the collapse. For the first one we assume a power law $`f(m_{prog})m_{prog}^\alpha `$, with $`\alpha =2.35`$ (Salpeter’s law) and $`\alpha =2.7`$ (Scalo’s law). Regarding the relation between the mass of progenitor stars, $`m_{prog}`$, and that of remnants, $`m`$, we use (after conversion from baryonic to gravitational mass) the results of Woosley $`\&`$ Weaver ? and of Fryer $`\&`$ Kalogera ?, who have performed systematic calculations of the evolution and explosion of non rotating massive stars. It must be stressed again that such evolutionary calculations are still subject to many uncertainties and this reflects in some differences between their results, obtained using different assumptions. The remnant masses they report are calculated a long time after the collapse, so that the possible fall-back is included. The distribution function we are searching for is simply given by

$$f(m)=f(m_{prog}(m))\frac{dm_{prog}}{dm}$$
(5)
In Fig.(5) we plot the function $`f(m)`$ obtained using a fit of the progenitor mass-remnant mass relation found by Woosley $`\&`$ Weaver ? (see their table 3), considering different values of $`m_{min}`$ and of the exponent $`\alpha `$. The maximum progenitor mass considered by Woosley $`\&`$ Weaver is $`m_{prog}=40M_{}`$. According to them, for masses greater than this, the effect of stellar wind becomes relevant and strongly modifies the star evolutionary path: a smooth convergence of remnant masses to $`1.5M_{}`$, corresponding to neutron stars or low mass black holes, is expected for progenitors of mass greater than about $`40M_{}`$. As a consequence, black holes with mass greater than $`10.3M_{}`$ \- which is what is predicted by Woosley $`\&`$ Weaver for a $`40M_{}`$ progenitor star if stellar winds are neglected - should never be formed, for progenitors of solar metallicity.$`^\text{a}`$footnotetext: $`^\text{a}`$Black holes of mass greater than that value could be produced if the metallicity is lower. Also Fryer $`\&`$ Kalogera ? show that low mass remnants are produced, for progenitors more massive than about $`40M_{}`$, if stellar winds are taken into account. On the other hand, if stellar winds are neglected, they find that, for masses above about $`40M_{}`$, the star collapse takes place in the prompt way, i.e. a black hole of mass nearly equal to the mass of the progenitor star is produced in a dynamical timescale. In this case the black hole mass distribution function is simply proportional to $`m^\alpha `$. The final mass distribution function, resulting from the progenitor mass-remnant mass relation found by Fryer $`\&`$ Kalogera, is plotted in Fig.(5) (we refer to their “most likely” model, in which $`50\%`$ of the explosion energy goes in unbinding the star).

There are clear differences between the distributions functions plotted in Figs.(5,5). First of all, the expected range of masses for black holes is larger in the second case, reaching $`14.7M_{}`$. Moreover, in this case the function $`f(m)`$ is monotonically decreasing, while in Fig.(5) there is a minimum around $`3M_{}`$ followed by a relative maximum at $`4.5M_{}`$. These differences are a consequence, obviously, of the different relations for progenitor vs. remnant masses$`^\text{b}`$footnotetext: $`^\text{b}`$In particular, the minimum in Fig.(5) is due to an increase of the slope of $`m_{prog}(m)`$ for $`m_{prog}30M_{}`$ while the following maximum is produced by the first term in the right-hand side of Eq.(5), which is dominant for large $`m_{prog}`$. and are due, mainly, to two reasons. First, explosion simulations are carried on by Woosley $`\&`$ Weaver in one-dimension, while Fryer $`\&`$ Kalogera simulations are two-dimensional. Second, Woosley $`\&`$ Weaver constrain the kinetic energy of ejected material at infinity$`^\text{c}`$footnotetext: $`^\text{c}`$Defined as the difference between the explosion energy and the binding energy of the ejecta to a constant value, about $`1.210^{51}erg`$, while in Fryer $`\&`$ Kalogera this quantity decreases for increasing progenitor mass. This explains why their remnant masses are greater.

The theoretical distribution functions we have found cannot be compared with observations, which refer to black candidates belonging to binary systems. It is widely accepted that the evolution of high mass stars in binary systems is quite different from that of single stars of comparable mass. A detailed discussion on massive binary systems evolution and on the possible mechanisms and bias effects which can explain the observed black hole candidates mass distribution can be found in refs. ? and ?.

## 6 Conclusions

In this paper we have considered the collapse of a massive star to a black hole, exploring the models of prompt and delayed collapse. According to recent simulations, the formation of stellar mass black holes should take place through a delayed collapse. The prompt collapse could happen for very massive progenitors (mass greater than $`40M_{}`$) only if stellar winds were negligible, an assumption which appears to be not very reliable. In the delayed collapse all black holes were born with a mass equal to its minimum value, or just a little greater, and then slowly accrete matter up to their final mass. We have shown that the main burst of gravitational radiation is emitted when the black hole forms so that the gravitational energy spectrum is peaked in the range $`4.5÷9kHz`$, depending on the initial mass and the angular momentum of the black hole. Such frequencies do not match very well with the sensitivity curve of ground-based interferometers which reach their best sensitivity in the band $`60÷1000Hz`$. We have estimated the detectability of the emitted gravitational radiation in the delayed case, the most reliable scenario, and in the prompt case, for progenitor stars more massive than $`40M_{}`$. For each collapse model, we have considered different values of black hole minimum mass and angular momentum and different laws for progenitors mass distribution. We have also derived the theoretical black hole final mass distribution function.

Delayed collapses are detectable only inside the Local Group of galaxies by interferometers of the first generation. Obviously, delayed collapses are less detectable than the prompt ones, due to the lower matching of their characteristic frequencies to the sensitivity curve of the detectors. On the other hand, the initial stage of the delayed collapse, with the formation of a neutron star following a strongly asymmetric explosion (if $`m_{min}=2M_{}`$), could be a promising source of gravitational waves. Detection perspectives of delayed collapses are not so much better for adavnced interferometers because, if distances up to the Virgo Cluster are considered (where we expect $`1÷5ev/yr`$), the signal-to-noise ratio is much lower than one, unless a very high degree of asimmetry is produced. On the contrary, prompt collapses could be detected with large enough $`SNR`$. Delayed collapses to a black hole belong to a class of high frequency sources of gravitational waves, which comprises various processes involving compact objects, as for instance the excitation of neutron star w-modes ?, the coalescence of two neutron stars with the formation of a light black hole ?, dynamical instabilities in neutron stars ?, and some kinds of secular instability in neutron stars ?. Such sources cannot be efficiently detected by present resonant detectors and forthcoming interferometers because they are expected to emit at frequencies higher than those at which both interferometric and resonant detectors have their best sensitivity. In past years local arrays of small resonant detectors, which are particularly suited for the detection of high frequency gravitational radiation, have been proposed?,?. A detailed study of their detection performances, considering different geometries, dimensions and materials, has been done ?.

The non-continuous background of gravitational waves produced by the ensemble of the star collapses to a black hole, which occurred at a higher rate in the early phases of the Universe, has recently been calculated ?. It would be interesting to repeat the calculation in the case of delayed collapse. In such a case we have the superposition of energy spectra nearly with the same shape and all peaked at nearly the same frequency, in the range $`4.5÷9kHz`$ (this holds also for star metallicity $`Z=0`$). Roughly speaking, as most of the collapses take place at redshift $`z2`$, we expect to have, at the detector, a background spectrum strongly peaked somewhere in the band $`1.5÷3kHz`$.

Aknowledgements

I want to thank S. Frasca and M. A. Papa for the useful discussions and suggestions, T. Piran for his encouragement, V. Ferrari for her careful reading of the manuscript and the anonymous referees for their comments aiming to an improvement of this paper.

References


text:
# On the strength of the Kerr singularity and cosmic censorship

## I Introduction

The cosmic censorship hypothesis of Penrose  states that a physically realistic gravitational collapse can never result in a “naked” singularity—that is, all singularities arising in such situations should always be enclosed within an event horizon and hence invisible to distant observers. This hypothesis plays a fundamental role in the theory of black holes. Unfortunately, in spite of much sustained effort, no complete proof (or convincing counterexample) to cosmic censorship has been found as yet. Many partial results, however, have been established which considerably restrict the class of possible naked singularities. A cosmic censorship theorem of this type has previously been proved by one of the present authors . The aim of this paper is to answer a certain essential question closely related to that result.

There exist many exact solutions of the Einstein equations in which naked singularities do occur. Among these solutions the most important one, due to its direct astrophysical applications, is undoubtedly that of Kerr . This solution depends on two parameters $`m`$ and $`a`$ (see below), and represents the exterior gravitational field of a rotating body with mass $`m`$ and angular momentum $`ma`$, as measured from infinity in geometrized units with $`c=G=1`$. As is well known , the maximal analytic extension of the Kerr solution with $`m>0`$ and $`a0`$ contains the ring curvature singularity that may be interpreted as the outcome of collapse of a rotating object. When $`|a|m`$, this singularity is always hidden behind an event horizon. But if $`|a|>m`$, there is no event horizon and the singularity is visible for all observers; moreover, there are closed timelike curves through every point of the spacetime . Clearly, the Kerr solution is highly idealized and the singularity with $`|a|>m`$ cannot be a counterexample to Penrose’s hypothesis. However, it is not unlikely that the pathologies occurring in the case $`|a|>m`$ could also arise in more general scenarios of the collapse of a rapidly rotating star.

The censorship theorem of Ref.  shows that, under certain reasonable assumptions, a physically realistic collapse developing from a regular initial state cannot lead to the formation of a final state resembling the Kerr solution with $`|a|>m`$ —i.e. of a naked singularity accompanied by closed timelike curves. An important role in this result plays a certain inextendibility condition, which is assumed to hold for all (achronal) null geodesics terminating at the Kerr-like naked singularity under consideration. This condition characterizes the curvature strength of the singularity; roughly speaking, it holds for a given null geodesic $`\lambda `$ if the curvature near the singularity is strong enough so that at least one irrotational congruence of Jacobi fields along $`\lambda `$ is forced to refocus as $`\lambda `$ approaches the singularity (for more details on this condition, see Refs. ).

It is well known  that the inextendibility condition will always hold for achronal null geodesics terminating at the so-called strong curvature singularities defined by Tipler  (see below). These singularities have the property that all objects approaching them are crushed to zero volume. One often assumes that all singularities arising in physically realistic collapse should be of the strong curvature type (see, e.g., Refs. \[9–11\]). However, Israel  has suggested that the Kerr singularity fails to be strong in the sense of Tipler’s definition, for it tends to cause repulsive effects. Since the inextendibility condition is similar in spirit to Tipler’s definition (in both cases some Jacobi fields are refocused), the question one immediately asks is whether this condition can still be expected to hold for null geodesics terminating at the Kerr singularity—i.e. whether the censorship theorem of Ref.  can be applied to singularities of this type. In this paper we will obtain a result that provides some positive answer to this question. Namely, we will show that, in general, all null geodesics reaching the Kerr singularity with $`m>0`$ and $`a0`$ do in fact terminate, contrary to the suggestion of Israel, at Tipler’s strong curvature singularity.

## II Preliminaries

To begin with, we need to recall some of the basic facts on the Kerr solution. In Boyer and Lindquist coordinates $`(r,\theta ,\varphi ,t)`$ the Kerr metric is given by (cf. Ref. , p. 161):

$$ds^2=\rho ^2\left(\frac{dr^2}{\mathrm{\Delta }}+d\theta ^2\right)+(r^2+a^2)\mathrm{sin}^2\theta d\varphi ^2dt^2+\frac{2mr}{\rho ^2}(a\mathrm{sin}^2\theta d\varphi dt)^2,$$
(1)
where $`\rho ^2r^2+a^2\mathrm{cos}^2\theta `$ and $`\mathrm{\Delta }r^22mr+a^2`$. As mentioned earlier, the constant $`m`$ represents the mass of the metric source while $`a`$ is its angular momentum per unit mass. The Kerr spacetime is stationary and axisymmetric, with Killing vector fields $`\xi ^a(/t)^a`$ and $`\omega ^a(/\varphi )^a`$. Moreover this spacetime is Ricci flat and is of Petrov type D. The ring singularity is located in the equatorial plane $`\theta =\pi /2`$ at points where $`r=0`$. As shown by Carter , the only null geodesics which can reach this singularity are those lying strictly in the equatorial plane on the positive $`r`$ side. The equations of motion for these geodesics are (cf. Ref. , p. 328):

$$u^r\frac{dr}{ds}=\pm \left[E^2+\frac{2m}{r^3}\left(LaE\right)^2\frac{1}{r^2}\left(L^2a^2E^2\right)\right]^{1/2},$$
(2)
$$u^\varphi \frac{d\varphi }{ds}=\frac{1}{\mathrm{\Delta }}\left[\frac{2ma}{r}E+\left(1\frac{2m}{r}\right)L\right],$$
(3)
$$u^t\frac{dt}{ds}=\frac{1}{\mathrm{\Delta }}\left[\left(r^2+a^2+\frac{2a^2m}{r}\right)E\frac{2ma}{r}L\right],$$
(4)
where $`u^r`$, $`u^\varphi `$ and $`u^t`$ are the coordinate basis components of the tangent vector, $`𝐮`$, to a given null geodesic $`\lambda (s)`$ parametrized by an affine parameter $`s`$. (Note that the corresponding component $`u^\theta d\theta /ds`$ of $`𝐮`$ identically vanishes as $`\lambda (s)`$ lies in the equatorial plane.) The quantities $`E`$ and $`L`$ are the constants of the motion associated with the Killing vectors $`\xi ^a`$ and $`\omega ^a`$; they are defined as follows: $`E\xi _au^a`$ and $`L\omega _au^a`$. Physically, $`E`$ and $`L`$ can be interpreted, respectively, as the energy at infinity and the angular momentum about the symmetry axis, $`\theta =0`$, of a photon moving along $`\lambda (s)`$. It is also worth noting here that if $`L=aE`$, then $`\lambda (s)`$ must belong to one of the two principal null congruences associated with the algebraic type D of the solution (see Ref. , p. 329). The $`\pm `$ signs in Eq. (2) correspond to outgoing and ingoing geodesics, respectively.

Let us also recall that, according to Tipler’s definition , an affinely parametrized null geodesic $`\lambda (s)`$ is said to terminate in a strong curvature singularity at affine parameter value $`s_0`$ if the following holds (cf. Ref. , p. 160): Let $`\mu (s)`$ be a 2-form on the normal space to the tangent vector to $`\lambda (s)`$ determined by two linearly independent vorticity-free Jacobi fields $`𝐙_1(s)`$ and $`𝐙_2(s)`$ along $`\lambda (s)`$, i.e. $`\mu (s)𝐙_1𝐙_2`$. For all $`\mu (s)`$ that vanish for at most finitely many $`s`$ in some neighbourhood $`(s_0,s_1]`$ of $`s_0`$, we have $`lim_{ss_0}\mu (s)=0`$. Very useful criteria for determining whether a given null geodesic terminates in a strong curvature singularity have been found by Clarke and Królak . One of these criteria, which will be used in proving our result, can be formulated as the following

Proposition 2.1. Let $`\lambda (s)`$ $`(0<ss_1)`$ be an affinely parametrized null geodesic and let $`\{𝐄_i\}`$ $`(i=1,2,3,4)`$ be a pseudo-orthonormal tetrad parallely propagated along $`\lambda (s)`$, with $`𝐄_1𝐄_1=𝐄_2𝐄_2=𝐄_3𝐄_4=𝐄_4𝐄_3=1`$, all other scalar products vanishing and $`𝐄_4=𝐮`$, where $`𝐮`$ is the tangent vector to $`\lambda (s)`$. Let $`C^m_{4n4}`$ $`(m,n\{1,2\})`$ be some component of the Weyl tensor with respect to the tetrad $`\{𝐄_i\}`$. If there exists some affine parameter value $`s_0(0,s_1]`$ such that $`C^m{}_{4n4}{}^{}Ks^2`$ on $`(0,s_0]`$, where $`K`$ is some positive constant, then $`\lambda (s)`$ terminates in Tipler’s strong curvature singularity as $`s0`$.

Proof. The proof follows immediately from Proposition 7 of Ref. .

## III The main result

We are now in a position to state and prove our main result.

Theorem 3.1. Let $`\lambda (s)`$ $`(s>0)`$ be an affinely parametrized null geodesic in the Kerr spacetime with $`m>0`$ and $`a0`$. Suppose that $`\lambda (s)`$ approaches the ring singularity as $`s0`$. If $`\lambda (s)`$ does not belong to either of the two principal null congruences, then $`\lambda (s)`$ terminates in Tipler’s strong curvature singularity as $`s0`$.

In very brief outline, the proof of our result runs as follows. We first construct a certain pseudo-orthonormal tetrad $`\{𝐄_i\}`$ parallely propagated along the geodesic $`\lambda (s)`$. We next find one of the components of the Weyl tensor with respect to the tetrad $`\{𝐄_i\}`$. It turns out that if $`\lambda (s)`$ is not a member of the principal null congruences, then this component must grow along $`\lambda (s)`$ as fast as $`r^5`$, where $`r`$ is the Boyer-Lindquist radial coordinate on $`\lambda (s)`$. Using Eq. (2), we then show that this tetrad component of the Weyl tensor must diverge along $`\lambda (s)`$ at least as fast as $`s^2`$, where $`s`$ is the affine parameter. This implies, by Proposition 2.1, that $`\lambda (s)`$ must then terminate in Tipler’s strong curvature singularity. The rigorous proof is as follows.

Proof. Let us assume that the spacetime is parametrized by the Boyer-Lindquist coordinates $`(r,\theta ,\varphi ,t)`$. Since $`\lambda (s)`$ reaches the ring singularity at $`r=0`$, it must lie entirely in the equatorial plane $`\theta =\pi /2`$ on the positive $`r`$ side. The coordinate $`r`$ will change along $`\lambda (s)`$ according to Eq. (2). From this equation it follows easily that there may exist at most two positive values of $`r`$ on $`\lambda (s)`$ for which $`dr/ds=0`$, and so we will have $`dr/ds0`$ along $`\lambda (s)`$ for all $`r`$ in a sufficiently small interval about 0. It is thus clear that there must exist some affine parameter value $`s_1>0`$ such that $`dr/ds>0`$ along $`\lambda (s)`$ for all $`s(0,s_1]`$, because $`r>0`$ and $`s>0`$ on $`\lambda (s)`$, and $`r0`$ along $`\lambda (s)`$ as $`s0`$. Let us now fix such a parameter value $`s_1`$. According to Proposition 2.1, in order to prove the theorem, it suffices to show that the rate of growth of the curvature along $`\lambda (s)`$ is strong enough for $`s`$ in some small interval about 0. For convenience, and without loss of generality, we can thus assume that the affine parameter on $`\lambda (s)`$ ranges over the interval $`(0,s_1]`$.

Let $`𝐮`$ be the tangent vector to $`\lambda (s)`$. The coordinate basis components $`u^r`$, $`u^\varphi `$ and $`u^t`$ of $`𝐮`$ are given by Eqs. (2)-(4), with the $`+`$ sign in Eq. (2) as $`dr/ds>0`$ along $`\lambda (s)`$. The corresponding component $`u^\theta `$ of $`𝐮`$ identically vanishes as $`\lambda (s)`$ lies entirely in the equatorial plane. Given any point $`p\lambda (s)`$, we define $`𝐤`$ to be the vector in the tangent space $`T_p`$ with the following coordinate basis components: $`k^r=k^\varphi =k^t=0`$ and $`k^\theta =r^1`$, where $`r`$ is the radial coordinate of $`p`$. One can readily verify that, with respect to the scalar product given by the Kerr metric (1), $`𝐤`$ is a unit spacelike vector orthogonal to $`𝐮`$, i.e. $`𝐤𝐤=1`$ and $`𝐤𝐮=0`$. In addition, $`𝐤`$ is parallely transported along $`\lambda (s)`$ (see Appendix A). The pair $`\{𝐤,𝐮\}`$ can easily be extended to the pseudo-orthonormal tetrad $`\{𝐄_i\}`$ mentioned in Proposition 2.1. To see this, let us first fix some point $`q\lambda (s)`$. Let us now take a unit spacelike vector orthogonal to the plane spanned by the vectors $`𝐮`$ and $`𝐤`$ in the tangent space $`T_q`$; we will denote this vector by $`𝐧`$. (Clearly, such a vector can always be found since $`𝐮`$ is null and $`𝐤`$ is spacelike.) Now let $`\{𝐞_i\}`$ $`(i=1,2,3,4)`$ be some orthonormal basis in the space $`T_q`$, with the spacelike vectors $`𝐞_2`$ and $`𝐞_3`$ chosen so that $`𝐞_2=𝐤`$ and $`𝐞_3=𝐧`$. Using the spacelike vector $`𝐞_1`$ and the timelike vector $`𝐞_4`$, we now define $`𝐯_1\alpha (𝐞_4+𝐞_1)`$ and $`𝐯_2\alpha (𝐞_4𝐞_1)`$, where $`\alpha 0`$ is some constant. Evidently, these new vectors are null and each of them is orthogonal to both $`𝐤`$ and $`𝐧`$. Moreover, as $`𝐯_1𝐯_2=2\alpha ^20`$, at least one of them, say $`𝐯_1`$, is not parallel to the vector $`𝐮`$, i.e. we must have $`𝐮𝐯_10`$. By suitable choice of the constant $`\alpha `$ in the definition of $`𝐯_1`$ one can always normalize $`𝐯_1`$ so that $`𝐮𝐯_1=1`$. The pseudo-orthonormal tetrad $`\{𝐄_i\}`$ in the space $`T_q`$ can now be chosen as follows: $`𝐄_1𝐧`$, $`𝐄_2𝐤`$, $`𝐄_3𝐯_\mathrm{𝟏}`$ and $`𝐄_4𝐮`$. By parallely transporting this tetrad along $`\lambda (s)`$ one obtains the pseudo-orthonormal basis at each point of $`\lambda (s)`$.

We shall now find one of the components $`C^i_{jkl}`$ of the Weyl tensor with respect to the tetrad $`\{𝐄_i\}`$. The components $`C_{mjkl}`$ of the Weyl tensor with respect to $`\{𝐄_i\}`$ can be obtained from the coordinate components $`C_{abcd}`$ of the Weyl tensor according to

$$C_{mjkl}=C_{abcd}E_m^aE_j^bE_k^cE_l^d,$$
(5)
where $`E_m^a`$ (resp., $`E_j^b`$, $`E_k^c`$ and $`E_l^d`$) is the $`a`$th (resp., $`b`$th, $`c`$th and $`d`$th) coordinate basis component of the vector $`𝐄_m`$ (resp., $`𝐄_j`$, $`𝐄_k`$ and $`𝐄_l`$) of the tetrad $`\{𝐄_i\}`$. Having the tetrad components $`C_{mjkl}`$, we can now easily find the tetrad components $`C^i_{jkl}`$ of the Weyl tensor:

$$C^i{}_{jkl}{}^{}=\eta ^{im}C_{mjkl},$$
(6)
where $`\eta ^{im}`$ is the inverse of the matrix $`\eta _{im}𝐄_i𝐄_m`$; that is, we have $`\eta ^{11}=\eta ^{22}=\eta ^{34}=\eta ^{43}=1`$ and $`\eta ^{im}=0`$ in all other cases. Since the tetrad $`\{𝐄_i\}`$ is chosen so that $`𝐄_2=𝐤`$ and $`𝐄_4=𝐮`$, and the vectors $`𝐤`$ and $`𝐮`$ are given in the explicit form, we can find, applying (5) and (6), an explicit expression for the component $`C^2_{424}`$ of the Weyl tensor with respect to $`\{𝐄_i\}`$. This is done in Appendix B; the result is

$$C^2{}_{424}{}^{}=\frac{3m(LaE)^2}{r^5}.$$
(7)
The task is now to show that if $`\lambda (s)`$ is not a member of the principal null congruences, then there exists some affine parameter interval $`(0,s_0]`$ of $`\lambda (s)`$ on which $`C^2{}_{424}{}^{}Ks^2`$, where $`K`$ is some positive constant. To do this, let us first recall that the affine parameter $`s`$ and the radial coordinate $`r`$ on $`\lambda (s)`$ are related by Eq. (2), with the + sign as $`dr/ds>0`$ on $`(0,s_1]`$. Let us now rewrite this equation in the form

$$\frac{ds}{dr}=\frac{r^{3/2}}{F(r)},$$
(8)
where $`F(r)\left[r^3E^2r(L^2a^2E^2)+2m(LaE)^2\right]^{1/2}`$. As $`dr/ds>0`$ on $`(0,s_1]`$, it is obvious that $`ds/dr>0`$ on $`(0,r_1]`$, where $`r_1`$ denotes the value of the coordinate $`r`$ of the point $`\lambda (s_1)`$. By (8) it is thus clear that $`F(r)`$ is strictly positive on $`(0,r_1]`$. It is also clear, as $`lim_{r0}F(r)=|LaE|\sqrt{2m}`$, that $`F(r)`$ is bounded on $`(0,r_1]`$. So there exists a positive number $`F_0sup\{F(r)|0<rr_1\}`$. Consider now the function $`y(r)2r^{5/2}(5F_0)^1`$ defined on $`(0,r_1]`$. Since $`dy/dr=r^{3/2}/F_0`$ and $`F_0F(r)>0`$ on $`(0,r_1]`$, by (8) we have $`ds/dr|_r^{}dy/dr|_r^{}>0`$ for all $`r^{}(0,r_1]`$. From this it follows easily, as $`lim_{r0}y(r)=lim_{r0}s(r)=0`$, that $`s(r)y(r)`$ for all $`r(0,r_1]`$. Evidently, both $`s(r)`$ and $`y(r)`$ are strictly increasing on $`(0,r_1]`$, and so there exist their inverse functions. Let $`r:(0,s_1](0,r_1]`$ be the inverse function of $`s(r)`$ and let $`z:(0,s_0](0,r_1]`$ be the inverse function of $`y(r)`$ (note that $`s_0s_1`$ since $`s(r)y(r)`$ on $`(0,r_1]`$). From the fact that $`s(r)y(r)`$ for all $`r(0,r_1]`$, it follows immediately that $`r(s)z(s)`$ for all $`s(0,s_0]`$. Note also that $`z(s)=(5F_0s/2)^{2/5}`$, which is clear from the definition of $`y(r)`$. We thus have $`r(s)(5F_0s/2)^{2/5}`$ on $`(0,s_0]`$. Combining this inequality with (7), and taking into account the fact that $`m>0`$, we obtain

$$C^2{}_{424}{}^{}=\frac{3m(LaE)^2}{r^5(s)}Ks^2$$
(9)
for all $`s(0,s_0]`$, where $`K12m(LaE)^2(5F_0)^2`$. Suppose now that $`\lambda (s)`$ does not belong to either of the two principal null congruences; then we have $`LaE`$ (see Sec. 2), and hence $`K>0`$. By (9) and Proposition 2.1 it is thus clear that $`\lambda (s)`$ must terminate in Tipler’s strong curvature singularity as $`s0`$, which is the desired conclusion.

## IV Concluding remarks

We have examined the curvature strength of the Kerr singularity with $`m>0`$ and $`a0`$. We have shown that every null geodesic reaching this singularity, with the exception of those belonging to the principal null congruences, must in fact terminate in Tipler’s strong curvature singularity. This is a typical property of null geodesics approaching the Kerr singularity because the principal null geodesics are very special and can be considered to form a set of “measure zero” in the family of all null geodesics reaching the singularity (Ref. , pp. 328, 329). The existence of these special geodesics is due to the high symmetry of the solution and one can expect that such geodesics will not occur in more general spacetimes. Thus if one attempts to define the curvature strength of a Kerr-like naked singularity, which could possibly arise in a generic collapse, one way of doing this is to assume that all null geodesics approaching such a singularity will behave much the same as a typical null geodesic approaching the Kerr singularity with $`|a|>m`$ —i.e. will also terminate in a singularity of the strong curvature type. Using this assumption, one may then attempt to formulate and prove a theorem which would constraint or prohibit the occurrence of Kerr-like naked singularities in generic collapse situations. A theorem of this type was established in Ref. .

In this context, it is worth recalling that the inextendibility condition assumed in the theorem of Ref.  may in fact hold for a much more general class of possible singularities than only those of the strong curvature type, for the curvature need not necessarily diverge along geodesics satisfying this condition . However, we have checked that this condition fails to hold for the principal null geodesics approaching the Kerr singularity (in fact, it will always fail to hold for principal null geodesics in any Ricci-flat spacetime). Thus the theorem of Ref.  does not exclude the possibility that a naked Kerr singularity accompanied by closed timelike curves could develop from some non-singular initial data. However, from the proof of that theorem it may be concluded that this can happen only if the inextendibility condition fails to hold for all null geodesics terminating in the past at the naked Kerr singularity and generating the future Cauchy horizon due to the formation of this singularity. According to our result, this is possible only if all these geodesics belong to the principal null congruences. This would be a very special case.

## Acknowledgments

We wish to thank the referee for helpful comments. This work was supported by the Polish Committee for Scientific Research (KBN) under grant No. 2 P03B 073 15.

## A

In this appendix, we demonstrate that the vector $`𝐤`$ is parallely transported along the null geodesic $`\lambda (s)`$ with the tangent vector $`𝐮`$. To this end we only need to show that

$$\left(\frac{k^c}{x^b}+\mathrm{\Gamma }^c{}_{ab}{}^{}k_{}^{a}\right)u^b=0,$$
(A1)
where $`\mathrm{\Gamma }^c_{ab}`$ are the connection coefficients, which can be obtained from the metric tensor $`g_{ab}`$ according to

$$\mathrm{\Gamma }^c{}_{ab}{}^{}=\frac{1}{2}g^{cd}(\frac{g_{bd}}{x^a}+\frac{g_{da}}{x^b}\frac{g_{ab}}{x^d}).$$
(A2)
Since $`u^\theta =k^r=k^\varphi =k^t=0`$ and $`k^\theta =r^1`$, the components of Eq. (A1) in the $`(r,\theta ,\varphi ,t)`$ coordinate system take the form
for $`c=r`$:

$$\left(\mathrm{\Gamma }^r{}_{\theta r}{}^{}u_{}^{r}+\mathrm{\Gamma }^r{}_{\theta \varphi }{}^{}u_{}^{\varphi }+\mathrm{\Gamma }^r{}_{\theta t}{}^{}u_{}^{t}\right)r^1=0,$$
(A3)
for $`c=\theta `$:

$$r^2u^r+\left(\mathrm{\Gamma }^\theta {}_{\theta r}{}^{}u_{}^{r}+\mathrm{\Gamma }^\theta {}_{\theta \varphi }{}^{}u_{}^{\varphi }+\mathrm{\Gamma }^\theta {}_{\theta t}{}^{}u_{}^{t}\right)r^1=0,$$
(A4)
for $`c=\varphi `$:

$$\left(\mathrm{\Gamma }^\varphi {}_{\theta r}{}^{}u_{}^{r}+\mathrm{\Gamma }^\varphi {}_{\theta \varphi }{}^{}u_{}^{\varphi }+\mathrm{\Gamma }^\varphi {}_{\theta t}{}^{}u_{}^{t}\right)r^1=0,$$
(A5)
for $`c=t`$:

$$\left(\mathrm{\Gamma }^t{}_{\theta t}{}^{}u_{}^{r}+\mathrm{\Gamma }^t{}_{\theta \varphi }{}^{}u_{}^{\varphi }+\mathrm{\Gamma }^t{}_{\theta t}{}^{}u_{}^{t}\right)r^1=0.$$
(A6)
Applying (A2) to the Kerr metric (1), we can now calculate the connection coefficients appearing in Eqs. (A3)-(A6); the result is

$$\mathrm{\Gamma }^r{}_{\theta r}{}^{}=\frac{a^2\mathrm{sin}\theta \mathrm{cos}\theta }{\rho ^2},\mathrm{\Gamma }^\theta {}_{\theta r}{}^{}=\frac{r}{\rho ^2},\mathrm{\Gamma }^\varphi {}_{\theta \varphi }{}^{}=\frac{(\rho ^2+2a^2mr\mathrm{sin}^2\theta )\mathrm{cos}\theta }{\rho ^4\mathrm{sin}\theta },$$

$$\mathrm{\Gamma }^\varphi {}_{\theta t}{}^{}=\frac{2mar\mathrm{cos}\theta }{\rho ^4\mathrm{sin}\theta },\mathrm{\Gamma }^t{}_{\theta \varphi }{}^{}=\frac{2ma^3r\mathrm{sin}^3\theta \mathrm{cos}\theta }{\rho ^4},\mathrm{\Gamma }^t{}_{\theta t}{}^{}=\frac{2ma^2r\mathrm{sin}\theta \mathrm{cos}\theta }{\rho ^4},$$

$$\mathrm{\Gamma }^r{}_{\theta \varphi }{}^{}=\mathrm{\Gamma }^r{}_{\theta t}{}^{}=\mathrm{\Gamma }^\theta {}_{\theta \varphi }{}^{}=\mathrm{\Gamma }^\theta {}_{\theta t}{}^{}=\mathrm{\Gamma }^t{}_{\theta r}{}^{}=\mathrm{\Gamma }^r{}_{\theta \varphi }{}^{}=0.$$

Substituting these coefficients into Eqs. (A3)-(A6), and putting $`\theta =\pi /2`$ as $`\lambda (s)`$ lies in the equatorial plane, we can now readily see that Eqs. (A3)-(A6) are satisfied, as desired.

## B

We give in this appendix some details of the computation of the component $`C^2_{424}`$ of the Weyl tensor with respect to the tetrad $`\{𝐄_i\}`$. Since $`𝐄_2=𝐤`$ and $`𝐄_4=𝐮`$, and $`u^\theta =k^r=k^\varphi =k^t=0`$ and $`k^\theta =r^1`$, the expression (5) for $`C_{2424}`$ takes the form

$$C_{2424}=r^2\left[C_{\theta r\theta r}(u^r)^2+C_{\theta \varphi \theta \varphi }(u^\varphi )^2+C_{\theta t\theta t}(u^t)^2\right]$$

$$+2r^2\left(C_{\theta r\theta \varphi }u^ru^\varphi +C_{\theta r\theta t}u^ru^t+C_{\theta \varphi \theta t}u^\varphi u^t\right),$$
(B1)
where the components $`u^r`$, $`u^\varphi `$ and $`u^t`$ of $`𝐮`$ are given by Eqs. (2)-(4), with the + sign in Eq. (2). We recall that the Weyl tensor $`C_{abcd}`$ is defined by

$$C_{abcd}=R_{abcd}+g_{a[d}R_{c]b}+g_{b[c}R_{d]a}+\frac{1}{3}Rg_{a[c}g_{d]b},$$

where $`R_{abcd}`$, $`R_{ab}`$ and $`R`$ denote the curvature tensor, the Ricci tensor and the curvature scalar, respectively. Since the Kerr spacetime is Ricci flat, $`R_{ab}`$ and $`R`$ will vanish, and so $`C_{abcd}=R_{abcd}`$. In this case the coordinate components of the Weyl tensor can be obtained from the metric components $`g_{ab}`$ according to

$$C_{abcd}=\frac{1}{2}(\frac{^2g_{ac}}{x^dx^b}+\frac{^2g_{bd}}{x^cx^a}\frac{^2g_{bc}}{x^dx^a}\frac{^2g_{ad}}{x^cx^b})+g_{ef}(\mathrm{\Gamma }^e{}_{ca}{}^{}\mathrm{\Gamma }_{}^{f}{}_{bd}{}^{}\mathrm{\Gamma }^e{}_{da}{}^{}\mathrm{\Gamma }_{}^{f}{}_{bc}{}^{}),$$

where the connection coefficients are given by (A2). Applying this formula to the Kerr metric (1), we can now calculate the coordinate components of the Weyl tensor appearing in (B1); the result is

$$C_{\theta r\theta r}=\frac{mr(3a^2\mathrm{cos}^2\theta r^2)}{\rho ^2\mathrm{\Delta }},$$

$$C_{\theta \varphi \theta \varphi }=\frac{mr(3a^2\mathrm{cos}^2\theta r^2)[a^2\mathrm{\Delta }\mathrm{sin}^2\theta +2(r^2+a^2)^2]\mathrm{sin}^2\theta }{\rho ^6},$$

$$C_{\theta t\theta t}=\frac{mr(3a^2\mathrm{cos}^2\theta r^2)(2a^2\mathrm{sin}^2\theta +\mathrm{\Delta })}{\rho ^6},$$

$$C_{\theta \varphi \theta t}=\frac{mra(3a^2\mathrm{cos}^2\theta r^2)[\mathrm{\Delta }+2(r^2+a^2)]\mathrm{sin}^2\theta }{\rho ^6},$$

$$C_{\theta r\theta \varphi }=C_{\theta r\theta t}=0.$$

Inserting these expressions in (B1), and putting $`\theta =\pi /2`$ as $`\lambda (s)`$ lies in the equatorial plane, we get

$$C_{2424}=\frac{3m(LaE)^2}{r^5}.$$

Finally, we note that $`C_{2424}=C^2_{424}`$, which is clear from (6).