# Multi-messenger observations of double neutron stars in Galactic disk
with gravitational and radio waves

Wen-Fan Feng  MOE Key Laboratory of Fundamental Physical Quantities
Measurements, Hubei Key Laboratory of Gravitation and Quantum Physics, PGMF,
Department of Astronomy and School of Physics, Huazhong University of Science
and Technology, Wuhan 430074, China  Jie-Wen Chen  MOE Key Laboratory of
Fundamental Physical Quantities Measurements, Hubei Key Laboratory of
Gravitation and Quantum Physics, PGMF, Department of Astronomy and School of
Physics, Huazhong University of Science and Technology, Wuhan 430074, China
Yan Wang<EMAIL_ADDRESS>MOE Key Laboratory of Fundamental Physical
Quantities Measurements, Hubei Key Laboratory of Gravitation and Quantum
Physics, PGMF, Department of Astronomy and School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, China  Soumya D. Mohanty
Department of Physics and Astronomy, University of Texas Rio Grande Valley,
Brownsville, Texas 78520, USA  Department of Physics, IIT Hyderabad, Kandai,
Telangana-502284, India  Yong Shao  School of Astronomy and Space Science,
Nanjing University, Nanjing, 210023, China  Key Laboratory of Modern Astronomy
and Astrophysics, Nanjing University, Ministry of Education, Nanjing, 210023,
China

###### Abstract

We evaluate the prospects for radio follow-up of the double neutron stars
(DNSs) in the Galactic disk that could be detected through future space-borne
gravitational wave (GW) detectors. We first simulate the DNS population in the
Galactic disk that is accessible to space-borne GW detectors according to the
merger rate from recent LIGO results. Using the inspiraling waveform for the
eccentric binary, the average number of the DNSs detectable by TianQin (TQ),
LISA, and TQ+LISA are 217, 368, and 429, respectively. For the joint GW
detection of TQ+LISA, the forecasted parameter estimation accuracies, based on
the Fisher information matrix, for the detectable sources can reach the levels
of $\Delta P_{\mathrm{b}}/P_{\mathrm{b}}\lesssim 10^{-6}$,
$\Delta\Omega\lesssim 100~{}{\mathrm{deg}}^{2}$, $\Delta e/e\lesssim 0.3$, and
$\Delta\dot{P}_{\mathrm{b}}/\dot{P}_{\mathrm{b}}\lesssim 0.02$. These
estimation accuracies are fitted in the form of power-law function of signal-
to-noise ratio. Next, we simulate the radio pulse emission from the possible
pulsars in these DNSs according to pulsar beam geometry and the empirical
distributions of spin period and luminosity. For the DNSs detectable by
TQ+LISA, the average number of DNSs detectable by the follow-up pulsar
searches using the Parkes, FAST, SKA1, and SKA are 8, 10, 43, and 87,
respectively. Depending on the radio telescope, the average distances of these
GW-detectable pulsar binaries vary from 1 to 7 kpc. Considering the dominant
radiometer noise and phase jitter noise, the timing accuracy of these GW-
detectable pulsars can be as low as 70 ${\rm ns}$ while the most probable
value is about 100 $\mu{\rm s}$.

##  I Introduction

Binary pulsars are emitters of both radio and gravitational waves (GWs) and,
as has been the case since PSR B1913+16 [1], short-period binary pulsars
provide excellent laboratories for constraining the theories of gravity [2].
An excellent example is the Double Pulsar PSR J0737-3039A/B [3]. Due to
various technical limitations described below, current radio surveys for such
systems can only find orbital periods longer than about an hour. However,
upcoming space-borne GW detectors, i.e., LISA [4, 5], TianQin (TQ) [6], and
Taiji [7], will be able to detect systems with shorter orbital periods. It is
therefore interesting to consider the prospects of radio follow-up
observations of short-period double neutron star (DNS) systems detected using
GWs.

Through radio surveys, a total of 20 DNS systems have been discovered up to
now [8, 9], in which the shortest orbital period is found in PSR J1946+2052
with a period of $1.88~{}{\rm hr}$ [10]. To detect pulsars in binary systems,
pulsar searching algorithms, in principle, need to consider the Doppler
smearing effect in arrival times of radio pulses due to orbital motion in
addition to the routine de-dispersion, folding, radio frequency interference
removal, etc. When integration time for pulsar searching, typically $10-30$
min, becomes a considerable fraction of the orbital period, the Doppler
smearing will be worsened and must be corrected in a searching algorithm.
Several methods, such as acceleration search [11] and jerk search [12], have
been proposed to compensate for the loss of search sensitivity due to the
binary motion, but they are only applicable for integration times that are
shorter than tens of percent of the orbital period. The other methods, such as
dynamic power spectrum [13] and phase modulation search [14], have also been
developed to extend the applicability to the case that the integration time
can be approximately equal to or even longer than the orbital period, however,
with a cost of the high demand of computational resources.

As a consequence of the above limitations, all 359 binary pulsars found so far
[15, 16] have orbital periods longer than $1.2~{}{\rm hr}$. Tight Galactic DNS
systems with orbital periods in the $3-60$ min range will emit GWs with
frequencies in the mHz band accessible to space-borne GW detectors. From the
discoveries by LIGO [17] and Virgo [18] of GWs from two DNS merger events,
GW170817 [19] and GW190425 [20], the inferred merger rate density is
$10-1700~{}\mathrm{Gpc^{-3}~{}yr^{-1}}$ for a 90% credible interval [21].
Taking $920~{}\mathrm{Gpc^{-3}~{}yr^{-1}}$ as the fiducial point estimate, as
in [22], we get a corresponding DNS merger rate of $210~{}\rm{Myr}^{-1}$ for a
typical Milky-Way type galaxy. From this inferred merger rate, one expects
$\sim 2100$ DNS systems with orbital periods range less than $1~{}{\rm hr}$ in
our Galaxy [23]. Therefore, despite the challenges of detecting DNSs directly
by radio surveys, the space-borne GW detectors can potentially detect DNS
systems with orbital periods less than $1~{}{\rm hr}$.

Several recent works have investigated the number of the DNSs in our Galaxy
(and possible nearby galaxies) that can be detected by LISA. Seto estimated
the frequency distribution of the DNSs in the Local Group of the galaxies
based on the higher merger rate $1540~{}\mathrm{{Gpc}^{-3}{yr}^{-1}}$ derived
from the combination of GW170817 and the LIGO-Virgo O1 data release [19]. This
work shows that LISA can detect about 5 DNSs in the Large Magellanic Cloud and
Andromeda with signal-to-noise ratios (SNRs) greater than 10 for a 10-year
observation [24]. Lau et al. simulated a set of DNSs in the LISA band using
the population synthesis code COMPAS [25] and found that for a 4-year
observation and the merger rate of $33~{}\mathrm{Myr}^{-1}$, LISA can detect
about 35 DNSs with SNRs greater than 8 (33 in our Galaxy, 2 in the Large
Magellanic Cloud and Small Magellanic Cloud) [26]. Andrews et al. used the
merger rate of $210~{}\mathrm{Myr}^{-1}$ to obtain the number of the Galactic
DNSs in the mHz band. They found that, for an observation of 4 (8) years and
the SNR threshold of 7, LISA can detect 240 (330) DNSs in our Galaxy, 1 (4) in
Andromeda [23].

The above studies have focused on the detection of DNSs using GWs alone. To
mitigate the problems, outlined earlier, associated with detecting short-
period DNSs using radio observations, Kyutoku et al. discussed a multi-
messenger strategy to detect Galactic pulsar-NS binaries with orbital periods
less than 10 min [27]. They showed that LISA can detect these DNSs with SNRs
greater than 100 and accurately determine their sky positions and orbital
parameters. In the scenario of synchronous observation of LISA and Square
Kilometre Array (SKA) [28] starting in the mid-2030s, the information provided
by LISA will enable a reduction of 2 to 6 orders of magnitude in the number of
telescopes pointing comparing with the all-sky search and a further reduction
of 9 orders of magnitude in the number of trials for orbital modulation in
order to correct the Doppler smearing. On the side of astrophysics, Thrane et
al. explored the prospect of using the Lense-Thirring precession effect to
constrain the equation of state of NSs based on LISA+SKA multi-messenger
observations [29]. They showed that using 10-year observation of SKA the
measurement accuracy of the mass-radius relation of NSs can reach $0.2\%$,
which is an order of magnitude better than the X-ray pulsation measurements
[30] and the GW observations [31].

Multi-messenger observations, including GW and radio (and possibly X-ray), can
be highly complementary and provide a more complete picture of the DNS
population and the structure of the NS. Regarding the DNS population, ground-
based GW detectors have currently detected coalescing DNSs, while radio
telescopes have detected DNSs with orbital periods of more than an hour. There
is a gap in the detection of the DNSs with periods of a few minutes to an
hour. The orbital periods of the DNSs that can be detected by space-borne GW
detectors (e.g., LISA and TQ) can fill this gap. Once the pulsar binaries are
detected from the GW-detectable DNSs, this will directly map the short-period
pulsar binaries, which will help to study the DNS formation and population.
Furthermore, with the detection of such a system we can study one of the most
important problems of the NS, namely the equation of state. Although both GW
detection and radio pulse detection can constrain the equation of state in
separate ways and have their own detection ranges and accuracies, the
information obtained from multi-messengers is expected to be able to constrain
it more tightly [32, 29].

In this work, we study the detectability of Galactic disk DNSs in mHz
frequency band using a network of space-borne GW detectors and the prospects
for radio follow-up of such detected systems. In addition to the number of
detectable sources, this work forecasts the accuracies of the estimated
parameters of the DNSs by taking into account the antenna responses of space-
borne GW detectors to the inspiraling waveforms of these systems. The follow-
up search for the possible radio pulsar components in the DNS systems is
studied based on both the intrinsic properties of the pulsars and the
specifications of the large radio telescopes including the Parkes radio
telescope [33], the Five-hundred-meter Aperture Spherical radio Telescope
(FAST) [34], SKA phase 1 (SKA1) [35] and SKA [36]. We estimate the timing
precision of the detectable pulsars by the radio telescopes, which are crucial
for the further investigations of the ultra-compact DNSs as laboratories for
astrophysics and fundamental physics.

The rest of this paper is organized as follows. Sec. II presents the
simulation of the DNS population in the Galactic disk. Sec. III.1 introduces
the GW signal model for the DNS in an eccentric orbit, the noise
characteristics for LISA and TQ, and the Fisher information approach to obtain
the parameter estimation accuracy. The numbers of detectable DNSs and
parameter estimation accuracies expected for TQ, LISA, and TQ+LISA network are
discussed in Sec. III.2. Sec. IV investigates the follow-up searches of the
pulsars based on the Parkes, FAST, SKA1, and SKA. The pulsar beam geometry and
the empirical distributions of the spin periods and luminosities are given in
Sec. IV.1. The synergy for the multi-messenger observations with the GW
detectors and the radio telescopes is given in Sec. IV.2. The paper concludes
in Sec. V.

