Reflecting boundary conditions in numerical relativity as a model for black hole echoes

Recently, there has been much interest in black hole echoes, based on the idea that there may be some mechanism (e.g., from quantum gravity) that waves/fields falling into a black hole could partially reflect off of an interface before reaching the horizon. There does not seem to be a good understanding of how to properly model a reflecting surface in numerical relativity, as the vast majority of the literature avoids the implementation of artificial boundaries, or applies transmitting boundary conditions. Here, we present a framework for reflecting a scalar field in a fully dynamical spherically symmetric spacetime, and implement it numerically. We study the evolution of a wave packet in this situation and its numerical convergence, including when the location of a reflecting boundary is very close to the horizon of a black hole. This opens the door to model exotic near-horizon physics within full numerical relativity.


I. INTRODUCTION
The recent discovery of gravitational waves (GWs) has led to the new field of GW astronomy, opening novel windows into some of the deepest physical mysteries in our universe. Chief amongst these is the true nature of black hole horizons, which appear to have both temperature and entropy, if one applies laws of quantum mechanics to the fields in their vicinity [1,2]. Yet, according to Einstein's theory of General Relativity (GR), horizons comprise nothing but empty space. Even more problematic is when the thermal nature of black holes (i.e. Hawking radiation) leads to their evaporation, leaving us to wonder what happens to the information behind the horizon that cannot escape due to relativistic causality. This is known as the black hole information paradox [3,4]. Given that GW astronomy can now probe into the vicinity of black hole horizons, it is natural to ask whether it can give us an empirical window into resolving these theoretical mysteries. In particular, one may wonder whether quantum effects near the horizon would lead to echoes in observations that would have been absent in GR with classical horizons [5].
Motivated by this, there has been a recent interest in modelling black hole echoes phenomenologically [5][6][7][8][9][10][11][12][13]. The idea behind black hole echoes is that there may be some mechanism, resulting from quantum phenomena for example, by which wave packets falling into a black hole could partially reflect off of an interface before reaching the horizon, thus resulting in a potentially detectable effect. In the context of GW astronomy, this has opened a novel window to study exotic near-horizon physics phenomenologically, from current and future observations [5,10,11,13]. However, interpreting this data requires modelling strong gravitational systems in numerical relativity, where there does not currently seem to be a good * Corresponding Author: cdailey@pitp.ca understanding of how to properly implement a reflecting surface (but see [14] for an approximate method and [15] for an adjacent approach to boundaries). There has also been some recent skepticism that black hole echoes can be detected at all due to the formation of an apparent horizon before the reflection occurs [16], and having a proper way to model reflecting surfaces allows for rigorous testing beyond order-of-magnitude estimates.
One way to implement a reflecting surface is to have an artificial boundary in the domain of the numerical problem where the surface is to be located, and impose a set of boundary conditions (BCs) there. This requires an understanding of the full initial boundary value problem (IBVP) for a particular framework of numerical relativity, which the majority of studies tend to avoid by imposing no boundaries or by having boundaries out of causal contact with the initial conditions. This work contains some examples of the full spherically symmetric IBVP in numerical relativity with a scalar field. Of course there are no GWs in spherical symmetry, but see Appendix A for how this can be done analagously in cylindrical symmetry.
The generic way of treating the BCs in hyperbolic systems of partial differential equations is to identify the characteristic modes and speeds of the system, and for each mode that has a characteristic speed leading it out of the domain, no BCs are applied, while all of the modes propagating into the domain require a boundary treatment. It is then advantageous to pick a formulation of Einstein's equations where the characteristic modes and speeds are simple and well-behaved. Here, the Einstein-Christoffel (EC) formulation is used due to its simple characteristic structure, symmetric hyperbolicity, and previously successful numerical implementation [17].
In Section II, we describe the EC system and the scalar field equations adapted from [17]. In Section III, we derive BCs based on conservation laws while Section IV dictates the initial conditions we consider. Section V describes our numerical implementation using summation arXiv:2301.05778v2 [gr-qc] 30 Sep 2023 by parts (SBP) operators and BCs implemented with simultaneous approximation terms (SATs). Finally, we present our results for several reflecting IBVPs in Section VI, discuss near-horizon boundaries, and discuss the evolution of the apparent horizon as a response to [16].

