General relativity reveals pulsar beams Pulsars are rotating neutron stars that emit beams of radio waves along their magnetic poles, seen as regular pulses if the beam points toward Earth. Desvignes et al. monitored a pulsar for more than a decade, observing how its radio pulses vary. General relativity causes precession of the rotation axis, because of the influence of a binary companion. In 2005, two pulses per rotation were visible, one from each magnetic pole, but by 2018 one had precessed out of our line of sight and disappeared. Mapping the radio emission across the magnetic pole determines the beaming angle, the angular region in which a radio observer can detect a pulsar. Science, this issue p. 1013

Abstract Binary pulsars are affected by general relativity (GR), causing the spin axis of each pulsar to precess. We present polarimetric radio observations of the pulsar PSR J1906+0746 that demonstrate the validity of the geometrical model of pulsar polarization. We reconstruct the (sky-projected) polarization emission map over the pulsar’s magnetic pole and predict the disappearance of the detectable emission by 2028. Two tests of GR are performed using this system, including the spin precession for strongly self-gravitating bodies. We constrain the relativistic treatment of the pulsar polarization model and measure the pulsar beaming fraction, with implications for the population of neutron stars and the expected rate of neutron star mergers.

Pulsars are fast-spinning neutron stars measuring ~1.2 to 2.2 solar masses ( M ☉ ) with strong magnetic fields that emit a beam of radio waves along their magnetic axes above each of their opposite magnetic poles. According to Einstein’s theory of general relativity (GR), space-time is curved by massive bodies. Predicted effects of this theory include relativistic spin precession in binary pulsars (1). This precession arises from any misalignment, by an angle δ, of the spin vector of each pulsar with respect to the total angular momentum vector of the binary, most likely caused by an asymmetric supernova (SN) explosion imparting a kick onto the neutron star(s) (2, 3). This precession causes the viewing geometry to vary, which can be tested observationally.

As a pulsar rotates, its radio beams sweep the sky. If one of the beams crosses our line of sight (LOS), its emission is perceived as being pulsed. When averaged over several hundreds of pulsar rotations, the pulses typically form a stable pulse profile. Evidence for a variable pulse profile attributed to changes in the viewing geometry caused by spin precession have been observed and modeled for the binary pulsar PSR B1913+16 (4, 5). Polarization information can provide an additional, independent tool to study relativistic spin precession (6, 7). We expect that the position angle (PA) sweep of a pulsar’s linearly polarized emission resulting from geometrical effects can be described by the rotating vector model (RVM) (8). This simple model, which assumes a magnetic dipole centered on the pulsar, relates the PA to the projection of the magnetic field line direction as the pulsar beam rotates and crosses our LOS. The resulting gradient of the PA sweep as a function of the pulsar rotational phase (9) depends only on the magnetic inclination angle α and the impact parameter β (i.e., the angle of closest approach between the observer direction and the magnetic axis). For LOSs crossing opposite sides of the same magnetic pole, the RVM predicts opposite slopes of the PA swing, which becomes steeper for smaller β (10). RVM has been extended to include rotational and relativistic effects between the pulsar and observer frame (11, 12), in principle allowing emission heights for the observed radio emission to be estimated. Although the RVM matches observations of young pulsars that present a smooth PA swing (13), deviations are also observed, so for large emission heights, the pure dipole approximation may not be valid for a rotating plasma-loaded magnetosphere (14). Evidence for the central assumption of all these models, i.e., the geometrical meaning of the PA sweep, has so far been missing.

