Spin evolution of a proto-neutron star Giovanni Camelio,1,∗ Leonardo Gualtieri,1,† Jos´e A. Pons,2,‡ and Valeria Ferrari1,§ 1Dipartimento di Fisica, “Sapienza” Universit`a di Roma & Sezione INFN Roma1, P.A. Moro 5, 00185, Roma, Italy. 2Departament de F´ısica Aplicada, Universitat dAlacant, Ap. Correus 99, 03080 Alacant, Spain. We study the evolution of the rotation rate of a proto-neutron star, born in a core-collapse supernova, in the first seconds of its life. During this phase, the star evolution can be described as a sequence of stationary configurations, which we determine by solving the neutrino transport and thestellarstructureequationsingeneralrelativity. Weincludeinourmodeltheangularmomentum loss due to neutrino emission. We find that the requirement of a rotation rate not exceeding the mass-shedding limit at the beginning of the evolution implies a strict bound on the rotation rate at later times. Moreover, assuming that the proto-neutron star is born with a finite ellipticity, we determinetheemittedgravitationalwavesignal,andestimateitsdetectabilitybypresentandfuture 6 ground-based interferometric detectors. 1 0 PACSnumbers: 04.40.Dg,97.60.Jd,26.60.Kp 2 l u I. INTRODUCTION rate of the progenitor is sufficiently high. On the other J hand,astrophysicalobservationsofyoungpulsarpopula- 2 tions (see [14] and references therein) show typical peri- 1 When a supernova explodes, it leaves a hot, lepton- ods (cid:38) 100 ms. richand(presumably)rapidlyrotatingremnant: aproto- ] neutron star (PNS). In the early stages of its evolution, Thequasi-stationaryevolutionofaPNSdrivenbyneu- E the PNS cools down and loses its high lepton content, trinotransportinasphericallysymmetricspace-timehas H whileitsradiusandrotationratedecrease. Inthisphase, been extensively studied in the past, quite often adopt- . ahugeamountofenergyandofangularmomentumisre- ing an equation of state (EoS) obtained within a finite- h p leased, mainlyduetoneutrinoemission[1–3]. Afraction temperature, field-theoretical model solved at the mean - of this energy is expected to be emitted in the gravi- field level [3, 15–17]. This approach yields a sequence o tational wave channel; indeed, as a consequence of the ofthermodynamicalprofilesdescribingthestructureand r t violent collapse non-radial oscillations can be excited, the early evolution of a non-rotating PNS. A different s making PNSs promising sources for present and future approach has been used in [6], where an EoS obtained a [ gravitational detectors [4–7]. within a finite-temperature many-body theory approach In the first tenths of seconds after its birth, the wasemployed,buttheneutrinotransportequationswere 2 PNS is turbulent and characterized by large instabil- not explicitly solved (a set of entropy and lepton frac- v 5 ities. During the next tens of seconds, it undergoes tion profiles were adopted, having the same qualitative 4 a more quiet, “quasi-stationary” evolution (the Kelvin- behaviourasthoseof[3]). Wealsomentionthatthenon- 9 Helmholtz phase), which can be described as a sequence radial oscillations of the quasi-stationary configurations 2 of equilibrium configurations [1, 3]. In this article, we obtainedwiththesedifferentapproacheshavebeenstud- 0 study the evolution of the rotation rate of a PNS during iedin[4,6,18],wherethequasi-normalmodefrequencies 1. this quasi-stationary, Kelvin-Helmholtz phase. An accu- of the gravitational waves emitted in the early PNS life 0 rate modeling of this phase is needed, for instance, to were computed. 6 compute the frequencies of the PNS gravitational wave The evolution of rotating PNSs has been studied 1 emission. Moreover, it provides a link bewteen super- in [19], where the thermodynamical profiles obtained : v nova explosions, a phenomenon which is still not fully in [3] for a non-rotating PNS were employed as effective i understood, and the properties of the observed popula- one-parameterEoSs;therotatingconfigurationswereob- X tion of young pulsars. Current models of the evolution tainedusingthenon-linearBGSMcode[20]tosolveEin- ar of progenitor stars [8], combined with numerical simula- stein’s equations. A similar approach has been followed tions of core collapse and explosion (see e.g. [9–13]), do in [21], which used the profiles of [3] and [6]. The main not allow sufficiently accurate estimates of the expected limitations of these works are the following. rotation rate of newly born PNSs; they only show that the minimum rotation period at the onset of the Kelvin- • The evolution of the PNS rotation rate is due not Helmoltz phase can be as small as few ms, if the spin onlytothechangeinthemomentofinertia(i.e.,to the contraction), but also to the angular momen- tumchangeduetoneutrinoemission[22]. Thiswas neglected in [19], and described with an heuristic formula in [21]. ∗ [email protected] † [email protected] ‡ [email protected] • As we shall discuss in this paper, when the PNS § [email protected] profiles describing a non-rotating star are treated 2 as an effective EoS, one can obtain configurations c = G = 1). The perfect fluid of the star is described which are unstable to radial perturbations. by the stress-energy tensor T = ((cid:15)+p)u u +pg , µν µ ν µν where uµ = (e−φ/2,0,0,0) is the fluid four-velocity and In this article, we study the quasi-stationary evolution (cid:15),p are the energy density and pressure of the fluid, re- of a spherically symmetric PNS, solving the relativistic spectively. The gravitational mass inside a radius r is equations of neutrino transport and of stellar structure. m(r) = r(1−e−λ)/2. At the surface of the star, r = R, The details of our code will be discussed in [23], where the pressure vanishes and the metric matches with the it will be applied to more recent EoSs. Here, we employ exterior Schwarzschild metric, with M = m(R) as the the same EoS used in [3] (i.e. GM3 [24]), to study the gravitationalmassofthestar. Wealsodefinethebaryon spin evolution of the PNS in its first tens of seconds of number inside a radius r, life. To model an evolving, rotating PNS, we use the profiles of entropy per baryon and lepton fraction s(a), (cid:90) r Y (a) (a is the number of baryons enclosed in a sphere a(r)=4π eλ/2n r(cid:48)2dr(cid:48), (2) L b passing through the point considered) obtained with our 0 evolution code which describes a non-rotating PNS. Our where n is the baryon number density. The position in- approach is different from that used in [19], as will be b sidethestarcanbedescribedeitherbythecoordinater, discussed in detail in Sec IIIB. In order to determine or by the enclosed baryon number a. We also define the the PNS spin evolution, we model the evolution of angu- rest-mass density ρ = m n (m is the neutron mass), larmomentum(duetoneutrinoemission)usingEpstein’s n b n and the baryon mass of the star M =m a(R). We will formula[22]. Wealsodiscussthegravitationalwaveemis- b n use M = 1.60 M , which corresponds, in the calcula- sion which could be associated with this process. b (cid:12) tionsofthispaper,toagravitationalmassof1.55M at Theplanofthepaperisasfollows. InSec.IIwebriefly (cid:12) 200msfromthecorebounce,whichreducesto∼1.4M describe our approach to model the PNS evolution in its (cid:12) in the first ten seconds of life of the PNS. quasi-stationaryphase. InSec.IIIwedescribeourmodel of a rotating PNS. In Sec. IV we show our results, and Since the PNS has a temperature of several MeV, its in Sec. V we draw our conclusions. The details of the EoS is non-barotropic, and can be written as slowly rotating model are described in Appendix A. (cid:15)=(cid:15)(p,s,{Y } ), (3) i i II. EARLY EVOLUTION OF A where s is the entropy per baryon, and Y = n /n i i b PROTO-NEUTRON STAR is the fraction of the i-th specie, with number density n . Assuming that the matter is in beta equilibrium, i The quasi-stationary, Kelvin-Helmholtz phase of a the dependence on the composition {Y } can be cast i i PNSstartsfewhundredsofmsafterthecorebounce[1,3]. into a dependence on only the electron-type lepton frac- This phase consists of two evolutionary stages. In the tion Y = n /n . Different choices of thermodynamical L L b first few tens of seconds, neutrinos diffuse from the low- variables are possible, for instance replacing the entropy entropy core to the high-entropy envelope, deleptonizing per baryon s with the temperature T. In this paper, the core and increasing its temperature. In the second we use the finite-temperature EoS GM3 of Glendenning phase,thestarisleptonpoorbuthot,theentropygradi- and Moszkowski [24], obtained within a field-theoretical ent is smoothed out, and thermally produced neutrinos model solved at the mean field level, where the interac- cooldownthePNS.Afteraboutoneminute,thestarbe- tions between baryons are mediated by the exchange of comestransparenttoneutrinosandcanbeconsideredas mesons; it contains only nucleonic degrees of freedom. a “mature” neutron star, with a radius of ∼10−15 km This is the same EoS employed in [3]; we consider the and a temperature <1 MeV. case of matter composed by electrons, protons and neu- ThePNSevolutionintheKelvin-Helmholtzphasecan trons. More recent EoSs, based on a many-body theory be considered as a sequence of quasi-stationary configu- approach, will be considered in a future work [23]. rations, because the hydrodynamical timescale is much InordertosolvetheTOVequations,weneedtoknow, smaller than the evolution timescale. Following [3], we at each point, the energy density as a function of the model this phase by solving the general relativistic neu- pressure; thus, we need to know the EoS and the ther- trino transport equations coupled with the structure modynamical profiles s(a), Y (a), which are obtained by L equations, assuming spherical symmetry. In each quasi- solving the transport equations stationary configuration, the spacetime metric has the form ∂Y ∂(4πeφ/2r2F ) L + ν =0, (4) ds2 =−eφ(r)dt2+eλ(r)dr2+r2(dθ2+sin2θdϕ2), (1) ∂t ∂a ∂s ∂Y ∂(4πeφr2H ) T +µ L +e−φ/2 ν =0, (5) where φ(r) and λ(r) are radial functions, obtained by ∂t ν ∂t ∂a solving the Tolman-Oppenheimer-Volkov (TOV) equa- tions (in this paper we use geometrized units, in which where F and H are, respectively, the neutrino number ν ν 3 and energy fluxes 6 0.2 s F = − e−λ+2φT2(cid:34)D ∂(Teφ/2) +(Teφ/2)D ∂η(cid:35), (6) 5 15..00 ss ν 6π2(cid:126)3 3 ∂r 2∂r 10. s H = − e−λ+2φT3(cid:34)D ∂(Teφ/2) +(Teφ/2)D ∂η(cid:35), (7) 4 ν 6π2(cid:126)3 4 ∂r 3∂r s 3 2 whereη =µ /T istheelectron-typeneutrinodegeneracy ν 1 parameterandD aretheneutrinodiffusioncoefficients, n which are computed assuming the diffusion approxima- tion [3]. In order to preserve causality and stabilize the 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 code in the semi-transparent regions near the PNS sur- m [M ] b O• face, we apply a flux-limiter [25]. Our code evolves the PNS by iteratively solving, at FIG. 1. Entropy per baryon as a function of the enclosed each time-step, (i) the transport equations (4) and (5) baryon mass m =m a at t=0.2,1,5,10 s (k =1). using an implicit scheme and (ii) the TOV equations by b n B relaxation method. The time evolution keeps the baryon mass M constant, and provides, at each time-step, a b quasi-stationaryconfigurationofthe(non-rotating)star, described by the profiles of all the thermodynamical quantities (p, (cid:15), nb, s, YL, etc.) as functions of a (or of 0.35 r). Westartourintegrationat200msfromcorebounce. 0.2 s 1.0 s The initial profiles (which are the same employed in [3]) 0.3 5.0 s are the result of core-collapse simulations [26]. In Fig- 10. s ures 1 and 2 we show the evolutionary profiles of the 0.25 entropy per baryon s and the electron-type lepton frac- 0.2 tion Y as functions of the enclosed baryon mass. We L L have checked that the total energy and lepton number Y 0.15 are conserved during the evolution within a few percent in the early stages of the evolution, and with more ac- 0.1 curacy in later stages. We remark that this error can be significantly reduced by reducing the timestep; however, 0.05 this accuracy is sufficient for the aims of this work. Our code will be described in detail in a future work [23]. 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 m [M ] b O• Recently,differentapproacheshavebeenappliedinthe study of the PNS evolution (see e.g. [17]), in which the neutrino spectrum is described with greater accuracy by means of multi-group codes. However, since in this work FIG. 2. Electron-type lepton fraction as a function of the enclosed baryon mass at t=0.2,1,5,10 s. we are not interested in the details of the neutrino emis- sion, we prefer to employ a simpler and faster energy- averaged approach (as in [3]). As mentioned above, our code also employs a flux-limiter [25], which makes III. A MODEL OF ROTATING it difficult to establish the precise location of the neu- PROTO-NEUTRON STARS trinosphere. Both the neutrinosphere and the neutrino spectra are better determined with more complex core- A. Slowly rotating stars in general relativity collapse codes, which however are far slower, while our PNS code is suitable to run for longer evolution times. Typical core-collapse codes run for at most 500 ms af- We model a rotating PNS using the perturbative ap- ter core bounce, whereas we can easily explore the first proach of Hartle and Thorne [27, 28] (see also [29]). The minute ofPNSlife, at theendof whichthe starbecomes rotating star is described as a stationary perturbation of neutrino-transparent. a spherically symmetric background, for small values of 4 the angular velocity Ω = 2πν, i.e., ν (cid:28) ν (ν is the Einstein’s equations, expanded in powers of Ω and in ms ms mass-shedding frequency, at which the star starts losing Legendre polynomials, can be written as a set of ordi- massattheequator,seeSec.IIID).Asshownin[21],this narydifferentialequationsfortheperturbationfunctions; “slowrotation”approximationisreasonablyaccuratefor theseequationsaresummarizedinAppendixA.Foreach rotationratesupto∼0.8ofthemass-sheddinglimit,pro- value of the central pressure p (or, equivalently, of the c viding values of mass, equatorial radius and moment of central energy density (cid:15) ) and of the rotation rate Ω, c inertia which differ by (cid:46)0.5% from those obtained with the numerical integration of the perturbation equations fully relativistic, nonlinear simulations. In our approach yieldstheperturbedfunctions,andthenthevaluesofthe we assume uniform rotation; PNSs are expected to have multipolemomentsofthestar(inparticular,themassM a significant amount of differential rotation at birth [30] and the angular momentum J), and of its baryonic mass which, however, is likely to be removed by viscous mech- M . ThesequantitiescanbewrittenasM =M(0)+δM, b anisms,suchas,forinstance,magnetorotationalinstabil- J =δJ,M =M(0)+δM ,etc.,wherethequantitieswith b b b ity [31], in a fraction of a second. superscipt (0) refer to the non-rotating star with central This work should be considered as a first step towards pressurep , andthequantitieswithδ arethecorrections c a more detailed description of rotating PNSs, in which due to rotation. we shall include differential rotation. Givenanon-rotatingstarwithcentralpressurep and c The spacetime metric, up to third order in Ω, can be baryonmassM ,therotatingstar(withspinΩ)withthe b written as same central pressure has a baryon mass M(0) +δM , b b ds2 =−eφ(r)[1+2h (r)+2h (r)P (µ)]dt2 which is generally larger than Mb. Therefore, a rotating 0 2 2 star with same baryon mass M as the non-rotating one, (cid:20) (cid:21) b 2m (r)+2m (r)P (µ) +eλ(r) 1+ 0 2 2 dr2 has necessarily a smaller value of the central pressure, r−2M(r) p +δp ,withδp <0(thisisnotsurprising: whenastar c c c +r2[1+2k (r)P (µ)] (8) is set into rotation, its central pressure decreases). 2 2 We mention that in [32] the neutrino transport equa- ×[dθ2+sin2θ{dφ−[ω(r)+w (r)+w (r)P(cid:48)(µ)]dt}2] 1 3 3 tions for a rotating star in general relativity have been solved by using an alternative approach. In this ap- whereµ=cosθ andP (µ)istheLegendrepolynomialof n proach (which is believed to be accurate for slowly ro- order n, the prime denoting the derivative with respect tating stars [32]) the structure and transport equations to µ. The perturbations of the non-rotating star are de- for a spherically symmetric star are modified by adding scribed by the functions ω (of O(Ω)), h , m and h , 0 0 2 a centrifugal force term, to include the effect of rotation. m , k (of O(Ω2)), and w , w (of O(Ω3)). The energy- 2 2 1 3 momentum tensor is Tµν =(E +P)uµuν +Pgµν (9) B. Including the thermodynamical profiles where g , uµ are the metric and four-velocity in the In order to integrate the structure equations of a cold µν rotating configuration, and we denote by calligraphic neutronstarweneedtoassignanequationofstatewhich, letters thermodynamical quantities (energy, density and inthecaseofPNSs,isnon-barotropic,i.e. (cid:15)=(cid:15)(p,s,Y ), L pressure) in the rotating star. An element of fluid, at thus we also need to know the profiles of entropy and position (r,θ) in the non-rotating star, is displaced by leptonfractionthroughoutthestar. AsdiscussedinSec- rotation to the position tion II, these profiles are obtained by our evolutionary code for spherical, non-rotating PNS at selected values r¯=r+ξ(r,θ), (10) of time. The non-rotating profiles can be used to compute the where ξ(r,θ) = ξ (r) + ξ (r)P (µ) + O(Ω4) is the La- 0 2 2 structure of a rotating PNS in different ways. A possible grangian displacement. approach is the following. In the Hartle-Thorne approach, one assumes that if Let us consider a spherical PNS with baryon mass M the fluid element of the non-rotating star has pressure P b at a given value of the evolution time t. The numeri- and energy density (cid:15), the displaced fluid element of the cal code discussed in Sec. II provides the functions p(a), rotating star has the same values of pressure and energy (cid:15)(a), s(a), Y (a), whereweremindthataistheenclosed density. In other words, the Lagrangian perturbations L baryon number. If we replace the inverse function of of the thermodynamical quantities (cid:15), P vanish (see [27], p(a) into the non-barotropic EoS, we obtain an “effec- eq.(6)); the modification of these quantities is only due tivebarotropicEoS”,(cid:15)˜(p)=(cid:15)(p,s(a(p)),Y (a(p)),which to the displacement (10): L canbeusedtosolvetheTOVequationsforthespherical d(cid:15) dP configuration to which we add the perturbations due to δ(cid:15)(r,θ)=− ξ(r,θ) , δP(r,θ)=− ξ(r,θ). (11) rotation, according to Hartle’s procedure. Since the ro- dr dr tatingstarmusthavethesamebaryonmassasthespher- We remark that as long as we neglect terms of O(Ω4), ical star, one can proceed as follows: (i) solve the TOV δ(cid:15)(r,θ)(cid:39)δ(cid:15)(r¯,θ). equations for a spherically symmetric star with central 5 pressure p +δp ; (ii) solve the perturbation equations ing PNS. Therefore, we are neglecting the effect of rota- c c for a chosen value of the rotation rate, to determine the tion on the time evolution of the PNS. To be consistent, actual baryon mass of the rotating star with same cen- we should have integrated the transport equations ap- tral pressure; (iii) iterate these two steps modifying δp propriate for a rotating star, which are much more com- c until the baryon mass coincides with the assigned value plicated. Sincetheseapproximationsaffectthetimescale M . This approach was used in [19], where the rotating of the stellar evolution, we would like to estimate how b star was modeled solving the fully non-linear Einstein faster, or slower, the rotating star looses its thermal and equations. lepton content with respect to the non-rotating one. However, this procedure has somerelevantdrawbacks. Sincetheevolutiontimescaleisgovernedbyneutrinodif- Indeed, during the first second after bounce the star is fusion processes, at each time step of the non-rotating very weakly bound, and it may happen that the proce- PNS evolution, we have computed and compared the dureaboveyieldsδpc >0,whichindicatesthatthesecon- neutrino diffusion coefficients Dn (see Eqns. (6), (7)) for figurations are in the unstable branch of the mass-radius non-rotatingandrotatingconfigurations. Thelatterhave diagram. We think that this is caused by the unphysical been obtained by replacing the profiles (p(a), (cid:15)(a), etc.) treatment of the thermodynamical profile (effectively, as ofanon-rotatingPNSwiththoseofarotatingPNS(com- a barotropic EoS). puted as discussed above in this Section). In the upper This problem did not occur in the simulations of [19] and middle panels of Fig. 3 we plot D2,D3 and D4 as becausetheauthorsconsideredadifferent, stablebranch functions of the enclosed baryon mass mB = mna, for ofthemass-radiuscurvecorrespondingtothe“effective” the non-rotating (solid line) and rotating (dashed line) EoS (cid:15)˜(p), at much lower densities. Indeed, for t (cid:46) 0.5 s, configurations,att=0.2s,t=1.2sandt=10s. Inthe at the center of the star they had n ∼ 10−2fm−3 (i.e., lower panels we plot the neutrino number density and b rest-mass density ρ∼1013g/cm3), which corresponds to the total energy density at the same times. We assume the outer region of the star modeled in [3]. When the Mb = 1.6M(cid:12) and that the initial angular momentum, central density is so low, only a small region of the star Jin, is equal to the maximum angular momentum Jmax, isdescribedbytheGM3EoS;therestisdescribedbythe abovewhichmass-sheddingsetsin(seeSec.IVAforfur- low-density EoS used to model the PNS envelope, which ther details). We see that the diffusion coefficients of does not yield unstable configurations. the rotating configurations are larger than those of the non-rotating star. For m (cid:46) 1 M the relative differ- Since wewantto modelthe PNSconsistently withthe B (cid:12) ence |Drot−Dnon rot|/|Dnon rot| is always smaller than evolutionary models given in [3], we decided to imple- n n n ∼ 10 − 20%, and becomes smaller than a few percent ment the non-rotating profiles in an alternative way. As after the first few seconds. in the previous approach, we consider the spherical con- figuration obtained by the evolution code at time t, with central density p and baryon mass (constant during the In the outer region m (cid:38) 1 M and early times, the c B (cid:12) evolution) M . To describe the rotating star, we use the relative difference seems larger, in particular for the co- b GM3 EoS (cid:15) = (cid:15)(p,s,Y ); since we are restricting our efficient D , but this has no effect for two reasons: first, L 3 analysis to slowly rotating stars, the entropy and lep- as shown in the two lower panels of Fig. 3, both the ton fraction profiles s(a) and Y (a) of the non-rotating neutrino number density and the total energy density L star are a good approximation for those of the rotating are much smaller than in the inner core; therefore, even star. We follow the steps discussed at the end of Sec- though the diffusion coefficients of the rotating star are tion IIIA: (i) solve the TOV equations for a star with larger than those of the non-rotating one, few neutrinos central pressure p +δp ; at each value of a, the energy are trapped in this region and transport effects do not c c density is (cid:15)(p,s(a),Y (a)); (ii) solve Hartle’s perturba- contribute significantly to the overall evolution; second, L tionequations,findingthebaryonmassofthestarrotat- the differences become large in the semi-transparent re- ingtoagivenratewiththisreducedcentralpressureand gion, when the mean free path becomes comparable to find the correction to the baryon mass due to rotation; (or larger than) the distance to the star surface. In this (iii) iterate the first two steps, finding δp such that the region the diffusion approximation breaks down and in c baryon mass of the rotating star is M . We remark that practice the diffusion coefficients are always numerically b the energy density of the rotating star in step (ii) is re- limited (a flux-limiter approach). lated to that of the non-rotating star in step (i) by the Hartle-Thorne prescription described above Eq. (11). From the above discussion we can conclude that the Sinceweareusinganappropriatenon-barotropicEoS, rotating star looses energy and lepton number through the instability discussed above disappears, and the cen- neutrino emission faster than the non-rotating one. This tral pressure of the rotating star is, as expected, smaller effect is larger at the beginning of the evolution, i.e. for than that of the non-rotating star with same baryon t (cid:46) 2 s, and is of the order of ∼ 10−20%, but becomes mass. negligible at later times. Consequently our rotating star We stress again that we are using the numerical so- coolsdownandcontractsoveratimescalewhich,initially, lution of the transport equations (5) for a non-rotating is∼10−20%shorterthanthatofthecorrespondingnon- PNS, to build quasi-stationary configurations of a rotat- rotating configuration. 6 100000 0.2 s m] 1000 1.2 s 10 s [2 D 10 0.1 100000 m] [3 1000 D 10 m] 100000 [4 1000 D 10 0.1 0.01 3] -m 0.0001 [fν n 1e-06 3] m V/f 100 e M [ε 10 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 m [M ] B O• FIG. 3. Neutrino diffusion coefficients D (n=2,3,4), as functions of the enclosed baryon mass, computed using the density n and thermodynamical profiles of the non-rotating (solid line) and rotating (dashed line) configurations, at t=0.2 s, t=1.2 s and t = 10 s (upper and middle panels). Profiles of neutrino number density and energy density (lower panels). We assume M = 1.6M and that the angular momentum is the maximum allowed J = J (for J > J , the PNS reaches the b (cid:12) in max in max mass-shedding limit during its evolution, see Sec. IVA). C. Evolution of the angular momentum and of the carries away ∼ 10% of the star gravitational mass [35], rotation rate and also a significant fraction of the total angular mo- mentum [36]. To our knowledge, the most sensible es- timate of the neutrino angular momentum loss in PNSs Once the equations describing the rotating configura- has been done by Epstein in [22] tion are solved for each value of the evolution time t and for an assigned value of the rotation rate Ω, the solution dJ 2 of these equations allows one to compute the multipole =− qL R2Ω (12) dt 5 ν moments of the rotating star, including the angular mo- mentum J. Conversely we can choose, at each value of where R is the radius of the star, L = −dM/dt is the ν t, the value of the angular momentum, and determine, neutrino energy flux, and q is an efficiency parameter, using a shooting method, the corresponding value of the which depends on the features of the neutrino transport rotation rate. andemission. Ifneutrinosescapewithoutscattering,q = Ifwewanttodescribetheearlyevolutionofarotating 1;if,instead,theyhaveaveryshortmeanfreepath,they PNS,weneedaphysicalprescriptionforthetimedepen- are diffused up to the surface, and then are emitted with denceofJ. Forinstance,wemayassumethattheangular q = 5/3. As discussed in [22] (see also [37–39]), q = 5/3 momentumisconstant,asin[19](seealso[33,34]). How- should be considered as an upper limit of the angular ever, in the first minute of a PNS life, neutrino emission momentum loss by neutrino emission. A more recent, 7 alternative study [40] indicates an angular momentum In the case of a newly born PNS the situation is dif- emission smaller than this limit. In the following, we ferent. As the star contracts, due to the processes re- shall consider Epstein’s formula with q = 5/3, and this lated to neutrino production and diffusion, its rotation hastobemeantasanupperlimit. Wealsomentionthat rate increases. If the PNS has a finite ellipticity, it emits asimplifiedexpressionbasedonEpstein’sformulaforthe gravitationalwaves, whoseamplitudeandfrequencyalso angularmomentumlossinPNSshasbeenderivedin[36] increase as the star spins up. The timescale of this pro- and used in [21]. cess is of the order of tens of seconds. In our model, for simplicity we shall assume that the PNS ellipticity remains constant over this short time interval. D. Mass-shedding frequency Unfortunately, the ellipticity of a PNS is unknown. In cold, old NSs it is expected to be, at most, as large as ∼10−5−10−4 [48, 49] (larger values are allowed for EoS As mentioned in Sec. IIIA, the perturbative approach including exotic matter phases [50, 51]). For newly born which we use to model a rotating star is accurate up to ν (cid:46) 0.8ν , where ν is the mass-shedding fre- PNSs, it may be larger, but we have no hint on its ac- ms ms tual value. To our knowledge, current numerical simu- quency. The only quantity which is poorly estimated lations of core-collapse do not provide estimates of the is, of course, the mass-shedding frequency itself. There- PNS ellipticity. We remark that although there is ob- fore,ν willbedeterminedusinganumericalfitderived ms servational evidence of large asymmetries in supernova in [41] from fully relativistic, non-linear integrations of explosions [52, 53], there is no evidence that they can be Einstein’s equations: inherited by the PNS. In the following, we shall assume (cid:115) (cid:15) = 10−4, but this should be considered as a fiducial M/M ν (Hz)=a (cid:12) +b (13) value: the gravitational wave amplitude (which is linear ms R/1km in(cid:15))canbeeasilyrescaledfordifferentvaluesofthePNS ellipticity. where a = 45862 Hz and b = −189 Hz. We remark that the coefficients a,b of this fit do not depend on the EoS. IV. RESULTS E. Gravitational wave emission A. Spin evolution of the proto-neutron star IftheevolvingPNSisbornwithsomedegreeofasym- In Figure 4 we show how the angular momentum metry, it emits gravitational waves. Assuming that the changes according to Epstein’s formula (12) as the PNS starrotatesaboutaprincipalaxisofthemomentofiner- evolves. We assume q = 5/3 and baryonic mass M = tia tensor, i.e., that there is no precession1, gravitational b 1.6M . We consider different values of the angular mo- waves are emitted at twice the orbital frequency ν, with (cid:12) mentum J at the beginning of the quasi-stationary amplitude [44–47] in phase(t=0.2safterthebounce): J =2.02×1048ergs, in J = 3.71×1048ergs and J = 8.08×1048ergs. We 4G(2πν)2I (cid:15) in in h (cid:39) 3 . (14) find that, in the first ten seconds after bounce, 13% of 0 c4r the initial angular momentum is carried away by neutri- nos if J = 2.02×1048ergs or J = 3.71×1048ergs; The deviation from axisymmetry is described by the el- in in 20% of the initial angular momentum is carried away if lipticity (cid:15), defined as J = 8.075×1048ergs. As mentioned above, q = 5/3 in I −I should be considered as an upper bound; for smaller val- (cid:15)= 1 2 (15) ues of q, the rate of angular momentum loss would be I 3 smaller. where I , I and I are the principal moments of inertia The corresponding evolution of the PNS rotation fre- 1 2 3 of the PNS and I is assumed to be aligned with the quency is shown in Figure 5. In the same Figure we 3 rotation axis. For old neutron stars, the loss of energy also show the mass-shedding frequency ν , computed ms throughgravitationalwavesiscompensatedbyadecrease using the fit (13). We see that if J =8.08×1048ergs, in of rotational energy, which contributes to the spin-down the curves of ν(t) and of ν (t) cross during the quasi- ms ofthestar(themaincontributiontothespin-downbeing stationary evolution; before the crossing, the PNS spin that of the magnetic field). is larger than the mass-shedding limit. This means that a PNS with such initial angular momentum would lose mass. If we require the initial rotation rate to be smaller than the mass-shedding limit, we must impose 1 Free precession requires the existence of a rigid crust [42], thus Jin ≤ Jmax ≡ 3.72 × 1048ergs. We remark that the it should not occur in the first tens of seconds of the PNS life, value of J is not affected by the efficiency of angu- max whenthecrusthasnotformedyet[43]. lar momentum loss q: if q < 5/3, J has the same max 8 Howeverthisspindowntimescaleismuchlongerthanthe timescaleofthequasi-stationaryevolutionweareconsid- Mb=1.6 MO• ering; therefore it is unlikely that after this early phase 2 the PNS rotation rate is larger than ∼300 Hz (i.e., that its period is smaller than ∼ 3 ms), unless some spin- 1.8 up mechanism (such as e.g. accretion) sets in. A less 1.6 efficient angular momentum loss (q < 5/3) would mod- 2J (km) 11..24 Jin=JmJJaiinnx===200...095020 kkkmmm222 ewroaIutteldilsyrweinomcratrheinansteohtteihnsigasmtfihena.atlmvoaldueel,sbofutprteh-esugpeenrenroavlapsitcetlularer 1 evolution [8] predict a similar range of the PNS rota- tion rate and angular momentum. Among the models 0.8 considered in [8], the only one with J > J (and ro- max 0.6 tation period smaller than 3 ms) is expected to collapse to a black hole. Other works [9, 10] have shown that if 0.4 2 4 6 8 10 the progenitor has a rotation rate sufficiently large, the t (s) PNS resulting from the core-collapse can have periods of few ms; our results suggest that this scenario is un- likely, unless there is a significant mass loss in the early FIG.4. Angularmomentumevolutionduetoneutrinolosses, Kelvin-Helmoltz phase. foraPNSwithbaryonicmassM =1.6M andinitialangu- b (cid:12) lar momentum J =(2.02, 3.71, 8.08)×1048erg s. in B. Gravitational wave emission If the PNS has a finite ellipticity (cid:15) (which we assume, for simplicity, to remain constant during the first ∼ 10 Mb=1.6 MO• 1000 s of the PNS life), it emits gravitational waves with fre- quency f(t)=2ν(t) and amplitude given by Eq. (14), 900 νms Jin=2.00 km2 800 Jin=Jmax 700 Jin=0.50 km2 h (cid:39) 4G(2πν(t))2I3(t)(cid:15). (16) 600 0 c4r z) (Hν 500 As the spin rate ν(t) increases, both the frequency and 400 the amplitude of the gravitational wave increase; there- 300 fore, the signal is a sort of “chirp”; this is different from the chirp emitted by neutron star binaries before co- 200 alescence, because the amplitude increases at a much 100 milder rate. In Figure 6 we show the strain amplitude 0 √ √ (cid:113) 2 4 6 8 10 h˜(f) f = f (h˜ (f)2+h˜ (f)2)/2, where h˜ (f) + × +,× t (s) are the Fourier transform of the two polarization of the gravitational wave signal FIG.5. EvolutionofthePNSrotationrate,correspondingto 1+cos2i h =h cos(2πf(t)t) (17) the angular momentum profiles shown in Figure 4. + 0 2 h =h cosisin(2πf(t)t), (18) × 0 andiistheanglebetweentherotationaxisandthelineof value, but the rotation rate grows more rapidly than in sight. In Figure 6 the signal strain amplitude, computed Figure 5. assuming optimal orientation, J =J , (cid:15)=10−4 and in max It is interesting to note that, since νms has a steeper a distance of r = 10 kpc, is compared with the sensi- increase than ν(t), even when the bound ν ≤νms is sat- tivitycurvesofAdvancedVirgo2, AdvancedLIGO3,and uratedatthebeginningofthequasi-stationaryphasethe of the third generation detector ET4. We see that the frequencybecomesmuchsmallerthanthemassshedding frequencyatlatertimes. Thisisana-posterioriconfirma- tion that the slow rotation approximation is appropriate to study newly born PNSs. For t>10 s, the PNS radius 2 https://inspirehep.net/record/889763/plots does not change significantly, and the star starts to spin 3 https://dcc.ligo.org/LIGO-T0900288/public down due to electromagnetic and gravitational emission. 4 http://www.et-gw.eu/etsensitivities 9 signal is marginally above noise for the advanced detec- up phase. For a PNS of M = 1.6M we find that one b (cid:12) tors, but it is definitely above the noise curve for ET. minuteafterbouncethestarwouldrotateatν (cid:46)300Hz, This signal would be seen by Advanced Virgo with a corresponding to a rotation period τ (cid:38)3.3×10−3s. min signal-to-noiseratioSNR=1.4,andbyAdvancedLIGO IfthePNSisbornwithafiniteellipticity(cid:15),whilespin- with SNR = 2.2, too low to extract it from the detec- ning up it emits gravitational waves at twice the rota- tor noise; however, since the signal-to-noise ratio scales tion frequency. This signal increases both in frequency linearly with the ellipticity, a star born with (cid:15) = 10−3 and amplitude. We find that for a galactic supernova, would be detected with SNR = 14 and SNR = 22 by if (cid:15) = 10−3 this signal could be detected by Advanced Advanced Virgo and LIGO, respectively. The third gen- LIGO/Virgo with a signal-to-noise ratio (cid:38)14. To detect eration detectors like ET would detect the signal coming farthersources,thirdgenerationdetectorslikeETwould from a galactic PNS born with (cid:15) = 10−4 with a very be needed. large signal-to-noise ratio, i.e. SNR = 22. If the source We remark that the actual value of PNS ellipticities isintheVirgocluster(d=15Mpc),theellipticityofthe is unknown, and depends on the details of the supenova PNSshouldbeaslargeas5·10−2 tobeseenbyETwith core collapse. Accurate numerical simulations of super- SNR=8. nova explosion addressing this issue are certainly needed to provide a quantitative estimate of the range of (cid:15). We also remark that in our approach the effects of the PNS rotation are consistently included in the structure ε=10-4 r=10 kpc equations, but they are neglected when solving the neu- 10-20 Advanced Virgo trino transport equations. We estimate that due to this ET approximation, we overestimate the evolution timescale Advanced LIGO 10-21 Jin=Jmax at early times of, at most, ∼ 10−20%. Moreover, since we are not interested in the details of the neutrino dy- 2 10-22 namics and we need a fast code to evolve the star for 1/ν tensofseconds,weperformenergyaveragestodetermine ) ν the neutrino diffusion coefficients, and we apply a flux- h(10-23 limiter; these approximations should not significantly af- fect the thermodynamical evolution of the PNS and its 10-24 gravitational wave emission. Thisworkisafirststepinthestudyoftheearlyevolu- tion of PNSs. A paper with a detailed description of our 10-25 10 100 1000 numericalcode,anditsextensiontomorerecentEoSs,is ν (Hz) in preparation [23]. Further developments shall include differentialrotation,convectionandgeneralizationofthe neutrino transport equations to rotating PNSs. √ FIG. 6. The strain amplitude h˜(f) f of the gravitational wave signal emitted by a PNS with (cid:15)=10−4, J =J , at in max a distance r = 10 kpc, is compared with the noise curves of ACKNOWLEDGMENTS Advanced Virgo, Advanced LIGO and ET. We thank O. Benhar and A. Lovato for useful dis- cussions on the EoS of the PNS. We also thank S. ReddyandL.F.Robertsforusefuldiscussionsontheneu- trino cross sections. This work was partially supported V. CONCLUDING REMARKS by “NewCompStar” (COST Action MP1304), and by the H2020-MSCA-RISE-2015 Grant No. StronGrHEP- In this paper we have studied the angular momentum 690904. J.A.P. acknowledges support by the MINECO loss, the time dependence of the rotation rate and the grants AYA2013-42184-P and AYA2015-66899-C2-2-P. gravitationalwaveemissionofanewlybornPNS,during thefirsttensofsecondsafterbounce. Theearlyevolution of the rotating PNS has been modeled using the entropy Appendix A: Hartle-Thorne equations and lepton fraction profiles consistently computed solv- ingthegeneralrelativistictransportequationsforanon- Here,webrieflydescribetheequationsoftheperturba- rotating star; angular momentum loss due to neutrino tive Hartle-Thorne approach discussed in Sec. IIIA. For emission has been modeled using Epstein’s formula [22]. further details we refer the reader to [27, 28, 54] and to During this early evolution, the star spins up due to the Appendix of [29]. contraction. By requiring that the initial rotation rate The spacetime metric (up to order O(Ω3)) is given by does not exceed the mass-shedding limit, we have esti- Eq. (8); it depends on the background functions φ(r), matedthemaximumrotationrateattheendofthespin- λ(r) = −log(1 − 2m(r)/r), and on the perturbations 10 functions h (r), m (r) (l = 0,2), k (r), ω(r), w (r) (l = and l l 2 l 1,3). The energy and pressure (Eulerian) perturbations are dm d(cid:15) u2 8πr5((cid:15)+P)χ2 0 =4πr2 [δp ((cid:15)+P)]+ + δP =((cid:15)(r)+P(r))(δp0(r)+δp2(r)P2(µ)) dr dP 0 12r4 3(r−2M) d(cid:15)/dr dδp u2 m (1+8πr2P) δ(cid:15)= δP , (A1) 0 = − 0 dP/dr dr 12r4(r−2M) (r−2M)2 anddependontheperturbationfunctionspl(r)(l=0,2). −4π((cid:15)+P)r2δp0 + 2r2χ (cid:104)u The background spacetime is described by the TOV r−2M 3(r−2M) r3 equations: (r−3M −4πr3P)χ(cid:21) + . (A6) dm r−2M =4πr2(cid:15) Matchingtheinteriorandtheexteriorsolutionsatr =R, dr itispossibletocomputethecorrectiontothemassdueto dφ m+4πr3P =2 stellarrotation,δM =m (R)+J2/R3,andthemonopo- 0 dr r(r−2m) lar stellar deformation. The baryonic mass correction dP (cid:15)+P dφ δM =δm (R) is given by solving the equation =− . (A2) b b dr 2 dr The mass of the non-rotating configuration is obtained (cid:20)(cid:18) (cid:19) dδm m 1 by matching at the stellar surface r = R the interior b =4πr2eλ/2 1+ 0 + r2(cid:36)2e−φ (cid:15) dr r−2m 3 solution with the exterior (Schwarzschild) solution, i.e., (cid:21) M =m(R). Moreover,thebaryonicmassM ofthenon- d(cid:15)/dr b + ((cid:15)+P)δp . (A7) rotating configuration is obtained integrating the equa- dP/dr 0 tion dm /dr =4πr2eλ/2ρ, and computing M =m (r). b b b The spacetime perturbation to first order in Ω is de- scribedbythefunctionω(r), whichisresponsibleforthe The l=2 perturbations satisfy the equations dragging of inertial frames; it satisfies the equations dχ u 4πr2((cid:15)+P)χ dv dφ 1 1dφ (cid:20)8πr5((cid:15)+P)χ2 u2 (cid:21) = − (A3) 2 =− h +( + ) + dr r4 r−2M dr dr 2 r 2dr 3(r−2M) 6r4 du 16πr5((cid:15)+P)χ dh (cid:20) dφ r dφ 4M (cid:21) = , (A4) 2 = − + ( )−1(8π((cid:15)+P)− ) h dr r−2M dr dr r−2M dr r3 2 where(cid:36) =Ω−ω,j(r)=e−φ/2(1−2M/r)1/2,χ=j(cid:36)and − 4v2 (dφ)−1+ u2 (cid:20)1dφr− 1 (dφ)−1(cid:21) u=r4jd(cid:36)/dr. TheangularmomentumJ isobtainedby r(r−2M) dr 6r5 2dr r−2M dr matching the interior with the exterior solution χ(r) = 8πr5((cid:15)+P)χ2 (cid:20)1dφ 1 dφ (cid:21) Ω−2J/r3, u(r)=6J at r =R. The moment of inertia, + r+ ( )−1 , (A8) 3(r−2M) 2dr r−2M dr at zero-th order in the rotation rate, is I =J/Ω. The perturbations to second order in Ω are described by the metric functions hl(r), ml(r) (l = 0,2), k2(r), where v2 = k2 + h2. Matching the interior and exte- and by the fluid pressure perturbations δpl. The l = 0 rior solutions, it is possible to determine the quadrupole perturbations satisfy the equations moment of the PNS and its quadrupolar deformation. d (cid:18) χ2r3 (cid:19) The equations for the peturbations at O(Ω3), wl(r) δp +h − =0 dr 0 0 3(r−2M) (l = 1,3), have a similar structure but they are longer and are not reported here; we refer the reader to [29, χ2r3 δp +h − =0 54]. They yield the octupole moment, the third-order 2 2 3(r−2M) corrections to the angular momentum and the second- (A5) order corrections to the moment of inertia. [1] A. Burrows and J. M. Lattimer, Astrophys. J. 307, 178 arXiv:astro-ph/9807040 [astro-ph]. (1986). [4] V.Ferrari,G.Miniutti, andJ.A.Pons,Mon.Not.Roy. [2] W. Keil and H. T. Janka, Astron. Astrophys. 296, 145 Astron. Soc. 342, 629 (2003), arXiv:astro-ph/0210581 (1995). [astro-ph]. [3] J. A. Pons, S. Reddy, M. Prakash, J. M. Lattimer, [5] C. D. Ott, Class. Quant. Grav. 26, 063001 (2009), and J. A. Miralles, Astrophys. J. 513, 780 (1999), arXiv:0809.0695 [astro-ph].