II. THEORETICAL FRAMEWORK FOR THE EINSTEIN-CHRISTOFFEL SYSTEM
In general, the physical metric g µν in spherical symmetry can be expressed with respect to the spherical coordinates x α = (t, r, θ, φ), with the line element: where α is the lapse, β r is the shift, and γ ij are elements of the 3-metric. The extrinsic curvature is defined as where £ β represents the Lie derivative with respect to the shift vector. In order to make the evolution equations first order in space, the functions are defined. Here, Γ k ij are the Christoffel symbols associated with the 3-metric. The EC framework leaves the densitized lapseα ≡ α/( √ γ rr γ θθ ) and the shift as freely specifiable gauge functions. It is understood that wherever α may appear in the evolution equations, it is set with α =α √ γ rr γ θθ . While we present a general formalism in what follows, we shall come back to our particular choice of gauge in Section (II C).

A. The Einstein-Christoffel vacuum equations
The evolution equations in the EC system in spherical symmetry were introduced in detail in [17], which we modify to the following form: ∂ t f rrr − β r ∂ r f rrr + α∂ r K rr = 12 There are a few differences in this form of the equations compared to those presented in [17]. First, we write the equations explicitly in terms of the operator D i , which acts on a 3-covector f i as sons for this will be given in Section V A). Also, instead of defining transverse elements like γ T = γ θθ /r 2 , we keep the evolution equations in a form that does not explicitly depend on the coordinates. The principal part of this system makes up the left hand side and source-like terms make up the right hand side, which suggests the simple characteristic structure of this formulation of Ein-stein's vacuum equations. We should note that there is only gauge dynamics possible in the above system until we add stress-energy terms, which we do next.

B. Scalar Field and Source Terms
Scalar fields of mass m are subject to the Klein-Gordon equation ∇ µ ∇ µ ϕ − m 2 ϕ = 0. This equation can be reduced to a first order system and keep the same characteristic structure as the EC system with the definition of two auxiliary variables: ψ r ≡ ∂ r ϕ and Π ≡ −(∂ t ϕ − β r ∂ r ϕ)/α. The evolution equations then become With the unit normal to the spatial foliation n µ , the stress-energy tensor T µν is decomposed into the energy density ρ = n µ n ν T µν , the momentum density S i = −γ iµ n ν T µν , the trace T = T µ µ , and the spatial stress S ij = γ iµ γ jν T µν . In spherical symmetry the nonzero elements for the Klein-Gordon field are explicitly These elements modify the right hand side Eqns. (6)(7)(8) in this framework as and this system is subject to a set of five constraints: Eqns. (21) and (22) are the Hamiltonian and Momentum constraints respectively, Eqns. (23) and (24) come from definition Eq. (3), and Eq. (25) was defined at the beginning of this section. The characteristic modes for the complete system of equations include a set of three modes that propagate along the timelike normal to the foliation with characteristic speed −β r : and a set of six modes that propagate along the light cone with characteristic speeds c ± ≡ −β r ± α/ √ γ rr :

C. Gauge Choices in the Bulk
In this work, we will consider one particular gauge choice, that is thatα and β r are time independent: Other types of gauge choices have been discussed in [17], but it is not clear to us that many of these choices do not change the characteristic structure of the evolution equations. If one considers a hyperbolic time evolution of the gauge variables, in general this will directly affect the characteristic structure, but it is not clear whether the characteristics are well-defined ifα and β r are prescribed with elliptic gauge conditions. Therefore, we restrict our analysis to our gauge choice Eq. (30), where it is clear that the characteristic structure is as presented in the previous section, and leave other gauge choices for future work.

