Gravitational collapse with quantum fields

Gravitational collapse into a black hole has been extensively studied with classical sources. We develop a new formalism to simulate quantum fields forming a black hole. This formalism utilizes well-established techniques used for classical collapse by choosing a convenient coherent state, and simulates the matter fields quantum mechanically. Divergences are regularized with the cosmological constant and Pauli-Villars fields. Using a massless spherically symmetric scalar field as an example, we demonstrate the effectiveness of the formalism by reproducing some classical results in gravitational collapse, and identifying the difference due to quantum effects.


Introduction
The simplest model of a black hole forming in 4D is the collapse of a massless scalar field. The classical aspects of scalar collapse have been extensively studied in spherical symmetry analytically by Christodoulou [1][2][3][4] and numerically in [5] and also [6] where Choptuik famously uncovered the phenomenon of critical collapse that exhibits a range of fascinating properties [7][8][9][10].
The semi-classical analysis by Hawking [11] shed new light on the quantum behavior of black holes, discovering the celebrated Hawking radiation. However, non-linearly solving scalar collapse semi-classically is a tremendous task, as the quantum evolution of matter fields with full gravitational effects is difficult to solve even in spherical symmetry. Early attempts towards investigating quantum scalar collapse were made in [12][13][14], as well as studying the quantum behaviour of a collapsing spherical shell [15,16], and the modelling of the quantum radiation with classical fields [17]. Numerical evidence for Hawking radiation has been found in analogue systems [18,19]. Simpler models such as the collapse of a dilaton in 2D gravity have been studied [20][21][22][23], but these do not contain the spherical harmonics that originate in the 4D system.
In this letter, we establish a new method to numerically study the collapse of a fully quantum source of matter into a black hole. This is done by coupling a quantum scalar field to classical gravity and numerically simulating its evolution towards gravitational collapse. By choosing a convenient coherent state, the gravitational source is naturally split into the classical contribution plus the quantum fluctuations, which allows us to modify the well-studied classical collapse to achieve a quantum simulation. This method has many applications. For example, it can be used to investigate how quantum effects change the universality and scaling of classical critical collapse, and also to simulate Hawking radiation in real time.
We have focused on the simplest model of a massless scalar field, but this formalism is applicable to generic free quantum fields, including fermionic fields, which are inherently quantum in nature. Some initial results will be presented in this letter, which are comparable to those of a classical collapse found in, for example, [24] and from which we can already identity quantum corrections in gravitational collapse. More details of the method, along with its applications, will be presented elsewhere [25].
Classical collapse The main idea is to promote the classical scalar involved in gravitational collapse to a quantum field operator. Before that, let us briefly review the classical collapse. Our notation and method follow that of [24]. The metric in spherical symmetry is chosen to be ds 2 = −α 2 (t, r)dt 2 + A(t, r)dr 2 + r 2 B(t, r)dΩ 2 , (1) where dΩ 2 = dθ 2 + sin 2 θdφ 2 . The stress-energy tensor of the classical massless scalar field Φ(t, r) is The Einstein equations and the scalar field equation, respectively, are where Λ is the cosmological constant, taken to be zero in the classical case but needed later to regularize divergences in the quantum case 1 . The Einstein equations can be used to find a set of well-behaved first order evolution equations for the metric fields α(t, r), A(t, r), B(t, r) , or rather re-formulated variables with these fields, as done in [24]. An important ingredient in the classical collapse is the gauge choice of the lapse function α(t, r). In [24], maximal slicing is used, but this requires integrating spatially an ODE at each time step, which is numerically expensive. We instead use the "1+log" gauge [26], which results in α(t, r) becoming a dynamical variable and is nowadays more commonly used due to its simplicity and robustness.
Quantization In the quantum simulations, the scalar field is promoted to be an operator and can be expanded in terms of spherical mode functions aŝ where Y m l (θ, ϕ) are the spherical harmonics and u k,l (t, r) = r l u k,l (t, r) are the quantum mode functions, which in Minkowski spacetime are just the spherical Bessel functionsũ k,l (t, r) = k √ πω e −iωt j l (kr). The operatorâ k,l,m , is a ladder operator that satisfies . Also, the quantum state is chosen to be a particular coherent state where the first exponential is to set χ|χ = 1. Coherent states are "minimal uncertainty states" (see e.g. [27]), which are useful to relate quantum systems to classical systems, and this will allow us to make a direct comparison with classical simulations. For this coherent state we have the standard property that it is an eigenstate of the lowering operator,â In addition, we define the expectation value of the field operatorΦ in the coherent state |χ which we note is spherically symmetric and may be written as Note that φ(t, r) and the classical field Φ(t, r) satisfy the same equation of motion. This crucial fact allows us to tap into the established method for classical scalar collapse, and simply add the quantum modes on top of the classical fields. Specifically, the evolution equations of the quantum mode functions can be found using the operator equation Φ = 0 and the operator expansion (5) so that: where π k,l = (A 1 2 B/α)∂ t u k,l and ψ k,l = ∂ r u k,l . Then the expectation values of the quantum stress-energy ten-sorT µν are nicely split into the "classical" coherent part plus the terms due to quantum fluctuations. For example, the two time derivative bilinear term appearing in the expectation value of the stress-energy tensor can be written as (11) We emphasize that this convenient separation is due to our choice of the coherent state, and the scalar field is treated fully quantum mechanically.
Note that the sums appearing in the stress-energy tensor such as in Eq. (11) are divergent. These need to be regularized, which can be done by introducing a non-zero cosmological constant and five Pauli-Villars fields [28][29][30] whose masses satisfy the relations Note that in our method these Pauli-Villars fields are also simulated quantum mechanically with the full machinery of mode function evolutions. The technical details of these will be explained in a future paper [25].
Numerical setup In the simulation we use a uniform spatial grid consisting of 200 points with dr = 0.025 and dt = dr/4. For spatial derivatives, tenth order finite difference methods are used, and a tenth order implicit Runge-Kutta method [31] is used for time integration. Such high order numerical methods are needed because the evolution equations such as Eq. (10) contain terms involving 1/r, 1/r 2 and other similar terms, which are numerically unstable if not regulated properly. An alternative route to increase the accuracy would be to make use of the spectral method [25]. One also needs to use some artificial dissipation, which acts as damping for the numerical errors. In our code, this is done by adding the Kreiss-Oliger terms [32] to each field's evolution equation. Following Alcubierre's notation [24], the Kreiss-Oliger term is fourth order with = 0.5 for the quantum mode functions and = 0.1 for the metric fields and φ(t, r).
A crucial choice in the simulation is the number of quantum modesũ k,l used, which essentially sets up the quantum field operator. In practice, we can only use a finite number of quantum modes and generally the more quantum modes one uses, the better the simulated quantum operator is. We find that the maximum number of quantum modes one can use is N l = N k = 150, where N l is the number of l-modes and N k is the number of k-modes. This limit arises from the amplitude of the modes becoming too small, saturating the double precision limit of the simulations. In addition, the evolution equations of the quantum modes are increasingly unstable as N l and N k increase. A good consistency check is that when the coherent state is the vacuum, i.e. z(k) = 0, all stress-energy tensor components need to vanish. Using our methods, we have managed to simulate a system with N l = N k = 70, i.e., almost five thousand modes, for a sufficiently long time in a sufficiently large spatial region to capture the collapse to a black hole.
Throughout our simulations M P is set to unity, that is, all the dimensionful quantities are expressed in Planck units.
Results In Fig. 1, the initial stress-energy tensor components are illustrated when the coherent state is the vacuum, i.e., φ(t, r) = 0. This is done with N l = N k = 70 as mentioned above. The cosmological constant is set so that T t t − Λ equals to zero. However, the stress-energy tensor expectation value is not precisely proportional to the metric, as would be expected in the quantum vacuum, due to the truncation to a finite number of quantum modes, and hence T r r − Λ and T θ θ − Λ are not exactly zero, but under control. Also, the flatness of the stress-energy tensor components breaks down at around r = 5 before the end of the spatial grid, meaning that the evolution of the fields can only be trusted inside that region.
Initial conditions for the metric fields are derived from the Hamiltonian and momentum constraints, as in [24]. The "classical" scalar configuration is initially taken to be φ(t = 0, r) = f (r) + f (−r) andφ(t = 0, r) = 0, where f (r) = a(r/D) 2 exp −((r − R)/D) 2 , with R = 0.5 and D = 0.5 and various amplitudes a. The quantum modes initially are such that they satisfy Eq. (10) in Minkowski spacetime. For time evolution, we find that the wavepacket splits into two smaller wave-packets travelling in opposite directions, one towards the centre at r = 0 and one towards the outer boundary. With a sufficiently large initial amplitude the former packet forms a black hole in the centre.
The presence of the black hole is seen in two ways, the first being the collapse of the lapse function α(t, r), and the second being the formation of an apparent horizon. The latter is calculated using the expansion of the outgoing null geodesics H, given by where K θ θ is an extrinsic curvature component. If H vanishes for some r, then there exists an apparent horizon at that location r = r AH , with an area given by To validate the code, we have performed convergence studies for both the number of points on the grid and the number of quantum modes used, by computing the L 2norm of the Hamiltonian constraint (see e.g. [33]) and the error in the ADM mass (defined by ∆m ADM (t) = (m ADM (t = 0)−m ADM (t))/m ADM (t = 0)), respectively; see Fig. 2. The convergence with respect to the number of grid points is shown for 4 different scenarios: evolutions involving black hole/no black hole formation and 4 th and 12 th order derivative Kreiss-Oliger terms. All the simulations converge, and as expected better convergence can be achieved with 12 th order Kreiss-Oliger terms for cases without black hole formations.
We present our results in Fig. 3 for both no black hole formation (left column) and a black hole formation (right column). These are achieved by using the same parameters (mentioned above) but with different initial amplitudes of φ: a = 0.1 and a = 1.0 respectively. The first two rows show the evolution of φ(t, r) = Φ and α(t, r = 0) respectively. In the top left of Fig. 3, one can see that when the amplitude is small, the in-going wave-packet simply gets reflected at the origin. When the amplitude is sufficiently large, a black hole forms. In the top right of Fig. 3, we see that soon after separating both the in-going and out-going wave-packet freeze, indicating a black hole forming, which contains all of the in-going bigger wave-packet and most of the out-going one. In the second row, the collapse of the lapse function α(t, r) can be observed on the right, where the area of the apparent horizon, A AH , is also illustrated. The T µ ν FIG. 4. Convergence study of quantum corrections for the stress-energy tensor component ∆T t t = χ|T t t |χ −T t c t (with T t c t from the classical collapse) at t ≈ 0.94.
components at t ≈ 0.94 can be seen in the last row, and we see that the stress-energy tensor is peaked inside the apparent horizon.
To compare to classical gravitational collapse, we have also performed the corresponding simulations in the purely classical theory, with the initial data of Φ in the classical simulations matching the initial data of Φ in the quantum ones. On the scale of the plots in Fig. 3 one is not able to distinguish the classical from quantum results. Given that, in Fig. 4, we extract the quantum corrections of T t t (i.e., the difference between T t t of the quantum simulation and the corresponding classical T t c t ) at t ≈ 0.94 in the case where a black hole forms. To assess the reliability of these corrections, we perform a convergence study for various numbers of quantum modes used. We see that the quantum corrections already converge at N mode = 4900.
Summary We have developed a formalism to simulate fully quantum fields collapsing to form a black hole. In this formalism, by choosing a convenient coherent quantum state, the equation of motion for the expectation value of the field operator coincides with that of the classical field in classical numerical relativity simulations, which means that all previous numerical relativity methods may be applied to the evolution of fully quantum fields by additionally attaching the evolution of the quantum modes. For free fields, including fermionic fields and U(1) gauge fields, all the mode functions can be evolved, as explicitly shown in this paper, while for interacting quantum fields some approximations are needed for the evolution of the quantum fields, such as the Hartree approximation [34,35]. We have focused on spherical symmetry in our simulations, which is computationally much less expensive, but 1/r n -type terms appear in the equations of motion, which require extra care to handle in numerical schemes. The formalism can be applied to full 4D simulations, whose equations of motion are even free of the difficult 1/r n terms, but they require considerably more computing resources. Detailed investigations of the quantum effects of gravitational collapse such as the presence of Hawking radiation and corrections to critical collapse will be discussed in a future work [25].