Modeling the Galactic Neutron Star Population for Use in Continuous Gravitational-wave Searches

Within an order of magnitude, the age of the Milky Way is ∼10¹⁰ yr and it has a Galactic supernovae rate of ∼1 per century (Diehl et al. 2006). We can therefore estimate that N₀ ∼ 10⁸ neutron stars (NS) have been born in our Galaxy to date. Of that number, a relatively small fraction are known through electromagnetic searches—a few thousand mostly radio pulsars—see, for example, the ATNF Pulsar Database (Manchester et al. 2005). Gravitational waves (GW) may be a means to discover some of the remaining unknown NSs and study the distribution of their shapes.

Any rotating NS with asymmetric deformations will produce continuous gravitational waves (CGWs) via quadrupole radiation (Zimmermann & Szedenits 1979; Lasky 2015) and the observed background of CGWs from GW detectors may reveal unknown NSs (Caride et al. 2019; Abbott et al. 2019a). A rotating NS radiates CGWs with strain amplitude h₀ according to

$$\begin{eqnarray}&&{h}_{0}=\displaystyle \frac{4{\pi }^{2}G}{{c}^{4}}\displaystyle \frac{{I}_{{zz}}{f}_{\mathrm{GW}}^{2}}{d}\epsilon ,\end{eqnarray} \tag{ 1 }$$

where d is the distance to the source and the gravitational-wave (GW) frequency is f_GW = 2ν for an NS rotating with spin frequency ν (Riles 2017; Abbott et al. 2019a). This relation is notable in that it is linearly dependent on the NS's ellipticity,

$$\begin{eqnarray}&&\epsilon =\sqrt{\displaystyle \frac{8\pi }{15}}\displaystyle \frac{{Q}_{22}}{{I}_{{zz}}}=\displaystyle \frac{{I}_{{xx}}-{I}_{{yy}}}{{I}_{{zz}}},\end{eqnarray} \tag{ 2 }$$

defined here in either terms of the quadrupole moment Q₂₂ or fractional difference in principle moments of inertia (Owen 2005; Riles 2017).

We expect a distribution of ellipticities across Galactic NSs. The maximum allowed ellipticity may be limited by the breaking strain of the NS crust to ≲ few × 10⁻⁶ (Ushomirsky et al. 2000; Horowitz & Kadau 2009; Gittins et al. 2020). In order for an NS to support a larger , there may need to be an exotic solid phase in the core, such as crystalline quark matter (Owen 2005; Johnson-McDaniel & Owen 2013). Constraining the ellipticity further, observations of millisecond pulsars (MSPs) suggest that the NS ellipticity reaches a minimum near ≈ 10⁻⁹ (Woan et al. 2018).

CGWs from Galactic NS are expected to be ∼10⁻⁴ times lower in amplitude than the GW signal from the binary mergers of compact objects (Riles 2017). There have been many CGW searches from known pulsars, see, for example, Abbott et al. (2019c). Furthermore, one can gain sensitivity to weak CGW signals by integrating for a long time. For known pulsars, however, radio or X-ray spin-down luminosity place an observational limit on the power in GW radiation. Alternatively, there are a number of all-sky searches for CGWs from unknown NSs (Abbott et al. 2019a; Steltner et al. 2021; Dergachev & Papa 2020, 2021; Dergachev & Alessandra Papa 2021). An unknown NS could be a strong CGW source with an unconstrained spin-down power. However, to date no CGW signals have been detected.

The physics implications of negative CGW searches are presently unclear, but we can use these results to infer some interesting limits on the NS population. Because the CGW strain amplitude depends inversely on the distance to the NS (h₀ ∼ 1/d), the lack of detected CGWs constrains the NS population within that distance from Earth. Of course, the strain amplitude also depends on the spin frequency of the NS producing the GWs as f_GW = 2ν and the ellipticity of the NS (${h}_{0}\sim {f}_{\mathrm{GW}}^{2}\epsilon$). Assuming a spatial distribution and spin distribution for Galactic NSs near Earth, we can therefore infer limits on the ellipticities of the NS population producing CGWs within a distance d from Earth. Doing so would allow us to estimate how many NSs are actually being probed by a given CGW search.

In this paper, we develop a very simple population model of the spatial and spin distribution of Galactic NSs. We then use this model to estimate the number of Galactic NSs probed by recent CGW searches. We detail our distribution choices in Section 2. These distributions can be used to help optimize future searches and to infer the physics implications of search results. In Section 3 we use our NS distribution, along with search limits on h₀, to infer a distribution of upper limits on NS ellipticities. Lastly, we discuss the implications and conclusions of our findings and possible future studies to improve these limits in Section 4. We also discuss our assumptions and their impact on the results; for example, the time-dependence of the CGW source, the effects of binary pairs of NS, and the NS birth history in the Milky Way.