III. BOUNDARY CONDITIONS
We now restrict the problem to a spatial domain where r ∈ [a, b]. One way to gain some intuition on BCs for the scalar wave equation is to study the case of a static background. In such a case we have a time-like Killing vector ξ ν and thus a conserved current J µ = ξ ν T µ ν . We then define the static energy and rewrite the conservation law ∇ µ J µ = 0 in terms of the characteristics The flux term on the right hand side suggests the BCs that can be applied to allow E s to enter/leave the domain: for some coefficients k a,b . While k a,b do not necessarily need to be constants, the values k a,b = ±1 keep E s exactly conserved in the domain for each boundary.
In this work, a conservative BC that reverses the sign of an incoming pulse is referred to as a Dirichlet type (k a,b = −1), while one that retains the sign is referred to as a Neumann type (k a,b = +1). The motivation for this naming scheme comes from the fact that these BCs reduce to the classical Dirichlet (∂ t ϕ = 0) and Neumann (∂ r ϕ = 0) BCs in Minkowski space where c ± = ±1. It is important to note that we are considering domains and gauge conditions where c + c − < 0 throughout so that both boundaries always require BCs.
Once the metric becomes time dependent, J µ is no longer conserved. Instead, the notion of what a reflection means needs to be rooted in a quantity that describes energy in a dynamical spacetime. One way of defining this is with the Misner-Sharp mass, which in spherical symmetry can be written as [17] M (r) = √ γ θθ To obtain the BCs on the incoming scalar modes U ± ϕ , we express the derivatives of the Misner-Sharp mass in terms of the state vector: These relationships can be derived by taking the appropriate derivative of Eq. (35) and substituting time derivatives with the evolution equations and spatial derivatives with the constraints. It is interesting that these relationships can be expressed solely in terms of the state vector, which enables us to impose BCs on the scalar characteristics to control M . To see this, let us integrate Eq. (37) to define the Misner-Sharp "energy" enclosed within our domain In the continuum, we clearly have We can then write a conservation law by taking a time derivative and substituting Eq. (36). Written here for a massless scalar (i.e., m = 0), we have The flux term on the right hand side suggests the general BCs that can be applied to allow E M to enter/leave the domain: where the coefficient values k a,b = ∓1 define the Dirichlet/Neumann type BCs on the scalar characteristics that keep E M exactly conserved at each boundary. Interestingly, in the case of ∂ t γ θθ = 0, these reduce to the static background BCs, Eqns. (33-34), which motivates their connection to the classical Dirichlet/Neumann BCs. The expansion of outgoing null geodesics is proportional to the value of U − θ , and thus the apparent horizon is the outermost surface where U − θ = 0. We call the surface where c + = 0 the characteristic horizon. Only when c + < 0 are BCs no longer required at r = a as both characteristics are then leaving the domain. If c + c − < 0 and −U − θ U + θ < 0, i.e. the boundary r = a is between the apparent horizon and the characteristic horizon, the BCs at r = a are not well-defined unless k a = 0. We will not consider any simulations where such cases are present in this work.
The BCs on the incoming angular characteristics U ± θ are given by Eqn. (35): where the values of M (a) and M (b) at each time step are given by integrating their time derivatives with the scalar BCs applied We thus see that, in the case of Dirichlet/Neumann type BCs, the Misner-Sharp mass at both boundaries remains constant.
As long as β r > 0, the boundary r = b has three additional incoming modes U 0 r , U 0 θ , and U 0 ϕ . As suggested in [18], these can all be fixed using the constraints (23-25) by replacing the r derivatives in their evolution equations: all evaluated at r = b.
While fixing the masses M (a) and M (b), and scalar characteristics (Eqs. 40-41), fixes all the physical degrees of freedom at the boundaries, the incoming radial characteristics U ± r are left arbitrary, and thus should be connected to the residual gauge freedoms. Although we have specified the gauge dynamics in the bulk by choosingα and β r , additional gauge freedom can propagate into the domain from outside, and thus needs to be specified by fixing the incoming radial modes.
We shall consider a few choices for the BCs of U ± r . First, one could impose ∂ t U ± r = 0 at the corresponding boundaries: which ensures there are no incoming radial characteristics into the domain. Next, we could impose ∂ t K rr = 0: and finally we could impose ∂ t f rrr = 0: It is not clear if different choices may be more advantageous, but we find that for long term stability, the condition ∂ t K rr = 0 at both boundaries seems to work the best.

