format: text

text:
# Energy transport, overshoot, and mixing in the atmospheres of very cool stars

## 1. Introduction

The increasing number of late M-type stars, brown dwarfs, and extrasolar planets discovered by infrared surveys and radial velocity searches has spawned a great deal of interest in the atmospheric physics of these objects. Their atmospheres are substantially cooler than e.g. the solar atmosphere, allowing the formation of molecules, or even liquid and solid condensates. Convection is a ubiquitous process in these atmospheres shaping the thermal structure and the distribution of material. Hydrodynamical simulations of solar and stellar granulation including a realistic description of radiative transfer have become an increasingly powerful and handy instrument for studying the influence of convective flows on the the structure of late-type stellar atmospheres as well as on the formation of their spectra. Here we report on efforts to construct hydrodynamical models for mid M-type atmospheres. We consider this as an intermediate step on the way to model brown dwarf and planetary atmospheres. The motivation was twofold: Firstly, pre-main-sequence evolutionary models of M-dwarfs and brown dwarfs based on mixing-length theory (MLT) depend sensitively on the mixing-length parameter $`\alpha _{\mathrm{MLT}}`$ (Baraffe et al. 2002). Secondly, the distribution of dust clouds in brown dwarfs depends on the efficiency of mixing by convective overshoot. Conventionally, convection is described within the simplistic picture of MLT dependent on free parameters. Our hydrodynamical models provide a description essentially from first principles free of the uncertainties of MLT, putting stellar models on a firmer footing.

## 2. Model overview

Figure 1 illustrates in the location of the hydrodynamical model atmospheres considered in this investigation in the $`\mathrm{T}_{\mathrm{eff}}`$-$`\mathrm{log}\mathrm{g}`$-plane: three M-type models are located at $`\mathrm{T}_{\mathrm{eff}}`$=2800 K with $`\mathrm{log}\mathrm{g}`$=3.0, 4.0, and 5.0 which form a $`\mathrm{log}\mathrm{g}`$-sequence. To assess the temperature effects a somewhat hotter M-type model is located at $`\mathrm{T}_{\mathrm{eff}}`$=3280 K and $`\mathrm{log}\mathrm{g}`$=4.0. For comparison we also considered a solar model at $`\mathrm{T}_{\mathrm{eff}}`$=5777 K, $`\mathrm{log}\mathrm{g}`$=4.44. All models possess solar chemical composition.

The RHD simulations were performed with a convection code developed by Å. Nordlund and. R.F. Stein (see Stein & Nordlund 1998, and references therein). The code implements a consistent treatment of compressible gas flows together with non-local radiative energy exchange in three spatial dimensions. The radiative transfer is treated in LTE approximation, the wavelength dependence of the radiation field is represented by a small number of wavelength bins. Open lower and upper boundaries, as well as periodic lateral boundaries are assumed. The effective temperature of a model (i.e. the average emergent radiative flux) is controlled indirectly by prescribing the entropy of inflowing material at the lower boundary. Magnetic fields are neglected. Opacities and equation-of-state have been adapted to the conditions encountered in the M-type atmospheres. In particular, the equation-of-state accounts for $`\mathrm{H}_2`$ molecule formation, the opacities include contributions of molecular lines. The opacities were extracted from the opacity data base of the PHOENIX model atmosphere code (for a description of the code see Hauschildt & Baron 1999). The application of the Nordlund & Stein convection code to M-type objects in described in more detail in Ludwig, Allard, & Hauschildt (2002, hereafter LAH).

Figure 2 visualises the flow structure (here via the temperature field) generically encountered in our hydrodynamical models for late-type stars. The computational domain is a Cartesian box representing a small volume at the stellar surface. It comprises the optically thin photospheric layers as well the uppermost layers of the convection zone which extends deeper into the star. A granulation pattern is visible in layers around optical depth unity. Concentrated, plume-like downdrafts dominate the sub-photospheric flow, the optically thin layers are harbouring a mixture of overshooting and wave motions with associated temperature fluctuations.

We want to derive quantitative estimates of the mixing by convective overshoot, as well as obtain a measure of the efficiency of the convective energy transport. For addressing these issues the RHD models have to give a reasonably accurate representation of the actual atmospheric conditions. Here we are particularly concerned about the radiative energy transport, which is complicated by the huge number of molecular absorption lines. In our hydrodynamical models we use a multigroup technique (dubbed Opacity Binning Method, hereafter OBM) for modelling the radiative energy exchange which employs four groups for representing the wavelength dependence of the radiation field (see LAH for details). The wavelength groups have been optimised for the $`\mathrm{T}_{\mathrm{eff}}`$=2800 K and $`\mathrm{log}\mathrm{g}`$=5.0 model. Figure 3 illustrates the accuracy which is achieved with the OBM. While there are differences between the atmospheric structure based on the OBM and PHOENIX models, which are based on direct opacity sampling, the OBM provides a significant improvement with respect to the simple grey approximation. Differences increase as one chooses atmospheric parameters away from the ones for which the OBM was optimised. We did not generate specifically optimised groupings for all atmospheres under consideration, since we were interested to study the systematic change of model properties with effective temperature and gravitational acceleration. I.e it appeared advantageous to leave the input physics the same in all models, only optimised for one temperature and gravity.

To mitigate the effects of the still present shortcomings in the radiative transfer we took a differential approach when measuring model properties: whenever possible we compared hydrodynamical and hydrostatic model atmospheres based on the same radiative transfer scheme, i.e. here the OBM.

## 3. Granulation in M-type objects

Our small sample of hydrodynamical M-type atmospheres gives some insight into the properties of stellar granulation at cooler temperatures. Figure 4 compares the granulation pattern of the $`\mathrm{T}_{\mathrm{eff}}`$=2800 K, $`\mathrm{log}\mathrm{g}`$=3.0 model with the one of a solar model. The first thing to note is that surface convection in M-type objects produces a granular pattern qualitatively resembling solar granulation: bright extended regions of upwelling material which are surrounded by dark concentrated lanes of downflowing material. The dark lanes form an interconnected network. Looking more closely, granules are less regularly delineated in M-type objects, the inter-granular lanes show a higher degree of variability in terms of their strength. A feature which is uncommon in the solar granulation pattern are the dark “knots” found in or attached to the inter-granular lanes. The knots are associated with strong downdrafts which carry a significant vertical component of angular momentum. The width of the inter-granular lanes to the typical granular size is smaller in our M-type objects. Inspecting the velocity field (not shown) in vicinity of the continuum forming layers shows less pronounced size differences. This indicates that the relatively broader lanes in the solar case are the result of a stronger smoothing of the temperature field due to a more intense radiative energy exchange, i.e. the effective Peclét number of the flow is larger around optical depth unity in the cooler objects.

## 4. Mixing by atmospheric overshoot

Stellar atmospheres around $`\mathrm{T}_{\mathrm{eff}}2800\mathrm{K}`$ are too hot to allow for a significant formation of dust grains. However, already at slightly cooler effective temperatures, grain formation sets in and dust grains become major opacity contributors, i.e. an important factor in determining the thermal structure of the atmosphere. In fact, the spectral energy distribution in the range of effective temperatures $`2500\mathrm{K}T1500\mathrm{K}`$ is crucially linked to the distribution of dust grains in the atmosphere (e.g. Allard, Hauschildt, & Schwenke 2000). The amount of dust which is present is determined by chemical condensation and evaporation processes as sources and sinks, as well as macroscopic transport processes which carry dust grains away from their sites of formation. In M-type and cooler atmospheres the transport is dominated by two opposing processes: gravitational settling of dust grains and their mixing due to the presence of velocity fields, either related to convection or global circulations induced by rotation. In our models, no dust formation takes place but we nevertheless find it worthwhile to give a characterisation of the atmospheric mixing due to convection and convective overshoot. This can give at least a first order approximation of how convective mixing might operate when dust is actually present.

Formally, we are interested in a statistical representation of the mean transport properties of the convective velocity field in vertical direction. The horizontal advection of dust grains by the convective velocity field probably produces horizontal inhomogeneities in the dust distribution. We neglect these here since i) we are targeting at the application of our results in standard 1D stellar atmosphere models, ii) the horizontal inhomogeneities are small scale (the size of a convective cell), i.e. hardly observable and iii) the uncertainties in our understanding of the dust formation process itself perhaps limits the achievable accuracy anyway. In view of the last point we present a proxy of the convective mixing only, without trying to derive a detailed statistical description of the transport properties of the convective velocity field, i.e. extracting its effective transport coefficients (c.f. Miesch, Brandenburg, & Zweibel 2000, for an example of a more stringent treatment).

As proxy measuring the mixing properties of the velocity field we introduce the mixing velocity $`v_{\mathrm{mix}}`$

$$v_{\mathrm{mix}}\frac{F}{\rho }$$
(1)
($`\rho `$ denotes the mass density, $`.`$ temporal and horizontal averaging) where $`F`$ is the temporally and horizontally averaged vertical gross mass flux

$$F\rho u_\mathrm{z}_{\mathrm{t},\mathrm{up}}$$
(2)
($`u_\mathrm{z}`$ denotes the vertical velocity). For $`F`$ the average is taken over surface areas with upwards directed flow only. Note, that the average total mass flux vanishes. $`F`$ measures the amount of mass flowing in and out a certain layer. If the in- and outflow happens sufficiently disordered, this mass exchange goes hand in hand with mixing of material.

As we shall see below $`v_{\mathrm{mix}}`$ exhibits an exponential height dependence. With such an dependence a characteristic mixing frequency $`f_{\mathrm{mix}}`$ can be written as

$$f_{\mathrm{mix}}\frac{1}{\rho }\frac{dF}{dz}=v_{\mathrm{mix}}\left(\frac{1}{H_\mathrm{v}}+\frac{1}{H_\rho }\right),$$
(3)
where $`H_\mathrm{v}`$ denotes the scale height of $`v_{\mathrm{mix}}`$, and $`H_\rho `$ the density scale height.

