Rotating Hyperonic Neutron Stars

ROTATING HYPERONIC NEUTRON STARS Anna Campoy Ordaz Master’s Thesis December call 2018/2019 Supervisors: Àngels Ramos Gómez Laura Tolós Rigueiro Abstract Neutron stars are among the most compact objects known in the Universe, with central densities several times that encountered in the center of atomic nuclei. The properties of neutron stars, such as their masses and radii, depend strongly on their internal composition, and hence the equation of state, which is still an open question. Most neutron stars are observed as pulsars, which rotate with periods as low as the millisecond. It is therefore of fundamental importance to study how rotation affects its structure and this is the main purpose of this project. To this aim, the TOV equations with the Hartle- Thorne approach for rotating neutron stars are solved for different angular velocities of rotation employing two different hadronic models: A pure nucleonic equation of state for nuclear matter, widely used in the literature, and a recently developed hyperonic equation of state, the FSU2H model, which allows two-solar-mass neutron stars with hyperons in their cores. It is found that, due to its stiffer behaviour, the FSU2H equation of state leads to smaller Keplerian frequencies than the pure nucleonic model (H&H) and, hence, the FSU2H rotating neutron stars suffer bigger deformations and reach higher moments of inertia. Particularly, the less massive neutron stars are the systems that experience the larger effects due to rotation. Contents 1 Introduction 1 2 Theoretical framework 3 2.1 Historicaloverview................................ 3 2.2 Structure of neutron stars . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 3 Equations of state of β-stable neutron star matter 8 3.1 Beta-stability conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2 General Equations-of-State considerations . . . . . . . . . . . . . . . . . . . 11 3.3 Nucleonic matter: H&H model . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.4 Nucleonic and Hyperonic matter: FSU2H model . . . . . . . . . . . . . . . 14 3.5 Equationsofstate ................................ 18 4 Non-rotating neutron stars: TOV equations 20 5 Rotating neutron stars: Hartle-Thorne approach 22 5.1 Hartle-Thorne approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5.2 Numericalcalculation .............................. 28 6 Results 30 6.1 Non-rotating neutron stars . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 6.2 Rotatingneutronstars.............................. 32 6.2.1 Deformation ............................... 33 6.2.2 Mass.................................... 38 6.2.3 Momentofinertia ............................ 42 7 Conclusions 45 i 8 Bibliography 47 A Appendix 50 A.1 Programme .................................... 50 ii Chapter 1 Introduction A neutron stars is a type of stellar compact remnant that results from the gravitational core-collapse of an ordinary giant star with a mass in the range between 8-25Mduring a Type-II, Ib or Ic supernova event [1]. The supernova explosion occurs when the red giant star has exhausted all its possibilities for energy production by nuclear fusion. Therefore, the pressure gradient provided by radiation is not sufficient to balance the gravitational attraction, consequently, the star becomes unstable and it eventually collapses. The inner dense regions of the star collapse first and gravitational energy is released and transferred to the outer layers, blowing them away as an expanding and fast-moving shock wave into the surrounding interstellar medium. This turns into an expanding shell of gas and dust which is observed as the supernova remnant. Once formed, neutron stars no longer actively generate heat and cool over time. They are in hydrostatic equilibrium as the gravitational collapse is compensated by neutron degeneracy pressure and repulsive nuclear forces. Neutron degeneracy pressure is based on the Pauli exclusion principlei; however, this mechanism can just hold up an object until 0.7M. For this reason, repulsive nuclear forces play a crucial role in supporting the more massive neutron stars observed. Nevertheless, if the remnant star exceeds a maximum mass Mmax ∼3M(Kalogera [2]), the neutron degeneracy pressure and the nuclear forces will not be enough to support the neutron star and it will collapse into a black hole. In this context, neutron stars are the most compact objects known without an event horizon and have masses in the range between 1-2M. The upper limit of 2Mhas been determined in very recent measurements of neutron star binary systems (Demorest et al. [3], Antoniadis et al. [4] and Cromartie et al. [5] works). Meanwhile, their estimated range for radii goes from 10.7 km to 13.1 km (Lattimer & Prakash [6]). iQuantum mechanical principle which states that two or more identical fermions cannot occupy the same quantum state simultaneously. 1 The properties of neutron stars, such as the mass and radius, strongly depend on their dense interior although its composition is still unknown and, hence, the equation of state (EoS) remains also undetermined. The knowledge of the EoS is a fact of great interest for many fields of Physics, ranging from nuclear physics, particle physics to astrophysics. As the true behaviour of the nuclear EoS at very high densities remains uncertain, there exist many hypothesis about the particles and phases of matter that could be possibly present in the core of neutron stars. Such models in the literature include hyperons, pion and kaon condensates, as well as deconfined quark matter. Furthermore, the recent detection of the first direct gravitational wave event (GW170817) due to the merging of two neutron star in a binary system by LIGO and Virgo collaboration on 17 August 2017 has opened a new window of possibilities into astrophysics, dense matter, gravitation, and cosmology. This event might unlock a new path through new observables for a more detailed knowledge of the EoS at high densities for neutron stars and, as for instance their tidal deformability [7]. For example, a maximum star limit at 2.17 Mhas been suggested by Margalit & Metzger [8] taking into account the measurement of the tidal deformation. It should be emphasized that neutron stars are mostly observed as pulsars, namely neutron stars in rotation that emit periodic electromagnetic radiation. They are typically classified in terms of their period of rotation into millisecond- or second-type pulsars [9]. In consequence, the main objective of this project is to study how rotation affects neutron stars by analysing their properties after solving the hydrostatic equilibrium equations within the Hartle-Thorne approach for rotating systems. Particularly, the behaviour of the dense neutron star interior is described by an hyperonic EoS recently released by the Barcelona group [10], the so-called FSU2H model. Results are compared with those obtained with a pure nucleonic EoS widely used in the literature [11]. This purpose is achieved by first solving the non-rotating neutron stars equations and later implementing the ones for rotating neutron stars via a numerical procedure, which implies the construction and development of a FORTRAN programme, a task that has represented a personal challenge. Nowadays, pulsars rotating with periods bigger than 700 Hz have not been observed yet, although a physics explanation preventing a more rapid rotation does not exist. For this reason, and motivated by the recently proposed sub-millisecond pulsars after the GW170817 event, we will explore different angular velocities (Ω=1000, 3000, 5000 Hz) into the neutron stars obtained by the different EoS. 2 Chapter 2 Theoretical framework In this chapter, a brief summary of the historical discoveries related to neutron stars is presented. Furthermore, the structure of neutron stars is also explained in a concisely way. 2.1 Historical overview It was not until 1932 when James Chadwick discovered the neutron. Only one year later, Walter Baade and Fritz Zwicky proposed the concept of neutron stars during their investigations of supernovae: a compact object formed by the collapse of the central core of the presupernova star [12]. This compact star made of "closely packed neutrons" would have to be degenerate and it could not generate energy after its formation; just only radiate its remaining heat through the slow process of photon diffusion. Nevertheless, due to the small size of neutron stars in comparison with white dwarfsi, they would be rarely observed by optical telescopes. Then, no searches for neutron stars were performed after Baade and Zick’s paper. A few years later, in 1939, Oppenheimer and Volkoff [14] and Tolman [15] derived the equations that describe the structure of a static star with spherical symmetry in General Relativity, and performed the first theoretical calculation of the equilibrium conditions of neutron stars and their properties assuming an ideal gas of neutrons at high density. In their calculation, they found that stable static neutron stars could not have masses larger than ∼0.75 M, a value much lower than the Chandrasekhar limit of white dwarfs ∼1.44 M. This low value is a consequence of the simple description they made of the state of matter of the neutron star interior in terms of non-interacting neutrons, and it is iStellar core remnants composed mostly of electron-degenerate matter, which represent the end point of stellar evolution for main-sequence stars with masses from 0.07 to 10M. The first white dwarf (40 Eridani B) was officially discovered in 1914 by Walter Adams [13]. 3 an indication that the role of the nuclear forces is fundamental to determine properly the structure of these objects. In the mid-1950s Wheeler and collaborators [16] built a more realistic neutron star matter EoS. They calculate the EoS for a non-interacting gas made of neutrons, protons and electrons under β-equilibrium conditions with respect to weak interaction processes. These authors developed a rigorous framework in which to understand compact stars, being their work of great importance even today. In 1959 the effect of the nucleon-nucleon interaction on the EoS and structure of neutron stars was studied for the first time, finding maximum mass values of Mmax ∼2M (Cameron [17]). At that time it started to be clear that, in addition to neutrons, protons and electrons, other particles such as muons, mesons or hyperons (baryons with strangeness content) could also be present in the interior of neutron stars. In spite of many efforts, the behaviour of the nuclear EoS at very high densities remained still uncertain. Another important theoretical step was the idea that neutron superfluidity could occur in neutron star interiors, that is, a singlet-state neutron pairing which operates at subnuclear densities in the inner crust of the neutron star but disappears in the core where the neutron-neutron interaction becomes repulsive [18]. The possibility of having in the core proton pairing in a singlet-state and neutron pairing in a triplet-state was understood later. It is well known that superfluidity in nuclear matter implies important consequences for cooling and other phenomena [19]. The main processes of a strict neutron star cooling theory seems to be the emission of neutrinos and photons in a later stage. Photon cooling consists on the photon diffusion mechanism as indicated before, but the neutrino emission processes have a stronger relevance on the cooling rates of neutron stars. Neutrinos are emitted in β-decay processes Eqs. (2.1), where thermally excited neutrons or protons undergo direct or inverse beta decay and the generated neutrinos carry away energy as they escape from the star. Consequently, the star cools and the number of thermally excited nucleons drops [20]. β−:n→p+e−+ ¯νe β+:p→n+e++νe (2.1) Neutrino cooling can be used to constraint the dense matter EoS by comparing the results of theoretical cooling models with the observation of the thermal radiation from neutron stars (Bahcall and Wolf [21,22]). It was not until 1967 that Pacini [23] showed that rapidly rotating neutron star with a strong dipole magnetic field could transform its rotational energy into electromagnetic radiation and accelerate particles to high energies, powering in this way a surrounding nebula, like the Crab nebula (named as Crab Pulsar later on in 1968 by Pacini). 4 After that, on August 1967, the first radio pulsar (named PSR B1919+21) was discovered by Bell and Hewish [24]; and the name of "pulsar" was selected due to their stable signals. However, it was not until the argument of Gold [25] that pulsars were thought to be highly magnetised rotating neutron stars. Furthermore, the simultaneous discoveries of Crab and Vela pulsars in 1968 (as mentioned before), both located in supernova remnants, also confirmed the prediction of Baade and Zwicky: neutron stars are formed in supernova explosions. Indeed, since 1968, understanding the properties of neutron stars was a matter of primordial importance and a lot of theoretical work was developed. In 1971, the discovery of pulsating compact X-ray sources (X-ray pulsars) stimulated this scenario even further. Those sources were believed to come from a neutron star in a close binary system which is accreting matter from its ordinary companion star [26]. Another important step in the history of neutron stars was the observation of the first binary pulsar (PSR J1913+16) by Hulse and Taylor in 1974 [27]. During the next years, several satellites were launched in order to observe neutron stars down to the last detail. A last historical surprise about neutron stars has been the very recent first direct detection of gravitational waves from the merger of two neutron stars in 2017 [7], which has shaken up the entire physics’ world anew. 2.2 Structure of neutron stars Neutron stars are supported against gravitational collapse mainly by the neutron degeneracy pressure, as mention in the previous section. Having masses of the order of 1-2M and radii within the range 10-13 km, an average density of the order of ∼1014 g/cm3is required. However, the internal structure of these objects can be described as an "onion" with a wide range of densities into their different shells (Figure 2.1). The atmosphere is the most external region and constitutes the thin layer, having a variation of its thickness from ∼10 cm in hot neutron stars to ∼mm in cold ones. The densities in the atmosphere are in the range n≤106g/cm3. It is a plasma layer where matter is formed by atoms of heavy nuclei, mainly nuclei around the iron mass number, distributed in a lattice in order to minimise the Coulomb energy. The observed thermal spectrum of the neutron stars comes from this region and, although the theoretical study of neutron stars atmospheres has been carried out by many authors, the current models are still incomplete due to some uncertainties. The main problems that need to be solved are the difficulty of the calculation of the EoS, the ionisation equilibrium and the spectral opacity of the atmospheric plasma (see ref. [28] for a detailed review). Then, we find the crust which is formed by the outer and inner crust, having a total thickness of ∼1 km and contributing only to a few percent of the total mass of the neutron star. 5 where the term d2E dx2 pis associated with the symmetry energy S(ρ)by the empirical mass formula and the higher-order derivatives are assumed to be small. From the definition of the chemical potential and imposing the known beta-equilibrium conditions, we are able to express proton and neutron chemical potentials in terms of the symmetry energy and proton fraction: µp=∂ε ∂ρp =∂(Eρ) ∂ρp =∂(Eρ) ∂ρ ∂ρ ∂ρp +∂(Eρ) ∂xp ∂xp ∂ρp =∂E ∂ρ ρ+E+ (1 −xp)∂E ∂xp µn=∂ε ∂ρn =∂(Eρ) ∂ρn =∂(Eρ) ∂ρ ∂ρ ∂ρn +∂(Eρ) ∂xp ∂xp ∂ρn =∂E ∂ρ ρ+E−xp ∂E ∂xp ⇒µn−µp=−∂E ∂xp = 4S(ρ)(1 −2xp).(3.17) Finally, other properties of interest in the characterisation of EoSs are the incompressibility modulus Kat non-zero pressure: K= 9∂P ∂ρ ,(3.18) and the sound speed vs, which depends also on the density of the nuclear medium thorugh the relation: vs c2 =dP dε =dP dρ dρ dε =K 9(mnc2+E+P/ρ).(3.19) The tracking of the dependence on density of vsis important, since a superluminal behaviour can occur at higher densities for some EoSs. 12 3.3 Nucleonic matter: H&H model Heiselberg and Hjorth-Jensen [11] detailed a simple parameterization of a non-relativistic nuclear matter model including three-body forces at temperatures much lower than the typical Fermi energies (T=0). The binding energy per nucleon in nuclear matter used by Heiselberg and Hjorth-Jensen (H&H) is adopted to be the sum of a compressional and a symmetry term: E=Ecomp(n) + S(n)(1 −2xp)2=E0uu−2−δ 1 + δu +S0uγ(1 −2xp)2,(3.20) being S(n)the symmetry energy and xpis the proton fraction, defined as xp=np/n. The ratio of the baryon density to nuclear saturation density corresponds to u=n/n0. The first term is parametrized in order to reproduce the saturation density, binding energy and compressibility of symmetric nuclear matter. The binding energy per nucleon at saturation density is taken to be E0=-15.8 MeV (without Coulomb energies). The introduction of δparameter is justified by applying the condition that the sound speed v2 s=∂P/∂ε does not exceed the speed of light. A value of δ=0.2 is chosen in order to avoid having a EoS that becomes superluminical at high densities. In agreement with the experimental value, the compressibility related to the first term is found to be K0=18E0/(1+δ)≃200 MeV. For the second term (the symmetry energy term), the best fit gives S0=32 MeV and γ=0.6. In neutron star matter, one needs to impose beta-stable conditions. In the H&H model, charge neutral uniform matter is made mainly of neutrons, protons, electrons and muons, and their composition is determined by the requirements of chemical and electrical equilibrium as exposed in Section 3.1. From Eqs. (3.3) and (3.11), the following conditions for matter in β-equilibrium with only nucleonic degrees of freedom are obtained: µn=µp+µenp=ne.(3.21) Notice that neutrinos are not included in the conditions for β-stable matter as already mentioned before. Moreover, assuming ultra-relativistic electrons (pe∼100 MeV), their chemical potential is found to be: µe=1 n ∂ε ∂xp =kFe= (3π2ne)1/3= (3π2np)1/3= 4S(n)(1 −2xp).(3.22) 13 The dependence of the electron chemical potential on the symmetry energy and proton fraction (last term of Eq. (3.22)) is found by using the expressions for the proton and neutron chemical potentials [Eq. (3.17)] from the previous section, keeping in mind the β-equilibrium conditions [Eq. (3.22)]. Therefore, from these results, one can conclude that the symmetry energy S(ρ)has a important role in studies of neutron star matter in beta-equilibrium. Therefore one is able to express the proton fraction in terms of the symmetry energy: nxp=(4S(n)(1 −2xp))3 3π2.(3.23) In the case of the H&H parametrization, the energy per particle for β-stable matter turns out to be given analitically by: E=E0uu−2−δ 1 + δu +S0uγ 2√a tan(2ψ)!2 ,(3.24) where a= 2(4S0uγ)3/π2nand tan ψ= (tan φ 2)1/3with tan φ=−2√a 3. 3.4 Nucleonic and Hyperonic matter: FSU2H model Tolós, Centelles and Ramos [10] started from the Relativistic Mean Field model of matter, where the baryons interact through the exchange of mesons, in order to provide a covariant description of the EoS and nuclear systems. The Lagrangian density can be expressed as the sum of the baryons (b), leptons (l) and mesons (m) Lagrangian densities: L=X bLb+Lm+X l=e,µ Ll,(3.25) with the lepton (electrons and muons) and mthe possible exchanged mesons (σ,ω,ρand φ). Their respective Lagrangians are given by 14 Lb=¯ Ψb(iγµ∂µ−qbγµAu−mb+gσbσ−gωbγµωµ−gφbγµφµ−gρbγµ~ Ib~ρ µ)Ψb, (3.26) Ll=¯ ψl(iγµ∂µ−qlγµAµ−ml)ψl, (3.27) Lm=1 2∂µσ∂µσ−1 2m2 σσ2−κ 3!(gσN σ)3−λ 4!(gσN σ)4−1 4ΩµνΩµν +1 2m2 ωωµωµ +ζ 4!(gωN ωµωµ)4−1 4~ Rµν ~ Rµν +1 2m2 ρ~ρµ~ρ µ+ Λωg2 ρN ~ρµ~ρ µg2 ωN ωµωµ −1 4PµνPµν +1 2m2 φφµφµ−1 4FµνFµν,(3.28) where Ψband Ψlare the baryon and lepton Dirac fields, respectively, while the mesonic and electromagnetic field strength tensors are Ωµν =∂µων−∂νωµ,~ Rµν =∂µ~ρν−∂ν~ρµ, Pµν =∂µφν−∂νφµand Fµν =∂µAν−∂νAµ. The meson couplings to a certain baryon are denoted by g(with Nindicating nucleon). The electromagnetic couplings are written as qand the baryon, meson and lepton masses by m. The isospin operator is denoted by the vector ~ Ib. The Dirac equations for baryons and leptons are presented below in which the effective baryon masses are set to m∗ b=mb−gσbσ, (iγµ∂µ−qbγ0A0−m∗ b−gωbγ0ω0−gφbγ0φ0−gρbI3bγ0ρ0 3)Ψb= 0, iγµ∂µ−qlγ0A0−mlψl= 0.(3.29) Now, following the Euler-Lagrange equations, the equations of motion for the meson fields in the mean-field approximation are expressed in terms of the respective meson mean-field expectation values (¯σ=< σ >,¯ρ=< ρ0 3>,¯w=< w0>and ¯ φ=< φ0>) as m2 σ¯σ+κ 2g3 σN ¯σ2+λ 3!g4 σN ¯σ3=X b gσbns b, m2 ω¯ω+ζ 3!g4 ωN ¯ω3+ 2Λωg2 ρN g2 ωN ¯ρ2¯ω=X b gωbnb, m2 ρ¯ρ+ 2Λωg2 ρN g2 ωN ¯ω2¯ρ=X b gρbI3bnb, m2 φ¯ φ=X b gφbnb,(3.30) where I3bis the third component of isospin of baryon bwith the convention I3p= +1/2. 15 The baryonic and leptonic Fermi energies in terms of their Fermi momenta, kF b and kFl, are described by EFb =qk2 Fb +m∗2 b, EFl =qk2 Fl +m2 l,(3.31) In this way, the scalar and vector densities for baryons band leptons llook like: ns b=m∗ b 2π2"EFbkF b −m∗2 bln kFb +EF b m∗ b#, nb=k3 Fb 3π2, nl=k3 Fl 3π2.(3.32) The energy density of the system is given by ε=X b εb+X l εl+1 2m2 σ¯σ2+1 2m2 ω¯ω2+1 2m2 ρ¯ρ2+1 2m2 φ¯ φ2+κ 3!(gσ¯σ)3 +λ 4!(gσ¯σ)4+ζ 8(gω¯ω)4+ 3Λω(gρgω¯ρ¯ω)2,(3.33) with the energy densities of baryons and leptons: εb=1 8π2kFbE3 Fb +k3 FbEF b −m∗4 bln kFb +EF b m∗ b, εl=1 8π2kFlE3 Fl +k3 FlEFl −m4 lln kFl +EF l ml.(3.34) The pressure of the system is computed using the thermodynamic relation P=X i µini−ε, (3.35) with the baryonic and leptonic chemical potentials given by 16 µb=EFb +gωb ¯ω+gρb I3b¯ρ+gφb ¯ φ , µl=EFl.(3.36) In this context, the parameters of the EoS are fitted with a view to fulfil the saturation properties of nuclear matter and finite nuclei, constraints of the high-density nuclear pressure from heavy-ion collisions and the 2Mobservations, while keeping small stellar radius in line with certain analysis [6]. The parameters are shown in the Table 3.1. mσ(MeV) 497.479 mw(MeV) 782.500 mρ(MeV) 763.000 g2 σN 102.7200 g2 wN 169.5315 g2 ρN 197.2692 κ4.0014 λ-0.013298 ξ0.008 Λw0.045 Table 3.1: Parameters of the FSU2H model from [10]. The mass of the nucleon is mN=939 MeV. In brief, the parametrization of FSU2H model leads to obtain the following properties at saturation density (n0=0.1505 fm−3). The energy per particle results in E/A=-16.28 MeV, the compressibility is K=238.0 MeV and the effective nucleon mass in symmetric nuclear matter leads to m∗ N/mN=0.593. The symmetry energy is set to Esym=30.5 MeV whose slope and curvature are defined by L=44.5 and Ksym=86.7 MeV, respectively. Finally, the pressure of pure neutron matter is PPNM =2.30 MeV fm−3at n0. As the interior of neutron stars stays globally neutral and in conditions of β-equilibrium [Eqs. (3.1) & (3.2)], therefore, the chemical potentials and the number densities of each particle in a neutron star are related to the beta-conditions already explained [see Eqs. (3.10) & (3.11)]. Finally, solving the Dirac equations (3.29) for the baryons and leptons with the mesonic field equations (3.30) for a given total baryon density nand obtaining the corresponding chemical potential from each species, one can determine the energy density and pressure of the β-stable neutron star as functions of density. 17 3.5 Equations of state In this section, we are going to represent the EoSs for the H&H and FSU2H models employed in this study, in order to analyse them and distinguish their characteristics. Indeed, Figures 3.1 and 3.2 show the energy density and pressure of the system, respectively, as functions of the baryon density for the H&H and FSU2H models. 0 2 4 6 8 10 12 14 0 1 2 3 4 5 6 7 ε[102MeV/fm3] ρ/ρ0 H&H FSU2H Figure 3.1: Energy density vs. baryon density for the different models presented in the text: H&H (Heiselberg & Hjorth-Jensen [11]) and FSU2H (Tolós, Centelles & Ramos [10]). 0 1 2 3 4 5 6 0 1 2 3 4 5 6 7 P[102MeV/fm3] ρ/ρ0 H&H FSU2H Figure 3.2: Pressure vs. baryon density for the different models presented in the text: H&H (Heiselberg & Hjorth-Jensen [11]) and FSU2H (Tolós, Centelles & Ramos [10]). 18 On the one hand, from Figure 3.2, one can extract some conclusions. The FSU2H model has a stiffer behaviour in the range n0≤n≤5n0than the H&H model, although the appearance of hyperons (∼2n0) softens the pressure of the FSU2H system so that the H&H EoS becomes stiffer from ≃5n0onwards. The change in the pressure conduct due to the introduction of hyperons can be explained in terms of the Fermi pressure, which reduces with the presence of this new species that starts to fill up its own Fermi sea. For that reason, the H&H model gives a harder EoS for larger densities. The lack of new species prevents the reduction of the pressure and the electrons, protons and neutrons continue occupying higher energy states of their own Fermi seas, rising the pressure of the neutron star. In this context, greater pressure values mean having a system harder to compress and, accordingly, it is able to resist more against gravity. In fact, a stronger bearing at saturation densities is translated into larger neutron star radii and a stiffer conduct at suprasaturation densities is marked by higher masses. On the other hand, the FSU2H energy density in function of the baryon density is larger than the H&H EoS in all the density range explored in Figure 3.1. Finally, Figure 3.3 represents the dependence of the pressure on the energy density for both models in study. As expected this plot shows the same behaviour sustained by Fig. 3.2 as both pressure and energy density increase monotonically with density. We show Figure 3.3 because pressure and energy density enter in the structure equations of neutron stars, presented in the next chapter. 0 1 2 3 4 5 6 0 2 4 6 8 10 12 14 P[102MeV/fm3] ε[102MeV/fm3] H&H FSU2H Figure 3.3: Pressure vs. energy density for the different models presented in the text: H&H (Heiselberg & Hjorth-Jensen [11]) and FSU2H (Tolós, Centelles & Ramos [10]). 19 Chapter 4 Non-rotating neutron stars: TOV equations In order to determine the structure of a neutron star, Einstein’s General Relativity Theory is needed as we are dealing with one of the densest objects in the Universe. For this reason, the equations that describe a non-rotating star in equilibrium (or Schwarzschild star) start from Einstein’s field equations. Using Gravitational Units (G=c=1)i, the Einstein’s field equations are given by Gµν =−8πTµν,(4.1) where Gµν is the so-called Einstein tensor and Tµν is the energy-momentum tensor. The solution of the Einstein’s equations for a static, spherically symmetric and relativistic star are known as the Tolman-Oppenheimer-Volkoff equations (TOV equations) [14, 15]. For a given EoS, and a given value of the central density, the stellar structure equations of hydrostatic equilibrium in General Relativity are: dp(r) dr =−(ε(r) + p(r))(M(r)+4πr3p(r)) r(r−2M(r)) (4.2) and dM(r) dr = 4πr2ε(r).(4.3) By re-writing Equation (4.2), we arrive to dp(r) dr =−M(r)ε(r) r21 + p(r) ε(r)1 + 4πr3p(r) M(r) 1−2M(r) r.(4.4) iWe are going to use Gravitational Units from now to the end of this project. 20 The equation can be interpreted as follows [1]. On the one hand, the left-hand side of Eq. (4.4) is the net force acting outward on the surface of the shell by the pressure difference dp(r). On the other hand, the first factor of the right-hand side of Eq. (4.4) is the attractive gravity Newtonian force acting on the shell by the mass enclosed. The second factor on the right-hand side of the equation refers to the corrections from General Relativity. Thus, these equations express the balance at each rbetween the internal pressure and the gravitational attraction of the mass enclosed at r. Moreover, the integration of the pressure p(r)and the enclosed mass to a given radius, M(r), must be performed outward; starting at the star’s centre (r=0) where the mass of the star is M=0 and the pressure is pC. We finish when the pressure becomes p= 0, which defines the edge of the star. The point at which the pressure vanishes defines the radius of a non-rotating star, R. The integration of equation (4.3) from the centre of the star up to its radius defines its gravitational mass: MG= 4πZR 0 drr2ε(r).(4.5) 21 It must be simultaneously integrated together with the following equation dv2(r) dr =−2dν(r) dr h2(r)+1 r+dν(r) dr  −r3 3 dj(r)2 dr ¯w(r)2+j(r)2 6r4d¯w(r) dr 2!,(5.21) with boundary conditions v2(0)=h2(0)=0 and v2(∞)=h2(∞)=0. To sum up, we have shown all necessary expressions that determining the structure of a rotating neutron star. 5.2 Numerical calculation In this section, we explain how to deal with the equations that permit us to determine the deformation of rotating neutron stars. We have developed a Fortran programme within the Hartle-Thorne approach that is shown in Appendix A.1. First of all, we read our data from the EoS and solve the TOV equations [Eqs. (4.2) & (4.3)] in order to know the gravitational mass and radius for the static solution. We start the integration from the origin having as initial conditions M(0)=0 and a particular value of the central pressure p(0) = pc, until it becomes zero with Rthe radius of the star and M(R)its gravitational mass. The integration of the metric function ν(r)[Eq. (5.5)] is done together with the TOV equations, selecting the convenient value of ν(0) = 0. As it is figured out, since M(r)and p(r)are the variables of the TOV equations, one can obtain ν(r)simultaneously. Once we have the result for the metric function ν(r), we find the correct condition of the metric function at Rby applying Eq. (5.7). We also work out the j(r) function [Eq. (5.3)] and the relative dragging frequency [Eq. (5.2)], being the arbitrary constant value fixed to ¯wc= 1 for convenience in order to accomplish the boundary conditions, and, finally, the total angular momentum [Eq. (5.9)]. Now, using Equation (5.8), we are able to find the value for the total angular velocity of the star, named as Ωarb. As we are interested to chose the total angular momentum of the star Ω, one must scale the relative angular velocity ¯w(r)as ¯wconst =Ω Ωarb .(5.22) Consequently, one should apply the rescale not only to the relative angular velocity but also to the respective derivative and total angular momentum of the star. 28 Furthermore, we are able to calculate the moment of inertia [Eq. (5.10)] and the Keplerian frequency ΩK. We use the approximation [12] ΩK= 0.65 M R31/2 .(5.23) Next, we begin the integration for the monopole mass [Eq. (5.14)], the monopole pressure [Eq. (5.15)], the h2function [Eq. (5.20)] and v2function [Eq. (5.21)] with the respective boundary conditions already mentioned. Once the integration has finished, we are able to get the increment of mass due to rotation given in Eq. (5.13). Lastly, we can do the integration for the quadrupole pressure [Eq. (5.19)] and the stretching functions [Eq. (5.18)], which allow us getting the polar and equatorial radii [Eq. (5.17)] and the star eccentricity [Eq. (5.16)]. All these steps are repeated for every star, that means for every central density given in the EoS data file. 29 Chapter 6 Results In this chapter, we are going to show the results obtained after solving the TOV equations for non-rotating neutron stars and, afterwards, the ones obtained after applying the Hartle- Thorne approach of rotating neutron stars, which is the main interest of this project. The results for rotating neutron stars will be classified according to the deformation in size and mass due to rotation and, the moment of inertia associated to those rotating neutron stars. The calculations have been performed for the nucleonic H&H and the hyperonic FSU2H EoSs, for different rotating velocities (Ω=1000, 3000 & 5000 Hz). Those values have been chosen taking into account that pulsars rotating with periods bigger than 700 Hz (∼4500 s−1) have not been observed yet, despite the fact that it does not exist an upper limit banning a more rapid rotation. Therefore, we will compare the results of the two models for the selected rotating velocities and try to understand their differences. 6.1 Non-rotating neutron stars First of all, it is important to understand the features of non-rotating neutron stars. For that reason, in Figure 6.1, we show the mass and radius relation obtained for both EoSs after solving the TOV equations (eqs. (4.2) & (4.3)). For stars with approximately the same mass, we denote a substantial difference between their radii. The FSU2H model gives rise to the larger radii as compared to the H&H EoS, a fact that is related to its stiffness (see Table 6.1). Even though, both EoSs give radii within the estimated range (10.7-13.1 km). Consequently, it turns out that only the FSU2H model achieves the maximum mass limit of 2Mgiven by observations. The maximum non-rotating neutron star masses are found to be of 1.94Mfor H&H EoS and 2.02Mfor FSU2H model, with the corresponding radii of R=10.21 km and R=12.13 km, as seen in Table 6.1. 30 In brief, a stiffer EoS leads to larger masses but also bigger stellar radii. M/MR(km) ρc(fm−3)εc(MeV/fm3)Ωk(s−1) H&H 1.04 12.06 0.40 400 5773 1.40 11.91 0.52 537 6815 1.90 10.77 0.97 1153 9239 1.94 10.21 1.21 1555 10092 FSU2H 1.02 13.04 0.28 276 5088 1.39 13.29 0.33 327 5761 1.90 13.14 0.53 573 6843 2.02 12.13 0.87 1040 7971 Table 6.1: Mass (M) and radius (R) obtained for H&H and FSU2H models, together with the corresponding central density and energy density. The Keplerian frequency is also shown for each case. 0 0.5 1 1.5 2 2.5 10 11 12 13 14 15 16 M/Ms R[km] H&H FSU2H Figure 6.1: Plot of the mass and radius relation for the FSU2H and H&H EoSs for non-rotating neutron stars. 31 6.2 Rotating neutron stars Before starting the analysing of the rotating neutron stars results, it is of great utility to understand first the relation of our non-rotating neutron stars parameters with the Keplerian frequency. The Keplerian frequency indicates the upper limit on the rotation of stable neutron stars. Therefore, calculating the Keplerian frequency of a star of mass Mand radius Rwill establish the maximum angular velocity that is capable to support before becoming unstable. Thus, the Keplerian frequency informs us about the angular velocities that our neutron stars can hold up in the rotation. Figure 6.2 illustrates the dependence of the Keplerian angular velocity ΩKon the gravitational neutron star mass [Eq. (5.23)], and consequently, on the underlying EoS also. The Keplerian frequency values increase as the gravitational mass of the stellar sequence does. 0 1000 2000 3000 4000 5000 6000 7000 8000 0 0.5 1 1.5 2 ΩK[s-1] M/Ms H&H FSU2H Figure 6.2: Keplerian angular velocity vs gravitational mass of non-rotating neutron stars for different EoSs. Referring to the EoSs in study, we oserve that the two different models have led to quite different values for ΩKas a function of M. This difference highlights with the rise of the gravitational mass. According to Glendenning [12], the dependence of Pon ε(Fig. 3.3) has an important impact on the Keplerian frequency. The EoS must be soft at low and intermediate densities and hard at high densities for a neutron star to withstand fast rotation. Such behaviour agrees with the nucleonic H&H model. 32 Consequently, the largest values for ΩK(M)result from the H&H EoS; the reason of its behaviour lies in the smaller values of the star radius that this model produces. Thus, the H&H model will allow higher rotational frequencies for a given mass value. Lastly, the three horizontal dashed lines represent the three angular velocities that will be employed in our calculations of rotating neutron stars (Ω=1000, 3000, 5000 s−1), setting a limit on the neutron stars that have smaller Keplerian frequency value. Therefore, the neutron stars with Keplerian angular velocities ΩKbelow these limiting lines will not be taken into account in the following, as they have arrived at their critical equatorial velocity at which mass-shedding sets in. 6.2.1 Deformation Firstly, we represent the dependence of the equatorial radius (Req) with the central energy density for different angular velocities, together with the non-rotating case, for both models in study. Figures 6.3 and 6.4 represent the results obtained with the H&H and FSU2H EoSs, respectively. For the non-rotating case, the relation of the radius with the central energy density gives us an idea about the compactness of the star. Not only the increase of the central energy density is linked with a decrease of the radius but also with bigger gravitational masses. In consequence, the size of the star is reduced with increasing central energy densities, leading to more compact neutron stars. Comparing the results of the H&H model with those of the FSU2H one for the same radius, we see that the nucleonic EoS gives a smaller central energy density and gravitational mass, outcoming into more diluted neutron stars than the FSU2H ones. Moreover, the FSU2H model shows a "wiggle" for ε∼2-4·102MeV/fm3which indicates that the same radius is obtained for different central energy densities, meaning different mass configurations and, on consequence, distinct stars. This is clearly seen in Figure 6.1, where the mass-radius plot for non-rotating neutron stars is shown. Once stars rotate, stars with the same central energy density but rotating faster become more extended due to the equatorial radius growth. Moreover, the deformation of the star is more evident for those stars with lower central energy densities, that is, more diluted stars. Lastly, comparing both studied models, we observe that FSU2H gives rise to larger equatorial radii values due to its harder EoS for densities around saturation, as previously discussed. 33 10 11 12 13 14 15 16 0 2 4 6 8 10 12 14 Req[km] εc[102MeV/fm3] H&H Ω=0s-1 Ω=103s-1 Ω=3*103s-1 Ω=5*103s-1 Figure 6.3: Equatorial radius vs. central energy density for different angular frequencies Ωfor the H&H EoS. 10 11 12 13 14 15 16 0 2 4 6 8 10 12 14 Req[km] εc[102MeV/fm3] FSU2H Ω=0s-1 Ω=103s-1 Ω=3*103s-1 Ω=5*103s-1 Figure 6.4: Equatorial radius vs. central energy density for different angular frequencies Ωfor the FSU2H EoS. 34 The deformation of the star from its spherical form due to rotation can be also shown by plotting the relation between the polar and equatorial radii. Figures 6.5 and 6.6 show this relation for the nucleonic H&H and the hyperonic FSU2H models, respectively, for different angular velocities. In accord with the previous conclusions for Figures 6.3 & 6.4, the higher the value of the equatorial radius, the more deviation of the star from its spherical shape. Thus, stars with smaller gravitational masses and, on consequence, bigger non-rotating radii (lower values of εc) suffer larger deformations on their aspect than the more massive stars. Moreover, the deviation from the spherical configuration increases as the angular velocity does, meaning having more oblate neutron stars for faster rotating configurations. By comparing the results for both models, the representation of polar and equatorial radii seems to give a more extensive deformed aspect for the star configurations of the FSU2H EoS, in concordance with the last conclusions of Figures 6.3 and 6.4. This behaviour will be illustrated more clearly by the observables represented in the following Figures 6.7 and 6.8. Notice that the number of stars decreases with the implementation of higher rotation. This is due to the fact that the ending point on the right of each angular velocity function is the last star which fulfils the Keplerian condition, and the higher the rotation, the faster this limit of instability is reached. Finally, in Figures 6.7 and 6.8, we show the eccentricity as a function of the rotational mass for the different angular velocities and both EoSs. This observable also illustrates the deformation of the star due to rotation but by the description of the star’s surface [Eq. (5.16)]. For a fixed frequency, the deformation is larger (higher eccentricity) for the less massive stars as this systems are characterised by a smaller Keplerian frequency. Thus, the less massive stars have a larger size (bigger radii and lower central energy density) because they are rotating closer to their mass-shedding limit. Consequently, those systems are more deformed from their original non-rotating spherical configuration. With the growth of the angular velocity Ω, the eccentricity increases its value approaching unity, which is in agreement with the rise of the degree of deformation of the equatorial and polar radii. In brief, the deformation increases as the system is getting closer to its Keplerian frequency. Now, by comparing both figures, since the FSU2H model is characterised by smaller Keplerian frequencies, it suffers larger deformations due to rotation than the H&H, thereby showing higher eccentricity values. This behaviour is also supported by our previous discussion of deformation. Note that eccentricity values very close to unity are not obtained as those values would come from rotating stars with frequencies above the Keplerian limit, and hence, unstable. 35 9 10 11 12 13 14 15 16 17 9 10 11 12 13 14 15 16 17 Rp[km] Req[km] H&H Ω=0s-1 Ω=103s-1 Ω=3*103s-1 Ω=5*103s-1 Figure 6.5: Polar vs. equatorial radius for different angular frequencies Ωfor the H&H EoS. 9 10 11 12 13 14 15 16 17 9 10 11 12 13 14 15 16 17 Rp[km] Req[km] FSU2H Ω=0s-1 Ω=103s-1 Ω=3*103s-1 Ω=5*103s-1 Figure 6.6: Polar vs. equatorial radius for different angular frequencies Ωfor the FSU2H EoS. 36 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 e=0.08 e=0.24 e=0.41 e Mrot/Ms H&H Ω=103s-1 Ω=3*103s-1 Ω=5*103s-1 Figure 6.7: Eccentricity vs. rotational mass for different angular frequencies Ωfor the H&H EoS. 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 e=0.10 e=0.30 e=0.49 e Mrot/Ms FSU2H Ω=103s-1 Ω=3*103s-1 Ω=5*103s-1 Figure 6.8: Eccentricity vs. rotational mass for different angular frequencies Ωfor the FSU2H EoS. 37 0 0.5 1 1.5 2 2.5 3 8 9 10 11 12 13 14 15 16 M=1.5Ms I[1045g*cm2] Req[km] H&H Ω=0s-1 Ω=103s-1 Ω=3*103s-1 Ω=5*103s-1 Figure 6.17: Moment of inertia vs. equatorial radius for different angular frequencies Ωfor H&H EoS. 0 0.5 1 1.5 2 2.5 3 8 9 10 11 12 13 14 15 16 M=1.5Ms I[1045g*cm2] Req[km] FSU2H Ω=0s-1 Ω=103s-1 Ω=3*103s-1 Ω=5*103s-1 Figure 6.18: Moment of inertia vs. equatorial radius for different angular frequencies Ωfor FSU2H EoS. 44 Chapter 7 Conclusions The main aim of this project has been to study the rotation of neutron stars and to understand how it could affect to their mass, radius and moment of inertia. We have analysed the results obtained by the recently developed hyperonic FSU2H EoS and compared them with those of the nucleonic H&H model, being able to distinguish the differences in behaviour depending on the presence or absence of hyperons in the neutron star composition and, thus, arriving to diverse conclusions. On the one hand, from the non-rotating neutron stars, one is able to foresee the behaviour due to rotation via the analysis of the EoS. A stiffer EoS is connected with characteristic higher gravitational masses and radii and, in consequence, smaller Keplerian frequencies. In this context, the Keplerian frequency sets the limit of angular velocity that a neutron star is capable to withstand. Thus, the FSU2H model due to its harder behaviour gives a sign of allowing lower rotational frequencies, leading of bigger affectations due to rotation. On the other hand, from the rotating neutron stars, one is interested in the specific changes to their properties, such as the mass and radius. Thereby, rotating neutron stars characterised by smaller masses suffer larger deformations in radii than the more massive ones due to their lower Keplerian frequencies, leading to a more oblate configuration. Moreover, the faster the rotation, the larger the deformation and the increment of mass. Referring to the moment of inertia of rotating neutron stars, the systems described by higher masses lead to larger values for the moment of inertia. When increasing the rotational frequency, we obtained a less intuitive result at first sight: the increment of the angular velocity leads to lower moment of inertia values, in spite of an increase of the equatorial radius. Therefore, by comparing the H&H and FSU2H models, we have found that the hyperonic FSU2H model reaches a more oblate configuration due to rotation than the H&H EoS. This fact was previously predicted by its stiffer behaviour close to saturation densities that leads to larger radii deformations due to its distinctive smaller Keplerian frequencies. Additionally, the FSU2H EoS gives rise to more massive non-rotating neutron stars, larger rotational masses and, hence, higher moments of inertia. 45 Finally, in future projects, it would be interesting to calculate the tidal deformability of the rotating hyperonic neutron stars obtained. Furthermore, one must be intrigued from the very recent developed FSU2H model including ∆particles [34] as its analysis could be of great importance in understanding how neutron stars are affected to rotation. Thereby, we would be able to compare the rotating ∆neutron stars with the hyperonic ones already obtained. 46 Chapter 8 Bibliography [1] Isaac Vidaña. A short walk through the physics of neutron stars, 2018. [2] V. Kalogera and G. Baym. The maximum mass of a neutron star. The Astrophysical Journal, 470(1):L61–L64, Oct 1996. [3] P. B. Demorest, T Pennucci, S. M. Ransom, M. S. E. Roberts, and J. W. T. Hessels. A two-solar-mass neutron star measured using Shapiro delay. Nature, 467:1081–1083, 10 2010. [4] J Antoniadis, P. C. C. Freire, N Wex, T. M. Tauris, R. S. Lynch, M. H. van Kerkwijk, M Kramer, C Bassa, V. S. Dhillon, T Driebe, J. W. T. Hessels, V. M. Kaspi, V. I. Kondratiev, N Langer, T. R. Marsh, M. A. McLaughlin, T. T. Pennucci, S. M. Ransom, I. H. Stairs, J van Leeuwen, J. P. W. Verbiest, and D. G. Whelan. A Massive Pulsar in a Compact Relativistic Binary. Science, 340:448, 4 2013. [5] H. Cromartie et al. A very massive neutron star: relativistic Shapiro delay measurements of PSR J0740+6620. 2019. [6] J. M. Lattimer and M. Prakash. The Equation of State of Hot, Dense Matter and Neutron Stars. Phys. Rept., 621:127–164, 2016. [7] B. P. Abbott et al. GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral. Phys. Rev. Lett., 119(16):161101, 2017. [8] Ben Margalit and Brian D. Metzger. Constraining the maximum mass of neutron stars from multi-messenger observations of gw170817. The Astrophysical Journal, 850(2):L19, Nov 2017. [9] S. Woosley and T. Janka. The physics of core-collapse supernovae. Nature Physics, 1(1):147–154, Dec 2005. [10] L. Tolos, M. Centelles, and A. Ramos. The Equation of State for the Nucleonic and Hyperonic Core of Neutron Stars. Publ. Astron. Soc. Austral., 34:e065, 2017. [11] M. Heiselberg, H.and Hjorth-Jensen. Phases of dense matter in neutron stars. Phys. Rept., 328:237–327, 2000. 47 [12] [Glendenning]_Compact_stars(BookFi). [13] Wikipedia. White Dwarf. [14] J. R. Oppenheimer and G. M. Volkoff. On massive neutron cores. Phys. Rev., 55:374– 381, Feb 1939. [15] R. C. Tolman. Static solutions of einstein’s field equations for spheres of fluid. Phys. Rev., 55:364–373, Feb 1939. [16] B.K. Harrison, M. Wakano and J.A. Wheeler. Matter-energy at high density: end point of thermonuclear evolution. La structure et évolution de l’univers, R. Stoops, Bruxelles, 1958. [17] A. G. W. Cameron. Pycnonuclear reations and nova explosions. The Astrophysical Journal, 130:916, Nov 1959. [18] R. A. Wolf. Some Effects of the Strong Interactions on the Properties of Neutron-Star Matter. The Astrophysical Journal, 145:834, Sep 1966. [19] D. Pines and M. A. Alpar. Superfluidity in neutron stars. Nature, 316:27–32, Jul 1985. [20] J. M. Lattimer. Viewpoint: A Rapidly Cooling Neutron Star. [21] J. N. Bahcall and R. A. Wolf. Neutron stars. ii. neutrino-cooling and observability. Phys. Rev., 140:B1452–B1466, Dec 1965. [22] J. N. Bahcall and R. A. Wolf. An Observational Test of Theories of Neutron-Star Cooling. The Astrophysical Journal, 142:1254–1256, October 1965. [23] F. Pacini. Energy Emission from a Neutron Star. Nature, 216:567–568, Nov 1967. [24] Bell S. J. Pilkington J. D. H. Scott P. F. Collins R. A. Hewish, A. Observation of a Rapidly Pulsating Radio Source. Nature, 217:709–713, Feb 1968. [25] T. Gold. Rotating Neutron Stars and the Nature of Pulsars. Nature, 221:25–27, Jan 1969. [26] S. Rappaport and P. C. Joss. Binary x-ray pulsars. Nature, 266:123–125, Mar 1977. [27] R. A. Hulse and J. H. Taylor. Discovery of a pulsar in a binary system. The Astrophysical Journal Letters, 195:L51–L53, January 1975. [28] D G Yakovlev P. Haensel, A.Y. Potekhin. Neutron Stars 1. Equation of State and Structure, volume 326. Springer, 2007. [29] J. A. Pons, D. Vigano, and N. Rea. Too much "pasta" for pulsars to spin down. Nature Phys., 9(arXiv:1304.6546):431–434, Apr 2013. [30] W. G. Newton. A taste of pasta?, Jun 2013. Nature Publishing Group, a division of Macmillan Publishers Limited. All Rights Reserved. 48 [31] J. B. Hartle. Slowly rotating relativistic stars. i. equations of structure. Astrophys. J., 150:1005, 1967. [32] J. B. Hartle and K. S. Thorne. Slowly Rotating Relativistic Stars. II. Models for Neutron Stars and Supermassive Stars. Astrophys. J., 153:807, September 1968. [33] M. Salgado, S. Bonazzola, E. Gourgoulhon, and P. Haensel. High precision rotating netron star models 1: Analysis of neutron star properties. , 291:155–170, November 1994. [34] P. Ribes, A. Ramos, L. Tolos, C. Gonzalez-Boquera, and M. Centelles. Interplay between delta particles and hyperons in neutron stars. The Astrophysical Journal, 883(2):168, Oct 2019. 49 Appendix A Appendix A.1 Programme c... Program that solves the TOV equations for a given EoS and the Hartle-Thorne approach for rotating neutron stars c... Determines M-R diagram c... Gives Mass profile for a couple of stars c... Gives rotation parameters and characterises the rotational deformation of the star c... (Adpated by Anna Campoy from Àngels Ramos and Laura Tolós code of non-rotating neutron stars) IMPLICIT NONE integer i,ncrust,ncore,ni,nf,nmax,nstep,k,ipre,ndiscard,imax integer column parameter (nmax=5000) real*8 pai,pai2,conv1,conv2,conv3,conv4,conv5,conv6,conv7 real*8 pmn,msolar,mmax,rmax real*8 rho(nmax),en(nmax),pre(nmax),eder(nmax) real*8 nur(100000),w(100000),mass(100000),press(100000) real*8 functh2(100000),functv2(100000) real*8 eps0(100000),eps2(100000),press2(100000) real*8 derj(100000),derw(100000),dernu(100000) real*8 m1,p1,en1,dens1,fntv1,dx,nu0,w0,w1,fntj1,mass1,press1 real*8 functv2_1, functh2_1,v2,h2,eder0,eder1 real*8 derw1,derj1,dernu1 real*8 ma,prener,ener,x,dens,fntv,fntvR,fntj,mass0,press0 real*8 dm,dp,dv,fac1,fac2,jintx,dmass0,dpress0,dh2,dv2 real*8 omega,wconst,omega_arb,Jsum,inmom,massinc real*8 polarR,equatR,eccent,kepler real*8 xx real*8 premin 50 c... input files open(unit=7,file=’EoS_crust_FSU2H_reduced.dat’) c... output files open(unit=8,file=’M_vs_R_npe.dat’) open(unit=9,file=’profile_npe1.dat’) open(unit=10,file=’profile_npe2.dat’) open(unit=12,file=’profile_npe3.dat’) open(unit=11,file=’derederp.dat’) write(11,’(a)’) ’#rho e p de/dp’ open(unit=13,file=’M_vs_R_rot.dat’) open(unit=14,file=’M_vs_R_rot-monopole.dat’) open(unit=15,file=’M_vs_R_rot-def.dat’) open(unit=17,file=’w(i).dat’) open(unit=18,file=’usefull_parameters.dat’) c... conversion factors pai=3.14159 pai2=pai*pai mmax=0.d0 rmax=0.d0 conv1=1.3234d-16 !(cm^(-2))*(fm^3/MeV) [de Mev/fm^3 a cm^(-2)] conv2=1.d-5 !km/cm [de cm a Km] conv3=1.d39 !fm^3/cm^3 [de fm^-3 a cm^-3] conv4=7.4237d-29 ! [g/cm^3 a cm^(-2)] (1g a cm) conv5=2.9979d10 !1s a cm conv6=4.0383d38 ! conv5/conv4 [cm a g.cm/s] conv7=1.d-15 !cm3 a km3 pmn=939.566d0 !MeV msolar=1.4766d5 !cm omega=0./conv5 !cm^(-1) column=0 c... write in file radius and mass write(8,’(2a)’) ’#r(km), M/Msol, rho_c(fm-3),’, & ’ E_c(MeV/fm^3), nu (adim), index’ c... write in file metric function, j function, dragging angular velocity, arbitrary angular velocity, total angular momentum and inertial momentum write(13,’(3a)’) ’#r(km), M/Msol, nu(adim),’, & ’ j(adim), w(s-1), Omega_arb(s-1),’, & ’ J (g*cm^2/s), I(g*cm^2), kepler(s^-1) index’ 51 c... write in file monopole mass and pressure, v2 and h2 functions, and the increase of mass due to rotation write(14,’(3a)’) ’#r(km), M/Msol, monopole mass(g)’, & ’ monopole pressure, v2 function,’, & ’ h2 function, massinc(g), index’ c... write in file quadrupole pressure, spherical and quadrupole stretching functions, polar and equatorial radius, and eccentricity write(15,’(3a)’) ’#r(km), M/Msol, quadr pressure,’, & ’ sph stretch(cm), quadr stretch(cm), polar R(km),’, & ’ equatorial R(km), eccentricity, index’ c... write in file rel. ang. velocity and Kepler frequency write(17,’(2a)’) ’#r(km), M/Msol, w(i)/omega,’, & ’ kepler(s^-1), index’ c... write in file rotational mass, polar and equatorial radius, Kepler freq., eccentricity, central energy density and moment of inertia write(18,’(3a)’) ’ #r(km), M/Msol, M/Msol rot,’, & ’ Req(km), Rpol(km), kepler(s^-1), eccent,’, & ’ E_c(MeV/fm^3), I(g*cm^2), index’ c... read files of density, energy density and presure (crust and core) ncrust=129 ncore=108 nf=ncrust+ncore !total number of points read(7,*) !read the points of the EoS file do i=1,nf read(7,*) xx,en(i),pre(i) !rho: fm^-3, en: MeV/fm^3, pre: MeV/fm^3 rho(i)=xx*conv3 !cm^(-3) c... write(*,*) rho(i)/conv3,en(i),pre(i) en(i)=en(i)*conv1 !cm^(-2) pre(i)=pre(i)*conv1 !cm^(-2) enddo c... computing derivative of energy density w.r.t pressure density do i=1,nf-1 eder(i)=(en(i+1)-en(i))/(pre(i+1)-pre(i)) if(pre(i+1).eq.pre(i)) eder(i)=eder(i-1) c if(i.eq.ncrust) eder(i)=eder(i-1) enddo eder(nf)=eder(nf-1) !set last value equal of the derivative to the previous one 52 write(11,’(4g20.12)’) rho(1)/conv3,en(1)/conv4,pre(1)/conv4, & eder(1) do i=2,nf-1 if(i.eq.83.or.i.eq.ncrust.or.eder(i).lt.0.d0) & eder(i)=eder(i-1)+(eder(i+1)-eder(i-1))*(rho(i)-rho(i-1))/ & (rho(i+1)-rho(i-1)) write(11,’(4g20.12)’) rho(i)/conv3,en(i)/conv4,pre(i)/conv4, & eder(i) enddo write(11,’(4g20.12)’) rho(nf)/conv3,en(nf)/conv4,pre(nf)/conv4, & eder(nf) c stop premin=pre(1) c... step for integration of Euler method nstep=100000 !with step of 1 meter, we go up to 100 Km dx=100 !cm (step of 1 meter) c... which initial/central pressure to use to calculate the corresponding mass-radius ni=ncrust+2 !where to start c... loop over different stars (with different central densities) do k=ni,nf c do k=213,213 c... m1: initial mass; p1: initial pressure; en1: initial energy density; dens1: initial density; fntv1: initial nu function m1=0.d0 !initial mass p1=pre(k) !initial pressure en1=en(k) !initial energy density dens1=rho(k) !initial density fntv1=0.d0 !initial metric function c... write(*,*) k, p1/conv1,en1/conv1,dens1/conv3 c... solve differential equations with Euler method c... starting values ma=m1 prener=p1 ener=en1 dens=dens1 fntv=fntv1 x=1 !cm 53 press0=press0+dpress0 h2=h2+dh2 v2=v2+dv2 mass(i)=mass0 ! defining monopole mass array press(i)=press0 ! defining monopole pressure array functv2(i)=v2 ! defining v2 function array functh2(i)=h2 ! defining h2 function array c... write(*,*) ma/msolar,prener/conv4,x*conv2,fntv c... determine index of new pressure within the pressure array do ipre=1,nf-1 if(((prener-pre(ipre))*(prener-pre(ipre+1))).lt.0.) goto 31 enddo 31 continue c... linear interpolation (to determine energy and density at the new position) ener=en(ipre)+(en(ipre+1)-en(ipre))/(pre(ipre+1)-pre(ipre))* & (prener-pre(ipre)) c... linear interpolation (to determine energy derivative w.r.t. pressure at the new position) eder0=eder(ipre)+ & (eder(ipre+1)-eder(ipre))/(pre(ipre+1)-pre(ipre))* & (prener-pre(ipre)) c... write interesting parameters in order to make some figures write(17,’(4g25.16,i5)’) x*conv2,ma/msolar,w(i)/omega, & kepler*conv5,k c print *, ’ mass0=’,mass0 c print *, ’press0=’,press0 enddo c... increase of mass of the star due to rotation massinc=mass(imax)+Jsum**2/x**3 ! cm write(14,’(7g14.6,i5)’) x*conv2,ma/msolar,mass(imax)/conv4, & press(imax),functv2(imax),functh2(imax), & massinc/msolar,k c... outside the star, relative angular velocity behaves as: do i=imax+1,imax+8000 x=x+dx w(i)=omega-2.*Jsum/x**3 write (17,*) x*conv2,column,w(i)/omega 60 enddo c ... initialising again the loop m1=0.d0 !initial mass p1=pre(k) !initial pressure en1=en(k) !initial energy density c... write(*,*) k, p1/conv1,en1/conv1,dens1/conv3 c... solve differential equations with Euler method c... initial conditions ma=m1 prener=p1 ener=en1 x=1 !cm c ... initialising integration loop in order to solve eqs 49 & 50 do i=1,imax dm=4.*pai*x**2.*ener*dx ! cm dp=-(ener+prener)*(ma+4.*pai*x**3.*prener)/(x*(x-2.*ma))*dx ! cm^(-2) c... new mass, pressure and radial position ma=ma+dm prener=prener+dp x=x+dx c... definition of quadrupole pressure perturbation array, and spherical and quadrupole stretching function arrays press2(i)=-functh2(i)-1./3.*x**2*w(i)**2*EXP(-2*nur(i)) ! dimentionless eps0(i)=-press(i)*(prener+ener)*dx/dp ! cm eps2(i)=-press2(i)*(prener+ener)*dx/dp ! cm c... determine index of new pressure within the pressure array do ipre=1,nf-1 if(((prener-pre(ipre))*(prener-pre(ipre+1))).lt.0.) goto 41 enddo 41 continue c... linear interpolation (to determine energy at the new position) ener=en(ipre)+(en(ipre+1)-en(ipre))/(pre(ipre+1)-pre(ipre))* & (prener-pre(ipre)) enddo 61 c... gives polar and equatorial radius, and eccentricity (description of the star’s shape at its surface) equatR=x+eps0(imax)-eps2(imax)/2. ! cm polarR=x+eps0(imax)+eps2(imax) ! cm eccent=(1-(polarR/equatR)**2)**0.5 ! dimentionless write(15,’(8g14.6,i5)’) x*conv2,ma/msolar,press2(imax), & eps0(imax),eps2(imax),polarR*conv2, & equatR*conv2,eccent,k write(18,’(9g14.6,i5)’) x*conv2,ma/msolar,(ma+massinc)/msolar, & equatR*conv2,polarR*conv2,kepler*conv5, & eccent,en1/conv1,inmom/conv4,k 12 continue enddo !k different central densities c... write maximum mass and the corresponding radius open(unit=22, file=’M_R_max_grav-pr.res’) print*, mmax, rmax write(22,*) ’Mmax Rmax’ write(22,’(2g14.6)’) mmax,rmax close(5) close(7) close(8) close(9) close(10) close(12) close(13) close(14) close(15) close(17) close(22) end 62