IV. INITIAL CONDITIONS
For the initial background, we consider a Schwarzschild black hole in Kerr-Schild coordinates, but with a rdependent mass function: The initial values of f rrr and f rθθ are formed from constraint Eqs. (23) and (24). The initial values of extrinsic curvature components are formed via Eqs. (4) and (5), with ∂ t γ θθ = 0 assumed, but ∂ t γ rr kept arbitrary. The yet unspecified functions M (r, 0) and ∂ t γ rr will be used to satisfy the constraints. For the scalar field, a spherical pulse of amplitude A, total width 2σ, and location µ is modeled as a section of a polynomial: and is set to zero outside of (µ − σ) ≤ r ≤ (µ + σ).
The degree of this polynomial assures that the first three derivatives are continuous. The initial condition on ψ r is obtained by taking the r derivative. For the initial condition on Π, a common choice is Π = 0, which describes a pulse that will break into two parts traveling at speeds c + and c − . Another choice is to specify that the pulse is initially traveling at speed c − , such that the initial time derivative is and is zero outside of (µ − σ) ≤ r ≤ (µ + σ). The initial value of Π is then completed with the definition Π = −(∂ t ϕ − β r ψ r )/α. Given initial conditions on the scalar field, the Hamiltonian and momentum constraints must be satisfied at the initial time. With the two functions M (r, 0) and ∂ t γ rr , one can show that these constraints are satisfied if These are easy conditions to satisfy, as M (r, 0) can be integrated with a standard solver starting at r = a with initial internal mass M 0 or from r = b with initial total mass M tot . The value of ∂ t γ rr can be replaced outright in the initial value of K rr .

V. NUMERICAL IMPLEMENTATION
In this section, we will introduce our numerical implementation of the IBVPs we are considering. We use Summation-By-Parts (SBP) derivative operators in an effort to preserve numerical stability by ensuring the energy of the system is bounded. In order to impose boundary conditions (BCs) in a stable fashion we use simultaneous approximation terms (SATs). Finally, we use numerical dissipation to prevent instability due to highfrequency noise.

A. Summation-By-Parts Derivative Operators
A popular way to maintain stability in numerical hyperbolic IBVPs is to use finite differencing operators that satisfy SBP, which is the discrete analog of integration by parts. Using SBP operators can be proven to be closely related to the well-posedness and the numerical stability of many conservative problems [19,20]. In order to introduce our covariant SBP scheme, we review some covariant integration properties. In a 3-dimensional domain V with coordinates x and the 2-dimensional boundary of that domain ∂V with coordinates ξ, the covariant Gauss's theorem states that for a vector field F i where the determinants of the metric on V and ∂V are γ and s respectively and n i is the unit normal to the boundary ∂V . If the vector field is represented by a scalar factor and a vector factor as F i = uv i , in spherical symmetry, this property can be written as where gence. This can be thought of as the covariant version of integration by parts in 1-dimension. Our goal will be to develop a numerical implementation of our IBVP that obeys a discrete analog of this property in an effort to maintain numerical stability. For our spatial domain r ∈ [a, b] we define a grid where all functions will be sampled with n points spaced by a distance h. An n × n finite differencing operator D paired with a symmetric and positive definite n × n norm operator Σ are said to satisfy SBP if [19] ΣD where the boundary operator B = diag(−1, 0, . . . , 0, 1). For this property to mimic the covariant property Eq. (61), we insert a matrix Γ with the values of √ γ injected along the diagonal: The operator ∇ ≡ Γ −1 DΓ approximates the covariant 3divergence of 3-vectors, the operator D approximates the covariant scalar gradient, and W ≡ ΣΓ is the covariant norm operator and approximates covariant integration. Then, in one spatial dimension, for a scalar grid function u and the one non-zero component of a 3-vector grid function v r , which directly mimics the continuous property Eq. (61).
In the case of the wave equation on a static background, one can write the energy conservation law in a directly analogous way to the covariant property Eq. (61) (see Appendix B), and thus when the system is discretized using the SBP operators, the system remains strictly stable in the sense that the energy is bounded [19,20]. It is no longer clear how to properly do this when the spacetime becomes dynamic. The state vector in the EC system is defined in terms of lower-index functions (f ijk ), which runs counter to the covariant SBP property that is defined in terms of an upper-index function. We have discovered one way that seems to be stable in at least all of the cases considered here, but a more elegant discretization paradigm may exist, perhaps one where the SBP properties can be made to mimic the dynamic conservation law Eq. (39), and the equations of motion can be written in terms of covariant divergences of upper-index functions. In the numerical scheme presented here, the evolution equations written in Section II are discretized by replacing scalar gradients with an SBP operator (i.e. ∂ r → D) and where present the operator D r → Γ −1 DΓ.
In this work, we use the D 6−3 SBP operator defined in [21] that minimizes the so-called average boundary truncation error. This operator is 6 th order accurate in the interior and 3 rd order accurate near the boundaries. It was reported in [21] that this operator achieves near 4 th order global convergence with a diagonal norm operator, which in principle will pair nicely with standard 4 th order Runge-Kutta time integration.