### 4.1. Subsonic filtering

As mentioned previously, the atmospheric velocity field is a superposition of advective, overshooting motions and acoustic waves generated by convection in deeper layers (see also Ludwig & Nordlund 2000). The wave motions contribute little if at all to the mixing due to their spatially coherent, oscillatory character. The overshooting, convective motions tend to decay with distance from the Schwarzschild stability boundary, while the wave motions tend to increase in amplitude due to the exponential decrease of the density with increasing height. This leads to the situation that beyond a certain height the atmospheric velocity field is even dominated by wave motions. In order to get a reliable estimate of the mixing it is therefore necessary to remove the wave contributions to the velocity field before evaluating the gross mass flux $`F`$.

We removed the wave contributions by subsonic filtering — a technique invented in the context of solar observations for cleaning images from “noise” stemming from the solar 5 minute oscillations (Title et al. 1989). Figure 5 schematically illustrates this filtering technique. In short, one considers a time sequence of images and removes features with horizontal phase speeds $`v_{\mathrm{phase}}`$ greater than a prescribed threshold. This is achieved by Fourier filtering of spatial-temporal data in the $`k`$-$`\omega `$ domain. For every depth layer in our data cubes we performed a 3D Fourier analysis (one temporal, two spatial dimensions) of the vertical mass flux retaining only contributions below a preset velocity threshold. In practice, acoustic and convective contributions are not as cleanly separated as shown in the Fig. 5, and one must find the right balance between removing as much acoustic components as possible while retaining as much as possible convective contributions. We always studied a sequence of phase speed thresholds in order to judge the success of the procedure.

### 4.2. Mixing in the solar atmosphere

Despite we are not primarily interested in the Sun here, it is worthwhile to look at the Sun as a reference object to compare with. Figure 6 shows $`v_{\mathrm{mix}}`$ in our solar model for various degrees of subsonic filtering. It is clearly visible that the subsonic filtering has the strongest impact on $`v_{\mathrm{mix}}`$ in the uppermost atmospheric layers. As advertised before $`v_{\mathrm{mix}}`$ exhibits an exponential decline with height ($`\mathrm{log}Pz`$) after appropriate subsonic filtering. This feature was put forward by Freytag, Ludwig, & Steffen (1996) as generic feature of convective overshoot. Moreover, we see that too low a velocity threshold removes also convective features, as visible by the reduction of the velocity in the deeper, convection dominated layers. Quantitatively, for the Sun we find that $`H_\mathrm{v}=(1.9\pm 0.2)H_\mathrm{p}`$, where $`H_\mathrm{p}`$ denotes the local pressure scale height.

The exponential decline of $`v_{\mathrm{mix}}`$ is also motivated from studying linear convective modes. In Fig. 6 we plotted the velocity profiles of to linear convective eigenmodes with horizontal wavelength of 1.5 and 5.0 Mm, respectively. We used the temporally and horizontally averaged hydrodynamical model as background mode. The absolute velocity amplitude of the modes has been scaled to match $`v_{\mathrm{mix}}`$ from the hydrodynamical model. We see that the velocity field of the modes exhibit an exponential “leaking” into the formally convective stable layers, and the decline is more rapid for the mode with shorter wavelength. We further see that the mode with 5 Mm wavelength shows a decline similar to the decline of $`v_{\mathrm{mix}}`$. Interestingly, the wavelength of the fitting mode is significantly larger than the horizontal scale of the dominant convective structures on the Sun — the granules with typical sizes of around 1.5 Mm.

### 4.3. Mixing in M-type atmospheres

Figure 7 shows $`v_{\mathrm{mix}}`$ for the M-type models at $`\mathrm{T}_{\mathrm{eff}}`$=2800 K. With decreasing gravity the zone of convective instability extends further and further into the optically thin atmosphere, leaving little room for overshoot in the lowest gravity model. Reading off an exponential decline rate is very uncertain here. However, we find a slow decline with $`H_\mathrm{v}6H_\mathrm{p}`$. For the $`\mathrm{log}\mathrm{g}`$=4.0 model we find a value close to solar of $`1.7H_\mathrm{p}`$, for the $`\mathrm{log}\mathrm{g}`$=5.0 model $`0.5H_\mathrm{p}`$. Qualitatively, overshooting is less pronounced in models of higher gravity. This is of course related to the fact that the buoyancy forces scale proportional to gravity, making buoyancy more effective in confining the convective motions to the formally unstable regions. Empirically we find that the scale height of $`v_{\mathrm{mix}}`$ roughly scales as $`g^{\frac{1}{2}}`$.

Trying to match the velocity profile in the overshooting regions worked only partially. For the lower gravity models the fit is not very good which might be related the situation that convection reaches very heigh atmospheric layers. We oriented the horizontal wavelength of the linear modes at the largest sizes of structures the computational box could accommodate in the respective models.

## 5. The efficiency of the convective energy transport

Convection is an important energy transport mechanism in M-type stars. In standard model atmospheres it is treated in the framework of MLT. The question is whether the simplistic MLT is actually capable to provide a reasonable description of the convective energy transport under conditions encountered in the M-type atmospheres. Figure 8 shows a comparison of the thermal structure of the hydrodynamical model atmospheres and standard 1D hydrostatic model atmospheres in radiative-convective equilibrium assuming different mixing-length parameters. The hydrodynamical models have been averaged temporally and horizontally on surfaces of constant optical depth. This procedure ensures a particularly good preservation of the energy transport properties of the hydrodynamical models (Steffen, Ludwig, & Freytag 1995). Besides the mixing-length parameter MLT contains a number of further “hidden” parameters intrinsic to the specific formulation of MLT which was chosen. We emphasise that a well-defined calibration of the mixing-length parameter must be always given with reference to the specific formulation. Here, we are using the formulation as given by Mihalas (1978).

The sensitivity of the structure of the standard models to the mixing-length parameter increases with decreasing gravity. This is in line with the earlier statement that stellar structure models depend sensitively on the mixing-length parameter at lower gravities. Towards the main sequence the sensitivity is largely reduced making a precise calibration of $`\alpha _{\mathrm{MLT}}`$ difficult — but of course also less important. The M-type atmospheres offer the remarkable opportunity to study convection also under optically thin conditions. We roughly separate the models into a optically thin and thick part and derive the corresponding mixing-length parameters suitable to fit the hydrodynamical structures in the different regions. The mixing-length parameter of the optically thick region is most relevant for global stellar structure models, while the one of the optically thin part is most relevant for stellar atmospheres.

The entropy gradients of hydrodynamical and hydrostatic models correspond closely in the highest, convectively stable layers (with $`\frac{ds}{d\tau }<0`$). By construction the standard models are in radiative equilibrium in these layers. This indicates that in the hydrodynamical models overshooting and wave heating are not particularly efficient, i.e. in the hydrodynamical models these layers are also close to radiative equilibrium. In the convectively unstable layers (with $`\frac{ds}{d\tau }>0`$) the situation is markedly different. The stratifications (except for the $`\mathrm{log}\mathrm{g}`$=5.0 model which is close to convective equilibrium, i.e. almost adiabatically stratified) are neither close to radiative nor convective equilibrium. The detailed balance between radiative and convective energy transport determines the thermal structure. This makes the stratification sensitive to the values of the mixing-length parameter which provides a measure of the efficiency of the convective energy transport. As indicated in Fig. 8 we find a mixing-length parameter between 2.1 and 2.8 in the optically thin regions, and between 1.5 and 1.9 in the optically thick region. The extreme values of 1.5 and 2.1 are found in the $`\mathrm{log}\mathrm{g}`$=5.0 model and are somewhat uncertain due to low sensitivity to changes of the mixing-length parameter. The $`\mathrm{T}_{\mathrm{eff}}`$=3280 K and $`\mathrm{log}\mathrm{g}`$=4.0 model (not plotted) shows values of 2.8 (optically thin part) and 2.2 (optically thick part).

Besides the detailed quantitative results two aspects are apparent: Firstly, it is not possible the specify a unique mixing-length parameter fitting the hydrodynamical stratification over the whole depth range. This illustrates principal shortcomings in the MLT. Secondly, the derived mixing-length parameters are falling in the range commonly considered. If one does not have high demands on the accuracy MLT provides a reasonable scaling of the energy transport efficiency with a fixed — perhaps solar calibrated — mixing-length parameter in the regime of M-type atmospheres.

## 6. Final remarks

While our results do not apply to the brown dwarf regime directly, it is tempting to extrapolate the mixing properties we found in our mid M-type objects. Our coolest and highest gravity model points towards a modest extra mixing by convective overshoot in brown dwarf atmospheres. The mixing would stay confined to layers in vicinity of the Schwarzschild stability boundary, and would not lead to a complete distribution of dust in the atmosphere (see LAH for more details).

Mixing-length theory does perform reasonably well in M-type atmospheres. To get a better quantitative description one might try to calibrate beside $`\alpha _{\mathrm{MLT}}`$ the “internal” parameters of MLT. M-type atmospheres appear particularly well suited for this undertaking since convection takes place under optically thick and thin conditions.

We reiterate that our results apply to atmospheres of solar metallicity. We expect a markedly different outcome for metal poor atmospheres.

#### Acknowledgments.

This project benefitted from financial support of the Walter Gyllenberg Foundation, and the Swedish Vetenskapsrådet.

## References

Allard, F., Hauschildt, P.H., & Schwenke, D., 2000, ApJ, 540, 1005

Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P.H., 2002, A&A, 382, 563

Freytag, B., Ludwig, H.-G., & Steffen, M., 1996, A&A, 313, 397

Hauschildt, P.H., & Baron, E., Journal of Computational and Applied Mathematics, 102, 41