In this section, we develop a simple model for the distribution of NSs in the galaxy to provide a simple first estimate of the number of NSs probed to a variety of ellipticities via CGW data. We explain our calculations of the maximum distance from Earth that is probed at a given frequency in Section 2.1. We then give our assumptions and calculations for the spatial distribution of NSs in Section 2.2 and detail our choice of spin-frequency distribution in Section 2.3. Finally, we show the calculation of the unknown NS population in Section 2.4.

2.1. Gravitational-wave Strain Data

Equation (2) gives the definition of , characterized by an asymmetric deformation on the surface of an NS. This asymmetry will cause a rapidly rotating NS to emit CGWs with a strain amplitude given by Equation (1). The strain amplitude h₀ is sensitive to the frequency of the GW signal and has a complicated behavior.

We use data from Abbott et al. (2019a), who presented a detailed analysis of their CGW search limits on h₀ as a function of f_GW at 95% confidence. Within this data are the constraints from the three pipelines SkyHough (Krishnan et al. 2004), Frequency-Hough (Astone et al. 2014), and TDFstat (Jaranowski et al. 1998), which have different sensitivities in the range of frequencies considered by Abbott et al. (2019a). Because there exists some overlap in the strain among the three pipelines, we define a grid of 20–1922 Hz and take the smallest of the three h₀ at each frequency to use in our calculations. We plot this and the data from the three pipelines in Figure 1.

Zoom In Zoom Out Reset image size

By simple inversion of Equation (1), we can solve for the maximum distance from Earth to which an NS with frequency f_GW and ellipticity has been excluded. Explicitly,

$$\begin{eqnarray}&&d({f}_{\mathrm{GW}},\epsilon )=\displaystyle \frac{4{\pi }^{2}G}{{c}^{4}}\displaystyle \frac{{I}_{{zz}}{f}_{\mathrm{GW}}^{2}\epsilon }{{h}_{0}({f}_{\mathrm{GW}})}.\end{eqnarray} \tag{ 3 }$$

The distance can be easily obtained by fixing a desired value for , using a canonical value for I_zz = 10⁴⁵ g cm², and then choosing a desired f_GW value. For illustration, we plot the distance versus f_GW for various values in Figure 2. Note that current GW interferometers are insensitive below 20 Hz.

Zoom In Zoom Out Reset image size

Because the maximum distance in Figure 2 goes as $d\sim {f}_{\mathrm{GW}}^{2}\epsilon$, the greatest distances probed are for f_GW ≳ 1000 Hz and large ellipticity ∼ 10⁻⁵. In particular, NSs with ∼ 10⁻⁵ are excluded to around d ≈ 20 kpc (likely the entire Galactic disk) for f_GW ≳ 1000 Hz and to d ≈ 500 pc for f_GW = 100 Hz. For NSs with ∼ 10⁻⁶, closer to the breaking strain of the crust, they are excluded to approximately d ≈ 2 kpc for f_GW ≳ 1000 Hz and around d ≈ 50 pc for f_GW = 100 Hz. For NSs with ∼ 10⁻⁹ they are excluded up to d ≈ 2 pc for f_GW ≳ 1000 Hz and d ≪ 1 pc for f_GW = 100 Hz.

2.2. Spatial Distribution

The NSs that fall within the range of the GW detector as defined by Equation (3) reside in the Galactic disk. The spatial distribution of Galactic NSs is believed to approximately follow an exponential distribution in the vertical direction above the disk and a Gaussian-like distribution in the radial direction (Binney & Merrifield 1998; Faucher-Giguère & Loeb 2010; Taani et al. 2012). We shall adopt the following equation for the 3D density of NS in the Galaxy

$$\begin{eqnarray}&&{\rho }_{c}({r}_{c},z)=\displaystyle \frac{{N}_{0}}{4\pi {\sigma }_{r}^{2}{z}_{0}}\exp \left[-\displaystyle \frac{{r}_{c}^{2}}{2{\sigma }_{r}^{2}}\right]\exp \left[-\displaystyle \frac{| z| }{{z}_{0}}\right],\end{eqnarray} \tag{ 4 }$$

where r_c is the cylindrical radius from the Galactic center, σ_r is a radius parameter, N₀ is the total number of NSs, and z₀ is the disk thickness. For σ_r we adopt a value of 5 kpc as in Faucher-Giguère & Loeb (2010) and we use N₀ = 10⁸ as discussed in Section 1. For normal stars, z₀ is often in the range of 0.5–1.0 kpc (Binney & Merrifield 1998; Binney & Tremaine 2008). However, NSs may follow a different distribution due to supernova kicks. Therefore, to probe a wider range of models for the distribution, we choose to vary z₀ for the values in Table 1.