PSR J1906+0746 (right ascension 19h06m48.86s, declination +07°46′25.9′′, J2000 equinox) is a young pulsar with spin period P s ~ 144 ms in a 4-hour orbit around another neutron star. When it was discovered in 2004 (15), PSR J1906+0746 showed two polarized emission components separated by nearly half a period (or ~180° of pulse longitude). The main pulse (MP) and interpulse (IP) indicated a nearly orthogonal geometry where emission from both magnetic poles is visible from Earth. Comparison with archival data from the Parkes Multibeam Pulsar Survey (PMPS) (16) revealed that only the stronger MP had been visible in 1998 (15), suggesting effects of relativistic spin precession. Timing analysis of this pulsar (17) has measured three relativistic corrections to the orbit, the post-Keplerian (PK) parameters. Two of the PK parameters, the periastron advance ω ˙ and the time dilation γ, allow the determination of the pulsar and companion masses, m p = 1.29 ± 0.01 M ☉ and m c = 1.32 ± 0.01 M ☉ , respectively, assuming that GR is correct. The third PK parameter, the orbital decay due to gravitational wave emission P ˙ b , is consistent with these measurements, providing a test of GR with 5% precision (17). GR also provides an estimate of the orbital inclination angle i = 43.7° ± 0.4° or 136.3° ± 0.4° (due to a 180° degeneracy) and the spin precession rate Ω p = 2.234° ± 0.014° per year. Independent measurements of these two terms are potentially observable in precessing systems, allowing these GR-predicted values to be tested. Observations of the double pulsar (18) and PSR B1534+12 (7, 19) match the GR predictions of spin precession with precisions of 13 and 20%, respectively.

We monitored PSR J1906+0746 from 2012 to 2018 with the 305-m William E. Gordon Arecibo radio telescope using the Puerto Rico Ultimate Pulsar Processing Instrument (PUPPI) tuned at a central frequency of 1.38 GHz. We supplement those observations with archival data from the Nançay and Arecibo telescopes recorded between 2005 and 2009 (17, 20). In total, our dataset comprises 47 epochs spanning from July 2005 to June 2018.

A previously observed trend of changing separation between MP and IP profiles with time (17) has continued. Our data show that the slope of the PA under the MP gradually flattens with time, whereas the MP becomes simultaneously weaker and subsequently undetectable toward the end of 2016. This indicates that our LOS has moved out of the MP emission beam.

The progressive steepening of the PA curve under the IP eventually leads to a flip of the PA around May 2014 (although it appears flat because of the 180° ambiguity in the PA), followed by a flattening. This behavior is in agreement with the RVM for a pole crossing our LOS. The observed change in the PA curve provides evidence that the PA swing has a geometrical origin linked to the magnetic field.

The presence of MP and IP emission, covering a wide range of pulse longitudes, allows a precise determination of the viewing geometry as a function of time. We performed a simultaneous modeling of the polarimetric profiles using the RVM including the effects of relativistic spin precession (21), hereafter referred to as the precessional RVM. A total of 53 parameters are included in this model (see table S1); the main parameters are the angle between the rotation axis and the magnetic axis of the MP, α MP ; the misalignment angle, δ; the inclination angle, i; and the precession rate, Ω p . The phase of the inflection point of the RVM under the MP, ϕ 0 M P , k , for each epoch k was also included in the analysis as a free parameter. The geometry of the system is illustrated in Fig. 1. We used the Bayesian sampling tool PolyChord (22) to explore the parameter space of the model (20). The 1σ uncertainty intervals were derived from the one-dimensional marginalized posterior distributions (shown in fig. S3), giving the 68% confidence level on each parameter.

Fig. 1 Geometry of PSR J1906+0746. The pulsar rotates with a spin period P s = 144 ms around its spin vector S shown with the vertical red arrow. S precesses with a period of P p = 360°/Ω p ~160 yr around the total angular momentum vector (misaligned by the angle δ from S) that can be approximated by the orbital momentum vector X, perpendicular to the Y–Z orbital plane. The orbital inclination angle is i. As the pulsar spins, its magnetic pole corresponding to the MP emission, B MP , and inclined with an angle α MP sweeps the sky along the dashed blue trajectory. The MP beam, with the extent pictured by the dotted blue circle, crosses our LOS represented by the K vector if | β M P | < 22 ° . This allows us to observe a cut, shown with the red curve, through the MP beam. After half a rotation of the pulsar, our LOS can also potentially cut through the IP beam if | β I P | < 22 ° . For clarity, only the MP beam is shown.

We measure α MP = 99°.41 ± 0°.17 (and therefore α IP = 180° − α MP = 80°.62 ± 0°.17), which combined with a separation between the MP and IP close to 180° confirms that PSR J1906+0746 is an orthogonal rotator and that the MP and IP emissions originate from opposite magnetic poles. The large δ =104° ± 9° is consistent with the pulsar having formed in an asymmetric SN explosion, producing a tilt of the spin axis of the younger pulsar in the binary pair (17) relative to the pre-SN orbit (3). As the spin vector contribution to the total angular momentum vector is negligible, we also interpret δ as the angle between the orbital angular momentum vector and the pulsar spin vector. With δ ~ 104°, this suggests that the pulsar spin axis lies close to the orbital plane of this binary system.

