Population synthesis and parameter estimation of neutron stars with continuous gravitational waves and third-generation detectors

Precise measurement of stellar properties through the observation of continuous gravitational waves from spinning non-axisymmetric neutron stars can shed light onto new physics beyond terrestrial laboratories. Although hitherto undetected, prospects for detecting continuous gravitational waves improve with longer observation periods and more sensitive gravitational wave detectors. We study the capability of the Advanced Laser Interferometer Gravitational-Wave Observatory, and the Einstein Telescope to measure the physical properties of neutron stars through continuous gravitational wave observations. We simulate a population of Galactic neutron stars, assume continuous gravitational waves from the stars have been detected, and perform parameter estimation of the detected signals. Using the estimated parameters, we infer the stars' moments of inertia, ellipticities, and the components of the magnetic dipole moment perpendicular to the rotation axis. The estimation of the braking index proved challenging and is responsible for the majority of the uncertainties in the inferred parameters. Using the Einstein Telescope with an observation period of 5 yrs, point estimates using median can be made on the moments of inertia with error of ~ 10 - 100% and on the ellipticities with error of ~ 5 - 50%, subject to the inference of the braking index. The perpendicular magnetic dipole moment could not be accurately inferred for neutron stars that emit mainly gravitational waves.


INTRODUCTION
One hundred years after Einstein (1915) published his paper on the general theory of relativity, the first gravitational wave detection (the binary black hole merger GW150914, Abbott et al. 2016) was made on September 14, 2015 by the Advanced Laser Interferometer Gravitational-Wave Observatory (LIGO, Aasi et al. 2015), a major milestone for gravitational wave astronomy.The observation of the binary neutron star merger GW170817 by LIGO and Virgo (Acernese et al. 2014) in both gravitational and electromagnetic radiation (Abbott et al. 2017a,b) is an example of multi-messenger astronomy and offers hope to further our understanding of neutron stars, which are known for their extreme densities and currently poorly understood physics.
Gravitational waves from mergers of binary neutron stars only convey information at the end of their life cycle, when they are under strong mutual gravitational perturbation.Alternatively, a neutron star with asymmetries about its axis of rotation can generate continuous gravitational waves (henceforth "continuous waves") that are longlived and quasi-monochromatic (Zimmermann & Szedenits 1979;Bonazzola & Gourgoulhon 1996;Riles 2017;Sieniawska & Bejger 2019).As continuous waves are generated by neutron stars in their equilibrium state, they could provide complementary insights into the stars' physical properties not conveyed through binary neutron star mergers.Continuous waves are weaker than gravitational waves from mergers, however, and have yet to be detected by the current second generation detectors LIGO and Virgo (Piccinni 2022;Riles 2023).
Efforts to detect continuous waves from neutron stars include development of, and improvements to, data analysis techniques for various types of searches (Tenorio et al. 2021;Wette 2023).These are typically specialised to a variety of sources, including but not limited to targeted searches and narrow-band searches for known pulsars (Abbott et al. 2022b;Zhang et al. 2021), directed searches for supernova remnants (Abbott et al. 2021b;Lindblom & Owen 2020), and all-sky surveys for undiscovered neutron stars (Abbott et al. 2019;Abbott et al. 2022a;Covas et al. 2022).Directly increasing the sensitivities of gravitational wave detectors by further suppressing the noise sources also improves the prospect of a first detection.Current upgrades to the second generation detectors are anticipated to reduce the noise amplitude spectral density by a factor of two (Miller et al. 2015).
The planned third generation detectors are expected to achieve a tenfold increase in strain sensitivity in a wide frequency range, offering hope to the detection of gravitational waves from other sources beyond binary compact object mergers (Sathyaprakash et al. 2012;Bailes et al. 2021;Hall 2022).Proposed third generation groundbased detector concepts include the Einstein Telescope (ET) and the Cosmic Explorer, both of which are planned to begin construction in the 2030s.The Einstein Telescope is expected to have three arms of 10 km length, arranged in an equilateral triangle formation and located underground.It will have three detectors at the vertices of the triangle and two interferometers at each detector, forming a total of six interferometers (Punturo et al. 2010;Sathyaprakash et al. 2012).Cosmic Explorer is planned to be constructed on the surface with a design similar to that of LIGO, but with arm length increased tenfold to 40 km (Reitze et al. 2019;Evans et al. 2021).
In this work we present a study of the potential for observations of continuous waves from isolated neutron stars by LIGO and the third generation detector Einstein Telescope, and the estimation of the continuous wave signal parameters.In addition, we also estimate the errors from the inference of neutron star properties, including the principal moment of inertia, ellipticity, and the component of magnetic dipole moment perpendicular to the rotation axis, using a theoretical framework developed by Lu et al. (2023).After reviewing relevant background in Section 2, we synthesise a population of neutron stars emitting continuous waves and electromagnetic radiation in Section 3, and infer the stellar properties of ten neutron stars with the largest continuous wave characteristic strain amplitude using Bayesian inference in Section 4. We conclude with a discussion in Section 5.

