{'fraction_non_alphanumeric': 0.024386325665081608, 'fraction_numerical': 0.01863513687186737, 'mean_word_length': 4.8926542976145395, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Application of Millisecond Pulsar Timing to the Long-Term Stability of Clock Ensembles

Terrestrial timescales show instabilities due to the physical limitations of the atomic clocks. Stricter environmental isolation and increased numbers of improved cesium clocks and cavity-tuned hydrogen masers have resulted in time scales more accurate by a factor of about five. The use of different clock ensembles results in measurable changes in some millisecond pulsar timing data. We investigate the possible application of millisecond pulsars to define a precise long-term time standard and positional reference system in a nearly inertial reference frame. Although possible quantitative contribution of the two longest studied millisecond pulsars to terrestrial timescales appears minimal, they may prove useful as independent standards in identifying error sources that are difficult to detect due to the finite lifetime and common reference frame of terrestrial clocks. New millisecond pulsars, perhaps some with even better timing properties, may be discovered as a result of the current global pulsar search efforts.


Introduction
The 1982 discovery of the 1.6 millisecond pulsar B1937+21 1 has provided an object with rotational stability comparable to the best atomic time standards over periods exceeding a few years. 2 The microwave beams of millisecond pulsars, which are rigidly anchored to rotating neutron stars, can act as precision celestial clocks owing to the combined effects of their large rotational energies, 10 50−52 erg, and low energy loss rates. Millisecond pulsars provide an opportunity for astronomical observations to define International Atomic Time (TAI) over long periods.
Because each pulsar has a different spin and spin-down rate, timing measurements of millisecond pulsars cannot provide a more accurate timescale than current atomic timescales. Since the best foreseeable daily-averaged pulsar timing measurements will be accurate to 0.1 microsecond, two observations separated by 5 × 10 7 seconds would be needed to reach the 2.5 × 10 −15 s/s TAI monthly precision 3 (defined as the internal rms of the data set, ignoring systematic errors), although more frequent pulsar observations would reduce this considerably.
The existence of a number of millisecond pulsars distributed across the sky leads to the possibility of timing several pulsars against each other, with terrestrial clocks providing merely the endpoint times and a means of interpolating observations. 4,5,6 Such a program can smooth out variations in TAI, better determine the masses of the outer planets, 7 place limits on perturbations produced by a primordial spectrum of gravitational radiation, 8,9 constrain the nature of interstellar turbulence, 10 and provide an astrometric tie for the planetary and radio reference frames. 11 The search for new millisecond pulsars is a matter of intense activity and steady progress. More than 30 millisecond pulsars have been found in globular clusters, and over 25 have been found in the galactic disk, over a wide range of galactic latitudes and longitudes. Unlike slow pulsars, more than 75% of the millisecond pulsars have stellar companions; millisecond pulsars are often termed recycled pulsars because it is believed they were spun-up by their companions. 12 It is estimated that over 30,000 in the galaxy are beaming in our direction. 13 Non-cluster millisecond pulsars are discovered at a rate of about 1 per 200-300 square degrees. 14 Assuming a factor of two improvement in search efficiency and ignoring serendipity, perhaps up to 400 field millisecond pulsars could be discovered within the next few decades.
In this work we summarize the problems and status of terrestrial timescales and pulsar observations, and explore the contributions pulsars can make to terrestrial timescales.

Terrestrial Clocks
According to the International System (SI), the second is defined so that the frequency of the (F, m f )=(3,0) to (4,0) hyperfine transition of cesium 133 is exactly 9.192631770 GHz, in the absence of a magnetic field and located on the geoid. 15 As recognized in the definition, the great majority of the timing data used to compute TAI are from cesium-based frequency standards. 16 Cesium atoms in the (3,0) state are created in an oven and isolated from cesium atoms in other states by their deflection in an inhomogeneous magnetic field. They are then allowed to travel through a magnetically shielded cavity, in which they are exposed to 9 GHz microwaves at each end. The probability that the 9 GHz microwaves will induce a transition to the (4,0) state has a Ramsey-pattern dependence upon frequency of the microwaves; 17 the number of ions which have made the transition to the (4,0) state is measured by the current created after deflection by an inhomogeneous field at the exit point of the cavity. The principal error sources are the magnetic field correction and the phase difference between the microwave fields at the two interrogation regions. 18 Recently, Hewlett-Packard introduced their model 5071 cesium standards, with a significantly improved electronics package. 19 This model cesium clock has demonstrated superior performance both in the laboratory 20 and in the field. 21 The world's most accurate laboratory standards are maintained by the Physikalisch-Technische Bundesanstalt (PTB) in Germany and the National Institute of Standards and Technology (NIST), which achieve long-term stability through use of special beam-path and field switching techniques and the use of a longer cesium travel path. 22 A second frequency standard used for high precision time keeping is the hydrogen maser. Hydrogen masers are based upon the 21-cm transition of atomic hydrogen and can achieve much greater precision in the short-term. 23 However, their long-term accuracy is limited by variations in the dimensions of the confining cavity and interactions with the cavity walls. 24 Special Teflon coatings have been developed which minimize the wall-shift. Recently the Sigma Tau Standards Corporation has begun marketing masers whose cavity capacitance is stabilized and optimized through an attached varactor diode and a tuning circuit which measures the maser amplitude at 7.5 kHz on each side of the microwave resonance. 25 The precision of the best types of frequency standards at the U.S. Naval Observatory (USNO) as a function of the averaging time are shown in Figure 1, which was adapted from Breakiron. 20 Note that masers have significantly better performance on short timescales, but on the longest scale, 60 days, their performance is similar to the cesiums. Following the methodology of Barnes et al., 26 we summarize in Table 1 the error performance of the clocks as a sum of white and random walk frequency components. Since these clocks have not been in existence long enough to gather data on still longer timescales, it is likely that future analyses will identify other components to the variation of the cesium standards. A long-term analysis of maser data will be more difficult to perform since they undergo discontinuous changes in frequency and frequency drift (the time derivative of the frequency) on timescales of a few months, which are removed by comparison with other clocks.
Within the next decade, we can expect further improvements in terrestrial standards. The closest to implementation are optically pumped mercury stored-ion devices. Their advantage is based upon the high frequency of the observed transition (40.5 GHz) and coherence times which allow line-widths as low as 17 mHz. Although the accuracy of prototypes was limited by vacuum contamination, 27 recent designs based upon a linear ion trap have achieved frequency stabilities better than 10 −15 s/s on timescales of a few hours 28 (frequency is defined hereafter as the time derivative of the timing error in units of time, not phase). Somewhat further away are clocks based upon trapped yttrium, 29 single trapped mercury ions, 30 and cesium atomic fountains, which achieve long coherence times by first cooling atoms in "optical molasses" to 700 nanokelvins and then allowing them to rise out of the measurement region until gravity brings them back. 31

Global Time Transfer
The Global Positioning System (GPS) has within the last few years become the chief instrument for comparing clocks of different institutions, as well as for determining position on the Earth. It consists of a constellation of 24 satellites, run by the U.S. Air Force Second Satellite Operations Squadron, each containing at least one cesium frequency standard and circling the Earth twice per sidereal day. 32 These satellites constantly broadcast their coded position and timing infor-mation, in both classified and unclassified channels, enabling anyone with the proper equipment to determine, by triangulation, the position and time at the receiver. Depending upon the kind of GPS receiver, the time can be communicated through an output signal and/or compared digitally to the time input from a standard at the site. The individual timing characteristics of the satellites are monitored by the USNO, which forwards the information to the Air Force for use in determining the frequency and time offset corrections and in measuring orbital variations. Although the specified GPS time precision is 1 microsecond, since January 1993 GPS time has kept within 260 nanoseconds of the USNO Master Clock #2 (UTC(USNO)), and the rms difference over that period was 60 nanoseconds. Satellites also broadcast on-line predictions of GPS−UTC(USNO), which are specified to be less than 90 nanoseconds, but are usually below 20 nanoseconds (which is close to the measurement error). Evaluations of GPS−UTC(USNO) are published as USNO Series 4, and are also available via modem, Internet, and the World Wide Web home page http://tycho.usno.navy.mil.
While the received GPS signal is distorted and delayed by unmodeled tropospheric, ionospheric, relativistic, and multi-path effects, the largest error contribution to unclassified use is the intentional distortion of up to 340 nanoseconds imposed by the military for security reasons and termed "selective availability" (SA). SA consists of both the insertion of a variable time delay and a misrepresentation of satellite orbital parameters, which are used to determine the propagation time for the signal to reach the observer.
If one is merely interested in time difference between remote clocks, it is possible to eliminate the error due to the variable time delay through the "Common View" method. In Common View, observations are scheduled so that both sites observe the same satellite at the same time. While the difference neatly cancels satellite timing errors, it cannot eliminate the errors due to the incorrectly broadcast orbital parameters. It is possible to correct for this after the fact, using satellite orbital elements routinely disseminated by the International GPS Service for Geodynamics, which relies upon multiple observations for its analysis. Common View techniques routinely provide time-transfer precision of 5 to 15 nanoseconds. 33,34 While it is possible, in principle, to correct any GPS measurement in real time for SA using Common View, such scheduling tends to result in fewer observations which are also conducted at lower elevation angles, so that tropospheric and ionospheric modeling errors are more significant. At the USNO, unclassified GPS transfers are achieved using the "melting-pot" method, which forsakes Common View in order to obtain a straight forward average of as many satellites as possible. Even in the presence of SA, a precision of 20 nanoseconds is attainable if one is able to average at least seventy 13-minute satellite passes over a 48-hour period. 35 The advantage of this method is even greater for receivers which can observe more than one satellite at the same time. Another form of time-transfer consists of using available non-GPS satellites to send signals to and from the two remote sites. Aside from achieving sub-nanosecond time transfer precision, such studies often reveal systematic and perhaps seasonal differences with conventional GPS measurements, of order 10 nanoseconds. 36

Terrestrial Timescales and the BIPM
The goal of timescale formulation is to approximate "true" coordinate time, from an average of differential clock timing measurements. In essence, the clocks serve as frequency standards to create a heterodyned tone of known or defined value. Individual timing events are realized either through summation of cycles or use of "zero-crossings" of the sinusoidal output voltage. In practice, timescales are usually generated in the frequency domain, in which a mean frequency is computed and used to define the time from an arbitrary offset.
Internationally, the official responsibility for timescale generation lies with the International Bureau of Weights and Measures (BIPM) in Paris, France as first established by the Convention du Metre in 1875. The BIPM produces TAI by independently evaluating the individual characteristics of approximately 200 clocks maintained by almost 50 institutions. 37 Evaluations are made at 10-day intervals, and reported monthly in their Circular T. The BIPM also evaluates timescales generated by many national institutions, including GPS time. Some of these timescales are free-running, such as the USNO "A.1" series, which is currently comparable to TAI in accuracy. Others, such as UTC(NIST) and UTC(USNO), are steered so as to track UTC closely. The UTC(USNO), which is the national standard for legal and military matters, is available on-line and is in practice based upon the output of a single maser that is steered daily through the addition of small frequency offsets so as to approach the computed mean of all contributing USNO clocks, which itself is adjusted towards extrapolations of the BIPM-computed TAI−UTC(USNO). As a result of this steering TAI and UTC(USNO) have since 1991 kept within 200 nanoseconds.
Once all the time comparisons have been assembled, the problem of timescale generation reduces to calculating weights and offsets for each standard. The BIPM first generates a free-running timescale, the Echelle Atomique Libre (EAL). The EAL is computed by evaluating clock weights each 60 days. The weights are determined from the variance of the frequency offsets over the past year, after allowing for any intentional rate adjustments. A maximum weight is chosen so that the timescale will not depend upon just a few of the clocks. While both the general and specific procedures are periodically under review, in 1994 about half of all clocks used had this maximum weight, while most of the rest were weighted less than 20% of the maximum. The BIPM does not subtract offsets from the clock frequencies, and this can lead to minor discontinuities as clocks are added, subtracted, or reweighted. The TAI timescale is created by gently steering the EAL so that its long-term behavior conforms to the "primary" frequency standards, which are currently the two "long-tube" primary standards of the PTB 22 and NIST-7. 37 This is a practical way for a real time system to take into account the fact that the primary standards have greater long-term accuracy than accounted for by the precision-based weighting system. Again, there is a minor problem in that no allowance is made for fact that the NIST-7 frequency is corrected by about 1.9 × 10 −14 s/s for the AC Stark effect due to the microwave background experienced by the cesium atoms, 38 whereas the PTB frequencies are not.
At the USNO, the A.1 timescale is currently determined using a weighting function which weights clocks by their type. Recent maser data are given a much higher weight than recent cesium data, but older maser data are progressively downweighted, with zero weight assigned to maser data more than 60 days old. 39 Also, data from the older model cesium standards are included at 0.65 the weight of model 5071 cesium clocks. Unlike the BIPM, each clock is introduced to the mean after subtracting an initial frequency offset. This has the effect of reducing discontinuities, and the drawback of subtly over-weighting initial clocks is soon lost in the randomwalk noise as long as no systematic errors are present. For historical reasons, the A.1 is intentionally steered a constant 19 ns/day. 40 Using data from 1993 and 1994, the BIPM estimates the 100-day error of the EAL as slightly less than the TAI, and both to be about 2.5×10 −15 s/s; 3 the errors in the 1980's were above 10 −14 s/s. The recent improvement in both TAI and the USNO free-running timescale A.1 is evident in Figures 2-4, which gives their difference as a function of time. Part of this apparent improvement is due to the fact that the contribution of the USNO clocks to TAI has varied (from 30% in the 1980's to a low near 15% in 1990 to almost 40% in 1994) but most is due to the world-wide increase in number, quality, and environmental isolation of the individual frequency standards. Similar increases in precision are shown in the comparisons with the independent PTB timescale.
For practical reasons, other timescales are often used which differ from TAI by specific and calculable amounts. Coordinated universal time (UTC) differs from TAI only by an integral number of leap seconds, which crudely take into account variations in the rate of Earth rotation. 41 These leap seconds are extra seconds included in UTC at UT midnight on either December 31 or June 30 so as to minimize UT1−UTC, where UT1 is time related to the rotational angle of the earth after allowing for polar motion. 41,42,43 They have lately been inserted at a rate approaching once a year, and on July 1, 1994, TAI−UTC= 29 seconds. Measurements and predictions of UT1, polar motion, and nutation parameters 44 are available on-line from the home page http://maia.usno.navy.mil. In order to allow for continuity with ephemeris time (ET), Terrestrial Time (TT, also called Terrestrial Dynamical Time, TDT) is defined to be offset by exactly 32.184 seconds from the "ideal" form of TAI, which is what TAI would be if no instrumental noise were present. Although in the future TT may be realized from an average of many sources, including pulsars, currently TT is realized only from TAI, by post-processing the raw clock data and improving upon the initial weighting and steering corrections. 45 Since 1982, the BIPM-computed numerical difference TT−TAI has varied by over 3 microseconds absolutely, and over 2 microseconds after removing the quadratic component of the variation. In order to provide continuity in the interpretation of old data and allow for the effects of general relativity in a consistent manner, the IAU in 1991 passed resolutions recognizing timescales which distinguish between centers of reference at the solar system barycenter (Barycentric Coordinate Time, TCB), the center of the Earth (Geocentric Coordinate Time, TCG), and the geoid (TT, TAI,  UTC, ET). 46,47,48