This analysis independently determines Ω p = 2.17° ± 0.11° per year and i = 45° ± 3° in addition to the three previously measured PK parameters (17). This allows us to perform two additional self-consistent tests of GR in this system (see Fig. 2). Both values agree with the GR predictions within 1σ. The constraint on Ω p tests GR at a 5% uncertainty level, which is tighter than the precession rate measurement in the double pulsar system (18). The parametrization of our model allows us to determine i without ambiguity, in contrast to information obtained solely from pulsar timing where only sin i can be derived.

Fig. 2 Mass-mass diagram. The black lines delimit the 1 − σ contours from the measurements of the orbital period decay P ˙ b , the periastron advance ω ˙ , and the time dilation γ (17). The red dotted and dashed lines delimit our two additional constraints, the measurements of spin-precession rate Ω p and inclination angle i, respectively. The gray parameter space is excluded by the mass function because sin i ≤ 1 . All parameters are consistent with m p = 1.29 ± 0.01 M ☉ and m c =1.32 ± 0.01 M ☉ .

From our fitted model, we can describe the evolution of the impact parameters for both the MP and IP, β MP (t) and β IP (t) (see fig. S4) and the latitudinal extent of the pulsar beam (i.e., the largest value of | β | for which emission is observed). Our model predicts that in 1998, β MP ~ 5° and β IP > 22°, coinciding with the strong detection of the MP and the nondetection of the IP in the PMPS data (15). The appearance of the IP between 1998 and 2004 corresponds to β IP (2001) ~ +20°. By 2018, our LOS to the IP had traversed to β IP (2018) ~ −6°, so we infer a latitudinal extent for the IP of ~20°. For the MP, we are only able to map the southern part of the emission beam, − 5 ° < β M P ≲ − 22 ° , after which the emission is no longer detected. The appearance of the IP and the disappearance of the MP at the same | β | suggests that the MP and IP share the same latitudinal extent of ~22°. We therefore predict that the IP will disappear from our LOS around 2028 and will reappear between 2070 and 2090 (fig. S4). The MP should reappear between 2085 and 2105.

With our determination of the geometry of the system, we can use the observed pulse profiles recorded at different epochs, including the reprocessed PMPS data (20), to reconstruct (20) the sky-projected beam maps of the radio emission from both magnetic poles (Fig. 3) and its polarization properties (Fig. 4). Previous attempts at mapping the pulsar emission have been published for other pulsars [e.g., (23)]. Such plots must not be viewed as strict maps of active field lines anchored on the polar cap but rather as projections of the emission on the sky.

Fig. 3 Beam maps of the radio emission. (A to E) show the PA measurements (in black) and the fit by the precessional RVM (the red curve) corresponding to different LOSs through the IP. (F) and (G) show the two-dimensional contours of the reconstructed Stokes I emission maps (on a logarithmic color scale), as projected on the sky, for the MP and the IP, respectively. The reference frame of these maps is fixed on the vertical spin vector S at longitude zero and 180° pointing toward negative latitude. The red crosses at (0,0) and (180,0) indicate the pulsar’s magnetic pole on the MP and IP maps, respectively, and the dashed circles show increments of 2° in the beam map. The dotted lines represent the LOSs at a given year, indicated on the left side, and, for the IP, refer to the PA measurements (A to E). The hatched areas correspond to the parts of the maps that were not observed.

Fig. 4 Polarization beam maps. Same as Fig. 3, F and G, but showing fractional linear polarization for the MP (A) and IP (B) and fractional circular polarization for the MP (C) and IP (D).

The map shows emission from both sides of the IP magnetic pole, ruling out models predicting that radio emission is restricted to one side of the pole [e.g., (24)]. The IP emission pattern is not symmetric in the latitudinal direction or with the MP. The IP weakens when our LOS crosses the magnetic pole. This finding matches theoretical predictions of the current density in the polar cap for an orthogonal rotator, albeit at low radio emission height (25, 26). Radio emission is produced at a height above the pulsar where any higher-order magnetic field components should have diminished. In a standard scenario where emission is produced in the entire open field-line region, we would expect the radio beam to be circular. Yet the observed beam is elongated and smaller than expected in the longitudinal direction (Fig. 3).

