Gravitational collapse of quantum fields and Choptuik scaling

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. By choosing a convenient coherent state, this formalism taps into well-established techniques used for classical collapse and adds on the evolution of the mode functions of the quantum field operator. 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 the quantum effects. We also find that Choptuik scaling in critical collapse survives in the semiclassical simulation, and furthermore the quantum deviation from the classical Choptuik scaling decreases when the system approaches the critical point.


Introduction
Black holes are fascinating objects in the universe that form as the end states of gravitational collapse, a unique manifestation of the strength and special nature of gravity. Since Hawking's famous calculation [1], we know that black holes have a temperature, as they can evaporate via Hawking quanta. This is an astonishing result, as it connects general relativity to quantum field theory, which is another pillar of modern physics that underlies all known matter sources in the universe. While stationary astronomical black holes are expected to exhibit little in the way of quantum effects outside the event horizons, the dynamical formation of black holes or other compact stars may involve extreme matter sources that are quantum in essence. Moreover, primordial black holes may form in the early universe [2] and can have a high temperature and large Hawking radiation. Purely

JHEP02(2022)183
theoretically, black holes along with their quantum properties have caused great puzzles but also great advances in our understanding of gravity and quantum field theory.
The dynamics of classical gravitational collapse with various matter sources have been studied over the years. The simple pressure-less homogenous dust model was studied as early as the late 1930s [3,4]. This has been generalized to spherically symmetric inhomogeneous dust models [5][6][7][8][9][10] and more general perfect fluid models [11][12][13][14][15][16]. While black holes are the natural end states of gravitational collapse, it is singularities that must form as the end states of gravitational collapse from reasonable matter sources, as required by the singularity theorems of Hawking and Penrose [17,18]. The cosmic censorship conjecture [19] suggests that naked singularities should be cloaked by horizons under appropriate assumptions so as to avoid loss of predictability at the singularities. However, naked singularities do appear in the above models of gravitational collapse. A major purpose of the gravitational collapse studies is to clarify the scopes and the assumptions of the cosmic censorship conjecture. See [20] for a review on gravitational collapse and spacetime singularities.
At a more fundamental level, one may model the matter source as a scalar field, which provides arguably the simplest arena to probe the nonlinearities of general relativity and its interaction with matter. Classical scalar collapse in spherical symmetry has been analytically investigated by Christodoulou [21][22][23][24][25][26], which also provides counter-examples to the naive cosmic censorship [25]. Indeed, as numerically uncovered by Choptuik [27], interesting critical phenomena exist in spherical scalar collapse that are characterized by scale invariance and universal scaling of the black hole mass with initial data. The critical collapse happens at the boundary between the parameter region of regular initial data that can form black holes and the region that can not, and the singularity at critical collapse is exactly a naked singularity, as explicitly shown in an analytical example [28]. [29,30] provided some further insights into the critical phenomena. Critical collapse has also been observed with axisymmetric gravitational waves [31], spherical collapse of radiation fluid [32] and a non-Abelian field [33]. A real analytical solution of the Choptuik spacetime has recently been constructed [34]. See [35] for a review of critical collapse.
In this paper, we study the semiclassical gravitational collapse of a massless scalar in spherical symmetry. That is, we shall keep the metric field at the classical level but fully quantise the scalar field. Such semiclassical scalar collapse is a demanding task, early attempts along this line of enquiry include [36][37][38]. The quantum properties of spherical shell collapse have also been investigated by [39][40][41], and [42][43][44][45] studied dilaton collapse in 2D gravity, which is drastically different from a 4D model for a number of reasons. First of all, by restricting to spherical symmetry or a 2D setup, these models essentially neglected all the effects of the spherical harmonics in the mode expansion, which are important for a thorough understanding of a dynamical gravitational collapse. On the other hand, we only restrict to spherical symmetry for gravity and for expectation values involving the scalar operators, and we are able to take into account as many as 100 2 spherical harmonic modes for the matter fields before our current numerical method becomes unstable. (With improved numerics, we expect to include even more spherical harmonic modes.) Because of this, the 4D setup is numerically challenging and computationally much more costly. Indeed, in order to achieve sufficient stability for our simulations, it is necessary to im-JHEP02(2022)183 plement some of the more developed numerical relativity methods for the gravitational evolution. Another complication that comes with our 4D model is the need to regularize various UV "divergences", which will be dealt by introducing several Pauli-Villars ghost fields in this paper.
We will present a formalism that allows us to simulate the gravitational collapse of a quantum scalar field into a black hole. This is the companion paper of [46] in which we reported the use of this formalism to obtain some initial results. (ref. [47], which appeared after [46], studied a similar system with a different method that only included the spherical mode for the quantum field.) In the present paper we provide a more thorough overview of the formalism with all essential technical details. Specifically, we couple a spherically symmetric massless quantum scalar field to classical gravity. Importantly, the notion of coherent states is utilised so that we can build the semiclassical system upon the previously well-studied classical collapse. Namely, the expectation value of the quantum scalar field operator in the coherent state is identified with the classical scalar field of classical collapse. This is a crucial point in the formalism because it allows us to split the stress-energy tensor components into a sum of contributions from the "classical field" and contributions from quantum modes.
It is also essential in our formalism to make use of the latest stable numerical relativity formulations of the Einstein equations, as numerically the quantum field simulation corresponds to a simulation of thousands of mode functions, which is computationally much more demanding in numerical accuracy and stability than the classical collapse. Also, to further sustain a sufficiently long period of simulation, we use as high as a tenth order finite difference method and introduce artificial dissipation terms in our equations of motion. Another major difference in simulations with products of quantum fields evaluated at the same spacetime point is that there will be UV divergences that needs cancelling, similar to loop calculations in quantum field theory. We find that the Pauli-Villars regularisation with five fields plus a cosmological constant is sufficient for our purposes.
We reproduce the results of classical scalar collapse to a black hole with our semiclassical simulations and extract the quantum deviations for the stress-energy tensor. We also simulate the Choptuik scaling for the semiclassical critical collapse and find that when approaching the critical point the quantum deviations from the classical Choptuik scaling actually decrease. We validate our results with various convergence studies.
The paper is organised as follows. In section 2 the classical collapse of a spherically symmetric scalar field is described, which will provide a basis for the quantum system. The evolution of both matter fields and metric fields are described in detail. In section 3 the aforementioned classical scalar field is quantised and expanded in terms of spherically symmetric mode functions. The chosen quantum state for the system is described, along with the expectation values of quantities relevant for calculation of the evolution equations of the geometry. We also describe how to regularise the system, since it has divergences that must be dealt with. In section 4 the numerical methods used to run the simulation and how the quantum vacuum is built up are discussed. Then, in section 5, we present the results about the evolution of the quantum collapse into a black hole, highlighting the quantum corrections to the classical collapse, and also the results on the quantum Choptuik JHEP02(2022)183 scaling near critical collapse. Lastly, we provide the convergence studies of our simulations in section 6, and conclude in section 7.

