Applying chiral field theory to neutron star moments of inertia

In this work we examine the applications of neutron stars and their unique structure toprovide limits on the range of possible equations of state which describe nuclear matter.We follow the history of developments in nuclear physics and discuss the current state ofknowledge of neutron stars, including why these objects are well suited to the study ofnuclear physics. We then move on to examine recent work to apply a chiral effective-fieldmodel for the nuclear equation of state to a calculation of the moment of inertia object Ain the binary pulsar system J0737-3039. The methods used in this approach appear usefulfor ruling out different classes of equations of state, as well as potentially putting limits ondifferent structural properties of neutron stars that are difficult to measure directly, such asthe radius. We also examine other recent developments in the study of neutron stars andnuclear physics, such as the application of gravitational wave astronomy to observe neutronstar mergers.

APPLYING CHIRAL FIELD THEORY TO NEUTRON STAR MOMENTS OF INERTIA by Robert J. Stahulak A Senior Honors Thesis Submitted to the Faculty of The University of Utah In Partial Fulfillment of the Requirements for the Honors Degree in Bachelor of Science In The Department of Physics and Astronomy Approved: TTTTTTTTTTTTTTTTT Dr. Yue Zhao Supervisor TTTTTTTTTTTTTTTTT Dr. Peter Trapa Chair, Department of Physics and Astronomy TTTTTTTTTTTTTTTTT Tamara Young Department Honors Advisor TTTTTTTTTTTTTTTTT Dr. Sylvia D. Torti Dean, Honors College May 2019 ABSTRACT In this work we examine the applications of neutron stars and their unique structure to provide limits on the range of possible equations of state which describe nuclear matter. We follow the history of developments in nuclear physics and discuss the current state of knowledge of neutron stars, including why these objects are well suited to the study of nuclear physics. We then move on to examine recent work to apply a chiral effective-field model for the nuclear equation of state to a calculation of the moment of inertia object A in the binary pulsar system J0737-3039. The methods used in this approach appear useful for ruling out different classes of equations of state, as well as potentially putting limits on different structural properties of neutron stars that are difficult to measure directly, such as the radius. We also examine other recent developments in the study of neutron stars and nuclear physics, such as the application of gravitational wave astronomy to observe neutron star mergers. ii TABLE OF CONTENTS ABSTRACT ii 1 INTRODUCTION 1 1.1 UNDERLYING PHYSICS . . . . . . . . . . . . . . . . . . . 1 1.2 MOTIVATION . . . . . . . . . . . . . . . . . . . . . . . . . 6 2 METHODS 9 2.1 APPROACH . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 PERTURBATION THEORY . . . . . . . . . . . . . . . . . . 10 2.3 CHIRAL EFFECTIVE-FIELD THEORY . . . . . . . . . . . 11 2.4 BAYESIAN SAMPLING OF THE EQUATION OF STATE . . 12 2.5 TOLMAN-OPPENHEIMER-VOLKOV EQUATIONS . . . . 13 3 RESULTS 16 4 CONCLUSION 18 5 ACKNOWLEDGEMENTS 19 6 REFERENCES 20 A CODE FOR SOLVING THE TOV EQUATIONS 24 iii 1 1. INTRODUCTION 1.1. UNDERLYING PHYSICS 1.1.1. HISTORY OF NEUTRON STAR PHYSICS In the early 20th century, advances in telescope technology led to many new types of stellar objects and events being categorized. While supernova had been observed throughout history, observations of their remnants led to a new interest in understanding the underlying processes of these stellar explosions. While little was known at the time, it was speculated that observed supernovae were part of the natural life cycle of stars, and that they could cause the star to transition into a new form of stellar object, a neutron star [1]. Neutrons were known to be produced in normal stars over time, and if a star were to expel its outer ionized layers of gasses, the neutrons would be left behind, being unaffected by the intense electrical forces driving the explosion. With the advent of X-ray telescopes, a whole new spectrum of observation was opened up. With this came the first observations of pulsars and neutron stars in the nebulous remains of supernovae. With these discoveries came an interest in understanding the structure of the dim, massive objects. The development of general relativity also led to an interest in creating dynamical models of stellar objects using the new equations of motion. In particular, astronomers sought to model the extremely high-energy processes of supernova; in this process, it became clear that some must leave a dense core behind, providing a theory for the lineage of the neutron stars [8]. 1.1.2. DEGENERATE MATTER While the results of quantum mechanics are known most often to produce effects at very small scales, quantum effects also play a role in some of the largest and densest objects in the known universe. Ordinary matter is composed of nuclei, discrete neutrons and protons surrounded by dense electron clouds, bound together by electromagnetic interactions between these outer shells. However, the extremes of pressure and energy found throughout 2 the universe are capable of blurring our usual discretization of matter. Degenerate matter is a state of leptonic/hadronic matter formed at extremes of density and pressure. Fermions, particle with non-integer spin, obey the Pauli exclusion principle. Thus, any two fermions are incapable of occupying the same quantum state. In very dense matter, when large numbers of fermions are spatially confined, they begin to fill the majority of low energy states available. This prevents the matter from collapsing further, as many particles are forced to remain in states of high potential. This effect produces an apparent pressure, termed the degeneracy pressure, which effectively prevents more thermal energy from escaping the system [19]. Dense stellar bodies such as neutron stars, supernovae, and other potentially more exotic objects can form this bizarre state of matter, and the laws which govern their behavior are dominated by the effects of the degeneracy pressure. 1.1.3. NEUTRON STAR STRUCTURE Neutron stars are the densest form of stable matter known in the universe. The leftover core of a supermassive star gone nova, the nuclear material which composes a neutron star is denser than an atomic nucleus. Four teaspoons of neutron star-matter have as much mass as Earth's moon. Conserving the original angular momentum of their parent star even at their reduced size, neutron stars rotate at extraordinary rates for a macroscopic object, up to hundreds of times a second. These objects also tend to keep the magnetization of their progenitor, and as such produce extremely strong magnetic fields and beams of x-ray radiation from their poles. The incredible density of this peculiar state of matter causes these stars to exhibit properties characterized by extreme physical theories. The form and structure of neutron stars represent an intersection of astronomy, general relativity, nuclear physics, and quantum mechanics [7]. Neutron stars evolve from the cores of massive stars which undergo a core-collapse supernova. In this process, the fusion reactions which push the stellar matter outwards against gravity fail, and the core of the star compresses to a state of quantum degeneracy. As the 3 core is composed entirely of fermions, its collapse is eventually stopped by the degeneracy pressure which emerges due to the Pauli exclusion principle. At this point, the rigidity of the core causes the energy of the collapse to propagate back outward through the upper layers of the star as a shockwave. This blasts the mantle away into space, leaving only the core, which maintains a high concentration of heavy elements produced in the star. The high density of the core alters the equilibrium state of the matter inside. Neutrons leak from the nucleus and most of the remaining protons transform into neutrons through either direct or indirect electron capture processes. This produces significant numbers of neutrinos and antineutrinos, which though normally intangible, are partially trapped by the extreme density and size of the core. As the proto-neutron star shrinks, it becomes transparent to these particles and loses significant amounts of energy through neutrino emission, cooling far below its original temperature [8]. After this point, depending on the strength of the degeneracy pressure and the mass of the core, the collapse of the neutron star may continue, leading the core to become a black hole. Because angular momentum in conserved in this evolution, as the remaining core shrinks it attains its extremely high rate of rotation. In addition, since all stars maintain a dynamic magnetic field due to convective currents in the core, the matter of the neutron star contains a very large frozen-in magnetization. This combined with the high rate of rotation produces large magnetic fields on the order of 107 Tesla at the magnetic poles of the star [22]. These fields lead neutron stars to emit strongly in the x-ray spectrum from their magnetic poles; neutron stars spinning about an axis different than their magnetic poles thus produce a lighthouse effect. These objects are termed pulsars. A typical neutron star consists of an inner and outer core, an upper crust, and beyond this the envelope, and atmosphere. The inner and outer core together contain the vast majority of the stars mass and are mostly composed of a continuous state of matter at incredible density. The crust is composed of more conventional matter, including elements such as iron. Progressing towards the core of the star, the energy required to remove a neutron 4 from the nucleus drops to zero, and nucleons leak from their confinement. As the gravitational deformation increases nearer the core, nucleons are gradually spaghettified and even inverted entirely as the density increases beyond that of the nucleus. This eventually transitions into the continuous state of dense matter found at the core [9]. However, while the rate of rotation of neutron stars for the most part remains highly predictable clocks, occasional changes in the rotation of these stellar bodies occur. A standing problem in the understanding of neutron stars, these glitches amount to an apparent increase in the rotary speed of the neutron star. This implies a reservoir of angular momentum that is able to quite suddenly alter the motion of the entire body. Just below the crust of a neutron star, the free neutrons are compressed to a state of continuous matter, with no clear distinction between neighboring particles. Similar to the phenomenon of superconductivity, the intense pressures of the core force the nucleons to form Cooper pairs, becoming capable of occupying the same state. This state of matter is understood to behave as a perfect superfluid, capable of flowing without any resistance [17]. It is theorized that stable vortices in this superfluid mantle may be capable of quickly transferring angular momentum between the crust and the rest of the star. In addition, the interactions of the neutron vortices may create secondary vortices in the few protons within this layer; it is hypothesized that these neutron-proton vortices may also be responsible for maintaining the extreme magnetic fields produced by neutron stars [23]. However, enough uncertainties remain in existing equations of state for the neutron star that the glitches could be explained by additional angular momentum stored in a thicker than expected stellar crust [18]. The glitches in pulsar periods would then be due to occasional interactions between the more conventional crust and the neutron fluid of the mantle. Because the structure of neutron stars is determined by established physical theories, measurements of the properties of these objects can serve to test those theories explaining the behavior of the dense matter of which they are formed. By comparing the predictions of thsese theories to measurements made of known neutron stars these theories can be 5 improved. One focus for current improvements in nuclear physics is the equation of state of dense matter, which can be used to predict the radius and density of neutron stars. 1.1.4. THE NUCLEAR EQUATION OF STATE The matter in the core of the star is typically modeled using general relativity as an incompressible fluid. While this type of solution is nonphysical in general relativity, it is very accurate in the case of a neutron star due to the extreme pressure of the degenerate matter. These models provide a basis for the equation of state of a neutron star. In core-collapse supernova, neutron stars, and the merging of neutron stars, highly dense and degenerate states of matter are achieved. The incompressibility of this matter heavily drives the outcomes of these events, and thus in understanding them it is important to model dense matter states accurately. Combining existing knowledge of quantum chromodynamics and nuclear theory, early equations of state were developed. However, the small scale of time and space in these experiments made precise measurements difficult. With the discovery of neutron stars, it became clear that these stellar objects could serve as a more accurate way of improving the models. The equation of state is an important tool in modelling complex processes in astrophysics. These equations relate the internal pressure of the star to its energy density. The degeneracy pressure, strong nuclear interactions, and weak thermal pressure all come together to prevent a neutron star from collapsing past the realm of dense matter into a black hole. Because of this, the equation of state makes predictions bounding the size and mass of neutron stars that can exist. Below the minimum mass threshold, a star would be incapable of creating the pressures necessary to maintain the state of continuous matter which defines it. If the mass of a star core is too high, it will collapse into a black hole and cease to radiate. Thus, the range of observed neutrons stars and pulsars can serve indirectly as a way of verifying current models of dense matter and general relativity, by evaluating the equations of state that they produce for these bodies [10]. 6 1.2. MOTIVATION My research will be comprised of examining different theoretical equations of state for neutron stars and deriving their associated moments of inertia. Since the moment of inertia of the star depends precisely on its structure, it can also serve as an indirect measure of the equation of state which governs this. In the future surveys of the moments of inertia of different known neutron stars will be made. By comparing these to theoretical predictions of different equations of state, limits can be placed on the physical parameters of these models based on the stellar observations. This in turn can lead to improvements of the nuclear theories on which these models are based. Specifically, my research will seek to survey existing models of neutron stars based on different theoretical equations of state. By calculating a model for the moment of inertia of a neutron star based on each equation of state, future measurements of these moments of inertia will serve as verification for each theory. The general approach followed in our research is well-supported; other studies have shown strong dependencies of the moment of inertia on the parameterization of the nuclear equation of state [11, 20]. It is therefore important to perform similar calculations for equations of state based on a chiral effective-field model, in order bring in perspective based on microscopic nuclear physics. 1.2.1. MEASURING THE MOMENT OF INERTIA Because many pulsars exist in binary systems, there is potential to extract information about the stars structure from the orbital dynamics of the system. Recent studies have suggested that a precise measurement of the moment of inertia is possible for one such binary pulsar, specifically the system PSR J037-3039 [3]. The binary system is highly relativistic, and as such it is anticipated that the advance of the periastron of the object A could be observed to a high enough precision that second order post-Newtonian effects would be evident [2]. The corrections due to the advance of the periastron are a purely relativistic effect, and depend on the spin-orbit coupling of the object [2]. This coupling is determined in part 7 by the moment of inertia of the object. It is anticipated that from the known dynamics of the system, these effects could be measured to a level of accuracy such that the moment of inertia could be determined up to a confidence of 10% [11]. With the possibility of this future measurement in mind, we seek to use current knowledge of the equation of state of nuclear matter to calculate predictions of this moment of inertia, with the hope of more tightly constraining these equations of state. Because the moment of inertia is proportional to the radius squared of the object in question, for a given mass the moment of inertia is highly sensitive to the stiffness of the composite matter. This makes the moment of inertia an ideal tool for testing different equations of state in the regime of extreme densities. 1.2.2. ALTERNATIVE CONSTRAINTS ON THE EQUATION OF STATE One reason motivating our exploration of neutron structure using nuclear theory is that bounds for these properties already exist based on astronomical observations and theoretical limits based on general relativity. This allows for a potential feedback between the study of nuclear theory and these other disciplines. Another property of neutron stars which has shown a strong dependence on the nuclear equation of state is the love number, or tidal deformability [14]. The tidal deformability shows a very strong linear correlation to the pressure predicted at two times nuclear density, a value which diverges widely across the range of existing equations of state. This implies that measurements of this property would be another method to tightly constrain the parameters of the equation of state. Preliminary studies indicate that the love number could be reasonably ascertained from observation of the final inspiral of a neutron star merger event, as the tidal deformability largely determines the dynamics of this process [6, 4]. With improvements to the LIGO observatory forthcoming, events such as the recent multi-messenger observation of a binary neutron star merger have the potential to significantly contribute to nuclear theory through constraining these observables [12, 21]. Future 8 observations of these events will provide more opportunities, and further efforts would do well to compare results between moment of inertia and love number constraints. Models of the neutron star radius have also been shown to have a strong dependence on the equation of state used. This would allow precise measurements of the radius of a neutron star of known mass to also tightly constrain the nuclear equation of state. Current methods which may produce these types of measurements include observing the pulsation and spectra of x-ray bursts on the neutron star surface [24]. The mass and radius of the star is encoded into the behavior of these hot spots through relativistic effects. 9 2. METHODS 2.1. APPROACH With the goal of predicting a range of moments of inertia for the object PSR J037-3039A, such that a precise measurement will exclude a large band of EOS, we begin by building a large set containing the structure of many neutron stars of varying masses. Using the TOV equations we will solve numerically to construct the neutron star structure. This set of coupled differential equations has a long history for describing neutron stars [16]. A set of equations of state is sampled from the posterior Bayesian distribution and used to construct the set of possible neutron stars; each equation of state consists of a discrete set of corresponding pressure and energy densities. In order to calculate energy densities at intermediate pressures, a form of polynomial interpolation is used, fitting the EOS as a whole and interpolating to the required values. We begin with a range of central energy densities for the neutron star. This is range is constructed so as to create stars which respect current limits on the maximum and minimum masses. By applying our knowledge of the nuclear equation of state, we can build the mass and energy density of the neutron star in differential layers from the center, using the density from the previous layer to calculate the pressure of the next. This process is repeated until the density reaches an anticipated value associated with the crust of the star, which is composed of more conventional matter. The actual numerical integration is done using a Runge-Kutta method, using the interpolated data and values from the previous steps. The Runge-Kutta method is a weighted mid-point method of numerical integration; in this work, we apply the fourth-order method, which includes four terms in the numerical integral. We perform this integration until the interpolated energy density reaches a regime associated with the more conventional matter of the neutron star crust. The bulk matter EOS is then stitched together with a fixed low density EOS in order to achieve realistic results. 10 We then continue the integration using this fixed EOS until the pressure reaches zero, representing the conventional edge of the star. While this procedure ignores the contribution of the atmosphere, this is considered to be negligible. By calculating the contribution to the moment of inertia at each differential layer, we can also calculate a total value for the moment of inertia, alongside the total mass and radius of the star. Though difficult to solve for directly, it is possible to obtain a simple differential equation for an intermediate variable φ , which is related to the moment of inertia through the simple relation shown in Equation 7. Thus we have a mass, radius, and moment of inertia value corresponding to a wide range of neutron stars, each with a unique central density and associated with a particular variation of the nuclear equation of state. In order to make predictions for the system PSR J0737-3039, we interpolate the data to the established mass of the object A, 1.338 times the mass of the sun [3]. We now examine the distribution of neutron stars with this mass, and analyze correlations between our observables. 2.2. PERTURBATION THEORY The nuclear equation of state is derived by applying perturbation theory to a noninteracting system of nucleons. Different interactions associated with the strong force are taken as different terms in the perturbative expansion for constructing the ground state Hamiltonian, each of which is constructed to respect the approximate chiral symmetries of QCD. Interaction diagrams for each of the terms in this expansion are shown in Figure 1. The leading order terms contain the simplest interactions, including only single pion exchanges between two nucleons and an effective contact interaction. The next to leading order term includes many two pion interactions. The next-to-next-to leading order term begins including three-nucleon interactions, and this is the term at which the expansion is stopped for our calculations. The N3 LO interactions can then be used to estimate the magnitude of the error in the expansion. 11 2.3. CHIRAL EFFECTIVE-FIELD THEORY Chiral effective-field theory (Chi-EFT) is an approach to constructing the Hamiltonian in low-energy nuclear interactions. This is done by creating an effective theory which respects the approximate chiral, parity, and charge conjugation symmetries of Quantum Chromodynamics. Because the goal is low-energy physics Chi-EFT does not consider quarks or gluons as degrees of freedom. In the Chi-EFT employed here, only pion exchanges between nucleons are treated precisely according to QCD [13]. The whole range of interactions considered are shown in Figure 1. Other interactions of higher energy are contained in contact terms of varying strength and proximity (crosses and dots in the figure) in the Hamiltonian, which act only over a very short range consistent with the typical energy scales of bulk nuclear matter. Given this approximation, it is possible that the resultant equations of state diverge significantly from reality in the case of extremely dense or hot matter, in which the higher energy strong interactions become more important. This is another reason motivating our predictions for neutron stars, as our calculations for neutron star masses and moments of inertia can be checked against limits and measured values established by the astronomical community. Finally, the equation of state is calculated from the ground state energy as the energy per particle, expressed as a function of the density and temperature. 12 Figure 1: Orders of expansion and associated interactions in perturbative Chi-EFT [15] 2.4. BAYESIAN SAMPLING OF THE EQUATION OF STATE The equations of state obtained purely from the chiral effective-field theory are taken as the prior distribution in applying Bayes' theorem. Bayes' theorem relates the probability of a hypothesis being true given a set of observations (the posterior likelihood function), and the probability of the observations occurring assuming the hypothesis to be true (the prior likelihood function). By taking each equation of state as a hypothesis, attempting to explain the data given by low-energy nuclear scattering experiments, Bayes' theorem can be used to increase the accuracy of the model. Formally, Bayes' theorem states P(H|E) = P(E|H) ∗ P(H) P(E) (Bayes' Theorem) Where P(H|E) is the probability that a hypothesis is true given the data (the posterior probability), and P(E|H) is the probability that the data is reproduced by the hypothesis (the prior distribution). Using Bayes' theorem, the set of equations which form the prior distribution are then modified based on the posterior probability that each individual equation 13 reproduces the data. These likelihood functions can be readily calculated using Bayes' theorem. The result of applying the weights to the ab initio calculations are shown in Figure 2. Figure 2: Prior and Posterior distributions of the Chi-EFT equation of state [14] 2.5. TOLMAN-OPPENHEIMER-VOLKOV EQUATIONS In order to construct the structure of the neutron star, we seek to solve the Tolman-OppenheimerVolkov (TOV) equations. This differential equation represents the hydrostatic equilibrium for a spherically symmetric body in general relativity, determining the radial pressure differential needed at every point on the interior in order to support the mass above. The general TOV equation is typically expressed dP Gmρ P 4πr3 P 2Gm −1 1+ 1− 2 =− 1+ 2 dr r ρc mc2 rc (Tolman-Oppenheimer-Volkov Equation) 14 The general TOV equation can be separated into differential equations for the mass and pressure. m + 4πr3 P dP = −(ρ + P) 2 dr r − 2mr dm = 4πr2 ρ dr (Separated TOV Equations) With an equation of state relating the energy density and pressure of the neutron matter, the separated TOV equations become readily solvable using numerical methods for linear differential equations. We approached the problem of solving these equations using a Runge-Kutta method of numerical integration. In order to calculate the moment of inertia, it is necessary to rearrange the TOV equations to find a differential relation for the moment of inertia in terms of the neutron star radius. In a slow rotating and spherically symmetric relativistic star, the moment of inertia is found to be I= 8π 3 Z R r4 (ρ + P)e λ −ν 2 (Moment of Inertia, Slow-Rotating [5]) ωdr 0 Where λ = −ln(1 − 2m/r), ν is the metric coefficent, and ω is the rotational drag function. Setting j = e−(λ +ν)/2 , ω must satisfy dj d 4 dω r j = −4r3 ω dr dr dr (1) with the boundary conditions 2I ωR = 1 − 3 R dω dr = 0. (2) R4 dω = 6 dr R (3) 0 Using Equation 1, the moment of inertia becomes 2 I=− 3 Z R 0 dj 1 r ω dr = dr 6 3 Z R 0 dj d r ω dr 4 15 Setting φ = dlnω dlnr , we can rewrite Equation 1 as φ dln j dφ = − (φ + 3) − (4 + φ ) dr r dr (4) 4πr2 dln j =− (ρ + p), dr r − 2m (5) where with the added boundary condition that φ (0) = 0. Using the boundary conditions in Equations 2, our definition of φ , and substituting into Equation 3, we find R3 φR I = φR ωR = (R3 − 2I) 6 6 (6) R3 φR I= 6 + 2φR (7) rearranging, we have Equation 4 is the differential equation which we solve numerically in the code found in appendix A. This allows for a simple substitution at the end of the calculation in order to find the moment of inertia in Equation 7. The main program is written in Fortran 90, and performs the Runge-Kutta integration by using an Aitken interpolation method to find the energy density at each step in the pressure, and calculating values for the derivatives in the equations above. The calculation of the pressure P, energy density ρ, and φ is repeated for each equation of state in a large set sampled from the prior distribution, and for a range of central densities within each equation of state. The resulting set of neutron star observables is discussed below. 16 3. RESULTS Using this dataset, we can construct a probability density for different combinations of values of observables of the neutron star, such as total mass, mass of the crust, radius, and moment of inertia. These histograms demonstrate the range of values for these observables predicted by the current nuclear theory, and the most likely values based on those predictions. Two of these histograms are shown in Figures 3, and 4. It is evident from Figures 4 that there is a strong, linear correlation in our dataset between the neutron star radius and moment of inertia. It is evident from the data that a measurement of the moment of inertia PSR J0737-3039A within 10% error would put bounds on the radius of the object with an error on the order of half a kilometer. With the mass of the object already established, this range could be compared with existing limits such as estimates based on x-ray observations. Figure 3 shows a very tight and mostly linear relationship between the calculated neutron star moment of inertia and the pressure at two times saturation density predicted by the equation of state. This pressure is highly characteristic of each particular equation of state, and differs wildly across the range of equations produced by the Bayesian modelling. This behavior is due to the fact that very little data exists which probes the behavior of nuclei in this region of extreme densities. Given this relationship, it is safe to conclude that a measurement of the neutron star moment of inertia within 10% would be sufficient to probe the behavior of nuclear matter at very high densities, and so tightly constrain the equation of state. 17 1.6 100 I (g cm2 ∗ 1045) 1.4 1.2 10−1 1.0 0.8 10−2 5 10 15 20 P2ρ0 (MeV fm−3) 25 30 Figure 3: Moment of Inertia v. Pressure at two times saturation density (M = 1.338 ∗ Msun ) 1.6 100 1.5 I (g cm2 ∗ 1045) 1.4 1.3 1.2 10−1 1.1 1.0 0.9 0.8 10−2 8 9 10 11 12 13 14 R (km) Figure 4: Moment of Inertia v. Neutron Star Radius (M = 1.338 ∗ Msun ) 18 4. CONCLUSION We began by constructing a distribution of nuclear equations of state based on a chiral effective-field theory. This was done by calculating the ground state energy of a system of nucleons through an effective interaction nuclear potential built on interactions which respect the approximate chiral symmetries of quantum chromodynamics. Translating the ground state energy into an equation of state based on density and temperature, a prior distribution was made by considering the higher order terms of the perturbative expansion. This distribution was further refined by applying Bayesian statistics to construct a posterior distribution based on the equations' ability to reproduce nuclear scattering data. Sampling from this distribution, we calculated ranges for the mass distribution, total mass, radius, and moment of inertia of a neutron star by solving the Tolman-OppenheimerVolkov equations numerically. By analyzing this set of neutron stars, we obtain probability distributions which show tight correlations between various neutron star observables, most importantly between the moment of inertia and the pressure at two times saturation density. While satisfactory measurements are still forthcoming, it is clear that the neutron star moment of inertia is a powerful tool for probing the behavior of nuclear matter at high densities. Though this method is not new, it has yet to be used in a way which accurately utilizes the QCD physics which underlies the behavior of nucleons. 19 5. ACKNOWLEDGEMENTS I would like to thank my advisors at the Texas A&M Cyclotron institute, Drs. Jeremy Holt and Yeunhwan Lim. Dr. Holt's patience and quiet knowledge allowed me to grow into this project, and Dr. Lim's constant willingness to help got me through to the end. Without the both of them, I would never have been able to complete this work. I would also like to thank my Honors Thesis advisor, Dr. Yue Zhao, for being willing to go out on a limb and help me through this process. Thank you to the NSF REU grants PHY 1659847 and PHY 1652199 which supported my research and made this work possible. 20 6. REFERENCES [1] W. Baade and F. Zwicky, Remarks on super-novae and cosmic rays [6], Physical Review 46 (1934), no. 1, 76-77. [2] Joaquim A Batlle and Rosario López, Revisiting the border between Newtonian mechanics and General Relativity : The periastron advance, 10 (2014), 65-72. [3] M. Bejger, T. Bulik, and P. Haensel, Constraints on the dense matter equation of state from the measurements of PSR J0737-3039A moment of inertia and PSR J0751+1807 mass, Monthly Notices of the Royal Astronomical Society 364 (2005), no. 2, 635- 639. [4] D. Gondek-Rosińska, M Bejger, T Bulik, E Gourgoulhon, Paweł Haensel, F Limousin, K Taniguchi, and L Zdunik, The final phase of inspiral of neutron stars: Realistic equations of state, Advances in Space Research 39 (2007), no. 2, 271-274. [5] James B Hartle, Slowly Rotating Relativistic Stars. I. Equations of Structure, Astrophysical Journal 150 (1967), no. December, p.1005. [6] Benjamin D Lackey, Koutarou Kyutoku, Masaru Shibata, Patrick R Brady, and John L Friedman, Extracting Equation of State Parameters from Black Hole-Neutron Star Mergers, Proceedings of the Wisconsin Space Conference (2015), 1-21. [7] J. M. Lattimer and M. Prakash, The Physics of Neutron Stars, Science 304 (2004), no. April, 1-7. [8] James M. Lattimer, Introduction to neutron stars, EXOTIC NUCLEI AND NUCLEAR/PARTICLE ASTROPHYSICS (V). FROM NUCLEI TO STARS: Carpathian Summer School of Physics 2014 1645 (2015), no. May, 61-78. [9] , Neutron Star Physics and EOS, EPJ Web of Conferences 109 (2016), 07001. 21 [10] James M. Lattimer and Madappa Prakash, The equation of state of hot, dense matter and neutron stars, Physics Reports 621 (2016), 127-164. [11] James M. Lattimer and Bernard F. Schutz, Constraining the Equation of State with Moment of Inertia Measurements, The Astrophysical Journal 629 (2005), no. 2, 979- 984. [12] LIGO Scientific Collaboration, Virgo Collaboration, Fermi GBM, INTEGRAL, Icecube Collaboration, AstroSat Cadmium Zinc Telluride Imager Team, IPN Collaboration, The Insight-Hxmt Collaboration, ANTARES Collaboration, The Swift Collaboration, AGILE Team, The 1M2H Team, The Dark Energy Camera GW-EM Collaboration, the DES Collaboration, The DLT40 Collaboration, GRAWITA, :, GRAvitational Wave Inaf TeAm, The Fermi Large Area Telescope Collaboration, ATCA, :, Australia Telescope Compact Array, ASKAP, :, Australian SKA Pathfinder, Las Cumbres Observatory Group, OzGrav, DWF, AST3, CAASTRO Collaborations, The VINROUGE Collaboration, MASTER Collaboration, J-GEM, GROWTH, JAGWAR, Caltech NRAO, TTU-NRAO, NuSTAR Collaborations, Pan-STARRS, The MAXI Team, TZAC Consortium, KU Collaboration, Nordic Optical Telescope, EPESSTO, GROND, Texas Tech University, SALT Group, TOROS, :, Transient Robotic Observatory of the South Collaboration, The BOOTES Collaboration, MWA, :, Murchison Widefield Array, The CALET Collaboration, IKI-GW Follow-up Collaboration, H. E. S. S. Collaboration, LOFAR Collaboration, LWA, :, Long Wavelength Array, HAWC Collaboration, The Pierre Auger Collaboration, ALMA Collaboration, Euro VLBI Team, Pi of the Sky Collaboration, The Chandra Team at McGill University, DFN, :, Desert Fireball Network, ATLAS, High Time Resolution Universe Survey, RIMAS, RATIR, and SKA South Africa/MeerKAT, Multi-messenger Observations of a Binary Neutron Star Merger, 12 (2017). [13] Yeunhwan Lim and Jeremy W. Holt, Structure of neutron star crusts from new Skyrme 22 effective interactions constrained by chiral effective field theory, Physical Review C 95 (2017), no. 6, 1-11. [14] , Neutron star tidal deformabilities constrained by chiral effective field theory, Phys. Rev. Lett. 121 (2018), 1-6. [15] R. Machleidt and D. R. Entem, Chiral effective field theory and nuclear forces, Physics Reports 503 (2011), no. 1, 1-75. [16] J. R. Oppenheimer and G. M. Volkoff, On massive neutron cores, Physical Review 55 (1939), no. 4, 374-381. [17] Dany Page, Madappa Prakash, James M. Lattimer, and Andrew W. Steiner, Rapid cooling of the neutron star in cassiopeia a triggered by neutron superfluidity in dense matter, Physical Review Letters 106 (2011), no. 8, 1-4. [18] J. Piekarewicz, F. J. Fattoyev, and C. J. Horowitz, Pulsar glitches: The crust may be enough, Physical Review C - Nuclear Physics 90 (2014), no. 1, 1-11. [19] A. Y. Potekhin, The Physics of Neutron Stars, Compute (2010), 1-7. [20] Carolyn A. Raithel, Feryal Özel, and Dimitrios Psaltis, Model-independent inference of neutron star radii from moment of inertia measurements, Physical Review C 93 (2016), no. 3, 1-5. [21] Jocelyn S. Read, Charalampos Markakis, Masaru Shibata, Koji Uryu, Jolien D E Creighton, and John L. Friedman, Measuring the neutron star equation of state with gravitational wave observations, Physical Review D - Particles, Fields, Gravitation and Cosmology (2009). [22] A. Reisenegger, Origin and evolution of neutron star magnetic fields, International Workshop on Strong Magnetic Fields and Neutron Stars (2002), 33-49. 23 [23] K. M. Shahabasyan and M. K. Shahabasyan, Vortex structure of neutron stars with triplet neutron superfluidity, Astrophysics 54 (2011), no. 3, 429-438. [24] Anna L. Watts, Nils Andersson, Deepto Chakrabarty, Marco Feroci, Kai Hebeler, Gianluca Israel, Frederick K. Lamb, M. Coleman Miller, Sharon Morsink, Feryal Özel, Alessandro Patruno, Juri Poutanen, Dimitrios Psaltis, Achim Schwenk, Andrew W. Steiner, Luigi Stella, Laura Tolos, and Michiel Van Der Klis, Colloquium: Measuring the neutron star equation of state using x-ray timing, Reviews of Modern Physics (2016). 24 A. CODE FOR SOLVING THE TOV EQUATIONS module tov CONTAINS SUBROUTINE RK42(pre, dp, rho, y, mass, edata, pdata, ndata, nmax) !! Runge-Kutta 4th order method with beta stable nuclear matter IMPLICIT NONE REAL(8) :: pre, dp, u, y, pre1, pre2 INTEGER :: NMAX REAL(8), DIMENSION(NMAX) :: edata, pdata, REAL(8) :: energy, mun, mup, mue, eb, ed, REAL(8) :: dydp, dmdp, mass, df INTEGER :: I, J, K, iter REAL(8) :: mass1, dm1, mass2, dm2, mass3, REAL(8) :: dy, y1, dy1, y2, dy2, y3, dy3, real(8) :: rho1, rho2, rho3, rho real(8) :: dedp, dedp1, dedp2, dedp3 pre1 = pre + 0.5*dp pre2 = pre + dp ndata ed1, ed2, ed3, ed4 dm3, mass4, dm4 y4, dy4 !!---------------------------------------------------!!Begin the first step of the Runge-Kutta integration call AITKEN_piece(NMAX, pdata, edata, pre, ed, dedp, DF) CALL deriv(ed, pre, mass, y, dydp, dmdp) dy1 = dp*dydp dm1 = dp*dmdp mass1 = mass + 0.5*dm1 y1 = y + 0.5*dy1 !!---------------------------------------------------!!Begin the second step !! interpolate the energy density call AITKEN_piece(NMAX, pdata, edata, pre1, ed1, dedp1, DF) CALL deriv(ed1, pre1, mass1, y1, dydp, dmdp) dy2 = dp*dydp 25 dm2 = dp*dmdp mass2 = mass + 0.5*dm2 y2 = y + 0.5*dy2 !!---------------------------------------------------!!Begin the third step CALL deriv(ed1, pre1, mass2, y2, dydp, dmdp) dy3 = dp*dydp dm3 = dp*dmdp mass3 = mass + dm3 y3 = y + dy3 !!---------------------------------------------------!!Begin the fourth step !!interpolate the energy density call AITKEN_piece(NMAX, pdata, edata, pre2, ed3, dedp3, DF) CALL deriv(ed3, pre2, mass3, y3, dydp, dmdp) dy4 = dp*dydp dm4 = dp*dmdp pre = pre2 y = y + (dy1 + 2*dy2 + 2*dy3 + dy4)/6. mass = mass + (dm1 + 2*dm2 + 2*dm3 + dm4)/6. CALL AITKEN(NMAX, pdata, ndata, pre2, rho3, df) !! interpolation of number density rho = rho3 Return End SUBROUTINE RK42 SUBROUTINE AITKEN(NNMAX, XXI, FFI, X, F, DF) ! ! Subroutine performing the Lagrange interpolation with the ! Aitken method. F: interpolated value. DF: error estimated. ! Copyright (c) Tao Pang 1997. ! use parameters, only: iprint IMPLICIT NONE 26 INTEGER, PARAMETER :: NMAX = 21, N = 5 INTEGER :: I, J, kmin, NNMAX, i1, i2 REAL(8), INTENT(IN) :: X REAL(8), INTENT(OUT) :: F, DF REAL(8) :: X1, X2, F1, F2 REAL(8), INTENT(IN), DIMENSION(NNMAX):: XXI, FFI REAL(8), DIMENSION(N) :: XI, FI REAL(8), DIMENSION(NMAX):: FT REAL(8) :: Mini, minn if (x .eq. xxi(1)) then f = ffi(1) df = 0.0 return endif if (x .eq. xxi(nnmax)) then f = ffi(nnmax) df = 0.0 return endif goto 1055 LOOP_DEC: DO I = 1, NNMAX - 1 IF ((XXI(I) - X)*(XXI(I + 1) - X) .LT. 0) THEN kmin = i exit loop_dec ENDIF ENDDO LOOP_DEC 1055 continue kmin = 1 call locate(xxi, nnmax, x, kmin) i1 = kmin if (kmin .lt. 3) then kmin = 3 endif if (kmin .gt. nnmax - 3) then kmin = nnmax - 3 27 endif DO I = 1, 5 XI(I) = XXI(KMIN - 3 + I) FI(I) = FFI(KMIN - 3 + I) ENDDO DO I = 1, N FT(I) = FI(I) END DO ! DO I = 1, N - 1 DO J = 1, N - I X1 = XI(J) X2 = XI(J + I) F1 = FT(J) F2 = FT(J + 1) FT(J) = (X - X1)/(X2 - X1)*F2 + (X - X2)/(X1 - X2)*F1 END DO END DO F = FT(1) !! !! !! !! Sometimes f does not in the range of f(1) and f(5) below function prevents it if (f .gt. fi(1) .and. f .gt. fi(5) .or. & f .lt. fi(1) .and. f .lt. fi(5)) then f = fi(3) + (fi(4) - fi(3))/(xi(4) - xi(3))*(x - xi(3)) endif if (f .gt. ffi(i1) .and. f .gt. ffi(i1 + 1)) then if (iprint .eq. 1) then print *, 'outside the two points I' endif i2 = i1 + 1 f = ffi(i1) & + (ffi(i2) - ffi(i1))/(xxi(i2) - xxi(i1))*(x - xxi(i1)) endif if (f .lt. ffi(i1) .and. f .lt. ffi(i1 + 1)) then if (iprint .eq. 1) then print *, 'outside the two points II' 28 endif i2 = i1 + 1 f = ffi(i1) & + (ffi(i2) - ffi(i1))/(xxi(i2) - xxi(i1))*(x - xxi(i1)) endif DF = (ABS(F - F1) + ABS(F - F2))/2.0 END SUBROUTINE AITKEN SUBROUTINE locate(xx, n, x, j) implicit none INTEGER :: j, n REAL(8) :: x, xx(n) !! !! !! !! !! Given an array xx(1:n), and given a value x, returns a value j such that x is between xx(j) and xx(j+1). xx(1:n) must be monotonic, either increasing or decreasing. j=0 or j=n is returned to indicate that x is out of range. INTEGER :: jl, jm, ju jl = 0 !!Initialize lower ju=n+1 and upper limits. ju = n + 1 10 continue if (ju - jl .gt. 1) then jm = (ju + jl)/2 !! If we are not yet done, compute a midpoint, if ((xx(n) .ge. xx(1)) .eqv. (x .ge. xx(jm))) then jl = jm else ju = jm endif goto 10 endif if (x .eq. xx(1)) then j = 1 elseif (x .eq. xx(n)) then j = n - 1 else j = jl endif 29 return END subroutine !! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !! !! Start of the moment of inertia calculation !! !! subroutine tovphi_p_gen(rhocrust, pcrust, infile, outfile, massmax) implicit none real(8) :: rhocrust character(len=80) :: infile, outfile, tempfile real(8), allocatable, dimension(:) :: pdata, ndata, edata real(8), real(8), real(8), real(8), dimension(1000) dimension(1000) dimension(1000) dimension(1000) :: :: :: :: ncdata, mdata, radata, crdata, radcrdata mdata1, rdata1, indata, incdata k2data, lbdata, bcdata, lsdata k2data1, lbdata1, bcdata1, lsdata1 real(8), dimension(1000) :: dim_icr1, dim_icr2 real(8), dimension(20000) :: dim_m, dim_r, dim_y, dim_n real(8), dimension(20000) :: dim_phi, dim_p integer :: nmax, eof real(8) :: a1, b1, c1, d1, e1, f1, h1, l1 real(8) :: pre, pre0, rhoc, rho, mst, y, y0, rado, rhoold real(8) :: phi, inertia real(8) real(8) real(8) real(8) real(8) :: :: :: :: :: pcrust radn, rad_crust, radius, radmax, cru_thick massmax, rhocmax, mass, df, dp, gue, crumin rho0, rhoci, z, z0, lp0 premin real(8) real(8) real(8) real(8) integer integer real(8) real(8) :: :: :: :: :: :: :: :: lovek2, bc, yr, lambt, lambs lovek2a, bc1, yr1, lambt1, lambs1 mass_crust, mst1, phi_crust, mass_core, phi_core inertia_core, inertia_crust, z1 i, j, k, imax(1), imax1, jmax, iinteg istr, nmax1, istr1, icnt dedpc !! dedp check iovermr2 30 integer, parameter :: id = 100 real(8) :: dphi(0:id), drad(0:id), dr, hd(0:id), rtemp real(8) :: xrc, xr, sum1, phitemp, omegr, omegrc real(8) :: inert, inertc real(8), parameter :: pi = 3.1415925635 real(8) :: ede, sum2 real(8), dimension(20000) :: dim_lnj real(8) :: dlnj(0:id), lnj, jrc real(8) :: radc, pt, et, ecrust, icr1, icr2, itot real(8), parameter :: ap = 1.3234d-6 !! 1MeV/fm^3 = 1.3234E-6 /km^2 !! for simpson rule hd = 1.0 hd(0) = 17./48.0 hd(1) = 59./48.0 hd(2) = 43./48.0 hd(3) = 49./48.0 hd(id) = hd(0) hd(id - 1) = hd(1) hd(id - 2) = hd(2) hd(id - 3) = hd(3) rho0 = 0.16 !! Set such that we reproduce established mass limits nmax = 0 !! don't forget initial set up for nmax=0 OPEN (unit=41, file=infile, status='old', recl=80, action='read', & IOSTAT=eof) DO WHILE (EOF >= 0) nmax = nmax + 1 READ (41, *, IOSTAT=EOF) a1, b1, c1 ENDDO nmax = nmax - 1 !!Allocate arrays to contain energy density (e), pressure (p), !!and number density (n) data ALLOCATE (edata(nmax), pdata(nmax), ndata(nmax)) REWIND (41) !!Read in the set of points describing the current EOS DO I = 1, NMAX READ (41, *) edata(i), pdata(i), ndata(i) 31 ENDDO CLOSE (41) imax = maxloc(edata) istr1 = imax(1) istr = max(istr1, 3) do i = 2, nmax/2 !! note that dedpc = 1/cs^2 dedpc = (edata(i + 1) - edata(i - 1))/(pdata(i + 1) - pdata(i - 1)) dedpc = abs(dedpc) if (dedpc .lt. 0.999) istr = i enddo istr = max(istr, 3) !! Procedure to check dedp = 1/cs^2 is reasonable nmax1 = nmax - istr + 1 rhoci = ndata(istr + 1) massmax = 0.0 loop_rho: do k = 0, 4000 rhoc = rhoci - 0.005*real(k) if (rhoc .lt. 0.08 .or. rhoc .lt. rhocrust - 0.01) then exit loop_rho endif gue = rhoc !! Interpolate the central pressure value !! from central density, EOS data CALL AITKEN(NMAX1, ndata(istr:nmax), pdata(istr:nmax), gue, pre, df) !! note z=r^2 pre0 = pre 32 rho mst y = z = phi = rhoc = 0.0 2.0 !! for the Love number 0.0 !! for the radius z= r^2 = 0.0 iinteg = 0 !!Perform step-wise numerical integration to calculate the mass, !!radius, moment of inertia for the current EOS, central density !!We step forward in small pressure intervals, and calculate the !!requisite radius, mass, and phi contribution from the pressure layer LOOP_MASS: DO J = 1, 20000 iinteg = iinteg + 1 z0 = z rado = z0**0.5 rhoold = rho !!Determine whether we are in the bulk matter or in the crust; if (pre .gt. 0.05*pre0) then dp = -0.01*abs(pre0) else dp = -0.1*abs(pre) endif 120 continue !! Runge-Kutta 4th order method with beta stable nuclear matter !! Do numerical integration, updating the values of rho, z (r^2), !! mst(mass), and phi(for the moment of inertia calculation) call RK42phi(pre, dp, rho, ede, z, mst, phi, & edata(istr:nmax), pdata(istr:nmax), & ndata(istr:nmax), nmax1) IF (RHO .LT. 0) THEN WRITE (*, *) 'DOES NOT CONVERGE AT rho', rho ENDIF radn = z**0.5 dim_m(j) = mst dim_r(j) = radn dim_y(j) = y 33 dim_n(j) = rho dim_phi(j) = phi dim_p(j) = pre !! be careful with the sig dim_lnj(j) = -4.*pi*z/(radn - 2.*mst)*(ede + pre)*ap !! !! !! IF We check if the integration loop as reached the edge of the star; if this is true, we have reached vacuum and the calculation is complete (rho .lt. 1.0e-12) then exit LOOP_MASS ENDIF ENDDO LOOP_MASS !! for the calculation of crust contribution CALL AITKEN(iinteg, dim_p(1:iinteg), dim_r(1:iinteg), & pcrust, rad_crust, df) !! interpolation of radius !! for the calculation of crust contribution CALL AITKEN(iinteg, dim_p(1:iinteg), dim_phi(1:iinteg), & pcrust, phi_core, df) !! interpolation of radius !! for the calculation of crust contribution CALL AITKEN(iinteg, dim_p(1:iinteg), dim_m(1:iinteg), & pcrust, mst1, df) !! interpolation of radius dr = (radn - rad_crust)/real(id) do i = 0, id - 1 rtemp = rad_crust + dr*real(i) drad(i) = 1./rtemp CALL AITKEN(iinteg, dim_r(1:iinteg), dim_phi(1:iinteg), & rtemp, phitemp, df) !! interpolation of phi dphi(i) = phitemp call AITKEN(iinteg, dim_r(1:iinteg), dim_lnj(1:iinteg), & rtemp, lnj, df) !! interpolation of log(j) dlnj(i) = lnj enddo drad(id) = 1./radn dphi(id) = phi 34 sum1 = dr*sum(drad*hd*dphi) sum2 = dr*sum(dlnj*hd) jrc = exp(-sum2) mass_core = mst1/1.4766 !! mass of core z1 = rad_crust**2 !! the radius of the core-crust transition !! For full NS configuration ! ! mass = mst/1.4766 !! 1Msun = 1.4766km radius = z**(0.5) inertia = phi*z**(1.5)/(6.+2.*phi)/1.4766 inertia = phi/(6. + 2.*phi)*z**0.5/mst itot = phi*z**(1.5)/(6. + 2.*phi) inert = phi*z**(1.5)/(6.+2.*phi) !! in unit of M_sun km^2 !! in unit of I/M R^2 !! in unit of km^3 !! km^3 unit omegr = 1.-2*inert/z**(1.5) omegrc = omegr*exp(-sum1) inertc = rad_crust**3/6.*omegrc*phi_core*jrc inertia_core = inertc/1.4766 !! in unit of M_sun km^2 !! crust contribution to moment of inertia inertia_crust = inertia - inertia_core itot = phi/(6.+2.*phi)*z**0.5/mst !! in unit of I/M R^2 radc = rad_crust pt = pcrust CALL AITKEN(NMAX1, pdata(istr:nmax), edata(istr:nmax), & pcrust, ecrust, df) !! interpolation of pressure et = ecrust !!Calculate the moment of inertia contribution from the crust call analytic_crust_moi(mst, radn, radc, pt, et, & itot, icr1, icr2) cru_thick = radius - rad_crust !!Track the maximum mass neutron star !!produced by this equation of state 35 if (mass .gt. massmax) then massmax = mass rhocmax = rhoc radmax = radius crumin = cru_thick endif icnt = icnt + 1 ncdata(icnt) = rhoc mdata(icnt) = mass radata(icnt) = radius crdata(icnt) = cru_thick indata(icnt) = inertia incdata(icnt) = inertia_crust dim_icr1(icnt) = icr1 dim_icr2(icnt) = icr2 radcrdata(icnt) = rad_crust enddo loop_rho nmax1 = icnt deallocate (edata, pdata, ndata) imax = maxloc(mdata(1:nmax1)) imax1 = imax(1) open (44, file=outfile, status='unknown') do i = imax1, nmax1 mass = mdata(i) radius = radata(i) inertia = indata(i) iovermr2 = inertia/(mass*radius**2) inertia_crust = incdata(i) icr1 = dim_icr1(i) icr2 = dim_icr2(i) write (44, '(14es14.5)') ncdata(i), mdata(i), radata(i), & inertia, mass/radius, iovermr2, & inertia_crust/inertia, & inertia_crust, &!!1.4766 , icr1,icr2, & 36 radcrdata(i) enddo close (44) return end subroutine tovphi_p_gen subroutine analytic_crust_moi(mst, radn, radt, pt, et, itot, icr1, icr2) !! note that this is km unit !! radn : radius of neutron star !! radt : radius of core !! pt : pressure at the transition !! et : energy density at the crust !! itot : total moment of inertia (km^3 unit) !! icr1 : 1st formula !! icr2 : 2nd formula implicit none real(8) :: mst, radn, radt, pt, et, itot, icr1, icr2 real(8) :: rads real(8), parameter :: ap = 1.3234d-6 !! 1MeV/fm^3 = 1.3234E-6 /km^2 real(8), parameter :: pi = 3.14159265358979 rads = 2.*mst icr1 = 16./3.*pi*(radt)**6/rads & *ap*pt*(1.-(rads/radn)*itot) & *(1.+48./5.*(radt/rads - 1)*(pt/et)) !! eq (14) !! km unit icr2 = 16./3.*pi*(radt)**6/rads & *ap*pt*(1.-0.21/(radn/rads - 1.0)) & *(1.+48./5.*(radt/rads - 1)*(pt/et)) return end subroutine SUBROUTINE deriv_phi(ed, pre, m, z, dzdp, dmdp, phi, dphidp) !! Calculat the numerical derivatives we need !! for the Runge-Kutta calculation !! derivative w.r.t. dp -> dpressure IMPLICIT NONE REAL(8) :: ed, pre, m, dedp, z, dzdp, dmdp, phi, dphidp REAL(8) :: ap, qr, elr, nur 37 real(8), parameter :: pi = 3.14159265358979 ap = 1.3234d-6 !! 1MeV/fm^3 = 1.3234E-6 /km^2 !! m -> km !! z -> km^2 if (z .lt. 1.0d-12) then dzdp = -2.0/((ed + pre))*(3.0 - 8.*pi*ed*ap*z) & /(4.*pi)/(ap*(ed + 3.*pre)) dmdp = 2.*pi*ap*ed*z**0.5*dzdp dphidp = -24./5./(ed + 3.*pre) else dzdp = -2.0/(ed + pre)*z*(z**0.5 - 2.*m) & /(m + 4.0*pi*pre*z**(1.5)*ap) dmdp = 2.0*pi*ap*ed*(z**0.5)*dzdp dphidp = (phi*(phi + 3.)*(z**0.5 - 2.*m) & - 4.0*pi*(phi + 4.0)*z**(1.5)*ap*(ed + pre)) & /(ed + pre)/(m + 4.*pi*ap*pre*z**(1.5)) endif Return End SUBROUTINE deriv_phi SUBROUTINE RK42phi(pre, dp, rho, ede, & z, mass, phi, edata, pdata, ndata, nmax) !! Runge-Kutta 4th order method with beta stable nuclear matter IMPLICIT NONE !! z = r^2 !! mass = m REAL(8) :: pre, dp, u, z, pre1, pre2 INTEGER :: NMAX REAL(8), DIMENSION(NMAX) :: edata, pdata, ndata REAL(8) :: energy, mun, mup, mue, eb, ed, ed1, ed2, ed3, ed4 38 real(8) :: ede REAL(8) :: dzdp, dmdp, mass, df INTEGER :: I, J, K, iter REAL(8) :: mass1, dm1, mass2, dm2, mass3, dm3, mass4, dm4 real(8) :: dz, z1, dz1, z2, dz2, z3, dz3, z4, dz4 real(8) :: dedp, dedp1, dedp2, dedp3, dedp4 real(8) :: rho1, rho2, rho3, rho real(8) :: phi, dphidp, dphi, dphi1, dphi2, dphi3, dphi4 real(8) :: phi1, phi2, phi3, phi4 real(8) :: df1, df2, df3, df4 real(8) :: av_ed real(8), parameter :: pi = 3.14159265358979 real(8) :: ap ! ap = 1.3234e-6 !! 1MeV/fm^3 = 1.3234E-6 /km^2 pre1 = pre + 0.5*dp pre2 = pre + dp CALL AITKEN(NMAX, pdata, edata, pre, ed, df) call deriv_phi(ed, pre, mass, z, dzdp, dmdp, phi, dphidp) dz1 = dp*dzdp dm1 = dp*dmdp dphi1 = dp*dphidp z1 = z + 0.5*dz1 mass1 = mass + 0.5*dm1 phi1 = phi + 0.5*dphi1 CALL AITKEN(NMAX, pdata, edata, pre1, ed1, df) call deriv_phi(ed1, pre1, mass1, z1, dzdp, dmdp, phi1, dphidp) dz2 = dp*dzdp dm2 = dp*dmdp dphi2 = dp*dphidp z2 = z + 0.5*dz2 mass2 = mass + 0.5*dm2 phi2 = phi + 0.5*dphi2 call deriv_phi(ed1, pre1, mass2, z2, dzdp, dmdp, phi2, dphidp) 39 dz3 = dp*dzdp dm3 = dp*dmdp dphi3 = dp*dphidp z3 = z + dz3 mass3 = mass + dm3 phi3 = phi + dphi3 CALL AITKEN(NMAX, pdata, edata, pre2, ed3, df) call deriv_phi(ed3, pre2, mass3, z3, dzdp, dmdp, phi3, dphidp) dz4 = dp*dzdp dm4 = dp*dmdp dphi4 = dp*dphidp !!-------------------------------------------------pre = pre2 ede = ed3 z = z + (dz1 + 2.*dz2 + 2.*dz3 + dz4)/6. mass = mass + (dm1 + 2.*dm2 + 2.*dm3 + dm4)/6. phi = phi + (dphi1 + 2.*dphi2 + 2.*dphi3 + dphi4)/6.0 call AITKEN(NMAX, pdata, ndata, pre2, rho3, DF) rho = rho3 Return End SUBROUTINE RK42phi Name of Candidate: Robert J. Stahulak Birth date: March 19th, 1997 Birth place: Salt Lake City, Utah Address: 1030 E 600 S Salt Lake City, UT 84102