Table 1. Adopted Parameter Values for the Population Distributions Used in This Study

Parameter	Symbol	Adopted Value(s)
Radius parameter	σr	5 kpc
Disk thickness	z0	0.1, 2.0, 4.0 kpc
Distance to GC	Re	8.25 kpc
Normalization	N0	108 stars

Download table as:  ASCIITypeset image

We then perform a coordinate transformation from the cylindrical (r_c , z) to the average 3D distance from Earth, d. This can be done by first transforming r_c to be centered on the Earth via $\vec{{r}_{c}}\longrightarrow \vec{r}+\vec{{R}_{e}}$, where R_e = 8.25 kpc is the distance from the Galactic center to Earth (Gravity Collaboration et al. 2019). Plugging this substitution into Equation (4) and integrating over the angular direction gives us

$$\begin{eqnarray}&&\rho ^{\prime} (r,z)=\displaystyle \frac{{N}_{0}{{\rm{e}}}^{-\tfrac{| z| }{{z}_{0}}}}{2{\sigma }_{r}^{2}{z}_{0}}{I}_{0}\left(\displaystyle \frac{{R}_{e}r}{{\sigma }_{r}^{2}}\right)\exp \left[-\displaystyle \frac{({r}^{2}+{R}_{e}^{2})}{2{\sigma }_{r}^{2}}\right],\end{eqnarray} \tag{ 5 }$$

where I₀ is the modified Bessel function. We note that this distribution is normalized to N₀,

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }{dz}{\int }_{0}^{\infty }{rdr}\rho ^{\prime} (r,z)={N}_{0}.\end{eqnarray} \tag{ 6 }$$

Now, using d for the 3D distance from Earth, $d=\sqrt{{r}^{2}+{z}^{2}}$, we can then arrive at an average 1D density ρ(d) by performing the following integral

$$\begin{eqnarray}&&\rho (d)={\int }_{-\infty }^{\infty }{dz}{\int }_{0}^{\infty }{rdr}\rho ^{\prime} (r,z)\delta (\sqrt{{r}^{2}+{z}^{2}}-d).\end{eqnarray} \tag{ 7 }$$

Integrating over the radial coordinate first, we then recast z in terms of a scaled variable x = z/d. With this, we have arrived at the probability density distribution

$$\begin{eqnarray}\begin{array}{rcl}\rho (d) & = & \displaystyle \frac{{N}_{0}{d}^{2}}{{\sigma }_{r}^{2}{z}_{0}}{\displaystyle \int }_{0}^{1}\exp \left[-\displaystyle \frac{{xd}}{{z}_{0}}\right]{I}_{0}\left[\displaystyle \frac{{R}_{e}d\sqrt{1-{x}^{2}}}{{\sigma }_{r}^{2}}\right]\\ & & \times \exp \left[-\displaystyle \frac{{R}_{e}^{2}+{d}^{2}(1-{x}^{2})}{2{\sigma }_{r}^{2}}\right]{dx},\end{array}\end{eqnarray} \tag{ 8 }$$

which gives the likelihood that an NS is a distance d from Earth. We plot the probability density distribution within 30 kpc of Earth in the left panel of Figure 3. Finally, we integrate Equation (8) to arrive at the cumulative distribution function N(d) at a distance d, defined here as

$$\begin{eqnarray}&&N(d)={\int }_{0}^{d}\rho (y){dy}.\end{eqnarray} \tag{ 9 }$$

We plot N(d) in the right panel of Figure 3.

Zoom In Zoom Out Reset image size

2.3. Spin-frequency Distribution

A neutron star spinning with frequency ν emits CGWs with frequency f_GW = 2ν according to Equation (1). The observed frequency of many MSPs is believed to be due to spin-up in an NS's low-mass X-ray binary phase (e.g., Radhakrishnan & Srinivasan 1982; Wijnands & van der Klis 1998; Papitto et al. 2013). Over the lifetime of this phase, asymmetric electron-capture reaction layers on the accreting neutron star may lead to an asymmetric deformation (Bildsten 1998; Ushomirsky et al. 2000) because of an asymmetry in the temperature distribution of the NS's magnetic field. The observed distribution of ν largely depends on the spin evolution in this phase (Bhattacharyya 2021). After this phase concludes, the spinning NS continues to emit CGWs, which affects the time evolution of both and ν.