##  II Galactic DNS population

In this work, we follow the method in Sec. 2.4 of [23] to synthesize the
Galactic population of tight DNSs. We assume the component masses of all DNSs
to be $m_{1}=m_{2}=1.4~{}M_{\odot}$. Hence, the chirp masses of the binaries
are $M_{\mathrm{c}}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}=1.22~{}M_{\odot}$.

Given the DNS merger rate and the Galactic disk model, we obtain the number
and locations of the DNSs at a specific evolution time. Here, we focus on
those DNSs that will merge within $10~{}{\rm Myr}$. Note that there is no need
to extend to longer merger times because the sensitivity of LISA or TQ is
limited to DNSs that are expected to merge within the next $\sim
10~{}\mathrm{Myr}$ (cf. Fig. 2). Using the merger rate
$R_{\mathrm{MW}}=210~{}\mathrm{Myr}^{-1}$ [23, 22], the expected number of
DNSs in our simulation is $N_{\mathrm{DNS}}=R_{\mathrm{MW}}\times
10~{}\mathrm{Myr}=2100$. We assume that these DNSs are randomly distributed in
the Galactic disk with a number density profile [37]

$f_{\mathrm{dis}}(r,z)\propto
e^{-r/L}\mathrm{sech\,}(z/\beta)^{2}~{}\mathrm{pc}^{-3}\,,$ (1)

where the disk scale length $L=2.5~{}\mathrm{kpc}$, the scale height
$\beta=0.2~{}\mathrm{kpc}$. The distribution for the azimuthal angle follows
$\mathcal{U}[0,2\pi]$. Here, $\mathcal{U}(a,b)$ denotes the continuous uniform
distribution within the interval $[a,b]$.

The differential equations of the binary orbit due to GW emission allow us to
obtain the evolution of these parameters with time. For the binary in an
eccentric orbit, the time evolution of its orbital semi-major axis $a$ and
eccentricity $e$ due to GW emission can be expressed as [38]:

$\displaystyle\left\langle\frac{{\rm d}a}{{\rm d}t}\right\rangle$
$\displaystyle=$ $\displaystyle-\frac{64}{5}\frac{G^{3}\mu
M^{2}}{c^{5}a^{3}}F(e)\,,$ (2) $\displaystyle\left\langle\frac{{\rm d}e}{{\rm
d}t}\right\rangle$ $\displaystyle=$
$\displaystyle-\frac{304}{15}\frac{G^{3}\mu
M^{2}e}{c^{5}a^{4}\left(1-e^{2}\right)^{5/2}}\left(1+\frac{121}{304}e^{2}\right)\,,$
(3)

where $\left\langle\cdot\right\rangle$ denotes the average over one orbital
period. $M$ and $\mu$ are the total mass and reduced mass of the system.

$F(e)={\left(1+\frac{73}{24}{{e}^{2}}+\frac{37}{96}{{e}^{4}}\right)}/{{{\left(1-{{e}^{2}}\right)}^{7/2}}}$
(4)

is the enhancement factor. Combining Eq. (2) and Eq. (3), one gets the
evolution of the semi-major as a function of the orbit eccentricity:

$a(e)=c_{0}\frac{e^{12/19}}{1-e^{2}}\left(1+\frac{121}{304}e^{2}\right)^{870/2299}\,.$
(5)

Let $g(e)=a(e)/c_{0}$, then $c_{0}$ satisfies $a(e_{0})=c_{0}g(e_{0})$.

One can solve the evolution of $e$ over time by inserting Eq. (5) in Eq. (3)
and numerically integrating the ordinary differential equation. Since Eq. (5)
contains $c_{0}$, the evolution of $e$ over time depends on the initial (or
reference) values $a_{0}$ and $e_{0}$, as well as the integration time. We
take the orbital parameters of PSR B1913+16 close to the merger as the
reference values: $P_{\mathrm{b0}}=0.2~{}\mathrm{s}$, $e_{0}=5\times 10^{-6}$
[23]. We draw the merger time of the DNSs in our simulations from a uniform
distribution over the past $10~{}\mathrm{Myr}$. Combining Kepler’s third law
with Eq. (5), one can obtain the relation between $e$ and orbital period
$P_{\mathrm{b}}$, then the evolution of $P_{\mathrm{b}}$ over time can also be
obtained.

The DNS population is generated by assigning these orbital parameters to the
aforementioned simulated sources.

##  III Detecting GW of DNS

###  III.1 DNS gravitational waveform, detector noise, and parameter
estimation accuracy

Under the transverse-traceless gauge and quadrupole moment approximation, we
can analytically deduce the gravitational waveforms for two polarizations of
an eccentric binary [39]

$\displaystyle h_{+}$ $\displaystyle=$
$\displaystyle-\frac{h_{0}}{2}\sum_{n=1}^{\infty}\left\\{\sin^{2}\iota
J_{n}(ne)\cos
nM_{\mathrm{A}}+\left(1+\cos^{2}\iota\right)\left[\frac{C_{n}-S_{n}}{2}\cos\left(nM_{\mathrm{A}}+2\Phi_{\mathrm{P}}\right)+\frac{S_{n}+C_{n}}{2}\cos\left(nM_{\mathrm{A}}-2\Phi_{\mathrm{P}}\right)\right]\right\\}\,,$
$\displaystyle h_{\times}$ $\displaystyle=$ $\displaystyle-
h_{0}\cos\iota\sum_{n=1}^{\infty}\left[{\frac{S_{n}-C_{n}}{2}\sin\left(nM_{\mathrm{A}}+2\Phi_{\mathrm{P}}\right)}+\frac{S_{n}+C_{n}}{2}\sin\left(nM_{\mathrm{A}}-2\Phi_{\mathrm{P}}\right)\right]\,,$
(6)

with

$\displaystyle h_{0}$ $\displaystyle=$
$\displaystyle\frac{4\left(GM_{\mathrm{c}}\right)^{5/3}\omega_{\mathrm{b}}^{2/3}}{c^{4}D}\,,$
(7) $\displaystyle{{S}_{n}}$ $\displaystyle=$
$\displaystyle-\frac{2{{\left(1-{{e}^{2}}\right)}^{{1}/{2}}}}{ne}\frac{{\rm
d}{{J}_{n}}(ne)}{{\rm
d}e}+\frac{2n{{\left(1-{{e}^{2}}\right)}^{{3}/{2}}}{{J}_{n}}(ne)}{{{e}^{2}}}\,,$
(8) $\displaystyle{{C}_{n}}$ $\displaystyle=$
$\displaystyle-\frac{\left(2-{{e}^{2}}\right){{J}_{n}}(ne)}{{{e}^{2}}}+\frac{2\left(1-{{e}^{2}}\right)}{e}\frac{{\rm
d}{{J}_{n}}(ne)}{{\rm d}e}\,.$ (9)

Here ${{\omega}_{\mathrm{b}}}={2\pi}/{{{P}_{\mathrm{b}}}}$ is the orbital
angular frequency, $D$ is the distance to the source, $\iota$ is the
inclination angle of the orbit, ${{J}_{n}}(ne)$ is the Bessel function of the
first kind $(n=1,2,3,...)$ with $n$ as the harmonic number,
${{M}_{\mathrm{A}}}={2\pi(t-{{t}_{0}})}/P_{\mathrm{b}}$ is the mean anomaly,
${{t}_{0}}$ is the time of periastron passage, and ${{\Phi}_{\mathrm{P}}}$ is
the periastron phase. The aforementioned waveforms can be transformed as
follows to include the effects of the orbital frequency derivative and Doppler
phase modulation in the GW phase:

$\displaystyle nM_{\mathrm{A}}\longrightarrow~{}$ $\displaystyle
nM_{\mathrm{A}}-n\frac{\pi\dot{P_{\mathrm{b}}}}{P_{\mathrm{b}}^{2}}(t-t_{0})^{2}+$
(10) $\displaystyle
2\pi\frac{n}{P_{\mathrm{b}}}\frac{R}{c}\sin\theta_{\mathrm{S}}\cos\left(2\pi
f_{m}t-\phi_{\mathrm{S}}\right)\,,$

where $R$ is the distance between the Sun and Earth and $f_{m}=1/(1~{}{\rm
yr})$ is the modulation frequency due to the annual revolution of the
detector’s constellation (for either LISA or TQ) around the Sun.
($\theta_{\mathrm{S}},\phi_{\mathrm{S}}$) are the colatitude and longitude of
the source in the heliocentric ecliptic coordinate system.

At the low frequency (i.e., long wavelength) limit, the strain signal induced
by GWs can be expressed as

$h_{\mathrm{I,II}}(t)={{F}^{+}_{\mathrm{I,II}}}(t)\cdot{{h}_{+}}(t)+{{F}^{\times}_{\mathrm{I,II}}}(t)\cdot{{h}_{\times}}(t)\,,$
(11)

where ${{F}^{+,\times}_{\mathrm{I,II}}}$ is the antenna response functions of
the detector for the two polarizations ($+,\times$) and two independent
Michelson-type interferometers (I, II) in the heliocentric ecliptic
coordinates. See [40] and [41] for their explicit expressions for LISA and TQ,
respectively.

The non-sky-averaged sensitivity curves for LISA [42] and TQ [43] can be
obtained by

$\displaystyle{{S}_{\rm L}}(f)$ $\displaystyle=$
$\displaystyle\frac{1}{{{L}_{\rm
L}^{2}}}\left[{{P}_{\mathrm{OMS}}}(f)+\frac{4{{P}_{\mathrm{acc}}}(f)}{{{\left(2\pi
f\right)}^{4}}}\right]\left[1+\frac{6}{10}{{\left(\frac{f}{{{f}_{*}^{\rm
L}}}\right)}^{2}}\right]$ (12) $\displaystyle+{{S}_{c}}(f)\,,$
$\displaystyle{{S}_{\rm T}}(f)$ $\displaystyle=$
$\displaystyle\frac{1}{{{L}_{\rm T}^{2}}}\left[\frac{4{{S}_{a}}}{{{\left(2\pi
f\right)}^{4}}}\left(1+\frac{{{10}^{-4}}~{}\mathrm{Hz}}{f}\right)+{{S}_{x}}\right]$
(13) $\displaystyle\times\left[1+0.6{{\left(\frac{f}{{{f}_{*}^{\rm
T}}}\right)}^{2}}\right]\,.$

