Docstoc

ABSTRACT

Document Sample
ABSTRACT Powered By Docstoc
					                                 ABSTRACT


The theory of General Relativity has been in existence for 90 years and has stood
up to all tests it has been subjected to in that time. The PPN parameter γ is
a measure of the accuracy of theories of gravity and assumes different values in
different theories. By measuring the Shapiro time delay of light it is possible to
constrain γ and thereby constrain gravitational theories. This Shapiro time delay
can be measured in our solar system but it is only in the vicinity of extremely
compact objects such as pulsars and black holes that it can be tested under the
immense gravitational fields that can only be found there. A pulsar in a binary
orbit about another compact object is the ideal system in which to test this
effect. In this work we have gone from Kepler’s laws of simple planetary motion
to deriving the equations that explain binary orbits to incorporating General
Relativity into these equations in order to obtain the equations for relativistic
particle orbits. We then evolved this theory even further so as to be able to
explain relativistic light ray orbits and then used this knowledge to model the
Shapiro delay in a binary system. With a working model it became possible to
simulate the Shapiro delay in a wide range of possible systems and then to use
these simulations to say something about what type of system should be focussed
on in future so as to measure the Shapiro delay and thereby constrain more tightly
the parameter γ.
                                  CONTENTS


1. Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                                          2

2. The   PPN Formalism . . . . . . . . . . . . . . . .                     .   .   .   .   .   .   .   .   .   .   .   .    8
   2.1   Background . . . . . . . . . . . . . . . . . .                    .   .   .   .   .   .   .   .   .   .   .   .    8
   2.2   The 10 PPN parameters and their origin . .                        .   .   .   .   .   .   .   .   .   .   .   .    9
   2.3   The Nordtvedt effect . . . . . . . . . . . . .                     .   .   .   .   .   .   .   .   .   .   .   .   12
   2.4   Variation of Newton’s constant . . . . . . .                      .   .   .   .   .   .   .   .   .   .   .   .   14
         2.4.1 Spin tests . . . . . . . . . . . . . . .                    .   .   .   .   .   .   .   .   .   .   .   .   14
         2.4.2 Orbital decay tests . . . . . . . . . .                     .   .   .   .   .   .   .   .   .   .   .   .   15
         2.4.3 Changes in the Chandrasekhar mass                           .   .   .   .   .   .   .   .   .   .   .   .   15

3. Orbital mechanics . . . . . . . . . .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   17
   3.1 Orbital geometry . . . . . . .      .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   17
        3.1.1 Orbital parameters . .       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   17
   3.2 Keplerian mechanics . . . . .       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   21
   3.3 Binary pulsar evolution . . . .     .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   26
        3.3.1 Evolutionary processes       .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   26

4. Topics from the Theory of Curved Space              .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   28
   4.1 Curvilinear coordinates . . . . . . .           .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   28
   4.2 The metric . . . . . . . . . . . . . .          .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   30
   4.3 Geodesics . . . . . . . . . . . . . .           .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   31

5. Relativistic particle orbits in a Schwarzschild             Geometry                .   .   .   .   .   .   .   .   .   34
   5.1 Physical and Mathematical concepts .                    . . . . . .             .   .   .   .   .   .   .   .   .   34
        5.1.1 Gravitational Redshift . . . . .                 . . . . . .             .   .   .   .   .   .   .   .   .   35
   5.2 Derivation of relativistic particle orbits              . . . . . .             .   .   .   .   .   .   .   .   .   36

6. Shapiro delay . . . . . . . . . . . . . . . . . . . . . . . . .                             . . .       .   .   .   .   40
   6.1 Background . . . . . . . . . . . . . . . . . . . . . . .                                . . .       .   .   .   .   40
   6.2 Mathematical derivation for the Shapiro time delay of                                   light       .   .   .   .   41
       6.2.1 Deflection angle . . . . . . . . . . . . . . . . .                                 . . .       .   .   .   .   44
       6.2.2 The Shapiro delay of light . . . . . . . . . . .                                  . . .       .   .   .   .   46

7. Strong-field Shapiro delay . . . . . . . . . . . . . . . . . . . . . . . . .                                             48
Contents                                                                                                                   1

8. Results of Simulations . . . .    . . . .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   54
   8.1 Procedure followed . . .      . . . .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   54
   8.2 Graphs and Results . . .      . . . .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   55
   8.3 Predicted Future Results      . . . .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   67
   8.4 3-Dimensional Parameter       spaces    .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   74

9. Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . .                                              82

10. References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                                         85


Appendix                                                                                                                   88

A. The connection or covariant derivative      . . . . . . . . . . . . . . . . .                                           89
   A.1 Torsion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                                           90
   A.2 Curvature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                                             91

B. MatLab code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .                                             94
                             1. INTRODUCTION


The process of stellar evolution is a long and complicated one with numerous
factors influencing the final state of a star. Generally the death of stars with a
birth-mass below about 8 Solar masses (M ) is a relatively gentle affair. They
will expand to many times their original diameter and go through a red giant
phase before the core, with a left-over mass below the Chandrasekhar mass of
1.4M , collapses to form a White Dwarf that is sustained against further collapse
by electron degeneracy pressure. This is the pressure created through the Pauli
exclusion principle preventing 2 electrons from occupying the same energy state
and the highest energy level that is filled is referred to as the Fermi energy.
    Those stars of birth-mass greater than 8M , come to a far more spectacular
end. These stars will explode as supernovae, releasing in an instant up to 1044 J
of energy and having a luminosity of billions of times that of the sun. Then,
as before, the core re-collapses and in the event of the core still having a mass
greater than 2.5M , although this limiting mass is still a matter of debate, there
is no degeneracy pressure that can retard the massive gravitational collapse and
it collapses to form a black hole.
    It is those stars whose cores have a mass of between 1.4M and 2.5M
that, upon collapse, form an extremely interesting class of astrophysical object
known as a neutron star. In 1932 the English physicist, Sir James Chadwick,
discovered the neutron, for which he won the Nobel Prize in 1935, and unwit-
tingly started a cascade of discoveries and theories in astrophysics that have
lasted to the present day. Two years after this discovery Walter Baade and Fritz
Zwicky1 surmised that a supernova explosion in a star that has insufficient mass
for complete gravitational collapse into a black hole, could be sustained by neu-
tron degeneracy pressure against further collapse. This is analogous to electron
degeneracy pressure that sustains white dwarf stars against collapse. Under or-
dinary circumstances the nuclear reactions of neutron and neutrino formation,
p+ + e−     n + ν, tend toward the production of proton-electron pairs but under
the massive self-gravitation of a forming neutron star, the forward reaction dom-
inates and neutron formation is the primary reaction. This is because, as central
density increases, the electron Fermi energy always increases up to the point
where inverse beta decay takes place and forces the electrons into the nuclei. The
  1
    Baade, W.,Zwicky, F., “Cosmic Rays from Super-novae”, Proc. Nat. Acad. Sci., 20, 254,
(1934)
                                                                                         3


reverse reaction barely takes place because all of the electron states into which
the electrons could move are already occupied, preventing decay of the neutrons
in the nucleus into protons and electrons. As a result these objects would consist
primarily of neutron matter and have a very small radius and enormous densities.
The likelihood of finding such a star, however, was deemed virtually nonexistent
due to the high probability that such a star would be extremely small and emit
very little visible light. The luminosity function is given by L ∝ R2 T 4 and the
predicted radius of such a star is approximately 10km. When this is compared
to a 1 solar mass star’s radius of approximately 700000km, the difference when
these values are squared is of the order of 109 . In mitigation of this the neutron
star’s temperature is also about 100 times greater than that of the sun and, as
the luminosity function contains a T 4 term, the loss in luminosity due to de-
creased radius would largely be compensated for by the gain due to increased
temperature. The fundamental difference is that, as a blackbody, the neutron
star’s radiation will peak short of the optical band and as such be scarcely visible
in the optical wavelengths. Theoretical debate on the properties of neutron stars
continued unabated as numerous great minds tried to elucidate their properties
from first principles.
    In 19672 , Antony Hewish and his research student, Jocelyn Bell, were un-
dertaking an interplanetary scintillation study using a large antenna that was
sensitive to relatively long, 3.7 m, wavelengths. Within a month Miss Bell had
detected a large periodic fluctuation that could not be scintillation. After further
observations it was realised that the signal had a 1.337 s period and thus could
not be a random effect. Initially it was feared/hoped to be a signal from an
extra-terrestrial source and was given the designation LGM, standing for Little
Green Men, but this thought was quickly laid to rest when the signal showed
no variation in period whatsoever. The link between this pulsating object and
neutron stars was not immediately made, however, and even once the association
had been made, it was not readily accepted. Although both Franco Pacini3 and
Thomas Gold4 published papers linking rapidly rotating neutron stars to observ-
able radio pulses, competing theories of rotating or pulsating white dwarf stars
remained. It took a year for the matter to largely be put to rest with the dis-
covery of the 89 ms Vela pulsar5 . This settled one side of the argument as it was
realised that a white dwarf star rotating at that speed would quickly be ripped
apart by the centrifugal forces created. Only a star with the small diameter of
the predicted neutron star could sustain such periodicities. With the discovery
  2
     Hewish, A., Bell, J., et al., “Observations of a Rapidly Pulsating Radio Source”, Na-
ture, 217, 709, (1968)
   3
     Pacini, F., “Energy Emission from a Neutron Star”, Nature, 216,567 , (1967)
   4
     Gold, T.,“Rotating Neutron Stars as the Origin of the Pulsating Source of Radio”, Na-
ture, 218,731 , (1968)
   5
     Large, M.I., Vaughan, A.E., Mills, B.Y.,“A Pulsar Supernova Association”, Nature, 220,
340-341, (1968)
                                                                                         4


that the periodicity was gradually decreasing, the last doubt was laid to rest as
only a rotating object would slow down as it lost energy while a vibrating one
would not. At last this class of object could be categorised and became known
to the world as pulsars with general designation of PSR, or pulsating source of
radio.
    The extremely large mass of neutron stars concentrated in such a small volume
means that the gravitational fields generated by these objects are just slightly
less than that of the super-dense black holes that general relativists have for so
long sought and studied. The obvious, and very large, advantage of neutron
stars is that they are visible in radio and occasionally in optical frequencies and
changes in behaviour can therefore be observed. A neutron star in a binary
system is a phenomenal testing ground for gravitational theories as the observed
electromagnetic beam generated just above the surface of the star will often pass
through the gravitational field generated by a binary companion. Should this
companion be another neutron star, or even better a black hole, the beam would
pass through extremely strong gravitational fields and thus strong field tests of
gravity could be performed that 30 years ago would have been mere speculation.
When the first double neutron star system was discovered by Russell Hulse and
Joseph Taylor (1975)6 , such tests became a reality and have continued unabated
with ever larger and more accurate data sets. The discovery of a double pulsar
system by Lyne, Burgay, Kramer et al. (2004)7 , was yet another step closer to an
understanding of matter under these extreme conditions. The ability to observe
the gravitational effects of two pulsars as their beams pass through one another’s
gravity is paramount to deepening our understanding of the manner in which
space-time is curved.
    It is clear that pulsars are observed through reception of the signal on earth
and the subsequent analysis of the signal and any changes contained therein.
In practice, however, the process is not quite so simple. Generally the pulse
of radiation received on earth is very small and is swallowed up in the noise of
the electronics and even the noise of the cosmic microwave background (CMB).
The method used to extract the information of the pulse is that of folding, or
integrating, the signal until the peak emerges from the noise. In a pulsar with a
known frequency, a long observation time of the pulsar is cut into segments of one
pulsar period’s length and all of these segments are added together so that the
pulse emerges from the noise. This is possible because the noise is stochastic and
does not add to any significant degree while the pulse intensity adds up over each
folding. In pulsar searches, the signal is folded over using a series of frequencies
until the peak emerges for a given frequency. The correct period can then be
found by refining this process to obtain the maximum signal to noise ratio for
  6
    Hulse, R., Taylor, J.,“Discovery of a Pulsar in a Binary System”, Ap. J., 195, L51-L53,
(1975)
  7
    Lyne, A.G., Burgay, M., Kramer, M., et al.,“”, Science, 303, 1153, (2004)
                                                                                          5


the signal.
    A complicating effect in this process is known as the Roemer delay. This effect
occurs in signals originating from pulsars that are in binary orbits. As the pulsar
itself moves through space, traversing periodically away from and then towards
the earth, the signal originating from the pulsar has to travel a varying distance
to our earth-based detectors. The result of this is a signal that varies with a
very similar periodicity to that of the Shapiro delay we are trying to measure.
Fortunately the Shapiro delay in the strong field regime has been measured so all
is not lost. In a system that is correctly oriented so as to make the Shapiro delay
a maximum, the delay can be measured and this has been done in the highly
relativistic binary system PSR J0737-30398 . So even though the effect is difficult
to detect, discovery of more relativistic systems and a better understanding of
the Shapiro delay should assist in measuring this effect more easily.
    This effect also occurs due to the earth’s orbit about the sun but fortunately
the orbital periods of the systems in which we are interested in this thesis are
so short that the motion of the earth during that time can be neglected. The
calculation of the magnitude of this effect is a fairly straightforward process that
requires only some simple orbital mechanics. This calculation will be dealt with
in a later chapter.
    In pulsar timing one also encounters the problem of dispersion. This is the
phenomenon where the higher frequency signals travel faster through the inter-
stellar medium and therefore arrive on earth earlier than lower frequency bands
resulting in a ’smearing out’ of the pulsar signal. For this reason de-dispersion
is necessary to obtain the true information from a signal by separating the in-
coming signal into frequency bands and compensating for the dispersive effect in
each band. There are a myriad of other effects that pulsar astronomers are aware
of and take into account when timing pulsars but these are not the topic of this
thesis and as such will not be discussed in detail.
    Suffice it to say that through the application of various tools it is possible to
measure pulse arrival times to an accuracy of up to 10−15 s and therefore even
minuscule variations in the arrival times can be observed and analysed as caused
by possible physical effects occurring in and around the neutron star. These
effects can have their source in the interior of the star, its crust, the surrounding
magnetosphere, the gravitational well generated by these strongly self-gravitating
bodies or a combination of these and only through continued observation and
study can an improved understanding of these dynamic objects be reached.
    Gravity is one of those aspects for which a better understanding is desired.
The theory of General Relativity (GR) formulated by Einstein in 19169 is still a
frontrunner in the race to understand gravity but further testing of GR and its
  8
    Burgay, M., D’Amico, N., Possenti, A., et al.,“The Highly Relativistic Binary Pulsar PSR
J0737-3039: Discovery and Implications”, ASP Conf. series, 328, 53, (2005)
  9
                                                           a
    Einstein, A., “Die Grundlage der Algemeinen Relativit¨tstheorie”, Ann. der Physik, 49,
769-822, (1916)
                                                                                                  6


competitors is essential in laying the matter to rest. The Parameterised Post-
Newtonian formalism was thus developed in the 1970’s10,11,12,13,14,15 as a frame-
work in which to compare theories of gravity through their effects on a set of
10 parameters. Although only a few of these parameters are applicable to pul-
sars, many tests were quickly formulated to allow for the testing of gravity in the
strong-field regime about pulsars. One of these tests is the Shapiro time delay16
of light; an effect manifested through a retardation of an electromagnetic wave as
it moves through the warping of spacetime caused by mass and energy. Within
the PPN framework, the Shapiro delay is a function of the PPN parameter γ and
as such, measurements of the effect should allow for a tighter constraint to be
placed upon that parameter. Through simulations of the timing effects occurring
in binary systems, such as Roemer delay and Shapiro delay, and the categorising
of which systems prove promising in detecting these effects, it becomes possible
to flag the type of system that should be used for the measurement of the Shapiro
delay and the resultant constraint on γ.


    To what extent has γ been constrained and what hope is there for further
constraint? The PPN parameter γ is constrained via two effects that are inti-
mately linked as we shall see later. These effects are the Shapiro time delay and
the deflection of light. The deflection of light is an effect that was predicted by
Einstein in 1916 and first tested by Eddington17 , unfortunately to low (around
30%) accuracy. Very Long Baseline Interferometry (VLBI) measurements made
in 1995 of the sources 3C273 and 3C279 as they approached superior conjunction
with our sun led to a constraint of 1 (1 + γ) = 0.9996 ± 0.001718 while an anal-
                                     2
 10
     Nordtvedt, K. Jr.,“Post Newtonian Metric for a General Class of Scalar-Tensor Gravita-
tional Theories and Observational Consequences”,Ap. J., 161, 1059, (1970)
  11
     Thorne, K. S., Will, C.M., “Theoretical Framework for Testing Relativistic Gravity I:
Foundations”,Ap. J., 163, 595, (1971)
  12
     Will, C.M.,“Theoretical Framework for Testing Relativistic Gravity II: Parametrized Post-
Newtonian Hydrodynamics and the Nordvedt Effect”,Ap. J., 163, 611, (1971)
  13
     Will, C.M.,“Theoretical Framework for Testing Relativistic Gravity III: Conservation Laws,
Lorentz Invariance and Values of the PPN parameters”,Ap. J., 169, 125, (1971)
  14
     Will, C.M., Nordtvedt, K. Jr.,“Conservation Laws and Preferred Frames in Relativistic
Gravity I: Preferred-Frame Theories and an Extended PPN Formalism”,Ap. J., 177, 757,
(1972)
  15
     Will, C.M., Nordtvedt, K. Jr.,“Conservation Laws and Preferred Frames in Relativistic
Gravity II: Experimental Evidence to Rule Out Preferred-Frame Theories of Gravity”,Ap.
J., 177, 775, (1972)
  16
     Shapiro, I.I., “The Fourth Test of General Relativity”,Phys. Rev. Lett., 13, 798, (1964)
  17
     Eddington, A.S., “The Total Eclipse of 1919 May 29 and the Influence of Gravity on Light”,
The Observatory, 42, 119-122, (1919).
  18
     Lebach, D. E., Corey, B. E., Shapiro, I. I., et al., Measurement of the solar gravitational de-
flection of radio waves using very-long-baseline interferometry, Phys. Rev. Lett., 75, 14391442,
(1995).
                                                                                                   7