BACKGROUND
Currently, neutron stars are observed primarily as pulsars.A pulsar emits electromagnetic radiation, typically directed along its magnetic dipole axis as a result of a strong stellar magnetic field (∼ 10 12 G).The open field lines of the magnetic field are strong enough to accelerate charges to a relativistic speed sufficient to escape the magnetosphere (Kramer 2005).If the dipole axis of the magnetic field is not aligned with the axis of rotation, the accelerated charges produce regular pulses of electromagnetic radiation at the same frequency as the rotation of the pulsar.A pulsar may be thought of as an extraterrestrial "lighthouse" that emits radiation at a fixed location for an observer.Pulsar evolution is well described using its period  and spindown , as both its characteristic age and magnetic field strength depend only on these two parameters (Kramer 2005).(The braking index  is defined and explained below in Eq. ( 5).) In addition to electromagnetic radiation, neutron stars may also emit continuous waves through several mechanisms, such as deformation away from axi-symmetry due to their magnetic fields (Zimmermann & Szedenits 1979;Bonazzola & Gourgoulhon 1996), quasi-normal fluid perturbations known as -modes (Owen et al. 1998;Andersson 1998;Friedman & Morsink 1998), or accretion of matter from a binary companion star (Bildsten 1998;Watts et al. 2008).Deformation of neutron stars about their rotation axes can be caused by cooling and cracking of the crust (Pandharipande et al. 1976), a non-axisymmetric magnetic field (Zimmermann & Szedenits 1979), magnetically-confined mountains (Melatos & Payne 2005), or electron capture gradients (Ushomirsky et al. 2000).The characteristic strain amplitude of a continuous wave signal from a deformed neutron star is given by (Jaranowski et al. 1998) where  is the gravitational constant;  is the speed of light;   is the principal moment of inertia aligned with the rotation axis;  is the equatorial ellipticity that characterises the extent of the neutron star's deformation;  is the distance to the detector; and  = 2 is the gravitational wave frequency, taken to be twice the rotational frequency .
As a neutron star rotates, it may emit energy through electromagnetic and gravitational wave radiation, which extracts rotational kinetic energy and in turn reduces its spin frequency.The spindown of a neutron star can be characterised as (Manchester et al. 1985) where  is a constant, and the braking index  is obtained by differentiating and rearranging Eq. ( 4): (5) It is expected from theory that neutron stars have a braking index ranging between 3 and 5, with 3 being the case when the neutron star emits only electromagnetic radiation (Ostriker & Gunn 1969), and 5 being the case for only continuous wave emission through timevarying mass quadrupole.Gravitational waves can additionally be emitted through a current-type quadrupole moment due to r-modes, which results in a breaking index of 7.There are measured braking indices from pulsars that span orders of magnitude outside this range (Johnston & Galloway 1999;Zhang & Xie 2012;Lower et al. 2021), however this likely reflects the difficulty of accurately measuring the second derivative  of the rotational frequency, or the impact of glitches (Ho et al. 2020) or timing noise (Vargas & Melatos 2023;Hobbs et al. 2004).We adopt the commonly accepted range of 3 ≤  ≤ 5.The amplitude of a continuous wave signal observed by a detector from a neutron star will be modulated by the detector's antenna pattern.We henceforth assume the neutron stars to be biaxial rotors with continuous wave frequency of  = 2 and amplitude given by (Jaranowski et al. 1998): where ℎ 0 is the characteristic strain amplitude; () is the phase of the signal with initial value  0 ; and cos  is the orientation angle, describing whether the neutron star is viewed face-on (spin-axis along the line of sight, cos  = ±1) or edge-on (spin-axis perpendicular to the line of sight, cos  = 0).The  +,× are the detector response functions for + or × polarised gravitational waves respectively, and are parameterised by the polarisation angle .The parameters ì A ≡ (ℎ 0 , cos , ,  0 ) are commonly referred to as the amplitude parameters.
The continuous wave signal of an isolated neutron star is expected to be nearly monochromatic in the neutron star reference frame; with respect to the detector reference frame, however, it is Doppler shifted by the rotation and the orbital motion of the Earth.On the time-scale of a day, the signal is modulated by the rotation of the Earth; while on the time-scale of a year, the signal is modulated by the orbit of the Earth around the Solar System Barycentre (SSB).(The velocity of the SSB relative to the neutron star could further affect the detected frequency, although this effect is approximately constant on the timescale of measurements.)In the neutron star frame of reference, the continuous wave signal phase is modelled, up to the second-order time derivative, as (Jaranowski et al. 1998) where  is time in the SSB frame; and  ,  ,  are the continuous wave frequency and its first-and second-order time derivatives, or spindowns, at a reference time  = 0.The neutron star time  can be converted to time  in the detector's frame of reference by where ì  is the position of the detector with respect to the SSB, and n is the unit vector pointing from the SSB to the neutron star.(Here we ignore relativistic effects, and assume the neutron star is at rest with respect to the SSB.)The parameters ì  ≡  ,  ,  , n are typically referred to as the phase evolution parameters.