The arm length for LISA and TQ are ${L}_{\rm L}=2.5\times 10^{9}~{}{\rm m}$
and ${L}_{\rm T}=1.7\times 10^{8}~{}{\rm m}$, respectively. The associated
transfer frequency for them are ${f}_{*}^{\rm L}=19.09~{}{\rm{mHz}}$ and
${f}_{*}^{\rm T}=0.28~{}{\rm{Hz}}$. The detailed budgets for single-link
optical metrology noise ${P}_{\mathrm{OMS}}$ and single test mass acceleration
noise ${P}_{\mathrm{acc}}$ for LISA and corresponding noise ${S}_{x}$ and
${S}_{a}$ for TQ can be found in [42] and [41, 43], respectively. $S_{c}$ is
the PSD of the stochastic GW signal associated with the unresolvable Galactic
binary population as a whole (dominated by WD binaries). Note that we do not
include $S_{c}$ for TQ since the instrumental noise for TQ is higher than it
in the concerned frequency range.

Assuming the noise of LISA and TQ are uncorrelated, the joint sensitivity
curve of LISA and TQ can be obtained from the following expression [44]

${S_{\rm{T+L}}}(f)\equiv\frac{{S_{\rm T}}(f){S_{\rm L}}(f)}{{S_{\rm
T}}(f)+{S_{\rm L}}(f)}\,.$ (14)

Given GW signal is sufficiently strong, the estimation accuracies of unknown
parameters in the GW signal can be approximately calculated by using the
Fisher information matrix (FIM),

${{\Gamma}^{ij}}\equiv\left(\frac{\partial
h}{\partial{{\lambda}_{i}}},\frac{\partial
h}{\partial{{\lambda}_{j}}}\right)\,,$ (15)

where (,) denotes the noise-weighted inner product that contains both the GW
signal and the detector noise [45]. Similarly, the optimal SNR for signal
detection is defined as ${\rm SNR}\equiv(h,h)^{1/2}$. The root-mean-square
(rms) error of parameter ${\lambda_{i}}$ is estimated as

$\Delta{\lambda_{i}}=\sqrt{{\Sigma_{ii}}}\,,$ (16)

where the variance-covariance matrix $\Sigma$ is the inverse of the FIM, i.e.,
$\Sigma={\Gamma^{-1}}$. The correlation coefficient between $\lambda_{i}$ and
$\lambda_{j}$ can be calculated by

$c_{ij}=\Sigma_{ij}/\sqrt{\Sigma_{ii}\Sigma_{jj}}\,.$ (17)

For a GW source with sky location $(\theta_{\mathrm{S}},\phi_{\mathrm{S}})$,
its sky localization error is defined as [46]

$\Delta\Omega=2\pi\sqrt{\Sigma_{\cos\theta_{\mathrm{S}}\cos\theta_{\mathrm{S}}}\Sigma_{\phi_{\mathrm{S}}\phi_{\mathrm{S}}}-\Sigma_{\phi_{\mathrm{S}}\cos\theta_{\mathrm{S}}}^{2}}\,.$
(18)

###  III.2 Results on detectable DNSs

We assume the observation time $T_{\mathrm{obs}}=4~{}\mathrm{yr}$. The
inclination angle $\cos\iota\sim\mathcal{U}(-1,1)$, the polarization angle
${\psi}\sim\mathcal{U}(0,\pi)$, and the periastron phase
$\Phi_{\mathrm{P}}\sim\mathcal{U}(0,2\pi)$. We draw $N_{\mathrm{DNS}}$ random
samples from the above distributions and assign them to the simulated DNSs in
the Galactic disk. The assumed maximum unknown parameter space for FIM is
$\bm{\lambda}=\left(\ln
h_{0},\cos\iota,\psi,\cos\theta_{\mathrm{S}},\phi_{\mathrm{S}},\ln
P_{\mathrm{b}},\Phi_{\mathrm{P}},\ln e,\ln\dot{P}_{\mathrm{b}}\right)$. We set
the SNR threshold to be 7, above which we regard the source as detectable.
According to the chirp frequency limit [37], i.e.,
$\dot{f}T_{\mathrm{obs}}\geqslant 1/T_{\mathrm{obs}}$, we get the frequency
range in which the chirp mass can be measured through GW signals

$\frac{f_{\mathrm{chirp}}}{\mathrm{~{}mHz}}\geqslant
1.75\left(\frac{{M}_{\mathrm{c}}}{1.22~{}M_{\odot}}\right)^{-5/11}\left(\frac{T_{\mathrm{obs}}}{4~{}\mathrm{yr}}\right)^{-6/11}F(e)^{-3/11}\,.$
(19)

Note that chirp mass and eccentricity are degenerate in the GWs of DNS when
the eccentricity is non-negligible. In this case, Eq. (19) is a necessary but
not sufficient condition for the measurement of chirp mass. However, as will
become clear later, in the current work, Eq. (19) only serves as a criterion
for selecting sources for subsequent evaluation of the estimation accuracy of
chirp mass.

For the DNSs with SNRs higher than the threshold but the GW frequency lower
than $f_{\mathrm{chirp}}$, their parameter space is $\bm{\lambda}$ with
$\ln\dot{P}_{\mathrm{b}}$ eliminated.

For a quasi-monochromatic source with frequency $f_{0}$ and observation time
$T_{\mathrm{obs}}$, the FIM can be simplified as [47, 48]

$\displaystyle{\Gamma}^{ij}$ $\displaystyle=$
$\displaystyle\frac{2}{S_{n}\left(f_{0}\right)}\sum_{\alpha=\mathrm{I,II}}\int_{-\infty}^{\infty}\partial_{i}\widetilde{h}_{\alpha}^{*}(f)\partial_{j}\widetilde{h}_{\alpha}(f)df\,,$
(20) $\displaystyle\simeq$
$\displaystyle\frac{2}{S_{n}\left(f_{0}\right)}\sum_{\alpha=\mathrm{I,II}}\int_{0}^{T_{\mathrm{obs}}}\partial_{i}h_{\alpha}(t)\partial_{j}h_{\alpha}(t)dt\,.$

Here $\widetilde{h}_{\alpha}^{*}(f)$ is the complex conjugate of the Fourier
transform of $h_{\alpha}(t)$. Note that $\mathrm{SNR}^{2}=\Gamma^{11}$ for the
parameter space we choose. We can replace the partial derivatives in FIM with
the central finite difference

$\frac{\partial
h}{\partial\lambda_{i}}\approx\frac{h\left(t,\lambda_{i}+\delta\lambda_{i}\right)-h\left(t,\lambda_{i}-\delta\lambda_{i}\right)}{2~{}\delta\lambda_{i}}\,,$
(21)

where $\delta\lambda_{i}/\lambda_{i}\lesssim 10^{-8}$, and their values can be
adjusted to obtain stable results [48]. In addition, the time integral in Eq.
(20) is transformed into a summation in which the sampling frequency is
selected to be $f_{\mathrm{s}}=4f_{0}$.

We perform 100 independent Monte Carlo simulations of the DNSs in the Galactic
disk and use Gaussian distributions to fit the numbers of the DNSs detectable
by TQ, LISA, and TQ+LISA network (ideally this indicates that TQ and LISA will
be operated synchronously and the data are analyzed jointly), respectively
(Fig. 1). The average numbers of detectable DNSs are 217, 368, and 429 for
those three cases.

Figure 1: The histograms for the number of detectable DNSs by TQ, LISA, and
TQ+LISA network from 100 independent Monte Carlo simulations. Black curves are
the best-fit Gaussian distributions.

####  III.2.1 Characteristic strain

In Fig. 2, we plot the characteristic strains for signal ($h_{\mathrm{c},n}$)
and noise ($h_{{\mathrm{noise},n}}$) which are defined by rewriting the
expression of SNR as follows,

$\displaystyle\mathrm{SNR}^{2}$ $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\mathrm{SNR}^{2}_{n}=\sum_{n=1}^{\infty}\frac{h^{2}_{\mathrm{c},n}}{h^{2}_{{\mathrm{noise},n}}}\,,$
(22) $\displaystyle=$
$\displaystyle\sum_{n=1}^{\infty}\frac{f_{n}\int_{0}^{T_{\mathrm{obs}}}2\left(h_{\mathrm{I},n}^{2}+h_{\mathrm{II},n}^{2}\right)dt}{f_{n}S(f_{n})}\,,$

where the subscript $n$ represents the $n$-th harmonics of the GW in Eq.
(III.1). Each source includes the influence of the antenna response functions.
Similarly, the joint characteristic strain of TQ and LISA can also be derived
from the combined SNR with the joint sensitivity in Eq. (14).

Given the orbital parameters of the reference source PSR B1913+16 and merger
time distribution $\mathcal{U}[0,10]~{}\mathrm{Myr}$), the maximum
eccentricity ($\approx 0.18$) in our simulated DNSs is relatively low. Thus,
the GW signal from the 2nd harmonic is dominant [39], so we only draw the
characteristic strain from the 2nd harmonic in Fig. 2. However, in our
calculation of the SNR and FIM, the first four harmonics ($n=1,2,3,4$) are
taken into account, while the cross terms between different harmonics are
neglected [49].

Since the benchmark used in the calculation of the orbital parameters of the
DNS comes from binary pulsar PSR B1913+16 at merger
($P_{\mathrm{b0}}=0.2~{}\mathrm{s},\,e_{0}=5\times 10^{-6}$), we add the
characteristic strain from the 2nd harmonic of this binary in Fig. 2 as a
reference [23].

Figure 2: The characteristic strain of the simulated DNSs from one of the
Monte Carlo simulations. The bottom abscissa is the GW frequency of the 2nd
harmonics $f_{n=2}$, and the top abscissa is the evolution of the orbital
eccentricity of the DNSs and the corresponding time before its merger based on
the parameters of the benchmark source PSR B1913+16. The red solid, blue solid
and green dashed curves represent the characteristic sensitivity of TQ, LISA,
and TQ+LISA, respectively. The red, blue, and green dots represent the
detectable DNSs by TQ, LISA, and TQ+LISA. The blue vertical dotted line is the
chirp frequency limit for a circular orbit, above which the evolution of the
frequency is detectable.

####  III.2.2 Parameter estimation accuracy

