# Revisiting the Galactic Double Neutron Star merger and LIGO detection rates

Kathrin Grunthal, Michael Kramer, Gregory Desvignes

► To cite this version:

Kathrin Grunthal, Michael Kramer, Gregory Desvignes. Revisiting the Galactic Double Neutron Star merger and LIGO detection rates. Monthly Notices of the Royal Astronomical Society, 2021, 507 (4), pp.5658-5670. 10.1093/mnras/stab2198. hal-03319388

HAL Id: hal-03319388

https://hal.science/hal-03319388

Submitted on 21 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Revisiting the Galactic Double Neutron Star merger and LIGO detection rates

K. Grunthal,<fnref n="1" /> M. Kramer and G. Desvignes<fnref n="1" />,<fnref n="3" />

<figure><img src="image_3.png" /></figure>

Accepted 2021 July 27. Received 2021 July 26; in original form 2021 June 5

## ABSTRACT

We revisit the merger rate for Galactic double neutron star (DNS) systems in light of recent observational insight into the longitudinal and latitudinal beam shape of the relativistic DNS PSR J1906 + 0746. Due to its young age and its relativistic orbit, the pulsar contributes significantly to the estimate of the joint Galactic merger rate. We follow previous analyses by modelling the underlying pulsar population of nine merging DNS systems and study the impact and resulting uncertainties when replacing simplifying assumptions made in the past with actual knowledge of the beam shape, its extent, and the viewing geometry. We find that the individual contribution of PSR J1906 + 0746 increases to $R = 6_{-5}^{+28}$ Myr$^{-1}$ although the values are still consistent with previous estimates, given the uncertainties. We also compute contributions to the merger rates from the other DNS systems by applying a generic beam shape derived from that of PSR J1906 + 0746, evaluating the impact of previous assumptions. We derive a joint Galactic DNS merger rate of $R_{\text{MW}}^{\text{gen}} = 32_{-9}^{+19}$ Myr$^{-1}$, leading to a LIGO detection rate of $R_{\text{LIGO}}^{\text{gen}} = 3.5_{-1.0}^{+2.1}$ yr$^{-1}$ (90 per cent conf. limit), considering the upcoming O3 sensitivity of LIGO. As these values are in good agreement with previous estimates, we conclude that the method of estimating the DNS merger and LIGO detection rates via the study of the radio pulsar DNS population is less prone to systematic uncertainties than previously thought.

**Key words:** stars: neutron – pulsars: general – neutron star mergers.

## 1 INTRODUCTION

Double neutron star (DNS) systems are of major interest for a number of reasons. First, their formation invokes a number of questions about the progenitor systems, the birth of neutron stars themselves via different possible mechanisms, and their size and origin of resulting properties like mass, spin, or orbital configuration (Tauris et al. 2017). Secondly, if (at least) one of the neutron stars is visible as a radio pulsar, DNS systems can act as tools for precision tests of strong-field gravity (e.g. Kramer et al. 2006). Indeed, the first DNS system known (PSR B1913 + 16), discovered in 1974 and now widely known as the Hulse–Taylor Pulsar (Hulse & Taylor 1975), provided the first evidence for gravitational waves (GW) (Taylor & Weisberg 1982). Such DNSs visible as binary pulsars are rare, and only about 20 DNS systems have been detected in radio surveys to date (e.g. Tauris et al. 2017; Haniewicz et al. 2021). This is not only due to their relative intrinsic rareness caused by the conditions during the formation process but also due to the difficulties in detecting the pulsar in fast binary orbits (e.g. Tauris et al. 2017).

The third reason, why these systems are of interest, is given by the fact that for compact DNS, the orbital decay of these systems inevitably leads to the merger of the two neutron stars, an event

that creates a copious amount of GWs during the merger, which can be picked up by terrestrial GW detectors. Indeed, following the discovery of the Hulse–Taylor pulsar, expected to merge with 300 Myr, it was anticipated that the first GW event seen by ground-based observatories, such as LIGO (Laser Interferometer Gravitational-Wave Observatory) or Virgo, would be a merger of two NSs (Barish 1999). Consequently, estimates of the DNS merger rate have been made frequently, based on population synthesis calculation or computations informed by the population of compact DNS observed as radio pulsars.

Interestingly, the first detected GW merger signal originated from a binary black hole, while GWs coming from a DNS merger were not detected until 2017 August. This latter event named GW170817 (Abbott et al. 2017a) revealed an electromagnetic counterpart, which was followed up across the electromagnetic spectrum, truly starting the era of multimessenger astronomy (Abbott et al. 2017b). Since then, however, only one additional DNS-merger event, GW190425, has been detected (Abbott et al. 2020). Hence, improving our knowledge of the DNS merger rate is still of great scientific interest. On one hand, it will help predicting the actual event rate, preparing also electromagnetic follow-up observations. On the other hand, once the actual event rate has been measured more precisely in upcoming detector runs, we can in turn use the comparison with theoretical predictions based on the current understanding of DNS populations in order to re-calibrate our knowledge about the formation and evolution of Galactic DNS systems.

<footnote n="1">¹Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany</footnote>

<footnote n="2">²Jodrell Bank Centre for Astrophysics, University of Manchester, Manchester M13 9PL, UK</footnote>

<footnote n="3">³LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Université de Paris, 5 Place Jules Janssen, F-92195, Meudon, France</footnote>

<footnote>* E-mail: mkramer@mpifr-bonn.mpg.de</footnote>
The detection rate for ground-based GW detectors can be estimated based on the assumption that the DNS merger rate in the Milky Way can be extrapolated to external galaxies in the GW detector range. Essentially, one can obtain the expected Galactic DNS merger rate via two routes. One is the computation by first principles, using ab initio population synthesis computations (e.g. Belczynski et al. 2008). Another way is to infer the merger rate from the known Galactic population of DNS systems, which are observed as binary radio pulsars. However, discovering the relevant DNS systems, i.e. those that merge due to the emission of GWs within a Hubble time, is typically subject to a number of selection effects (Lorimer & Kramer 2005). Therefore, it is necessary to use methods that take these into account.

Kim, Kalogera & Lorimer (2003) developed a method to derive a merger rate probability distribution using the known DNS population and Bayesian statistics. The combination of all individual DNS contribution then yields the merger rate of the Milky Way $R_{MW}$. By extrapolating that rate to the observable volume of the LIGO/Virgo system, the number of DNS mergers that the observatories will detect can be predicted. The method was also recently applied by Pol, McLaughlin & Lorimer (2019) and Pol, McLaughlin & Lorimer (2020). In the latter update, the authors presented a Milky Way merger rate of $R_{MW}^{Pol} = 37(+24, -11)$ Myr$^{-1}$ and an inferred LIGO detection rate of $R_{LIGO}^{Pol} = 4.2(+2.6, -1.3)$ yr$^{-1}$ (90 per cent confidence limit, c.l.). The extrapolation from $R_{MW}^{Pol}$ to $R_{LIGO}^{Pol}$ was done for the LIGO O3 range distance of 130 Mpc.

There are three factors that determine the rate constraints for pulsar binaries, namely the number of suitable binary pulsars, their effective lifetime, and their beaming fraction, i.e. the fraction of the sky potentially illuminated by a given pulsar. The number of pulsars can be estimated via survey simulations as explained further below, which takes the sensitivity of past and current observations (including selection effects due to binary motion), as well as the luminosity distribution of the sources into account.

Following O'Shaughnessy & Kim (2010) and Pol et al. (2019), the lifetime of a pulsar binary can be estimated from the sum of the current age of the system and the remaining detectable lifetime. The age can be estimated from the current spin-properties. For a non-recycled pulsar, one uses the characteristic age (e.g. Lorimer & Kramer 2005). For a recycled pulsar, one compares the rotation period and spin-down rate with their combination given by the so-called ‘spin-up line’. The remaining lifetime is the time that the pulsar needs to spin-down either until it crosses the pulsar death line (see e.g. Lorimer & Kramer 2005) or until it merges with its companion, whatever is smaller.

The beaming fraction takes into account that radio surveys only detect those systems, where the pulsar beam is directed towards Earth, whereas GW detectors are unaffected by such selection effects. Following Kim et al. (2003) or O'Shaughnessy & Kim (2010), one can introduce the pulsar beaming correction factor, $f_b$, as the inverse of the fraction of the sky that is illuminated by a given pulsar's beam (accounting for both magnetic poles). The beaming fraction depends on the geometry of the pulsar, i.e. the magnetic inclination angle between the spin and magnetic axes, $\alpha$, and the angular size of the emission beam, $\rho$. These two quantities also affect the pulse width, $W$, of the observed pulse profile, which obviously is only a one-dimensional cut through an emission beam with both longitudinal and latitudinal dimensions. In order to estimate the overall illuminated fraction of the sky, one usually makes the assumption that the inferred longitudinal dimension also applies to the latitudinal direction, i.e. a circular beam shape.

But even the inference of the longitudinal dimension, based on the pulse width, is often problematic as it requires, apart from $\alpha$ and the assumption that emission completely fills the pulsar beam area, also knowledge of the impact angle, $\beta$, i.e. the smallest angular separation of our line of sight with the magnetic axis. In this case,

$$ \cos \rho = \cos \alpha \cos(\alpha + \beta) + \sin \alpha \sin(\alpha + \beta) \cos \left( \frac{W}{2} \right) \quad (1) $$

Gil, Gronkowski & Rudnicki (1984). In principle, it is possible to determine $\alpha$ and $\beta$ from polarization information by applying the so-called ‘Rotating Vector Model’ (RVM) (Radhakrishnan & Cooke 1969). However, uncertainties due to covariances in the parameters can be large (e.g. Lorimer & Kramer 2005), while it is not clear how well the RVM is applicable to recycled pulsars (see e.g. Kramer et al. 2021).