The Techniques of Pulsar Timing
The precision of an individual pulsar time-of-arrival (TOA) measurement depends upon its pulse period and the signal-to-noise ratio as determined from the particular instrumental setups. 49,50 The fastest pulsar known today is still PSR B1937+21, whose 1.6 ms period is close to the theoretical limit determined from balancing the centripetal and gravitational forces, assuming a typical neutron star equation of state at densities of > ∼ 10 14 g cm −3 . 51 The simplest rotation model for a solitary neutron star can be given in terms of phase residuals as a power series in time: where φ is the pulse phase, Ω is the rotation frequency,Ω is the spin-down rate. The spin frequency is related to the pulse period by Ω = 2π/P (Ω = 2πν), whilė Ω = −2πṖ /P 2 . If the spin-down of the pulsar is driven by magnetic dipole radiation, due to misalignment of the magnetic and rotational axes, then the braking index n, defined by ΩΩ/Ω 2 ≡ 2 − PP /Ṗ 2 , is 3. If the initial pulsar period was P 0 , the pulsar age can be computed from Most pulsarsP measurements yield only upper limits, but in those pulsars for whicḧ P can be measured above the noise, n is found to be between 2 and 3. Inverting the equations above, we derive a theoretical valuë In the case of PSR B1937+21 the theoretical slowdown would result in a 0.1 microsecond residual from a constantṖ model after 11 years. The relativistic time-offlight effects and acceleration due to the galactic gravitational field can also make small, but non-negligible contributions toΩ and the theoreticalΩ. 52 Also, for pulsars in globular clusters, contributions toΩ from Doppler acceleration can be important where the gravitational potential of the cluster accelerates the pulsar, influencing the rotation period over several-year timescales. 53,54 Pulsar timing requires an explicit definition of both a steady atomic time standard and accurate monitoring of the motion of the Earth with respect to the background stars. This task is accomplished by using a clock whose offset is frequently calibrated against a real time steered UTC service, such as UTC(NIST) or UTC(USNO). This real time UTC realization is later corrected to a more accurate timescale using after-the-fact determinations published by the BIPM. Although TAI is often used, TT is recommended as the most accurate terrestrial scale available today. 5 Once a terrestrial timescale has been chosen, it must be converted to the barycentric frame (TCB). The omission of this time conversion would, among other things, result in a sinusoidal monthly error of about a microsecond as the Moon pulls the the Earth back and forth through the solar gravitational potential well. In most software, the time conversion happens automatically as part of the four-dimensional general-relativistic transformation which accounts for the Earth's orbit. The IAU has recommended that all users avoid transformations to Barycentric Dynamical Time (TDB), 46 which is scaled to remove secular differences with TT; as a result it diverges from TCB at a constant rate of 1.55 × 10 −8 , which would be absorbed into the fitted parameters.
The geometric, general relativistic, and astrometric effects due to the position and velocity of the antenna with respect to the solar system barycenter must also be removed. In most cases this is done using the IERS values for UT1−UTC and the DE200 ephemeris. 55 Since most observations used for computing the Earth's orbit are Earth-based they are relatively insensitive to the uncertainties in the masses of the outer planets. The most uncertain long-term component of the ephemeris for pulsar work is the motion of the Sun and known planets about the solar system barycenter. The largest error contribution is due to the uncertainty in the mass of Saturn, which could result in a TOA error of order 7 microseconds with a 30-year period. 49 This and other such errors would be identifiable through their different and systematic projections upon the pulsars observed. 4 Fukushima has examined the differences between the DE200, DE102, and DE245 ephemerides and found the differences between them (including their effects on the relativistic time correction TCB−TT) to be negligible for pulsar work. 48 Pulsar signals must be corrected for dispersion caused by the intervening plasma using the frequency dependence of the delay. For a homogeneous and isotropic medium the group velocity is given by v g = c(1 − ω 2 p /ω 2 ) 1/2 , where ω p is the plasma frequency and ω is the wave frequency. The plasma frequency in Gaussian units is ω 2 p = 4πn e e 2 /m, where n e is the mean electron density, and e and m are the electron charge and mass. The delay is proportional to the column density of electrons in the line of sight towards the pulsar. This column density is called the dispersion measure and is defined as DM = z 0 n e dl, where z is the pulsar distance. The dispersion measure is a directly measurable quantity determined from the differential pulse arrival time between two frequencies. To first order, observing v g gives a time delay between two frequencies as Observations of dispersion measures toward 706 radio pulsars give values that range from 1.8 to ∼ 1000 cm −3 pc. 56 For some lines of sight the turbulent electron population may cause different wavefronts to travel along slightly different paths, which can lead to delays that are not removable by the two-frequency formula, although some improvement can occur with multi-frequency observations. The unmodeled delay, which depends upon the total electron column density and the turbulent properties on the interstellar electron population, has been estimated in the worst cases to be at the sub-microsecond level for radio frequency observations at 1.4 GHz and above. 10 Another important consideration is that millisecond pulsars have steep radio spectra and are generally weak enough that observations at the largest telescopes are required in order to achieve an adequate signal-to-noise ratio. An estimate of the SNR achievable at current and future sites is provided in Table 2.
Once the TOAs from a pulsar are assembled and iteratively corrected for known effects, several pulsar parameters are typically solved for. Uncertainty in the a priori pulsar period and spin-down rate are identified through their quadratic dependence. The pulsar parallax is identified from the TOA by its semiannual periodicity, position errors through their annual periodicity, and proper motion by a growing annual periodicity. Pulsar parameters associated with orbits about companions are more complex. 57,58 More exotic error sources, such as gravitational radiation and unknown planets, can be identified through comparison of delays with other members of the pulsar ensemble. If the TOA errors are modeled as a multi-pole expansion over the sky, timing errors are related to the monopole term, ephemeris errors result in a dipole dependence, and gravitational radiation would cause a quadrupole and higher order effect. 4

