UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Critical collapse of Newtonian fluids Aguilar-Martinez, Silvestre 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2015_november_aguilarmartinez_silvestre.pdf [ 3.07MB ]
JSON: 24-1.0165797.json
JSON-LD: 24-1.0165797-ld.json
RDF/XML (Pretty): 24-1.0165797-rdf.xml
RDF/JSON: 24-1.0165797-rdf.json
Turtle: 24-1.0165797-turtle.txt
N-Triples: 24-1.0165797-rdf-ntriples.txt
Original Record: 24-1.0165797-source.json
Full Text

Full Text

Critical Collapse of Newtonian FluidsbySilvestre Aguilar-MartinezB.Sc., California State University, Fresno, 2002M.Sc., The University of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2015c© Silvestre Aguilar-Martinez 2015AbstractThis thesis constitutes a numerical study concerning the dynamics of an inviscidfluid subject to Newtonian gravity. Type-II critical phenomena has been previ-ously measured in gravitational collapse simulations of isothermal-gas-spheres inNewtonian gravity. Our first objective was to extend this work by applying themore general polytropic-gas equation-of-state to the spherically symmetric fluid. Weshowed that under generic conditions of critical collapse, the polytropic gas allowsfor scale-invariant solutions. These solutions display self-similarity of the first kindwith non-linear scaling between the space and time variables. One of these solutionswas identified as the critical solution in critical collapse simulations. Such solutionwas found to have a single unstable mode with a Lyapunov exponent whose valuedepends on the polytropic index (Γ) from the equation of state. We argued that thisbehavior constitutes evidence of type-II critical phenomena with a transition fromtype-II to type-I behavior occurring at Γ ≥ 6/5. Thus, the polytropic gas exhibitsboth types of critical behavior. These phenomena was investigated dynamically andalso from perturbation analysis.In the second phase of the project we extended the hydrodynamic model to treataxi-symmetric gravitational collapse. This allowed us to study the effect of angularmomentum on the critical solution. As previously predicted, infinitesimal initialrotation introduces a non-spherical, unstable axial mode into the dynamics. Themeasured scaling behavior of the specific angular momentum of the collapsed coreagrees with the predicted growth rate (Lyapunov exponent) of the axial perturba-tion. This two-mode linear regime modifies the scaling laws via the introduction ofiiAbstractuniversal functions that depend on the two-parameter family of initial data. Thepredicted universality of these functions was confirmed through careful measure-ments of the collapsed mass and its angular momentum near the collapse threshold.A two-parameter space survey reveals a universal behavior of the order-parameters,with no mass-gap even in the presence of finite initial rotation. The behavior changesslightly beyond some initial rotation threshold. The results then, can be interpretedas an intermediate convergence to a non-spherical self-similar critical solution witha single unstable mode.iiiPrefaceThe work presented in this manuscript contains original research conducted by theauthor in collaboration with the research supervisor, Professor Matthew W. Chop-tuik. This contribution is summarized and highlighted in Chap. 4 and Chap. 5.Only the discussion of the spin-up mode in Sec. 5.1 had previously been published.It should be pointed out that in the exposition of the formalism displayed in Chap. 2certain equations had to be rederived to tailor therefor the purposes of this inves-tigation. In particular, the discussion on type-II critical phenomena Sec. 2.8.1 wasrevised from what was done in Phys. Rev. D, 65(064019):1-10, (2002) to treat thepolytropic gas in Newtonian gravity. We repeated Gundlach’s perturbation analysisand derived the Newtonian-specific scaling laws for the collapsed mass and specificangular momentum of the collapsed core namely, Eqs. (2.118) and (2.120). Thesepreviously unpublished equations are analogous to the scaling laws for the black holemass and its specific angular momentum derived in Phys. Rev. D, 65(064019):1-10,(2002).ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xNotation Conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xivAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xviii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.1.1 Historical Context . . . . . . . . . . . . . . . . . . . . . . . . 51.2 Critical Phenomena in Newtonian Collapse . . . . . . . . . . . . . . 111.2.1 Newtonian Analogue of the Black Hole Mass . . . . . . . . . 121.3 Goals/Objectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141.3.1 Layout of the Rest of Thesis . . . . . . . . . . . . . . . . . . 152 Formalism and Equations of Motion . . . . . . . . . . . . . . . . . . 182.1 Fluid Dynamics in Newtonian Gravity . . . . . . . . . . . . . . . . . 19vTable of Contents2.1.1 Equation of State . . . . . . . . . . . . . . . . . . . . . . . . 212.1.2 Restriction to Spherical Symmetry . . . . . . . . . . . . . . . 242.1.3 Restriction to Axial Symmetry . . . . . . . . . . . . . . . . . 272.2 Self-similarity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 282.3 Self-similar Ansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . . 302.3.1 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . 332.4 Spherically Symmetric Linear Perturbations . . . . . . . . . . . . . 352.5 Non-spherical Linear Perturbations . . . . . . . . . . . . . . . . . . 382.6 Correspondence with General Relativity . . . . . . . . . . . . . . . . 442.6.1 Regular Newtonian Self-similar Solutions . . . . . . . . . . . 452.6.2 Regular GR Self-similar Solutions . . . . . . . . . . . . . . . 482.6.3 Non-spherical Perturbations . . . . . . . . . . . . . . . . . . 492.7 Initial Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 512.8 Critical Phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . 522.8.1 Type II Critical Behavior and Scaling Laws . . . . . . . . . . 532.8.2 Type I Critical Behavior . . . . . . . . . . . . . . . . . . . . 583 Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.1 Finite Differences . . . . . . . . . . . . . . . . . . . . . . . . . . . . 603.1.1 Independent Residual . . . . . . . . . . . . . . . . . . . . . . 663.2 Finite Volumes: Conservative Methods . . . . . . . . . . . . . . . . 683.2.1 Spherical Coordinates: Spherical Symmetry . . . . . . . . . 693.2.2 Cylindrical Coordinates: Axial Symmetry . . . . . . . . . . . 703.3 The Riemann Problem . . . . . . . . . . . . . . . . . . . . . . . . . 723.4 Approximate Riemann Solvers . . . . . . . . . . . . . . . . . . . . . 753.4.1 The HLLC Approximate Riemann Solver . . . . . . . . . . . 773.4.2 The Roe Solver . . . . . . . . . . . . . . . . . . . . . . . . . 793.4.3 Cell Boundary Variable Reconstruction . . . . . . . . . . . . 803.5 Time Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81viTable of Contents3.6 The Grid Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . 833.6.1 Spherical Symmetry Non-uniform Grid . . . . . . . . . . . . 833.6.2 Axial Symmetry Grid . . . . . . . . . . . . . . . . . . . . . . 853.7 Solution to Poisson’s Equation . . . . . . . . . . . . . . . . . . . . . 863.8 Fluid Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 883.8.1 Fluid Regularity Conditions in Spherical Symmetry . . . . . 893.8.2 Fluid Regularity Conditions in Axial Symmetry . . . . . . . 903.9 Conserved Quantities and Error Diagnostics . . . . . . . . . . . . . 933.9.1 Code Validation: Spherical Symmetry . . . . . . . . . . . . . 933.9.2 Code Validation: Axial Symmetry . . . . . . . . . . . . . . . 964 Results: Spherical Symmetry . . . . . . . . . . . . . . . . . . . . . . 1024.1 Self-similar Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . 1024.1.1 Solutions, 1 ≤ Γ < 6/5 . . . . . . . . . . . . . . . . . . . . . 1064.1.2 Growing Modes . . . . . . . . . . . . . . . . . . . . . . . . . 1134.2 Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . . . 1164.2.1 Simulations Γ ≈ 1 . . . . . . . . . . . . . . . . . . . . . . . . 1174.2.2 Connection with General Relativity . . . . . . . . . . . . . . 1254.2.3 Simulations 1 < Γ < 6/5 . . . . . . . . . . . . . . . . . . . . 1274.2.4 Simulations 6/5 ≤ Γ < 4/3 . . . . . . . . . . . . . . . . . . . 1345 Results: Axial Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . 1405.1 The Unstable Axial (Spin-up) Mode . . . . . . . . . . . . . . . . . . 1425.2 Numerical Solutions in Axial Symmetry . . . . . . . . . . . . . . . . 1465.2.1 Slow Rotation (q → 0) . . . . . . . . . . . . . . . . . . . . . 1475.2.2 Finite Initial Rotation . . . . . . . . . . . . . . . . . . . . . . 1545.2.3 Rapid Initial Rotation (Large q Regime) . . . . . . . . . . . 1645.2.4 Analogy with Statistical Mechanics . . . . . . . . . . . . . . 1685.3 Numerical Solution at Larger Values of Γ . . . . . . . . . . . . . . . 171viiTable of Contents5.3.1 Critical Solutions Γ < 6/5 . . . . . . . . . . . . . . . . . . . 1716 Conclusion and Further Work . . . . . . . . . . . . . . . . . . . . . . 1746.1 Extensions and Further Work . . . . . . . . . . . . . . . . . . . . . 177Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179AppendicesA Solution’s Behavior Near Sonic Points for Γ = 1 . . . . . . . . . . 187B Polar Perturbations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188viiiList of Tables2.1 Newtonian-GR similarity solutions. . . . . . . . . . . . . . . . . . . . 494.1 Values of the parameters Q0 and xs for the similarity solutions cor-responding to Γ = 1.12. . . . . . . . . . . . . . . . . . . . . . . . . . 1064.2 The Hunter-A parameters Q0 and xs for some choices of Γ ∈ [1, 6/5). 1094.3 Results of the stability analysis performed particularly on the Hunter-A solution for the range 1 ≤ Γ < 6/5. . . . . . . . . . . . . . . . . . . 1144.4 Spherically symmetric initial data. . . . . . . . . . . . . . . . . . . . 1175.1 Initial data profiles for the primitive variables ρ, vs, vφ, vz and Pused in the axisymmetric evolutions. . . . . . . . . . . . . . . . . . . 142ixList of Figures1.1 Behavior of the order-parameter MBH near the threshold of gravita-tional collapse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61.2 The Newtonian analogue of the black hole mass . . . . . . . . . . . . 172.1 This plot illustrates the spectrum of axial perturbation modes andits dependence on Γ . . . . . . . . . . . . . . . . . . . . . . . . . . . 442.2 Sketch of the regular self-similar solutions for the isothermal gas (Γ = 1) 473.1 Typical wave pattern solution to the vector Riemann problem (this isthe generic solution for a three-component vector Riemann problem,i.e. q = (q1 q2 q3)⊺) . . . . . . . . . . . . . . . . . . . . . . . . . . . 743.2 Plots illustrating the convergence of total mass (Mtotal) and energy(Etotal) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 953.3 The independent residual convergence test for the spherically sym-metric fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 973.4 Plots of convergence of the conserved quantities for the axi-symmetricfluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 983.5 Independent residual convergence tests for the axi-symmetric fluidmodel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1004.1 Plot of u(x)/x for various choices of integration parameters Q0 and xs 1044.2 Phase space plot for the computed solutions α(x) and u(x) at x = 0.4 1054.3 Plot of self-similar solutions α(x) corresponding to Γ = 1.12 . . . . . 107xList of Figures4.4 Plot of self-similar solutions u(x)/x corresponding to Γ = 1.12 . . . . 1084.5 Hunter-A solution for α(x) in the range 1 ≤ Γ < 6/5 . . . . . . . . . 1094.6 Hunter-A solution for u(x)/x in the range 1 ≤ Γ < 6/5 . . . . . . . . 1104.7 Larson-Penston solution for α(x) in the range 1 ≤ Γ ≤ 6/5 . . . . . . 1114.8 Multiple plots of the Hunter-A solutions at various values of Γ closeto the critical value of Γ = 6/5 . . . . . . . . . . . . . . . . . . . . . 1124.9 Linear perturbation functions δα(x) and δu(x)/x for the Hunter-Asolution where Γ = 1.1999 . . . . . . . . . . . . . . . . . . . . . . . . 1154.10 Plot of the similarity parameter Q0 versus central density ρ(t, 0) forvarying degrees of fine-tuning . . . . . . . . . . . . . . . . . . . . . . 1184.11 Snapshots of the evolution of α(x) for critical initial data using model-A, at Γ = 1.00001 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1214.12 Snapshots of the evolution of u(x)/x for critical initial data usingmodel-A, at Γ = 1.00001 . . . . . . . . . . . . . . . . . . . . . . . . . 1224.13 Scaling behavior of the collapsed mass (M) and maximum centraldensity ρmaxc at Γ = 1.00001 . . . . . . . . . . . . . . . . . . . . . . . 1254.14 Evolution of Q0 versus the central density ρ(t, 0) at multiple valuesof Γ < 6/5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1284.15 Plots of the dimensionless density variable α(x) emerging from crit-ical initial data (p ≈ p⋆) for Γ = 1.04, 1.08, 1.12, 1.16 taken atintermediate times . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1304.16 Plots of the dimensionless density variable α(x) emerging from criticalinitial data (p ≈ p⋆) for Γ = 1.04, 1.08, 1.12, 1.16 taken at late times 1314.17 Scaling behavior for the collapsed mass at Γ = 1.04, 1.08, 1.12, 1.16,and 1.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1324.18 Plots of the collapsed massM with initial data resembling the Hunter-A solutions at Γ = 1.04, 1.08, 1.12, 1.16, respectively . . . . . . . . . 1334.19 Plot of various measurements of Q0 versus ρc for Γ = 1.28 . . . . . . 135xiList of Figures4.20 Snapshots of the critical evolution of model-A for Γ = 1.28 . . . . . . 1364.21 Scaling of the critical solution’s lifetime (T0) for various choices ofΓ > 6/5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1384.22 Plots of M(collapsed) and the solution’s lifetime divided by the totalmass for two distinct 1-parameter families of initial data at Γ = 1.28 1395.1 Radial function plot of the spin-up mode (δuΦ(x)/x) for the Hunter-Asolution at five values of Γ . . . . . . . . . . . . . . . . . . . . . . . . 1445.2 Plots of the spin-up mode’s azimuthal-component velocity field vφ forΓ = 1.00001 and Γ = 1.12 . . . . . . . . . . . . . . . . . . . . . . . . 1455.3 Plot of Q0 versus central density (ρ(t, 0, 0)) for critical initial data . 1485.4 Critical-evolution measurements of vφ for the spin-up mode, withΓ = 1.00001 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1515.5 The scaling behavior of the collapsed mass (M) and its specific an-gular momentum (a) for supercritical initial data near the collapsethreshold p⋆ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1525.6 Scaling behavior of the collapsed mass M , and the specific angularmomentum a, for initial data that closely resembles the Hunter-Asolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1535.7 Plot of function F (δ) for three different 2-parameter families of initialdata . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1565.8 Plot of function G(δ) for three different 2-parameter families of initialdata . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1575.9 Universal functions F (δ) and G(δ) near the collapse threshold δmax . 1595.10 Curve of the gravitational collapse threshold on the p− q parameterspace for initial data model-A at Γ = 1.00001 . . . . . . . . . . . . . 1615.11 Survey of collapsed core masses M in the p− q parameter space . . 1625.12 Survey of the collapsed core’s specific angular momentum a in thep− q parameter space . . . . . . . . . . . . . . . . . . . . . . . . . . 163xiiList of Figures5.13 Displays of vφ corresponding to the critical solutions for the case of“large” initial rotation (q = 0.5) with Γ = 1.00001 . . . . . . . . . . 1665.14 Measurements ofM and a near the threshold of gravitational collapsepcr for data with “large” initial rotation (q = 0.5) with Γ = 1.00001 . 1675.15 Calculations in axi-symmetric critical fluid collapse with Γ = 1.12and slow initial rotation (q = 10−14) . . . . . . . . . . . . . . . . . . 173xiiiNotation ConventionsDuring the writing of this thesis we tried to be as consistent as possible with thenotation and definitions previously used in related works. In particular, the em-ployed notation and convention follows from the work of the most repeatedly citedauthors. As it will be made clear, this practice will assist in the comparison of thephysical results. We provide a partial list of the symbols which we anticipate couldpotentially generate confusion while reading this thesis.G Newton’s constant, set to unity in all numerical calculations.x Similarity variable, r/tnn Scaling exponent, 2− Γτ Zooming time coordinate,− ln(1− t/t0)Q Convergence factorQ0 Similarity variable related to central density, ln(4πGρc(t0 − t)2)q 1×N array of conservative variables (continuum case), denotes(q1, . . . , qN )⊺. All bold-faced symbols represent 1-dimensional arrays.p 1×N array of primitive variables (continuum case), denotes(p1, . . . , pN )⊺Q Array of conservative variables (discrete case)P Array of primitive variables (discrete case)~x Position vector (x1, x2, x3). All quantities with~. denote spatial vectors.xivNotation Conventionsq Magnitude of vector factor on the azimuthal velocity used in setting initialangular momentum |~q|δ Magnitude of small vector quantity |~δ| related to the initial angular momentuma Magnitude of the specific angular momentum of the collapsed core, |~a|M Mass of the collapsed coreF (~δ) Universal function related to the mass of the collapsed core~G(~δ) Universal function related to the specific angular momentumof the collapsed core. In axisymmetry ~G(~δ) = G(δ)zˆxvAcknowledgementsI would like to express my sincere gratitude to my research advisor Professor MatthewChoptuik. I am truly indebted to him for the guidance he provided throughout thisresearch project. I specially appreciate his highly rigorous and critical style of teach-ing while granting his students a great degree of freedom in the pursuit of researchideas. I feel truly honored for being a part of his research group during these years.I also would like to extend special thanks to the members of my research committeeProfessors, Jeremy Heyl, William Unruh and Hirohisa Tanaka.I am obliged to professor Frans Pretorius; the programs PAMR and AMRDthat he developed, played a significant role in expediting our code development andcalculations. I also owe many thanks to Scott Noble; many of the ideas used inthe development of our spherically symmetric hydrodynamic code were inspired byand built upon his work of the subject. He [Scott] was also very generous with histime in answering many of my hydro questions via email. I also owe a great dealto Professor Carsten Gundlach; much of our work is motivated by his studies oncritical phenomena. I would like to specially thank him for our Skype conversationwhich help me clarify important concepts pertaining to this research.Moreover, I would like to take this opportunity to thank my friends and col-leagues at UBC: Daoyan Wang, Arman Akbarian, Graham Reid, Jason Penner,Bruno Mundim, and Ben Gutierrez. Great gratitude is extended to Lynn Bagloleand Shanda Walters. I thank them for their encouragement and for reading part ofthis manuscript. I also want to use this chance to thank my high school counselorDennis Conner and my friend Robert Bergen. Their support and generosity arrivedxviAcknowledgementsat a critical time in my life and help set into motion the events which led to thewriting of this dissertation. Furthermore, I thank my family for their patience andunwavering support.I am profoundly grateful to The Natural Sciences and Engineering ResearchCouncil of Canada (NSERC) which provided the financial support that made thiswork possible.xviiDedicationEste trabajo esta´ dedicado con gran afecto a mi padre, Alvaro y a la memoria demi madre, Delfina.xviiiChapter 1IntroductionEinstein’s theory of General Relativity (GR) as well as Newton’s law of gravityprovide a mathematical framework that successfully model the effect of gravitationwhere matter-energy act as a source. An important consequence of both theories isthe production of compact states of matter through the mechanism of gravitationalcollapse [51, 70, 104]. Some of the most exciting areas of research in astronomyand astrophysics involve just such a scenario. Examples of these are neutron stars,black holes, active galactic nuclei, supernovae, collapsars, etc. Much effort has beenallocated towards the investigation of these systems using facilities such as γ-ray andX-ray telescopes [86, 110]. Parallel to this development a more theoretical approachhas emerged with the goal of solving the dynamical equations (GR/Newtonian)using numerical techniques in an attempt to simulate/model these phenomena. Theliterature on the subject is too immense to list, therefore, we direct the reader to ashort compendium of reviews and applications of the methods to various physicalsystems, [7, 20, 23, 25, 26, 29, 67, 92, 100, 107]. Numerical methods have provento be instrumental in advancing our understanding of the process of gravitationalcollapse.The perfect fluid model is a subclass of the matter models often used in numer-ical simulations of astrophysical and cosmological phenomena. This model ignoresnon-adiabatic effects such as viscosity and radiative transfer [3, 25, 26, 67]. In in-flationary cosmology the perfect fluid is used to model the contribution of radiationenergy to the stress-energy tensor, see [62, 79, 106] and references therein. It alsoserves as a basic model of stellar structure, see discussion on barotropic fluids in11.1. Overview[58, 80]. The perfect fluid is also widely used in investigations of a more fundamen-tal nature associated with the interaction of gravity and matter. Its ramificationsare far reaching in fields such as critical phenomena, naked singularity formation,cosmic censorship, and quantum gravity. The project presented herein is of thelatter theoretical orientation. Our aim was to consider an idealized fluid model cou-pled to gravity in order to investigate critical gravitational collapse using numericaltechniques. This was done in an attempt to shed more light on the surprising criticalphenomena which is known to emerge at the collapse threshold [10].1.1 OverviewIn numerical calculations of the collapse of a massless real scalar field in GR, Chop-tuik (1993) discovered intricate, small-scaled, periodic structures on the scalar fieldprofile that developed at the threshold of gravitational collapse [15]. This behav-ior would later be called Critical Phenomena in Gravitational Collapse. To date,a myriad of other matter models have been studied, with some showing a similarcritical behavior. An up to date list of these can be found in Gundlach’s review ofthe subject [10]. Following Choptuik’s discovery, Evans and Coleman (1994) identi-fied similar behavior in spherically symmetric critical collapse of a photon gas [22].The Evans and Coleman critical solution, like Choptuik’s counterpart was found tobe scale-invariant (self-similar). This however, displayed a continuous symmetry,unlike the “periodic echoing”, discrete symmetry found in Choptuik’s solution. Inspite of these differences, a general description of critical phenomena can be made ina model independent way. In what follows, we provide a brief description of criticalphenomena in gravitational collapse, focusing on the perfect fluid model. The morecomplete formalism is given in Sec. 2.8.Initial conditions are imposed on the matter fields via a set of parameters e.g.temperature, width of matter distribution, amplitude of the energy density, etc.Numerical integration of Einstein’s equations and the equations of motion for the21.1. Overviewmatter fields yield the following empirical facts; the matter-energy initially “im-plodes” under the influence of gravity–the spatial profile of the energy distributionnarrows effectively strengthening the gravitational interaction, which in some caseswill lead to singularity formation, e.g. a black hole in General Relativity [104]. Ifwe allow control of the initial data through a single control parameter p, i.e. a1-parameter family of data, one discovers that the outcome of the evolution can bedivided into two, sometimes three qualitatively different states.1. For a certain range of the parameter p the initial implosion is followed by an“explosion”, an effect that leads to the dispersion of the matter-energy awayfrom the center. We call these data subcritical.2. If the implosion continues leading to singularity/black hole formation; we callthis the supercritical regime.3. After some time the evolution settles into a static configuration, a star-likestate.If generic choices of p yield end-states described only by 1 or 2 (collapse/dispersal)for a particular matter type under investigation then, we find that these results areassociated with two regimes in the domain of p. For example, if p > p⋆ where p⋆is some threshold, then the end-state corresponds to complete collapse (black holeformation) regardless of the difference p − p⋆. Similarly, for p < p⋆ the evolutionleads unambiguously to dispersal of the matter. We can the tune the initial data tothe threshold of singularity formation, by setting p = p⋆. The spacetime and matterfields resulting from this fine-tuning is known as the critical solution. This solutionhas been found to possess qualities which are reminiscent of phase transitions instatistical mechanics [10, 34].In analogy with phase transitions in statistical mechanics the critical solutionis said to belong to two distinct types. Type-I critical behavior refers to staticmetastable configurations of the matter-energy. The solution is time-translation31.1. Overviewsymmetric. This symmetry can be continuous, that is the solutions to the equationsof motions Z(r, t) is invariant under Z(r, t) → Z(r, t + ∆t) for any ∆t ∈ R, ordiscrete i.e. invariant under Z(r, t) → Z(r, t + ν∆t′) where ∆t′ is a specific periodcharacteristic of the solution, and ν ∈ Z. However, this is only true if we have fine-tuned the, initial condition to precisely the threshold parameter p⋆. Experimentally,it is impossible to fine tune the initial state to p = p⋆ due to limitations in numericalprecision. The critical solution is extrapolated from the behavior near the collapsethreshold. In particular, the time the computed solution converges to the criticalsolution as a result of fine tuning p to p⋆ follows the scaling lawt ∼ −σ ln |p− p⋆|, (1.1)where t is the life-time of an almost static state which resembles the critical solutionand σ, called the scaling exponent, is a constant whose significance is discussedin Sec. 2.8. Furthermore, if we explore the supercritical regime (p close to p⋆ forwhich singularity formation is the outcome) the mass of the black hole (MBH) has aminimum, finite value, thus there exists a “mass-gap” in the size of the possible blackholes. This is the reason for the label Type-I critical phenomenon. MBH, playingthe role of the order parameter, changes discontinuously at the critical value p = p⋆,and thus reminds us of first-order phase transitions in statistical mechanics. See topplot in Fig. 1.1 for a typical example.For the case of Type-II critical phenomena, arguably the more interesting case,the solution is scale symmetric (self-similar), where again, the symmetry can beeither continuous or discrete. The scaling law now involves the order parameter,MBH such that,lnMBH ∼ γ ln |p− p⋆|, (1.2)for values of p in the supercritical regime, γ is again called the scaling exponentdiscussed in Sec. 2.8. In this case there is no mass gap in the spectrum of MBH,41.1. Overviewsee bottom plot in Fig. 1.1. The transition between flat (dispersal) and singularfinal spacetimes can be made continuously via the critical solution, and hence it isanalogous to second-order phase transitions in statistical mechanics [10, 34].Another property of the critical solution is universality. It can be verified (bynumerical integration of the governing equations) that the critical solution and thescaling exponent (γ), are independent of the details of, and choice of 1-parameterinitial data family. Evidence of universality has been found in many models thatdisplay type-II critical behavior [9, 15, 22, 52, 55, 73, 103]. This property of thecritical solution can be understood in terms of linear perturbation theory as out-lined by [99]. The single-unstable-mode structure of the critical solution makes itan intermediate attractor of the critical evolution. The unstable mode grows ac-cording to 1/(t0 − t)Re(λ0), where t0 is the time of singularity formation and λ0 isthe Lyapunov exponent of the perturbation mode. It can be shown that γ = 1/λ0.This formalism will be discussed in more detail in Sec. Historical ContextThis brief historical review is focused on those studies involving the perfect fluidmodel. Its presentation is aimed at establishing the relevance of our project. Type-II critical phenomena was observed in critical collapse simulations of a perfect fluidthat obeys an ultrarelativistic equation of state (EoS),P = kρ, (1.3)where the variables P and ρ represent the pressure and energy density, respectively,with k being a constant of proportionality [22, 47, 73, 74, 96]. Evans and Colemanfound the critical solution for k = 1/3 (photon gas). Calculation of the black holemass resulting from supercritical evolutions of the radiation fluid found that the massfollowed a power-law similar to that of Choptuik’s scalar field. The black hole masswas shown to obey Eq. (1.2), with scaling exponent γ ≈ 0.36. Koike et al. (1994)51.1. OverviewFigure 1.1: Behavior of the order-parameterMBH near the threshold of gravitationalcollapse. This plot illustrates the typical behavior of MBH in the vicinity of thecritical parameter (p = p⋆). A gap exists in the possible values of MBH (top) acrossthe collapse threshold parameter p⋆ whereas for type-II transitions MBH “turns on”at arbitrarily small values (bottom) as p → p⋆− (p approaches p⋆ from the left inthe given example).61.1. Overviewstudied linear perturbations about the Evans-Coleman solution and discovered ithad a single growing mode. From their result, a more precise calculation of thescaling exponent was possible. They found, γ ≡ 1/λ0 ≈ 0.3558019. It was soonestablished that universality did not extend to including different values of k fromEq. (1.3) [66]. Nielsen and Choptuik (2000) [73, 74] found continuously self-similarcritical solutions for selected values of k ∈ (0, 1]. They measured, for instance, mass-scaling exponents γ ≈ 0.15 at k = 0.05 and γ ≈ 1.0 at k = 1. Their work includedspherically symmetric collapse simulations and solutions to the autonomous systemof equations obtained by assuming a self-similar ansatz.This work was extended by Noble (2008) [75] in his studies of the ideal gas EoS,P = kρ0ǫ. (1.4)The variables ρ0 and ǫ are respectively the rest-mass energy density and specificinternal energy, and k is again a constant parameter. Noble (2008) carried out sim-ulations of critical collapse of an ideal gas in spherical symmetry [75]. It was revealedthat the critical solution was identical to that previously found by Nielsen et. al.[73, 74] using the simpler ultrarelativistic EoS (1.3). These results were interpretedheuristically as the tendency of the collapsing ideal gas towards an ultrarelativisticlimit, modeled by P = kρ. This follows from the fact that for rapidly collapsingmatter, ρ = ρ0 + ǫρ0 ≈ ǫρ0.Ori and Piran (1990) studied self-similar solutions of a spherically symmetric per-fect fluid with EoS (1.3) in GR. They showed that the limit k → 0 (Eq. (1.3)) yieldsthe weak-field approximation—the solutions reduce to those describing, isothermalself-similar flows in Newtonian gravity [77]. This Newtonian system will be discussedin detail in Sec. 2.3. Investigations involving the Newtonian isothermal self-similarfluid model date back the work of Larson (1969)[60] and independently Penston(1969) [81]. The regular (analytic) self-similar solution found simultaneously byLarson and Penston, later named the Larson-Penston (LP) solution, was proposed71.1. Overviewas a candidate solution to describe the early stages of gravitational collapse of molec-ular clouds. Shu (1977) [95] found a static self-similar solution which is singular atthe origin and used it as initial data in collapse simulations. The work of Hunter(1977) [56] provided a more complete picture of the spectrum of analytic self-similarsolutions to the isothermal gas model. He [Hunter] applied a graphical techniquewhich allowed him to determine in principle all of the regular self-similar solutionsto the isothermal gas. He found the solution spectrum to be discrete and infinite.This technique was also applied in this work, details of its application are given inSec. 4.1.1. The solutions were shown to form a hierarchical structure Sec. 2.6.2.The first member in the hierarchy is the previously known Larson-Penston solution,with subsequent solutions labeled Hunter-A,B,C,... Extrapolation of this Hunter-spectrum hints at a convergence to the static singular solution found by Shu [95].More details of these solution are provided in Sec. 2.6.2 and Chap. 4.Similar to the work done by Nielsen and Choptuik [73, 74], Harada et al. per-formed spherically-symmetric dynamical simulations of critical collapse of a perfectfluid in GR with a focus on the near-Newtonian regime, specifically k ∈ (0, 0.036] inEq. (1.3) [47]. They showed that the type-II critical phenomena previously measuredcontinued to k → 0 (the Newtonian limit). The results of the numerical experimentswere interpreted as follows. The dynamical solution of generic initial data was shownto have late time (t0 − t → 0) convergence to a general relativistic self-similar so-lution. This solution had been previously studied by Ori and Piran by solving theautonomous system (one obtained by assuming a self-similar ansatz for the met-ric and fluid fields in spherically-symmetric, Schwarzschild-like spacetime) [77]. Werefer to this as the Ori-Piran solution. Fined-tuned initial data (critical solution)showed intermediate convergence to a different self-similar solution. This criticalsolution had similar properties to the critical solution found by Evans and Coleman[22] for the radiation fluid. Here, we use the label Evans-Coleman solution to refer tothe critical solution for arbitrary k ∈ (0, 1]. Like Koike et. al.’s stability analysis of81.1. Overviewthe radiation gas [99], Harada (2001) [44] carried out a similar linear stability workon the Ori-Piran and Evans-Coleman similarity solutions for k ∈ (0, 0.036]. Theirresults failed to find any unstable modes for the Ori-Piran solution, thus providedfurther support towards its role as a “global attractor solution” during gravitationalcollapse. This was clearly in agreement with the dynamical calculations and the latetime behavior of the solution. On the other hand, the Evans-Coleman solution (thecritical collapse solution for the General Relativistic perfect fluid calculated in [22]),according to Harada et. al.’s analysis [44] indicated the presence of a single unstablemode. The Lyapunov exponent of this mode was related to the scaling laws of theorder parameters e.g. the black hole mass (MBH). From these results a value ofthe scaling exponent γ ≡ 1/λ0 was computed for different values of k near k = 0.The value of the unstable mode at k = 0 was extrapolated to be 1/λ0 ∼ 0.11 [44].Based on these results Harada et. al. [44, 47] predicted that a similar type-II criti-cal phenomena would emerge for a strictly Newtonian calculation with appropriatere-definition of the order-parameter given that there are no black holes of the New-tonian gravity. This analogous quantity is the collapsed mass M(collapsed), whichlike the black hole mass, MBH is proportional to a length scale. Thus, the collapsedmass, is the order-parameter to be measured in Newtonian critical collapse. Thedefinition of M(collapsed) or simply M and further details are given in Sec. 1.2.1.Linear stability analysis of the Newtonian isothermal gas system conducted in[64] confirmed that the first member of the Hunter series [56], namely, the Hunter-Asolution contained a single unstable mode. The reciprocal of the growth rate forthe unstable mode was calculated to be 1/λ0 ≈ 0.10567, in agreement with theestimated value of 0.11, [44, 47]. Given the single growing mode property of theHunter-A solution Maeda et. al. (2001) [64] predicted its role as a critical solutionin Newtonian critical collapse of isothermal spheres. In contrast, calculations relatedto the stability of the Larson-Penston solution indicated an absence of any unstablemodes, suggesting its potential role as a universal attractor during the evolution91.1. Overviewof collapsing isothermal gas spheres [64]. Dynamical simulation of critical collapseof an isothermal gas in Newtonian gravity confirmed the Hunter-A solution as theNewtonian analogue of the Evans-Coleman solution [48]. Harada and Maeda (2003)also showed that the Larson-Penston solution is a universal attractor and thus isthe Newtonian analogue of the Ori-Piran solution [48]. The meaning of “criticalphenomena in Newtonian collapse” is made more concrete in Secs. 1.2 and 1.2.1.The connection with the GR self-similar solutions from the perspective of dynamicalsimulations in the Newtonian limit (k → 0 in Eq. (1.3)) was made clear thanks tothe work of Snajdr (2006) [96]. Quadruple precision numerics were used to calculatespherically symmetric critical collapse evolutions of a perfect ultra-relativistic fluidin GR. The high level of precision in the calculation of the numerical solution resolvedthe convergence to the Newtonian solutions as k → 0, in accordance with the initialpredictions of Ori and Piran [77]. The numerical calculation of the scaling lawsshowed convergence to the Newtonian results of [44, 47, 48, 64].Beyond spherical symmetry Gundlach has argued that the critical collapse so-lution of a perfect fluid that behaves in accordance with EoS (1.3) and is endowedwith initial infinitesimal angular momentum contains a second unstable mode inthe vanishing k regime [10, 24, 30–34]. The second, non-spherical1 growing mode isattributed to an ℓ = 1 axial perturbation. Gundlach (2002) [32] calculates modifi-cations to the scaling laws (e.g. MBH) taking into account a non-spherical, unstableaxial perturbation mode using the formalism of Koike et. al. Specifically, he calcu-lated the effect that the axial mode would have on the scaling of the BH mass (MBH)and its angular momentum (~LBH). These calculations constitute predictions for thecritical behavior of the perfect fluid (with EoS (1.3)) near the collapse threshold.Thus, the modified scaling laws are a test for the presence of this extra unstablemode. A discussion of the non-spherical perturbations related to our system is foundin Sec. 2.5. The consequences of any non-spherical unstable mode(s) in our model1The second mode can only come via non-spherical perturbations, since by definition of thecritical solution there can only be one unstable spherical mode101.2. Critical Phenomena in Newtonian Collapseare examined in Sec. 2.8.1.Similar non-spherical linear perturbation analyses about the Larson-Penston so-lution has been carried out in Newtonian fluids [38–42]. Simulations of gravitationalcollapse of an isothermal gas in Newtonian gravity have identified the growth of ax-ial and polar (the bar-mode) non-spherical perturbations [68, 69]. The bar-modeperturbation (polar perturbations with ℓ = 2, Appx. B) is believed to be responsi-ble for the process of fragmentation and formation of binary star systems [68]. Thespectrum of all axial perturbations can be computed explicitly for the isothermalideal gas [41], where it is clear that all self-similar solutions are subject to axialinstabilities (Sec. 2.5). Axisymmetric evolutions of the Newtonian fluid would re-press any ℓ = 2, m 6= 0 bar mode instabilities since these are φ-dependent functions,(φ is the azimuthal coordinate in polar spherical coordinates (r, θ, φ)) and axisym-metry requires that all of the fields are independent of φ. The only potentiallyobservable ℓ = 2 polar perturbation corresponds to m = 0, however, for reasonsdiscussed in Sec. 2.7 this mode is also eliminated given the symmetry restrictionsimposed on the initial data. Therefore, it is expected that axisymmetric evolutionsof critical collapse in Newtonian gravity have two unstable perturbation modes, theknown spherical one found in [64] and the axial mode [40, 41]. The effect of thisnon-spherical axial perturbation on the critical phenomena measured in [48, 64], weargue, is analogous to that predicted by Gundlach in [32]. A major part of thisproject involves this investigation.1.2 Critical Phenomena in Newtonian CollapseBefore we can properly discuss critical phenomena in Newtonian collapse some con-cepts must be clarified. We considered the system described by a spherically sym-metric isothermal gas. For concreteness, it is assumed that Temperature T (equiv-alently the internal energy) is the control parameter as was done in [48]. Knownempirically from numerical calculations performed by Harada et. al. [48] that for111.2. Critical Phenomena in Newtonian CollapseT < T ⋆ the end result is a singular solution (the supercritical regime). This is char-acterized by an exponential growth of the density and pressure at the origin. Theparameter T ⋆ represents a critical temperature. Close examination of this singularsolution reveals a tendency towards self-similarity. It becomes apparent that thisself-similar end-state is described by the Larson-Penston solution. Such end-statesare analogous to the General Relativistic self-similar black hole solution. More issaid about this association throughout Sec. 2.6. The quantity analogous to the blackhole mass in Newtonian gravity is the collapsed mass M . The mass of the collapsedregion in Newtonian gravity is defined in Sec. 1.2.1. In the case of dispersion, i.eevolutions for which T > T ⋆, (the subcritical regime) after the initial implosion thegas then disperses out to infinity, leaving a rarefied region at the origin where thedensity and pressure tend to zero. This solution is analogous to flat spacetime inGR. The critical case (T ≈ T ⋆), as described in the previous section the solutiontends towards the one-mode unstable Hunter-A solution at intermediate times, laterit either converges to the LP solution, or disperses. This is interpreted as a conse-quence of the growth of the unstable mode just as in the general relativistic case.The properties of the critical solution give rise to the features of the system whichare characteristic of type-II critical phenomena, e.g. scaling laws and universality.1.2.1 Newtonian Analogue of the Black Hole MassParticular singular solutions in Newtonian gravity play an analogous role to theblack hole solution of General Relativity. In the specific case of an isothermal gasan example of these singular solutions is the self-similar Larson-Penston solution.Similar to General Relativity, certain regions of the parameter space of initial con-ditions will lead to gravitational collapse in purely Newtonian gravity. The fieldsdescribing the fluid diverge exponentially at the origin. Like black hole formation inGR, the Newtonian collapse results in the formation of collapsed cores or compactobjects. Naturally, just like the black hole mass, we should be able to associate a121.2. Critical Phenomena in Newtonian Collapsemass quantity to these objects. Borrowing the definition used in Harada (2003) [48],the mass of the core is defined as the integral of fluid elements over the innermostregion with “in-falling” radial velocity. In spherical coordinates, with the origin setat the center of collapse, the fluid elements in the collapsing region have negative ra-dial velocity. The radius at which the velocity field changes sign defines the boundsfor the integration, and the physical size of the collapsing core. This is illustrated isFig. 1.2(a). This innermost region of negative radial fluid velocity forms at the earlystages of collapse, its evolution can be tracked over time as shown in Fig.4.13(a).The core’s mass stops evolving as the collapsing fluid approaches the time of sin-gularity formation t0 (Fig. 4.13(a)). The mass of the core approaches a fixed valuethat depends on the initial conditions, this is what we called the collapsed mass Mthe quantity which is analogous to the black hole mass MBH. More is said aboutthis tendency to a fixed value of the collapsed mass as t→ t0 in Sec. 4.2.1.This situation is similar in axisymmetry except that the velocity in the core’sregion is not purely radial, however, the core can still be defined by a negativeradial component of the velocity field. In cylindrical coordinates (s, φ, z) where, sis the distance from the axis of symmetry (the z-axis), the fluid velocity is given by~v(t, s, φ, z) = (vs, vφ, vz). Explicitly in terms of the coordinates, vs ≡ s˙, vφ ≡ sφ˙and vz ≡ z˙, where ˙≡ ddt . The radial velocity in these coordinates is then,vr(s, z) =svs + zvz√s2 + z2. (1.5)So, the sign of vs and the coordinate z determine whether or not the velocity is“in-falling”. Once we have defined the core the measurement of the mass M andspecific angular momentum ~a follows by integrating over the core’s region.131.3. Goals/Objectives1.3 Goals/ObjectivesThe general goal of this research project was to generate numerical solutions ofthe Euler equations of fluid dynamics (Chap. 2) coupled to Newtonian gravity inspherical and axial symmetry. The purpose of this undertaking is to study crit-ical gravitational collapse. The first objective that we set out to achieve was toextend the spherically-symmetric-fluid work in [48, 64] to accommodate for a moregeneric/realistic EoS i.e. the polytropic ideal gas law (2.18). This was done by nu-merical computation of the solutions to the dynamical equations in spherical symme-try under conditions of critical collapse. Furthermore, assuming a self-similar ansatzand the simpler barotropic EoS (2.22) the fluid equations can be transformed intoa system of ordinary differential equations (ODEs). We also looked for analyticself-similar solutions of this system, and studied their linear stability. The relation,between the dynamical and self-similar solutions was analyzed. The role of theadiabatic index (Γ) in both the dynamical solutions and the spectrum of similaritysolutions was also investigated.The second phase of the project involved extending our project to treat axisym-metric fluids. Given the significant increase in computation, very little work hasbeen done in the form of dynamical simulations of critical collapse beyond spher-ical symmetry. One notable exception is the critical collapse of gravity waves inaxisymmetry [1]. To date, there are no calculations of critical collapse involvinga fluid dynamical model beyond spherical symmetry. Our goal was to provide thefirst such calculations, albeit in Newtonian gravity. We investigated non-sphericaleffects i.e. rotation (initial angular momentum) on the critical phenomena in New-tonian collapse first measured in [48, 64]. Given the analogy between the criticalbehaviors of the isothermal gas in Newton’s gravity and the ultrarelativistic perfectfluid in GR, we argue that this analogy extends to the addition of angular mo-mentum. Both systems are affected by the presence of an additional unstable butnon-spherical mode. We measured this effect in the Newtonian system and tested141.3. Goals/Objectivesthe predicted behavior using the modified scaling laws derived in [34], and discussedin Sec. 2.8.1. These results allow us to reflect on the parallels with the GR perfectfluid system. In particular, we considered the implications that can be drawn fromour Newtonian calculation to the full GR+fluid model in the limit where k → 0with EoS (1.3). Furthermore, the presence of a second growing mode in the criticalcollapse solution permits an extension of the analogy with statistical mechanics.This is discussed in Chap. 5. We also briefly discussed the effect of the adiabaticindex Γ on the critical behavior.1.3.1 Layout of the Rest of ThesisThe equations of motion are presented in Chap. 2, along with the theoretical basisfor our numerical experiments. Following a presentation of the general fluid model,we imposed symmetry restrictions that lead us to the spherically symmetric Eulerfluid equations. Two sets of such systems are presented, one pertaining to the moregeneric polytropic gas and the other to a constant entropy barotropic EoS. Strictlyspeaking only the latter simpler system, as discussed in Sec. 2.1.2, allows for scale-invariant solutions. The axisymmetric fluid equations are also introduced followedby a discussion of the spherically symmetric linear perturbation to the isentropicfluid equations. A brief discussion of the non-spherical linear perturbations relevantto the axisymmetric evolution is introduced. Lastly, in Chap. 2, the modified,rotation inclusive formalism for type-I and type-II “Newtonian” critical phenomenaobserved in the numerical experiments is presented.In Chapter 3 we provide a brief overview of all the numerical techniques used inthis project. This begins with presentation of finite difference techniques, followedby a derivation of the finite volume equations necessary for proper treatment ofhydrodynamics and conservation laws. We included a brief discussion of the typeof shock-capturing methods that were used. We then listed the structure of thediscrete equations in both spherical and axial symmetry, time-integration technique,151.3. Goals/Objectivesregularity and boundary conditions. Finally, code-validation data is given thatinclude the independent residual test, and evidence of the code’s compliance withthe conservation of mass, angular momentum and total energy.Chapter 4 contains all of our the spherically symmetric calculations. The dy-namical solutions obtained using the spherically symmetric code are listed and dis-cussed here. First, we presented the calculation of the self-similar solutions to theisentropic fluid equations, followed by the calculation of their unstable linear per-turbations modes. We analyzed the dependence of the results on the chosen valueof the so-called adiabatic index Γ in the EoS. Notice that both, the polytropic idealgas, P = (Γ− 1)ρǫ and the isentropic EoS P = KρΓ (with K being constant) havedependencies on the adiabatic index Γ. This is followed by calculation of the explicittime-dependent equations. These are used to simulate critical gravitational collapse.These results are compared and with the calculated self-similar solutions. Again,the role of Γ is documented.Chapter 5 is devoted to the presentation of the axisymmetric calculations. Itbegins with the presentation of the relevant non-spherical perturbation mode, i.e.the axial mode. Its dependence on Γ is also noted. Calculations of the evolution ofcritical collapse data with rotation are given. The results are compared the thoseobtained in spherical symmetry by turning “off” the angular momentum at theinitial state as a check for consistency. The axisymmetric code is used to investigatethe effect of angular momentum on the critical solution by calculation of the orderparameter e.g. the collapsed mass. The connection of these results with GR andstatistical mechanics is also discussed.Finally in Chapter 6 we include some concluding remarks regarding the impor-tance and significance of the accomplished work. We mentioned possible extensionsof the current project and further related work.161.3. Goals/Objectives(a) (b)(c)Figure 1.2: The Newtonian analogue of the black hole mass is defined in terms ofthe total mass contained in the inner-most region of the density distribution withnegative radial fluid velocity measured at t near t0, the time of singularity formation.Panel (a) sketches this situation in spherical symmetry. Early in the evolution sucha region is formed, the radial distance measured from the origin (r = 0) to thefirst point where the radial velocity changes sign (from negative to positive) definesthe size of the collapsing core rcore as seen in (a). Such a region can similarlybe defined in axial symmetry be calculating the radial fluid velocity according toEq. (1.5). Panel (b) sketches this region for a typical axisymmetric evolution, herewe highlighted the region where vr(t, s, z) < 0. The core’s region corresponds to thearea in (b) depicted in red. Again, this region can be unambiguously defined evenat the early in the evolution. Once this region is defined, calculation of its mass andangular momentum is straightforward. Panel (c) shows the evolution of the core’smass. Notice that the core’s size and mass changes with time but as its evolutionnears t0 it approaches a fixed value. This value is what is referred to as the collapsedmass M(collapsed) or simply M . This is the quantity that is analogous to the blackhole mass MBH.17Chapter 2Formalism and Equations ofMotionThis project involves a particular fluid dynamical model. The equations governingthe dynamics of fluid flows belong to a class of time-dependent partial differentialequations (PDEs) known as conservation laws. In their most general representationthey constitute a non-linear, vector-valued, multi-dimensional system of hyperbolicequations. In a coordinate system given by {~x, t} such system in index notationcan be written as,∂qk(~x, t)∂t+3∑j=1∂f jk(q)∂xj= ψk(q; ~x, t), (2.1)where q = (q1, . . . , qk, . . . , qN )⊺ is an N × 1 matrix of so-called conservative vari-ables, N denotes the dimension, i.e the number of equations. The function f j(q)is also a N × 1 matrix referred to as the physical flux along the jth direction; inthe general case this is a non-linear function of the conservative variables. Finally,ψ(q; ~x, t) is the source function.A specially simple example of a conservation law is given by the advection equa-tion in one spatial dimension,∂q0∂t+∂∂y(νq0) = 0. (2.2)In its simplest form q0is a scalar, thus Eq. (2.2) is a scalar conservation law, and νis a constant. Note that we used the variable y to denote a generic spatial dimension182.1. Fluid Dynamics in Newtonian Gravityin order to reserve x for the definition of the similarity variable given by Eq. 2.44.Equation (2.2) describes the ‘advection’ of quantity q0with constant speed ν. TheBurger’s equation,∂q0∂t+∂∂y(q202)= 0, (2.3)provides a less trivial example of a scalar conservation law, it has applications inmany areas including fluid dynamics e.g. [61]. In this case, the flux is a non-linearfunction of q0; this leads to nonuniform characteristic flow speed equivalent to q0since,∂q0∂t+ q0∂q0∂y= 0, (2.4)this in turn can lead to discontinuities in the profile of q0even if it is initially smooth.This arises from the non-linearity of the flux function, a fact that complicates thesearch for solutions to conservation laws. However, given the ubiquity of these typesof systems in many areas of physics and engineering, much effort has been devotedto solving them. This effort consists of developing adequate numerical methods todeal with possibly discontinuous solutions. The numerical approach used in thisproject is discussed in Chap. 3. In this chapter we present the particular systemof conservation laws which comprises our physical model. The particular geometryunderlying our symmetry assumptions is also presented here.2.1 Fluid Dynamics in Newtonian GravityWe considered the non-relativistic, self-gravitating, inviscid fluid model given byEuler’s equations of fluid dynamics coupled to Newtonian gravity. They correspondto the inviscid form of the more general Navier-Stokes equations. The literaturecontains many derivations of the Euler equations, one approach being the applica-tion of Newtons laws, principles of conservation and thermodynamics, the other,the slow flow, weak-field limit of general relativistic conservation laws. We willnot provide an extra derivation here, instead we refer the reader to the [28, 57] for192.1. Fluid Dynamics in Newtonian Gravitythese details. The Euler equations constitute a vector-value system of inhomoge-neous conservation laws for the mass, momentum and energy of infinitesimal fluidelements. The gravitational interaction enters the formalism as a source-term inEuler’s equations. Newton’s law of gravity is described by Poisson’s equation. Thedifferential, coordinate independent form of these is given by,∂ρ∂t+∇ · (ρ~v) =0, (2.5)∂ρ~v∂t+∇ · (ρ~v ⊗ ~v) =−∇P − ρ∇ϕ, (2.6)∂E∂t+∇ · ((E + P )~v) =− ρ~v · ∇ϕ, (2.7)with,∇2ϕ = 4πGρ. (2.8)The variables ρ(~x, t), ~v(~x, t), P (~x, t), represent the density, velocity and pressurerespectively, at location ~x and time t. These are sometimes referred to as theprimitive variables. The set formed by {ρ, ρ~v, E} are the conservative variables.The quantity ρ~v is the momentum per unit mass, and E(~x, t), defined byE = ρǫ+12ρv2, (2.9)is the internal and kinetic contributions to the total energy density of the fluid; thevariable, ǫ(~x, t) represents the internal energy per unit mass. This function is relatedto the other state variables, ρ and P via an equation of state such as,ǫ = ǫ(ρ, P ). (2.10)Finally, the system is coupled to Newton’s law of gravity, equation (2.8), whereϕ(~x, t) is the Newtonian potential. The contribution to the energy missing inEq. (2.9) is the gravitational potential energy −12ρϕ. Note that we can incorporate202.1. Fluid Dynamics in Newtonian Gravitythis term into Eq. (2.9) leading to an alternative expression for Eq. (2.7), namely,∂ET∂t+∇ · ((ET + P )~v) = −ρ∂ϕ∂t, (2.11)where, ET = ρǫ+12ρv2 + 12ρϕ.2.1.1 Equation of StateThe fluid model under consideration is comprised by the six equations (2.5)–(2.8)that relate a total of seven unknown quantities, namely ρ, ~v (three components),P , and ǫ. Clearly, an extra condition is necessary in order to uniquely specify thestate of the system. This necessary condition is supplied by the equation of state.It gives a relationship among the thermodynamic quantities, bulk properties of thesystem. In this project we considered the equations of state (EoS) corresponding toan ideal gas. For such EoS the pressure P in the gas is related to the density (ρ)and the temperature (T ) according to,P =kBmρT. (2.12)Where kB is the Boltzmann constant and m is the mass of the constituent particles.The ideal gas law can be derived from first principles from the kinetic theory ofgases [91].We further restrict our attention to polytropic gases, in which the internal energyin a fixed volume is proportional to its temperature, and,ǫ = cvT, (2.13)the constant of proportionality cv is known as the specific heat at constant volume.A process involving a small change in temperature (dT ) leads to a corresponding212.1. Fluid Dynamics in Newtonian Gravitychange in internal energy (dǫ),dǫ = cvdT. (2.14)For a process that includes expansion or contraction but is held at constant pressure,the change is internal energy contributes to the change in temperature of the gasand the work done according to the first law of thermodynamics,dǫ = cpdT − Pd(1/ρ). (2.15)The constant cp is called the specific heat at constant pressure. Note that the changein volume per unit mass is written as d(1/ρ). If we define,h ≡ ǫ+ Pρ, (2.16)where h is known as the enthalpy and since the process is held at constant pressurewe can integrate Eq. (2.15) to obtain,h = cpT. (2.17)Solving for the temperature in Eq. (2.13) and Eq. (2.17) then combining these resultsusing the ideal gas law (2.12) yields,(cp − cv)ǫ = cv Pρ.With the definition Γ ≡ cp/cv we obtain the familiarP = (Γ− 1)ρǫ, (2.18)EoS for a polytropic ideal gas.During a thermodynamic process a quantity of typical interest is the entropy.Often it is more convenient to speak of the entropy per unit mass (specific entropy);222.1. Fluid Dynamics in Newtonian Gravitythis is commonly denoted by s. Roughly speaking s is a measure of the degree of“disorder” in a physical system. From the first law of thermodynamics, a smallchange in entropy ds is related to the other thermodynamics quantities accordingto,Tds = dǫ+ Pd(1/ρ). (2.19)Using the polytropic EoS (2.18) we can eliminate ǫ and obtain,dscv=dPP− Γdρρ. (2.20)Integrating this equation and solving for the pressure leads to,P = K0es/cvρΓ, (2.21)where K0 is a constant of integration2. Eq. (2.21) is an alternative expression forthe polytropic EoS, where the specific entropy is included instead of the internalenergy. If the evolution of the fluid takes place at constant entropy, these type offluid solutions are referred to as isentropic flows. Under these conditions (s = s0 =constant) it is clear that the EoS is simply,P = KρΓ, (2.22)where K = K0es0/cv is also a constant. Thus the pressure is solely a function of thedensity; the EoS in this case is said to be barotropic. Therefore, constant entropyconditions for a polytropic gas lead to a barotropic EoS. As we will see in the nextsection isentropic flows provide a significant simplification to the dynamical model.A further special circumstance is obtained by setting Γ = 1 in the barotropic EoS2Physically, K0 = P0ρ−Γ0 e−s/cv represents a reference point from which the changes in the stateof the gas are measured. P0, ρ0 and s0 are the reference values of pressure, density and entropy.232.1. Fluid Dynamics in Newtonian Gravity(2.22). Doing this we get,P = Kρ, (2.23)which is the EoS for an isothermal gas.2.1.2 Restriction to Spherical SymmetryFor spherically symmetric flows, the coordinate free system (2.5)-(2.8), when writtenin spherical polar coordinates (r, θ, φ) becomes,∂ρ∂t+1r2∂∂r(r2ρv) = 0, (2.24)∂ρv∂t+1r2∂∂r(r2(ρv2 + P )) =2Pr− ρ∂ϕ∂r, (2.25)∂E∂t+1r2∂∂r(r2(E + P )v) = −ρv∂ϕ∂r, (2.26)1r2∂∂r(r2∂ϕ∂r) = 4πGρ. (2.27)Here, the fields ρ(r, t), v(r, t), P (r, t), and ϕ(r, t) have no angular dependence. There-fore, this system is essentially one dimensional. Nonetheless, we would like to elim-inate the factors of 1/r2 on the radial derivatives and turn Eqs. (2.24)–(2.26) into aCartesian-like conservation law with sources similar to Eq. (2.1). This is done so thatthe application of the numerical methods discussed in Chap. 3 is straightforward.The system can be cast into the form of Eq. (2.1) in at least two ways:• Variable substitution r → r3• Expansion, i.e. 1r2∂∂r(r2f)→ 2fr+∂f∂rThe variable substitutions r → r3 effectively, regularizes the gradient operator,1r2∂∂rby 3∂∂r3242.1. Fluid Dynamics in Newtonian Gravitywhich follows from application of the chain rule. Under this transformation, thesystem now has the form,∂∂tρρvE+ 3 ∂∂r3r2ρvr2(ρv2 + P )r2(E + P )v =02Pr− ρ∂ϕ∂r−ρv∂ϕ∂r , (2.28)which is that of Eq. (2.1) with the qualification that now the spacial derivative of theflux is with respect to r3. The so-called expansion of the gradient operator generatesa term that appears to be singular at r = 0, namely 2f/r. It should be noted thatf ∼ r near the origin, so the term 2f/r is regular at r = 0. This term is collectedas a geometric source term in the right-hand-side. Note that the mass conservationEq. (2.24) acquires a source term. Explicitly, the system then becomes,∂∂tρρvE+ ∂∂rρv(ρv2 + P )(E + P )v =−2ρvr−2ρv2r− ρ∂ϕ∂r−2(E + P )vr− ρv∂ϕ∂r . (2.29)Again, this model has the Cartesian-like form of Eq. (2.1). Since finite differencemethods (Chap. 3) were used to solve Poisson’s Equation (2.27) we are not concernedwith the particular form of this operator in spherical symmetry. In our studiesrestricted to spherical symmetry both expressions Eq. (2.28) and Eq. (2.29) weresolved numerically. We found no discernible difference in the results of the numericalexperiments (Chap. 4), following the use of either Eq. (2.28) or Eq. (2.29). Thissystem is supplemented with the polytropic ideal gas EoS (2.18).Fluid Equation for Isentropic FlowsFurther simplification of the equations of motion can be achieved by assuming isen-tropic conditions. In such a case the ideal gas law EoS reduces to the barotropic form252.1. Fluid Dynamics in Newtonian Gravitygiven by Eq. (2.22). This assumption reduces the number of equations needed todescribe the flow dynamics. It can be easily checked that the energy equation (2.7)becomes redundant. In addition, given spherical symmetry we can use Gauss’ lawto integrate Poisson’s equation and introduced the variable M(r, t) to denote thetotal mass contained within a sphere of radius r, thus we arrive at the system,∂M∂t+ 4πr2ρv = 0, (2.30)∂ρv∂t+1r2∂∂r(r2(ρv2 + P )) =2Pr− GMρr2, (2.31)∂M∂r= 4πr2ρ. (2.32)Notice that we have chosen to work with the variable M(r, t) instead of the Newto-nian potential ϕ(r, t). This was done in order to ease comparison with previous workwhere the gravitational interaction in the isentropic fluid model is given in termsof M(r, t). In particular, our notation matches that used in [98]. The Newtonianpotential however, can it be easily computed from,∂ϕ∂r=GMr2. (2.33)As it will be explicitly evident in Sec. 2.2 the system of fluid equations (2.30)–(2.32) which is only valid for a barotropic EoS, i.e. (2.22) allows for particularscale-invariant (self-similar) solutions.262.1. Fluid Dynamics in Newtonian Gravity2.1.3 Restriction to Axial SymmetryUsing a cylindrical coordinate system (s, φ, z) and assuming axial symmetry, Eu-ler+Poisson equations take the form,∂ρ∂t+1s∂∂ss(ρvs) +∂∂z(ρvz) = 0, (2.34)∂∂t(ρvs) +1s∂∂ss(ρv2s + P ) +∂∂z(ρvsvz) =ρv2φ + Ps− ρ∂ϕ∂s, (2.35)∂∂t(ρvφ) +1s∂∂ss(ρvsvφ) +∂∂z(ρvφvz) = −ρvsvφs, (2.36)∂∂t(ρvz) +1s∂∂ss(ρvsvz) +∂∂z(ρv2z + P ) = −ρ∂ϕ∂z, (2.37)∂E∂t+1s∂∂ss(E + P )vs +∂∂z(E + P )vz = −ρ(vs∂ϕ∂s+ vz∂ϕ∂z), (2.38)1s∂∂ss∂∂sϕ+∂ϕ∂z2= 4πρ. (2.39)All of the fields are functions of s, z and t. Here, we use the variable s to representthe distance from the axis of symmetry, namely the z-axis. Normally, ρ is used, butin our case that variable is reserved for the fluid density. Similar to the sphericallysymmetric case we proceed to write the axisymmetric system in the Cartesian-like form of Eq. (2.1) with sources. The left-hand-side has the form of a generictwo-dimensional conservation law in Cartesian coordinates while geometric termsgenerated by expanding ∂s(sf) are collected on the right-hand-side as source terms.272.2. Self-similarityIn vector notation, the axisymmetric fluid equations are,∂∂tρρvsρvφρvzE+∂∂sρvsρv2s + Pρvsvφρvsvz(E + P )vs+∂∂zρvzρvzvsρvzvφρv2z + P(E + P )vz=−ρvss−ρ(v2s − v2φ)s− ρ∂ϕ∂s−2ρvsvφs−ρvsvzs− ρ∂ϕ∂z−(E + P )vss− ρ(vs∂ϕ∂s+ vz∂ϕ∂z). (2.40)Finite volume, and multigrid numerical methods were used to obtain approximatesolutions to the above axisymmetric system, Eq. (2.40), and Poisson’s Eq. (2.39).Just as was done in spherical symmetry the polytropic EoS (2.18) was used tocomplement this fluid-dynamical model.2.2 Self-similarityThe property of self-similarity can often be found in physical systems that ex-tend across different space and time scales. The concept has a well defined math-ematical meaning (see Buckingham’s π theorem [4]). Heuristically, the property ofself-similarity refers to a system (physical, mathematical e.g. geometrical 3) thatresembles itself at different scales. In mathematics, solutions to time-dependentpartial differential equations may display self-similarity. This is understood as aconsequence of the scale invariance of the dynamical equations, i.e. the system is3Examples of geometrical self-similarity are fractal objects such as the Sierpinski triangle or theKoch snowflake.282.2. Self-similarityinvariant under a rescaling transformation. Two important examples of systemswhich accept scale-invariant solutions are the Euler’s equations of gas dynamics inNewtonian gravity and the general relativistic perfect fluid. We must qualify thatin both cases only with an appropriate choice of EoS for the fluid can the solutionsto be rendered scale-invariant. In the case of the General Relativistic perfect fluid,only the barotropic EoS (1.3) admits self-similar solutions [65].We focused on the type of self-similarity that is involved in the solutions to ahydrodynamic system coupled to Newtonian gravity. Under certain circumstancesthe solution to the dynamical equations (e.g. the density field) can be written asa product of a dimensionless variable, a function of time and some combination ofother invariant parameters and constants of the physical system. Furthermore, ifthe functional dependence of the dimensionless variable is given by r/l(t), r beingthe spatial coordinate and l(t) some time-dependent length scale, the solution isthen said to be self-similar. Considering only one spatial dimension (e.g. sphericalsymmetry) this can be written as,pj(r, t) = cj(t)θj(x), (2.41)where x = r/c0(t). The explicit time dependence is completely contained in thefunctions cj(t), along with all other parameters and constants; these have beensuppressed in the writing of Eq. (2.41). The index j denotes the j-th field variable.If the form of {cj(t)} can be obtained directly from dimensional analysis then thisis referred to as a self-similar solution of the first kind [5]. Alternatively, if thiscannot be obtained from dimensional considerations alone, but requires solving aneigenvalue problem, then it is called incomplete self-similarity, or self-similarity ofthe second kind [5]. This type of self-similar solutions are not relevant in our project,thus they will not be discussed any further.Self-similar solutions of the first kind have been identified in the sphericallysymmetric fluid-dynamic+Newtonian gravity models [56, 60, 81, 98, 108]. For some292.3. Self-similar Ansatzmodels the dimensionless variable is x ∝ r/t. This implies a linear scaling betweenthe space and time variable i.e. the dimensionless functions {θj(x)} remain invariantunder (r, t) → (ηr, ηt), where η ≡ constant. Other models [98, 108], the scalingbetween the space and time variables is non-linear, r ∝ tσ, where σ is a functionof the model parameters. All of these concepts can be formalized by introducingthe machinery of group theory. This is done specifically in the context of Galilean“spacetime” in [14].Even though this project is strictly Newtonian, important parallels exist betweenself-similarity in Newtonian and General Relativistic hydrodynamic systems. Forsimplicity, we only consider self-similarity in spherical symmetry. Self-similarity ofthe first kind in General Relativity implies the existence of a homothetic vectorfield ξ such that the solution (metric (gµν), and matter fields) in the appropriatecoordinates (τ ≡ − ln t, ζ ≡ r/t) is re-scaled under Lie differentiation along ξ =−∂/∂τ , e.g. Lξgµν = 2gµν , since gµν(τ, ζ) = e−2τ g˜µν(ζ), where g˜µν(ζ) is scaleinvariant. Note that the solution’s spatial dependence enters only through the ratioζ = r/t. These solutions were first studied in the pioneering work of Cahill andTaub [11]. The connection between this homothetic solution and its Newtoniancounterpart is discussed in Sec. 2.6. A generalization of self-similarity of the secondkind exists in General Relativity, called kinematic self-similarity [12, 13]. At thepresent time kinematic self-similarity does not seem to occur in critical phenomenain gravitational collapse in GR and for this reason it will not be discussed here.2.3 Self-similar AnsatzA spherically symmetric ideal gas at constant entropy is modeled by Eqs. (2.30)–(2.32). This system admits a self-similar solution. Our interest in this isentropicmodel is to investigate the self-similar solutions that might describe the fluid evo-lution in the ‘epoch’ prior to gravitational collapse. If we define the time variable tsuch that this event (collapse) happens at t = 0, implies we are concerned with the302.3. Self-similar Ansatzinterval ∞ < t < 0. Borrowing from the approach and notation used by Suto andSilk (1988) [98], the transformationt→ −t, v → −v, ρ→ ρ P → P, and M →M, (2.42)leaves the system of Eqs. (2.30)–(2.32) invariant. Following this transformation weconsidered only 0 < t <∞. In some parts of this thesis we have used another shiftof the time coordinate, namely, t → (t0 − t) where t0 denotes the time of collapse.However, unless otherwise specified we use t and the event of interest happens att = 0.The nature of the self-similar solution can be determined by dimensional analysis.We began by assuming a factorization of the time dependence in all of the dynamicalfields according to,r = c0(t)x,ρ(r, t) = c1(t)α(x), v(r, t) = c2(t)u(x),P (r, t) = c3(t)β(x), M(r, t) = c4(t)m(x). (2.43)The fields u(x), α(x), β(x), m(x) are the scale-invariant variables related to v, ρ, Pand M , respectively. Assuming a power-law scaling c0(t) ∝ tn the other coefficientsc1(t), c2(t), c3(t), and c4(t) can be computed following substitution into Eqs. (2.30)–(2.32). With the correct coefficients (2.43) becomes,x =r√κtn, (2.44)v(r, t) =√κtn−1u(x), (2.45)ρ(r, t) =α(x)4πGt2, (2.46)P (r, t) =κt2n−44πG(α(x))Γ, (2.47)M(r, t) =κ3/2t3n−2(3n − 2)Gm(x). (2.48)312.3. Self-similar AnsatzAs expected the coefficients involve the parameters and constants G, K (throughκ = K(4πG)1−Γ) and Γ. The exponent, n = 2 − Γ [98]. If suitable solutions forthe similarity variables u(x), α(x), β(x), m(x) can be found, the general case Γ 6= 1corresponds to self-similar solutions of the first kind with non-linear scaling r ∝ t2−Γ.The case where Γ = 1 (n = 1) corresponds to an isothermal gas, and linear scalingr ∝ t.The autonomous system that results from applying this self-similar ansatz, Eqs. (2.44)–(2.48) into the isentropic gas model, Eqs. (2.30)–(2.32) is,(nx− u)dαdx− αdudx= −2αx(x− u) (2.49)ΓαΓ−2dαdx− (nx− u)dudx= −nx− u3n− 2α− (n− 1)u, (2.50)dmdx= (3n − 2)x2α. (2.51)After some algebraic manipulations we arrive at,dαdx=α(nx− u)2 − ΓαΓ−1[(n− 1)u+ nx− u3n− 2α− 2(x− u)(nx− u)x], (2.52)dudx=1(nx− u)2 − ΓαΓ−1[(nx− u)(n − 1)u+ (nx− u)23n − 2 α− 2Γx− uxαΓ−1],(2.53)together with the algebraic relation,m(x) = αx2(nx− u). (2.54)Thus solving equations (2.52) and (2.53) for the similarity variables α(x) and u(x)with appropriate boundary conditions plus the algebraic condition (2.54) completelydetermines the self-similar flow.322.3. Self-similar Ansatz2.3.1 Boundary ConditionsIn considering solutions to Eq. (2.52), and Eq. (2.53) we restricted our attention tothose which are regular (analytic) at x = 0. This condition is imposed by assuminga Taylor series expansion of α(x) and u(x) in x near the origin. The expansionsare inserted into Eq. (2.52) and Eq. (2.53) in order to determine the expansioncoefficients. Up to O(x4) the asymptotic solutions are given by,u(x) =23x− α1−Γ∗15Γ(α∗ − 23)(n− 23)x3 + . . . , (2.55)α(x) = α∗ − α2−Γ∗6Γ(α∗ − 23)x2 + . . . (2.56)where α∗ ≡ α(x = 0). Here, as was done in [48, 64] we will refer to the parameterQ0 = lnα∗ when specifying the above boundary conditions. Studies involving otherasymptotic solutions which are not regular at at x = 0 exist in the literature.However, they were not considered in this project. These are discussed in [63, 98,105]. Furthermore, we seek solutions to Eq. (2.52) and Eq. (2.53) that exists overan arbitrary range in x i.e., x ∈ (0,∞). In doing so we must be cautious near thesingular point,nx− u =√Γα(Γ−1)/2, (2.57)where the denominator of equations (2.52) and (2.53) vanishes. Known also as thesonic curve, this set of points consists of the all points for which the flow speed v(r, t)relative to the curve r/t equals the local sound speed, cs =√κtn−1(ΓαΓ−1)1/2. Inour treatment we required that α and u be analytic everywhere including across thesonic critical curve. Thus, in the vicinity of the sonic point x = xs+δx it is assumedthat the solution can be Taylor expanded,α(x) = α0 + α1(x− xs) + α2(x− xs)2 + . . . , (2.58)u(x) = u0 + u1(x− xs) + u2(x− xs)2 + . . . . (2.59)332.3. Self-similar AnsatzThe parameter xs labels the sonic point and becomes a new parameter in the problemwhose value is yet to be determined. As a necessary condition for continuity is thatthe numerator of Eq. (2.52) must vanish at the sonic point;0 = (n− 1)u0 + nxs − u03n− 2 α0 − 2(xs − u0)(nxs − u0)xs. (2.60)This, together with Eq. (2.57) allows us to find the coefficients {α0, u0} in terms ofxs. If n = 1 (isothermal gas), the solution for α0 and u0 can be obtained in closedform in terms of xs. However, if n 6= 1 then {α0, u0} can only be computed numer-ically by first determining a numerical value for xs. We can generate conditions forthe higher order coefficients namely, α1, u1, α2, u2, . . . by examining the analyticbehavior of Eqs. (2.53), (2.52) near the sonic point xs. This is accomplished throughthe substitution of the expansions (2.58), (2.59) into Eqs. (2.53), (2.52). Keepingterms up to linear order in x− xs yields,F1 + F2(x− xs) +O((x− xs)2) =0, (2.61)G1 +G2(x− xs) +O((x− xs)2) =0, (2.62)where, F1, G1, F2, G2 are complicated expressions involving the coefficients, α1, u1,α2, u2 and xs. Clearly, from the linear independence (x − xs)k where k ∈ Z, thecoefficients F1, G1, F2, G2 must all vanish, and this gives a set of equations thatallows us to determine the coefficients α1, u1, α2, u2. One important feature whichaids computation of the solution is the fact that each set of conditions involves onlycoefficients of the corresponding and lower orders. The nonlinear set of conditions342.4. Spherically Symmetric Linear Perturbationshas the following general form,F0(α0, u0,Γ, xs) = 0, (2.63)G0(α0, u0,Γ, xs) = 0, (2.64)F1(α0, u0, α1, u1,Γ, xs) = 0, (2.65)G1(α0, u0, α1, u1,Γ, xs) = 0, (2.66)F2(α0, u0, α1, u1, α2, u2,Γ, xs) = 0, (2.67)G2(α0, u0, α1, u1, α2, u2,Γ, xs) = 0. (2.68)Clearly,F0 =nxs − u0 −√Γα(Γ−1)/20 , (2.69)G0 =(n− 1)u0 + nxs − u03n − 2 α0 − 2(xs − u0)(nxs − u0)xs. (2.70)We solved Eqs. (2.63)–(2.68) in sequence (via an iterative method such as Newton’smethod), obtaining the lowest order solution(s) α0, u0 and using them in the higherorder ones. A consistency check can be performed on equations Eqs. (2.63)–(2.68),where setting n = 1 (isothermal case) should yield a system that is soluble in closedform. In this case, the solutions for the coefficients given in terms of xs shouldmatch those obtained in [47, 56, 60, 81, 95]. These expressions are also included inAppx. A.2.4 Spherically Symmetric Linear PerturbationsIncluded with our study of spherically symmetric self-similar solutions to Eq. (2.52),and Eq. (2.53) is an investigation of their linear stability. We considered sphericallysymmetric linear perturbations about the self-similar solutions, α¯(x), u¯(x), m¯(x).352.4. Spherically Symmetric Linear PerturbationsThese can be introduced in the following form,α(r, t) = α¯(x) + εα1(x, t),u(r, t) = u¯(x) + εu1(x, t), (2.71)m(r, t) = m¯(x) + εm1(x, t),for ε << 1. The explicit time dependence is contained in the quantities α1(x, t), u1(x, t)and m1(x, t). Plugging the above ansatz, (2.71)–(2.72) into the dynamical systemgoverning isentropic flows, i.e. Eqs. (2.30)-(2.32) will yield, up to linear order in εthe system,0 =m′1 − (3n − 2)x2α1, (2.72)0 =(nx− u¯)x2α1 − x2α¯u1 −m1 − m˙1t3n− 2 , (2.73)0 =− tm˙1′3n− 2 + (nx− u¯)x2α′1 − (u¯′ − 2)x2α1 − (α¯′u1 + α¯u′1)x2− 2x(α¯u1 + α1u¯), (2.74)0 =− tm˙1(3n− 2)2x2 −α¯u13n− 2 − (nx− u¯)u′1 + (n− 1)u1− (α1α¯′ − α¯α′1)Γα¯Γ−3 + u1u¯′ + tu˙1 +(nx− u¯)α13n− 2 . (2.75)The notation “prime” and/or “dot” (′ and/or )˙ indicate partial differentiation withrespect to x and/or t. The time dependence of the perturbations is assumed tofollow,α1(x, t) = δα(x)eλτ(t) ,u1(x, t) = δu(x)eλτ(t), (2.76)m1(x, t) = δm(x)eλτ(t),362.4. Spherically Symmetric Linear Perturbationswhere,τ = − ln t. (2.77)Using these ansatz in Eqs. (2.72)–(2.75) and dropping the “bar” notation on theself-similar solutions by simply writing α and u together with some algebraic ma-nipulations, we obtained a system of ODEs for the perturbations functions δα(x)and δu(x). These are,δα′ =1(nx− u)2 − ΓαΓ−1{[− ΓαΓ−2α′ + 2u(nx− u)x+(nx− u)α3n− 2− λ+ (u′ − 2− λ)(nx− u)]δα+[(u′ + 3n− λ− 1)α+(nx− u)α′ − 2uαx− α23n− 2− λ]δu}, (2.78)δu′ =1α((nx− u)2 − ΓαΓ−1){[(u′ − λ− 2)ΓαΓ−1 − 2uxΓαΓ−1−(nx− u)ΓαΓ−2 − (nx− u)ΓαΓ−2α′ + (nx− u)2α3n − 2− λ]δα+[2xΓαΓ + ΓαΓ−1α′ + (nx− u)(u′ + n− λ− 1)α−(nx− u)α23n− 2− λ]δu}, (2.79)δm′ =− (3n − 2)x2δα, (2.80)δm =(3n − 2)(nx− u)x23n − 2− λ δα −(3n− 2)x2α3n− 2− λ δu. (2.81)Again, we required that the perturbations also be analytic everywhere. In particular,we impose regularity at x = 0 and at the sonic point (x = xs). Previously wedemanded analyticity of the self similar solutions α and u at x = 0 leading us tothe Eqs. (2.55) and (2.56). Applying these asymptotic solutions to the perturbationequations we find that their leading order behavior near the origin (x << 1) is given372.5. Non-spherical Linear Perturbationsby,δα =3eQ0λδu0, (2.82)δu = xδu0. (2.83)Where δu0 is a free parameter which can be used to normalize the perturbationmode. At the sonic point a necessary condition for analyticity is that the numeratorin Eq. (2.78) vanishes, i.e.(nx− u)[− (nx− u)α′α+2ux+ u′ − 2− λ+ α3n− 2− λ]δα+[u′ + 3n− λ− 1 + (nx− u)α′α− 2ux− α3n− 2− λ]αδu = 0. (2.84)Given the above condition (2.84) it can be easily checked that the numerator inEq. (2.79) also vanishes. Once the expansions (2.58) and (2.59) have been insertedinto Eq. (2.84) we have a boundary condition at the sonic point. Note that sincewe only have one condition the perturbations can only be determined up to a con-stant factor, again, this constant factor is a fundamental freedom that can use it tonormalize the perturbation function. As it will be discussed in more detail later on,the boundary conditions (2.82)–(2.84) can only be satisfied for certain values of theLyapunov exponent, λ. We will further restrict our search for growing perturbationmodes λ > 0 which would tell us the stability of the solutions α and u. Our sta-bility analysis is a generalization of the work of [47] for a polytropic EoS. Again, aconsistency check can be performed in the limit Γ → 1 where their results shouldbe recovered.2.5 Non-spherical Linear PerturbationsClearly, a more complete study concerning the stability of the spherically symmet-ric self-similar solutions mentioned in Sec. 2.3 must include perturbations arising382.5. Non-spherical Linear Perturbationsfrom small departures from spherical symmetry. Previously, several papers havebeen devoted to this topic, notably [38–41]. These investigations ware conductedconsidering general non-spherical perturbations in terms of scalar and vector spher-ical harmonics as outlined in [6]. These are most conveniently treated in sphericalcoordinates r, φ, θ, r is the radial distance from a fixed origin and φ and θ are theazimuthal and polar angles, respectively. Since the background solution is the spher-ically symmetric self-similar solution the formulation of the system of perturbationequations also includes the variables x and τ as defined in Eqs. (2.44) and (2.77).Non-spherical perturbations about the self-similar solutions for the scalar quantitiesin our model e.g. the dimensionless density variable α, generally have the form,α(τ, x, φ, θ) = α¯(x) +∑ℓ,mδαmℓ (x)Ymℓ (φ, θ)eλmℓ τ . (2.85)δαmℓ is a small radial functions, λmℓ are the corresponding growth rates and Ymℓ (φ, θ)are the spherical harmonics. The linear independence and completeness propertiesof the spherical harmonics allows for the representation of any angular dependentfunctions in terms of their sum. Similar deviations from spherical symmetry can beconsidered for the other scalar quantities such as m(τ, x, φ, θ).General perturbations for the vector-valued quantities have a more compli-cated representation. Fortunately, such general representations are possible throughvector-spherical-harmonics. As explained in [6] the general vector perturbation canbe done by taking sums of the vector quantities,~Y mℓ (φ, θ) ≡Y mℓ (φ, θ)rˆ, (2.86)~Ψmℓ (φ, θ) ≡∇Y mℓ (φ, θ), (2.87)~Φmℓ (φ, θ) ≡rˆ ×∇Y mℓ (φ, θ). (2.88)Therefore, perturbations about the self-similar velocity field ~¯u(x) = (u¯(x), 0, 0) have392.5. Non-spherical Linear Perturbationsthe form,~u(τ, x, φ, θ) = ~¯u(x) +∑ℓ,mδuΨ(x)~Ψmℓ (φ, θ)eλτ +∑ℓ,mδuΦ(x)~Φmℓ (φ, θ)eλτ (2.89)Indices {ℓ, m} are implied on the Lyapunov exponent λ and the small quantitiesδuΨ, δuΦ; they have been suppressed to keep the notation from getting too cluttered.Terms which are proportional to ~Ψmℓ are labeled as “polar” perturbations, whereasthose proportional to ~Φmℓ are called “axial” perturbations. This naming conventionis also used in [30, 33].Linear perturbations in the small δ-quantities for the Euler+Poisson systemwith EoS (2.22) decouple into either, purely polar or, purely axial perturbations.This means that to linear order, the equations resulting from considering smallperturbations proportional to Eqs. (2.86)–(2.88) of the spherically symmetric fieldsdo not have a coupled mixtures of ~Ψmℓ and~Φmℓ terms, but an independent sum ofthe two. According to the work done in [39] the lowest-order polar perturbation ofany possible self-similar solutions of the polytropic gas is given by an ℓ = 2 mode.This is called the bar-mode instability. With our choice of initial data (Sec. 2.7)this mode can only be observed in full three dimensional simulations, see [38]. Dueto its limited relevance in our dynamical axisymmetric model, and the fact thatno evidence of unstable polar modes were detected in our simulations (Chap. 5)we did not include further analysis of polar perturbations. Nevertheless, explicitcalculation of all possible unstable polar modes about the spherically-symmetricself-similar solutions for the polytropic gas is left as extension to this thesis. Otherreasons for our omission of these polar modes are discussed in Sec. 2.6. Nevertheless,the polar perturbation equations about the self-similar solutions for the polytropicgas are provided in Appx. B. The axial perturbations about the self-similar velocity402.5. Non-spherical Linear Perturbationsfield ~¯u(x) are given by,urδuθδuφ =u¯(x)δuΦ(x)x sin θ∂∂φY mℓ (θ, φ)eλτ−δuΦ(x)x∂∂θY mℓ (θ, φ)eλτ . (2.90)Note that the dimensionful fluid velocity is obtained through~v(r, φ, θ) =√κtn−1urδuθδuφ (2.91)from Eq. (2.45) and recalling that κ = K(4πG)1−Γ and n = 2 − Γ. Inserting thisexpression for the velocity into the continuity Eq. (2.5) and keeping only the linear-order terms in the perturbations it is easy to show thatδαmℓ (x) = 0. (2.92)Recall that δαmℓ (x) represent non-spherical perturbations about the spherically sym-metric dimensionless density variable α¯(x). From this and Poisson’s equation followsthatδϕmℓ (x) = 0 (2.93)for the perturbations of the Newtonian potential. We have already mentioned thatthe energy Eq. (2.7) becomes redundant for a barotropic EoS, therefore, the axialperturbations in this case are completely governed by the momentum equation (2.6).412.5. Non-spherical Linear PerturbationsWhen this is written in spherical coordinates we get,∂ρvr∂t+1r2∂∂r(r2ρv2r ) +1r sin θ∂∂θ(sin θρvrvθ) +1r sin θ∂∂φ(ρvrvφ)−ρ(v2θ + v2φ)r= −∂P∂r− ρ∂ϕ∂r, (2.94)∂ρvθ∂t+1r2∂∂r(r2ρvrvθ) +1r sin θ∂∂θ(sin θρv2θ) +1r sin θ∂∂φ(ρvθvφ)+ρvrvθr− cot θρv2φr= −1r∂P∂θ− ρr∂ϕ∂θ, (2.95)∂ρvφ∂t+1r2∂∂r(r2ρvrvφ) +1r sin θ∂∂θ(sin θρvθvφ) +1r sin θ∂∂φ(ρv2φ)+ρvrvφr+ cot θρvθvφr= − 1r sin θ∂P∂φ− ρr sin θ∂ϕ∂φ. (2.96)The linear perturbations are obtained by inserting equations (2.90), (2.92), and(2.93) into the above system; we find that,∂∂t(ρδvθ) +1r2∂∂r(r2ρvrδvθ) =0, (2.97)∂∂t(ρδvφ) +1r2∂∂r(r2ρvrδvφ) =− ρvrδvφr. (2.98)In term of the similarity variables this becomes,(λ+ 3− n)αδuΦx+ nxddx(αδuΦx)− 1x2ddx(xαuδuΦ) =αuδuΦx2(2.99)Combining this equation with relation (2.49) and some algebraic manipulations weget the following simple equation for the axial perturbations,(λ− 2n + 1)δuΦ + (nx− u)dδuΦdx= 0. (2.100)422.5. Non-spherical Linear PerturbationsUsing the fact that near the origin u ≈ 2x/3 we can solve for λ,λ = 3− 2Γ +(Γ− 43)limx→0xdδuΦdxδuΦ=13+(Γ− 43) limx→0xdδuΦdxδuΦ− 2 . (2.101)We only consider perturbations which are regular at the origin, in particular vφ ∼ rnear r = 0. This implies that δuΦ(x) can be Taylor expanded in powers or x thus,δuΦ(x) =∞∑ν=0xℓ+1aℓνx2ν where, x≪ 1, ℓ = 1, 2, . . . (2.102)and aℓν are scalar coefficients. The above expression Eq. (2.102) ensures regularityof the vφ near the origin. Inserting expression (2.102) into Eq. (2.101) we get,λℓ,ν =13+(Γ− 43)(ℓ− 1 + 2ν) . (2.103)In the case where the overall perturbation grows linearly in x the fastest growingmode corresponds to (ℓ = 1, ν = 0) with growth rate λ = 1/3, provided Γ < 4/3.Therefore, all regular self-similar solutions of the polytropic gas have an (ℓ = 1)unstable axial mode. For the an isothermal gas (Γ = 1) only the ℓ = 1 mode isunstable. If we consider Γ ∈ [1, 4/3) higher ℓ-values can be unstable, however, theirgrowth rates are always less than 1/3, see Fig. 2.1. For example, in the case where1 < Γ / 1.17 the polytropic gas contains two unstable axial modes, namely ℓ = 1 anℓ = 2. However, if we restrict our attention to axial symmetry, along with equatorialsymmetry of the initial data (Sec. 2.7), the ℓ = 2 mode is suppressed. The ℓ = 1axial perturbation mode is often called the spin-up mode, [38, 41]. A formal solution432.6. Correspondence with General Relativityfor the spin-up mode can be obtained by direct integration of Eq. (2.100)δuΦ(x) = δuΦ0 exp(−∫λ− 2n+ 1nx− u dx). (2.104)δuΦ0 is a constant of integration. From this expression it is clear that all regularself-similar solutions are unstable against the spin-up mode.Figure 2.1: This plot illustrates the spectrum of axial perturbation modes and itsdependence on Γ. For the isothermal gas (Γ = 1) only the ℓ = 1 perturbation modeis unstable with growth rate λ1,0 = 1/3. The ℓ = 1 mode has the same growthrate (λ1,0 = 1/3) independent of the value of Γ, thus all self-similar solutions of thepolytropic gas are unstable against this axial perturbation. Known as the spin-upmode, it is the fastest growing axial mode for 1 ≤ Γ < 4/3. Values of Γ > 1 lead tohigher unstable ℓ modes, from these, the ℓ = 2 mode is first to “turns on” (becomesunstable) at Γ > 1 but its growth rate is less than that of the spin-up mode for thevalues of Γ that were considered. Given, axial symmetry and the symmetry of theinitial data (Sec. 2.7), this ℓ = 2 mode is not observable. The ℓ = 3 mode becomesunstable at Γ ' Correspondence with General RelativityThis investigation was initially motivated by previous studies involving sphericallysymmetric self-similar solutions of a perfect fluid in General relativity (GR) [10, 11,442.6. Correspondence with General Relativity22, 45–47, 73, 74, 77, 96]. All these studies considered a barotropic EoS (1.3), whichdescribes an isothermal gas. However, in GR the interpretation typically given toEq. (1.3) is that of an ultra-relativistic perfect fluid. This can be derived from theideal gas law P = kρ0ǫ, where ǫ is the specific energy density, and ρ0 is the densitymeasured in the rest frame of the fluid element. The product ρ0c2 (where c is thespeed of light) is the rest-mass energy density of the fluid. The total energy densityof the fluid element is given by ρc2 = ρ0c2(ǫ/c2 + 1) (internal plus rest-mass energydensity). In the ultrarelativistic limit, ρ0ǫ≫ ρ0c2, so ρ ≈ ρ0ǫ and we get Eq. (1.3).To distinguish the ultrarelativistic EoS from the isothermal gas, Eq. (2.23) we usea lower-case k in the former case.All the calculations we have performed were strictly Newtonian. The relevantresults from GR will simply be reported without repeating any of the calculations.The interested reader will be provided with the appropriate references for completedetails and discussions. Ori and Piran calculated general relativistic self-similarsolutions of a perfect fluid with EoS (1.3) in spherical symmetry [77]. They showedthat by taking a particular limit of k → 0 (from the EoS) the self-similar solutions ofthe general relativistic autonomous system reduced to the Newtonian isothermal-gasself-similar solutions. Therefore, the limit k → 0 in EoS (1.3) can be considered theNewtonian limit of the general relativistic perfect fluid. Through this transition theNewtonian solutions inherit the properties of their general relativistic counterpartsas discovered by [44, 47, 77, 96]. The Newtonian solutions are briefly discussed insection Regular Newtonian Self-similar SolutionsRegular Newtonian self-similar solutions for an isothermal gas were first studied in[56, 60, 81]. One of these solutions simultaneously discovered by Larson and Penston(called the Larson-Penston (LP) solution) is regular in the similarity variable butsingular at the origin in the dimensionful variables such as the fluid density ρ(r, t).452.6. Correspondence with General RelativityThe blow-up4 occurs at a finite time, and has finite mass. The solution can alsobe studied after the formation of the singularity [64, 77]. After collapse, the corecontinues to grow by accreting mass thus this solution is analogous to black holesolutions in GR, see [77]. Hunter (1977) [56] reported an infinite branch of self-similar solutions to the Newtonian-isothermal-gas-system different from the Larson-Penston solution. These solutions, again, contained a singularity in the regularcoordinates (r, t) at a finite future time t0. However, in this case the “collapsing core”shirks to zero as t → t0 and thus has zero mass. After collapse, the fluid expandsoutward, and the central density tends to zero as t →∞ [64, 77]. When expressedin terms of the similarity variable x ∝ r/t all these solutions are finite and regular atthe center (x = 0). Their properties at the origin (x = 0) can be used to order thesolutions into a hierarchy. In the literature the solutions are labeled as Hunter-A,Hunter-B, Hunter-C, etc, with Hunter-A (HA) having the lowest amplitude in thedimensionless density similarity variable, α(x) in our notation. A sketch of thesesolutions for t0 − t > 0 is provided in Fig. 2.2, the hierarchical structure formedby the amplitude of α(x) can be appreciated in Fig. 2.2(a). The Hunter-solutionsfor the self-similar fluid velocity u(x) as seen in Fig. 2.2(b) all contain oscillationswhich also follow a hierarchical paradigm. The Larson-Penston does not containoscillations in the profile of u(x), thus it is a pure collapse solution. A third type ofself-similar-solution reported by Maeda et. al. (2001) is the so-called homogeneoussolution. This solution is particularly simple and it can be obtained in closed form,see Table 4.1. It represents either an expanding or a collapsing solution. The lattercorresponds to a “big crunch” singularity; the collapse happens everywhere at t0.This solution may have some relevance in Newtonian cosmology [82]. This solutionis also plotted in Fig. 2.2.Spherical linear perturbations of the Newtonian isothermal gas self-similar so-lutions were studied in [64]. The calculations reveal that the LP solution has no4Mathematicians often call the singular behavior in the governing PDEs as a “blow-up”.462.6. Correspondence with General Relativity(a) (b)Figure 2.2: Sketch of the regular self-similar solutions for the isothermal gas (Γ = 1).Plots of the dimensionless density α(x) are provided in panel (a). Aside from thehomogeneous solution, the Larson-Penston solution has the lowest amplitude withα(0) ≈ 1.6658 (a). Members of the Hunter-family have increasingly higher ampli-tudes arranged in hierarchical order. As seen in panel (b), the Hunter-solutionscontains oscillations in the self-similar velocity profile (u(x)). The number of theseoscillations also follows the same hierarchical structure, with the Larson-Penstonand homogeneous solutions having 0 oscillations, Hunter-A 1, Hunter-B 2, Hunter-C 3, and so on. The absence of any oscillation in the profile of the Larson-Penstonsolution means it describes homologous collapse of the core; the homogeneous solu-tions describes simultaneous collapse. All of these plots correspond to the regimewhere t < t0, (t0 being the time of collapse) thus they illustrate the behavior of thefields prior to the time of singularity formation. These solutions are also calculatedand discussed in [64].unstable spherically symmetric linear perturbations modes. In contrast, the Huntersolutions contain growing (unstable) perturbation modes. The Hunter-A solutionhas a single unstable mode. The Hunter-B solution has two unstable modes, Hunter-C has three and so on. All the Hunter solutions investigated in [64] follow thispattern. The homogeneous solution was found to be unstable against fluctuationsemanating from the sonic point. The fluctuations grow into a shock thus destroyingthe regularity of the solution. This is known as “kink”-mode instability.472.6. Correspondence with General Relativity2.6.2 Regular GR Self-similar SolutionsThe pioneering work of Cahill and Taub (1971) [11] on spherical similarity solutionsof Einstein’s equations for a perfect fluid promoted interest in the investigation ofgravitational collapse. This is due to the fact that self-similar solutions provide asimpler avenue for the analysis of the gravitational collapse process. An exampleof this was the work Ori and Piran (1990) [77]; they studied and classified themultiple similarity solutions of the GR+perfect fluid system. A summary of theirclassification and description is given below.(1) Black-hole solutions: Singular solutions with finite mass. The singularity(Black hole) forms at a finite proper time. The collapse happens homologously.After black hole formation, the mass continues to grow by accretion. Thesesolutions are the general relativistic versions of the Larson-Penston blow-upsolution. Some authors refer to this as the GRLP (general relativistic Larson-Penston) others called it the Ori-Piran solution(s) [46].(2) Repulsive solutions: These solutions contain oscillations in their velocity pro-files. Again, the singularity forms at a finite proper time t0 but immediatelyvanishes. The central density tends to zero for t ≫ t0, the gas expands out-wardly. These solutions are analogous to the Newtonian Hunter solutions.The first member of this family is the general relativistic Hunter-A solution(GRHA). This is also known as the Evans-Coleman5solution [22, 46].(3) Asymptotically Friedmann solutions This solutions describe a “big crunch”singularity. The entire spacetime collapses at t0. This is the general relativisticversion of Newtonian homogeneous solutions previously described.The spectrum of similarity solutions which are analytic everywhere is discrete justas its Newtonian counterpart.5In most cases, the Evans-Coleman solutions refers to the 1-mode unstable self-similar solutionfor the radiation gas k = 1/3.482.6. Correspondence with General RelativityLinear stability analysis of the Ori-Piran (GRLP) solution reveal an absenceof spherically symmetric growing mode perturbations [45, 59]. The Evans-Coleman(GRHA) solution was shown to contain a single such perturbation mode [22, 66, 99].The specific value of the Lyapunov exponent (growing mode) has been found todepend on the parameter k of the EoS (1.3). This single analytic mode structureof the Evans-Coleman solution persists over the entire range k ∈ [0, 1] [73, 74].Solving the general relativistic spherical-linear-perturbation equations reveal thatas k → 0 the value of the unstable mode (Lyapunov exponent) converges to theknown value for the Hunter-A solution [45]. Further evidence of this convergencecan also be found in dynamical simulations of critical collapse of an ultrarelativisticperfect fluid in the limit (k → 0) [96]. We will talk more about this result in ourdiscussion of critical phenomena, Chap. 4.Table 2.1: Newtonian-GR similarity solutions.Newtonian General RelativityLarson-Penston Ori-Piran (GRLP)Hunter-A Evans-Coleman (GRHA)Hunter-B GRHBHomogeneous solutionFriedmann Solution†α(x) = 2/3, u(x) = −2x/3Static Singular SolutionStatic Singular Solution†α(x) = 2/x2, u(x) = 0† Their explicit form can be found in Ref. [77]2.6.3 Non-spherical PerturbationsGundlach (2002) calculated non-spherical linear perturbations about the similaritysolutions for the perfect fluid in GR with EoS (1.3) [32, 33]. Similar to the Newtoniansystem the linear-order axial and polar perturbations are decoupled. The growth492.6. Correspondence with General Relativityrate spectrum of all axial perturbations can be computed in closed form, this isgiven by,λnℓ (k) =2(1 − 3k)− (1 + 3k)(ℓ+ 2n)3(1 + k), (2.105)where k is the EoS proportionality constant in Eq. (1.3). The parameter n ∈{0, 1, 2, . . . }, it denotes a particular solution of the radial equations. Naturally,Eq. (2.105) is the General Relativistic version of Eq. (2.103). Clearly, setting, n = 0in Eq. (2.105) yields the least negative, slowest decaying perturbation for any k andℓ. For k ∈ [0, 1] with n = 0, all ℓ ≥ 2 are negative and thus decaying, but forℓ = 1 within 0 < k ≤ 1/9 there is exactly one growing (λ01 > 0) axial perturbation,namely,λ01 =1− 9k3(1 + k). (2.106)Similar to the Newtonian fluid the axial mode is common to all similarity solutions.In particular, the Evans-Coleman solution, which we know has an unstable sphericalmode is also unstable against the axial mode. Gundlach showed that in the range0 < k ≤ 1/9 the Evans-Coleman solution is stable against all polar perturbations.This was shown by numerical analysis of the perturbations equations [33]. Wemust point out that the k-range over which no unstable polar perturbations werefound includes the Newtonian limit (k → 0), hence, the Newtonian equivalent tothe Evans-Coleman solution i.e. the Hunter-A solution is predicted to be absentof any unstable polar modes. Note that in the Newtonian limit the growth rate ofthe unstable axial mode λ01 = 1/3 which is the known spin-up mode discussed inSec. 2.5. The stability structure of the Evan-Coleman solution is exactly analogousto the Newtonian Hunter-A solution, at least for the range where k ∈ [0, 1/9]. Thisis a key property of this unstable solution that motivated this project. Both theEvans-Coleman solution (also known as GRHA) and the Hunter-A solutions havethe same mode structure. Therefore, we expect that their role as critical collapsesolutions result in very similar type-II critical phenomena. This has already been502.7. Initial Datashown in spherical symmetry. Our goal here is to test this in axisymmetry with theaddition of rotation.2.7 Initial DataA large portion of our research involve investigating the role of angular momentumin gravitational collapse beginning with generic initial data in axial symmetry. Inthe case of initial infinitesimal rotation, it is expected that the angular momentumof the final compact object be related to the initial axial perturbations [32]. Theseaxial perturbations in general can be described by a vector-value (three functions)function of the spatial coordinates i.e. ~g(r, θ, φ). Following Gundlach’s suggestionswe can explore this 3∞-dimensional parameter space by means of vector-value fac-tor ~q [32]. The initially rotating data, or more specifically the vφ (in cylindricalcoordinates (s, φ, z)) velocity field is specified by the product ~qg0(r, θ, φ). The restof the initial data is constructed such that ~q = 0 reduces to the spherical symmetry.In addition to the three parameters (components) of ~q, there is an extra controlparameter labeled as p which controls other spherically symmetric aspect of theinitial conditions e.g. the amplitude of the initial density profile. Thus, the initialcondition is controlled by a 4-family of initial data.The quantities of interest during the evolution of the fluid are the mass of thecompact object (M) and its specific angular momentum (~a). Symmetry requiresthat these quantities be related to the parameters p and ~q by way of,M(p,−~q) =M(p, ~q), ~a(p,−~q) = −~a(p, ~q). (2.107)Various forms for the profiles g0(r, θ, φ) which physically specify differential rotationwere tested, all with the reflection symmetry about the equatorial plane.512.8. Critical Phenomena2.8 Critical PhenomenaNumerical calculations of the evolution of a spherically symmetric real scalar fieldin General Relativity at the threshold of black hole formation led to the discov-ery of critical phenomena in gravitational collapse [15]. The numerical experimentsrevealed that the final fate of this data belonged to one of two categories. One,characterized by complete collapse of the matter-energy, i.e. black hole formation,i.e. the supercritical regime. Two, the subcritical regime is characterized by com-plete dispersal of the matter-energy out to infinity leaving behind flat, Minkowskispace. A bisection search was carried out between these two regimes to locate thethreshold, critical value of a control parameter p, this parameter is labeled p⋆. Themass of the black hole that formed during supercritical evolutions was shown todepend of on the “distance” from the critical parameter p⋆ according the followingpower-law,MBH ∝ |p− p⋆|γ . (2.108)This is the now famous scaling-law for the black hole mass. The black hole massesresulting from evolutions across the critical parameter p⋆ appear to “turns-on” atarbitrarily small values. The exponent γ is independent of the particular choice of1-parameter family of initial data, thus it is universal. Furthermore, the criticalsolution (scalar field and metric functions) display discrete self-similarity, see [15]for more details. This solution is also described in [76, 85]. The self-similarity ofthe critical solution, power-law scaling and universality are reminiscent of 2th orderphase-transition in statistical mechanics. This is the reason for the label type-IIcritical phenomena.Type-II critical phenomena was also identified in critical collapse simulations of aradiation fluid by Evans and Coleman [22]. The critical solution is continuously self-similar–it is the single-unstable-mode Evans-Coleman solution previously known.The reciprocal of the scaling exponent γ is the Lyapunov exponent (growth rate)522.8. Critical Phenomenaof the unstable mode. The work Koike et. al. [99] and Maison [66] based on therenormalization group provided a coherent framework which associates the salientfeatures of critical phenomena i.e. universality and power-law scaling with thesingle-mode instability property of the critical solution.2.8.1 Type II Critical Behavior and Scaling LawsThe formalism provided Koike et. al. applies to any autonomous system of equationswhose solution under linear spherically symmetric perturbations has a single growingmode. Here, we would like to summarize and generalize their analysis for the case ofa spherically symmetric self-similar solution containing two unstable modes. One ofthese modes is spherically symmetric, while the other is a non-spherical axial mode.We borrow from Gundlach’s analysis [32], where he describes an intermediate linearregime governed by two unstable modes. Finally, we will apply this formalism to theNewtonian fluid model to derive scaling laws for the collapsed mass and maximumdensity.Let Z(τ, x, θ, φ) represent a solution to the dynamical system (PDEs) in sphericalcoordinates. The variables,x =rl(t),τ = − ln t, (2.109)where l(t) is some time dependent length scale. If the solution is solely a functionx i.e. Z(x) clearly then, the solution is spherically symmetric and self-similar. Fur-thermore, if Z(x) has a single spherically symmetric unstable mode we label thissolution by Z⋆(x). In addition, we assume that Z⋆(x) solution also has an axialnon-spherical unstable mode. Our aim is to construct the 4-parameter initial data,controlled by p and ~q as described in Sec. 2.7 that goes through the “intermedi-ate linear regime”. This is the region in phase space parametrized by (p, ~q) that532.8. Critical Phenomenayields solutions that are well approximated by Z⋆(x) plus the two unstable modes(spherical and axial-non-spherical). This is represented by,Z(τ, x, θ, φ; p, ~q) ≃ Z⋆(x) +A(p, ~q)Z0(x)eλ0τ + ~B(p, ~q) · ~Z1(x, θ, φ)eλ1τ . (2.110)Z0(x) and ~Z1(x, θ, φ) are the spherical and non-spherical mode functions respec-tively. Their corresponding growth rates are λ0 and λ1. We have explicitly includedthe dependence on the initial data parameters p and ~q through the constant A(p, ~q),and constant vector ~B(p, ~q). Recall that we also require that ~q = 0 reduces to spher-ical symmetry. This implies ~B(p, 0) = 0. The solution then only has one unstablemode. In order to properly inherit the spherically symmetric results discussed in[10, 99] we must require the existence of initial critical value parameter p = p⋆ suchthat A(p⋆, 0) = 0. We refer to (p = p⋆, ~q = 0) as the critical point since both A and~B vanish there.For simplicity, we assumed the functions A(p, ~q), and ~B(p, ~q) are analytic at thecritical point as in [10, 32]. Their leading order behavior around the critical pointis given by,A(p, ~q) ≃ ∂A(p⋆, 0)∂p(p− p⋆) + second order, (2.111)~B(p, ~q) ≃ ∂~B(p⋆, 0)∂~q~q + second order. (2.112)By definition of a critical solution A(p⋆, 0) = 0 = ~B(p⋆, 0). With this more explicitform of the constants, the intermediate linear regime is now given by,Z(τ, x, θ, φ; p, ~q) ≃ Z⋆(x) + C0(p − p⋆)Z0(x)eλ0τ + ~qC1 · ~Z1(x, θ, φ)eλ1τ . (2.113)where the scalar C0 ≡ ∂pA(p⋆, 0) and the 3× 3 matrix C1 ≡ ∂~q ~B(p⋆, 0).At some intermediate time (τ∗) our solution is well approximated by (2.113).542.8. Critical PhenomenaWe can relabel the small deviations from Z⋆(x) by introducing,ǫ ≡ C0(p− p⋆)eλ0τ∗ , (2.114)~δ ≡ ~qC1eλ1τ∗ , (2.115)where ǫ and ~δ are now the small fixed quantities. Solving for τ∗ in (2.114) andreplacing it in (2.115) yields,~δ =~qC1(C0ǫ−1)λ1/λ0(p− p⋆)λ1/λ0. (2.116)Following Gundlach’s argument [32] we takeZ(τ∗, x, θ, φ; p, ~q) ≃ Z⋆(x) + ǫZ0(x) + ~δ · ~Z1(x, θ, φ). (2.117)as Cauchy data. Later on, when the perturbation(s) have grown sufficiently large,ǫ ∼ O(1), the solution will either disperse or collapse. Nevertheless, the only scalein our system is set by e−τ∗ [99], which from Eq. (2.114) we get that, e−τ∗ =(C0ǫ−1|p− p⋆|)1/λ0 .At this point, we would like to apply this formalism to the our Newtonian the axi-symmetric polytropic gas, Eqs. (2.34)–(2.39). We assume that dynamical solutionto this Newtonian model with the “appropriate initial conditions” controlled byp and ~q is able to approach a 2-mode unstable spherically symmetric self-similarsolution. So, the solution at some intermediate time is described by the linearregime given by Eq. (2.117). The total mass enclosed inside radius r, at time t forany similarity solution is given by Eq. (2.48). Its dimensions are given by the explicittime dependent part i.e. M(r, t) ∝ t3n−2. Therefore, the mass of the compact objectthat may form following the solution’s departure from the linear regime scales ast3n−2∗ , where t∗ = e−τ∗ , by definition. The remaining dependence can only be on the552.8. Critical Phenomenadimensionless parameter ~δ and the sign of (p − p⋆). We then have,M(p, ~q) = (C¯0|p − p⋆|)(4−3Γ)/λ0F (~δ) (2.118)We have set C0ǫ−1 ≡ C¯0, and substituted n = 2−Γ. The function F (~δ), along withλ0 are universal as stipulated in [32]. Similarly, the scale of the specific angular mo-mentum, denoted by ~a is set by t2n−1∗ . This is determined from dimensional analysissince the units of ~a are related the units of M(r, t) (mass) and ~v(r, t) (velocity)according to[~a] =[G][M ][~v], (2.119)where [ ] means “units of” and G is the Gravitational constant. Therefore,~a(p, ~q) = (C¯0|p − p⋆|)(3−2Γ)/λ0 ~G(~δ). (2.120)The vector function ~G(~δ) is also universal. The scaling laws for an isothermal gasare obtained by setting Γ = 1, the solutions have self-similarity of the first kind. Inthis case Eq. (2.118) and Eq. (2.120) are formally identical to equations (17) and(19) derived in [32].In the limit ~q → 0 we should recover the spherically symmetric results. FromEq. (2.107) we conclude that F (~δ) must be even in ~δ, whereas ~G(~δ) must be odd.Thus, in the limit ~q → 0 the scaling law for the mass is,M =[C¯0(p − p⋆)](4−3Γ)/λ0F (0)=[C¯0(p − p⋆)](4−3Γ)/λ0 , (2.121)where we have recovered Choptuik’s scaling law for the mass of the compact ob-ject as a function of the collapse threshold [15]. The self-similar critical solution(Z⋆(x ∝ r/t2−Γ)) in this case yields a scaling exponent that depends explicitly onthe adiabatic index Γ, as well as λ0 due to the non-linear scaling r ∝ t2−Γ. From562.8. Critical PhenomenaEq. (2.107) we conclude that the leading order behavior of ~G(~δ) in the limit ~q → 0is proportional to ~δ. Then, the scaling law for ~a follows,~a =[C¯0(p− p⋆)](3−2Γ)/λ0~δ=C1~q[C¯0(p − p⋆)](3−2Γ−λ1)/λ0 . (2.122)Since λ1 = 1/3 and ignoring all the family-dependent constants we can write,M ∝ (p − p⋆)(4−3Γ)/λ0 as ~q → 0, (2.123)~a ∝ ~q(p − p⋆)(2/3)(4−3Γ)/λ0 as ~q → 0. (2.124)In axial symmetry, the vector quantities can be treated as scalars once the axis ofsymmetry is aligned with the z-axis. With these conditions and using cylindricalcoordinates (s, φ, z) we have,~a = azˆ, ~δ = δzˆ, ~G(~δ) = G(δ)zˆ. (2.125)zˆ is the unit-vector along the z-axis. The above discussion is a summary of [32].Notice that the scaling laws for M and ~a highlight the importance of Γ = 4/3, (n =2/3), at this value, the autonomous system given by Eq. (2.52) and Eq. (2.53)becomes singular and no self-similar solutions exists.In a similar way we can generate the scaling law for the maximum central densityattained during subcritical evolutions ρ(0, t) ≡ ρmax. Since ρmax ∝ t−2∗ we find thatρmax ∝ |p − p⋆|−2/λ0 as ~q → 0. (2.126)If critical collapse initial data of the ideal gas model in Newtonian gravity reaches aHunter-A-type similarity solutions then we expect the characteristic variables suchas the mass and specific angular momentum obey the scaling laws just described.They represent the key predictions to be tested against dynamical evolutions of572.8. Critical Phenomenacritical collapse. These tests are presented in Chap. 4 and Chap. Type I Critical BehaviorA possible result for the evolution of our fluid+gravity model endowed with criticalinitial data is a static, star-like solution, Z⋆(r) instead of a Hunter-A-type of scale-invariant solution. Similar to type-II critical behavior, this solution is only observedfollowing the fine tuning of a 1-parameter family of initial data to the brink ofgravitational collapse. Again, this solution is by definition a critical solution. Here,were discuss the general properties of this solution. The formalism of Koike et.al. [99] applies equally to the static critical solution. The only difference is thatthe static case requires the use of regular coordinates (r, t) as oppose to zoomingcoordinates (ln(x), τ). The formalism is nearly identical to the self-similar case, sowe will not repeat it here. This is done in detail in [72, 76]. In this case, the massof static solution does not change; proximity to the threshold value p⋆ does not leadto smaller black holes but instead, a mass-gap exists set by the length scale of thestatic solution. Nevertheless, a scaling law for the lifetime T0 of the critical solutioncan be deduced. As before, this is a result of the solution being funneled throughthe linear regime. This is given by,T0 ∝ − 1λ0ln |p− p⋆|. (2.127)Here, λ0 is the growth rate of the spherically symmetric unstable mode. This kindof behavior is called type-I critical phenomenaAt the current time there is no formalism that introduces slow rotation of thestatic critical solution as non-spherical linear perturbations in a way that is anal-ogous to what is done for the similarity solutions, Sec. 2.8.1. Heuristically, onepossibility is to incorporate the effect of rotation on the critical solution by adding582.8. Critical Phenomenaa non-spherical function, i.e.Z(t, r, θ, φ; p, ~q) ≃ Z⋆(r) +A(p, ~q)Z0(r)eλ0t + f(t, r, θ, φ; p, ~q). (2.128)Clearly, we should require that f(t, r, θ, φ; p, ~q = 0) ≡ 0. Then, at ~q = 0 an expansionof A(p, ~q) near p = p⋆ yields the scaling law (2.127), recovering the sphericallysymmetric result. We could write the function f(t, r, θ, φ; p, ~q) in terms of the linearnon-spherical modes of the critical solution Z⋆(r). However, as discussed at the endof Chap. 4 we do not have access to the solution Z⋆(r) nor to its linear sphericalperturbations making it impossible to check Eq. (2.128) directly. This is a resultof not having enough equations to solve the hydrostatic system, see Sec. 4.2.4.Therefore, we will not say much more about the linear regime in the case of type-Icritical phenomena. The only available tool for checking whether the quasi-staticsolutions encountered in this project exhibit type-I critical behavior is Eq. (2.127).Due to the lack of a theoretical basis for treating rotation in critical collapse of ourNewtonian model the topic is only marginally explored in this project. It is left forfuture work.59Chapter 3Numerical MethodsThis chapter contains a brief presentation of the methods and techniques used inthe numerical analysis of the equations of motion introduced in Chap. 2. The mainobjective of these methods is to provide approximate solutions to the time-dependentnonlinear system of partial differential equations (PDEs). Due to the nature ofthis project a combination of numerical techniques were used to simultaneouslyprovide solutions to the gravitational field determined by the Newtonian potentialϕ(~x, t), and the dynamical fluid variables ρ(~x, t), ~v(~x, t) and P (~x, t) (density, velocityand pressure). Considerably more effort was allocated to finding solutions to thehydrodynamic component which is given by a system of conservation laws. Thenonlinear nature of these system is known to give rise to discontinuities (shocks)even from smooth initial conditions.3.1 Finite DifferencesArguably the most common of all numerical techniques used in the production ofapproximate solutions to PDEs are finite differencing algorithms. These approacheshinge on generating a finite-difference representation of a system of PDEs. Thisis done by Taylor expansions of the continuum solution at nearby points in orderto approximate the derivatives in the equations through finite differences. The re-sult of the procedure is a system of equations that relates the solution at variouspoints on a discrete domain. The finite-difference equations comprise an algebraicsystem of equations whose solutions represent approximations to the continuum sys-603.1. Finite Differencestem (PDEs). The finite difference approximation can in general be made arbitrarythough some control parameter, e.g. h, the grid spacing, at the cost of a greaternumber of algebraic manipulations. We provide a brief introduction to the method,a more complete discussion can be found in [16, 35, 71, 83]. Here, we borrow thenotation used in [16].A system of PDEs can be written in the following abstract notation,Lu = f, (3.1)the solution u, can be a scalar or vector-valued function. The function f is sometimescalled the “source”, and L is a type of differential operator. For what follows it isassumed that both u and f are continuous and smooth. In general u = u(~y, t), andf = f(u; ~y, t) where ~y denotes the position vector in some coordinate system and t isthe time. For simplicity we assume that u is a scalar function and is dependent upontime and only one spatial variable i.e. (y, t). Similarly, we assume that f = f(u, y, t)and consequently it is also a scalar function. A commonly given example is the linearwave equation,(∂tt − ∂yy)u(y, t) = 0, (3.2)clearly in this case, L = ∂tt − ∂yy, and f(u, y, t) = 0.As previously stated, the process of approximating the continuum solution u(y, t)follows from the construction of a discrete system analogous to (3.1). First, we con-sider a discrete domain, a subset of points selected from the continuum. To simplifythe discussion a set of N uniformly spaced points is selected {y1, y2, . . . , yj, . . . , yN},where yj+1−yj = h. Thus, h is the fundamental control parameter of the discretiza-tion. Again, adopting the notation used in [16], the discrete system has the form,Lhuh = fh. (3.3)Where now, uh represents the solution of the discrete system, the other quantities613.1. Finite Differenceswith superscript h are the discrete version of the quantities in (3.1). Equation (3.3)is called the finite difference approximation (FDA) of the PDE, Eq. (3.1). If ourconstruction of the difference operator, Lh, is correct then the discrete solutionuh should approximate u. As we said before, the degree of the approximation isdetermined by the grid spacing parameter h. By design, the FDA (3.3) shouldreduce to Eq. (3.1) as h→ 0.Our goal is to find the discrete solution uh by inverting the difference operatorLh. For nonlinear algebraic equations a close-form discrete solution is in generalunattainable. Nevertheless, an approximate solution to the algebraic system, u˜h,can still be provided through an iterative process such as Newton’s method. Thisleads to the introduction of a residual function, defined asrh = Lhu˜h − fh. (3.4)It measures the degree to which the difference equation (3.3) is satisfied.Now, we address the method for generating discrete representation(s) of thedifferential operator L. This can be done by way of multiple Taylor series expansionsof the continuum system. This allows us to represent any combination of derivatives,or higher derivatives as combinations of finite differences between neighboring points.The operator can be approximated to the desired order in the power of the distancebetween points, i.e. O(hp). The number of terms we choose to keep in the Taylorseries expansion of L determines the order of our finite difference approximation.Naturally, this truncation of the series introduces another source of error. A Taylorseries expansion of the original PDE has the form,τh = Lhu− fh. (3.5)The quantity τh is called the truncation error, it contains higher order terms notincluded in the definition of Lh. Note that u is the solution to the continuum system.623.1. Finite DifferencesThe solution to the difference equation uh is said to converge if it satisfies, uh → uas h→ 0. For consistency we must add,limh→0τh = 0, (3.6)This ensures that the finite difference approximation reduces to the continuum equa-tion (PDE).Now that the basics concepts have been laid out, we discuss how we quantify thesolution error involved in the numerical approximation. Naturally, this is definedas,eh = u− uh. (3.7)Similar to the other error quantities, rh and τh we would like to express eh as afunction of h. More specifically, we want to express eh as a sum of linear powersof h. This will allow us to specify the order of approximation of the numericalalgorithm. However, the definition of the solution error involves the exact solutionu, which is what we want to determine in the first place. Fortunately, the differenceu − uh ≡ eh should be a smooth function and thus expressible as an infinite seriesin powers of h. Then,eh = u− uh =∞∑l=1hlel. (3.8)This is known as the Richardson’s expansion [87]. The coefficients el = el(y, t).Let us assume that the truncation error for some discretization scheme is O(h2).ThereforeLhu− fh = τh = τ2(x, t)h2 +O(h3).Furthermore, let us assume that we could solve the FDA equation exactly i.e. Lhuh−633.1. Finite Differencesfh = 0. Then,Lh(u− uh) =∞∑l=1hl(Lhel) = τ2(x, t)h2 +O(h3).The following ansatz for eh trivially satisfies the above condition.eh =∞∑l=1hl+1el+1. (3.9)Note that the leading term in the series has a factor of h2. We use the leading orderterm to define the order of the solution error, which for this discretization schemeis O(h2). Assuming linearity of the operator L the truncation and solution errorsare related via,Lhu− fh = τh = Lh(uh + eh)− fh= Lheh. (3.10)We can analyze the principle of convergence using the Richardson’s expansion.Suppose a numerical solution uh has been obtained for a difference equation of thetype given by Eq. (3.3). If a second-order discretization scheme has been employed,the truncation error is given by,τh = τ2h2 + τ3h3 + . . . . (3.11)Based on Richardson’s analysis the solution error in such a case should follow,uh = u+ e2h2 + e3h3 + . . . . (3.12)643.1. Finite DifferencesAt two higher levels of resolution, namely 2h and 4h the solution errors are,u2h =u+ e2(2h)2 + e3(2h)3 + . . . , (3.13)u4h =u+ e2(4h)2 + e3(4h)3 + . . . . (3.14)For what follows it is crucial that we assumed the error calculation is performed atthe same point (y, t) in all three cases. The difference between two levels e.g. h and2h to leading order is,u2h − uh ≈ (22 − 1)e2h2. (3.15)Similarly, for levels 2h and 4hu4h − u2h ≈ 22(22 − 1)e2h2. (3.16)The ratio of these two results defines a quantity called the convergence factorQ(y, t) =u4h − u2hu2h − uh . (3.17)Clearly, in this case Q(y, t) = 4. This defines the order of convergence of thenumerical scheme, for the present example a convergence factor (Q(y, t)) of 4 denotes2nd-order convergence. For illustration purposes it often useful to synthesize theinformation contained in Q(y, t) by computing a spatial norm, typically the l2-norm.The notation ‖·‖l2 denotes the l2-norm, which is defined as,∥∥∥uh∥∥∥2= 1NN∑j=1(uhj )212. (3.18)Applying this to the convergence factor, we obtain a time series,Q(t) =∥∥u4h − u2h∥∥‖u2h − uh‖ . (3.19)653.1. Finite DifferencesIf the solution converges i.e. limh→0 uh = u and it is properly described by theRichardson’s expansion this quantity should remain constant. In the above caseQ(t) = 4. Similarly, we could calculate convergence factors for other orders of thediscretization. It is easy to show that for 1st-order schemes Q(t) = Independent ResidualEven if the numerical algorithm passes the consistency and convergence tests out-lined above, it does not guarantee that uh is converging to the desired continuumsolution u pertaining to system (3.1). Since the above tests only indicate a narrow-ing of the gap between uh and some continuum function as h→ 0, it is possible thatuh approaches to the wrong continuum solution. The independent residual test isdesigned to avoid this pitfall.To illustrate this point we consider the following. Suppose that a second-orderscheme (τh ∼ O(h2)) has been successfully implemented such that a numericalsolution uh to the FDA (3.3) has been obtained. This implies,• consistency: limh→0(Lhu− fh) = 0,• convergence: limh→0(u− uh) = 0.There is considerable freedom in the formulation of the FDA. This follows from thevariety of ways that we can use Taylor series expansions to construct Lh. Therefore,we can easily generate an alternative difference operator L˜h. Like Lh, this new oper-ator must also satisfy the consistency requirement, limh→0 L˜hu−fh = 0. Therefore,a 2nd-order approximation to L in terms of the mesh-spacing h follows,L˜h = L+E2h2 +O(h3), (3.20)where the E2(y, t)s are the coefficients generated by the Taylor series expansion.663.1. Finite DifferencesLikewise, convergence to the continuum solution allows,uh = u+ e2h2 +O(h3). (3.21)Then,L˜huh =(L+ E2h2 + . . .) (u+ e2h2 + . . .)= Lu+ (E2u+ Le2)h2 +O(h3). (3.22)Thus, we see that the operator L˜h should follow the same quadratic convergenceas Lh. Suppose that an error has been introduced in our representation of Lh, butwe have correctly solved Lhuh − fh = 0, satisfying the consistency requirement.Therefore, we still should measured the desired convergence but the solution that isapproached follows,uh = u+ e0 +∞∑l=1elhl, (3.23)where, e0 is an O(1) error in the solution introduced by not having implementedthe correct/desired Lh operator. In the limit h → 0, uh approaches u + e0 insteadof u. The action of L˜h on this numerical solution is,L˜huh = Lu+ L˜he0 +O(h2). (3.24)The O(1) term, L˜he0 clearly violates the consistency requirement, unless L˜he0 = 0which is highly unlikely. Obviously, errors could still creep in via the constructionof L˜h. However, since this operator is only used once after the numerical solutionshas been obtained it is easier to ensure its correctness. Consequently, if uh passesthe independent residual test we can be highly confident that we have the numericalsolution to the differential equation we set out to solve. This test was the primarydiagnostic tool we used to establish the correctness of our numerical results. Formore information, and other examples on the independent residual test see [37, 71,76, 85].673.2. Finite Volumes: Conservative Methods3.2 Finite Volumes: Conservative MethodsThe finite difference techniques described in Sec. 3.1 were used in order to calculatethe Newtonian potential ϕ(~x, t) by solving Poisson’s equation numerically. In thissection we present the numerical methods used in solving the equations of fluiddynamics with sources. As previously mentioned, these equations are susceptible tothe formation of shocks and discontinuities in the fluid variables even from smoothinitial conditions. This issue leads to undefined gradients in the equations of motion.Nevertheless, this problem can be confronted by considering the integral form of thefluid equations. These can be obtained directly by integrating the differential formof these equations over a finite volume Vj . We apply this to a general vector-valued differential conservation law Eq. (2.1), more specifically to the k-th equation.Therefore, ∫Vj∂qk∂tdV +∫Vj∇ · ~fk(q)dV =∫Vjψk(q, ~x, t)dV. (3.25)The quantities qk(~x, t), ~fk(q), ψk(q, ~x, t) introduced in Chap. 2 represent the k-th conserved variable, flux and source functions respectively. Application of thedivergence theorem yields,∂∂t∫VjqkdV +∮Sj~fk · nˆdS =∫Vjψk(q, ~x, t)dV. (3.26)Sj is the surface which encloses the volume Vj, and nˆ is a vector normal to thissurface. The derived integral equation (3.26) does not involve gradients of thedynamical variable qk, and so discontinuities in qk do not lead to ill-defined terms.The above equation (3.26) is a typical mathematical statement of conservation. Thechange in quantity qk inside a cell of volume Vj equals the flux over the boundaryof such cell Sj plus a source or sink of qk inside Vj.We can formally integrate equation (3.26) over a time interval ∆t = tn+1 − tn683.2. Finite Volumes: Conservative Methodsto produce,∫Vjqk(~x, tn+1)dV −∫Vjqk(~x, tn)dV +∫ tn+1tn∮Sj~fk · nˆdSdt =∫ tn+1tn∫VjψkdV dt. (3.27)Dividing Eq. (3.27) by the volume element Vj allows us to identify the quantity,Qk,j(t) ≡ 1Vj∫Vjqk(~x, t)dV, (3.28)which is the average value of qk over the finite cell Vj . To make this discussionmore concrete in the following section we discuss the integral conservation law in agiven coordinates system. This allows us to generate a discrete equation for the cellaverages Qk,j.3.2.1 Spherical Coordinates: Spherical SymmetryWe considered writing Eq. (3.27) in spherical polar coordinates. For simplicity weassume spherical symmetry i.e. all the physical fields e.g. qk are independent ofthe coordinates θ and φ. Under this assumption it is convenient to define the cell’svolume as Vj = 4π(r3j+1/2 − r3j−1/2). Which is a spherical shell of inner and outerboundaries given by rj−1/2 ≡ rj −∆r/2 and rj+1/2 ≡ rj +∆r/2, respectively. Notethat ∆r is considered small. For this geometry, the surface integral in Eq. (3.27)becomes,∫ tn+1tn∮Sj~fk · nˆdSdt =∫ tn+1tn(∫Sij+1/2fkdS −∫Sij−1/2fkdS)dt.=∫ tn+1tn4πr2j+1/2fk(t, rj+1/2)dt−∫ tn+1tn4πr2j−1/2fk(t, rj−1/2)dt, (3.29)693.2. Finite Volumes: Conservative Methodswhere Sj±1/2 denotes the outer/inner boundaries of the cell, and fk is the radial com-ponent of the flux vector corresponding to the kth conservation law. The notationcan be simplified through the following definition,Fn+1/2k,j+1/2 =1∆t∫ tn+1tnfk(t, rj+1/2)dt. (3.30)This quantity is known as the numerical flux. With the additional definitionΨn+1/2j =1∆tVj∫ tn+1tn∫VjψdV dt, (3.31)we finally arrive at a discrete equation for the conservation laws where the dynamicalvariables are the cell averages Q(x, t),Qn+1j −Qnj +3∆tr3j+1/2− r3j−1/2(r2j+1/2Fn+1/2j+1/2 − r2j−1/2Fn+1/2j−1/2)= ∆tΨn+1/2j . (3.32)We have suppressed the index k to keep the notation from becoming unmanageable.It is implied that Eq. 3.32 represents the kth component of a vector-valued nonlinearconservation law. The components are coupled through Fn+1/2j+1/2 and Ψn+1/2j .3.2.2 Cylindrical Coordinates: Axial SymmetryDue to its relevance to the calculations performed in this thesis we expressed theintegral conservation law, Eq. (3.26) in cylindrical coordinates (s, φ, z). The phys-ical system is assumed to be axisymmetric, so all the fields e.g. qk are independentof the azimuthal angle φ (if the axis of symmetry coincides with the z-axis). Inthis situation there are two non-trivial flux vector components for the kth equation,namely, f sk and fzk (the other component fφk is trivial in the sense that it does notappear in the equations of motion since, ∂fφk /∂φ = 0).The cell boundaries are given by [si−1/2, si+1/2]× [zj−1/2, zj+1/2], while its vol-703.2. Finite Volumes: Conservative Methodsume, Vij = 2π(s2i+1/2 − s2i−1/2)(zj+1/2 − zj−1/2). Once again, si±1/2 = si ± ∆s/2,where ∆s is a small parameter. The same holds for zj∓1/2 = zi ± ∆z/2 with ∆zalso being a small parameter, typically ∆s = ∆z. The surface integral in Eq. 3.26becomes∫ tn+1tn∮S~fk · nˆdSdt =∫ tn+1tn(∫ zj+1/2zj−1/22πsi+1/2fsk(t, si+1/2, z)dz−∫ zj+1/2zj−1/22πsi−1/2fsk(t, si+1/2, z)dz)dt+∫ tn+1tn(∫ si+1/2si−1/22πf sk(t, s, zj+1/2)sds−∫ zi+1/2zi−1/22πf sk(t, s, zj−1/2)sds)dt. (3.33)Similarly, following the definitions of the numerical fluxes, we obtainF si±1/2,j =1∆t∆z∫ tn+1tn∫ zj+1/2zj−1/2f s(t, si±1/2, z)dzdt, (3.34)F zi,j±1/2 =2∆t(s2i+1/2 − s2i−1/2)∫ tn+1tn∫ si+1/2si−1/2f z(t, s, zj±1/2)sdsdt. (3.35)Notice, we have suppressed the superscript n + 1/2 which denotes a temporal av-erage, as well as the index k which indicates our discussion applies to an arbitrarycomponent of a vector of conservation laws. The average of the source function overthe cell Vij is,Ψn+1/2ij =1∆tVij∫ tn+1tn∫VijψdV dt. (3.36)Combining, Eqs. (3.33)–(3.36) along with Eq. (3.28), we arrive at equation,Qn+1ij −Qnij +2∆t(r2i+1/2− r2i−1/2)(si+1/2Fsi+1/2,j − si−1/2F si−1/2,j)+∆t∆z(F zi,j+1/2 − F zi,j−1/2)= ∆tΨij. (3.37)713.3. The Riemann ProblemThis is the discrete version of the conservation law, Eq. 2.1 subject to axial symme-try. The dynamical variables are the cell-averages Qij.Eqs. 3.32 and 3.37 are finite volume discrete formulations of the conservationlaws whose solutions are the cell-averages. A significant impediment in solving forthe cell averages is the fact that the numerical fluxes e.g. Fj±1/2 are functions ofthe solutions to the continuous problem qk, which are unknown. Effective numericalmethods provide adequate representation of the numerical fluxes so that the finitevolume equations can be used to solve for the cell averages. A brief description ofthese methods is given below, Sec. The Riemann ProblemThe source-free, non-linear conservation-law,∂tq+ ∂xf(q) = 0, (3.38)where q = (q1, . . . , qN )⊺, f = (f1, . . . , fN )⊺, and the initial data,q(0, x) = qL, x < 0qR, x > 0(3.39)represents a general statement of the Riemann problem. By definition the Riemannproblem is one-dimensional; it is extensively discussed in [61, 101]. Here we presentonly the features which have relevance to our hydrodynamical model. The conditionsentailed by the Riemann problem are typical of certain problems in physics andengineering i.e. those involving fluid dynamics.The solution to the Riemann problem clearly depends on the form of the fluxfunction f(q). For the scalar case (q = q) with initial data qL > qR (Fig. 3.1) andf(qL), f(qR) > 0, a formal solution can be provided. The Rankine-Hugoniot jump723.3. The Riemann Problemcondition,f(qR)− f(qL) = s · (qR − qL), (3.40)gives the shock speed, s (see Ref. [61] for more details). At later times the solutionis described by,q(t, x) = qL, x < stqR, x > st.(3.41)Alternatively, if qL < qR with f(qL), f(qR) > 0 the solution to the scalar Riemannproblem is given by,q(t, x) =qL, x <∂f(qL)∂qtqˆ(x/t),∂f(qL)∂qt < x <∂f(qR)∂qtqR, x >∂f(qR)∂qt.(3.42)The discontinuity disappears after the initial time, the two state qL and qR are con-nected by a continuous self-similar state qˆ(x/t). This is known as a rarefaction fan[101].The problem of direct relevance to this project is the Riemann problem for avector-valued, nonlinear, coupled system given by Eq. 3.38. Discontinuous initialdata can be given through any combination of the components for the left andright vectors qL and qR. The solution of the vector-valued Riemann problem cansimultaneously yield shocks, and rarefaction waves, as well as contact discontinuities.A possible wave pattern for such initial data is given by Fig. 3.1 In general, theevolution gives rise to intermediate states q∗L and q∗R. The specific values forthese states as well as the characteristic waves (shocks, rarefaction fan, contactdiscontinuities) displayed in Fig. 3.1 together represents the solution to the vectorRiemann problem.Our numerical problem requires that we solve for the cell averages e.g. Qn+1j733.3. The Riemann ProblemyRarefactionFanContact DiscontinuityShocktq∗L q∗RqRqLFigure 3.1: Typical wave pattern solution to the vector Riemann problem (this is thegeneric solution for a three-component vector Riemann problem, i.e. q = (q1 q2 q3)⊺[101]). The straight lines represent the various waves e.g. shocks waves. Theyseparate the different states; knowledge about these states constitute the completesolution to the Riemann problem.6. Implicit in our finite-volume formulations given by Eqs. 3.32 and 3.37, is therequirement that the solution at the interfaces could be accessed in order to constructFj+1/2. In general, adjacent cell averages will not be equal i.e. Qj 6= Qj+1. Azeroth-order extrapolation at the boundary according to Qj = QˆLj+1/2, and Qj+1 =QˆRj+1/2 leads to a situation that is identical to the Riemann problem with QˆLj+1/2 andQˆRj+1/2 playing the role of the discontinuous initial data. Thus the solution to theRiemann problem could be used to obtained the solution at the interface. Recall thatthis is a requirement for the construction of Fj+1/2. Generating an approximationto the solution at the cell borders via this (zeroth-order) extrapolation procedure inconjunction with the Riemann problem solution is known as the Godunov method6The bold-faced symbol represents vector notation. Alternatively, we could use index notationand refer to the k-th component Qn+1k,j as was done in Sec. 3.2743.4. Approximate Riemann Solvers[61]. With this procedure for computing Fj+1/2 the cell averages can be integratedforward in time according to the finite volume formulas Eqs. 3.32 and 3.37. Thespecific way in which the Riemann solution is used to construct the numerical fluxesis discussed in the following sections. In this work, we have used two commonlyused methods, Roe’s approximate Riemann solver and the Harten-Lax-Leer-contact(HLLC) approximate Riemann solver. These are discussed extensively in [61, 88,101] and reviewed below.3.4 Approximate Riemann SolversBefore we introduce the approximate Riemann solvers Roe and HLLC, we shouldaddress the use of the qualifier “approximate” attached to their name. This followsfrom the fact that the numerical fluxes, Fj+1/2 are constructed by solving the Rie-mann problem corresponding to a “linearized” conservation law. The discontinuousdata at the cell boundary is applied to a linear set of conservation laws obtainedthrough a linearization procedure of the original equations of motion.For concreteness, let us write the general one-dimensional vector conservationlaw in quasi-linear form. Starting with Eq. (2.1) in index notation we can get,∂qk∂t+N∑m=1∂fk∂qm∂qm∂y= ψk. (3.43)The quantities involved are defined at the start of Chap. 2. Defining the componentsof the Jacobian matrix is the expression,Ak,m ≡ ∂fk(q)∂qm. (3.44)So that,∂qk∂t+N∑m=1Ak,m∂qm∂y= ψk. (3.45)753.4. Approximate Riemann SolversIn vector notation this equation becomes,∂q∂t+A∂q∂y= ψ. (3.46)Clearly the matrixA is a function of the vector of conservative variables q. However,over very small changes in time, e.g. tn+1 − tn = ∆t we may assume that q isapproximately constant which implies thatA is also approximately constant. If ∆t issufficiently small we may approximate the solution to the non-linear problem by thesolving the linear Riemann problem (Eq. (3.46) with constant A). This solution, ormore specifically the characteristic structure of Jacobian matrix (A(q = constant))is then used to reconstruct the flux function at the cell borders via the HLLCalgorithm (described in Sec. 3.4.1). In this sense the HLLC and Roe are approximateRiemann solvers.In order to solve the linear Riemann problem posed by Eq. (3.46) an appropri-ate representation of A from the states qL, qR is required. This representation isdenoted by Aˆ(qL,qR). The following conditions must be satisfied by Aˆ(qL,qR)[61],1. All the eigenvalues of Aˆ are real i.e. Aˆ is non-singular (diagonalizable).2. Aˆ(qL,qR) −→ A(q) in the limit qL, qR → q3. Aˆ(qL,qR)(qL − qR) = f(qL)− f(qR)The first condition ensures that the system of linear equations is hyperbolic, i.e.finite, real characteristics speeds. The second condition guarantees consistency withthe nonlinear conservation law. Finally, the third condition states that the shockspeed is given by the Rankine-Hugoniot condition, Eq. (3.40). These requirementsare discussed in [61] in extensive detail.763.4. Approximate Riemann Solvers3.4.1 The HLLC Approximate Riemann SolverThe HLLC solver was proposed by Harten, Lax and van Leer [50]. One beginsby assuming a three-wave pattern solution of the vector Riemann problem similarto that displayed in Fig. 3.1. These waves-speeds are labeled as {SL, SR, S∗}.After some short time period ∆t the initial two-state system given by Eq. (3.39)may in general evolve into a four-state solution such as {qL, q∗L, q∗R, qR}. Thespeeds SL and SR are eigenvalues of the Jacobian Aˆ and correspond to the fastestcharacteristics speeds emanating from the discontinuity. Toro et. al. [102], in orderto restore the full characteristic structure of the Euler equations, introduced thethird speed S∗ (another eigenvalue of Aˆ), which is the speed of propagation of thecontact discontinuity separating states q∗L and q∗R (Fig. 3.1). States separatedby a contact discontinuity have characteristic speeds that run parallel the q∗L–q∗Rboundary. See [102].For the specific case of the source-free Euler equations of fluid dynamics, theinitial data for the Riemann problem is given by, qL(pL) and qR(pR) where thestates pL = (ρL, vL, PL)⊺ and pR = (ρR, vR, PR)⊺ are the vectors of primitivevariables. In this case the fastest speeds can be defined according to,SL = vˆ − cˆs, SR = vˆ + cˆs. (3.47)The variables, vˆ and cˆs are the Roe-average7 velocity, and local sound speed definedby,vˆ =√ρLvL +√ρRvR√ρL +√ρR, (3.48)cˆs =((Γ− 1)(Hˆ − 12vˆ2))1/2, (3.49)7This particular average ensures condition 2 is satisfied i.e. Aˆ as qL, qR → q for the Eulerequations773.4. Approximate Riemann Solverswhere H = (E + P )/ρ is the enthalpy, with its Roe-average given by,Hˆ =√ρLHL +√ρRHR√ρL +√ρR. (3.50)Next, the speed of propagation of the contact discontinuity S∗ is given by,S∗ =PR − PL + ρLvL(SL − vL)− ρRvR(SR − vR)ρL(SL − vL)− ρR(SR − vR) . (3.51)With these definitions the numerical flux vector at the interface between cells j andj + 1 is constructed according to,FHLLCj+1/2 =FL if 0 ≤ SL,FL∗ if SL ≤ 0 ≤ S∗,FR∗ if S∗ ≤ 0 ≤ SR,FR if 0 ≥ SR,(3.52)Where,FX∗ = FX + SX(qX∗ − qX), (3.53)and X ∈ {L, R}. Finally, the intermediate ∗ states are given by the formula,qX∗ = ρX(SX − vXSX − S∗)1S∗EXρX+ (S∗ − vX)(S∗ +PXρX(SX − vX)) . (3.54)The generalization to more than one spatial dimension is easily done by solvingthe linearized Riemann problem along each direction independently. Thus, expres-sions (3.47)–(3.54) apply for all the independent directions. For more details on theHLLC solver see [101].783.4. Approximate Riemann Solvers3.4.2 The Roe SolverThe Roe solver relies on a full decomposition of the linearized Riemann problem intoits characteristic structure. Again, for the special case of the Euler equations in onedimension this is done through diagonalization of the Jacobian matrix Aˆ evaluatedat the Roe averages e.g. Eqs. (3.48)– (3.50). The characteristic structure is thenused to construct the numerical flux according to,FRoej+1/2 =12(f(qLj+1/2) + f(qRj+1/2)−∑k|λk|ωkχk), (3.55)where, χk and λk are kth eigenvectors and eigenvalues of the Jacobian. The quan-tities ωks are the solutions to equation,qR − qL =∑kωkχk, (3.56)and are called the “jumps”. They are the coefficients involved in the decompositionof the discontinuity qR−qL into the linear characteristic waves, χk. Lastly, f(qXj+1/2)is the physical flux evaluated at the L/R initial state. In the case of the Eulerequations this is,f(q) =ρvρv2 + P(E + P )v . (3.57)The extension to higher dimensions is again straightforward. Complete informationon the Roe solver can be found in [61]. It is easy to note that the Roe solver iscomputationally more expensive than the HLLC solver since it involves completediagonalization of the Jacobian matrix Aˆ at every cell interface.793.4. Approximate Riemann Solvers3.4.3 Cell Boundary Variable ReconstructionThe Godunov method is 1st-order accurate in the local grid spacing. This is dueto the fact that the L/R-states involved in setting up the Riemann problem areobtained through a zero-order extrapolation of the cell averages {Qj}. For example,in the interface between the jth and the (j + 1)th cells, the L/R states which definethe Riemann problem are Qj and Qj+1 respectively. A higher level of accuracy canbe achieved by providing a better approximation to the L/R states from the cellaverages. The goal of high resolution methods is to accomplish this through a specialextrapolation procedure. Such procedure should take into account the discontinuousnature of the data i.e. the presence of shocks, and avoid the introduction of spuriousoscillations.A linear piece-wise reconstruction of the fluid variables at the cell boundariesleads to a 2nd order accurate (in the spatial grid spacing) computational scheme.In order to avoid spurious oscillation near shocks the scheme is only 1st-order ac-curate at such points, and also at local extrema. These techniques are known asessentially-non-oscillatory (ENO) methods [94]. We should clarify that in practicethe extrapolation is performed on the cell-averaged primitive variables ({Pj}). Fromthese we then compute the L/R conservative variables of the Riemann problem. Thisprocess can be written symbolically as follows,PLj+1/2 = Pj + σj(yj+1/2 − yj), (3.58)PRj+1/2 = Pj + σj(yj+1/2 − yj+1). (3.59)Here, Pj8 are the cell-averaged primitive variables, with yj being its coordinatelocation in the one-dimensional discrete domain, and σj is called the slope limiter,8Whenever we speak of the discrete conservation laws we use the symbols Pj and Qj to referto the primitive and conservative variables. Lower-case variables p and q are their correspondingcontinuum representation. We also use the lower-case symbols when referring to the generic solutionof the Riemann problem.803.5. Time IntegrationIt is chosen from the slopes,sj−1/2 ≡Pj −Pj−1yj − yj−1 , sj+1/2 ≡Pj+1 −Pjyj+1 − yj , (3.60)according to,σj = minmod(sj−1/2, sj+1/2), (3.61)where, applied to each component separately,minmod(a, b) =0 if ab < 0a if |a| < |b| ab > 0b if |b| < |a| ab > 0.(3.62)The above, Eq. (3.62) defines the “minmod” slope limiter. Other ENO such limitersexists, see [72, 76]. Throughout, this project we have exclusively used the minmodlimiter to perform the variable reconstruction. Finally, the conservative variablesQLj+1/2, QRj+1/2 can be computed from PLj+1/2, PRj+1/2 via the definition, Q(P).From the definition of the minmod limiter it follows that at shocks and extremaσj = 0, thus the extrapolation is piece-wise constant and locally the scheme is1st-order accurate.3.5 Time IntegrationThe time integration of the equations of motion (Euler+Poisson equations) was donethrough a two-step explicit method. This procedure is called Huen’s method and itis 2nd-order accurate in the time step ∆t. In one spatial dimension, e.g. sphericalsymmetry, it can be simply stated as,Q˜j =Qnj +∆tLˆ(Qn), (3.63)Qn+1j =12(Q˜j +Qnj +∆tLˆ(Q˜)), (3.64)813.5. Time Integrationwhere we have lumped together all the complicated operation involved in the spatialintegration and source terms into the operator Lˆ(Qn).Since we also performed calculations in axial symmetry, the methods we neededto apply are effectively two-dimensional. In this case the time integration throughHuen’s two step predictor-corrector method becomes,Q˜ij =Qnij +∆tLs(Qn) + ∆tLz(Qn), (3.65)Qn+1ij =12(Q˜ij +Qnij +∆tLs(Q˜) + ∆tLz(Q˜)), (3.66)where Ls(Qn), Lz(Qn) are the spatial integration i.e. solution to the linearizedRiemann problem along the two independent directions s and z (in cylindrical co-ordinates (s, φ, z)). The indices i and j refer to the location (si, zj). The sourceterms are included in either or both of the terms Ls(Qn) and Lz(Qn).An important property of Huen’s two-step method is that it is known to beTotal-Variation-Diminishing (TVD). It is a statement regarding the stability of thealgorithm. The concept was first introduced by Harten [49], with the total variationbeing defined as,TV(Qnj ) =∑j|Qnj+1 −Qnj |. (3.67)or in two-dimensions9 [27, 61],TV(Qnij) =∑i,j∆s|Qni+1,j −Qni,j|+∆z|Qni,j+1 −Qni,j|. (3.68)The method is said to be TVD if TV(Qn+1ij ) ≤ TV(Qnij) and TV(Qnij) is monotoni-cally decreasing in time.9Again, we assume axial symmetry and we use cylindrical coordinates (s, φ, z).823.6. The Grid Structure3.6 The Grid StructureHere, we present the grid structure used in approximating the solution to Eqs. (2.5)–(2.8) in spherical and axial symmetry. The spherically symmetric grid along with itsadaptive-mesh-refinement (AMR) facility closely follows the algorithm used in [76],where a single dynamical, non-uniform grid was used. We provide a brief summaryof this algorithm in Sec. 3.6.1. For the calculations in axial symmetry we used amultitude of uniform grids. The dynamically chosen grids ensure uniform truncationerror over the entire computational domain. This AMR algorithm was developedby Berger and Collela [8].3.6.1 Spherical Symmetry Non-uniform GridA nonuniform grid structure was used to discretized the spherically symmetric com-putational domain. This algorithm was developed and previously used in [76],therefore we will not describe it here in great detail. The grid domain (distancer from the origin) is divided into three regions, Ωa, Ωb, Ωc. The innermost regionΩa : 0 < r ≤ ra consists of a uniform grid of Na points and uniform spacingri+1 − ri = ∆ra. There are Nb points in Ωb : ra < r ≤ rb spaced according toln(ri+1) − ln(ri) = ∆R, with uniform ∆R. The transition from the uniform regionto the logarithmic, nonuniform region can be made smooth by choosingeRa+∆R − eRa = ∆ra, (3.69)where Ra = ln(ra). This is required so that we can use the same routines, e.g. cell-boundary variable reconstruction (Sec. 3.4.3) across the Ωa − Ωb domain interface.Lastly, Ωc : rb ≤ r ≤ rc is again a uniform grid of Nc points and uniform spacingri+1 − ri = ∆rc. The transition between the nonuniform (logarithmic) grid and the833.6. The Grid Structureuniform grid is again smooth provided we impose,eRb − eRb−∆R = ∆rc, (3.70)where, Rc = ln(rc).To completely specify the grid it is sufficient to provide the parameters {Na, Nb,Nc, ∆ra}. While the other grid parameters are given by,∆R = ln(Na + 1Na), ∆rc = ∆ra(Na + 1Na)Nb−1. (3.71)These follow from conditions (3.69) and (3.70). Trivially,ra = Na∆ra, rb = (Na + 1)∆rc, rc = Nc∆rc + rb. (3.72)The AMR procedureDuring gravitational collapse of an ideal gas, the solution, as characterized by thefields ρ, ~v and P evolves over a wide range of spatial and temporal scales. Properresolution of these requires dynamical refinement of the mesh-spacing. For thespherically symmetric code, this requirement was accomplished by utilizing essen-tially non-oscillatory (ENO) methods in order to construct a 3rd order interpolatingpolynomial that was used to add points in regions of high gradients, see [93] formore details regarding ENO-methods.During the evolution, if the normalized gradient i.e.1||ρ||l2∂ρ∂r, (||ρ||l2 is the l-2norm of the density field) exceeded a certain threshold, then the AMR procedure wasapplied, otherwise the evolution was continued. The application of the techniquecan be outlined as follows:• Reduce the grid spacing by a factor of 1/2 in the uniform grid such thatra → ra/2, keeping the number of grid points in this region constant, butreducing ∆ra → ∆ra/2. The new points are interpolated using the 3rd order843.6. The Grid StructureENO polynomial.• The region left behind in the reduction ra → ra/2 is interpolated and incor-porated into the logarithmic grid. Therefore, Nb → Nb + N ′b, where N ′b =NINT(ln(2)/∆R). The function NINT() operates on a real number to providethe nearest integer. The interpolation in this region is also done via the 3rdorder ENO interpolant.• Repeat process after a set number of evolution time steps.This process doubles the resolution over a collapsing volume near the origin in orderto adequately resolve the form the fields ρ, v, P and ϕ. This AMR algorithm is verysimilar to the one used in [76].3.6.2 Axial Symmetry GridWe used cylindrical coordinates (s, φ, z) in computing the evolution of the axisym-metric fluid. The axis of symmetry was chosen to coincide with the z-axis. Thedomain is defined by [0, smax]× [−zmax, zmax]. This is known as the base grid. Thetotal number of points in the base grid is given by, N0s × N0z. These integer pa-rameters define the uniform base grid spacing through,∆s0 =smaxN0s, ∆z0 =2zmaxN0z(3.73)The physical boundaries are set by smax and zmax, these are chosen such thatϕ ∼ 1√s2max + z2max<< 1. (3.74)Recall that ϕ(s, z) is the Newtonian potential in axisymmetry.The AMR in this case is accomplished through application of Berger-Collela’salgorithm [8]. Unlike the single grid used in the spherically symmetric calculations,here, a nested hierarchy of grids is used in order to resolve the high gradient regions.853.7. Solution to Poisson’s EquationThe hierarchy of grids is evolved. At the highest level of resolution, a new grid canbe added or deleted based on a prescribed truncation error criterion. We will notexplain the algorithm here, instead we refer the reader to the original work of Bergerand Collela, [8].Given the large dynamical range involved in simulations of gravitational col-lapse and the localized nature of the phenomena, the use of AMR greatly expeditesthe computations. For example, the density amplitude evolves through the range∼ 1–107, for a typical critical collapse simulation. Nevertheless, the search for acritical solution requires dozens of runs, each taking longer as the critical solutionis approached. On a single CPU the computation time would be an insurmountableobstacle even with the use of AMR. As a result, we were also required to carry outour calculations in parallel. We used a computational infrastructure called PAMR(Parallel Adaptive Mesh Refinement) developed by Frans Pretorius [21, 84]. Thisfacility allows for each grid in the hierarchy to be partitioned into smaller rectan-gular grids and distributed among a set number of processing units on a computercluster. Synchronization and interpolation of the boundary points after every stepin the evolution is handled by the PAMR program [84].3.7 Solution to Poisson’s EquationIt is simple to show that a 2nd-order discretization of the spherically symmetricPoisson’s equation, Eq. (2.27) can be put into the linear form,MΦ = 4π̺, (3.75)whereM turns out to be a tridiagonal matrix, and the vectors Φ = (ϕ1, . . . , ϕj , . . . , ϕN ),̺ = (ρ1, . . . , ρj , . . . , ρN ) store the discrete value of the Newtonian potential and thefluid density respectively. Note that N ≡ Na + Nb + Nc. Solution to Eq. (3.75)is easily obtained by inverting the matrix M. This was done using the standard863.7. Solution to Poisson’s EquationLinear-Algebra package (LAPACK) routines.The boundary and regularity conditions are incorporated into the tridiagonalstructure ofM. We require that ϕ(t, r) be an even function of r, i.e. ∂ϕ(t, 0)/∂r = 0.The discrete form of the 2nd derivative which shows up in M at the r = 0 becomes,1r2∂∂rr2∂ϕ∂r∣∣∣∣r=0→ 6(ϕ2 − ϕ1)∆r2, (3.76)with ϕ1 ≡ ϕ(0, t) and ϕ2 ≡ ϕ(∆r, t). At the outer boundary we impose the conditionthat the potential should fall asrϕ→ constant at large r.In terms of the discrete system we get that,rNϕN = rN−1ϕN−1. (3.77)Where N ≡ Na + Nb + Nc. The boundary conditions on the fluid variables e.g. ρare discussed in the Sec. 3.8.The 2-dimensional axisymmetric Poisson’s equation (2.39) was solved using themultigrid method. This methods employs a hierarchy of fine and coarse grids. Theelliptic equations are thereby solved on every grid using Gauss-Seidel relaxation.The solution is either injected from the fine grid to the coarse grid or interpolatedfrom the coarse to the fine grid. The advantage of the methods over straightforwardrelaxation is expediency in acquisition of the numerical solution. Shifting betweendifferent grid levels turns out to reduces the number of relaxation sweeps needed forsolution convergence. The method and its convergence properties are discussed in[16].We required evenness of the axisymmetric Newtonian potential with respect tothe axis of symmetry, coordinate s. This means that the discrete approximation to873.8. Fluid Boundary Conditionsthe 2nd derivative with respect to s at s = 0 follows,1s∂∂ss∂ϕ∂s∣∣∣∣s=0→ 4(ϕ2,j − ϕ1,j)∆s2. (3.78)Once again, at the outer boundary we assume the Newtonian potential follows√s2 + z2ϕ(t, s, z) = constant at large√s2 + z2.In terms of the discrete rectangular domain we have three outer boundaries, so weused the conditions,√s2Ns−1 + z2jϕNs−1,j =√s2Ns + z2jϕNs,j, (3.79)√s2i + z2Nz−1ϕi,Nz−1 =√s2i + z2Nzϕi,Nz , (3.80)√s2i + z22ϕi,2 =√s2i + z21ϕi,1 (3.81)In all of the discrete expressions we have suppressed the time dependence.3.8 Fluid Boundary ConditionsThe numerical solution for the fluid variables correspond to the cell-averaged values.Naturally, these averages are located at the cell centers of the discretized domain.Notice that given either spherical or axisymmetric geometry the origin, r = 0 or(s = 0, z = 0) cannot correspond to a cell center. Thus, the regularity conditionsfor the fluid variables which are valid at the origin must be somehow extrapolatedto a distance that is half of cell away from the origin (∆r/2 in spherical symmetry),i.e. the location of the first cell average.883.8. Fluid Boundary Conditions3.8.1 Fluid Regularity Conditions in Spherical SymmetrySo-called conservative regularity conditions for this type of cell-centered data wereproposed by [76]. In spherical symmetry this is done by first assuming the fullsolution (not the cell-averages) follows a Taylor series expansion near the origin.The solution is thus approximated by a polynomial expansion of degree M− 1 inaccordance with,qj(r) =M−1∑m=0am(r − rj)m. (3.82)The j index denotes the discrete cell locations near r, and am is the vector or unde-termined coefficients. There areM such vectors given Eq. (3.82). However, since weare concerned with extrapolating cell-averages, these are obtained by integration,via,Qj(ri) =1Vi∫ViqjdV=3r3i+1/2 − r3i−1/2M−1∑m=0am(∫ ri+1/2ri−1/2(r − rj)mr2dr). (3.83)The index-j denotes the point at which the solution is Taylor expanded, whereasthe index-i is the location of the cell-average being extrapolated. We demand thatthese averages obey (in a local sense) the conservation properties of the numericalsolution. Therefore, the cell-averages calculated via expansion (3.83) must equal thenumerical solution obtained by solving the finite volume equations e.g. Eq. (3.32).These equivalence allows us to determine the vector coefficients am. If the expansionis done about the origin, i.e. rj = 0, and we keep terms in the expansion up toorder M = 4 then determining the ams lead to the following relationships among893.8. Fluid Boundary Conditionsthe cell averages.Even: Qi =11627(3311Qi+1 − 2413Qi+2 + 851Qi+3 − 122Qi+4), (3.84)Odd: Qi =136883(35819Qi+1 − 16777Qi+2 + 4329Qi+3 − 488Qi+4). (3.85)The above relations are obtained following the behavior of the fields near the origin.Even fields such as the fluid density ρ are finite at r = 0 and have vanishing firstderivatives i.e. ∂ρ(t, 0)/∂r. Whereas odd fields such as the fluid velocity v vanishat the origin but have finite spatial first derivative. Relations (3.84) and (3.85) wereused to extrapolate the cell-averages for the even and odd fields near r = 0. Notethat ith cell-average is the one being extrapolated, thus it corresponds to the firstphysical cell in the domain; in general, there could be additional ghost cells to theinterior of ri. In our numerical work i = 1 and so there are no additional ghostcells. Therefore special care is taken in the computation of the numerical fluxes at1 + 1/2-interface 10. These techniques were first used in [76].3.8.2 Fluid Regularity Conditions in Axial SymmetryEquivalent conservative regularity conditions were applied to the axisymmetric fluid.To our knowledge this constitutes the first time this is done for a 2-dimensionalaxisymmetric fluid code. The condition of regularity of the fluid fields in this caseextends along the axis of symmetry, the z-axis in cylindrical coordinates (s, φ, z).Again, we assume a series expansion in s such as,qj(s, z) =M−1∑m=0am(z)(s − sj)m. (3.86)10The minmod limiter requires 2 cells to the left and right of the cell border. In the case of1 + 1/2 in spherical symmetry the assumptions of evenness or oddness of the fluid variable aredirectly incorporated into the cell-border reconstruction and thus into the numerical flux903.8. Fluid Boundary ConditionsClearly, the solution depends now on the z-axis. This dependence is included in theundetermined coefficients am(z)s. The corresponding cell-averages for this expan-sion can be computed through the integral,Qj(si, zk) =1Vi,k∫Vi,kqjdV=2s2i+1/2 − s2i−1/2M−1∑m=0am(zk)(∫ i+1/2i−1/2(s− sj)mrdr). (3.87)We clarify that Eq. (3.87) computes the i-kth cell average (located at (sj , zk)) fromexpansion of the solution at (sj , zk).Similarly, we impose local fluid conservation of the extrapolated values. Consis-tency requires an equivalence between these cell averages computed using Eq. (3.87)and the numerical solution of the axisymmetric finite volume Eqs. (3.37). If theexpansion is performed on the axis of symmetry (sj = 0) then calculation of thecoefficients ams yield the following relationships between the cell averages near thez-axis,Even: Qi,k =2Qi+1,k − 107Qi+2,k +12Qi+3,k − 114Qi+4,k, (3.88)Odd: Qi,k =11627(1419Qi+1,k − 635Qi+2,k + 161Qi+3,k − 18Qi+4,k). (3.89)The set of cell averages closest to the z-axis corresponds to index i = 1, this is thefirst physical cell. Note that we can then use the next four cell centers to extrapolatethe value of the i = 1 cell average from Eqs. (3.88) and (3.89). Eq. (3.88) imposesthe s = 0 regularity condition for the even variables such as ρ(t, s, z), P (t, s, z) andvz(t, s, z), whereas Eq. (3.89) is used for the odd fields, vs(t, s, z) and vφ(t, s, z). Allof the mathematical expression involved in the discussion of regularity conditions,namely, Eqs. (3.82)–(3.89) are time dependent, this however, was suppressed in theabove equations.Lastly, at the outer boundary we used “outflow” boundary conditions. These913.8. Fluid Boundary Conditionsare simple to implement and are non-reflective [19, 36, 54, 78, 89, 90]. Following theupdate of the interior cells through the prescribed Godunov-type method describedin secs. 3.4 and 3.5 the cell(s) at the outer boundary are updated (copied) usingthe value of the last (outermost) physical cell. In spherical symmetry this is easilyrepresented by,QN ≡ QN−1, (3.90)where N denotes the number of cells such that rN = rmax. The above condition,Eq. (3.90) can be thought of as an approximation to,∂q(r, t)∂r∣∣∣∣r=rN= 0, (3.91)where q(r, t) is the continuous fluid variable. It is easy see that the solution atthe boundary satisfying Eq. (3.91) is q(r − vt); if the radial fluid velocity is alwayspositive, then, q(r − vt) describes a function which is advected in the positive r-direction. In all of our simulations the radial velocity is always positive or zero atthe boundary, thus corresponding to outflow (it is always set to zero at the initialtime). In practice, the boundary is set far away from the center, so, the dynamicsthere do not have enough time to significantly disturb the fluid at the boundary.This minimizes fluid loss due to outflow, maintaining approximate conservation ofmass. For the axisymmetric code these condition is given by,QNx,k ≡QNx−1,k for k = 1, 2, . . . , Nz,Qi,1 ≡Qi,2 for i = 1, 2, . . . , Nx − 1, (3.92)Qi,Nz ≡Qi,Nz−1 for i = 1, 2, . . . , Nx − 1,with the identification that [s1, sNx ] ≡ [smin, smax], and [z1, zNz ] ≡ [zmin, zmax].923.9. Conserved Quantities and Error Diagnostics3.9 Conserved Quantities and Error DiagnosticsTo ascertain the correctness of the numerical solutions we applied the independentresidual test, compliance with conservation laws as well as various consistency checksbetween the spherically symmetric and axisymmetric calculations. First, we checkedthat our numerical algorithm conserves total energy Etotal, total mass Mtotal, andtotal angular momentum Ltotal. The results of our calculation involving criticalcollapse are also checked for consistency with previous, related calculations. Theseare discussed in more detail in Chaps. 4 and Code Validation: Spherical SymmetryWe first checked that our spherically symmetric fluid model conserves total massand energy. In spherical symmetry, these quantities are defined by the integrals,Mtotal(t) = 4π∫ rmax0ρ(r, t)r2dr, (3.93)Etotal(t) = 4π∫ rmax0E(r, t)r2dr, (3.94)where,E(r, t) =PΓ− 1 +12ρv2 +12ρϕ, (3.95)is the total energy density of the fluid element. These integrals were approximatednumerically at every time step. Fig. 3.2 illustrates the conservation property of thesolution. This is shown via convergence of the measured quantities Mtotal and Etotalat four resolution levels. As the resolution is incremented, the conserved quantitiesconverged to a constant value. More specifically, the deviation from the measuredtime averages, 〈Mtotal〉 and 〈Etotal〉 tends to zero as the resolution increases, Fig 3.2.The resolution is adjusted and monitored according to the following procedure.At the lowest resolution level L0 we start with an initial spacing ∆ra0 overthe uniform grid Ωa. Higher levels L1, L2, L3, . . . are initialized according to933.9. Conserved Quantities and Error Diagnostics∆ra0/2, ∆ra0/4, ∆ra0/8. During the evolution this relative resolution is maintainedover the Ωa domain by scripting/recording the times and locations the AMR pro-cedure is activated at L0. Then, this script is read off at levels L1, L2, L3, . . .to appropriately double the resolution in both ∆ra and ∆t at the correct ra andt values. Therefore, over the uniform mesh r ∈ (0, ra] the resolution is preciselydouble at each higher level. However, over the logarithmic region (r ∈ [ra, rb]) theresolution cannot be precisely doubled while satisfying condition (3.69). This is thereason for stating that the resolution in the entire domain is only approximatelydoubled. Nevertheless, in simulations involving gravitational collapse most of theinteresting dynamics occur over region covered by the uniform grid (Ωa) where theresolution is most easily monitored.To clinch the correctness of the numerical solution we checked for convergence ofthe independent residual introduced in Sec. 3.1.1. If the solution obtained is correctthen applying an alternative discretization of the dynamical model to this solutiondefines the independent residual. The independent residuals for the Euler+PoissonEqs. (2.29) and (2.27) in spherical symmetry, are labeled by, Iρ, Iρv, IE and Iϕ areexpected to converge to zero as the resolution is incremented. Note that here werefer to the Cartesian-like, vector form of Euler’s equation in order use the morecompact matrix notation in our definition of the independent residual given by,I =Qn+1i −Qni∆t+f(Qni+1)− f(Qni−1)4∆r+f(Qn+1i+1 )− f(Qn+1i−1 )4∆r+−Sni + Sn+1i2, (3.96)where I ≡ (Iρ, Iρv, IE)⊺ and f(Qni ) are the physical fluxes computed from the nu-merical solution i.e. the cell-averages Qni at time tn and position ri. The aboveEq. (3.96), is a 2nd-order finite difference approximation to the Euler equations andis independent of the methods used calculate the numerical solutions (finite volumeand Godunov-type methods) therefore, it comprises an adequate definition of the in-943.9. Conserved Quantities and Error DiagnosticsFigure 3.2: Plots illustrating the convergence of total mass (Mtotal) and energy(Etotal) towards a constant value as the resolution is increased. Therefore, indi-cating the conservation of total mass and energy for the spherically symmetric hy-drodynamic model. Plotted here are the deviation from the mass and energy timeaverages, 〈Mtotal〉 and 〈Etotal〉 respectively. These are given at four levels of resolu-tion, L0, L1, L2, L3, each successive level has approximately doubled the resolutionof the previous one. These levels correspond to initial (t = 0) resolution over theuniform interior grid (Ωa) of ∆ra0, ∆ra0/2, ∆ra0/4, ∆ra0/8, where our control pa-rameter is ∆ra0. The approximate doubling of resolution (the resolution can only beexactly doubled over interior domain Ωa) was achieved by scripting the regriddingprocess, so that at each level Lj, the mesh refinement is performed at the same loca-tions and at the same time. This accounts for the coincidence in the ‘step’ structureof the plotted datadependent residual. The independent residual for Poisson’s equation (2.27) is givenby the following 1st-order discretization,Iϕ =ϕni+2 − 2ϕni+1 + ϕni∆r2+2riϕni+1 − ϕni∆r− 4πρni . (3.97)These residuals were evaluated at many resolutions levels.953.9. Conserved Quantities and Error DiagnosticsThe resolution levels Lns are again set by scripting i.e. recording the times andlocations the AMR facility is activated on the lowest level, L0 during the evolution.This script is then read-off on L1 to ensure that the resolution is always approxi-mately doubled that of L0, on L2 it is quadrupled, and so on. For our applied HRSCmethod, the expected convergence of the independent residuals’ l2-norm is 2nd-orderexcept at shocks and extrema. At these locations only 1st-order convergence is ex-pected. If we scale the l2-norms at levels L1, L2, L3. . . . by a factor of 2, 4, 8, . . .respectively, then the l2-norms should overlap if the convergence is of order 1. For2nd-order the scaling factors should be 4, 16, 64, . . . . Since our dynamical modeldoes involve shocks and extrema we expect less than 2nd-order convergence. Wescaled the l2-norms of the independent residuals according to order 1 and plottedthe results in Fig. 3.3. We clearly see higher than 1st order convergence. The initialdata corresponding to these convergence test is discussed in Table 4.4, Model-A.3.9.2 Code Validation: Axial SymmetryIn the absence of any external torques the axially symmetric fluid conserves totalangular momentum ~Ltotal, in addition to total mass (Mtotal) and energy (Etotal). Incylindrical coordinates these conserved quantities are given by the following inte-grals,Mtotal(t) =2π∫ zmax−zmax∫ smax0ρ(t, s, z)sdsdz, (3.98)|~Ltotal(t)| =Lz = 2π∫ zmax−zmax∫ smax0ρ(t, s, z)vφ(t, s, z)s2dsdz, (3.99)Etotal(t) =2π∫ zmax−zmax∫ smax0E(t, s, z)sdsdz. (3.100)These integrals were numerically approximated at every time step, this served toconfirm the conservation of mass, angular momentum and energy. The plots pro-vided in Fig. 3.4 clearly illustrate this for a sample run corresponding to Γ = 1.00001.The measured deviation from the time average quantities, i.e. |Mtotal − 〈Mtotal〉 |,963.9. Conserved Quantities and Error DiagnosticsFigure 3.3: The independent residual convergence test for the spherically symmetricfluid. The l2-norms of the independent residuals corresponding to the three conser-vation laws (fluid model) and Poisson’s equation (Newtonian gravity) are presentedhere at four levels of resolution L0, L1, L2, L3. At each subsequent level the reso-lution is doubled. These panels illustrate the l2-norms of the independent residualsIρ, Iρv, IE and Iϕ. These l2-norms at levels L0, L1, L2, L3 were scaled by theirrespective factors 1, 2, 4, 8, such that 1st-order convergence of the solution wouldmake the scaled l2-norm data were overlap. Since the HRSC methods we used arespatially 1st-order at shocks and extrema and 2nd-order everywhere else, we expectthat the l2-norms to indicate higher than 1st-order convergence but less than secondorder. This fact is shown in the above panels|~Ltotal−〈|~Ltotal〉|, and |Etotal−〈Etotal〉 | converges to zero with increasing resolution.These quantities are plotted in Fig. 3.4 at four levels of resolution (L0, L1, L2, L3).The mesh spacing (∆x, ∆z) and time step (∆t) is multiplied by a factor of 1/2 at973.9. Conserved Quantities and Error Diagnosticseach successive higher level.Figure 3.4: Plots of convergence of the conserved quantities for the axi-symmetricfluid. These plots indicate conservation of total mass (Mtotal), angular momentum(~Ltotal) and energy (Etotal). Once again, the deviation from their respective aver-ages, 〈Mtotal〉,〈|~Ltotal|〉, and 〈Etotal〉 converge to zero as resolution improves. Fourresolution levels, labeled by L0, L1, L2, L3 were tested where level Lj has twice theresolution of level Lj−1. Again, the regridding procedure was scripted so that theAMR activation is coincident in both space and time in all the refinement levelstested. This accounts for the observed ‘step’ structure in the data.As was done for the spherically symmetric fluid, the independent residual testwas applied to the axisymmetric fluid dynamical equations. If applied correctly,98Floorthe test allows us to claim, with very high degree of certainly that the employednumerical algorithm has indeed provided a solution to the equations of motion [16].The corresponding 2nd-order discretization used to define the independent residualfor the axi-symmetric Euler equations (2.40) isI =Qn+1i,k −Qni,k∆t+f s(Qni+1,k)− f s(Qni−1,k)4∆s+f s(Qn+1i+1,k)− f s(Qn+1i−1,k)4∆s+f z(Qni,k+1)− f z(Qni,k−1)4∆z+f z(Qn+1i,k+1)− f z(Qn+1i,k−1)4∆z− Sni,k + Sn+1i,k2. (3.101)Similarly, I ≡ (Iρ, Iρvs , Iρvφ , Iρvz IE)⊺. Note that f s and f z are the physical fluxesalong the s and z directions, respectively. The independent residual for the axi-symmetric Poisson’s equations is,Iϕ =ϕni+2,k − 2ϕni+1,k + ϕni,k∆s2+ϕni,k+2 − 2ϕni,k+1 + ϕni,k∆z2+1siϕni+1 − ϕni∆s− 4πρni . (3.102)Measurements of the l2-norms for the independent residuals of Euler’s and Poisson’sEqs. (2.40) and (2.39) in axial symmetry are provided in Fig. 3.5. Again, we observedhigher than 1st-order convergence of the independent residual in accordance withHRSC methods. These sample convergence test pertain to octant-symmetric initialdata, more specifically we used model-A in Table 5.1, with Γ = 1.00001.FloorTo conclude the chapter on the numerical methods we mention briefly the treatmentgiven to ‘vacuum’ regions in the computational domain. We used the minimum-atmosphere floor technique to ensure the fluid velocities are defined everywhere99FloorFigure 3.5: Independent residual convergence tests for the axi-symmetric fluidmodel. Plotted above are the l2-norms of the independent residuals for the Eu-ler equations (hydrodynamics) and Poisson’s equation (Newtonian gravity) in axialsymmetry. These residuals are labeled, Iρvs , Iρvφ , Iρvz , Iρ, IE and Iϕ. The inde-pendent residual data is presented at four resolution levels, L0, L1, L2, L3. At eachsubsequent level the resolution is doubled. The l2-norms at levels L0, L1, L2, L3were scaled by their respective factors 1, 2, 4, 8, so if the convergence were 1st or-der everywhere and at all times the rescaled norms would approximately overlap.Since the methods used to solve the equation of hydrodynamics motion are spatiallysecond-order except at shocks and extrema (HRSC methods) we expect less thansecond order convergence but higher than first order.at all times. The underlying assumption is that these region(s) are dynamicallyinsignificant due to their low energy content. In all of our simulations the minimum-atmosphere, or “floor” for the fluid density was set to ρatm = 10−12 and Patm =100Floor10−12 for pressure. In the spherically symmetric model, the atmosphere has zeroinitial velocity. In axial-symmetry, only the axial component vφ is initially non-zero in the rarefied (floor) regions. With these conditions, the fluid evolves withoutthe generation of undefined, negative density or pressure regions, even after theatmosphere acquires a finite velocity.101Chapter 4Results: Spherical SymmetryIn this chapter we present the results for all the calculations conducted in sphericalsymmetry. We begin with the spectrum of self-similar solutions for the polytropicgas. The nature of these solutions is discussed along with the role of the adiabaticindex Γ. This discussion is followed by an analysis of the spherically symmetriclinear perturbations, where again, the role of Γ is examined. We comment on theisothermal gas limit Γ→ 1 for which the results are well known [56, 60, 64, 81, 108].4.1 Self-similar SolutionsContinuous self-similar solutions satisfying equations (2.52), and (2.53), along withregularity conditions (2.55)-(2.59) were determined via numerical integration. Recallthat this model contains two undetermined parameter α∗ and xs, the amplitude ofthe dimensionless density and location of the sonic point respectively. We found thatanalytic solutions are only possible for specific values of the parameter α∗ = exp(Q0)and xs, thus, the spectrum of self-similar solutions is discrete. The correct valuesfor Q0 and xs were found using a graphical technique initially used in [56] to findthe analytic self-similar solution for an isothermal gas. This technique consists inintegrating Eqs. (2.52), and (2.53) forward i.e. from x ≈ 0 up to a ‘matchingpoint’ xM < xs, as well as backwards from xs to xM . Since Q0 is associated withthe boundary conditions of α(x) and u(x) at the origin it is reasonable that theforward integrations should be parametrized by different choices of Q0. Similarlyxs determines the analytic behavior of the fields at the sonic point therefore, the1024.1. Self-similar Solutionsbackward integrations were parametrized by xs. An example of these integrations interms of the dimensionless velocity over the similarity variable x (u(x)/x) is givenin Fig. 4.1. Curves represented by blue-long-dashed lines correspond to forwardintegration at different values of Q0. The red-short-dashed lines are the backwardsintegrations. Note that these solutions have different values at the matching pointxM (terminating point for both integrations). For a specific choice of Q0 and xs thesolutions obtained by the forwards and backwards integrations join smoothly at xM ,as shown by the black dotted curve in Fig. 4.1. A systematic method to determineall possible smooth matchings of the solution at a prescribed xM is described bellow.The solutions at the matching point are stored in the arrays fQ0 ≡ {α(xM ),u(xM ) Q0}, bxs ≡ {α(xM ), u(xM ), xs}, where clearly fQ0 and bxs contain theforward and backward integration data respectively. For purposes of illustrationwe consider Γ = 1.12 (non-isothermal gas) in our model. For this case, we inte-grated Eqs. (2.52) and (2.53) forwards with choices of Q0 ranging from 0 to about50, whereas for the backwards integrations xs was varied over an interval rangingfrom xM to about 4xM . Hundreds of runs were carried out for each case, the re-sults/solutions at xM (α(xM ), u(xM )) were plotted on a phase diagram given byFig. 4.2. The intersections between the curves fQ0 and bxs correspond to smoothmatchings of the solutions at x = xM = 0.4. These determine the correct choicesof parameters Q0 and xs that yields regular (analytic) solutions. The same analysiswas applied for other values of Γ.The spectrum of solutions can be ordered according to the value of the parame-ter Q0 (from here onward we refer to Q0 instead of α∗ to specify the value α at theorigin). We labeled the solution with the lowest allowed value of Q0 as the “Larson-Penston-type” solution. This solution has some relevance in astrophysical scenariosinvolving core collapse, and has been discussed in [63, 98, 105, 108]. Subsequentmatchings occur at higher values of Q0. This is referred to as the “Hunter-type”branch of solutions, see Fig. 4.2. The limiting (Q0 → ∞) case is a singular, static1034.1. Self-similar SolutionsFigure 4.1: Plot of u(x)/x for various choices of integration parameters Q0 andxs. The blue-long-dashed line represents the numerical solution from x = 0 to aso-called matching point xM (xM = 0.4 in this case). The matching point is chosenso that the autonomous system (Eqs. (2.52) and (2.53)) is well-behaved over theinterval [0, xM ]. Specifically, the denominator in Eq. (2.52) is non-zero. It is notedthat different choices of Q0 yield different solutions at xM . Similarly, the red-short-dashed curves represent the solutions obtained from integrating these equations fordifferent choices of xs (the choice for xs represents a guess for the location of thesonic point). These integrations start at the chosen value of xs and terminate atxM . Note that xM < xs and so the integrations from xs are xM are done in thereverse direction. Again, the solutions at xM vary with xs. Nevertheless, for aspecific choice of Q0 and xs the solutions obtained by the forwards and backwardsintegrations match smoothly at xM , a fact that is illustrated by the black-dottedline. This particular choice of parameters represent an analytic solution to theautonomous system.(u = 0) solution. The structure of the spectrum is essentially identical to that of theisothermal gas [56]. These solutions were labeled “Hunter-type-A, B, C, D, . . . ” in1044.1. Self-similar SolutionsFigure 4.2: Phase space plot for the computed solutions α(x) and u(x) at x = 0.4,where the adiabatic index Γ = 1.12. The curve fQ0 are the results of the forward in-tegrations of equations (2.52), and (2.53) from x ≃ 0 to x = xM = 0.4 parametrizedby Q0; whereas the curve bxs are the results of backward integrations parametrizedby xs starting from the sonic point xs to x = xM = 0.4. The points where the curvescross correspond to the correct preregistration of the solutions. Only a discrete setof parameters {Q0, xs} results in allowed solutions. The Larson-Penston (LP) so-lution has the lowest value for the parameter Q0. As shown the Hunter solutions(HA-HD) are generated as the fQ0 begins to spiral around a common point in thebxs curve. This graphical technique was used in [56].accordance with their increasing central value parametrized by Q0. To our knowl-edge this is the first time these Hunter-type solutions are explicitly determined forΓ > 1. Since the spectrum structure for higher values of Γ is essentially identical tothe isothermal case from here onward we shall call these set of regular solutions assimply the Larson-Penston (LP), and Hunter-A, B, C, D, . . . .1054.1. Self-similar SolutionsThe specific values of Q0 and xs for the LP solution and the first four Huntersolutions are found in Table 4.1. Notice that the Hunter solutions have much highercentral values since α(0) = exp(Q0). The profiles for the self-similar density vari-able α(x), corresponding to these solutions are plotted in Figure 4.3. The velocityvariable u(x) is plotted in Figure 4.4. A distinguishing feature of the Hunter so-lutions are the oscillations in the velocity profiles. The number of nodes in theirprofiles follows from their placement in the spectrum hierarchy, just as was found in[64] for Γ = 1. The Larson-Penston solutions lacks any of these oscillations whichimplies the fluid is collapsing into the origin for all times, i.e. this solution de-scribes coherent core collapse. In contrast, the Hunter-A solutions contains a regionadjacent to the collapsing core where the fluid is not in-falling. Since these plotsrepresent self-similar solutions, this rarefied region “chases” the collapsing core evac-uating fluid away from the origin. Higher order Hunter solutions contain more ofthese increasingly compact regions where the fluid velocity changes direction. Allof these however, contain a rarefaction wave adjacent to a collapsing core as shownin Fig. 4.4.Table 4.1: Values of the parameters Q0 and xs for the similaritysolutions corresponding to Γ = 1.12.Solution Q0 xsLarson-Penston 1.32732279 2.713336217Hunter-A 10.62019209 6.110589148×10−1Hunter-B 14.63542969 1.649524387×10−1Hunter-C 21.75502169 1.012505889Hunter-D 27.18143020 1.3431697064.1.1 Solutions, 1 ≤ Γ < 6/5In what follows we considered the effect of changing the polytropic index Γ on theHunter and LP solutions. We found that the solutions’ central value parameter Q0becomes larger as Γ is increased. The entire spectrum of analytic solutions shifts to1064.1. Self-similar SolutionsFigure 4.3: Plot of self-similar solutions α(x) corresponding to Γ = 1.12. Thesesolutions plotted here on a logarithmic scale are the first four Hunter solutions(A −D) as well as the Larson-Penston solution which has the lowest central value(x = 0). The Hunter family appears to have an infinite and discrete structure withits members having ever-increasing central values. This is similar to the resultspresented in [64] for the isothermal gas.higher Q0 values, whereas the sonic point xs migrates towards the origin. This factmake computation of the solutions increasingly more difficult. The matching point,xM for the Hunter solutions is “squeezed” into a narrow region where the gradientsin Eqs. (2.52), (2.53) are very high.Aside from the LP solution, the Hunter-A solution has the lowest Q0 amongthe Hunter family, so it is the least expensive to integrate since the gradients areless steep than the other family members Fig. 4.2. We computed the parameters1074.1. Self-similar SolutionsFigure 4.4: Plot of self-similar solutions u(x)/x corresponding to Γ = 1.12 An oscil-latory patterns is revealed near the origin for the Hunter solutions. Such oscillationsare absent in the Larson-Penston solution. Furthermore, the number of nodes intheir profiles is ordered. This is similar to the ordering of allowed central densityparameter Q0. This means that the Hunter-A solution has the lowest parameter Q0of the Hunter-family and displays one node, whereas the Hunter-D solution containsfour nodes and has the fourth largest Q0 value. The results shown here are verysimilar to those found in [64] Fig. (2) for Γ = 1.Q0, xs for the Hunter-A solution for various values of Γ starting at the known Γ = 1(the isothermal gas). We found that the center of the dimensionless density variableα(x) increases nonlinearly with linear increments in Γ. The results of the numericalintegration are tabulated in 4.2. The density and velocity profiles for the Hunter-A solutions are presented in Fig. 4.5 and Fig. 4.6. We noticed that the node inthe velocity field also changes non-linearly. This behavior is not shared by the LP1084.1. Self-similar SolutionsTable 4.2: The Hunter-A parameters Q0 and xs for somechoices of Γ ∈ [1, 6/5).Γ Q0 xs1.0 7.45615862 7.390727573×10−11.00001 7.45631897 7.390753939×10−11.04 8.18305221 7.385798698×10−11.12 10.62019209 6.110589148×10−11.16 13.42242634 3.953918340×10−11.18 16.67083350 2.067449815×10−11.19 20.36959616 9.053763430×10−2Figure 4.5: Hunter-A solution for α(x) in the range 1 ≤ Γ < 6/5. In contrast tothe Larson-Penston solution, the Hunter-A solution does not behave linearly withrespect to linear changes in Γ. As illustrated in this plot. This solution (Hunter-A)does not seem to exists for Γ ≥ 6/5, its behavior near this critical value is displayedseparately, Fig. 4.8.1094.1. Self-similar SolutionsFigure 4.6: Hunter-A solution for u(x)/x in the range 1 ≤ Γ < 6/5. The plottedu(x)/x profiles are the Hunter-A solutions chosen at uniformly distributed Γ valuesover this range. We observed the single node which characterizes this solution growsnonuniformly as Γ is varied from 1 to 6/5.solution which as seen in Fig. 4.7 changes linearly over the range 1 ≤ Γ < 6/5.We discovered that as we approach Γ = 6/5, the central value pertaining to theHunter-A solution diverges. This limit was calculated by choosing Γ according to,Γ(i+ 1) =65− 6/5− Γ(i)2i, for i = 0, 1, 2, . . . ,∞. (4.1)Where Γ(i) was used in the ith integration, i.e. Γ is chosen such that interval6/5− Γ(i) is reduced by one half at every trial, as a result each new increment getsexponentially smaller. The results of the integrations were plotted in Fig. 4.8. Thisplot exhibits the exponential growth of the density profiles (α(x)) for the Hunter-A1104.1. Self-similar Solutionssolution in this limit. We concluded that Hunter-A disappears at Γ = 6/5, and giventhat it has the lowest central value, the entire Hunter family vanishes at this point.Our results also indicate that the Hunter solutions do not reappear in the range6/5 < Γ < 4/3. Our investigation looked for matchings in fQ0 and bxs (Fig. 4.2) inthe latter range but failed to find any. A single matching corresponding to the LPsolution persists in the range 6/5 ≤ Γ < 4/3.Figure 4.7: Larson-Penston solution for α(x) in the range 1 ≤ Γ ≤ 6/5. The Larson-Penston solutions was calculated at uniformly spaced values of Γ in this interval.Notice that the dimensionless density profile (α(x)) near the center increases linearlywith respect of a linear increase in Γ. Notice also that we are able to compute thissolution (LP) across the seemingly special value of Γ = 6/5.1114.1. Self-similar SolutionsFigure 4.8: Multiple plots of the Hunter-A solutions at various values of Γ close tothe critical value of Γ = 6/5. Γ is chosen at every trial such that the interval Γ-6/5of the previous trial is reduced by a factor of 1/2. The amplitude of the Hunter-Asolution as shown here increases exponentially.1124.1. Self-similar Solutions4.1.2 Growing ModesWe looked exclusively for growing mode solutions (λ > 0) to Eqs. (2.78)–(2.84).This was done systematically by subjecting the spectrum of regular self-similarsolutions for a given Γ to a linear stability analysis. Equations (2.78)–(2.80) wereintegrated using the LSODE facility in Maple-12. The relative tolerances were setto a maximum of 10−18 for all of our calculations, and we use equation (2.81) tocheck for consistency. The correct values for λ were determined by the shootingmethod-technique. First, the equations were integrated “outward” from x ∼ 0towards the sonic point xs. Near xs, the solutions δα(x) and δu(x) go to ±∞ forsome non-negative choice of λ. Since regularity implies that the solution is finite atthe sonic point, we looked for λ > 0 such that δα(x) and δu(x) are finite and satisfyEq. (2.84) for x ≈ xs.As was expected the possible non-negative solutions for λ depend on the specificself-similar solution being tested. There were no analytic growing modes found forthe Larson-Penston solution at any Γ in the interval [1, 4/3). Varying the growthrate, λ, from 0 to a vary large value (∼ 1032) did not change the blow-up behaviorof δα, δu near the sonic point for the LP solution. This implies that there are nonon-negative values of λ that can satisfy Eq. (2.84) and thus no unstable modes.The Hunter-A solution has precisely one growing mode for Γ ∈ [1, 6/5), where theHunter-A exists. The Hunter-B,C,D,. . . solutions have 2, 3, 4, . . . growing modesrespectively. In the limit Γ → 1 the growth rate, λ, reduces to that computedin [64] for the isothermal gas. We discovered that the overall effect of Γ on thestability of the Hunter solutions is to render them more unstable, i.e. the growthrates becomes larger as Γ increases. Nevertheless, the overall mode structure andhierarchy of the solutions remains unchanged, until we reach Γ = 6/5.Table 4.3 provides a list of the computed numerical values of the unstable modesfor selected self-similar solutions. Notice that the unstable mode for the Hunter-Asolution also seem to diverge at the special value of Γ = 6/5.1134.1. Self-similar SolutionsTable 4.3: Results of the stability analysis performed particu-larly on the Hunter-A solution for the range 1 ≤ Γ < 6/5.Γ Solution Mode λ1.0 Hunter-A 1 9.46371.00001 Hunter-A 1 9.46431.001 Hunter-A 1 9.52471.04 Hunter-A 1 1.2729×1011.12 Hunter-A 1 3.6734×101Hunter-B 1 1.0337×1012 2.6841×102Hunter-C 1 2.1588×1012 3.3496×1023 9.4086×1031.16 Hunter-A 1 1.3456×1021.18 Hunter-A 1 6.4143×1021.19 Hunter-A 1 3.9349×1031.1996 Hunter-A 1 1.228×1081.1998 Hunter-A 1 1.344×1091.1999 Hunter-A 1 1.49×10101.19995 Hunter-A 1 1.67×1011The solution profile, and specifically its amplitude has a similar dependence onΓ, becoming singular at the origin as Γ approaches 6/5 (we fixed the free parameterδu0 = 1 in Eqs. (2.82) and (2.83) in all our calculations). This behavior is expectedsince the self-similar solutions themselves, as was argued in Sec. 4.1.1 are singularthere. A plot of the calculated perturbation functions, δα and δu for the Hunter-Asolutions are presented in Fig. 4.9 in the extreme case where Γ = 1.1999. Again, weperformed a consistency check by setting Γ = 1 and confirming that our solutionprofiles reduce to those provided in [64].At this point we would like to summarize the results presented thus far. Wefound the spectrum of analytic similarity solutions for the polytropic gas model byintegrating (numerically) the autonomous system Eqs. (2.52)–(2.54) and imposingregularity (analyticity) at the origin and the sonic point. The spectrum of solutionsis discrete; it follows a hierarchical structure identical to the regular self-similar1144.1. Self-similar Solutions(a) (b)Figure 4.9: Linear perturbation functions δα(x) and δu(x)/x for the Hunter-A so-lution where Γ = 1.1999. The amplitude of δα(x) provided in panel (a) divergesas Γ → 6/5. However, for Γ < 6/5 this amplitude can be rescaled (normalized),this freedom is entailed by the parameter δu0 in Eqs (2.82) and (2.83). Panel (b)presents the profile for δu(x)/x. The features of mode functions are “squeezed” nearthe originsolutions for an isothermal gas found previously in [56, 64]. In the limit Γ→ 1 oursolutions reduce to the isothermal spectrum. This structure persists for Γ ∈ [1, 6/5).In this range the first member of the Hunter-branch of solutions (Hunter-A) containsa single unstable mode, suggesting its potential role as a critical solution [64, 66, 99].We found that the Lyapunov exponent corresponding to the unstable mode(s) of theall the Hunter solutions increases with increasing Γ. Interestingly, the entire Hunter-branch of self-similar solutions disappears for Γ ≥ 6/5, hinting of a possible changein critical behavior of the polytropic gas at this value of Γ. This was investigatedfrom dynamical simulations; the results are presented in the next section.1154.2. Numerical Simulations4.2 Numerical SimulationsIn this section we present the numerical solutions of the spherically symmetric dy-namical model given by Eqs. (2.24)–(2.27) for 1 ≤ Γ < 4/3. Previously, type-II critical phenomena had been identified in dynamical simulations of Newtoniangravitational collapse of an isothermal gas [48]. The work of Harada et. al. beginsby considering an isentropic hydrodynamic model (equivalent to equations (2.30)–(2.32), with Γ = 1) from the outset. We begin with a more general approachby including the energy density conservation law Eq. (2.26) (Euler equations) andadopting the more generic polytropic ideal gas law, given by EoS (2.18).Initial DataIn most of our calculations, the density and pressure profiles were initialized accord-ing to a Gaussian function. Nevertheless, a second 1-parameter family of initial datawas used in some cases. We labeled these sets as, Models A and B. Their functionalforms are provided in Table 4.4. In both cases, the radial velocity field of the initialdata was set to zero. Note that the control parameter p, modulates the amplitudeof the pressure profile.For dynamical evolutions of Eqs (2.24)–(2.27) involving larger values of the pa-rameter Γ, it was necessary to select initial data that “resembled” the correspondingHunter-A solution. The reason for this selection is made explicit later on. The gen-eral form of this initial data is then given by,Z(0, x) = Z⋆(x) + fp(x). (4.2)Where Z⋆(x) is the Hunter-A solution, and fp(x) is a function that describes thedeviation from Hunter-A. We allow the parameter p to control one aspect of this1164.2. Numerical Simulationsfunction, e.g. the amplitude of a Gaussian profile. Our specific choice for fp(x) was,fp(x) ≡pe−r2+ ερ(r)εv(r)εP (r) . (4.3)The entries in Eq. (4.3) expresses the deviation from Hunter-A for the primitivevariables ρ(0, r), v(0, r) and P (0, r). The variables ερ, εv, and εP are unknownfunctions, representing all other numerical errors associated with approximating theHunter-A solution. These functions are small compared to the manually introducedGaussian perturbation.Table 4.4: Spherically symmetric initial data.Model A Model Bρ(0, r) = e−r2ρ(0, r) =1(1 + r2)2if r < 110−12 if r ≥ 1v(0, r) = 0 v(0, r) = 0P (0, r) = pe−r2P (0, r) ={ p(1 + r2)2if r < 110−12 if r ≥ 14.2.1 Simulations Γ ≈ 1The numerical experiment consisted of evolving 1-parameter families of initial dataof the type given in Table 4.4, until we could unambiguously, identify its final state.The observed initial stages of the evolution, regardless of its final fate are verysimilar. Qualitatively, we observed a deepening of the Newtonian potential well(ϕ(r, t)) with subsequent contraction of the matter at the origin. The central den-sity increases as a core of collapsing matter forms around the center. Following thistransient behavior, we found that for values of our control parameter p less thansome threshold p⋆, the central density continues to grow; the central core collapses1174.2. Numerical SimulationsFigure 4.10: Plot of the similarity parameter Q0 versus central density ρ(t, 0) forvarying degrees of fine-tuning. This plot shows the effect of fine tuning the controlparameter p to the critical value p⋆, for Γ = 1.00001. Note that |p − p⋆| ≡ δp.The value of Q0 characterizes the self-similar solutions (Table 4.1). Q0 = 4πGρct2(with G set to 1 and ρc ≡ ρ(0, t)) computed from the dynamical solutions evolvestowards the self-similar Larson-Penston (LP) value, irrespective of the fine tuning.The fine-tuned evolution approaches the numerical value of Q0 corresponding to theHunter-A (HA) solution at intermediate times. At late times, the spherical unstableperturbation mode interrupts this convergence, and the dynamically computed valueof Q0 again approaches the Larson-Penston value.homologously. The central density grows exponentially, and over exponentially de-creasing time scales. We followed this evolution until the central density had grownover 14 orders of magnitude; at which point the simulation is halted, and we de-clare this final state as the Newtonian analogue of black hole (BH) formation. Incontrast, the subcritical regime, p > p⋆ the central density ceases to grow at somemaximum value ρmax, this is followed by a gradual decrease in the central densityleading, ultimately to complete evacuation of the matter in the central region. Wecall this case, “dispersal”. Once these two regimes have been identified we focusedon the behavior of the solution as we fine tuned our initial data to the threshold ofsingularity “BH” formation, i.e. we fine tune p→ p⋆. We call the solution at p = p⋆the critical solution.1184.2. Numerical SimulationsWe first studied the near isothermal gas Γ = 1.00001, choosing this value of Γallows us to compare our results with those obtained in [48, 64] where an strictlyisothermal gas was considered. A bisection search which partitions the intervalseparating “collapse” and “dispersal” solutions was utilized to zero-in on the criticalvalue p⋆. The search was performed until |p − p⋆|/p⋆ → 10−14. We investigatedthe nature of the collapse at the origin (r = 0) by computing the dimensionlessdensity α(x = 0) from the dynamical solution. For comparison purposes it is moreconvenient to work with Q0 = eα(0). Based on Eq. (2.43) we can computeQ0(t) = ln(4πGρ(t, 0)(t0 − t)2). (4.4)This requires that we know the time of singularity formation t0. In order to identifythis parameter in our simulations we applied a similar procedure to that used byHarada et. al. [48]. First, we assumed the center ρ(t, 0) is at all times collapsing ina self-similar fashion such that ρ(t, 0)(t0 − t)2 is a constant. From this assumptionwe can compute tn0 , the tn-prediction for the time of collapse. Assuming the solutionis self-similar, Q0 is a constant and we have the relation,ρn+1c (tn0 − tn+1)2 = ρnc (tn0 − tn)2, (4.5)where ρnc ≡ ρ(tn, 0). Then, tn0 is used in Eq. (4.4) to determineQ0. If this assumptionis correct, Q0 should be constant or approaching a constant value. More importantly,since Q0 parametrizes the hierarchy of self-similar solutions as discussed in 4.1.1 wecan monitor its possible convergence to any of these values. Another approach tocompute t0 also used in [48] involved directly measuring the time required for thefluid to effectively collapse. Numerically, this was accomplished by measuring theelapsed time during which the initial central density increases over 14 orders ofmagnitude.By fine tuning p in our initial data close to the critical value of p⋆ we found that1194.2. Numerical Simulationsthe computed quantity of Q0(t) converges to the value corresponding to the Hunter-A solution. In fact, the Hunter-A solution is the critical solution as suggested byits characteristic single-unstable-mode. Evidence for this is seen in Fig. 4.10, wherewe have plotted our computed value of Q0(t) against the central density, ρ(t, 0),along with the known values of Q0 for the Hunter-A and LP solutions. As thecentral density increases during critical collapse, the Hunter-A solution becomes an“intermediate attractor” of the evolution. As the density increases the growing modeeventually disrupts this convergence. The growing perturbation effectively “pushes”the intermediate state away from Hunter-A and towards the stable LP similaritysolution. The LP solution as revealed by our mode analysis in Sec. 4.1.2 lacks anygrowing modes, so it is expected that any set of collapsing data will converge to it.This prediction is corroborated by Fig. 4.10. This result is essentially identical theone presented by Harada et. al. in [48]. However, we obtained this result via a moregeneral EoS (2.18).So far the presented results apply only to the central density. The convergence tothe self-similar Hunter-A solutions extends over a region surrounding the center. Aseries of snapshots for critical evolution are presented in Fig. 4.11. The plots pertainto the dimensionless density variable α(x) defined in Eq. (2.43). Generic initial dataconvergences to the Hunter-A solution at intermediate times (t1 − t5) in Fig. 4.11(blue curve). We expect the same local convergence of the other dimensionlessvariables defined in Eq. (2.43). We also looked at the velocity profile u(x) forcritical evolution and found that in fact at intermediate times it matches the profileof the Hunter-A solution. A rarefaction wave envelopes the collapsing core, thus itssize diminishes. This feature of the Hunter-A solution is replicated by the criticalsolution as seen in Fig. 4.12 panels t1 − t5. Over time the unstable mode beginsto dominate the evolution and the solution deviates from Hunter-A. This is againverified by continuing the evolution of critical data past its intermediate convergenceto the 1-mode unstable solution in Fig. 4.11 and Fig. 4.12 (t6 − t9), the profiles of1204.2. Numerical Simulationsboth α(x) and u(x) are described by the LP solution. This late-time behavior isobserved in supercritical evolutions irrespective of fine-tuning.Figure 4.11: Snapshots of the evolution of α(x) for critical initial data using model-A, at Γ = 1.00001. This plot presents the evolution of the dimensionless densityvariable computed from the numerical solution ρ(r, t) according to Eq. (2.46) forcritical initial data (|p− p⋆|/p⋆ → 10−14). The evolution data is represented by theblue curve. The top-left plot t1 represents the initial data whereas the bottom-rightt9 represents the final state. Maximum approach to the 1-mode unstable Hunter-Asolution (dotted line), is observed at intermediate times (t5 in the above panel) Aswe expected as the unstable mode grows the solution finally settles into a profilethat matches the LP solution (dashed plotted), t9.1214.2. Numerical SimulationsFigure 4.12: Snapshots of the evolution of u(x)/x for critical initial data usingmodel-A, at Γ = 1.00001. The plotted u(x) profiles which correspond the criticalsolution (|p−p⋆|/p⋆ → 10−14) are calculated from Eq. (2.45) at each time step. Theflow near the origin replicates the features of the Hunter-A solution after some time(dotted line). In particular the node in the velocity profile which is responsible forthe shrinking of collapsing core. Again the top-left plot is a snapshot of evolution aninstant after the initial time (at t = 0 u(x) = 0), the central panel (t5) correspondsto the maximum attained convergence to Hunter-A solution. Again, we see that atthe final stage (t9) u(t, r) converges to the profile of the LP solution (dashed line).1224.2. Numerical SimulationsPhysical Interpretation of the ResultsThe Hunter-A solutions contains a rarefaction wave that travels towards the origin.This is entailed by the node on the profile of u(x). This wave “catches” up with thecollapsing center precisely at the time of singularity formation t0. Loosely speak-ing the collapsing region (core) becomes infinitely small, even though the densitydiverges at this time (t = t0). The fine tuning of initial data to the critical valuep⋆ is not ideal due to numerical limitations. The data can be above or below thetrue threshold value. For subcritical data the rarefaction reaches the center beforethe singularity forms leading to dispersal of the gas away from the origin. Undersupercritical conditions the rarefaction wave never reaches the center and the coreceases to shrink as the core collapses coherently.The critical data goes through the linear regime discussed in Sec. 2.8.1 (withoutrotation, ~q = 0). The final fate of the data is sealed by the sign of ǫ ≡ (p −p⋆) exp(λ0τ). Collapsing data will have the form, Z(x, τ) = Z⋆(x) + ǫδZ0(x) atthe linear regime (intermediate times). This occurs late enough in the evolutionthat the stable perturbations have died out, but early enough that ǫ is still small.The dispersal situation has the opposite sign, i.e. Z(x, τ) = Z⋆(x) − ǫδZ0(x). Forcollapsed data, the time of departure from the critical solution sets the mass ofthe collapsed core. The longer the initial data converges to the Hunter-A solution,the smaller the core. If our formalism is correct, the scaling of the collapsed massshould follow Eq. (2.121). We selected a series of supercritical initial data evolutionsand measured the mass of the collapsed core. This was done by simply integratingover the spherical shells which have in-falling velocity, over the innermost regionof the cloud. The point at which velocity changes direction defines the radius ofthe core. The integration is performed at late times, once the solution has beganits convergence to the LP solution and the core is collapsing coherently. At thispoint, neither the mass of the core nor its radius change prior to the blow up of thecentral density. This property of the LP solution is what allowed us to define the1234.2. Numerical Simulationscollapsed core. This is the same method used by Harada et. al. in their simulationsof spherically symmetric isothermal gas collapsed [48].A similar test can be conducted for the subcritical regime by measuring themaximum central density attained by the fluid before it began to disperse. Again,perturbation analysis and convergence to the 1-mode unstable, Hunter-A solutionpredicts a scaling of the maximum central density according to Eq. (2.126). Theresults of these calculations are plotted in Fig. 4.13. These plots show the scalingof both the collapsed mass and maximum density. Both plots indicate convergenceto the predicted scaling law based on the existence of a single unstable mode. Theslopes of the lines on the log-log scale of Fig. 4.13 are related to the reciprocal of theLyapunov exponent of the unstable mode and the adiabatic index Γ, as dictated byEqs. (2.121) and (2.126).The results presented so far are consistent with the existence of type-II criti-cal phenomena in the gravitational collapse of a soft fluid (Γ = 1.00001). As wehave stated before, our critical collapse simulations pertained to a a more generalNewtonian hydrodynamic model than previously considered. This is entailed by theinclusion of the energy conservation law, Eq. (2.7) and non-barotropic EoS (2.18).This model allows for generic entropy distributions, however, only under uniformentropy conditions can the self-similar solutions discussed in Sec. 4.1.1 be realized.We found that since the critical solution has vanishingly small mass (a feature oftype-II critical phenomena), then the critical solution originates from a vanishinglysmall region of the initial conditions. Over this region, the entropy is arbitrarilyclose to a constant, unless shocks develop. No shocks seem to arise in any of theinitial conditions we considered. Therefore, the critical solution is isentropic, thusit is reasonable to expect that one to the self-similar solutions of the autonomoussystem Eqs. (2.52)–(2.53) is realized as the critical solutions (Hunter-A).1244.2. Numerical Simulations(a) (b)Figure 4.13: Scaling behavior of the collapsed mass (M) and maximum centraldensity ρmaxc at Γ = 1.00001. The collapsed mass (a) and maximum density (b)were calculated from supercritical and subcritical evolutions, respectively. The solidline is the predicted scaling law behavior near the critical point (|p−p⋆| → 0) derivedfrom Eqs. (2.121) and (2.126). Agreement between the numerical simulation results(blue and red) and the results from perturbation theory (black) improve near thecritical value p⋆, this is expected since the scaling laws only apply near criticality,where all other stable modes have died out. The bisection search ultimately fails aswe approach the numerical precision limit (10−16), and thus we are unable zero-inon the exact value of the critical parameter p⋆. At this point, we cannot reliablymeasure the difference |p−p⋆|. This failure begins to be evident in both sets of dataas |p− p⋆| ∼ 10− Connection with General RelativitySo far our studies of critical collapse of an ideal gas in Newtonian gravity have re-vealed the existence of type-II critical phenomena. This was made evident by theconvergence to a 1-mode unstable critical solutions (Hunter-A) and the particularscaling properties of the collapsed mass and maximum density. The results pre-sented so far pertain to the nearly isothermal gas, Γ = 1.00001. It is clear thatour results are in agreement to those of presented in [48, 64] for the isothermal gasmodel. Previously, Snajdr and Choptuik [96] obtained these results from a General1254.2. Numerical SimulationsRelativistic critical collapse of a perfect fluid with an ultra-relativistic EoS (1.3).This was done by considering the limit k → 0, which, according to the work of Oriand Piran (1990) [77], certain regular similarity solutions of the GR perfect fluidwith EoS (1.3) reduce to the weak-field (Newtonian limit) self-similar solutions ofthe isothermal gas. Snajdr and Choptuik found that the critical solution resembledthe Newtonian Hunter-A solution as k tended to zero. Their normal mode analysisconfirmed that the value of the growing mode also converged to that of Hunter-A.This was previously estimated in [45]. Further confirmation was obtained throughcalculations of the scaling behavior of the black hole mass. Again, the scaling expo-nent agree with our purely Newtonian calculations and those of [48, 64]. We noticedthat the “transition” into the Newtonian regime maintains the solution’s structurebut leads to larger values of the unstable modes. This is evident in the data col-lected in [96]. This trend continues as Γ becomes greater than one, the unstablemode increases as evident in Table 4.3. The Hunter-A solution becomes infinitelyunstable at Γ = 6/5 where this solution structure seems to break down.1264.2. Numerical Simulations4.2.3 Simulations 1 < Γ < 6/5We proceeded with our investigations of critical collapse by considering a stiffer EoS(Γ > 1). We recall at this point the discovered absence of 1-mode unstable regularsolutions for Γ ≥ 6/5. In Sections 4.1.1 and 4.1.2 we found the Hunter-A, 1-modeunstable solutions exist in the range, 1 < Γ < 6/5. The value of the unstable mode,like the amplitude of the Hunter-A solution was found to be dependent on Γ. Both,this amplitude and its unstable mode become singular as Γ→ 6/5 (see Table 4.3).The presence of this 1-mode unstable solution suggests that similar critical behaviorfound at Γ ≈ 1, persists up to, but not including Γ = 6/5.We found similar convergence to the corresponding Hunter-A solution at theselarger Γs during critical collapse. As Γ becomes larger, however, convergence to thecorresponding Hunter-A solution becomes more difficult to ascertain. We attributethis in part to numerical precision limitations which hindered our ability to suffi-ciently fine tune the initial data to the critical parameter 11. More importantly, thevalue of the unstable mode grows significantly as we approach Γ = 6/5; which meansmuch faster growth of the perturbation. For example, in the case of Γ = 1.00001 theunstable mode grows approximately as 1/(t0 − t)9.46, compared to 1/(t0 − t)36.7 forΓ = 1.12. Calculations of Q0 for critical data were carried out at various values of Γ.These results are plotted in Fig. 4.14. We see that the initial stages of the evolutionproceed with the convergence to the Hunter-A solutions. The higher the value ofΓ, the higher its instability as implied by the growth rate (Lyapunov exponent) ofthe unstable mode. This means that generic initial data, fine-tuned as closely asour numerical precision permits (|p− p⋆|/p⋆ ∼ 10−15) does not have enough time to“shed away” the stable mode dependence that defined the initial state. Therefore,we see only a partial convergence to its Hunter-A solution.We computed the profiles for α(x) and u(x) at various values of Γ from the11Using double precision arithmetic, the minimum difference (p− p⋆) which could be resolved is10−15 for a control parameter p that is of order 1. If quadruple precision were employed we couldresolve differences down to 10−31.1274.2. Numerical SimulationsFigure 4.14: Evolution of Q0 versus the central density ρc for critical collapse atmultiple values of Γ < 6/5. All the presented evolutions in this plot correspond tocritically collapsing data (|p− p⋆|/p⋆ → 10−14). As Γ is increased towards the valueof 6/5 the Hunter-A critical solutions becomes more unstable due to the growthof its unstable mode. The calculations are in agreement with this prediction. Ofthe presented cases Γ ∈ {1.04, 1.08, 1.12}, Γ = 1.12 shows the least amount ofconvergence to the Hunter-A solution (dotted line). All cases show convergence tothe value of Q0 corresponding the Larson-Penston (dashed line) solution at the laterstages of the evolution.numerical solutions, ρ(r, t) and v(r, t). We present the results of the solutions’maximum approach to the Hunter-A solution as indicated by plots of the variableα(x) in Fig. 4.15. Again, we observed a decrease in the degree of convergence at1284.2. Numerical Simulationshigher Γ values. As already suggested by our calculations of Q0(t) presented inFig. 4.14 the late stages of critical-initial-data evolution are well described by theircorresponding LP solutions. This fact is also evident in Fig. 4.16 where we observedthat the evolution data for α(x) closely matches the profile of the LP solution. Latetime convergence to the LP solution is unaffected at Γ = 6/5, this is expected sinceit continues to be stable until Γ reaches 4/3.The scaling exponent for the collapsed mass is predicted to be ((4− 3Γ)/λ0) asderived in Eq. (2.121). As Γ approaches 6/5, we know that λ0 grows nonlinearly,leading eventually to the vanishing of the scaling exponent. We measured the col-lapsed masses for Γ ∈ (1, 6/5] near their respective collapsed threshold (p⋆). Themass scaling is consistent with the expected behavior of the 1-mode unstable criticalsolution, Fig. 4.17. However, as we mentioned, the large value of the growth rateλ0, along with our precision limitation prevents the critical solution from a closeapproach to the Hunter-A solution. This however, does not change the fact thatit is an intermediate attractor of the evolution. The larger value of the Lyapunovexponent means the solution does not have sufficient time to drive away the stablemodes and the solution only marginally enters the linear regime.We can manually remove the stable modes by choosing initial data close to theHunter-A solution. In this case, the initial data is represented by the Hunter-Asolution plus some generic background perturbation. Again, the perturbation iscontrolled by a single parameter (labeled, once again by p) which can be fine-tunedto the threshold of collapse (labeled by p⋆). Notice that this p⋆ is not equal to zerosince that would mean we have provided the exact Hunter-A solution. Instead p⋆is a nontrivial value found as before through fine-tuning; therefore, there is still ageneric aspect to the solutions deviations from the Hunter-A solution at the initialtime. Under these circumstances, the initial data should be strongly attracted tothe Hunter-A solution, which should be well described the linear regime, Eq. (2.113)(with ~q = 0). This is particularly evident in the scaling of the collapsed mass when1294.2. Numerical SimulationsFigure 4.15: Plots of the dimensionless density variable α(x) emerging from criticalinitial data (p ≈ p⋆) for Γ = 1.04, 1.08, 1.12, 1.16 taken at intermediate times.The dotted line represents the respective Hunter-A (HA) solutions obtained viaintegration of the autonomous system, Eqs. (2.52) and (2.53). While the blue cure isthe dynamical solution for critical evolution i.e. |p−p⋆|/p⋆ → 10−15. The evolutionof α(x) approaches the Hunter-A solution, near the center, at intermediate times.The convergence becomes more tenuous at larger values of Γ.subjected to this type of initial data. We performed experiments for the same valuesof Γ as previously studied and we found very good agreement with perturbationtheory. These results are found in Fig. 4.18. The mass follows very closely thescaling law, given in Eq. 2.121 (F (~δ = 0) = 1) derived by assuming the critical1304.2. Numerical SimulationsFigure 4.16: Plots of the dimensionless density variable α(x) emerging from crit-ical initial data (p ≈ p⋆) for Γ = 1.04, 1.08, 1.12, 1.16 taken at late times. Thedashed and blue line represent the Larson-Penston (LP) and dynamical solutionsrespectively. The evolution data is taken late in the evolution when the unstableperturbation has grown and the solution has been driven away from the Hunter-Asolution. We clearly see that in all cases the solution resembles the Γ-dependentLarson-Penston solution which lacks any spherically symmetric growing modes andtherefore is stable.solution reaches the linear regime.At Γ = 6/5 supercritical evolutions of generic initial data (Model-A) near thethreshold yielded approximately the same value of the collapsed mass, see Fig. 4.17.This indicates the development of a mass gap, a fact which is consistent with the1314.2. Numerical Simulationsscaling exponent going to zero. From this and the fact that there are no 1-modeunstable solutions for Γ ≥ 6/5 we conclude that type-II critical behavior for thisfluid model ends at this value of Γ. The transition into type-I behavior is discussedin the next section.Figure 4.17: Scaling laws for the collapsed mass at Γ = 1.04, 1.08, 1.12, 1.16 and1.2. The data is in agreement with the scaling law for the mass derived from thesingle unstable mode of the critical solution Eq. (2.121), also plotted. As the Γapproaches 6/5, greater precision is required to get close enough to the 1-modeunstable solution. The slope of the line is related to ∼ 1/λ0, the reciprocal of theunstable mode, thus as λ0 diverges, the slope tends to zero. Consequently, a massgap develops in the spectrum of possible collapsed mass for Γ ≥ 6/5. The dynamicalresults agree with the prediction that this gap is generated at Γ = 6/5.1324.2. Numerical SimulationsFigure 4.18: Plots of the collapsed mass M with initial data resembling the Hunter-A solutions at Γ = 1.04, 1.08, 1.12, 1.16, respectively. Additional fine-tuning of theinitial data lead to the fading away of any deviations from Hunter-A. At late timesthe only the unstable mode grows, so the solutions is ensured to reach the linearregime where the scaling law for the mass of the core holds; this is evident in thisplot. Calculations of the core mass are in close agreement with the predicted scalinglaw and the respective value of the unstable mode’s growth rate. In this log-log plot,the slope formed by the series of collapsed cores becomes shallower at large Γ; thisimplies we need larger degree of fine-tuning to the critical parameter p⋆ to achieveinfinitesimal “collapsed cores”.1334.2. Numerical Simulations4.2.4 Simulations 6/5 ≤ Γ < 4/3Calculations of critical collapse solutions were also carried out at Γ ∈ [6/5, 4/3).The character of the critical collapse solution does in fact change for this range ofthe parameter Γ. Once more, we computed the quantity Q0(t) from the dynamicalsolutions. The late time evolution of this parameter indicates converge to its appro-priate Larson-Penston value, irrespective of fine-tuning. This is shown in Fig. 4.19.For fine-tuned data at intermediate times, the central density stops growing andbecomes static while Q0 continues to vary according to ∼ (t0 − t)2. The amount oftime the central density hovers around this temporary maximum depends on howclose the control parameter p in the initial data is a critical value p⋆. It is self-evidentthat in this case the central density is not describe by a self-similar solution i.e. onewhich “blows up” as (t0 − t)2 leading to a constant Q0. Instead, the intermediatecritical solutions seems to approach a static, star-like configuration, also shown inFig. 4.19.Snapshots of density, fluid velocity and specific entropy density provided inFig. 4.20 illustrate a convergence to a static solution. The density is plotted ona log-log scale, the other two fields namely, v and 3s/cv are on a linear-log graph.This evolution belongs to critical data (|p−p⋆|/p⋆ ∼ 10−14). The density at the cen-ter quickly grows as a result of the initial implosion. Once it has reached a certainvalue, the fluid velocity tends to zero over a region near the origin. This configu-ration remains virtually unchanged for a long time relative to its initial transientbehavior, the fluid then either collapses in a manner described by the LP solution(Fig. 4.19), or it disperses away from the center. The critical data presented inFig. 4.20 also shows the presence of shocks. The first outgoing shock forms almostinstantly at the start of the simulation. The second shock is generated followingthe initial compression of the fluid after which a shock wave is formed that travelsradially outward. Notice that as the fluid crosses the shock front by falling intothe compact object the specific entropy given in the form of 3s/cv increases discon-1344.2. Numerical Simulationstinuously. In regions where the flow is smooth the entropy field is simply advected(Fig. 4.20). This behavior is common to all critical initial data where 6/5 ≤ Γ < 4/3.Like the self-similar case the static solution is increasingly more difficult to resolvein the limit Γ → 6/5+. This parameter value implies a transition in the type ofcritical solution, a “boundary” separating static–similarity critical solutions.Figure 4.19: Plot of various measurements of Q0 versus ρc for Γ = 1.28. The fine-tuning (|p−p⋆|/p⋆ → 10−14) of the initial data do not yield in this case convergenceof Q0 to a constant value during intermediate times. Recall, Q0 is computed fromthe dynamical data (Equations (4.4)). Instead Q0 varies while the central densityreaches an intermediate maximum. The fluid evolution becomes nearly static at thismaximum density (ρ(0, t) ∼ 26). In the case of supercritical data, this is followedby rapid collapse. At late times Q0 converges to the value of the Larson-Penstonself-similar solution, regardless of fine-tuning.We cannot be certain that the static critical solution is characterized by a sin-gle unstable mode. In fact, the static solution cannot be determined from an au-tonomous system, that is by assuming a static, spherically symmetric ansatz. Thesetype of static solutions must obey the spherically symmetric Euler+Poisson equa-1354.2. Numerical SimulationsFigure 4.20: Snapshots of the critical evolution of model-A for Γ = 1.28. The profilesof density (ρ, black), fluid velocity (v, red) and entropy (3s/cv , blue, cv is the specificheat at constant volume) demonstrate convergence to a static intermediate state.The central density grows to a value of ∼ 26 where it remains for some time beforeultimately collapsing. Also evident in these plots is the presence of two outgoingshocks. The first shock develops almost instantly at the start of the evolution, whilethe second one forms following the contraction and subsequent expansion of the gasat the core’s surface. Notice that the fluid which crosses the shock by falling ontothe core acquires a jump in the entropy. In regions where the flow is smooth, theentropy is simply advected.tions (2.24)–(2.27) with v = 0,∂P∂r= −ρ∂ϕ∂r, (4.6)1r2∂∂r(r2∂ϕ∂r) = 4πGρ. (4.7)1364.2. Numerical SimulationsTogether with the polytropic ideal gas-law Eq. (2.18) we have three equations andfour unknowns, namely ρ(r), P (r), ϕ(r) and ǫ(r). The static state of the fluidis the result of the nonlinear evolution of the equations of motion and it cannotbe determined independently. This issue carries over to the linear perturbationanalysis, whereby the solution to the eigenvalue problem requires that we know thetime-independent solution of one of the field variables e.g. ǫ(r) and its perturbation.Nevertheless, the stability properties of the critical solution can be inferred fromthe numerical simulations. In particular, the life-time of the critical solution mustfollow Eq. (2.127) if it in fact contains a single growing mode. We measured this life-time by calculating the elapsed time the central maximum hovers about a commonvalue before collapse or dispersal. The results of these measurements are plottedagainst p − p⋆ in Fig. 4.21 for multiple values of Γ. We can clearly see the linearrelationship, from which the value of the unstable mode can be estimated. We notedthat the unstable mode decreases as we approach Γ = 4/3. Furthermore, in contrastto the self-similar critical solution, no evidence of universality was found in this case(Γ ∈ [6/5, 4/3)). The value of the unstable mode as well as the critical solutionsshowed dependence on the particular 1-parameter family being used. Initial datamodels A and B from Table 4.4 with Γ = 1.28 were tested. The results, plotted inFig. 4.22(b) support the absence of universality. The slope of the linear relationshipstated in Eq. (2.127) equals the reciprocal of the unstable mode λ0, so a difference inthe slopes imply different Lyapunov exponents (λ0). Also plotted in Fig. 4.22(a) arethe collapsed masses of the critical solutions. As stated earlier the critical solutionassociated with type-I critical phenomena is a finite, scale dependent solution (thescale set by the static solution). This scale is determine by the initial data and thecollapsed mass is proportional to this scale. As evident in Fig. 4.22(a), both sets ofdata exhibit a mass gap in the profile of the core’s mass that is clearly not universal.1374.2. Numerical SimulationsFigure 4.21: Scaling of the critical solution’s lifetime (T0) for various choices ofΓ > 6/5. The lifetime of this intermediate solution scales in a manner that isconsistent with the existence of a single unstable mode perturbing the intermediatestate. This is evident from its linear relation with log |p − p⋆|, as predicted inEq. (2.127). The Lyapunov exponent of the unstable mode determines the slope ofthe linear relationship. Steeper slopes imply smaller Lyapunov exponents which inturn imply more stable solutions. We can then note that as we Γ approaches the4/3 (the value corresponding to a gas of photons) the critical solution becomes morestable.1384.2. Numerical Simulations(a) (b)Figure 4.22: Plots of M(collapsed) (a) and the solution’s lifetime divided by thetotal mass (b) for two distinct 1-parameter families of initial data at Γ = 1.28. Onthe left panel we see the scale dependence of the critical solution is made evident bythe mass-gap, a characteristic of Type-I critical phenomena. The critical solution isa quasi-static configuration of finite size set by the choice of 1-parameter family ofinitial data. The lifetime of the critical solutions follows a scaling law, Eq. (2.127)as shown in (b). This scaling is a necessary property of type-I critical phenomena.The two distinct families of data presented follow slightly different scaling. The two1-parameter families, Models A and B yield different slopes i.e. different λ0s forthe linear relationship given by Eq. (2.127). This appears to support the absence ofuniversality for the unstable mode, λ0.139Chapter 5Results: Axial SymmetryThe number of numerical studies concerning critical phenomena in gravitationalcollapse beyond spherical symmetry remain limited. This situation is rather unfor-tunate given that angular momentum is expected to play a significant role duringthe late stages of collapse. As the length scale of the gravitation interactions short-ens, slow initial rotation will lead to high tangential velocities of the final compactobject. In order to investigate the role of angular momentum during critical collapsewe adopted axial symmetry for the geometry of the physical system. At the presenttime the only works which have produced numerical simulations of critical collapsebeyond spherical symmetry can be found in, [1, 2, 17, 18, 53, 97]. These pertain toaxisymmetric collapse of gravity waves and the scalar field (real and complex).We have implemented a numerical algorithm to solve the Euler equations of fluiddynamics coupled to Newtonian gravity in axial symmetry. Our aim was to investi-gate the role of angular momentum near the threshold of gravitational collapse. Thiswork was built upon the results of the spherically symmetric calculations (Chap. 4).We are particularly interested in the addition of slow (infinitesimal) initial rotationto otherwise spherically symmetric initial data near this threshold. The results ofthe numerical experiments also allowed us to verify the scaling laws predicted in[32], which we have rederived for our Newtonian system in Sec. DataInitial DataThree distinct 2-parameter families of initial data were chosen for the fluid evolutionwith axial symmetry. Rotation about the z-axis breaks spherical symmetry, whichcan be recovered by setting the velocity component vφ = 0. Otherwise, the threesets of data posses reflection symmetry about the equatorial plane (x − y plane),these are explicitly given in Table 5.1. We have chosen models A and B to coincidewith the spherically symmetric 1-parameter sets given in Table 4.4 if q = 0 (notethat this effectively sets vφ to zero). As we did in spherical symmetry, model A willbe the most commonly used set.In certain cases a fourth 2-parameter family of initial data was used. Similar towhat was done in spherical symmetry, this data represents the linear regime nearthe Hunter-A solution. The initial state is therefore represented by,Z(0, x, θ) = Z⋆(x) + fp(x, θ) + qZ1(x, θ). (5.1)As before, Z⋆(x) is the Hunter-A solution. The p-controlled perturbation expressedin cylindrical coordinates (s, φ, z) is given by,fp(s, z) ≡pe−s2−z2 + ερ(s, z)εvs(s, z)εvφ(s, z)εvz (s, z)εP (s, z). (5.2)Once more, the variables such as ερ(r, θ) are the cumulative errors generated inapproximating the Hunter-A solution. Note that we have introduced a Gaussianperturbation about the background Hunter-A solution of the density profile. Ro-tation is introduced by initializing the azimuthal component of the velocity fieldvφ according to the ℓ = 1 spin-up mode Z1(x, θ). In cylindrical coordinates and1415.1. The Unstable Axial (Spin-up) Modeassuming axial symmetry (m = 0), we haveqZ1(s, z) = q00δuΦxs√s2 + z200. (5.3)This expression is equivalent to Eq. (2.90) for ℓ = 1 andm = 0. The function δuΦ(x)give the radial dependence of the spin-up mode whose formal solution is given inEq. (2.104). Notice that x ∝ r/t2−Γ and so, x ∝ √s2 + z2/t2−Γ.Table 5.1: Initial data profiles for the primitive variables ρ, vs, vφ, vzand P used in the axisymmetric evolutions.Variable Model A Model B Model Cρ(0, s, z) e−s2−z2 1(1 + (s2 + z2)2)2cos3(π√s2 + z24)if√s2 + z2 < 210−12 if√s2 + z2 ≥ 2vs(0, s, z) 0 0 0vφ(0, s, z) qse−s2 qs√1 + s2{q sin(πs4)if√s2 + z2 < 20 if√s2 + z2 ≥ 2vz(0, s, z) 0 0 0P (0, s, z) pe−s2−z2 p(1 + (s2 + z2)2)2p cos3(π√s2 + z24)if√s2 + z2 < 210−12 if√s2 + z2 ≥ 25.1 The Unstable Axial (Spin-up) ModeSection 2.5 highlighted the relevance of the axial perturbation modes. Originallydetermined in the analysis of Hanawa [38, 40] the axial mode spectrum has a growthrate given by expression (2.103) provided regularity (analyticity) of the similarity1425.1. The Unstable Axial (Spin-up) Modesolution. As discussed in Sec. 2.5 the ℓ = 1 axial mode is the only relevant (growing)axial perturbation mode which can be observed through axisymmetric evolution ofour fluid model. Recall that for 1 < Γ / 1.17 the ℓ = 2 mode is also unstable,in axial symmetry only the m = 0 mode is observable but since our initial data(Table 5.1) is by construction symmetric about the equatorial plane the ℓ = 2 modeis suppressed during the evolution. Therefore, provided that 1 ≤ Γ / 1.17 thereis only one unstable axial mode that we predicted to be observable, namely thespin-up mode. This mode has a growth rate given by λ1 = 1/3. Given its unstablenature and our axisymmetric setup, the spin-up mode is expected to determine theangular momentum of the collapsed core when slowly rotating, supercritical initialconditions are considered. In general, there is no a priori reason to rule out polarperturbations, nevertheless these have been left out of our analysis due to indirectsuggestions of their limited relevance, these are mentioned in Sec. 2.5 and Sec. 2.6.3.We now present our calculations of the spin-up mode for the Hunter-A solutionat various values of Γ. Recall that a formal solution is provided by Eq. (2.104) fromwhich it is clear that the spin-up mode couples only to the radial velocity profile ofthe similarity solution u(x). So we proceeded by solving Eq. (2.100) in simultaneitywith Eqs. (2.52) and (2.53). These results are plotted in Fig. 5.1. We plotted theradial part of the spin-up mode given by δuΦ/x. These profiles contain a self-similar“hump”, its amplitude decreases with increasing Γ as shown in Fig. 5.1. All of theseconverge to the same behavior at small x, as illustrated in the inset of Fig. 5.1.To see the full structure of the spin-up mode in three dimensions we must com-plement the radial function just described with the mode’s angular dependence.This is given by,δuφ(τ, x, φ, θ) ∼ δuΦx∂∂θY 01 (θ, φ) =12√3πδuΦxcos θ.Note that only the m = 0 term is relevant in the context of our axisymmetric fluidmodel with the symmetry axis aligned with the z-axis. Fig. 5.2 presents the angular1435.1. The Unstable Axial (Spin-up) ModeFigure 5.1: Radial function plot of the spin-up mode (δuΦ(x)/x) for the Hunter-Asolution at five values of Γ. The angular part can be obtained by the vector sphericalharmonics according to Eq. (2.90). All of these functions have a similar profile, withthe peak of δuΦ/x that slightly decreases with increasing Γ. In all these calculationswe chose the free parameter δuΦ0 = 1. A log-log plot of these profiles is providedin the inset, where we can see that they all have the same behavior near the origin,which follows δuΦ(x)/x = δuΦ0x.dependence of the spin-up mode for the Hunter-A solution in the cases, Γ = 1.00001and Γ = 1.12. A similar spatial domain was chosen in both cases. The time prior tosingularity formation (t0− t) in Eq. (2.91) was chosen independently in each case sothat the profiles’ maxima are of similar order (∼ 10−2). The full ℓ = 1 axial-modespatial solutions for Γ = 1.00001, 1.12 are plotted in Fig. 5.2(a) and Fig. 5.2(b),respectively. The angular dependence is identical for all values of Γ, however due1445.1. The Unstable Axial (Spin-up) Modeto the small difference in the radial dependence illustrated in Fig. 5.1 the profilesdiffer slightly with increasing Γ.(a) (b)Figure 5.2: Plots of the spin-up mode’s azimuthal-component velocity field vφ forΓ = 1.00001 and Γ = 1.12. Panels (a) and (b) present only the ℓ = 1 contributionof the vector harmonic function ~Φ0ℓ(θ, φ) at Γ = 1.12 and Γ = 1.00001 respectively.Since this is the only observable axial mode which has a positive growth rate asindicated by Eq. (2.103) and the symmetry of our fluid model, slowly rotating initialdata near the threshold of gravitational collapse vφ is predicted to resemble thisprofile at late times.1455.2. Numerical Solutions in Axial Symmetry5.2 Numerical Solutions in Axial SymmetryThe bulk of our calculations in axial symmetry were conducted at Γ = 1.00001. Thisvalue of the adiabatic index (Γ) yielded results which closely resemble previous cal-culations obtained using the spherically symmetric isothermal gas model [44–48, 64].The results of the isothermal-gas coupled to Newtonian-gravity are significant sincethey represent a particular limit (k → 0 in Eq. (1.3)) of the spherically symmetricGeneral Relativistic self-similar perfect fluid solutions [77]. This correspondence be-tween the GR perfect fluid and the Newtonian isothermal gas motivated our choice(Γ = 1.00001) for the polytropic index. It seems reasonable to suppose that Gund-lach’s work on non-spherical linear perturbations of the same General Relativisticsystem [10, 24, 30, 32, 33] also includes the Newtonian isothermal gas limit. Fromthis work we have made the observation that the non-spherical perturbation modestructure is preserved as we look at the Newtonian limit discussed in [77]. Morespecifically, the unstable axial mode found by [33] whose growth is given by (2.105)becomes the spin-up mode under the Newtonian limit. Based on this observationwe predicted that the introduction of nonspherical initial data characterized by slowrotations into our Newtonian axisymmetric model will be analogous to the predictedbehavior of the General Relativistic system.We argue that the situation is exactly analogous in the case of critical collapse.We know of this exact analogy in the cases of spherical symmetry, Chap. 4. In thenonspherical case (slow rotation) the critical solution or more specifically its depar-ture from the Hunter-A solution is governed by two growing modes, thus maintainingits exact analogy with the general relativistic system [30, 32, 33]. In particular therole that non-zero angular momentum plays in both systems should likewise beanalogous. It is precisely the effect of non-zero angular momentum on the criticalsolution that we set out to explore.1465.2. Numerical Solutions in Axial Symmetry5.2.1 Slow Rotation (q → 0)We solved Euler’s equations (2.34)–(2.39) using the finite volume numerical tech-niques outlined in Chap. 3. We then conducted numerical experiments using the2-parameter families of initial data given in Table 5.1. The degree of initial rotationis controlled by the parameter q. In the first set of experiments q is set to a verysmall relative value. Specifically, q = 10−14, and since q enters the initial conditionsas a factor in vφ, this azimuthal velocity field is many order of magnitude smallerthan vs and vz after the first time step in the evolution. This scenario is chosen inorder to test the prediction that the spherically symmetric results are obtained inthe limit q → 0. We prepared critical initial data with parameter p near p⋆, wherep⋆ is a factor of the pressure profile that sets the threshold of gravitational collapsein the absence of rotation (q = 0). Turning on the initial rotation by a very smallamount for instance by setting q = 10−14 does not change the threshold value p⋆12.In this sense we consider this the addition of infinitesimal rotation.The first quantity we checked for was the convergence of the collapsing center tothe Hunter-A solution. Following what was done in spherical symmetry, we computethe self-similar quantity Q0 from the dynamical solution for the density at the originρ(t, 0, 0). We found precisely the same convergence of Q0(t) to the Hunter-A solu-tion that we observed in spherical symmetry. All three families exhibit convergenceto the same value of Q0. This fact is presented in Fig. 5.3. As was expected, thediscrepancies which characterizes the differences in the initial data sets are “washedaway” by their common convergence to the intermediate attractor (the Hunter-Asolution), hence the observed universality of the critical solution. Finally, these crit-ical fluid evolutions proceed with their convergence to the Larson-Penston solution.Due to the high computational costs we are unable to follow this convergence to thesame degree possible in spherical symmetry. Eventually, this convergence will also12There is an inherent uncertainty in the computed value of p⋆ due to limitations in numericalprecision. Therefore, the change introduced by adding small (infinitesimal) q = 10−14 is smallerthan this uncertainty.1475.2. Numerical Solutions in Axial Symmetrybe interrupted by the spin-up mode of the Larson-Penston solution. Also containedin Fig. 5.3 is the spherically symmetric critical data corresponding to Model A withΓ = 1.00001.Figure 5.3: Plot of Q0 versus central density (ρ(t, 0, 0)) for critical initial data. Thevalue of Q0 was computed from the critical solution using Q0(t) = 4πGρ(0, 0, t)(t0−t)2 for the three distinct 2-parameter families of initial data presented in Table 5.1.The degree of initial rotation was controlled by the parameter q which was setto 10−14 in all of the three cases. We note that as the central density grows thecomputed value of Q0 converges to the value corresponding to the Hunter-A solution,given in Table. 4.2. After some time the convergence to Hunter-A breaks down andthe solution begins to approach the Larson-Penston (L-P) solution. The results ofthe spherically symmetric critical solution are also provided for comparison purposes.The initial degree of rotation can be made small enough that we reproduced thesame spherically symmetry results. Nevertheless, as soon as q is nonzero the linearregime (equation (2.110)) contains an extra unstable mode. From the discussionof Sec. 2.8.1 we should observe an unstable nonspherical axial mode during theevolution of the rotating fluid. This mode enters our formalism through the velocitycomponent vφ, and so it is through monitoring of this velocity field component thatwe were able to detect it. In essence, if this 2-unstable-mode picture is correct thevelocity field vφ should resemble the profiles corresponding to the spin-up modegiven in Fig. 5.2.Indeed, snapshots of vφ at its closest approach to the Hunter-A solution for1485.2. Numerical Solutions in Axial Symmetryall three families A-C displayed in Fig. 5.4(a)–Fig. 5.4(c) reproduce the featuresof the ℓ = 1 spin-up mode plotted in Fig. 5.4(d) (this is a plot of the azimuthalvelocity field computed from the spin-up mode). The time of closest approach isdetermined by the collapsing center, Fig.5.3. Since the spin-up mode is presumablythe only nonspherical growing structure (with all others decaying), the initiallydistinct profiles are expected to become similar. The differences in the observedscales between the calculation of vφ from the critical evolutions and that computedfrom the spin-up mode can be accounted by the self-similar radial part, namelyδuΦ(x)/x. Fig. 5.4(d) represents a plot of vφ taken at a significantly time before thecollapse occurs, so the structures of the spin-up mode are not as “compressed” asthey are in panels (a)–(c) of Fig. 5.4. The peak in these profiles is determined byδuΦ(x)/x and an overall family-dependent scale set by κ as well as t0 − t (t0 is thetime of collapse), see Eq. (2.91).The data presented thus far is consistent with the prediction that a growingaxial perturbation is solely responsible for the solution’s final angular dependence.Critical initial data (p ≈ p⋆) leads to the known funneling of the solution to the two-mode linear regime described in Sec. 2.5 and given by Eq. (2.110). The interplay ofthese two modes yields a modified version of Choptuik’s mass scaling law Eq. (2.118).The specific angular momentum is also expected to follow similar scaling behaviorthis is given by Eq. (2.120). The case under current discussion is the limit ofinfinitesimal initial rotation. The mass and specific angular momentum scalinglaws in this limit are respectively given by M ∝ |p − p⋆|(4−3Γ)/λ0 and a ∝ q|p −p⋆|(3−2Γ−λ1)/λ0 (Eqs. (2.123) and (2.124)).The scaling law for the specific angular momentum in the limit as q → 0 has anexplicit dependence on the growth rate of the axial perturbations λ1, Eq. (2.124).We computed the collapsed massM and its specific angular momentum a for a seriesof supercritical runs. First of all, we confirmed that the results are in agreement withspherically symmetric data. This is expected since the mass scaling is independent1495.2. Numerical Solutions in Axial Symmetryof q, to leading order, i.e. ∂F (δ)/∂q|q=0 = 0. All three 2-parameter families ofinitial data yielded the same mass scaling behavior as evident in Fig. 5.5. The massscaling independently computed in spherical symmetry produced the same results.The angular momentum scaling behavior also displays universality for increasinglyfine-tuned data (p → p⋆). More importantly, the scaling behavior is consistent thepredictions of Eqs. (2.123) and (2.124). Indeed, the measured angular momentumof the intermediate state has the imprint of the axial growing mode. Supplying theknown values, λ0 ≈ 9.4643 (from Table 4.3) and λ1 = 1/3, we can compute theperturbation theory prediction. The computed data converges to this prediction asp→ p⋆ (Fig. 5.5).Choosing initial data near the Hunter-A solution allowed us to identify the in-fluence of the two growing modes on the scaling behavior of the collapsed mass andits angular momentum. Plotted in Fig. 5.6 are the measurements of quantities M(collapsed mass) and a (collapsed mass’s specific angular momentum) for a seriesof supercritical evolutions near the collapsed threshold with initial data that ap-proximates the Hunter-A solution, Eqs. (5.1)–(5.3). Notice that with this type ofinitial data we have already, from the outset “trimmed out” most of the decayingperturbations. We conjectured that a very clear signal of the two growing modesshould be observed through the scaling of M and a. Again, we also plotted thepredictions from perturbation theory (Fig. 5.6). Clearly, our measurements of Mand a closely match the predictions coming from perturbation theory.1505.2. Numerical Solutions in Axial Symmetry(a) (b) (c)(d)Figure 5.4: Critical-evolution measurements of vφ for the spin-up mode, with Γ =1.00001. Panels (a) through (c) are the result of critical evolution of initial datacorresponding to models A, B and C from Table 5.1 taken at its closest approachthe to the Hunter-A solution. The initial rotation is controlled by the parameter qwhich was set to 10−14 in the presented cases. Panel (d) is the result of our explicitcomputation of the spin-up mode. Recall that the spin-up mode and its angulardependence can be computed explicitly, the angular part is given by the ℓ = 1vector harmonic, i.e. by ∂Y 01 (θ, φ)/∂θ. All three initial data sets (A-C) converge tothe same profile; one where the only growing structure is described by the spin-upmode (d), thus indicating universality. The scales differences among panels (a-c)and (d) is accounted by the self-similar nature of the radial profile i.e. δuΦ(x)/x.Panel (d) is generated at a significant time before the collapse happens. The size andlocation of the peak is determined by the radial profile (δuΦ(x)/x) and an overallscale set by the local speed of sound, see Eq.(2.91).1515.2. Numerical Solutions in Axial SymmetryFigure 5.5: The scaling behavior of the collapsed mass (M) and its specific angularmomentum (a) for supercritical initial data near the collapse threshold p⋆. The casespresented here belong to the three distinct 2-parameter families of initial data givenin Table 5.1 with q = 10−14. The perturbation analysis developed by Gundlach[32], discussed in Sec. 2.8.1 provided predictions for the scaling behavior of the massand angular momentum of the compact object near the collapse threshold. Ourcalculations are in agreement with Gundlach’s perturbation theory predictions inthis case of slow, initial rotation. All three families converge to the predicted linearrelationship suggesting universality. This is explained by convergence to a commonscale-invariant state (the Hunter-A solution). Note that in general a should be avector, but in axisymmetry only the z-component is non-zero.1525.2. Numerical Solutions in Axial SymmetryFigure 5.6: Scaling behavior of the collapsed mass M , and the specific angular mo-mentum a, for initial data that closely resembles the Hunter-A solution, Eqs (5.1)–(5.3). The initial data was endowed with an infinitesimal amount of rotation, setby q = 10−14. The infinitesimal rotation introduces a non-spherical unstable axialmode (spin-up) with growth rate λ1 = 1/3. The existence of this non-sphericalgrowing mode leads to a scaling-law for a, this is given by Eq. (2.124). The scalingof M follows Eq. (2.123). Our calculations of M and a during supercritical evolu-tions which are now very close to the Hunter-A solutions behave in accordance withthe predictions of perturbation theory. More specifically, the slopes of the linearrelationships are consistent with the values λ0 = 9.46430101 and λ1 = 1/3, thegrowth rates of the spherical and axial modes respectively.1535.2. Numerical Solutions in Axial Symmetry5.2.2 Finite Initial RotationIn the previous section we considered the evolution of our fluid model subject toinfinitesimal initial rotation. The next step in our investigation was to add a finiteamount of initial angular momentum to the initial state in order to measure thefunctions F (δ) and G(δ) in Eqs. (2.118) and (2.120), which are presumably universal.The degree of rotation is controlled by q; in this case it is a finite, yet small quantitysuch that the critical solution still goes through the linear regime given by theexpansion (2.113). Prior to presenting our calculations of F (δ) and G(δ), we wouldlike to address the question as to what happens to critical initial data when smallbut finite rotation is included.We wished to know the fate of initial data near the threshold p⋆ upon the additionof a finite amount of initial rotation. In [32] Gundlach discussed two possible results.Initial data, Z⋆ + δZ1 which teeter on the brink of gravitational collapse, eitherformed a compact object followed by complete collapse, or it dispersed, leavingbehind empty space. In his article [32], he labels these as follows,• Possibility 1: Critical initial data collapses even with the addition of angularmomentum through δZ1.• Possibility 2: Intuitively, the addition of angular momentum should gener-ate an outward “centrifugal” pressure to the critical state that leads to thedispersal of the gas.The experiment was conducted using our Newtonian model and the results obtainedcorrespond to Possibility 2 discussed in Sec. VII of [32]. Possibility-2 as Gundlachpointed out is the more physically intuitive possibility, we expect that the additionof angular momentum promotes the outward dispersal of the fluid. This is indeedwhat we found in our numerical experiments, adding rotation to already criticaldata leads unambiguously to fluid dispersal.The universal functions depend on the initial parameters via δ ∝ q/|p−p⋆|λ1/λ0 .1545.2. Numerical Solutions in Axial SymmetryWe varied δ by fixing p and choosing a series of values for q starting from zero. Wechose p in the supercritical regime, so that at q = 0, our choice of p yields collapsingdata. Nevertheless, p is near the threshold value p⋆. Note that p⋆ represents thethreshold of collapse only if q = 0. As we have just discussed, finite q shifts thisthreshold. We can represent this as pcr(q), with the condition that pcr(0) ≡ p⋆13. The calculated values of F (δ) and G(δ) are plotted in Fig. 5.7 and Fig. 5.8,respectively. These were calculated from the measured values ofM , a and the scalinglaws (2.118) 2.120. Three 2-parameter families were tested, with the normalizationconditionsF (0) = 1,∂G(0)∂δ= 1, (5.4)which follows from the freedom to normalize the mode functions Z0(x) and Z1(x, θ, φ).In essence, this amounts to dividing out the family-dependent constants C¯0(4−3Γ)/λ0 ,C1C¯0(3−2Γ−λ1)/λ0 in Eqs. (2.118) (2.120) for q → 0 data. We found evidence to sup-port the universality of these function since all three families yielded nearly identicalresults.More detailed measurements of F (δ) and G(δ) were generated at the large δregime. We know already that the addition of initial angular momentum shiftsthe threshold of gravitational collapse (Possibility 2). We began with supercriticaldata at fixed p < p⋆, then we added just enough amount of initial rotation bymodulating q until we find the threshold of collapse, we labeled this qmax. Thethreshold can be defined in terms of the quantity δmax. Where δmax = qmax/|p −p⋆|λ1/λ0 . Alternatively, we could start with fixed q and modulate p until p = pcr, sothat, δmax = q/|pcr − p⋆|λ1/λ0 . By definition, both F (δ) and G(δ) must vanish for|δ| > δmax. The universal functions are non-trivial over the region of the p-q planedefined by,C¯0(p− p⋆) < 0, and C1q < δmax(C¯0|p− p⋆|)λ1/λ0 . (5.5)13This definition was also borrowed from Ref. [32]1555.2. Numerical Solutions in Axial SymmetryFigure 5.7: Plot of function F (δ) for three different 2-parameter families of initialdata. These 2-parameter families are specified in Table 5.1. This plot confirmsspeculations concerning the universality of F (δ) stated in the perturbation theoryanalysis of [32]. As the initial rotation is turned off (δ → 0), the function F (0)→ 1for the three families (after appropriate rescaling), a fact which is consistent withthe spherically symmetric results. This data belongs to the case where Γ = 1.00001in the polytropic EoS (2.18).Gundlach [32] speculated on the possible nature of these functions near the thresholdδmax. Assuming that the situation described by possibility 2 is realized (as indeedit was in our case), Gundlach provides two subcategories for the behavior of theorder-parameters M and a near the collapse threshold given by δmax or pcr. Theseare,1565.2. Numerical Solutions in Axial SymmetryFigure 5.8: Plot of function G(δ) for three different 2-parameter families of initialdata. The three 2-parameter families of initial data can be found in Table 5.1. Ourthree measurements of G(δ) corresponding to the initial data families (A-C) arenearly identical; this is in agreement with the predicted universality of G(δ) [32].We obtained the expected linear behavior of G(δ) near δ = 0 (inset).• Possibility 2a: A mass-gap develops in the collapsed mass spectrum across theδmax boundary. This implies a discontinuity in the universal functions.• Possibility 2b: There is no mass-gap; for this possibility Gundlach proposed apower law behavior for the universal functions which must vanish across the1575.2. Numerical Solutions in Axial Symmetrythreshold boundary, e.g.,F (δ) = KM (δmax − δ)βM δ . δmax0, δ > δmax(5.6)where βM is again a universal exponent and KM is a family dependent con-stant. Similarly, for G(δ),G(δ) = Ka(δmax − δ)βa δ . δmax0. δ > δmax(5.7)Our results indicate continuity of F (δ) and G(δ) across the δmax boundary, in agree-ment with Gundlach’s Possibility 2b also discussed in Sec.VII of Ref. [32]. Thecalculations agree with the proposed behavior near δmax as dictated by Eqs. (5.6)and (5.7). These results are plotted in Fig. 5.9. Both F (δ) and G(δ) seem to vanishas δ → δmax according to an power law. The results for the three families appear tobe governed by a scaling law whereby the scaling exponents, βM and βa are givenby,βM = (4− 3Γ)/λ0 and βa = (3− 2Γ− λ1)/λ0, (5.8)the same as those of scaling laws for M and a given in Eqs. (2.118) and (2.120).This fact is illustrated in Fig. 5.9. If the exponents for the assumed power-laws ofF (δ) and G(δ) Eqs. (5.6) and (5.7) near δmax are indeed those given in Eq. (5.8)then, it is easy to show that the mass M and specific angular momentum a of thecollapsed core can be written as,M ∝ (qmax − q)(4−3Γ)/λ0 and a ∝ (qmax − q)(3−2Γ−λ1)/λ0 , (5.9)provided that δmax is found by fine tuning q to qmax. If δmax is obtained by varyingp → pcr then in order to see the behavior of M and a near pcr we need to performan expansion of F (δ) and G(δ) in the small quantity |pcr − p|. These requires a bit1585.2. Numerical Solutions in Axial Symmetryof algebra to show that,M ∝∼ (pcr − p)(4−3Γ)/λ0 and a ∝∼ (pcr − p)(3−2Γ−λ1)/λ0 . (5.10)Where the symbol ∝∼ indicates the presence of higher order terms of the quantity|pcr − p| that are omitted in Eq. (5.10).(a) F (δ) (b) G(δ)Figure 5.9: Universal functions F (δ) and G(δ) near the collapse threshold δmax.Calculations were performed using three distinct 2-parameter families of initial data.The results indicate a vanishing of F (δ) and G(δ) as δ → δmax according to a powerlaw ∼ (δmax − δ)β . The scaling exponents for F (δ) and G(δ), respectively, βM andβa seem to be the same as those which govern the scaling of M and a in equations(2.118) and (2.120)From the above discussion we have learned that, whether we approach the col-lapse threshold by varying p→ pcr at fixed q or q → qmax at fixed p, M and a havethe same scaling behavior that depend only on the “distance” from the collapsethreshold. The numerically calculated form of the universal functions F and G in-dicate continuity of the spectrum of collapsed cores across the threshold boundary.There appears to be no mass gap in the spectrum even with finite initial angular1595.2. Numerical Solutions in Axial Symmetrymomentum Fig. 5.9.A parameter survey was conducted over the p − q space. For this survey weselected Model-A from Table 5.1 to define the initial state. The parameter domainwas chosen to be a rectangular grid given by (p, q) ∈ [0.99p⋆, 1.01p⋆] × [−0.1, 0.1].We then proceeded to find all the pcrs at every q on the grid by way of multiplebinary searches 14. Similarly, we found the qmaxs at every p value. This procedureallowed us to trace the critical collapse curve (pcr, qmax). We can clearly establishthat the critical curve follows a parabola in the p− q space as displayed in Fig. 5.10.The critical parameters pcr and qmax are modeled by pcr(q) = p⋆ − Aqb, qmax(p) =A˜(p−p⋆)b˜ respectively. The scalar parameters A or A˜ are dependent on the family ofinitial data, whereas b or b˜ are universal parameters set by the relative dimensionalityof p and q. A least-square curve fit was applied to this data in order to approximatethese curve parameters. We discovered that for this data A ≈ 0.114 and b ≈ 1.99, avalue which is suspiciously close to 2.Calculations for M and a near the critical curve were generated by approachingthe curve along the two independent direction on the p − q plane. From thesecalculations we were able to ascertain that the behavior of M and a is indeedthat given by equations (5.9) and (5.10). Much like the spherically symmetric casearbitrarily small cores can be produced with increased proximity to the criticalcurve. The survey results are plotted in Figs. 5.11 and 5.12. The quantities M anda seem to be continuous across the critical curve. Furthermore, the scaling of Mand a is identical (universal) anywhere near this curve and only seems to dependon its distance (in the p − q plane) from it. The results are also consistent withthe symmetry requirements M(p,−q) = M(p, q) and a(p,−q) = −a(p, q) as seen inFigs. 5.11 and 5.12.We used the results generated by our parameter survey to compare with Gund-14A binary search refers to the process of finding the threshold of gravitational collapse by nu-merically evolving the equations of motion for subcritical and supercritical values of the controlparameter. The subcritical-supercritical gap is narrowed until the ‘critical’ value is identified.1605.2. Numerical Solutions in Axial SymmetryFigure 5.10: Curve of the gravitational collapse threshold on the p − q parame-ter space for initial data Model-A at Γ = 1.00001. Binary searches were carriedout along the two independent directions over this space to determine the collapsethreshold curve. This was done by fixing p < p⋆ and then fine tuning (via a binarysearch) q to qmax. As expected, we found that for the same value of p anothercollapse threshold could be identified at q = qmin, the value of qmin matches oursymmetry expectations i.e. qmin = −qmax. We plotted the experimentally deter-mined data (p, qmax) and (p, −qmax). Similarly, at fixed q the binary search (in p)yields pcr for the threshold parameter, naturally, pcr = p⋆ in the case where q = 0which corresponds to spherical symmetry. Again, we plotted (pcr, q). A power-lawcurve of the form, pcr(q) = p⋆ + Aqb was fitted to this data using the method ofleast-squares and determine the parameters A and b. For this data A ≈ 0.114 andb ≈ 1.99. The first parameter, A is family dependent, but b is a universal parameterthat depends on the relative dimensionality of p and q. The value of b (1.99) issuspiciously close to 2, thus the relation between p and q is a parabola.lach’s schematic plots of M and a found in Fig. 2 of Ref. [32]. The plots presentedthere are the consequence of considering the form of F (δ) and G(δ) to be that of1615.2. Numerical Solutions in Axial Symmetryequations (5.6) and (5.7) (Possibility 2b). These plots have the same general “topol-ogy” as our Figs. 5.11 and 5.12. This means continuity of the spectrum of collapsedmasses across p(qmax). In Gundlach’s paper the surface plots of M and a over thep − q space are constructed by assuming the modes’ growing rates are λ0 = 2 andλ1 = 1. However, for our Newtonian model the growing rates are, λ0 ≃ 9.4643, andλ1 = 1/3.Figure 5.11: Survey of collapsed core masses (M) in the p − q parameter space.These calculations were performed using initial data model A from Table 5.1, andthe adiabatic index was set at Γ = 1.00001. The threshold of gravitational collapsedwas found by determining all qmaxs and pcrs; this allowed us to calculate cores’masses near the threshold along the two independent directions in the p− q plane.The evidence indicates that arbitrarily small collapsed cores can be generated acrossthe threshold boundary, i.e. the is no mass-gap in the core’s spectrum. The behaviornear this boundary is the same over the presented parameter ranges and it is givenby (5.9) and (5.10). In the p − q plane the collapsed threshold seems to follow aparabolic curve with vertex set (p = p⋆, q = 0). These results display the expectedsymmetry about q = 0 and they look similar to those plotted in Ref. [32] Fig. 2.1625.2. Numerical Solutions in Axial SymmetryFigure 5.12: Survey of the collapsed core’s specific angular momentum a in thep− q parameter space. All the computations were conducted using the initial dataset model-A found in Table 5.1, with Γ = 1.00001. Careful binary searches werecarried along the two independent directions in the p − q plane, in order to findthe threshold values pcr at every value of q, and qmax at every value of p. Oncethe set {pcr, qmax} over the selected part of the p − q plane has been determinedwe proceeded to compute a near this threshold values. The angular momentum ofthe collapsed core seems to go continuously to zero in going from the supercriticalto the subcritical region. Furthermore, we confirmed that the behavior of a nearthe curve (pcr, qmax) is given by equations (5.9) and (5.10). The collected data fora also displays the expected anti-symmetry a(p,−q) = −a(p, q). These results aresimilar to those predicted in Ref. [32], Fig. 2.1635.2. Numerical Solutions in Axial Symmetry5.2.3 Rapid Initial Rotation (Large q Regime)Some of what has been said so far concerning the behavior of M and a near thecritical collapse solutions eventually breaks down in the large q regime. For initialdata described by Model-A in Table 5.1 the parameter q ∈ (0, 0.1] yield the resultsdescribed above. For values q & 0.1 the critical solutions looks qualitatively differ-ent. Concretely, the velocity field component vφ contains features that cannot beaccounted for, solely by the ℓ = 1 spin-up mode. Snapshots of vφ for the criticalevolutions of 2-parameter families A, B and C Table. 5.1 at q = 0.5 are presented inFig. 5.13. Clearly, we need to include higher order vectors harmonics to adequatelyaccount account for the angular dependence displayed by vφ, as evident in Fig. 5.13.The solution, however, does seem to continue to have features of self-similarity. Al-beit, the self-similarity displayed by the velocity field component vφ seems to be ofa discrete nature, at least with regard to the field vφ. This is evident in Fig. 5.13.This solution also seems to be a universal critical solution.Initial data with large values of q does not necessarily go through the lin-ear regime described in Sec. 2.8.1 from which the scaling laws for collapsed massEq. (2.118), and specific angular momentum Eq. (2.120) were derived. In fact, thefeatures of the profile for vφ (Fig. 5.13) suggest a departure from the linear regime.Nevertheless, we can calculate the scaling behavior of M and a near the thresholdof gravitational collapse. This results in a series of collapsed masses whose scalingbehavior near (p→ pcr) remains unchanged, even over this “extreme” regime. How-ever, the specific angular momentum clearly diverts from its slow rotation behavior.The scaling exponent of a near the threshold of collapsed pcr seems to be identicalto that of the collapsed mass. These results are plotted in Fig. 5.14. The scaling ofa is clearly different from what is predicted by Eq. (5.10). On a speculative note,the observed results are consistent with the scaling exponent βa = 3 − 2Γ insteadof βa = 3− 2Γ− λ1 (Eq. (5.8)) for the G(δ)-ansatz in Eq. (5.7). Interestingly, theseresults can be derived by assuming a different “linear regime”, where the critical so-1645.2. Numerical Solutions in Axial Symmetrylution approaches a one-mode unstable, nonspherical, self-similar solutions (clearlydifferent from Hunter-A). This regime could be represented by,Z(τ, x, θ, φ; p˜) ≃ Z⋆(x, θ, φ) + C0(p˜− p˜cr)Z0(x, θ, φ)eλ0τ . (5.11)Here Z⋆(x, θ, φ) is one-mode-unstable similarity solution to the axisymmetric isen-tropic gas15 with Z0(x, θ, φ) being its linear unstable mode function. In this case p˜is the generic parameter with non-trivial critical value p˜cr. Notice that we could usep, q, δ or any other parameter to represent the 1-parameter family of initial data.Repeating the formalism of Sec. 2.8.1 it is easy to show that eτ∗ ∝ (p˜ − p˜cr)1/λ0 ,where τ∗ represents the sole length scale in the system. From dimensional analysiswe know |~a| ∝ t2n−1∗ and so,|~a| ≡ a ∝ |p˜− p˜cr|(3−2Γ)/λ0 . (5.12)Our measurements of a near p˜cr are consistent with Eq. (5.12). This correspondsto the case where Γ = 1.00001 and λ0 = 9.4643 in which case (4 − 3Γ)/λ0 ≈(3− 2Γ)/λ0 ∼ 1/λ0. Thus the scaling of M and a are essentially indistinguishable,see Fig. 5.14.15In order to test this assumption we could proceed by assuming a similarity, axisymmetric ansatzto the Euler equations and then show that one of the self-similar solutions has a single unstablemode with growth rate ∼ λ0. Then, we would need to show that this one-mode solution is theintermediate attractor in critical collapse. This however, is beyond the scope of our current project.1655.2. Numerical Solutions in Axial Symmetry(a) (b) (c)Figure 5.13: Displays of vφ corresponding to the critical solutions for the case of“large” initial rotation (q = 0.5) with Γ = 1.00001. Panels (a)–(c) correspond toinitial data sets A–C respectively. The critical solution appears to be universal asthe different initial data sets A, B and C converge to a common solution. Evidentin panels (a)–(c) is the appearance of a type of discrete self-similarity in the profileof vφ.1665.2. Numerical Solutions in Axial SymmetryFigure 5.14: Measurements of M and a near the threshold of gravitational collapsepcr for data with “large” initial rotation (q = 0.5) with Γ = 1.00001. Panels (a)–(c)correspond to initial data sets A–C, respectively. The scaling of the massM remainsvirtually unchanged compared to the result obtained using slow rotating data. Thespecific angular momentum a scales differently at large q, its scaling law matches thebehavior of the collapsed mass, with scaling exponent, ∼ 1/λ0. One interpretationof the results is the existence of nonspherical one-mode unstable critical solutionwith Lyapunov exponent ∼ λ0. Then, by dimensional analysis, both the specificangular momentum a and the collapsed mass M have very similar scaling whenΓ = 1.00001, since, (4 − 3Γ)/λ0 ≈ (3 − 2Γ)/λ0 ∼ 1/λ0. This clearly differs fromthe data presented in Fig. 5.5. The critical solution is clearly universal as all threedata sets, A, B and C result in the same scaling near the threshold of gravitationalcollapse.1675.2. Numerical Solutions in Axial Symmetry5.2.4 Analogy with Statistical MechanicsSome of the most surprising and interesting elements of critical phenomena in grav-itational collapse stem from its analogy with elements of phase transitions in Statis-tical Mechanics. We will not delve into the nature of statistical mechanical systems;our purpose here is to simply outline the recognized parallels that it shares withthe self-gravitating systems. We also would like to describe the extension of thisanalogy resulting from the present work. For an introduction to phase transitionsin Statistical Mechanics see [109]. As pointed out in [10, 32] the formal calculationof critical exponents corresponding to a phase transition in a thermodynamic sys-tem is mathematically identical to the linear perturbation analysis that we used todetermine the growing mode perturbations. The critical solution is a fixed point ofa dynamical system. In statistical mechanics the fixed points contain two growingperturbation modes. One of them is related to temperature T . While the other isconnected to a generalized external force e.g. an external magnetic field ~Bext.One example commonly given is the liquid-gas phase transition. Provided otherthermodynamic variables are held constant this transition occurs at a critical tem-perature T ⋆. Allowing control over one other thermodynamic variable such as thepressure P provides the generalized external force, and thus a second growing mode.For this system the order parameter is the difference in density between the twostates i.e. ρliquid − ρgas. For T near T ⋆ the order parameter follows a power-law,ρliquid − ρgas ∝ |T − T ⋆|γ , (5.13)where γ is the scaling exponent, related to the growth rates of the perturbationmodes. The density interface changes continuously at the critical temperature T ⋆.This is an example of a second order phase transition, thus analogous to type-IIcritical phenomena in gravitational collapse.The other example which we have borrowed from [10, 32] follows from consider-1685.2. Numerical Solutions in Axial Symmetryation of a ferromagnetic material at high temperature. The order parameter is themagnetization ~m. The generalized force is the applied/external magnetic field ~Bext.At ~Bext = 0 and T < TC , where TC is the Curie Temperature, the material assumesspontaneous magnetization ~m. As temperature is increased to TC , |~m| vanishesaccording to|~m| ∝ (TC − T )γ . (5.14)However, a finite external magnetic field breaks rotational symmetry above the CurieTemperature; the magnetic moments align themselves with ~Bext. Hence the materialbecomes paramagnetic. This effect is attributed the to second growing mode aroundthe scale-invariant solution at T = TC , ~Bext = 0.With the confirmed existence of a second growing mode coming from the additionof angular momentum to the initial state in our gravitation collapse scenario, a nearexact analogy can be drawn with the ferromagnet plus ~Bext at high temperature.Gundlach and Goldenfeld [32] pointed out that the second growing mode wouldmake the two systems analogous provided the following associations are made16.The vector ~q which specifies the initial angular momentum plays the role of ~Bext.While the specific angular momentum of the collapsed core ~a plays to role of ~m. Thecollapsed massM is analogous to the correlation length ξ given that ξ ∝ (TC−T )−γ .Clearly, p is associated with T as is p⋆ with TC .Qualitatively different outcomes follow from initial data that straddle the col-lapse threshold. In spherical symmetry, Z⋆+εZ0, and Z⋆−εZ0 respectively generatethe dispersal and collapse outcomes. Beyond spherical symmetry the addition of an-gular momentum through nonzero ~q adds a second mode to the linear regime i.e.Z⋆±εZ0+~δ · ~Z1. In this case, the outcome (collapse/dispersal) depends solely on themagnitude of the ~δ. Evidence for this can seen in Fig. 5.11 and Fig. 5.12. Similarly,for the ferromagnet T = TC − ε and T = TC + ε correspond to two qualitatively16Their discussion pertained to the general relativistic fluid collapse, however, for reasons madeexplicit throughout this thesis their conclusions regarding the ferromagnet analogy applies likewiseto Newtonian collapse1695.2. Numerical Solutions in Axial Symmetrydistinct states. Again, the introduction of the applied magnetic field to the stateT = TC ± ε yields an outcome that is dependent on | ~Bext| but independent of itsdirection. The parameters ~q, and ~Bext both set a preferred direction in their respec-tive systems. The final angular momentum of the compact object (~a) is aligned withthe direction of ~q from the initial state. Similarly, the spontaneous magnetization ~mwhich occurs below the Curie temperature TC has no preferred direction until the~Bext break this symmetry and ~m always aligns with ~Bext. The symmetry breakingoccurs at the values ~q = 0 and ~Bext = 0, hence these are the critical values requiredto eliminate the second growing mode. The other critical values are, off course, p⋆and TC . In the latter case, the values cannot be trivially found but depend on otherdetails of the physical model.1705.3. Numerical Solution at Larger Values of Γ5.3 Numerical Solution at Larger Values of ΓFinally, we present the results for other selected values of the parameter Γ. We choseΓ within the two regimes known from our investigations in spherical symmetry togenerate both type-II and type-I critical behavior (Chap. 4). Once more, we inves-tigated the effect of adding an infinitesimal amount of initial rotation on evolutionof the fluid at the threshold of gravitational collapse.5.3.1 Critical Solutions Γ < 6/5In spherical symmetry 1-mode unstable solutions (Hunter-A) exist within the rangeΓ ∈ [1, 6/5). We know the mode’s Lyapunov exponent grows exponentially as Γ→6/5 (Table 4.3). Due to the this fact it becomes increasingly difficult for the evolutionof generic initial data to reach the linear regime. As a result, the initial data usedfor the dynamical evolutions at Γ = 1.12 was carefully constructed in accordance toEqs. (5.1)–(5.3). With infinitesimal amount of rotation (q = 10−14) we fine tunedp to the collapse threshold. We observed the same “enhanced” convergence to theHunter-A solution that we witnessed with Γ ≈ 1. The deviations from the Hunter-A solution are dominated by the two growing modes (spherical and axial). Weconcluded that for the range of Γ where a Hunter-A solution exists (Sec. 4.1.1) thecritical behavior is dominated by the two unstable modes. Thus, for Γ ∈ [1, 6/5)the critical solution is described by linear regime Eq. (2.117).This linear regime solution dictates the behavior of physical quantities such asM and a. As was done in the case of Γ = 1.00001 these quantities can be used toidentify the unstable modes through calculation of their scaling behavior. Shown inFig. 5.15(b), the scaling ofM and a follow precisely as dictated by the linear regime.Since our solution is already placed in the linear regime at the initial time, the scalinglaws follow closely the expectations from perturbation theory. The fine-tuning ofthis data further shows the “attractor” nature of the Hunter-A solution with onlythe growing structures of the two unstable modes. The growth of the central density1715.3. Numerical Solution at Larger Values of Γconverges to that predicted by the Hunter-A solution. This is evident in Fig. 5.15(a).Lastly, we computed the corresponding universal function F (δ) and G(δ) using aseries of supercritical data near p⋆ with q = 10−14. These are plotted in Figs. 5.15(c)and 5.15(d). We found no difference in the behavior of these functions which furthersupports our conclusion that this 2-mode structure of the critical solutions continuesup to but excluding the value of Γ = 6/5.1725.3. Numerical Solution at Larger Values of Γ(a) Q0(t) vs ρ(t, 0, 0) (b) Scaling of M and a(c) F (δ) (d) G(δ)Figure 5.15: Calculations in axi-symmetric critical fluid collapse with Γ = 1.12 andslow initial rotation (q = 10−14). The initial data used here is given by Eqs. (5.2)and (5.3). This data places the initial state near the Hunter-A solution. Therefore,we can see close convergence to appropriate Q0 value of Hunter-A solution as shownin 5.15(a). Series of calculations of M and a indicate their scaling behavior isgoverned by the same 2-unstable-mode linear regime. The spin-up mode has thesame growth rate λ1 = 1/3, but the spherical mode has a much larger growth rateof λ0 ∼ 37. The scaling laws of M and a are consistent with these exponents asshown in panel 5.15(b). The universal functions F (δ) and G(δ) were also computed(panels 5.15(c) and 5.15(d)). The results are formally identical to those obtained inthe near isothermal case (Γ = 1.00001).173Chapter 6Conclusion and Further WorkThis thesis was devoted to the investigation of the critical phenomena generatedby a Newtonian self-gravitating fluid at the threshold of gravitational collapse. Weextended the previous work on the spherically symmetric isothermal gas model byconsidering the more general polytropic gas. In the second phase of the our in-vestigation we have relaxed the symmetry restrictions and considered axisymmetricfluid collapse. This allowed us to analyze initially rotating data. We have combinedperturbation theory and full nonlinear hydrodynamics in spherical and axial sym-metric to construct a more complete picture of the critical phenomena emerging atthe threshold of gravitational collapse.First, we summarize what we’ve learned from studying our spherically symmetricmodels. The Euler equation of fluid dynamics with barotropic EoS (2.18) describe aisentropic fluid, which allows self-similar solutions of the first kind. Applying a self-similar ansatz we looked for analytic similarity solutions. We found the spectrum ofsolutions to be discrete and infinite, with the adiabatic index Γ in EoS (2.18) havingan effect on the spectrum of solutions. As a check of consistency we found thatthe previously known isothermal gas self-similar spectrum is recovered by settingΓ = 1. Increasing values of the adiabatic index was found to shift the amplitudefor the density similarity variable. This pattern persists for Γ ∈ [1, 6/5), where thespectrum of similarity solutions follow the same structure as that of the isothermalgas. This structure is identified as the Hunter-branch of solutions and the Larson-Penston solution. It was discovered that for Γ ≥ 6/5 the discrete but infiniteHunter-branch of solutions vanishes leaving the LP solution as the only similarity174Chapter 6. Conclusion and Further Worksolution in this regime. The LP solution persists for Γ ≤ 4/3.Calculation of the analytic, spherically symmetric linear perturbations aboutthe computed self-similar solutions reveal that the stability properties of the solu-tions mirror those of the isothermal gas. Again, by setting Γ = 1 we recoveredthe known results. More importantly, the first member of the Hunter branch, theso-called Hunter-A solution contains a single unstable mode. All members of theHunter family contain unstable modes. The Hunter-A is significant given that theformalism of Koike et. al. [99] states that the critical solution responsible for criti-cal phenomena must have a single relevant (unstable) mode. As we approached theunexpected transition value of Γ = 6/5, the Lyapunov exponent corresponding tothe unstable mode of the Hunter-A solution diverges. This seems to imply that theHunter-A solution becomes increasingly unstable as we approach this value for Γ.On the other hand, the LP solution is stable under linear perturbations. Therefore,we expect it to describe generic features of gravitational collapse.Under conditions of critical gravitational collapse we calculated (numerically)solutions to Euler’s equations of fluid dynamics subject to a polytropic EoS (1.4).We found that, similar to the isothermal gas the critical solution converges at inter-mediate times to the now Γ-dependent Hunter-A solution. At the late stages of theevolution the solution showed converges to the LP solution. This late time conver-gence to the LP solution was independent of the initial data. The calculation resultsfor Γ ∈ [1, 6/5) are consistent with type-II critical phenomena. However, the non-linear scaling between the time t and radial coordinate r required a modification tothe collapsed mass M and specific angular momentum scaling laws which acquiredan explicit dependence on the adiabatic index Γ. These results are consistent withthose presented in [48] (isothermal gas critical collapse).Interestingly, in the regime set by Γ ∈ [6/5, 4/3) no Hunter-A solutions existsfor the isentropic model. The critical behavior measured is consistent with type-Icritical phenomena. The critical solution evolves towards a metastable state whose175Chapter 6. Conclusion and Further Worklifetime T0 follows a liner scaling law. Following the departure from this state, pre-sumably due to an unstable mode, the solution evolves towards the correspondingLP solution. Therefore, we concluded that the polytropic gas model exhibits bothtypes of critical behavior depending on the parameter Γ. However, the type-I be-havior cannot be verified due to the inaccessibility the static solutions resulting fromnot having enough physical condition to solve the hydrostatic system, see Sec. 4.2.4.Our axisymmetric work constitutes the first study of critical fluid collapse wherethe effects of angular momentum were studied dynamically. In spherical symmetrythe fluid model plus Newtonian gravity is completely analogous the more realisticGR fluid model. Furthermore, it can be interpreted as an degenerate case of theGR system for k → 0 in EoS (1.3). Therefore, the results obtained hereby directlyapply to the GR system for vanishing k17. Alternatively, we can view the Newtoniangravity+fluid system as a toy model for studying the effect of slow rotation on thecritical solution.The analogy between the Newtonian and GR fluid models extends beyond spher-ical symmetry. The GR critical fluid collapse solution, the Evans-Coleman solutionis unstable against a nonspherical axial mode. Likewise, we confirmed the presenceof an axial unstable mode, the so-called spin-up mode on the Newtonian criticalcollapse solution (Hunter-A). In this case, the linear regime contains two growingmodes, the spherical and the spin-up modes. This linear regime is identical to thatdiscussed by Gundlach in [32]. Grounded on this analogy we borrowed Gundlach’slinear-regime-derived scaling behavior of the black-hole mass and its angular mo-mentum. These were used to generated our own predictions for the scaling of thecollapsed mass and its specific angular momentum. The results are in agreementwith the presence of an extra, axial growing mode.Our results indicate that all spherically symmetric self-similar solutions in ourmodel are unstable against the axial mode. Of particular interest is the Ori-Piran17This can only be correct for the isothermal gas, or for Γ→ 1 in the polytropic gas.1766.1. Extensions and Further Worksolution found in [77] which was shown to have a naked singularity in the regimewhere 0 < k < 0.0105. The solution being the GR equivalent of the LP solutioncontains no unstable spherical modes [43, 47]. Thus, it suggests that such a solutionshould be a global attractor of the evolution. However, according to [31, 32] all suchsolutions are unstable against nonspherical axial perturbations for k < 1/9. Ourresults are consistent with the instability of the Ori-Piran solution in the limit ask → 0 when rotation effects are included. Therefore, it cannot be a global attractorof the evolution.The detection of the growing axial mode about the critical collapse solutionsupports Gundlach’s [10, 33] suggested extension to the analogy with statistical me-chanics. A ferromagnet at high temperature T subject to an external magnetic field~Bext is analogous to the fluid under conditions of critical collapse with initial angu-lar angular momentum. The analogous parameters are p → T and ~q → ~Bext. Thevector quantity, the specific angular momentum ~a plays the role of the magnetiza-tion ~m. Both vanish according to a power-law at the critical point (p = p⋆, ~q = 0)for the gravity+fluid model and (T = T ⋆, ~Bext = 0) for the ferromagnet. The col-lapsed mass, or Newtonian black hole is analogous to the correlation length in theferromagnet system. Both systems manifest the properties of universality and scaleinvariance.6.1 Extensions and Further WorkAn obvious extension to the present project is the elimination of all symmetryrestrictions in the fluid model and carry out 3-dimensional critical collapse calcula-tions. This would allow us to observe in full generality all the growing nonsphericalmodes of the critical solution. It is unclear whether the Hunter-A solution is un-stable against a polar mode, the so-called bar mode. No evidence of it was foundin our axisymmetric simulations. Perturbation theory of the GR critical solution inthe Newtonian limit (k → 0) predicts that no other nonspherical modes (besides the1776.1. Extensions and Further Workaxial mode) should be present [33]. Three-dimensional dynamical calculation of crit-ical collapse as well as computation of the linear nonspherical polar perturbationsabout the Hunter-A solution would allow us to test this prediction.The polytropic gas model displayed behavior upon critical gravitational collapsethat is consistent with both, type-I and type-II critical phenomena. Large values ofΓ led to the onset of type-I critical phenomena. These solutions were first discoveredin this project. Their role as critical solutions suggests a one-mode linear stabilitystructure, however, an alternative method of computing them needs to be developedwhich would allow us to study their stability properties. In particular, calculationsof the Lyapunov exponent for these solutions are predicted to match their scalingexponent of the T0 scaling law. Equipped with this knowledge we can generate amore in depth analysis regarding the effects of initial rotation on these solutions. Itwould be useful to investigate these effects in three dimensions.Ultimately, we would like to construct a GR hydro code subject EoS (1.3) inaxial symmetry. This would allow us to vary the parameter k and test the otherperturbation-theory-based predictions made in [31–33]. This is a much more difficulttask due to the inherent complexities of solving simultaneously Einstein’s equationsand the conservation laws for mass, energy and momentum of the fluid in axisym-metry. Such difficulties are exacerbated when one considers extreme conditions ofcritical gravitational collapse. The Newtonian model presented in this project, al-though simpler than the GR system provides much of the same insights regardingthe nature of critical phenomena in gravitational collapse.178Bibliography[1] Andrew M. Abrahams and Charles R. Evans. Phys. Rev. Lett., 70:2980–2983,May 1993.[2] Andrew M. Abrahams and Charles R. Evans. Phys. Rev. D., 49:3998–4003,April 1994.[3] Nils Andersson. Relativistic fluid dynamics: Physics for many different scales.Living Rev. Relativity, 10(1), 2007.[4] G. I. Barenblatt. Scaling, self-similarity, and intermediate asymptotics, vol-ume 14 of Cambridge Texts in Applied Mathematics. Cambridge UniversityPress, 1996.[5] G. I. Barenblatt and B. Ya Zeldovich. Ann. Rev. Fluid Mech., 4:285–312,1972.[6] R. G. Barrera, G. A. Estevez, and J. Giraldo. Eur. J. Phys., 6:287–294, 1985.[7] Beverly K. Berger. Numerical approaches to spacetime singularities. LivingRev. Relativity, 5(1), 2002.[8] M. J. Berger and P. Colella. Journal of Computational Physics, 82:64–84,1989.[9] M. Birukou, V. Husain, G. Kunstatter, E. Vaz, and M. Olivier. Phys. Rev. D,65(104036):1–7, 2002.[10] Gundlach C. Living Rev. Relativity, 10(5), 2007.179Bibliography[11] A. H. Cahill and M. E. Tubb. Commun. Math. Phys., 21:1, 1971.[12] B. J. Carr and A. A. Coley. Gen. Relativ. Gravit., 16:31–71, 1999.[13] B. Carter. and R. N. Henriksen. Ann. Physique. Supp. No. 6, 14:47, 1989.[14] B. Carter. and R. N. Henriksen. J. Math. Phys., 32:2580–2597, 1991.[15] M. W. Choptuik. Phys. Rev. Lett., 70:9, 1993.[16] Matthew W. Choptuik. Mexican school on gravitation and mathematicalphysics. VII Mexican School: Lectures. 2006.[17] M.W. Choptuik, E.W. Hirschmann, S.L. Liebling, and F. Pretorius. Phys.Rev. D, 68(044007):1–9, 2003.[18] M.W. Choptuik, E.W. Hirschmann, S.L. Liebling, and F. Pretorius. Phys.Rev. Lett., 93(131101):1–4, 2004.[19] C. K. Chu and A. Sereny. J. Comput. Phys., 1974.[20] Gregory B. Cook. Initial data for numerical relativity. Living Rev. Relativity,3(5), 2000.[21] W. E. East, F. Pretorius, and B. C. Stephens. Phys. Rev. D, 85(124010), 2012.[22] C. R. Evans and J. S. Coleman. Phys. Rev. Lett., 72:1782, 1994.[23] J. A. Faber and F. A. Rasio. Binary neutron star merger. Living Rev. Rela-tivity, 15(8), 2012.[24] L. Ferna´ndez-Jambrina and L. M. Gonza´lez. Current Trends in RelativisticAstrophysics. Springer, Madrid, 2003.[25] Jose´ Font. Numerical hydrodynamics in general relativity. Living Rev. Rela-tivity, 6(4), 2003.180Bibliography[26] Jose´ Font. Numerical hydrodynamics and magnetohydrodynamics in generalrelativity. Living Rev. Relativity, 11(7), 2008.[27] J. B. Goodman and R. J. LeVeque. Math. Comp., 45(171):15–21, 1985.[28] E. Gourgoulhon. An introduction to relativistic hydrodynamics. EAS Publ.Ser., 21:43–79, 2006.[29] P. Grandcle´ment and J. Novak. Spectral methods for numerical relativity.Living Rev. Relativity, 12(1), 2009.[30] C. Gundlach and J. M. Martin-Garcia. Phys. Rev. D, 61:084024, 2000.[31] Carsten Gundlach. Phys. Rev. D, 57:7080–7083, 1998.[32] Carsten Gundlach. Phys. Rev. D, 65(064019):1–10, 2002.[33] Carsten Gundlach. Phys. Rev. D, 65(084021):1–22, Mar 2002.[34] Carsten Gundlach. Phys. Rev., 376:339–405, Mar 2003.[35] B. Gustaffsson, H-O Kreiss, and J. Oliger. Time Dependent Problems andDifference Methods. John Wiley & Sons, Inc., New York, 1995.[36] B. Gustafsson and H. O. Kreiss. J. Comput. Phys., 1979.[37] Benjamin Gutierrez. Relativistic Scattering of Solitons in Nonlinear FieldTheory. PhD thesis, The University of British Columbia, Vancouver, 2010.[38] T. Hanawa and T. Matsumoto. Astrophys. J., 521:703–707, 1999.[39] T. Hanawa and T. Matsumoto. Astrophys. J., 540:962–968, 2000.[40] T. Hanawa and T. Matsumoto. Publ. Astron. Soc. Japan, 52:241–247, 2000.[41] T. Hanawa and K. Nakayama. Astrophys. J., 484:238–244, 1997.[42] T. Hanawa, K. Saigo, and T. Matsumoto. Astrophys. J., 558:753–760, 2001.181Bibliography[43] T. Harada. Phys. Rev. D, 58(104015):1–10, 1998.[44] T. Harada. Class. Quantum Grav., 18:4549–4567, 2001.[45] T. Harada. Stability criterion for self-similar solutions with perfect fluids ingeneral relativity, 2001.[46] T. Harada. Self-Similar Solutions, Critical Behaviour and Convergence toAttractor in Gravitational Collapse, 2003.[47] T. Harada and H. Maeda. Phys. Rev. D, 63:084022, 2001.[48] T. Harada, H. Maeda, and B. Semelin. Phys. Rev. D, 67:084003, 2003.[49] A. Harten. J. Comput. Phys., 135:260–278, 1997.[50] A. Harten, P. Lax, and B. van Leer. SIAM review, 25(1):35–61, 1983.[51] J. B. Hartle. GRAVITY An Introduction to Einstein’s General Relativity.Addison Wesley, San Francisco, 2003.[52] S. H. Hawley and M. W. Choptuik. Phys. Rev. D, 62(104024):1–19, 2000.[53] James Healy and Pablo Laguna. Critical Collapse of Scalar Fields BeyondAxisymmetry. Gen. Relativ. Gravit., 46:1722, 2014.[54] G. W. Hedstrom. J. Comput. Phys., 1979.[55] S. Hod and T. Piran. Phys. Rev. D, 55(6):3485–3496, 1996.[56] C. Hunter. The Astronomical Journal, 218:434–845, 1977.[57] J. D. Anderson Jr. Modern Compressible Flow: with historical perspective.New York: McGraw-Hill, second edition, 1990.[58] R. Kippenhahn, A. Weigert, and A. Weiss. Stellar Structure and Evolution.Springer, New York, 2012.182Bibliography[59] T. Koike, T. Hara, and S. Adachi. Phys. Rev. D, 59:104008, 1999.[60] R. B. Larson. Mon. Not. R. astr. Soc., 145:271–295, 1969.[61] R. J. LeVeque. Finite-Volume Methods for Hyperbolic Problems. CambridgeUniversity Press, Cambridge, 2004.[62] A. R. Linddle. An introduction to cosmological inflation, 1999.[63] Y-Q. Lou and W-G Wang. Mon. Not. R. Astron. Soc, 372:885–900, 2006.[64] H. Maeda and T. Harada. Phys. Rev. D, 64:124024, 2001.[65] H. Maeda, T. Harada, H. Iguchi, and N. Okuyama. Phys. Rev. D, 66:027501,2002.[66] D. Maison. Phys. Lett. B, 366:82–84, 1996.[67] J. M. Mart´ı and E. Mu¨ler. Numerical hydrodynamics in special relativity.Living Rev. Relativity, 6(7), 2003.[68] T. Matsumoto and T. Hanawa. Astrophys. J., 521:659–670, 1999.[69] T. Matsumoto, T. Hanawa, and K. Nakamura. Astrophys. J., 478:569–584,1997.[70] C. W. Misner, K. S. Thorne, and J. A. Wheeler. Gravitation. W. H. Freeman,San Francisco, 1973.[71] Bruno Coutinho Mundim. A Numerical Study of Boson Star Binaries. PhDthesis, The University of British Columbia, Vancouver, 2010.[72] David Wayne Neilsen. Extremely Relativistic Fluids in Strong-Field Gravity.PhD thesis, The University of Texas, Austin, 1999.[73] D. H. Nielsen and M. W. Choptuik. Class. Quantum. Gravit., 17:761–782,2000.183Bibliography[74] D. H. Nielsen and M. W. Choptuik. Class. Quantum. Gravit., 17:733–759,2000.[75] S. C. Noble and M. W. Choptuik. Phys. Rev. D, 78:064059, 2008.[76] Scott Charles Noble. A Numerical Study of Relativistic Fluid Collapse. PhDthesis, The University of Texas, Austin, 2003.[77] A. Ori and T. Piran. Phys. Rev. D, 42:1068–1090, 1990.[78] I. Orlanski. J. Comput. Phys., 1976.[79] Anninos P. Living Rev. Relativity, 4(2), 2001.[80] T. Padmanabhan. Theoretical Astrophysics Volume II: Star and Stellar Sys-tems. Cambridge University Press, Cambridge, 2001.[81] M. V. Penston. Mon. Not. R. astr. Soc., 144:425–448, 1969.[82] Max Pettini. Physical cosmology. Lectures. 2010.[83] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. NumericalRecipes in FORTRAN. Cambridge University Press, New York, second edition,1993.[84] F. Pretorius and M. W. Choptuik. J. Comput. Phys., 218:246–274, 2006.[85] Frans Pretorius. Numerical Simulations of Gravitational Collapse. PhD thesis,University of British Columbia, Vancouver, 2002.[86] Dimitrios Psaltis. Probes and tests of strong-field gravity with observationsin the electromagnetic spectrum. Living Rev. Relativity, 11(9), 2008.[87] L. F. Richardson. The approximate arithmetical solution by finite differencesof physical problems involving differential equations, with and application tothe stresses in a masonry dam. Philosophical Transactions of the Royal Societyof London. Series A, 210:307–357, 1911.184Bibliography[88] P. L. Roe. J. Comp. Phys., 43:357–372, 1981.[89] D. H. Rudy and J. C. Strikwerda. J. Comput. Phys., 1980.[90] D. H. Rudy and J. C. Strikwerda. J. Comput. Phys., 1981.[91] F. W. Sears and G. L. Salinger. Thermodynamics, Kinetic Theroy, and Sta-tistical Thermodynamics. Addison-Wesley Publishing Company, third edition,1986.[92] M. Shibata and K. Taniguchi. Coalescence of black hole-neutron star binaries.Living Rev. Relativity, 14(6), 2011.[93] C-W. Shu. Leture Notes in Mathematics, 1697:325–432, 1998.[94] Chi-Wang Shu. Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory Schemes for Hyperbolic Conservation Laws. Technical report,NASA, November 1997.[95] F. H. Shu. The Astronomical Journal, 214:488–497, 1977.[96] M. Snajdr. Class. Quantum. Gravit., 23:3333–3352, 2006.[97] Evgeny Sorkin. On critical collapse of gravitational waves. Class.Quant.Grav.,28:025011, 2011.[98] Y. Suto and J. Silk. The Astronomical Journal, 326:527–538, 1988.[99] T. Hara T. Koike and S. Adachi. Phys. Rev. Lett., 74:5170, 1995.[100] Jonathan Thornburg. Event and apparent horizon finders for 3+1 numericalrelativity. Living Rev. Relativity, 10(3), 2007.[101] E. F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics.Springer, New York, 2009.185[102] E. F. Toro, M. Spruce, and W. Speares. Restoration of the contact surface inthe hll-riemann solver. Technical Report CoA-9204, Department of AerospaceScience, College of Aeronautics, Cranfield Institute of Technology, UK, March1993.[103] J. F. Ventrella and M. W. Choptuik. Phys. Rev. D, 68(044020):1–10, 2003.[104] R. M. Wald. General Relativity. The University of Chicago Press, Chicagoand London, 1984.[105] W-G Wang and Y-Q. Lou. Astrophys. Space Sci., 311:363–400, 2007.[106] S. Weinberg. Cosmology. Oxford University Press, Oxford New York, 2008.[107] Jeffrey Winicour. Characteristic evolution and matching. Living Rev. Rela-tivity, 15(2), 2012.[108] A. Yahil. The Astronomical Journal, 265:1047–1055, 1983.[109] J. M. Yeomans. Statistical Mechanics of Phase Transitions. Clarendon Press;Oxford University Press, New York; Oxford, 1992.[110] N. Yunes and X. Siemens. Living Rev. Relativity, 16(9), 2013.186Appendix ASolution’s Behavior Near SonicPoints for Γ = 1We assume the solution for the similarity variables α(x) and u(x) in Eqs. (2.52) and(2.53) is analytic. Therefore, at the sonic point the solution can be Taylor expanded.For the isothermal case, Γ = 1 the coefficients in the expansion can be computed inclosed-form in terms xs. There are two possible expansions at the sonic point, welabel these as type 1,α(x) =2xs+2xs − 6x2s(x− xs) + x2s − 8xs + 13x3s(x− xs)2 + . . . , (A.1)u(x) =(xs − 1) + xs − 1xs(x− xs)− xs − 12x2s(x− xs)2 + . . . (A.2)and type 2,α(x) =2xs− 2x2s(x− xs)− x2s − 6xs + 7x3s(2xs − 3)(x− xs)2 + . . . , (A.3)u(x) =(xs − 1) + 1xs(x− xs) + x2s − 5xs + 52x2s(2xs − 3)(x− xs)2 + . . . (A.4)The free parameter is xs. Both expansions can be used to set the starting point, (xs)for the “backward” integration to the matching point xM , discussed in Sec. BPolar PerturbationsGeneric nonspherical linear perturbations about any spherically symmetric solutionto the Eqs. (2.5)–(2.8) have the following general form,ρ(t, r, θ, φ)vr(t, r, θ, φ)vθ(t, r, θ, φ)vφ(t, r, θ, φ)P (t, r, θ, φ)ϕ(t, r, θ, φ)=ρ¯(t, r) + δρ(t, r, θ, φ)v¯r(t, r) + δvr(t, r, θ, φ)v¯θ(t, r) + δvθ(t, r, θ, φ)v¯φ(t, r) + δvφ(t, r, θ, φ)P¯ (t, r) + δP (t, r, θ, φ)ϕ¯(t, r) + δϕ(t, r, θ, φ), (B.1)when written in spherical coordinates (r, θ, φ). The “bar” quantities represent thesolution to the spherically symmetric model, Eqs. (2.24)–(2.27). Inserting theseansatz into Eqs. (2.5)–(2.8) we obtain the following system for the linear perturba-188Appendix B. Polar Perturbationstions,∂δρ∂t+1r2∂∂rr2(vrδρ+ ρδvr) +1r sin θ∂∂θ(sin θρδvθ) +1r sin θ∂∂φ(ρδvφ) = 0, (B.2)∂∂t(ρδvr + vrδρ) +1r2∂∂rr2(2ρvrδvr + v2rδρ) +ρvrr(1sin θ∂∂θ(sin θδvθ)+1sin θ∂δvφ∂φ)= −KΓρΓ−1[(Γ− 1)δρρ∂ρ∂r+∂δρ∂r]− ρ∂δϕ∂r− δρ∂ϕ∂r, (B.3)∂∂t(ρδvθ) +1r2∂∂r(r2ρvrδvθ) +ρvrδvθr= −1rKΓρΓ−1∂δρ∂θ− ρr∂δϕ∂θ, (B.4)∂∂t(ρδvφ) +1r2∂∂r(r2ρvrδvφ) +ρvrδvφr= − 1r sin θKΓρΓ−1∂δρ∂φ− ρr sin θ∂δϕ∂φ, (B.5)1r2∂∂rr2∂δϕ∂r+1r2 sin θ∂∂θsin θ∂δϕ∂θ+1r2 sin2 θ∂2δϕ∂φ2= 4πGδρ. (B.6)We have assume that the fluid obeys EoS (2.18), thus we can write the δP in termsof δρ. It clear that the fluid quantities represent the spherically symmetric solutions,therefore, we omit the bar-notation.Complementary to the discussion found in Sec. (2.5), we consider polar pertur-bations about the spherically symmetric self-similar solutionsρ(t, r, θ, φ) =14πGt2(α(x) + δα(x)eλτY mℓ (θ, φ)), (B.7)~v(r, θ, φ) =√κtn−1u(x) + δu(x)eλτY mℓ (θ, φ)δuΨ(x)eλτℓ+ 1∂∂θY mℓ (θ, φ)δuΨ(x)eλτ(ℓ+ 1) sin θ∂∂φY mℓ (θ, φ) , (B.8)ϕ(t, r, θ, φ) =κt2n−2(ϕ˜(x) + δϕ˜(x)eλτY mℓ (θ, φ)). (B.9)Where ϕ(t, r) ≡ κt2n−2ϕ˜(x) is the similarity variable for the Newtonian potentialdefined in terms of m(x) via Eq. (2.33) such that,∂ϕ˜(x)∂x=13n− 2m(x)x2. (B.10)189Appendix B. Polar PerturbationsInserting the polar perturbation ansatz Eqs. (B.7)–(B.9) into the equations for gen-eral linear perturbations (B.2)–(B.6) we obtain the following autonomous system,−(λ+ 2)δα − nx∂δα∂x+1x2∂∂xx2(αδu + uδα)− ℓxαδuΨ = 0, (B.11)(n− 3− λ)(αδu + uδα) − nx ∂∂x(αδu+ uδα) +1x2∂∂xx2(2αuδu + u2δα)− ℓxαuδuΨ = −Γ ∂∂x(αΓ−1δα) − α∂δϕ˜∂x− δα∂ϕ˜∂x, (B.12)(n− 3− λ)αδuΨ − nx ∂∂x(αδuΨ) +1x2∂∂xx2αuδuΨ +αuδuΨx=−ℓ+ 1x(ΓαΓ−1δα + αδϕ˜), (B.13)1x2∂∂xx2∂δϕ˜∂x− ℓ(ℓ+ 1)x2δϕ˜ = δα. (B.14)Where the parameters κ = K(4πG)1−Γ and n = 2− Γ.Again we only consider perturbations which are analytic. Applying the ex-pansions given by Eqs. (2.55) and (2.56) into the autonomous system for the polarperturbations, Eqs. (B.11)–(B.14) yield that the perturbations mush vanish at x = 0according to,δα(x) =δα0α2−Γ∗ xℓ, (B.15)δu(x) =δu0ℓxℓ−1, (B.16)δuΨ(x) =δu0(ℓ+ 1)xℓ−1, (B.17)δϕ˜(x) =−[Γδα0 −(λ+ 2Γ− 3 + ℓ(4− 3Γ)3)δu0]xℓ (B.18)Where δα0, and δu0 are free constant parameters.It can be shown that at the sonic point (nx − u) = √Γα(Γ−1)/2 the necessarycondition for regularity is,[(Γ− 1)α′α(nx− u)2 +(2ux− λ− 2)(nx− u) + (n− 1)u+ ϕ˜′]δα+[3(n − 1)− λ+ 2u′]αδu− ℓ(nx− u)xαδuΨ + αδϕ˜′ = 0 (B.19)190Appendix B. Polar PerturbationsWhere ∂(·)/∂x ≡ (·)′. The above condition, Eq. (B.19) can be used to solve for eitherδα or δu to be subsequently substituted into Eqs. (B.13) and (B.14) evaluated atthe sonic point. The solutions of these equations, i.e. δϕ˜′ and δuΨ are determinedup to and overall constant.191


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items