ysis made in 2004 of almost 2 million VLBI measurements made of 541 sources
using 87 VLBI sites improved this constraint to 0.99992 ± 0.0002319 , the tightest
constraint on this parameter using this technique.
    The Shapiro time delay of light is the second means of constraining γ. First
tested in 197720 with the Viking mission to Mars, it has been most accurately
constrained with the recent Cassini mission to Saturn21 where Doppler tracking
constrained γ − 1 to (2.1 ± 2.3) × 10−5 meaning that γ must be within 0.0012%
of unity, an even tighter constraint than light deflection tests can offer.
    Within the weak field regime it is hoped that NASA’s Laser Astrometric Test
Of Relativity (LATOR) mission will constrain γ to better than 1 part in 108 but
as launch is expected no sooner than 2009-2010, much work can still be done in
the interim.
    The testing of gravity in the strong field regime is still regarded as a fairly
fledgling science in that systems suitable for this testing are scarce; only 8 dou-
ble neutron star binaries are known22 and it is these systems that are required
for such investigations. As the population of compact binaries increases, and
existence/non-existence of a pulsar-black hole binary is determined, we can only
improve our knowledge of the universe and the gravity that holds it together.




 19
     Shapiro, S. S., Davis, J. L., Lebach, D. E., and Gregory, J. S., Measurement of the solar
gravitational deflection of radio waves using geodetic very-long-baseline interferometry data,
1979-1999, Phys. Rev. Lett., 92, 121101, (2004).
  20
     Shapiro, I.I., Reasenberg, R.D., et al., “The Viking Relativity Experiment”,J. Geophys.
Res., 82, 4329-4334, (1977).
  21
     Bertotti, B., Iess, L., and Tortora, P., A test of general relativity using radio links with the
Cassini spacecraft, Nature, 425, 374376, (2003).
  22
     Stairs, I.H., “Pulsars in Binary Systems: Probing Binary Stellar Evolution and General
Relativity”, Science, 304, 547, (2004).
                          2. THE PPN FORMALISM


                                  2.1 Background
The testing of General Relativity through the timing of radio pulsars falls into two
broad categories: setting limits on the magnitudes of parameters that describe vi-
olation of equivalence principles, and verifying that the measured post-Keplerian
timing parameters of a given binary system match the predictions of strong field
GR better than competing theories. This thesis will primarily focus on the latter
but it would be negligent not to mention the former and give some explanation
of what is meant in the following few pages.
    The Equivalence Principle basically states that there exists no experiment by
which one could tell the difference between a uniform gravitational field and a
uniform acceleration in the absence of a gravitational field and that these two
frames may thus be regarded as equivalent. Over time, it was realised that three
different forms of the Equivalence Principle needed to be distinguished. These
are the Weak Equivalence Principle (WEP), the Einstein Equivalence Principle
(EEP) and the Strong Equivalence Principle (SEP) and differ in the constraints
they impose on the system. The WEP states that objects of varying composi-
tions and masses will experience the same acceleration in an external gravitational
field. The EEP goes a little further and adds Lorentz invariance (non-existence of
preferred reference frames) and positional invariance (non-existence of preferred
locations) for non-gravitational experiments. This means that experiments will
have the same outcomes in inertial and freely falling reference frames. Lastly the
SEP adds Lorentz invariance and positional invariance for gravitational experi-
ments, thereby including experiments on objects with strong self-gravitation.
    For some years now, the WEP and EEP have been experimentally verified
and there is little doubt with regards to their accuracy in explaining observed
phenomena. The Strong Equivalence Principle, however, still needs to stand up
to experimental scrutiny. In order to more easily determine whether a violation of
the SEP had occurred, a formalism was developed to facilitate comparisons with
other gravitational theories in the strong gravitational limit. This parameterised
post-Newtonian (PPN) formalism, mentioned in the previous chapter, uses ten
parameters (γP P N , β, ξ, α1 , α2 , α3 , ζ1 , ζ2 , ζ3 , ζ4 ,) to compare competing theories
with one another. Pulsar timing allows for strict limits to be placed on α1 , α3 and
ζ2 as well as placing limits on other SEP violations that constrain combinations
of PPN parameters. Possible violations are the Nordtvedt effect, the existence of
2.2 The 10 PPN parameters and their origin                                            9


dipolar gravitational radiation and possible time variation of Newton’s constant.
The parameters and SEP violations will now be discussed.

                2.2 The 10 PPN parameters and their origin
The PPN metric is given by1 ,


      g00   = −1 + 2U + 2βU 2 − 2ξΦW + (2 + 2γ + α3 + ζ1 − 2ξ2 )Φ1
            + 2(3γ − 2β + 1 + ζ2 + ξ)Φ2 + 2(1 + ζ3 )Φ3 − 2(3γ + 3ζ4 − 2ξ)Φ4
            − (ζ1 − 2ξ)A − (α1 − α2 − α3 )ω 2 U − α2 ω i ω j Uij + (2α3 − α1 )ω i Vi
            + O( 3 ),                                                              (2.1)
                1                                    1
      g0i   = − (3 + 4γ + α1 − α2 + ζ1 − 2ξ2 )Vi − (1 + α2 − ζ1 + 2ξ)Wi
                2                                    2
              1             i        j          5/2
            −   (α1 − 2α2 )w U − α2 w Uij + O( ),                                  (2.2)
              2
      gij   = (1 + 2γU + O( 2 ))δij                                                (2.3)
                                                                                   (2.4)

where the various potential terms are given by
                                             ρ
                                   U=              d3 x ,                         (2.5)
                                          |x − x |
                                     ρ(x − x )i (x − x )j 3
                           Uij =                         d x,                     (2.6)
                                          |x − x |3
                        ρ ρ (x − x )   (x − x ) (x − x ) 3 3
                ΦW =              3
                                     ·           −          d xd x ,              (2.7)
                          |x − x |      |x − x |   |x − x |
                                      ρ [v · (x − x )]2 3
                              A=                       d x,                       (2.8)
                                          |x − x |3
                                            ρv2 3
                                   Φ1 =            d x,                           (2.9)
                                          |x − x |
                                            ρU
                                   Φ2 =            d3 x ,                        (2.10)
                                          |x − x |
                                            ρΠ
                                   Φ3 =            d3 x ,                        (2.11)
                                          |x − x |
                                             ρ
                                   Φ4 =            d3 x ,                        (2.12)
                                          |x − x |
  1
    Will, C.M., “The Confrontation Between General Relativity and Experiment”,Liv. Rev.
Rel., submitted (astro-ph:gr-qc/0510072),(2005)
2.2 The 10 PPN parameters and their origin                                           10

                                            ρ vi 3
                                 Vi =              d x,                          (2.13)
                                          |x − x |
                                  ρ [v · (x − x )](x − x )i 3
                        WI =                               d x,                   (2.14)
                                          |x − x |3
    It is within these equations that the 10 PPN parameters appear and below is
a brief description of each.
    The first parameter, γP P N , measures the space curvature produced by a unit
rest mass. In GR it takes the value 1 and is measured through observing light
deflection due to the presence of mass. In theories that predict some violation of
the SEP, known as semi-conservative theories, as well as some fully conservative
theories, that don’t violate the SEP, this parameter takes on a value other than
1 and so is a prime candidate for testing GR.
    The parameter β measures any non-linearity in the superposition law for
gravity and has been accurately measured through measurement of the precession
of Mercury’s perihelion. This solar system test was one of the first big successes
of general relativity and over the years has not been found to contradict GR in
any way. The parameter, usually measured as |β − 1| has been constrained to
3 × 10−3 , in agreement with the prediction of GR that β = 1.
    ξ is the parameter that measures any violation of positional invariance or
the existence of preferred-location effects. In GR the value is zero as GR does
not support the supposition that the outcome of local gravitational experiments
depends on the location of the laboratory with respect to the source of gravitation.
The basic result of such an effect is the inconstancy of the gravitational constant
as locally measured. The result of this, in terms of the parameter ξ, is the
existence of anomalous torques that cause a random alignment of the sun’s spin
axis with the ecliptic. ξ has been constrained to a few parts in 107 .
    The parameters α1 , α2 and α3 all measure various aspects of preferred-frame
effects. As with ξ, GR predicts that there is no effect of the velocity of the
laboratory with respect to the mean rest frame of the universe on the outcome of
local experiments and so these parameters should all be zero. An orbital evolution
of a binary system due to a relative velocity between it and the “universal”
background reference frame given by the CMB, will result in a non-zero α1 ,
however. Similarly a random alignment of the sun’s spin axis with the ecliptic,
as with a non-zero ξ, will result in a non-zero α2 . In the event of a self-acceleration
of a rotating body orthogonal to both its spin and its absolute velocity, the term
α3 would be non zero and violate both local Lorentz invariance and conservation
of momentum.
    The last four PPN parameters (ζ1 , ζ2 , ζ3 and ζ4 ) all indicate a non-
conservation of momentum. ζ2 along with α3 predicts an acceleration of the
center of mass of a binary system according to
                                          πm1 m2 (m1 − m2 )
                  acm = (α3 + ζ2 )                                enp .          (2.15)
                                     Pb [(m1 + m2 )a(1 + e2 )]3/2
2.2 The 10 PPN parameters and their origin                                            11




Fig. 2.1: A brief description of the 10 PPN parameters with their values in some of
          the competing theories (Will, 2005). It should be noted that there is no
          discrepancy in the values of the last 5 parameters. This was not always the
          case and certain theories predicted non-zero values for these parameters but
          over time it was realised that these theories were flawed and only the surviving
          theories are shown here.




Fig. 2.2: An overview of a few of the PPN parameters with their values in some of
          the more plausible competing theories (Will, 2005). Again it should be noted
          that there is no discrepancy in the values of the last 5 parameters.
2.3 The Nordtvedt effect                                                              12


In this equation np is a unit vector from the center of mass to the periastron of
body 1 with a mass m1 . m2 is the mass of the companion object, e the eccentricity
of the system, a the semi-major axis length and Pb is the binary orbital period.
                                                                             ˙
This acceleration leads to a Doppler type contribution to the pulsar’s P , which is
the first time derivative of its spin period, in that, in the case of higher eccentricity
orbits2 , the projection of acm along the line of sight to the earth will change with
time producing a second derivative to the pulse period given by
                                               2
               ¨ P                    2π           X(1 − X) eω cos ω
                                                               ˙
               P = (α3 + ζ2 )m2 sin i                     2 (1 − e2 )3/2
                                                                                 (2.16)
                  2                   Pb           (1 + X)
                                              ˙
where X = m1 /m2 , ω is the periastron angle, ω its rate of change and i is the
orbital inclination. The PPN parameters are summarised in Figures (2.1) and
(2.2).

                          2.3 The Nordtvedt effect
The SEP states that differing masses will experience the same acceleration in an
external gravitational field even when taking their self-gravitation into account.
If the SEP were violated, it would imply that the earth and the moon, containing
different fractional contributions from self-gravitation, would fall differently in
the sun’s gravitational field. This would result in a “polarisation” of the orbit
in the direction of the sun. This was first pointed out by Nordvedt3 when he
suggested testing the effect using laser ranging to the moon. The Nordtvedt
effect does not constrain any one parameter but rather a combination of seven of
them given by
                                     10         2    2    1
                  η = 4β − γ − 3 −      ξ − α1 + α2 − ζ1 − ζ2 .                  (2.17)
                                     3          3    3    2
                                                                   1
In General Relativity η = 0 while in tensor-scalar theory η = 2+ω + 4Λ. In
the case of η = 0 the earth will fall towards the sun with a slightly different
acceleration to that of the moon resulting in a perturbation of the earth-moon
orbit causing a polarisation of the orbit. This results in a perturbation of the
earth moon distance given by:
                          δr = 13.1η cos (ωm − ωs )t     [m],                    (2.18)
where ωm , ωs are the angular frequency of the moon and sun about the earth,
respectively. Lunar laser-ranging allows for measurement of the earth-moon dis-
tance to an accuracy of about 1 cm. This in turn allows for a limit to be placed
of η < 0.001.
  2
    Will, C.M., “Is momentum conserved?                A test in the binary system
PSR 1913+16”,Astrophys. J., 393. L59-L61. (1992)
  3
    Nordtvedt, K., “Testing Relativity with Laser Ranging to the Moon”, Phys. Rev., 170.
1186-1187. (1968)
2.3 The Nordtvedt effect                                                                  13


   In the strong-field formalism, however, η is not used directly but rather as
part of a larger function ∆i , that contains higher order η terms and which relates
the gravitational and inertial masses for various objects, i, through
                  mgrav
                                 = 1 + ∆i                                             (2.19)
                 minertial   i
                                                                         2
                                              E grav            E grav
                                 = 1+η                     +η                + ....   (2.20)
                                              mc2      i        mc2      i

where η’ is the first order correction to η. This equation perhaps highlights the
difference between strong-field and weak-field gravity better than any other. In
the case of the earth, the quantity E grav /mc2 takes on a value of the order of 10−10
and even for the entire solar system the quantity takes on a value no bigger than
10−5 . For neutron stars and black-holes, however, this quantity takes on values
of 0.2 and 0.5, respectively, and therefore the higher order terms in equation
(2.20) cannot be ignored. It is this factor that distinguishes weak from strong
field gravity.
    In the event of a SEP violation, the two components of a binary system will fall
differently in the gravitational field of the galaxy, denoted g, and the equations of
motion will contain an extra acceleration ∆net g, where ∆net = ∆pulsar −∆companion
in the case of a binary pulsar system. This ∆net g term will affect the evolution
of the system, in particular the evolution of the eccentricity. This implies that
there will be competition between the general relativistic evolution of eccentricity
through advance of periastron and this new forcing of the eccentricity towards
alignment with the projection of g, g⊥ , onto the orbital plane, and this term will
be constant. The time-dependent eccentricity vector can thus be written as

                                    e(t) = eF + eR (t).                               (2.21)

Here eR (t) is the relativistic evolution of the eccentricity vector due to advance
                ˙
of periastron (ω) and eF is the forced component. The equation for this forced
component for low eccentricity is given by4

                                           1 ∆net g⊥ c2
                                 |eF | =                    ,                         (2.22)
                                           2 F GM (2π/Pb )2

where M is the total mass of the system, G is Newton’s constant, Pb is the
orbital period and F = 1 in the case of GR. There are a few considerations to be
taken into account when selecting binary pulsar systems for useful testing of this
                                             ˙
effect. The rate of advance of periastron (ωR ) must be greater than the rate of
galactic rotation in order for g to be constant. Furthermore, the pulsar must be
                          ˙
sufficiently old that the ω-induced rotation of the eccentricity e has completed
  4
    Wex, N., “New limits on the violation of the Strong Equivalence Principle in strong field
regimes”. Astron. Astrophys., 317, 976-980. (1997).
2.4 Variation of Newton’s constant                                                     14




Fig. 2.3: Illustration of the polarisation of a low eccentricity orbit due to the presence
          of a forcing vector g. The diagram shows the forced eccentricity eF and
          the eccentricity evolving under general relativistic advance of periastron eR
          through an angle θ


enough turns that the component eR (t) can be assumed to be randomly oriented
in equation (2.21). One can only assume this randomness has been achieved in
pulsars with a characteristic age of τc     ˙
                                        2π/ω. For this reason young pulsars can
be ruled out as valid systems in which to perform this test.


                    2.4 Variation of Newton’s constant
As mentioned in the discussion of the PPN parameters themselves, some theo-
ries that allow violations of the SEP also allow for a time varying gravitational
constant, G. Although this variation is expected to occur on a Hubble timescale,
 ˙
G/G ∼ H◦ ∼ 0.7 × 10−10 yr−1 , it should be testable with the aid of neutron stars,
white dwarfs and binary systems.

                                      2.4.1 Spin tests
It is fairly obvious that any variation in the gravitational constant would lead to
a variation in the binding strength between atoms in a neutron star and thereby
affect its moment of inertia. This would lead to a variation in the spin period on
the same timescale as the variation in G. This change has been written as5

                                 P˙            ∂ ln I       G˙
                                           =                                       (2.23)
                                 P     ˙
                                       G
                                               ∂ ln G   N
                                                            G,
  5
    Goldman, I., “Upper Limit on G Variability Derived from the Spin-Down of PSR
0655+64”,Mon. Not. R. Astron. Soc., 244,184-187, (1990).
2.4 Variation of Newton’s constant                                                     15


for I the moment of inertia at a constant number of baryons N . The neutron star
equation of state will give the baryon density for the star and, since there is no
agreement as to what the correct equation of state is, there is still much debate
                                                         ˙
over what the baryon density and therefore the ratio G/G is. The value has,
however, been calculated for the pulsar PSR B0655 + 64 by Goldman (1990) and
               ˙
found to be |G/G| ≤ (2.2 − 5.5) × 10−11 yr−1 . This indicates that only pulsars
with fairly large characteristic ages, of the order of tens of Gyr, can be used
to perform this test as in younger pulsars the effect would simply not have had
enough time to become manifest.

                              2.4.2 Orbital decay tests
Since the orbital period of a binary system is dependent on G, it is obvious that
a variation therein would lead to an orbital variation with time. This problem
was first considered in 1988 by Damour, Gibbons and Taylor6 who found that
                                      ˙
                                     Pb           ˙
                                                  G
                                              = −2 .                               (2.24)
                                     Pb   ˙
                                          G
                                                  G

It was realised, however7 , that a changing G would also affect the internal struc-
ture of the neutron star and thus the total mass and angular momentum. Taking
these factors into consideration it was shown that the change should be given by
            ˙
           Pb                 m 1 c1 + m 2 c2   3 m1 c2 + m2 c1        ˙
                                                                       G
                    =− 2−                     −                                    (2.25)
           Pb   ˙
                G
                               m1 + m2          2  m1 + m2             G

where ci is the compactness of body “i”. Compactness is defined by Nordtvedt
(1990) to be
                                       G δmi
                                  ci =         ,                         (2.26)
                                       mi δG
where δmi is the variation in mass of body “i” and δG is the variation in grav-
itational constant. Equation (2.25) shows that with improved timing and mea-
              ˙
surement of Pb , and with tighter constraints on the component masses, tighter
constraints can be placed on a variation of G.

                     2.4.3 Changes in the Chandrasekhar mass