NEUTRON STAR POPULATION SYNTHESIS
Since no continuous waves have been detected from the Galactic neutron star population, a simulated neutron star population must realistically contain stars that emit electromagnetic and/or continuous gravitational waves.Such a population has been previously investigated by a number of authors (Palomba 2005;Knispel & Allen 2008;Wade et al. 2012;Pitkin 2011;Woan et al. 2018;Cieślar et al. 2021;Reed et al. 2021).In this work we utilise Monte-Carlo methods to simulate a neutron star population.We assign neutron star parameters drawn from theoretical probability distributions, evolve the neutron stars temporally, and use their final spin frequency to simulate continuous wave signals.

Simulation method and parameters
The initial spin frequency  0 can be obtained through the inverse of the initial period  0 .Three models for the initial period distribution have been proposed (Palomba 2005;Knispel & Allen 2008).The first model is a lognormal distribution with P0 = 5 ms,  = 0.69, and excludes all  0 < 0.5 ms.The second model uses the same lognormal distribution but with all  0 < 10 ms set to 10 ms; this mimics the potential existence of -modes in young neutron stars, which increases the spin period to 10 ms1 once -modes become saturated.The third model simply uses a uniform distribution of  0 ∈ [2, 15] ms; this is intended to accommodate both -modes, and the fall-back of matter after the supernova which decreases the period through additional angular momentum (Watts & Andersson 2002).A uniform initial period distribution was favoured by simulations of Galactic binary neutron stars performed in Sgalletta et al. (2023).We follow Palomba (2005) and Knispel & Allen (2008) and adopt the third model for  0 .
We assume an average birth rate of one Galactic neutron star every 100 years, as evidenced by the observation of core-collapse supernovae in the Galaxy (Diehl et al. 2006;Rozwadowska et al. 2021).An array of random variables drawn from a Poisson distribution with rate  = 100 yrs is generated and cumulatively summed, and the values of the summed array are taken as the ages of the simulated stars.Similar to Palomba (2005), an upper bound of  max = 10 8 yrs is used, as a compromise between the computer memory required to store a simulated population and realistic expectations of Galactic neutron star ages.
Instead of evolving neutron stars through the Galactic potential, we populate the Galaxy with stationary neutron stars according to a theoretical spatial distribution, similar to the approach used by Reed et al. (2021).Current understanding suggests that the spatial distribution of Galactic neutron stars follows a Gaussian distribution in the radial direction and a double-sided exponential (Laplace) distribution in the vertical direction (Faucher-Giguère & Loeb 2010;Binney & Tremaine 2008).In the Galacto-centric cylindrical frame (, , ), the distribution for the radial distance  is given by p with  > 0 and  = 5 kpc.The vertical distribution is given by with  0 ∈ [0.5, 1] kpc for conventional stars (Binney & Tremaine 2008).In their work, Reed et al. (2021), however, chose  0 = 0.1, 2, 4 kpc to investigate the effect of supernova kicks on the distribution of neutron stars.We use  0 = 2 kpc as a compromise between a clustered vertical distribution ( 0 = 0.1 kpc) and a spread-out distribution ( 0 = 4 kpc).We use a uniform distribution of angle  between each neutron star and the Sun with respect to the Galactic centre, with  ≡ 0 for the Sun.With the location of the Earth assumed to be coincident with the Sun ( e = 8.25 kpc,  e = 0,  e = 0.02 kpc; Reed et al. 2021;Humphreys & Larsen 1995), the distance  between a neutron star at (, , ) and the Earth can then be calculated.One can determine the spindown  of a neutron star from the conservation of energy, as follows.Assuming the total kinetic energy budget from the star's rotation is expended either in electromagnetic or continuous wave emission gives (Lu et al. 2023) Due to the expected scale of the ellipticity  ≪ 1, neutron stars can be assumed to be nearly spherical, and hence the rotational energy can be approximated as that of a rotating sphere (Wette et al. 2008): The energy emitted by a rotating magnetic dipole is (Ostriker & Gunn 1969) where  0 is the vacuum permeability and  is the radius of the neutron star.The component of the magnetic dipole moment perpendicular to the rotational axis   can be related to the magnetic field strength at the surface  by (Condon & Ransom 2016) 2 The simulation assumes, for simplicity, that the magnetic dipole moment is entirely perpendicular to the rotation axis, i.e. sin  = 1.The power emitted in continuous waves is given by (Ostriker & Gunn 1969;Riles 2017) Substituting Eqs. ( 13) and ( 15) into Eq.( 12) and rearranging yields the equation for the spindown : Wade et al. ( 2012) solve Eq. ( 16) and obtain an analytical expression for the neutron star age  as a function of its final and initial spin frequencies,  and  0 respectively: where and  =  GW / EM are used for brevity.
To obtain the spin frequency evolution () of a neutron star, we need to invert Eq. ( 17).One approach is to numerically solve Eq. ( 17) through root-finding; this method suffers from poor convergence, however, due to the rapid decrease of  at small .Alternatively, one evaluates Eq. ( 17) at an array of sample frequencies {  } to obtain the corresponding time array {  }, then interpolate to obtain a continuous function of ().This approach, however, is limited by the range of sample frequencies  lower ,  upper , which corresponds to a range for {  } of [ lower =  ( upper ),  upper =  ( lower )].The upper bound for {  } can easily be set as  upper =  0 , giving a lower bound  lower =  ( upper ) = 0.The lower bound for {  } may, however, be higher than the true value, i.e.  upper =  ( lower ) <  age , and hence interpolation cannot be performed.This occurs for very old stars with large  age and very small .While one can simply lower  lower , doing so requires significantly more points in {  } (assuming linear sampling) and is generally inefficient.
We therefore employ a hybrid approach of Piece-wise Cubic Hermite Interpolation (PCHIP) at small  and root-finding using the Levenberg-Marquardt algorithm at large  to solve Eq. ( 17). Figure 1 illustrates this process.Interpolation is tractable when ( age ) ∈ [ 0 , 0.005 0 ], which is the case for most neutron stars.For older neutron stars with ( age ) < 0.005 0 , root-finding is then used, starting with the initial guess  =  lower ,  = 0.005 0 .The true age is plotted as the green line.