Classical collapse
In this section the theory for the simulation of the collapse of a spherically symmetric classical massless scalar field is discussed using the methods of [48] but with a different gauge choice. First, the geometrical side of the Einstein equation will be presented using variables that make the evolution equations better suited for lattice simulations, and then we show how to include the matter fields. Before connecting the geometry to the matter equations, the gauge choices are reviewed briefly. Lastly, the choices for the initial conditions for the simulation are shown. Note that in the simulation all the fundamental constants are set to one, = c = M P = 1. If displayed, their purpose is simply illustration.

Geometry
We will restrict our simulations to have spherical symmetry in 3+1D in this paper. For classical simulations, it is certainly viable to do full 3 + 1D simulations with modern hardware, but real time simulations with quantum fields involve solving partial differential equations for many quantum modes simultaneously, which is computationally very expensive. For our purposes, the following spherically symmetric line element is chosen: where α(t, r) is the lapse function and A(t, r) and B(t, r) are functions governing the spatial part of the metric. The lapse function α(t, r) has a crucial role to play in the collapse -not only will it be used to choose a Bona-Masso type gauge condition, but also one can use it to identify the formation of a black hole, due to the "collapse of the lapse" phenomenon [48]. However, we shall not directly use α(t, r), A(t, r) and B(t, r) as our dynamical variables, with which the Einstein equations are second order differential equations. To cast the Einstein equations as first order equations, we can define the following variables where the dot represents differentiation with respect to t, and the prime represents differentiation with respect to r. Note that λ is defined in order to guarantee the local flatness of the metric (at r = 0). Local flatness in these coordinates manifests itself as A(t, r = 0) = B(t, r = 0) for all t. Many 1/r terms in the evolution equations discussed below are only regular (not divergent) numerically, if this local flatness condition is exactly met. Small numerical errors in A(t, r = 0) and B(t, r = 0) could mean exponential divergences in the numerical JHEP02(2022)183 system. Hence, using the definition of λ(t, r), and further requiring it to be antisymmetric, the local flatness condition is exactly met at each time step.
In addition, to ensure that the system of equations is strongly hyperbolic for any gauge choice of α(t, r), we must perform some further changes of variables. Namely, instead of (K A , D A ), we will use (K,Ũ ) respectively, The detailed motivation for these choices are explained more in depth in [48], but essentially for some gauge choices for the lapse function α(t, r) the original system is not strongly hyperbolic. Hence, using the above change of variables, our system is strongly hyperbolic for any gauge choice of α(t, r), which are discussed at the end of this subsection. Then the non-zero components of the Einstein tensor can be written as The G r t and G t t components of the Einstein equations are actually constraint equations, which do not contain time derivatives with the new dynamical variables above. They are called the momentum and Hamiltonian constraint respectively where G t r and G t t are understood to be replaced with the corresponding components of the stress-energy tensor via the Einstein equations. These constraints provide a sanity check on whether numerical errors are under control in the simulation.

JHEP02(2022)183
Our metric ansatz, eq. (2.1), contains gauge/coordinate degrees of freedom. Particularly, we may freely choose α as a gauge choice, which slices the spacetime into hypersurfaces with the same times. Since we are dealing with gravitational collapse, we need to choose a singularity avoiding gauge choice. This ensures that even when a singularity is present inside an apparent horizon, the spatial grid outside of the horizon may still be evolved. The two established families of singularity avoiding gauge choices are maximal slicing conditions and Bona-Masso type slicing conditions [49].
The maximal slicing conditions are used in [48], which is equivalent to demanding that the extrinsic curvature of the equal time hypersurface is constant in time. This leads to a constraint equation involving α that needs to be integrated with respect to r in spherical symmetry.
The Bona-Masso type slicing conditions [49], on the other hand, treat α as a dynamical variable, which is artificially specified by an evolution equation along with the evolution of other physical variables. These are also called hyperbolic slicing conditions, since α(t, r) obeys a wave equation with a source containing K.
Although maximal slicing is a robust gauge choice, it is computationally expensive and slow, since one needs to integrate spatially at every time step. The Bona-Masso type gauges possess similar robustness, but requires less computational power, since one just has an additional dynamical field. Therefore, in our simulations, we use the Bona-Masso slicing conditions. Specifically, the lapse function α(t, r) is chosen to satisfy the following evolution equationα where f (α) is a suitable function of α, to be tuned in a particular problem. The α evolution also leads to the evolution of D α , which is given bẏ The most used gauge choice within the Bona-Masso family is the 1+log gauge, which is given by This choice ensures strong singularity avoidance and is the one we use.

Matter
Let us now introduce the scalar field involved in the simulation. We shall for simplicity work with a free massless scalar Φ, which couples minimally to gravity. Classically, the scalar field acts as a stress-energy source to the Einstein equation, and gravity also influences the matter via the Klein-Gordon equation. We have the following stress-energy tensor for the scalar We shall introduce the following new field variables for the scalar field 14)

JHEP02(2022)183
which are needed to formulate the scalar evolution equations with only first order derivatives. With these variables, the stress-energy tensor can be decomposed into Here n µ is the unit timelike vector orthogonal to the equal time hypersurface, and γ ij is the inverse of the spatial part of the metric. Note that in our coordinates these quantities are simply ρ = −T t t , j A = T t r , S A = T r r and S B = T θ θ .

