ournal of C osmology and A stroparticle hysics J Evolution of black holes through a nonsingular cosmological bounce

. We study the classical dynamics of black holes during a nonsingular cosmological bounce. Taking a simple model of a nonsingular bouncing cosmology driven by the combination of a ghost and ordinary scalar ﬁeld, we use nonlinear evolutions of the Einstein equations to follow rotating and non-rotating black holes of diﬀerent sizes through the bounce. The violation of the null energy condition allows for a shrinking black hole event horizon and we ﬁnd that for suﬃciently large black holes (relative to the minimum Hubble radius) the black hole apparent horizon can disappear during the contraction phase. Despite this, we show that most of the local cosmological evolution remains largely unaﬀected by the presence of the black hole. We ﬁnd that, independently of the black hole’s initial mass, the black hole’s event horizon persists throughout the bounce, and the late time dynamics consists of an expanding universe with a black hole of mass comparable to its initial value.


Introduction
A proposed alternative to cosmic inflation is the idea that the universe underwent a bounce: a transition from a stage of contraction to expansion [1][2][3][4][5]. In a singular bounce, the universe passes through a classical singularity where the cosmological scale factor becomes small, curvature invariants blow up, and quantum gravity effects presumably become highly relevant to determining the future dynamics of the universe [6][7][8][9]. An alternative, which we focus on here, is a nonsingular bounce. For such cosmologies, so long as the spacetime curvature does not become Planckian, there is the possibility that quantum gravity effects could be subdominant to classical effects, in which case one may be able to describe the dynamics of the bounce using classical physics. Nonsingular bouncing cosmologies require violating the null convergence condition (NCC), which states that for all null vectors k µ , R µν k µ k ν ≥ 0 [4,5,[10][11][12]. In Einstein gravity, the NCC is equivalent to the null energy condition which is satisfied by most standard classical field theories [12]. Nonsingular bouncing cosmologies hence require non-standard matter terms, or modifications to Einstein gravity, for example, Horndeski theories including ghost condensation [13,14] or (cubic) Galileon/Horndeski models [15][16][17][18][19][20]. While perturbative studies of these theories suggest they may be free of ghost or gradient -1 -

JCAP09(2022)063
instabilities [16,19], less is known about which models will remain (strongly) hyperbolic through a bounce, when the solution is presumably not in the weakly coupled regime [21][22][23]. 1 An important open question is what happens in bouncing cosmologies in the inhomogeneous and non-perturbative regime. While there are several analytical and numerical studies of the dynamics of bouncing cosmologies during their contraction phase [25][26][27][28], there are relatively few studies of the dynamics of the bounce [19,[29][30][31][32], and none that consider the dynamics of black holes beyond the restriction to spherical symmetry [33]. Previous studies of black hole-cosmological bounces have either constructed initial data for black hole bouncing solutions [34], worked in a perturbative limit [35][36][37][38], or made use of analytic solutions (e.g. generalizations of the McVittie solutions [39][40][41]), that are limited by the fact that the metric evolution is prescribed ad-hoc, and from that the implied matter type and evolution is derived. The question of what happens to a black hole in a nonsingular cosmological bounce is particular salient for several reasons. On the one hand, the bounce necessarily requires a violation of the assumptions made in black hole singularity theorems and results on black hole horizons (namely the NCC) [42], so there is a question of whether the black hole will survive the bounce, or if the bounce mechanism will also reverse gravitational collapse, and if this will possibly lead to a naked singularity. On the other hand, one might also worry what the backreaction of the black hole's gravity will be on the bounce in the neighbourhood of the black hole. An extreme scenario would be if the bounce failed to happen in the vicinity of the black hole, possibly leading to a patch of contraction that grows into the expanding spacetime, as happens, e.g., in scenarios where the Higgs boson is destabilized during inflation, and goes to its true vacuum at negative energy densities [43,44].
Here, we address these questions by studying the nonlinear dynamics and evolution of black holes in a particular nonsingular bouncing cosmology (details of which are described below). Black holes can be expected to form during the contraction of matter and radiation dominated universes [36,45,46], and will generally be present from previous eras in cyclic cosmologies [47,48]. However, it is common to invoke a smoothing phase during contraction (e.g. ekpyrosis [47,49,50]), and argue that Hubble patches containing a black hole will be rare. Regardless, we view our work as serving two main purposes: (1) to study the dynamics and robustness of a nonsingular bouncing model when a very large perturbation, namely a black hole, is introduced, and (2) to explore the dynamics of the black hole and cosmological horizons during the bounce.
To avoid the difficulties related to finding a motivated theory that can give rise to bouncing solutions while also having well-posed evolution equations in the inhomogeneous regime, and thus being suitable for describing black hole dynamics, we will work with a bouncing cosmology model that incorporates a minimally coupled scalar field with a ghost field (i.e. a field which contributes to a negative cosmological energy density), to drive the bounce. While ghost fields are known to give problematic quantum mechanical theories (for a discussion of this in the context of cosmology see [51,52]; see also [53]), we take the point of view of [29,30,32] and treat the ghost field as an effective model for NCC violation. Quantum stability and unitarity is a distinct issue requiring a separate analysis (see, e.g., [54]). Unlike earlier work with this model, we do not restrict ourselves to cosmological spacetimes that have planar symmetry [32], or to small linear perturbations about a background bouncing spacetime [29,30]. Instead, we consider contracting cosmological initial data that contains a black hole, and work in an axisymmetric spacetime. This allows us to examine the effect that a large inhomogeneity has on the dynamics of the spacetime near and during the bounce. 1 We note that the model proposed in [19] is known to break down shortly after the bounce has ended [24].

JCAP09(2022)063
Following the growing number of studies making use of techniques from numerical relativity to study cosmological phenomena involving black hole dynamics [44,[55][56][57][58][59][60][61][62], we use numerical solutions to follow the evolution of different size black holes, both non-spinning and spinning, through a bounce, considering those both bigger and smaller than the minimum Hubble radius. Our main results are that the black holes persist to the expanding phase, and that the nonsingular bouncing model under study is fairly robust under large perturbations, in the sense that the local spacetime expansion around the black hole successfully bounces for all of the cases we explored. For large enough black holes, we find the black hole apparent horizon collides with the cosmological horizon, and temporarily disappears during the contraction phase. Nevertheless, the black hole apparent horizon eventually reappears (with finite radius event horizon throughout) and this does not disrupt the bounce at late times.
In principle a nonsingular, classical bounce could occur at any characteristic length scale that is larger than the scale at which quantum gravity effects become important (presumably the Planck scale: l P ∼ 10 −33 cm in geometric units). Given this, the length scale of a classical nonsingular bounce can still be extremely small compared to the typical length scale of say, an astrophysical black hole (e.g. in [48] the bounce happens at a typical length scale of ∼ 10 −25 cm ∼ 10 8 l P ). One may expect then that if any Hubble patch were to contain a black hole, that the black hole would be much larger than the minimum size of the Hubble patch. For example, even a black hole with a mass of m BH ∼ 10 15 g 2 at the bounce would still have a size of ∼ 10 20 l P ; this is orders of magnitude larger than the example bounce scale mentioned above. For this reason, we will be more interested in considering black holes whose size is comparable or larger than the bounce scale (which we take to be 1/|H min |, where H min < 0 is the maximum contraction rate).
The remainder of this paper is as follows. We discuss the nonsingular bouncing model we use in section 2. Our numerical methods and diagnostics for evolving the nonsingular bounce are outlined in section 3. Our numerical results are described in section 4, and we conclude in section 5. In appendix A, we discuss our numerical methodology in more detail, in appendix B, we define various quasi-local notions of black hole and cosmological horizons, and in appendix C, we provide an overview of the McVittie spacetime, an analytic solution to the Einstein equations of a black hole embedded in a cosmology, of which our numerical simulations can be seen as a generalization.

Ghost field model
We consider a theory that has two scalar fields φ and χ coupled to gravity:

JCAP09(2022)063
This model has a canonically normalized scalar field φ with a potential V (φ) = V 0 e −cφ , and a massless ghost field χ. The covariant equations of motion for (2.1) are Nonlinear, inhomogeneous cosmological solutions to the model (2.1) were studied in [32]. There, the authors considered a toroidal universe with a planar perturbation in one of the spatial directions. In this work, we consider an asymptotically bouncing FLRW universe with an initial black hole; see section 3 and appendix A for more details on our numerical methodology.
Strictly speaking, the ghost field should be stabilized by some mechanism at the quantum level. We choose to ignore this and treat (2.1) as a purely classical theory. As the equations of motion (2.2) have a well-posed initial value problem, 3 we expect the model should admit at least short-time classical solutions from generic initial data.

Homogeneous bouncing cosmology
Here we briefly review homogeneous, isotropic bouncing solutions for the system (2.2) (see also [30,32]), and discuss the values used for our asymptotic initial data. We work with harmonic coordinates (g µν Γ α µν = 0), so that the metric line element is The scalar field equations and Friedmann equations are then where the is the derivative with respect to the harmonic time coordinate t related to proper time by dτ ≡ a 3 dt, and H is the harmonic Hubble parameter where H is the Hubble parameter defined with respect to the proper time H ≡ (da/dτ )/a. We define effective energy densities ρ and pressures P for the two scalar fields: A requirement for having a nonsingular bounce is that w < −1, which coincides with violation of the NCC. 4 For example, if we consider the null vector When the NCC holds, we see thatḢ < 0, so that we have cosmic deceleration during expansion (H > 0), or cosmic acceleration during contraction (H < 0). When the NCC is violated, H > 0, and cosmic contraction can be slowed down, and even reversed to make a bounce. We also define effective equations of state for the fields φ and χ: From the Friedmann equations (2.4), one can determine that the energy density of the field f scales as ρ f ∝ a −3(1+w f ) .

Initial conditions
For our initial conditions, we first set the free initial data by superimposing the homogeneous initial conditions for the cosmological scalar fields and metric with the metric of a (rotating) black hole spacetime. We then solve the constraint equations for the full metric using a conformal thin sandwich solver [65] (see appendix A for more details). For the cosmological free initial data, we consider an initially contracting FLRW universe dominated by the canonical scalar field φ (that is, with the initial condition ρ χ ρ φ ). In this limit, with the potential V = V 0 e −cφ , φ can obey a scaling solution such that the effective equation of state is roughly constant and equal to w φ = c 2 3 √ 16π − 1 [30] (see more generally [66][67][68]). For c > √ 96π and V 0 < 0, the scaling solution in a contracting universe is ekpyrotic: the contracting solution is a dynamical attractor, and density perturbations are smoothed out in each Hubble patch during contraction [25-28, 49, 50]. In this limit w φ ≥ 1 = w χ , so if ρ φ > ρ χ initially during contraction, it remains so for all remaining time (recall ρ f ∝ a −3(1+w f ) ), and there cannot be a bounce. We instead consider the scenario where c < √ 96π and V 0 > 0 so that w φ < 1 = w χ , which is required in order to obtain a nonsingular bounce with the massless ghost field we consider [30,32]. As a result, the asymptotic, contracting, solution is not an attractor and the initial condition must be fine-tuned in order to keep w φ constant during the contracting phase. We justify this by noting that our main goal is to just explore the bouncing phase, and not to give a completely realistic description of a bouncing cosmology.
Setting c < √ 96π implies w χ > w φ , so the negative energy density of χ -which we choose to be initially negligible -grows faster than the positive energy density of the JCAP09(2022)063 canonical field during the contraction. Because of this, the total scalar field energy density ρ φ + ρ χ eventually goes through zero, and the sign ofȧ switches from being negative to being positive. At this point, the universe goes from contraction to expansion. From the Friedmann equations (2.4), we see that once expansion has begun, the ghost field energy quickly diminishes and becomes negligible again compared to the energy density of φ [30,32,69].
In (2.10), we present our choice of asymptotic FLRW initial data, which, as discussed above, is fine-tuned to allow for the asymptotic cosmological value of w φ to remain roughly constant during contraction up until the bouncing phase. The initial values for φ, φ , χ, χ , a, and a are: Here η 0 ≡ η(0) is the initial value of the ratio between the energy densities of the two scalar fields We compute ρ φ , ρ χ in the code using formulas (3.1) and (3.3).
In a similar fashion to [32], we choose c = √ 48π so that φ initially behaves like matter with w φ = 0. Such a matter-like contracting phase can generate scale invariant adiabatic perturbations that would seed structure formation in the early expansion phase.

Overview of numerical method and diagnostics
We evolve the system (2.2) nonlinearly using the harmonic formulation, and work with an axisymmetric spacetime. We spatially compactify our numerical domain, and evolve the boundary using the homogeneous FLRW equations of motion (2.4). See appendix A for a more thorough discussion on our numerical methods.
In order to characterize our results, we make use of several diagnostic quantities. We define the following stress-energy tensors so that the Einstein equations read From T µν (φ) and T µν (χ) we define the corresponding energy densities where n µ is the time-like unit normal vector to hypersurfaces of constant time. We additionally compute the local expansion rate where K is the trace of the extrinsic curvature on each constant t time slice (as in, e.g., [44,57]). We note that H K asymptotes to H at the boundary of our domain where r is the proper circumferential radius (see equation (3.7)). We define an effective scale factor on each time slice where γ 3 is the determinant of the (three-dimensional) metric intrinsic to each constant time hypersurface. We are mainly interested in computing (3.3) to (3.6) on the black hole surface, and at different coordinate radii far away from the black hole. For non-rotating black holes, we track their values as a function of the distance from the center of the black hole. We compare the values to their homogeneous counterpart given by (2.6) and (2.4d). In axisymmetric spacetimes, the coordinate radius on the equator r co is related to the proper circumferential radius r through the relation where γ zz is the value of the spatial metric along the symmetry axis. In spherical symmetry, eq. (3.7) reduces to the areal radius. To characterize the boundaries of black holes in our dynamical setting, we will consider two surfaces: event horizons and apparent horizons. The black hole event horizon is the boundary behind which null rays no longer escape to the asymptotic region. We compute its approximate location by integrating null surfaces backwards in time [70][71][72] (we restrict this to spherically symmetric cases, where it is sufficient to consider spherical null surfaces). We define the apparent horizon of the black hole, on the other hand, on each time slice, as the outermost marginally outer trapped surface, i.e. the surface for which the outgoing null expansion θ (l) vanishes and the inward null expansion θ (n) is negative and such that θ (l) > 0 immediately outside the black hole (and θ (l) < 0 immediately inside). 5 In analogy to black hole apparent horizons, we will also use marginally trapped surfaces to define the location of the cosmological apparent horizon. We will refer to this simply as the cosmological horizon, but we note that this is not to be confused with the event horizon or the particle horizon commonly used in cosmology. During the contracting phase, the cosmological horizon is defined as the surface for which the outgoing null expansion θ (l) vanishes, and the inward null expansion θ (n) is negative, but θ (l) > 0 immediately inside the cosmological horizon. During the expanding phase, the cosmological horizon is defined as the surface for which the ingoing null expansion θ (n) vanishes and the outward null expansion θ (l) is positive and such that θ (n) > 0 outside the cosmological horizon.
In a homogeneous spacetime, the cosmological apparent horizon is simply the sphere with coordinate radius equal to the comoving Hubble radius, R H = (aH) −1 , yielding an area JCAP09(2022)063 For our black hole spacetimes, we will always take the cosmological horizon to be centered on the black hole. This is because we are interested in the dynamics in the vicinity of the black hole, and it is this surface that is the most relevant to understanding the behavior of the black hole horizon. See appendix A for more details on our numerical implementation and appendix B for more details on the various definitions of horizons we use.
From the area of the black hole apparent horizon A B , we define an areal mass M A ≡ A B /(16π). The spacetime we study here violates the NCC, and thus we expect to find instances where M A decreases. Similarly, the second law of black hole thermodynamics states that so long as the NCC is satisfied, the area of a black hole event horizon must increase into the future [73]. This can be extended to the cosmological setting assuming that the universe does not again collapse, and a notion of infinity can be defined [74]. However, here we are evolving a black hole in a spacetime that violates the NCC, and find that the event horizon does decrease in area.
The cosmological and black hole apparent horizons that we find on each time slice can also be thought of as foliations of three dimensional surfaces called holographic screens [75][76][77] or Marginally Trapped Tubes (MTTs) [78] in general, and dynamical horizons [79][80][81] if they obey certain extra conditions (we review the definitions of these concepts in appendix B). Though one can formulate area laws for these surfaces, in spherical symmetry they do not place any constraints on whether the area increases to the future. We keep track of the MTTs corresponding to the cosmological and black hole apparent horizons, and in particular, compute when they are spacelike or timelike in nature.
For the black holes, we compute the equatorial circumference of the horizons c eq , and define their corresponding equatorial radii r eq = c eq /2π, which in the case of spherical symmetry is also equal to the areal radius, r A = A B 4π . When studying rotating black holes we can also associate an angular momentum to the apparent horizons whereφ i is the axisymmetric Killing vector, and, using the Christodoulou formula, we can define a mass Since the scalar fields do not carry any angular momentum in axisymmetry, the total angular momentum of the black hole remains constant throughout the evolution of our spacetime. Thus, we will only be interested in the total mass and circumferential radius of the black hole.

Results
We begin by studying the evolution of non-spinning black holes in an asymptotically bouncing universe (section 4.1-4.3) using the method described in section 3. Though we do not explicitly enforce spherical symmetry, we find no evidence of any instabilities that break that symmetry if our initial data respects it. We consider non-spinning black holes in section 4.4. We find that the qualitative behavior of our solutions can be divided into two regimes, which can be distinguished by the ratio of the areal radius of the initial black hole horizon, r BH,0 and the minimum size of the Hubble radius of the background cosmology R H,min ≡ min t |1/H| = −1/H min (where H min < 0 is the maximum contraction rate). When -8 -

JCAP09(2022)063
R H,min /r BH,0 3.5, the black holes pass through the bounce freely. When R H,min /r BH,0 < 3.5, we find that the locally defined cosmological and black hole apparent horizons merge, and cease to exist for a period of time during the contracting and bouncing phase. We note that the horizons merger at R H,min /r BH,0 > 1, as the black hole grows in size during the contraction phase (see figures 3, 6; we will discuss this more in the following subsections).
For every initial data setup we considered, we find that the black hole continues to exist after the bounce phase ends: the late-time evolution always consists of a black hole in an expanding universe with the ghost field energy density decreasing at a faster rate than the canonical scalar field energy. Moreover, we find that the late time black hole mass remains similar to the initial black hole mass, regardless of the ratio of the initial black hole radius and minimum Hubble patch radius. In the following sections, we quantify these observations and extrapolate our findings to the regime where the Hubble radius shrinks to a much smaller size compared to the radius of the black hole.

Small black hole regime
We first consider solutions where R H,min /r BH,0 3.5 (see above for definitions). In figure 1, we show the Hubble parameter (left panel) computed from (3.4) and the ratio of scalar fields (right panel) computed from (2.11), (3.1) and (3.3) as a function of harmonic time for different coordinate radii. We also plot the value these quantities take at spatial infinity, where we assume homogeneous FLRW boundary conditions (see section 2.1). While the bounce seems to be pushed to slightly earlier harmonic times when the black hole is present, most of the local cosmological evolution remains unaffected by the presence of the black hole and follows the same qualitative evolution as the background cosmology (section 2.2). To determine how the cosmology is affected in a region close to the black hole, in figure 2 we plot the spatial dependence of η and H K /|H min | as a function of distance again along the equator at different times. Although the local expansion rate and the ratio of the energy densities can differ from their background values by up to 15 − 60% and 9 − 16% respectively, beyond r ∼ 10 − 25r BH,0 both quantities quickly asymptote to their respective background values. Note that the coordinate radius differs from the proper radius by the local scale factor; see eq. (3.7). The effective scale factor computed from (3.6) at different coordinate radii is plotted in figure 12 (see appendix A). Again we find that far enough from the black hole, the value of scale factor remains largely unaffected by the presence of the black hole. We caution that these quantities will also be subject to gauge effects -in particular from our choice of the lapse function (see section A). As we describe below, towards the end of the simulations we find strong variation in the rate of which time advances at different spatial points.
We next present several results regarding the behavior of the area of the black hole, as measured by either the event or apparent horizon. Naively, one expects the accretion of the canonical/ghost field to result in an increase/decrease in mass of the black hole [82]. That being said, it is less clear how a black hole embedded in a cosmology driven by a canonical/ghost field may behave [83,84]. Figure 3 depicts the evolution of the black hole's areal radius. We find that during the contracting phase prior to the bouncing/NCC violation phase, the canonical scalar field energy density exceeds that of the ghost field; see figure 1 and the solid purple curve in figure 2. The black hole's proper area increases during this time; (first region in figure 3 where H < 0,Ḣ < 0). Once the bouncing phase starts (t|H min | ∼ 120 in figure 1), the black hole starts to shrink as one may expect since the ghost field energy in this regime is comparable to the canonical scalar field energy density (second region where H < 0 andḢ > 0 in figures 1 Notice that the black hole reaches its maximum size slightly before the universe at large scales bounces, as the ghost field begins to dominate at an earlier time the closer one gets to the black hole horizon. The slight difference in the maximum absolute value of the FLRW value of H K /|H min | at t|H min | ∼ 120, 400 is due to numerical error in our integration.   and 2). Near the end of the bouncing phase the universe is expanding (region where H > 0 andḢ < 0), yet the black hole's size is still shrinking in this region, as the ghost field energy density still dominates over the canonical scalar field energy density in the region near the black hole (in other words, η < 1 in the region close to the black hole, see figure 2), although at an increasingly slower rate as the ghost field energy density quickly diminishes in time.
After the end of the bouncing phase, the universe continues to expand, the ghost field decays to dynamically irrelevant values, and the black hole begins growing in size (fourth region where H > 0 andḢ < 0 and the dotted purple curve). The left panel of figure 3 also shows the areal radius of the cosmological horizon. We see that during the contracting phase, the cosmological horizon shrinks from r C,0 = 75r BH,0 to a minimum radius of r C,min = 4.34r BH,0 at t ∼ 50r BH,0 . This is similar to the value the Hubble radius (R H ≡ |1/H|), would shrink to in the absence of a black hole. This value is indicated by the diamond in figure 3. From this we conclude that -at least in this regime -the presence of the black hole does not qualitatively change the dynamics of the spacetime. Past this point of closest encounter, the cosmological horizon tends to r C → +∞ which defines the location of the bounce (lim H→0 1/H = ∞). Once the universe switches from contraction to expansion, the cosmological horizon is defined as the location where the ingoing null expansion vanishes and outgoing null expansion is positive. After the bounce, the cosmological horizon at first shrinks to a minimum size before re-expanding to +∞. We note that the areal radius of the cosmological horizon is no longer symmetric about the bounce once a black hole is present.
We also compute the signature of the MTTs associated with the horizons (see appendix B for definitions), which we plot in figure 3.
First we study the properties of the black hole MTT in more detail. Using the terminology of appendix B, the black hole is a future marginally trapped tube foliated by future marginally outer trapped surfaces (alternatively called a future holographic screen). The area law of dynamical horizons states that if the MTT is spacelike (i.e. if it is a dynamical horizon), then -11 -JCAP09(2022)063 the area of the black hole should increase in the outward radial direction, while if the MTT is timelike (i.e. we have a timelike membrane with Θ (n) < 0), then the area should increase into the past. Looking at figure 3 we find that (as expected) these laws are obeyed at all times, even during the bouncing phase.
We next look at the cosmological horizon. We consider the contracting and expanding phases separately. During the contracting phase, the cosmological horizon is a MTT foliated by future marginally inner trapped surfaces (alternatively, it is a future holographic screen). From the area law of future holographic screens [76,77], we expect the cosmological horizon to obey the same area law as the black hole during the contracting phase. Our findings agree with this expectation: we find that the cosmological horizon is timelike when it decreases in time and spacelike when it increases in the outward direction. During the expanding phase, however, the cosmological horizon ceases to be a MTT. Instead, we find that it satisfies the definition of a past holographic screen (as the ingoing null expansion now vanishes). From [76,77], we still expect its area to increase in the future on timelike portions and in the outward direction on spacelike portions. Again we find that this is satisfied at all times during the expanding phase.
We conclude by looking at the event horizon shown in the right panel of figure 3. Our main finding here is that the event horizon no longer lies outside the apparent horizon at all times. This is a result of the violation of the NCC [42]. Interestingly, this behavior begins not during the bouncing phase of cosmological evolution (between the two dashed grey lines) when the NCC is violated, but before the bouncing phase has begun. This is because the event horizon is not a quasi-local quantity, so it can "anticipate" the bouncing/NCC violation phase. In general, we find that the event horizon always increase until it crosses the apparent horizon of the black hole, after which it decreases. Once the bouncing phase ends, the event horizon crosses the apparent horizon again, after which it starts increasing and remains larger than it for all future times.
We were not able to evolve the spacetime to arbitrarily large proper times. We ascribe this to gauge artefacts which impede the stable numerical evolution of the solution. In particular, the lapse function appears to become distorted in the spacetime region between the black hole and the asymptotically homogeneous regime, which causes that interior region to advance in time much faster compared to elsewhere in the simulation. (This is evident in the rightmost panels of figures 1, 4, and 12.) That being said, based on the simulations we have run, we conjecture that the black hole asymptotes to close to its initial mass at t → ∞ with no significant gain or loss of energy. That is, the end state is described by a black hole embedded in an expanding matter like FLRW universe with a negligible amount of ghost field and matter energy density. This is illustrated in figure 13 of appendix A where we consider a black hole with half the mass of the one depicted in figure 3, i.e. we consider a black hole such that the ratio of the minimum Hubble radius to the initial radius of the black hole is R H,min /r BH,0 = 8.69. Figure 13 shows that, overall, the black hole's size changes by a negligible amount. In this particular case, the final size of the apparent horizon of the black hole is ∼ 6% larger that its initial value, the small difference being an artefact of the initial data. More importantly, figure 13 also shows that the event horizon asymptotes to the apparent horizon at late times.

Large black hole regime
We next consider solutions where R H,min /r BH,0 < 3.5. The nonlinear evolution of one particular case is shown in figure 4. As is the case for the lower initial mass evolutions (figure 1), we see  that the cosmological evolution remains unaffected far away from the black hole. The bounce is pushed to even earlier times, as one may expect since a large black hole could presumably accelerate the rate of cosmological contraction. Figure 5 shows that in the region near the black hole apparent horizon, the local expansion rate and the ratio of the energy densities now differ from their background value by up to 15-75% and 13-60%. Beyond r ∼ 2-12r BH,0 , both quantities asymptote to their respective background values. The behavior of the black hole and cosmological apparent horizons, which is shown in figure 6, is qualitatively different for the large black hole initial data as compared to the small black hole initial data (R H,min /r BH,0 3.5). Similar to the cases studied in the section 4.1, Spacelike Timelike BH EH FLRW Figure 6. Same as figure 3 but for a black with initial mass such that the Hubble radius of the background cosmology R H ≡ |H −1 | shrinks from an initial value of R H,0 = 75r BH,0 to 2.17r BH,0 (here r BH,0 is the initial black hole radius). Notice that the location where H = 0 (that is, where the Hubble radius diverges) does not exactly coincide to where the cosmological horizon blows up, as the cosmological horizon is measured locally (in the interior of the computational domain), while H = 0 is determined by the asymptotic cosmological evolution. For more discussion on how we define the cosmological horizon, see section 3.
the cosmological horizon shrinks at first. Unlike those earlier cases though, it eventually merges with the expanding black hole apparent horizon. Following the merger, the spacetime has no apparent horizons for some time until they re-emerge. After that, the cosmological and black hole apparent horizons follow a similar trajectory to the horizons studied in section 4.1 during the cosmological expansion phase.
The merging of black hole and cosmological apparent horizons has been observed in McVittie spacetimes [85,86] (see also appendix C) and can be interpreted the following way. As the apparent horizon of the black hole grows and the cosmological horizon shrinks during the contraction of the universe, we reach a point in time at which the black hole horizon coincides with the cosmological horizon. At this point, one cannot distinguish between the black hole and the cosmological horizon (recall that during the contraction the outward null expansion is negative outside of the cosmological horizon). A finite time later, before the bounce, but after the background Hubble radius reaches its minimum size, the effective Hubble radius has increased to a sufficiently large value so that the black hole solution again fits within the cosmological horizon. At this point, the cosmological and black hole apparent horizons reappear. We note that the black hole event horizon persists throughout the evolution of the spacetime, so in this sense the black hole never disappears; see figure 6. We next investigate the physical properties of this process in more detail.
We first address the question of whether a naked singularity forms after the black hole and cosmological horizons collide [83,85,86]. The formation of a naked singularity would signal a breakdown of the theory -either through the formation of a blowup in curvature, or through necessitating new boundary conditions to be set at the singularity boundary [87]. Our simulations suggest no naked singularity is formed. More concretely, the outward null expansion during this period is negative everywhere, so the entire spacetime is essentially trapped, and no new boundary conditions need to be specified. In particular, we can continue to excise a central region corresponding to the inside of the black hole. Additionally, considering the event horizon shown in figure 6, we see that it remains finite at all times. Note that just like in the case studied earlier in section 4.1, the event horizon is smaller than the apparent horizon before, and during the bouncing phase, and turns around when it crosses the apparent horizon. We next consider the behavior of the marginally (anti-)trapped tubes and their signature, shown in figure 6. Note that while the black hole MTT is spacelike and increasing in time before it merges with the cosmological horizon, when it reappears from the merger, its signature remains spacelike even though its area continues to decrease in time. Since the area of the black hole always increases in the outward radial direction, this implies that while the outward direction points into the future before the merger, it points into the past when it reappears. The black hole apparent horizon undergoes another signature change at the bounce (indicated by the grey vertical solid line) after which it behaves like the case studied above (i.e. the signature of the horizon becomes timelike, and decreases as we evolve forwards in time). Similarly, we find that the cosmological horizon follows the same trend as the case in section 4.1, except for a brief period of time just before it merges with the black hole apparent horizon: here the horizon signature becomes spacelike. We see that the cosmological horizon and black hole apparent horizons have the same signature when they annihilate and re-emerge.
A natural question to ask is whether the collision of the apparent horizons during the contraction phase is an artefact of the particular matter model we use, or is a more general consequence of a contracting universe. To explore this, we consider the same initial conditions as the ones used in figure 4, but now evolve only with the canonically normalized scalar field. The results of this are plotted in figure 7. We find that during contraction, the apparent horizons, with and without the presence of a ghost scalar field, behave in a similar fashion. In both cases, the black hole apparent horizon merges with the cosmological horizon at the same areal radius. This is in line with our earlier observation that the black hole horizon's size exceeds the cosmological horizon before the bouncing phase starts, (i.e. before the ghost field has a significant impact on the evolution of the system). The black hole and cosmological apparent horizon merge earlier by around t ∼ 2|H min | in the case of contraction without the ghost scalar field. This is consistent with the notion that the accretion of the ghost field

JCAP09(2022)063
should slow down the rate at which the black hole can grow in size, which would delay the time of merger of the two horizon.
Finally, we note that the signature of the cosmological horizon becomes spacelike in this setup just before merging with the black hole horizon for both cases. During this phase of evolution, the cosmological horizon is a dynamical horizon whose area decreases with time.

Dependence on black hole size
In this section, we explore in more detail how the properties of the spacetime during the bounce change as a function of R H,min /r BH,0 .
As described in section 4.1, for initial data where R H,min ≈ 4.34r BH,0 the black hole apparent horizon persists through the whole bounce, and the spacetime evolution near the black hole qualitatively resembles the asymptotic cosmological evolution. As this behavior will hold to an even greater degree for smaller black holes (relative to R H,min ), we are more interested in the opposite regime, considering larger black holes. As mentioned in section 1, for astrophysical black holes we expect R H,min r BH . We find that, when R H,min 3.5r BH,0 , (see section 4.2), the black hole apparent horizon collides with the cosmological horizon while the universe is still contracting. In this section, we therefore explore how this behaviour changes as one increases the initial mass of the black hole. We note that for numerical reasons, 6 we will restrict to evolutions where R H,min > 0.86r BH,0 . However, as we argue below, we already see some consistent trends as r BH is varied within this regime.
In the left panel of figure 8, we plot the radius of the black hole apparent horizon normalized by its initial value as a function of time. For evolutions where the black hole and cosmological MTTs do not collide, we find that although the area of the apparent horizon always reaches it maximum and minimum values at around the same harmonic time (t ∼ 110R H,min for the maximum value, and t ∼ 440R H,min for the minimum value), the value the maximum and minimum take does change as a function of initial black hole area. Independently of the black hole's initial size, the area of the apparent horizon is close to one around the bounce or in other words around the time where the total energy density of the background cosmology is zero. In the low mass regime, the variation in the black hole's size increases with increasing initial black hole area.
However as the initial size of the black hole increases, the maximum change in the radius of the apparent horizon eventually peaks at a value of r AH,max /r AH,0 ∼ 2.6. In this case, the ratio of the minimum Hubble radius of background cosmology to initial radius of the black hole corresponds to the threshold beyond which the horizons merge. Beyond this peak, although the horizons merge at successively earlier times (and always before the bouncing phase starts), with increasing initial black hole radius, the relative increase in the radius of the apparent horizon when the horizons merge saturates at a value of r AH,max /r AH,0 ∼ 2.5. Within the range of masses we were able to evolve, the apparent horizons always reappear, from which we conjecture that the presence of black holes in bouncing cosmologies do not disrupt the bounce. We were not able to evolve the space time to arbitrarily later proper time but based on all the simulations we have run, we conjecture that the black hole asymptotes to close to its initial radius as t → ∞.

JCAP09(2022)063
In the right panel of figure 8, we plot the radius of the black hole event horizon normalized by its initial value as a function of time. We do not compute the evolution of the event horizon past the bounce for black holes with initial black hole radius such that the minimum Hubble radius is smaller than R H,min < 2.90r BH,0 , as for those cases the event horizon cannot be located to the desired accuracy (see appendix A for more details on the computation of the event horizon). For the set of initial radii we do compute, we find that the area of the event horizon reaches a maximum at successively earlier times with increasing initial black hole radius, always before the bouncing phase starts and always when the event horizon crosses the apparent horizon of the black hole. Beyond this point, the event horizon decreases in size, until it crosses the apparent horizon again, after which it starts increasing. This minimum happens at successively earlier times with increasing initial black hole radius. While the maximum size of the event horizon throughout the evolution increases with increasing initial black hole radius, the minimum decreases.
We next argue that the behavior of the event horizon in the region leading up to the bounce (where H < 0), can be at least qualitatively captured by studying null rays in the background FLRW spacetime. The reasoning is as follows: it is reasonable to assume that in the regime where R H /r BH,0 1, the evolution of null rays near the black hole horizon will not be greatly influenced by the background cosmological evolution. Likewise, we assume that in the regime where the black hole is "large" (R H /r BH,0 1), the trajectories of null rays exterior to the black hole are more influenced by the cosmological evolution, 7 and in the background FLRW spacetime, during the contraction phase, outward radial null rays have decreasing proper radius when they are inside the Hubble radius. Following this line of thought, we integrate null rays backward in time in the background FLRW spacetime given by eq. (2.3), starting from the latest time for which H < 0 and R H = r BH,0 . Figure 9 shows the trajectories of a few such null rays for different ratios of R H,min /r BH,0 . We find that the proper radius of the null rays increases (as we go backwards in time) until the ray crosses R H , after which it decreases. This is consistent with the behavior of the event horizon in the right panel of figure 8 and suggests that, at least for this part of the evolution, the size of the black hole is determined by the evolution of the background cosmology. This simple calculation also shows that as R H,min /r BH,0 decreases, the maximum radius of the null ray increases. This agrees with what we see in our full numerical simulations. Extrapolating this trend to arbitrarily small R H,min /r BH,0 suggests that for arbitrarily large black holes the peak of the event horizon will diverge. However, we are working with a cosmological solution that has undergone an infinite number of e-folds of contraction to the past (see section 2.2). If one were to consider a bouncing model that had only had a finite period of contraction (for example if we considered a cyclic cosmology [47]), then the maximum of the event horizon would always be finite.
Evolving forward in time, into the region where the universe is expanding (H > 0), we find that the event horizon continues to decrease until it crosses the apparent horizon, at which point it begins to increase in size. However, this behavior can not be captured by integrating the null geodesics in the background spacetime, which suggests that the influence of the black hole on the geometry is more relevant when H > 0, and for radii less than r BH,0 . Due to numerical issues, we are unable to evolve far enough in time to determine if the minimum of the event horizon keeps decreasing and eventually reaches a point where the event horizon ceases to exist as R H,min /r BH,0 → 0.  . Proper radius of outward, radial null rays in the background FLRW spacetime (described in section 2.1). The left panel shows an example null ray that begins at a specified radius r null,0 , increases until it crosses the Hubble radius (during contraction), and then decreases until it reaches the Hubble radius again at r null,0 . The right panel shows the same thing for different ratios of the minimum Hubble radius to r null,0 .
Finally, we note that (as is shown in figure 13) one expects the apparent and event horizons to converge to the same value at late times, but for reasons mentioned earlier in this section, we are not able to evolve long enough in time to show this happens for initial data with R H,min < 4.34r BH,0 .

Spinning black holes
Up to this point, we have only considered non-spinning black holes (spherically symmetric spacetimes). However, our methods can be applied equally well to spinning black holes spacetimes. We have considered several such cases, finding the same qualitative behavior as for non-spinning black hole initial data. We illustrate this with a representative example case: initial data where the black hole is initially spinning with a dimensionless spin value of a 0 = 0 a 0 = 0.5 Figure 10. The circumferential radius of the black hole in figure 6 with zero spin (a = 0) compared to dimensionless spin of a = 0.5. a 0 = 0.5. As we find little difference compared to the spacetimes with non-spinning black holes, here we only present the results for a black hole with R H,min = 2.17r BH,0 , and the same mass as in section 4.2. Figure 10 compares the circumferential radius of the black hole along the equator for the spinning and non-spinning cases. We find that the addition of spin causes the horizons to merge at a slightly later time as compared to a comparable non-spinning case. We do not plot the behavior of the asymptotic background cosmology as it is the same regardless of whether the black hole is spinning or not. As was mentioned in section 3, the angular momentum of the black hole is constant since the scalar field does not carry angular momentum.

Discussion and conclusion
We have considered the first numerical evolution of black holes through a nonsingular bouncing cosmology. As in [30,32], we worked with a model that has two scalar fields: a canonically normalized field with an exponential potential and a ghost field. We have additionally considered asymptotically cosmological initial data that is tuned to allow for a matter-like (effective equation of state w = 0) contraction which is then followed by a bounce that ends with cosmological expansion. In [32], translational symmetries were assumed in two spatial directions, which precludes the formation of black holes. By contrast, in this work we considered axisymmetric spacetimes which allowed us to study the behavior of black holes through a bounce. While only a small fraction of Hubble patches are expected to have a black hole during the late stages of ekpyrotic contraction [4,36,46], our setup allows us to examine the robustness of the ghost-field bounce, which in turn serves as an effective classical model of NCC violation. We found two qualitatively different kinds of spacetime evolution, which depended on the ratio of the minimum Hubble radius of the background cosmology to the initial radius of the black hole. For black holes with initial radius smaller than ∼ 3.5 times the minimum size of the Hubble radius of the background cosmology, the black hole passes through the bounce freely and the background cosmology remains largely unaffected (see section 4.1). Beyond this limit, we found that while regions far away from the black hole still bounce freely, regions close to the black hole evolve differently (see section 4.2). In particular, we found -19 -

JCAP09(2022)063
that during the contracting phase, the cosmological horizon and the black hole apparent horizon merge and cease to exist for a brief period of time. Some finite time later, before the bounce but after the background Hubble radius reached its minimum size, the cosmological and black hole apparent horizons separate. Within the range of masses we considered, we found that the black hole size (as measured by its horizon radius), varies significantly during its evolution. However, regardless of the initial mass of the black hole, we found that the late time evolution consists of a black hole in an expanding universe with a mass similar to its initial value. Although we were not able to evolve spacetimes where the Hubble radius shrinks to a much smaller size compared to the radius of the black hole, we conjecture that the black hole always survives through the bounce. This means that black holes created (or already present) in the contraction phase [36,46] can persist to have observational consequences in the post-bounce era.
We found instances where the event and apparent horizons decrease as a result of our spacetime violating the NCC. Independently of the NCC being violated, we found that in the regime where the black hole and cosmological apparent horizon collide, the latter becomes spacelike shortly before merging with the black hole. This is consistent with the observation that the signature of the marginally (anti-)trapped tubes changes such that any merging/reappearing pair of horizons always has the same signature.
Finally, we point out a few directions for future research. One would be to study the dynamics in a setup where the asymptotic cosmology is not prescribed. For example, this could be accomplished by considering a toroidal/periodic setup, and then considering a "lattice" of black holes [34,55]. This setting would allow for the study of the impact of black holes, as well as other perturbations on the overall dynamics of the bounce. While small perturbations have not been found to appreciably change the dynamics of a nonlinear bounce when translational symmetries are assumed [32], it would be interesting to see if perturbations could be more disruptive in the presence of a black hole, and in a less-symmetric spacetime that does not preclude large-scale anisotropies.
Another direction would be to consider other models of cosmological bounces. While we believe that the main conclusions we find here do not depend strongly on the details of the bounce model, it would be interesting to determine what differences would result from potentially more realistic models of a bounce. As we mention in the Introduction, the cosmological bounce scale may be many orders of magnitude smaller than the initial size of a primordial black hole. Due to the numerical instabilities (as described in section 4.2), we were unable to carry out evolutions in the regime R H,min /R BH,0 1. It would be interesting to see if our results still hold in this limit.
Another interesting question is the degree to which a ghost field, which can reverse cosmic contraction, may similarly affect gravitational collapse and singularity formation in a black hole interior. NCC violating fields such as ghost fields have been used to construct singularity free black hole-like solutions, such as wormholes [88][89][90][91], so it is not entirely implausible that there could be nontrivial dynamics near the center of a black hole that accretes a ghost field. In this study, we ignored the dynamics deep inside the black hole, excising that region from our domain. Exploring this would require coordinates better adapted to studying the interior of black hole spacetimes, such as null coordinates [92].

A Numerical methodology
We solve the equations of motion (2.2) using the generalized harmonic formulation as described in [93]. The numerical scheme we use follows that of [94], which we briefly summarize here. We discretize the partial differential equations in space, using standard fourth-order finite difference stencils, and in time, using fourth-order Runge-Kutta integration. We control high frequency numerical noise using Kreiss-Oliger dissipation [95]. We use constraint damping to control the constraint violating modes sourced by truncation error, with damping parameter values similar to those used in black hole evolutions using the generalized harmonic formulation [93]. We fix the gauge freedom by working in harmonic coordinates, x α = 0. During the expansion phase, we dynamically adjust the time step size in proportion to the decreasing global minimum of 1/α where α is the lapse (this would be α = 1/a 3 in a homogeneous FLRW universe, see eq. (2.3)) in order to avoid violating the Courant-Friedrichs-Lewy condition [57,96].
Following [97], we dynamically track the outer apparent horizon of the black hole, and excise an ellipsoid-shaped region interior the horizon. We typically set the ratio of the maximum ellipsoid axis to the maximum black hole radial value to be 0.6.
We compute the event horizon by integrating null surfaces backwards in time [70][71][72] (we restrict this to spherically symmetric cases, where it is sufficient to consider spherical null surfaces). Since we are not able to evolve the spacetime to infinite proper time (at which point the event and apparent horizon would coincide), we cannot precisely determine the final position of the event horizon. Instead, we use the apparent horizon as the approximate location of the event horizon and choose a range of initial guesses around this value. For two surfaces initially separated by 2.5r BH,0 , we find that their separation decreases to 0.1r BH,0 within ∼ 4 × 10 −3 |H min | −1 when evolving the null surfaces backwards in time, after which we consider the location of the event horizon to be accurate to the desired accuracy. Note that the separation rapidly decreases when integrating backwards in time, a direct consequence of the divergence of the null geodesics going forward in time [71].

JCAP09(2022)063
We additionally make the use of compactified coordinates so that physical boundary conditions can be placed at spatial infinity [97]: so thatx i = 1 corresponds to x i = ∞. Unlike in [97] though, we work in an asymptotically FLRW spacetime instead of an asymptotically flat spacetime, similar to what is done in [44].
That is, at our spatial boundary we set where the lapse α(t) = a(t) 3 ; and the scale factor a(t) satisfies the Friedmann equations, eq. (2.4).
We use Berger-Oliger [98] style adaptive mesh refinement (AMR) supported by the PAMR/AMRD library [99,100]. Typically our simulations have 9-12 AMR levels (using a 2 : 1 refinement ratio), with each nested box centered on the initial black hole and between 128 and 256 points across the x-direction on the coarsest AMR level. The interpolation in time for the AMR boundaries is only third-order accurate, which can reduce the overall convergence to this order in some instances. As we restrict to axisymmetric spacetimes, we use the modified Cartoon method to reduce our computational domain to a two-dimensional Cartesian half-plane [97].
We construct initial data describing a black hole of mass M (t = 0) = M 0 in an initially contracting FLRW spacetime described in section 2.2. We solve the constraint equations using the conformal thin sandwich formalism, as described in [65]. More precisely, we choose the initial time slice to have constant extrinsic curvature is given by (2.4d), and the initial values for {φ, φ , χ, χ , a, a } are fixed by (2.10) (a similar approach was employed in [44,57]). Without loss of generality, we choose the initial value of the ratio between the energy density of the φ and χ fields and V 0 to be such that during the contraction phase, the Hubble radius of the background cosmology R H ≡ |H −1 | shrinks from an initial value of R H (t = 0) = 75r BH,0 to 4.34r BH,0 (here r BH,0 is the initial black hole radius). We considered a range of initial black hole masses, keeping the initial ratio of Hubble to black hole radius to be 75, but changing the minimum Hubble to initial black hole radius ratio from 4.34r BH,0 all the way to 0.87r BH,0 . We also study some black holes with an initial dimensionless spin of a 0 = 0.5.
Finally, we present a convergence test of our code and setup. In figure 11, we present the time evolution of the apparent horizon of the black hole and the norm of the constraint violations C α ≡ x α integrated over the coordinate radius r ≤ 265M 0 , for a non-spinning black hole with initial mass such that R H,min = 1.45r BH,0 , for different numerical resolutions. For this case, the lowest resolution is 128 points across the x-direction on the coarsest AMR level with 10 levels of mesh refinement and a spatial resolution of dx/M 0 ≈ 0.004 on the finest level. The medium and high resolutions correspond, respectively, to an increased resolution of 3/2 and 2× that of the lowest resolution run. We find that the constraints converge to zero at roughly third order. This is because the convergence is dominated by the third order time interpolation on the AMR boundaries. The medium resolution in the convergence study is equivalent to the resolution we use for all the other cases studied here. We place the mesh refinement such that the radius of the black hole resides inside the finest AMR level initially. During the evolution, the mesh refinement is adjusted according to truncation error estimates to maintain roughly the same level of error. The effective scale factor |γ 3 | 1/6 computed from (3.6) for a black hole with initial mass such that the Hubble radius of the background cosmology R H ≡ |H −1 | shrinks from an initial value of R H,0 = 75r BH,0 to 4.34r BH,0 (left)/2.17r BH,0 (right). The solid line shows the corresponding background solution and the dashed and dash-dotted lines the values at different coordinate radii. The vertical grey line is the time at which the black hole reaches its maximum mass as observed by the apparent horizon.

B.1 General definitions and properties
Nonsingular classically bouncing cosmologies require the violation of the NCC [4,5,10,11,30]. The NCC plays a fundamental role in the classical area law for black holes [42,101]. Given this, we pay particular attention to the dynamics of the black hole horizon in our simulations. In addition to the event horizon (which can only be computed once the whole spacetime is known [72]), there are several other quasi-local definitions of black hole horizons which we measure: dynamical horizons [80,81,102], apparent horizons [33,[79][80][81][102][103][104], and holographic screens [75][76][77] (also called marginally trapped tubes [78]). For completeness, we collect the definitions and some of the basic properties of these horizons in this appendix. Wherever applicable, we also discuss how these definitions can be extended to define cosmological horizons. We refer the reader to [33, 75-77, 79-81, 102-104] for more thorough reviews on this subject.
Trapped surfaces and apparent horizons. Let S be a smooth, closed, orientable spacelike two-dimensional submanifold in a four-dimensional spacetime (M, g ab ). We then define two linearly independent, future-directed, null vectors normal to S, normalized 8 such that where by convention l α and n β are respectively the outgoing and ingoing 9 null normals. The two-metric induced on S isq αβ = g αβ + l α n β + n α l β , (B.2) and the null expansions are defined as The

JCAP09(2022)063
A black hole apparent horizon is a future marginally outer trapped surface. Within the context of cosmology, the cosmological apparent horizon of an expanding FLRW spacetime is a past marginally inner anti-trapped surface. For a contracting FLRW spacetime, the cosmological apparent horizon is a future marginally inner trapped surface.
Dynamical horizons and holographic screens. We now have all the ingredients to introduce the concept of a Marginally Trapped Tube (MTT): a MTT is a smooth, threedimensional submanifold that is foliated by MTSs. If a MTT is everywhere spacelike, it is referred to as a dynamical horizon [80,81]. If it is everywhere timelike, it is called a timelike membrane (TLM). 10 Finally, if it is everywhere null then we have an isolated horizon.
We next outline the various ingredients that go into the area law of dynamical horizons (the quasi-local horizon that appears the most frequently in our numerical solutions). It is straightforward to derive an area law for purely spacelike or purely timelike dynamical horizons. Consider first the spacelike case. Let H be a dynamical horizon and S be a member of the foliation of future marginally trapped surfaces. Since H is spacelike, we can define a future-directed unit timelike vector normal to H,τ a and a unit outward pointing spacelike vector tangent to H and normal to the cross-sections of H,r α . A suitable set of null normals is then Then since (by the definition of a dynamical horizon) Θ (l) = 0 and Θ (n) < 0, it follows that the extrinsic curvature scalar of S is where D α is the covariant derivative operator on H. This shows that the area of the crosssections of a spacelike dynamical horizon increases alongr α . We emphasize that this does not necessarily imply that the area increases in time. In spherical symmetry, we explicitly show below (section B.2) that the outward vector points in the future when the area increases in time and in the past when the area decreases in time. For a timelike dynamical horizon, the roles ofτ α andr α are interchanged. In this case,r α is no longer tangential to H, and is instead the unit spacelike vector normal to H. Additionally,τ α is instead the unit timelike vector tangent to H and orthogonal to the cross-sections of H. The area law then becomes i.e. the area of a timelike dynamical horizons decreases alongτ α . Note this law does not rely on any energy conditions, such as the NCC. Finally, we note that the area law defined in [80,81] only applies to dynamical horizons and timelike membranes. The definition does not include marginally anti-trapped tubes, which are often present in cosmological settings, or marginally trapped tubes which may not have a definite signature at a given time. To remedy this, Bousso and Engelhardt [76,77] formulated and proved a new area theorem applicable to an entire hypersurface H of indefinite 10 More recently this surface has also been called a timelike dynamical horizon; see appendix B of [81]. signature. The area theorem is based on a few technical assumptions but should be applicable to most hypersurfaces foliated by marginally trapped or anti-trapped surfaces S, called leaves. In this context, marginally (anti-)trapped tubes are referred to as future (or past) holographic screens. More precisely, the Bousso-Engelhardt area theorem is: the area of the leaves of any future (past) holographic screen, H, increases monotonically along H. The direction of increase along a future (past) holographic screen is the past (future) on timelike portions of H or exterior on spacelike portions of H. Thus H only evolves into the past (future) and/or exterior of each leaf.

B.2 Dynamical horizons and timelike membranes in spherical symmetry
As most of our simulations are performed in an essentially spherically symmetric spacetime, here we consider the properties of dynamical horizons for these spacetimes in more detail. The main purpose of this section is to illustrate how the area law for dynamical horizons [80,81] reduces to an essentially tautological statement about the dynamics of the horizon area.
We use r to denote the areal radius, and we will work with a gauge such that r is also a coordinate of the spacetime, that is we will consider a metric of the form ds 2 = α ab dx a dx b + r 2 dϑ 2 + sin 2 θdϕ 2 , (B.7) where α ab is a two-dimensional metric that is function of (t, r) (here t is the timelike coordinate). We recall that in spherical symmetry the expansion for a null vector v µ is [105] The last expression follows from our imposing a gauge such that the areal radius is also a coordinate of our spacetime. We consider the level sets of a function where we have defined Nτ ≡ √ −∇ α F ∇ α F . We next find the unit spacelike vector orthogonal toτ α ,r αr α = 1,r ατ α = 0. We writer α aŝ where Nr is the normalization. Defining the null vectors according to (B.4), a surface r(t) is trapped if Θ (l) = 0, Θ (n) < 0, and it is anti-trapped if Θ (l) > 0, Θ (n) = 0. The area law for dynamical horizons states that the area of the dynamical horizon must increase in the direction ofr α as we evolve alongr α [80,81]. From the form ofr α , we see that this reduces to: ifṙ > 0, then the dynamical horizon area increases in the direction of increasing time, and ifṙ < 0, then the dynamical horizon areas increases in the direction of decreasing time.

JCAP09(2022)063
Case 2: the level sets of F are timelike. Analogous to the case when the level set is spacelike, we define a unit spacelike vector orthogonal to the level set of F : Nr (ṙ, −1, 0, 0) , (B.12) and a unit timelike vector orthogonal tor α , where Nt is the normalization. Again, a surface r(t) is trapped if Θ (l) = 0, Θ (n) < 0, and it anti-trapped if Θ (l) > 0, Θ (n) = 0. The area law for timelike membranes states that the area of the timelike membrane must decrease in the direction ofτ α as we evolve alongτ α [80,81]. From the form ofτ α , we see that this statement then reduces to: ifṙ > 0, then the membrane area increases in the direction of increasing time, and ifṙ < 0, then the membrane area increases in the direction of decreasing time.

C The McVittie spacetime
Here we briefly review the McVittie spacetime [106] (see also [85,[107][108][109][110][111][112]), which is an analytic solution to the Einstein equations that describes a spherically symmetric black hole embedded in an asymptotically cosmological spacetime provided the cosmology asymptotes (in time: t → ∞) to a de-Sitter cosmology -for more discussion on this point, see [112]. 11 The two most salient properties of the McVittie spacetime are that the spacetime is spherically symmetric and satisfies the no-accretion condition, G r t = 0, which in turn implies that the stress-energy component T r t = 0. Thus, there is no radial flow of cosmic fluid in the McVittie solution (this assumption can be dropped for some generalizations of the McVittie spacetime [111]). We relax all of these assumptions in our numerical simulations, in addition to working in a set of coordinates that allows us to extend our spacetime past the black hole horizon, which to our knowledge has not yet been accomplished for the McVittie spacetime or its generalizations. While our numerical solutions differ in many of their properties from the McVittie spacetime, the McVittie spacetime serves as a useful analytic example to understand some of the properties of dynamical, apparent, and event horizons in spacetimes that have a black hole and an asymptotic cosmological expansion (see section B).
We consider only spatially flat McVittie solutions. The spacetime metric in isotropic coordinates is