{'fraction_non_alphanumeric': 0.04413977595719779, 'fraction_numerical': 0.034370745455848284, 'mean_word_length': 4.2584777694046725, 'pattern_counts': {'":': 0, '<': 24, '<?xml version=': 0, '>': 5, '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}}
Evidence for a Massive Neutron Star from a Radial-Velocity Study of the Companion to the Black Widow Pulsar PSR B1957+20

The most massive neutron stars constrain the behavior of ultra-dense matter, with larger masses possible only for increasingly stiff equations of state. Here, we present evidence that the black widow pulsar, PSR B1957+20, has a high mass. We took spectra of its strongly irradiated companion and found an observed radial-velocity amplitude of K_obs=324+/-3 km/s. Correcting this for the fact that, due to the irradiation, the center of light lies inward relative to the center of mass, we infer a true radial-velocity amplitude of K_2=353+/-4 km/s and a mass ratio q=M_PSR/M_2=69.2+/-0.8. Combined with the inclination i=65+/-2 deg inferred from models of the lightcurve, our best-fit pulsar mass is M_PSR=2.40+/-0.12 M_sun. We discuss possible systematic uncertainties, in particular in the lightcurve modeling. Taking an upper limit of i<85 deg based on the absence of radio eclipses at high frequency, combined with a conservative lower-limit to the motion of the center of mass, K_2>343 km/s (q>67.3), we infer a lower limit to the pulsar mass of M_PSR>1.66 M_sun.


INTRODUCTION
One of the outstanding problems in physics is the behavior of matter at extreme densities. This behavior is modelled using quantum-chromodynamics calculations, but these cannot yet determine reliably the densities at which, e.g., meson condensation and the hadron to quark-gluon phase transition occur. At densities slightly above nuclear and high temperature, models can be tested with heavy-nuclei collision experiments. For higher densities and low temperatures, only comparison with neutron-star parameters is possible (for a review, e.g., Lattimer & Prakash 2007).
The different models lead to different equations of state, which predict different mass-radius relations for neutron stars. Unfortunately, most attempts at observational tests have been frustrated by susceptibility to systematic errors and modeling uncertainties. The most robust tests have involved measurements of extrema. For instance, the fastest measured spin period, 1.4 ms (Ter 5ad; Hessels et al. 2006), excludes the stiffest equations of state, for which neutron stars would be too large to spin so fast.
A problem with measurements of the extrema, is that whether they occur in Nature depends not only on whether they are allowed physically, but also whether they are expected astronomically. From models of stellar evolution, Timmes et al. (1996) find neutron-star mass distributions at birth with two narrow peaks, at 1.3 and 1.8 M ⊙ , containing remnants of stars with initial masses smaller and larger than ∼ 19 M ⊙ , respectively. In binaries, however, where much of the stellar envelope is removed during the evolution, they expect only to form neutron stars in the lower mass bin. The above may explain why until recently most accurate masses were all close to 1.4 M ⊙ : most were measured for binaries containing pulsars with neutron star companions, where the preceding evolution predicts relatively little mass has been accreted (for a review, Stairs 2004). The one exception was the X-ray binary Vela X-1, for which a higher mass of 1.86 ± 0.16 M ⊙ was inferred (Barziv et al. 2001;Quaintrell et al. 2003). Such a high mass would imply, e.g., that meson condensation is not important. However, the large uncertainty prevented a definitive conclusion.
However, neutron stars can become more massive after birth by accretion. Accretion also leads to "recycling" of radio pulsars: it increases their spin frequency and decreases their magnetic fields. Most recycled pulsars are accompanied by low-mass, 0.2 M ⊙ companions. For these systems to form in the age of the Universe, the companion must originally have been a star of 0.8 M ⊙ . Thus, the companion lost almost all of its mass, and some of it should have landed on the neutron star, leading to a concomitant increase in the neutronstar mass. Initial tests, however, were suggestive but not conclusive: for pulsars with low-mass white dwarf companions, masses around 1.4-1.7 M ⊙ were found (van Kerkwijk et al. 1996;Jacoby et al. 2005;Bassa et al. 2006;Verbiest et al. 2008), somewhat larger than typical, but not yet very constraining. Presumably, in these systems much of the mass actually left the binary.
More recently, higher masses were found in different types of binaries, starting with a number of pulsar binaries in globular clusters (Ransom et al. 2005;Freire et al. 2008a,b). In these cases, however, the masses rely on observations of periastron advance, which is assumed to be due to general relativistic effects only (rather than classical ones such as due to rotationally and tidally induced quadrupoles), and statistical arguments that the inclinations are unlikely to be very low. Thus, it was still possible to doubt that very massive neutron stars could exist.
While we were writing and revising this work, however, such doubts disappeared, with accurate mass determinations for PSR J1614-2230 (1.97 ± 0.04 M ⊙ , Demorest et al. 2010) and PSR J1903+0327 (1.67 ± 0.02 M ⊙ ), both relying on measurements of Shapiro delay, which is not easily mimicked by other processes. These masses exclude many of the soft equa-tions of state, such as Kaon condensation as envisaged by Brown & Bethe (1994).
Intriguingly, both of the above systems do not have low-mass white-dwarf companions like most binary pulsars, but rather a more massive, carbon-oxygen white dwarf (PSR J1614-2230) and a solar-mass main sequence star (PSR J1903+0327). Thus, both systems also had different evolutionary histories (puzzling for PSR J1903+0327; see Freire et al. 2010). In this sense, our approach is similar, in that we specifically target another group of binary pulsars with different properties and evolutionary histories, and, therefore, perhaps different masses, the so-called "black widow pulsars." In black-widow systems, a millisecond pulsar is accompanied by a low-mass, few 0.01 M ⊙ companion, which is bloated and strongly irradiated by the pulsar, leading to outflows strong enough to eclipse the pulsar signal for significant fractions of the orbit. The irradiation causes strong heating on the side of the companion facing the pulsar, and, as a result, strong orbital brightness variations of the optical counterparts van Paradijs et al. 1988;Stappers et al. 1996). From detailed modelling of the lightcurves, the inclinations can be constrained (Callanan et al. 1995;Stappers et al. 1999;Reynolds et al. 2007), which, when combined with velocity information, can be used to derive masses.
Here, we present a radial-velocity study of PSR B1957+20 , the proto-type black-widow system, which has the brightest and best-studied counterpart. In § 2, we describe our observations and data reduction, and in § 3 we constrain the properties of the companion from the spectra, determine radial velocities, and fit an observed radial-velocity amplitude. In § 4, we discuss the available constraints on the radius and inclination, and in § 5 the corrections we need to make because our velocity amplitude is that of the center of light, which, because of the irradiation, is shifted towards the pulsar relative to the center of mass. We present our final constraints on the masses in § 6 and discuss how these may be made more secure in the future.

OBSERVATIONS & REDUCTION
We obtained pilot observations at the Keck telescope of the companion of PSR B1957+20 on the night of 15 June 2007, and a larger set of spectra on the nights of 4 and 5 August 2008 (see Table 1). For all observations, the seeing was good, ranging from 0.6 to ∼ 1 ′′ , but the sky was not photometric. For relative flux calibration, we obtained spectra of a number of spectrophotometric standards from Bohlin et al. (1995) (see below). We obtained exposures of internal flat fields and Hg/Kr/Ar and Ne/Ar arc spectra interspersed with the observations.
The spectra were obtained using the two-armed Low Resolution Imaging Spectrometer (LRIS, Oke et al. 1995;McCarthy et al. 1998). We employed the atmospheric dispersion corrector, used a 0. ′′ 7 slit, set to position angle 35 • to cover both the companion and a nearby star (at ∼ 1. ′′ 3, hereafter the contaminator 1 ), and split the light with a dichroic at 6800 Å. In the blue arm, we used the 600 line mm −1 grism, blazed at 4000 Å, which covers 3100-5600Å at a resolution ∆λ ≃ 3.2 Å or ∆v ≃ 220 km s −1 (for the 0. ′′ 7 slit). The detector is a mosaic of two Marconi CCDs, each with 4096 × 2048 pix-  (Arzoumanian, Fruchter, & Taylor 1994; we ignored orbital-period derivatives, see text). Col. (6): Flux ratio relative to the contaminator in the 5100-5300 Å range. Col. (7): Flux ratio in the 8450-8650 Å range. Col. (8): Radial velocity relative to the contaminator inferred from the blue spectra (we estimate a barycentric velocity of the contaminator of −25.6 ± 1.3 km s −1 ). els of 15 µm on the side (0. ′′ 135 on the sky), which we binned by two in the dispersion direction (the only direction it can be binned). On the red side, we used the 1200 line mm −1 grating, blazed at 7500 Å, set to cover 7600-8900 Å at ∆λ = 2.1 Å or ∆v ≃ 75 km s −1 . Here, the detector was a Tektronix CCD with 2048 × 2048 pixels of 24 µm on the side (0. ′′ 215 on the sky), which we read out unbinned.
For the reduction, we used the Munich Image Data Analysis System ESO-MIDAS, and routines running in the MIDAS environment. For all images, we subtracted bias as determined from the overscan regions. For the blue images, we subsequently corrected for small-scale variations in efficiency by dividing by a spatially averaged flat field, normalized using a third-degree polynomial, and with the bluest, poorly exposed part shortward of 4000 Å replaced by unity. For the red images, we simply divided by the flat field, normalized using a bi-linear fit.
The spectra and their uncertainties were extracted by fitting, at each dispersion position, the sum of three stellar profiles (for the companion, contaminator, and another star at ∼ 3. ′′ 2 on the other side of the target) and a constant sky. For the profiles, we used Moffat functions of the form P = A/(1+(x/w) 2 ) δ , with power δ = 5 for the blue images and δ = 6 for the red ones. At each dispersion position, only the three amplitudes A and the sky were fitted; the other parameters were determined globally (with the central position and width w allowed to vary quadratically with dispersion position, and the relative positions fitted as constants). The fits to the images were generally good for the pulsar fields, with reduced χ 2 near unity, and somewhat poorer for the higher signalto-noise images of the spectrophotometric standards (where slight mismatches between the Moffat function and the true point spread function are more apparent).
Wavelength calibration was done using arc spectra. For the blue arm, we used well-exposed images of the arc lamps taken at the start of the night to define an overall solution, which re-quired a fourth-order polynomial to give an adequate dispersion solution, with root-mean-square residuals of 0.14 Å for 24 lines. Next, we found offsets relative to this solution from less well-exposed arc frames taken interspersed between the observations. For the red arm, the individual arc frames were well exposed, and a third-degree polynomial sufficed to give solutions with residuals of ∼ 0.04 Å (for typically 31 lines).
For flux calibration, we first corrected all spectra approximately for atmospheric extinction using a curve made by combining the CFHT values (Bèland, Boulade, & Davidge 1988) shortward of 5200 Å with the better sampled La Silla values longward of 5200 Å (ESO users manual 1993; see also Tüg 1977). Next, for the blue spectra, we calculated response curves by comparing our observed spectra for Feige 110 with the calibrated STIS spectra given by Bohlin et al. (1995): we slightly smoothed our spectra to match the STIS resolution, divided the two, and smoothly interpolated the ratio. For the 2007 red spectra, we proceeded similarly, except that we compared our observation of BD +28 4211 with a third-degree polynomial fit to the STIS spectrum of Bohlin et al. (1995), since the STIS spectrum was relatively noisy and the intrinsic spectrum should be smooth. Furthermore, we took care to fit the ratio spectrum avoiding telluric absorption, so that we could use the deviations -scaled with airmass -to correct to first order for telluric absorption in our target spectra. For the 2008 red spectra, our procedure was similar, except that we compared our observed spectrum of Hz 43A with the calibrated model spectrum (which is again smooth over the relevant wavelength range; Bohlin et al. 1995). Since the conditions were not photometric, the above gives good relative fluxes over the observed wavelength range, but only approximate absolute fluxes. To place the companion spectra on the same flux scale, we scaled all sets of spectra to the same average contaminator flux in the 5100-5300 Å and 8450-8650 Å ranges (for blue and red, respectively). After this scaling, integrating over B filter band passes, we find that our inferred companion magnitudes are roughly consistent with the ones observed by Reynolds et al. (2007).

SPECTRA AND VELOCITIES
In Fig. 1, we show averaged blue and red spectra for phases within 0.15 from superior conjunction of the companion, when the irradiated side is seen face on, and for phases between 0.15 and 0.3 from superior conjunction, when the view is more sideways.
Below, we first determine the spectral type for the two phases. While this will be only an "average" spectral type over the surface, it gives a way to determine the influence of irradiation. We will find that the spectrum appears rather normal, suggesting that the irradiated energy is dissipated well below the photosphere (in contrast to the case of, e.g., the pre-cataclysmic variable NN Ser, where the companion is irradiated with ultraviolet light from a hot white dwarf; see Parsons et al. 2010).
Given the normal appearance of the spectra, we use them to help determine reddening and distance, and to constrain the surface gravity and radius of the companion. In § 5, we will find that these provide additional evidence that the companion is close to filling its Roche lobe. Next, we describe how we determined the velocities from the individual spectra, and use these to fit an orbit. We discuss the required corrections to this orbit in § 5.
3.1. Spectral Type, Reddening, and Distance We compared our blue spectra with classification spectra shown in the on-line atlas by R. O. Gray. 2 We find that the contaminator has spectral type G1-2 V, while that of the companion is slightly earlier at its brightest phase (F9-G0, as seen from, e.g., the weaker G band near 4300 Å), and similar when seen from the side (about G1). Its luminosity class is slightly higher, intermediate between IV and III, as can be seen from the stronger Sr II λ4077 line, and also supported by the lower strength of the Fe I λ4046 line relative to the Mn I λ4030 line. Consistently, in Fig. 1, one sees that β Hyi, which has spectral type G2 IV, is intermediate in luminosity class between the contaminator and the companion. For β Hyi, the mean density and radius have been measured using asteroseismology and interferometry, leading to a surface gravity log g = 3.952 ± 0.005 (North et al. 2007). For the companion of PSR B1957+20, we thus infer log g 4.
The blue spectra of the companion show emission cores in the Ca II H and K lines, suggestive of an active chromosphere. Furthermore, Hβ seems rather weak, possibly due to being filled in by poorly subtracted emission from the Balmerdominated bow shock nebula (Kulkarni & Hester 1988; the contaminator might also be affected by this). Finally, the Mg Ib triplet is a bit weaker than expected.
In the red spectra, we find that the spectrum of the contaminator is as expected for a G1-2 V star, but that the spectra of the companion cannot be classified as easily. In particular, the Ca II IR triplet nearly absent, perhaps being filled in by chromospheric emission. Furthermore, the Na I λλ8183,8195 lines are much stronger than seen in the contaminator: combined equivalenth widths of ∼ 0.9 and 1.5 Å for the face-on and more sideways view of the companion, respectively, and ∼0.7 Å for the contaminator. Probably, this reflects that in the red, a large contribution to the light arises from cooler regions of the companion, which will have much stronger Na I absorption (e.g., Zhou 1991 finds equivalent widths of ∼ 0.7 Å and ∼ 1.5 Å for early G and early M dwarfs, respectively). The strength of the Na I line also shows that the surface gravity of the companion is similar to that of a main-sequence or sub-giant star; for giants, the Na I equivalent width does not exceed 1 Å (Zhou 1991) for any temperature.
We can use the spectral information to constrain the reddening to the source. Starting with the contaminator, from archival Hubble Space Telescope observations with WFPC2, which we analyzed with HSTphot (Dolphin 2000), we measure m F675W = 19.968 ± 0.010 and m F814W = 19.458 ± 0.015. Using the transformations of Holtzman et al. (1995), this corresponds to R = 20.06 ± 0.03 and I = 19.40 ± 0.02. Since a G1-2 V star has (R − I) 0 = 0.34 ± 0.02 and M R ≃ 4.3 (Cox 2000), we infer an extinction A V = 1.42 ± 0.18 mag and distance of ∼ 8 kpc (where we use the extinction curve of Schlegel, Finkbeiner, & Davis 1998).
The reddening to the companion to PSR B1957+20 should be similar to that to the contaminator, since its observed colors at maximum are about as much bluer (from our blue spectra, ∆(B − V ) = −0.07 ± 0.02) as the intrinsic color difference expected from the spectral type (∆(B −V ) 0 = −0.05 ± 0.02). The reddening can be compared with the run of reddening with distance measured using red-clump stars. Using the technique of Durant & van Kerkwijk (2006), we find that the reddening FIG. 1.-Spectra of the companion to PSR B1957+20. We show averages both around superior conjunction (blue trace), when the illuminated side is most in view, and nearer quadrature (red trace). Also shown are the spectrum of the nearby "contaminator" (green trace, offset by 5 units), as well as a UVES spectrum of the G2 IV star β Hyi (Bagnulo et al. 2003) (black trace; convolved to the same resolution, scaled to match the brightness of the conjunction spectrum, and offset by 10 units; the red part of that spectrum is not shown, as it has residual ripples and does not cover the Ca II triplet). The right-hand enlargement shows the Ca II H and K lines, which have cores in emission for the companion, and lines of Sr II at 4077 Å, Fe I at 4046 Å, and Mn I at 4030 Å, which we use to estimate the luminosity class and surface gravity.
is A V ≃ 0.5 at d ≃ 1.5 kpc, increases first slowly and then more sharply to A V ≃ 1.2 at d ≃ 2 kpc, and remains constant thereafter. Given that the contaminator is certainly well beyond 2 kpc, and that the reddening of PSR B1957+20 and the contaminator are similar, we conclude that PSR B1957+20 is at a distance d 2 kpc. This is consistent with the distance in the range of 1.5 to 2.5 kpc inferred from the pulsar dispersion measure of 29 cm −3 pc using models of the Galactic electron distribution (Taylor & Cordes 1993;Cordes & Lazio 2002). From the observations of Callanan et al. (1995), the companion at maximum brightness is about 0.4 mag brighter than the contaminator in the R band, which, using the R-band magnitude above, 3 implies R = 19.7. This corresponds to an absolute magnitude M R ≃ 7.2 − 5 logd 2 and an effective radius

Radial Velocities
We determined velocities by fitting the flux-calibrated blue spectra of both the pulsar companion and the contaminator with a template based on the high-resolution spectrum of the G2 IV star β Hyi (HD 2151) from the UVESPOP library (Bagnulo et al. 2003). For our fitting, we convolved the UVES spectrum with a truncated Gaussian to match the resolution of the observed spectra corresponding to the ∼ 0. ′′ 7 seeing and 0. ′′ 7 slit. We ignored rotational and orbital broad-3 From Fig. 1 of Callanan et al. (1995), we read off R = 19.89 ± 0.03 for the contaminator, with a quoted systematic uncertainty of 0.05 mag. This is somewhat brighter than what we measured from the HST data. We use the fainter magnitude to obtain a conservative limit on the radius .   FIG. 2.-Radial-velocity measurements of the companion of PSR B1957+20 as a function of orbital phase. The colored circles mark velocity measurements from different nights, with errors as indicated. They are repeated in black for clarity. The drawn curve is the best-fit circular orbit, and the dotted one the pulsar orbit inferred from radio timing. ening, which are well below our resolution for the blue spectra and just comparable to that of the red spectra (for companion radius R c ≃ 0.25 R ⊙ , one finds v sin i = 2πR c sin i/P ≃ 30 km s −1 sin i; for integration time t int = 1800 s, the maximum smearing is K obs sin(2πt int /P) ≃ 80 km s −1 , where K obs = 324 km s −1 is the radial-velocity amplitude [see below]).
The template was fitted for a grid of velocities between −600 and +600 km s −1 with a step size of 5 km s −1 , at each velocity fitting for the normalization and possible variation with wavelength using a quadratic function. Typical reduced χ 2 values for the best fits were χ 2 red ≃ 1.2. The best fit velocity was determined using a quadratic fit to the χ 2 values within 40 km s −1 of the minimum. Looking at the results for the contaminator, it is clear that systematic variations are present, with root-mean-square of 13 km s −1 , likely due to small shifts in placement in the slit and/or uncorrected atmospheric dispersion. As our velocities, we thus take the difference between the velocities inferred for the companion and the contaminator; these velocities and their corresponding uncertainties are listed in Table 1.
To see how sensitive our results are to our choices, we tried fitting only part of the blue spectrum, and using different UVESPOP stars or model atmospheres as templates. We found that the velocities were consistent to within ∼1σ, which we will use as an estimate of the associated systematic uncertainty below. We also tried fitting the red spectra, but found that for most, the absence of the expected lines combined with the poor signal-to-noise did not allow us to obtain a reliable velocity.
We fitted the velocities with a circular orbit using epoch of ascending node T 0 = MJD 48196.0635242 and orbital period P = 33001.91484 s (Arzoumanian, Fruchter, & Taylor 1994), but found we could obtain a good fit only if we left the phase free (Fig. 2). Compared to the prediction, ascending node occurs ∼ 350 s later (phase offset 0.011 ± 0.002). However, the orbital period is known to vary quasi-periodically (Arzoumanian et al. 1994;Nice, Arzoumanian, & Thorsett 2000); indeed, from the first and second orbital period derivatives measured by Arzoumanian et al. (1994), the expected phase offset is 0.0406. Keeping the phase free, our fit has χ 2 red = 0.92 for 12 degrees of freedom (15 measurements, 3 parameters). The observed radial-velocity amplitude is K obs = 324 ± 3 km s −1 and the systemic velocity, measured relative to the contaminator, is γ obs = 110 ± 5 km s −1 (here, we multiplied the formal errors with √ 2 to account for the systematic uncertainty in the velocities related the choice of template and fitting region discussed above).
The above systemic velocity is relative to the velocity of the contaminator. We tried to measure the latter in two different ways. First, a simple average of the velocity measurements from the blue spectra yields −39 ± 4 km s −1 (where we corrected for the velocity of 22.7 ± 0.9 km s −1 of β Hyi, Evans 1967; we verified that we obtained the same result within 2 km s −1 using other UVESPOP stars and model atmospheres). Possible evidence that this is not reliable, however, comes from the large scatter in the velocities from individual spectra (see above). A similar scatter is seen in the O I 5577 sky emission line. The latter also shows an average offset, of 6 km s −1 , but it is not clear whether one can apply this, since the line is very close to the red edge of the blue spectra, where the wavelength solution may be less reliable.
As an alternative, therefore, we measured the velocities of the contaminator using the red spectra, where, given the higher resolution, any shifts due to offsets in the slit and/or atmospheric dispersion should be smaller. Here, we unfortunately could not use the UVESPOP spectra, since these have a gap over the Ca II IR triplet, and thus we determined velocities relative to a T eff = 6000 K, log g = 4.5 model spectrum from Zwitter et al. (2004). We find that the scatter in these velocities is substantially smaller (root-mean-square of 6 km s −1 ), and that the (weighted) average velocity is −25.1 ± 1.5 km s −1 . From the OH sky line at 8344.602 Å, we infer that any offset due to wavelength calibration errors are small; including these, we find an average velocity of −25.6 ± 1.3 km s −1 .
As a check on these numbers, we can compare the velocities with what is expected for a star that is following simple Galactic rotation. Assuming a flat rotation curve with Θ = 220 km s −1 and a distance to the Galactic center of 8.5 kpc (Cox 2000), as well as a peculiar velocity of the Sun relative to the local standard of rest of (U,V,W ) = (10.00, 5.25, 7.17) km s −1 (Dehnen & Binney 1998), we find that for distances of 6, 8, and 10 kpc for the contaminator, the expected radial velocities are +16.7, −1, and −24 km s −1 , consistent with our measurements. (Arguably, one should consider asymmetric drift in our estimates, in which case the velocities would be more negative by ∼ 15 km s −1 .) Overall, we believe the contaminator velocity from the red spectra is more reliable. For PSR B1957+20, the implied systemic velocity is 84 ± 5 km s −1 .

RADIUS AND INCLINATION
The companion is irradiated by the pulsar and presents a hot and a cold side. Hence, the center of light does not coincide with the center of mass, but is shifted somewhat towards the pulsar. As a result, the spectroscopic observations underestimate the true radial-velocity amplitude. The correction depends on the stellar radius, the temperature distribution, and the inclination. Here, we discuss the available constraints on these properties.
The lightcurve provides strong constraints on the system parameters, especially if simplifying assumptions can be made for the temperature distribution. For PSR B1957+20, high-quality lightcurves were presented by Callanan et al. (1995) and Reynolds et al. (2007). These authors also fit their lightcurves with models.
In these lightcurve synthesis models, it is assumed that the companion's shape is that of an equipotential surface, and that its temperature distribution is given by some background temperature (modified suitably by gravity darkening) that is increased by irradiation by an isotropic pulsar wind such that the outgoing flux equals the sum of the background and irradiation fluxes (for a detailed description, see Orosz & Hauschildt (2000), who wrote the ELC code used by Reynolds et al. (2007)). It is assumed the irradiation does not affect the temperature structure of the atmosphere (the "deep heating" approximation), such that each surface element can be taken to radiate as predicted by a model atmosphere for a single star. The main free parameters are the background temperature and the irradiating flux (which set the temperature distribution), the extent to which the companion fills its Roche lobe (which determines its shape), and the inclination of the orbit. Reynolds et al. (2007) find that such models reproduce the lightcurves in detail, and they infer that the companion nearly fills its Roche lobe, up to a filling factor R nose /R L1 in the range 0.81 < R nose /R L1 < 0.87 (where R nose is the radius of the star in the direction of the pulsar, and R L1 the distance to the inner Lagrangian point); this corresponds to a volume-equivalent radius R 2 in the range 0.946 < R 2 /R RL < 0.974 (where R RL is the volume-equivalent radius of the Roche lobe). They also infer an inclination i in the range 63 • < i < 67 • .
Formally, the above ranges are at the 3σ level (Reynolds et al. 2007). This, however, does not take into account uncertainties in the models or the extent to which the underlying assumptions hold. For instance, the models may have temperature-dependent missing opacities, the metallicity or hydrogen abundance may not be solar, the pulsar wind that irradiates the companion may not be isotropic (though the lightcurve is remarkably symmetric), or the deep heating approximation may not be valid (although our blue spectra suggest it is not bad). Given these issues, we will treat the ranges as 1σ uncertainties below.
In addition to the above, an uncertainty that is more difficult to constrain relates to the extent to which heat is redistributed. Any redistribution would lead to a smoother temperature distribution, and thus it would require a larger inclination to obtain the same observed modulation amplitude. In principle, it should be possible to constrain the temperature distribution directly by fitting multi-band lightcurves and spectra simultaneously, similar to what has been done for NN Ser by Parsons et al. (2010), which is also strongly irradiated (and for which no evidence for heat transport was found). We hope to pursue this in the future.
Here, we will try to use the observations to set limits. For our purposes, the most important quantities are the companion radius, which determines the correction to the radial-velocity amplitude (see below), and the inclination, on which the final masses depend as 1/ sin 3 i.
For the radius, our spectra yield independent clues. First, the surface gravity is limited to log g 2 4. Since the minimum companion mass is M 2,min = 0.022 M ⊙ (from the pulsar mass function and our observed radial-velocity amplitude), one infers R 2,g = (g 2 /GM 2 ) 1/2 0.25 R ⊙ . Similarly, from the distance limit inferred from the reddening, we found R 2,d > R eff 0.25 R ⊙ .
These radii can be compared with the radius of the Roche lobe. For a companion much less massive than the pulsar, R RL ≃ 0.46a(M 2 /[M 1 + M 2 ]) 1/3 (Paczyński 1971). The scaling yields the well-known result that the mean density of the Roche lobe is determined just by the period; numerically, ρ RL = 0.185 g cm −3 (P/1 d) −2 (Eggleton 1983). For PSR B1957+20, ρ RL = 1.27 g cm −3 , and hence the size of the Roche lobe is R RL = (3M 2 /4πρ RL ) 1/3 = 0.29 R ⊙ (where the numerical value is for the minimum mass). Thus, we conclude that R 2,g /R RL 0.86. This is an overall lower limit, since R 2,g /R RL ∝ M 1/6 2 . Our limit is consistent with what was inferred from the lightcurve fit (as well as from theoretical considerations of the cause of orbital period variations; Applegate & Shaham 1994). It implies a lower limit R nose /R L1 0.7, which we will use below.
For the inclination, it is more difficult to set stringent limits. The fact that at low radio frequencies, a symmetric eclipse is seen that lasts ∼ 8% of the orbit (Fruchter et al. 1990;Ryba & Taylor 1991), while no eclipse is seen at high frequencies (Fruchter & Goss 1992), suggests that the conditions along the line of sight do not change too strongly during the eclipse. This seems easier to understand if the line of sight does not pass close to the companion, i.e., if the inclination is intermediate, i ≃ arccos(R E /a) ≃ 75 • (for fractional eclipse radius R E /a ≃ 0.08π), consistent with the inference from the lightcurve. Strictly, however, the eclipses only set a weak upper limit: the absence of eclipses at high frequencies implies i < arccos(R 2 /a) 85 • .
A lower limit to the inclination can be set from the large brightness contrast between superior and inferior conjunction (a factor 100 in the R band, Reynolds et al. 2007), which shows that at inferior conjunction at most a small part of the irradiated hemisphere is visible. Conservatively, we estimate that this requires i > 50 • . A similar lower limit seems reasonable from the fact that the radio eclipses last long; for lower inclinations, it is difficult to envisage a physical eclipse region that does not extend all the way to the pulsar. Overall, we conclude that all observations are consistent with the model inferences of Reynolds et al. (2007), with R nose /R L1 = 0.84 ± 0.03 and i = 65 ± 2 deg. Considering possible systematic uncertainties, a secure constraint on the radius seems to be 0.7 < R nose /R L1 < 1, while for the inclination the constraint is 50 < i < 85 deg.

MOTION OF THE CENTER OF MASS
The ratio between the observed radial-velocity amplitude of the center of light and the actual one of the center of mass can be written as K obs /K 2 = 1 − f eff R nose /a 2 , where the effective normalized emission radius f eff is constrained to be between 0 (uniform emission) and 1 (emission from the tip of the star facing the pulsar only). It should depend primarily on the surface brightness and line strength distributions in the observed band, with minor contributions arising from the exact shape of the star (determined by the mass ratio, filling factor, and the degree of co-rotation) and the orbital inclination. Indeed, for the somewhat similar situation of irradiation-induced Bowen emission lines on a Roche-lobe filling companion, Muñoz-Darias et al. (2005) found that the "K-correction" factor depends mainly on mass ratio (which determines the radius) and is nearly independent of inclination.
We investigated this hypothesis by adapting a lightcurve synthesis code to deal with spectra (the code was used by Stappers et al. (1999) to model pulsar irradiation, and is similar to that of Orosz & Hauschildt (2000)). The code produces not just fluxes, but also synthetic spectra, by summing Doppler-shifted spectra over all surface elements. While we do not yet have a suitable set of atmosphere models in hand to model our observations reliably, we have used the code to estimate the effect on the radial-velocity amplitude, by generating mock spectra for different sets of binary parameters. We generated spectra at the same orbital phases as our observations, determined radial velocities using the same method as done for the real observations ( § 3.2), and fitted circular orbits using the same weights. 4 At a fixed mass ratio, we found that the orbital inclination had negligible effect on f eff , 1% over the range 50 < i < 90 deg. The effects of the filling factor, for the range 0.7 < R nose /R L1 < 0.95, were a bit larger, though still small, at ∼ 7%. The influence of the precise values of strength of the irradiation is ∼ 6%, for the range of irradiation strengths that give front-side temperatures consistent with our spectral type (6000-6500 K); the backside temperature does not matter much, since it is constrained to be ∼ 2900 K (Reynolds et al. 2007), much too low to contribute B-band flux. As expected, with our choice of scaling with R nose , the effect of the mass ratio is very small, 1% for the range of values that are able to yield the observed K obs . Overall, we infer f eff ≃ 0.60 ± 0.04.
For comparison, looking from the side at a spherical star that is dark on one side and emits isotropically at a uniform temperature on the other, simple integration yields f eff = 4/3π = 0.42. Assuming instead a temperature distribution T = T 0 cos 1/4 θ, as expected for irradiation by a parallel beam, we find f eff = 0.62 for T 0 = 6400 K (assuming black-body emission and a linear limb darkening law with u = 0.6). This is very similar to what we find from our model, showing that the increase in brightness dominates over effects such as an increase in line strength with decreasing temperature. Before using the above to estimate K 2 , it is useful first to consider the limits. For uniform emission, K 2 = K obs = 324 ± 3 km s −1 , which yields a lower limit to the mass ratio for the system q min = M 1 /M 2 = K obs /K 1 = 63.6 ± 0.6 (where the pulsar's radial-velocity amplitude K 1 = 2πa 1 sin i/P orb = 5.09272 ± 0.00004 km s −1 is known from radio timing; Arzoumanian et al. 1994). At the other extreme, for light emitted by the nose of a Roche-lobe filling star, i.e., f eff R nose = R L1 , one finds K obs = 0.843K 2 , i.e., K 2 = 384 ± 4 km s −1 and q max = 75.4 ± 0.7. We will take these values as limits in § 6 (accounting for the fact that neither can be close to realistic by ignoring the uncertainty on K obs ).
For our more general case, combining our estimate of f eff = 0.60 ± 0.04 above with R nose /R L1 = 0.84 ± 0.03, and solving for K 2 (taking into account that R L1 /a depends on q), one finds K 2 = 353 ± 4 km s −1 . For the conservative range in companion radius, 0.7 < R nose /R L1 < 1, we find 348 ± 4 < K 2 < 358 ± 4 km s −1 , i.e., it corresponds to a 5 km s −1 uncertainty in K 2 . Adding this in quadrature to the 2σ uncertainty of 8 km s −1 arising from the uncertainties in K obs and f eff , we infer a conservative range in radial-velocity amplitude 343 < K 2 < 363 km s −1 . From the above, we conclude that the mass ratio is q = 69.2 ± 0.8, and that a conservative range is 67.3 < q < 71.3.

INFERRED MASSES AND CONCLUSIONS
In Fig. 3, we show our constraints on the masses. One sees that most likely, PSR B1957+20 is massive, with M PSR = 2.40 M ⊙ . Taking the inferences from Reynolds et al. (2007) on the lightcurve and the corresponding inclination, and our correction to the radial-velocity amplitude at face value, the formal uncertainty is small, ∼ 0.12 M ⊙ .
As discussed in § 4, however, the lightcurve modelling relies on a number of assumptions, especially that there is no heat transport over the face of the star. From our conservative constraints on both the inclination and the mass ratio, we find a lower limit to the mass of 1.66 M ⊙ .
Thus, from our work we conclude that PSR B1957+20 certainly is more massive than the canonical 1.35 M ⊙ , and likely substantially more massive. Indeed, it may well be more massive even than PSR J1614-2230 (1.97 ± 0.04; Demorest et al. 2010), and thus allow even more stringent constraints on the equation of state. The large mass also suggests a large amount of mass was transferred in the preceding phase as an X-ray binary, although this conclusion depends on the initial mass. However, even if that were as high as the mass found for Vela X-1 (∼ 1.9 M ⊙ , Barziv et al. 2001), our measurements suggest the pulsar has accreted about half a solar mass. 5 To confirm the high inferred mass will require more secure constraints on the orbital inclination and, to a lesser extent, 5 Since Vela X-1 has a massive companion, its mass must still be close to the one it formed with. -Mass-mass diagram for PSR B1957+20 and its companion. The cross indicates our best-fit solution with 1σ uncertainties, and the surrounding blue parallelogram the convervative region including our best estimate of possible systematic uncertainties (see text). The physically allowed region (light blue) is limited by the constraints sin i ≤ 1 (thick, blue line), and K L1 ≤ K obs ≤ K 2 (thick, red lines). Here, the second constraint arises because the observed radial-velocity amplitude K 2 is measured using light emitted from the side facing the pulsar, and the center of light cannot be further away from the pulsar than the center of mass, or closer to the pulsar than the first Lagrangian point L1 (with velocity amplitudes K 2 and K L1 , respectively). For reference, we show contours of constant inclination i (dotted), calculated using the pulsar mass function. the correction factor for the radial-velocity amplitude. For this purpose, most important would be to model the lightcurve in more bands, and to check explicitly what inclinations are possible for less-constrained, perhaps even arbitrary temperature distributions. It may be especially valuable to model the spectra at the same time, thus avoiding the indirect calculation of a correction factor. We are currently working on rewriting our code for this purpose.
An improved, nearly model-independent constraint on the mass ratio could be obtained from a near-infrared radialvelocity curve. Since in the near-infrared the contribution of the cold side to the light budget is more important, the radialvelocity amplitude of the center of light at infrared wavelengths should be much closer to that of the companion's center of mass, and hence the uncertainty in the corrections much less important.
Furthermore, one could determine the projected rotational velocity v sin i from high-resolution spectra. This would allow one to check the predictions from the models, in particular for the mass ratio and the filling factor, on which v sin i depends most strongly.
Finally, radio observations could provide a complementary improved constraint on the inclination, from mapping the eclipse at a larger range of frequencies than was done by Fruchter & Goss (1992), and making use of the large increases in sensitivity, especially at high frequency, that have been made over the last decades.
The data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. We also used data from the UVES Paranal Observatory Project UVESPOP (ESO DDT Program ID 266.D-5655). We made extensive use of SIMBAD and ADS.