Evolution equations
The Einstein equations and the scalar equation of motion are respectively given by where M P = 1/ √ 8πG is the reduced Planck mass and we have also added the (bare) cosmological constant Λ as a stress-energy source. As we are interested in an asymptotically flat spacetime, for classical simulations, we can set the cosmological constant to zero, but the presence of a cosmological constant is required to regularise the quantum field in the semiclassical system, as will be explained in section 3.4 in detail.
As mentioned above, we want to cast these field equations in a first order form, which can be achieved by evolving the following variables introduced previously where the first three fields are matter fields, the others are metric fields, and the (temporarily added) superscript indicates whether the field is symmetric or antisymmetric around r = 0, as we shall be using initial data with definite parity to help with numerical stability. The Klein-Gordon equation (2.20) can be cast in the following forṁ (2.24)

JHEP02(2022)183
As for the metric fields, the relevant components of the Einstein equations are the following: where the first two equations are constraint equations, i.e., eq. (2.8) and eq. (2.9), and the latter two are evolution equations of our chosen dynamical fields. For A(t, r), B(t, r) and D B (t, r), their evolution equations can be found straightforwardly by simply using their definitions:Ȧ The evolution of the gauge degrees of freedom α(t, r) and D α (t, r) are to be chosen by hand, that is, by imposing the following slicing conditionṡ We will work in 1+log gauge throughout, for which we choose f (α) = 2/α. To find the evolution equations of K(t, r) and K B (t, r), the Einstein equations have to be used. The equation of motion for K B (t, r) is derived by combining the G r r and G t t components, which givesK For the evolution of K, the equations for G r r , G t t and G θ θ can be used to finḋ For λ, using its definition, we can geṫ