During the formation of white dwarf stars, complete core collapse is prevented
through electron degeneracy pressure. The mass at which this electron degener-
  6
     Damour, T., Gibbons, G.W., Taylor, J.H., “Limits on the Variability of G Using Binary
Pulsar Data”, Phys. Rev. Lett., 61, 1151-1154, (1988)
   7                   ˙
     Nordtvedt, K., “G/G and a Cosmological Acceleration of Gravitationally Compact Bod-
ies”,Phys. Rev. Lett., 65, 953-956, (1990).
2.4 Variation of Newton’s constant                                                    16


acy pressure can no longer oppose gravitational collapse is known as the Chan-
drasekhar mass, MCh , and is expressed as
                                                3/2
                                           ¯
                                           hc         1
                                 MCh ∼                   .                        (2.27)
                                           G          m2
                                                       n

      ¯
Here h is Planck’s constant and mn is the neutron mass. Clearly a change in the
value of G would affect the Chandrasekhar mass and so by comparing the masses
of very old to very young pulsars, it should be possible to identify any change
in mass over time taking into account the obvious mass-loss that older pulsars
would have experienced over that time. Such an analysis was done by Thorsett8
                                             ˙
who found a limit in the change of G to be G/G = (−0.6 ± 4.2)−12 yr−1 at a 95%
confidence level, the strongest limit found thus-far.


    In summary, it has been shown that the PPN formalism has been used for
many years to distinguish between competing theories of gravity, and in some
cases even to lead physicists to discard a theory for not meeting experimental
criteria. Although constraining a parameter perfectly is never a scientifically
possible outcome, the placing of tighter limits through both strong and weak
field tests of gravity allows for the elimination of certain aspects of gravitational
theories, bringing us closer to an understanding of this most fundamental of
forces. Pulsars stand in the forefront of this testing process by allowing for the
direct testing of three of the parameters and three further combinations thereof.
Only through improved and more extensive timing and analysis of these most
dynamic of objects, will the scientific community have sufficient data available
to them to put gravitational theories under the microscope and finally put the
theory first elucidated by Sir Isaac Newton in his Philosophæ Naturalis Principia
Mathematica in 1686 to rest.




  8
   Thorsett, S.E., “The Gravitational Constant, the Chandrasekhar Limit, and Neutron Star
Masses”,Phys. Rev. Lett., 77, 1432-1435, (1996)
                       3. ORBITAL MECHANICS


In a study of the effects of General Relativity on signals emitted by a pulsar in
a binary system, it is important to understand the mechanics of binary orbits.
In order to do this efficiently, it is necessary to develop the theory from simple
elliptical orbits to that of slightly more complex binaries. In the next sections we
will evolve this theory, along with a fairly detailed discussion of how the binary
orbits under discussion actually formed.

                           3.1 Orbital geometry
The analysis of the geometry describing orbits leads to the understanding of
the parameters that are used in calculating the many effects occurring in such
systems. It has long been known that all orbits can be described via the family of
conic sections known as ellipses. The analysis of ellipses illustrates relationships
between parameters leading to equations for position and velocity as functions of
time. An analysis of Keplerian orbits now follows;

                            3.1.1 Orbital parameters
There are seven parameters needed to describe a satellite in orbit. These seven
parameters, known as the Keplerian parameters, define the ellipse of the orbit, the
position of the ellipse about the parent body and the position of the orbiting body
within the ellipse. In classical Keplerian mechanics, the ellipse has constant shape
and orientation and for this analysis we will use this approximation, generalising
the concepts in later sections. The Keplerian parameters are: the Epoch, the
orbital inclination, the right ascension of the ascending node, the argument of
perigee, the eccentricity, the mean motion and lastly the mean anomaly.
  1. The Epoch is simply a number that specifies the time at which all other
     parameters were measured. It is usually measured in Mean Julian Days
     and is often corrected to the solar system barycenter in order to have a
     standard time measurement on the earth’s globe.

  2. The orbital inclination is the angle between the plane of the orbit and the
     plane of the sky, that is the perpendicular plane to the line of sight from
     earth. An inclination angle of 0◦ is a face-on orbit while an inclination angle
     of 90◦ represents an edge-on orbit.
3.1 Orbital geometry                                                             18




       Fig. 3.1: Basic vectors defining the parameters for a Keplerian orbit


  3. The right ascension (RA) of the ascending node is the celestial equivalent
     of terrestrial longitude. The RA is measured as an increasing angle in the
     east-west direction, along the equator from the zero point at the vernal
     equinox of the component of interest’s orbit. The ascending node is the
     point at which the orbiting body crosses through the ecliptic of the binary
     system moving from south to north where north is defined as the direction
     to the perigee.

  4. The argument of perigee is the angle that determines the rotation of the
     ellipse within its orbital plane. It is the angle between the line of nodes and
     the semi-major axis measured at the center of the ellipse.

  5. Eccentricity is a measure of how distorted from a circle the orbit is. An
     eccentricity of 0 is a circular orbit while an eccentricity of close to 1 is a
     long and narrow orbit. Eccentricities of 1 are physically impossible as that
     would imply a linear oscillatory motion from perigee to apogee.

  6. The mean motion, n, is a concept used to indicate the approximate size
     of the orbit. It is the inverse of the orbital period and so is measured in
     number of revolutions per unit time. It is also related to the semi-major
     axis a via
                                             GM
                                      n=                                   (3.1)
                                              a3
     Where G is the gravitational constant and M the mass of the parent body.
3.1 Orbital geometry                                                           19


  7. Lastly the mean anomaly is used to describe the position of the orbiting
     body at any instant in time. It is an angle that traces out the 360◦ of
     an orbit in equal time increments and is measured from the semi-major
     axis at the center of the ellipse. For circular orbits (eccentricity = 0) the
     line to which the angle is measured would point directly to the orbiting
     body. For eccentric orbits, however, the velocity of the orbiting body is
     a function of position within the orbit and so the line to which the mean
     anomaly is measured is not commensurate with the position of the body.
     Consider figure (3.1). For a system with primary focus f and origin o, the
     mean anomaly M is defined such that the area of triangle boz is always
     equal to the area of triangle af z. In this scenario a is the point on the
     auxiliary circle created by extending the perpendicular to the semi-major
     axis through the position of the orbiting body, p, and b is the point on the
     auxiliary circle creating the line which defines the angle that is the mean
     anomaly. Thus only at perigee and apogee does the line to which the mean
     anomaly is measured, point towards the orbiting body. At this point it is
     also perhaps useful to indicate the other two anomalies often referred to
     when dealing with orbits. The eccentric anomaly, E, is measured from the
     semi-major axis to the line ao, is a measure of the eccentricity of the orbit
     and is related to the mean anomaly via M = E − e sin E. Secondly the true
     anomaly, T , is measured at the focus of the ellipse and is the angle between
     the semi-major axis and the line from the orbiting body to the prime focus
     of the ellipse pf .



Continuing with our geometrical analysis: For an ellipse with semi-major axis a,
eccentricity e and true anomaly θ such as the one illustrated in figure (3.2) it is
always true that




                 Fig. 3.2: Diagram illustrating ellipse parameters
3.1 Orbital geometry                                                          20



                                  r + r = constant                          (3.2)
but when θ = θ = 0 as at perigee we can write the distances r and r as
                                 r = 2ae + (a − ae)                         (3.3)
and
                                     r = a − ae,                            (3.4)
and so
                                     r + r = 2a.                            (3.5)
But, from geometrical considerations
                               r cos θ = r cos θ − 2ae                      (3.6)
and
                                  r sin θ = r sin θ .                       (3.7)
Squaring these equations we have
                     r 2 cos2 θ = r2 cos2 θ + 4aer cos θ + 4a2 e2           (3.8)
and
                                 r 2 sin2 θ = r2 sin2 θ                     (3.9)
which then gives the relationship
                   r 2 = r2 cos2 θ + 4aer cos θ + 4a2 e2 + r2 sin2 θ.      (3.10)
Using equation (3.5) then yields
                        (2a − r)2 = r2 + 4aer cos θ + 4a2 e2 ,             (3.11)
such that
                             r(1 + e cos θ) = a(1 − e2 ),                  (3.12)
or in final form
                                      a(1 − e2 )
                             r=                       ;                   (3.13)
                                  1 + e cos (θ + θ◦ )
where the θ◦ takes into account any possible rotation of the ellipse about the
origin. This equation then relates the position of the orbiting body to the fun-
damental parameters explained above. Furthermore, this equation is used for
the determination of the Roemer delay in a binary orbit as discussed in Chapter
1. With some knowledge of the orientation of the orbit in space and the above
equation, the variation of distance through which the pulsar signal propagates
can be determined. Using r and the inclination angle of the orbit, i, one finds
that the distance to the pulsar as it moves in it’s orbit is given by
         R = [(r cos(π/2 − i) sin(θ))2 + (r cos(π/2 − i) cos(θ))2 ]1/2 .   (3.14)
3.2 Keplerian mechanics                                                         21


                           3.2 Keplerian mechanics
In a system with a mass m orbiting a mass M there are no preferred directions
or points other than the vector r going from mass M to mass m and possibly
the vector v which is the rate of change of r with time. Therefore it is irrelevant
how the basis vectors are chosen and as such we may as well choose them at M
      ˆ                       ˆ
with r pointing outward and θ in the direction of increasing θ, as shown in figure
(3.3).




                   Fig. 3.3: Basis vectors for an orbiting system

     Newton II states
                                      ¯    a
                                      F = m¯                                (3.15)
and
                                    ¯    GM m ˆ
                                    F = − 2 R,                              (3.16)
                                          r
          ¯    r
and since r = rˆ we have
                                     r
                                    d¯   dr        r
                                                  dˆ
                               ¯
                               v=      =    ˆ
                                            r+r                             (3.17)
                                    dt   dt       dt
                                         ˙r      ˆ
                                               ˙θ;
                                       = rˆ + rθ                            (3.18)
so
                         v
                        d¯   dr˙       dˆ dr ˙ ˆ
                                        r              ˙
                                                     dθ ˆ        ˆ
                 a=
                 ¯         =     r + r + θθ + r θ + rθ
                                 ˆ ˙                          ˙ dθ          (3.19)
                        dt   dt        dt   dt        dt        dt
                           = rr + rθ
                             ¨ˆ ˙     ˆ ˙ ˙ˆ     ¨ˆ      ˙ ˙r
                                     ˙θ + rθθ + rθθ + rθ(θˆ)                (3.20)
                                    ˙ r       ¨     ˙ ˆ
                           = (¨ − rθ2 )ˆ + (rθ + 2rθ)θ.
                              r                                             (3.21)
And so it is seen that
                           ¯ ¨ˆ      ¨ˆ
                           a = rr + rθθ +    extra terms.                   (3.22)
It is these extra terms that indicate how the basis is changing as the object moves
through the space and, in combination, are often referred to as a centrifugal force
despite being only an ”imaginary” force. One therefore finds
                        GM m            ˙ r        ¨
                    −                                   ˙˙ ˆ
                             r = m(¨ − rθ2 )ˆ + m(rθ + 2rθ)θ,
                             ˆ     r                                        (3.23)
                         r2
3.2 Keplerian mechanics                                                       22


thus
                                     GM         ˙
                                 −       = r − rθ2
                                           ¨                               (3.24)
                                      r2
and
                                          ¨    ˙˙
                                     0 = rθ + 2rθ.                         (3.25)
The second of these ordinary differential equations can be multiplied out by r to
show
                                  ˙˙     ¨
                           0 = 2rrθ + r2 θ                                 (3.26)
                                 d 2 ˙         d ˙
                             =      (r ) θ + r2 θ                          (3.27)
                                 dt            dt
                                d 2˙
                             =    (r θ)                                    (3.28)
                               dt
and so
                                  ˙
                               r2 θ = k, a constant                        (3.29)
which implies
                                       k 2 GM
                               0=r−
                                 ¨         + 2 .                          (3.30)
                                        r3    r
    Solving these differential equations to find r = r(t) and θ = θ(t) is done by
letting r = u−n . In this case
                                 d −n
                            r=
                            ˙       u = −nu−n−1 u;
                                                ˙                          (3.31)
                                 dt
but
                                       du dθ  ˙ du
                               ˙
                               u=            =θ ,                          (3.32)
                                       dθ dt    dθ
and
                                 ˙   k
                                 θ = 2 = ku2n ,                            (3.33)
                                    r
therefore
                                                du
                                     u = ku2n
                                     ˙             .                       (3.34)
                                                dθ
Equation (3.31) can then be written as
                                                        du
                            r = −nu−n−1 ku2n
                            ˙                                              (3.35)
                                                        dθ
                                                 du
                                 = −nun−1 k         ,                      (3.36)
                                                 dθ
such that the second derivative is
                      d          du       dθ
                r =
                ¨        −nkun−1                                           (3.37)
                      dθ         dθ       dt
                                                                
                                                 2
                                           du                d2 u
                 = −nk(n − 1)un−2                   − nkun−1 2  ku2n .   (3.38)
                                           dθ                dθ
3.2 Keplerian mechanics                                                                     23


                                                       ¨
But equation (3.30) already gave us a relationship for r that can be rewritten as

                                   r = k 2 u3n − Gmu2n
                                   ¨                                                     (3.39)

from which we finally get
                                                             2
                                                        du                      d2 u
         k 2 u3n − Gmu2n = −nk 2 (n − 1)u3n−2                    − nk 2 u3n−1        .   (3.40)
                                                        dθ                      dθ2

Choosing n = 0, 1 the du term is eliminated thus simplifying our equation but
                         dθ
the n = 0 solution is trivial, therefore we choose n to be 1. In this case
                                         d2 u
                        0 = −k 2 u2           − k 2 u3 + GM u2 , or                      (3.41)
                                         dθ2
                                     d2 u       GM
                                  0=    2
                                          +u− 2 ;                                        (3.42)
                                     dθ          k
         GM
   but    k
              is a constant and this allows us to write

                              d2          GM                      GM
                         0 =     2
                                     u− 2                +u−                             (3.43)
                              dθ           k                       k2
                                2
                              d
                            =      (s) + s.                                              (3.44)
                              dθ2
An equation of this form has the solution

                                    s = A cos (θ + θ0 )                                  (3.45)

as it is analogous to the analysis of the harmonic oscillator. We thus have
                                                         GM     1
                           u = A cos (θ + θ0 ) +            2
                                                              =                          (3.46)
                                                          k     r
and therefore
                                              1
                        r =       GM                                                     (3.47)
                                   k2
                                         + A cos (θ + θ0 )
                                     2
                                   k                    1
                           =                     Ak2
                                                                     .                   (3.48)
                                  GM        1+   GM
                                                       cos (θ + θ0 )
This is analogous to equation (3.13) and so we must have that
                           Ak 2                                      k2
                      e=                  and      a(1 − e2 ) =         .                (3.49)
                           GM                                       GM
   At perigee/apogee the separation of the two bodies will be a min/max,
                                                                     ˆ
rmin /rmax but more importantly the velocity will be entirely in the θ direction,

                                            ¯    ˙ˆ
                                            v = rθθ.                                     (3.50)
3.2 Keplerian mechanics                                                         24


This implies
                           ¯     1       1  ˙
                           Ek = m¯2 = m(rθ)2 .
                                   v                                        (3.51)
                                 2       2
But from equation (3.29) we can deduce that
                                       ˙
                                   (rm θ)2   k2
                                           = 2,                             (3.52)
                                     2      2rm
the kinetic energy per unit mass. The total energy per unit mass will then be
the combination of kinetic and potential energies
                                        k2   GM
                                   =     2
                                           −    ,                           (3.53)
                                       2rm   rm
and therefore
                               2     GM          k
                              rm +        rm −     =0                       (3.54)
                                                 2
which obviously leads to
                                                       −1
                              MG            M 2 G2 2 
                        rm =  2 ±                + 2         .             (3.55)
                              k              k4    k

Here we have used equation (3.53). Upon comparison with equation (3.13) and
(3.49) and remembering that at perigee/apogee cos (θ + θ0 ) = ±1, we see that

                                         M 2 G2 2
                               A=              + 2.                         (3.56)
                                          k4    k
Now by combining equations (3.13) and (3.49) one gets
                                            k2
                                  rm =              ,                       (3.57)
                                         M G(1 + e)
such that
                            M 2 G2 (1 + e)2 M 2 G2 (1 + e)
                          =                 −                               (3.58)
                                 2k 2            k2
                              2 2
                            M G (1 + e)[1 + e − 2]
                          =                                                 (3.59)
                                       2k 2
                              2 2 2
                            M G (e − 1)
                          =                                                 (3.60)
                                 2k 2
                              MG
                          = −      ,                                        (3.61)
                              2a
                    2
                    k
since a(1 − e2 ) = 2M . From this it is then easily seen that the total energy per
unit mass is given by
                                MG       v2 M G
                              −        =    −       ,                        (3.62)
                                 2a       2     r
3.2 Keplerian mechanics                                                    25


and upon rewriting we obtain the orbital speed for a mass m orbiting a mass M
                                        2 1
                             v2 = M G    −  .                          (3.63)
                                        r a
3.3 Binary pulsar evolution                                                            26


                         3.3 Binary pulsar evolution
It is fairly obvious that a binary pulsar is a system of two, mutually orbiting
stellar bodies, one of which is a detectable pulsar. In the case of strong field
General Relativistic tests, however, what is meant by a binary pulsar must be
stated more precisely. It means a binary system wherein there exists a pulsar
with a compact binary companion such as a white dwarf star, another neutron
star or in the most extreme case, a black hole. The reason that an extrasolar
planet or main-sequence star are excluded is that the gravitational field generated
by such bodies is not sufficient to affect the signal from the pulsar sufficiently to
allow tests of General Relativity. Stellar evolutionary models have indicated that
only one in a hundred stars will evolve to become a neutron star or black hole
and furthermore have shown that around fifty percent of stars occur in binary
systems. This implies that only one out of every two hundred stellar systems will
form a compact binary of some sort. Of course the chances that one member of
the binary system is a pulsar beaming towards earth is again many times smaller
and therefore the short list of a few hundred detected compact binary systems is
understandable. The list of only 8 systems where both components are neutron
stars is thus also understandable in light of the rarity of such systems.

                            3.3.1 Evolutionary processes