B. Application of Boundary Conditions
The BCs are applied using SATs. This method applies BCs in a weak fashion, adding an exponential decay term to the evolution equations at the boundary. For two arbitrary characteristic grid functions U ± (r, t) with characteristic speeds c ± , the application of the general BCs U + (a, t) = g a (t) and with strengths s a and s b which dictate the exponential decay scale [19]. One can show that for proper SBP energy conservation, in the context of the wave equation on a static background, the strengths must be equal to the magnitude of the incoming characteristic speeds [20]. At the left boundary we have s a → c + and at the right boundary we have s b → −c − . In the case of the wave equation around a static background black hole in Kerr-Schild coordinates, we have s b = 1 and 0 < s a < 1 where s a = 0 when a = 2M 0 as the boundary is then placed exactly on the horizon, where BCs are not to be imposed anyway since there are no longer any incoming characteristics into the domain at r = a.

C. Numerical Dissipation
To stabilize nonlinear evolution it is common to add numerical dissipation, ensuring that high frequency noise does not destabilize the evolution. We use the A 6 numerical dissipation operator that applies to the D 6−3 operator defined in [21], and it is added to the right hand side of the evolution equations as where ϵ > 0 is the amount of dissipation, which is usually chosen to be of order unity. The coefficients of both operators D 6−3 and A 6 can be found in [21].

VI. RESULTS
First, to demonstrate long term stability and energy conservation of this framework, an IBVP with two distinct reflection boundaries, one of Dirichlet type and one of Neumann type, will be demonstrated in several situations. Then we will investigate situations where the boundary r = a is brought very close to the horizon of the black hole. Finally, we discuss the movement of the apparent horizon and situations where the initial pulse can reflect without the apparent horizon moving into the domain.

A. Enclosed Reflecting Boundary Problems
Here we consider several enclosed reflecting boundary problems. These problems involve an initial pulse and BCs of either Dirichlet or Neumann type so that energy is conserved in the continuum limit. We consider three situations for this type of problem, first a Minkowski spherical wave equation, second a static background Schwarzschild black hole in Kerr-Schild coordinates, and third the fully dynamic problem around a black hole. In all of these examples, the initial conditions and evolutions are kept as similar as possible to aid in a direct comparison.
Here, we consider a spatial domain r ∈ [3,8]M 0 and a temporal domain t ∈ [0, 200]M 0 . As the mass M 0 is not relevant in Minkowski space, we simply take r ∈ [3,8] and t ∈ [0, 200] in this case. We use a Dirichlet type BC for r = a (i.e. k a = −1) and a Neumann type BC for r = b (i.e. k b = 1). For the problems with a black hole, the internal mass is set to M 0 = 1. The initial pulse follows the model Eq. (56) with parameters µ = 6M 0 , 2σ = M 0 , and A = 0.05/M 0 . We use initial condition Eq. (57) so that the pulse is initially traveling with speed c − = −1.
The time integration is done using 4 th order Runge-Kutta time stepping, with time step size h t = h/4. We also apply the numerical dissipation operator with strength ϵ = 1 in all cases so that the energy loss is directly comparable to the fully dynamic case.

Static Minkowski Background
where the ϕ (i) are the solutions obtained with a certain grid spacing h (i) and where the ℓ 2 norm is defined for a grid function f as We use h (i) = h (i−1) /2 as the resolution is doubled each run. Five runs were conducted starting at h (1) = 0.02M 0 . This leaves three subsets of resolutions to calculate the solution convergence order p order as the resolution is increased, deviating when a reflection occurs where the error is dominated by the boundary region accurate to 3 rd order. At the bottom we see the wave solution in the r-t plane, where it is apparent that a pulse propagating in the −r direction has a constant characteristic speed while one traveling in the +r direction appears to accelerate along curves. Figure 2 depicts the numerical solution to the IBVP with a static Schwarzschild background in Kerr-Schild coordinates. The same five resolutions are used to obtain the Low, Mid, and High convergence order plots, which shows once again consistent 4 th order convergence as the resolution is increased. Since the characteristic speeds here are c − = −1 and c + = (r − 2M 0 )/(r + 2M 0 ), pulses traveling in the −r direction follow straight lines in the r-t plane and stay the same width, but pulses traveling in the +r direction appear to accelerate along curves and have a characteristic speed that varies across the pulse, so the width increases as the pulse propagates in the +r direction.