Ludwig, H.-G., & Nordlund, Å., 2000: in: Stellar Astrophysics, Proceedings of the Pacific Rim Conference 1999, eds. K.S. Cheng, H.F. Chau, K.L. Chan, & K.C.Leung (Kluwer Academic Publishers), 37

Ludwig, H.-G., Allard, F., & Hauschildt, P.H., 2002, A&A, in press

Miesch, M.S., Brandenburg, A., & Zweibel, E.G., 2000, Phys. Rev. E, 61, 457

Mihalas, D., 1978, Stellar Atmospheres (Freeman and Company)

Steffen, M., Ludwig, H.-G., & Freytag, B., 1995, A&A, 300, 473

Stein, R.F., & Nordlund, Å., 1998, ApJ, 499, 914

Title, A.M., Tarbell, T.D., Topka, K.P., Ferguson, S.H., & Shine, R.A., 1989, ApJ, 336, 475


text:
# Monopoles and Solitons

## 1 Introduction

Non linear wave equations and solitons have been studied for a long time. Let us consider some less well known application of solitons. We start with monopoles. Ever since Dirac deduced theoretically the existence of the monopole in 1931, it has eluded physicists . At the same time the possibility of realising huge amounts of energy using monopoles has been an exciting prospect. In 1980 when the fiftieth Anniversary of the monopole was being commemorated, Dirac himself expressed his belief that the monopole did not exist . Some scholars have indeed dismissed the monopole , while in a model based on quantized vortices in the hydrodynamical formulation, the monopole field can be mathematically identified with the momentum vector . Monopoles had also been identified with solitons .
In any case, it has been noted that the existence of free monopoles would lead to an unacceptably high density of the universe , which in the light of latest observations of an ever expanding universe  would be difficult to reconcile.
We will now show that monopoles arise due to the non commutative structure of space time being ignored, and this would also provide an explanation for their being undetected.

## 2 The Monopole

Let us start by reviewing Dirac’s original derivation of the Monopole (Cf.ref.). He started with the wave function

$$\psi =Ae^{ı\gamma },$$
(1)
He then considered the case where the phase $`\gamma `$ in (1) is non integrable. In this case (1) can be rewritten as

$$\psi =\psi _1e^{ıS},$$
(2)
where $`\psi _1`$ is an ordinary wave function with integrable phase, and further, while the phase $`S`$ does not have a definite value at each point, its four gradient viz.,

$$\kappa ^\mu =^\mu S$$
(3)
is well defined. We use natural units, $`\mathrm{}=c=1`$. Dirac then goes on to identify $`\kappa `$ in (3) (except for the numerical factor $`hc/e`$) with the electromagnetic field potential, as in the Weyl gauge invariant theory.
Next Dirac considered the case of a nodal singularity, which is closely related to what was later called a quantized vortex (Cf. for example ref.). In this case a circuit integral of a vector as in (3) gives, in addition to the electromagnetic term, a term like $`2\pi n`$, so that we have for a change in phase for a small closed curve around this nodal singularity,

$$2\pi n+e\stackrel{}{B}𝑑\stackrel{}{S}$$
(4)
In (4) $`\stackrel{}{B}`$ is the magnetic flux across a surface element $`d\stackrel{}{S}`$ and $`n`$ is the number of nodes within the circuit. The expression (4) directly leads to the Monopole.
Let us now reconsider the above arguments in terms of recent developments.
The Dirac equation for a spin half particle throws up a complex or non Hermitian position coordinate. Dirac identified the imaginary part with Zitterbewegung effects and argued that this would be eliminated once it is realized that in Quantum Mechanics, space time points are not meaningful and that on the contrary averages over intervals of the order of the Compton scale have to be taken to recover meaningful physics . Over the decades the significance of such cut off space time intervals has been stressed by T.D. Lee and several other scholars . Indeed with a minimum cut off length $`l`$, it was shown by Snyder  that there would be a non commutative space time structure, and infact at the Compton scale we would have (Cf.ref.)

$$[x,y]=0(l^2)$$
(5)
and similar relations. The Planck scale ofcourse, is the Compton scale for a Planck mass.
Infact starting from the Dirac equation itself, we can deduce directly the non commutativity (5) as recently shown . This non commutative feature has also been recently stressed, both in Quantum Gravity and in Quantum SuperStrings theory .
Let us now return to Dirac’s formulation of the monopole in the light of the above comments. As noted above, the non integrability of the phase $`S`$ in (2) gives rise to the electromagnetic field, while the nodal singularity gives rise to a term which is an integral multiple of $`2\pi `$. As is well known  we have

$$\stackrel{}{}S=\stackrel{}{p}$$
(6)
where $`\stackrel{}{p}`$ is the momentum vector. When there is a nodal singularity, as noted above the integral over a closed circuit of $`\stackrel{}{p}`$ does not vanish. Infact in this case we have a circulation given by

$$\mathrm{\Gamma }=\stackrel{}{}S𝑑\stackrel{}{r}=\mathrm{}𝑑S=2\pi n$$
(7)
It is because of the nodal singularity that though the $`\stackrel{}{p}`$ field is irrotational, there is a vortex - the singularity at the central point associated with the vortex makes the region multiply connected, or alternatively, in this region we cannot shrink a closed smooth curve about the point to that point. Infact if we use the fact as seen above that the Compton wavelength is a minimum cut off, then we get from (7) using (6), and on taking $`n=1`$,

$$\stackrel{}{}S𝑑\stackrel{}{r}=\stackrel{}{p}𝑑\stackrel{}{r}=2\pi mc\frac{l}{2mc}=\frac{h}{2}$$
(8)
$`(l=\frac{\mathrm{}}{2mc}`$ is the radius of the circuit and $`h=2\pi `$ in the above natural units).
In other words the nodal singularity or quantized vortex gives us the mysterious Quantum Mechanical spin half (and other higher spins for other values of $`n`$). In the case of the Quantum Mechanical spin, there are $`2\times n/2+1=n+1`$ multiply connected regions, exactly as in the case of nodal singularities. Indeed in the case of the Dirac wave function, which is a bi-spinor $`\left(\begin{array}{cc}\mathrm{\Theta }\hfill & \\ \varphi \hfill & \end{array}\right)`$, it is well known that far outside the Compton wavelength, it is the usual spinor $`\mathrm{\Theta }`$, preserving parity under reflections that predominates, whereas at and near the Compton scale it is the spinor $`\varphi `$ which predominates, where under a reflection $`\varphi `$ goes over to $`\varphi `$. This double connectivity of the Dirac spinor was shown to lead immediately to the same electromagnetic potential we had obtained from the nonintegrability of the phase above, which again was identical to that from Weyl’s gauge invariant theory (Cf.ref. for details).
Let us see all this in a little greater detail . We start with a non integrable infinitessimal parallel displacement of a four vector,

$$\delta a^\sigma =\mathrm{\Gamma }_{\mu \nu }^\sigma a^\mu dx^\nu $$
(9)
The $`\mathrm{\Gamma }`$’s are the Christoffel symbols. This represents the extra effect in displacements, due to curvature. In a flat space, all the $`\mathrm{\Gamma }`$’s on the right side would vanish. Considering partial derivatives with respect to the $`\mu `$-th coordinate, this would mean that, due to (9),

$$\frac{a^\sigma }{x^\mu }\frac{a^\sigma }{x^\mu }\mathrm{\Gamma }_{\mu \nu }^\sigma a^\nu ,$$
(10)
The second term on the right side of (10) can be written as

$$\mathrm{\Gamma }_{\mu \nu }^\lambda g_\lambda ^\nu a^\sigma =\mathrm{\Gamma }_{\mu \nu }^\nu a^\sigma ,$$

where we have linearised the metric,

$$g_{\mu \nu }=\eta _{\mu \nu }+h_{\mu \nu },$$

$`\eta _{\mu \nu }`$ being the Minkowski metric and $`h_{\mu \nu }`$ a small correction whose square is neglected. From (10) we conclude that,

$$\frac{}{x^\mu }\frac{}{x^\mu }\mathrm{\Gamma }_{\mu \nu }^\nu $$
(11)
We can identify

$$A_\mu =\mathrm{\Gamma }_{\mu \nu }^\nu $$
(12)
from the above using minimum electromagnetic coupling exactly as in Dirac’s monopole theory.
If we use (11), we will get the commutator relation,

$$\frac{}{x^\lambda }\frac{}{x^\mu }\frac{}{x^\mu }\frac{}{x^\lambda }\frac{}{x^\lambda }\mathrm{\Gamma }_{\mu \nu }^\nu \frac{}{x^\mu }\mathrm{\Gamma }_{\lambda \nu }^\nu $$
(13)
Let us now use (12) in (13): The right side does not vanish due to the electromagnetic field (12) and we have a non-commutativity of the momentum components of quantum theory. Indeed the left side of (13) can be written as

$$[p_\lambda ,p_\mu ]\frac{0(1)}{l^2},$$
(14)
$`l`$ being the Compton wavelength. In (14) we have utilised the fact that at the extreme scale of the Compton wavelength, the Planck scale being a special case, the momentum is $`mc`$.
From (12), (13) and (14), we have,

$$Bl^2\frac{1}{e}=\left(\frac{\mathrm{}c}{e}\right),$$
(15)
where $`B`$ is the magnetic field.
Equation (15) is the well-known equation for the magnetic monopole. Indeed it has been shown by Saito and the author  that a non commutative spacetime at the extreme scale shows up as a powerful magnetic field.
To recapitulate, the Monopole was shown by Dirac to arise because of two separate issues. The first was the non integrability of the phase $`S`$ given in (2), which gave rise to the electromagnetic potential (3) which was equivalent to the Weyl potential (12) (which latter was dismissed because it was adhoc). The other issue was that of nodal singularities or alternatively the multiply connected nature of space which gave rise to a term like $`2\pi n`$ as in (4). In effect there would be free monopoles. However all this was considered in the context of the usual commutative Minkowski spacetime. Effectively this means that terms $`0(l^2)`$ as in (5) are neglected.
However once such terms are included, in other words once the non commutative structure of spacetime to this order is recognised, firstly the previously supposedly adhoc Weyl electromagnetic formulation automatically follows as in (12) and furthermore the first term in the monopole expression (4) immediately gives the Quantum Mechanical spin, and the elusive monopole appears as the magnetic effect at the Compton (or Planck scale). Indeed in recent times the fact that non commutative spacetime gives rise to spin has been recognized .