There are believed to be four evolutionary scenarios1 that can give rise to a com-
pact binary system, with the development of a neutron star-black hole binary
being ignored due to the non-detection of such a system at the time of writing
of this thesis. The first scenario results in the formation of a long-period binary
system with a millisecond pulsar and a low-mass white dwarf companion. One
starts with a main-sequence (> 8M ) star and a low mass (∼ 1M ) companion.
The primary star undergoes standard stellar evolution and explodes as a super-
nova with the core collapsing to form a neutron star. The lower mass companion
undergoes its usual, much slower evolution and, during its expansion phase, over-
flows its Roche-lobe allowing accretion of the matter onto the companion neutron
star. This is a very slow steady process and results in the system being visible in
the x-ray regime due to production of hard and soft x-rays when the incident ac-
creted matter is pulled onto the pulsar along the magnetic field lines and strikes
the pulsar magnetosphere. The final result is a low-mass white dwarf in slow
orbit about a recycled millisecond pulsar.
    In the case of a (> 8M ) main-sequence star with an intermediate-mass
(∼ 5M ) companion, one again has neutron star formation through a super-
nova explosion of the primary star. However, due to the much larger size of the
companion, one does not have an accretion driven evolution but rather common-
  1
    Stairs, I.H., “Pulsars in Binary Systems: Probing Binary Stellar Evolution and General
Relativity”, Science, 304, 547, (2004).
3.3 Binary pulsar evolution                                                     27


envelope evolution where the neutron star spirals into and eventually expels the
envelope of the companion. This results in a mildly-recycled pulsar in a close
orbit with a massive (∼ M ) white dwarf.
    The third scenario is where the evolutionary process begins to get very in-
teresting. In the event that neither star is massive (> 8M ), but rather there
exists one star at (∼ 7M ) while the other is slightly smaller at (∼ 5M ), mass
transfer will occur from the primary star to its lighter companion due to the
primary’s faster evolution. Through this process the primary loses enough mass
to evolve normally into a white dwarf star. The companion is now much larger
than it initially was and so with its standard evolution it expands to envelop the
white dwarf primary until the white dwarf spirals in and expels the envelope.
The result of this is a white dwarf - helium star binary that, upon supernova
explosion of the helium star, forms a young pulsar in orbit about a relatively
massive white dwarf companion.
    The last possibility for this sort of binary formation is certainly the rarest.
For two main-sequence stars, both of which have mass (> 8M ), the primary
explodes as a supernova to form a neutron star. There then comes a period of
mass transfer from the companion to the neutron star in the form of a companion
wind. Again common-envelope evolution occurs with the neutron star expelling
the envelope of the secondary to leave a neutron star - helium star binary system.
At this stage Roche-lobe overflow can occur from the helium star to the neutron
star, resulting in a high-mass x-ray binary, but the eventual supernova explosion
of the helium star to form a second neutron star is what makes this such an
interesting system. The first such system, PSR 1913+16, was detected by Russell
Hulse and Joseph Taylor (1975) and their study of this system resulted in their
being awarded the 1993 Nobel prize in physics. The first double pulsar system,
PSR J0737-3039, in which both neutron stars are beaming towards earth was
only discovered in 2004 by Lyne, Burgay, Kramer et al. (2004).
    A neutron star is formed during the final stages of stellar evolution when
a star has converted almost all the hydrogen in its core to heavier elements.
When electron degeneracy pressure can no longer support the the star, the inner
layers collapse under the influence of gravity while the outer layers are blasted
off into space in a supernova explosion. Under this free-fall acceleration, the
matter in the core of the star collapses into itself reducing its volume by up to
109 times. At these volumes with the large mass still present in the core of the
star, the density exceeds nuclear density and electrons are forced into the nuclei
to form neutrons. This then allows for neutron degeneracy pressure to sustain
the star against further collapse. In this process a star is formed with a density
in excess of 1013 g/cm3 and magnetic fields of the order of 1010 Gauss or more. As
a result, any free charged particles near the surface of the star would be seized
by the magnetic field and accelerated along the field lines. By the process of
bremsstrahlung, these particles would radiate electromagnetic radiation and this
energy could be detected over many parsecs of space.
      4. TOPICS FROM THE THEORY OF CURVED SPACE


In any attempt to test a theory, one must have a firm enough grasp of the subject
matter so as to understand any behaviour that deviates from prediction. An
overview of some of the basic principles of General Relativity is mandatory and
will help in the full understanding of the concepts that will be introduced later.
    The understanding of General Relativity requires the ability to switch coor-
dinate systems at will, knowing that events can be described within an infinity
of different coordinate systems.

                       4.1 Curvilinear coordinates
Consider figure (4.1). Construct a rectilinear, or cartesian, 2-dimensional coor-
dinate system according to X i = X i (x) with basis vectors Ei . Pick any point
within that space such that it has coordinates x1 and x2 . The x1 grid line go-
                                                     o   o
ing through that point will then be given, in terms of the rectilinear coordinate
system, by
                                X 1 = X 1 (x1 , x2 )
                                                 0                          (4.1)
and
                                X 2 = X 2 (x1 , x2 ).
                                                 0                           (4.2)
Just as the components of the x2 grid line will be

                                X 1 = X 1 (x1 , x2 )
                                            0                                (4.3)

and
                                X 2 = X 2 (x1 , x2 ).
                                            0                                (4.4)
In this case the components of the basis e1 in the (X 1 , X 2 ) coordinate system
will be given by
                                   ∂X 1 1 2
                                       (x , x )                              (4.5)
                                    ∂x1 0 0
and
                                   ∂X 2 1 2
                                       (x , x );                             (4.6)
                                   ∂x1 0 0
and therefore the basis vectors are
                                          ∂X i
                                   e1 =        Ei                            (4.7)
                                          ∂x1
4.1 Curvilinear coordinates                                                         29




      Fig. 4.1: Curvilinear coordinates embedded in a cartesian coordinate system


and
                                             ∂X i
                                    e2 =          Ei ;                          (4.8)
                                             ∂x2
or in index notation
                                        ∂X i
                                    eα =     Ei .                         (4.9)
                                        ∂xα
So within the coordinate system X i = X i (x) we have a new basis given by
equation (4.9). Both sets of bases must be expressible in terms of one another,
else they would not be basis vectors, and therefore the relation
                                             ∂xα
                                     Ei =         eα                           (4.10)
                                             ∂X i
must hold. Taking the derivative of eα will determine how this new basis changes
as we move through the curvilinear coordinate system. This is done as follows:
                         ∂eα      ∂2X i          ∂X i ∂Ei
                               =         (x)Ei +                               (4.11)
                         ∂xβ     ∂xβ ∂xα         ∂xα ∂xβ
                                   2 i
                                  ∂ X
                               =         (x)Ei                                 (4.12)
                                 ∂xβ ∂xα
                                       ∂Ei
where the second term is zero since    ∂xβ
                                             = 0 because the Euclidian basis does not
change. Thus
                                  ∂X i
                   eα (x + δx) =       (x + δx)Ei                              (4.13)
                                  ∂xα
                                   ∂X i        ∂2X i
                                =       (x) + β α (x)δxβ Ei                    (4.14)
                                   ∂xα        ∂x ∂x
4.2 The metric                                                                   30

                                     ∂2X i
                              = eα +   β ∂xα
                                             (x)δxβ Ei                       (4.15)
                                    ∂x
                                     ∂2X i          ∂xσ
                              = eα + β α (x)δxβ          eσ .                (4.16)
                                    ∂x ∂x           ∂X i
We now define
                               ∂ 2 X i ∂xσ
                                           ≡ Γσ ;
                                              αβ                             (4.17)
                              ∂xβ ∂xα ∂X i
and as such equation (4.16) becomes

                         eα (x + δx) − eα (x) = Γσ δxβ eσ .
                                                 αβ                          (4.18)

If we then divide through by δxβ and take the limit as it goes to zero, we get
                                   ∂eα
                                       = Γσ eσ .
                                          αβ                                 (4.19)
                                   ∂xβ
This equation tells us the way in which the bases, in the context of flat spaces,
change while one moves through a curvilinear coordinate system. By determining
the Γ’s for a given coordinate system, the basis can be found for every point within
that coordinate system. This is an essential piece of information needed for any
further calculations within a given spacetime.


                               4.2 The metric
How does one measure distance on a 2-dimensional surface embedded within
a 3-dimensional Euclidian space? Well, consider figure (4.2). In this system




Fig. 4.2: Diagram highlighting measurement on a 2-dimensional plane embedded in
          3-dimensions.

the question becomes, how do we measure the distance A→ B? Ordinarily on
4.3 Geodesics                                                                                            31


a flat space, such as can be approximated over an infinitesimal region of the
curved space, we would sum over all the infinitesimal lengths δs by using the
pythagorean relationship δs2 = δu2 + δv 2 but on curved surfaces this relationship
does not hold. The relationship δs2 = δx2 + δy 2 + δz 2 does, however. But since
x = x(u, v) we can write
                                   ∂x        ∂x
                              δx =    δu +      δv,                          (4.20)
                                   ∂u        ∂v
and squaring this yields
                                        2                     2
                           2       ∂x        2   ∂x                          ∂x ∂x
                   δx =                     δu +                  δv 2 + 2         δuδv.              (4.21)
                                   ∂u            ∂v                          ∂u ∂v

The line element δs2 can then be written as
                   2                    2                     2                   2             2
      2       ∂x            ∂y 2              2  ∂z                   2∂x              2   ∂y
 δs       =            δu +                 δu +                  δu +                δv +          δv 2 +
              ∂u            ∂u                   ∂u                    ∂v                  ∂v
                   2
              ∂z                    ∂x ∂x ∂y ∂y ∂z ∂z
                       δv 2 + 2          +     +      δuδv;                                           (4.22)
              ∂v                    ∂u ∂v ∂u ∂v ∂u ∂v
which is written more concisely in index notation as
                                                 ∂xi ∂xj α β
                                    δs2 = δij              δu δu                                      (4.23)
                                                 ∂uα ∂uβ
                                            = gαβ δuα δuβ .                                           (4.24)
                ∂x ∂x  i       j
Here gαβ = δij ∂uα ∂uβ and the δij term is the Euclidian metric in a Cartesian
system of coordinates. It is this gαβ term that is defined as the metric for a
particular curved space.

                                            4.3 Geodesics
A straight line may be described as the shortest distance between two points.
Light paths are also often used to describe straight lines because they naturally
follow the straightest route between two points. It is this property of light that is
of fundamental importance to this thesis. Mathematically, how is this indicated?
Consider 2 points connected by a curve C with the property (by Hamilton’s
variational principle) that
                                            S=             ˙
                                                      L(x, x, λ)dλ                                    (4.25)
                                                  C

where the curve is parameterised by xi = xi (λ). In this case the curve C for
which S is an extremum is the solution to the Euler-Lagrange equations
                                               d       ∂L             ∂L
                                        0=                        −                                   (4.26)
                                              dλ       ∂ xi
                                                         ˙            ∂xi
4.3 Geodesics                                                                        32


Now consider curves of shortest length in the surface xi = xi (u1 , u2 ).
                                                   2
                                    S=                 ds,                   (4.27)
                                               1

but from equation (4.24) we know that ds2 = gαβ δuα δuβ and so
                                         2
                               S=             gαβ δuα δuβ .                  (4.28)
                                     1

But the u’s are also parameterised by λ and as such
                                             duα
                                 duα =           (λ)dλ,                      (4.29)
                                             dλ
and hence we can write
                                                   ∂uα       ∂uβ
                   gαβ δuα δuβ = gαβ (u(λ))            (λ)dλ     (λ)dλ.      (4.30)
                                                   ∂λ        ∂λ
Finally we can rewrite equation (4.28) as
                                         2
                               S=            gαβ uα uβ dλ.
                                                 ˙ ˙                         (4.31)
                                     1

                                                                            ∂L
Now in order to apply the Euler-Lagrange solution it is necessary to find    ∂ uα
                                                                              ˙
                                                                                 .   In
order to avoid confusion of indices this is written as
         ∂L       ∂
            α
              =      (gρσ uρ uσ )1/2
                           ˙ ˙                                               (4.32)
         ∂u
          ˙     ∂ uα
                  ˙
                1                    ∂gρσ ρ σ         ∂ uρ
                                                        ˙         ∂ uσ
                                                                    ˙
              =   (gρσ uρ uσ )−1/2
                       ˙ ˙                 u u + gρσ α uσ + gρσ uρ α
                                           ˙ ˙              ˙   ˙            (4.33)
                2                    ∂ uα
                                       ˙              ∂u˙         ∂u˙
                1
              =   (gρσ uρ uσ )−1/2 (0 + gρσ δα uσ + gρσ uρ δα )
                       ˙ ˙                     ρ
                                                 ˙      ˙ σ                  (4.34)
                2
                1
              =   (gρσ uρ uσ )−1/2 (2gασ uσ );
                       ˙ ˙               ˙                                   (4.35)
                2
where we have used that facts that the metric gρσ is not a function of uα and
                                                                       ˙
          ρ
that gρσ δα = gασ . Furthermore we have

             ∂                 1                 ∂gρσ ρ σ
                (g uρ uσ )1/2 = (gρσ uρ uσ )−1/2
               α ρσ
                    ˙ ˙              ˙ ˙             ˙ ˙
                                                     u u +0+0 .              (4.36)
            ∂u                 2                 ∂uα

The Euler-Lagrange equation then becomes
                   d                          1                ∂gρσ
             0=      (gρσ uρ uσ )−1/2 gασ uσ − (gρσ uρ uσ )−1/2 α uρ uσ
                          ˙ ˙             ˙         ˙ ˙             ˙ ˙      (4.37)
                  dλ                          2                ∂u
4.3 Geodesics                                                                         33


Writing out the first term of this equation we get

 d                                 d                                            ∂gασ ∂uβ σ
   (gρσ uρ uσ )−1/2 gασ uσ
        ˙ ˙             ˙     =      [gρσ uρ uσ ]−1/2 gασ uσ + (gρσ uρ uσ )−1/2
                                          ˙ ˙             ˙         ˙ ˙                 ˙
                                                                                        u
dλ                                dλ                                            ∂uβ ∂λ
                                +(gρσ uρ uσ )−1/2 gασ uσ
                                      ˙ ˙             ¨                              (4.38)
                                                    ∂gασ
                              = 0 + (gρσ uρ uσ )−1/2 β uβ uσ
                                         ˙ ˙             ˙ ˙
                                                     ∂u
                                +(gρσ uρ uσ )−1/2 gασ uσ
                                      ˙ ˙             ¨                              (4.39)

    This is because (gρσ uρ uσ )1/2 = L and since L = L(x, x) is not a function
                         ˙ ˙                               ˙
       d                         d −1
of λ, dλ L = 0 and therefore dλ L = 0. This implies that L is constant along
the Euler-Lagrange curve and therefore equation (4.37) finally becomes (using
(4.39)):

                                     ∂gασ β σ
             0 = 0 + (gρσ uρ uσ )−1/2
                          ˙ ˙             u u + (gρσ uρ uσ )−1/2 gασ uσ
                                           ˙ ˙       ˙ ˙             ¨
                                      ∂uβ
                    1                ∂gµν
                  − (gρσ uρ uσ )−1/2 α uµ uν
                           ˙ ˙            ˙ ˙                                     (4.40)
                    2                ∂u
                                   ∂gασ β σ             1 ∂gρσ ρ σ
                = (gρσ uρ uσ )−1/2
                       ˙ ˙            β
                                        u u + gασ uσ −
                                        ˙ ˙       ¨             ˙ ˙
                                                               u u .              (4.41)
                                   ∂u                   2 ∂uα

Now with a little mathematical wizardry we are able to change some summation
indices, without influencing the content of the equation in the slightest, and use
the symmetry condition of the metric to show
                ∂gασ β σ ∂gαβ σ β 1 ∂gαβ σ β 1 ∂gασ β σ
                    ˙ ˙
                    u u =     ˙ ˙
                              u u =       u u +
                                          ˙ ˙         ˙ ˙
                                                      u u                         (4.42)
                ∂uβ       ∂uσ       2 ∂uσ       2 ∂uβ
and with this identity equation (4.41) becomes

                            1 ∂gρσ ρ σ 1 ∂gαρ σ ρ 1 ∂gασ ρ σ
              0 = gασ uσ −
                      ¨             ˙ ˙
                                   u u +       u u +
                                               ˙ ˙         ˙ ˙
                                                           u u                    (4.43)
                            2 ∂uα        2 ∂uσ       2 ∂uρ
                 = gασ uσ + Γαρσ uρ uσ
                       ¨         ˙ ˙                                              (4.44)

(see equation (4.17) above). But we can use the properties of metrics and the
                     µ
identity g µν gνλ = δλ and write

                             0 = g βα gασ uσ + g βα Γαρσ uρ uσ
                                          ¨              ˙ ˙                      (4.45)
                                   β      β
                                 ¨
                               = u + Γσρ                                          (4.46)

    These concepts are sufficient to understand the details of the derivations done
in later sections of this work. Other general relativistic concepts can be found in
the appendices.
            5. RELATIVISTIC PARTICLE ORBITS IN A
                 SCHWARZSCHILD GEOMETRY


                 5.1 Physical and Mathematical concepts
The motion of particles under the influence of a gravitational field is an impor-
tant aspect in the understanding of how light behaves under the same conditions.
Taking the relativistic corrections into account is essential to developing an ac-
curate picture of space-time curvature. For a star represented as an isolated,
non-rotating source of gravity the Einstein field equations can be solved to de-
scribe the geometry of the spacetime surrounding the star. The solution to this
problem gives rise to the Schwarzschild metric1 whose line element is given by
                                                  −1
                  2GM               2GM
  ds2 = − 1 −        2r
                        (cdt)2 + 1 − 2                 dr2 + r2 (dθ2 + sin2 (θ)dφ2 ) (5.1)
                   c                 cr
When this equation is compared with the static, weak-field metric

                       2Φ(xi )              2Φ(xi )
      ds2 = − 1 +              (cdt)2 + 1 −         (dx2 + dy 2 + dz 2 )             (5.2)
                         c2                   c2