Dynamic Schwarzschild Background
Here we demonstrate a fully dynamic IBVP, where the scalar pulse is coupled to gravity. Along with the solution convergence order, we also define the constraint convergence order where the norm of the constraints is defined as This somewhat simpler definition for the constraint convergence order is possible because we know the exact solution (i.e. C = C r = 0) whereas we don't necessarily know the exact solution for the scalar field ϕ.
We use total reflection BCs on the scalar field Eqns. The results are plotted in Figure 3. At the initial time slice, c + = [r − 2M (r)]/[r + 2M (r)], so the outgoing characteristic speed varies across the pulse as in the static background case, but also as M (r) changes across the pulse width, so there is an additional dispersion due to the coupling to gravity. In this particular gauge, ∂ t c ± = ±α∂ t γ θθ , so the characteristic speeds change proportional to the change in the areal radius, R = √ γ θθ . We notice that our choice of BCs causes the boundaries to move inward in R , making the reflected pulse at the inner boundary move more slowly as time progresses.
In the continuum limit, each of the above IBVPs conserves the energy Eq. (39). Figure 4 shows that the energy loss at t = 200M 0 in these discrete problems converges to zero at about 6 th order. Since the strength of the SAT boundary at r = a depends on the characteristic speed c + , it is applied more weakly in the Kerr-Schild and dynamic cases, leading to significantly more energy dissipation at the boundary. This is apparent in Figure 4 as the energy decreases about three orders of magnitude more in the Kerr-Schild case than the flat spherical case at t = 200M 0 , although still converges to zero at 6 th order. shows the solution convergence order that starts near 4 th order but changes dramatically as the solution becomes more dependent on boundary interactions and the solution errors move out of phase with each other. Plot B) shows the constraint convergence order slowly settles to about 4 th order as the resolution is increased. Plot C) shows the scalar field solution shown as a function of the areal radius R = √ γ θθ over coordinate time t. It is apparent that over the course of the evolution, the inner boundary falls toward the horizon by ∼ 0.5M0 and the outer boundary by ∼ 2.7M0. The pulse also appears to disperse much more than in Fig. 2, due to the variation of c± as the Misner-Sharp mass M (r) changes across the pulse width. Plot D) shows the Misner-Sharp mass in the same R-t plane, demonstrating the coupling to gravity.

B. Reflections close to the horizon
Black hole echoes are usually discussed in the context of quantum gravity, suggesting that quantum phenomena may result in a reflecting surface. If this is the case, one might expect such a surface to be a distance on the order of the Planck scale away from the horizon (e.g., [5,6,22]). To simulate this, one might be interested in applying reflecting BCs very close to the horizon. Although this problem is well-posed in principle, in practice it can be difficult to achieve numerically.
In the context of the static black hole background, if a pulse of width w is incident on a reflecting surface, the reflected pulse will have width ∼ w|c − /c + |, so one should be sure that h is at least several times smaller than this. There is also usually a significant amount of high frequency noise that propagates at high speed off of the boundary r = a, which may need to be controlled with numerical dissipation. When c + is very small at r = a, the SATs become very weak and thus lose much more energy, so in order for the discrete problem to conserve energy to a certain tolerance, much more resolution might be required than the naive estimate h < w|c − /c + |.