JHEP02(2022)183
which is well defined in the continuum limit. However, the 1/r factor on the right hand side makes this evolution equation numerically highly unstable when discretising the equation and putting it on a lattice. To overcome this problem, the momentum constraint (2.8) can be used to remove this factor and we arrive aṫ (2.37) Note that other 1/r terms in the evolution equations (such as in eq. (2.24)) are automatically regular by demanding the dynamical variables to be symmetric or antisymmetric. However here, in eq. (2.36), the combination of (K − 3K B )/r must be regular. Analytically, this is true since the local flatness condition A(t, . Numerically, on the other hand, small errors can spoil this cancellation hence making the term unstable. Thus the need to make it regular by using the momentum constraint. Similarly, using its definition and the Hamiltonian constraint eq. (2.9), the evolution forŨ can be found to bė This completes the derivation of the evolution equations for our dynamical variables, which, unlike the original Einstein equations or even the ADM decomposed version, are generally well behaved numerically after discretisation.

Initial conditions
To evolve the system, some initial conditions must be chosen. Some care must be taken for a system with gauge degrees of freedom, in which case some of the equations of motion are actually constraint equations. We follow a free evolution scheme in which the constraints are solved once initially when preparing the initial data. Then the keeping of the constraint equations can be used to monitor how well the simulation works. For the matter fields, we choose the initial conditions to be where here the superscript 0 denotes the initial condition of the corresponding quantity at t = 0. A natural choice for our purposes is to take the function f (r) to be a Gaussian wave packet. In our simulations, we use the following family of functions where a is the amplitude and D the width of the wave packet. Note that the above definition of f (r) ensures that the initial scalar field is an even function, and thus all initial JHEP02(2022)183 1/r terms in the equations of motion are regular. The initial conditions for the gravity fields are chosen to be such that they solve the constraint equations. Note that the choice Π 0 = 0 is convenient, since then we have j 0 A = 0, which means that the momentum constraint in eq. (2.8) is automatically solved by the choices K 0 = K 0 B = 0. The initial conditions for the other gravity fields require some calculations. To do this, we first take A 0 (r = 0) = 1, and integrate out the following equation to find the A 0 (r) function, which comes from the Hamiltonian constraint along with the initial condition ρ 0 = (Ψ 0 ) 2 /(2A 0 ). Then to find λ 0 andŨ 0 , we use their definitions along with A 0 : Thus we have initial conditions for all classical dynamical variable fields. This completes the setup for the classical gravitational collapse in spherical symmetry, which provides the base framework for the quantum/semi-classical simulations. We will run the classical simulations along with the semi-classical ones for a comparison later. Before that, in the following, we shall first set up our formalism to simulate quantum fields on curved space.

Quantisation
We want to develop a semiclassical simulation scheme where the metric fields are still treated classically but we can consistently treat the matter fields as quantum fields. This is a valid approach in situations where the quantum effects of the matter fields have become important but the quantum gravity effects have yet to kick in. For our case, it is only the classical scalar field that needs to be promoted to be a quantum scalar field. To this end, the c-number Φ needs to be upgraded to an operatorΦ, and its equation of motion is treated as an operator equation. We can not directly simulate the operator equation, but once we expandΦ in terms of creation and annihilation operators, theΦ operator equation can be decomposed into the evolution equations for the functions in front of the creation and annihilation operators, the mode functions, which are c-numbers and can be put in a lattice. Of course, the difficulty is that there are infinitely many of these mode functions. Since the classical simulation is spherically symmetric, it is also convenient to further expand the mode functions in terms of the spherical harmonics Y (θ, ϕ). The feedback of the quantum scalar field to the classical geometry is obtained via the semiclassical Einstein equations,

JHEP02(2022)183
by taking the quantum expectation value of the stress-energy tensor of the scalar field. Explicitly, we will work with the following semiclassical system of equations of motion whereT µν = ∂ µΦ ∂ νΦ − 1 2 g µν g ρσ ∂ ρΦ ∂ σΦ and a coherent state will be chosen for the quantum state |χ .
The key element of our formalism is that, due to a convenient choice of the quantum coherent state, the sources for the Einstein equation can be nicely separated into a classical background plus some quantum fluctuations. This fact allows us to utilise all the existing classical numerical relativity setup, reviewed in the previous section, and simply add the quantum mode functions of the matter fields on top of the simulation of the classical fields. Crucially, though, the scalar field remains fully quantum mechanical.

Quantum scalar field
As mentioned, the real scalar field is now promoted to a Hermitian quantum field operator Φ, which can be expanded in spherical mode functions aŝ where Y m l (θ, ϕ) are the spherical harmonics and in the following we shall mostly use the mode functions with a factor of r l stripped off For notational simplicity, we have suppressed the k, l indices for u andũ, and we will also suppress these indices for the π and ψ quantities defined shortly, when not confusing. The u(t, r) = u k,l (t, r)'s are the unknown quantities to be solved in our simulations. The operatorsâ k,l,m andâ † k,l,m are the corresponding ladder operators for each quantum mode denoted by the indices k, l, m. The standard route to quantisation is then to impose the canonical commutation relation between the scalar field and its conjugate momentum If we require the mode functions satisfy this is equivalent to requiring the ladder operators to obey

JHEP02(2022)183
We want to emphasize that although the classical geometry is restricted to spherical symmetry, the quantum operators are allowed to fluctuate non-spherical symmetrically. The spherical symmetry is reflected in the quantum state we will choose below. Similar to the classical scalar field, we also define the following variables for the quantum mode functions in order to formulate the evolution equations in a first order form We then note that the evolution equations for the mode functions u and ψ are jusṫ For the π evolution equation, one can use eq. (3.3) along with eq. (2.24) to find thaṫ which actually carries the real dynamics of the quantum scalar field. The last term in eq. (3.11) is only strictly regular if A(t, 0) = B(t, 0). Similar regularity conditions were encountered in the case of the classical collapse as well. There, the use of λ ensured regularity and hence one can use it here as well. Thus, the final dynamical equation to be used for the scalar mode functions is given bẏ We emphasize that here we have suppressed the k, l mode indices for u, ψ and π, that is, the above equation represents an infinite set of c-number equations, which however are decoupled as we work with a free scalar field. To get a good approximation for the quantum scalar field, one should evolve as many of these mode functions as possible, which makes the semiclassical simulations much more expensive computationally than the classical ones. In order to find suitable initial conditions for the quantum mode functions, let us recall that the equation of motion in Minkowski space for with j l (kr)'s being the spherical Bessel functions. This means that the solution for u in Minkowski space is given by (3.14)

JHEP02(2022)183
For simplicity, we shall choose this as our basis of functions at t = 0. This entails the following initial conditions at t = 0

Coherent state
After setting up the evolution of the field operators, one also needs to specify the quantum state of a system. We choose the quantum state to be a spherically symmetric coherent state which is an eigenstate of the lowering operator,â k,0,0 |χ = z(k) |χ , andâ k,l,m |0 = 0. The first exponential in the |χ definition above is fixed by normalisation χ|χ = 1. Coherent states are known to be useful to connect quantum systems to classical ones, and they are "minimal uncertainty states" (see e.g. [50]), which are in a sense "the closest to the classical states". This specific choice of the quantum state will be crucial in this formalism, and it will allow for a direct comparison between the semiclassical collapse and classical collapse.
It is worth pointing out that the state we are considering is an eigenstate of the initial lowering operator, and that in the Heisenberg picture, where operators evolve, the state need not be an eigenstate of the time-evolved lowering operator. This, however, is not a factor in the simulations we perform. Let us define the expectation value of the quantum scalar field operator in the chosen coherent state to be φ(t, r) ≡ χ|Φ(t, r) |χ . (3.17) which is often called the "classical field". Indeed, we may identify it with the scalar field in the classical gravitational collapse. To see this, note that, making use of the coherent state definitionâ k,l,m |χ = δ l,0 δ m,0 z(k) |χ , we can get where h.c. stands for the Hermitian conjugate of the previous terms and the fraction is due to the fact that Y 0 0 = 1/2 √ π. By the evolution equation of the mode functionũ k,0 (t, r), we see that this expectation value φ(t, r) satisfies the same equation of motion as the classical scalar field in section 2. Hence the coherent state choice allows us to separate the classical and quantum contributions in the stress-energy tensor, which is crucial to tap into the established methods to simulate classical gravitational collapse. This is discussed in the next subsection.

Quantum stress-energy tensor
By virtue of the chosen coherent state above, the objects appearing in the stress-energy tensor automatically separate into a sum of contributions from the coherent state expectation value φ(t, r), which only involves the spherically symmetric mode functionũ k,0 (t, r), and contributions from the mode functionsũ k,l (t, r). Specifically, using eq. (3.3), the choice of the quantum coherent state eq. (3.16) and some identities for the spherical harmonics, it is straightforward to find that the expectation values for the bilinears appearing in the relevant stress-energy tensor components are where we have explicitly put back the reduced Planck constant and the speed of light c. Note that here φ(t, r) is taken to be the same as the scalar field in the classical gravitational collapse. The fact that the quantum expectation values of these operators separate into the classical coherent part and the extra quantum part is merely due to the choice of the coherent state. Note that eqs. (3.22)-(3.23) contain only quantum fluctuations, as expected, since the classical pieces of those terms do not depend on (θ, ψ). In addition, by setting to zero, one can readily switch off quantum contributions. Therefore, in our semiclassical simulations, we shall take the decomposed stress-energy tensor quantities ρ, j A , S A and S B to be where we have defined Ô ≡ χ|Ô |χ and used the following relations for the dynami-

Regularisation
The semiclassical eq. (3.1) needs regularisation, since the quantum expectation values on the right hand side diverge, which manifests itself as the infinite sums over the mode functions appearing in eqs. (3.19)-(3.23). In our lattice simulations, these divergences appear as large, unphysical numbers in the evolution equations that need to be regularised, as is to be expected for a quantum field theory. In a lattice simulation the popular regularisation scheme used in analytic perturbative computations, dimensional regularisation, can not be used. In general, dimensional regularisation is useful because it maintains all the symmetries of a theory, including Lorentz symmetry. So if dimensional regularisation were to be used, the introduced bare cosmological constant Λ would be sufficient to regularise the UV divergences in eq. (3.1). However, the finite number of modes (or lattice points) used in a numerical scheme necessarily breaks Lorentz invariance, which means that the UV divergences can not be cancelled by a simple cosmological constant counterterm. (In other words, Lorentz invariance is not a symmetry of a lattice field theory, which has a built-in hard cutoff, and in a regularization scheme with a hard cutoff, the UV divergences in eq. (3.1) are not proportional to η µν .) Instead, we shall introduce Pauli-Villars fields [51] to cancel the various divergences that appear. As we shall see, we need to introduce several Pauli-Villars fields to achieve our goal. This is because while one Pauli-Villars ghost field would eliminate the original divergences in eq. (3.1), it also introduces some less severe divergences, which need to be balanced/canceled by introducing extra Pauli-Villars fields. The Pauli-Villars regularisation scheme consists of adding some extra auxiliary (ghost) fields to the physical one, which by construction cancel the UV divergences. This happens by giving the ghost terms in the Lagrangian the appropriate signs, so that the divergences disappear in the stress-energy tensor. These ghost fields must have larger masses than the physical field, in order to only influence the UV physics. The ghost fields are chosen also to cancel each others' contributions to the various divergences.
The Pauli-Villars regularisation scheme consists of adding some extra auxiliary (ghost) fields to the physical one, which by construction cancel the UV divergences. This happens by giving the ghost terms in the Lagrangian the appropriate signs, so that the divergences disappear in the stress-energy tensor. These ghost fields must have larger masses than the physical field, in order to only influence the UV physics. The ghost fields are chosen also to cancel each others' contributions to the various divergences.
Note that even though our quantum field will live on a curved spacetime, it is sufficient to fix the masses of the ghosts by first regularising in Minkowski spacetime. This is because the divergences are UV in nature, hence short distance. In short distances, due to local flatness, the Minkowski approximation still holds. Thus, once the UV divergences are fixed in Minkowski, they are cured on curved spacetimes as well.

JHEP02(2022)183
To illustrate the divergences involved in the system, let us take a massless scalar field in Minkowski spacetime. The Lagrangian is just This scalar field may be expanded in terms of mode functions similarly as in eq. (3.3). Then, the energy density and pressure for the field can be found to be [52] where k is just the wave number and a cut-off M has been introduced. These integrals are straightforwardly evaluated to be (3.32) One can see that these terms are divergent, and need to be regularised. The introduction of a bare cosmological constant Λ, which is equivalent to the normal ordering for the quantum fields, can only cancel divergences of the form vac = −P vac , and so is not able to cure the divergence in eqs. (3.31)-(3.32), which does not satisfy the equation of state of vacuum energy. Thus, one needs to introduce Pauli-Villars ghost fields to cancel the divergent density and pressure contributions in eqs. (3.31)-(3.32). Let us see how introducing one ghost field G 1 with mass m 1 (which is non-zero) changes these expressions. The ghost field is introduced in the Lagrangian with opposite sign kinetic and mass terms, so that Then, the energy density and pressure for the system is Hence the large M expansion of these integrals is modified to be

JHEP02(2022)183
The original divergence coming from the physical field is cancelled, however, due to the nonzero mass of the ghost field, there are new, quadratic, divergences introduced. In addition, now there are finite terms in the expansions from the ghost field, which are unphysical. Therefore, both the new divergences and the finite ghost contributions need to be cancelled. Notice that the third terms in eqs. (3.35)-(3.36) can be cancelled by the cosmological constant Λ, so that However, we are still left with the quadratic divergences (the second terms) and the unphysical finite ghost contributions (the fourth terms). To cancel these, we need to introduce more Pauli-Villars fields. On the other hand, to make sure that the first terms in eqs. (3.37)-(3.38) keeps being cancelled, one needs a total number of fields that is even and therefore two Pauli-Villars fields will not suffice. Note that the third terms will be cancelled by the cosmological constant no matter the number of Pauli-Villars fields. Hence, one needs to introduce at least three Pauli-Villars fields G 1 , G 2 and G 3 , with masses m 1 , m 2 and m 3 . To cancel the M 4 terms, one simply uses the opposite sign for the kinetic and mass terms in the Lagrangian for two of them, e.g. (3.39) The solutions to this set of equations always involve m 1 or m 3 being zero. However, the Pauli-Villars masses must be larger than the physical mass, and therefore these are not suitable solutions. In order to obtain a set of masses that are nonzero, two more Pauli-Villars fields have to be introduced, G 4 and G 5 , with G 5 having the opposite sign kinetic and mass terms. Thus, the masses must obey (3.40) There are many solutions to this set of equations, and we choose the following: Therefore, using five Pauli-Villars fields with the above particular masses, the system is regularised. Note that, together with the physical field, this means one has to evolve six

JHEP02(2022)183
dynamical quantum fields now. The Lagrangian for the full regularised semiclassical system is then: (3.42) where now three ghost fields (with i = 1, 3, 5) have opposite sign kinetic and mass terms. Hence, the system with the Lagrangian in eq. (3.42) produces regularised stress-energy tensor components and so the expressions in eqs. (3.24)-(3.27) will also be finite.

Simulation setup
We shall solve the coupled, highly nonlinear differential equations of the quantum field plus gravity system numerically by discretising the spacetime and putting the equations on a lattice. In this section, the numerical scheme of the spacetime grid is presented along with techniques to solve the discretised system. A subsection is also devoted to the initialisation of the quantum field, which involves choosing a specific number of quantum modes that build up the full quantum operator. It is illustrated how the stress-energy tensor components change for different numbers of modes.

Numerical methods
In the simulation we use a uniform spatial grid consisting of 500 points with dr = 0.025 and dt = dr/4. For spatial derivatives, tenth order finite difference methods are used, and similarly a tenth order implicit Runge-Kutta method [53] is used for time integration. Such high order numerical methods are necessary due to the evolution equations, such as eq. (3.11), containing terms involving the 1/r and/or 1/r 2 factor, which are numerically unstable due to the coordinate singularity at r = 0 if not treated carefully.
To this end, one also needs to use some artificial dissipation, which acts as damping for the numerical errors. In our code, this is done by adding Kreiss-Oliger terms [54] to the evolution equation of each field. For example, for φ(t, r), at time step n with the original evolution step denoted schematically by G(φ n ), φ at the next time step would be calculated by Then, adding the extra dissipation term, the evolution equation G(φ n ) is modified as where is a positive constant which is smaller than one, and 2N is the order of dissipation, and is an integer. This added term essentially damps the modes with wavelength close to the grid spacing dr.
In order to preserve the tenth order accuracy of the system, the dissipation term must be at least twelfth order. Even though the accuracy of the system is maintained to be high, in our experience the simulation is more robust when lower order dissipation is used. Thus, following Alcubierre's notation [48], the Kreiss-Oliger term is fourth order with = 0.5 for JHEP02(2022)183 the quantum mode functions and = 0.1 for the metric fields and φ(t, r). Thus, we use the tenth order methods to stabilise the evolution around r = 0, but our overall accuracy is fourth order, which is sufficient for our purposes.
The semiclassical simulation starts with the initial conditions described in subsection 2.4 with the coherent state expectation value of the quantum field replacing the classical scalar field. In addition, the initial quantum mode functions are described in eqs. (3.15). We also run the classical simulation with the same initial conditions for a comparison. We solve the constraint equations only once at t = 0 and then use the Hamiltonian constraint to monitor the accuracy of the our simulations.
Note that the accuracy and stability of the system are dependent on the initial coherent state expectation value of the quantum field. This is parameterised by the Gaussian configuration parameters a, D. In terms of stability, the most sensitive part of the simulation is the quantum mode functions, specifically the mode functions with high k and l values. These dynamical variables have the highest tendency to develop instabilities at the origin. Even though both the spatial derivatives and time iteration method are tenth order and dissipation is added artificially to the system, instabilities still arise for high k and l modes, which limits the number of k and l modes that can be added to achieve better accuracy. The artificial dissipation, on the other hand, greatly delays the appearance of the instabilities.
One also needs to choose numerical parameters for the mode functions and the integrals, appearing in the stress-energy tensor components, that become discrete sums. Firstly, the mass m 1 of the ghost field G 1 must be chosen, which is simply taken to be m 1 = 1.0 in Planck units through our simulations. Crucially, this needs to be larger than the energy of the initial scalar field, which is satisfied by the above mentioned initial conditions. In addition, the continuous integrals of the mode functions over k, in e.g. eqs. (3.19), must be numerically approximated. To this end, a minimum wave number k min is chosen, which will then be the minimal step between wave numbers dk = k min , and the integrals over k then become summations multiplied by dk. In our simulations we choose k min = π/15 throughout. Lastly, one needs to choose a finite number of mode functions to involve in the simulation, namely the number of k-modes and l-modes, which will be denoted by N l and N k respectively. This will be discussed in the next subsection.

Vacuum
As mentioned, a crucial choice in the system is the number of quantum modesũ k,l used, which essentially defines 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 simulated the quantum operator is. We find that the maximum number of quantum modes one can use is N l = N k = 150 in our current numerical setup. This limit arises from the amplitude of the modes becoming too small, saturating the double precision limit of the simulations.
Due to the finite amount of mode functions, the quantum field is well-defined only in a limited spatial region around r = 0. In addition, thanks to instabilities from the high k and l valued mode functions, one can simulate the well-defined quantum field only for JHEP02(2022)183 a finite amount of time also. Hence, there are both spatial and temporal limitations of the simulation. More specifically, for a particular set of initial conditions for the coherent state expectation value, there is a characteristic spatial region where the quantum field is accurately modelled, as well as a characteristic time frame during which all the relevant evolution takes place and through which all evolved mode functions must stay stable. Thus the initial coherent state expectation value must be chosen so that it vanishes outside the region of well-defined quantum field operator, as well as the evolution taking place before the mode functions become unstable. This ensures that everything stays physical in the chosen spatial region and time frame.
Let us discuss the case when the coherent state expectation value is just the vacuum, i.e., z(k) = 0. Since the coherent state expectation value does not contain divergences, but only the quantum mode contributions do, if the quantum mode contributions are well-defined, then the full quantum operator is as well. The field operator in the vacuum contains exactly these quantum mode contributions only. Since all stress-energy tensor components need to vanish in a well-defined quantum vacuum, this is a good consistency check. When this is true, we expect that the full quantum field operator is well-defined.
As mentioned before, using a larger number of mode functions, one achieves a welldefined quantum vacuum in a larger spatial region. To illustrate this, the stress-energy tensor components, for the case of zero coherent state expectation value, are presented in figure 1 for various numbers of mode functions with N k = N l . The stress-energy tensor components have been regularised in the previously presented manner. The cosmological JHEP02(2022)183 constant has been chosen such that the T t t component is zero at r = 0. This is merely an ad-hoc choice, which then introduces a small deviation of the T r r and T θ θ components from zero, as seen in figure 1. The relative difference between the components is unchanged by the choice of the cosmological constant.
From figure 1, it is evident that more mode functions make a better quantum vacuum. The spatial region in which all stress tensor components (approximately) vanish, is more or less constant, namely when r < 5 for the chosen set of initial conditions. The deviation from zero of the T r r and T θ θ components, on the other hand, greatly decreases as the number of mode functions is increased. This deviation introduces a small systematic error in our simulation, but generally the errors coming from other sources of the full evolution are actually greater than this error.
Another set of plots are shown in figure 2, where the number of l-modes and k-modes are not equal (N k = N l ), in order to demonstrate how varying one and holding the other number constant influences the vacuum. Note that in these plots, the spatial region is increased in preferable cases, compared to figure 1, now r stretching out to r max = 20. In the first row, N k is held constant, whilst N l is increased. The region where the vacuum is welldefined becomes larger as this happens. The relative difference between the components stays constant throughout. In the second row, we hold N l constant but increase N k . We see now that the relative difference between the components decreases, as one adds more k mode functions. Unfortunately, the flat region where the quantum field is well-defined is actually shrinking as well. This is the reason why, in figure 1, even though both the number of l-modes and k-modes are increased, the flat region stays relatively constant.

JHEP02(2022)183
For our purposes, we have managed to simulate a quantum field with N l = N k = 50, i.e., 2500 modes, for a sufficiently long time in a sufficiently large spatial region to capture the gravitational collapse to a black hole as well as Choptuik scaling.
Having established the physical region of the quantum field operator, a "good" set of initial conditions for the coherent state expectation value is presented in the next section which vanishes outside the aforementioned spatial region, and its relevant evolution finishes before the mode functions spoil the simulation.

Gravitational collapse and Choptuik scaling
Once the initial conditions are determined for both coherent state expectation value φ(t, r) and quantum modesũ k,l , we are ready to simulate the gravitational dynamics with a quantum field. In this section, results are shown for initial conditions ending up collapsing into a black hole and also dispersing to infinity. In addition, the system is studied close to the critical amplitude in order to demonstrate Choptuik scaling in critical collapse. We also compare with the classical collapse and identify the quantum deviations in these phenomena.

Black hole formation
After establishing that the quantum field is well-defined in a sufficient region around r = 0, one can add a nonzero initial coherent state expectation value, and evolve it. The evolution can have two distinct final spacetime structures: Minkowski or Schwarzschild, depending on the initial conditions. We will show the evolution of a chosen initial coherent state expectation value for both cases of black hole formation and no black hole formation.
As discussed, choosing the coherent state expectation value is equivalent to choosing the initial classical scalar configuration. We focus on an initial classical scalar configuration described by eq. (2.39) and eq. (2.40) with D = 1.0. This initial amplitude can be tuned to achieve an evolution involving a black hole formation or otherwise. For the former case, we can choose a = 5.0 and for the latter a = 1.0 is a possible choice. Note that these results were obtained using N l = N k = 50, in other words, 2500 quantum modes. The classical simulations without adding the quantum modes are also performed for a comparison. As expected, for our choices of initial conditions, the classical simulations are a good approximation of the semiclassical ones. More specifically, the plots presented in figure 3 would be almost indistinguishable from the ones obtained using classical matter. Crucially, though, they are not exactly the same, we will come back to this at the end of this section.
The results for the two separate simulations are presented in figure 3. The left column corresponds to the case of no black hole formation (subcritical), and the right column to when the quantum field collapses and a black hole forms (supercritical). (The critical case is presented in the next subsection.) The different rows show various quantities relevant to the evolutions. Let us go through these row by row.
The first row shows plots of the coherent state expectation value of the quantum field, namely φ(t, r) = Φ , at three times: at the very beginning, in the middle of the  simulation, and at the end of it. In the case of no black hole formation, one can observe that the initial Gaussian wave-packet immediately starts to drop, after which it propagates outward. In the case of the black hole formation (right column), the wave-packet starts a similar evolution, however, it quickly freezes around r = 0 by virtue of the singularity avoiding gauge choice for the lapse function α(r, t). One can observe a small fraction of the wave-packet still escaping, albeit extremely slowly.
The second row shows the evolution of the lapse function at the very centre, α(0, t). This is significant, since the lapse function vanishes if a black hole forms, due to our gauge JHEP02(2022)183 choice [48]. Thus α(r, t) essentially signals the strength of the curvature at any given point.
It is plotted at the centre of the spherically symmetric grid, since that is where α(r, t) first reaches zero in case of gravitational collapse. On the left hand side α(0, t) is plotted when no black hole forms. One can see that at its minimal value it reaches about 0.5-0.6, signalling quite a significant curvature. However, it starts increasing and asymptotically approaches one, in other words, Minkowski spacetime is recovering. The plot on the right hand side, corresponding to the supercritical case, shows the evolution of the central value of the lapse function as well as the area of the appearing apparent horizon. The apparent horizon is found using the expansion of the outgoing null geodesics Θ, given by where the physical significance of K B is that it is the (θ, θ) component of the extrinsic curvature. If Θ vanishes for any given r, then there is an apparent horizon located there with radius r = r AH and area A AH = 4πr 2 AH B| r AH . Thus, in figure 5.1, the central value of α(r, t) is plotted along with A AH /4π. It can be seen that α(0, t) collapses to zero rather quickly, during which an apparent horizon appears with area A AH /4π ≈ 0.8. The latter then asymptotes to a constant value of A AH /4π ≈ 1.25. The presence of a black hole is unequivocally signalled by these two phenomena; the collapse of the lapse function, and the appearance of an apparent horizon.
The last row in figure 3 shows the expectation value of the stress-energy tensor components in the coherent state at the end of the simulation, at t = 3.0. When no black hole forms, as expected, these values are vanishing at r = 0 and follow the wave-packets propagating outward. On the other hand, in the supercritical case (right figure), all the components of the stress-energy tensor are centered around r = 0, which is again expected, since most of the energy is located in the black hole around the center. Additionally, note that the stress-energy tensor components are larger by three orders of magnitude in the black hole case compared to the no black hole case.
Therefore we have seen that the simulation goes as expected for both subcritical and supercritical cases, with clear signs of black hole formation in the supercritical case. These plots, however, are very similar when the matter is just classical as well. The obvious step then is to compare some values in the semiclassical simulation to the classical one. As a proxy to measure the quantum effects, we shall plot the differences between the quantum expectation values of the stress-energy components and the classical stress-energy components: where T µ (c) ν corresponds to the purely classical simulation, without any quantum modes added. This is plotted in figure 4. Interestingly, as it can be seen in the figure, the largest quantum effects are cloaked by the horizon of the black hole. Nevertheless, there are some quantum deviations spilling out of the horizon. Also note that even though the quantum field is said to be well defined only until r = 5 (visible in figure 1), the quantum mode contributions look flat until around r = 10. This is just because the quantum mode contributions diverge after the well-defined region (0 < r < 5) rather slowly. JHEP02(2022)183 Figure 5. Black hole mass M BH against initial amplitude a for both classical and semiclassical cases. In the left graph, both the M BH and a axis are linear, whereas in the right graph they are plotted (natural) logarithmically to illustrate the proportionality. In addition, the bottom part of the left graph shows the deviation of the semiclassical black hole masses from the classical black hole masses, in percentages. As the system is tuned closer to the critical point, the deviation decreases.
Due to the gradients of functions at the centre reaching arbitrarily large values as the critical amplitude is approached, since we are closer to the singularity, the simulation becomes more unstable and the data points are limited in the vicinity of the critical point. Nevertheless, the value of the critical amplitude and the critical exponent can be still extracted using curve fitting. This has been done with the classical points, using a fitting function of the form: F (a, γ, a crit , c) = c(a − a crit ) γ , where a is the data array of the amplitudes and the other three entries, the critical exponent γ, the critical amplitude a crit and the constant of proportionality c, are parameters of the fit. This function is the expected critical behaviour of the mass of the formed black holes, classically well-known. Using the non-linear least squares fit, the best values of the three parameters are found. Note that in the fitting, not all classical data points in figure 5 have been used, only the ones with smaller amplitudes, as those are the ones expected to obey the critical behaviour. Indeed, one can see the significant departure from the critical behaviour in case of the last few points (with largest amplitudes). The specific values of the fit plotted in figure 5 is the following: γ = 0.379 ± 0.006, a crit = 1.87 ± 0.01, c = 0.442 ± 0.003.
The accepted value of the classical scaling exponent in the literature [30,55] is γ = 0.374 ± 0.001, which is in close agreement with the value found here. As mentioned above, the difference between the semiclassical and classical black hole masses decreases for smaller black holes as the system flows to the universal scaling limit.

JHEP02(2022)183
On the other hand, the curvature at the centre at the critical fixed point is theoretically infinite, so one would naively expect that quantum effects also get larger close to the critical point. Additionally, the strength of Hawking radiation is believed to scale as M −2 BH , which also suggests this. Whether the semiclassical gravitational collapse approaches the classical one at the critical point is beyond the scope of this paper.

Convergence analysis
Code validation is important to ensure that the approximate numerical system does indeed converge toward a true exact solution, ruling out the presence of systematic numerical errors that could be falsely identified with physical effects. In this section, the convergence properties of our semiclassical system are analysed in order to validate the simulation numerically. Our simulations have two significant parameters: the number of points used on the physical grid (N grid ) and the number of mode functions (N l , N k ) used in the semiclassical simulation. These will be analysed separately, since they each individually must be convergent to a true solution for the convergence of the whole semiclassical system.
Let us first discuss the number of grid points, N grid . Specifically, by this, we mean the number of discrete spatial points and temporal steps that our 2D computational grid is divided into. Thus N grid is inversely proportional with the separation between grid points spatially dr and also temporally (through the relation dt = dr/4). Hence, as one uses more grid points for a given physical region, the separation between the points decreases with the grid approaching the smooth continuous spacetime. To check convergence for this parameter, the L 2 -norm of the Hamiltonian constraint is analysed, which is a staple object in numerical relativity used for such purposes. This is plotted on the left hand side of figure 6, for four different cases: supercritical evolution with low and high order artificial damping and subcritical evolution with low and high order damping. As the purpose of this convergence study is to calibrate the numerical accuracy of our finite difference method for the grid and the damping term, these simulations are done using classical matter, in other words, without any quantum mode functions added. From the figure, convergence for all cases is clear. It is also indicated on the graph which specific configurations are used to plot figures 3 and 4.
In order to check convergence with respect to the number of quantum modes added, another staple numerical relativity object is chosen; the ADM mass. This is calculated using the so-called "Schwarzschild-like" mass [48], which converges to the correct ADM mass very rapidly, sufficiently far away from the center of the spacetime where all the mass is contained. The "Schwarzschild-like" mass is defined by where r is an arbitrary radial coordinate, A sphere = 4πB(r, t)r 2 is the area of spheres at radial coordinate r and A(r, t) is the rr component in metric, eq. (2.1). Then the ADM mass is just  where R is some sufficiently large value of radial distance. Then, to measure the convergence of the system with respect to the number of added quantum modes, we define the deviation from the initially calculated ADM mass as time evolves This is plotted on the right hand side of figure 6 with various different numbers of quantum mode functions added, for the black hole formation case at t = 1.5 and r ≈ 8.0. The convergence is clear as the number of modes, N mode , is increased. Again, the specific amount of quantum modes used to obtain the results presented in figure 3 and 4 is indicated as well. While monitoring the convergence of the Hamiltonian constraints and the black hole mass provides a useful sanity check for the whole numerical setup, to ensure convergence of figure 4 specifically, we plot ∆T t t for various different numbers of quantum modes in figure 7. From this, we see that adding even more quantum modes would not qualitatively change the observed quantum deviations for the stress-energy tensor.

JHEP02(2022)183
We also want to check the convergence of the Choptuik scaling. For this, we shall check the convergence of the black hole mass differences defined in eq. (5.4) for figure 5. This can be seen in figure 8, where we plot the percentages of the black hole mass differences for various numbers of quantum modes. We see that increasing the number of quantum modes makes the black hole mass corrections converge rather quickly. There is little difference between the curves with N mode = 1600 and N mode = 2500 (N mode = 2500 being used in figure 5) for small initial amplitude a, and good convergence has also been achieved for the larger a in figure 5.

Conclusions
In this paper a new formalism has been introduced to collapse quantum fields into black holes. Coherent states have been utilised in order to relate the semiclassical system to the purely classical one. This allows one to replicate the already well-studied classical simulations in a semiclassical setting, hence giving rise to a convenient framework to study quantum effects in gravitational collapse. The simulations have been validated by convergence analysis for the L 2 -norm of the Hamiltonian constraint and the deviation in the ADM mass during the evolution. In addition, the convergence of the quantum effects themselves in the results has been illustrated.
The formalism has been implemented in numerical simulations to solve the coupled differential equations of the Einstein field equations along with the Klein-Gordon (opera-JHEP02(2022)183 tor) equation of a fully quantum mechanical massless scalar field in a spherically symmetric spacetime. It has been shown that using a finite number of quantum modes the quantum field operator is well-defined in a finite but sufficiently large spatial grid where interesting dynamics can take place. Results have been presented for initial conditions of both the subcritical and supercritical gravitational collapse, which are in agreement with the expected classical evolution. The presence of quantum effects have been explicitly demonstrated in the case of supercritical collapse by comparing the stress-energy tensor components in the semiclassical simulation to the ones in the corresponding classical simulation. The emerging quantum effects are found to be located around the apparent horizon.
The black hole masses in the semiclassical and classical systems have also been compared. The semiclassical black hole mass seems to be always smaller than its classical counterpart, and also the difference between them grows for increasing initial energy input in the simulation. The scaling behaviour of the black hole mass with the initial amplitude of the input wave-packet was investigated for both classical and semiclassical systems. It is found that they both obey the Choptuik scaling near the critical amplitude. The scaling exponent and critical amplitude have been computed using a least squares fitting, with the results in agreement with the accepted classical values.
There are various avenues of possible future work. Among others, one option is to use this formalism to find numerical evidence for the appearance of Hawking radiation in the supercritical case. In terms of the numerical sector, it has been mentioned that the higher order quantum modes introduce some instabilities to the simulations, which means that the semiclassical evolution can only be maintained for a certain time period before the accuracy is spoiled. While this has not posed existential difficulties for the current work, this might be crucial, for instance, to study the Hawking radiation. Various improvement can potentially be made for the numerical schemes, for example by employing the spectral methods or adaptive mesh refinement. These techniques might sustain stability for the semiclassical simulation for a longer period of time, which would enable a more extensive study of the emerging quantum effects in gravitational collapse. Another future direction is to consider self-interacting quantum matter fields, which, however, is more complicated as we then need to solve a coupled system between 2-point and higher correlation functions. Nevertheless, a Hartree approximation can be made, which often produces reasonable results [56].