As a simple starting point, we assume that the true distribution of Galactic NS spin frequencies is the same as the observed spin-frequency distribution of 2811 pulsars from the ATNF Pulsar Database (Manchester et al. 2005). In the left panel of Figure 4, we show a histogram of the f_GW expected from pulsars in the database. Within the sample, there are 489 pulsars with a spin frequency above 10 Hz that produce f_GW > 20 Hz and fall within GW detector sensitivity. We note the maximally rotating pulsar in the catalog is PSR J1748-2446ad (Hessels et al. 2006) and has ν ≈ 716 Hz (f_GW ≈ 1432 Hz), which is the fastest rotating pulsar yet observed.

Zoom In Zoom Out Reset image size

Using a kernel-density estimator (KDE) from Virtanen et al. (2020; package documentation in SciPy Community 2021) we calculate a probability distribution function (PDF) of f_GW from the pulsar distribution, Φ. This method uses a Gaussian with Scott's rule (Scott 2015) of n^−1/5 and is convenient as it produces a continuous function of f_GW, which makes solving integrals with it much easier. We then normalize Φ to unity

$$\begin{eqnarray}&&1={\int }_{-\infty }^{\infty }{\rm{\Phi }}(\mathrm{log}f)d(\mathrm{log}f)\end{eqnarray} \tag{ 10 }$$

giving us a normalized PDF of f_GW. We note that Φ(f) is continuous, but drops off quickly after f_GW ≳ 1430 Hz because there are no observed pulsars with spin frequencies above ν = 716 Hz. Above 2000 Hz, we set Φ(f > 2000 Hz) = 0 since the data from Section 2.1 does not go above f_GW = 2000 Hz. We show the PDF in the right panel of Figure 4. Note that the KDE was fit in ${\mathrm{log}}_{10}({f}_{\mathrm{GW}})$ to achieve higher accuracy at low values of f_GW. Henceforth, all log₁₀ shall be shortened to just log.

2.4. Unknown NS Population Estimate

Ground-based GW detectors are insensitive for f_GW < 20 Hz and based on our catalog, ≳83% of known pulsars are therefore spinning too slowly to be detectable via CGWs. As a first-order estimate we can say that CGW searches are at most sensitive to the remaining 17%, that is, approximately ∼17 million NSs throughout the Galaxy,

$$\begin{eqnarray}&&0.17\approx {\int }_{\mathrm{log}(20{Hz})}^{\infty }{\rm{\Phi }}(\mathrm{log}f)d(\mathrm{log}f).\end{eqnarray} \tag{ 11 }$$

In practice, ground-based GW detector sensitivity rapidly declines below f_GW < 100 Hz, so it may be difficult to detect slowly spinning stars.

We can now estimate the number of unknown NSs probed at a given , defined as N_⋆. We define our grid of values to range from $\mathrm{log}(\epsilon )=[-9,-5]$. Then, we perform the following integral to calculate N_⋆

$$\begin{eqnarray}\begin{array}{rcl}{N}_{\star }(\epsilon ) & = & {\displaystyle \int }_{\mathrm{log}({f}_{1})}^{\mathrm{log}({f}_{2})}N(d(\mathrm{log}f,\epsilon ))\\ & & \times {\rm{\Phi }}(\mathrm{log}f){\rm{d}}(\mathrm{log}f),\end{array}\end{eqnarray} \tag{ 12 }$$

where f₁ = 20 Hz is the minimum sensitivity of GW detectors and f₂ is the spin-down limited frequency that an NS could produce a detectable signal in GWs. This value is dependent on the particular search and will mostly affect the number of fast-spinning highly elliptical NSs. In Abbott et al. (2019a), the maximum spin-down considered in their search was $| \dot{f}| =1\times {10}^{-8}\,\mathrm{Hz}\ {{\rm{s}}}^{-1}$ for the SkyHough and TDFStat pipelines and $| \dot{f}| =2\times {10}^{-9}\,\mathrm{Hz}\ {{\rm{s}}}^{-1}$ for Frequency-Hough. It should be noted that the limit on $| \dot{f}|$ for Frequency−Hough quoted here is the value considered in the range of 512-1024 Hz, whereas in the range of 20–512 Hz $| \dot{f}| =1\times {10}^{-8}\,\mathrm{Hz}\ {{\rm{s}}}^{-1}$. For an NS with ellipticity

$$\begin{eqnarray}&&{f}_{\max }=225\,\ \mathrm{Hz}\left(\displaystyle \frac{| \dot{{\text{}}f}{| }^{1/5}}{{\epsilon }^{2/5}}\right),\end{eqnarray} \tag{ 13 }$$