The passage of the LOS over the IP magnetic pole coincides with a drop in fractional linear polarization. At the same time, the circular polarization (Stokes V) appears to change sign when crossing the magnetic pole (Fig. 4 and fig. S5) (20). The fractional linear polarization of the MP decreases as our LOS moves away from the center of the MP beam, whereas some unpolarized pulse components are shown to appear and disappear on time scales of months (20).

If we assume that the emission region is symmetric around the magnetic meridian and that the radiation is beamed in the forward direction of the relativistic charge flow along the magnetic field lines, then the rotation of the pulsar magnetosphere in the observer frame causes the pulse profile to precede the magnetic meridian and the PA to lag behind it.

Models predict that the total phase shift Δ ϕ S between the midpoint of the pulse profile and the inflection point of the PA swing is given by Δ ϕ S ≈ 4 h e m / R L C , where h em is the emission height and the light cylinder radius R LC = cP s /2π, with c being the speed of light (11, 12). This relationship is expected to be valid for small emission heights, h e m ≲ 0.1 R L C . But these models predict different contributions to the aforementioned shifts in the pulse profile and the PA curve. There is also a shift in absolute value of the PA, but this can be neglected here as α ~ 90° (27).

Measuring the phase shifts Δ ϕ S M P , k for the MP at a given epoch k [and its impact parameter β(k)] gives an estimate of the MP emission heights h e m M P , k , assuming all shifts are attributable to rotational effects. The emission heights for the IP h e m I P , k can be derived in the same way, also using the RVM inflection points for the IP given by ϕ 0 I P , k = ϕ 0 M P , k + 180 ° . These results show that the emission heights range from close to the surface of the pulsar, for low impact parameter β (when our LOS is atop the magnetic pole), to ~320 km for large β, matching the theoretical prediction (28) (Fig. 5), although this theoretical prediction only depends on α and β. Deviations are observed for LOSs close to the beam edge of the MP, and large β for the IP, where we may observe a different active patch of the IP beam as suggested by the polarization properties (e.g., the sign of Stokes V under the IP has flipped between 2009 and 2012) (see Fig. 4 and fig. S5).

Fig. 5 Radio emission heights. Phase shift Δ ϕ S of the PA from Stokes I (left axis) with the derived emission height (right axis) for the MP (blue data points) and IP (red data points) as a function of β. The blue and red dashed curves represent the theoretical emission height as a function of β.

The latitudinal extent of a pulsar beam determines the pulsar beaming fraction (i.e., the portion of sky illuminated by the pulsar beam). This parameter affects the predicted number of Galactic double neutron stars (DNSs) and, hence, the expected gravitational wave detection rate for neutron star mergers (29). Usually, the latitudinal extent of the beam cannot be independently determined and so the expected beam radius ρ for a conal emission beam is used. This beam shape is typically expected, although its internal structure and/or filling factor is still a matter of debate (10). On the basis of beam cuts given by pulsar profiles, various studies [e.g., (30)] have suggested ρ ∼ 6.5 ° P s − 0.5 , where the dependence on P s is consistent with the expectation for dipolar field lines.

Previous studies that estimated the DNS population number and merger rate (29) from a sample of pulsars dominated by only two pulsars (PSR J0737−3039A and PSR J1906+0746) have assumed for PSR J1906+0746 a uniform radio luminosity across a conal beam of ρ ~ 17°. This value is smaller than our observed latitudinal extent of 22° that, combined with α MP = 99°.4, gives a large beaming fraction of 0.52 but with substantial luminosity variations. The simple conal beam model with uniform luminosity cannot be used to reliably constrain the DNS merger rate for PSR J1906+0746.

Supplementary Materials science.sciencemag.org/content/365/6457/1013/suppl/DC1 Materials and Methods Supplementary Text Figs. S1 to S6 Tables S1 and S2 References (34–46)

http://www.sciencemag.org/about/science-licenses-journal-article-reuse This is an article distributed under the terms of the Science Journals Default License.