Synthesised population
method, as described in Section 3.1, and can be understood through Fig. 2b.The population may be partitioned into four regions in the - space.The purple region in Fig. 2b contains mainly young neutron stars; they do not have the time to spin down significantly, and hence have higher final spin frequencies (smaller periods) and higher period derivatives than older stars.This positions them in the upper left region of the - space.Conversely, older stars have more time to spin down, so they have higher periods and lower period derivatives.At a constant birth rate, there are more older neutron stars than young ones, e.g. 10 times more neutron stars with age 10 7 yrs than those with age 10 6 yrs.This results in a lack of simulated neutron stars in the purple region and a high concentration of older neutron stars.
For a neutron star with negligible ellipticity, its spindown is dominated by electromagnetic radiation, and it follows a trajectory similar to path (1) in Fig. 2b.It has a slower spindown of  ∝ −3 , corresponding to  ∝  −1 , and will eventually occupy the green region in Fig. 2b.A neutron star with both significant ellipticity and a strong magnetic field initially radiates energy mostly through continuous waves, then switches to mostly electromagnetic radiation as the spin frequency  becomes smaller.Its braking index evolves from  ∼ 5 to  ∼ 3, resulting in a non-linear trajectory exemplified by path (2).
The energy emission of a neutron star with significant ellipticity and low magnetic field strength is dominated by continuous wave radiation.Its trajectory is therefore steeper, with spindown  ∝ − 5 , or  ∝  −3 , as shown by path (3) in Fig. 2b.It can reach a slower  than stars dominated by electromagnetic radiation with the same characteristic age.This results in the discontinuity in the line of maximum  age seen in Fig. 2b.This can also be seen from Eq. (1), which yields different characteristic ages for stars with different braking indices.Given their steeper gradient in the - plane, these neutron stars are typically located in the red region.
For a neutron star with both insignificant ellipticity and weak magnetic field strength, its spindown will be small, and there will be no significant difference between its initial and final spin frequencies, resulting in a short or even unobservable trajectory demonstrated by path (4) and the stationary star next to it in Fig. 2b.The blue region typically occupied by these neutron stars is defined by the minimum and maximum initial period, and the minimum ellipticity.
To be precise, in Fig. 2b, the upper bound of the population (green  solid and dashed lines) is defined by the maximum magnetic field strength  max through Eq. ( 2); the right boundaries of the population (two purple solid lines) are defined by the maximum age through Eq. ( 1) and the maximum ellipticity; the two solid blue lines are defined by the minimum and maximum initial period; and the lower bound of the population (solid red line) is defined by the minimum ellipticity through Eq. ( 16).Inside the population envelope, the red dashed line used to bound the blue region is defined by the intersection of the blue line  0,max and  max ; the intersection of that line with  0,min is then used to define the purple region; and the red and green regions are then separated by the red dashed line defined by  max .While these four regions provide a basic intuition regarding the distribution of the population, they do not map one-to-one to all neutron stars.For example, the majority of the continuous wavedominated (red) region, except near the boundary of  max , can still contain neutron stars with predominantly electromagnetic radiation.
The Australian Telescope Network Facility (ATNF) hosts a catalogue of observed pulsars (Manchester et al. 2005), which are superimposed on the synthesised population in Fig. 2a.The majority of the pulsars are clustered at the centre of the diagram, occupying the electromagnetic radiation-dominated (green) region described in Fig. 2b.In addition, the majority of neutron stars are observed at low frequencies  = 2 = 2/ ≲ 10 Hz, where present-day gravitational wave detectors have limited sensitivity.This is consistent with the current lack of detection of continuous waves through targeted searches of existing pulsars (Abbott et al. 2022b,c).Note that, while most of the observed pulsars coincide with the simulated population, some pulsars have characteristic ages > 10 8 yrs.This is a direct consequence of  max = 10 8 yrs set for the simulation.Regardless, such discrepancies are not likely to drastically affect the usefulness of the simulation in respect of the detectability of continuous waves from young neutron stars.
A subset of pulsars is located in the lower left region of Fig. 2a, with  < 0.01 s.These pulsars are millisecond pulsars and are typically in binary systems.Even though they coincide with the low  region of the simulated population, the underlying physics of the two are different as we did not consider spin-ups due to accretion in our simulation.However, as argued in Wade et al. (2012), while we did not explicitly account for millisecond pulsars, we did not exclude them either.An old neutron star that has spun-up through recycling can be thought of as a young neutron star born with a high spin frequency.
Similar to Fig. 2a, Cieślar et al. (2021) plotted simulated neutron stars with continuous waves detectable by the Einstein Telescope in their Fig. 7. Examining the  axis, we notice differences between the two sets of results.The detectable neutron stars simulated by Cieślar et al. have higher , younger characteristic ages and therefore fall in the young (purple) region in Fig. 2b.This may be due to their assumption of a decay model for the ellipticity , which means that for older neutron stars their ellipticities, and consequently ℎ 0 , are lower.In their model, therefore, only young neutron stars are detectable.