where $| \dot{f}|$ is in units of [Hz s⁻¹] (Abbott et al. 2019b). Using the strain data from Abbott et al. (2019a), we will adopt ${f}_{2}={f}_{\max }$ for $| \dot{f}| =2\times {10}^{-9}\,\mathrm{Hz}\ {{\rm{s}}}^{-1}$ for simplicity since the overall search is limited by the smallest maximum value of $| \dot{f}|$.

3.1. Neutron Star Estimates

In Figure 5 are model populations from solving Equation (12) for the chosen values of z₀ and these show the number of NSs probed to a given ellipticity. We also record the characteristic values for each of the models in Figure 5 in Table 2. We find that between ≈3–5.3 × 10⁵ stars of the 1.7 × 10⁷ NSs with f_GW > 20 Hz have been probed to ∼ 10⁻⁵–corresponding to between ≈1.7 and 3.1% of the model population. By contrast, only between ≈1 and 5.3 × 10⁴ stars are probed to ∼ 10⁻⁶, corresponding to only between ≈0.06 and 0.31% of the visible population. For comparison, we show the limits on from the analysis of Abbott et al. (2017) in Figure 5 using CGW searches of known pulsars at known f_GW. This data probes significantly fewer NSs at ellipticities above ≳ 10⁻⁷ than using our method. This is because most of the NSs in the galaxy are unknown. We discuss this further in Section 4.

Zoom In Zoom Out Reset image size

Table 2. Estimates for Number of NSs Probed at the Value of Ellipticity for Different Values of Disk Thickness z₀

${\mathrm{log}}_{10}(\epsilon )$	z0 = 0.1 kpc	z0 = 2.0 kpc	z0 = 4.0 kpc
−5.00	5.3 × 105	4.1 × 105	3.0 × 105
−5.25	3.8 × 105	2.7 × 105	1.8 × 105
−5.50	2.1 × 105	1.2 × 105	7.6 × 104
−5.75	1.2 × 105	5.4 × 104	3.2 × 104
−6.00	5.3 × 104	1.8 × 104	1.0 × 104
−6.25	2.0 × 104	4.8 × 103	2.6 × 103
−6.50	6.5 × 103	1.0 × 103	540
−6.75	1.8 × 103	190	99
−7.00	430	35	18
−7.25	95	6	3
−7.50	19	1	1
−7.75	4	0	0
−8.00	1	0	0
−8.25	0	0	0
−8.50	0	0	0
−8.75	0	0	0
−9.00	0	0	0

Download table as:  ASCIITypeset image

3.2. Effects of Improved Strain Sensitivity

The strain amplitude h₀ used in this study is limited by the sensitivity of GW interferometers and the parameters of the search. We can see the effects of improving the sensitivity of h₀ directly in Equation (3) in that the distance we can be sensitive to will increase with decreasing h₀. This will then increase our estimate for the total number of NSs probed at a given ellipticity. As an example of this effect, we test what would happen to N_⋆ should a new search reduce h₀ in either the high-frequency (f_GW ≥ 1000Hz) or low-frequency (f_GW ≤ 100Hz) regimes by a factor of 2.

We present the predicted values of N_⋆ for improved detector sensitivity in Figure 6 for a disk model, which has z₀ = 2.0 kpc and we also tabulate characteristic values in Table 3. In this figure, we show the original prediction for N_⋆ in Figure 5 along with the number of new NSs probed, ΔN_⋆, when h₀ is decreased by a factor of 2 in the high- or low-frequency regimes. We find that improving the high-frequency regime by a factor of 2 has a much larger effect on the number estimates compared to improving the low frequency by a factor of 2. This is due to the larger number of MS-pulsars from our catalog compared to those with f_GW ≲ 100 Hz and because one is sensitive to greater distances at higher frequencies. In this simple example, we see that lowering the value of h₀ by a factor of 2 in the high-frequency regime can add nearly three times as many new NSs as the current estimates for ≲ 10⁻⁶.

Zoom In Zoom Out Reset image size

Table 3. Estimates for Number of New NSs Probed with z₀ = 2.0 kpc When the Strain Amplitude h₀ Sensitivity in Different Frequency Regimes Is Increased by a Factor of 2

log10()	ΔN⋆(≤100 Hz)	ΔN⋆(≥1000 Hz)	ΔN3g
−5.00	390	0	4.1 × 106
−5.25	74	0	5.5 × 106
−5.50	14	0	6.4 × 106
−5.75	3	0	6.2 × 106
−6.00	0	920	4.2 × 106
−6.25	0	6.9 × 103	1.9 × 106
−6.50	0	2.2 × 103	5.6 × 105
−6.75	0	450	1.3 × 105
−7.00	0	84	2.8 × 104
−7.25	0	15	5.5 × 103
−7.50	0	3	1.0 × 103
−7.75	0	1	190
−8.00	0	0	35
−8.25	0	0	6
−8.50	0	0	1
−8.75	0	0	0
−9.00	0	0	0

