Open Collections

UBC Undergraduate Research

Critical Collapse of a Massive Scalar Field in Newtonian Gravity - Numerical Investigations of the Schodinger-Poisson.. Inwood, Andrew 2009

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
Inwood_Andrew_UBC_2009_Physics_449_Honours_Thesis.pdf [ 606.37kB ]
Metadata
JSON: 1.0085959.json
JSON-LD: 1.0085959+ld.json
RDF/XML (Pretty): 1.0085959.xml
RDF/JSON: 1.0085959+rdf.json
Turtle: 1.0085959+rdf-turtle.txt
N-Triples: 1.0085959+rdf-ntriples.txt
Original Record: 1.0085959 +original-record.json
Full Text
1.0085959.txt
Citation
1.0085959.ris

Full Text

Critical Collapse of a Massive ScalarField in Newtonian GravityNumerical Investigations of the Sch¨odinger-PoissonSystembyAndrew InwoodA THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFBACHELOR OF SCIENCE (HONOURS)inThe Faculty of Science(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)May 2009c© Andrew Inwood 2009AbstractCritical collapse in General Relativity (GR) has been studied extensively, both numerically andanalytically, for a wide variety of matter models including massive and massless scalar fieldscoupled to gravity, as well as a variety of gas models. Comparatively less work has been donetowards the investigation of critical collapse in Newtonian gravity, although some progress hasbeen made using an isothermal gas model. In this paper we numerically investigate criticalcollapse of a massive scalar field in spherical symmetry, coupled to gravity in the Newtonianlimit. The evolution of such a field is described by the Schr¨odinger-Poisson (SP) system.A Crank-Nicolson finite-differencing approximation is coded using the Rapid Numerical Pro-totyping Language (RNPL), and the updates are coded in FORTRAN. The SP system is investi-gated using Gaussian initial data, centred at the origin. By varying the amplitude of the Gaussiandata, no evidence for critical behaviour is found, and so the SP system is modified through thegravitational coupling in an attempt to induce critical behaviour. This modification is realizedby making the exponent in the Poisson equation variable.Type I critical behaviour is found using the Gaussian initial data centred at the origin for anexponent epsilon1 = 3. By performing a bisection search the critical value of the amplitude, A, is foundto 15 digits of precision. The critical solution is found to be periodic in time, characteristic ofType I phenomena, with a period of T ≈ 0.02. By surveying the 1-dimensional parameter spaceof A, the Lyapunov exponent associated with the growing mode in Type I phenomena is foundto 3 significant digits to be λ ≈ 163.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.1 Critical Collapse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21.2 Self-Similarity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.3 Cosmic Censorship . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.4 Newtonian Collapse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.5 Blow-up in Newtonian Self-Gravitating Systems . . . . . . . . . . . . . . . . . . . 62 Schr¨odinger-Poisson System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.2 Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2.1 Spherical Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.2.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.3 Finite Difference Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112.3.2 Tridiagonal Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.3.3 Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.3.4 Norm Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152.3.5 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.3.6 Change of Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18iiiTable of Contents3 Results: Schr¨odinger-Poisson System . . . . . . . . . . . . . . . . . . . . . . . . . 194 Modified Schr¨odinger-Poisson System . . . . . . . . . . . . . . . . . . . . . . . . . 254.1 Parameter Space Survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 264.2 Critical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 295 Summary and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 325.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34AppendicesA Perturbation Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35A.1 Type I Critical Phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35A.2 Type II Critical Phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36B Independent Residual Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . 38ivList of Figures1.1 The phase space portrait near a critical point. The red line represents a 1-parameter family of initial gravity configurations, and the blue lines extendingfrom the initial data represent trajectories in the phase space of possible Z. Thered line is drawn to cross the singularity threshold at p = pstar, and for this value ofp the gravity configuration tends towards the critical solution. Image taken fromGundlach[8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.1 The first bound-state wavefunction (S is a scaled wavefunction, Ψ). The solutionwas calculated by Moroz, Penrose, and Tod[13] numerically using a Runge-KuttaNAG routine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.2 The first excited state wavefunction. Each excited state has one more zero thanthe previous. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.3 Discretization of space-time domain. Grid functions (ie. V and Ψ) are labelledwith subscripts to indicate to the spatial coordinate and superscripts to indicatethe time coordinate. The spatial spacing is uniform and given by h, and the timespacing is set to be proportional to h. The proportionality constant λ is knownas the Courant factor. In a Crank-Nicolson scheme, the equations of motion arecentred at the half time level, which is coloured in red. . . . . . . . . . . . . . . . . 122.4 Time evolution of the norm for the Schr¨odinger equation, calculated at multipleresolutions. The difference between the value of the norm at the initial time ascalculated at levels 0 and 2 is about 1 part in 100,000. . . . . . . . . . . . . . . . . 183.1 For a small value of A, the self-gravitation is weak and the solutions immediatelydisperse. This means the probability density, |Ψ|2, becomes uniform as t →∞. . . 20vList of Figures3.2 The time derivative of |Ψ|2 at the origin is initially negative and tends to zero ast → ∞. This indicates that the solution is dispersing as t → ∞. The resolutionfor these calculations is 1025 points, unless indicated otherwise. . . . . . . . . . . . 213.3 When A = 800, the wavefunction initially grows at the origin before ultimatelydispersing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.4 For A = 800, the growth rate at the origin is initially positive, indicating that massis accumulating at the origin. Ultimately the growth rate becomes negative andtends to zero indicating dispersion. . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.5 By increasing A, it is observed that the growth rate at the origin exhibits similarbehaviour on shorter time scales. This is evidence that the SP system scales intime with A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.6 The growth rate for A = 5000 exhibits the same sort of behaviour as for A = 1500and A = 1000. For this calculation the resolution was increased to 2049 points. . . 244.1 For values of A below the critical value, the growth rate is initially positive butultimately becomes negative and tends to zero as the solution disperses, just as forthe SP system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 264.2 For A = 80, the growth rate makes a transition in behaviour, and appears to growwithout bound. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 274.3 By varying A through the critical value, the transition is clearly seen through thegrowth rate at the origin. For A = 80 the growth rate diverges from the other linesvery quickly. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284.4 The same picture as Fig.4.3, but on a larger time scale. Here it can be seen thatthe growth rate also increases without bound for A = 75.5, and that the time atwhich blow-up occurs is almost twice as late for A = 80. . . . . . . . . . . . . . . . 284.5 The solution exhibits periodic behaviour as is typical with Type I critical phenom-ena. The period of oscillation is found to be approximately, T ≈ 0.02. Resolutionis 32800 points on the range 0 < r < 40. . . . . . . . . . . . . . . . . . . . . . . . . 294.6 Lifetimes of subcritical solutions calculated at resolution of 32800 points. Thescaling exponent σ is determined from the slope, σ ≈ 0.0061. Supercritical datashould give similar results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31viList of FiguresB.1 The real part of the alternative residual calculated for the “modified SP system”ie. epsilon1 = 3, for subcritical data. The residual is calculated at the base level ofdiscretization, as well as the first and second level, with 1025 points at base level.Clearly as the level of discretization increases, the residual is tending to zero whichis the desired behaviour. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39viiAcknowledgementsI would like to thank my supervisor, Dr. Matthew Choptuik, for his support in this project. Myengaging conversations with Dr. Choptuik have caused me to develop a strong interest in thefields of Gravitation and Computational Physics, and I am certain that the knowledge I havegained from him will serve me well in my future endeavours.1Chapter 1IntroductionIt has long been known that solutions to Einstein’s equations with physically reasonable matterin the initial state can develop singularities [10]. In fact, if one can construct a so called “trappingregion” then the Hawking-Penrose singularity theorems guarantee that a singularity will form.The fact that singularity formation in General Relativity (GR) has been proven is interestingbecause it is a sign that the validity of the theory has broken down in this region. In order tofully understand the physics in this region, a theory of quantum gravity will be required.The saving grace of GR is that in almost all situations, singularities are enclosed by an eventhorizon—a so called “clothed singularity”. This means that the singularity would not be causallyconnected to an observer lying outside of the event horizon, and so the predictions of GR areunaffected in most practical applications. However, beyond the singularity theorems, the detailsof singularity formation are not well understood, and it is not clear that generically singularitieswill be clothed, and might in fact form beyond the confines of an event horizon. Such a singularityis called a “naked singularity”. This is related to the topic of Cosmic Censorship, which will bediscussed further.1.1 Critical CollapseBy 1973 it was clear that singularities must form in GR from physically reasonable and sufficientlymassive initial states, and it was also clear that for sufficiently light initial states, singularityformation would be impossible. However, due to the complexity of the Einstein equations in GR,it is difficult to say whether or not a singularity will form in the “middle ground” where the massof the initial state is neither light nor heavy. This, as well as the question of whether or notsingularities are all clothed, motivated Choptuik[3] to investigate the middle ground numericallyin the early 1990’s.Choptuik investigated a spherically symmetric, massless scalar field coupled to gravity in GR,21.2. Self-Similarityand by varying 1-parameter families of initial data, investigated the transition between solutionsin which all of the mass disperses away from the origin, and those in which the mass accumulatesat the origin and forms a singularity. For a phase space picture of this, see figure 1.1. Remarkably,it was discovered that for initial data precisely tuned between the dispersive and singular types,the solution would evolve toward a critical solution which was independent of the family of initialdata—this property is known as the Universality of the critical solution. This behaviour ofthe solutions is analogous to phase transitions in statistical mechanics. There are two types oftransitions, named Type I and Type II, which are discussed in the appendices. Furthermore,the critical solutions were found to exhibit a discrete self-similarity, and also contained a nakedsingularity at the origin.1.2 Self-SimilarityFurther numerical investigations on the critical behaviour of self-gravitating systems[6][5] revealedtwo different kinds of critical solutions: ones that are time-independent, and ones that exhibiteither continuous self-similarity (CSS) or discrete self-similarity (DSS) (for a review of self-similarsolutions in GR, see [2]), and are labelled Type I and Type II critical phenomena, respectively.They are so labelled because of their resemblance to Type I and Type II phase transition instatistical mechanics. For Type I critical phenomena the mass of the end state singular solutionis a certain fraction of the mass of the critical solution, but for Type II the mass of the end statesingular solution obeys a power scaling law, and initial conditions can be set to make the resultingblack hole mass infinitesimally small.1.3 Cosmic CensorshipCosmic Censorship encompasses the idea that even though singularities exist in GR, they may nothave any influence on the physics outside of the event horizon which surrounds it. At present, anystatement of cosmic censorship is at the level of a conjecture, because the available mathematicaltools are insufficient to provide anything as far as a mathematical proof or theorem. There aretwo basic forms of the cosmic censorship conjecture, the weak and strong versions.The weak cosmic censorship conjecture[14] claims that any singularities which form will beenclosed by an event horizon, and therefore the physics beyond the confines of the event horizon31.3. Cosmic CensorshipFigure 1.1: The phase space portrait near a critical point. The red line represents a 1-parameterfamily of initial gravity configurations, and the blue lines extending from the initial data representtrajectories in the phase space of possible Z. The red line is drawn to cross the singularitythreshold at p = pstar, and for this value of p the gravity configuration tends towards the criticalsolution. Image taken from Gundlach[8].41.4. Newtonian Collapsewill be unaffected. This is the “weak” version because it still allows for the singularity to affectobservers who travel through the event horizon, into the black hole. The strong cosmic censorshipconjecture[11] claims that time-like singularities do not form, and therefore even observers whotravel into the black hole will not be affected by the singularity.The “middle ground” self-gravitating regime is an excellent laboratory to investigate cosmiccensorship, since naked singularities have been found in numerical studies of spherical criticalcollapse. However, since most evidence for naked singularity formation is empirical, it is notalways clear whether or not the singularities are stable under perturbations to the model, suchas evolution beyond spherical symmetry. The precise form of cosmic censorship is one of thebiggest open questions in classical GR. Determining the precise form of cosmic censorship andunderstanding under what conditions self-similarity arises are two of the main goals of researchin the field of gravitational collapse.1.4 Newtonian CollapseInvestigations of critical collapse in purely Newtonian self-gravitating systems have been car-ried out for a Newtonian isothermal gas[9], and also scalar fields with or without mass or self-interaction potential. These systems also exhibit self-similar solutions which have naked singular-ities, and so provide another window for investigating the range of validity of cosmic censorship.In the field of critical collapse, much less work has been done in Newtonian gravity comparedto GR; however, due to the reduced complexity associated with Newtonian gravity, it offers agood laboratory for investigating critical collapse under deviations from spherical symmetry andother simplifications that are usually commonplace. For these reasons, working with Newtoniangravity could offer insight into cosmic censorship and self-similarity. It should be noted that sincethese analyses are purely Newtonian, they could potentially have been carried out long beforethe development of GR. There is presently no explanation for the fact that solutions exhibit self-similarity in terms of dynamics, so an understanding of the Newtonian case may provide insightinto the role that the dynamics play in creating self-similar solutions.51.5. Blow-up in Newtonian Self-Gravitating Systems1.5 Blow-up in Newtonian Self-Gravitating SystemsHere we are interested in investigating critical collapse in Newtonian self-gravitating systemswhich exhibit “blow-up”, ie. solutions which become singular in finite time. Such solutions areinteresting for a number of reasons. As mentioned above, investigating the evolution of solutionsas they become singular will provide insight into how singular solutions form out of the dynamics.Blow-up solutions are also interesting in the context of Newtonian gravity because it is not clearthat critical phenomena will occur for any matter model. This is because in GR, the universalityof coupling of mass, energy, and momentum guarantees that there will be black hole formationin any matter model; however, Newtonian gravity only sets a universality of coupling of massto momentum, and it is not clear that all matter models will admit initial data that result inblow-up.6Chapter 2Schr¨odinger-Poisson System2.1 BackgroundIn GR, the evolution of a self-gravitating, complex massive scalar field satisfies the Einstein-Klein-Gordon (EKG) equation. In the Newtonian limit of weak gravity and slow velocities comparedto the speed of light, the EKG equation can be shown to reduce to the Schr¨odinger-Poisson (SP)system[7]:i∂Ψ∂t = −12∇2Ψ+gVΨ (2.1)∇2V = |Ψ|2 (2.2)where we are working in units with G = c = ¯h = 1, g is a coupling constant, and the physicalquantities are given byxphy = ¯hmcxcomp (2.3)tphy = ¯hmc2tcomp (2.4)Ψphy = c2√4piGΨcomp (2.5)The physical interpretation of the SP system is that the wavefunction of a single particleevolves according to the Schr¨odinger equation, and feels a potential that is produced by the par-ticle itself through |Ψ|2. Therefore the self-gravitation is produced through the Poisson equation.The SP system is known to have an infinite set of spherically symmetric static solutionscalled Newtonian Boson Stars (whereas the static solutions to the EKG equations are calledRelativistic Boson Stars), see figure 2.1. These solutions are smooth at the origin, and can beordered according to the number of zeroes they have.It is not known as far as a mathematical proof whether or not this system permits solutionsthat exhibit blow-up, and so the first order of business in this research is to solve the SP system72.1. BackgroundFigure 2.1: The first bound-state wavefunction (S is a scaled wavefunction, Ψ). The solution wascalculated by Moroz, Penrose, and Tod[13] numerically using a Runge-Kutta NAG routine82.1. BackgroundFigure 2.2: The first excited state wavefunction. Each excited state has one more zero than theprevious.92.2. Constraintsnumerically using finite-differencing techniques, and vary a 1-parameter family of initial data inan attempt to find blow-up solutions.2.2 Constraints2.2.1 Spherical SymmetryWorking in spherical symmetry allows us to write the laplacian as:∆f = 1r2 ∂∂r r2∂f∂r¶(2.6)However, smoothness at the origin will be better enforced by the numerical scheme if the laplacianis rewritten in terms of ∂(r3),∆f = 3 ∂∂(r3) r2∂f∂r¶(2.7)2.2.2 Boundary ConditionsAt the inner boundary (near the origin) we imposed that V and Ψ are smooth:dVdr |r=0 = 0 (2.8)dΨdr |r=0 = 0 (2.9)As for the outer boundary, the natural boundary condition is that both V and Ψ go to zero asr goes to infinity. Although it is possible to write an equivalent boundary condition at a finiter [1], these are difficult to implement numerically. There are a number of options available fordealing with the outer boundary numerically, such as a so called “absorbing boundary condition”[12]; however, these sort of advanced techniques should not be required to establish whether ornot blow-up occurs in the SP system, so to simplify things we impose that all functions are zeroat a finite outer boundary r = 1,V|r=1 = 0 (2.10)Ψ|r=1 = 0 (2.11)The main disadvantage to this approach is that it limits how long a time over which theequations can be integrated, because if the wavefunction propagates out to the outer boundary itwill reflect and the solution becomes unphysical. But again, if blow-up is occurring in the model,102.3. Finite Difference Approximationthen all the mass falls into the origin very quickly, faster than the wavefunction can propagate tothe outer boundary.2.3 Finite Difference Approximation2.3.1 DiscretizationTo discretize the spatial domain into Nr grid points, FORTRAN notation is used, ie:r1 = 0r2 = hr3 = 2h...rNr = 1Figure 2.3 displays a schematic diagram of the discretized domain. h is the standard notation forthe spacing of the spatial coordinate, and we use h and dr synonymously. Note that subscriptsindicate the spatial point, and superscripts indicate the time level.In order to discretize the laplacian as written in equation 2.7, and maintain accuracy to O(h2)it is necessary to introduce “half spatial points”:rj+12≡ rj + dr2 (2.12)rj−12≡ rj − dr2this slight abuse of notation makes the equations more readable. Now the laplacian takes theform:∆f = 3r3j+12 −r3j−12 r2j+12∗ fj+1 −fjdr¶−r2j−12∗ fj −fj−1dr¶¶(2.13)We choose to implement a Crank-Nicolson FDA. The essence of a Crank-Nicolson scheme isthat the equations of motion are centred at the so called half time level. This is simple to do forthe time-derivative of Ψ, by writing:Ψn+1j −Ψnjdt = ∂Ψ∂t¶n+12j+O(h2) (2.14)112.3. Finite Difference ApproximationFigure 2.3: Discretization of space-time domain. Grid functions (ie. V and Ψ) are labelled withsubscripts to indicate to the spatial coordinate and superscripts to indicate the time coordinate.The spatial spacing is uniform and given by h, and the time spacing is set to be proportional toh. The proportionality constant λ is known as the Courant factor. In a Crank-Nicolson scheme,the equations of motion are centred at the half time level, which is coloured in red.122.3. Finite Difference ApproximationTo centre the right-hand side of the Schr¨odinger equation at the half time level and maintainO(h2) accuracy, it is necessary to average the Ψn and Ψn+1 time levels. Taking into account allof these considerations, and using the discretized laplacian, the Schr¨odinger equation becomes:i‡Ψn+1j −Ψnj·dt = −12 ·3r3j+12−r3j−12"r2j+12· 12‡Ψn+1j+1 −Ψn+1j +Ψnj+1 −Ψnj·dr (2.15)−r2j−12· 12‡Ψn+1j −Ψn+1j−1 +Ψnj −Ψnj−1·dr#+ g2‡V n+1j Ψn+1j +V nj Ψnj·and the FDA of the Poisson equation becomes|Ψnj|2 = 3r3j+12 −r3j−12"r2j+12· V nj+1 −V njdr¶−r2j−12· V nj −V nj−1dr¶#(2.16)Next the boundary conditions need to be treated numerically. The outer boundary is trivial,and can be satisfied by setting the outermost values of V and Ψ equal to zero. The inner boundarycan be treated by introducing the O(h2) accurate forward-differencing derivative operator:(fr)nj = −32fnj +2fnj+1 − 12fnj+2h +O(h2) (2.17)Using this operator, the inner boundary condition can be rewritten as:fn1 = −13fn3 + 43fn2 (2.18)where f is either V or Ψ. This essentially allows for all references to V1 and Ψ1 to be removedfrom the finite-differencing scheme, since we can now solve a smaller set of equations to getΨ2,...,ΨNr, and use eq. 2.18 to get Ψ1 trivially.2.3.2 Tridiagonal FormAs it turns out, both eqs. 2.15 and 2.16 can be brought into a canonical form called tridiagonalform. The Poisson equation in tridiagonal form looks like:c−V n+1j−1 +c0V n+1j +c+V n+1j+1 = ¡ΨjΨ∗j¢n+1 (2.19)132.3. Finite Difference Approximationwherec− = 3dr0@r2j−12r3j+12−r3j−121A (2.20)c0 = − 3dr0@r2j−12+r2j+12r3j+12−r3j−121A (2.21)c+ = 3dr0@r2j+12r3j+12−r3j−121A (2.22)At the inner boundary, eq. 2.19 will have the form:c−V1 +c0V2 +c+V3 = Ψ(2)Ψ∗(2) (2.23)However, the reference to V1 can be removed using the boundary condition, eq. 2.18:V1 = 43V2 − 13V3 (2.24)Then the equations at the boundary become: 43c− +c0¶V2 + c+ − 13c−¶V3 = Ψ(2)Ψ∗(2) (2.25)NowallreferencestoV1 havebeenremoved, andwearesolvingasmallersystemforV2,...,VNr.Note that having removed the reference to V1, the Poisson equation, eq. 2.19, can be written inmatrix form as:0BBBBBBBBB@43c− +c0 c+ −13c− 0 0 0 0c− c0 c+ 0 0 0 ...0 c− c0 c+ 0 00 0 c− c0 c+ 0...1CCCCCCCCCA0BBBBBBBBBBBBBBBB@V2V3V4V5V6V7...1CCCCCCCCCCCCCCCCA=0BBBBBBBBBBBBBBBB@|Ψ2|2|Ψ3|2|Ψ4|2|Ψ5|2|Ψ6|2|Ψ7|2...1CCCCCCCCCCCCCCCCA(2.26)Where it is understood that c−, c0, and c+ are all dependent on r. Only the three main diagonalsof the matrix operating on V are non-zero, and this is why this form is called tridiagonal form.The Schr¨odinger equation can also be brought into tridiagonal form, although this equationwill depend on values of V n+1 and also the nth time level V n and Ψn. We just wrap the extra142.3. Finite Difference Approximationdependence in a function F, and then the tridiagonal form is:d−Ψn+1j−1 +d0Ψn+1j +d+Ψn+1j+1 = F(j) (2.27)whered− = 12 3r3j+12 −r3j−12r2j−12· 12 dtdr (2.28)d0 = i− 12 3r3j+12 −r3j−12"r2j+12· 12 dtdr +r2j−12· 12 dtdr#− g2V n+1j (2.29)d+ = 12 3r3j+12 −r3j−12r2j+12· 12 dtdr (2.30)F(j) = iΨnj − 12 3r3j+12 −r3j−12"r2j+12· 12 dtdr ¡Ψnj+1 −Ψnj¢−r2j−12· 12 dtdr ¡Ψnj −Ψnj−1¢#+ g2V nj Ψnj(2.31)And again the equation near the inner boundary can be written as: 43d− +d0¶Ψ2 + d+ − 13d−¶Ψ3 = F(2) (2.32)2.3.3 CodeWe have now developed all of the machinery to take an initial state Ψ1 and calculate the stateat the next time level Ψ2 satisfying the SP system and boundary conditions. This calculationhas been reduced to inverting tridiagonal matrices. However, since the Schr¨odinger and Poissonequations are coupled, an initial guess at V n+1 is required to invert the matrix for Ψ. Then V n+1and Ψn+1 are updated iteratively until both are within tolerance.Routines have been developed to quickly and efficiently invert tridiagonal matrices, and weuse the routines found in LAPACK1 for tridiagonal matrices: DGTSV and ZGTSV. A skeletoncode for solving the SP system is written in RNPL (Rapid Numerical Prototyping Language),and the updates are written in FORTRAN 77.2.3.4 Norm CalculationAs usual, the physical interpretation of |Ψ|2 is that it is the spatial probability distributionfunction of the particle. This means that integrating over it tells us the probability of finding the1LAPACK (Linear Algebra PACKage): http://www.netlib.org/lapack/152.3. Finite Difference Approximationparticle in a region of space. However, since we are working in spherical symmetry, there is anr2 weighting in this probability, so even though |Ψ|2 may be large at the origin, the probabilityof finding the particle there may be small because of this weighting.Therefore, to get a good impression of where the mass is, it is important to look at theintegrated probability:I(r) = 4piZ r0ξ2|Ψ(ξ)|2 dξ (2.33)If the equations of motion have been set up such that Ψ is exhibiting blow-up behaviour, then onecould tell that most of the mass is concentrated near the origin if eq. 2.33 goes to 1 for a smallvalue of r. To calculate this integral numerically we employ a Riemann sum using the midpointrule,I(rj) = I(rj−1)+4pidr r(j)− dr2¶2 (Ψj +Ψj−1)2‡Ψ∗j +Ψ∗j−1·2 (2.34)and I(r1) ≡ 0.2.3.5 Error AnalysisThe ansatz underlying all of our error analysis, originally observed by Richardson, is that thesolution, Ψ, can be expanded in a power series in h, the mesh spacing:Ψ(t,r;h) = Ψ(t,r)+he1(t,r)+h2e2(t,r)+... (2.35)The method of convergence testing we employ requires us to calculate the solution of the FDAusing three different values of h, in the ratio 1 : 2 : 4. Also noting that for centred FDA, like theCrank-Nicolson scheme, the odd powers of h in the Richardson expansion (equation 2.35) dropout, then the expansions at different h are:Ψ(t,r;h) = Ψ(t,r)+h2e2(t,r)+h4e4(t,r)+... (2.36)Ψ(t,r;h) = Ψ(t,r)+(2h)2e2(t,r)+(2h)4e4(t,r)+... (2.37)Ψ(t,r;h) = Ψ(t,r)+(4h)2e2(t,r)+(4h)4e4(t,r)+... (2.38)We can now compute the convergence factor, Q(t):Q(t) ≡ ||Ψ4h −Ψ2h||||Ψ2h −Ψh|| (2.39)162.3. Finite Difference ApproximationWhere ||·|| is any suitable spatial norm. For the case of the SP system in spherical symmetry,Q(t) takes the form:Q(t) =Rrmax0 r2(Ψ4h −Ψ2h)drRrmax0 r2(Ψ2h −Ψh)dr(2.40)It can be easily shown that if the FDA is converging then:limh→0Q(t) = 4 (2.41)For a more detailed discussion of convergence testing based on the Richardson expansion forFDAs, see the notes by Choptuik[4].To check that the PDE had been coded properly, it is useful to use the norm of Ψ as adiagnostic. For the Schr¨odinger equation, not coupled to gravity, the norm is explicitly conserved,which results from the fact that the total probability of finding a particle somewhere in spaceis always unity. Checking for norm conservation is a way to check that the code is producingsensible results.Figure 2.4 displays the time evolution of the norm of Ψ for a free space Schr¨odinger equation,calculated at multiple resolutions. The plots are converging to a constant, which is the expectedbehaviour for the Schr¨odinger equation. The norm calculated at different resolutions do not agreeat the initial time due to the fact that the finite mesh spacing introduces uncertainty into thenorm calculation; however, at the initial time the norm as calculated at level 0 and level 2 varyby only about 1 part in 100,000.The norm is so closely conserved for the Schr¨odinger equation because the Crank-NicolsonFDA itself also conserves the norm, and so the norm is conserved to greater precision than mightbe expected by an analysis involving only the Richardson expansion. This poses a problem inusing the PDE as a diagnostic, since in the case that the FDA has not been coded properly,we may be convinced of its correctness based on norm conservation attributable to the Crank-Nicolson scheme. To test that the FDA has been coded properly, the results produced by ithave to be compared to those produced by an alternative FDA (ideally one which is not normpreserving). Such a comparison is called an alternative residual calculation, the implementationof which is outlined in appendix B.172.3. Finite Difference ApproximationFigure 2.4: Time evolution of the norm for the Schr¨odinger equation, calculated at multipleresolutions. The difference between the value of the norm at the initial time as calculated atlevels 0 and 2 is about 1 part in 100,000.2.3.6 Change of VariablesIt should be noted that if self-similarity is found in the process of these investigations, it willbecome necessary to solve the equations under arbitrary change of variables, depending on theform of the self-similarity. This has mostly been implemented into the code, although at the timeof writing this feature has not been required, as self-similarity has not been found.18Chapter 3Results: Schr¨odinger-Poisson SystemInvestigations were carried out for a 1-parameter family of initial data in hopes of finding criticalbehaviour in the SP system. We choose to initialize the wavefunction, Ψ, as being real, Gaussian,and centred at the origin. Although any parameter would do (either the coupling constant g,the amplitude of the Gaussian data A, or even the width of the Gaussian), in this research weare interested in investigating which matter models permit blow-up, and so we choose to fix thecoupling constant g equal to 1, and vary the amplitude of the Gaussian, since this correspondsto increasing the total mass of the initial state.For A = 1, figure 3.1, the mass of the initial state is light enough that the self-gravity of thesystem is not significant, and the solutions behave like solutions to the free space Schr¨odingerequation. The mass disperses from the origin, ultimately evolving towards a state of equalprobability on the interval 0 ≤ r ≤ 1.Figure 3.2 displays the time derivative of |Ψ|2 evaluated at the origin. For A = 1, this rate isobserved to be negative at the initial time, and goes to zero in the limit as t → ∞. If blow-upwere to occur, this rate should be observed to grow without bound, which is why it is useful tolook at these plots.If A is set much larger than 1, a different behaviour is observed. For A = 800, figure 3.3, theself-gravity of the system is now strong enough that initially the mass starts to accumulate at theorigin, but as t → ∞ the mass still disperses. Again referring to the growth rate at the origin,figure 3.4, now the rate is positive at the initial time, but ultimately becomes negative as themass starts to disperse away from the origin. Based on the observed change of behaviour as Ais increased, it might be expected that for sufficiently large A the self-gravitation will be strongenough to cause blow-up in the solutions.Figures 3.5 and 3.6 display the growth rate at the origin for increasing values of A. What isseen is that the overall character of the growth rate remains the same with some scaling in time.This evidence suggests that there is no blow-up in the SP system as it stands, since the growth19Chapter 3. Results: Schr¨odinger-Poisson SystemFigure 3.1: For a small value of A, the self-gravitation is weak and the solutions immediatelydisperse. This means the probability density, |Ψ|2, becomes uniform as t →∞.rate seems to scale with A.20Chapter 3. Results: Schr¨odinger-Poisson SystemFigure 3.2: The time derivative of |Ψ|2 at the origin is initially negative and tends to zero ast → ∞. This indicates that the solution is dispersing as t → ∞. The resolution for thesecalculations is 1025 points, unless indicated otherwise.21Chapter 3. Results: Schr¨odinger-Poisson SystemFigure 3.3: When A = 800, the wavefunction initially grows at the origin before ultimatelydispersing.22Chapter 3. Results: Schr¨odinger-Poisson SystemFigure 3.4: For A = 800, the growth rate at the origin is initially positive, indicating that massis accumulating at the origin. Ultimately the growth rate becomes negative and tends to zeroindicating dispersion.Figure 3.5: By increasing A, it is observed that the growth rate at the origin exhibits similarbehaviour on shorter time scales. This is evidence that the SP system scales in time with A.23Chapter 3. Results: Schr¨odinger-Poisson SystemFigure 3.6: The growth rate for A = 5000 exhibits the same sort of behaviour as for A = 1500and A = 1000. For this calculation the resolution was increased to 2049 points.24Chapter 4Modified Schr¨odinger-Poisson SystemWhile it has not been proven that blow-up cannot occur in the SP system, we should keep inmind that the main goal of this research is to investigate blow-up in Newtonian self-gravitatingsystems. Since there is currently no evidence that the SP system will permit blow-up, in order toproceed with the research we modify the SP system in a way that is expected to permit blow-upsolutions. It should be noted that the modification is only motivated by our interest in producingblow-up solutions, and not motivated by an attempt to model a known physical system. In thisresearch we hope to find what kind of matter models will permit blow-up, and whether or notthese models are physically relevant is to be decided by an as of yet undiscovered theory.In an attempt to induce blow-up, one could introduce non-linearities into the SP system.There are many ways one might consider doing this, but not all are relevant to gravitationalcollapse. For example, one could try introducing a Ψ|Ψ|2 term into the Schr¨odinger equation toproduce the usual non-linear Schr¨odinger (NLS) equation, but now coupled to a Poisson equation;however, this is guaranteed to produce blow-up solutions for the focussing NLS equation, and doesnot result from the gravitational coupling.In this research, it is desired to investigate blow-up which occurs as a result of the gravitationalcoupling, and so we choose to make the following modification to the SP system:i∂Ψ∂t = −12∆Ψ+gVΨ (4.1)∆V = |Ψ|epsilon1where epsilon1 ≥ 2. The only difference here is the epsilon1 in the Poisson equation. All of the machinerydeveloped in the earlier chapters can be recycled, including the FDA (with trivial modification ofthe discretized Poisson equation) and independent residual checker, and so will not be reproducedhere.254.1. Parameter Space SurveyFigure 4.1: For values of A below the critical value, the growth rate is initially positive butultimately becomes negative and tends to zero as the solution disperses, just as for the SPsystem.4.1 Parameter Space SurveyHere we present results for when epsilon1 = 3. Again we set g = 1 and investigate the same 1-parameterfamily of initial data, the Gaussian centred at the origin.For A = 70, the solution has a similar behaviour to solutions to the regular SP system, wherethe growth rate at the origin is initially positive as the mass accumulates at the origin, butultimately becomes negative and tends to zero as the mass disperses, see figure 4.1. However, forA = 80, figure 4.2, a very different behaviour arises, where the growth rate remains positive forall time and seems to increase without bound.Figures 4.3 and 4.4 show the growth rate at the origin for the values A = {70,75.3,75.5,80}.The change in behaviour of the solutions as A is increased can be readily seen in these graphs.For sufficiently small A the growth rate tends to zero whereas for large A the growth rate growswithout bound. In fact, there is a change in behaviour between A = 75.3 and A = 75.5. ForA = 75.5, the growth rate tends to infinity, although this is not apparent in figure 4.3. In figure4.4 it is clear that for A = 75.5, the growth rate eventually increases rapidly, although this is at264.1. Parameter Space SurveyFigure 4.2: For A = 80, the growth rate makes a transition in behaviour, and appears to growwithout bound.a much later time than for A = 80.This is clear evidence of critical behaviour, for presumably if A is tuned arbitrarily close tothe critical value, then the growth rate will start to increase at arbitrarily large times. Thisagrees with the standard picture of Type I critical collapse as described in the appendices, wherethe solutions tend towards a so called critical solution before “deciding” between a singular ordispersive end state.274.1. Parameter Space SurveyFigure 4.3: By varying A through the critical value, the transition is clearly seen through thegrowth rate at the origin. For A = 80 the growth rate diverges from the other lines very quickly.Figure 4.4: The same picture as Fig.4.3, but on a larger time scale. Here it can be seen thatthe growth rate also increases without bound for A = 75.5, and that the time at which blow-upoccurs is almost twice as late for A = 80.284.2. Critical Solution4.2 Critical SolutionTo get the tightest bound on the critical value of A, a bisection search was performed to acquireAstar to 15 digits of precision. The calculations were performed at higher resolutions to check theconvergence of Astar, and at the time of this writing, based on the convergence as we increase theresolution, Astar has been determined to 3 significant digits. The closer we tune A to the criticalvalue, the longer the solution will persist before dispersing or becoming singular. This poses aproblem because reflections off the outer boundary limits the time over which the solutions canevolve reliably, and so to deal with this we increase the spatial domain to 0 < r < 40.Figure 4.5: The solution exhibits periodic behaviour as is typical with Type I critical phenomena.The period of oscillation is found to be approximately, T ≈ 0.02. Resolution is 32800 points onthe range 0 < r < 40.Figure 4.5 shows the time evolution of the point at r = 0, using a resolution of 32800 points onthe interval 0 < r < 40. The solution is oscillatory, with a period of T ≈ 0.02. The periodic natureof the solution is also typical of Type I critical behaviour. As is described in the appendices, the294.2. Critical Solutionnatural parameterization, pi, of Type I phenomena is given by:pi = ln|A−Astar| (4.2)This form is natural because it is found that the “lifetime”, τ, of the solution, the time up towhich the solution is “close” to the critical solution, is proportional to pi in the limit that themesh spacing goes to zero.τ → σpi (4.3)To find σ we calculate the slope of a τ vs pi plot, and check the convergence of σ with increasingresolution. In order to do this we need a working definition of τ. We define τ as the time at whichthere is a 1% fractional difference between the solution and the critical solution at the origin, ie:||Ψ|2(A)−|Ψ|2(Astar)||Ψ|2(Astar) |r=0 > 0.01 (4.4)Of course we only know Astar to a certain precision set by our bisection search, and so when lookingat subcritical data we use the lower bound on Astar produced by the bisection search, and whenlooking at supercritical data we use the upper bound. The fact that we do not know Astar to infiniteprecision means that pi has more uncertainty when pi →−∞, and so when calculating σ we haveto be careful to look at values of pi such that the effects of approximating Astar are small.Figure 4.6 shows a plot of τ vs pi, using subcritical data. By calculating the slope using dataproduced at different resolutions, we can see that σ has been found accurately to 2 significantdigits, with σ ≈ 0.0061. Referring to the perturbation analysis in the appendices, σ is related tothe Lyapunov exponent, λ, associated with the growing mode in Type I solutions, as λ = −1/σ.From this we see that λ ≈ 163.304.2. Critical SolutionFigure 4.6: Lifetimes of subcritical solutions calculated at resolution of 32800 points. The scalingexponent σ is determined from the slope, σ ≈ 0.0061. Supercritical data should give similarresults.31Chapter 5Summary and Future Work5.1 SummaryPartially motivated by work on a Newtonian isothermal gas model [9], the Schr¨odinger-Poisson(SP) system was investigated to establish whether or not it permits blow-up, ie. solutions which,starting from physically reasonable initial data, become singular in finite time. The SP systemwas investigated numerically using a Crank-Nicolson finite-difference approximation (FDA), anda family of 1-parameter initial data was surveyed to look for blow-up. As the parameter wasincreased, the general behaviour of solutions appeared to be scaled in time, and this is evidencethat the SP system does not permit blow-up for this family of data.The SP system was modified by changing the coupling to gravity, namely by replacing |Ψ|2 inthe Poisson equation with |Ψ|epsilon1. Investigations were carried out with epsilon1 = 3, using the same familyof initial data. Critical behaviour was readily established in this case by varying the amplitudeof initial data. The critical value of the amplitude occurred between A = 75.3 and A = 75.4.The critical value was found to greater precision using a bisection search, which provided about15 digits of certainty. For values of A tuned close to the critical value, Astar, it was found thatsolutions quickly tended towards a critical solution before ultimately dispersing for subcriticaldata or becoming singular for supercritical data. The critical solution was found to be periodicin time, and that the lifetime of the solutions scaled with ln|A−Astar|, as is characteristic of TypeI critical phenomena.At this stage, we have not ruled out the possibility of Type II critical phenomena, and a moreextensive parameter space survey would need to be carried out, using different families of initialdata, to establish a possible phase transition between Type I and Type II solutions. This wouldbe of great interest since Type II phenomena permit mass scaling rather than just scaling of thelifetime.325.2. Future Work5.2 Future WorkAt the time of this writing, there remains many possibilities for future work. The most importanttask would be to determine the minimum value of epsilon1 for which blow-up occurs. This could presentsomedifficultiesbecausewehavenotbeenabletoprovethatblow-updoesnotoccurforsufficientlysmall epsilon1. Also, it would be straightforward to acquire tighter bounds on Astar and the Lyapunovexponent λ by doing more expensive computations at higher resolutions.Since it has been established that a critical solution exists for the case when epsilon1 = 3, it would beinteresting to attempt to construct the critical solution numerically. Then it would be possible toget a better prediction of the period of oscillation. Also, if other critical solutions are found whichexhibit some form of self-similarity, it would be necessary to implement a change of variables toresolve the solution over different scales. This feature has already been mostly implemented intothe code, and would not be difficult to implement fully.It would also be of interest to determine if the singularities are stable beyond spherical sym-metry, and thus one future possibility is to solve these systems imposing cylindrical symmetry orother to see if results similar to the spherical case are found.In order to observe the evolution of the blow-up solutions in more detail, it will be requiredto implement an adaptive mesh refinement scheme. Furthermore, the way we treat the outerboundary can be computationally taxing, because it requires more points as we increase theevolution time. A more sophisticated treatment at the outer boundary may prove useful infuture calculations.33Bibliography[1] V. A. Baskakov and A. V. Popov. Wave Motion, 14(123), 1991.[2] B. J. Carr and A. A. Coley. Class. Quantum Grav., 16(R31), 1999.[3] Matthew W. Choptuik. Phys. Rev. Lett., 70(9), 1993.[4] Matthew W. Choptuik. Lectures for vii mexican school on gravitation and mathematicalphysics: Relativistic astrophysics and numerical relativity: Numerical analysis for numericalrelativists, 2006.[5] Matthew W. Choptuik and David W. Neilsen. Class. Quantum Grav., 17:761–782, 1993.[6] C. R. Evans and J. S. Coleman. Phys. Rev. Lett., 72(1728), 1994.[7] R. L. Guenther, 1995.[8] Carsten Gundlach and Jos´e M. Mart´ın-Garc´ıa. Living Rev. Rel., 10(5), 2007.[9] T. Harada, H. Maeda, and B. Semelin. Phys. Rev. Lett. D, 67(084003), 2003.[10] S. W. Hawking and G. F. R. Ellis. The large-scale structure of spacetime. Cambrige Uni-versity Press, 1973.[11] Steven W. Hawking and W. Israel. Cambridge University Press, 1979.[12] M. Israeli and S. A. Orszag. J. Comput. Phy., 41(115), 1981.[13] Irene M. Moroz, Roger Penrose, and Paul Tod. Class. Quantum Grav., 15:2733 – 2742, 1998.[14] Roger Penrose. Revistas del Nuovo Cimento, 1(252), 1969.34Appendix APerturbation AnalysisIn this appendix we presenta perturbation analysis of critical phenomena in gravitational collapse,and distinguish between Type I and Type II transitions. The mathematical ideas presented inthis section are based on those from Gundlach[8].A.1 Type I Critical PhenomenaIn investigating critical collapse, one is concerned with a function Z which represents the strengthof the gravitational field, and is in general a function of both space and time, Z = Z(t,r). Zcan also depend on any number of parameters depending on the physical system in question. Inthe phase space of parameters, there will be a surface defining initial Z which will tend towardsthe critical solution. Varying a single parameter, p, holding all others fixed will produce initialconfigurations Z which generically will have an intersection with this surface. Fig.1.1 shows aone parameter family of initial configurations, and their resulting “trajectories” in phase space,ultimately resulting in one of three possible outcomes: singularity formation, total dispersion, ora critical solution.The critical solution in a Type I critical phenomenon is observed to be either independent oftime or periodic in time. Here we focus on the case where the solution is independent of time,and so we can write Zstar = Zstar(r), and its general linear perturbation can be written asZ(r,t) = Zstar(r)+XiCi(p)eλitZi(r) (A.1)Numerical simulations suggest that when looking at a 1-parameter family of initial configurations,there is only one growing mode (λ0 > 0), and so for Z close to the critical solution we have theapproximationZ similarequal Zstar(r)+ dC0dp (pstar)(p−pstar)eλ0tZ0(r) (A.2)Initial configurations Z will all approach the critical solution before the growing mode forces35A.2. Type II Critical Phenomenathem into one of the two stable configurations. The time spent near the critical solution will belonger for initial configurations with |p − pstar| closer to zero. This can be made quantitative indefining a time, tp bydC0dp |p−p∗|eλ0tp ≡ epsilon1 (A.3)where epsilon1 is an arbitrarily small positive constant. By setting epsilon1 to a particular small value thiscauses tp to be interpreted as the time at which Z satisfiesZ(r,tp) similarequal Zstar(r)±epsilon1Z0(r) (A.4)where Z0(r) is the initial configuration of gravity, and the sign in front of epsilon1 is the sign of |p−pstar|.From its definition we can see that tp scales astp = − 1λ0ln|p−pstar|+const (A.5)tp can be thought of as the time up to which the solution Z(t,r) is significantly dependent on|p−pstar|. After t = tp the solution will quickly converge to a stable solution and the informationof the magnitude of |p−pstar| is lost.It is also a matter of observation through numerical simulations that the mass of the stablesingular solution is a certain fraction of the mass of the critical solution, independent of initialdata. Therefore, Type I critical collapse usually occurs in physical systems where a mass scale isset.A.2 Type II Critical PhenomenaType II critical phenomena are observed to have critical solutions which exhibit scale-invariance,or self-similarity. This symmetry comes in a continuous and discrete form, as it did in Type Iphenomena, and to facilitate this symmetry we need to introduce new variables, for examplex = −rt, τ = −ln(−tl), t < 0 (A.6)where the origin of t has been moved so t = 0 and τ → ∞ corresponds to the time of criticalcollapse in the critical case, and l is a dimensionful parameter with units length which dependson the 1-parameter family of initial data.If we now define gravity configurations by Z = Z(τ,x), and the critical solution is Zstar(τ,x),then we can see that Zstar(τ,x) is continuously scale-invariant if it is independent of τ, and this36A.2. Type II Critical Phenomenaproperty is called continuous self-similarity (CSS). The discrete version is whenZ∗(τ,x) is periodicin τ withZstar(τ,x) = Zstar(τ +∆,x) (A.7)and this is called discrete self-similarity (DSS), because going in from τ to τ +∆, the solution isthe same but the space and time scales (x and t) are smaller by a factor of e−∆.Working through the same perturbation analysis as with Type I, we arrive at an analogousrelation for time scalingτp = 1λ0ln|p−pstar|+const (A.8)Where τ represents a measure of the time spent near the critical solution. In analogy with Eq.A.4,and shifting the origin of t to correspond with τ = τp, the solution can be written asZ(0,r) = Zstar(− rLp)±epsilon1Z0(− rLp) (A.9)where Lp = le−τp, which now represents an overall scale.Now, in the units we are working in, c = G = 1, mass also has units of length, which meansthe mass of the singularity is proportional to LpM ∝ Lp ∝ (p−pstar)1/λ0 (A.10)where the last proportionality follows from the definition of τ, Eq.A.8, and we have found the socalled critical exponent γ = 1/λ037Appendix BIndependent Residual CalculationThe independent residual calculation is a method to confirm that the solutions to the FDA thatthe code is producing converge to the continuum solutions in the limit that h → 0. This is doneby devising an alternative FDA defined on the same mesh, figure 2.3. This can be written as:Alternative Residual = ˜L[˜Ψ] (B.1)where ˜L represents the finite difference operators used in the new FDA and ˜Ψ is the solution tothe alternative FDA.The solutions to both FDA’s should be approximations to the continuum solutions of the SPsystem, and in the limit as h → 0 it is expected that solutions to both FDA will approach thecontinuum solutions. This means that if we take the solutions to the first FDA, and operate onit with the alternative FDA operator ˜L, then in the limit as h → 0 the result should tend to zero,ie:limh→0˜L[Ψ] = 0 (B.2)The idea is to choose an ˜L which is simpler in form than the original FDA (the Crank-Nicolsonscheme in our calculations). Of course the increased simplicity of the FDA usually corresponds toreduced accuracy, which is why the alternative FDA is not used to solve the SP system originally.However, allweareinterestedinisobservingthebehaviourofeq. B.2ash → 0, andthisbehaviourcan be observed by simply calculating the solution Ψ at several resolutions. An example of thisbehaviour is shown in figure B.1.For the alternative FDA, we choose to centre the equations at the full time level (as opposedto the half time level as is done with the Crank-Nicolson scheme). This requires the use of a FDAoperator for the time derivative which is O(dt2) accurate at the full time level. This is given by:Ψn+1j −Ψn−1j2dt =∂Ψnj∂t +O(dt2) (B.3)38Appendix B. Independent Residual CalculationFigure B.1: The real part of the alternative residual calculated for the “modified SP system” ie.epsilon1 = 3, for subcritical data. The residual is calculated at the base level of discretization, as wellas the first and second level, with 1025 points at base level. Clearly as the level of discretizationincreases, the residual is tending to zero which is the desired behaviour.39Appendix B. Independent Residual CalculationTherefore the operator ˜L is given by˜L[Ψnj] = i‡Ψn+1j −Ψn−1j·2dt+ 32drr2j+12‡Ψnj+1 −Ψnj·−r2j−12‡Ψnj −Ψnj−1·r3j+12−r3j−12−gV nj Ψnj (B.4)Notice that this choice of alternative FDA couples data at three time levels, whereas the originalFDA only coupled data at two time levels. This complicates things slightly, but is manageable.40

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
United States 1 0
City Views Downloads
Ashburn 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

Share

Embed

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"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.52966.1-0085959/manifest

Comment

Related Items