As shown in Fig. 3, the rms errors of all parameters in the FIM decrease, in
general, with the increase of SNR. Among them, the relative estimation
accuracy of the orbital period (or frequency) is the highest, $\Delta
P_{\mathrm{b}}/P_{\mathrm{b}}\lesssim 10^{-6}$. It is clear that there are
strong correlations between the amplitude and the inclination angle as well as
the polarization angle and the periastron phase. For a source with SNR around
100, the sky localization $\Delta\Omega\lesssim 1~{}{\mathrm{deg}}^{2}$, the
relative error of orbital eccentricity $\Delta e/e\lesssim 0.1$, the relative
error of the orbital period derivative
$\Delta\dot{P}_{\mathrm{b}}/\dot{P}_{\mathrm{b}}\lesssim 0.004$. As
represented by the orange dots, the orbital period derivative
$\dot{P}_{\mathrm{b}}$ can be measured in about $17\%$ of all sources. As
shown below, the measurement of $\dot{P}_{\mathrm{b}}$ is critical for the
subsequent determination of the chirp mass and source distance.

Figure 3: The rms errors of the parameters for the DNSs from one of the Monte-
Carlo simulations. All sources have an SNR of the TQ+LISA network larger than
the threshold (vertical dashed lines near the left edges). The green dots
correspond to relatively low-frequency systems with $\dot{P}_{\mathrm{b}}$ set
to zero and excluded from the parameter estimation space. Overall, the
parameter estimation accuracies improve as the SNR increases. There is a
strong anti-correlation between the amplitude and the inclination as well as
the polarization angle and the periastron phase, hence their distributions
show larger scattering.

Following the analysis of the correlation between parameters in [48], the
strain waveforms in Eq. (11) with inclination $0^{\circ}\leq\iota\lesssim
45^{\circ}$ have the same waveform structure, so by changing the amplitude or
the inclination, we can obtain a very similar signal, i.e. there is a strong
anti-correlation between the $h_{0}$ and $\cos{\iota}$. However, for
inclination $45^{\circ}\lesssim\iota\leq 90^{\circ}$, the magnitude and
structure of the signal change for different inclinations, so they are no
longer correlated. Secondly, the polarization angle $\psi$ is measured in the
plane perpendicular to the GW propagation direction, while the periastron
phase angle $\Phi_{\rm p}$ is measured in the orbital plane. When the sources
are getting close to being face-on to the detector, $\psi$ and $\Phi_{\rm p}$
will appear to have correlations due to partial overlap between those two
planes. This explains the poor parameter accuracy for some detectable sources
(see Fig. 3).

Figure 4: Same as Fig. 3, but remove the inclination $\iota$ and polarization
angle $\psi$ in the two pairs of strongly correlated parameters.

In the case that the electromagnetic observations can provide an accurate
measurement of the inclination angle (so we remove the inclination angle
$\iota$ from the FIM), the amplitude $h_{0}$ can then be well constrained
through GW observation (see Fig. 4). Likewise, this is also the case for the
periastron phase and polarization angle. By comparing corresponding panels in
Fig. 3 and Fig. 4, we can see that, when one of a pair of strongly correlated
parameters is removed from the FIM, the error of the other parameter will
reduce to the lower boundary of its widely dispersed distribution in Fig. 3.

The fitting result for the parameter estimation accuracy that does not include
the frequency derivative (see the green dot in Fig. 4) is

$\displaystyle\Delta\Omega$ $\displaystyle\sim$ $\displaystyle
90\left(\frac{10}{\mathrm{SNR}}\right)^{2}\mathrm{deg^{2}}\,,$ (23)
$\displaystyle\frac{\Delta h_{0}}{h_{0}}$ $\displaystyle\sim$ $\displaystyle
0.1\left(\frac{10}{\mathrm{SNR}}\right)\,,$ (24) $\displaystyle\frac{\Delta
e}{e}$ $\displaystyle\sim$ $\displaystyle
0.2\left(\frac{10}{\mathrm{SNR}}\right)\,,$ (25) $\displaystyle\frac{\Delta
P_{\mathrm{b}}}{P_{\mathrm{b}}}$ $\displaystyle\sim$ $\displaystyle 3.8\times
10^{-7}\left(\frac{10}{\mathrm{SNR}}\right)\,,$ (26)
$\displaystyle\Delta\Phi_{\mathrm{P}}$ $\displaystyle\sim$ $\displaystyle
0.1\left(\frac{10}{\mathrm{SNR}}\right)\mathrm{rad}\,.$ (27)

When the orbital period (or frequency) derivative is included in the FIM,
except for the eccentricity, the other parameters still maintain a power law
of the SNR and hold the same as before. For $30\lesssim\mathrm{SNR}\lesssim
100$, the relative errors of eccentricity and the orbital period derivative
are roughly inversely proportional to the SNR. But when $\mathrm{SNR}\gtrsim
100$, their results tend to have larger dispersion. The significant dispersion
of the relative error of $e$ is caused by the interplay between the detector’s
sensitivities at different frequencies and the evolution of $e$, which results
in a larger dispersion of $e$ and subsequently a larger dispersion of $\Delta
e/e$ for larger SNRs than for lower SNRs.