Note. The rightmost column is the result of decreasing the strain at all frequencies by a factor of 10.

Download table as:  ASCIITypeset image

While improved sensitivity in the high-frequency regime will increase N_⋆, it is also worth examining the sensitivity of future third-generation GW detectors—for example, the Einstein Telescope (Punturo et al. 2010) and the Cosmic Explorer (Dwyer et al. 2015). This generation of detectors at present is estimated to be a factor of 10 times more sensitive than present detectors. To explore this possibility we revisit the model of N_⋆ with z₀ = 2 kpc, but with h₀ reduced by a factor of 10 at all frequencies.

The resulting improvement to N_⋆—which we define as ΔN_3g —is shown in Figure 6. We see that there is an increase in N_⋆ for all by at least an order of magnitude and for ∼ 10⁻⁷ the improvement is approximately 3 orders of magnitude. While for the original study we were well below the observational limit of ∼17 million NSs (probing ≲3%), with the sensitivity of third-generation GW detectors this limit is much closer to being reached (probing ≲40%), see Table 3. We note that this assumes the same search parameters as in Abbott et al. (2019a). An improved search could further improve these limits as well.

3.3. Alternative Searches

We present now an example of our methodology using new data for the strain sensitivity. Here, we follow the same process for both the determination of Φ(f) and N(d) as described in Section 2. For our strain sensitivity, we use the results of Steltner et al. (2021) for frequencies 20–500 Hz, Dergachev & Papa (2020) for frequencies 500–1700 Hz, and Dergachev & Papa (2021) for frequencies 1700–2000 Hz. We show this data in Figure 7. These latter two searches were intended to search for low ellipticity NSs. As such, the maximum spin-down allowed for a given ellipticity is considerably lower, $| \dot{f}| =2.5\times {10}^{-12}\,\mathrm{Hz}\ {{\rm{s}}}^{-1}$. The search performed by Steltner et al. (2021) did consider a higher spin-down, $| \dot{f}| =2.6\times {10}^{-9}\,\mathrm{Hz}\ {{\rm{s}}}^{-1}$. For this reason, we have now split up the integral in Equation (12) with one integral being ${f}_{1}=20\,\mathrm{Hz}\,\&\,{f}_{2}={f}_{\max }(| \dot{f}| =2.6\times {10}^{-9}\,\mathrm{Hz}\ {{\rm{s}}}^{-1})$ and the second being ${f}_{1}=500\,\mathrm{Hz}\,\&\,{f}_{2}={f}_{\max }(| \dot{f}| =2.5\,\times {10}^{-12}\,\mathrm{Hz}\ {{\rm{s}}}^{-1})$. This ensures that we do not hinder the breadth of the search of Steltner et al. (2021) with the lower $| \dot{f}|$ value.

Zoom In Zoom Out Reset image size

Indeed, our results show that using this improved data does yield more NSs probed at small , see Figure 8. This can largely be attributed to the overall decrease in h₀ at f_GW ≥ 500 Hz when compared to the strain from Abbott et al. (2019a). However, at ≳ 10⁻⁷ this data set probes significantly fewer NSs than for Abbott et al. (2019a). This is expected since NSs with these high were not prioritized in the study of h₀(f_GW) cited above. We note that the low-frequency regime has been considered by Dergachev & Alessandra Papa (2021) since this work was submitted.

Zoom In Zoom Out Reset image size

Prior searches for CGWs from known pulsars involve searching a well-defined number of NSs near an expected f_GW for each source (Abbott et al. 2017). Limits on in Figure 5 show the number of unknown NSs where GW (assuming a source with a given ) have been searched for and not found. We see that the models begin to result in similar estimates near N_⋆ ∼ 10⁶ above ≳ 10⁻⁶. This may be the maximum allowed by the NS's crust (Ushomirsky et al. 2000; Horowitz & Kadau 2009; Gittins et al. 2020), which is only slightly disfavored by our results. We predict that only ≳10⁵, or 0.1%, of Galactic NSs have been probed above = 10^−5.5. This places the limit that one in ten million NSs may have such an ellipticity or we would have detected a signal in gravitational waves.

The largest ellipticity in our tested range = 10⁻⁵, though heavily disfavored from studies of the breaking strain of an NS's crust, cannot be ruled out entirely using current CGW data. From our results, we only rule out this ellipticity for ≈1.6% of all Galactic NSs. This may be somewhat unrealistic, for example, a millisecond pulsar would produce a very large strain amplitude in CGW signal with such a large ellipticity.