C. Apparent Horizons
Finally, as a response to [16], we can study when an apparent horizon moves into the domain of the simulation, making BC Eq. (40) ill-defined. In [16], Guo and Mathur argue that for merging equal mass black holes, one should never expect echoes of GWs due to the formation of an apparent horizon that envelopes the waves. They make the assumption that the waves have a wavelength on the same order of the black hole masses, but any localized wavepacket should contain shorter, as well as longer, wavelengths. If a black hole absorbs a sufficiently localized infalling wavepacket with a total mass M ϕ = M tot −M 0 , then the apparent horizon should move from R = 2M 0 to R = 2M tot as t → ∞. This suggests  6. Plot of the maximum pulse mass M ϕ of an initial pulse before the apparent horizon moves into the domain for various widths of the initial pulse. Each case is asymptotic to M ϕ = 0.05M0 as the width goes to zero (σ → 0), which corresponds to an infinitely sharp pulse with a = 2Mtot. that a reflecting surface should be located at R > 2M tot to ensure that the wavepacket does not become enveloped by a trapped surface. We demonstrate, however, that this estimate can be relaxed for sufficiently wide pulses, and thus the frequency content of a wavepacket plays an important role in whether a reflection can occur.
We consider a domain r ∈ [2.1, 203]M 0 with the initial pulse location µ = 103M 0 and vary its width from σ = 10M 0 to σ = 100M 0 . We consider two situations, one where the boundary r = a has Neumann type reflective condition k a = +1, and one where it has transmitting (or absorbing) condition k a = 0. If the mass of the initial pulse M ϕ is too large, the apparent horizon will move into the domain from its initial position of r = 2M 0 which is initially outside of the domain. We then maximize the mass of the initial pulse by adjusting the amplitude A to just before this occurs, and plot the results of this maximum mass for the varying initial widths. The transmitting case k a = 0 demonstrates reflection just from the Schwarzchild potential, a feature that occurs regardless of BCs we apply. Figure 6 shows that the case k a = +1 (reflective BC) allows a significantly higher maximum mass M ϕ than the transmitting case k a = 0 (absorptive BC), which suggests that there is a possibility for low frequency constituents of the incident scalar pulse to reflect near the horizon in a way discernible from reflection due just to the Schwarzchild potential. Indeed, low frequency (or wide) pulses, comparable to the Hawking frequency 1/(8πM 0 ) are also expected to have the largest reflection (from the horizon) based on quantum mechanical arguments [6].
It is important to note that as we have observed previously, our gauge conditions tend to move the boundaries inward, and so the apparent horizon can move into the domain either because it is growing in areal radius, or due to the boundary r = a shrinking in areal radius. The orange curve in Figure 6 should then be thought of as a lower bound that could be higher if the boundary was fixed in areal radius (through a different choice of gauge/BCs).

VII. CONCLUSION AND FUTURE PROSPECTS
In this work, we developed and demonstrated a full initial boundary value problem in numerical relativity in spherical symmetry. In particular, we derived boundary conditions (BCs) based on conservation laws involving a quasi-local mass-energy measure, which is the main result of this work. These BCs are general, allowing one to dictate how much mass-energy enters or leaves the domain. We demonstrated how this framework can model black hole echoes, which has mostly been discussed previously either within linear perturbation theory that ignored backreaction or order of magnitude estimates.
However, we only considered a gauge condition whereα and β r were both time independent. The BC framework presented in Section III can be applied in the case of a different gauge choice, although the characteristic structure of the evolution equations may change in a non-trivial way. To properly model a black hole echo event, one may want some way of controlling the position of the boundary r = a relative to the apparent horizon. For example, one might require that the reflection surface remain a constant proper distance from the apparent horizon, so that a pulse can be reflected without the concern that the apparent horizon would move into the domain. Similar conditions are considered in [23,24], where gauge conditions that ensure the inner boundary of the domain stays just inside the apparent horizon are used, which could be adapted to our case where we want the boundary just outside the apparent horizon. We leave a better understanding of how to control the boundary r = a to future work, that may be prescribed by yet-unknown quantum gravitational processes.
This BC framework may also be applied in a full 3D simulation. We considered here in spherical symmetry the Misner-Sharp mass as a way to measure the quasilocal mass-energy of the spacetime, but this is not a unique choice. In 3 dimensions, one would need to choose a quasi-local mass-energy measure (e.g. those found in [25]) and derive BCs such that it can be conserved in the same sense as what we have done here. Future work may include a 3D generalization of this framework.

ACKNOWLEDGMENTS
We are thankful for valuable discussions with David Brown, Luis Lehner, and the Strong Gravity group at the Perimeter Institute. Numerical implementation for this project was accomplished with the Julia programming language [26] and the packages [27]. This research was funded thanks in part to the Canada First Research Excellence Fund through the Arthur B. McDonald Canadian Astroparticle Physics Research Institute, the Natural Sciences and Engineering Research Council of Canada, and the Perimeter Institute. Research at Perimeter Institute is supported in part by the Government of Canada through the Department of Innovation, Science and Economic Development and by the Province of Ontario through the Ministry of Colleges and Universities.
Dirichlet/Neumann conditions ∂ t ψ = 0 and ∂ ρ ψ = 0. We can then define a reflecting IBVP as we did in spherical symmetry. The BCs are simple due to the gauge fixed assumption of the line element considered here, but this may not be the case if a more general line element where the "perimeter radius" ρ is allowed to be dynamic, as it was in our dynamic examples in spherical symmetry.