## 3 Non linear Equations

Let us now come to non linear wave equation. In diverse areas of physical application a non-linearity leads to what may be called auto catalysis or auto production. For example in the well known trimolecular model (Brusselator) a cubic non-linearity is required in the rate equations. In fact the non-linearity has to be at least of the order three. The rate equations are then of the form

$$\frac{X}{t}=K_1A(K_2B+K_4)X+K_3X^2Y+D_1_r^2X,$$

and a companion equation, where the symbols have the usual significance.
Such a process is also well known in the upper atmosphere, in the formation of the triple oxygen molecule, ozone:

$$O+O_2+mO_3+m.$$

It is also known that the rate equations in a number of biochemical reactions involving enzyme catalysis also exhibit in some limiting cases, cubic terms, for example the Glycolytic pathway (Cf.).
In fact, in most problems of cooperative phenomena physics, for example plasmas or lasers, at least cubic non-linearities are required for such cooperative behaviour.
In other areas, such as zoology or biology or sociology too non-linearities cause auto production (Cf. for a non technical discussion). A very well-known example is the logistic equation.
There is a similar situation in quantum field theory also. It is not surprising that auto catalysis or auto production should be a non-linear feature.
In fact for a linear system not only $`\psi `$ but also $`\alpha \psi `$ is a solution.
However if the system is non-linear, it can be written as

$$M(\psi )=L(\psi )+N(\psi )=0,$$

where $`L`$ denotes the linear part and $`N`$ the non-linear part. In a first approximation, we can take

$$L(\psi _0)=0,$$

and $`N`$ can be linearized by suitably substituting $`\psi _0`$ for $`\psi `$ to get, say, $`L^{(0)}(\psi )`$.
The system would now be approximately described by the linear equation

$$[L+L^{(0)}](\psi )=0.$$

This process, if convergent can be continued. At each stage, the coefficients of $`\psi `$ would depend on the linear approximation of up to that stage. This precisely is the characteristic of auto production.
A similar technique has been used recently for a Ricatti equation derived from the Schrodinger equation of non-relativistic quantum mechanics.

## 4 The Non-linear Schrodinger Equation

In the light of the above comments, a non-linear Schrodinger equation was deduced and used to argue that the origin of inertial mass lies in self interacting amplitudes within, typically the Compton wavelength. A cubic non-linearity associated with auto catalysis is responsible for the generation of the inertial mass.
In this case the equation is given by

$$ı\mathrm{}\frac{\psi }{t}=\frac{\mathrm{}^2^2\psi }{x^2}+H(x,x^{})\psi (x^{})𝑑x^{}$$
(16)
where

