Critical Collapse of a Massive Scalar Field in Newtonian Gravity Numerical Investigations of the Schödinger-Poisson System by Andrew Inwood A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF BACHELOR OF SCIENCE (HONOURS) in The Faculty of Science (Physics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) May 2009 c© Andrew Inwood 2009 Abstract Critical collapse in General Relativity (GR) has been studied extensively, both numerically and analytically, for a wide variety of matter models including massive and massless scalar fields coupled to gravity, as well as a variety of gas models. Comparatively less work has been done towards the investigation of critical collapse in Newtonian gravity, although some progress has been made using an isothermal gas model. In this paper we numerically investigate critical collapse of a massive scalar field in spherical symmetry, coupled to gravity in the Newtonian limit. The evolution of such a field is described by the Schrödinger-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 Gaussian data, no evidence for critical behaviour is found, and so the SP system is modified through the gravitational coupling in an attempt to induce critical behaviour. This modification is realized by making the exponent in the Poisson equation variable. Type I critical behaviour is found using the Gaussian initial data centred at the origin for an exponent ² = 3. By performing a bisection search the critical value of the amplitude, A, is found to 15 digits of precision. The critical solution is found to be periodic in time, characteristic of Type I phenomena, with a period of T ≈ 0.02. By surveying the 1-dimensional parameter space of A, the Lyapunov exponent associated with the growing mode in Type I phenomena is found to 3 significant digits to be λ ≈ 163. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1 Critical Collapse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Self-Similarity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.3 Cosmic Censorship . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.4 Newtonian Collapse . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 1.5 Blow-up in Newtonian Self-Gravitating Systems . . . . . . . . . . . . . . . . . . . 6 2 Schrödinger-Poisson System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1 Spherical Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 Finite Difference Approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.2 Tridiagonal Form . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3.3 Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.4 Norm Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.3.5 Error Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.6 Change of Variables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 iii Table of Contents 3 Results: Schrödinger-Poisson System . . . . . . . . . . . . . . . . . . . . . . . . . 19 4 Modified Schrödinger-Poisson System . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.1 Parameter Space Survey . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.2 Critical Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 5 Summary and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5.1 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5.2 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 Appendices A Perturbation Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 A.1 Type I Critical Phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 A.2 Type II Critical Phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 B Independent Residual Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 iv List of Figures 1.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 extending from the initial data represent trajectories in the phase space of possible Z. The red line is drawn to cross the singularity threshold at p = p?, and for this value of p the gravity configuration tends towards the critical solution. Image taken from Gundlach[8]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2.1 The first bound-state wavefunction (S is a scaled wavefunction, Ψ). The solution was calculated by Moroz, Penrose, and Tod[13] numerically using a Runge-Kutta NAG routine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2 The first excited state wavefunction. Each excited state has one more zero than the previous. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.3 Discretization of space-time domain. Grid functions (ie. V and Ψ) are labelled with subscripts 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 to h. 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. . . . . . . . . . . . . . . . . 12 2.4 Time evolution of the norm for the Schrödinger equation, calculated at multiple resolutions. The difference between the value of the norm at the initial time as calculated at levels 0 and 2 is about 1 part in 100,000. . . . . . . . . . . . . . . . . 18 3.1 For a small value of A, the self-gravitation is weak and the solutions immediately disperse. This means the probability density, |Ψ|2, becomes uniform as t→∞. . . 20 v List of Figures 3.2 The time derivative of |Ψ|2 at the origin is initially negative and tends to zero as t → ∞. This indicates that the solution is dispersing as t → ∞. The resolution for these calculations is 1025 points, unless indicated otherwise. . . . . . . . . . . . 21 3.3 When A = 800, the wavefunction initially grows at the origin before ultimately dispersing. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.4 For A = 800, the growth rate at the origin is initially positive, indicating that mass is accumulating at the origin. Ultimately the growth rate becomes negative and tends to zero indicating dispersion. . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.5 By increasing A, it is observed that the growth rate at the origin exhibits similar behaviour on shorter time scales. This is evidence that the SP system scales in time with A. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.6 The growth rate for A = 5000 exhibits the same sort of behaviour as for A = 1500 and A = 1000. For this calculation the resolution was increased to 2049 points. . . 24 4.1 For values of A below the critical value, the growth rate is initially positive but ultimately becomes negative and tends to zero as the solution disperses, just as for the SP system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 4.2 For A = 80, the growth rate makes a transition in behaviour, and appears to grow without bound. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.3 By varying A through the critical value, the transition is clearly seen through the growth rate at the origin. For A = 80 the growth rate diverges from the other lines very quickly. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.4 The same picture as Fig.4.3, but on a larger time scale. Here it can be seen that the growth rate also increases without bound for A = 75.5, and that the time at which blow-up occurs is almost twice as late for A = 80. . . . . . . . . . . . . . . . 28 4.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. Resolution is 32800 points on the range 0 < r < 40. . . . . . . . . . . . . . . . . . . . . . . . . 29 4.6 Lifetimes of subcritical solutions calculated at resolution of 32800 points. The scaling exponent σ is determined from the slope, σ ≈ 0.0061. Supercritical data should give similar results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 vi List of Figures B.1 The real part of the alternative residual calculated for the “modified SP system” ie. ² = 3, for subcritical data. The residual is calculated at the base level of discretization, 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 which is the desired behaviour. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39 vii Acknowledgements I would like to thank my supervisor, Dr. Matthew Choptuik, for his support in this project. My engaging conversations with Dr. Choptuik have caused me to develop a strong interest in the fields of Gravitation and Computational Physics, and I am certain that the knowledge I have gained from him will serve me well in my future endeavours. 1 Chapter 1 Introduction It has long been known that solutions to Einstein’s equations with physically reasonable matter in the initial state can develop singularities [10]. In fact, if one can construct a so called “trapping region” 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 interesting because it is a sign that the validity of the theory has broken down in this region. In order to fully 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 event horizon—a so called “clothed singularity”. This means that the singularity would not be causally connected to an observer lying outside of the event horizon, and so the predictions of GR are unaffected in most practical applications. However, beyond the singularity theorems, the details of singularity formation are not well understood, and it is not clear that generically singularities will be clothed, and might in fact form beyond the confines of an event horizon. Such a singularity is called a “naked singularity”. This is related to the topic of Cosmic Censorship, which will be discussed further. 1.1 Critical Collapse By 1973 it was clear that singularities must form in GR from physically reasonable and sufficiently massive initial states, and it was also clear that for sufficiently light initial states, singularity formation 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 mass of the initial state is neither light nor heavy. This, as well as the question of whether or not singularities are all clothed, motivated Choptuik[3] to investigate the middle ground numerically in the early 1990’s. Choptuik investigated a spherically symmetric, massless scalar field coupled to gravity in GR, 2 1.2. Self-Similarity and by varying 1-parameter families of initial data, investigated the transition between solutions in which all of the mass disperses away from the origin, and those in which the mass accumulates at 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 initial data—this property is known as the Universality of the critical solution. This behaviour of the solutions is analogous to phase transitions in statistical mechanics. There are two types of transitions, 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 naked singularity at the origin. 1.2 Self-Similarity Further numerical investigations on the critical behaviour of self-gravitating systems[6][5] revealed two different kinds of critical solutions: ones that are time-independent, and ones that exhibit either continuous self-similarity (CSS) or discrete self-similarity (DSS) (for a review of self-similar solutions 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 in statistical mechanics. For Type I critical phenomena the mass of the end state singular solution is a certain fraction of the mass of the critical solution, but for Type II the mass of the end state singular solution obeys a power scaling law, and initial conditions can be set to make the resulting black hole mass infinitesimally small. 1.3 Cosmic Censorship Cosmic Censorship encompasses the idea that even though singularities exist in GR, they may not have any influence on the physics outside of the event horizon which surrounds it. At present, any statement of cosmic censorship is at the level of a conjecture, because the available mathematical tools are insufficient to provide anything as far as a mathematical proof or theorem. There are two 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 be enclosed by an event horizon, and therefore the physics beyond the confines of the event horizon 3 1.3. Cosmic Censorship Figure 1.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 extending from the initial data represent trajectories in the phase space of possible Z. The red line is drawn to cross the singularity threshold at p = p?, and for this value of p the gravity configuration tends towards the critical solution. Image taken from Gundlach[8]. 4 1.4. Newtonian Collapse will be unaffected. This is the “weak” version because it still allows for the singularity to affect observers who travel through the event horizon, into the black hole. The strong cosmic censorship conjecture[11] claims that time-like singularities do not form, and therefore even observers who travel into the black hole will not be affected by the singularity. The “middle ground” self-gravitating regime is an excellent laboratory to investigate cosmic censorship, since naked singularities have been found in numerical studies of spherical critical collapse. However, since most evidence for naked singularity formation is empirical, it is not always clear whether or not the singularities are stable under perturbations to the model, such as evolution beyond spherical symmetry. The precise form of cosmic censorship is one of the biggest open questions in classical GR. Determining the precise form of cosmic censorship and understanding under what conditions self-similarity arises are two of the main goals of research in the field of gravitational collapse. 1.4 Newtonian Collapse Investigations 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 compared to GR; however, due to the reduced complexity associated with Newtonian gravity, it offers a good laboratory for investigating critical collapse under deviations from spherical symmetry and other simplifications that are usually commonplace. For these reasons, working with Newtonian gravity could offer insight into cosmic censorship and self-similarity. It should be noted that since these analyses are purely Newtonian, they could potentially have been carried out long before the 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 insight into the role that the dynamics play in creating self-similar solutions. 5 1.5. Blow-up in Newtonian Self-Gravitating Systems 1.5 Blow-up in Newtonian Self-Gravitating Systems Here we are interested in investigating critical collapse in Newtonian self-gravitating systems which exhibit “blow-up”, ie. solutions which become singular in finite time. Such solutions are interesting for a number of reasons. As mentioned above, investigating the evolution of solutions as 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 clear that critical phenomena will occur for any matter model. This is because in GR, the universality of coupling of mass, energy, and momentum guarantees that there will be black hole formation in any matter model; however, Newtonian gravity only sets a universality of coupling of mass to momentum, and it is not clear that all matter models will admit initial data that result in blow-up. 6 Chapter 2 Schrödinger-Poisson System 2.1 Background In 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 compared to the speed of light, the EKG equation can be shown to reduce to the Schrödinger-Poisson (SP) system[7]: i ∂Ψ ∂t = −1 2 ∇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 physical quantities are given by xphy = h̄ mc xcomp (2.3) tphy = h̄ mc2 tcomp (2.4) Ψphy = c2√ 4piG Ψcomp (2.5) The physical interpretation of the SP system is that the wavefunction of a single particle evolves according to the Schrödinger 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 solutions called Newtonian Boson Stars (whereas the static solutions to the EKG equations are called Relativistic Boson Stars), see figure 2.1. These solutions are smooth at the origin, and can be ordered according to the number of zeroes they have. It is not known as far as a mathematical proof whether or not this system permits solutions that exhibit blow-up, and so the first order of business in this research is to solve the SP system 7 2.1. Background Figure 2.1: The first bound-state wavefunction (S is a scaled wavefunction, Ψ). The solution was calculated by Moroz, Penrose, and Tod[13] numerically using a Runge-Kutta NAG routine 8 2.1. Background Figure 2.2: The first excited state wavefunction. Each excited state has one more zero than the previous. 9 2.2. Constraints numerically using finite-differencing techniques, and vary a 1-parameter family of initial data in an attempt to find blow-up solutions. 2.2 Constraints 2.2.1 Spherical Symmetry Working in spherical symmetry allows us to write the laplacian as: ∆f = 1 r2 ∂ ∂r ( r2 ∂f ∂r ) (2.6) However, smoothness at the origin will be better enforced by the numerical scheme if the laplacian is rewritten in terms of ∂(r3), ∆f = 3 ∂ ∂(r3) ( r2 ∂f ∂r ) (2.7) 2.2.2 Boundary Conditions At the inner boundary (near the origin) we imposed that V and Ψ are smooth: dV dr |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 as r goes to infinity. Although it is possible to write an equivalent boundary condition at a finite r [1], these are difficult to implement numerically. There are a number of options available for dealing 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 or not blow-up occurs in the SP system, so to simplify things we impose that all functions are zero at 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 the equations can be integrated, because if the wavefunction propagates out to the outer boundary it will reflect and the solution becomes unphysical. But again, if blow-up is occurring in the model, 10 2.3. Finite Difference Approximation then all the mass falls into the origin very quickly, faster than the wavefunction can propagate to the outer boundary. 2.3 Finite Difference Approximation 2.3.1 Discretization To discretize the spatial domain into Nr grid points, FORTRAN notation is used, ie: r1 = 0 r2 = h r3 = 2h ... rNr = 1 Figure 2.3 displays a schematic diagram of the discretized domain. h is the standard notation for the spacing of the spatial coordinate, and we use h and dr synonymously. Note that subscripts indicate 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+ 1 2 ≡ rj + dr2 (2.12) rj− 1 2 ≡ rj − dr2 this slight abuse of notation makes the equations more readable. Now the laplacian takes the form: ∆f = 3 r3 j+ 1 2 − r3 j− 1 2 ( r2 j+ 1 2 ∗ ( fj+1 − fj dr ) − r2 j− 1 2 ∗ ( fj − fj−1 dr )) (2.13) We choose to implement a Crank-Nicolson FDA. The essence of a Crank-Nicolson scheme is that the equations of motion are centred at the so called half time level. This is simple to do for the time-derivative of Ψ, by writing: Ψn+1j −Ψnj dt = ( ∂Ψ ∂t )n+ 1 2 j +O(h2) (2.14) 11 2.3. Finite Difference Approximation Figure 2.3: Discretization of space-time domain. Grid functions (ie. V and Ψ) are labelled with subscripts 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 to h. 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. 12 2.3. Finite Difference Approximation To centre the right-hand side of the Schrödinger equation at the half time level and maintain O(h2) accuracy, it is necessary to average the Ψn and Ψn+1 time levels. Taking into account all of these considerations, and using the discretized laplacian, the Schrödinger equation becomes: i ( Ψn+1j −Ψnj ) dt = −1 2 · 3 r3 j+ 1 2 − r3 j− 1 2 [ r2 j+ 1 2 · 1 2 ( Ψn+1j+1 −Ψn+1j +Ψnj+1 −Ψnj ) dr (2.15) − r2 j− 1 2 · 1 2 ( Ψn+1j −Ψn+1j−1 +Ψnj −Ψnj−1 ) dr ] + g 2 ( V n+1j Ψ n+1 j + V n j Ψ n j ) and the FDA of the Poisson equation becomes |Ψnj |2 = 3 r3 j+ 1 2 − r3 j− 1 2 [ r2 j+ 1 2 · ( V nj+1 − V nj dr ) − r2 j− 1 2 · ( V nj − V nj−1 dr )] (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 boundary can be treated by introducing the O(h2) accurate forward-differencing derivative operator: (fr)nj = −32fnj + 2fnj+1 − 12fnj+2 h +O(h2) (2.17) Using this operator, the inner boundary condition can be rewritten as: fn1 = − 1 3 fn3 + 4 3 fn2 (2.18) where f is either V or Ψ. This essentially allows for all references to V1 and Ψ1 to be removed from 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 Form As it turns out, both eqs. 2.15 and 2.16 can be brought into a canonical form called tridiagonal form. The Poisson equation in tridiagonal form looks like: c−V n+1j−1 + c0V n+1 j + c+V n+1 j+1 = ( ΨjΨ∗j )n+1 (2.19) 13 2.3. Finite Difference Approximation where c− = 3 dr  r2j− 12 r3 j+ 1 2 − r3 j− 1 2  (2.20) c0 = − 3 dr r2j− 12 + r2j+ 12 r3 j+ 1 2 − r3 j− 1 2  (2.21) c+ = 3 dr  r2j+ 12 r3 j+ 1 2 − r3 j− 1 2  (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 = 4 3 V2 − 13V3 (2.24) Then the equations at the boundary become:( 4 3 c− + c0 ) V2 + ( c+ − 13c− ) V3 = Ψ(2)Ψ∗(2) (2.25) Now all references to V1 have been removed, and we are solving a smaller system for V2, . . . , VNr. Note that having removed the reference to V1, the Poisson equation, eq. 2.19, can be written in matrix form as:  4 3c− + c0 c+ − 13c− 0 0 0 0 c− c0 c+ 0 0 0 . . . 0 c− c0 c+ 0 0 0 0 c− c0 c+ 0 ...   V2 V3 V4 V5 V6 V7 ...  =  |Ψ2|2 |Ψ3|2 |Ψ4|2 |Ψ5|2 |Ψ6|2 |Ψ7|2 ...  (2.26) Where it is understood that c−, c0, and c+ are all dependent on r. Only the three main diagonals of the matrix operating on V are non-zero, and this is why this form is called tridiagonal form. The Schrödinger equation can also be brought into tridiagonal form, although this equation will depend on values of V n+1 and also the nth time level V n and Ψn. We just wrap the extra 14 2.3. Finite Difference Approximation dependence in a function F , and then the tridiagonal form is: d−Ψn+1j−1 + d0Ψ n+1 j + d+Ψ n+1 j+1 = F (j) (2.27) where d− = 1 2 3 r3 j+ 1 2 − r3 j− 1 2 r2 j− 1 2 · 1 2 dt dr (2.28) d0 = i− 12 3 r3 j+ 1 2 − r3 j− 1 2 [ r2 j+ 1 2 · 1 2 dt dr + r2 j− 1 2 · 1 2 dt dr ] − g 2 V n+1j (2.29) d+ = 1 2 3 r3 j+ 1 2 − r3 j− 1 2 r2 j+ 1 2 · 1 2 dt dr (2.30) F (j) = iΨnj − 1 2 3 r3 j+ 1 2 − r3 j− 1 2 [ r2 j+ 1 2 · 1 2 dt dr ( Ψnj+1 −Ψnj )− r2 j− 1 2 · 1 2 dt dr ( Ψnj −Ψnj−1 ) ] + g 2 V nj Ψ n j (2.31) And again the equation near the inner boundary can be written as:( 4 3 d− + d0 ) Ψ2 + ( d+ − 13d− ) Ψ3 = F (2) (2.32) 2.3.3 Code We have now developed all of the machinery to take an initial state Ψ1 and calculate the state at the next time level Ψ2 satisfying the SP system and boundary conditions. This calculation has been reduced to inverting tridiagonal matrices. However, since the Schrödinger and Poisson equations are coupled, an initial guess at V n+1 is required to invert the matrix for Ψ. Then V n+1 and Ψn+1 are updated iteratively until both are within tolerance. Routines have been developed to quickly and efficiently invert tridiagonal matrices, and we use the routines found in LAPACK1 for tridiagonal matrices: DGTSV and ZGTSV. A skeleton code 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 Calculation As usual, the physical interpretation of |Ψ|2 is that it is the spatial probability distribution function of the particle. This means that integrating over it tells us the probability of finding the 1LAPACK (Linear Algebra PACKage): http://www.netlib.org/lapack/ 15 2.3. Finite Difference Approximation particle in a region of space. However, since we are working in spherical symmetry, there is an r2 weighting in this probability, so even though |Ψ|2 may be large at the origin, the probability of 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 the integrated probability: I(r) = 4pi ∫ r 0 ξ2|Ψ(ξ)|2 dξ (2.33) If the equations of motion have been set up such that Ψ is exhibiting blow-up behaviour, then one could tell that most of the mass is concentrated near the origin if eq. 2.33 goes to 1 for a small value of r. To calculate this integral numerically we employ a Riemann sum using the midpoint rule, I(rj) = I(rj−1) + 4pidr ( r(j)− dr 2 )2 (Ψj +Ψj−1) 2 ( Ψ∗j +Ψ ∗ j−1 ) 2 (2.34) and I(r1) ≡ 0. 2.3.5 Error Analysis The ansatz underlying all of our error analysis, originally observed by Richardson, is that the solution, Ψ, 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 FDA using three different values of h, in the ratio 1 : 2 : 4. Also noting that for centred FDA, like the Crank-Nicolson scheme, the odd powers of h in the Richardson expansion (equation 2.35) drop out, 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) 16 2.3. Finite Difference Approximation Where || · || is any suitable spatial norm. For the case of the SP system in spherical symmetry, Q(t) takes the form: Q(t) = ∫ rmax 0 r 2(Ψ4h −Ψ2h) dr∫ rmax 0 r 2(Ψ2h −Ψh) dr (2.40) It can be easily shown that if the FDA is converging then: lim h→0 Q(t) = 4 (2.41) For a more detailed discussion of convergence testing based on the Richardson expansion for FDAs, see the notes by Choptuik[4]. To check that the PDE had been coded properly, it is useful to use the norm of Ψ as a diagnostic. For the Schrödinger 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 space is always unity. Checking for norm conservation is a way to check that the code is producing sensible results. Figure 2.4 displays the time evolution of the norm of Ψ for a free space Schrödinger equation, calculated at multiple resolutions. The plots are converging to a constant, which is the expected behaviour for the Schrödinger equation. The norm calculated at different resolutions do not agree at the initial time due to the fact that the finite mesh spacing introduces uncertainty into the norm calculation; however, at the initial time the norm as calculated at level 0 and level 2 vary by only about 1 part in 100,000. The norm is so closely conserved for the Schrödinger equation because the Crank-Nicolson FDA itself also conserves the norm, and so the norm is conserved to greater precision than might be expected by an analysis involving only the Richardson expansion. This poses a problem in using 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 it have to be compared to those produced by an alternative FDA (ideally one which is not norm preserving). Such a comparison is called an alternative residual calculation, the implementation of which is outlined in appendix B. 17 2.3. Finite Difference Approximation Figure 2.4: Time evolution of the norm for the Schrödinger equation, calculated at multiple resolutions. The difference between the value of the norm at the initial time as calculated at levels 0 and 2 is about 1 part in 100,000. 2.3.6 Change of Variables It should be noted that if self-similarity is found in the process of these investigations, it will become necessary to solve the equations under arbitrary change of variables, depending on the form of the self-similarity. This has mostly been implemented into the code, although at the time of writing this feature has not been required, as self-similarity has not been found. 18 Chapter 3 Results: Schrödinger-Poisson System Investigations were carried out for a 1-parameter family of initial data in hopes of finding critical behaviour 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 we are interested in investigating which matter models permit blow-up, and so we choose to fix the coupling constant g equal to 1, and vary the amplitude of the Gaussian, since this corresponds to 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 the system is not significant, and the solutions behave like solutions to the free space Schrödinger equation. The mass disperses from the origin, ultimately evolving towards a state of equal probability on the interval 0 ≤ r ≤ 1. Figure 3.2 displays the time derivative of |Ψ|2 evaluated at the origin. For A = 1, this rate is observed to be negative at the initial time, and goes to zero in the limit as t → ∞. If blow-up were to occur, this rate should be observed to grow without bound, which is why it is useful to look at these plots. If A is set much larger than 1, a different behaviour is observed. For A = 800, figure 3.3, the self-gravity of the system is now strong enough that initially the mass starts to accumulate at the origin, 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 the mass starts to disperse away from the origin. Based on the observed change of behaviour as A is increased, it might be expected that for sufficiently large A the self-gravitation will be strong enough 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 is seen 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 growth 19 Chapter 3. Results: Schrödinger-Poisson System Figure 3.1: For a small value of A, the self-gravitation is weak and the solutions immediately disperse. This means the probability density, |Ψ|2, becomes uniform as t→∞. rate seems to scale with A. 20 Chapter 3. Results: Schrödinger-Poisson System Figure 3.2: The time derivative of |Ψ|2 at the origin is initially negative and tends to zero as t → ∞. This indicates that the solution is dispersing as t → ∞. The resolution for these calculations is 1025 points, unless indicated otherwise. 21 Chapter 3. Results: Schrödinger-Poisson System Figure 3.3: When A = 800, the wavefunction initially grows at the origin before ultimately dispersing. 22 Chapter 3. Results: Schrödinger-Poisson System Figure 3.4: For A = 800, the growth rate at the origin is initially positive, indicating that mass is accumulating at the origin. Ultimately the growth rate becomes negative and tends to zero indicating dispersion. Figure 3.5: By increasing A, it is observed that the growth rate at the origin exhibits similar behaviour on shorter time scales. This is evidence that the SP system scales in time with A. 23 Chapter 3. Results: Schrödinger-Poisson System Figure 3.6: The growth rate for A = 5000 exhibits the same sort of behaviour as for A = 1500 and A = 1000. For this calculation the resolution was increased to 2049 points. 24 Chapter 4 Modified Schrödinger-Poisson System While it has not been proven that blow-up cannot occur in the SP system, we should keep in mind that the main goal of this research is to investigate blow-up in Newtonian self-gravitating systems. Since there is currently no evidence that the SP system will permit blow-up, in order to proceed with the research we modify the SP system in a way that is expected to permit blow-up solutions. It should be noted that the modification is only motivated by our interest in producing blow-up solutions, and not motivated by an attempt to model a known physical system. In this research we hope to find what kind of matter models will permit blow-up, and whether or not these 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 gravitational collapse. For example, one could try introducing a Ψ|Ψ|2 term into the Schrödinger equation to produce the usual non-linear Schrödinger (NLS) equation, but now coupled to a Poisson equation; however, this is guaranteed to produce blow-up solutions for the focussing NLS equation, and does not result from the gravitational coupling. In this research, it is desired to investigate blow-up which occurs as a result of the gravitational coupling, and so we choose to make the following modification to the SP system: i ∂Ψ ∂t = −1 2 ∆Ψ+ gVΨ (4.1) ∆V = |Ψ|² where ² ≥ 2. The only difference here is the ² in the Poisson equation. All of the machinery developed in the earlier chapters can be recycled, including the FDA (with trivial modification of the discretized Poisson equation) and independent residual checker, and so will not be reproduced here. 25 4.1. Parameter Space Survey Figure 4.1: For values of A below the critical value, the growth rate is initially positive but ultimately becomes negative and tends to zero as the solution disperses, just as for the SP system. 4.1 Parameter Space Survey Here we present results for when ² = 3. Again we set g = 1 and investigate the same 1-parameter family 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, where the growth rate at the origin is initially positive as the mass accumulates at the origin, but ultimately becomes negative and tends to zero as the mass disperses, see figure 4.1. However, for A = 80, figure 4.2, a very different behaviour arises, where the growth rate remains positive for all 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 grows without bound. In fact, there is a change in behaviour between A = 75.3 and A = 75.5. For A = 75.5, the growth rate tends to infinity, although this is not apparent in figure 4.3. In figure 4.4 it is clear that for A = 75.5, the growth rate eventually increases rapidly, although this is at 26 4.1. Parameter Space Survey Figure 4.2: For A = 80, the growth rate makes a transition in behaviour, and appears to grow without bound. a much later time than for A = 80. This is clear evidence of critical behaviour, for presumably if A is tuned arbitrarily close to the critical value, then the growth rate will start to increase at arbitrarily large times. This agrees with the standard picture of Type I critical collapse as described in the appendices, where the solutions tend towards a so called critical solution before “deciding” between a singular or dispersive end state. 27 4.1. Parameter Space Survey Figure 4.3: By varying A through the critical value, the transition is clearly seen through the growth 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 that the growth rate also increases without bound for A = 75.5, and that the time at which blow-up occurs is almost twice as late for A = 80. 28 4.2. Critical Solution 4.2 Critical Solution To get the tightest bound on the critical value of A, a bisection search was performed to acquire A? to 15 digits of precision. The calculations were performed at higher resolutions to check the convergence of A?, and at the time of this writing, based on the convergence as we increase the resolution, A? has been determined to 3 significant digits. The closer we tune A to the critical value, the longer the solution will persist before dispersing or becoming singular. This poses a problem because reflections off the outer boundary limits the time over which the solutions can evolve 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 on the range 0 < r < 40. Figure 4.5 shows the time evolution of the point at r = 0, using a resolution of 32800 points on the interval 0 < r < 40. The solution is oscillatory, with a period of T ≈ 0.02. The periodic nature of the solution is also typical of Type I critical behaviour. As is described in the appendices, the 29 4.2. Critical Solution natural parameterization, pi, of Type I phenomena is given by: pi = ln |A−A?| (4.2) This form is natural because it is found that the “lifetime”, τ , of the solution, the time up to which the solution is “close” to the critical solution, is proportional to pi in the limit that the mesh spacing goes to zero. τ → σpi (4.3) To find σ we calculate the slope of a τ vs pi plot, and check the convergence of σ with increasing resolution. In order to do this we need a working definition of τ . We define τ as the time at which there is a 1% fractional difference between the solution and the critical solution at the origin, ie: ||Ψ|2(A)− |Ψ|2(A?)| |Ψ|2(A?) |r=0 > 0.01 (4.4) Of course we only know A? to a certain precision set by our bisection search, and so when looking at subcritical data we use the lower bound on A? produced by the bisection search, and when looking at supercritical data we use the upper bound. The fact that we do not know A? to infinite precision means that pi has more uncertainty when pi → −∞, and so when calculating σ we have to be careful to look at values of pi such that the effects of approximating A? are small. Figure 4.6 shows a plot of τ vs pi, using subcritical data. By calculating the slope using data produced at different resolutions, we can see that σ has been found accurately to 2 significant digits, with σ ≈ 0.0061. Referring to the perturbation analysis in the appendices, σ is related to the Lyapunov exponent, λ, associated with the growing mode in Type I solutions, as λ = −1/σ. From this we see that λ ≈ 163. 30 4.2. Critical Solution Figure 4.6: Lifetimes of subcritical solutions calculated at resolution of 32800 points. The scaling exponent σ is determined from the slope, σ ≈ 0.0061. Supercritical data should give similar results. 31 Chapter 5 Summary and Future Work 5.1 Summary Partially motivated by work on a Newtonian isothermal gas model [9], the Schrödinger-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 system was investigated numerically using a Crank-Nicolson finite-difference approximation (FDA), and a family of 1-parameter initial data was surveyed to look for blow-up. As the parameter was increased, the general behaviour of solutions appeared to be scaled in time, and this is evidence that 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 in the Poisson equation with |Ψ|². Investigations were carried out with ² = 3, using the same family of initial data. Critical behaviour was readily established in this case by varying the amplitude of 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 about 15 digits of certainty. For values of A tuned close to the critical value, A?, it was found that solutions quickly tended towards a critical solution before ultimately dispersing for subcritical data or becoming singular for supercritical data. The critical solution was found to be periodic in time, and that the lifetime of the solutions scaled with ln |A−A?|, as is characteristic of Type I critical phenomena. At this stage, we have not ruled out the possibility of Type II critical phenomena, and a more extensive parameter space survey would need to be carried out, using different families of initial data, to establish a possible phase transition between Type I and Type II solutions. This would be of great interest since Type II phenomena permit mass scaling rather than just scaling of the lifetime. 32 5.2. Future Work 5.2 Future Work At the time of this writing, there remains many possibilities for future work. The most important task would be to determine the minimum value of ² for which blow-up occurs. This could present some difficulties because we have not been able to prove that blow-up does not occur for sufficiently small ². Also, it would be straightforward to acquire tighter bounds on A? and the Lyapunov exponent λ by doing more expensive computations at higher resolutions. Since it has been established that a critical solution exists for the case when ² = 3, it would be interesting to attempt to construct the critical solution numerically. Then it would be possible to get a better prediction of the period of oscillation. Also, if other critical solutions are found which exhibit some form of self-similarity, it would be necessary to implement a change of variables to resolve the solution over different scales. This feature has already been mostly implemented into the 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 or other 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 required to implement an adaptive mesh refinement scheme. Furthermore, the way we treat the outer boundary can be computationally taxing, because it requires more points as we increase the evolution time. A more sophisticated treatment at the outer boundary may prove useful in future calculations. 33 Bibliography [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 mathematical physics: Relativistic astrophysics and numerical relativity: Numerical analysis for numerical relativists, 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é M. Mart́ın-Garćı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. 34 Appendix A Perturbation Analysis In this appendix we present a perturbation analysis of critical phenomena in gravitational collapse, and distinguish between Type I and Type II transitions. The mathematical ideas presented in this section are based on those from Gundlach[8]. A.1 Type I Critical Phenomena In investigating critical collapse, one is concerned with a function Z which represents the strength of the gravitational field, and is in general a function of both space and time, Z = Z(t, r). Z can also depend on any number of parameters depending on the physical system in question. In the phase space of parameters, there will be a surface defining initial Z which will tend towards the critical solution. Varying a single parameter, p, holding all others fixed will produce initial configurations Z which generically will have an intersection with this surface. Fig.1.1 shows a one parameter family of initial configurations, and their resulting “trajectories” in phase space, ultimately resulting in one of three possible outcomes: singularity formation, total dispersion, or a critical solution. The critical solution in a Type I critical phenomenon is observed to be either independent of time or periodic in time. Here we focus on the case where the solution is independent of time, and so we can write Z? = Z?(r), and its general linear perturbation can be written as Z(r, t) = Z?(r) + ∑ i Ci(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 the approximation Z ' Z?(r) + dC0 dp (p?)(p− p?)eλ0tZ0(r) (A.2) Initial configurations Z will all approach the critical solution before the growing mode forces 35 A.2. Type II Critical Phenomena them into one of the two stable configurations. The time spent near the critical solution will be longer for initial configurations with |p − p?| closer to zero. This can be made quantitative in defining a time, tp by dC0 dp |p− p∗|eλ0tp ≡ ² (A.3) where ² is an arbitrarily small positive constant. By setting ² to a particular small value this causes tp to be interpreted as the time at which Z satisfies Z(r, tp) ' Z?(r)± ²Z0(r) (A.4) where Z0(r) is the initial configuration of gravity, and the sign in front of ² is the sign of |p− p?|. From its definition we can see that tp scales as tp = − 1 λ0 ln |p− p?|+ const (A.5) tp can be thought of as the time up to which the solution Z(t, r) is significantly dependent on |p− p?|. After t = tp the solution will quickly converge to a stable solution and the information of the magnitude of |p− p?| is lost. It is also a matter of observation through numerical simulations that the mass of the stable singular solution is a certain fraction of the mass of the critical solution, independent of initial data. Therefore, Type I critical collapse usually occurs in physical systems where a mass scale is set. A.2 Type II Critical Phenomena Type 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 I phenomena, and to facilitate this symmetry we need to introduce new variables, for example x = −r t , τ = − ln(− t l ), t < 0 (A.6) where the origin of t has been moved so t = 0 and τ → ∞ corresponds to the time of critical collapse in the critical case, and l is a dimensionful parameter with units length which depends on the 1-parameter family of initial data. If we now define gravity configurations by Z = Z(τ, x), and the critical solution is Z?(τ, x), then we can see that Z?(τ, x) is continuously scale-invariant if it is independent of τ , and this 36 A.2. Type II Critical Phenomena property is called continuous self-similarity (CSS). The discrete version is when Z∗(τ, x) is periodic in τ with Z?(τ, x) = Z?(τ +∆, x) (A.7) and this is called discrete self-similarity (DSS), because going in from τ to τ +∆, the solution is the 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 analogous relation for time scaling τp = 1 λ0 ln |p− p?|+ 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 as Z(0, r) = Z?(− r Lp )± ²Z0(− r Lp ) (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 means the mass of the singularity is proportional to Lp M ∝ Lp ∝ (p− p?)1/λ0 (A.10) where the last proportionality follows from the definition of τ , Eq.A.8, and we have found the so called critical exponent γ = 1/λ0 37 Appendix B Independent Residual Calculation The independent residual calculation is a method to confirm that the solutions to the FDA that the code is producing converge to the continuum solutions in the limit that h→ 0. This is done by 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 to the alternative FDA. The solutions to both FDA’s should be approximations to the continuum solutions of the SP system, and in the limit as h → 0 it is expected that solutions to both FDA will approach the continuum solutions. This means that if we take the solutions to the first FDA, and operate on it with the alternative FDA operator L̃, then in the limit as h→ 0 the result should tend to zero, ie: lim h→0 L̃[Ψ] = 0 (B.2) The idea is to choose an L̃ which is simpler in form than the original FDA (the Crank-Nicolson scheme in our calculations). Of course the increased simplicity of the FDA usually corresponds to reduced accuracy, which is why the alternative FDA is not used to solve the SP system originally. However, all we are interested in is observing the behaviour of eq. B.2 as h→ 0, and this behaviour can be observed by simply calculating the solution Ψ at several resolutions. An example of this behaviour is shown in figure B.1. For the alternative FDA, we choose to centre the equations at the full time level (as opposed to the half time level as is done with the Crank-Nicolson scheme). This requires the use of a FDA operator for the time derivative which is O(dt2) accurate at the full time level. This is given by: Ψn+1j −Ψn−1j 2dt = ∂Ψnj ∂t +O(dt2) (B.3) 38 Appendix B. Independent Residual Calculation Figure B.1: The real part of the alternative residual calculated for the “modified SP system” ie. ² = 3, for subcritical data. The residual is calculated at the base level of discretization, 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 which is the desired behaviour. 39 Appendix B. Independent Residual Calculation Therefore the operator L̃ is given by L̃[Ψnj ] = i ( Ψn+1j −Ψn−1j ) 2dt + 3 2dr r2 j+ 1 2 ( Ψnj+1 −Ψnj ) − r2 j− 1 2 ( Ψnj −Ψnj−1 ) r3 j+ 1 2 − r3 j− 1 2 − gV nj Ψnj (B.4) Notice that this choice of alternative FDA couples data at three time levels, whereas the original FDA only coupled data at two time levels. This complicates things slightly, but is manageable. 40