The theoretical upper limit on of ∼few × 10⁻⁶ has likewise not been ruled out. In fact, our results suggest that ≲0.1% of all Galactic NSs have been probed at this ellipticity for all values of the disk thickness. Therefore there is great need to continue searching for CGWs arising from NSs with ellipticities near this value. With further studies of NS ellipticities and CGW searches, this limit may become more apparent.

Our methodology used in this study attempted to keep things as simple as possible. Several additional complications to the study could be introduced to further constrain N_⋆. First, our choice of Equation (3) as the distribution of Galactic NSs is a simple model that largely follows the star formation pattern in the Galactic disk. In reality, NSs may have a much different distribution, in part resulting from large transverse velocity kicks during their birth. We have attempted to mitigate this by picking different values for z₀, which either condense (z₀ = 0.1 kpc) or expand (z₀ = 4.0 kpc) the density distribution of NSs as seen from Earth.

Most NSs are expected to be born with a transverse space velocity from a supernovae kick (Shklovskii 1970). Analyses of NS orbits suggest that fewer than ≲20% are retained in the disk and a greater fraction remain in bound orbits in the Galactic halo (Sartore et al. 2010). Furthermore, some NSs have a sufficient space-velocity to escape the Galactic potential entirely (Arzoumanian et al. 2002; Katsuda et al. 2018; Nakamura et al. 2019). As a result, kicked NSs leaving the disk will spread the density distribution in Equation (4) to larger z₀ than is typical for other stellar populations. Additionally, our calculated estimate for the total number of NSs probed in Table 2 could be reduced by more than a factor of 2 depending on the real distribution of supernovae kick velocities. A clear next step with this type of estimate would be to self-consistently include an empirical density distribution of NSs that can account for a kicked NS population.

Additionally, we have chosen to neglect CGWs arising from NSs in binaries because of the added complication it would cause on the GW signal and on the search parameters. However, in future work, it would be useful if the GW search treated binary NSs and isolated NSs separately. The newly developed BinarySkyHough (Covas & Sintes 2019) pipeline is much better equipped to search for CGWs in binaries than its predecessor SkyHough (Krishnan et al. 2004), which was used in Abbott et al. (2019a). By better constraining the values of , this can also further improve the search parameter computation time.

We find that the disk thickness parameter z₀ from Equation (4) has a significant impact on the estimated number of nearby NS. These nearby sources of CGWs would be vital in constraining ellipticities ≲10⁻⁷. In the thin disk approximation (${z}_{0}\longrightarrow 0$), Equation (3) behaves like ρ(d) ≈ d for small values of d. However, for other values of z₀, where d ≪ z₀, the distribution is instead ρ(d) ≈ d². This has significant effects on nearby number estimates as any stars lying above the plane of the disk are then condensed, thereby increasing the total number of stars estimated. This is easily seen in the right-hand panel of Figure 3, and the effects on the estimated numbers of NSs are seen in Figure 5 and in Table 2. Figure 5 is the result of Equation (12) for the values of z₀ used in Table 1. We see that N_⋆ is very sensitive to z₀ for small values of , decreasing for increasing z₀. Better determination of the disk thickness of the Galactic NS population is important for constraining ≲ 10^−5.5, where an order of magnitude difference exists between the models used here.

We explore the implications of a new CGW search with improved h₀(f_GW) sensitivity using the current generation of GW detectors. Our results show that improving sensitivity in the high-frequency regime (f_GW ≥ 1000 Hz) can have the greatest impact on the search for CGWs. From Figure 6, we can see that the improvements in h₀ sensitivity can have much higher returns on the total number of new NSs probed. At ≲ 10⁻⁶, for example, this results in finding approximately three times N_⋆ new NSs.

Interestingly, this is not true for very high values of . We see that the high-frequency regime has a turnoff point at ∼ 10⁻⁶ which occurs for two reasons. First, the original search already probed a significant fraction of visible NSs in the disk for > 10⁻⁶, and so fewer new NSs would become visible. Second, for > 10^−5.5, ${f}_{\max }\lt 1000\,\mathrm{Hz}$ and so the improvement is no longer limiting the contribution of h₀(f_GW) to the integral in Equation (12). From these two points, we see that the strain sensitivity at high frequency has a significant impact on searches for NSs with moderate ellipticity. Conversely, improving the low-frequency regime (f_GW ≤ 100 Hz) certainly increases the number probed, it is approximately 4 orders of magnitude smaller in effect than improving the high-frequency regime for moderate ellipticity. This is because there is a much larger number of MS-pulsars with spin frequencies in excess of ν > 50 Hz, as discussed in Section 2.3.

