Starquakes in millisecond pulsars and gravitational waves emission
MMNRAS , 1–12 (2020) Preprint 5 February 2021 Compiled using MNRAS L A TEX style file v3.0
Starquakes in millisecond pulsars and gravitational waves emission
E. Giliberti , (cid:63) , G. Cambiotti Dipartimento di Fisica, Universit`a degli Studi di Milano, Via Celoria 16, 20133,Milano, Italy Dipartimento di Scienze della Terra, Universit`a degli Studi di Milano, Via Cicognara 7, Milano, 20129, Italy Istituto Nazionale di Fisica Nucleare, sezione di Milano, Via Celoria 16, 20133 Milano, Italy
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
So far, only transient Gravitational waves (GWs) produced by catastrophic eventsof extra-galactic origin have been detected. However, it is generally believed thatthere should be also continuous sources of GWs within our galaxy, such as accretingneutron stars (NSs). In fact, in accreting NSs, centrifugal forces can be so strongto break the neutron star crust (causing a starquake), thus producing a quadrupolemoment responsible for the continuous emission of GWs. At equilibrium, the angularmomentum gained by accretion and lost via GWs emission should balance each other,stopping the stellar spin-up.We hereinafter investigate the above physical picture within the framework of aNewtonian model describing compressible, non-magnetized and self-gravitating NSs.In particular, we calculate the rotational frequency need to break the stellar crust ofan accreting pulsar and we estimate the upper limit for the ellipticity due to this event.Depending on the equation of state (EoS) and on the mass of the star, we calculatedthat the starquake-induced ellipticity ranges from 10 − to 10 − . The correspondingequilibrium frequency that we find is in good agreement with observations and, for allthe scenarios, it is below the observational limit frequency of 716 .
36 Hz. Finally, wealso discuss possible observational constraints on the ellipticity upper limit of accretingpulsars.
Key words: star: neutron – gravitational waves –
Gravitational waves (GWs) detections have widen ourknowledge of astrophysical events. First the discover of blackholes merger (Abbott 2016), and then the neutron stars(NSs) coalescence (Abbott 2017) have opened new possiblewindows for the study of extreme compact objects. How-ever, until now direct detections of GWs came only fromcatastrophic, transient events, with an extra-galactic origin.We expect, nonetheless, that continuous signals should comealso from our Galaxy, emitted by fast rotating, accreting pul-sars.Observation of Low-Mass X-ray Binaries (LMXBs) hasshown a paucity of stars rotating near the centrifugal break-up frequency (Lattimer & Prakash 2007; Chakrabarty2008), opening the outstanding question of why these ob- (cid:63) [email protected] This very rough estimate of the frequency beyond whichthe centrifugal force tears apart the star is given by the Ke-plerian rotational frequency ν limit = (cid:112) GM/ πa , where M (cid:39) . M odot , a (cid:39) km are the canonical stellar mass and radius,respectively, and G is the gravitational constant. jects seems to spin well under that limit. In fact, prelim-inary estimates by Cook et al. (1994) suggested that thespin-up timescale for NSs in LMXBs should be large enoughto make them reach at least the rotational frequency of 1kHz. On the contrary, the actual fastest spinning accretingNS has ν = 599 Hz (Galloway et al. 2005), and also mil-lisecond pulsars, that are thought to be the ultimate fate ofLMXBs (Bhattacharya & van den Heuvel 1991), rotate witha frequency lower than the break-up frequency. In particu-lar, Chakrabarty et al. (2003) using Bayesian statistics haveshown that the actual distribution of the Accreting Millisec-ond X-ray Pulsars (AMXPs) gives a theoretical maximumspin frequency of 760 Hz.One possible explanation for this behaviour is that thesekind of pulsars emit GWs that make them slow down. Inparticular, many works (Bildsten 1998; Ushomirsky et al.2000; Watts et al. 2008) suggest that accreting millisecondpulsars can reach an equilibrium configuration when the an-gular momentum gained from infalling material is lost byGWs emission. Papaloizou & Pringle (1978) and Wagoner(1984) firstly suggested that GWs emission can explain thedynamical equilibrium of NSs, but only more recently Bild-sten (1998) proposed the possibility of mountains forming © a r X i v : . [ a s t r o - ph . H E ] F e b Giliberti et al. on accreting objects as a concrete mechanism for generatinga non-zero ellipticity and, as a consequence, GWs.It is also interesting to remind that multi-million molec-ular dynamic simulations have shown that the crustal break-ing strain can be quite large (Horowitz & Kadau 2009; Baiko& Chugunov 2018) and, therefore, that the crust can sus-tain a maximum ellipticity large enough to generate GWsdetectable from Earth by the current generation of inter-ferometers (Haskell et al. 2006; Johnson-McDaniel & Owen2013).In the literature, two main mechanisms able to producea static ellipticity have been studied: thermal mountains(Bildsten 1998; Ushomirsky et al. 2000; Haskell et al. 2015),and magnetically confined mountains (Cutler 2002; Melatos& Payne 2005; Haskell et al. 2008; Vigelius & Melatos 2009;Priymak et al. 2011). The former are due to pycno-nuclearreactions that heat the accreted material deep into the crust.The latter, instead, are caused by a local enhancements ofthe magnetic field structure (related to accretion) that cansustain mountains.Using the higher breaking strain threshold, and model-ing a NS as a homogeneous, incompressible object, Fattoyevet al. (2018) claimed that starquakes (Ruderman (1969),Baym et al. (1969)) can happen only on accreting, rapidlyrotating star, where the centrifugal force is large enough tomake the crust reach the failure threshold. They also intro-duced the hypothesis that the breaking of the crust mightproduce a quadrupolar deformation sufficient to emit enoughenergy through GWs to prevent the stellar spin-up. How-ever, in Fattoyev et al. (2018) the ellipticity produced bystarquakes is not self-consistently calculated as well as theresultant evolutionary path of the NS, i.e. the reaching ofthe equilibrium angular velocity.In this work we use the model described in Gilibertiet al. (2020) to study the deformation of a rotating, com-pressible, non-magnetized, self-gravitating NS to explore theidea that a sequence of starquakes can act as a trigger forGWs emission.In particular, we will study the following physical pic-ture. A NS is accreting mass from a companion: the infallingmaterial creates a disk that transfers angular momentum tothe star and spins-up the central object. The NS will thusrotate faster and faster til the breaking condition is reached:in that moment a first starquake occurs, altering the stel-lar axial symmetry by creating a non-null ellipticity. Fromthen on, the star radiates GWs. The balance between theangular momentum gained from accretion and the one lostby emission will bring the star through a sequence of breaksand finally to a dynamical equilibrium frequency.The paper’s organization reflects the different stepsneeded to explore the NS’s evolutionary path: first of allin Section 2 we summarize the model described in Gilibertiet al. (2020), focusing on the stellar configuration; in section3 we study the problem of crust failure and find an estima-tion of what are the typical frequencies necessary to breakthe crust. Once it is shown that fast rotating stars can reachthe breakup frequency, we move forward in section 4 by in-troducing an upper limit for the ellipticity (cid:15) due to a seriesof starquakes on the NS. In this way, in section 5 we will beready to study the dynamical equilibrium frequency. Finally,in section 6 we make a comparison of the ellipticity upper
Figure 1.
Sketch of the stellar structure for the model consideredin this work (not in scale). The star is divided into a fluid core,extending from the centre to r = r c , and an elastic crust, thatgoes from the core-crust boundary up to the stellar surface r = a . limit predicted by our model with the one estimated usingobservations of fast rotating pulsars. We want to analyze the stressing effect of material accret-ing from a companion on a NS, accelerating the object. Ouraim is to study the deformation of a self-gravitating, non-magnetized , compressible NS, under the effect of the cen-trifugal force, in order to find what is the maximum rota-tional frequency before crust-breaking. For this purpose, weuse a general Newtonian model (Giliberti et al. 2020), wherethe star consists of a fluid core, extending from the originto the radius r c , and an elastic crust, that covers the regionfrom r c to the stellar surface r = a (Fig. 1). The outer-crustboundary is placed at the density 1 × g/cm in order toguarantee the numerical stability of the solution against thecomputational problems due to the very rapid variation ofthe density in the outermost layers (Ushomirsky et al. 2000).Following Ushomirsky et al. (2000), the crust-core transitionis set at the fiducial density 1 . × g/cm , that impliesa core-crust transition at r c ≈ . a for a standard neu-tron star with M = 1 . M (cid:12) . We use MATHEMATICA 11 toperform numerical computation.For describing the stellar matter, we choose the equa-tion of state (EoS) of a polytrope of index n = 1, since it In this work we are interested mainly on LMXB and millisecondpulsars, that have both typically very low magnetic field at thesurface B ≈ − G (Manchester et al. 2005): we expect that inthis condition B has a very small impact on crustal deformation(Franco et al. 2000). In a spherical coordinate system r is the radial distance fromthe star’s center, θ is the colatitude and ϕ is the longitude.MNRAS000
Sketch of the stellar structure for the model consideredin this work (not in scale). The star is divided into a fluid core,extending from the centre to r = r c , and an elastic crust, thatgoes from the core-crust boundary up to the stellar surface r = a . limit predicted by our model with the one estimated usingobservations of fast rotating pulsars. We want to analyze the stressing effect of material accret-ing from a companion on a NS, accelerating the object. Ouraim is to study the deformation of a self-gravitating, non-magnetized , compressible NS, under the effect of the cen-trifugal force, in order to find what is the maximum rota-tional frequency before crust-breaking. For this purpose, weuse a general Newtonian model (Giliberti et al. 2020), wherethe star consists of a fluid core, extending from the originto the radius r c , and an elastic crust, that covers the regionfrom r c to the stellar surface r = a (Fig. 1). The outer-crustboundary is placed at the density 1 × g/cm in order toguarantee the numerical stability of the solution against thecomputational problems due to the very rapid variation ofthe density in the outermost layers (Ushomirsky et al. 2000).Following Ushomirsky et al. (2000), the crust-core transitionis set at the fiducial density 1 . × g/cm , that impliesa core-crust transition at r c ≈ . a for a standard neu-tron star with M = 1 . M (cid:12) . We use MATHEMATICA 11 toperform numerical computation.For describing the stellar matter, we choose the equa-tion of state (EoS) of a polytrope of index n = 1, since it In this work we are interested mainly on LMXB and millisecondpulsars, that have both typically very low magnetic field at thesurface B ≈ − G (Manchester et al. 2005): we expect that inthis condition B has a very small impact on crustal deformation(Franco et al. 2000). In a spherical coordinate system r is the radial distance fromthe star’s center, θ is the colatitude and ϕ is the longitude.MNRAS000 , 1–12 (2020) tarquakes and gravitational waves allows us to study in a somewhat realistic way the star’sphysical characteristic as a function of its mass (for a fixedmass, a Newtonian approach gives larger stars with respectto Relativistic ones). With this choice, the stellar mass andradius are independent of each other, giving us the oppor-tunity to associate the Relativistic radius a with a givenmass M . In other words, we can use the realistic mass-radiusrelation of a chosen EoS, obtained from the integration ofthe Tolman-Oppenheimer-Volkoff equations, for fixing thestar’s radius a once the stellar mass M is chosen. In par-ticular, in this paper we will use two different EoSs andtheir mass-radius relation: SLy (Douchin & Haensel 2001)and the stiffer BSk21 (Goriely et al. 2010). Both are unified EoSs, covering consistently both the core and the crust; fur-thermore, the SLy EoS is a standard choice for making pre-dictions about crustal quadrupoles, see e.g. (Horowitz 2010;Johnson-McDaniel & Owen 2013). The comparison betweenEoSs is extremely important, since it makes possible to seethe impact of different stars configuration on the deforma-tion of NSs.Furthermore, following the analysis of Giliberti et al.(2020), we study the behaviour of the same star in two sce-narios, parametrized by different values of the adiabatic in-dex. In the first, the response of the stellar matter to pertur-bations is characterized by its equilibrium adiabatic index,namely γ ∗ = 2. This is the situation in which the typical dy-namical timescale are small compared to the reactions ones.In the second scenario, instead, the contrary is true and theadiabatic index differ from its equilibrium value. We referto it as the frozen adiabatic index γ f and consider the twocases of γ f = 2 . γ f = ∞ . We choose this two val-ues to give a preliminary insight of the spectrum of all thepossible allowed frozen indices: γ f = 2 . γ f = ∞ representsan incompressible response of an initial compressible starto external forces. We refer to Giliberti et al. (2020) for adeep discussion on the difference between these two scenar-ios and their impact on displacements and stresses of thestar’s crust.In our model, once the EoS and the value of the adi-abatic index have been chosen, there are only other twoparameters needed to complete the description of the NS’sconfiguration. They are the stellar elastic moduli: the bulkmodulus κ and the shear modulus µ . The first is given bythe relation κ ( r ) = γ P ( r ) , (1)where P is the local pressure. In particular, the initial nonrotating configuration will be always characterized by κ = γ ∗ P , while the response of the star to the centrifugal forcewill be modeled with κ = γ ∗ P or κ = γ f P , according to thedynamical timescale of the spin-up. Concerning the shearmodulus we use the same prescription as Cutler et al. (2003) µ = 10 − P. (2) Starting from an unstressed, spherical, symmetric, non-rotating configuration , the spin-up caused by the infallingmaterial coming by the companion will cause a deformationof the star, that will become oblate. The stress in the crustwill grow as the rotational velocity increases, till a breakingcondition is reached. We use the Tresca criterion to estab-lish when the crust will break, i.e. when we have a starquake.This criterion states that the crust will fail when the strainangle α , defined as the difference between the maximum andthe minimum eigenvalues of the strain tensor, is half of thebreaking strain σ max (Christensen 2013) α = σ max . (3)Since the deformation due to rotation is proportional to thefrequency squared α = ˜ αν , (4)the Tresca criterion will be satisfied at a threshold frequency,that we call breaking frequency ν b . ˜ α is a term dependingonly on the structure of the NS, namely its mass and EoS.From Eqs (3) and (4) we immediately get ν b = (cid:114) σ max α max , (5)where ˜ α max is the maximum values of ˜ α ( r, θ ) on the wholeNS crust, i.e. in the range r c ≤ r ≤ a , 0 ≤ θ ≤ π ( θ is thecolatitude angle). Therefore, the larger σ max , the larger thebreaking frequency, as expected.Unfortunately the value of σ max is very uncertain, rang-ing from 10 − of the first theoretical estimation of Ruderman(1991), to the more recent values obtained with moleculardynamic simulations of 10 − by Horowitz & Kadau (2009) orof 0 .
04 proposed with semi-analytical approaches by Baiko& Chugunov (2018). In the present work we use the largerbreaking strain 10 − for two reasons. The first is to compareour results with the ones obtained by Fattoyev et al. (2018),that used the same threshold. The second is to give an up-per limit for the breaking frequency ν b . In fact, in the caseof the lowest estimation of σ max = 10 − , we can see (cf Eq(5)) that the breaking frequency will decrease of about twoorders of magnitude compared to our choice.The results of this first analysis, coming from the modelbriefly introduced in section 2, are shown in Fig 2, wherethe curves for ν b are plotted for different EoSs and differentvalues of the adiabatic index. The main features of Fig 2 arethe following:(i) For a given NS mass, a softer EoS produces less com-pact stars and, so, gives larger breaking frequency values. In-deed, as shown by Giliberti et al. (2019) and Giliberti et al.(2020), the star’s compactness, M/a , is a key parameterscontrolling the deformations. We are interested in calculating the displacement field betweena configuration rotating with velocity Ω and one rotating at Ω + δ Ω, where δ Ω > δ Ω) − Ω .MNRAS , 1–12 (2020) Giliberti et al. (ii) The larger the adiabatic index value, the lower thebreaking frequency. We know (Giliberti et al. 2020) that in-compressible star ( γ f = ∞ ) will develop larger strains ˜ α max with respect to a compressible one: using Eq (5) it is clearthat this means that the larger γ the smaller the breakingfrequency.(iii) A small change in the adiabatic index value giveslarge changes in the breaking frequency curve. As an ex-ample, for the BSk21 EoS, we can see that the ν b curvefor γ f = 2 . γ ∗ = 2 and the other for γ f = ∞ . This is a typical fea-tures of compressible, self-gravitating NSs, due to the small-ness of the shear modulus with respect to the bulk modulus(Chamel & Haensel 2008). In fact, if the ratio µ/κ is small,the star, despite its elastic crust, behaves essentially like afluid. This means that the incompressible limit is reach evenwith a small departure of the adiabatic index from its equi-librium (Giliberti et al. 2020).(iv) For a typical M = 1 . M (cid:12) NS, we can say that thebreaking frequency is in the range 200 −
600 Hz, well belowthe maximum observed rotational frequency (Hessels et al.2006) ν o = 716 .
36 Hz . (6)In this respect, our analysis refines the results obtained byFattoyev et al. (2018), which predict larger breaking frequen-cies in the range 400 − ν o . The reasonswhy so far we have not observed any NS with this frequencycould be many and different, and are all compatible andunderstandable within our model. Indeed we expect veryfew NSs rotate with a frequency larger than ν o . In fact, themajority of NSs show a breaking frequency that is smallerthan 700 Hz and only massive stars with a softer EoS canreach a breaking frequency of 800 −
900 Hz.Moreover, we can think that the fast NSs we can ob-serve are actually near their equilibrium frequency , i.e. thefrequency at which the angular momentum gained from ac-cretion is equal to the one lost by GWs emission, that mustnot be confused with the breaking frequency , which is thetypical frequency at which the crust starts to fail. And, as wewill see in the following sections, the equilibrium frequencyis typically smaller than the breaking one.
Now that we know that the sufficiently fast rotating NSs canreach the condition for crust breaks, we can study the possi-ble consequence of starquake events and, as a consequence,the creation of a non-null ellipticity. In fact, a break of thecrust cause a local deformation, that brings the star awayfrom the pure axyal simmetry, see Fig. 3. This deformationcan be evaluated by the ellipticity (cid:15) , defined as (cid:15) = A − BC , (7)where A and B are the principal moments of inertia alongthe equatorial axes, while C is the one along the rotationaxis. The exact value of (cid:15) depends on the crustal propertiesof the star and on its seismic history. As a first approach, one M / M ⊙ ν b ( H z ) Figure 2.
The breaking frequency ν b as a function of the stellarmass for different EoS and adiabatic indices. The solid lines areobtained with the SLy EoS, while the dashed ones for the BSk21one. Different colours indicates different adiabatic indices values:red for γ ∗ = 2, blue for γ f = 2 . γ f = ∞ . Theblack, dotted line indicates the rotation frequency of the fastestrotating pulsar observed ν o = 716 .
36 Hz (Hessels et al. 2006). would be led to estimate the effect of every single starquake;however this is a too difficult task. In fact, our knowledgeof the NS crust physics is extremely poor, and a reasonabledescription of a quake involves a very large number of un-known parameters (dip and strike angles, displacement dis-continuity, fault area etc.). For these reason it seems morereasonable to follow a different approach.We look at starquakes as the attempts of a stressed starto achieve the equilibrium fluid shape despite the constrain-ing action of its elastic crust. As Eq (7) shows, in the simplecase of an uniform rotation the star is axially symmetric,i.e. A = B and thus (cid:15) = 0. Whatever the angular velocity,the crust of the star is stressed, since it cannot achieve thecorresponding equilibrium configuration that it would haveif it was completely fluid . The elastic crust, in fact, con-strains the star to have a more prolate shape with respectto the fluid one. However, the axial symmetry of centrifu-gal deformation can be broken by starquakes, that can cre-ate a mountain on the NS surface. Starting from an initial(pre-starquake) configuration, through a sufficient numberof breaking events, the star will, therefore, tend towards itsfluid configuration. In other words, the cumulative effect ofa sequence of many crust failures is to give a more oblateshape to the star. Clearly, it sound physically reasonable tostate that the dynamical equilibrium will be achieved witha sequence of events, since we do not expect that a singlequake could release all the stresses of the crust and bringinstantaneously the star to its fluid configuration. On theother hand, the breaking of the crust leads the NS to get anellipticity different from zero and thus, to emit GWs.Thus, in order to calculate the maximum ellipticity (i.e.the one reached after a “complete” sequence of starquakes)that a NS could have at a given angular velocity, we comparethe principal moment of inertia of two different configura-tions, rotating at the same frequency. The first is the one ofa NS with a solid crust, while the second is the one of a purefluid star. The difference between these two configuration, in MNRAS000
36 Hz (Hessels et al. 2006). would be led to estimate the effect of every single starquake;however this is a too difficult task. In fact, our knowledgeof the NS crust physics is extremely poor, and a reasonabledescription of a quake involves a very large number of un-known parameters (dip and strike angles, displacement dis-continuity, fault area etc.). For these reason it seems morereasonable to follow a different approach.We look at starquakes as the attempts of a stressed starto achieve the equilibrium fluid shape despite the constrain-ing action of its elastic crust. As Eq (7) shows, in the simplecase of an uniform rotation the star is axially symmetric,i.e. A = B and thus (cid:15) = 0. Whatever the angular velocity,the crust of the star is stressed, since it cannot achieve thecorresponding equilibrium configuration that it would haveif it was completely fluid . The elastic crust, in fact, con-strains the star to have a more prolate shape with respectto the fluid one. However, the axial symmetry of centrifu-gal deformation can be broken by starquakes, that can cre-ate a mountain on the NS surface. Starting from an initial(pre-starquake) configuration, through a sufficient numberof breaking events, the star will, therefore, tend towards itsfluid configuration. In other words, the cumulative effect ofa sequence of many crust failures is to give a more oblateshape to the star. Clearly, it sound physically reasonable tostate that the dynamical equilibrium will be achieved witha sequence of events, since we do not expect that a singlequake could release all the stresses of the crust and bringinstantaneously the star to its fluid configuration. On theother hand, the breaking of the crust leads the NS to get anellipticity different from zero and thus, to emit GWs.Thus, in order to calculate the maximum ellipticity (i.e.the one reached after a “complete” sequence of starquakes)that a NS could have at a given angular velocity, we comparethe principal moment of inertia of two different configura-tions, rotating at the same frequency. The first is the one ofa NS with a solid crust, while the second is the one of a purefluid star. The difference between these two configuration, in MNRAS000 , 1–12 (2020) tarquakes and gravitational waves terms of moment of inertia, will give us the maximum valueof (cid:15) . Let us introduce how to use our model for calculating theNS’s inertia tensor. In the following, we adopt the notationof Sabadini et al. (2016); for a brief summary of the no-tation used, main quantities and equations introduced, seeAppendix A.Consider an initially non-rotating star that is spun upto a given frequency ν . The centrifugal perturbation actingon the object involves both the (cid:96) = 0 and (cid:96) = 2 sphericalharmonics; however, since we are interested into the calcula-tion of the stellar ellipticity, we can focus only on the latter.In fact, the (cid:96) = 0 term would give the same contributions toall the principal moments of inertia and these contributionscancel each other in the difference A − B . Furthermore, asfar as the difference A − B is small, the contribution to C can be neglected within a first-order perturbation theory.Therefore, at the first order approximation Eq (7) can bewritten as: (cid:15) = ∆ A − ∆ BI , (8)where ∆ A and ∆ B are the changes of the principal momentsof inertia A and B due to (cid:96) = 2 spherical harmonic pertur-bations, and I is the unstressed stellar moment of inertia.In the Cartesian reference frame, the changes of the in-ertia tensor I ij due to (cid:96) = 2 spherical harmonic perturbationcan be obtained according to the following expression∆ I ij = (cid:88) m = − Q mij π (cid:90) a ρ ∆2 m ( r ) r dr, (9)where ρ ∆2 m and Q mij are, respectively, the spherical harmoniccoefficients of degree (cid:96) = 2 and order m of the density dis-tribution and of the matrix Q ij defined by Q ij = 13 δ ij − ˆ r i ˆ r j . (10)Here δ ij is the Kronecker delta and ˆ r i are the Cartesiancomponents of the radial unit vector. The deformation ofthe star is described by two Poisson equations, one for theperturbed gravitational potential φ ∆ ∇ φ ∆ = 4 πGρ ∆ , (11)and the other for the centrifugal potential φ C φ C = − π ν . (12)It is quite natural to introduce the total perturbed potentialΦ ∆ , as the sum of φ ∆ and φ C . By expanding also the to-tal perturbed potential in spherical harmonics we can write(Chao & Gross 1987)Φ ∆2 m ( a ) = − πG a (cid:90) a ρ ∆2 m ( r ) r dr. (13)Therefore, from Eqs (9) and (13), it follows that the per-turbed tensor of inertia can be written as:∆ I ij = − (cid:88) m = − Q mij a G Φ ∆2 m ( a ) . (14) Figure 3.
Sketch of the stellar deformation (not in scale). The staris divided into a fluid core (yellow) and an elastic crust (blue).
Top : Vision of the star from the three axis x, y, z in the uniformrotation configuration.
Bottom : Vision of the star from the threeprincipal axis x, y, z after a starquake that has created a mountainin the x direction. The above expression is extremely useful for the calcula-tion of the inertia changes, since it involves only the valueof the total potential at the star’s surface, that can be eas-ily obtained within the model presented in Giliberti et al.(2020).We also observe that in the case of a uniform rota-tion the inertia tensor can be expressed in a diagonal form,namely ∆ I = Diag[∆ A, ∆ A, ∆ C ] , (15)where ∆ C is the change of the moments of inertia alongthe rotational axis C . Finally, we note that by choosing acoordinate system in which the rotational axis z coincideswith θ = 0, we can restrict ourselves to the only (cid:96) = 2, m = 0 spherical harmonic, and write Q = Diag [1 / , / , / . (16)In this case ∆ A and ∆ C satisfy the relation∆ C = − A = − a G Φ ∆20 ( a ) . (17)In fact, the perturbation terms due to the (cid:96) = 2 , m = 0spherical harmonic contributes only to the deviatoric partof the inertia tensor.With the Equations (14)-(17) we are now equipped todeduce an upper limit for the ellipticity caused by star-quakes, as it will be shown in the following subsection. The pre-starquakes (rotating, stressed, elastic) and final (ro-tating, fluid) configurations will be characterized by two
MNRAS , 1–12 (2020)
Giliberti et al. slightly different inertia tensors, that can both be writtenin a diagonal form, ∆ I E = Diag[∆ A E , ∆ A E , ∆ C E ] and∆ I F = Diag[∆ A F , ∆ A F , ∆ C F ], where E stands for elas-tic and F for fluid , respectively. The explicit calculation ofthese two tensors is done by using Eqs (14), (15) and (16).If we assume, as said above, that a pure fluid star willbe more oblate with respect to an elastic one at the samerate of rotation, we can state that∆ C E ≤ ∆ C F (18)∆ A F ≤ ∆ A E . (19)In between the initial and the final fluid-like configurations,also the tensor of inertia - that is not necessarily symmetricdue to the intrinsic nature of the rupture process - can begiven in the diagonal form. Considering again only the de-viatoric (cid:96) = 2 harmonic term, we can therefore write, usingthe apex Q for the post-quakes configuration: I Q = Diag[∆ A Q , ∆ B Q , ∆ C Q ] , (20)where ∆ B is the change of the moments of inertia along oneof the principal axis. In this case ∆ A Q (cid:54) = ∆ B Q . Followingabove considerations, we expect that the post-quakes con-figuration will be “something” between the elastic and thefluid ones; therefore we require that∆ C E ≤ ∆ C Q ≤ ∆ C F ∆ A F ≤ ∆ A Q , ∆ B Q ≤ ∆ A E . (21)Using Eq (21) we see that the maximum difference ∆ B Q − ∆ A Q can be expressed as ∆ A E − ∆ A F , and thus we are ableto obtain an upper limit for the ellipticity due to a series ofstarquakes (see Eq 17): (cid:15) max = a I G (cid:104) Φ ∆ F ( a ) − Φ ∆ E ( a ) (cid:105) . (22)Therefore, to compute (cid:15) max , we have to build two differentrotating configurations for the same star: the first has anelastic crust, as sketched in Fig 1, while in the second theobject is completely fluid. For each of these two configu-rations we can extract the perturbed total potential valueat the stellar surface, calculate the corresponding change inthe inertia tensor through Eq (14) and, finally, get the valueof the maximum ellipticity using Eq (22). In order to modelthe the elastic configuration, we used a typical shear mod-ulus for cold catalyzed matter (Eq (2)) since the fact thataccreting NSs can reach very high temperature in the crust(10 K) should not affect much the shear modulus shape(Chamel & Haensel 2008). Moreover, high temperature canlower the µ value, which means that the crust is more simi-lar to a fluid than it is in the cold configuration (Hoffman &Heyl 2012). Since we are looking for for the maximum differ-ence between the elastic and fluid configuration, our choicefollows straightforwardly.In Table 1 are reported the values of (cid:15) max for a M =1 . M (cid:12) NS rotating at ν = ν o (Eq (6)) both with SLy andBSk21 EoSs. As we can see, also in this case a small changein the value of the adiabatic index value causes a large dif-ference in the ellipticity’s value. This results can be com-pared with the ones of Johnson-McDaniel & Owen (2013), As explained briefly in section 2, our configuration is fixed withthe choice of the NS’s mass, EoS and adiabatic index value. γ ∗ = 2 γ f = 2 . γ f = ∞ SLy 5 . × − . × − . × − BSk21 1 . × − . × − . × − Table 1.
Maximum ellipticity (22) calculated with SLy and BSk21EoSs for a standard NS with mass M = 1 . M (cid:12) . In all the cases ν = ν o (Eq (6)). γ ∗ = 2 γ f = 2 . γ f = ∞ SLy 39 30 25BSk21 11 8 6
Table 2.
Ratio of the maximum ellipticity of Eq (22) for a M =1 M (cid:12) NS and M = 2 M (cid:12) NS, i.e. (cid:15) max (1 M (cid:12) ) (cid:15) max (2 M (cid:12) ) , calculated with SLyand BSk21 EoSs. which are typically used in the GWs literature as bench-mark. Note that if the star is in its equilibrium configuration(i.e. γ ∗ = 2), the upper limit value that we find is lower thanthe maximum value of ellipticity ( (cid:15) ≈ − ) that a standard1 . M (cid:12) NS can sustain (Johnson-McDaniel & Owen 2013).But, as soon as we depart from this condition and γ = 2 . − . However, we underline that thesetwo estimated values of ellipticities comes from very differentperspective. In fact, the one of Johnson-McDaniel & Owen(2013) is the maximum elastic deformation that the starcan sustain before breaking, while the ellipticity given byour Eq (22) is the upper limit of the deformation that canbe reached due to the breaking process. Therefore, from ourmodel we expect that starquakes could produce somewhatlarge ellipticities in accreting (or fast rotating) NSs and inturn GWs.Our model shows also a strong dependence of the star’sellipticity on the stellar mass, see Table 2. For both theEoSs a M = 1 M (cid:12) object can produce an ellipticity aboutone order of magnitude larger than an heavier M = 2 M (cid:12) star. In section 3 we have shown that typical 1 . M (cid:12) NSs rotatingwith frequencies in the range 200 −
700 Hz may undergoa series of starquakes, and consequently, emit GWs. Nowwe will focus our attention on the consequences that thisemission might have on the NSs’ dynamical equilibrium.Our NS is tearing some material from its companionand so it is gaining angular momentum, with a rate N acc that is roughly given by (Ushomirsky et al. 2000) N acc = ˙ M √ GMa. (23)At the same time the star, that reached the crust breakingcondition and that has a non-null ellipticity, is losing angularmomentum through GWs emission at a rate N GW . This lat-ter can be written, using standard symbols, as (Ushomirskyet al. 2000) N GW = 128 π GI ν (cid:15) c . (24) MNRAS000
700 Hz may undergoa series of starquakes, and consequently, emit GWs. Nowwe will focus our attention on the consequences that thisemission might have on the NSs’ dynamical equilibrium.Our NS is tearing some material from its companionand so it is gaining angular momentum, with a rate N acc that is roughly given by (Ushomirsky et al. 2000) N acc = ˙ M √ GMa. (23)At the same time the star, that reached the crust breakingcondition and that has a non-null ellipticity, is losing angularmomentum through GWs emission at a rate N GW . This lat-ter can be written, using standard symbols, as (Ushomirskyet al. 2000) N GW = 128 π GI ν (cid:15) c . (24) MNRAS000 , 1–12 (2020) tarquakes and gravitational waves We can now use (cid:15) max to get a lower limit for the equi-librium frequency value. In fact, since the crust ruptures arecaused by fast rotation, we can express (cid:15) max as (cid:15) max = ˜ (cid:15)ν , (25)where ˜ (cid:15) is a function of the EoS and of the stellar mass. Byequating Eq (23) and Eq (24), and using the above expres-sion (25) we get the equilibrium frequency ν eq = K ˙ M / ( Ma ) / I / ˜ (cid:15) / , (26)where we have explicitly written all the terms depending onthe stellar mass and EoS, and C is the constant defined by K = (cid:113) c / π / √ G . (27)From Eq (26) we can obtain ν eq , both for SLy and BSk21EoSs, as a function of the mass and of the adiabatic indexvalue. The dynamical equilibrium clearly depends on therate ˙ M at which the star is accreting. We consider two dif-ferent thresholds that roughly constraints the region wherethe astrophysical values for this kind of objects can be found.In particular, they are the same values given by Ushomirskyet al. (2000): an upper limit of ˙ M = 2 × − M (cid:12) / yr and alower one of ˙ M = 10 − M (cid:12) / yr.The results for γ ∗ = 2 are shown in Fig 4, while thestudy of the effect of different adiabatic indeces on the equi-librium frequency is exemplified for a M = 1 . M (cid:12) NS inTables 3 and 4. In fact, we expect that a more compressiblestar has a smaller maximum ellipticity and thus a largerequilibrium frequency, if compared with an incompressibleone, and this is exactly what happens. The value of γ ∗ = 2in Fig 4 has been chosen since the curves plotted in this caseare the largest ones between the equilibrium and the frozenscenario.Let us note that, as said in section 3, the expected equi-librium frequency is smaller than ν o . Furthermore, ν eq (Ta-bles 3 and 4) is always lower than the breaking frequency(cf. Fig 2 with Fig 4), i.e. ν eq < ν b . (28)Despite the different values obtained for different adiabaticindices, the above relation remains valid. This relation sug-gest the following physical picture. An old star accretes somematerial from a companion, increasing its angular velocity.Stresses develop into the crust, till the breaking strain isreached: crust fails, the star loses its axial symmetry andstarts to emit GWs. The rate of angular momentum lostby this emission is greater than the one gained from accre-tion, and the star spins-down, till the equilibrium is reached.However, we remind that the estimated values of ν eq calcu-lated in this way are the lower value, since in Eq (26) wehave used our upper limit for the ellipticity. Furthermore,since the breaking frequency is calculated for σ max = 10 − ,when the breaking strain is smaller, we get a lower breakingfrequency.Thus, from our analysis we expect the breaking of thecrust for rapidly rotating pulsars, even in the case of veryhigh (10 − ) breaking strain; moreover, our model, throughthe development of a large ellipticity, suggests also why wedon’t observe any NS with rotational frequency above 700Hz. M / M ⊙ ν e q ( H z ) Figure 4.
Equilibrium frequency ν eq as function of the stellar massfor SLy (solid) and BSk21 (dashed) EoSs, fixed γ ∗ = 2. The curvesare calculated for two different mass accretion rates: ˙ M = 2 × − M (cid:12) / yr (orange) and ˙ M = 10 − M (cid:12) / yr (purple). The black,dashed line represents the actual maximum observed rotationalfrequency ν o . γ ∗ = 2 γ f = 2 . γ f = ∞ ν b (Hz) 585 478 265 ν Maxeq (Hz) 331 157 140 ν Mineq (Hz) 183 87 78
Table 3.
Breaking frequency ( ν b ) and equilibrium frequencies( ν eq ) calculated with SLy EoS for a M = 1 . M (cid:12) NS and dif-ferent mass accretion rates: ˙ M = 2 × − ˙ M (cid:12) / yr ( ν Maxeq ) and˙ M = 1 × − ˙ M (cid:12) / yr ( ν mineq ). γ ∗ = 2 γ f = 2 . γ f = ∞ ν b (Hz) 499 377 240 ν Maxeq (Hz) 262 128 114 ν Mineq (Hz) 145 80 63
Table 4.
Breaking frequency ( ν b ) and equilibrium frequencies( ν eq ) calculated with BSk21 EoS for a M = 1 . M (cid:12) NS and dif-ferent mass accretion rates: ˙ M = 2 × − ˙ M (cid:12) / yr ( ν Maxeq ) and˙ M = 1 × − ˙ M (cid:12) / yr ( ν mineq ). In the previous sections, we used our upper limit value (cid:15) max to calculate the NSs’ equilibrium rotational frequency. How-ever, just why it is an upper limit , we do not expect all thesources to reach the maximum possible deformation mea-sured by (cid:15) . In this section we will follow a somewhat reversepath. Starting from observational data coming from elec-tromagnetic and GWs observations we will constraint theellipticity of observed NSs.
The O O O MNRAS , 1–12 (2020)
Giliberti et al. tating NS. Searches focused both on wide-parameter sources(Abbott et al. 2005; Abbott et al. 2007; Abbott et al. 2008,2009, 2016, 2017, 2018) and signal coming from specific tar-get (Abbott et al. 2019c,a; Abbott et al. 2019d). In partic-ular, the recent paper (Abbott et al. 2019d) put constraintson the fiducial ellipticity for a selection of rapidly rotating( ν r >
100 Hz) pulsars. These estimations can be very usefulif compared with our maximum ellipticity value, Eq (22): infact, we can assume that during their life real pulsars reachonly a fraction β ≤ (cid:15) = β(cid:15) max . (29)If we state that the ellipticity of NSs is due only to thestarquakes mechanism, we can extract the value of β simplyas the ratio between LIGO/Virgo fiducial ellipticities (cid:15) L/V and our upper limit β GW = (cid:15) L/V (cid:15) max . (30)In this way, β GW > (cid:15) value; if, on the contrary, β GW < (cid:15) max currently developed on the NS.In Fig 5 we show β GW , obtained using the definition(30) and calculated for 1 . M (cid:12) , with both the EoSs. Usingthe SLy EoS we get a larger value of β GW = 0 . β GW = 0 . β GW >
1, but estimations for pul-sars close to Earth could be relevant in the next future. Forexample, the recent constraints (Abbott et al. 2019b) forJ043-4715 ( (cid:15) < . × − ) and J071-6830 ( (cid:15) < . × − )are comparable with our upper limit of about (cid:15) max (cid:39) − for both pulsars. The latest LIGO-Virgo observational run(O3) is very important: the complete analysis of its signal,integrated over many months, could show the first directdetection of continuous GWs (and thus a measurement of (cid:15) ); if this were not the case, we could use the new data tolower the estimated value of (cid:15) (i.e., the estimation of β ). Inthe debate whether to search for GWs emitted by slowly ro-tating or fast rotating objects our model suggests to searchGWs emission from high rotating pulsars that we expect tohave larger ellipticities (and therefore greater gravitationalemission power) than the slowly rotating ones. Observation of LMXBs can also be used to extract a value of (cid:15) , giving an useful benchmark to compare it with our upperlimit (cid:15) max . In fact, assuming that the measured rotationalfrequency of an observed LMXB is its equilibrium frequency ,one can obtain the corresponding ellipticity (cid:15) acc = C / ν / r (cid:68) ˙ M (cid:69) / ( Ma ) / I , (31) ●● ● ● ● ●● ●● ●● ●●● ●● ● ●● ● ●●● ●●● ●●●●●● ● ●●● ● ●●● ●● ●●●● ●●● ● ●● ●● ●●● ●● ● ●● ●●●● ●● ●●● ● ● ●● ●● ● ●●● ●● ● ●● ● ●● ●●●● ● ●●●● ● ●● ●●●● ●●●● ●● ● ●●● ● ●●● ●●● ●● ●●● ●●●● ▲▲ ▲ ▲ ▲ ▲▲ ▲▲ ▲▲ ▲▲▲ ▲▲ ▲ ▲▲ ▲ ▲▲▲ ▲▲▲ ▲▲▲▲▲▲ ▲ ▲▲▲ ▲ ▲▲▲ ▲▲ ▲▲▲▲ ▲▲▲▲ ▲▲ ▲▲ ▲▲▲ ▲▲ ▲ ▲▲ ▲▲▲▲▲▲ ▲▲▲ ▲ ▲ ▲▲ ▲▲ ▲ ▲▲▲ ▲▲ ▲ ▲▲ ▲ ▲▲ ▲▲▲▲ ▲ ▲▲▲▲ ▲ ▲▲ ▲▲▲▲ ▲▲▲▲ ▲▲ ▲ ▲▲▲ ▲ ▲▲▲ ▲▲▲ ▲▲▲▲▲ ▲▲▲▲ ● SLy ▲ BSk21
100 200 300 400 500 600 7000.1110100 ν rot ( Hz ) β G W Figure 5.
Ratio of the observational ellipticity given by non-detections of continuous GWs and our upper limit (cid:15) max . β GW (30) is calculated for M = 1 . M (cid:12) with two different EoS: SLy(green circles) and BSk21 (red triangles). The dashed, black curveindicates β GW = 1. where acc stands for accretion and (cid:68) ˙ M (cid:69) is the average massaccretion rate during outburst . In the following, we use thedata elaborated by Haskell et al. (2015) (see Table 1 therein),that give, for each star, its rotational frequency, distanceand mass accretion rate. It is convenient to introduce theparameter β acc = (cid:15) acc (cid:15) max . (32)The meaning of β acc is straightforward: every star with β acc ≤ β acc > M = 1 . M (cid:12) while the adiabatic in-dex value at γ ∗ = 2: the values of (cid:15) max for γ f = 2 . , ∞ arealways larger than the one obtained with γ ∗ , giving a lower β acc .The value of β acc as a function of observed rotationalfrequency is shown shown in Fig 6. About 95% of the starsfall in the category β acc ≤ β acc found with this method is 0 . .
09 for SLy EoS. In our sample only J1756.9-2508, in the case of Sly EoS, has a value of β acc larger than1; this is due to the fact that this star has a rotating fre-quency of “only” 182 Hz, and thus has an (cid:15) max value smallerthan that of other stars. These results confirms that star-quakes mechanism could explain why these stars have all afrequency smaller than 700 Hz. In fact, even a small fractionof our maximum value (cid:15) max is enough for LMXBs to reacha dynamical equilibrium at frequency smaller than ν o . β We can also obtain an alternative estimation of β basedon the actual observational upper value of frequency. If weassume that for a given NS (cid:15) is only a fraction β o of (cid:15) max , Typically accreting NSs show short bursts, lasting from daysto months, with a corresponding high accretion rate, and verylong period of recovery, during which the accretion is orders ofmagnitude smaller than in active phase (Watts et al. 2008).MNRAS000
09 for SLy EoS. In our sample only J1756.9-2508, in the case of Sly EoS, has a value of β acc larger than1; this is due to the fact that this star has a rotating fre-quency of “only” 182 Hz, and thus has an (cid:15) max value smallerthan that of other stars. These results confirms that star-quakes mechanism could explain why these stars have all afrequency smaller than 700 Hz. In fact, even a small fractionof our maximum value (cid:15) max is enough for LMXBs to reacha dynamical equilibrium at frequency smaller than ν o . β We can also obtain an alternative estimation of β basedon the actual observational upper value of frequency. If weassume that for a given NS (cid:15) is only a fraction β o of (cid:15) max , Typically accreting NSs show short bursts, lasting from daysto months, with a corresponding high accretion rate, and verylong period of recovery, during which the accretion is orders ofmagnitude smaller than in active phase (Watts et al. 2008).MNRAS000 , 1–12 (2020) tarquakes and gravitational waves ● ●● ●● ●● ● ● ● ● ● ●● ●●●● ● ▲ ▲▲ ▲▲ ▲▲ ▲ ▲ ▲ ▲ ▲ ▲▲ ▲▲▲▲ ▲ ● SLy ▲ BSk21
200 300 400 500 6000.0050.0100.0500.1000.5001 ν r ( Hz ) β a cc Figure 6.
Ratio of the observational ellipticity assuming dynamicalequilibrium in LMXB objects (cid:15) acc (Eq (32)) and our upper limit (cid:15) max (Eq (22)) for a M = 1 . M (cid:12) NS with γ ∗ = 2. Green dotsare calculated with SLy EoS while red triangles with BSk21. Thedashed, black curve indicates β acc = 1. i.e. (cid:15) = β o (cid:15) max ; β ≤ , (33)one can express the equilibrium frequency as a function of β o , as in Eq (26), namely ν eq ( β o ) = C ˙ M / ( Ma ) / I / β / o ˜ (cid:15) / , . (34)Then, given an EoS and fixed the stellar mass, we state thatthe minimum feasible value of β o is the one that satisfies thecondition ν eq ( β o ) = ν o . (35)Using this selection criterion for different masses M and dif-ferent mass accretion rates ˙ M we construct the curves shownin Fig. 7. In all the cases β o <
1, and for the smaller massaccretion rate ˙ M = 10 − M (cid:12) / yr (with M = 1 . M (cid:12) ), wefind, for the SLy EoS, β o = 0 . β o = 0 . In the previous two sub-sections we focused only on accret-ing NSs, while in the following the analysis is extended alsoto spinning-down objects. Millisecond pulsars are thought tobe the evolutionary descendant of accreting objects (Alparet al. 1982; Bhattacharya & van den Heuvel 1991), i.e. oldNSs that have been spun up to high rotational frequenciesvia accretion. If we assume that LMXB objects can developlarge ellipticity due to starquakes, we can expect that mil-lisecond stars too have a non-zero (cid:15) , i.e. a residual part oftheir initially larger quadrupolar deformation. In this case For these objects it could be useful to take into accounts alsothe long-time evolution of elastic layers due to a non-null viscos-ity. However, at the present time, the viscosity ν of NSs’ crust isessentially unknown, even if in the last years some first estima-tions have appeared (Kwang-Hua 2018; Lander & Gourgouliatos2019). The inclusion also of this parameter in a consistent model M / M ⊙ β o Figure 7.
Minimum value of β o , defined implicitly in Eq (35), forSLy (solid) and BSk21 (dashed) EoSs. The curves are calculatedfor a NS with a mass of M = 1 . M (cid:12) , γ ∗ = 2 for two differentmass accretion rates: ˙ M = 2 × − M (cid:12) / yr (orange) and ˙ M =10 − M (cid:12) / yr (purple). The black, dotted line represents β = 1. we can introduce a simple model to explain the millisecondactual decreasing period. In fact, these objects lose energyvia electromagnetic and GWs emission. Therefore, follow-ing Woan et al. (2018), that assume I = 10 g cm and avacuum dipole radiation, we can write˙ P − = 0 . (cid:18) P (cid:19) (cid:18) I g cm (cid:19) − (cid:18) B s G (cid:19) + (36)+2 . (cid:18) P (cid:19) (cid:18) I g cm (cid:19) (cid:16) (cid:15) − (cid:17) , where B s is the surface magnetic field of the star. In orderto to calculate the ellipticity needed to explain the observedstellar spin-down in the case of pure gravitational wave emis-sion, i.e. in the gravitar limit (Palomba 2005), we can neglectthe magnetic field and put B s = 0 in Eq (35). Solving for (cid:15) (the subscript gr stands for gravitar) we get (cid:15) gr = 10 − (cid:115) . (cid:18) I g cm (cid:19) (cid:18) ˙ P − (cid:19) (cid:18) P (cid:19) . (37)We can compare this value with our upper limit, Eq (22),assuming a standard star configuration (with M = 1 . M (cid:12) and an adiabatic index γ ∗ = 2). As representatives of mil-lisecond pulsars we select all the non-accreting stars with ν r >
100 Hz in the ATNF database ( ). In Fig. 8 we showthe ratio β gr = (cid:15) gr (cid:15) max . (38)The spin-down of all the stars with β gr ≤ (cid:15) max correspondingto the NS’s rotational frequency.On the contrary, we observe that for all the objects with β gr > can help the understanding of the global dynamics of a realisticNS, and is clearly a very interesting field for future research.MNRAS , 1–12 (2020) Giliberti et al. ▲ ▲▲ ▲▲ ▲ ▲ ▲▲ ▲ ▲▲▲ ▲▲ ▲▲▲ ▲ ▲▲ ▲ ▲ ▲▲ ▲▲ ▲▲ ▲ ▲ ▲▲▲ ▲ ▲▲ ▲ ▲▲ ▲▲ ▲▲ ▲▲▲▲ ▲ ▲▲ ▲▲ ▲▲ ▲ ▲▲▲▲ ▲ ▲▲ ▲▲ ▲▲ ▲ ▲ ▲▲ ▲ ▲▲▲ ▲▲▲ ▲▲▲ ▲ ▲▲▲ ▲▲▲ ▲▲▲ ▲▲ ▲▲▲▲ ▲▲ ▲▲ ▲▲ ▲ ▲▲ ▲▲ ▲ ▲ ▲▲▲▲▲ ▲ ▲▲▲ ▲▲ ▲▲ ▲▲▲▲ ▲ ▲▲ ▲ ▲ ▲▲ ▲ ▲ ▲▲ ▲ ▲▲ ▲▲▲ ▲▲ ▲▲ ▲▲ ▲▲ ▲▲ ▲▲ ▲▲ ▲ ▲▲ ▲▲ ▲▲ ▲▲▲▲ ▲ ▲▲▲ ▲▲ ▲▲ ▲ ▲▲ ▲▲▲ ▲ ▲ ▲ ▲▲▲ ▲ ▲▲▲▲ ▲ ▲▲ ▲▲▲ ▲▲▲▲▲ ▲▲ ▲ ▲ ▲▲▲▲▲ ▲▲ ▲▲ ▲ ▲ ● ●● ●● ● ● ●● ● ●●● ●● ●●● ● ●● ● ● ●● ●● ●● ● ● ●●● ● ●● ● ●● ●● ●● ●●●● ● ●● ●● ●● ● ●●●● ● ●● ●● ●● ● ● ●● ● ●●● ●●● ●●● ● ●●● ●●● ●●● ●● ●●●● ●● ●● ●● ● ●● ●● ● ● ●●●●● ● ●●● ●● ●● ●●●● ● ●● ● ● ●● ● ● ●● ● ●● ●●● ●● ●● ●● ●● ●● ●● ●● ● ●● ●● ●● ●●●● ● ●●● ●● ●● ● ●● ●●● ● ● ● ●●● ● ●●●● ● ●● ●●● ●●●●● ●● ● ● ●●●●● ●● ●● ● ● ▲ BSk21 ● SLy ν r ( Hz ) β g r Figure 8.
Ratio of the observational ellipticity in the gravitar limit (cid:15) gr (Eq (37)) and our upper limit (cid:15) max for a M = 1 . M (cid:12) NSwith γ ∗ = 2. Green dots are calculated with SLy EoS while redtriangles with BSk21. The dashed, black line indicates β gr = 1. gravitar hypothesis, using the full expression of Eq (35). Theminimum value of β gr in our sample is β gr = 0 . . (39)In the case of SLy roughly 90% of the selected star have β gr <
1, while for the stiff BSk21 EoS the percentage rise to97%. These results show that the starquakes mechanism canproduce very large ellipticity, and thus that actual value of (cid:15) for millisecond pulsars can in principle be produced by crustrupture on their progenitors accreting stars. Fig 8 shows alsoanother interesting aspect. The slowest NSs in our cataloguehave a β gr value that is larger than 1, which means that forthese objects our upper limit is smaller than the gravitarellipticity value.Finally we note that also the recent estimation ofJ1023+0038 pulsar’s ellipticity (Bhattacharyya 2020) con-firm the range of expected β form LMXBs. In fact, for a1 . ± . M (cid:12) NS with γ ∗ = 2, the ratio between Bhat-tacharyya’s ellipticity and our maximum values is about10 − both for SLy and BSk21 EoSs. GWs emission from rapidly rotating NSs is a very actual andinteresting field of research, and so the study of the maxi-mum mountains on pulsars. Many scenarios as been invokedto produce a non-null ellipticity, such as thermal mountainsor particular magnetic fields configurations; other studies fo-cused on the maximum quadrupole moment that a NS cansustain. In this paper it is presented for the first time (atleast at the best of our knowledge) a realistic and consistentcalculation of a new mechanism to produce GWs emission:starquakes. This kind of mechanism has been already re-cently proposed by Fattoyev et al. (2018), but here it is thefirst time that the full problem is consistently studied. Infact, we construct a model to study the reaching of crustbreaking, the maximum ellipticity produced in a sequenceof crust ruptures and the equilibrium frequency reached byaccreting object. Our calculations shows that NSs crust canfail due to centrifugal stresses when the frequency of the ro-tating star is in the range 200 −
900 Hz, depending on the EoS and the mass of the star. In general, the equilibriumfrequency is found to be smaller than the breaking one and,therefore, it is also below the actual observable threshold of716 .
36 Hz. The study of how large the ellipticity due to star-quakes can be has given an upper limit for (cid:15) lying between10 − and 10 − , depending on the EoS and the mass of theobject. The comparison between different EoS showed thatthe stiffer ones produce larger ellipticities and, consequently,bring the star towards a lower equilibrium frequency.In the case of highly spinning pulsars, and for γ differ-ent from γ ∗ , we also showed that crust failure can producean ellipticity comparable with the maximum theoreticallyexpected value given by Johnson-McDaniel & Owen (2013).But with a great difference; in fact, the mechanism produc-ing the quadrupolar deformation is due to crust breaking inour model, instead of be sustained by the elastic crust.We found that the stellar mass affects the star’s re-sponse and that a 2 M (cid:12) objects creates an ellipticity aboutone order of magnitude smaller than the lighter 1 M (cid:12) ones.Moreover, (cid:15) depends strongly also on the adiabatic index:even just a small difference from the adiabatic equilibriumvalue leads to large differences on the ellipticity.Last, but not least, comparing our upper limit estima-tion with the observational data coming from both accretingobjects and millisecond pulsars, and calculating within ourmodel the ellipticity deduced from the data, we found a valueof (cid:15) that is even a very small fraction of our model upperlimit ( (cid:15)/(cid:15) max = 10 − ) can in principle explain the vast ma-jority of the available data (90% to 97% depending on theEoS used).Our model explains why rotating NSs with frequenciesgreater than about 700 Hz are not observed and, at the sametime, it expects that accreting NSs stars could generate, withthe right combination of distance, mass and frequency, GWsof intensity detectable by LIGO-Virgo O3 run ,which com-plete data analysis is still ongoing.It also shows that the evolutionary scenario depictedin this paper seems to be sufficiently robust to be a com-petitive candidate in describing NSs achieve a dynamicalequilibrium.All these results are obtained in a Newtonian frame-work. A natural and important step for further improve-ment of this model is its generalisation to General Relativity,which permits also to use realistic EoSs for the descriptionof the whole star. ACKNOWLEDGMENTS
The authors are grateful for useful discussions with L. Per-otti, B. Haskell, M. Antonelli and P. M. Pizzochero.
DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding author.
REFERENCES
Abbott B. e. a., 2016, Phys. Rev. Lett., 116, 061102Abbott B. e. a., 2017, Phys. Rev. Lett., 119, 161101MNRAS000
Abbott B. e. a., 2016, Phys. Rev. Lett., 116, 061102Abbott B. e. a., 2017, Phys. Rev. Lett., 119, 161101MNRAS000 , 1–12 (2020) tarquakes and gravitational waves Abbott B., et al., 2005, Phys. Rev. D, 72, 102004Abbott B., et al., 2007, Phys. Rev. D, 76, 082001Abbott B., et al., 2008, Phys. Rev. D, 77, 022001Abbott B., et al., 2009, Phys. Rev. D, 79, 022001Abbott B. P., et al., 2016, Phys. Rev. D, 94, 042002Abbott B. P., et al., 2017, Phys. Rev. D, 96, 062002Abbott B. P., et al., 2018, Phys. Rev. D, 97, 102003Abbott B. P., et al., 2019a, Phys. Rev. D, 99, 122002Abbott B. P., et al., 2019b, Phys. Rev. D, 100, 024004Abbott B. P., et al., 2019c, Phys. Rev. D, 100, 122002Abbott B. P., et al., 2019d, ApJ, 879, 10Alpar M. A., Cheng A. F., Ruderman M. A., Shaham J., 1982,Nature, 300, 728Baiko D. A., Chugunov A. I., 2018, MNRAS, 480, 5511Baym G., Pethick C., Pines D., Ruderman M., 1969, Nature, 224,872Bhattacharya D., van den Heuvel E. P. J., 1991, Phys. Rep., 203,1Bhattacharyya S., 2020, MNRAS, 498Bildsten L., 1998, ApJ, 501, L89Chakrabarty D., 2008, AIP Conference ProceedingsChakrabarty D., Morgan E., Muno M., Galloway D., WijnandsR., Van Der Klis M., Markwardt C., 2003, Nature, 424, 42Chamel N., Haensel P., 2008, Living Reviews in Relativity, 11, 10Chao B. F., Gross R., 1987, Geophysic R. Astro. Soc., 91, 569Christensen R., 2013, The Theory of Materials Failure. OxfordUniversity PressCook G. B., Shapiro S. L., Teukolsky S. A., 1994, ApJ, 423, 117Cutler C., 2002, Phys. Rev. D, 66, 084025Cutler C., Ushomirsky G., Link B., 2003, ApJ, 588, 975Douchin F., Haensel P., 2001, A&A, 380, 151Fattoyev F. J., Horowitz C. J., Lu H., 2018, preprint,( arXiv:1804.04952 )Franco L. M., Link B., Epstein R. I., 2000, ApJ, 543, 987Galloway D. K., Markwardt C. B., Morgan E. H., ChakrabartyD., Strohmayer T. E., 2005, ApJ, 622, L45Giliberti E., Antonelli M., Cambiotti G., Pizzochero P. M., 2019,Publ. Astron. Soc. Australia, 36, e036Giliberti E., Cambiotti G., Antonelli M., Pizzochero P., 2020,MNRAS, 491, 1064Goriely S., Chamel N., J.M. P., 2010, Phys. Review C, 82, 035804Haskell B., Jones D. I., Andersson N., 2006, MNRAS, 373, 1423Haskell B., Samuelsson L., Glampedakis K., Andersson N., 2008,MNRAS, 385, 531Haskell B., Priymak M., Patruno A., Oppenoorth M., Melatos A.,Lasky P. D., 2015, MNRAS, 450, 2393Hessels J. W. T., Ransom S. M., Stairs I. H., Freire P. C. C.,Kaspi V. M., Camilo F., 2006, Science, 311, 1901Hoffman K., Heyl J., 2012, MNRAS, 426, 2404Horowitz C. J., 2010, Phys. Rev. D, 81, 103001Horowitz C. J., Kadau K., 2009, Phys. Rev. Lett., 102, 191102Johnson-McDaniel N. K., Owen B. J., 2013, Phys. Rev. D, 88,044004Kwang-Hua C. W., 2018, Ap&SS, 363, 184Lander S. K., Gourgouliatos K. N., 2019, MNRAS, 486, 4130Lattimer J. M., Prakash M., 2007, Phys. Rep., 442, 109Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, AJ,129, 1993Melatos A., Payne D. J. B., 2005, ApJ, 623, 1044Palomba C., 2005, MNRAS, 359, 1150Papaloizou J., Pringle J. E., 1978, MNRAS, 184, 501Priymak M., Melatos A., Payne D. J. B., 2011, MNRAS, 417,2696Ruderman M., 1969, Nature, 223, 597Ruderman M., 1991, ApJ, 382, 576Sabadini R., Vermeersen B., Cambiotti G., 2016, Global Dy-namics of the Earth: Applications of Viscoelastic Relax-ation Theory to Solid-Earth and Planetary Geophysics. Springer Netherlands, https://books.google.pl/books?id=33xBDAAAQBAJ
Ushomirsky G., Cutler C., Bildsten L., 2000, MNRAS, 319, 902Vigelius M., Melatos A., 2009, MNRAS, 395, 1972Wagoner R. V., 1984, Annals of Physics, 278, 345Watts A. L., Krishnan B., Bildsten L., Schutz B. F., 2008, MN-RAS, 389, 839Woan G., Pitkin M. D., Haskell B., Jones D. I., Lasky P. D., 2018,preprint, ( arXiv:1806.02822 ) APPENDIX A: INERTIA AND SPHERICALHARMONICS EXPANSION
The inertia tensor is defined as I ij = (cid:90) ρ ( r ) (cid:0) r δ ij − r i r j (cid:1) dV. (A1)The initial, non-rotating stellar configuration will be de-formed by centrifugal force, that gives rise to a change ofthe star’s density profile, can be written as ρ ( r ) = ρ ( r ) + ρ ∆ ( r, θ, φ ) , (A2)where we have put in evidence the initial, unstressed profile ρ and the local perturbation ρ ∆ . We can use the sphericalsymmetry of the problem to recast Eq (A1). To that purposelet us first introduce the spherical harmonics, defined as Y (cid:96)m ( θ, ϕ ) = P (cid:96)m cos ( θ ) e imϕ , (A3)where P (cid:96)m are the associated Legendre polynomials P (cid:96)m ( x ) = (cid:96) (cid:96) ! (cid:0) − x (cid:1) m/ d (cid:96) + m ( x − ) (cid:96) dx (cid:96) + m m > − m ( (cid:96) − m )!( (cid:96) + m )! P (cid:96)m ( x ) m < ρ ∆ ( r, θ, ϕ ) = ∞ (cid:88) (cid:96) =0 m = (cid:96) (cid:88) m = − (cid:96) ρ ∆ (cid:96)m ( r ) Y (cid:96)m ( θ, ϕ ) (A5)and the total perturbed potential asΦ ∆ ( r, θ, ϕ ) = ∞ (cid:88) (cid:96) =0 (cid:96) (cid:88) m = − (cid:96) Φ ∆ (cid:96)m ( r ) Y (cid:96)m ( θ, ϕ ) . (A6)In this work, all the perturbed terms are the ones due torotation, thus we can focus only on the centrifugal potential φ C . It can be expanded as a sum of the only (cid:96) = 0 and (cid:96) = 2terms, i.e. φ C ( r, θ, ϕ ) = φ C ( r ) Y ( θ, ϕ ) + m = (cid:96) (cid:88) m = − (cid:96) φ C m ( r ) Y m ( θ, ϕ ) . (A7)If we substitute the density expansion in Eq (A1), we canexpress the inertia tensor I as I = I + I ∆ , (A8)where we highlighted the unperturbed tensor of inertia I .With some straightforward algebra the perturbed inertiatensor can be divided into two terms I ∆ = I ∆00 + I ∆20 . (A9)Note that by choosing a coordinate system in which the rota-tional axis z coincides with the one at θ = 0, the centrifugal MNRAS , 1–12 (2020) Giliberti et al. ▲ ▲ ▲ ▲ ▲ ▲ ▲ - - - - M / M ⊙ ϵ m ax Figure B1.
Maximum elliptiticy given by a sequence of starquakes,calculated as the difference between the fluid and the deformedelastic configuration as in Eq (22), showed as blue dots, and asthe difference between the fluid and the spherical elastic con-figuration, Eq (B1), reported as red triangles. The ellipticity iscalculated for different stellar masses ranging from 1 M (cid:12) to 2 M (cid:12) ,keeping ν = 1 Hz and γ ∗ = 2 fixed. potential contains only the m = 0 order of the (cid:96) = 2 har-monic term. The contributions of these spherical harmonicterms perturbed to the inertia tensor are, respectively, I ∆00 = 8 π δ ij (cid:90) a ρ ∆00 ( r ) r dr, (A10)and I ∆20 = (cid:18) δ ij − ˆ r i ˆ r j (cid:19) π (cid:90) a ρ ∆20 ( r ) r dr. (A11)For simplicity, in the main text we use the notation∆ I = I ∆20 . (A12) APPENDIX B: ALTERNATIVE ESTIMATION OFELLIPTICITY
Another heuristic estimation of ellipticity can be obtainedby assuming that the elastic crust can keep the NS in aspherical configuration, despite the centrifugal forces due torotation (Fattoyev et al. 2018). In this case Φ ∆ E ( a ) = 0 andone gets (see Eq 22) (cid:15) max = a I G Φ ∆ F ( a ) . (B1)We observe this rough approximation to be larger than our (cid:15) max , since the rotating configuration is squeezed towardsthe equatorial plane by the fast rotation. In Fig B1 we com-pare, for different NSs masses and ν = 1 Hz, the ellipticitygiven by our upper limit (cid:15) max and the one given by Eq (B1)for the case γ ∗ = 2. The first of the two ellipticity valueshas clearly a different dependence on the stellar mass and,furthermore, it is 5 orders of magnitude smaller than thesecond. This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000