CONTINUOUS WAVE PARAMETER ESTIMATION
We now investigate using continuous waves to infer properties of neutron stars in the simulated population.Continuous waves may be detected in various ways, such as through a targeted search of known pulsars, or from a blind all-sky search.In this work we assume that a continuous wave signal has already been detected and sufficiently localised to allow a targeted search with maximum sensitivity. 3revious work has investigated the detectability of continuous waves from a population of neutron stars.Wade et al. (2012) compared the continuous wave amplitude with the estimated noise curve of the detectors; a neutron star is considered detected if its strain amplitude is above the noise curve.Wade et al. also derived a theoretical framework to describe the detectability of continuous waves from neutron stars, which holds for young neutron stars aged ≲ 10 7 years.Cieślar et al. (2021) took into account the spatial position of the neutron star relative to the detector, calculated a signal-to-noise (/) ratio, and assumed detection only for stars with / > 11.4 (based on Abbott et al. 2004).Reed et al. (2021) used the results from existing continuous wave searches (Abbott et al. 2019;Steltner et al. 2021;Dergachev & Papa 2020, 2021) to obtain a function of strain amplitude and frequency, which is then applied to the simulated population and the ATNF catalogue to determine detectability.Pitkin (2011) used Markov Chain Monte Carlo (MCMC) methods to obtain the posterior distributions of signal amplitude parameters, and constraints on the possible range of gravitational quadrupole moment  22 and magnetic field strength , assuming LIGO, Virgo, and Einstein Telescope data.We simulate continuous wave signals from the synthesised population and use Bayesian inference to estimate the continuous wave signal parameters.From the posteriors of the parameters we then infer the physical properties of the neutron stars and estimate the errors from the inference.