Third-generation detectors may dramatically increase the number of NSs probed. Given the blue dashed–dotted curve in Figure 6 we can see that improving h₀ by a factor of 10 increases N_⋆ by more than a factor of 100–1000 times. We note that this assumes for Equation (13) $| \dot{f}| =2\times {10}^{-9}\,\mathrm{Hz}\,\ {{\rm{s}}}^{-1}$, which may not be the true limit considered when searches with these instruments take place. Despite this, however, just the improvements to h₀ we estimate will probe almost 40% of all the NSs in the Galaxy at large . In addition, improvements in search techniques and computer resources may further increase the number of NSs probed.

We conclude our discussion with the analysis of Section 3.3. The data used here is comprised of several additional analyses of the data from Abbott et al. (2019a), however, now with improved strain sensitivity. Both searches use different techniques and have different goals for performing their respective searches. For example, Dergachev & Papa (2020) were primarily interested in finding low NS, which allows for a smaller maximum value for $| \dot{f}|$. If one sets f_max = 2000 Hz, then for = 10⁻⁸ $| \dot{f}| \approx 5.55\times {10}^{-12}\,\mathrm{Hz}\,\ {{\rm{s}}}^{-1}$. This allows for more restrictive limits on h₀, as seen in Figure 7. While the value of $| \dot{f}| =2.5\times {10}^{-12}\,\mathrm{Hz}\,\ {{\rm{s}}}^{-1}$ is consistent with pulsar data, this constraint on the analysis has two effects on the overall results. First, it reduces the search parameter space considerably and therefore allows for a better determination of h₀(f_GW), as seen in Figure 7. As we have shown in Section 3.2, reducing the strain amplitude does increase N_⋆.

However, the second, and most important effect for this work, is that it limits the amount of detectable NSs. For example, taking $| \dot{f}| =2\times {10}^{-9}\,\mathrm{Hz}\ {{\rm{s}}}^{-1}$ as we did in Section 3.1, for = 10⁻⁶ this means ${f}_{\max }\,\approx \,1029\,\mathrm{Hz}$. Using $| \dot{f}| =2.5\times {10}^{-12}\,\mathrm{Hz}\ {{\rm{s}}}^{-1}$ instead, ${f}_{\max }\,\approx \,220\,\mathrm{Hz}$. This means that this data set may be inefficient when looking for highly elliptical NSs because none of the MSPs are being probed. Interestingly though, the improved data set does probe more NSs at > 10^−5.75. This is because the values of f_max for both sets exclude the highest frequencies from the search, meaning only the low-frequency sources contribute to N_⋆. Since we see in Figure 7 that the strain used in this analysis is smaller than from Abbott et al. (2019a), slightly more NSs are probed.

Note that there are two possible approaches when selecting the optimal choice for $| \dot{f}|$ because the ellipticity distribution of NSs is unknown. Should a future CGW search occur with the intent of probing the highest ellipticities near ∼ 10⁻⁶ one should consider using a higher $| \dot{f}|$ limit. Possibly the most promising value of $| \dot{f}|$ is slightly higher than considered by Abbott et al. (2019a), 5.55 × 10⁻⁸ Hz s⁻¹. This value will ensure that ${f}_{\max }=2000\,\mathrm{Hz}$ for = 10⁻⁶ so that any fast-spinning NSs will not be excluded. On the other hand, if the intent is to find lower ellipticity NSs—for instance, near ∼ 10⁻⁹—one should consider a deeper search with a lower $| \dot{f}|$ limit.

We have detailed estimates on the total number of NSs probed with GW detectors. In doing so, we have shown that continuous GW searches suggest that fewer than about 1 in 10,000 NSs have an ellipticity ≳10⁻⁶. Additionally, we have shown that the disk thickness strongly affects the number counts of nearby neutron stars while leaving more distant stars largely unaffected. We have explored the effects of improving strain amplitude sensitivity at higher frequencies, which can increase the amount of NSs probed to a given ellipticity. These estimates are important for setting upper limits on the ellipticity of an NS as well as detecting radio-quiet NS that may be nearby, yet unobserved. Finally, we discuss the impact of third-generation detectors and find that they may probe 100–1000 times more NSs than have presently been probed.

We would like to thank the anonymous referee for helpful comments and kind words. We would also like to extend our gratitude to M. Papa for useful comments and for providing the data used in Section 3.3. We also thank N. Andersson, F. Gittins, V. Dergachev, F. De Lillo, and P. Covas for their helpful comments.

This material is based upon work supported by the U.S. Department of Energy Office of Science, Office of Nuclear Physics under Awards DE-FG02-87ER40365 (Indiana University) and Number DE-SC0008808 (NUCLEI SciDAC Collaboration).