$\frac{\Delta
e}{e}\sim\left\\{\begin{array}[]{ll}0.08\left(\frac{50}{\mathrm{SNR}}\right),\quad\mathrm{SNR}\lesssim
100\\\
0.04\left(\frac{200}{\mathrm{SNR}}\right)^{0.3},\quad\mathrm{SNR}>100\end{array}\right.$
(28)
$\frac{\Delta\dot{P}_{\mathrm{b}}}{\dot{P}_{\mathrm{b}}}\sim\left\\{\begin{array}[]{l}0.01\left(\frac{50}{\mathrm{SNR}}\right),\quad\mathrm{SNR}\lesssim
100\\\
0.001\left(\frac{200}{\mathrm{SNR}}\right)^{2},\quad\mathrm{SNR}>100\end{array}\right.$
(29)

Using the standard deviation propagation similar to that in [50], for the
chirping sources, the relative error of the chirp mass and the distance of the
source can be derived by

$\frac{\Delta M_{\rm c}}{M_{\rm c}}=\left[\left(\frac{11}{5}\frac{\Delta
P_{\rm b}}{P_{\rm b}}\right)^{2}+\left(\frac{3}{5}\frac{\Delta\dot{P}_{\rm
b}}{\dot{P}_{\rm b}}\right)^{2}+\left(\frac{3}{5}\frac{dF(e)}{de}\frac{\Delta
e}{F(e)}\right)^{2}\right]^{1/2}\,,$ (30)

and

$\frac{\Delta D}{D}=\left[\left(\frac{5}{3}\frac{\Delta M_{\rm c}}{M_{\rm
c}}\right)^{2}+\left(\frac{2}{3}\frac{\Delta P_{\rm b}}{P_{\rm
b}}\right)^{2}+\left(\frac{\Delta h_{0}}{h_{0}}\right)^{2}\right]^{1/2}\,.$
(31)

The fitting equations for the above parameters are

$\displaystyle\frac{\Delta M_{\rm c}}{M_{\rm c}}$ $\displaystyle\sim$
$\displaystyle\left\\{\begin{array}[]{l}0.007\left(\frac{50}{\mathrm{SNR}}\right),\quad\mathrm{SNR}\lesssim
100\\\
0.001\left(\frac{200}{\mathrm{SNR}}\right)^{2},\quad\mathrm{SNR}>100\end{array}\right.$
(34) $\displaystyle\frac{\Delta D}{D}$ $\displaystyle\approx$
$\displaystyle\frac{\Delta h_{0}}{h_{0}}\sim
0.1\left(\frac{10}{\mathrm{SNR}}\right)\,.$ (35)

With the small relative error in the estimation of chirp mass ($\Delta M_{\rm
c}/M_{\rm c}\lesssim 1\%$), one may distinguish DNS ($M_{\rm
c}=1.22~{}M_{\odot}$) from typical WD-NS binary ($M_{\rm
c}=0.78~{}M_{\odot}$). But for a few extremely WDs with $M_{\rm c}\lesssim
1.2~{}M_{\odot}$, as discussed by [26], these WD binaries still can not be
distinguished from DNSs. A general constraint can be made from the binary
evolution models, such as the eccentric orbit and the disc binaries favor the
DNSs. The electromagnetic follow-up around the estimated sky localization can
provide further evidence to discern these binaries. As an example, we discuss
the follow-up by radio telescopes in the next section.

##  IV Radio follow-up of pulsars

In this section, we evaluate the numbers and distances of GW-detectable DNSs
in Sec. III that can be detected as radio pulsar binaries based on the pulsar
beam geometry and the detection conditions. The measurement accuracies of the
time of arrivals (TOAs) for the detectable pulsars are estimated by
considering the dominant sources of noise, namely, radiometer noise and phase
jitter noise.

###  IV.1 Pulsar beam geometry and empirical distributions

The detectability of a pulsar depends on the orientation of the pulsar’s spin
axis, magnetic axis, and the observer’s line of sight. In the pulsar beam
geometry, the intrinsic pulse width $W$ measured in longitude of rotation
satisfies the following relationship [51, 52]

$\sin^{2}\left(\frac{W}{4}\right)=\frac{\sin^{2}(\rho/2)-\sin^{2}((\xi-\alpha)/2)}{\sin\alpha\sin\xi}\,,$
(36)

where $\alpha$ is the angle between the spin axis and the magnetic axis, $\xi$
is the angle between the spin axis and the line of sight, and $\rho$ is the
angular radius of the cone-shaped beam that is centered on the magnetic axis.
Note that $\xi-\alpha$ is usually defined as the impact angle which is the
closest approach between the magnetic axis and the line of sight when the
former is rotating around the spin axis [53].

According to [54], the angle $\xi$ follows a probability distribution,
$\cos\xi\sim\mathcal{U}(0,1)$, and $\alpha\sim\mathcal{U}(0,\pi/2)$. For a
cone in a circular shape, the condition under which the radio beam can sweep
across the line of sight is

$|\xi-\alpha|<\rho\,.$ (37)

Otherwise, the beam will miss the observer. In the case that Eq. (37) is
satisfied while the right-hand side of Eq. (36) is in the range of $[0,1]$,
the intrinsic pulse width $W$ can be obtained from Eq. (36).

It has been found that the distribution of $\rho$ is strongly correlated with
pulsar’s spin period $P$, which follows an empirical formula [55]

$\rho=\left\\{\begin{array}[]{l}5.4^{\circ}(P/{\rm{s}})^{-1/2},\quad
P>30\mathrm{~{}ms}\,,\\\ 31.2^{\circ},\quad P\leq
30\mathrm{~{}ms}\,.\end{array}\right.$ (38)

To account for the fluctuations in the actual situation, in addition to the
fiducial value $\bar{\rho}$ obtained from the above equation, the observed
value of $\rho$ needs to add a random component $p$,

$\log\rho=\log\bar{\rho}+p\,.$ (39)

Here, $p\sim\mathcal{U}(-0.15,0.15)$[56].

By fitting the pulsars found in the Parkes Multibeam Pulsar Survey (PMPS) and
the Parkes High-latitude Pulsar Survey (PHPS), the probability distribution of
the spin periods for the isolated normal pulsars (NPs; $P>30\mathrm{~{}ms}$)
in our Galaxy can be expressed as a lognormal function [57],

$f_{\mathrm{dis}}(P)\propto\exp\left[-\frac{(\log
P-2.7)^{2}}{2\times(-0.34)^{2}}\right]\,.$ (40)

Here, $\log$ is the base 10 logarithm. In addition, based on the current
observations of about 20 DNSs [8, 9], we assumed that the distribution of the
spin periods for the millisecond pulsars (MSPs; $P\leq 30\mathrm{~{}ms}$) in
the DNS systems follows the uniform distribution

$P\sim\mathcal{U}(10,30)~{}\rm{ms}\,.$ (41)

The probability distribution of pulsar’s pseudo luminosity at 1400 MHz
$L_{1400}$ can be fitted as follows [58]:

$f_{\mathrm{dis}}(L_{1400})\propto\exp\left[-\frac{\left(\log
L_{1400}+1.1\right)^{2}}{2\times 0.9^{2}}\right]\,.$ (42)

Note that this formula ignores geometric and beaming factors of the radio
beam, thus $L_{1400}$ is not the true physical luminosity of the radio
emission from the pulsar [59]. On the other hand, although Eq. (42) is
obtained from the data of NPs, it has been shown that this relationship also
holds for MSPs [60]. Thus, in the following simulations, we apply this
luminosity distribution to both NPs and MSPs in the DNS systems.

Given the pulsar distance $D$, the corresponding flux density around 1400 MHz
on the Earth is

$S_{1400}=\frac{L_{1400}}{D^{2}}\,.$ (43)

###  IV.2 SNR of integrated pulse profile

According to the standard evolution scenario of the DNSs formed in isolated
binaries, their progenitor systems have experienced multiple stages of mass
transfer and two supernova explosions (e.g., [61, 62]). It is believed that
the first-born NS can be recycled to be an MSP during mass-transfer stages,
while the second-born NS is expected to appear as an NP. These are consistent
with observed features of pulsars in Galactic DNSs. Thus we assume that each
DNS in our simulation is composed of one NP and one MSP. Their spin periods
can be assigned according to the distributions given in Eq. (40) and Eq. (41),
respectively. The angular radius of the radio beam $\rho$ can be calculated
from Eq. (38) and Eq. (39). The flux density can be drawn from Eq. (42) and
Eq. (43).

Using the radiometer equation [53], the SNR of the integrated pulse profile of
a pulsar can be calculated by [57]:

$\mathrm{SNR}=\frac{S_{1400}\mathcal{G}\sqrt{n_{\mathrm{p}}\Delta\nu\tau}}{\beta
T}\sqrt{\frac{P-W_{\mathrm{eff}}}{W_{\mathrm{eff}}}}\,,$ (44)

where $S_{1400}$ is in mJy, $\mathcal{G}$ is the effective telescope gain
(K/Jy), $n_{\mathrm{p}}$ is the number of polarizations, $\Delta\nu$ is the
observation bandwidth (MHz), $\tau$ is the integration time (s), $\beta$ is
the SNR loss factor, $T$ is the system temperature (K), and $P$ is the spin
period (ms). The effective pulse width $W_{\mathrm{eff}}$ (ms) can be obtained
by

$W_{\text{eff }}=\sqrt{W^{2}+t_{\text{sam }}^{2}+\Delta
t^{2}+\tau_{\mathrm{s}}^{2}}\,,$ (45)

where $t_{\text{sam}}$ is the data sampling time in the telescope backend
system. The pulse profile broadening in a single filter channel caused by
interstellar medium dispersion $\Delta t$ and scattering $\tau_{\mathrm{s}}$
can be evaluated as follows [63, 64]

$\displaystyle\frac{\Delta t}{\mathrm{ms}}$ $\displaystyle=$
$\displaystyle~{}8.3\times 10^{6}\times\mathrm{DM}\times\frac{\Delta
f_{\mathrm{MHz}}}{f_{\mathrm{MHz}}^{3}}\,,$ (46)
$\displaystyle\log\left(\frac{\tau_{\mathrm{s}}}{\mathrm{ms}}\right)$
$\displaystyle=$ $\displaystyle-6.46+0.154\log(\mathrm{DM})$ (47)
$\displaystyle+1.07(\log\mathrm{DM})^{2}-3.86\log f_{\mathrm{GHz}}\,.$

Here, $\Delta f_{\mathrm{MHz}}$ is the frequency channel width (MHz),
$f_{\mathrm{MHz}}=1400~{}{\mathrm{MHz}}$, and
$f_{\mathrm{GHz}}=1.4~{}{\mathrm{GHz}}$. DM is the dispersion measure
($\mathrm{pc/cm^{3}}$) for which we adopt the YMW16 electron density model
[65].

We set the SNR threshold of the integrated pulse profile to be 9 as in [63],
above which the pulsar is regarded as detectable. Meanwhile, as one can see
from Eq. (44), the SNR is also determined by the technical specifications of
radio telescopes. In order to see the potential to detect the pulsars in our
simulations, in the following we consider mock searches of these pulsars based
on the Parkes [66], FAST [56, 67], SKA1-Mid [68, 69], and SKA-Mid [70], of
which the necessary technical specifications and observation parameters are
listed in Table 1. Note that the declination ranges of visible sky are
$[-90^{\circ},25^{\circ}]$ for Parkes [71], $[-14.4^{\circ},65.6^{\circ}]$ for
FAST [72], and $[-90^{\circ},45^{\circ}]$ for SKA1-Mid and SKA-Mid [73].

Table 1: The technical specifications and observation parameters for the pulsar searches based on the Parkes, FAST, SKA1-Mid, and SKA-Mid. The SKA-Mid is assumed to have the same parameters as the SKA1-Mid except for $\mathcal{G}=110~{}\mathrm{KJy^{-1}}$.  | Parkes | FAST | SKA1-Mid
---|---|---|---
Correction factor $\beta$ | 1.2 | 1.2 | 1.2
Gain $\mathcal{G}~{}(\mathrm{KJy^{-1}})$ | 0.7 | 16.5 | 15
Integration time $\tau~{}(\mathrm{s})$ | 2100 | 600 | 2100
Sampling time $t_{\mathrm{sam}}~{}(\mathrm{\mu s})$ | 250 | 50 | 64
System temperature $T~{}(\mathrm{K})$ | 25 | 25 | 30
Bandwidth $\Delta\nu~{}(\mathrm{MHz})$ | 288 | 400 | 300
Channel width $\Delta f~{}(\mathrm{MHz})$ | 3 | 0.39 | 0.02
Number of polarizations $n_{\mathrm{p}}$ | 2 | 2 | 2

###  IV.3 Number of detectable pulsar binaries

We obtain the numbers of the DNSs in Sec. III.2 that are jointly detectable by
the space-borne GW detectors and the radio telescopes. The histograms of the
numbers from 100 independent Monte Carlo simulations are shown in Fig. 5 (for
TQ), Fig. 6 (for LISA), and Fig. 7 (for TQ+LISA). The corresponding mean
values are listed in Table 2.

(a) TQ+Parkes

(b) TQ+FAST

(c) TQ+SKA1-Mid

(d) TQ+SKA-Mid

Figure 5: Histograms of the numbers of the pulsar binaries jointly detectable
by TQ and the radio telescopes (Parkes, FAST, SKA1-Mid, and SKA-Mid). MSP+NS
represents the DNSs that contain an MSP and a NS with no visible radio pulses.

(a) LISA+Parkes

(b) LISA+FAST

(c) LISA+SKA1-Mid

(d) LISA+SKA-Mid

Figure 6: Histograms as in Fig. 5, but for LISA.

(a) TQ+LISA+Parkes

(b) TQ+LISA+FAST

(c) TQ+LISA+SKA1-Mid

(d) TQ+LISA+SKA-Mid

Figure 7: Histograms as in Fig. 5, but for TQ+LISA.  Table 2: The mean numbers of the pulsar binaries from the 100 Monte Carlo simulations that can be jointly detectable by the GW detectors and the radio telescopes with the specifications in Table 1. In each parenthesis, the numbers from left to right are for TQ, LISA, TQ+LISA, respectively. | MSP-NP | MSP-NS | NP-NS
---|---|---|---
Parkes | (0.2, 0.3, 0.3) | (2.7, 4.5, 4.7) | (1.5, 2.3, 2.6)
FAST | (0.6, 0.8, 0.8) | (4.1, 6.0, 6.6) | (1.2, 2.1, 2.4)
SKA1-Mid | (1.6, 2.7, 3.0) | (14, 23, 26) | (7.3, 12, 14)
SKA-Mid | (4.4, 7.1, 8.0) | (28, 47, 55) | (13, 21, 24)

Since LISA (red curve in Fig. 2) is more sensitive than TQ (blue curve) around
$0.1-10~{}{\rm mHz}$, the mean numbers of jointly detectable pulsar binaries
by LISA and each radio telescope are larger than the corresponding ones for TQ
by $50-70\%$ in most of the cases. When operating TQ and LISA together
(TQ+LISA), these numbers will be increased by up to $17\%$ than LISA alone.

In a similar fashion, the numbers of observed pulsar binaries increase with
improved sensitivity of the pulsar searches based on the Parkes, FAST,
SKA1-Mid, and SKA-Mid. Not surprisingly, due to lower sensitivity, the
possibility of detecting short-period double pulsars with an MSP and a NP by
the Parkes is quite low; While FAST and SKA1-Mid will have good chances and
SKA-Mid will guarantee the detection of several double pulsars. All radio
telescopes have the potential to detect, at least, several MSP-NS and NP-NS
systems jointly with the GW detectors. The mean numbers for MSP-NS are larger
than the ones for NP-NS binaries in all cases, this is due to the fact that
the angular radius of the radio beam (cf., Eq. (38)) for MSP is typically
larger than NP. FAST is expected to discover 5-10 pulsar binaries detected by
TQ, LISA, and TQ+LISA; While SKA1-Mid (SKA-Mid) will further increase the
numbers for FAST by a factor of 4 (8).

###  IV.4 Distances of detectable pulsar binaries

For those DNSs in Sec. III.2 that are jointly detectable by the space-borne GW
detectors and the radio telescopes, the histograms of their distances from 100
independent Monte Carlo simulations are shown in Fig. 8, the corresponding
mean values are listed in Table 3.

The average distances of pulsar binaries increase with the increasing
sensitivity of the radio telescopes. MSP-NP binaries tend to be in the closer
region of 1-5 kpc, while NP-NS (MSP-NS) binaries are in the range of 3-7 kpc
(2-6 kpc).

Table 3:  The mean distances (in unit of kpc) of the pulsar binaries from the 100 Monte Carlo simulations that can be jointly detectable by the GW detectors and the radio telescopes. In each parenthesis, the values from left to right are for TQ, LISA, TQ+LISA, respectively. | MSP-NP | MSP-NS | NP-NS
---|---|---|---
Parkes | (0.9, 1.1, 1.2) | (1.8, 2.1, 2.2) | (3.4, 3.5, 3.6)
FAST | (1.7, 2.0, 2.1) | (2.5, 2.7, 2.8) | (4.4, 4.4, 4.5)
SKA1-Mid | (2.8, 2.8, 3.0) | (4.3, 4.4, 4.5) | (5.8, 5.8, 5.9)
SKA-Mid | (4.6, 4.6, 4.7) | (5.4, 5.5, 5.6) | (6.6, 6.7, 6.8)
Figure 8:  Histograms of the distances for the pulsars in Sec. IV.3 that are
simultaneously detectable by the GW detectors and the radio telescopes.

###  IV.5 Timing accuracies of detectable pulsars

The total rms error of TOA measurement ($\sigma_{\mathrm{t}}$) is mainly
determined by the radiometer noise ($\sigma_{\mathrm{r}}$) and the pulse phase
jitter noise ($\sigma_{\mathrm{j}}$) [74, 75],

$\sigma_{\mathrm{t}}=\sqrt{\sigma_{\mathrm{j}}^{2}+\sigma_{\mathrm{r}}^{2}}\,.$
(48)

Here,

$\sigma_{\mathrm{j}}\approx 0.28\,W_{\mathrm{eff}}\sqrt{\frac{P}{\tau}}\,,$
(49)

and

$\sigma_{\mathrm{r}}\approx\frac{W_{\mathrm{eff}}T}{S_{1400}\mathcal{G}\sqrt{2\Delta\nu\tau}}\sqrt{\frac{W_{\mathrm{eff}}}{P-W_{\mathrm{eff}}}}\,.$
(50)

Based on the equations above, we calculate $\sigma_{\mathrm{t}}$ for the
detectable MSPs and NPs, respectively, combined from 100 Monte-Carlo
simulations in Sec. IV.3. Each column of the histograms in Fig. 9 shows the
$\sigma_{\mathrm{t}}$ for the corresponding cases in Figs. 5, 6, and 7. As one
can see, the most probable values of the total rms error distributions for the
MSPs are a factor of $\sim 1.3-8.2$ smaller than that for the NPs. This is
mainly due to the narrower pulse widths of the MSPs than that of the NPs. In
addition, for SKA-Mid with the highest sensitivity, the total rms errors for
the MSPs can be as low as $70~{}\mathrm{ns}$ with most ones concentrated
around $100~{}\mathrm{\mu s}$; while $\sigma_{\mathrm{t}}$ for NPs can be as
low as $5~{}\mathrm{\mu s}$ with most ones concentrated around
$130~{}\mathrm{\mu s}$.

Figure 9: Histograms of the total rms errors of the TOA measurements for the
pulsars in Sec. IV.3 that are simultaneously detectable by the GW detectors
and the radio telescopes. In each panel, the pulsars from 100 Monte-Carlo
simulations are combined. The red and blue histograms represent the results
for the MSPs and NPs, respectively.

##  V Conclusions

In this work, we simulate the DNS population in the Galactic disk that will
merge within the next $10~{}{\rm Myr}$. The DNS merger rate is adopted from
the recent LIGO results. The numbers of the detectable DNSs and the associated
parameter estimation accuracies obtained by TQ, LISA, and TQ+LISA are
estimated by FIM in 100 Monte Carlo simulations. Based on the beam geometry
and empirical distributions of pulsar spin period and luminosity, we get the
numbers and distances of these DNSs comprising radio pulsars that can be
potentially detected by the Parkes, FAST, SKA1-Mid, and SKA-Mid. In addition,
we estimate the timing accuracies of the TOAs for these pulsars by evaluating
their dominating radiometer and jitter noises.

There are caveats caused by the assumptions used in our simulations. Firstly,
we follow [23] in setting the orbital parameters of B1913+16 at the merger as
the benchmark ($P_{\mathrm{b0}}=0.2~{}\mathrm{s}$, $e_{0}=5\times 10^{-6}$) in
the calculation of the orbital parameters of the simulated DNSs. However, when
we set J0737-3039 [76] as the benchmark, i.e.,
$P_{\mathrm{b0}}=0.18~{}\mathrm{s}$ and $e_{0}=1\times 10^{-6}$ at the merger,
the eccentricities of the simulated DNSs are generally smaller ($e\lesssim
0.04$) than the ones in Sec. III, according to the orbital evolution equations
in Sec. II. By repeating the simulation procedure in Sec. III, we find that
the distributions of the estimated parameters for the latter case remain
unchanged except for $\Delta e/e$, the overall amplitude of which increases by
a factor of 5 while the power law form in Eq. (28) holds the same. Secondly,
there are only five MSPs in DNSs are discovered to date [8], due to the lack
of prior knowledge, we assume that the distribution of their spin periods
obeys a uniform distribution (cf. Eq. 41). This will be improved when more
MSPs are discovered in future radio pulsar surveys. We also assume that the
luminosity distribution of the pulsars in tight binaries is the same as that
of the normal pulsars. If a long period of time has elapsed from the wide
binaries (at the beginning of their formation) to the present tight binaries,
the luminosity of pulsars in the latter population may be dimmer because of
the decay of the luminosity with time [63].

In this work, we only consider the SKA1-Mid and SKA-Mid in pulsar
observations. However, the SKA1-Low and SKA-Low and their precursors or
pathfinders working mainly in $50-350$ MHz may have an important contribution
in the search for pulsars in DNSs, especially when a pulsar has a steep power
spectrum in flux density [77]. Furthermore, observation of pulsars in low
frequencies can effectively improve the timing accuracy through the
measurement and mitigation of the dispersion effects, and for the pulsars near
the ecliptic plane, the solar wind effects [78].

For the Galactic DNSs, both LISA and TQ have the potential to detect their 3rd
harmonic signals (cf. Eq. (III.1)). Combining with the 2nd harmonic signals
(cf. Fig. 2), we can determine the total masses of these eccentric DNSs
through the GW frequency shift caused by the periastron advance [79]. Thus, we
can obtain the component masses of the binary by their chirp mass and total
mass. With the aid of multi-messenger observation realized by space-borne GW
detectors and radio telescopes, these ultra-compact binary systems with short
orbital periods ($\lesssim 10$ min) are expected to probe the highly
relativistic regime of astrophysics, such as constrain the equation of state
of NS with remarkable fidelity [29] or tight constraint on the parameter space
for specific alternative theories of gravity [80]. These considerations are
currently under our investigation.

###### Acknowledgements.

Y.W. gratefully acknowledges support from the National Key Research and
Development Program of China (Nos. 2022YFC2205201 and 2020YFC2201400), the
National Natural Science Foundation of China (NSFC) under Grants No. 11973024,
Major Science and Technology Program of Xinjiang Uygur Autonomous Region (No.
2022A03013-4), and Guangdong Major Project of Basic and Applied Basic Research
(Grant No. 2019B030302001). J.W.C acknowledges the support from China
Postdoctoral Science Foundation under Grant No. 2021M691146. We thank the
anonymous referees for helpful comments and suggestions.

## References

  * Hulse and Taylor [1975] R. A. Hulse and J. H. Taylor, Discovery of a pulsar in a binary system., The Astrophysical Journal Letters 195, L51 (1975).
  * Stairs [2003] I. H. Stairs, Testing General Relativity with Pulsar Timing, Living Reviews in Relativity 6, 5 (2003), arXiv:astro-ph/0307536 [astro-ph] .
  * Kramer  _et al._ [2021] M. Kramer, I. H. Stairs, _et al._ , Strong-Field Gravity Tests with the Double Pulsar, Physical Review X 11, 041050 (2021), arXiv:2112.06795 [astro-ph.HE] .
  * Amaro-Seoane  _et al._ [2017] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, _et al._ , Laser interferometer space antenna, arxiv e-prints (2017), arXiv preprint arXiv:1702.00786 (2017).
  * Amaro-Seoane  _et al._ [2022] P. Amaro-Seoane, J. Andrews, _et al._ , Astrophysics with the Laser Interferometer Space Antenna, arXiv e-prints , arXiv:2203.06016 (2022), arXiv:2203.06016 [gr-qc] .
  * Luo  _et al._ [2016] J. Luo, L.-S. Chen, H.-Z. Duan, Y.-G. Gong, S. Hu, J. Ji, Q. Liu, J. Mei, V. Milyukov, M. Sazhin, C.-G. Shao, V. T. Toth, H.-B. Tu, Y. Wang, Y. Wang, H.-C. Yeh, M.-S. Zhan, Y. Zhang, V. Zharov, and Z.-B. Zhou, TianQin: a space-borne gravitational wave detector, Classical and Quantum Gravity 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Hu and Wu [2017] W.-R. Hu and Y.-L. Wu, The Taiji Program in Space for gravitational wave physics and the nature of gravity, National Science Review 4, 685 (2017), https://academic.oup.com/nsr/article-pdf/4/5/685/31566708/nwx116.pdf .
  * Andrews and Mandel [2019] J. J. Andrews and I. Mandel, Double Neutron Star Populations and Formation Channels, The Astrophysical Journal Letters 880, L8 (2019), arXiv:1904.12745 [astro-ph.HE] .
  * Agazie  _et al._ [2021] G. Y. Agazie, M. G. Mingyar, _et al._ , The green bank northern celestial cap pulsar survey. VI. discovery and timing of PSR j1759+5036: A double neutron star binary pulsar, The Astrophysical Journal 922, 35 (2021).
  * Stovall  _et al._ [2018] K. Stovall, P. C. C. Freire, _et al._ , PALFA discovery of a highly relativistic double neutron star binary, The Astrophysical Journal 854, L22 (2018).
  * Faulkner  _et al._ [2004] A. J. Faulkner, I. H. Stairs, M. Kramer, A. G. Lyne, G. Hobbs, A. Possenti, D. R. Lorimer, R. N. Manchester, M. A. McLaughlin, N. D’Amico, F. Camilo, and M. Burgay, The Parkes Multibeam Pulsar Survey - V. Finding binary and millisecond pulsars, Monthly Notices of the Royal Astronomical Society 355, 147 (2004), arXiv:astro-ph/0408228 [astro-ph] .
  * Bagchi  _et al._ [2013] M. Bagchi, D. R. Lorimer, and S. Wolfe, On the detectability of eccentric binary pulsars, Monthly Notices of the Royal Astronomical Society 432, 1303 (2013), arXiv:1302.4914 [astro-ph.SR] .
  * Chandler [2003] A. M. Chandler,  _Pulsar searches: From radio to gamma-rays_ , Ph.D. thesis, CALIFORNIA INSTITUTE OF TECHNOLOGY (2003).
  * Ransom  _et al._ [2003] S. M. Ransom, J. M. Cordes, and S. S. Eikenberry, A New Search Technique for Short Orbital Period Binary Pulsars, Astrophys. J. 589, 911 (2003), arXiv:astro-ph/0210010 [astro-ph] .
  * ATN [2021] ATNF Pulsar Catalogue, https://www.atnf.csiro.au/research/pulsar/psrcat/ (2021), [last accessed 25-Nov-2022].
  * Manchester  _et al._ [2005] R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, The Australia Telescope National Facility Pulsar Catalogue, The Astronomical Journal 129, 1993 (2005), arXiv:astro-ph/0412641 [astro-ph] .
  * Abbott  _et al._ [2009] B. P. Abbott, R. Abbott, R. Adhikari, P. Ajith, B. Allen, G. Allen, R. S. Amin, S. B. Anderson, W. G. Anderson, and et al, LIGO: the laser interferometer gravitational-wave observatory, Reports on Progress in Physics 72, 076901 (2009).
  * Acernese  _et al._ [2008] F. Acernese, M. Alshourbagy, P. Amico, F. Antonucci, S. Aoudia, P. Astone, S. Avino, L. Baggio, G. Ballardin, and et al, Status of virgo, Classical and Quantum Gravity 25, 114045 (2008).
  * Abbott  _et al._ [2017] B. P. Abbott, R. Abbott, T. D. Abbott, F. Acernese, K. Ackley, J. Adams, LIGO Scientific Collaboration, and Virgo Collaboration, GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral, Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2020] B. Abbott, R. Abbott, T. Abbott, S. Abraham, F. Acernese, K. Ackley, C. Adams, R. Adhikari, V. Adya, C. Affeldt, _et al._ , GW190425: Observation of a compact binary coalescence with total mass $\sim$ 3.4 m⊙, The Astrophysical Journal Letters 892, L3 (2020).
  * The LIGO Scientific Collaboration  _et al._ [2021] The LIGO Scientific Collaboration, the Virgo Collaboration, and the KAGRA Collaboration, The population of merging compact binaries inferred using gravitational waves through GWTC-3, arXiv e-prints , arXiv:2111.03634 (2021), arXiv:2111.03634 [astro-ph.HE] .
  * Abbott  _et al._ [2019] B. Abbott, R. Abbott, _et al._ , GWTC-1: a gravitational-wave transient catalog of compact binary mergers observed by ligo and virgo during the first and second observing runs, Physical Review X 9, 031040 (2019).
  * Andrews  _et al._ [2020] J. J. Andrews, K. Breivik, C. Pankow, D. J. D’Orazio, and M. Safarzadeh, LISA and the Existence of a Fast-merging Double Neutron Star Formation Channel, The Astrophysical Journal Letters 892, L9 (2020), arXiv:1910.13436 [astro-ph.HE] .
  * Seto [2019] N. Seto, Search for neutron star binaries in the Local Group galaxies using LISA, Monthly Notices of the Royal Astronomical Society 489, 4513 (2019), arXiv:1909.01471 [astro-ph.HE] .
  * Riley  _et al._ [2022] J. Riley, P. Agrawal, _et al._ , Rapid stellar and binary population synthesis with COMPAS, The Astrophysical Journal Supplement Series 258, 34 (2022).
  * Lau  _et al._ [2020] M. Y. M. Lau, I. Mandel, A. Vigna-Gómez, C. J. Neijssel, S. Stevenson, and A. Sesana, Detecting double neutron stars with LISA, Monthly Notices of the Royal Astronomical Society 492, 3061 (2020), arXiv:1910.12422 [astro-ph.HE] .
  * Kyutoku  _et al._ [2019] K. Kyutoku, Y. Nishino, and N. Seto, How to detect the shortest period binary pulsars in the era of LISA, Monthly Notices of the Royal Astronomical Society 483, 2615 (2019), arXiv:1812.02177 [astro-ph.HE] .
  * Dewdney  _et al._ [2009] P. E. Dewdney, P. J. Hall, R. T. Schilizzi, and T. J. L. W. Lazio, The Square Kilometre Array, IEEE Proceedings 97, 1482 (2009).
  * Thrane  _et al._ [2020] E. Thrane, S. Osłowski, and P. D. Lasky, Ultrarelativistic astrophysics using multimessenger observations of double neutron stars with LISA and the SKA, Monthly Notices of the Royal Astronomical Society 493, 5408 (2020), arXiv:1910.12330 [astro-ph.HE] .
  * Riley  _et al._ [2019] T. E. Riley, A. L. Watts, S. Bogdanov, P. S. Ray, R. M. Ludlam, S. Guillot, Z. Arzoumanian, C. L. Baker, A. V. Bilous, D. Chakrabarty, K. C. Gendreau, A. K. Harding, W. C. G. Ho, J. M. Lattimer, S. M. Morsink, and T. E. Strohmayer, A nicer view of psr j0030+0451: Millisecond pulsar parameter estimation, The Astrophysical Journal 887, L21 (2019).
  * Hernandez Vivanco  _et al._ [2019] F. Hernandez Vivanco, R. Smith, E. Thrane, P. D. Lasky, C. Talbot, and V. Raymond, Measuring the neutron star equation of state with gravitational waves: The first forty binary neutron star merger observations, Phys. Rev. D 100, 103009 (2019).
  * Miller  _et al._ [2020] M. C. Miller, C. Chirenti, and F. K. Lamb, Constraining the Equation of State of High-density Cold Matter Using Nuclear and Astronomical Measurements, Astrophys. J. 888, 12 (2020), arXiv:1904.08907 [astro-ph.HE] .
  * Hobbs  _et al._ [2020] G. Hobbs, R. N. Manchester, _et al._ , An ultra-wide bandwidth (704 to 4 032 MHz) receiver for the Parkes radio telescope, Publications of the Astronomical Society of Australia 37, e012 (2020), arXiv:1911.00656 [astro-ph.IM] .
  * Nan  _et al._ [2011] R. Nan, D. Li, C. Jin, Q. Wang, L. Zhu, W. Zhu, H. Zhang, Y. Yue, and L. Qian, The Five-Hundred Aperture Spherical Radio Telescope (fast) Project, International Journal of Modern Physics D 20, 989 (2011), arXiv:1105.3794 [astro-ph.IM] .
  * Dewdney  _et al._ [2016] P. Dewdney, W. Turner, R. Braun, J. Santander-Vela, M. Waterson, and G.-H. Tan,  _SKA1 System Baseline Design_ (Square Kilometre Array Organisation, Jodrell Bank, UK, 2016).
  * Weltman  _et al._ [2020] A. Weltman, P. Bull, _et al._ , Fundamental physics with the Square Kilometre Array, Publications of the Astronomical Society of Australia 37, e002 (2020), arXiv:1810.02680 [astro-ph.CO] .
  * Nelemans  _et al._ [2001] G. Nelemans, L. R. Yungelson, and S. F. Portegies Zwart, The gravitational wave signal from the Galactic disk population of binaries containing two compact objects, Astronomy and Astrophysics 375, 890 (2001), arXiv:astro-ph/0105221 [astro-ph] .
  * Peters [1964] P. C. Peters, Gravitational Radiation and the Motion of Two Point Masses, Physical Review 136, 1224 (1964).
  * Moreno-Garrido  _et al._ [1995] C. Moreno-Garrido, E. Mediavilla, and J. Buitrago, Gravitational radiation from point masses in elliptical orbits: spectral analysis and orbital parameters, Monthly Notices of the Royal Astronomical Society 274, 115 (1995).
  * Cornish and Rubbo [2003] N. J. Cornish and L. J. Rubbo, LISA response function, Phys. Rev. D 67, 022001 (2003).
  * Hu  _et al._ [2018] X.-C. Hu, X.-H. Li, Y. Wang, W.-F. Feng, M.-Y. Zhou, Y.-M. Hu, S.-C. Hu, J.-W. Mei, and C.-G. Shao, Fundamentals of the orbit and response for TianQin, Classical and Quantum Gravity 35, 095008 (2018), arXiv:1803.03368 [gr-qc] .
  * Robson  _et al._ [2019] T. Robson, N. J. Cornish, and C. Liu, The construction and use of LISA sensitivity curves, Classical and Quantum Gravity 36, 105011 (2019), arXiv:1803.01944 [astro-ph.HE] .
  * Huang  _et al._ [2020] S.-J. Huang, Y.-M. Hu, V. Korol, P.-C. Li, Z.-C. Liang, Y. Lu, H.-T. Wang, S. Yu, and J. Mei, Science with the TianQin Observatory: Preliminary results on Galactic double white dwarf binaries, Phys. Rev. D 102, 063021 (2020), arXiv:2005.07889 [astro-ph.HE] .
  * Tinto and de Araujo [2016] M. Tinto and J. C. N. de Araujo, Coherent observations of gravitational radiation with LISA and gLISA, Phys. Rev. D 94, 081101 (2016), arXiv:1608.04790 [astro-ph.IM] .
  * Finn [1992] L. S. Finn, Detection, measurement, and gravitational radiation, Phys. Rev. D 46, 5236 (1992), arXiv:gr-qc/9209010 [gr-qc] .
  * Cutler and Vecchio [1998] C. Cutler and A. Vecchio, Lisa’s angular resolution for monochromatic sources, in  _AIP Conference Proceedings_ , Vol. 456 (American Institute of Physics, 1998) pp. 95–100.
  * Cutler [1998] C. Cutler, Angular resolution of the LISA gravitational wave detector, Phys. Rev. D 57, 7089 (1998), arXiv:gr-qc/9703068 [gr-qc] .
  * Shah  _et al._ [2012] S. Shah, M. van der Sluys, and G. Nelemans, Using electromagnetic observations to aid gravitational-wave parameter estimation of compact binaries observed with LISA, Astronomy and Astrophysics 544, A153 (2012), arXiv:1207.6770 [astro-ph.IM] .
  * Mikóczi  _et al._ [2012] B. Mikóczi, B. Kocsis, P. Forgács, and M. Vasúth, Parameter estimation for inspiraling eccentric compact binaries including pericenter precession, Phys. Rev. D 86, 104027 (2012).
  * Berti  _et al._ [2005] E. Berti, A. Buonanno, and C. M. Will, Estimating spinning binary parameters and testing alternative theories of gravity with LISA, Phys. Rev. D 71, 084025 (2005), arXiv:gr-qc/0411129 [gr-qc] .
  * Gil [1981] J. Gil, The frequency dependence of component separation in pulsar integrated profiles. part 1. low frequencies, Acta Physica Polonica. Series B 12, 1081 (1981).
  * Gil  _et al._ [1984] J. Gil, P. Gronkowski, and W. Rudnicki, Geometry of the emission region of PSR 0950+08., Astronomy and Astrophysics 132, 312 (1984).
  * Lorimer and Kramer [2004] D. R. Lorimer and M. Kramer,  _Handbook of Pulsar Astronomy_ , Vol. 4 (Cambridge University Press, 2004).
  * Gil and Han [1996] J. A. Gil and J. L. Han, Geometry of Pulsar Emission and Pulse Width Distribution, Astrophys. J. 458, 265 (1996).
  * Kramer  _et al._ [1998] M. Kramer, K. M. Xilouris, D. R. Lorimer, O. Doroshenko, A. Jessner, R. Wielebinski, A. Wolszczan, and F. Camilo, The Characteristics of Millisecond Pulsar Emission. I. Spectra, Pulse Shapes, and the Beaming Fraction, Astrophys. J. 501, 270 (1998), arXiv:astro-ph/9801177 [astro-ph] .
  * Smits  _et al._ [2009a] R. Smits, D. R. Lorimer, M. Kramer, R. Manchester, B. Stappers, C. J. Jin, R. D. Nan, and D. Li, Pulsar science with the Five hundred metre Aperture Spherical Telescope, Astronomy and Astrophysics 505, 919 (2009a), arXiv:0908.1689 [astro-ph.IM] .
  * Lorimer  _et al._ [2006] D. R. Lorimer, A. J. Faulkner, A. G. Lyne, R. N. Manchester, M. Kramer, M. A. McLaughlin, G. Hobbs, A. Possenti, I. H. Stairs, F. Camilo, M. Burgay, N. D’Amico, A. Corongiu, and F. Crawford, The Parkes Multibeam Pulsar Survey - VI. Discovery and timing of 142 pulsars and a Galactic population analysis, Monthly Notices of the Royal Astronomical Society 372, 777 (2006), arXiv:astro-ph/0607640 [astro-ph] .
  * Faucher-Giguère and Kaspi [2006] C.-A. Faucher-Giguère and V. M. Kaspi, Birth and Evolution of Isolated Radio Pulsars, Astrophys. J. 643, 332 (2006), arXiv:astro-ph/0512585 [astro-ph] .
  * Arzoumanian  _et al._ [2002] Z. Arzoumanian, D. F. Chernoff, and J. M. Cordes, The Velocity Distribution of Isolated Radio Pulsars, Astrophys. J. 568, 289 (2002), arXiv:astro-ph/0106159 [astro-ph] .
  * Bagchi  _et al._ [2011] M. Bagchi, D. R. Lorimer, and J. Chennamangalam, Luminosities of recycled radio pulsars in globular clusters, Monthly Notices of the Royal Astronomical Society 418, 477 (2011), arXiv:1107.4521 [astro-ph.SR] .
  * Tauris  _et al._ [2017] T. M. Tauris, M. Kramer, P. C. C. Freire, N. Wex, H. T. Janka, N. Langer, P. Podsiadlowski, E. Bozzo, S. Chaty, M. U. Kruckow, E. P. J. van den Heuvel, J. Antoniadis, R. P. Breton, and D. J. Champion, Formation of Double Neutron Star Systems, Astrophys. J. 846, 170 (2017), arXiv:1706.09438 [astro-ph.HE] .
  * Shao and Li [2018] Y. Shao and X.-D. Li, On the Role of Supernova Kicks in the Formation of Galactic Double Neutron Star Systems, Astrophys. J. 867, 124 (2018), arXiv:1810.03324 [astro-ph.HE] .
  * Bates  _et al._ [2014] S. D. Bates, D. R. Lorimer, A. Rane, and J. Swiggum, PSRPOPPy: an open-source package for pulsar population simulations, Monthly Notices of the Royal Astronomical Society 439, 2893 (2014), arXiv:1311.3427 [astro-ph.IM] .
  * Bhat  _et al._ [2004] N. D. R. Bhat, J. M. Cordes, F. Camilo, D. J. Nice, and D. R. Lorimer, Multifrequency Observations of Radio Pulse Broadening and Constraints on Interstellar Electron Density Microstructure, Astrophys. J. 605, 759 (2004), arXiv:astro-ph/0401067 [astro-ph] .
  * Yao  _et al._ [2017] J. M. Yao, R. N. Manchester, and N. Wang, A New Electron-density Model for Estimation of Pulsar and FRB Distances, Astrophys. J. 835, 29 (2017), arXiv:1610.09448 [astro-ph.GA] .
  * Manchester  _et al._ [2001] R. N. Manchester, A. G. Lyne, F. Camilo, J. F. Bell, V. M. Kaspi, N. D’Amico, N. P. F. McKay, F. Crawford, I. H. Stairs, A. Possenti, M. Kramer, and D. C. Sheppard, The Parkes multi-beam pulsar survey - I. Observing and data analysis systems, discovery and timing of 100 pulsars, Monthly Notices of the Royal Astronomical Society 328, 17 (2001), arXiv:astro-ph/0106522 [astro-ph] .
  * Huang and Wang [2020] W. J. Huang and H. G. Wang, Role of the Fan Beam Model in Population Synthesis of Isolated Radio Pulsars, Astrophys. J. 905, 144 (2020), arXiv:2010.13127 [astro-ph.HE] .
  * Swart  _et al._ [2022] G. P. Swart, P. E. Dewdney, and A. Cremonini, Highlights of the SKA1-Mid telescope architecture, Journal of Astronomical Telescopes, Instruments, and Systems 8, 011021 (2022).
  * Cieślar  _et al._ [2020] M. Cieślar, T. Bulik, and S. Osłowski, Markov Chain Monte Carlo population synthesis of single radio pulsars in the Galaxy, Monthly Notices of the Royal Astronomical Society 492, 4043 (2020).
  * Smits  _et al._ [2009b] R. Smits, M. Kramer, B. Stappers, D. R. Lorimer, J. Cordes, and A. Faulkner, Pulsar searches and timing with the square kilometre array, Astronomy and Astrophysics 493, 1161 (2009b), arXiv:0811.0211 [astro-ph] .
  * Manchester  _et al._ [2013] R. N. Manchester, G. Hobbs, _et al._ , The Parkes Pulsar Timing Array Project, Publications of the Astronomical Society of Australia 30, e017 (2013), arXiv:1210.6130 [astro-ph.IM] .
  * Qian  _et al._ [2019] L. Qian, Z. Pan, D. Li, G. Hobbs, W. Zhu, P. Wang, Z. Liu, Y. Yue, Y. Zhu, H. Liu, D. Yu, J. Sun, P. Jiang, G. Pan, H. Li, H. Gan, R. Yao, X. Xie, F. Camilo, A. Cameron, L. Zhang, and S. Wang, The first pulsar discovered by FAST, Science China Physics, Mechanics, and Astronomy 62, 959508 (2019), arXiv:1903.06318 [astro-ph.IM] .
  * Janssen  _et al._ [2015] G. Janssen, G. Hobbs, M. McLaughlin, C. Bassa, A. Deller, M. Kramer, K. Lee, C. Mingarelli, P. Rosado, S. Sanidas, A. Sesana, L. Shao, I. Stairs, B. Stappers, and J. P. W. Verbiest, Gravitational Wave Astronomy with the SKA, in  _Advancing Astrophysics with the Square Kilometre Array (AASKA14)_ (2015) p. 37, arXiv:1501.00127 [astro-ph.IM] .
  * Cordes and Shannon [2010] J. M. Cordes and R. M. Shannon, A Measurement Model for Precision Pulsar Timing, arXiv e-prints , arXiv:1010.3785 (2010), arXiv:1010.3785 [astro-ph.IM] .
  * Wang [2015] Y. Wang, Noise in pulsar timing arrays, in  _Journal of Physics Conference Series_, Journal of Physics Conference Series, Vol. 610 (2015) p. 012019, arXiv:1505.00402 [astro-ph.IM] .
  * Kramer  _et al._ [2006] M. Kramer, I. H. Stairs, R. N. Manchester, M. A. McLaughlin, A. G. Lyne, R. D. Ferdman, M. Burgay, D. R. Lorimer, A. Possenti, N. D’Amico, J. M. Sarkissian, G. B. Hobbs, J. E. Reynolds, P. C. C. Freire, and F. Camilo, Tests of General Relativity from Timing the Double Pulsar, Science 314, 97 (2006), arXiv:astro-ph/0609417 [astro-ph] .
  * Swainston  _et al._ [2021] N. A. Swainston, N. D. R. Bhat, M. Sokolowski, S. J. McSweeney, S. Kudale, S. Dai, K. R. Smith, I. S. Morrison, R. M. Shannon, W. van Straten, M. Xue, S. M. Ord, S. E. Tremblay, B. W. Meyers, A. Williams, G. Sleap, M. Johnston-Hollitt, D. L. Kaplan, S. J. Tingay, and R. B. Wayth, Discovery of a Steep-spectrum Low-luminosity Pulsar with the Murchison Widefield Array, The Astrophysical Journal Letters 911, L26 (2021), arXiv:2104.10018 [astro-ph.HE] .
  * Kumar  _et al._ [2022] P. Kumar, S. M. White, K. Stovall, J. Dowell, and G. B. Taylor, Pulsar observations at low frequencies: applications to pulsar timing and solar wind models, Monthly Notices of the Royal Astronomical Society 511, 3937 (2022), arXiv:2202.01251 [astro-ph.HE] .
  * Seto [2001] N. Seto, Proposal for Determining the Total Masses of Eccentric Binaries Using Signature of Periastron Advance in Gravitational Waves, Phys. Rev. Lett. 87, 251101 (2001), arXiv:astro-ph/0111107 [astro-ph] .
  * Miao  _et al._ [2021] X. Miao, H. Xu, L. Shao, C. Liu, and B.-Q. Ma, Stringent Tests of Gravity with Highly Relativistic Binary Pulsars in the Era of LISA and SKA, Astrophys. J. 921, 114 (2021), arXiv:2107.05812 [astro-ph.HE] .