one sees that they are the same when one realises that Cartesian coordinates
(x, y, z) have been converted to spherical (r, θ, φ) and the Newtonian potential
Φ is given as − GM . So M must be the total mass of the source of curvature.
                   r
This means that M is comprised of any source of mass as well as any source of
energy such as electromagnetic fields, nuclear interaction energy and even the
energy in the spacetime curvature itself. In this regime the spacetime curvature
for spherical symmetry is dependent only on the total mass M and not on how
it is distributed inside the source.
     Two fundamental properties of the Schwarzschild solution are time indepen-
dence and spherical symmetry. These properties give rise to, amongst others,
two vectors, known as Killing vectors, that highlight the symmetries in a system.
Named after the German physicist Wilhelm Killing, Killing vectors are vectors
that allow for the isolation of certain symmetries within a solution of Einstein
field equations and these two Killing vectors will be used extensively in later
sections.
  1                    ¨
   Schwarzschild, A., “Uber das Gravitationsfeld eines Massanpunktes nach der Einsteinschen
Theorie”, Proc. Royal Prus. Ac. Sci., 1, 189-196, (1916).
5.1 Physical and Mathematical concepts                                           35


                             5.1.1 Gravitational Redshift
Consider a stationary observer emitting a light signal with frequency ω at time
t◦ from a fixed Schwarzschild coordinate radius R. The photon will lose some
energy in climbing out of the gravitational well of the central mass. Since the
energy of a photon is related to its frequency, as measured by an observer through
E = hω; another stationary observer situated at r
      ¯                                                  R (or ≈ ∞), will measure
a lower frequency for the emitted photon. This is the gravitational redshift.
The magnitude of this effect is most easily shown using the time independence
property of the Schwarzschild geometry and the conserved quantity arising as a
result of this. When considering the two symmetries of spatial, (φ), independence
and time, (t), independence, the Killing vectors are ξ and η (not to be confused
with the η introduced in the Nordtvedt effect) and the conserved quantity due
to the spatial symmetry is given by
                                        ξ · u = const                          (5.3)
for a four-velocity u. It is therefore obvious that ξ · p must also be conserved for
the four-momentum p because p = mu. The energy of a photon as measured by
an observer is given by
                                    E = −p · uobs ,                            (5.4)
and thus
                                    hω = −p · uobs .
                                    ¯                                          (5.5)
But the observer is stationary and so the spatial components of the four-velocity
are zero and only the time component contributes anything to the frequency.
Normalisation of the four velocity yields
                       uobs(r) · uobs (r) = gαβ uα (r)uβ (r) = −1,
                                                 obs   obs                     (5.6)
and since all but the time component are zero, we have
                                  gtt [ut (r)]2 = −1.
                                        obs                                    (5.7)
The coefficient of the tt component of the metric, from equation (5.1), is
           −1/2
      2M
 1−    r
                  and thus in general

              uα (r) = [(1 − 2M/r)−1/2 , 0, 0, 0] = (1 − 2M/r)−1/2 ξ α ,
               obs                                                             (5.8)
while for a stationary observer at a radius r this becomes
                              uobs (r) = (1 − 2M/r)−1/2 ξ.                     (5.9)
If we combine this with equation (5.5) we can see that the frequency measured
by the stationary observer at R, the Schwarzschild radius, is given by
                                                −1/2
                                           2M
                           hωR = 1 −
                           ¯                            (−ξ · p)R .          (5.10)
                                            R
5.2 Derivation of relativistic particle orbits                                    36


Clearly for the observer at infinite radius

                                 hω∞ = (−ξ · p)∞ ,
                                 ¯                                            (5.11)

but it has already been noted that ξ · u is conserved and so must be the same at
infinity as at R. Therefore it is clear that the frequency redshift due to gravity
must be given by
                                            2M −1/2
                           ω∞ = ωR 1 −               .                     (5.12)
                                             R

              5.2 Derivation of relativistic particle orbits
In this analysis it has been seen that the two Killing vectors for spatial and time
independence are ξ and η, whose components are (1, 0, 0, 0) and (0, 0, 0, 1)
in the coordinate basis associated with the Schwarzschild solution. Consider the
4-velocity of the particle, u. Because of the conservative nature of Killing vectors,
the quantities formed by taking the products of the the 4-velocity and the two
Killing vectors are also conserved. They only return the dt and dz components
respectively. This is because the metric is both t and φ independent. These
quantities are given the names
                                                 2M    dt
                            e = −ξ · u = 1 −                                  (5.13)
                                                  r    dτ
and
                                                  dφ
                               l = η · u = r2 sin2 θ                    (5.14)
                                                  dτ
and represent the conservation of energy per unit mass and the conservation of
angular momentum per unit rest mass, respectively. This is seen from the fact
that for flat space the energy of the particle is given by

                                                 dt
                                E = mut = m                                   (5.15)
                                                 dτ

and for sufficiently large r, e approximately reduces to this. Similarly l becomes
the angular momentum per unit mass when dealing with sufficiently small veloc-
ities.
    For simplicity it is convenient to restrict one’s attention to particles moving
in a meridional ”plane” about the gravitating body. Consider the 4-velocity of
a particle at an instant in time. If the coordinates are oriented in such a way
that dφ/dτ = 0 at the instant that φ = 0, then l is zero and by the conservation
of angular momentum and equation (5.3), dφ/dτ must remain zero. As a result
the particle must remain within the meridional plane of φ = 0. Changing our
viewpoint slightly we could simply remain within the plane of θ = π/2 in which
case uθ = 0.
5.2 Derivation of relativistic particle orbits                                            37


   The 4-velocity is normalised such that

                            u · u = gαβ uα uβ                                          (5.16)
                                        dxα dxβ
                                  = gαβ         = −1,                                  (5.17)
                                         dτ dτ
and when the Schwarzschild metric is substituted in, one gets
                                                       −1
                 2M              2M
            − 1−    (ut )2 + 1 −                            (ur )2 + r2 (uφ )2 = −1.   (5.18)
                  r               r
Here it must be noted that geometrised units have been used such that G = c = 1.
Now, in general, uα = dα/dτ ; using that knowledge and equations (5.14) and
(5.15) to eliminate dt/dτ and dφ/dτ one gets
                            −1                          −1          2
                       2M                       2M             dr           l2
              − 1−               e2 + 1 −                               +      = −1.   (5.19)
                        r                        r             dτ           r2

Solving for e2 − 1 and dividing by two:
                                   2
               e2 − 1   1    dr            1           2M               l2
                      =                +        1−                 1+         −1 ,     (5.20)
                  2     2    dτ            2            r               r2

which illustrates the similarity to the Newtonian energy integral if one identifies
the term on the left as , the total energy per unit mass, and the second term
on the right as the effective potential Vef f (r). Thus orbits in the Schwarzschild
geometry can also be described using an equivalent of the Newtonian effective
potential, with the difference lying in the existence of a third, relativistic term.
This is seen by multiplying out the right-hand side term of equation (5.20) to get

                      2M            l2               2M   l2 2M l2
                   1−            1+ 2          −1 =−    + 2− 3                         (5.21)
                       r           r                  r  r    r

This third, cubic term prevents runaway to infinity for very small r as happens
in the Newtonian case but rather brings it back down to a finite value. This is
clearly shown in figure (5.1). What we thus have can be written as:
                                                 2
                                 1         dr
                              ε=                     + Vef f (r)                       (5.22)
                                 2         dτ

    What does the effective potential actually mean? In chapter 3 we saw that the
energy per unit mass for a body in orbit about another was given by equation
(3.62). That equation was the classical form that indicated a potential going
to infinity for small r. By implication no particle can approach a gravitating
source beyond a certain point in the Keplerian theory due to the runaway of the
5.2 Derivation of relativistic particle orbits                                         38




Fig. 5.1: Illustration of the Effective Potential for a particle a distance r from a source
          mass M , comparing the Newtonian and Relativistic cases.


potential to infinity. The Relativistic correction applied in this section introduces
a cubic term that serves to bring the potential back down for small values of r, see
figure (5.1). As a result it becomes possible to push a particle over that increased
potential until it falls into the well that captures the particle. By analysing this
potential for various energies we see in figure (5.1) that there exist four possible
types of orbit, depending on energy, labeled a, b, c and d.

   1. In (a) we have a particle approaching a gravitating source with sufficient
      energy to overcome the potential ’hump’ of the relativistic scenario but an
      impact parameter low enough so as to get caught. Such a particle will be
      pulled into the gravitating source and be absorbed in what is known as a
      plunge orbit.

   2. In scenario (b) we find two positions of equilibrium where a particle will
      orbit the parent body in a perfectly circular orbit. The upper position,
      however, is very unstable and a small perturbation in energy will send it
      either crashing into the parent body or flying off to infinity and so only
      in the case where the incident particle’s energy lies at the minimum will a
      stable circular orbit be found.

   3. Case (c) is known as a scattering orbit as the incident particle will approach
      from infinity, encounter the potential ’hump’, which it cannot overcome,
      and thus slide back out to infinity. The trajectory of the incident particle is
      merely bent about the parent body to continue on an altered, but otherwise
5.2 Derivation of relativistic particle orbits                                    39


      ordinary path. In this case the there exists an angle δ through which the
      trajectory has been bent and it is this angle that is fundamental to the
      Shapiro delay that will be discussed later. This is the type of orbit many
      comets follow as they are only temporarily captured by our sun on their
      journey through space.

   4. The last scenario is one where the particle rolls back and forth along between
      the points making up the line (d) in an orbit described as precessing. In
      such an orbit the particle’s eccentric orbit gradually moves about the parent
      body tracing out a large circle over time. This orbit is perhaps the strangest,
      but the existence of which came as the first proof of Einstein’s theory of
      Relativity when it explained the precession of Mercury’s perihelion. The
      precessing orbit is one where the positions of the perihelion (point of closest
      approach) and aphelion (point of farthest distance) gradually rotate about
      the parent body and trace out a circle with time. The precession of perihelia
      is another interesting test of gravitational theories and much effort has gone
      into understanding and predicting this effect.

   In following chapters we will see that for the case of light rays, not all these
orbits exist and that the scattering orbit is the one of most interest in the testing
of GR using signals from radio pulsars.
                                 6. SHAPIRO DELAY


                                     6.1 Background
The Shapiro delay (Shapiro, 1964) is a general relativistic effect occurring every-
where that an electromagnetic wave passes through the gravitational field gener-
ated by mass or energy density. The effect does, however, require a sufficiently
large mass to become detectable using current instruments. An electromagnetic
wave will always follow a path that in the geometrical approximation is the null
geodesic about a gravitating body. Keplerian geometry obviously does not take
into account the bending of spacetime and predicts a linear path of propagation
for the electromagnetic wave. There is thus a difference between the predicted
paths taken by the beam, for Keplerian and Einsteinian geometry respectively,
and therefore a corresponding difference in path length. Since, in the absence of
interfering media, the wave propagates at velocity c, there will be a difference
in propagation time for the beam. This difference in arrival time computed for
Keplerian and Einsteinian geometries is known as the Shapiro delay.
    This effect was accurately measured for the first time in the weak field regime
with the Viking mission to Mars in 19761 . Both of the landers and both of
the orbiters carried transponders that transmitted signals to earth at or close to
grazing incidence with the sun. The landers transmitted signals at the 10 cm
wavelength while the orbiters transmitted at both the 10 cm wavelength and at
3 cm. This was to compensate for the dispersive effect of the solar corona on
the propagating signals. The experiment was run when the Earth and Mars were
moving in opposing directions relative to one another i.e. on opposite sides of the
sun (See figure (6.1)). In this way the effect could be observed to grow and then
fade as the signal approached and then receded from grazing incidence with the
sun. That experiment measured the delay to be a maximum of approximately
250 µsec, a remarkable achievement given the fairly primitive equipment used
and the complexity of the problem. Using this result it was possible to constrain
the PPN parameter, γ, to 0.2% within Einstein’s theory. Further solar system
experiments2 constrained this parameter even further and a point was reached
  1
    Shapiro, I.I., Reasenberg, R.D., et. al., “The Viking Relativity Experiment”,J. Geophys.
Res., 82, 4329-4334, (1977).
  2
    Lebach, D. E., Corey, B. E., Shapiro, I. I., et al., “Measurement of the solar gravitational de-
flection of radio waves using very-long-baseline interferometry”, Phys. Rev. Lett., 75, 14391442,
(1995)
6.2 Mathematical derivation for the Shapiro time delay of light                           41




Fig. 6.1: Method for measuring the Shapiro time delay of light as performed during
          the Viking mission to Mars in 1976/7. (Shapiro, et al., 1977)


where only strong field tests could possibly allow for further constraints on the
theory.
    The discovery of the first binary pulsar system in 1974 (Hulse & Taylor,
1975) promised to be the system where the time delay could be measured in
the strong field3,4 but unfortunately the orbital elements weren’t favourable to
accurate measuring of the delay and more than thirty years would pass before
the challenge could be taken up again. In 2004, the first double pulsar binary
system was discovered (Lyne, et al., 2004) and allowed for the best strong field
testing of Relativity to date. The orbital elements are such that measurement of
the Shapiro time delay of light is measurable in the strong field regime and has
so far yielded excellent results. A detailed discussion of this effect follows below.


 6.2 Mathematical derivation for the Shapiro time delay of light
This derivation largely parallels that which was done to calculate the expected
delay in the signal during the Viking mission. A beam is emitted from Earth at
a certain time, reflected off a reflector the far side of the sun (in the case of the
Viking mission, at Mars) such that the beam passes as close to the sun as possible.
  3
    van Straten, W., et al., “A test of General Relativity from the Three-Dimensional Orbital
Geometry of a Binary Pulsar”, Nature, 412, 158-160, (2001)
  4
    Camilo, F., Foster, R.S., Wolszczan, A., “High Precision Timing of J1713+0747: Shapiro
Delay”, Ap. J., 437, L39-L42, (1994)
6.2 Mathematical derivation for the Shapiro time delay of light                                42


The round-trip travel time is then measured here on earth using an atomic clock.
In the case of pulsars we do not have the luxury of having reflectors in place but
fortunately pulsars are such stable clocks that, if our theories on the composition
of the pulsar are correct, the emission time of every pulse is extremely accurately
known and so the predicted arrival time equally well known. Any deviation
from this arrival time can therefore be attributed to some physical effect such as
the Shapiro delay and the highly predictable arrival times fulfil the role of the
reflector. Consequently, a similar derivation is sufficient to determine the Shapiro
delay for the Viking experiment and for binary pulsars.
    In the case of orbits of light rays, a very similar treatment is followed to that
done for particles in the previous chapter. The fundamental difference is that
light cannot be treated as a beam of particles with rest mass and so the concepts
of energy and momentum per unit mass have no meaning. This is not a problem
if one describes the world line of the light rays in terms of the family of affine
parameters that parameterise the path. Additionally, it is useful to deal with
the ratio between e and l as opposed to these parameters independently, in order
to eliminate the mass dependence of these parameters. As before, because the
Schwarzschild metric is independent in both t and φ, one gets
                                                    2M        dt
                            e ≡ −ξ · u = 1 −                                                 (6.1)
                                                     r        dλ
and
                                                 dφ
                                l ≡ η · u = r2 sin2 θ
                                                    ,                       (6.2)
                                                 dλ
where λ is the affine parameter that parameterises the path. The tangent vectors
that are described by this parameter must thus also be null because they are the
derivatives of points lying along the light path which is described by the affine
parameter of the curve that traces out the geodesics and therefore

                                             dxα dxβ
                               u · u = gαβ           = 0,                                    (6.3)
                                             dλ dλ
as opposed to the particle case where this equation equated to -1. Upon sub-
stituting in the Schwarzschild metric and focussing on the equatorial plane of
θ = π/2, one gets
                                2                −1            2                  2
               2M         dt             2M              dr             2    dφ
          − 1−                      + 1−                           +r                 = 0.   (6.4)
                r         dλ              r              dλ                  dλ

This reduces to
                               −1                   −1             2
                       2M                      2M         dr                l2
               − 1−                 e2 + 1 −                           +       = 0,          (6.5)
                        r                       r         dλ                r2
6.2 Mathematical derivation for the Shapiro time delay of light                    43


when equations (6.1) and (6.2) are used to eliminate the differentials in favour of
the e and l terms. Multiplying by (1 − 2M/r)/l2 and rewriting yields
                                           2
                         e2    1    dr             1      2M
                           2
                             = 2               +     2
                                                       1−    .                  (6.6)
                         l     l    dλ             r       r
This again has the form of a Newtonian energy integral where the second term
on the right is the effective potential and the ’energy’ is given by (l2 /e2 )−1 which,
for reasons that will become obvious, is designated 1/b2 . What is this b2 term? If
the light ray is sufficiently far away from the source of gravitation then one can
approximate the surrounding space-time to be flat and the Schwarzschild polar
coordinates can be replaced by ordinary Cartesian coordinates. For a light ray
traveling parallel to the x-axis a distance d from it;
                                   l   r2 dφ/dλ     dφ
                            b≡       ≈          = r2 .                          (6.7)
                                   e    dt/dλ       dt
Geometrically it is obvious that for very large r, φ ≈ d/r and dr/dt ≈ −1.
Therefore
                                 dφ    dφ dr     d
                                    =        = 2,                              (6.8)
                                 dt    dr dt     r
which shows that b = d. Therefore b is just the impact parameter of the incident
light ray approaching a source of gravitation from infinity. This impact parameter
b is in geometrised units and determines the various types of orbit a light ray can
take about a gravitating body. Unlike for particle orbits, a light ray can only
follow circular, scattering or plunge orbits although it must be mentioned that
the circular orbit is only possible about a black hole; see below. For neutron
stars, that leaves the scattering and plunge orbits that depend upon the impact
parameter, b. As a result, equation (6.6) can be written as:
                                                   2
                             1     1       dr
                               2
                                 = 2                   + Vef f (r),             (6.9)
                             b     l       dλ
for
                                         1       2M
                             Vef f (r) =     1−         .                      (6.10)
                                         r2        r