$$H(x,x^{})=(\psi (x^{})|\psi (x)>,$$
(17)
All this can be generalized immediately to the three dimensional case.

$$E\psi =\frac{\mathrm{}^2}{2m}^2\psi +H(r,r^{})\psi (r^{})d^3r^{}$$
(18)
If next in an equation like (17), the amplitude for a particle at $`r^{}`$ to be at $`r`$ vanishes outside a small interval, so that a $`\delta `$ function can be introduced in (17), then we have the equation

$$E\psi =\frac{\mathrm{}^2}{2m}^2\psi +g\psi ^3$$
(19)
Let us now consider partial wave decomposition equation (18) in spherical symmetry. This gives

$$[\frac{d^2}{dr^2}+k^2\frac{l(l+1}{r^2}],u=gu^3,k^2=\frac{2m}{\mathrm{}^2}E$$
(20)
Further specializing to the case $`l=0`$ and $`k^20`$ and in the spirit of the considerations in Section 1, neglecting the cubic term in (19), we have in the zeroeth approximation for $`u,ur`$. Now writing in (20), $`u^3`$ as $`u^2u`$, that is linearizing equation (20) and using the zeroeth approximation we get

$$[\frac{d^2}{dr^2}+(k^2r^2)]u=0(k^20)$$
(21)
Equation (21) is the well known Harmonic Oscillator equation with degenerate energy levels. As is well known, a set of Harmonic oscillators as above represents an assembly of Bosons. Thus we have a collection of closely packed nearly zero energy Bosons similar to the Bose-Einstein condensation.
Interestingly the link between Solitons arising from the non-linear equations and Bose-Einstein condensation is being investigated, for example by Khaykovich and co-workers at the ENS Laboratory in Paris and also at the European Laboratory for Non-linear Spectrascopy in Italy. In these experiments a Bose-Einstein condensate of a dilute atomic gas of Lithium atoms is used, the inter atomic interaction providing the non-linearity. This in turn ensures the Solitonic propagation.
All this could have varied applications in fields ranging from Particle Physics to Non-linear Optics, including the possibility of high speed fibre optic communication.


text:
# New Limits on Gravitational Radiation using Pulsars

## 1 Introduction

There are two categories of sources that could produce a detectable level of stochastic gravitational radiation from the distant Universe. One is of cosmological origin and the other involves galaxy evolution. First, low-frequency ($`10^{20\mathrm{to}0}`$ Hz) relic gravitational radiation may be generated during inflation, according to string theory models of the early Universe, via both early evolution of the extra dimensions and decay of cosmic strings. Cosmic strings, while now disfavored as the driver behind structure formation in the Universe, still arise in many grand unified theories of particle physics after the epoch of inflation (Caldwell, Kamionkowski, & Wadley 1998; Hogan 2000; see also Maggiore 2000). Second, the coalescence of massive black holes (MBHs) during galaxy evolution could also produce a detectable level of Gravitation Waves (GWs). (See Lommen & Backer 2001 for more on this.) Three cosmology experiments are being conducted or planned by a number of investigators around the globe to either detect or place new limits on the stochastic background of gravitational radiation: polarization of the Cosmic Microwave Background Radiation (pCMBR), Pulsar Timing Array (PTA), and Laser Interferometer Space Array (LISA).

The PTA consists of a regular observing schedule on a handful of the millisecond pulsars that are most accurate as clocks. So far, we have used the Arecibo Telescope in conjunction with the Green Bank 140ft telescope, but we plan on including the new Green Bank 100-m telescope that has just been commissioned and a collection of other European and Australian telescopes to maximize the observation time and the range of observations available.

Pulsar timing residuals are sensitive to gravitational radiation of wavelengths of about 1 yr, making the PTA complementary to LISA which will measure much shorter wavelengths, and the pCMBR which will measure much longer wavelengths.

Using timing residuals from Arecibo Observatory observations of PSRs B1937+21 and B1855+09 over 8 yr, Kaspi,Taylor, & Ryba (1994, hereafter KTR94) placed a limit of

$$\mathrm{\Omega }_gh^2<6\times 10^8(95\%\mathrm{confidence})$$

where $`\mathrm{\Omega }_g`$ is the fractional energy density in gravitational waves per logarithmic frequency interval and the Hubble constant H$`{}_{o}{}^{}=100h`$ km s-1 Mpc-1. We have connected the KTR94 data to 10 years of Green Bank 140ft telescope data and 3 years of Arecibo Observatory data. The overlapping data sets span a total of nearly 17 years, doubling the baseline of Kaspi et al’s results. The sensitivity of the PTA is proportional to $`1/(timespan)^4`$, and therefore we expect an increase in sensitivity of about 16 over the KTR94 results.

First, in §2 we describe the 3 different sets of observations that have gone into calculating the limits presented in this article. In §4 we connect the three data sets. In §5 we show the calculation of the new gravitational wave limit. §6 compares the precision of pulsar clocks to those of atomic time standards. §7 is an update of the relationship between $`\dot{P}`$ and timing noise presented by Arzoumanian et al. (1994). In §8 we discuss the possibility that a planet orbits pulsar PSR B1937+21. Finally in §10 we present our conclusions.

## 2 Observations

At the NRAO111The National Radio Astronomy Observatory (NRAO) is operated by Associated Universities, Inc., under contract with the National Science Foundation. 140-foot (42.7 m) telescope we observed PSR J1713+0747, PSR B1855+09, and PSR B1937+21 in addition to others 4-6 times per year at radio frequencies near 800 and 1400 MHz from 1989 October to 1999 July. For details on the analysis of these data see the report by Backer et al. (1993) on results from the first half of the data set. Observations of PSR J1713+0747 started after its discovery in 1994 (Camilo, Thorsett, & Kulkarni, 1994).

We have been conducting monthly observations at 0.43 GHz, 1.4 GHz and 2.4 GHz of an array of MSPs using the NAIC222The National Astronomy and Ionosphere Center Arecibo Observatory is operated by Cornell University under contract with the National Science Foundation. Arecibo Observatory 300 m telescope since December 1997. These data are used to make precision arrival time measurements for a variety of astrophysical goals. We used the Arecibo-Berkeley Pulsar Processor (ABPP), which is a multi-channel, coherent dispersion removal processor333‘coherent’ means that the dispersion is removed in the voltage domain prior to power detection. with 112-MHz total bandwidth capability (For a detailed technical description of the hardware see Backer et al.(1997.) For PSRs J1713+0747 and B1855+09 we use 56-MHz bandwidth for observations at 1.4 GHz, and 112-MHz for 2.4 GHz, while for PSR B1937+21 we use 45-MHz bandwidth at 1.4 GHz and 56-MHz at 2.4 GHz.

Calibrated total intensity profiles were formed from signals with orthogonal circular polarization. The profiles were then cross correlated with a template to measure times of arrival (TOAs) relative to the observatory atomic clock. Small errors in the observatory UTC clock, of order 1 $`\mu `$s, were corrected based on comparison of local time to transmissions from the Global Positioning System of satellites (GPS). The templates used for cross correlations were constructed by fitting a set of Gaussian components to long-term averages of observations, using the software described in Kramer et al. (1994) and Kramer (1994). This model fitting scheme is described extensively in Lommen (2001). We use the model to generate noise-free templates with a specific common fiducial point at all frequencies. Lommen (2001) shows that the templates generated in this manner have at most 2 $`\mu `$s of error in the fiducial point between 1400 and 2380 MHz.

In addition to the two data sets mentioned above, we also have incorporated the archival Arecibo data from KTR94 on PSR B1855+09 and B1937+21 into our analysis. See the original paper for details on these data.

## 3 Analysis

As we just mentioned, TOAs are calculated calculated via cross-correlation with a template. Before these TOAs can be analyzed several corrections are required:

$$TOA^{}=TOA+\mathrm{\Delta }_{clock}+\mathrm{\Delta }_{DM}+\mathrm{\Delta }_{backend}$$

$`\mathrm{\Delta }_{clock}`$ is comprised of a number of clock corrections which in our case transform a TOA which is referenced against an observatory maser, to one which is referenced against UTC. $`\mathrm{\Delta }_{DM}`$ causes TOA to be the time that the pulse would have arrived had it been at infinite frequency (lower frequencies are delayed more by the ISM). $`\mathrm{\Delta }_{backend}`$ corrects for any delays that are dependent on the observing backend used. In our case we have a digital latency which is dependent on channel bandwidth and a mid-scan correction which is a result of our on-line folding. After all of these corrections, TOA is independent of observatory clock, frequency, and observing backend. TOA marks the space-time event of the arrival of a defined fiducial point of the pulse at the telescope on a defined time scale.

## 4 Connection of Data Sets

In the TEMPO analysis arrival time of pulses at the observatory is referred to that at the barycenter of the solar system. The solar time corresponding to the rotation of the earth is known as UT1. UT1 is tabulated by measuring the difference between UT1 and UTC. UTC is a solar time scale that runs at the rate of TAI, international atomic time. BIPM tables of (UT1-UTC) are used to orient the observatory in inertial space and make the translation from observatory to earth-center. (UT1-UTC) is between 7 and 32 seconds over the course of our observations, and is known to 0.1 ms at all times. This translates to an uncertainty in TOA of only $``$ 0.1 ns. To calculate the earth’s position TEMPO uses standard JPL ephemerides DE405.1950.2050. We know the geodetic position of each observatory to within meters. The delays that remain, that are most difficult to account for, are the delays between the arrival of the pulse at the fiducial point of the antenna as defined by VLBI coordinates and the completion of its transmission through the ‘back-end’ of the receiving system. These are typically in the range of 100-3500 ns. At Arecibo, for example, we know that the travel time for a signal from the control room to the Gregorian receivers and back is approximately 7 $`\mu `$s 444Mike Nolan, private communication, 2001 Aug 19 largely resulting from the $``$1200 ft of fiber-optic cabling through which the signal must travel from the platform to the control room. This would give us a $`3.5\mu s`$ delay in the TOAs. Typically, we ignore these delays because we take data from a single observatory, with a single operating system, and any constant delay such as this is seemingly irrelevant.

However, these delays become important for connecting multi-observatory data sets. If all such delays at all observatories engaged in pulsar timing were known, it would be possible to have a “Universal Pulsar Timing Array (UPTA)” in which all TOAs from all pulsar data from all over the world could be combined a priori in single meaningful data base. This would necessitate the need for uniform template fiducial point processing. For example, if the world’s pulsar astronomers could piece together 17-year data sets on 10 different MSPs with comparable precision, the sensitivity of the UPTA would go down by a factor of 1/$`\sqrt{10}`$.

In our comparatively small experiment, we attempt this, but can only account for everything but the final 15-60 $`\mu `$s, and so in the end we end up doing a bootstrap connection. We have 3 overlapping data sets. For each pulsar, the average separation between two overlapping sections is used as the offset between the two data sets. These offsets, which will be included in the ITOA format that we will make publicly available, are shown for each pulsar and for each adjacent data set in Table 1. In the case of the Green Bank and ABPP data sets, the TOAs were generated using the same series of Gaussians with the same fiducial point, so all template discrepancies are removed to within 2 $`\mu `$s (see Lommen 2001 for more on this subject). 2 $`\mu `$s in fact, represents the maximum error in going between 1400 and 2380 MHz. The data we have connected is centered from 1330 MHz to 1420 MHz, and the error roughly scales with the size of the bandwidth over which one must construct templates (see Lommen 2001). Therefore the maximum error in connecting our range of data, due to our template cross-correlation is roughly a tenth of 2 $`\mu `$s, or 0.2 $`\mu `$s. In the case of the offset between the KTR94 data sets and the Green Bank data sets, we added, a priori, an amount of offset corresponding to the difference between our fiducial point and the fiducial point of KTR94, who used the peak of the template. These a priori numbers were a phase of 0.443259 for B1855+09 and 0.353460 for B1937+21. The numbers shown in the table, are the offsets in addition to these template offsets. The sign of the offsets are meaningful, e.g. the numbers under the column “GB-ABPP” are the quantities acquired from subtracting the average ABPP residual from the average Green Bank residual during the epoch while they overlap. The offsets between data sets that we processed ourselves (GB - ABPP) are all less than 8 $`\mu `$s.

In order to perform a weighted fit across multiple data sets, it is necessary to be very careful of the average weighting of one data set relative to another. Said another way, with a single set of data, all taken with a single observing system and analyzed uniformly, one need not worry about, for example, overestimating all the error bars by a factor of 3. However, when one is combining two such data sets, if one has overestimated error bars, and a second has underestimated error bars, the fit will tremendously favor the latter. In order to account for this problem, we normalized the error bars in each data set by the total variance per day in those data, i.e.

$$\sigma _{\mathrm{normalized}}=\sigma _{\mathrm{original}}\times \frac{\sqrt{_i1/\mathrm{err}_i^2}}{L}$$

where $`L`$ is the length of the data set, in days, and $`\sigma _{original}`$ are the error bars that were published by KTR or calculated by us, using the method created by Christophe Lange (private communication) which was adapted from Downs & Reichley (1983).

We used the standard TEMPO package555http://pulsar.princeton.edu/tempo to perform weighted fits to $`\alpha `$, $`\delta `$, $`\mu _\alpha `$, $`\mu _\delta `$, $`P`$, and $`\dot{P}`$, in each of the 3 pulsars. Additionally in J1713+0747 and B1855+09 we fit for the 5 Keplerian parameters and the Shapiro delay parameters, m2 and $`\mathrm{sin}i`$ (See Stairs, 1998 for an excellent summary of the process of fitting for the Shapiro delay).

We show the resulting residuals in Figure 1. In (a), (b), and (c), we show J1713+0747, B1855+09, and B1937+21, respectively, fitted for the parameters not marked with footnote “c” in the table. Additionally, in panel (d), we show B1937+21 fitted for $`\ddot{P}`$($`\ddot{\nu }`$).

In the next three sections, we discuss various implications of the residuals resulting from these fits. First, we assume the residuals shown in Figure 1 are the result of GWs and use them to place new limits on the GWB. Second, we assume the same residuals are the result of timing noise intrinsic to the pulsars, and determine whether these pulsars are more or less noisy than expected based on the characteristics of the population at large. Finally, we consider the possibility that the residuals in B1937+21 are the result of a planet orbiting that body.

## 5 Limit on Gravitational Radiation

These data, representing a continuous 17-y set, provide a new limit on the level of background gravitational radiation present throughout the galaxy, and ostensibly the universe. In short, we assume that all the observed fluctuations in the timing residuals are due to the stochastic background of gravitational radiation. This yields the maximum possible gravitational wave spectrum present.

To determine the limit they place, we first need to measure the spectrum of fluctuations that we observe in the timing residuals. The process of computing an observed noise spectrum from irregularly sampled data has been given much consideration by previous authors (Groth 1975; Deeter & Boynton 1982; Deeter 1984; Cordes & Downs 1985; KTR94). These authors are particularly concerned with “red” spectrum where imperfect sampling leads to coupling of spectral estimates at low frequency. Deeter (1984) compares the various methods and determines that the method of fitting the data to orthonormal polynomials produces the most meaningful results. We therefore use this method, as did both Stinebring et al. (1990) and KTR94 in previous generations of this experiment.

Stinebring et al. (1990) and KTR94 start with a calculation which shows that the energy density of background gravitational radiation per frequency octave would have a frequency dependency $`\nu ^5`$ and they go about putting limits on that assumed spectrum. We do the same. However, Rajagopal & Romani (1995), Phinney (2001), and Jaffe & Backer (2002) recently show that the energy density per frequency octave would be proportional to $`\nu ^{4.3}`$ rather than $`\nu ^5`$. Either way, this is steep compared to the intrinsic spectrum observed in young pulsars, and would not change our results significantly.

A series of polynomials, $`p_i(t)`$ where $`i`$ is of the set $`0,1,2,3`$ of corresponding order, is generated orthogonally over the sampling of the timing residuals. We started with the following set:

$$p_{j=0,1,2,3}(t)=\underset{i=0}{\overset{j}{}}t^i$$

The process of orthogonalization is done via standard Gram Schmidt reduction. A linear combination of the $`p_i`$’s is found which minimizes the quantity

$$\chi ^2=\frac{1}{\sigma ^2}\underset{n=1}{\overset{N}{}}\left(r(t_n)\underset{j=0}{\overset{3}{}}C_jp_j(t_n)\right)^2$$

where $`r(t_n)`$ is the measured residual at time $`t_n`$. The first three $`C_j`$’s are covariant with the fitting of the phase, period, and period derivative in the pulsar model. $`C_3`$, however, is a measure of the amplitude of the spectral density of the variance in the timing residuals at a frequency corresponding to 1/T where T is the length of the data set. $`S_m=<C_3^2>`$ is the corresponding estimate of the power(variance) contained near frequency 1/T. In order to sample the spectrum of $`S_m`$ over multiple frequencies, we divide the data into smaller and smaller subsets in time, by factors of 2. We refer to the divisor as the frequency index. A frequency index of 2 implies the data set was divided in half, and that therefore we are measuring frequencies of 2/T, or periods of about 8.5 years.

The results of the spectral measurements are shown in the solid line Figure 2 for PSRs J1713+0747, B1855+09, and B1937+21, using the fits corresponding to parts a, b, and c of Figure 1, i.e. fits through $`\dot{\nu }`$ only. Though this method works well for irregularly sampled data, the value of the $`S_m`$s are still dependent on sampling. In order to gauge the extent of this dependence we have simulated data sets of various spectral indices, with the same sampling as the actual data. Randomly distributed Gaussian noise was transformed into the Fourier domain, multiplied by a function with spectral index 0,2,3, or 5, normalized to 1 $`\mu `$s2 y at a frequency of 1 y-1, and transformed back to a time series. Once in the time domain, the data were sampled identically to the actual data. 10,000 such realizations were created for each of the 4 spectral indices. The power spectra of each was measured using the method of orthonormal polynomials described above. The results of this analysis are shown in dotted lines in Figure 2 for each of the 4 spectral indices.

The average value of the spectral estimator, $`S_m`$, with spectral index 5 represents the estimate of $`S_m`$ in the presence of a gravitational radiation, $`<S_m(\mathrm{\Omega }_gh^2)>`$, where $`\mathrm{\Omega }_g`$ is the energy density in GWs, $`\rho `$, expressed as a fraction of the closure density, $`\rho _c`$ of the universe:

$$\mathrm{\Omega }_g=\frac{\rho }{\rho _c}=\frac{8\pi G\rho }{3H_o^2}$$

where $`H_o`$ is the Hubble constant. We generally express $`\mathrm{\Omega }_g`$ in terms of $`h`$, i.e. $`\mathrm{\Omega }_gh^2`$ where $`H_o=100h\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$. We need to translate the normalizing amplitude, (1 $`\mu `$s2 y at a frequency of y-1) into a particular value of $`\mathrm{\Omega }_gh^2`$. We use the relationship derived from the definition of $`\mathrm{\Omega }_g`$ (from Stinebring et al. 1990)

$$P_g(f)=\frac{H_0^2}{8\pi ^4}\mathrm{\Omega }_gf^5.$$

where $`P_g(f)`$ is spectral density, i.e. the same quantity we are attempting to measure with $`S_m`$. This means that values of $`<S_m(\mathrm{\Omega }_gh^2)>`$ from data with 1 $`\mu `$s y at a frequency of y-1 need to be divided by 752 in order to represent the spectral density when $`\mathrm{\Omega }_gh^2=10^7`$ which is the value we represent in Table 2.

Figure 2 shows that the measured spectra of both J1713+0747 and B1855+09 are both quite flat compared to the spectral index 5 model. We would expect J1713+0747 to be somewhat flatter than the others simply by virtue of the plot representing a higher frequency region of the spectrum. PSR B1937+21 demonstrates a steep spectrum, comparable to the spectral index 5 model spectrum, but this may be due to other sources as we will discuss later (e.g. see §8).

In order to calculate the expected influence of gravitational waves on the spectral measurements, $`S_m`$, we needed to know the contribution from purely white noise. To this end, we simulated residuals that were only influenced by our known quantities of instrumental noise. We created white noise at the same sampling as that of the three pulsars we are considering. The white data from each observatory were normally distributed random numbers normalized such that their RMS matched the median of instrumental noise from that pulsar at that observatory. The value used for instrumental noise was the standard deviation of the mean of the timing residuals about the mean for a particular day, i.e.

$$stdev=\frac{\sqrt{_i(r_iu)^2}}{N}$$

where $`r_i`$ is a single residual of which there are $`N`$ averaged to form a daily average and $`u`$ is the mean of the residuals on that day. We created 10,000 random realizations for each of the 3 pulsars, and measured their spectrum using the same technique shown above. The average value of these data we call $`<S_m>_w`$.

The values of $`S_m`$, $`<S_m>_w`$, and $`<S_m>_g`$ for each value of the frequency index, m, are given in Table 2.

In order to attempt to rigorously calculate an upper limit to background gravitational radiation present in our residuals, we duplicate the analysis of Thorsett & Dewey (1996, hereafter TD96). There is significant controversy around this method, which we address following the results. TD96 create a statistic, $`Stat_1=mS_m/[<S_m>_w+<S_m(\mathrm{\Omega }_gh^2)>_g]`$, which is normally distributed with $`m`$ degrees of freedom. The statistic is essentially the measured value of the spectral estimator divided by the predicted value, based on Monte Carlo simulations.

We use a Neyman-Pearson test to determine whether the distribution of $`Stat_1`$, for various values of $`\mathrm{\Omega }_gh^2`$ is significantly different than the corresponding statistic in the absence of gravitational radiation. The validity of the Neyman-Pearson test is complicated, but it has a simple construction. Consider a particular measurable, A, with probability distribution S1. Correspondingly, consider another measurable, B, with probability distribution S2. The Neyman-Pearson test essentially tests the extent to which the distributions overlap. If they have considerable overlap, they are indistinguishable, if not, they aren’t. The quantitative test of ‘considerable’ involves computing the line which delineates 5% of the area under the curve, representing the least likely outcomes of measuring A. Consider Figure 3 which shows two probability distributions. The area to the left of the dotted line is 5% of the total area of the solid curve. The mean value of the dotted curve, falls just to the left of the dotted line. Thus, we can say that these two distributions are different from each other at the 95% confidence level.

One of the subtleties of the Neyman-Pearson test is that one of the two hypothesis is considered the ‘null’ hypothesis, and this is the hypothesis one accepts unless the data compels one to choose otherwise. Thus, the Neyman-Pearson test, by conjecture has a preferred answer, i.e. whatever you choose as the null hypothesis. Our null hypothesis will be that there is no gravitational radiation in the universe. This is of course incorrect. This is the first issue McHugh et al. (1996) have with this method. We are inclined to agree that this is not a good construction.

TD96 used a likelihood ratio, which is an unnecessarily complicated way to construct the Neyman-Pearson test. We merely need to compare the two distributions that are shown in the numerator and the denominator of the ratio TD96 use.

$$S1=\mathrm{\Pi }_{m=1,2,4,8}\chi _m^2\left(\frac{mS_m}{<S_m>_w+<S_m(\mathrm{\Omega }_gh^2)>_g}\right)$$

$$S0=\mathrm{\Pi }_{m=1,2,4,8}\chi _m^2\left(\frac{mS_m}{<S_m>_w}\right).$$

S0 represents the null statistic, and is shown in the solid line of Figure 4. S1 represents the statistic when a certain amount of gravitational radiation $`\mathrm{\Omega }_gh^2`$ is present. We have plotted S1 in Figure 4 when its value is equal to the 5% limit that we calculate, namely, $`\mathrm{\Omega }_gh^22.8\times 10^6`$ at a 95% confidence level.

TD96 calculated $`\mathrm{\Omega }_gh^21.0\times 10^8`$ at a 95% confidence level using only the KTR94 B1855+09 data. Using the technique of TD96 our expanded set of B1855+09 data actually place a weaker limit than the KTR94 data alone, because of the large scatter in the Green Bank data points. The technique separates the data into smaller and smaller sections in order to measure shorter frequencies. The average of these shorter frequency sections is therefore falsely enlarged in our calculations. Using only the $`m=1`$ term of our data we have: $`\mathrm{\Omega }_gh^29.8\times 10^6`$ at a 90% confidence level.

McHugh et al. (1996) rejected the analysis of TD96 because of their use of hypothesis testing rather than parameter estimation. They completed a full Bayesian analysis of the same data which resulted in a much less restrictive limit, $`\mathrm{\Omega }_gh^2<9.3\times 10^8`$, than the analysis of TD96 using the same data.

These opaque calculations cover up the essential nature of the results. After KTR94 we can make a much more simple estimate of the limit these data place on the GWB. We use KTR94 equation 6 which relates the energy density, $`\rho `$, in gm cm-3 of a gravitational wave to its affect on TOA in $`\mu `$s, $`A`$, and a frequency, f, in Hz.

$$\rho =\left(\frac{243\pi ^3}{416G}\right)A^2f^4.$$

The largest amplitude sinusoid that one could conceivably fit to the B1855+09 data has $`A=3\mu `$s and $`f=1/17yr=1.9\times 10^9`$ Hz. This gives $`\rho =3.2\times 10^{38}`$ gm cm-3, or using

$$\frac{\rho }{\rho _c}=\frac{8\pi G\rho }{3H_o^2}=2\times 10^9h^2,$$

which is more than an order of magnitude smaller than the limit found by KTR94. We assume $`H_o=100h`$ km s-1 Mpc-1.

As explained in Lommen & Backer (2001), a gravity wave propagating through the galaxy would have an effect both on the emitting site (the pulsar) and the receiving site (the Earth). Neglecting geometrical considerations for a moment, this could imply that the limit we have just derived will on average be a factor of 2 too high, i.e., that the response that we see in the pulsar timing, since it results from effects at both sites, is essentially doubled. The factor is not so simple, but depends on (a) the direction of propagation and polarization of the impinging gravitational wave, and (b) the geometrical relationship of the pulsars to the Earth. For a stochastic background of gravitational waves, however, we can calculate the expected effect. If the GW at the emission and reception sites constructively interfere then $`C_3`$ will increase by a factor of 2, and $`S_m`$ will increase by a factor of 4. If emission and reception destructively interfere then $`S_m`$ will be close to 0 (it will certainly be non-zero due to irregular sampling). The expectation value, i.e. the ensemble average over the 4$`\pi `$ directions of propagation, and over the wavelengths in question is a factor of 2, so the limit is actually half what we, or any of these other groups, measure. This of course leaves aside the issue of the variance in the GWB. The background value we measure using these data place a limit on the GWB in a specific place and epoch: in our galaxy and during the 2 decades in which this experiment took place. Jaffe & Backer (2002) are working on a realistic version of the GWB which takes these specifics into account.

We have assumed up until now, that PSR B1855+09 produces the best limit on the GWB since it demonstrates the most stable timing and a long baseline. However, what if there was a geometrical situation in which B1937+21 was experiencing a large perturbation due to GWs while a null occurred at PSRs J1713+0747 and B1855+09? We first need to consider that the limited range of data on J1713+0747 was preventing our detection of any cubic therein. In Figure 5 we performed the same fit that we did for part (c) of Figure 1 but using only the limited range of the J1713+0747 data. The obvious cubic is gone but there is a ‘sawtooth’ present at the level of 2 $`\mu `$s that we do not see in the J1713+0747, so we conclude that we are seeing no comparable gravitational wave in J1713+0747.

To attempt to find a geometrical situation which produces an approximate null in J1713+0747 and B1855+09 we inspected the range of possible geometrical multiplicative factors, $`(1\gamma )/2`$ over the 4$`\pi `$ sphere of possible gravitational wave directions for each of the three pulsars. See Lommen & Backer (2001) for an explanation of this factor. Figure 6 shows the factor for 2$`\pi `$ values of RA while DEC=0$`^{}`$for each of the 3 pulsars. B1937+21 is shown by the solid line. J1713+0747 is the dotted line, and B1855+09 is the dashed line. The three pulsars are close enough together in the sky that $`(1\gamma )/2`$ is highly correlated over the sphere. We chose the DEC=0 line to show because it essentially produces the most favorable spot on the sphere for a B1937+21 gravity wave enhancement with a near-null in J1713+0747 and B1855+09. At DEC=80$`^{}`$J1713+0747 has $`(1\gamma )/2=0.005`$, B1855+09 has $`(1\gamma )/2=0.05`$, while B1937+21 has $`(1\gamma )/2=0.12`$. That gives B1937+21 a 2.4$`\times `$ enhancement over B1855+09 and a 24$`\times `$ enhancement over J1713+0747. The relative strengths of J1713+0747 and B1855+09 can be exchanged somewhat as you can see in Figure 6. We conclude that it would be very difficult to arrange for the ‘deviant’ residuals we see in B1937+21 to be produced by gravitational radiation without producing a similar, but slightly smaller, $`10\mu `$s effect, in either J1713+0747 or B1855+09. In addition, a gravitational wave disturbance, which, when reduced by a factor $`(1\gamma )/2=0.12`$ still yields a 30 $`\mu `$s residual deviation is unrealizable using any known gravitational wave producer in the universe (see Lommen & Backer 2001).

Finally, we consider whether emission and reception site disturbances could be destructively interfering in the case of B1855+09 and constructively interfering in the case of PSR B1937+21. The stochastic GWB can be thought of as “crinkled” space-time. We approximate the crinkling as a sine-wave of equal amplitude but arbitrary phase near each of the 2 pulsars and also near the earth. We begin by computing the enhancement of the amplitude of a sine wave by adding two equal sine-waves together, with a variable displacement, $`\varphi `$, between the two. This enhancement we call $`E(\varphi )`$. Thinking of B1937+21 and the earth as emitter and receiver, we may, for example, calculate the probability that $`E(\varphi )`$ will be 5 or greater: 0.001. To obtain the probability that B1937+21’s perturbation will be enhanced over B1855+09’s perturbation by a factor of 25 (roughly the ratio of their peak deviations) or greater we compute the quotient:

$$Q(\mathrm{\Delta }\varphi )=\frac{E(\varphi )}{E(\varphi \mathrm{\Delta }\varphi )}$$

for $`\mathrm{\Delta }\varphi =[0,2\pi ]`$. The probability that $`Q(\mathrm{\Delta }\varphi )>25`$, 0.016, is the probability that the disturbance in B1937+21 is 25 times that in B1855+09. Thus, the residuals we observe in B1855+09 render it highly unlikely that the large amplitude quasi-sinusoid observed in B1937+21 is the result of a GW. Consideration of PSR J1713+0747 would suppress the probability further but not by an equal factor.

## 6 Fractional Stability of Terrestrial Atomic Time Standards vs. Fractional Stability of Neutron Star Rotations

At the baselines considered in this paper, namely longer than 10 years, the fractional stability of neutron star rotations rivals that of terrestrial atomic time standards. To make the comparison quantitative we used the statistic $`\sigma _z`$ proposed by Matsakis, Taylor, & Eubanks (1997) for describing pulsar and clock stabilities. The recipe for computing $`\sigma _z`$ given in Matsakis, Taylor, & Eubanks (1997) is very complete. We only give a bare outline of the computation process here. We divide the timing residual or clock comparison, $`X(t)`$, into smaller and smaller subsections of time, starting with the full length, T, and going to T/2, T/4, T/8, etc. To each subset we fit the function

$$X(t)=c_o+c_1(tt_0)+c_2(tt_0)^2+c_3(tt_0)^3.$$

$`\sigma _z`$ is related to the average of $`c_3^2`$ by the following formula

$$\sigma _z(\tau )=\frac{\tau ^2}{2\sqrt{5}}<c_3^2>^{1/2}$$

where the angular brackets denote the average, and $`\tau `$ is the length of the subset of data. Figure 7 shows the fraction stability, $`\sigma _z`$ vs. dataspan for \[UTC-GPS\], \[TAI-PTB\], \[TAI-USNO\], and \[PSR-TAI\] for PSRs B1855+09, B1937+21. To compute the error we used the estimate provided by Matsakis, Taylor, & Eubanks (1997), which underestimates the error in the case of data with significant red noise, as is evident in the pulsar data.

Figure 7 shows that PSR B1855+09 beats the \[TAI-USNO\] time scale at the longest time periods, and is essentially equally as accurate as \[TAI-PTB\] at timescales of 19 years. Since the fractional stability of B1855+09 improves steadily with increased dataspan, we expect that as additional data on B1855+09 is acquired, this pulsar will easily beat \[TAI-PTB\] at $``$25 year timescales unless TAI and PTB improve or the pulsar becomes less stable. Note that TAI is significantly better than USNO alone even though TAI is an aggregate that contains USNO, and similarly, for TT(BIPM) which contains TAI. It is also interesting to note that the PTA allows the possibility of ignoring the earth clock altogether; one degree of freedom in the PTA can be used to solve for time, i.e., PSR vs PSR.

## 7 Update on Timing Noise as a Function of Period Derivative

Pulsar astronomers have long wondered whether the noise they see in timing residual plots, such as we show in Figure 1 is due to intrinsic instability in the rotation of the pulsar, or is an error in the model, such as a missing planet or unmodeled DM variations. Arzoumanian et al. (1994) published a compendium of measurements of “timing noise” in a number of slow and millisecond pulsars, and suggested that there is a relationship between timing noise and $`\dot{P}`$, i.e., the larger the $`\dot{P}`$ the larger the timing noise. This implies first that the noise is intrinsic to the pulsar, and second that for each pulsar there is a minimum RMS in the residuals that no amount of accounting for systematics will lessen. This has obvious implications for the PTA, and our effort to achieve sub-$`\mu `$s RMSs on a collection of MSPs. The data presented here along with recently published results on MSPs allow us to update the picture of timing noise in MSPs as a function of period derivative. This was most recently done by Arzoumanian et al. (1994) who quantified the noisiness of a pulsar using the noise parameter, $`\mathrm{\Delta }_8`$, defined by

$$\mathrm{\Delta }_8=\mathrm{log}\left(\frac{1}{6\nu }|\ddot{\nu }|t^3\right),$$

where $`t=10^8`$ s. Arzoumanian et al. (1994) quote the best fit result to the scatter plot

$$\mathrm{\Delta }_8=6.6+0.6\mathrm{log}\dot{P}$$

which brought forth the rather depressing notion that we would never achieve better than $`1\mu s`$ accuracy for the MSPs owing to the fact that $`\dot{P}`$ for the population is $`10^{20}`$ or higher.

We are happy to report that some of the MSPs used to create the low $`\dot{P}`$ end of the spectrum have been determined to have lower $`\mathrm{\Delta }_8`$’s than were presented in Arzoumanian et al. (1994). We present the new MSP data in Table 3 and the corresponding plot in Figure 8.

The value used for PSR B1937+21 is the best fit to $`\ddot{\nu }`$ letting all the parameters vary. Arzoumanian et al. (1994) were only able to place an upper limit on $`\mathrm{\Delta }_8`$ for PSR B1855+09 and so are we, but ours is significantly lower, -6.9, down from -6. This comes from the best-fit value of $`\ddot{\nu }`$ of $`8\pm 4\times 10^{29}`$, where the uncertainty quoted is that given by TEMPO. This is significantly smaller than the value previously obtained by KTR94, $`\ddot{\nu }=1.0\pm 0.9\times 10^{27}`$. The value for $`\ddot{\nu }`$ that we use to calculate $`\mathrm{\Delta }_8`$ is different from the value one obtains by only allowing $`\ddot{\nu }`$ to vary. With only 1 degree of freedom, the value is not an appropriate upper limit. For J1713+0747 the fit refines nicely to give $`\ddot{\nu }=2.6\pm 0.2\times 10^{27}`$.

PSR J0437-4715 was discovered by Johnston et al. (1993). This 5.8-ms pulsar is very bright and has been shown to be very good for timing. The RMS of the residuals, folded and averaged at the binary orbital period, is only 35 ns.

Since a third body was discovered around PSR B1257+12, a new upper limit on $`\ddot{\nu }<1.35\times 10^{25}`$ has been published by Wolszczan et al. (2000). The value is an upper limit since Wolszczan suggests the distinct possibility of a fourth body. While these new measurements do not change the value of $`\mathrm{\Delta }_8`$ significantly from the value reported by Arzoumanian et al. (1994), it is now clear that the value represents an upper limit.

PSR B1534+12 has an updated $`\ddot{\nu }<6\times 10^{28}`$ (I.Stairs, private communication) which yields $`\mathrm{\Delta }_8<5.4`$.

PSR B1957+20 is possibly influenced by DM variations, and also by an instable orbital system, as shown by Applegate & Shaham (1994) and Arzoumanian, Fruchter, & Taylor (1994). On the basis of the relationship between DM and RMS DM shown by Backer et al. (1993) we would expect PSR B1957+20, with a DM of 29, to have an RMS DM of about 0.0002 cm-3pc. We can calculate the expected influence of these DM variations on the estimation of $`\ddot{\nu }`$ which was done on single-frequency data at 430 MHz. A DM of 0.0002 cm-3pc corresponds to a timing residual of 5 $`\mu `$s. A 5 $`\mu `$s trend over the course of the observations performed by Arzoumanian, Fruchter, & Taylor (1994) for example, could produce a false $`\ddot{\nu }`$ of $`1\times 10^{25}`$ given the 5-year observation length. This $`\ddot{\nu }`$ corresponds $`\mathrm{\Delta }_8=4.6`$, which is the value attributed to B1957+20 by Arzoumanian et al. (1994). In other words, the timing noise measured by Arzoumanian et al. (1994) could be entirely due to unmodeled DM variations. We therefore make this point an upper limit and emphasize that the true value is probably much smaller.

It seems, then, that a significant fraction of MSPs fall below the line fitted by Arzoumanian et al. (1994). Either the relationship is not well described by a line, or the line has a steeper slope than was published. Either way, this bodes well for the PTA, which relies on submicrosecond accuracy in a handful of pulsars.

Figure 8 is somewhat misleading in that it makes the $`\ddot{\nu }`$ measured in B1937+21 look average. Many of the $`\ddot{\nu }`$’s measured are from single-frequency timing, so DM variations may be contributing to what we are calling ‘timing noise.’ There is, in fact, no MSP that displays such a convincing, significant $`\ddot{\nu }`$. We measure $`\ddot{\nu }`$ in B1937+21 to 4 significant digits, where as most other MSP measurements are upper limits or at best known to 20%. It is possible that B1937+21 is displaying noise whose source is internal to the neutron star crust. Alternatively, perhaps B1937+21 has a magnetic field structure that is evolving with a 20-y timescale, but we see such evolution at this magnitude in no other object. As an alternative we consider the possibility that B1937+21 is hosting a planet.

## 8 Planet Around PSR B1937+21?

The residuals shown in part (c) of Figure 1 are the result of a best-fit model $`\alpha `$, $`\delta `$, $`\mu _\alpha `$, $`\mu _\delta `$, $`P`$, and $`\dot{P}`$, and includes daily corrections for DM. One possible source of their cubic structure is the existence of a planet of binary period, $`P_b=17.6\pm 0.2y`$ which yields an orbital separation of 8 AU, using Kepler’s Law with a central mass of 1.4 $`\mathrm{M}_{}`$. The residuals after including the planet in the model are shown in Figure 9. A$`\mathrm{sin}i`$ is $`27.1\pm 0.7\mu `$s, yielding a companion mass of $`0.08/\mathrm{sin}i`$ $`\mathrm{M}_{}`$, where $`i`$ is the inclination of the orbit to the line of sight. If real, this planet would be the smallest known extrasolar planet besides the lunar-mass planet around PSR B1257+12 (Wolszczan et al., 2000). A planet of similar mass ($`0.15/\mathrm{sin}i`$$`\mathrm{M}_{}`$) has been found around PSR B1257+12 but is much closer to the pulsar (0.20 AU vs 8 AU) Wolszczan (1994).

The planetary model adds 3 free parameters to the standard fit, $`T_o`$, $`P_b`$, and A$`\mathrm{sin}i`$. We compare the goodness of this model to adding the three parameters, $`\ddot{\nu }`$, $`\stackrel{\mathrm{}}{\nu }`$, and $`\stackrel{\mathrm{}.}{\nu }`$. The best-fit planetary model fits the data minutely better than the best $`\ddot{\nu }`$ through $`\stackrel{\mathrm{}.}{\nu }`$ fit (RMS of 1.14 $`\mu `$s vs 1.18 $`\mu `$s). The best fit parameters are $`\ddot{\nu }=1.07(2)\times 10^{26}`$s-3, $`\stackrel{\mathrm{}}{\nu }=1.68(5)\times 10^{34}`$s-4, and $`\stackrel{\mathrm{}.}{\nu }=1.42(3)\times 10^{42}`$s-5. We show in Figure 10 the $`\chi ^2`$ as a function of orbital period, and RMS as a function of orbital period, both of which show a minimum near 17 years, which is the length of our data set. The dotted line in both parts of the figure shows the value of the quantity when the fit to $`\ddot{\nu }`$ is turned on instead of the addition of the planet. Only $`\nu `$, $`\dot{\nu }`$, T0, and a$`\mathrm{sin}i`$ were allowed to vary in the fits. We did not allow $`\alpha `$, $`\delta `$, $`\mu _\alpha `$ and $`\mu _\beta `$ to vary for consistency with the 2-step fitting process described for B1937+21 in §4.

It is possible that the residuals you see in Figure 1 are the result of some intrinsic property of the pulsar (see §7). For years this pulsar was thought to have significant “timing noise” causing the phase to wander as it does (KTR94). We were concerned that there may always be a well in the figure of $`\chi ^2`$ vs $`P_b`$ (e.g. Figure 10) where $`P_b`$ is roughly equal to the length of the data set. We tested this possibility by doing the following. We simulated an extended data set, N years into the future, by assuming that $`\ddot{\nu }`$ is fixed. We then fit a planet to each data set and looked at the position of the well (in other words, the best-fit $`P_b`$) vs N. What we found was that a ‘planet’ will produce a constant $`\ddot{\nu }`$ up through the year 2020, which is as far as we tested. There is a correlation between orbital period and the length of the data set, as shown in Figure 11. Some of the fits are not stable, such as the point at data span=18 y which yields a 3-y period. Notice that in general, the best-fit planetary orbit is much longer than the data set. This is easily understood if you consider that the planetary model needs to imitate the $`\ddot{\nu }`$ and can do so by extracting a fraction of an orbit from a putative planet. The non-continuous jumps to a higher $`P_b`$ shown in Figure 11 are doublings of the orbital period. The $`\chi ^2`$ and RMS for these simulated fits remains very similar to what we observe. We conclude that we could be fooled by a constant $`\ddot{\nu }`$ into supposing the presence of a planet. However, the previous section, §7, shows that the large and highly significant $`\ddot{\nu }`$ in B1937+21 is strange relative to the population of MSPs. The only thing that will resolve this is additional $``$ 10 years of data, assuming the orbit is about 20 years.

PSR B1937+21 is regarded as having too large a $`\dot{P}`$ for being a recycled pulsar, i.e. its large $`\dot{P}`$ alone implies it is a ’young’ pulsar although it is generally regarded as old Backer et al. (1993). We wondered if perhaps the unknown existence of an orbiting planet could have been skewing the measurement of $`\dot{P}`$ since the discovery of the pulsar. In fact, it could not. Our best fit $`\dot{P}`$ without the planet, and without assuming a $`\ddot{\nu }`$ is slightly smaller, by 1 part in $`10^6`$ than our best fit $`\dot{P}`$ with the planet.

## 9 Implications of Planet on Possible Evolutionary Scenarios of PSR B1937+21

Millisecond pulsars are thought to be created in supernovae and spun up by accretion from a companion star during the red-giant phase of the companion. Forming planets around such objects is complex and is the subject of much debate (Wolszczan 1998 and references therein). In fact, finding planets, such as this one, in similar millisecond pulsars helps to solve the mystery in the creation of isolated millisecond pulsars. At the present, of the $``$70 known millisecond pulsars, 9 are solitary systems, meaning the presence of a companion has not been detected. Given that we suspect these objects have been spun-up by a companion, this presents an inconsistency. However, if planets are somehow the remains of an ablated or otherwise significantly reduced companion, then we would expect to find similar planets around the other 8 known millisecond pulsars. A planet with a 17-y orbit or longer, as we describe would not have been found in any other system, simply by lack of dataspan.

## 10 Conclusion

We have demonstrated that connection of data sets with sub-$`\mu `$s accuracy across multiple telescopes and many years is possible and yields high precision results for pulsar model parameters. The new limit that our data place on the energy density in background gravitational radiation, $`\frac{\rho }{\rho _c}=2\times 10^9h^2`$, is below that yielded by the work of KTR94 by more than an order of magnitude.

The data suggest the existence of a small ($`<1`$ $`\mathrm{M}_{}`$) planet around B1937+21. The best-fit orbital period is currently 17.6 y but that may change with the addition of more data. We have ruled out the possibility that the cubic present in B1937+21’s residuals are due to a GW that has geometrically escaped detection in both J1713+0747 and B1855+09. We also show that if the cubic is due to a $`\ddot{\nu }`$ intrinsic to the pulsar, that B1937+21 is unique. We find this alternative therefore unlikely.