Bayesian inference of simulated signals
From the simulated neutron star population, we select the 10 neutron stars with the largest ℎ 0 for the parameter estimation study, which are summarised in Table 1.We note that the first three neutron stars in our population have ℎ 0 that exceeds the best upper limit of ℎ 0 ∼ 1.1 × 10 −25 ruled out by existing targeted searches using O3 data (Abbott et al. 2022a).This could happen as we did not incorporate such restrictions during the simulation.However, the continuous wave frequency range in which the best upper limit is obtained is  = 100 − 200 Hz.None of the top three stars have frequencies in this range, so there is no strong contradiction with the result of existing searches.
We utilise the CWInPy (CW Inference in Python, Pitkin 2022) package and perform software injections to simulate continuous wave signals.Currently, the only third generation detector supported by CWInPy is the D configuration of the Einstein Telescope (ET-D).Hence, we limit the scope of this work to LIGO and ET-D.We assume Gaussian noise with a standard deviation defined by the amplitude spectral density specific to each detector at the continuous wave signal frequency  = 2.Comparing the sensitivities of LIGO (Kissel 2021;Betzwieser 2021) and ET-D (Hild et al. 2011;Evans et al. 2016), the latter will have noise amplitude spectral density ∼ 10 times smaller than that of LIGO, resulting in higher signal-to-noise ratios.
For LIGO, we simulate detector data at times identical to the real data recorded by each of the two LIGO detectors (Hanford and Livingston) across all observation runs, O1-O3 (Abbott et al. 2021a(Abbott et al. , 2023)).For the Einstein Telescope, we perform two separate analyses with gap-less detector data spanning 2 years and 5 years respectively.The 2-year analysis is approximately the combined effective observation time of LIGO from O1 to O3b, while the 5-year analysis is similar to the total amount of time spanned by these observation runs (2015-2020), during which the continuous wave signal can evolve.
The amplitude parameters ì A and frequency parameters ì  of the simulated neutron stars are used to generate continuous wave signals in CWInPy.Each simulated signal is then heterodyned and downsampled to a frequency of 1/60 Hz (Dupuis & Woan 2005).We perform parameter estimation of both amplitude and frequency parameters, namely ℎ 0 , cos ,  0 , ,  ,  ,  .To obtain the posterior of a parameter, one requires the likelihood, the evidence, and the prior.The likelihood comes from the Gaussian noise assumption, and we use nested sampling (Skilling 2004, 2006, see Ashton et al. 2022 for an intuitive illustration) to evaluate the evidence and obtain the posterior.The nested sampling scheme is implemented in CWInPy via the Bayesian Inference Library (BILBY, Ashton et al. 2019) which uses dynesty (Speagle 2020;Higson et al. 2019) for nested sampling.
The choice of priors represents the initial assumptions made for the parameters, which in turn affects the performance of the sampler and consequently the estimated posteriors.The prior on ℎ 0 is chosen to be a uniform distribution [0, 10 −24 ] to include any likely continuous wave strain amplitude; in the simulated population, ℎ 0 ≲ 2 × 10 −25 .The value of cos  is unknown a priori and is thus chosen uniformly over [−1, 1].The prior on  0 is chosen only between [0, ], as in CWInPy  0 is the initial rotational phase offset, and will cover the full phase range when converted to a gravitational wave frequency via  = 2.The prior on  is set to [0, /2] to account for the degeneracy in  given by the transformation of  →  + /2 (Jones 2015).
For targeted searches, the phase parameters  = 2,  ,  are assumed to be roughly determined based on previous observations, e.g. from electromagnetic detection of a known pulsar, or from follow-up of a continuous wave candidate found in an all-sky search.To ensure a physical braking index, we substitute  for  using Eq. ( 5) and impose the constraint 3 ≤  ≤ 5. We use uniform prior distributions on  ,  ,  centred at the true values with single-side widths 3, where  is obtained from the inverse Fisher information matrix (Jaranowski & Królak 2010;Lu et al. 2023).The priors are chosen to be broad enough to not restrict the posteriors produced by the Bayesian sampler, but narrow enough to allow convergence of the posteriors.Explicitly,  for each parameter is given by: Here, is the sensitivity depth, which is related to the signal-to-noise ratio  2 by (Behnke et al. 2015;Dreissigacker et al. 2018) is the observation time; and  ℎ is the power spectral density of the strain noise in the detector.
Figure 3 shows a corner plot of the prior and posterior distributions of ℎ 0 ,  ,  ,  for continuous waves from a simulated neutron star in ET-D data with an observation time of 2 years.The 1D distributions for each parameter are shown along the diagonal; the off-diagonal plots show the 2D distributions over pairs of parameters.We plot relative errors E defined by i.e. the true values of each parameter are located at E = 0.The convergence of the posteriors is consistent with our prior assumption of the detectability of continuous waves, in particular by the Einstein Telescope.Figure 3 also highlights the difficulty in accurately estimating ; we will see in Section 4.2 that this limitation dominates the uncertainty in inferring the stellar properties.

Inferring physical properties
Starting with Eq. ( 16), Lu et al. (2023) developed a theoretical framework where the continuous wave signal parameters ℎ 0 ,  ,  ,  and the distance to the neutron star  can be used to infer the physical properties of the star: the principal moment of inertia   , the