This effective potential clearly contains a cubic as well as a quadratic term but
not a linear one. Runaway to positive infinite potential still does not happen
although the shape of the effective potential graph is slightly different to that of
the particle case as is seen in figure (6.2). The light ray orbits possible due to the
Schwarzschild geometry parallel those available to particles in the Schwarzschild
geometry, with the difference that precessing orbits are not available to light rays.
This means that circular, scattering and plunge orbits are all possible for light
rays although the circular orbit is energetically unstable as it only occurs for the
maximum effective potential and so will quickly decay into a scattering or plunge
orbit with any perturbation in energy.
6.2 Mathematical derivation for the Shapiro time delay of light                     44




Fig. 6.2: Effective potential for a light ray approaching a gravitating source of mass M
          at a distance r

                               6.2.1 Deflection angle
It is clear that all material bodies bend light to some degree and this bending
results only in scattering and plunge orbits. For the case of the Shapiro delay,
the plunge orbit can be ignored as an absorbed photon cannot be detected by an
earthbound observer and so we will only consider the case of the scattering orbit.
When dealing with a scattering orbit, the light ray will be deflected through an
angle δφdef as shown in figure (6.3). In order to determine the deflection of the
light one needs to know how the angle φ changes as a function of distance from
the gravitational source, or dφ/dr. Once this differential equation is obtained it
can be easily solved to find the desired relationship between r and the deflection
angle φ. To do this, equations (6.2) and (6.6) are solved for dφ/dλ and dr/dλ,
respectively, and one divided into the other. The result is:
                                                              −1/2
                     dφ                   l2
                        = ±l(r2 sin2 θ)−1 2 − l2 Vef f (r)                      (6.11)
                     dr                   b

Simplifying yields
                                                   −1/2
                        dφ       1 1
                            = ± 2 2 − Vef f (r)         .                  (6.12)
                        dr       r b
where the sign gives the direction of the orbit. Defining the total angle swept
out by the beam as it approaches from infinity, is bent about the source and
continues onwards as ∆φ with the turning point being the bisector of the angle.
Thus ∆φ is the same as twice the angle swept out from the turning point, or point
6.2 Mathematical derivation for the Shapiro time delay of light                                 45


of closest approach r1 , to infinity. This can is seen in figure (6.3). Substituting




Fig. 6.3: Deflection angle for an incident electromagnetic wave upon a source of gravi-
          tation

in the definition of Vef f (r) this is then:

                                        ∞                              −1/2
                                            dr 1    1    2M
                          ∆φ = 2                  − 2 1−                                     (6.13)
                                       r1   r 2 b2 r      r

The point r1 is where 1/b2 = Vef f (r1 ), so by introducing a new variable r ≡ b/w,
the previous expression becomes:
                                                                              −1/2
                            b dw w2 1
                              ∞       w2  2M
               ∆φ = 2      − 2 2 2 − 2 1−    w                                       ;       (6.14)
                      b/w1   w b b    b    b

that easily reduces to:
                                       w1                              −1/2
                                                            2M
                          ∆φ = 2            dw 1 − w2 1 −      w              .              (6.15)
                                   o                         b
In this equation w1 is the point where 1/b2 = Vef f (r1 ) as before. It should be
clear that the value 2M/b determines the amount through which the beam is
bent. In the case of our sun, b ≥ R (= 6.96 × 105 km) and M = M (= 1.47
km) which means that 2M/b ≈ 10−6 , a very small number. As a result we can
expand equation (6.15) in powers of 2M/b. We first re-write (6.15) as:
                                               −1/2               −1              −1/2
                     w1            2M                      2M
       ∆φ = 2             dw 1 −      w               1−      w        − w2              .   (6.16)
                 o                  b                       b
6.2 Mathematical derivation for the Shapiro time delay of light                      46


Expand both inverse powers in terms of 2M/b to yield
                                                                                      −1/2
             w1         M    3 M2 2                       2M      M2
∆φ = 2            dw 1 + w +      w + ...              1+    w + 4 2 w2 + ... − w2           .
         o              b    2 b2                          b       b
                                                                               (6.17)
To first order this is
                                                          M
                                      w1             1+   b
                                                            w
                         ∆φ = 2            dw                       1/2
                                                                          .       (6.18)
                                  o                  2M
                                                1+    b
                                                        w   − w2

This integral can be now be looked up in integral tables to show that the final
solution is simply
                                          4M
                                ∆φ ≈ π +     ,                           (6.19)
                                           b
for small M/b. The relationship between ∆φ and δφdef is clearly shown in figure
(6.3) and is simply
                               δφdef = ∆φ − π.                           (6.20)
Finally the relativistic deflection of light can be written out; Remembering that
we used geometrised units for the derivation and replacing the correct factors for
G and c, this is
                               4GM
                      δφdef =             (for small GM/bc2 ).             (6.21)
                                bc2

                          6.2.2 The Shapiro delay of light
The derivation of the equation giving the amount by which a source will delay a
light beam is now shown. This derivation proceeds in a similar manner to that
done for the deflection angle. We require the differential equation dt/dr in terms
of the parameters e and l and this is obtained by solving equation (6.1) for dt/dλ,
equation (6.6) for dr/dλ and dividing the second into the first to obtain
                                                −1                  −1/2
                        dt   1   2M                   1
                           =± 1−                         − Vef f              ,   (6.22)
                        dr   b    r                   b2
where the ± indicates increasing or decreasing radius. In this derivation we are in-
terested in the variation of the total travel time for the beam in the Schwarzschild
coordinates centered on the sun. The total travel time for the beam is equal to
the travel time taken from the source of the signal at rp to the point of closest
approach indicated by r1 , plus the time taken from that point of closest approach
to the detector on earth at rR . In other words we measure time both forwards and
backwards from the point of closest approach to the source of curvature defined
as r1 . This is written as:

                            (∆t)total = t(rp , r1 ) + t(r1 , rR )                 (6.23)
6.2 Mathematical derivation for the Shapiro time delay of light                                   47


where in general
                                   r                    −1                    −1/2
                                            1    2M          1
                    t(r, r1 ) =        dr     1−                − Vef f (r)          .         (6.24)
                                  r1        b     r          b2
The expansion of this is then
                    r             2M  4M 2           1               3
 t(r, r1 ) =            dr 1 +       + 2 + ...                           2
                                                 1 + b2 Vef f (r) + b4 Vef f (r) + ...
                    r1             r   r             2               8
                                                                                  (6.25)
                     r                               2
                              2M  1               M b Vef f (r)
            =          dr 1 +    + b2 Vef f (r) +               + ... .           (6.26)
                    r1         r  2                    r
Also we know that
                                  1                   1     2M
                                    2
                                      = Vef f (r1 ) = 2 1 −                                    (6.27)
                                  b                  r1      r1
therefore
                                                             −1/2
                                                     2M
                                        b = r1    1−                                           (6.28)
                                                      r1