In cases where the geometry cannot be determined, one can resort to apply knowledge inferred for population properties to estimate $\rho$. Based on the analysis of non-recycled pulsars, a variety of authors (e.g. Gil, Kijak & Seiradakis 1993; Rankin 1993; Kramer et al. 1994; Gould & Lyne 1998) have determined that the value of $\rho$ scales with the pulse period as $\rho = A \times P^{-0.5}$, whereas A ranges from 4.9 to 6.5 deg $s^{0.5}$, depending on the study and the intensity level that the width, $W$, was measured for (see the discussion by Venkatraman Krishnan et al. 2019). This scaling relationship correctly reflects the period dependency as expected for a dipolar field line structure. Nevertheless, empirically the relationship breaks down for recycled pulsars (at about $P \sim 10-30$ ms), which typically appear to show a smaller beam size than expected from period scaling (Kramer et al. 1998).

O'Shaughnessy & Kim (2010) provided a careful analysis, combining knowledge from the known population of recycled and non-recycled pulsars and all possible information for the visible pulsars in DNS systems known at the time, in order to derive an effective beaming correction factor, $f_b$, in a recent variation of this work, Pol et al. (2019) also adopted overall the approach by O'Shaughnessy & Kim (2010) (and earlier by Kim et al. 2003) to assume a fixed $f_b$, $\text{eff}$ value (set to 4.6), while adopting the actually inferred beaming correction factor $f_b$ for those three DNS with a determined geometry, i.e. PSRs J0737 - 3039A, B1534 + 12, and B1913 + 16. Given the crucial nature of the beaming correction factor, it is desirable to estimate the impact of the simplifying assumptions made.

Recently, long-term observations of the relativistic binary pulsar PSR J1906 + 0746 allowed to map the emission beam of the visible, non-recycled pulsar of this DNS also in latitudinal direction. This is possible due to the effects of relativistic spin-precession, which slowly changes the impact angle and hence the location of the one-dimensional cut through the beam (Desvignes et al. 2019). Spanning observations of more than 10 yr, Desvignes et al. (2019) were able to determine the beam shape, even obtaining information from both magnetic poles (via the ‘main pulse’ and ‘interpulse’ emission separated by half a period) that can be combined to estimate for the first time reliably the previously unknown $f_b$ and the exact viewing geometry of the pulsar (see Section 2).