Stability of Pulsar Time Standards
Using data on the millisecond pulsars B1937+21 and B1855+09 made available by anonymous FTP from Princeton University, 2 the program TEMPO 59 was used to solve for the pulsar periods, spin-down rates, position, and proper motion. The residuals to the solutions are shown in Figure 5. As Kaspi et al. 2 found with the same data, the residuals of PSR B1937+21 appear to display red noise (long-period correlations, in this case well-modeled as a non-zeroP ), whereas the residuals of PSR B1855+09 do not. Since these two pulsars are close together in the sky, the red noise is considered to be intrinsic to the pulsar B1937+21 and not due to the timescale or ephemeris errors. The timing noise for millisecond pulsars seems to follow the correlation withṖ found for slower pulsars by Arzoumanian et al.: 2, 60, 61 where the "stability parameter" ∆ 8 = log 1 6ν |ν|(10 8 s) 3 with a scatter of ± one decade. Pulsar B1937+21 has a ∆ 8 = −5.5. If this inverse correlation is valid for millisecond pulsars, then the PSR J1713+0747 would be a particularly promising candidate to have a very low ∆ 8 stability parameter, 62 with a red noise less than that of PSR B1937+21 by a factor of about 5.
The effects of different input timescales are shown in Figure 6, which plots the difference between the residuals of parameter fits to PSR B1937+21, which differ only in the input timescale assumed. For these and subsequent solutions, the parameters were the pulsar period, period derivative, position, parallax, and proper motion. In Figure 7 these difference solutions are double-differenced with the differences between the timescales, as made available by the BIPM, with a secondorder term removed. The double-differences indicate how much timescale variation is soaked up into the fitted parameters, and hence is indeterminable from pulsar TOA data. The annual signature evident in the plots, for example, shows how pulsar timing data alone will never be able to identify strictly periodic annual variations in terrestrial standards. In the presence of an independent means of determining pulsar positions, proper motion, or parallax (such as astrometric Very Long Baseline Interferometry, or VLBI), such variations would be recognizable as displacements in the fitted parameters. Table 3 shows the current positional precision of a number of millisecond pulsars, derived from timing data. The coupling of the dynamic and radio reference frames will, along with VLBI observations of millisecond pulsars, aid in removing two degrees of freedom from the pulsar fitting models and improve our ability to constrain other sources of errors.
Slightly extending the analysis of Blandford et al., 63 it is possible to apply a Weiner filter to the residuals by Fourier-transforming the residuals, and then transforming back to the time domain after multiplication by a "transfer function," which corrects for the distortion caused by fitting to pulsar parameters. The exact correction factor would also include a factor incorporating the clock and pulsar stabilities as a function of frequency. Since pulsar data are much less accurate on short time scales, compared to terrestrial standards, pulsar data would be given insignificant weights at short periods. For the generation of the plots in Figure 8, a simpler analysis was used in which all periods less than 400 days were discarded. Figure 9 shows the "transfer function" for two different time ranges, which indicate what frequencies in the raw data are absorbed by the fitted parameters. 63 The low values at periods of one year and six months indicate that any noise in the data with those frequencies will be absorbed into the fitted pulsar parameters. It is evident that the accumulation of another decade of data will bring about a considerable improvement in the ability of pulsar timing data to discern long-term variations.
The contribution millisecond pulsar data can make to a timescale is entirely dependent upon the quantity and quality of the pulsar TOAs compared to the terrestrial standards. Included along with the simple clock model in Table 1 are parameters for the timing stability of PSR B1937+21, which is characterizable as a random walk in frequency, or perhaps an integrated random walk, which is even redder. 64 These numbers were chosen to be consistent with the residuals plotted in Figure 5, and they show that PSR B1937+21 is less accurate than a single model 5071 cesium standard. If the relation between timing stability andṖ is valid, then the pulsar J1713+0747 would be five times more accurate, and an ensemble consisting of several such pulsars could approach the accuracy of today's ensemble of terrestrial frequency standards. More important than the numerical evaluations is the fact that pulsars can provide an independent means of evaluating clock performance, and can serve to constrain the amount of extremely red noise in terrestrial clocks, which would be undetectable on the short timescales available so