and, to first order in M,
                                         M
                                             + ...).
                                        b = r1 (1 +                        (6.29)
                                          r1
To first order we can use this to eliminate b from the integral and integrate it to
get
                                                             
                                                             2
                              2
                                                 r+   r 2 − r1     r − r1            1/2
      t(r, r1 ) =       r2 − r1 + 2M ln                       +M                         .   (6.30)
                                                      r1           r + r1

In this equation we find that the first term is the Newtonian contribution to the
propagation time that would be present even in the absence of the gravitating
source. The other terms are the relativistic corrections due to the curvature of
the spacetime in the vicinity of the source. The total travel time for the electro-
magnetic signal is then given by substituting equation (6.30) into equation (6.23)
for the various radial parameters of interest. These parameters are rR which is
the distance from the source of curvature to the detector i.e. the earth, rp is the
distance from the source of curvature to the source of the electromagnetic ray i.e.
the pulsar; and r1 is the distance of closest approach between the electromagnetic
ray and the gravitating body as shown in figure (6.3). Using these parameters
the excess time of flight of the signal, which is the Shapiro delay, in the weak
field is written as
                                              2GM     4rR rp
                            ∆tShapiro ≈          3
                                                   ln    2
                                                                    +1 ,                       (6.31)
                                               c        r1
upon un-geometrising.
                  7. STRONG-FIELD SHAPIRO DELAY


The necessary extension to the previous chapter is the elucidation of the Shapiro
delay in the strong field regime. Such a regime is found in the vicinity of a
double neutron star system consisting of a pulsar and a neutron star, a rarer
pulsar-pulsar binary or possibly even the theorised pulsar-black hole binary. One
of the first analyses of this type of compact binary system was done by Blandford
and Teukolsky (1976)1 who corrected earlier work by Wheeler2 and although
refinements have been made since,3,4 the original analysis has been used here.
    In order to simplify the problem somewhat, the relativistic effects are regarded
as small perturbations of the classical orbit and it is assumed that this orbit is a
slowly precessing Keplerian ellipse. The reason for this approximation is that the
pulsars are so small compared to the orbit that they may be regarded as point
particles in a rotating frame centred on the centre of mass of the binary system.
When dealing with rotating frames, the Kerr metric is to be used but due to
the approximation made, a different version is used to solve the problem. This
version is1 :

 ds2 = −[1+2Φ+O(v 4 )]dt2 +O(v 3 )dxi dt+[1−2Φ+O(v 4 )](dx2 +dy 2 +dz 2 ). (7.1)

where O(v 4 ) represents all the higher order terms. To first approximation, the
Newtonian potential Φ can be regarded as the sum of two Newtonian potentials
centred on each of the neutron stars in the binary, hence
                                                −M1           M2
                     Φ(r, t) = Φ1 + Φ2 =                 −             ,                  (7.2)
                                             |r − r1 (t)| |r − r2 (t)|

using geometrised units and with M1 and M2 representing the masses of the
pulsar and it’s companion, respectively. One can then construct the equatorial
plane of a polar coordinate system such that it coincides with the orbital plane
  1
     Blandford, R., Teukolsky, S.A., “Arrival-time analysis for a pulsar in a binary system”,Ap.
J., 205, 580-591, (1976).
   2
     Wheeler, J.C., “Timing effects in Pulsed Binary Systems”, Ap. J. (Letters), 196, L67,
(1975).
   3
     Haugan, M., “Post-Newtonian Arrival-Time Analysis for a Pulsar in a Binary System”,
Ap. J., 296, 1-12, (1985).
   4
     Backer, D., Hellings, R., “Pulsar Timing and General Relativity”, Ann. Rev. Astron.
Astrophys., 24, 537-575, (1986).
                                                                                        49


of the binary system and in doing so the (r, θ, φ) coordinates of the pulsar and
companion are simply,

                      r1 = (r1 , π/2, φ),    r2 = (r2 , π/2, φ + π).                  (7.3)

In these equations, φ is the true anomaly as described in Chapter 3 and r1 and
r2 relate the positions of each object within their own ellipses and given by,
                       M2                     M1                 a(1 − e2 )
              r1 =            r,    r2 =            r,     r=               .         (7.4)
                     M1 + M 2               M1 + M2             1 + e cos φ
The time that would be measured by an ideal clock on the pulsar, the proper
time Tp , is determined in terms of the coordinate time, t, and is is found from
the Kerr metric defined above. The relationship is,
                                                           1
        (dT )2 = −ds2 = dt2 1 + 2Φ + O(v 4 ) −                (dx2 + dy 2 + dz 2 )
                                                          dt2
                                           2
                = dt2 [1 + 2Φ + O(v 4 ) − v1 ].                                       (7.5)

This can be re-written as,
                           dTp                        2
                               = [1 + 2Φ + O(v 4 ) − v1 ]1/2                          (7.6)
                            dt
which upon performing the Taylor expansion and incorporating higher order
terms into O(v 4 ) yields
                          dTp               1 2
                              = 1 + Φ(r1 ) − v1 + O(v 4 ).                            (7.7)
                           dt               2
From (7.2) it may appear that
                                                  −M1                M2
          Φ(r1 ) = Φ1 (r1 ) + Φ2 (r1 ) =                     −                  .     (7.8)
                                            |r1 (t) − r1 (t)| |r1 (t) − r2 (t)|
This would imply that Φ1 (r1 ) is infinite. However, Φ should be evaluated at the
surface of the pulsar i.e. at r = r1 + ε. It is thus a constant gravitational redshift
term. Also, − 1 v1 can be regarded as the transverse Doppler shift and is defined
               2
                 2

by
                                         2
                               2      M2         2 1
                              v1 =                 −     ,                       (7.9)
                                   M1 + M2 r a
as was seen in chapter 3(cf. the end of section 3.2). Upon dropping the constant
terms in equation (7.7) one obtains
                                              2
                           dTp      M2      M2     1
                               = 1−    −                                             (7.10)
                            dt      r    M1 + M2 r
                                    M2 (M1 + 2M2 )
                               = 1−                                                  (7.11)
                                    r M1 + M2
                                                                             50


The reason we can drop the constant terms becomes apparent after looking at
the timing equation for time of emission of the N -th pulse from a pulsar with a
rotational frequency ν,
                                      1 2 1 3
                        N = N0 + νTp + νTp + ν Tp .
                                        ˙     ¨                           (7.12)
                                      2     6
We want to find Tp and so after integration, any overall multiplicative constants
affecting Tp emerging from the integration can be absorbed into the ν of the
timing equation while any additive constants that come out can be absorbed into
N0 without affecting the final form of the result. From the mechanics introduced
in chapter 3, we know that

                                 r = a(1 − e cos E),                      (7.13)

for E the eccentric anomaly. Equation (7.11) can thus be written as

                     dTp     M2 M1 + 2M2     1
                         =1−                                              (7.14)
                      dt     a M1 + M2 (1 − e cos E)

since the constant term can be dropped. It has been stated in chapter 3 that the
mean anomaly is related to the eccentric anomaly via:

                                 M = E − e sin E                          (7.15)

But we know that for a binary with orbital period Pb ,
                                        2π
                                  M=       (t − t0 )                      (7.16)
                                        Pb
and since (2π/Pb )t0 can be dropped because it is a constant that relates to the
time of periastron passage, we have
                                   Pb    Pb
                             t=       E−    e sin E.                      (7.17)
                                   2π    2π
Therefore:
                                 Pb      Pb
                          dt =      dE −    e cos EdE;                    (7.18)
                                 2π      2π
such that
                                  M2 M1 + 2M2     1
                  dTp = dt 1 −                            ,               (7.19)
                                  a M1 + M2 (1 − e cos E)
which can finally be written as
                                      M2 M1 + 2M2 Pb
                     dTp = dt 1 −                    dE .                 (7.20)
                                      a M1 + M2 2π
                                                                                       51


Equation (7.14) can therefore be integrated to yield
                                        M2 M1 + 2M2 Pb E
                             Tp = t −                                              (7.21)
                                        a M1 + M2 2π
after one has transformed the integration variable from dt to dE as was shown
above.


   The propagating signal has to cross the large distances of space and in doing
so will be affected by the interstellar medium’s electron content. It will therefore
travel at a group velocity of less than the speed of light, which in the geometrised
units in which we are working is unity. Define from the line element (7.1)
                         ˆ
                        dt = (1 + Φ)dt,        |dˆ| = (1 − Φ)|dx|,
                                                 x                                 (7.22)
in order to accurately measure coordinate time at a constant point and to accu-
rately measure distances at a constant coordinate time. This lower velocity can
then be described by
                                    x
                                   dˆ
                                       =1− .                                 (7.23)
                                    ˆ
                                   dt
And so
                                dx
                                    = 1 − + 2Φ.                              (7.24)
                                dt
We now integrate along a straight line joining re to r1 by approximating the
world line of the photon as a straight line. This is because light rays follow null
geodesics and in the absence of a large gravitating source to warp the spacetime
between the pulsar and the earth, this path will be approximately straight. We
further integrate between the time of emission and the time of arrival to obtain:
                                         re (tarr )
                        tarr − tem =                  (1 + − 2Φ)|dx|.              (7.25)
                                        r1 (tem )

In this equation there are three terms to be integrated. The constant is a trivial
matter. The will result in a term that contains the frequency of the signal as
well as the dispersion constant, D, that is a function of the electron content, ne ,
along the path to be integrated and is thus itself an integral given by:
                                         e2           re
                                 D=                        ne |dx|.                (7.26)
                                        2πm         r1

Electron content models are notoriously difficult to obtain due to the variations
constantly occurring in the interstellar medium but the most common model
currently used is that of Cordes and Lazio 5 .
  5
    Cordes, J.M., “NE2001: A New Model for the Galactic Electron Density and it’s Fluctua-
tions” in Milky Way Surveys: The Structure and Evolution of our Galaxy, Proceedings of ASP
Conference 317, Eds, Clemens, D., Shah, R., Brainerd, T., p211, (2004)
                                                                                                   52


    It is however, the integral of Φ that is responsible for the relativistic time
delay across the orbit and it is this value that varies over the orbit while the
others depend only on the beam path taken. If one focuses only on this integral
and remembers that Φ1 can be ignored because it only adds a constant term, we
have only the contribution from Φ2 and it becomes possible to write the integral
as,
                         re               tarr      dt
                    −2      Φ|dx| ≈ 2M2                         .            (7.27)
                        r1               tem |x(t) − r2 (tem )|

Since we are integrating along a straight path along the line joining re and r1 ,
we know that the tangent vector to the curve has constant direction and so we
can write

               [x(t) − r1 (tem )][tarr − tem ] = [re (tarr ) − r1 (tem )][t − tem ].           (7.28)

Therefore
                                     t − tem
                     x(t) = r1 (tem ) +       [re (tarr ) − r1 (tem )],                        (7.29)
                                   tarr − tem
which is substituted into the integral to obtain
       tarr          dt                    tarr                              dt
2M2                              = 2M2                             t−tem
      tem     |x(t) − r2 (tem )|          tem     |r1 (tem ) +            [r (t )
                                                                 tarr −tem e arr
                                                                                    − r1 (tem )] − r2 (tem )|
                                                                                               (7.30)
The solution to this integral is then

            tarr − tem |re − r1 ||r + re − r1 | + |re − r1 |2 + r · (re − r1 )
      2M2               ln                                                     .               (7.31)
             |re − r1 |              |re − r1 |r + r · (re − r1 )

Looking at equation (7.29), to first approximation tarr − tem = |re − r1 | and the
expression outside the ln is simplified. Furthermore re          r1 and so the final
solution for the relativistic time delay in the strong field regime is,
                                                        2re
                                 ∆trel = 2M2 ln                                                (7.32)
                                                       r+r·n
where n = re and represents the unit vector pointing to the solar system barycen-
           re
ter. The distance between the binary system and the earth is practically constant
and so the numerator in the above expression is effectively constant and can be
ignored. The dot-product term is expressible as −r sin (ω + φ) sin i using the
Keplerian mechanics of chapter 3 and since r = 1 + e cos φ, we have
                                                    1 + e cos φ
                           ∆trel = 2M2 ln                             .                        (7.33)
                                                1 − sin i sin (ω + φ)

It must be noted that, from the start, this analysis was done assuming the cor-
rectness of General Relativity. Fortunately the analysis varies very little within
                                                                             53


the PPN framework, only introducing the parameter γ that replaces the number
2 in equation (7.33)with the expression (1 + γ). In GR γ takes on the value 1
while in certain other theories it is other than unity. The final equation in the
PPN formalism for the relativistic time delay of a pulsar signal in a compact
binary system is thus
                                    GM2        1 + e cos φ
                  ∆trel = (1 + γ)       ln                       ,        (7.34)
                                     c2    1 − sin i sin (ω + φ)
once we un-geometrise.
                    8. RESULTS OF SIMULATIONS


                          8.1 Procedure followed
The results displayed in the following sections were done with a program writ-
ten in the MATLAB programming language, for which the code can be found in
appendix B. A program was written to accommodate those orbital parameters
that are relatively easy to determine through simple observation of the pulsar in
its orbit. These parameters, such as the orbital period, eccentricity and semi-
major axis length, were then used to determine the relativistic effect of interest,
the Shapiro time delay. Essential parameters in this calculation were the in-
clination of the orbit and the mass of the companion body. These, although
extremely difficult to determine experimentally, were also included and allowed
to vary through a wide range of values so as to see the effects on the Shapiro
time delay. As mentioned previously, it has proven experimentally very difficult
to measure the Shapiro delay due to the fact that it varies quasi-sinusoidally with
orbital phase. This in itself is not problematic but the effect virtually disappears
into the Roemer time-of-flight delay caused by the finite size of the orbit and
finite speed of light. Because the pulsar is in orbit, its signal has to propagate
periodically through larger and then shorter distances to earth. This results in
a quasi-sinusoidal variation in the arrival time of the signal as it is received on
earth. Furthermore the magnitude of this signal is much greater than that of
the Shapiro delay and through this combination of similar periodicity and differ-
ing amplitude, the Shapiro delay signal is easily buried within the Roemer delay
signal and lost. For this reason the Roemer delay was also modeled so that the
contrast in signal could be observed.
    Ordinarily the Roemer delay refers to the annual sinusoidal variation in arrival
time due to the orbit of the earth about the sun but in the case of binary pulsars,
the timescale of variation of the signal of the pulsar in its orbit is much smaller
than that of the solar system Roemer delay. For this reason the Roemer delay
referred to in this thesis represents the time of flight delay of only the binary
system signal. In order to understand all the possible effects a large number of
simulations were run in which the various orbital parameters were varied so as
to note any trends in the signal that might be observed using current and future
equipment, and identify promising systems types.
8.2 Graphs and Results                                                         55


                         8.2 Graphs and Results
The first set of 6 results shows the comparison of the Roemer and Shapiro delays
in a simulated, yet typical, compact binary system where both components have
a mass of 1.4M . The semi-major axis of the pulsar’s orbit is set at a typical
value of 1 000 000 km. Each figure of the 6 figures (8.1-8.6) below contains 6
graphs. Successive graphs within a figure have higher inclination angles and show
both the Roemer and Shapiro time delays on the same set of axes. The Shapiro
delay has been scaled up 105 times to make it visible and to show the similarities
in phase of the two curves. Successive figures within the set (8.1-8.6) have a
higher eccentricity to show how the variation of this parameter affects the curves
by inducing a phase shift between the two competing curves.
8.2 Graphs and Results                                                                  56




Fig. 8.1: The eccentricity of the system in this set is 0. Please note that the two delays
          are not on the same scale. The Shapiro delay is 105 times smaller than the
          Roemer delay. The top left graph is at an inclination angle of 15◦ with the
          top right being at an inclination of 30◦ . The middle left graph is then at 45◦ ,
          the middle right at 60◦ , the bottom left at 75◦ and finally the bottom right is
          at an inclination angle of 90◦ .
8.2 Graphs and Results                                                                57




Fig. 8.2: The eccentricity of the system in this set is 0.1. Please note that the two
          delays are not on the same scale. The Shapiro delay is 105 times smaller than
          the Roemer delay. The top left graph is at an inclination angle of 15◦ with
          the top right being at an inclination of 30◦ . The middle left graph is then at
          45◦ , the middle right at 60◦ , the bottom left at 75◦ and finally the bottom
          right is at an inclination angle of 90◦ .
8.2 Graphs and Results                                                                58




Fig. 8.3: The eccentricity of the system in this set is 0.2. Please note that the two
          delays are not on the same scale. The Shapiro delay is 105 times smaller than
          the Roemer delay. The top left graph is at an inclination angle of 15◦ with
          the top right being at an inclination of 30◦ . The middle left graph is then at
          45◦ , the middle right at 60◦ , the bottom left at 75◦ and finally the bottom
          right is at an inclination angle of 90◦ .
8.2 Graphs and Results                                                                59




Fig. 8.4: The eccentricity of the system in this set is 0.4. Please note that the two
          delays are not on the same scale. The Shapiro delay is 105 times smaller than
          the Roemer delay. The top left graph is at an inclination angle of 15◦ with
          the top right being at an inclination of 30◦ . The middle left graph is then at
          45◦ , the middle right at 60◦ , the bottom left at 75◦ and finally the bottom
          right is at an inclination angle of 90◦ .
8.2 Graphs and Results                                                                60




Fig. 8.5: The eccentricity of the system in this set is 0.6. Please note that the two
          delays are not on the same scale. The Shapiro delay is 105 times smaller than
          the Roemer delay. The top left graph is at an inclination angle of 15◦ with
          the top right being at an inclination of 30◦ . The middle left graph is then at
          45◦ , the middle right at 60◦ , the bottom left at 75◦ and finally the bottom
          right is at an inclination angle of 90◦ .
8.2 Graphs and Results                                                                61




Fig. 8.6: The eccentricity of the system in this set is 0.8. Please note that the two
          delays are not on the same scale. The Shapiro delay is 105 times smaller than
          the Roemer delay. The top left graph is at an inclination angle of 15◦ with
          the top right being at an inclination of 30◦ . The middle left graph is then at
          45◦ , the middle right at 60◦ , the bottom left at 75◦ and finally the bottom
          right is at an inclination angle of 90◦ .
8.2 Graphs and Results                                                               62


    The figures (8.1 - 8.6) can be summarised best by combining the Shapiro
delays for varying eccentricities in a single plot. In figures (8.7) and (8.8) we
see how the Shapiro delay varies in a system with inclination 90◦ and one with
inclination 45◦ . It has been noted that the Shapiro delay is most easily measured
in systems with high inclination angles because of the sharp peak in the delay
curve and in this instance it can be seen that increasing eccentricity has little
effect. Only in the lower inclination systems does an increase in eccentricity lead
to a marked phase shift that should make it separable from the Roemer delay. In
figures (8.7) and (8.8) a simulated system was used with a 10M companion in
order to maximise the effect of the Shapiro delay and a semi-major axis length
of 3.3 ls in order to minimise the Roemer delay for ease of identification.




Fig. 8.7: Variation in Shapiro delay with an increase in the eccentricity of a system at
          an inclination angle of 90◦ and keeping all other orbital parameters constant.
8.2 Graphs and Results                                                               63




Fig. 8.8: Variation in Shapiro delay with an increase in the eccentricity of a system at
          an inclination angle of 45◦ and keeping all other orbital parameters constant.
8.2 Graphs and Results                                                             64


    These figures have illustrated the fact that the Roemer and Shapiro delays
both vary sinusoidally with phase and as a result it is extremely problematic to
distinguish the Shapiro delay from the Roemer delay owing to its much smaller
(10−5 s) size. It is in systems with a high inclination angle that it becomes
more favourable to measure the Shapiro delay due to the sharp peak that forms
when nearing an eclipsing orbit. It is also apparent that a high eccentricity
causes a distinct phase shift between the Roemer and Shapiro delays, making
the Shapiro delay distinguishable under certain circumstances. It is exactly these
circumstances that are important to understand so that future measurements
of the delay can be made in a wider range of binary system. It can be seen
from figure (8.10) that a decreasing semi-major axis length, i.e. an increasingly
compact binary system, leads to a decrease in the magnitude of the Roemer delay.
Fortunately a decrease in a does not affect the Shapiro delay and thus extremely
compact binary systems are ideal for separating these two effects. Furthermore
an increasing companion mass leads to an increasing Shapiro delay, figure (8.9),
without affecting the Roemer delay directly. Of course an increasing companion
mass would alter the orbital dynamics and thereby affect the Roemer delay; but
in the absence of any evolution of the orbital parameters due to the increase
in companion mass, say through the process of accretion, an increasingly heavy
companion will greatly assist in determination of the Shapiro delay.




Fig. 8.9: Figure illustrating the increase in the Shapiro delay as the companion mass,
          mc increases. Note that mc is given in units of solar masses.
8.2 Graphs and Results                                                               65




Fig. 8.10: Figure illustrating the increase in the Roemer delay as the semi-major axis
           length, a, increases. Note that a is given in units of light-seconds.




Fig. 8.11: Figure illustrating the variation in the Roemer delay as the eccentricity in-
           creases in a system with an inclination of 45◦ . This effect is shown for a
           binary system with a companion mass of 1.4M .
8.2 Graphs and Results                                                                 66




Fig. 8.12: Figure illustrating the variation in the Shapiro delay as the position of peri-
           astron (θ◦ ) is rotated about the center of mass from 90◦ to the right of the
           line of sight to 90◦ to the left of the line of sight. The system used was one
           with an eccentricity of 0.9, a companion mass of 1.4M and an inclination
           of the orbit of 45◦ .




Fig. 8.13: Figure illustrating the variation in the Roemer delay as the position of pe-
           riastron θ◦ is rotated about the center of mass. The system used was one
           with an eccentricity of 0.9, a companion mass of 1.4M and an inclination
           of the orbit of 45◦ .
8.3 Predicted Future Results                                                     67


                      8.3 Predicted Future Results
The following set of graphs parallels the first 6 shown in this chapter, in com-
paring Shapiro and Roemer delays in a binary system with all parameters except
inclination and eccentricity fixed. Each set of the 6 figures is plotted at a constant
eccentricity with the inclination varying from 15o to 90o as one moves through
the graphs going from left to right and top to bottom. Each consecutive figure
increases the eccentricity. The major difference in these graphs is a companion
mass, mc , equal to 10M . This means our companion is a low mass black hole
and the system is the very exotic and as yet unfound pulsar-black hole binary.
Furthermore the simulation was done using a compact, but not impossible semi-
major axis length of 2 000 000 km so as to minimise the Roemer delay of the
system. The reason for including this set of data is that a pulsar-black hole
binary is widely regarded as the holy grail of relativistic astrophysics and some
indication of the expected results would assist in the identification of the Shapiro
delay in such systems.
8.3 Predicted Future Results                                                         68




Fig. 8.14: The following set of 6 results shows the expected Roemer and Shapiro delays
           for a hypothetical system containing a black-hole as the companion object.
           Successive figures have increasing eccentricity while sequential graphs within
           each figure have increasing inclination angles. The eccentricity of the system
           in this set is 0. Please note that the two delays are not on the same scale.
           The Shapiro delay is 104 times smaller than the Roemer delay.
8.3 Predicted Future Results                                                       69




Fig. 8.15: The eccentricity of the system in this set is 0.2. Please note that the two
           delays are not on the same scale. The Shapiro delay is 104 times smaller
           than the Roemer delay.
8.3 Predicted Future Results                                                       70




Fig. 8.16: The eccentricity of the system in this set is 0.4. Please note that the two
           delays are not on the same scale. The Shapiro delay is 104 times smaller
           than the Roemer delay.
8.3 Predicted Future Results                                                       71




Fig. 8.17: The eccentricity of the system in this set is 0.6. Please note that the two
           delays are not on the same scale. The Shapiro delay is 104 times smaller
           than the Roemer delay.
8.3 Predicted Future Results                                                       72




Fig. 8.18: The eccentricity of the system in this set is 0.8. Please note that the two
           delays are not on the same scale. The Shapiro delay is 104 times smaller
           than the Roemer delay.
8.3 Predicted Future Results                                                       73




Fig. 8.19: The eccentricity of the system in this set is 0.9. Please note that the two
           delays are not on the same scale. The Shapiro delay is 104 times smaller
           than the Roemer delay.
8.4 3-Dimensional Parameter spaces                                          74


                 8.4 3-Dimensional Parameter spaces
In the following series of surface plots, the various parameters of the Shapiro
and Roemer delays have been plotted against one another to generate surfaces of
maximal delay. In this manner it is possible to see in 3-dimensions the regions
where the Shapiro delay is maximal and contrast it with where the Roemer delay
is minimal.
8.4 3-Dimensional Parameter spaces                                                         75




Fig. 8.20: This 2-dimensional plot shows the combined effect of eccentricity and incli-
           nation on the Shapiro delay for a standard compact binary system. What
           is meant by a standard binary system is one in which the companion is a
           1.4M neutron star, the semi-major axis length is 6.6 ls or 2 000 000 km, the
           periastron is positioned at 45◦ to the right of the line of sight, the eccentricity
           is 0.3 and the inclination angle is 45◦ . In all further simulations these values
           are used except for where the two parameters of interest are varied over the
           full range of possible values such as the eccentricity and inclination angle in
           this figure. In this plot it is seen that increasing inclination is the dominant
           parameter in the Shapiro delay but it must be noted that there does exist
           an increase due to increasing eccentricity, particularly for lower inclination
           angles as evidenced by the slope of the projected contour lines.
8.4 3-Dimensional Parameter spaces                                                 76




Fig. 8.21: This 2-dimensional plot shows the combined effect of eccentricity and incli-
           nation on the Roemer delay for a standard compact binary system.
8.4 3-Dimensional Parameter spaces                                                   77




Fig. 8.22: This 2-dimensional plot shows the combined effect of eccentricity and the
           position of the periastron on the Shapiro delay for a standard compact binary
           system. There is a clear symmetry about 1.6 rad (or 90◦ ) as at this angle
           the periastron is directly behind the companion thus maximising the effect.
8.4 3-Dimensional Parameter spaces                                                  78




Fig. 8.23: This 2-dimensional plot shows the combined effect of eccentricity and the
           position of the periastron on the Roemer delay for a standard compact bi-
           nary system. The symmetry noticed in the previous graph is apparent here
           again but in this case the gradient of the increase in delay is much greater
           and implies that at periastron positions of 0◦ and 180◦ there is a maximal
           difference between Shapiro and Roemer delay.
8.4 3-Dimensional Parameter spaces                                               79




Fig. 8.24: This 2-dimensional plot shows the combined effect of companion mass and
           semi-major axis length on both the Shapiro and Roemer delay for a stan-
           dard compact binary system. This figure illustrates the earlier statement
           that Shapiro delay is not dependent on orbit size and Roemer delay is not
           dependent on companion mass.
8.4 3-Dimensional Parameter spaces                                                      80




Fig. 8.25: This 2-dimensional plot shows the combined effect of inclination and the
           position of the periastron on the Shapiro delay for a standard compact binary
           system. There is a slight increasing effect due to the position of the periastron
           but this is clearly buried in the increase due to inclination angle
8.4 3-Dimensional Parameter spaces                                                     81




Fig. 8.26: This 2-dimensional plot shows the combined effect of inclination and the
           position of the periastron on the Roemer delay for a standard compact binary
           system. In this figure there is also a bulge in the Roemer delay at a periastron
           position of around 1.6 rad and a decrease at positions around 0 and 3.2 rad.
           This would indicate a difficulty in separating Shapiro and Roemer delays
           purely through systems with the periastron of the pulsar in the line of sight
           of the center of mass of the system.
                      9. CONCLUDING REMARKS


General Relativity is still, after 90 years of tests, the front-runner in the race to
understand gravity. If it has stood up so well to experiment, why continue to go
to such pains to develop new methods to test it? There are two reasons for doing
this: firstly gravity is one of the fundamental forces of nature and if we do not
understand the fundamental forces we cannot hope to understand the stranger
aspects of the universe. Secondly we are fairly confident that General Relativity
as it stands is not sufficient to allow us to unify the four fundamental forces of
nature into one Grand Unified Theory. Failures to quantise gravity have led us
to believe that to some degree it must be adjusted before unifying it with the
other forces. And thus we continue to test its framework and predictions in the
hope of finding, or not-finding, some discrepancy.
    General Relativity predicts a Shapiro time delay in an electromagnetic wave
propagating through the gravitational field of any object. The magnitude of this
delay is dependent on the region of the gravitational field that the electromag-
netic wave passes through. General Relativity takes the view that gravity is not
an exchange force but simply a geometry where all objects bend, or warp, the
spacetime fabric that makes up the universe. The greater the warp that the elec-
tromagnetic signal travels through, the greater the delay it experiences. It should
be immediately obvious that the heavier the body responsible for the bending of
spacetime or the closer that the signal passes to that body, the greater the delay
experienced by the signal. We have seen how this effect has been thoroughly
tested in solar systems tests from the Viking mission to Mars in the 70’s to the
Cassini mission to Saturn a few years ago. The PPN parameter testing the valid-
ity of various theories with respect to the Shapiro delay has been constrained to
within 0.012% of unity. Although this does not eliminate all competing theories
as can be seen from figure (2.2), it does indicate that in the Brans-Dicke the-
ory, for example, the parameter ωBD must be greater than 40000 to meet those
experimental results.
    These facts highlight the success of General Relativity in the weak-fields of
the solar system but more needs to be done in the strong-field regime found only
in the near vicinity of pulsars and black holes. It is for this reason that this
thesis involves the simulation of the Shapiro delay in compact binary systems
where one component is a pulsar and the other a neutron star or black hole. We
have seen that the Shapiro delay is occulted by the Roemer time of flight delay
of the signal as it traverses the finite size of the orbit. The greater the size of
                                                                                83


the orbit, measured through the semi-major axis and described by the parameter
a, the greater the size of the Roemer delay and the more difficult it becomes to
separate the very small Shapiro delay from it. Therefore it is essential for future
tests that we search for and time binary pulsar systems that are as compact as
possible with as heavy a companion as can be found. These two properties, if
found in the same system, will maximise the Shapiro delay while minimising the
Roemer delay and accordingly increase our chances of constraining γ.
    Simulations done here have shown, and it has been discussed in theory, that
high inclination systems are the most promising in which to measure the Shapiro
delay. High inclination means the pulsar will move in an almost eclipsing orbit
forcing its beam to pass at close to grazing incidence to the companion and
thereby maximise the delay. Furthermore it can be argued that only in these high
inclination systems does the beam pass sufficiently close to the companion so as to
pass through the strong field we wish to test. This may be true but there are many
factors to consider. For a binary system with a semi-major axis of 1 000 000 km
at an inclination angle of 89◦ , the impact parameter is approximately 20 000 km,
1000 times the diameter of the companion, assuming a neutron star companion.
At an inclination of 80◦ , which is still regarded as high inclination, the impact
parameter has increased by an order of magnitude and so drawing a line beyond
which the beam no longer experiences strong field is difficult. Additionally there
is the position of periastron to take into account. The maximum delay occurs
when the pulsar is behind the companion within the line of sight. Should the
semi-major axis line be at right angles to the line of sight, the pulsar will be
in that maximal position i.e. behind the companion, at a distance much closer
than the semi-major axis length that is used to estimate the impact parameter.
So should periastron lie in the line of sight, the impact parameter is increased
while if it lies at 90◦ to the line of sight, it will be minimised. Furthermore
the companion could be a black hole in which case the region of strong field
gravity is greatly enlarged. Clearly many considerations need to be taken into
account when determining what is meant by strong field. For the remainder of
this discussion we will assume that the entire pulsar binary system falls within
the category of strong field.
    Varying the eccentricity of the orbital system indicates that this could be a
very useful parameter to consider when testing for the Shapiro delay. In the
absence of a sufficiently high inclination, the Shapiro delay cannot be separated
from the Roemer delay as the two signals are in phase. As the eccentricity
increases, however, there is an increase in the phase shift between the two signals
with the result that subtracting out the Roemer delay should not eliminate the
Shapiro delay along with it. Additionally there is the advantage that, as the
inclination of the orbit decreases, so too does the Roemer delay. Unfortunately
there are some problems with this line of investigation. The very nature of
binary pulsar searches makes finding high eccentricity orbits extremely difficult.
The problem arises with the searching of binary systems in the 5-dimensional
                                                                                84


parameter space of parameters affecting binary orbits. Because of the difficulty
in this, one method of searching is to assume low eccentricity and longitude,
or position, of periastron resulting in a 3-dimensional parameter space that is
easier to search. A simpler method frequently used is known as an acceleration
search which assumes a constant acceleration over each orbital period. Obviously
this approximation will deviate further from reality as eccentricity increases and
so this indicates that most binary pulsar searches are only looking for close to
circular orbits. New methods of pulsar searching could assist in finding more high
eccentricity systems that could be used for the testing of our earlier supposition.
    Although this thesis has not attempted to constrain γ to any degree, it has
shown how the measurement of arrival times of pulses emitted by a pulsar in a
binary system can be used to constrain γ through determining the Shapiro delay.
We have seen how the parameter appears in the governing equations and how it
could be determined using pulsar timing data. By timing a greater number of
systems and obtaining sufficiently large data sets, it will be possible to obtain
Shapiro delay data for a greater number of pulsar binary systems and to fit the
data so as to place limits on the PPN parameter. By simulating the effect of
a black hole companion we have preempted the hoped for discovery of such a
system and shown the expected results that such a system would offer. The
measurement of the Shapiro delay of light in a binary pulsar system is but one
test in the quest to understand the universe but every step we take is one closer
to achieving that most noble of goals.
                            10. REFERENCES


• Baade, W.,Zwicky, F., “Cosmic Rays from Super-novae”, Proc. Nat. Acad.
Sci., 20, 254, (1934).
• Backer, D., Hellings, R., “Pulsar Timing and General Relativity”, Ann. Rev.
Astron. Astrophys., 24, 537-575, (1986).
• Blandford, R., Teukolsky, S.A., “Arrival-time analysis for a pulsar in a binary
system”,Ap. J., 205, 580-591, (1976).
• Burgay, M., D’Amico, N., Possenti, A., et al.,“The Highly Relativistic Binary
Pulsar PSR J0737-3039: Discovery and Implications”, ASP Conf. series, 328,
53, (2005).
• Camilo, F., Foster, R.S., Wolszczan, A., “High Precision Timing of J1713+0747:
Shapiro Delay”, Ap. J., 437, L39-L42, (1994)
• Cordes, J.M., “NE2001: A New Model for the Galactic Electron Density and it’s
Fluctuations” in Milky Way Surveys: The Structure and Evolution of our Galaxy,
Proceedings of ASP Conference 317, Eds, Clemens, D., Shah, R., Brainerd, T.,
p211, (2004).
• Damour, T., Gibbons, G.W., Taylor, J.H., “Limits on the Variability of G
Using Binary Pulsar Data”, Phys. Rev. Lett., 61, 1151-1154, (1988).
• Eddington, A.S., “The Total Eclipse of 1919 May 29 and the Influence of
Gravity on Light”, The Observatory, 42, 119-122, (1919).
• Einstein, A., “Die Grundlage der Algemeinen Relativit¨tstheorie”, Ann. der
                                                             a
Physik, 49, 769-822, (1916).
• Gold, T.,“Rotating Neutron Stars as the Origin of the Pulsating Source of
Radio”, Nature, 218,731 , (1968).
• Goldman, I., “Upper Limit on G Variability Derived from the Spin-Down of
PSR 0655+64”,Mon. Not. R. Astron. Soc., 244,184-187, (1990).
• Goldstein, H., Poole, C., Safko, J., Classical Mechanics 3rd edition, Addison
Wesley, San Francisco, (2002).
• Hartle, J.B., Gravity: An Introduction to Einstein’s General Relativity, Addison
Wesley, San Francisco, (2003).
• Haugan, M., “Post-Newtonian Arrival-Time Analysis for a Pulsar in a Binary
System”, Ap. J., 296, 1-12, (1985).
• Hewish, A., Bell, J., et al., “Observations of a Rapidly Pulsating Radio Source”,
Nature, 217, 709, (1968).
• Hulse, R., Taylor, J.,“Discovery of a Pulsar in a Binary System”, Ap. J., 195,
L51-L53, (1975).
                                                                               86


• Large, M.I., Vaughan, A.E., Mills, B.Y.,“A Pulsar Supernova Association”,
Nature, 220, 340-341, (1968).
• Lebach, D. E., Corey, B. E., Shapiro, I. I., et al., “Measurement of the solar
gravitational deflection of radio waves using very-long-baseline interferometry”,
Phys. Rev. Lett., 75, 14391442, (1995).
• Lorimer, D., Kramer, M., Handbook of Pulsar Astronomy, Cambridge Univer-
sity Press, Cambridge, (2005).
• Lyne, A.G., Burgay, M., Kramer, M., et al.,“A Double-Pulsar System: A Rare
Laboratory for Relativistic Gravity and Plasma Physics”, Science, 303, 1153,
(2004)
• Lyne, A., Graham-Smith, F., Pulsar Astronomy 3rd edition, Cambridge Uni-
versity Press, Cambridge, (2005).
• Nordtvedt, K., “Testing Relativity with Laser Ranging to the Moon”, Phys.
Rev., 170. 1186-1187. (1968).
• Nordtvedt, K. Jr.,“Post Newtonian Metric for a General Class of Scalar-
Tensor Gravitational Theories and Observational Consequences”,Ap. J., 161,
1059, (1970).
                    ˙
• Nordtvedt, K., “G/G and a Cosmological Acceleration of Gravitationally Com-
pact Bodies”,Phys. Rev. Lett., 65, 953-956, (1990).
• Pacini, F., “Energy Emission from a Neutron Star”, Nature, 216,567 , (1967).
                        ¨
• Schwarzschild, A., “Uber das Gravitationsfeld eines Massanpunktes nach der
Einsteinschen Theorie”, Proc. Royal Prus. Ac. Sci., 1, 189-196, (1916).
• Shapiro, I.I., “The Fourth Test of General Relativity”,Phys. Rev. Lett., 13,
798, (1964).
• Shapiro, I.I., Reasenberg, R.D., et. al., “The Viking Relativity Experiment”,J.
Geophys. Res., 82, 4329-4334, (1977).
• Shapiro, S. S., Davis, J. L., Lebach, D. E., and Gregory, J. S., “Measurement of
the solar gravitational deflection of radio waves using geodetic very-long-baseline
interferometry data, 1979-1999”, Phys. Rev. Lett., 92, 121101, (2004).
• Stairs, I.H., “Pulsars in Binary Systems: Probing Binary Stellar Evolution and
General Relativity”, Science, 304, 547, (2004).
• Stairs, I.H., “Testing General Relativity with Pulsar Timing”, Liv. Rev. Rel.,
6, 5, (2005).
• Thorne, K. S., Will, C.M., “Theoretical Framework for Testing Relativistic
Gravity I: Foundations”,Ap. J., 163, 595, (1971).
• Thorsett, S.E., “The Gravitational Constant, the Chandrasekhar Limit, and
Neutron Star Masses”,Phys. Rev. Lett., 77, 1432-1435, (1996).
• van Straten, W., et al., “A test of General Relativity from the Three-
Dimensional Orbital Geometry of a Binary Pulsar”, Nature, 412, 158-160, (2001).
• Weinberg, S., Gravitation and Cosmology: Principles and Applications of the
General Theory of Relativity, John Wiley and Sons, (1972).
• Wex, N., “New limits on the violation of the Strong Equivalence Principle in
strong field regimes”. Astron. Astrophys., 317, 976-980. (1997).
                                                                              87


• Wheeler, J.C., “Timing effects in Pulsed Binary Systems”, Ap. J. (Let-
ters), 196, L67, (1975).
• Will, C.M.,“Theoretical Framework for Testing Relativistic Gravity II:
Parametrized Post-Newtonian Hydrodynamics and the Nordvedt Effect”,Ap.
J., 163, 611, (1971).
• Will, C.M.,“Theoretical Framework for Testing Relativistic Gravity III: Conser-
vation Laws, Lorentz Invariance and Values of the PPN parameters”,Ap. J., 169,
125, (1971).
• Will, C.M., Nordtvedt, K. Jr.,“Conservation Laws and Preferred Frames in
Relativistic Gravity I: Preferred-Frame Theories and an Extended PPN Formal-
ism”,Ap. J., 177, 757, (1972).
• Will, C.M., Nordtvedt, K. Jr.,“Conservation Laws and Preferred Frames in Rel-
ativistic Gravity II: Experimental Evidence to Rule Out Preferred-Frame Theo-
ries of Gravity”,Ap. J., 177, 775, (1972).
• Will, C.M., “Is momentum conserved?           A test in the binary system
PSR 1913+16”,Astrophys. J., 393. L59-L61. (1992).
• Will, C.M., “The Confrontation Between General Relativity and Experi-
ment”,Liv. Rev. Rel., submitted (astro-ph:gr-qc/0510072),(2005).
APPENDIX
  A. THE CONNECTION OR COVARIANT DERIVATIVE


This is a mathematical function fundamental to GR with some very interesting
and special properties. The quantity X Y is a function that, on a differentiable
manifold M , takes in two vector field inputs and outputs a third vector field.
Practically what it does is measure the rate of change of the Y-vector field in the
direction of the X-vector field. Mathematically this is written as:

                               : χ(M ) × χ(M ) → χ(M )
                               : (X, Y ) −→ X Y

    As previously mentioned,      has some very important properties that have
deep consequences on the subject matter to be dealt with. In order to highlight
these properties we shall look at the function in terms of its subscript and main
arguments. It turns out that is linear in its subscript argument but not in the
second argument and as such is not a tensor. To see this one must first note that
  X Y can be re-written as

                                        δ i
                          XY   =Yj         Y Ei = (LX )Y i Ei               (A.1)
                                       δxj
where LY is the differential operator acting on the vector field Y such that Ly =
     δ
Y i δyi and Ei is the basis of the space.
     Using this we can now show that is linear in its subscript argument. So

                       X+X   Y = (LX+X Y i )Ei
                               = (LX Y i + LX Y i )Ei
                               = (LX Y i )Ei + (LX Y i )Ei
                               =   XY +      X Y;                           (A.2)

and

                                f XY    =   (Lf X Y i )Ei
                                        =   (fLX Y i )Ei
                                        =   f (LX Y i )Ei
                                        =   f X Y.                          (A.3)

   In the second argument, however, it is not linear
A.1 Torsion                                                                     90




                        X (Y   + Y ) = (LX [Y i + Y i )Ei
                                     = (LX Y i + LX Y i )Ei
                                     = (LX Y i )Ei + (LX Y i )Ei
                                     =   XY +      XY ;                       (A.4)

and

                       X (f   Y) =    (LX [fY i ])Ei
                                 =    ([LX f ]Y i + f [LX Yi ])Ei
                                 =    [LX f ](Y i Ei ) + f ([LX Y i ]Ei )
                                 =    (LX Y) + f X Y.                         (A.5)

This shows that, although the first criterion for linearity is met, the second is
clearly a Leibniz multiplication law and not linear. For to be a tensor it must
be linear in all its arguments. Rather, is known as a connection.


                                     A.1 Torsion
Is there a difference between        XY   and    Y X?

                                              i
                              XY    =   X (Y Ei )
                                            i        i
                                    =   X Y (Ei ) + Y (       X Ei )
                                    = (LX Y i )Ei                             (A.6)

since   X Ei   = 0 because the basis does not change. Similarly we find that

                                     YX   = (LY X i )Ei ,                     (A.7)

and therefore the difference is

                          XY    −    YX   = (LX Y i + LY X i )Ei .            (A.8)

This is simply the commutator bracket [X, Y]. This property of the connection is
known as symmetry. A connection is said to be symmetric if, on a differentiable
manifold,
                             XY −     Y X = [X, Y]                        (A.9)
and thus, in a Euclidean space

                                XY   −    YX   − [X, Y] = 0                 (A.10)
A.2 Curvature                                                                    91


The next step is to generalise this statement so that it applies to any case and
not just to an Euclidean space. The general format of this equation introduces
us to the concept of Torsion and the Torsion tensor T . In general

                         XY   −   YX   − [X, Y] = T (X, Y)                    (A.11)



In order to understand the significance of this it is necessary to under-
stand the significance and geometrical interpretation of the commutator bracket
[X, Y].
   Consider the field lines X and Y parameterised by σ and λ respectively as
shown in figure (3.1). We then pick a point (xo , yo ) and follow the X field line
that passes through that point by a set parameter value σ to the point P1 . Next
one moves along the Y field line passing through that point by a set parameter
value λ to the point P2 . If we then move the same parameter ’distance’, λ along
the Y field line passing through the point (xo , yo ) to a point P3 and then the set
parameter value σ along the X field line that passes through the point P3 to a
point P4 , the vector from point P2 to point P4 is the commutator [X, Y].




        Fig. A.1: Geometric interpretation of the commutator bracket [X, Y]



                               A.2 Curvature
In many cases one would be expected to use the connection twice in order to
determine the second covariant derivative of a vector field. The order of this
differentiation is very important to the outcome of the problem; by analysing the
difference between the orders in which differentiation can be applied, valuable
insight can be obtained. In mathematical notation:
A.2 Curvature                                                                       92


    If X, Y and Z are vector fields on a manifold M, then both X Z and Y Z
will be vector fields on M. We are interested in the second derivatives and so
consider first X Y Z. We know from our analysis of torsion that
                                                       i
                                        YZ    =   Y (Z Ei )
                                              = ( Y Z i )Ei
                                              = (LY Z i )Ei ;                   (A.12)

Therefore:
                                                          i
                                X       YZ   =   X ((LY Z )Ei )
                                             = ( X (LY Z i ))Ei
                                             = (LX (LY Z i ))Ei .               (A.13)

Similarly:
                                Y       XZ   = (LY (LX Z i ))Ei .               (A.14)
If we construct the difference between these two terms we have

                   X   YZ   −       Y    XZ       = (LX LY − LY LX )Z i )Ei
                                                  = (L[X,Y] Z i )Ei
                                                  = ( [X,Y] Z i )Ei
                                                               i
                                                  =   [X,Y] (Z Ei )
                                                  =   [X,Y] Z.                  (A.15)

   Therefore it can be seen that, for a Euclidean space:

                            X   YZ      −     Y    XZ   −   [X,Y] Z   =0        (A.16)

    A space in which this last statement holds is said to be flat or to have zero
curvature. If we now, as before, generalise this statement to all spaces, we find
that the above statement is usually not equal to zero but is equal to a rank 3
tensor that takes 3 vector field inputs and gives one vector field output. This
is because there is nothing in the definition of     to guarantee that this tensor
must be zero. This tensor is known as the Riemannian curvature tensor, R and
is defined by
                           ˙
                  R(X, Y)Z = X Y Z − Y X Z − [X,Y] Z                       (A.17)
The geometrical interpretation of the [X,Y] Z illustrates its meaning far more
clearly.
    First we consider a point (xo , yo ) with field lines X and Y going through that
point and parameterised by λ and σ respectively. At the point (xo , yo ) we have
                  ¯ ¯
tangent vectors (xo , y o ) to the field lines. We then parallel transport these tangent
vectors along the field lines to the points (xo +Xλ, yo ) and (xo , yo +Y σ) to obtain
A.2 Curvature                                                                     93




                   Fig. A.2: Geometric interpretation of             [X,Y] Z


          ¯ ¯
vectors (x1 , y 1 ). At these new points there are yet again tangent vectors to the
                                      ¯ ¯
field lines and these we shall call (x2 , y 2 ). In this scheme

                             [X,Y] Z
                                          ¯ ¯          ¯ ¯
                                       = (x2 , x1 ) − (y 2 , y 1 )             (A.18)

which equals zero in the case of a coordinate system because the commutator
[X, Y] equals zero in a coordinate system.
B. MATLAB CODE