TY - JOUR AU1 - Gittins,, Fabian AU2 - Andersson,, Nils AB - ABSTRACT The fastest-spinning neutron stars in low-mass X-ray binaries, despite having undergone millions of years of accretion, have been observed to spin well below the Keplerian break-up frequency. We simulate the spin evolution of synthetic populations of accreting neutron stars in order to assess whether gravitational waves can explain this behaviour and provide the distribution of spins that is observed. We model both persistent and transient accretion and consider two gravitational-wave-production mechanisms that could be present in these systems: thermal mountains and unstable rmodes. We consider the case of no gravitational-wave emission and observe that this does not match well with observation. We find evidence for gravitational waves being able to provide the observed spin distribution; the most promising mechanisms being a permanent quadrupole, thermal mountains, and unstable r modes. However, based on the resultant distributions alone, it is difficult to distinguish between the competing mechanisms. accretion, accretion discs, gravitational waves, stars: neutron, stars: rotation, X-rays: binaries 1 INTRODUCTION Neutron stars (NSs) are among the most compact objects that we observe in our Universe. For this reason, they are able to spin up to extremely high frequencies. The classical picture for the evolution of rapidly spinning NSs begins with an NS accreting gas via a circumstellar accretion disc from their companion in a low-mass X-ray binary (LMXB; Alpar et al. 1982; Radhakrishnan & Srinivasan 1982). This process causes the NS to spin up and, eventually, the NS accretes all the gas from its companion such that all that is left of the binary is a radio millisecond pulsar (RMSP). This scenario should, in theory, have no difficulty in spinning NSs up to their centrifugal break-up frequency (Cook, Shapiro & Teukolsky 1994), which is generally above |$\sim {1}{\, \mathrm{kHz}}$| for most equations of state (Lattimer & Prakash 2007). However, the fastest-spinning pulsar that has been observed to date is PSR J1748-2446ad that spins at 716 Hz (Hessels et al. 2006), well below the limit set by the break-up frequency. In fact, the distribution of spins for both LMXBs and RMSPs has been shown to have a statistically significant cut-off at 730 Hz (Chakrabarty et al. 2003; Patruno 2010). Accreting millisecond X-ray pulsars (AMXPs) are a sub-class of LMXBs that have been spun up to millisecond periods through accretion (see Patruno & Watts 2012 for a review on AMXPs). They are characterized by accretion rates |$\gtrsim 10^{-11}{\, \mathrm{M}_{\odot }\, \mathrm{yr}^{-1}}$| and comparatively weak magnetic fields (⁠|$\sim 10^{8}{\, \mathrm{G}}$|⁠). Another important sub-class of LMXBs are the nuclear-powered X-ray pulsars (NXPs). These pulsars show short-lived burst oscillations during thermonuclear burning on their surfaces and are distinct from AMXPs due to being powered by nuclear burning rather than accretion. 19 AMXPs and 11 NXPs have been observed to date. Haskell et al. (2018) show that the observed rotation rate limit on NSs does not correspond to centrifugal break-up and argue that additional spin-down torques are required to explain this effect. It is unclear what physical process prevents these NSs from spinning up to sub-millisecond periods. One candidate is the interaction between the magnetic field and the accretion disc (Ghosh & Lamb 1978; White & Zhang 1997; Andersson et al. 2005). Patruno, Haskell & D’Angelo (2012) demonstrated that a magnetic-field strength at the magnetosphere of |$\sim 10^{8}{\, \mathrm{G}}$| could be enough to explain the deficiency in accreting NSs above |$\sim {700}{\, \mathrm{Hz}}$|⁠. More recently, it has been shown that transient accretion can have a significant impact on the spin evolution of an accreting NS (Bhattacharyya & Chakrabarty 2017; D’Angelo 2017). In fact, Bhattacharyya & Chakrabarty (2017) noted that in the case of transient accretion, a magnetosphere of |$\sim 10^{8}{\, \mathrm{G}}$| would no longer be sufficient to explain the spin limit. It was suggested by Bildsten (1998) and Andersson, Kokkotas & Stergioulas (1999) that one would observe a spin-frequency limit to accreting NSs if they were emitting gravitational waves (GWs), thus providing a torque to balance the accretion torques (Papaloizou & Pringle 1978; Wagoner 1984). Rapidly rotating NSs have been the subject of many GW searches. These include searches of known radio pulsars (Abbott et al. 2004, 2005b, 2007a, 2008b, 2010, 2017c,d, 2018b; Abadie et al. 2011a,b; Aasi et al. 2014, 2015a), along with wide-parameter surveys for unknown pulsars (Abbott et al. 2005a, 2007b, 2008a, 2009, 2016, 2017b, 2018a; Abadie et al. 2012; Aasi et al. 2013). Assuming torque balance, then the brightest X-ray sources should thus be the loudest in GW emission. Hence, Scorpius X-1 has long been considered a potential candidate for GW emission and has been the focus of a number of targeted GW searches (Abbott et al. 2007b, 2017a,e; Aasi et al. 2015b). These searches, while only yielding upper limits so far, have helped develop and improve data analysis techniques that can be used in the future, when GW detectors further increase their sensitivities. Detecting GWs from rotating NSs will be a challenge, since only a few of the most rapidly accreting NSs are detectable with the current generation of GW detectors (Watts et al. 2008). The main limiting factors in detecting GWs from these systems are the precision with which the spin and the orbital parameters are measured. If these are not well known, then it becomes computationally very expensive to run GW searches, since the searches need to run over a large parameter space. Patruno, Haskell & Andersson (2017) performed a statistical analysis of accreting NSs in LMXBs and found that there is statistical evidence for the distribution to comprise two sub-populations: one at relatively low spin frequencies with an average of |$\approx {300}{\, \mathrm{Hz}}$| and the other at higher frequencies with an average of |$\approx {575}{\, \mathrm{Hz}}$|⁠. The two sub-populations were shown to transition from one to the other in the region around |$\approx {540}{\, \mathrm{Hz}}$|⁠. Patruno et al. (2017) suggested that GWs would provide a physical meaning to this transition region – above this point, GWs become significant. However, Patruno et al. (2017) noted that some of the sources that belong to the fast sub-population have a spin behaviour that is not straightforwardly reconciled with GW scenarios [see Hessels (2008) and Papitto et al. (2014) for additional work on the spin-frequency distribution of millisecond pulsars]. In this paper, we investigate whether an additional component is needed in order to explain the spin evolution of accreting NSs. We consider GW emission as the source of such a component and explore what might be the dominant GW-production mechanism. This paper is structured as follows. To begin with, in Section 2, we introduce the basic theory regarding accretion in LMXBs and discuss how transient accretion can affect the spin evolution of accreting NSs. In Section 3, we provide a brief review of gravitational radiation in these systems and the different mechanisms that can give rise to such radiation. We describe our model for the spin evolution of accreting NSs in Section 4. In Section 5, we summarize the results of our NS population simulations, which include the different GW-production mechanisms. Finally, we conclude and suggest future work in Section 6. 2 ACCRETION IN LMXBs In an LMXB, the companion star has overfilled its Roche lobe and is donating matter to the NS through the inner Lagrange point. Given this donated gas has some large specific angular momentum, it cannot be transferred directly to the surface of the NS and instead forms a circumstellar accretion disc around it. The gas from the accretion disc is channelled on to the magnetic poles of the NS along the magnetic-field lines. This channelling occurs at what is known as the magnetospheric radius, rm, which is the characteristic radius where the magnetic field dominates interactions. At this boundary, the magnetic field can, in turn, be distorted by the accreting gas. This coupling, between the field lines and the disc, results in a torque acting on the star that can spin it up or down depending on the relative difference between rm and the co-rotation radius, \begin{eqnarray*} r_\text{c} \equiv \left(\frac{G M}{\Omega ^2} \right)^{1 / 3}, \end{eqnarray*} (1) which is the location of a Keplerian disc that rotates with the same frequency of the NS, where G is the gravitational constant, M is the mass of the NS, Ω = 2|$\pi$|ν = 2|$\pi$|/P is its angular frequency, ν is the spin frequency, and P is the spin period. If rm < rc, the NS spins slower than the accretion disc and so the gas that is channelled on to the NS has greater specific angular momentum than it, thus acting to spin it up. Conversely, if rm > rc, the NS spins faster than the disc, which spins it down. The magnetospheric radius is a somewhat poorly understood quantity. It is generally defined as the point where the kinetic energy of in-falling gas becomes comparable to the magnetic energy of the magnetosphere. For the straightforward case where the gas is radially accreted on to the NS, one can calculate the Alfvén radius, rA, from \begin{eqnarray*} \frac{1}{2} \rho (r_\text{A}) v(r_\text{A})^2 = \frac{B(r_\text{A})^2}{8 \pi }, \end{eqnarray*} (2) where ρ(rA), v(rA), and B(rA) are the gas density, gas velocity, and magnetic-field strength, respectively, at rA. This calculation gives the standard expression for the magnetospheric radius (Pringle & Rees 1972): \begin{eqnarray*} r_\text{A} = \left(\frac{\mu ^4}{2 G M \dot{M}^2} \right)^{1/7}, \end{eqnarray*} (3) where μ = BR3 is the magnetic moment of the NS, |$\dot{M}$| is the mass accretion rate from the disc to the NS surface, and R is the radius of the NS. This picture becomes more complicated when considering accretion from a circumstellar disc. A factor ξ of order unity is introduced to correct for the non-spherical geometry of the problem and also account for the extended transition region between where mass is accreted on to the star and where mass is ejected in an outflow. This gives the magnetospheric radius as \begin{eqnarray*} r_\text{m} = \xi r_\text{A} = \xi \left(\frac{\mu ^4}{2 G M \dot{M}^2} \right)^{1/7}. \end{eqnarray*} (4) Typically, ξ is assumed to fall in the range of 0.5–1.4 (Wang 1996). This correction demonstrates that rm is sensitive to the coupling between the field lines and the accretion disc that plays a key role in understanding the magnetospheric radius. 2.1 Accretion-torque models The spin evolution of an NS is dictated by the torques acting upon the star. By measuring the time derivative of the spin period, |$\dot{P}$|⁠, one can gain insight into the physics of accretion processes, as well as other aspects such as magnetic-field strengths and GW emission. The variation in spin is related to the torque exerted on to the NS by the standard expression \begin{eqnarray*} N = -\frac{2 \pi I \dot{P}}{P^2}, \end{eqnarray*} (5) where I is the stellar moment of inertia. For an NS accreting from an accretion disc truncated at the magnetospheric transition, with rm < rc, the standard torque is (neglecting magnetic-field effects) \begin{eqnarray*} N_\text{acc} = \dot{M} r_\text{m}^2 \Omega _\text{K}(r_\text{m}) = \dot{M} \sqrt{G M r_\text{m}}, \end{eqnarray*} (6) where ΩK(rm) is the Keplerian angular velocity at rm. For rm > rc, the coupling between the magnetosphere and the accretion disc becomes important as the magnetic-field lines are threaded through the disc so extra torques due to magnetic stresses come into play. Ghosh & Lamb (1979) developed an accretion model based on detailed calculations of this coupling. The torque from this model predicts \begin{eqnarray*} \dot{P} \approx -5.0 \times 10^{-5} M_{1.4}^{3/7} I_{45}^{-1} \mu _{30}^{2/7} \left[ \left(\frac{P}{{1}{\, \mathrm{s}}} \right) \dot{M}_{-9}^{3/7} \right]^2 n(\omega _\text{s}) \, {\, \mathrm{s}\, \mathrm{yr}^{-1}}, \end{eqnarray*} (7) where |$M_{1.4} = M / {1.4}{\, \mathrm{M}_{\odot }}$|⁠, |$R_6 = R / {10}{\, \mathrm{km}}$|⁠, |$I_{45} = I / 10^{45}{\, \mathrm{g\, cm^{2}}}$|⁠, |$\mu _{30} = \mu / 10^{30}{\, \mathrm{G}\, \mathrm{cm^{3}}}$|⁠, |$\dot{M}_{-9} = \dot{M} / 10^{-9}{\, \mathrm{M}_{\odot }\, \mathrm{yr}^{-1}}$|⁠, and n(ωs) is the dimensionless torque that accounts for the magnetic field-accretion disc coupling and is a function of the fastness parameter ωs. The fastness parameter is defined as the ratio of the NS spin frequency to the Keplerian orbital frequency at the magnetospheric boundary, \begin{eqnarray*} \omega _\text{s} \equiv \frac{\Omega }{\Omega _\text{K}(r_\text{m})} \approx {3.1} \, \xi ^{3/2} M_{1.4}^{-5/7} \mu _{30}^{6/7} \left[ \left(\frac{P}{{1}{\, \mathrm{s}}} \right) \dot{M}_{-9}^{3/7} \right]^{-1}. \end{eqnarray*} (8) The sign of n(ωs) depends on whether the NS accretes the gas and spins up (the ‘slow rotator’ regime, ωs < 1) or ejects the gas and spins down (the ‘fast rotator’ regime, ωs > 1; Wang 1995). It is interesting to note that an NS can still be spun down at long spin periods (⁠|$P \gg {1}{\, \mathrm{s}}$|⁠) as the magnetic field can be strong enough to mean it would still be classified as a fast rotator. Ho et al. (2014) introduced a simple approximation to the Ghosh & Lamb model by considering angular momentum changes on the NS. Matter accreting at the magnetosphere, rm, has specific angular momentum \begin{eqnarray*} l_\text{acc} = \pm r_\text{m}^2 \Omega _\text{K}(r_\text{m}), \end{eqnarray*} (9) where the sign of lacc depends on whether there is prograde rotation between the accretion disc and the NS (lacc > 0) or retrograde rotation (lacc < 0). Prograde rotation will be assumed. What must also be accounted for is matter that is ejected from the NS, which will carry specific angular momentum \begin{eqnarray*} l_\text{m} = r_\text{m}^2 \Omega . \end{eqnarray*} (10) Both of these effects produce a torque on the NS. The net torque is obtained by summing these contributions: \begin{eqnarray*} N = \dot{M} (l_\text{acc} - l_\text{m}) = \dot{M} r_\text{m}^2 \Omega _\text{K}(r_\text{m}) (1 - \omega _\text{s}). \end{eqnarray*} (11) In this relation, one can see the standard spin-up torque due to disc accretion (equation 6) which is corrected by the fastness parameter to account for interactions spinning down the NS. (This model phenomenologically accounts for effects such as accretion disc-magnetic field coupling and outflows.) This expression can be related to the change in spin period using equation (5) to obtain \begin{eqnarray*} \dot{P} &\approx& -8.1\times 10^{-5} \, \xi ^{1/2} M_{1.4}^{3/7} I_{45}^{-1} \mu _{30}^{2/7} \left[ \left(\frac{P}{{1}{\, \mathrm{s}}} \right) \dot{M}_{-9}^{3/7} \right]^2\nonumber\\ &&\times \,(1 - \omega _\text{s}) \, {\, \mathrm{s}\, \mathrm{yr}^{-1}}. \end{eqnarray*} (12) It is clear from equation (12) that the fastness parameter dictates whether the NS spins up or down. A commonly considered aspect of an NS’s spin evolution is the spin equilibrium. This occurs when the spin rate is gradually adjusted until the net torque on the star is approximately zero and the accretion flow is truncated at the magnetospheric radius, rm ≃ rc. When an NS reaches spin equilibrium, it is straightforward to estimate its magnetic field, assuming that the accretion rate and spin are known. One can estimate the spin-equilibrium period, Peq, from equation (12) by setting |$\dot{P} = 0$|⁠, when ωs = 1: \begin{eqnarray*} P_\text{eq} &\approx& {8.2} \, \xi ^{3/2} M_{1.4}^{-5/7} \left(\frac{\mu }{10^{26}{\, \mathrm{G}\, \mathrm{cm^{3}}}}\right)^{6/7}\nonumber\\ &&\times \,\left(\frac{\dot{M}}{10^{-11}{\, \mathrm{M}_{\odot }\, \mathrm{yr}^{-1}}}\right)^{-3/7} \, {\mathrm{ms}}, \end{eqnarray*} (13) where we have scaled the period to characteristic AMXP values. The magnetic-field lines rotate with the NS. This produces magnetic-dipole radiation that causes the NS to spin down. The torque due to this is described by \begin{eqnarray*} N_\text{EM} = - \frac{2 \mu ^2 \Omega ^3}{3 c^3}, \end{eqnarray*} (14) where c is the speed of light in a vacuum. The change in spin due to magnetic-dipole radiation is \begin{eqnarray*} \dot{P}_\text{EM} \approx 3.1\times 10^{-8} \mu _{30}^2 I_{45}^{-1} \left(\frac{P}{{1}{\, \mathrm{s}}} \right)^{-1} {\, \mathrm{s}\, \mathrm{yr}^{-1}}. \end{eqnarray*} (15) The vast majority of pulsars are isolated and their spin evolution can be generally described by magnetic-dipole radiation. However, in the case of rapidly accreting NSs this effect can be essentially negligible. There are more accurate numerical models one can use to describe these torques (e.g. see Spitkovsky 2006). 2.2 Transient accretion Up until now, the accretion rate has been implicitly assumed to be steady. However, many LMXBs exhibit long periods of quiescence, which can be of the order of months to years, and short transient outbursts, which can last from days to weeks. These outbursts are believed to be caused by instabilities in the accretion disc and occur when the mass-accretion rate rises above a certain threshold (see e.g. Lasota 1997). As the companion star donates a steady flow of gas to the accretion disc, the disc gets larger and eventually reaches a critical mass to trigger an instability. This causes the accretion rate from the disc to the surface of the NS to increase by several orders of magnitude, giving rise to a transient outburst. Once the accretion disc has donated a sufficient amount of gas the system returns to a quiescent state until a new outburst occurs when the disc has accumulated enough mass from the companion and the cycle repeats (Done, Gierliński & Kubota 2007). D’Angelo (2017) and Bhattacharyya & Chakrabarty (2017) have shown that transient accretion with a varying accretion rate has a significant impact on the spin evolution of an NS. Both found that for a given long-term average accretion rate, these transients can spin up NSs to rates several times higher than that of persistent accretors; however, it takes approximately an order of magnitude longer to reach these spin-equilibrium periods. This demonstrates that for transient systems, like most LMXBs, it is not accurate to assume a time-averaged accretion rate but instead one must consider the outburst/quiescence phases. D’Angelo (2017) noted that the two key changes when considering transient accretion are that the torque over an outburst is significantly smaller than for the persistent case at a given accretion rate and that the equilibrium accretion rate is shifted to a lower value. This has the combined effect to increase the time it takes for a transient source to reach spin equilibrium and decrease its spin-equilibrium period. D’Angelo (2017) and Bhattacharyya & Chakrabarty (2017) found that the spin-equilibrium period and time to reach spin equilibrium are sensitive to the features of the accretion profile. They found that by increasing the duration of an outburst by a factor of 10 the spin-equilibrium period can decrease by up to a factor of 2. For their analysis, D’Angelo (2017) used a fast-rise, exponential-decay function to model the accretion profile (whereas Bhattacharyya & Chakrabarty 2017, used a simple sawtooth function): \begin{eqnarray*} f(t) = \exp \left(\sqrt{\frac{2}{F_\text{t}}} \right) \exp \left(- \frac{1}{10 t} - \frac{10 t}{F_\text{t}} \right) + f_\text{min}, \end{eqnarray*} (16) where t denotes the time from the beginning of the outburst, Ft is an approximate measure of the duration of the outburst and fmin is the minimum. It should be noted that this function models a single outburst/quiescence cycle, so in order to model multiple cycles one repeats this after a given recurrence time, Trecurrence. Time has arbitrary units in this model. The ratio of the maximum to the minimum is \begin{eqnarray*} \frac{f_\text{max}}{f_\text{min}} = \frac{1}{f_\text{min}} \exp \left(\frac{\sqrt{2} - 2}{\sqrt{F_\text{t}}} \right) + 1. \end{eqnarray*} (17) This accretion profile requires two normalizations. The first normalization chooses fmax/fmin to obtain fmin for a fixed Ft using equation (17). The second normalization is to demand that 〈f(t)〉 = 1. This normalization depends on Trecurrence and results in fmin no longer corresponding precisely to the minimum value. These normalizations allow one to choose the magnitude of the accretion outburst, with respect to the quiescent accretion rate, and also mean that one can simply choose an average accretion rate over one cycle by multiplying equation (16) by the chosen average. Thus, the time-dependant accretion rate is given by \begin{eqnarray*} \dot{M}(t) = \langle \dot{M} \rangle f(t), \end{eqnarray*} (18) where f(t) has been appropriately normalized. The canonical profile used by D’Angelo (2017) had Ft = 10, fmax/fmin = 693.97 and Trecurrence = 100 and is shown in Fig. 1. Figure 1. View largeDownload slide Accretion outburst profile, f(t), as a function of time, t, where Ft = 10, fmax/fmin = 693.97, and Trecurrence = 100. The accretion rate and time have arbitrary units. Figure 1. View largeDownload slide Accretion outburst profile, f(t), as a function of time, t, where Ft = 10, fmax/fmin = 693.97, and Trecurrence = 100. The accretion rate and time have arbitrary units. 3 GRAVITATIONAL RADIATION FROM ACCRETING NSs Rotating NSs emit GWs if they are asymmetric about the axis of rotation. Such mass asymmetries are referred to as mountains. For a spinning NS with a mass-quadrupole moment Q22, the braking torque due to GWs is given by \begin{eqnarray*} N_\text{GW} = - \frac{256 \pi }{75} \frac{G \Omega ^5 Q_{2 2}^2}{c^5}. \end{eqnarray*} (19) This corresponds to a spin-down rate of \begin{eqnarray*} \dot{P}_\text{GW} \approx 1.4\times 10^{-19} I_{45}^{-1} \left(\frac{Q_{2 2}}{10^{37}{\, \mathrm{g\, cm^{2}}}} \right)^2 \left(\frac{P}{{1}{\, \mathrm{s}}} \right)^{-3} \, {\, \mathrm{s}\, \mathrm{yr}^{-1}}. \end{eqnarray*} (20) In order to estimate how strong a quadrupole is needed in order to considerably influence the spin evolution of the NS it is useful to balance equation (19) with the accretion-magnetosphere torque from equation (11). This leads to \begin{eqnarray*} Q_{2 2} \approx 4.2&\times& 10^{37} \xi ^{1 / 4} M_{1.4}^{3 / 14} \mu _{30}^{1 / 7} \dot{M}_{-9}^{3 / 7} \left(\frac{\nu }{{500}{\, \mathrm{Hz}}}\right)^{- 5 / 2} \nonumber \\ &\times& \,(1 - \omega _\text{s}) \, {\, \mathrm{g\, cm^{2}}}. \end{eqnarray*} (21) For a typical AMXP with |$B \sim 10^{8}{\, \mathrm{G}}$| and |$\dot{M} \sim 10^{-11}{\, \mathrm{M}_{\odot }\, \mathrm{yr}^{-1}}$|⁠, this gives a quadrupole moment of |$Q_{2 2} \sim 10^{36}{\, \mathrm{g\, cm^{2}}}$| in order to achieve spin equilibrium at |$\nu \sim {500}{\, \mathrm{Hz}}$|⁠. One can express a mass-quadrupole in terms of a moment of inertia ellipticity, defined as (see e.g. Owen 2005) \begin{eqnarray*} \epsilon \equiv \sqrt{\frac{8 \pi }{15}} \frac{Q_{2 2}}{I}. \end{eqnarray*} (22) Therefore, in order to balance the accretion torque with GW spin-down, one requires ϵ ∼ 10−9. This is far smaller than the maximum deformation an NS crust can sustain for most reasonable equations of state (Johnson-McDaniel & Owen 2013). Recent population-based analysis has suggested that ϵ ≈ 10−9 is the minimum ellipticity of millisecond pulsars (Woan et al. 2018). An outstanding problem in understanding rapidly spinning accreting NSs is their peculiar spin distribution (see Fig. 2). It is this unusual shape and, in particular, the sharp cut-off at |$\sim {600}{\, \mathrm{Hz}}$| that has motivated the search for GWs from these systems. This is an appealing explanation since the braking torque due to GWs scales as the fifth power of the spin frequency for deformed, rotating NSs. Patruno et al. (2017) have shown that among AMXPs and NXPs, there appear to be two sub-populations. One sub-population is at a relatively low spin frequency, with a mean spin of |$\approx {300}{\, \mathrm{Hz}}$|⁠. The second sub-population has a higher peak and a mean of |$\approx {575}{\, \mathrm{Hz}}$|⁠. This faster sub-population has a very narrow range and is composed of a mixture of AMXPs and NXPs. The two sub-populations are separated by a transition region around |$\approx {540}{\, \mathrm{Hz}}$|⁠. Patruno et al. (2017) argued that, when considering various accretion torque models, no model naturally explains the presence of a fast sub-population and postulated that whatever mechanism that causes this clustering it must set in quickly – as soon as the pulsars reach a certain spin threshold. It was noted by Patruno et al. (2017) that this is a subtly different problem to the one of accreting NSs not spinning close to their break-up frequency. These two problems make GWs a promising avenue to explore. GWs can help justify the transition region between the two sub-populations and provide a physical meaning to it (the region in which GW emission starts to become significant), and naturally explain the cut-off at |$\sim {600}{\, \mathrm{Hz}}$|⁠. Figure 2. View largeDownload slide Distribution of spin frequencies for accreting NSs with millisecond periods. The accreting NS population comprises AMXPs, coloured in blue, and NXPs, coloured in orange. Figure 2. View largeDownload slide Distribution of spin frequencies for accreting NSs with millisecond periods. The accreting NS population comprises AMXPs, coloured in blue, and NXPs, coloured in orange. There are a number of different ways a mass asymmetry could arise in an accreting NS. Bildsten (1998) originally proposed that interior temperature asymmetries misaligned with respect to the spin-axis of the NS could produce a significant quadrupole through temperature-sensitive electron captures. Hotter regions of the crust would have electron captures at lower pressures and so the density drop would occur at higher altitudes in the hotter parts of the crust. This is known as a thermal mountain (Bildsten 1998; Ushomirsky, Cutler & Bildsten 2000; Melatos & Payne 2005; Haskell, Jones & Andersson 2006; Payne & Melatos 2006; Johnson-McDaniel & Owen 2013). Another mechanism through which mass-quadrupoles can be built are through mountains sustained by magnetic stresses, called magnetic mountains (Cutler 2002; Haskell et al. 2008). These can occur when an NS has a sufficiently large toroidal or poloidal magnetic field that will act to distort the NS into an oblate or prolate shape and will naturally produce a quadrupole if the spin- and magnetic-axes are misaligned. A third way through which GWs can arise is through internal r mode instabilities (Andersson 1998; Andersson et al. 1999; Levin 1999; Andersson et al. 2000; Andersson, Jones & Kokkotas 2002; Heyl 2002; Wagoner 2002; Nayyar & Owen 2006; Bondarescu, Teukolsky & Wasserman 2007). In a perfect fluid, these modes are unstable for all rates of rotation due to GW emission. We explored whether GWs could explain the observed distribution and, if so, whether there is a preference for any of the GW-production mechanisms. For our analysis, we did not consider mountains solely created by the magnetic field, nor did we consider magnetic mountains built through accretion. For these cases, the magnetic fields are not strong enough to sustain sufficiently large mountains for the spin evolution of these systems to be noticeably affected (see Haskell et al. 2008; Priymak, Melatos & Payne 2011; Haskell & Patruno 2017). 4 SPIN-EVOLUTION MODEL For this work, we constructed a model for the spin evolution of an accreting NS. We incorporated the accretion–magnetosphere coupling by using the model of Ho et al. (2014, equation 12) and included a torque due to GW spin-down (equation 20). The spin rate is a first-order time derivative and so can be evolved numerically. Our assumed canonical values for an LMXB are shown in Table 1. For our canonical NS, we did not include GW effects. We assumed constant density for our NSs, which affected the moment of inertia. For simplicity, we did not model the magnetic-field evolution. The time an NS is evolved for is denoted as the evolution time. Table 1. Canonical values for an LMXB. Mass / M⊙ Radius / km B / G Initial spin period / s ξ |$\langle \dot{M} \rangle$| / |$\mathrm{M}_{\odot }\, \mathrm{yr}^{-1}$| Q22/ |$\mathrm{g} \, \mathrm{cm}^2$| 1.4 10 108 0.1 0.5 5 × 10−11 0 Mass / M⊙ Radius / km B / G Initial spin period / s ξ |$\langle \dot{M} \rangle$| / |$\mathrm{M}_{\odot }\, \mathrm{yr}^{-1}$| Q22/ |$\mathrm{g} \, \mathrm{cm}^2$| 1.4 10 108 0.1 0.5 5 × 10−11 0 View Large Table 1. Canonical values for an LMXB. Mass / M⊙ Radius / km B / G Initial spin period / s ξ |$\langle \dot{M} \rangle$| / |$\mathrm{M}_{\odot }\, \mathrm{yr}^{-1}$| Q22/ |$\mathrm{g} \, \mathrm{cm}^2$| 1.4 10 108 0.1 0.5 5 × 10−11 0 Mass / M⊙ Radius / km B / G Initial spin period / s ξ |$\langle \dot{M} \rangle$| / |$\mathrm{M}_{\odot }\, \mathrm{yr}^{-1}$| Q22/ |$\mathrm{g} \, \mathrm{cm}^2$| 1.4 10 108 0.1 0.5 5 × 10−11 0 View Large Our model can evolve both persistent and transient accretors. For transient accretors, we used a fast-rise, exponential-decay function, described in Section 2.2 by equations (16)–(18), and evolved the time-averaged spin-derivative, |$\langle \dot{P}(P, \dot{M}) \rangle$|⁠, which, for a given NS, is a function of the spin and accretion rate. This average was obtained by averaging the spin derivative over one outburst/quiescence cycle. The time-average was evolved rather than the instantaneous spin rate, |$\dot{P}(P, \dot{M})$|⁠, to simplify the integration procedure. Otherwise the integration procedure would have needed to take into account the full fast-rise, exponential-decay features of the accretion profile. For persistent accretors this was not a problem and so we could simply evolve |$\dot{P}(P, \dot{M})$|⁠. Unless specified otherwise we used the following values for the transient accretion profile: |$F_\text{t} = {10}{\, \mathrm{yr}}$|⁠, |$T_\text{recurrence} = {100}{\, \mathrm{yr}}$|⁠, and fmax/fmin = 104. This was chosen for simplicity and to limit the explorable parameter space. Most of our simulations turned out to be relatively insensitive to the exact values of these parameters. Of course, should one be interested in modelling individual systems with this profile then particular care would need to be taken when tuning these parameters. Fig. 3 shows the spin evolution of the canonical accreting NS with persistent and transient accretion. As was found by D’Angelo (2017) and Bhattacharyya & Chakrabarty (2017), we see that the persistently accreting NS initially spins up faster and reaches a final spin of |$\nu = {678}{\, \mathrm{Hz}}$|⁠. The transient system spins up slower but obtains a faster final spin of |$\nu = {1055}{\, \mathrm{Hz}}$|⁠. However, neither of the systems were evolved long enough to reach spin equilibrium. The upper limit of |$10^{10}{\, \mathrm{yr}}$| for the evolution time was chosen since no system can evolve for longer than the age of the Universe. Figure 3. View largeDownload slide The spin evolution of the canonical accreting NS with persistent accretion (blue continuous line) and transient accretion (orange dashed line) with initial values from Table 1. The persistent accretor initially spins up faster than the transient accretor. However, towards the end of its evolution the persistent accretor begins to slow down and the transient accretor overtakes and reaches a faster final spin. Figure 3. View largeDownload slide The spin evolution of the canonical accreting NS with persistent accretion (blue continuous line) and transient accretion (orange dashed line) with initial values from Table 1. The persistent accretor initially spins up faster than the transient accretor. However, towards the end of its evolution the persistent accretor begins to slow down and the transient accretor overtakes and reaches a faster final spin. 5 SIMULATED POPULATIONS In order to obtain a distribution of spins with which to compare to the observed distribution, we used a Monte Carlo population synthesis method to draw the initial parameters from a given set of distributions and evolve each NS (see Possenti et al. 1998, for another NS population synthesis study). Each NS was assigned a mass, radius, magnetic-field strength, initial spin period, average accretion rate, mass-quadrupole moment, and evolution time. We evolved 1000 NSs in each simulation. The first simulations were evolved using the distributions shown in Table 2. We fixed the masses and radii at |${1.4}{\, \mathrm{M}_{\odot }}$| and |${10}{\, \mathrm{km}}$|⁠, respectively, to match the canonical values for NSs. Typically, AMXPs are measured to have magnetic fields of |$\sim 10^{8}{\, \mathrm{G}}$| and so the field strength was taken from a log-Gaussian distribution with mean μ = 8.0 and standard deviation σ = 0.1. The initial spin period was drawn from a flat distribution between |$0.01\,\hbox{and}\,0.1{\, \mathrm{s}}$|⁠, which our simulations turned out to be relatively insensitive to. The correction factor ξ was chosen to be 0.5. The average accretion rate was motivated by observations of LMXBs and was given by a log-Gaussian with μ = −11.0 + log10(5) and σ = 0.1. For the initial simulations we assumed there was no GW component. We found that for evolution times much less than |$10^{9}{\, \mathrm{yr}}$| the NSs would not have enough time to spin up to frequencies above |${100}{\, \mathrm{Hz}}$| and so the evolution time was taken from a flat-in-the-log distribution between |$10^{9}\,\hbox{and}\,10^{10}{\, \mathrm{yr}}$|⁠. The distribution was chosen to be flat-in-the-log in order for it to be scale-invariant (as was used by Possenti et al. 1998), thus parametrizing our uncertainty in the value of the evolution time. Table 2. Initial values and evolution parameters for population synthesis. Parameter Distribution Values Mass / M⊙ Single value 1.4 Radius / km Single value 10 log10(B / G) Gaussian μ = 8.0, σ = 0.1 Initial spin period / s Flat 0.01–0.1 ξ Single value 0.5 |$\log _{10}(\langle \dot{M} \rangle \, /\, {\, \mathrm{M}_{\odot }\, \mathrm{yr}^{-1}})$| Gaussian μ = −11.0 + log10(5), σ = 0.1 Q22 / |$\mathrm{g\, cm^{2}}$| Single value 0 Evolution time / yr Flat-in-the-log 109–1010 Parameter Distribution Values Mass / M⊙ Single value 1.4 Radius / km Single value 10 log10(B / G) Gaussian μ = 8.0, σ = 0.1 Initial spin period / s Flat 0.01–0.1 ξ Single value 0.5 |$\log _{10}(\langle \dot{M} \rangle \, /\, {\, \mathrm{M}_{\odot }\, \mathrm{yr}^{-1}})$| Gaussian μ = −11.0 + log10(5), σ = 0.1 Q22 / |$\mathrm{g\, cm^{2}}$| Single value 0 Evolution time / yr Flat-in-the-log 109–1010 View Large Table 2. Initial values and evolution parameters for population synthesis. Parameter Distribution Values Mass / M⊙ Single value 1.4 Radius / km Single value 10 log10(B / G) Gaussian μ = 8.0, σ = 0.1 Initial spin period / s Flat 0.01–0.1 ξ Single value 0.5 |$\log _{10}(\langle \dot{M} \rangle \, /\, {\, \mathrm{M}_{\odot }\, \mathrm{yr}^{-1}})$| Gaussian μ = −11.0 + log10(5), σ = 0.1 Q22 / |$\mathrm{g\, cm^{2}}$| Single value 0 Evolution time / yr Flat-in-the-log 109–1010 Parameter Distribution Values Mass / M⊙ Single value 1.4 Radius / km Single value 10 log10(B / G) Gaussian μ = 8.0, σ = 0.1 Initial spin period / s Flat 0.01–0.1 ξ Single value 0.5 |$\log _{10}(\langle \dot{M} \rangle \, /\, {\, \mathrm{M}_{\odot }\, \mathrm{yr}^{-1}})$| Gaussian μ = −11.0 + log10(5), σ = 0.1 Q22 / |$\mathrm{g\, cm^{2}}$| Single value 0 Evolution time / yr Flat-in-the-log 109–1010 View Large The resultant spin-frequency distributions for persistent and transient accretors are shown in Fig. 4. One can see for this simple case that for both persistent and transient accretion we do indeed obtain NSs that spin in excess of |${1}{\, \mathrm{kHz}}$|⁠. More generally, we observe that we get many NSs that spin faster than the observed spin-frequency limit of |$\sim {600}{\, \mathrm{Hz}}$| and, as one might expect, we find more high-frequency NSs for the transient case. This is because transient accretion enables these stars to spin to higher frequencies than with persistent accretion, provided they evolve for long enough. For both simulations, we have not obtained the characteristic behaviour of the observed distribution, since there is no evidence for a pile-up of NSs at high frequencies. Figure 4. View largeDownload slide Distributions of spin frequencies for simulated persistently-accreting NSs (left-hand panel) and transiently accreting NSs (right-hand panel) with initial distributions from Table 2. Figure 4. View largeDownload slide Distributions of spin frequencies for simulated persistently-accreting NSs (left-hand panel) and transiently accreting NSs (right-hand panel) with initial distributions from Table 2. In order to quantify how different our simulated populations are to the observed population we applied a Kolmogorov–Smirnov test to the distributions. This test allows one to compare two distributions by testing the null hypothesis that the two distributions are the same. We chose a significance level of α = 0.10 for this work, which meant that should we have found p-values less than this value we could reject the null hypothesis with 90 per cent certainty for that case. For the persistent accretors, we obtained a p-value of p = 7.7 × 10−2 and for the transient accretors we obtained p = 9.9 × 10−4. This meant that we could reject the null hypothesis at the 10 per cent significance level that the observed distribution is drawn from the persistent population or the transient population. We explored the effect that magnetic-dipole radiation (equation 15) has on transient accretors using the initial distributions in Table 2. The results are displayed in Fig. 5. The inclusion of this additional torque stops many of the systems from spinning up to sub-millisecond periods. We obtain p = 0.20, so we cannot rule out the null hypothesis with any statistical certainty. However, in regards to the shape of the distribution, we do not obtain a sharp peak at the observed spin-frequency limit. Instead, we find a broad peak in the range |$200\rm {-}600{\, \mathrm{Hz}}$|⁠. Figure 5. View largeDownload slide Distribution of spin frequencies for simulated transiently accreting NSs with initial distributions from Table 2 including the magnetic-dipole torque. Figure 5. View largeDownload slide Distribution of spin frequencies for simulated transiently accreting NSs with initial distributions from Table 2 including the magnetic-dipole torque. This demonstrates that a simple model for accretion is not sufficient to explain the observations of accreting NSs and also suggests that the inclusion of magnetic-dipole torques do not resolve this tension either. Therefore, an additional component needs to be included into the model. 5.1 Including gravitational-wave torques We explored whether including a GW component to the spin evolution of accreting NSs could give us the observed spin distribution. Motivated by the necessary quadrupole in order to achieve torque balance (equation 21), we repeated the same simulations but with a fixed |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$| for all NSs. (Physically, this can be interpreted as a permanent crustal mountain.) Fig. 6 shows the final-spin distributions for these simulations. This quadrupole has notably stopped the NSs from spinning up to sub-millisecond periods and has resulted in a pile-up centred on the |$500\rm {-}550{\, \mathrm{Hz}}$| bin for the persistent accretors and at |$550\rm {-}600{\, \mathrm{Hz}}$| for the transient accretors. This has appeared since the GW torque imposes a spin-frequency limit on the NSs. The peak for transient accretors is promising as this is where the peak lies for the spin distribution that we observe (cf. Fig. 2). Interestingly, there is also a broader peak at lower frequencies. We found that |$\approx 19{{\ \rm per\ cent}}$| of persistently-accreting NSs and |$\approx 14{{\ \rm per\ cent}}$| of transiently accreting NSs reached spin equilibrium by the end of the simulation. The systems that had reached spin equilibrium were clustered around the high-frequency peaks. We obtained p = 0.19 and p = 0.80 for the persistent and transient cases, respectively, and therefore were unable to reject the null hypothesis for both populations. Figure 6. View largeDownload slide Distributions of spin frequencies for simulated persistently accreting NSs (left-hand panel) and transiently accreting NSs (right-hand panel) with initial distributions from Table 2 and a fixed quadrupole of |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$|⁠. Figure 6. View largeDownload slide Distributions of spin frequencies for simulated persistently accreting NSs (left-hand panel) and transiently accreting NSs (right-hand panel) with initial distributions from Table 2 and a fixed quadrupole of |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$|⁠. We considered how magnetic-dipole radiation affects this picture for systems undergoing transient accretion. We used the same quadrupole and obtained the results shown in Fig. 7. Interestingly, this distribution is qualitatively similar to the results without magnetic-dipole torques (the right-hand panel of Fig. 6). We recover a pronounced peak at higher frequencies |$500\rm {-}550{\, \mathrm{Hz}}$|⁠, which by comparison has shifted down by |${50}{\, \mathrm{Hz}}$|⁠. Since the features of the distribution remain the same we argue that one could obtain a distribution with a peak that matches the observed distribution through slight adjustment of the initial values, e.g. the quadrupole. Such an adjustment would be justifiable since there is significant uncertainty in many of these parameters. We could not reject the null hypothesis for these results with p = 0.39. Figure 7. View largeDownload slide Distribution of spin frequencies for simulated transiently accreting NSs with initial distributions from Table 2 with a fixed quadrupole of |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$| and including the magnetic-dipole torque. Figure 7. View largeDownload slide Distribution of spin frequencies for simulated transiently accreting NSs with initial distributions from Table 2 with a fixed quadrupole of |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$| and including the magnetic-dipole torque. For accreting NS systems the magnetic-dipole torque is expected to be negligible during outbursts, but it could play an important role during the quiescent phases. We explored a range of outburst durations, Ft distributed uniformly between |$1\,\hbox{and}\,100{\, \mathrm{yr}}$|⁠, to assess the impact this made on the resultant spin distribution (Fig. 8). For this wide range of outburst lengths we find a broad peak between |$450\,\hbox{and}\,600{\, \mathrm{Hz}}$|⁠. This contrasts the narrow peak in Fig. 2. We obtain a p-value of p = 0.38. Figure 8. View largeDownload slide Distribution of spin frequencies for simulated transiently accreting NSs with initial distributions from Table 2 with a fixed quadrupole of |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$|⁠, including the magnetic-dipole torque and Ff distributed flat between |$1\,\hbox{and}\,100{\, \mathrm{yr}}$|⁠. Figure 8. View largeDownload slide Distribution of spin frequencies for simulated transiently accreting NSs with initial distributions from Table 2 with a fixed quadrupole of |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$|⁠, including the magnetic-dipole torque and Ff distributed flat between |$1\,\hbox{and}\,100{\, \mathrm{yr}}$|⁠. We investigated how sensitive the simulated populations were on the distribution of the evolution time. We ran a simulation with the same quadrupole and the evolution time distributed flat between |$10^{8}\,\hbox{and}\,10^{10}{\, \mathrm{yr}}$| for transient accretors. This was motivated by assuming that NSs are born at a uniform rate, which is intuitively what one might expect, and that there are no selection effects to suggest that we are more likely to observe younger systems. The resultant spin-distribution is shown in Fig. 9. As was observed in the case when the time was distributed flat-in-the-log, there exists a pronounced peak towards the higher spin frequencies. However, there are far fewer systems spinning at frequencies below this peak. For this simulation, we obtained a p-value of p = 1.2 × 10−2 and thus could reject the null hypothesis. This distribution does not match what we observe. Figure 9. View largeDownload slide Distribution of spin frequencies for simulated transiently accreting NSs with initial distributions from Table 2 except for a fixed quadrupole of |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$| and an evolution time distributed flat between |$10^{8}\,\hbox{and}\,10^{10}{\, \mathrm{yr}}$|⁠. Figure 9. View largeDownload slide Distribution of spin frequencies for simulated transiently accreting NSs with initial distributions from Table 2 except for a fixed quadrupole of |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$| and an evolution time distributed flat between |$10^{8}\,\hbox{and}\,10^{10}{\, \mathrm{yr}}$|⁠. 5.2 Thermal mountains One of the most promising avenues for producing a mass-quadrupole on a fast-spinning, accreting NS is through thermal mountains built during accretion phases through asymmetries in pycnonuclear reaction rates (Haskell & Patruno 2017). As an NS accretes matter composed of light elements, the matter becomes buried by accretion and is then compressed to higher densities. This causes the matter to undergo nuclear reactions such as electron captures, neutron emission, and pycnonuclear reactions (Haensel & Zdunik 1990). If the accretion flow is asymmetric this can cause asymmetries in density and heating which can produce a quadrupole moment. The quadrupole moment due to asymmetric crustal heating from nuclear reactions is approximated by (Ushomirsky et al. 2000) \begin{eqnarray*} Q_{2 2} \approx 1.3\times 10^{37} R_6^4 \left(\frac{\delta T_\text{q}}{10^{7}{\, \mathrm{K}}} \right) \left(\frac{E_\text{th}}{{30}{\, \mathrm{MeV}}} \right)^3 \, {\, \mathrm{g\, cm^{2}}}, \end{eqnarray*} (23) where δTq is the quadrupolar temperature increase due to the nuclear reactions and Eth is the threshold energy for the reactions to occur. The value δTq will be a fraction of the total heating (Ushomirsky & Rutledge 2001), \begin{eqnarray*} \delta T \approx 2&\times& 10^{5} \left(\frac{C}{k_\text{B} \, \text{baryon}} \right)^{-1} \left(\frac{p}{10^{30}{\, \mathrm{erg}\, \mathrm{cm^{-3}}}} \right)^{-1} \nonumber \\ &\times& \,\left(\frac{Q}{{1}{\, \mathrm{MeV}}} \right) \left(\frac{\Delta M}{10^{-9}{\, \mathrm{M}_{\odot }}} \right) \, {\, \mathrm{K}}, \end{eqnarray*} (24) where kB is the Boltzmann constant, C is the heat capacity, p is the pressure, Q is the heat released locally due the reactions, and ΔM is the accreted mass. Some of this heating will be converted into the quadrupolar temperature increase; however, it is quite unclear how much will be converted. Ushomirsky et al. (2000) estimate that δTq/δT ≲ 0.1, but, in reality, this ratio is poorly understood. These thermal mountains are built during accretion outbursts. During quiescence phases, the deformations are washed away on a thermal timescale (Brown, Bildsten & Rutledge 1998), \begin{eqnarray*} \tau _\text{th} \approx 0.2 \left(\frac{p}{10^{30}{\, \mathrm{erg}\, \mathrm{cm^{-3}}}} \right)^{3/4} \, \, \mathrm{yr}. \end{eqnarray*} (25) If the system is in quiescence for longer than this time-scale, then the thermal mountain will be washed away and a new mountain will be built during the next outburst. We implemented the expression for a quadrupole moment due to these reactions from equation (23) and assumed the following values for our NSs [as estimated by Haskell & Patruno (2017) for the pulsar J1023+0038]: C/baryon ≈ 10−6kB, |$E_\text{th} = {30}{\, \mathrm{MeV}}$|⁠, |$p = 10^{30}{\, \mathrm{erg}\, \mathrm{cm^{-3}}}$|⁠. and |$Q = {0.5}{\, \mathrm{MeV}}$|⁠. This meant that the quadrupole due to these reactions was dependent only on the accreted mass ΔM and the fraction δTq/δT. For this mechanism, we only considered transient accretion since in persistent accretion these mountains will not wash away, but, instead, will get progressively larger until the crust can no longer sustain them. This is effectively modelled through a fixed quadrupole that represents the largest mountain that can be built. In our model, we calculated ΔM by numerically integrating the accretion profile from the beginning of the outburst up to Ft. Our quiescence phases were long enough for the mountain to wash away during them. Unlike for our other prescriptions, we found that for thermal mountains, the specific values that parametrize the outburst features were very important in dictating the final-spin distribution. Based on observations of X-ray transients, we chose |$F_\text{t} = {0.1}{\, \mathrm{yr}}$|⁠, |$T_\text{recurrence} = {2.0}{\, \mathrm{yr}}$| and fmax/fmin = 104. We simulated NSs that built thermal mountains with δTq/δT = 4 × 10−4. The left-hand panel of Fig. 10 shows the resultant distribution of final spins. We can see qualitatively that this distribution has similar features to the fixed quadrupole case (see the right-hand panel of Fig. 6). Another promising aspect of the final-spin distribution is the prominence of the high-frequency peak. Like what is observed, this peak is narrow and much larger than the values in other frequency bins. We found a p-value of p = 0.60 for this distribution. Figure 10. View largeDownload slide Distributions of spin frequencies for simulated transiently accreting NSs that built thermal mountains during outburst phases with initial distributions from Table 2. The left-hand panel has a fixed δTq/δT = 4 × 10−4 and the right-hand panel has δTq/δT distributed flat-in-the-log between 10−4 and 10−2. Figure 10. View largeDownload slide Distributions of spin frequencies for simulated transiently accreting NSs that built thermal mountains during outburst phases with initial distributions from Table 2. The left-hand panel has a fixed δTq/δT = 4 × 10−4 and the right-hand panel has δTq/δT distributed flat-in-the-log between 10−4 and 10−2. The only constraint on the ratio of quadrupolar to total heating comes from the non-detection of X-ray emission of quadrupolar flux perturbations during quiescence phases in LMXBs, which gives δTq/δT ≲ 0.1 (Ushomirsky et al. 2000). Currently, there is no reason to believe that this fraction should be constant for all NSs. To account for our uncertainty in this fraction, we distributed δTq/δT flat-in-the-log between 10−4 and 10−2. The result of this simulation is shown in the right-hand panel of Fig. 10. The distribution peaks at low frequencies and then falls off towards higher frequencies. We obtained p = 7.0 × 10−2 for this case which meant that we could reject the null hypothesis. This shows how this prescription favours δTq/δT being a fixed value. Seeking a physical explanation for this preference of δTq/δT being a fixed value as opposed to being distributed is beyond the scope of this paper and has been left for future work. 5.3 Unstable r modes An r mode is a fluid mode of oscillation for which the restoring force is the Coriolis force. Andersson (1998) demonstrated that gravitational radiation destabilizes the r modes of rotating stars. These modes are generically unstable to GW emission (Friedman & Morsink 1998) and satisfy the Chandrasekhar–Friedman–Schutz instability, which facilitates the star finding lower energy and angular momentum configurations that allow the mode amplitude to grow (Chandrasekhar 1970; Friedman & Schutz 1978). The r-mode instability has long been considered a potential mechanism for imposing a spin limit on NSs in LMXBs (Andersson et al. 1999). The typical picture involves an NS being spun up through accretion until it enters the r-mode instability window. This instability region depends primarily on the spin of the NS and its core temperature. Once an NS has entered this region, it will emit gravitational radiation and begin to spin down until it reaches stability. This is expected to occur on a time-scale much shorter than the age of the system and should result in most LMXBs being stable. However, theoretical models for the r-mode instability demonstrate that many of the observed accreting NSs, in fact, lie inside the instability window (Ho, Andersson & Haskell 2011). This result would be consistent if the saturation amplitude for these systems was small, α ≈ 10−8–10−7, but this is at odds with predictions which suggest that the amplitude should be several orders of magnitude higher than this (Bondarescu et al. 2007). Owen et al. (1998) described a phenomenological model for the evolution of r modes and the spin of the star. In this model, the quadrupole moment for a constant-density NS that is unstable due to r modes is given by \begin{eqnarray*} Q_{2 2} \approx 1.67\times 10^{33} \left(\frac{\alpha }{10^{-7}} \right) M_{1.4} R_6^3 \left(\frac{P}{{1}{\, \mathrm{s}}}\right)^{-1} \, \mathrm{g\, cm^{2}}. \end{eqnarray*} (26) An interesting feature of this expression is its dependence on the spin of the NS. As a NS spins faster the quadrupole moment grows. This is different to what is expected from mountains. In fact, r modes and mountains could be differentiated from one another through the scaling of the quadrupoles as well as the frequency of the emitted GWs; for mountains the GW frequency is 2ν, whereas for r modes the frequency is 4ν/3. In order to simulate accreting NSs with unstable r modes, we implemented equation (26) into our model. We assumed that the mode-amplitude α remained constant for each NS. We repeated the previous simulations for persistent and transient accretors with unstable r modes and α = 10−7. Fig. 11 shows the final-spin distributions for those simulations. The unstable r modes were sufficient in both cases to give a peak at high spin frequencies. For the persistently accreting NSs, the peak was in the |$500\rm {-}550{\, \mathrm{Hz}}$| frequency bin, and for the transient accretors, the peak was in the |$550\rm {-}600{\, \mathrm{Hz}}$| bin. These distributions are similar to the case of a permanent quadrupole |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$| (see Fig. 6). For transient accretion with unstable r modes, the peak is narrower and more pronounced indicating that the magnitude of the GW torque sets in quickly. This is due to the scaling of the quadrupole in equation (26), since it depends linearly on the spin. For the persistent case, we found a p-value of p = 0.28 and for the transient case p = 0.57. Figure 11. View largeDownload slide Distributions of spin frequencies for simulated persistently accreting NSs (left-hand panel) and transiently accreting NSs (right-hand panel) with unstable r modes with initial distributions from Table 2 and α = 10−7. Figure 11. View largeDownload slide Distributions of spin frequencies for simulated persistently accreting NSs (left-hand panel) and transiently accreting NSs (right-hand panel) with unstable r modes with initial distributions from Table 2 and α = 10−7. We also conducted a simulation where α was distributed flat-in-the-log between 10−8 and 10−4. The result is shown in Fig. 12. We can see, similar to the thermal mountain distribution, that both distributions follow an exponentially decreasing behaviour. From these distributions, we found p = 1.2 × 10−2 when the NSs were persistently accreting and p = 1.5 × 10−3 when they were transiently accreting. From these p-values, we can reject the null hypothesis and note that the unstable-r-modes prescription produces more promising results when α is fixed, which is in agreement with current theoretical expectations (Arras et al. 2003; Bondarescu et al. 2007). Figure 12. View largeDownload slide Distributions of spin frequencies for simulated persistently accreting NSs (left-hand panel) and transiently accreting NSs (right-hand panel) with unstable r modes with initial distributions from Table 2 and α distributed flat-in-the-log between 10−8 and 10−4. Figure 12. View largeDownload slide Distributions of spin frequencies for simulated persistently accreting NSs (left-hand panel) and transiently accreting NSs (right-hand panel) with unstable r modes with initial distributions from Table 2 and α distributed flat-in-the-log between 10−8 and 10−4. 6 CONCLUSIONS An unresolved problem in the study of LMXBs is the unusual spin distribution of rapidly accreting NSs and, in particular, why no NS has been observed to spin close to the centrifugal break-up frequency. A potential explanation to this problem comes from GWs. Theoretically, GWs could be able to spin down these systems away from the break-up frequency. However, there are a number of different mechanisms that could give rise to gravitational radiation and it is unclear which are the most probable. It is also unclear whether GWs are the only way to explain the observed distribution of accreting NSs. For example, it has recently been suggested by Parfrey, Spitkovsky & Beloborodov (2016) that spin-down torques from an enhanced pulsar wind due to a disc-induced opening of the magnetic field could have a meaningful effect on the spin evolution of an accreting neutron star. Such a torque is not phenomenologically accounted for in our method and could be a direction for future work. In this paper, we have explored, within the context of our current understanding of accretion torques, whether an additional component is required in order to describe the spin evolution of accreting NSs. We investigated whether GW emission could be one such explanation and have compared competing GW mechanisms. We presented our model for the spin evolution of an accreting NS that accounts for accretion and magnetic-field effects, and also includes a GW spin-down component. Our model is able to simulate persistent and transient accretors. In our simulations with no GW torques we obtained NSs with much higher spins than what is observed. We did not obtain any of the characteristic behaviour of the observed spin distribution. In particular, there was no evidence of a pile-up at high frequencies. However, by adding a permanent quadrupole moment, motivated by torque balance, of |$Q_{2 2} = 10^{36}{\, \mathrm{g\, cm^{2}}}$| we obtained qualitatively similar behaviour to the observed distribution for the transiently accreting NS population. We considered the impact of magnetic-dipole radiation on the above results. We found that in the case of no GW emission one does not obtain the observed distribution. With the inclusion of GW emission, the resultant distribution is qualitatively similar to the case of no magnetic-dipole radiation. By varying the outburst duration with GW and magnetic-dipole torques, we obtained a distribution with a broad high-frequency peak. We investigated two GW-production prescriptions. For thermal mountains produced by asymmetric nuclear reactions in the crust, our model was sensitive to the precise features of the outburst profile, as well as the ratio of quadrupolar to total heating, δTq/δT. We found that a value of δTq/δT = 4 × 10−4 produced a similar distribution to what is observed. Promisingly, this gave the characteristic pile-up at high frequencies with a narrow, pronounced peak. We examined whether the distributions had a preference for δTq/δT being a single value or being distributed and found strong evidence arguing that it should be a fixed value. Accreting NSs with unstable rmodes and α = 10−7 produced similar results to the case with a fixed quadrupole moment and the thermal mountain prescription. This prescription favoured α being fixed as opposed to being distributed. The three cases that produced the most promising distributions that were qualitatively similar to the observed spin distribution – permanent quadrupole, thermal mountains, and unstable r modes – are almost indistinguishable from one another. Although, the r-mode instability case could, in theory, be differentiated from the other prescriptions. This distinction could come from the fact that the quadrupole moment due to unstable r modes scales linearly with the spin frequency of the star. Another key difference comes from the frequency of the GWs that are emitted through this channel. Unstable r modes emit GWs with a frequency of 4ν/3, whereas GWs due to deformations on an NS have a frequency of 2ν. For the r-mode scenario, the value for the saturation amplitude that we found to agree well with observation (α = 10−7) is many orders of magnitude below what is currently predicted. Theory would need to explain why this is so, or why the instability window is smaller than what is usually assumed. We have not addressed the spin distribution of RMSPs in this work. Future work could explore how the LMXB population evolves into the RMSP population and consider whether GWs are relevant in this process and can explain the observed distribution. In our modelling of transient accretion, we considered a simple fast-rise, exponential-decay function with a constant average accretion rate. However, in these systems it is expected that binary evolution will play a key role in the accretion rates and result in a long-term modulation of the accretion rate. This, of course, could have a significant effect on the resultant spin distribution. Such long-term variations could be explored in a future study. ACKNOWLEDGEMENTS The authors are grateful for useful discussions with B. Haskell and D. I. Jones. NA acknowledges support from the STFC via grant numbers ST/M000931/1 and ST/R00045X/1. REFERENCES Aasi J. et al. . , 2013 , Phys. Rev. D , 87 , 042001 https://doi.org/10.1103/PhysRevD.87.042001 Aasi J. et al. . , 2014 , ApJ , 785 , 119 https://doi.org/10.1088/0004-637X/785/2/119 Crossref Search ADS Aasi J. et al. . , 2015a , Phys. Rev. D , 91 , 022004 https://doi.org/10.1103/PhysRevD.91.022004 Aasi J. et al. . , 2015b , Phys. Rev. D , 91 , 062008 https://doi.org/10.1103/PhysRevD.91.062008 Abadie J. et al. . , 2011a , Phys. Rev. D , 83 , 042001 https://doi.org/10.1103/PhysRevD.83.042001 Abadie J. et al. . , 2011b , ApJ , 737 , 93 https://doi.org/10.1088/0004-637X/737/2/93 Crossref Search ADS Abadie J. et al. . , 2012 , Phys. Rev. D , 85 , 022001 https://doi.org/10.1103/PhysRevD.85.022001 Abbott B. et al. . , 2004 , Phys. Rev. D , 69 , 082004 https://doi.org/10.1103/PhysRevD.69.082004 Abbott B. et al. . , 2005a , Phys. Rev. D , 72 , 102004 https://doi.org/10.1103/PhysRevD.72.102004 Abbott B. et al. . , 2005b , Phys. Rev. Lett. , 94 , 181103 https://doi.org/10.1103/PhysRevLett.94.181103 Crossref Search ADS Abbott B. et al. . , 2007a , Phys. Rev. D , 76 , 042001 https://doi.org/10.1103/PhysRevD.76.042001 Abbott B. et al. . , 2007b , Phys. Rev. D , 76 , 082001 https://doi.org/10.1103/PhysRevD.76.082001 Abbott B. et al. . , 2008a , Phys. Rev. D , 77 , 022001 https://doi.org/10.1103/PhysRevD.77.022001 Abbott B. et al. . , 2008b , ApJ , 683 , L45 https://doi.org/10.1086/591526 Crossref Search ADS Abbott B. et al. . , 2009 , Phys. Rev. D , 79 , 022001 https://doi.org/10.1103/PhysRevD.79.022001 Abbott B. P. et al. . , 2010 , ApJ , 713 , 671 https://doi.org/10.1088/0004-637X/713/1/671 Crossref Search ADS Abbott B. P. et al. . , 2016 , Phys. Rev. D , 94 , 042002 https://doi.org/10.1103/PhysRevD.94.042002 Abbott B. P. et al. . , 2017a , Phys. Rev. D , 95 , 122003 https://doi.org/10.1103/PhysRevD.95.122003 Abbott B. P. et al. . , 2017b , Phys. Rev. D , 96 , 062002 https://doi.org/10.1103/PhysRevD.96.062002 Abbott B. P. et al. . , 2017c , Phys. Rev. D , 96 , 122006 https://doi.org/10.1103/PhysRevD.96.122006 Abbott B. P. et al. . , 2017d , ApJ , 839 , 12 https://doi.org/10.3847/1538-4357/aa677f Crossref Search ADS Abbott B. P. et al. . , 2017e , ApJ , 847 , 47 https://doi.org/10.3847/1538-4357/aa86f0 Crossref Search ADS Abbott B. P. et al. . , 2018a , Phys. Rev. D , 97 , 102003 https://doi.org/10.1103/PhysRevD.97.102003 Abbott B. P. et al. . , 2018b , Phys. Rev. Lett. , 120 , 031104 https://doi.org/10.1103/PhysRevLett.120.031104 Crossref Search ADS Alpar M. A. , Cheng A. F. , Ruderman M. A. , Shaham J. , 1982 , Nature , 300 , 728 https://doi.org/10.1038/300728a0 Crossref Search ADS Andersson N. , 1998 , ApJ , 502 , 708 https://doi.org/10.1086/305919 Crossref Search ADS Andersson N. , Kokkotas K. D. , Stergioulas N. , 1999 , ApJ , 516 , 307 https://doi.org/10.1086/307082 Crossref Search ADS Andersson N. , Jones D. I. , Kokkotas K. D. , Stergioulas N. , 2000 , ApJ , 534 , L75 https://doi.org/10.1086/312643 Crossref Search ADS Andersson N. , Jones D. I. , Kokkotas K. D. , 2002 , MNRAS , 337 , 1224 https://doi.org/10.1046/j.1365-8711.2002.05837.x Crossref Search ADS Andersson N. , Glampedakis K. , Haskell B. , Watts A. L. , 2005 , MNRAS , 361 , 1153 https://doi.org/10.1111/j.1365-2966.2005.09167.x Crossref Search ADS Arras P. , Flanagan E. E. , Morsink S. M. , Schenk A. K. , Teukolsky S. A. , Wasserman I. , 2003 , ApJ , 591 , 1129 https://doi.org/10.1086/374657 Crossref Search ADS Bhattacharyya S. , Chakrabarty D. , 2017 , ApJ , 835 , 4 https://doi.org/10.3847/1538-4357/835/1/4 Crossref Search ADS Bildsten L. , 1998 , ApJ , 501 , L89 https://doi.org/10.1086/311440 Crossref Search ADS Bondarescu R. , Teukolsky S. A. , Wasserman I. , 2007 , Phys. Rev. D , 76 , 064019 https://doi.org/10.1103/PhysRevD.76.064019 Brown E. F. , Bildsten L. , Rutledge R. E. , 1998 , ApJ , 504 , L95 https://doi.org/10.1086/311578 Crossref Search ADS Chakrabarty D. , Morgan E. H. , Muno M. P. , Galloway D. K. , Wijnands R. , van der Klis M. , Markwardt C. B. , 2003 , Nature , 424 , 42 https://doi.org/10.1038/nature01732 Crossref Search ADS PubMed Chandrasekhar S. , 1970 , Phys. Rev. Lett. , 24 , 611 https://doi.org/10.1103/PhysRevLett.24.611 Crossref Search ADS Cook G. B. , Shapiro S. L. , Teukolsky S. A. , 1994 , ApJ , 424 , 823 https://doi.org/10.1086/173934 Crossref Search ADS Cutler C. , 2002 , Phys. Rev. D , 66 , 084025 https://doi.org/10.1103/PhysRevD.66.084025 D’Angelo C. R. , 2017 , MNRAS , 470 , 3316 https://doi.org/10.1093/mnras/stx1306 Crossref Search ADS Done C. , Gierliński M. , Kubota A. , 2007 , A&A Rev. , 15 , 1 https://doi.org/10.1007/s00159-007-0006-1 Friedman J. L. , Morsink S. M. , 1998 , ApJ , 502 , 714 https://doi.org/10.1086/305920 Crossref Search ADS Friedman J. L. , Schutz B. F. , 1978 , ApJ , 222 , 281 https://doi.org/10.1086/156143 Crossref Search ADS Ghosh P. , Lamb F. K. , 1978 , ApJ , 223 , L83 https://doi.org/10.1086/182734 Crossref Search ADS Ghosh P. , Lamb F. K. , 1979 , ApJ , 234 , 296 https://doi.org/10.1086/157498 Crossref Search ADS Haensel P. , Zdunik J. L. , 1990 , A&A , 227 , 431 Haskell B. , Patruno A. , 2017 , Phys. Rev. Lett. , 119 , 161103 https://doi.org/10.1103/PhysRevLett.119.161103 Crossref Search ADS PubMed Haskell B. , Jones D. I. , Andersson N. , 2006 , MNRAS , 373 , 1423 https://doi.org/10.1111/j.1365-2966.2006.10998.x Crossref Search ADS Haskell B. , Samuelsson L. , Glampedakis K. , Andersson N. , 2008 , MNRAS , 385 , 531 https://doi.org/10.1111/j.1365-2966.2008.12861.x Crossref Search ADS Haskell B. , Zdunik J. L. , Fortin M. , Bejger M. , Wijnands R. , Patruno A. , 2018 , A&A , 620 , A69 https://doi.org/10.1051/0004-6361/201833521 Crossref Search ADS Hessels J. W. T. , 2008 , in Wijnands R. , Altamirano D. , Soleri P. , Degenaar N. , Rea N. , Casella P. , Patruno A. , Linares M. , eds, AIP Conf. Proc. Vol. 1068, A Decade of Accreting Millisecond X-ray Pulsars . Am. Inst. Phys , New York , p. 130 Google Preview WorldCat Hessels J. W. T. , Ransom S. M. , Stairs I. H. , Freire P. C. C. , Kaspi V. M. , Camilo F. , 2006 , Science , 311 , 1901 https://doi.org/10.1126/science.1123430 Crossref Search ADS PubMed Heyl J. S. , 2002 , ApJ , 574 , L57 https://doi.org/10.1086/342263 Crossref Search ADS Ho W. C. G. , Andersson N. , Haskell B. , 2011 , Phys. Rev. Lett. , 107 , 101101 https://doi.org/10.1103/PhysRevLett.107.101101 Crossref Search ADS PubMed Ho W. C. G. , Klus H. , Coe M. J. , Andersson N. , 2014 , MNRAS , 437 , 3664 https://doi.org/10.1093/mnras/stt2193 Crossref Search ADS Johnson-McDaniel N. K. , Owen B. J. , 2013 , Phys. Rev. D , 88 , 044004 https://doi.org/10.1103/PhysRevD.88.044004 Crossref Search ADS Lasota J. P. , 1997 , in Wickramasinghe D. T. , Bicknell G. V. , Ferrario L. , eds, ASP Conf. Ser. Vol. 121, IAU Colloq. 163: Accretion Phenomena and Related Outflows . Astron. Soc. Pac , San Francisco , p. 351 Google Preview WorldCat Lattimer J. M. , Prakash M. , 2007 , Phys. Rep. , 442 , 109 https://doi.org/10.1016/j.physrep.2007.02.003 Levin Y. , 1999 , ApJ , 517 , 328 https://doi.org/10.1086/307196 Crossref Search ADS Melatos A. , Payne D. J. B. , 2005 , ApJ , 623 , 1044 https://doi.org/10.1086/428600 Crossref Search ADS Nayyar M. , Owen B. J. , 2006 , Phys. Rev. D , 73 , 084001 https://doi.org/10.1103/PhysRevD.73.084001 Owen B. J. , 2005 , Phys. Rev. Lett. , 95 , 211101 https://doi.org/10.1103/PhysRevLett.95.211101 Crossref Search ADS PubMed Owen B. J. , Lindblom L. , Cutler C. , Schutz B. F. , Vecchio A. , Andersson N. , 1998 , Phys. Rev. D , 58 , 084020 https://doi.org/10.1103/PhysRevD.58.084020 Papaloizou J. , Pringle J. E. , 1978 , MNRAS , 184 , 501 https://doi.org/10.1093/mnras/184.3.501 Crossref Search ADS Papitto A. , Torres D. F. , Rea N. , Tauris T. M. , 2014 , A&A , 566 , A64 https://doi.org/10.1051/0004-6361/201321724 Crossref Search ADS Parfrey K. , Spitkovsky A. , Beloborodov A. M. , 2016 , ApJ , 822 , 33 https://doi.org/10.3847/0004-637X/822/1/33 Crossref Search ADS Patruno A. , 2010 , ApJ , 722 , 909 https://doi.org/10.1088/0004-637X/722/1/909 Crossref Search ADS Patruno A. , Watts A. L. , 2012 , preprint (arXiv:1206.2727) Patruno A. , Haskell B. , D’Angelo C. , 2012 , ApJ , 746 , 9 https://doi.org/10.1088/0004-637X/746/1/9 Crossref Search ADS Patruno A. , Haskell B. , Andersson N. , 2017 , ApJ , 850 , 106 https://doi.org/10.3847/1538-4357/aa927a Crossref Search ADS Payne D. J. B. , Melatos A. , 2006 , ApJ , 641 , 471 https://doi.org/10.1086/498855 Crossref Search ADS Possenti A. , Colpi M. , D’Amico N. , Burderi L. , 1998 , ApJ , 497 , L97 https://doi.org/10.1086/311286 Crossref Search ADS Pringle J. E. , Rees M. J. , 1972 , A&A , 21 , 1 Priymak M. , Melatos A. , Payne D. J. B. , 2011 , MNRAS , 417 , 2696 https://doi.org/10.1111/j.1365-2966.2011.19431.x Crossref Search ADS Radhakrishnan V. , Srinivasan G. , 1982 , Curr. Sci , 51 , 1096 Spitkovsky A. , 2006 , ApJ , 648 , L51 https://doi.org/10.1086/507518 Crossref Search ADS Ushomirsky G. , Rutledge R. E. , 2001 , MNRAS , 325 , 1157 https://doi.org/10.1046/j.1365-8711.2001.04515.x Crossref Search ADS Ushomirsky G. , Cutler C. , Bildsten L. , 2000 , MNRAS , 319 , 902 https://doi.org/10.1046/j.1365-8711.2000.03938.x Crossref Search ADS Wagoner R. V. , 1984 , ApJ , 278 , 345 https://doi.org/10.1086/161798 Crossref Search ADS Wagoner R. V. , 2002 , ApJ , 578 , L63 https://doi.org/10.1086/344502 Crossref Search ADS Wang Y.-M. , 1995 , ApJ , 449 , L153 https://doi.org/10.1086/309649 Wang Y.-M. , 1996 , ApJ , 465 , L111 https://doi.org/10.1086/310150 Crossref Search ADS Watts A. L. , Krishnan B. , Bildsten L. , Schutz B. F. , 2008 , MNRAS , 389 , 839 https://doi.org/10.1111/j.1365-2966.2008.13594.x Crossref Search ADS White N. E. , Zhang W. , 1997 , ApJ , 490 , L87 https://doi.org/10.1086/311018 Crossref Search ADS Woan G. , Pitkin M. D. , Haskell B. , Jones D. I. , Lasky P. D. , 2018 , ApJ , 863 , L40 https://doi.org/10.3847/2041-8213/aad86a Crossref Search ADS © 2019 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) TI - Population synthesis of accreting neutron stars emitting gravitational waves JF - Monthly Notices of the Royal Astronomical Society DO - 10.1093/mnras/stz1719 DA - 2019-09-01 UR - https://www.deepdyve.com/lp/oxford-university-press/population-synthesis-of-accreting-neutron-stars-emitting-gravitational-1DwyfmQz3c SP - 99 VL - 488 IS - 1 DP - DeepDyve ER -