3.
Corner plot of the (blue) and posterior (red) distributions of parameter estimation results for an example simulated neutron star (labelled J1704-0203).The ET-D detector with observation time of  = 2 yrs is used.
Relative errors E of the posteriors are shown, with true values at E = 0.The vertical dashed lines in the 1D histograms represent 1- confidence, and the three contours in the 2D histograms represent 3-.Note that the distributions of E (ℎ 0 ) is zoomed in to make the posterior distribution visible.
ellipticity , and the perpendicular magnetic dipole moment   : where  EM = 2 2 /3 3  0 and  GW = 32 4 /5 5 .We assume a 20% error in the measurement of  (i.e. ∼  truth (1 + 0.2) where  is drawn from the standard normal distribution) which is expected to be achievable by the upcoming radio telescope Square Kilometre Array (Smits et al. 2011).
Figure 4 shows example corner plots of the inferred physical properties, with the posterior distribution of the braking index  positioned at the top right.Comparing Fig. 4a and Fig. 4b, there is little difference between the posteriors inferred using LIGO versus ET-D with  = 2 yrs of data.As seen in Fig. 3,  and  are very well constrained, whereas  ∝  is not: the posteriors of  are uniform over the range 3 − 5. (In the 1D histogram of ,  truth = 5, and E ranges from E ( = 3) = −0.4 to E ( = 5) = 0.0).This means that the majority of the uncertainty in the inferred parameters stems from the uncertainty in  .Since  in both Figs.4a and 4b is badly estimated, the uncertainties in   , ,   are similar for both LIGO and ET-D at  = 2 yrs.With an observation time of  = 5 yrs, however (Fig. 4c),   [Eq.( 21)] is small enough so that the prior on  is a narrow subset of [3,5], and consequently the uncertainty on   and  is improved.
The inference of   is challenging due to the behaviour of Eq. ( 27). Figure 5 plots the theoretical model for   as a function  of  (or equivalently ).The posterior of  is shown in red, and the true   is plotted in green.Note the asymptotic behaviours of   at  = 3 and 5; consequently, when a neutron star has  ≈ 5, any small deviation of the estimated  will drastically affect the estimation of   .Therefore, while the posterior of   is supposed to contain the true value since the posterior of  contains the true value, in reality, the fact that we obtain the posterior through sampling, coupled with the asymptotic behaviour of   , makes it unlikely that any sampled  will be close enough to  truth to produce a good estimate of   .This is unsurprising: given that we are inferring neutron star properties from continuous wave detections, we expect most of the detected neutron stars to have a continuous wave-dominated spindown with  ≈ 5.As the magnetic dipole moment is related to the strength of the electromagnetic radiation, which is not the dominant form of radiation for these stars, the amount of energy emitted via electromagnetic radiation, and hence the ability to infer   through this method, is for neutron stars  ≈ 5. Similarly, the skewed distributions of   and  in Figs.4a and 4b are also a consequence of the theoretical framework [Eqs.( 25) and ( 26)] and the posterior on .Both equations have asymptotes at  = 3, with lim →3 +   = +∞ and lim →3 +  = 0.As the posterior on  approaches 3, therefore,   and  will increasingly diverge from their true values.Unlike   , however, a much higher fraction of posterior samples (when  > 3) produce good estimates of   and ; as seen in Fig. 4, the highest densities in the   - distributions coincide with the true values for both LIGO and ET-D.In Fig. 6, we show the inferred physical properties of the 10 neutron stars with the largest ℎ 0 in the simulated population, in descending order of ℎ 0 .The posterior distributions are shown in the form of violin plots, in which the widths of the "violins" give the posterior probability densities.The red lines represent the highest density intervals (HDI) of the distributions, defined as the smallest possible 90% credible intervals, i.e. the smallest bounds we can place on the physical properties.The varying sizes of the distributions are due to the different extents to which  is constrained.For neutron stars where  is well constrained, the uncertainties in the inferred properties are also smaller.This is confirmed in Fig. 7: for neutron stars with higher  , the uncertainties in the inferred parameters, as quantified by the relative error E, are smaller.This can be understood by considering the continuous wave signal searched for by CWInPy.The contribution of  to the frequency evolution of a continuous wave is more pronounced if  is larger, and hence CWInPy can more readily constrain larger  .From Fig. 7 we see that, for five years of continuous wave observations using ET-D, the physical properties of the 10 neutron stars in the synthesised population have errors that depend on the estimation of , with  being a contributing factor.For neutron stars with small  ,   can be inferred with an error of ∼ 100% for point estimates using the median.The 90-th percentile credible interval yields a maximum error of ∼ 1000%.For the ellipticity , a point estimate with relative error ∼ 50% can be made, and the 90-th percentile credible interval giving a maximum error of ∼ 100%.For neutron stars with more significant  , the median of the posterior of   is accurate enough to be within ∼ 10% of the true value, while the 90-th percentile credible interval gives a maximum error of ∼ 100%.The ellipticity  can be estimated using median with a relative error of ∼ 5%, and a maximum error from the credible interval of ∼ 50%.Few alternative methods exist to measure   .Separate measurement of neutron star mass and radius can provide a measurement of   , but measuring these two properties simultaneously is challenging (Steiner et al. 2015;Miller et al. 2019).Another technique measures   through the higher-order relativistic correction to the periastron advance of rapidly spinning binary pulsars (Damour & Schäeer 1988).To date, this approach can only be applied to the double pulsar PSR J0737-3039, with the measurement yielding errors of ∼ 10-20% (Miao et al. 2022), comparable with our point estimates of ∼ 10% for neutron stars with large  .So far, no alternative approach to measuring  is available other than through a continuous wave detection (Lu et al. 2023).We are unable to accurately infer the perpendicular magnetic dipole moment, but it may be inferred alternative methods, such from pulsars directly through the measurement of frequency and spindown (Kramer 2005).

