{'fraction_non_alphanumeric': 0.03669968164131588, 'fraction_numerical': 0.027178398773729514, 'mean_word_length': 4.367879746835443, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Triple Pulsar System PSR B1620-26 in M4

The millisecond pulsar PSR B1620-26, in the globular cluster M4, has a white dwarf companion in a half-year orbit. Anomalously large variations in the pulsar's apparent spin-down rate have suggested the presence of a second companion in a much wider orbit. Using timing observations made on more than seven hundred days spanning eleven years, we confirm this anomalous timing behavior. We explicitly demonstrate, for the first time, that a timing model consisting of the sum of two non-interacting Keplerian orbits can account for the observed signal. Both circular and elliptical orbits are allowed, although highly eccentric orbits require improbable orbital geometries. The motion of the pulsar in the inner orbit is very nearly a Keplerian ellipse, but the tidal effects of the outer companion cause variations in the orbital elements. We have measured the change in the projected semi-major axis of the orbit, which is dominated by precession-driven changes in the orbital inclination. This measurement, along with limits on the rate of change of other orbital elements, can be used to significantly restrict the properties of the outer orbit. We find that the second companion most likely has a mass m~0.01 Msun --- it is almost certainly below the hydrogen burning limit (m<0.036 Msun, 95% confidence) --- and has a current distance from the binary of ~35 AU and orbital period of order one hundred years. Circular (and near-circular) orbits are allowed only if the pulsar magnetic field is ~3x10^9 G, an order of magnitude higher than a typical millisecond pulsar field strength. In this case, the companion has mass m~1.2x10^-3 Msun and orbital period ~62 years.


INTRODUCTION
It is now a quarter century since the discovery of the first binary radio pulsar. Since 1974, about fifty other binaries have been discovered, including systems with a wide variety of companions: planets, "brown dwarfs" below the hydrogen burning limit, white dwarfs, neutron stars, and massive main se-quence stars.
In 1993, the first identification of a pulsar in a triple system was proposed. PSR B1620−26 is an 11 ms pulsar associated with the nearby globular cluster M4. It is in a 191 day, low-eccentricity orbit with a ∼ 0.3M ⊙ companion presumed to be a white dwarf   However, the early pulse timing data were not well described by a simple Keplerian model (Thorsett 1991). Backer (1993) was the first to realize that the data could be well fitted with the addition of a substantial cubic term to the pulsar timing model: a change in the apparent spin-down rate of the pulsar that was so large that it predicted a change in sign of the pulsar frequency derivative-from spindown to spin-up-on a timescale of ∼ 10 yrs. Such a cubic could be induced by a changing acceleration of the pulsar binary under an external gravitational force, which Backer proposed could be due to a second companion in a wide orbit. Further timing observations (Thorsett, Arzoumanian, and Taylor 1993, Backer, Foster, and Sallmen 1993, Backer and Thorsett 1995, Arzoumanian et al. 1996 led to measurements of the next two terms in the Taylor expansion of the pulsar frequency as well as to the first measurements of perturbations in the orbital elements of the inner binary caused by the tidal effects of the outer body. These observations in turn led to tighter constraints on the properties of the second orbit; the second companion most likely has a mass typical of a brown dwarf or planet, ∼ 0.01M ⊙ , and is in a ∼ 40 AU orbit (Rasio 1994, Sigurdsson 1995, Arzoumanian et al. 1996, Joshi and Rasio 1997.
We have now been observing PSR B1620−26 for over a decade. In §2, we describe the various observing systems that have been used. In §3, we describe the analysis of the data. For the first time, we present results from a full, two-orbit analysis, rather than a single Keplerian orbit combined with a Taylor series expansion in pulsar frequency. We describe the constraints that can be placed on the elements of both the inner and outer orbits from the timing data. Finally, in §4, we discuss the implications of our measurements for the understanding of the triple system.

OBSERVATIONS AND ANALYSIS SUMMARY
Since shortly after its discovery, PSR B1620−26 has been observed regularly at multiple radio frequencies using several different telescopes. We report on observations made using the Very Large Array (VLA), near Socorro, New Mexico; the 43 m telescope at the National Radio Astronomy Observatory in Green Bank, West Virginia; and the 76 m Lovell Telescope at Jodrell Bank, England.
At the VLA, observations were made on one or two days each two or three months from 1990 November 30 to 1998 September 21, excluding the period from 28 June 1996 to 28 September 1997-a to-tal of 49 days. A filter bank and the VLA's "High Time-Resolution Processor" were used to divide a 50 MHz bandpass at 1.66 GHz into 14 slightly overlapping 4 MHz channels in each of two orthogonal circular polarizations. The signal in each frequency channel was sampled and averaged synchronously with the predicted topocentric pulsar period, then the channels were shifted to account for dispersion and added, and the polarizations were summed, to produce a single integrated pulse profile every 5 minutes, using a Princeton Mark III pulsar timing system (Stinebring et al. 1992). The start time of each integration was referenced to external time standards using GPS. After eliminating data contaminated with radio-frequency interference (RFI), 486 individual arrival time measurements were available, comprising just over 40 hrs of observing time. Timing precision varied, with a weighted, root-meansquare precision of 31 µs.
At Green Bank, observations were made in campaigns that lasted several days each quarter from 1989 August 20 to 1998 August 7-a total of 165 days. At each epoch observations were made at two frequencies, varying between 400, 575, 800, and 1330 MHz. The "Spectral Processor" fast-fourier transform spectrometer was used to synthesize 512 channels across a 40 MHz passband (256 channels across 20 MHz before 1991 February) in each of two orthogonal polarizations, and to fold the channels synchronously with the topocentric pulsar period. The channels and polarizations were summed, as at the VLA, to produce a single integrated pulse profile for each integration period. Each integration was begun on a time signal from the site maser, and GPS was used to reference observatory time to external standards. After eliminating data contaminated with RFI, 763 arrival time measurements were available, comprising just over 60 hrs of observing time, with a weighted, r.m.s. precision of 40 µs.
At Jodrell Bank, observations were made on 490 days between 1987 October 15 and 21 October 1998. On each day, observations were made at either 400, 600, 1400, or 1600 MHz. At each frequency, two circular polarizations were observed, using a 2 × 64 × 0.125 MHz filterbank at 400 and 600 MHz, and a 2 × 32 × 1 MHz filterbank at higher frequencies.
The signals were detected, filtered, and the polarizations added, then were sampled and averaged synchronously with the predicted topocentric pulsar period. Shifting and adding of channels to account for interstellar dispersion was done in hardware. There were 537 arrival time measurements available, with weighted, r.m.s. precision of 45 µs.
In total, the data set spans exactly 31,391,908,721 pulsar rotations.
Analysis of the pulse arrival times was carried out with the standard software package TEMPO (Taylor and Weisberg 1989). The integrated profiles were cross-correlated with a high-signal-to-noise-ratio average profile to measure the offsets between the start of each integration and the arrival time of a pulse near the center of the integration (Taylor 1992). The arrival times were fitted with a model that included the pulsar phase at a reference epoch, the pulsar frequency f and frequency derivativeḟ , the pulsar position (α, δ) and proper motion (α cos δ,δ), the dispersion measure DM and a linear rate of change of the dispersion measure, and the parameters of the known binary orbit: the period P a , projected semimajor axis x 1a = a 1a sin i a /c, eccentricity e a , argument of periastron ω 1a (measured from the ascending node to periastron in the orbital plane), and a time of periastron T a , as well as other parameters discussed below. 3 The parameters were varied, and the differences between the model and observed arrival times were minimized in a least-squares sense. The data are weighted to account for large variations in their signal-to-noise ratios. Quoted uncertainties are, unless otherwise noted, intended to be approximately 68% confidence limits; in general, they are obtained by testing the robustness of the fitted parameter estimates and formal uncertainties when subsections of data are omitted, and through bootstrap Monte Carlo techniques. Typically, quoted error regions are between one and three times as broad as the formal uncertainties reported by the TEMPO fitting procedure. Timing "residuals"-the observed pulse arrival times minus the model predictions-are calculated; examples of the observed residuals will be shown below.

TIMING RESULTS AND DISCUSSION
In Fig. 1, we display the timing residuals for a model containing only a fit for f ,ḟ , the astrometric parameters, and a single Keplerian orbit. Clearly, the model poorly fits the data; the residuals are dominated by the cubic term reported by Backer (1993). That the residuals are cubic in form-rather than linear or quadratic-should not surprise us, since the model that has already been removed from the data contains a second order polynomial: where T 0 is an arbitrary time chosen near the center of the data span.
FIG. 1.-In the top panel, the best fit timing model containing only a fit for f ,ḟ , the astrometric parameters, and a single Keplerian orbit has been removed from the data. A large cubic term dominates the postfit timing residuals. In the bottom panel, the timing model given in Table 1 has been removed from the data. Note the difference in scales between the panels.
Qualitatively similar timing behavior is familiar from studies of young pulsars (e.g., Downs and Reichley 1983). So-called "red" timing noise is thought to arise from stochastic interactions between the crust of the neutron star and the superfluid vortices that carry angular momentum in the interior of the star. There are strong arguments against such an interpretation in the case of PSR B1620−26. First, our understanding of timing noise suggests that it is inversely correlated with the pulsar age: old, relatively cold millisecond pulsars have relatively little phase wander (Stinebring et al. 1990, Cordes 1993, Arzoumanian et al. 1994. Second, the magnitude of the observed cubic is very large for any pulsar, independent of age. During the eleven years of observations, the pulsar's apparent spin-down rate has varied from −8.1 × 10 −15 s −2 to −1.4 × 10 −15 s −2 , changing the apparent torque by a factor of more than five; at the current rate of change, the spin-down of the pulsar will change to an apparent spin-up in November 2000. While we cannot rule out intrinsic spin instabilities as the cause of the observed signal, such instabilities are, in magnitude, unlike any seen in any other pulsar, and have no known theoretical basis. Turning to external causes, the most promising explanation of the largef is a changing gravitational acceleration ("jerk") by a massive body external to the binary. As we have previously discussed (Thorsett, Arzoumanian, and Taylor 1993), the most likely candidate is a second gravitationallybound companion. Although several pulsars in globular clusters are observed to have positive apparent frequency derivatives due to acceleration in the mean cluster field (Wolszczan et al. 1989, Anderson 1992, Nice and Thorsett 1992, the central mass density and velocity dispersion in M4 is far too low to produce the observed jerk on the PSR B1620−26 system (Phinney 1993). A close passage by a cluster star on a hyperbolic orbit cannot be ruled out, but the chance of discovering the binary during such a transient event is very small (∼ 2 × 10 −5 , Phinney 1993). Timing measurements could eventually distinguish between ellipsoidal and hyperbolic orbits for a second companion, but we will not further consider the possibility here.

A Polynomial Model
Obviously, a theory-independent way to characterize the deviations of the observations from the standard timing model is to add additional terms to the polynomial for φ(t), essentially Taylor-expanding the residuals around the epoch where we have introduced the convenient notation f (n) ≡ d n f /dt n . Given a model for the system, such as acceleration of the binary in the gravitational field of a second companion, the measured f (n) can be related to the model parameters. This is the approach that has previously been used in analysis of PSR B1620−26 data Taylor 1993, Joshi andRasio 1997).
The binary pulsar model can be further generalized to allow secular or periodic variations in the orbital elements. The simplest extension is to allow the elements to vary linearly with time: for example, ω = ω 0 +ω(t − T 0 ), where both ω 0 andω are free parameters in the least-squares fit. In the limit of a Keplerian orbit, of course, we expect the time derivatives of the elements to vanish.
In Table 1, we present the timing parameters obtained from a fit to the PSR B1620−26 data using a model including frequency derivatives to f (5) and allowing for linear variations in the orbital elements. In general, the results are in excellent agreement with those we have previously published, with substantially smaller uncertainties. (Note that, in order to reduce parameter covariances, the epoch of the frequency expansion has been changed from that of Thorsett et al. (1993), to a point nearer the center of the current data set.) A significant measurement of the pulsar proper motion has been made for the first time; it can be compared with optical measurements of the cluster proper motion (Cudworth and Hansen 1993): µ RA = −11.6(7) and µ dec = −15.7(7)mas yr −1 . The measurements disagree at ∼ 95% confidence-the difference is −1.8 ± 1.2 mas yr −1 in right ascension and −9.3 ± 5 mas yr −1 in declination. At d = 1.7 kpc (Peterson, Rees, and Cudworth 1995), the difference corresponds to a pulsar velocity relative to the cluster of 78 ± 40 km s −1 -far above the cluster escape velocity. Although such a velocity could be imparted to the system either during the formation of the triple or in an interaction with another cluster star, the probability of discovering such a system on an escape trajectory is negligible, since it would spend only ∼ 10 4 yr as close to the cluster core as PSR B1620−26 is now observed. More likely, either the pulsar or cluster proper motion is in error. Further timing observations are needed to test this claim. (A chance angular coincidence of an unassociated pulsar and cluster can be ruled out by the projected position very near the cluster core, the similar distance estimates, and the relatively good proper motion agreement.) The observed rate of change of the dispersion measure, (dDM/dt)/DM ≈ 8 × 10 −6 yr −1 , is comparable to the 1.4 ×10 −5 yr −1 change observed in the millisecond pulsar PSR B1937+21, which is at a similar dispersion measure (Ryba 1991).
The only significant deviation from simple Keplerian orbital motion is a very large derivative of the projected semi-major axis x, with timescale x 1a /ẋ 1a ≈ 3 Myr. The measurement is significant at > 20σ. Only lower limits are available for the timescales of changes in the other elements, and the limits are relatively weak in comparison with the measured timescale of change in x a1 : P a /Ṗ a ∼ > 2.2 Myr, e a /ė a ∼ > 0.6 Myr, and |1 rad/ω 1a | ∼ > 0.4 Myr. We defer discussion of the orbital perturbations until §3.3. We wish to relate the observed frequency derivatives to the properties of the triple system. To a good first approximation, we can treat the motion in the system as the sum of two non-interacting Keplerian  (14)  NOTE.-Timing parameters relative to a model including a Keplerian orbit with elements allowed to vary linearly with time, and a Taylor expansion in the barycentric pulsar frequency extended through f (5) . Position is relative to the JPL DE405 solar system ephemeris. Numbers in parentheses are uncertainties in the last digits shown. See §2 for discussion of notation. ellipses. The pulsar (with mass m 1 ) orbits with the inner companion (mass m 2 ) around their common center of mass. The elements of the pulsar's orbit are x 1a , P a , e a , ω 1a , and T a . The inner binary (mass m 1 + m 2 ) then orbits with the outer companion (mass m 3 ) around the center of mass of the triple. The elements of the pulsar's motion in the second orbit are x 1b , P b , e b , ω 1b , and T b . These latter parameters can be related to the line-of-sight acceleration a ·n wherê n is a unit vector along the line of sight, which can in turn be related to the measured frequency derivatives:ḟ (1) where in each case we have neglected terms of order v/c smaller than the leading contribution. (As will be seen below, the reflex motion of the binary in its orbit with the third companion has an amplitude ∼ 6 m s −1 ∼ 2 × 10 −8 c.) An explicit expansion of a = a(x 1b , P b , e b , ω 1b , T b ) is tedious, but has been carried out in some special cases by Joshi and Rasio (1997). Substitution of the result into eqn. 1 yields a series of non-linear equations in five unknown parameters. Measurement of the five frequency derivativesḟ ,. . . , f (5) then gives a set of equations that can be inverted (at least numerically) to determine the orbital elements.
In fact, onlyḟ ,. . . , f (4) have so far been measured in the case of PSR B1620−26, with only a limit for f (5) ( Table 1). Hence inversion of eqn. 1 will yield a one parameter family of solutions. The situation is further complicated by the fact thatḟ includes an unknown contribution from the intrinsic pulsar spindown rate, soḟ =ḟ int +ḟ acc . Becauseḟ int is nearly constant, whileḟ has changed by a factor of five during the time period described, it is likely thatḟ ∼ḟ acc . It is, on the other hand, unlikely thatḟ ∼ḟ int , since such a large intrinsic torque implies a dipole field strength ∼ 3 × 10 9 G, an order of magnitude larger than the typical field for millisecond pulsars in the plane of the Galaxy (Camilo, Thorsett, and Kulkarni 1994). Although we recognize the possibility that fields of pulsars in the plane and in clusters might differ systematically, we assume a magnetic field strength of 3 × 10 8 G, implying thatḟ int ∼ 10 −2ḟ , orḟ ≈ḟ acc . Intrinsic contributions to higher order derivatives of f are expected to be negligible.
Using the method outlined by Joshi and Rasio (1997), we have inverted eqn. 1, letting the eccentricity e b range over discrete values between 0 and 1. For each value of e b , we have used the resulting orbital parameters to calculate f (5) , and compared it with the measured limit. The results are shown in Fig. 2. Orbits with eccentricity below e b = 0.17 are excluded because they predict a larger f (5) than observed. To test the sensitivity of our results to the assumptioṅ f acc =ḟ obs , we have repeated the calculation assumingḟ acc = 0.1ḟ obs . As seen in the figure, relaxing this assumption makes virtually no change in the inferred properties of the companion or its orbit (except reducing the mass m 3 by about 10%). However, orbits with smaller eccentricity (including circular orbits) are allowed. We will return to this point below.  Joshi and Rasio (1997). The inner binary mass was assumed to be m 1 + m 2 = 1.7M ⊙ . The bold line is the solution assuming that the observedḟ = −5.5 × 10 −15 s −2 on JD2448725.5 has negligible contribution from the intrinsic spin-down. The lighter line assumes that 90% of the spin-down rate on that date was intrinsic. These curves bound the likely region.

Double Keplerian Models
Within a simple Keplerian model, there are only five orbital elements accessible from line-of-sight velocity measurements, so only five coefficients in the polynomial expansion are independent. Continuation of the expansion beyond five terms could thus, in principle, provide a test of the triple hypothesis. Unfortunately, the current data span is inadequate to provide high-significance estimates of the terms beyond f (5) in the expansion. Even with the current data set, a polynomial model is not the optimal approach to estimating orbital parameters. For example, high order polynomial terms are significantly covariant (particularly odd terms with odd, and even terms with even). Also, the non-uniformity of the data, both in sampling intervals and signal-to-noise ratio achieved, makes it difficult to interpret the fit over the extended data set as a simple Taylor expansion around a particular epoch. To avoid these complexities, we introduce a timing model in which the pulsar orbits in a hierarchical triple system: it is assumed to move in a Keplerian orbit around its common center of mass with the inner companion (possibly with elements that vary linearly with time), and the inner binary then orbits about the center of mass of the triple in a second Keplerian elipse. Not surprisingly, the second orbit is underdetermined (ignoring for the moment the evolution of the elements of the inner orbit). We approach this limitation in several ways.
First, we assume the outer orbit is circular. The resulting orbital parameters are given in Table 2. If we assume that the mass of the inner binary is 1.7M ⊙ , then the mass function gives a "projected" mass of the second companion m 3 sin i = 1.2 ×10 −3 M ⊙ (about 20% more massive than Jupiter) and a separation of 19 AU (comparable to the size of the orbit of Uranus). Comparing Tables 1 and 2, we note that in the model with a circular outer orbit, acceleration contributes just over 10% of the observed spin-down at the epoch. Our result can thus be directly compared to the lighter curve in Fig. 2. We note the excellent agreement (within the published errors) between predictions from the polynomial model and the direct fit to the double-Keplerian orbit.
The circular orbit solution yields a relatively large intrinsic spin-down rate, and hence a large inferred dipole field strength of 2.6 × 10 9 G. As discussed above, this is relatively high for a millisecond pulsar, though well below the ∼ 8 × 10 9 G upper limit for an 11 ms pulsar on the spin-up line (Lyne and Smith 1998). The age of the pulsar, assuming it was born on the spin-up line and slowed with a braking index n = 3, is about 240 Myr-far smaller than the cluster age. If the neutron star was formed in the collapse of a massive star during the early life of the cluster, its spin-up to form a millisecond pulsar would be a relatively recent event. At first glance, it might be surprising to note that in this model the ascending node passage occurred during our limited data span, in April 1993. However, with eleven years of data and a 62 yr orbit the chances were actually better than one in three that our observations would include either an ascending or descending node passage.
As an alternate approach to the underdetermination of the outer orbital parameters, we could assume that the magnetic field of the pulsar is small (e.g., ∼ 3 × 10 8 , a typical value for fast pulsars (Camilo, Thorsett, and Kulkarni 1994)). Then the intrinsicḟ contributes negligibly to the observed value, so we take it to be zero. As with the polynomial fits described above, we are then able to produce a one parameter family of solutions, where it is convenient to take the free parameter as the eccentricity. We present orbital parameters for two values of e in Table 3. In each case, the results from the timing fits agree to better than 10% (and generally within 5%) with the parameters obtained from the polynomial fitting procedure in §3.1.

Orbital perturbations
To this point, we have neglected three-body effects. In a hierarchical triple like the PSR B1620−26 system, the orbits are at all times very nearly Keplerian ellipses. It is customary to consider the osculating orbital elements which are, at any moment, the Keplerian elements of the orbit tangent to the real orbit with the same velocity. Various algebraic and numeric techniques have been developed to determine the time evolution of these osculating elements (e.g., Brouwer and Clemence 1961). We note that variations in the projected orbital elements can also arise from proper motion, but such effects are expected to be small in this case (Arzoumanian et al. 1996).
The perturbations can be usefully divided into short-period (∼ P a ) terms, long-period (∼ P b ) terms, and "apse-node" terms (∼ P 2 b /P a ) terms (Brown 1936, Söderhjelm 1975. The amplitudes of the short-period terms are too small to be detected in the current data. Because the data span is much shorter than either the long-period or apse-node timescales, the observed "secular" perturbations are a sum of contributions from both terms. The general solution is quite complex, but can be greatly simplified if we assume that P b is much longer than the data span, so we can take the companion to be at a fixed position during the observations. Given the relatively crude perturbation measurements now available, this assumption is adequate for even the smallest allowed values of P b . Projected semi-major axis x (s) 6.4(2) Orbital period P b (yrs) 61.8(7) Time of ascending node T 0 (JD) 2449104(5) Mass function (M ⊙ ) 5.6(4) × 10 −10 a Intrinsic spin frequency derivatives beyond the first are assumed to vanish. Shown are spin parameters and parameters of the outer orbit; other parameters are consistent with those given in Table 1.
The orbital perturbations were first calculated in this approximation by Rasio (1994). They can be written (Rasio 1994, Joshi andRasio 1997) ω 1a = 3πη P a sin 2 θ 3 5 cos 2 φ 3 − 1 − 1 (2) e a = − 15πη 2P a e a sin 2 θ 3 sin 2φ 3 (3) where η = [m 3 /(m 1 + m 2 )](a a /r 3 ) 3 , a a = a 1a + a 2a is the semimajor axis of the inner binary, and r 3 , θ 3 , and φ 3 are the fixed spherical polar coordinates of the second companion with respect to the center of mass of the inner binary, choosing φ 3 measured from pericenter in the orbital plane and θ 3 such that sin θ 3 = 1 in the coplanar case. The "secular" perturbation of the projected semimajor axis is thusẋ 1a = x 1a cot i aia . To use our measurementsẋ 1a and constraints oṅ ω 1a andė a to further constrain the triple system parameters, we have performed Monte Carlo simulations following the procedure outlined in Joshi and Rasio (1997). We have assumed a uniform prior probability distribution for cos i a , cos i b , and α -the angle between the lines of nodes of the two orbits. Whereas Joshi and Rasio assumed a thermal distribution in e b (prior probability proportional to e b ), without any reason to expect that the system is in thermal equilibrium with the cluster, we have instead adopted a uniform distribution in e b , though we note that this choice has very little effect on our results.
Our procedure is straightforward. For each trial, we select values for cos i a , cos i b , α, and e b as described. Using e b , we calculate x 1b , P b , ω 1b , and T b as described in §3.1 (assumingḟ ≈ḟ acc ). Using cos i a , we find m 2 , and using cos i b and α we calculate r 3 , θ 3 , φ 3 , and η. (We assume a neutron star mass m 1 = 1.4M ⊙ .) Using eqns. 2-4, we calculate the expected orbital perturbations, which we compare with the measured values ofẋ 1a ,ω 1a , andė a and their uncertainties, rejecting trials with a probability determined from the appropriate three-dimensional gaussian distribution. 129 (2)  Some results of our analysis are shown in Figs. 3-6. Our results are in good agreement with the earlier analysis of Joshi and Rasio (1997), which were based on a preliminary version of the timing results published here. We find that the mass of the outer companion is quite tightly constrained: m 3 = 0.0118 +0.0087 −0.0048 (68%) or m 3 = 0.0118 +0.0373 −0.0073 (95%). The current distance of the outer companion from the center of mass of the binary is also well constrained, r 3 = 35 ± 6 AU (68%) or r 3 = 35 ± 10 AU (95%). Less well constrained is the orbital period P b . The most favored values are just above the minimum allowed from the f (n) measurements, or a few hundred years, but there is a long tail to very long orbital periods. These solutions correspond to very high eccentricity orbits in which the second companion is currently very near periastron. The 68% confidence upper limit to the orbital period is 1200 yrs. The inclination of the outer orbit is not well constrained by the data, with a posterior likelihood very close to the assumed cos i prior likelihood. However, the inner inclination is constrained to 40 ± 12 • (68%) or 40 ± 24 • (95%). This leads to a most likely inner companion mass of m 2 = 0.46M ⊙ , somewhat higher that the 0.3M ⊙ normally assumed.  6.-Results of Monte Carlo estimation of the posterior likelihood distribution for the inclinations of the inner and outer "binary" orbits, as discussed in the text. The inclination of the outer orbit is not strongly constrained (the prior distribution was uniform in cos i b ), but the inner orbit is constrained to a fairly low inclination, i a ∼ 40 • .
We have also confirmed that the orbital perturbation measurements and limits are consistent with the solution in whichḟ acc = 0.1ḟ obs . In particular, solutions with a circular outer orbit (Table 2) are acceptable if i a ∼ 40 • and α ∼ 100 • or α ∼ 280 • . Again, i b is poorly constrained, and m 3 = 1.55 ± 0.30 × 10 −3 M ⊙ , or about half again the mass of Jupiter.

CONCLUSIONS
Since the initial suggestion that pulsar PSR B1620−26 is the member of a triple system, constraints on the properties of the system components have gradually improved. The detection of a secular change in the projected semimajor axis of the inner binary, which has been discussed previously at conferences (e.g., Arzoumanian et al. 1996) but is presented in detail here for the first time, is an important confirmation of the triple nature of the system, since there is no known way for spin instabilities of the pulsar to produce timing fluctuations on the orbital timescale. The only simple way to understand the system is as a binary accelerating in an external tidal field, most likely due to a bound companion.
As we have shown, the orbital parameters of the outer orbit are still poorly constrained. The reason is not hard to understand. Because the span of available data is much shorter than the period of the outer orbit, the relative position of the inner binary and outer companion has not changed substantially since the pulsar was discovered. Although the position itself is known rather accurately, the relative velocity is not. This leads to a family of allowed solutions, all with periastron distance ∼ 35 AU (the current separation). Very high eccentricity orbits with very large semimajor axis can account for the observed acceleration and orbital perturbations, but only if the system is currently observed near periastron. Because a body moving in a high eccentricity orbit spends only a small fraction of the time near periastron, such solutions require significant fine tuning.
Another remaining source of uncertainty is the extent to which the observed frequency derivativeḟ obs is dominated by acceleration, rather than intrinsic pulsar spin-down. We have argued that it is most likely thatḟ acc ≈ḟ obs , in which case m 3 ∼ 0.01M ⊙ . Allowing a significant contribution fromḟ int reduces the magnitude of the acceleration by the third body, allowing masses as small as m 3 ∼ 10 −3 M ⊙ .
In either case, the second companion is almost certain substellar: below the ∼ 0.08M ⊙ hydrogen burning limit, and is most likely below or near the deuterium burning limit, ∼ 0.015M ⊙ (Leibert and Probst 1987). Whether it is called a "brown dwarf" or a "planet" is probably not important. The possibility that such objects might be found in globular clusters was first suggested just prior to the discovery of the triple nature of PSR B1620−26 (Sigurdsson 1992), and Sigurdsson (1993 has further suggested models for the formation of such a triple involving exchange interactions in the dense environment of M4. However, developing a complete model that can explain not only the system formation and stability and the pulsar spin-up, but also the non-zero eccentricity of the inner binary remains an open problem (Joshi and Rasio 1997).
The prospects for continued improvement of the timing constraints on the system parameters are good. We expect (eqns. 2-4) that the timescales for perturbation of ω 1a and e a should be comparable to that of x 1a , and as noted in §3.1 the current measurement uncertainties have nearly reached that level. Furthermore, the limits on f (5) already significantly constrain the allowed solutions, and we expect the uncertainty to improve rapidly with observing span.
We thank D. Backer, R. Foster, and J. Taylor for their substantial contributions to the early stages of this project, and an anonymous referee for helpful comments. M. McKinnon, T. Hankins, and M. Goss provided important support for pulsar timing at the VLA. We particularly thank K. Joshi, F. Rasio, and S. Sigurdsson for many interesting and helpful discussions over the years. Green Bank and the VLA are facilities of the National Radio Astronomy Observatory, operated by Associated Universities, Inc., under contract from the NSF, which has also provided direct support for this work-as have Caltech, Princeton, and the Sloan Foundation.