These measurements show additionally that the beam intensity profile can neither be modelled by a box-shaped function with the opening angle as the radius nor that the main pulse and the interpulse have the same intensity profile and intensity maxima (for a given magnetic latitude), as it is usually assumed. As PSR J1906 + 0746 significantly contributes to $R_{MW}^{Pol}$ (Lorimer et al. 2006; O'Shaughnessy & Kim 2010), it is therefore an ideal test case to study the impact of making particular assumptions on the robustness of the derived DNS merger rate estimates. One example that we can study
in detail is indeed the importance of the applied certain beaming correction factors, such as assuming a uniform value rather than using one that is better reflecting the intrinsic properties of a given pulsar. Such a study does not only allow us to evaluate the robustness of the simulations and to better understand the uncertainties but also provides us with an updated merger rate.

The plan for the paper is as follows. First, we briefly summarize the information about PSR J1906 + 0746, which plays a central role in our study by providing direct and accurate measurements of the extent of an emission beam and hence the beaming fraction. We then summarize the statistical approach adopted in this work, following Kim et al. (2003) and Pol et al. (2019). This is followed by a description how we use the new information on PSR J1906 + 0746 to improve on the previous studies, which we extend then also to other DNS systems, before we discuss the results and draw conclusions.

## 2 PSR J1906 + 0746

Pulsar J1906 + 0746 was discovered in 2006 (Lorimer et al. 2006). The 144-ms pulsar with a small characteristic age of only 113 kyr has an unseen companion in an eccentric ($e = 0.085$) 4-hr orbit. Even though the spin properties of the pulsar suggested that it is young and unrecycled, the exact nature of the binary companion was initially unclear. The matter was settled, when van Leeuwen et al. (2015) presented an improved timing solution including the measurement of three post-Keplerian orbital parameters (e.g. Lorimer & Kramer 2005), resulting in mass measurements for the pulsar of $1.291 \pm 0.01M_\odot$, and for the companion of $1.32 \pm 0.01M_\odot$, respectively. Combined with the orbital parameters, the masses identify the companion as another, recycled but unseen NS, while making PSR J1906 + 0746 the second-born NS in the system.

After the discovery, the pulsar was identified also in archival data, revealing that by the time of the discovery, emission from the second magnetic pole had become visible as an 'interpulse', which was missing earlier (Lorimer et al. 2006). The origin of this change in the observed emission is caused by a change in viewing geometry due to relativistic spin precession that was studied in great detail by Desvignes et al. (2019). Desvignes et al. were not only able to provide the best test of relativistic spin precession as predicted by General Relativity thus far, but they were also able to determine the viewing geometry accurately. They found the inclination angle of the magnetic axis related to the 'main pulse' to be $\alpha = 99.4 \pm 0.2^\circ$ and a large tilt angle of the spin axis of the pulsar to the total angular momentum vector of $\delta = 104 \pm 9^\circ$. The observed pulsed emission evolves with time, with the main pulse now having almost disappeared from sight, while the interpulse has remained to be observable for the moment. The slowly changing line of sight through the emission beams provides a 'tomographic' view of the emission regions as it is projected on the plane of the sky. In contrast to existing phenomenological models (see Lorimer & Kramer 2005 for an overview), the emission is not symmetric about the magnetic axis, neither in structure nor in extent. While the longitudinal extent of the emission beams varies between $10^\circ$ and $20^\circ$, the latitudinal extent is observed to be about $20^\circ$. Most importantly, however, the luminosity across the beam is highly non-uniform, decreasing noticeably before the beam 'edge' is reached. This is in sharp contrast to the simple conal beam model with uniform luminosity that has been applied so far in studies to infer the DNS merger rate.

Addressing the deviation from the simple model used so far is, actually, of particular importance for PSR J1906 + 0746. With its young age, scale factors inferred to estimate the number of similar pulsars existing in the Milky Way are unusually large, making the

contribution of PSR J1906 + 0746 to the overall estimated merger rate comparable to that of the Double Pulsar,<fnref n="1" /> as already pointed out by O'Shaughnessy & Kim (2010). It is therefore especially important to revisit the previously inferred merger rate in light of the new recent information provided by Desvignes et al. (2019).

## 3 CALCULATING DNS MERGER RATES FROM SIMULATIONS WITH THE PSRPOPPY PACKAGE

Since all known pulsars represent only a small fraction of the total pulsar population in the Milky Way Galaxy, it is necessary to use models and simulations to account for the various selection effects described earlier in order to estimate the underlying population. We follow previous works, in particular that of Pol et al. (2019), and use the simulation code PsrPopPy2² as well as the separate, available analysis code used by Pol et al. (2019).<fnref n="3" /> While we modify the code as described later to enable our updated analysis, we maintain the basic approach applied by Pol et al. (2019) to obtain an estimate for the size of the Galactic population of pulsars, $N_{tot}$, and to subsequently infer the rate of DNS mergers to be detected by LIGO. We describe this method in the following to a detail that is sufficient to motivate and clarify our modifications that we explain in Section 4.

### 3.1 Statistical analysis towards the LIGO NS-NS merger detection rate

In order to estimate the total Galactic DNS population, the basic idea presented by Kim et al. (2003) and also applied by Pol et al. (2019) is simple. It is assumed that each currently detected DNS system originates from a separate part of a total DNS population that represents the individual current properties of the system and its visible pulsar. The systems are therefore treated independently and not, for example, as an evolved version of another observed DNS system. While this snapshot analysis circumvents the need for sophisticated binary evolution codes (which may bring uncertainties by themselves), the method still needs to reflect the selection effects, such as probability of the pulsar beam pointing towards Earth. Thus, the population size behind any detected DNS pulsar must be of a particular size, so that statistically exactly one pulsar is seen with the sensitivities of the conducted radio surveys. The statistical considerations leading to the population size that returns one detected pulsar are derived by Kim et al. (2003).

As $N_{tot} \gg N_{obs}$, it is expected that the probability of observing $N_{obs}$ pulsars out of an $N_{tot}$-sized population follows the Poisson distribution

$$P_{\text{Pois}}(N_{\text{obs}}; \lambda) = \frac{\lambda^{N_{\text{obs}}} e^{-\lambda}}{N_{\text{obs}}!}, \quad (2)$$

where, by definition, $\lambda := \langle N_{obs} \rangle$. By varying $N_{tot}$, Kim et al. (2003) also found that

$$\lambda = \gamma N_{\text{tot}} \quad (3)$$

with the constant $\gamma$ depending on the properties of the pulsar population. By knowing $\gamma$ one can derive $N_{tot}$ from $N_{obs}$.

The final goal is estimating the Galactic DNS merger rate from simulated populations, where in the first step, the likelihood function

<footnote n="1">¹The Double Pulsar has a merger time-scale of only 85 Myr compared to 300 Myr for J1906 + 0746.</footnote>

<footnote n="2">²https://github.com/devanshkv/PsrPopPy2</footnote>

<footnote n="3">³https://github.com/NihanPol/2018-DNS-merger-rate</footnote>
$P(\lambda|DX)$ is derived, i.e. the probability to obtain a model hypothesis, here the expectation value $\lambda$, given a data set of pulsars D and the model priors X (physical parameters of each pulsar, e.g. luminosity, pulse period etc.). Since the only probability known so far is the probability to observe a sample of pulsars D with given $\lambda$ and X, i.e.

$$P(D|\lambda X) = P_{\text{Pois}}(1; \lambda(N_{\text{tot}}), X) = \lambda(N_{\text{tot}})e^{-\lambda(N_{\text{tot}})}, \quad (4)$$

the Bayes' theorem is applied to rewrite $P(\lambda|DX)$ in terms of $P(D|\lambda X)$ as

$$P(\lambda|DX) = P(\lambda|X) \cdot \frac{P(D|\lambda X)}{P(D|X)} \quad (5)$$

Following the arguments given by Kim et al. (2003), the likelihood function of $\lambda$ then is given as

$$P(\lambda) = P(\lambda|DX) = P_{\text{Pois}}(1; \lambda(N_{\text{tot}}), X) = \lambda(N_{\text{tot}})e^{-\lambda(N_{\text{tot}})}. \quad (6)$$

From there, the probability distribution for $N_{\text{tot}}$ can easily be obtained as

$$P(N_{\text{tot}}) = P(\lambda) \left| \frac{d\lambda}{dN_{\text{tot}}} \right| = \gamma^2 N_{\text{tot}} \cdot e^{-\gamma N_{\text{tot}}}, \quad (7)$$

where equation (3) is used.

Furthermore, the merger rate $\mathcal{R}$ can be calculated via

$$\mathcal{R} = \frac{N_{\text{tot}}}{\tau_{\text{life}}} f_{\text{b}} \quad (8)$$

for a given population size $N_{\text{tot}}$, with the beaming correction factor $f_{\text{b}}$ and their lifetime $\tau_{\text{life}}$. Using this equation, the desired merger rate probability is given as

$$P(\mathcal{R}) = \left( \frac{\gamma \tau_{\text{life}}}{f_{\text{b}}} \right)^2 \cdot \mathcal{R} \cdot e^{-(\gamma \tau_{\text{life}}/f_{\text{b}})\mathcal{R}}. \quad (9)$$

In order to determine an estimate of the number of DNS merger events to be detected by the LIGO/Virgo network, at first the Milky Way merger rate $\mathcal{R}_{\text{MW}}$ is determined via fitting equation (9) to a given data set. Then, this rate is extrapolated to the observable volume of LIGO. To this end, one must estimate the formation rate of DNS systems in other galaxies.

In the nearby Universe, the formation rate of binary compact objects is expected to be proportional to the star formation rate, where a measure for the same is the B-band luminosity of the given galaxy<fnref n="4" /> (Phinney 1991; Kopparapu et al. 2008). Inside a sphere with radius $r$, the DNS merger rate follows

$$\mathcal{R}_{\text{LIGO}} = \mathcal{R}_{\text{MW}} \frac{L_{\text{total}}(r)}{L_{\text{MW}}} \quad (10)$$

with the total blue luminosity $L_{\text{total}}(r)$ within the distance $r$ and the B-band luminosity of the Milky Way $L_{\text{MW}} = 1.7 \times 10^{10} L_{\odot}$ denoted in terms of the solar B-band luminosity $L_{\odot} = 2.16 \times 10^{33} \text{ erg s}^{-1}$ (Kopparapu et al. 2008). The actual (Advanced) LIGO range<fnref n="5" /> in run O3 after the last upgrade is estimated to be $D_{\text{r}}^{\text{LIGO}} = 130 \text{ Mpc}$ (Abbott et al. 2017a; Abbott et al. 2018). Following the arguments laid out by Kopparapu et al. (2008), the rate of DNS mergers seen by LIGO is given by (cf. equation 19 in Kopparapu et al. 2008)

$$\mathcal{R}_{\text{LIGO}} = 7.4 \cdot 10^{-3} \left( \frac{\mathcal{R}_{\text{MW}}/L_{\text{MW}}}{(10^{10} L_{\odot})^{-1} \text{ Mpc}} \right) \times \left( \frac{D_{\text{h}}^{\text{LIGO}}}{100 \text{ Mpc}} \right) \text{yr}^{-1} \quad (11)$$

where $D_{\text{h}}^{\text{LIGO}} = 2.28 \cdot D_{\text{r}}^{\text{LIGO}}$ is the horizon distance<fnref n="6" /> of LIGO in run O3 (Finn & Chernoff 1993; Chen et al. 2017). In the derivation of the formula, the sensitivity decrease of LIGO depending on the orientation of the GW source with respect to the ground-based detector is already taken into account, as pointed out recently by Pol et al. (2020). Therefore, this correction to equation (15) in Pol et al. (2019) is also applied here.

## 3.2 PsrPopPy2

In order to infer the underlying Galactic radio pulsar population in the presence of selection effects, we follow the example of Pol et al. (2019) and deploy the package PsrPopPy2. The package provides a wide range of tools for the simulation and modelling of the pulsar populations and their evolution as described in Bates et al. (2014). The simulations are done within PsrPopPy2 creating a pulsar population with given pulsar parameters and then analyzing their detection during one or many chosen model surveys that reflect the real-life observing campaigns that led to the observed population. In the 'snapshot method', the aim is to determine $N_{\text{obs}}$, the number of pulsars detected by a survey, depending on the total number of pulsars in the population $N_{\text{tot}}$. The analysis is based on a Monte Carlo simulation: A number $N_{\text{tot}}$ of pulsars in the Galaxy is created, where each pulsar is given random parameters in luminosity and physical properties, drawn from a number of models distributions each (see Bates et al. 2014). After that, the population is run through a simulated radio telescope survey. Based on the specific survey parameters encoded in the PsrPopPy2 framework, it is checked how many pulsars could have been detected on Earth, which is then used as follows.

The goal is to obtain $\gamma$ of equation (3) for each kind of pulsar. Populations with numbers of pulsars, $N_{\text{tot}}$, varying between 10 and a few thousands are created and subjected to the simulated surveys, saving the number of observed pulsars, $N_{\text{obs}}$. That yields a data set, to which a Poisson curve (cf. equation 2) is fitted to obtain $\lambda$ for each $N_{\text{tot}}$. With this information, $\gamma$ is calculated for each pulsar by a linear fit to all pairs of $N_{\text{tot}}$, $\lambda$, following equation (3). The size of the expected galactic DNS (pulsar) population $N_{\text{real}}$, i.e. the population size $N_{\text{tot}}$, that yields one detection, is given by rearranging equation (3) and inserting $\lambda = \langle N_{\text{obs}} \rangle = 1$. Therefore, $N_{\text{real}}$ is given by $1/\gamma$. Then, $N_{\text{real}}$ and $\gamma$ are used to calculate the probability density function (PDF) for the population size and the merger rate following equations (7) and (9).

# 4 UTILIZING INFORMATION ON THE BEAM SHAPE OF PSR J1906 + 0746

## 4.1 Beam tomography

As relativistic spin-precession causes the viewing geometry of PSR J1906 + 0746 to change, our line of sight (LOS) cuts through the emission beam in different ways. This makes it possible to use the observed pulse profiles to construct a latitudinal cross-section of the pulsar beam. In the case of PSR J1906 + 0746, the geometry is such that the angle between the magnetic and the spin axes is measured to be $99.4 \pm 0.2^\circ$. As a result, for a period of about 20 yr, both magnetic poles were visible to an observer on Earth as two pulses separated from each other by half a spin period. Due to an inhomogeneity in the

<footnote n="4">⁴The B-band luminosity is the blue luminosity $L_B$ extracted from the galaxy's absolute blue magnitude $M_B$.</footnote>

<footnote n="5">⁵Radius of a Euclidean sphere containing the same volume as the true redshifted volume (Chen et al. 2017).</footnote>

<footnote n="6">⁶The horizon distance is the farthest luminosity distance a source could be detected above a certain threshold (Chen et al. 2017).</footnote>
<figure><img src="image_4.png" /></figure>

**Figure 1.** Latitudinal beam shape of the main pulse and interpulse of J1906 + 0476. Each point shows the maximum intensity of each beam in one of the 47 pulse profiles. The intensity is given in arbitrary units, as we are interested only in the relative intensities. $\beta_{MP}$ denotes the value of $\beta$ measured with respect to the magnetic axis of the main pulse; $\beta_{IP}$ denotes the value of $\beta$ measured with respect to the magnetic axis of the interpulse. Due to the precession of the pulsar, $\beta$ also can be read as a timeline from the most right point of every pulse to the left. Uncertainties in the shown values can result from calibration errors or interstellar scintillation. We estimate them for each data point to be of the order of 10–20 per cent, which we do not mark on this logarithmic scale.

pulsar beams, the pulses observed from each pole are not of equal or uniform intensity but in fact differ and change with time and relative intensity. The pole, the emission of which was stronger at discovery, is producing the ‘main pulse’ (MP), whereas the other pole causes the observed ‘interpulse’ (IP). As the LOS moved out of the ‘main’ beam, its observed intensity faded progressively, making the second pulse, the interpulse, the currently stronger one observed.

The angular separation of the LOS to the corresponding magnetic pole is designated as $\beta$, which is usually constant for ordinary pulsars but here becomes a function of time due to precession. The various observed profiles consequently provide a cross-section through the emission beams, similar to a tomography of the beam, across both magnetic poles. In order to describe the latitudinal structure, we measure the maximum intensity of each pulse (main pulse and interpulse, for each pole, separately) in every pulse profile and record these values for each of the 47 observations, freely available online<fnref n="7" /> (see Desvignes et al. 2019), from 2005 July 10 to 2018 June 21. We do this by splitting the pulse profile recorded in 2048 phase bins into two parts for main and interpulse and treat them independently hereafter. The geometry derived by Desvignes et al. (2019) from the Kramer & Wex (2009) precession model allows us to assign each observation an impact angle $\beta(t)$. Hereby, we introduce $\beta_{MP}$ and $\beta_{IP}$ to indicate that the corresponding values of $\beta$ are measured with respect to the magnetic axis of the main pulse (MP) or interpulse (IP), respectively. As shown in Fig. 1, the beams do not have a uniform shape but, in fact, we identify a number of sub-maxima in the main pulse. The apparent rapid variation in some of the intensity values may be caused by calibration errors or may be due to scintillation caused by the interstellar medium, although we cannot rule out that they are real, reflecting a fine structure of the beam. We will later

smooth over these apparent short-term features, as they are irrelevant for the further study.

Overall, these unique data allow us to conclude already that the general assumption of a uniform, box-shaped intensity distribution is not the best depiction of the beam shape. We also note that the interpulse shows a clear intensity drop around $\beta_{IP} = 0^{\circ}$, i.e. above the pole, which is expected from magnetospheric theories as pointed out by Desvignes et al. (2019). It is also worth emphasizing that the bright region of the main pulse is less wide than expected (see Introduction). Combined, with this new perception about the intensity distribution of the pulsar beam, we conclude that we can improve on the modelling done by Pol et al. (2019) by taking these observations into account.

## 4.2 Reconstructing the latitudinal beam shape

Even though relativistic spin-precession provides access to the latitudinal beam structure, the information is not complete. Due to gaps in the observations between 1998 and 2005, as well as between 2005 and 2009 (Desvignes et al. 2019), a reconstruction requires a limited set of assumptions and extrapolation from the available data. This is done taking two conditions into account: First, the pulsar has been detected with a strong signal in 1998 (Lorimer et al. 2006), where $\beta_{MP}$ was greater than $0^{\circ}$ (Desvignes et al. 2019). Secondly, assuming that for an almost orthogonal rotator as we have here ($\alpha \sim 90^{\circ}$), we can expect to first order an axial symmetry with respect to the magnetic axis. Based on these assumptions, we extrapolate the beam in the form of mirroring intensity profiles of both main pulse and interpulse along the axis with $\beta_{MP/IP} = 0^{\circ}$, the results of which are shown in Fig. 2.

Finally, in order to obtain a continuous intensity distribution as a function of $\beta$ for our later modelling, we perform a spline interpolation. The lack of information about the main pulse’s beam profile in the range of $\beta_{MP} \in [-5^{\circ}, 5^{\circ}]$ poses a major problem. We apply the LSQUnivariateSpline function from the PYTHON package numpy.interpolate, because it returns a spline function with the option of setting explicit internal knots, so the edges can be shaped manually. The different shapes of the main pulse and the interpulse impose a slightly different interpolation process for both beams.

Starting with the main pulse, the values for $|\beta_{MP}| > 15^{\circ}$ show short-term oscillations as discussed before (cf. Fig. 1). Hence, for this range, we set grid points manually. It turns out that grid points at $\pm 16^{\circ}$, $\pm 17^{\circ}$, $\pm 18^{\circ}$, $\pm 19^{\circ}$, and $\pm 20^{\circ}$ yield a stable spline that represents the shape of the profile. Obviously, the distinct structures in the outer region of the profile can cause an instability of the interpolation, strongly depending on the position of the grid points. This restricts the number and the position of grid points. To model the pulse around $\beta_{MP} \in [-10^{\circ}, 10^{\circ}]$, we identify the grid points in this interval with the actual data points to ensure that the spline follows the explicit shape of the measured beam as well as ensures a smooth interpolation of the profile in $\beta_{MP} \in [-5^{\circ}, 5^{\circ}]$. We are uncertain, as to whether the intensity peak at $\beta_{MP} \sim -10^{\circ}$ is a real feature of the pulse profiles. Judging from the slopes of the intensity profile on either side, we are inclined to consider it as real and will therefore include it in our modelling of the beam shape.

We applied the same scheme to the interpulse. The edge of the observed beam shape was modelled using manually chosen grid points (set to $-2^{\circ}$ to $3^{\circ}$ on the left side and $35^{\circ}$ to $40^{\circ}$ on the right side, both in steps of $1^{\circ}$.) We show the results in Fig. 2, where the composition of both interpolated latitudinal beam shapes is depicted.

<footnote n="7">⁷https://zenodo.org/record/3358819</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Similar to Fig. 1, showing the latitudinal beam intensity profile for both the main (blue) and the interpulse (green) but performing a spline interpolation and mirroring the observed data at the location of the magnetic axis at each pole. For each pole, the darker shaded data points show the maximum intensities of the measured profiles. The lighter shaded data points indicate the mirrored values. The interpolating function used for each pulse is a cubic spline. $\beta_{MP}$ is measured with respect to the main pulse magnetic axis, the intensity is given in arbitrary units, as we are interested only in relative intensities. The location of both magnetic axes is marked by yellow lines. The red line marks the pulsar's (rotational) equator, i.e. $\alpha = 90^{\circ}$.</figcaption></figure>

## 4.3 Implementation in the code framework

As shown in Fig. 2, the measured latitudinal main pulse and interpulse beam shapes differ significantly in their shape and relative intensities. This is in sharp contrast to the assumptions in previous studies, where one would compute the beaming correction factor based on (measured or estimated) values for the angles $\alpha$ and $\rho$ according to

$$ \frac{4\pi}{f_b} = 2 \times \int_{0}^{2\pi} d\phi \int_{\max(0, \alpha-\rho)}^{\min(\pi/2, \alpha+\rho)} \sin\theta d\theta, \quad (12) $$

which imposes the very simple beaming model of two equal beams with hard-edged cones (cf. O'Shaughnessy & Kim 2010).

In the case of PSR J1906 + 0746, we can use the determined beam information directly. The overlap of both profiles (main and interpulse) as well as the relative intensity to each other pulse needs to be handled with care, however, as it does not matter for detection statistics, if the pulsar is detected by its main pulse, its interpulse, or both. For this reason, we construct an effective beam shape that reflects the illuminated sky, using the higher intensity of a given beam if the same spot is illuminated twice.

As we have modelled the intensity in a logarithmic scale to trace the beam edge more accurately, the intensity values obtained from the splines are exponentiated and afterwards normalized with respect to the maximum value of both poles, i.e. the maximum value of the main pulse spline.

In our simulations, PSR J1906 + 0746-like pulsars will be assigned a luminosity that is drawn from an assumed distribution describing the whole population (see Pol et al. 2019). We also draw a random value of $\zeta \in [0^{\circ}, 180^{\circ}]$, which is uniformly distributed in $\cos(\zeta)$, from which we derive $\beta_{MP}$ as $\zeta = \alpha + \beta$. The observable radio intensity is that of the model luminosity scaled to a new value according to the value of $\beta_{MP}$ and the inferred relative latitudinal intensity.

By choosing a value of $\zeta$ along the whole polar angle, the fraction of the sky that is illuminated by the pulsar's beam is already intrinsically considered within our simulations. Hence, when

later analyzing the results from the PSRPOPPY simulations (see Section 3.2), the beaming correction factor in equation (9) is set to $f_b = 1$.

## 5 A GENERIC LATITUDINAL BEAM SHAPE FOR OTHER PULSARS

PSR J1906 + 0746 is not the only DNS system where relativistic spin precession has been observed. Information on the beam structure is also available for the DNSs B1913 + 16 (Kramer 1998; Weisberg & Taylor 2002) and B1534 + 12 (Stairs, Thorsett & Arzoumanian 2004) and the relativistic NS-white dwarf system PSR J1141 – 6545 (Manchester et al. 2010; Venkatraman Krishnan et al. 2019). In none of these cases, however, has the beam been traversed completely yet, so that the important latitudinal extent cannot be determined as for PSR J1906 + 0746. [Only PSR J0737 – 3039B, the slow companion in the Double Pulsar (Lyne et al. 2004) has precessed out of our LOS (Perera et al. 2010), but here the beam is severely distorted by wind of pulsar A (McLaughlin et al. 2004)]. What is common to all beam shapes observed so far [including that of the non-recycled pulsar J1141 – 6545 (Manchester et al. 2010; Venkatraman Krishnan et al. 2019)] is that they do not seem to have a uniform structure, as it is assumed in all previous studies to infer DNS merger rates from the known population. We therefore improve on this assumption by utilizing the observed beam structure of PSR J1906 + 0746 to derive a more realistic template for a beam shape of other DNSs to provide a new estimate on the DNS merger rate. Even though the individual beam shapes may in reality still differ from this template, a study of the variations in the resulting merger rates allows us to quantify the corresponding systematic uncertainty related to the beam shape for the first time.

In order to derive a generic pulse profile, we consider that we have observed different parts of the beam in both main pulse and interpulse of PSR J1906 + 0746. A suitable combination of these well-measured portions in both beams describes the full latitudinal extent
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Exemplary composition of the generic latitudinal beam profile for one pulsar beam based on the observation of PSR J1906 + 0746 together with a polynomial fit function. Points marked in green are taken from the main pulse measurements and those marked blue belong to the interpulse measurements, shifted by a random offset of 0.45 in intensity. The fitted and normalized polynomial representing the corresponding implemented profile is shown. The shaded region indicates the range of all possible other fitting functions for all different offsets. It is limited by the fit functions with the minimal and maximal offset of 0.0 and 0.9 for the interpulse points.</figcaption></figure>

that is crucial for a realistic simulation of the pulsar detectability in a mock survey. We choose the edge of the main pulse as the edge of the generic beam shape. We use the centre of the interpulse to model the central part of the generic pulse shape. In order to account for the gaps in the beam coverage as observations were not regularly spaced (cf. Fig. 1), we make use of the spline-fitting discussed in Section 4.2 and shown in Fig. 2.

The pulsar beam size, and hence its latitudinal extent, increases with decreasing period, *P*. As discussed in the introduction, for dipolar field lines, one expects a scaling of $\rho = k \times P^{-0.5}$, which is confirmed statistically from observations. Given that the factor *k* has some uncertainty (see introduction and Venkatraman Krishnan et al. 2019), we adopt the value of *k* = 5.4 deg $s^{0.5}$ to be consistent with O'Shaughnessy & Kim (2010). Computing the corresponding $\rho$ for PSR J1906 + 0746, this implies a value of $\rho \sim 14^{\circ}$. This appears to be somewhat too small when inspecting Fig. 1, but the intensity of the main pulse has dropped by a factor of $\sim 100$ at this value of $\beta$, so that we consider this as a sufficiently good approximation, especially since all DNSs studied here have $P > 20$ ms. We show the resulting latitudinal beam shape after scaling the $\beta$-axis by a corresponding $\rho$ value in Fig. 3.

Before we apply the generic beam shape to simulated pulsars, we scale the profile in Fig. 3 to the appropriate pulse period. We also vary it in a random fashion in two aspects. One is resulting from the uncertainty in the intensity level of the clear minimum in the profile's centre upon the magnetic axis. Here, for each simulated pulsar, we vary the relative intensity at the centre of the beam by a random value in the range [0.,0.9]. The maximum value of 0.9 is chosen in order to preserve at least a small intensity drop, as this is also expected theoretically (Timokhin & Arons 2013; Gralla, Lupsasca & Philippov 2017).

The other aspect addresses the relative intensity of the beams above two magnetic poles. From observations of interpulse pulsars (e.g. Johnston & Kramer 2019) and from an apparent lack of interpulse pulsars in the population of normal pulsars (e.g. Weltevrede & Johnston 2008), we deduce that the relative intensity between the beams may be remarkably different and may depend

<table><caption>**Table 1.** Values for the magnetic inclination angles applied in the implementation of the generic pulse profiles.</caption><thead><tr><td>Pulsar</td><td>α (deg)</td><td>Reference</td></tr></thead><tbody><tr><td>J0737 − 3039A</td><td>79</td><td>Kramer et al. (2021)</td></tr><tr><td>B1534 + 12</td><td>103</td><td>Stairs et al. (2004)</td></tr><tr><td>J1756 − 2251</td><td>74</td><td>Ferdman et al. (2014)</td></tr><tr><td>J1906 + 0746</td><td>81</td><td>Desvignes et al. (2019)</td></tr><tr><td>B1913 + 16</td><td>153</td><td>Kramer (1998)</td></tr></tbody></table>

on some a priori unknown physical parameter. We therefore assign a random intensity ratio determined from two random draws from the lognormal luminosity distribution that is used in the simulations. By definition, we call the resulting less luminous pulse the interpulse.

Finally, we assign every simulated pulsar a geometry, which consists of the angles $\alpha$ and $\beta$. For pulsars representing the realizations of DNSs with known values of $\alpha$, we adopt those listed in Table 1. As $\alpha$ is yet unknown for the pulsars J0509 + 3801, J1757 − 1854, J1913 + 1102, and J1946 + 2052, for each simulated pulsar in the associated population we choose a random value distributed uniformly in $\cos(\alpha)$, consistent with $\zeta$ being uniformly distributed in $\cos(\zeta)$. In order to automatically account for the fraction of the sky illuminated by either main pulse or interpulse, we construct a combined profile as we had demonstrated for PSR J1906 + 0746 in Fig. 2, i.e. we locate the main pulse at a position $\beta_{MP} = 0^{\circ}$ and centre the weaker interpulse accordingly at $\beta = 2 \times (90^{\circ} - \alpha)$.

# 6 RESULTS

In the following, we apply our derived latitudinal beam shapes to simulations of the Galactic DNS population and computations of the inferred DNS merger rate as an enhancement of the work by Pol et al. (2019) and Pol et al. (2020). We do this in three steps. After first confirming that we can reproduce the results by Pol et al., we then modify the simulations by using only the new information on the beam shape for PSR J1906 + 0476 to gauge the impact of this change. In a third step, we then apply the newly derived generic latitudinal beam shape to all DNS systems, which allows us to compare the different resulting merger rates and to assess the importance of these changes.

## 6.1 Validating the simulation scheme

Before applying the derived beam shape model, we use the information provided by Pol et al.<fnref n="8" /> to compare our obtained results with those derived by Pol et al. (2019, 2020) (see Table 2). For most pulsars, we could use the actual simulation and parameter files used by Pol et al. for a detailed comparison. Only for PSR J0509 + 3801 added in Pol et al. (2020) this was not possible, and we derived our own input data. As a result, we derive a Galactic merger rate estimate of $R_{\text{MW,recr}}^{\text{Pol}} = 37_{-11}^{+24} \text{ Myr}^{-1}$ and an Advanced LIGO detection rate of $R_{\text{LIGO,recr}}^{\text{Pol}} = 3.9_{-1.2}^{+2.4} \text{ yr}^{-1}$, which are a little smaller but in very good agreement with the values presented in Pol et al. (2020). The small deviation in the result is likely caused by statistical fluctuations in the simulation framework, the limited number of simulation runs in both cases, or some differences in the input data (e.g. the estimated lifetime of PSR J0509 + 3801). In the following, we use our derived value as the reference to compare to the changes when introducing improved beam shape modelling.

<footnote n="8">⁸https://github.com/NihanPol/2018-DNS-merger-rate/tree/master/python</footnote>
<table><caption>**Table 2.** Parameters and results from the simulation using the values of $f_b$ adopted by Pol et al. (2019).</caption><thead><tr><th>Pulsar</th><th>f<sub>b</sub></th><th>δ</th><th>τ<sub>age</sub><br>(Myr)</th><th>t<sub>merger</sub><br>(Gyr)</th><th>γ</th><th>N<sub>tot</sub></th><th>R<br>(Myr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>J0509 + 3801</td><td>4.59</td><td>0.18</td><td></td><td></td><td>0.007437</td><td>614<sup>+2823</sup><sub>-422</sub></td><td>0.8<sup>+3.7</sup><sub>-0.6</sub></td></tr><tr><td>J0737 − 3039A</td><td>2.0</td><td>0.27</td><td>159</td><td>0.085</td><td>0.003183</td><td>627<sup>+2868</sup><sub>-463</sub></td><td>2.6<sup>+11.9</sup><sub>-1.9</sub></td></tr><tr><td>J0737 − 3039B<sup>9</sup></td><td></td><td></td><td></td><td>0.085</td><td></td><td></td><td></td></tr><tr><td>B1534 + 12</td><td>6.0</td><td>0.04</td><td>208</td><td>2.70</td><td>0.010687</td><td>567<sup>+2581</sup><sub>-386</sub></td><td>0.2<sup>+0.9</sup><sub>-0.1</sub></td></tr><tr><td>J1756 − 2251</td><td>4.59</td><td>0.03</td><td>396</td><td>1.69</td><td>0.009133</td><td>504<sup>+2292</sup><sub>-348</sub></td><td>0.3<sup>+1.4</sup><sub>-0.2</sub></td></tr><tr><td>J1757 − 1854</td><td>4.59</td><td>0.06</td><td>159</td><td>0.076</td><td>0.006522</td><td>706<sup>+3208</sup><sub>-495</sub></td><td>4.4<sup>+19.8</sup><sub>-3.1</sub></td></tr><tr><td>J1906 + 0746 (old)</td><td>4.59</td><td>0.01</td><td>0.11</td><td>0.30</td><td>0.019134</td><td>248<sup>+1082</sup><sub>-147</sub></td><td>4.1<sup>+18.0</sup><sub>-2.4</sub></td></tr><tr><td>J1906 + 0746 (new)</td><td>1</td><td>0.01</td><td>0.11</td><td>0.30</td><td>0.002781</td><td>361<sup>+1638</sup><sub>-266</sub></td><td>6.0<sup>+27.3</sup><sub>-4.5</sub></td></tr><tr><td>J1913 + 16</td><td>5.7</td><td>0.169</td><td>77</td><td>0.50</td><td>0.006681</td><td>651<sup>+2998</sup><sub>-458</sub></td><td>2.3<sup>+10.6</sup><sub>-1.7</sub></td></tr><tr><td>J1913 + 1102</td><td>4.59</td><td>0.06</td><td>2625</td><td>0.30</td><td>0.007019</td><td>816<sup>+3777</sup><sub>-586</sub></td><td>0.2<sup>+1.0</sup><sub>-0.1</sub></td></tr><tr><td>J1946 + 2052</td><td>4.59</td><td>0.06</td><td>247</td><td>0.046</td><td>0.004800</td><td>963<sup>+4363</sup><sub>-697</sub></td><td>3.3<sup>+14.9</sup><sub>-2.4</sub></td></tr></tbody></table>

Notes. Note 1: $f_b$ is the beaming correction factor, $\delta$ is the pulse duty cycle, and $\tau_{age}$ is the effective age of the pulsar. The numbers given for $N_{det}$, $N_{tot}$, and $R$ denote the peak values resulting from the probability distributions 7 and 9. $N_{det}$ is the number of pulsars beaming towards the earth, $N_{tot}$ is the total population number, and $R$ is the individual merger rate. $N_{tot}$ is gained by scaling $N_{det}$ with $f_b$.

Note 2: To remain consistency with the work by Pol et al. (2019), all errors are quoted on the 95 % confidence intervals.

Note 3: PSR J1906 + 0746 (old) denotes the results of the simulation based on the concept by Pol et al.; J1906 + 0746 (new) denotes the results of the simulation using the spline interpolated latitudinal intensity profile.

## 6.2 Applying beam shape modelling for PSR J1906 + 0746

Following the prescription described in Section 4.3, the simulations are repeated for all DNS sources but with modifications for J1906 + 0746. We also apply the correction discussed in context of equation (11) and a horizon distance of 130 Mpc. Following Pol et al. (2020), we also include the recently discovered DNS J0509 + 3801 in our sample and include the corresponding Green Bank North Celestial Cap Survey in our simulations (Lynch et al. 2018).

The results of our computations are shown in Figs 4 and 5, where we can study the impact of introducing a more realistic beam shape model. The impact of the new method is clearly seen in the increase of $N_{det}$ from $56_{-31}^{+236}$ to $N_{det}^{new} = 361_{-266}^{+1638}$. But note that in our new scheme, the selection effect due to beaming is now fully addressed in our survey simulations, so that for PSR J1906 + 0746, we set $f_b = 1$ in equation (8); therefore, $N_{det}^{new} = N_{tot}^{new}$. Thus, the peak value of the total population number increased very significantly by 45 per cent. The effect of the beam shape correction on the individual merger rate of J1906 + 0746 is shown in Fig. 4: it shifts from $4.1_{-2.4}^{+18.0}$ Myr$^{-1}$ to $R^{new} = 6.0_{-4.5}^{+27.3}$ Myr$^{-1}$. This increase also affects the total Milky Way merger rate and the calculated LIGO merger rate. The new rate distributions are shown in both plots in Fig. 5. The updated estimate on the Milky Way merger rate is $R_{MW}^{new} = 38_{-12}^{+27}$ Myr$^{-1}$. The new DNS merger rate estimate for LIGO resulting from the changes applied to PSR J1906 + 0746 is

$$R_{\text{LIGO}}^{\text{new,1906}} = 4.1_{-1.32}^{+2.98} \text{yr}^{-1}, \qquad (13)$$

which imposes indeed a significant increase. The errors of both values are quoted at the 90 per cent confidence level.

## 6.3 Effective beaming correction factor

Using the PDFs of the merger rates generated from the generic latitudinal profile simulation, it is possible to determine an effective beaming correction factor $\tilde{f}_{b,eff}$. We define it as the value for $f_b$, at which the merger rate PDF obtained from the simulation following the procedure from Pol et al. (2019) coincides with the merger rate PDF from the simulation using generic pulse profile, i.e. $\tilde{f}_{b,eff} = \gamma_{old}/\gamma_{new}$. $\tilde{f}_{b,eff}$ is iteratively determined by imposing a maximal difference between the maxima of the two PDFs. The individual results are shown in the last column of Table 3. For PSR J1906 + 0746, we apply an alternative way to compute $\tilde{f}_{b,eff}$ by using our spline interpolation of the beam shape directly. This results in $\tilde{f}_{b,eff} = 6.98$. This is significantly larger than the value adopted by O'Shaughnessy & Kim (2010) ($f_b,eff = 3.37$) due to the lack of additional information, or the uniform value adopted by Pol et al. (2019) ($f_b,eff = 4.60$). The impact of adopting an appropriate beam shape is clearly demonstrated.

## 6.4 General application of the generic latitudinal beam shape

In the last step, we adopt the generic beam shape mode for all pulsars listed in Table 3 following the procedure described in Section 5. Again, for the extrapolation to the LIGO DNS detection rate, we use a horizon distance of 130 Mpc. The numerical results are listed in Table 3. Comparing the values of $N_{det}$ from the previous simulation to the simulation with the generic profile, we first see a significant increase in the number of all DNS systems. However, when computing the merger rates [hence using $f_b = 1$ in equation (8) for all pulsars], most contributions to the overall rates decreased. The individual merger rates following from the simulation are shown in Fig. 6. Inspecting the location of the peaks in the DNS merger rates, compared to the earlier reproduction of the Pol et al. results (Fig. 4), we can see that most of the curves moved to slightly lower values. Not surprisingly, we encounter the same grouping of curves: The three pulsars B1534 + 12, J1913 + 1102, and J1756 - 2251 still contribute

<footnote n="9">⁹The J0737-3039 system was discovered with pulsar A; also, pulsar B will cross the death line significantly earlier than pulsar A. In addition to that, pulsar B introduces large uncertainties into the total merger rate (Pol et al. 2019). Therefore, it will not be taken into account for the analysis.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Galactic merger rate PDFs for the individual DNS systems. J1906 + 0746 is simulated using the extrapolated latitudinal beam profile (see Fig. 1). As result, the PDF for the J1906 + 0746 population extends to the largest rates in comparison.</figcaption></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Estimates on the Milky Way DNS merger rate (top plot) and the DNS merger detection rate by LIGO (bottom plot) using the latitudinal beam profile (Fig. 1) for PSR J1906 + 0746. The errors are quoted at the 90 per cent confidence interval, with 45 per cent limits to the left and the right.</figcaption></figure>

least to the overall DNS merger rate, but yet the other PDF curves appear closer to each other. The PDF from the PSR J0509 + 3801 population merger rate moved into the right group of curves and also the position of the individual PDF curves changed within that group. Importantly, PSR J1906 + 0746 represents the DNS population with the biggest contribution to the galactic DNS merger rate, which is indeed to be expected, given its young age. The contributions of PSRs J0737 – 3039A and J1913 + 1102 increased noticeably in relation to the other curves. The resulting estimate on the Milky Way DNS merger rate is $R_{\text{MW}}^{\text{gen}} = 32_{-9}^{+19} \text{ Myr}^{-1}$. The estimate for the LIGO detection rate becomes

$$R_{\text{LIGO}}^{\text{gen}} = 3.5_{-1.0}^{+2.1} \text{ Myr}^{-1}, \qquad (14)$$

which represents a decrease in comparison to the results of Pol et al. (2020), but is still consistent with their values. Again, the errors are quoted at the 90 per cent confidence level.

## 7 DISCUSSION

There are a number of caveats and systematic uncertainties in the applied method and framework. The caveats concerning the applied lognormal luminosity distribution used in the PsrPopPy2 simulation, the beaming correction factors, the effective lifetime of J1906 + 0746, and the extrapolation of the merger rate to the observable volume of LIGO outlined in Section 4.2 of the paper by Pol et al. (2019) are equally valid in this work. However, since they have been discussed in detail in the cited work, we mention them only for completeness at this point and refer to Pol et al. (2019) for more details.

One may also consider whether all systems studied here represent ‘true’ DNS systems, or whether an unseen companion may be indeed a heavy white dwarf. In principle, such possibility cannot be excluded. In the Double Pulsar, we have seen both components of the system. In the case of PSR J1906 + 0746, we observe the non-recycled pulsar. It is possible to find a young pulsar around a white
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Galactic merger rate PDFs for the individual DNS systems based on the simulation using the generic latitudinal beam profile. We still encounter the accumulation of two distinct groups in terms of the PDF's relative position, yet overall the curves shifted to the left, i.e. to smaller peak positions.</figcaption></figure>

<table><caption>**Table 3.** Parameters and results from the simulation using the generic latitudinal pulse profile. The errors are quoted at the 90% confidence interval.</caption><thead><tr><th>Pulsar</th><th>γ</th><th>N<sub>det, tot</sub></th><th>R (Myr<sup>-1</sup>)</th><th>f<sub>b, eff</sub></th></tr></thead><tbody><tr><td>J0509 + 3801</td><td>0.00138568</td><td>720<sup>+2700</sup><sub>-543</sub></td><td>0.9<sup>+3.6</sup><sub>-0.7</sub></td><td>5.24</td></tr><tr><td>J0737 – 3039A</td><td>0.0012233</td><td>817<sup>+2752</sup><sub>-619</sub></td><td>3.4<sup>+11.5</sup><sub>-2.6</sub></td><td>2.47</td></tr><tr><td>B1534 + 12</td><td>0.00334478</td><td>298<sup>+1366</sup><sub>-220</sub></td><td>0.1<sup>+0.5</sup><sub>-0.1</sub></td><td>3.20</td></tr><tr><td>J1756-2251</td><td>0.00332935</td><td>302<sup>+1370</sup><sub>-224</sub></td><td>0.2<sup>+0.8</sup><sub>-0.1</sub></td><td>2.76</td></tr><tr><td>J1757 – 1854</td><td>0.0022562</td><td>441<sup>+2001</sup><sub>-328</sub></td><td>2.7<sup>+12.4</sup><sub>-2.0</sub></td><td>2.87</td></tr><tr><td>J1906 + 0746</td><td>0.00350737</td><td>286<sup>+1302</sup><sub>-208</sub></td><td>4.8<sup>+21.7</sup><sub>-3.5</sub></td><td>5.50</td></tr><tr><td>J1913 + 16</td><td>0.00084191</td><td>1188<sup>+2604</sup><sub>-919</sub></td><td>3.2<sup>+7.0</sup><sub>-2.5</sub></td><td>6.71</td></tr><tr><td>J1913 + 1102</td><td>0.00217771</td><td>457<sup>+2065</sup><sub>-339</sub></td><td>0.1<sup>+0.7</sup><sub>-0.1</sub></td><td>3.21</td></tr><tr><td>J1946 + 2052</td><td>0.00187952</td><td>533<sup>+2325</sup><sub>-399</sub></td><td>1.8<sup>+7.9</sup><sub>-1.4</sub></td><td>2.55</td></tr></tbody></table>

*Notes.* Note 1: *f*<sub>b, eff</sub> is the effective beaming correction factor calculated as described in Section 6.3. The numbers given for *N*<sub>det, tot</sub> and *R* denote the peak values resulting from the probability distributions 7 and 9. *N*<sub>det</sub> is the number of pulsars beaming towards the earth, *N*<sub>tot</sub> is the total population number, and *R* is the individual merger rate. *N*<sub>det</sub> with *f*<sub>b</sub> = 1 for all pulsars, the distributions do not change from *N*<sub>det</sub> to *N*<sub>tot</sub> and thus the values are given in the same column.

*Note 2:* To remain consistency with the work by Pol et al. (2019), the errors on the population numbers are quoted on the 95% confidence intervals, whereas the errors on the merger rates are given with the 90% confidence interval.

dwarf companion, as in the case of PSR J1141 – 6545 (Antoniadis et al. 2011), but such evolution requires finely tuned initial system parameters (e.g. Tauris & Sennels 2000). In all other cases, the observed pulsar is recycled, and all systems have a significant eccentricity, which indicates a second supernova explosion with an asymmetric kick and/or significant mass loss, strongly suggesting that these systems are DNSs.

It has recently been suggested that a population of massive radio-quiet neutron stars in compact binary systems could exist (Vigna-

Gómez et al. 2021). In such a case, our derived LIGO detection rates will represent a lower limit. If there is indeed a large difference to the observed detection rates, compared to our estimates here, it would lend credibility to such conclusion. However, one should note that so far there is no evidence for a correlation between pulsar luminosity and pulsar mass (Lorimer & Kramer 2005).

In the following, we discuss further caveats associated with the specifics of this work.

## 7.1 Caveats on the latitudinal beam profile

### 7.1.1 Extrapolation from the PSR J1906 + 0746 measurements

The main uncertainty in deriving the latitudinal beam profile arises from the extrapolation of the main pulse. A fraction less than half of both latitudinal profiles is covered by the available observations and only the detection in 1998 gives information back in time to tie the extrapolation to. As no further measurements are available, the method of mirroring the beam shape is the only way to treat the pulse intensity profile. Therefore, the actual shape of the latitudinal profile could vary significantly from the one assumed in this work. During our simulations, we have seen that a variation of the opening angle strongly affects the merger rate estimate, so that a different pulse shape could lead to overall different results in the individual contributions of the populations to the combined galactic merger rate. In addition, the lack of data points for both main pulse and interpulse from 2009 to 2012 forces us to interpolate the resulting wide gap nearly linearly as first-order approximation. The real pulse shape may be different in this gap. As discussed previously, also the shape of the intensity distribution in the centre of each pulse poses an uncertainty to our merger rate estimate for the population of PSR J1906 + 0746. A more accurate profile could lead to a larger illuminated fraction of the sky and thus to a small decrease of the DNS merger rate estimate.
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Comparison of the merger rates derived for PSR J1906 + 0746 using the variation of methods discussed in this work. The curves correspond in their order as they appear in the labelling to the different simulation strategies (i)–(iv) enumerated in Section 7.2.</figcaption></figure>

### 7.1.2 Validity of the generic beam shape

The derivation of the generic profile from the profile of PSR J1906 + 0746 also bears some uncertainties, which can affect the resulting populations and thus the DNS merger rate estimate significantly.

We chose the contribution of both magnetic poles of PSR J1906 + 0746 to the generic profile in such a way that the resulting polynomial depicts the observed shape the most suitable way. But surely we could approach this in many different ways, especially the choice of the individual contributions. For example, one could add up both profiles and fit a suitable function accordingly. But taking a closer look at Fig. 1, this would lead to a smeared shape, as the edges of both profiles do not coincide. Compensating for that by rescaling parts of either the main pulse or the interpulse along the x-axis, one would have had to make alternative assumptions about the rescaling magnitude and also which part to rescale.

Yet another way of gaining a generic profile would have been working only with the interpulse data, as via mirroring the values, ending up with an almost complete profile. This option also hides two peculiarities. At first, the value of $\beta_{IP}$ corresponding to the maximum of the interpulse is larger than the value of $\beta_{MP}$ corresponding to the maximum of the main pulse. Therefore, taking the two main peaks from the latitudinal interpulse profile (cf. Fig. 1) as maxima for the generic pulse profile would lead to a slightly larger latitudinal extent of the same. This in turn leads to smaller population sizes, probably resulting in an underestimation of the contribution to the DNS merger rate. Especially, pulsar populations with large pulse widths are prone to that. Secondly, in order to take the centre variations of the profile into account, one would have to decide arbitrarily, which data points define the 'centre' of the profile.

Taking all this into account, we decided following the approach pointed out in Section 5. The most uncertain aspect of that strategy is the negligence of the plateaus near the middle of the interpulse, motivated by the circumstance that they caused numerical instabilities in the fitting process. The intensity level of the central values of the intensity profile above the pole, as well as the relative intensity of interpulse and main pulse, also bears some uncertainty. However, given that PSR J1906 + 0746 is the only pulsar for which such a measurement has been possible so far, we believe that our approach serves as a useful first-order approximation, improving upon the assumption of a uniformly filled beam, which clearly contradicts the observations.

## 7.2 Comparison of methods for PSR J1906 + 0746

Our main focus in this work is the implications of the recent observational results for PSR J1906 + 0746. In addition, for PSR J1906 + 0746, we have the unique opportunity to study the impact corresponding to the additional degrees of freedom introduced during the creation of the generic profile. This leads to a total of four different simulations on PSR J1906 + 0746 that we conducted, i.e. those

(i) following the strategy of Pol et al. (2019);

(ii) applying the profile extracted from the observations;

(iii) applying the generic profile<fnref n="10" /> using the known (fixed) $\alpha$; and

(iv) applying the generic profile<fnref n="11" /> and assuming $\alpha$ to be unknown, i.e. allowing it to vary.

The results of all four simulations are shown in Fig. 7; the values of $\gamma$ and $N_{obs}$ are collected in Table 4. We see that the rates resulting from simulations (ii)–(iv) are shifted to the right with respect to the merger rate obtained from the initial simulation by Pol et al. (2020), meaning that in all three cases the peak merger rate increased. Nevertheless, we also find that the usage of the generic profile leads to a noticeably smaller increase than the usage of the observed profile.

Obviously, allowing the actually observed profile to vary in a random fashion when applying the generic profile has a noticeable impact. Due to the variation in the intensity ratio between main and interpulse, as well as varying the central intensity levels, we create, on average, a pulsar beam that is brighter than the original one. As a result, a corresponding population of pulsars is detected more often in our virtual pulsar surveys than otherwise, leading to a decrease in the merger rate estimate. Overall, however, we can conclude that it is crucial for the simulation to take the latitudinal beam shape into consideration.

## 7.3 Comparison with other DNS merger rate estimates

The estimate for the rate of DNS merger detections by LIGO derived in this work can be compared to the ones obtained through different theoretical or phenomenological models. In comparison to the method used in this work, i.e. extrapolating the merger rate on the basis of the detected DNS systems, it is also possible to generate a population of Galactic DNS systems ab initio. Here, the different (stellar and binary evolution) formation processes towards a DNS system are considered. There is a rich literature on this topic. Here, as a recent example, we compare our results to the rate predictions by Chruslinska et al. (2017) and by Kruckow et al. (2018) shown in Fig. 8. We also consider the DNS merger rate calculated by the LIGO collaboration based on the observational input from GW170817 and GW190425.

### 7.3.1 LIGO DNS merger rate

Due to the detection of two DNS mergers, GW170817 and GW190425, the LIGO collaboration released a new estimate on the DNS merger rate based on both these events (Abbott et al. 2020). After the unit conversion as in Pol et al. (2020), we find

$$R_{\text{LIGO}} = 4.6^{+7.1}_{-3.4} \times \left( \frac{D_r}{100 \text{ Mpc}} \right)^3 \text{yr}^{-1}.$$

<footnote n="10">¹⁰Variation of the centre intensity values and the intensity ratio between the MP and IP allowed.</footnote>

<footnote n="11">¹¹See footnote 10.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Comparison of the DNS merger rate estimates calculated in this work to previous work using different methods. These are the ab initio estimates by Chruslinska et al. (2017) and Kruckow et al. (2018), the LIGO estimate by Abbott et al. (2020), and the estimate from the simulation following the Pol et al. (2020) scheme (which is slightly smaller here than the value derived by Pol et al. 2020 due to statistical fluctuations). Both estimates derived in this work are consistent with each other, as well as with the rate from Pol et al. (2020). All rates are scaled to an observational distance of LIGO of 130 Mpc.</figcaption></figure>

<table><caption>**Table 4.** Collection of all numerical results concerning the simulations on J1906 + 0746. These values are used to estimate the errors on the simulation of the other pulsars with unknown physical beam shape.</caption><thead><tr><th></th><th>N<sub>tot</sub></th><th>R (Myr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Pol et al. (2019, 2020)</td><td>248<sup>+1082</sup><sub>-147</sub></td><td>4.13<sup>+18.03</sup><sub>-2.44</sub></td></tr><tr><td>Measured profile</td><td>365<sup>+1669</sup><sub>-272</sub></td><td>6.09<sup>+27.82</sup><sub>-4.53</sub></td></tr><tr><td>Generic profile, fixed α</td><td>286<sup>+1306</sup><sub>-208</sub></td><td>4.76<sup>+21.77</sup><sub>-3.46</sub></td></tr><tr><td>Generic profile, random α</td><td>357<sup>+1630</sup><sub>-264</sub></td><td>5.96<sup>+27.16</sup><sub>-4.39</sub></td></tr></tbody></table>

Applying a luminosity distance of $D_r = 130$ Mpc, this gives a rate of $R_{LIGO} = 10.11_{-7.47}^{+15.60} \text{ yr}^{-1}$, where the errors are quoted at the 90 per cent confidence interval. This rate is also plotted in Fig. 8.

### 7.3.2 Ab initio simulations

Chruslinska et al. (2017) predicted a DNS merger rate density of $48.4 \text{ Gpc}^{-3} \text{ yr}^{-1}$, which translates to a rate of $0.45 \text{ yr}^{-1}$ using a LIGO range distance of 130 Mpc. This is clearly lower than the range of merger rates presented in this work. They also offer a variety of different models, where the most optimistic one derives a rate density of $600_{-300}^{+600} \text{ Gpc}^{-3} \text{ yr}^{-1}$, which corresponds to a merger rate range of $5.52_{-2.76}^{+5.52} \text{ yr}^{-1}$ using the previously introduced range distance.

In comparison, Kruckow et al. (2018) synthesized binary populations at different metallicities ($Z_{MW} = 0.0088$; $Z_{IzW18} = 0.0002$), taking into account their development towards DNS systems. Amongst others, they present three different Galactic merger rate densities relevant in the context of our analysis. The smallest one was derived using a default model, an intermediate one using an optimistic model to enhance the NS-NS merger rate, as well as the maximal merger rate density that can be realistically produced. Scaling all three to

the LIGO range distance of 130 Mpc, we arrive at a prediction of $0.82_{-0.64}^{+2.86} \text{ yr}^{-1}$.

It is remarkable, that for both these ab-initio simulations, only the most optimistic, respectively maximally realistic estimates are compatible with the results from the validation of the results by Pol et al. (2020) and our enhanced simulations using the generic profile, as can be seen in Fig. 8. The results of their reference models are both clearly lower than the estimates resulting from LIGO or derived here.

## 8 SUMMARY AND CONCLUSION

Improving on important assumptions made in previous estimates of the Galactic DNS merger rate and the resulting LIGO detection rate, we obtain results that are in good agreement with previous estimates when considering the whole sample of DNS included in this study. This is in contrast to the individual rate derived for PSR J1906 + 0746, which shows the expected large contribution to the overall rates and which increases using the measured beam shape in comparison to previous assumptions or generic profiles. Applying generic profiles derived from PSR J1906 + 0746's observation to the other DNSs makes, however, little overall difference in the total rate when combining the results. We therefore conclude that the method of estimating the DNS merger and LIGO detection rates via the study of the known Galactic radio pulsar DNS population is less prone to systematic uncertainties than previously thought and argue that the derived estimates should be considered as robust. Consequently, following the work presented here, we predict a detection of three to nine DNS mergers per year by Advanced LIGO (within 90 per cent confidence intervals). However, there are still a number of other systematic uncertainties present in this and other methods, so that actual results from the third operating run O3 of LIGO will provide
significant insight in the correctness of the different estimates and
may allow us to address the various caveats discussed earlier.

In the case of PSR J1906 + 0746, further measurements to
improve the implemented intensity profile can be made until the
pulsar disappears from view, constraining our extrapolations due to
even better knowledge of the beam shape. The current measurements
(Fig. 1) show significant differences in the shape of the main
pulse and the interpulse, imposing important constraints on pulsar
emission models. Hence, future work on these aspects might result
in interesting insights on the nature of pulsar beams and the beaming
process itself.

Moreover, the differences between ab initio simulations and the
simulations based on the population of detected DNS systems could
be resolved by including detected HMXB (high mass X-Ray binaries)
systems or radio-quiet NS-star binaries into the detection-based
simulations, as they depict preliminary stages to DNS systems.

Finally, for a number of purposes, it would be tremendously helpful
to expand the tracking of pulsars showing relativistic spin precession,
allowing us to have detailed studies such as for PSR J1906 + 0746
also for other DNS systems. This would not only help our under-
standing of the pulsar emission processes and the general structure
of pulsar emission beams but also help to decrease the overall
uncertainties of DNS merger rate predictions.

ACKNOWLEDGEMENTS

We thank Nihan Pol for helpful discussions and provision of his
simulation framework. MK and GD are supported by the European
Research Council for the ERC Synergy Grant BlackHoleCam under
contract no. 610058. We used the code to determine the SNR
degradation factor of PSR J0508 + 3801 developed and written
by Bagchi, Lorimer & Wolfe (2013).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request
to the corresponding author.

REFERENCES

Abbott B. P. et al., 2017a, Phys. Rev. Lett., 119, 161101

Abbott B. P. et al., 2017b, ApJ, 848, L12

Abbott B. P., Abbott R., Abbott T. D., Abernathy M. R., Acernese F., Ackley
K., 2018, Living Rev. Relativ., 21, 3

Abbott B. P. et al., 2020, ApJ, 892, L3

Antoniadis J., Bassa C. G., Wex N., Kramer M., Napiwotzki R., 2011,
MNRAS, 412, 580

Bagchi M., Lorimer D. R., Wolfe S., 2013, MNRAS, 432, 1303

Barish B. C., 1999, in American Physical Society (APS) Meeting of the
Division of Particles and Fields (DPF 99). preprint(arXiv:gr-qc/9905026)

Bates S. D., Lorimer D. R., Rane A., Swiggum J., 2014, MNRAS, 439, 2893

Belczynski K., Kalogera V., Rasio F. A., Taam R. E., Zezas A., Bulik T.,
Maccarone T. J., Ivanova N., 2008, ApJS, 174, 223

Chen H.-Y., Holz D. E., Miller J., Evans M., Vitale S., Creighton J., 2017,
Class. Quantum Gravity, 38, 055010

Chruslinska M., Belczynski K., Klencki J., Benacquista M., 2017, MNRAS,
474, 2937

Desvignes G. et al., 2019, Science, 365, 1013

Ferdman R. D. et al., 2014, MNRAS, 443, 2183

Finn L. S., Chernoff D. F., 1993, Phys. Rev. D, 47, 2198

Gil J., Gronkowski P., Rudnicki W., 1984, A&A, 132, 312

Gil J. A., Kijak J., Seiradakis J. H., 1993, A&A, 272, 268

Gould D. M., Lyne A. G., 1998, MNRAS, 301, 235

Gralla S. E., Lupsasca A., Philippov A., 2017, ApJ, 851, 137

Haniewicz H. T., Ferdman R. D., Freire P. C. C., Champion D. J., Bunting K.
A., Lorimer D. R., McLaughlin M. A., 2021, MNRAS, 500, 4620

Hulse R. A., Taylor J. H., 1975, ApJ, 195, L51

Johnston S., Kramer M., 2019, MNRAS, 490, 4565

Kim C., Kalogera V., Lorimer D. R., 2003, ApJ, 584, 985

Kopparapu R. K., Hanna C., Kalogera V., O'Shaughnessy R., González G.,
Brady P. R., Fairhurst S., 2008, ApJ, 675, 1459

Kramer M., 1998, ApJ, 509, 856

Kramer M., Wex N., 2009, Class. Quantum Gravity, 26, 073001

Kramer M., Wielebinski R., Jessner A., Gil J. A., Seiradakis J. H., 1994,
A&AS, 107, 515

Kramer M., Xilouris K. M., Lorimer D. R., Doroshenko O., Jessner A.,
Wielebinski R., Wolszczan A., Camilo F., 1998, ApJ, 501, 270

Kramer M. et al., 2006, Science, 314, 97

Kramer M. et al., 2021, MNRAS, 504, 2094

Kruckow M. U., Tauris T. M., Langer N., Kramer M., Izzard R. G., 2018,
MNRAS, 481, 1908

Lorimer D. R., Kramer M., 2005, Handbook of Pulsar Astronomy. Cambridge
Univ. Press, Cambridge, England

Lorimer D. R. et al., 2006, ApJ, 640, 428

Lynch R. S. et al., 2018, ApJ, 859, 93

Lyne A. G. et al., 2004, Science, 303, 1153

Manchester R. N. et al., 2010, ApJ, 710, 1694

McLaughlin M. A. et al., 2004, ApJ, 613, L57

O'Shaughnessy R., Kim C., 2010, ApJ, 715, 230

Perera B. B. P. et al., 2010, ApJ, 721, 1193

Phinney E. S., 1991, ApJ, 380, L17

Pol N., McLaughlin M., Lorimer D. R., 2019, ApJ, 870, 71

Pol N., McLaughlin M., Lorimer D. R., 2020, Res. Notes AAS, 4, 22

Radhakrishnan V., Cooke D. J., 1969, ApJ, 3, 225

Rankin J. M., 1993, ApJ, 405, 285

Stairs I. H., Thorsett S. E., Arzoumanian Z., 2004, Phys. Rev. Lett., 93,
141101

Tauris T. M., Sennels T., 2000, A&A, 355, 236

Tauris T. M. et al., 2017, ApJ, 846, 170

Taylor J. H., Weisberg J. M., 1982, ApJ, 253, 908

Timokhin A. N., Arons J., 2013, MNRAS, 429, 20

van Leeuwen J. et al., 2015, ApJ, 798, 118

Venkatraman Krishnan V., Bailes M., van Straten W., Keane E. F., Kramer
M., Bhat N. D. R., Flynn C., Osłowski S., 2019, ApJ, 873, L15

Vigna-Gómez A., Schröder S. L., Ramirez-Ruiz E., Aguilera-Dena D. R.,
Batta A., Langer N., Wilcox R., 2021, preprint (arXiv:2106.12381)

Weisberg J. M., Taylor J. H., 2002, ApJ, 576, 942

Weltevrede P., Johnston S., 2008, MNRAS, 387, 1755

This paper has been typeset from a *TEX/LATex* file prepared by the author.