DISCUSSION
This paper presents an analysis of the capability of LIGO and ET to measure physical properties of neutron stars using continuous gravitational waves.We first synthesised a population of neutron stars using Monte-Carlo techniques, performed parameter estimation using Bayesian inference on the ten gravitationally loudest simulated continuous wave signals, and finally inferred their physical properties, namely the stellar moment of inertia   , equatorial ellipticity , and the perpendicular magnetic dipole moment   .Targeted searches for continuous waves from the synthesised neutron stars produced well-constrained posteriors for the continuous wave strain amplitude ℎ 0 , frequency  , and its first order time derivative  , for both LIGO and ET-D.The inference of the braking index  proved challenging, and was the cause of the majority of the uncertainty in the inference of the physical properties   , ,   , but this can be improved with longer observation periods.Using the ET-D configuration with five years of observations, depending on estimation of , which is related to the size of  , 90-th percentile credible intervals can be placed on   and  with errors of ∼ 100 − 1000% and ∼ 50 − 100%, respectively; and point estimates using median can be made with errors of ∼ 10 − 100% and ∼ 5 − 50%, respectively.The perpendicular magnetic dipole moment could not be properly inferred for neutron stars with  ≈ 5 due to the asymptotic behaviour of Eq. ( 27).Lu et al. (2023) found that measurement of   , , and   , with relative errors of ≲ 27%, might be achievable with continuous waves.In comparison, our headline results (summarised above) are somewhat less promising for the accurate estimation of these parameters.The reason for this discrepancy is readily apparent from Table 1; the 10 neutron stars used for parameter estimation all have relatively small spindowns (|  | ≲ 1 × 10 −10 Hz/s, |  | ≲ 1 × 10 −22 Hz/s 2 ) compared to what is typically assumed for continuous wave searches.Indeed, as seen in Fig. 4 of Lu et al. (2023), the |  | of the 10 neutron stars are similar to (or even smaller than) the smallest |  | ∼ 1 × 10 −11 Hz/s covered by Lu et al. (2023), where errors were found to be worst (E ≳ 30).Our results, therefore, reaffirm the intuition that the detection, not only of continuous waves, but from a relatively luminous neutron star with high spindown, would be required for any useful estimation of its physical parameters.
Future work could improve upon some of the simplifying assumptions used in the neutron star population synthesis, such as performing -body simulations of the spatial evolution of the neutron stars through the Galactic potential, or adopting evolution models for the ellipticity and the magnetic field strength.In addition, directed or allsky searches could be performed on the synthesised population to provide a more comprehensive assessment of the detection prospects for LIGO and the Einstein Telescope.

FigureFigure 1 .
Figure2asummarises the simulated population, containing 10 6 neutron stars, in a - diagram.The shape of the simulated population is a result of the various assumptions and limitations of the simulation

Figure 2 .
Figure 2. (a) - diagram of a representative synthesised population (coloured pixels) with observational data from the ATNF superimposed (crosses).The top 10 neutron stars with the largest ℎ 0 are plotted as blue dots, from which we then perform parameter estimations and will be detailed in Section 4. A line of constant magnetic field strength  is drawn in green; lines of constant ellipticity  are drawn in red; lines of constant characteristic age  are drawn in purple.(b) A schematic illustration of four different regions in the - diagram; see the text for details.Example evolution trajectories from  = 0 to  =  age for four neutron stars are labelled (1) to (4).

Figure 4 .
Figure 4. Corner plots for the distributions of the logarithm of the principal moment of inertia   , ellipticity  , and perpendicular magnetic dipole moment   .Results are for the example simulated neutron star J1704-0203, and (a) Advanced LIGO for the duration of O1-O3, (b) ET-D with  = 2 yrs, and (c) ET-D with  = 5 yrs.

Figure 5 .
Figure 5. Theoretical model by Lu et al. (2023) of   as a function of  on top of the posterior distribution of  for ET-D at  = 5 yrs.The true value is plotted in green.

Figure 6 .
Figure 6.Violin plots showing the distributions of (a)   , (b)  , and (c)   for the parameter estimations with CWInPy on ten neutron stars, arranged in descending order of ℎ 0 .The 90% highest density intervals from each run are shown in red; the medians are shown in blue; and true values are shown as green crosses.ET-D is used with  = 5 yrs.

Figure 7 .
Figure 7. Magnitudes of relative errors in   (filled markers) and  (hollow markers) plotted against  for ET-D with  = 5 yrs.The maximum errors of the 90% HDI are shown in red, and that of medians are shown in blue.

Table 1 .
Continuous wave parameters and physical properties of the ten simulated neutron stars with the largest continuous wave strain amplitude ℎ 0 .